From e9f443528fae75ccb971eba313e244e7dde85062 Mon Sep 17 00:00:00 2001 From: ndeak Date: Mon, 18 Sep 2023 13:36:55 -0600 Subject: [PATCH] Added cvode functionality to SUNDIALS integrator (#3436) ## Summary ## Additional background ## Checklist The proposed changes: - [ ] fix a bug or incorrect behavior in AMReX - [x] add new capabilities to AMReX - [ ] changes answers in the test suite to more than roundoff level - [ ] are likely to significantly affect the results of downstream AMReX users - [ ] include documentation in the code and/or rst files, if appropriate --------- Co-authored-by: Nicholas Deak Co-authored-by: Ann Almgren --- .../SUNDIALS/AMReX_SundialsIntegrator.H | 133 +++++++++++++++++- 1 file changed, 131 insertions(+), 2 deletions(-) diff --git a/Src/Extern/SUNDIALS/AMReX_SundialsIntegrator.H b/Src/Extern/SUNDIALS/AMReX_SundialsIntegrator.H index 305da3af4ee..9f8cccedaa7 100644 --- a/Src/Extern/SUNDIALS/AMReX_SundialsIntegrator.H +++ b/Src/Extern/SUNDIALS/AMReX_SundialsIntegrator.H @@ -8,10 +8,12 @@ #include /* prototypes for ERKStep fcts., consts */ #include /* prototypes for ARKStep fcts., consts */ #include /* prototypes for MRIStep fcts., consts */ +#include /* access to CVODE solver */ #include /* manyvector N_Vector types, fcts. etc */ #include /* MultiFab N_Vector types, fcts., macros */ #include /* MultiFab N_Vector types, fcts., macros */ #include /* access to SPGMR SUNLinearSolver */ +#include /* access to SPGMR SUNLinearSolver */ #include /* access to FixedPoint SUNNonlinearSolver */ #include /* defs. of realtype, sunindextype, etc */ @@ -70,11 +72,13 @@ private: bool use_erk_strategy; bool use_mri_strategy; bool use_mri_strategy_test; + bool use_cvode_strategy; bool use_implicit_inner; SUNNonlinearSolver NLS; /* empty nonlinear solver object */ SUNLinearSolver LS; /* empty linear solver object */ void *arkode_mem; /* empty ARKode memory structure */ + void *cvode_mem; /* empty CVODE memory structure */ SUNNonlinearSolver NLSf; /* empty nonlinear solver object */ SUNLinearSolver LSf; /* empty linear solver object */ void *inner_mem; /* empty ARKode memory structure */ @@ -101,6 +105,7 @@ private: use_erk_strategy=false; use_mri_strategy=false; use_mri_strategy_test=false; + use_cvode_strategy=false; amrex::ParmParse pp("integration.sundials"); @@ -124,6 +129,10 @@ private: use_mri_strategy=true; use_mri_strategy_test=true; } + else if (theStrategy == "CVODE") + { + use_cvode_strategy=true; + } else { std::string msg("Unknown strategy: "); @@ -146,6 +155,7 @@ private: NLS = nullptr; /* empty nonlinear solver object */ LS = nullptr; /* empty linear solver object */ arkode_mem = nullptr; /* empty ARKode memory structure */ + cvode_mem = nullptr; /* empty CVODE memory structure */ NLSf = nullptr; /* empty nonlinear solver object */ LSf = nullptr; /* empty linear solver object */ inner_mem = nullptr; /* empty ARKode memory structure */ @@ -188,8 +198,10 @@ public: return advance_mri(S_old, S_new, time, time_step); } else if (use_erk_strategy) { return advance_erk(S_old, S_new, time, time_step); - } else { - Error("SUNDIALS integrator backend not specified (ERK or MRI)."); + } else if (use_cvode_strategy) { + return advance_cvode(S_old, S_new, time, time_step); + }else { + Error("SUNDIALS integrator backend not specified (ERK, MRI, or CVODE)."); } return 0; @@ -643,6 +655,123 @@ public: return timestep; } + amrex::Real advance_cvode (T& S_old, T& S_new, amrex::Real time, const amrex::Real time_step) + { + t = time; + tout = time+time_step; + hfixed = time_step; + timestep = time_step; + + // We use S_new as our working space, so first copy S_old to S_new + IntegratorOps::Copy(S_new, S_old); + + // Create an N_Vector wrapper for the solution MultiFab + auto get_length = [&](int index) -> sunindextype { + auto* p_mf = &S_new[index]; + return p_mf->nComp() * (p_mf->boxArray()).numPts(); + }; + + /* Create manyvector for solution using S_new */ + NVar = S_new.size(); // NOTE: expects S_new to be a Vector + nv_many_arr = new N_Vector[NVar]; // vector array composed of cons, xmom, ymom, zmom component vectors */ + + for (int i = 0; i < NVar; ++i) { + sunindextype length = get_length(i); + N_Vector nvi = amrex::sundials::N_VMake_MultiFab(length, &S_new[i]); + nv_many_arr[i] = nvi; + } + + nv_S = N_VNew_ManyVector(NVar, nv_many_arr, sunctx); + nv_stage_data = N_VClone(nv_S); + + /* Create a temporary storage space for MRI */ + Vector > temp_storage; + IntegratorOps::CreateLike(temp_storage, S_old); + T& state_store = *temp_storage.back(); + + SundialsUserData udata; + + /* Begin Section: SUNDIALS FUNCTION HOOKS */ + /* f routine to compute the ODE RHS function f(t,y). */ + udata.f = [&](realtype rhs_time, N_Vector y_data, N_Vector y_rhs, void * /* user_data */) -> int { + amrex::Vector S_data; + amrex::Vector S_rhs; + + const int num_vecs = N_VGetNumSubvectors_ManyVector(y_data); + S_data.resize(num_vecs); + S_rhs.resize(num_vecs); + + for(int i=0; inComp()); + S_rhs.at(i)=amrex::MultiFab(*amrex::sundials::getMFptr(N_VGetSubvector_ManyVector(y_rhs, i)),amrex::make_alias,0,amrex::sundials::getMFptr(N_VGetSubvector_ManyVector(y_rhs, i))->nComp()); + } + + BaseT::post_update(S_data, rhs_time); + BaseT::rhs(S_rhs, S_data, rhs_time); + + return 0; + }; + + udata.ProcessStage = [&](realtype rhs_time, N_Vector y_data, void * /* user_data */) -> int { + amrex::Vector S_data; + + const int num_vecs = N_VGetNumSubvectors_ManyVector(y_data); + S_data.resize(num_vecs); + + for (int i=0; inComp()); + } + + BaseT::post_update(S_data, rhs_time); + + return 0; + }; + /* End Section: SUNDIALS FUNCTION HOOKS */ + + /* Set up CVODE BDF solver */ + cvode_mem = CVodeCreate(CV_BDF, sunctx); + CVodeSetUserData(cvode_mem, &udata); + CVodeInit(cvode_mem, SundialsUserFun::f, time, nv_S); + CVodeSStolerances(cvode_mem, reltol, abstol); + CVodeSetMaxNumSteps(cvode_mem, 100000); + + for(int i=0; i= 0); + + // Copy the result stored in nv_S to state_new + for(int i=0; i /* Map */) override {}