diff --git a/.gitignore b/.gitignore index 3e9185713d..763405865e 100644 --- a/.gitignore +++ b/.gitignore @@ -36,6 +36,8 @@ compile_commands.json *.orig *.nvvp *.ptp-sync* +*.project +*.cproject *~ \#* diff --git a/CHANGELOG.md b/CHANGELOG.md index 2aff004841..40eb177e37 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,11 @@ Volta GPUs while the automatically selected value will vary across compilers and compiler versions. As such, users are encouraged to override this value with the architecture for their system. +Added a time-stepping module to ARKODE for low storage Runge--Kutta methods, LSRKStep. +This currently supports five explicit low-storage methods: the second-order Runge--Kutta--Chebyshev +and Runge--Kutta--Legendre methods, and the second- through fourth-order optimal strong stability +preserving Runge--Kutta methods. All methods include embeddings for temporal adaptivity. + The Trilinos Teptra NVector interface has been updated to utilize CMake imported targets added in Trilinos 14 to improve support for different Kokkos backends with Trilinos. As such, Trilinos 14 or newer is required and the diff --git a/benchmarks/diffusion_2D/main_arkode.cpp b/benchmarks/diffusion_2D/main_arkode.cpp index 39c6274a6e..a29a3bcf6a 100644 --- a/benchmarks/diffusion_2D/main_arkode.cpp +++ b/benchmarks/diffusion_2D/main_arkode.cpp @@ -15,7 +15,10 @@ * ---------------------------------------------------------------------------*/ #include "arkode/arkode_arkstep.h" +#include "arkode/arkode_lsrkstep.h" #include "diffusion_2D.hpp" +#include "sunadaptcontroller/sunadaptcontroller_imexgus.h" +#include "sunadaptcontroller/sunadaptcontroller_soderlind.h" struct UserOptions { @@ -28,6 +31,8 @@ struct UserOptions int maxsteps = 0; // max steps between outputs int onestep = 0; // one step mode, number of steps bool linear = true; // linearly implicit RHS + bool implicit = true; // implicit (ARKStep) vs explicit STS (LSRKStep) + ARKODE_LSRKMethodType lsrkmethod = ARKODE_LSRK_RKC_2; // LSRK method type // Linear solver and preconditioner settings std::string ls = "cg"; // linear solver to use @@ -43,6 +48,14 @@ struct UserOptions void print(); }; +// ----------------------------------------------------------------------------- +// LSRKStep-specific dominant eigenvalue function prototype +// ----------------------------------------------------------------------------- + +static int dom_eig(sunrealtype t, N_Vector y, N_Vector fn, sunrealtype* lambdaR, + sunrealtype* lambdaI, void* user_data, N_Vector temp1, + N_Vector temp2, N_Vector temp3); + // ----------------------------------------------------------------------------- // Main Program // ----------------------------------------------------------------------------- @@ -58,8 +71,8 @@ int main(int argc, char* argv[]) // Create SUNDIALS context MPI_Comm comm = MPI_COMM_WORLD; - SUNContext ctx = NULL; - SUNProfiler prof = NULL; + SUNContext ctx = nullptr; + SUNProfiler prof = nullptr; flag = SUNContext_Create(comm, &ctx); if (check_flag(&flag, "SUNContextCreate", 1)) { return 1; } @@ -109,6 +122,17 @@ int main(int argc, char* argv[]) return 1; } + // Return with error on unsupported LSRK method type + if (!uopts.implicit) + { + if ((uopts.lsrkmethod != ARKODE_LSRK_RKC_2) && + (uopts.lsrkmethod != ARKODE_LSRK_RKL_2)) + { + cerr << "ERROR: illegal lsrkmethod" << endl; + return 1; + } + } + // ----------------------------- // Setup parallel decomposition // ----------------------------- @@ -152,12 +176,8 @@ int main(int argc, char* argv[]) if (check_flag((void*)(uout.error), "N_VClone", 0)) { return 1; } } - // --------------------- - // Create linear solver - // --------------------- - - // Create linear solver - SUNLinearSolver LS = NULL; + // Set up implicit solver, if applicable + SUNLinearSolver LS = nullptr; SUNMatrix A = nullptr; #if defined(USE_SUPERLU_DIST) // SuperLU-DIST objects @@ -172,77 +192,96 @@ int main(int argc, char* argv[]) sunindextype* A_col_idxs = nullptr; sunindextype* A_row_ptrs = nullptr; #endif - - int prectype = (uopts.preconditioning) ? SUN_PREC_RIGHT : SUN_PREC_NONE; - - if (uopts.ls == "cg") - { - LS = SUNLinSol_PCG(u, prectype, uopts.liniters, ctx); - if (check_flag((void*)LS, "SUNLinSol_PCG", 0)) { return 1; } - } - else if (uopts.ls == "gmres") + if (uopts.implicit) { - LS = SUNLinSol_SPGMR(u, prectype, uopts.liniters, ctx); - if (check_flag((void*)LS, "SUNLinSol_SPGMR", 0)) { return 1; } - } - else - { -#if defined(USE_SUPERLU_DIST) - // Initialize SuperLU-DIST grid - superlu_gridinit(udata.comm_c, udata.npx, udata.npy, &grid); - - // Create arrays for CSR matrix: data, column indices, and row pointers - sunindextype nnz_loc = 5 * udata.nodes_loc; - - A_data = (sunrealtype*)malloc(nnz_loc * sizeof(sunrealtype)); - if (check_flag((void*)A_data, "malloc Adata", 0)) return 1; - - A_col_idxs = (sunindextype*)malloc(nnz_loc * sizeof(sunindextype)); - if (check_flag((void*)A_col_idxs, "malloc Acolind", 0)) return 1; + // --------------------- + // Create linear solver + // --------------------- - A_row_ptrs = - (sunindextype*)malloc((udata.nodes_loc + 1) * sizeof(sunindextype)); - if (check_flag((void*)A_row_ptrs, "malloc Arowptr", 0)) return 1; + // Create linear solver - // Create and initialize SuperLU_DIST structures - dCreate_CompRowLoc_Matrix_dist(&A_super, udata.nodes, udata.nodes, nnz_loc, - udata.nodes_loc, 0, A_data, A_col_idxs, - A_row_ptrs, SLU_NR_loc, SLU_D, SLU_GE); - dScalePermstructInit(udata.nodes, udata.nodes, &A_scaleperm); - dLUstructInit(udata.nodes, &A_lu); - PStatInit(&A_stat); - set_default_options_dist(&A_opts); - A_opts.PrintStat = NO; + int prectype = (uopts.preconditioning) ? SUN_PREC_RIGHT : SUN_PREC_NONE; - // SUNDIALS structures - A = SUNMatrix_SLUNRloc(&A_super, &grid, ctx); - if (check_flag((void*)A, "SUNMatrix_SLUNRloc", 0)) return 1; - - LS = SUNLinSol_SuperLUDIST(u, A, &grid, &A_lu, &A_scaleperm, &A_solve, - &A_stat, &A_opts, ctx); - if (check_flag((void*)LS, "SUNLinSol_SuperLUDIST", 0)) return 1; - - uopts.preconditioning = false; + if (uopts.ls == "cg") + { + LS = SUNLinSol_PCG(u, prectype, uopts.liniters, ctx); + if (check_flag((void*)LS, "SUNLinSol_PCG", 0)) { return 1; } + } + else if (uopts.ls == "gmres") + { + LS = SUNLinSol_SPGMR(u, prectype, uopts.liniters, ctx); + if (check_flag((void*)LS, "SUNLinSol_SPGMR", 0)) { return 1; } + } + else + { +#if defined(USE_SUPERLU_DIST) + // Initialize SuperLU-DIST grid + superlu_gridinit(udata.comm_c, udata.npx, udata.npy, &grid); + + // Create arrays for CSR matrix: data, column indices, and row pointers + sunindextype nnz_loc = 5 * udata.nodes_loc; + + A_data = (sunrealtype*)malloc(nnz_loc * sizeof(sunrealtype)); + if (check_flag((void*)A_data, "malloc Adata", 0)) return 1; + + A_col_idxs = (sunindextype*)malloc(nnz_loc * sizeof(sunindextype)); + if (check_flag((void*)A_col_idxs, "malloc Acolind", 0)) return 1; + + A_row_ptrs = + (sunindextype*)malloc((udata.nodes_loc + 1) * sizeof(sunindextype)); + if (check_flag((void*)A_row_ptrs, "malloc Arowptr", 0)) return 1; + + // Create and initialize SuperLU_DIST structures + dCreate_CompRowLoc_Matrix_dist(&A_super, udata.nodes, udata.nodes, + nnz_loc, udata.nodes_loc, 0, A_data, + A_col_idxs, A_row_ptrs, SLU_NR_loc, + SLU_D, SLU_GE); + dScalePermstructInit(udata.nodes, udata.nodes, &A_scaleperm); + dLUstructInit(udata.nodes, &A_lu); + PStatInit(&A_stat); + set_default_options_dist(&A_opts); + A_opts.PrintStat = NO; + + // SUNDIALS structures + A = SUNMatrix_SLUNRloc(&A_super, &grid, ctx); + if (check_flag((void*)A, "SUNMatrix_SLUNRloc", 0)) return 1; + + LS = SUNLinSol_SuperLUDIST(u, A, &grid, &A_lu, &A_scaleperm, &A_solve, + &A_stat, &A_opts, ctx); + if (check_flag((void*)LS, "SUNLinSol_SuperLUDIST", 0)) return 1; + + uopts.preconditioning = false; #else - std::cerr << "ERROR: Benchmark was not built with SuperLU_DIST enabled\n"; - return 1; + std::cerr + << "ERROR: Benchmark was not built with SuperLU_DIST enabled\n"; + return 1; #endif - } + } - // Allocate preconditioner workspace - if (uopts.preconditioning) - { - udata.diag = N_VClone(u); - if (check_flag((void*)(udata.diag), "N_VClone", 0)) { return 1; } + // Allocate preconditioner workspace + if (uopts.preconditioning) + { + udata.diag = N_VClone(u); + if (check_flag((void*)(udata.diag), "N_VClone", 0)) { return 1; } + } } - // -------------- - // Setup ARKStep - // -------------- + // ---------------------- + // Setup ARKStep/LSRKStep + // ---------------------- // Create integrator - void* arkode_mem = ARKStepCreate(NULL, diffusion, ZERO, u, ctx); - if (check_flag((void*)arkode_mem, "ARKStepCreate", 0)) { return 1; } + void* arkode_mem = nullptr; + if (uopts.implicit) + { + arkode_mem = ARKStepCreate(nullptr, diffusion, ZERO, u, ctx); + if (check_flag((void*)arkode_mem, "ARKStepCreate", 0)) { return 1; } + } + else + { + arkode_mem = LSRKStepCreateSTS(diffusion, ZERO, u, ctx); + if (check_flag((void*)arkode_mem, "LSRKStepCreateSTS", 0)) { return 1; } + } // Specify tolerances flag = ARKodeSStolerances(arkode_mem, uopts.rtol, uopts.atol); @@ -252,38 +291,60 @@ int main(int argc, char* argv[]) flag = ARKodeSetUserData(arkode_mem, (void*)&udata); if (check_flag(&flag, "ARKodeSetUserData", 1)) { return 1; } - // Attach linear solver - flag = ARKodeSetLinearSolver(arkode_mem, LS, A); - if (check_flag(&flag, "ARKodeSetLinearSolver", 1)) { return 1; } + // Configure implicit solver + if (uopts.implicit) + { + // Attach linear solver + flag = ARKodeSetLinearSolver(arkode_mem, LS, A); + if (check_flag(&flag, "ARKodeSetLinearSolver", 1)) { return 1; } #if defined(USE_SUPERLU_DIST) - if (uopts.ls == "sludist") - { - ARKodeSetJacFn(arkode_mem, diffusion_jac); - if (check_flag(&flag, "ARKodeSetJacFn", 1)) return 1; - } + if (uopts.ls == "sludist") + { + ARKodeSetJacFn(arkode_mem, diffusion_jac); + if (check_flag(&flag, "ARKodeSetJacFn", 1)) return 1; + } #endif - if (uopts.preconditioning) - { - // Attach preconditioner - flag = ARKodeSetPreconditioner(arkode_mem, PSetup, PSolve); - if (check_flag(&flag, "ARKodeSetPreconditioner", 1)) { return 1; } + if (uopts.preconditioning) + { + // Attach preconditioner + flag = ARKodeSetPreconditioner(arkode_mem, PSetup, PSolve); + if (check_flag(&flag, "ARKodeSetPreconditioner", 1)) { return 1; } - // Set linear solver setup frequency (update preconditioner) - flag = ARKodeSetLSetupFrequency(arkode_mem, uopts.msbp); - if (check_flag(&flag, "ARKodeSetLSetupFrequency", 1)) { return 1; } - } + // Set linear solver setup frequency (update preconditioner) + flag = ARKodeSetLSetupFrequency(arkode_mem, uopts.msbp); + if (check_flag(&flag, "ARKodeSetLSetupFrequency", 1)) { return 1; } + } + + // Set linear solver tolerance factor + flag = ARKodeSetEpsLin(arkode_mem, uopts.epslin); + if (check_flag(&flag, "ARKodeSetEpsLin", 1)) { return 1; } + + // Select method order + flag = ARKodeSetOrder(arkode_mem, uopts.order); + if (check_flag(&flag, "ARKodeSetOrder", 1)) { return 1; } - // Set linear solver tolerance factor - flag = ARKodeSetEpsLin(arkode_mem, uopts.epslin); - if (check_flag(&flag, "ARKodeSetEpsLin", 1)) { return 1; } + // Specify linearly implicit non-time-dependent RHS + if (uopts.linear) + { + flag = ARKodeSetLinear(arkode_mem, 0); + if (check_flag(&flag, "ARKodeSetLinear", 1)) { return 1; } + } + } + else // Configure explicit STS solver + { + // Select LSRK method + flag = LSRKStepSetMethod(arkode_mem, uopts.lsrkmethod); + if (check_flag(&flag, "LSRKStepSetMethod", 1)) { return 1; } - // Select method order - flag = ARKodeSetOrder(arkode_mem, uopts.order); - if (check_flag(&flag, "ARKodeSetOrder", 1)) { return 1; } + // Provide dominant eigenvalue function + flag = LSRKStepSetDomEigFn(arkode_mem, dom_eig); + if (check_flag(&flag, "LSRKStepSetDomEigFn", 1)) { return 1; } + } // Set fixed step size or adaptivity method + SUNAdaptController C = nullptr; if (uopts.hfixed > ZERO) { flag = ARKodeSetFixedStep(arkode_mem, uopts.hfixed); @@ -291,16 +352,17 @@ int main(int argc, char* argv[]) } else { - flag = ARKStepSetAdaptivityMethod(arkode_mem, uopts.controller, SUNTRUE, - SUNFALSE, NULL); - if (check_flag(&flag, "ARKStepSetAdaptivityMethod", 1)) { return 1; } - } - - // Specify linearly implicit non-time-dependent RHS - if (uopts.linear) - { - flag = ARKodeSetLinear(arkode_mem, 0); - if (check_flag(&flag, "ARKodeSetLinear", 1)) { return 1; } + switch (uopts.controller) + { + case (ARK_ADAPT_PID): C = SUNAdaptController_PID(ctx); break; + case (ARK_ADAPT_PI): C = SUNAdaptController_PI(ctx); break; + case (ARK_ADAPT_I): C = SUNAdaptController_I(ctx); break; + case (ARK_ADAPT_EXP_GUS): C = SUNAdaptController_ExpGus(ctx); break; + case (ARK_ADAPT_IMP_GUS): C = SUNAdaptController_ImpGus(ctx); break; + case (ARK_ADAPT_IMEX_GUS): C = SUNAdaptController_ImExGus(ctx); break; + } + flag = ARKodeSetAdaptController(arkode_mem, C); + if (check_flag(&flag, "ARKodeSetAdaptController", 1)) { return 1; } } // Set max steps between outputs @@ -377,22 +439,26 @@ int main(int argc, char* argv[]) // Free MPI Cartesian communicator MPI_Comm_free(&(udata.comm_c)); + (void)SUNAdaptController_Destroy(C); // Free timestep adaptivity controller ARKodeFree(&arkode_mem); - SUNLinSolFree(LS); + if (uopts.implicit) + { + SUNLinSolFree(LS); - // Free the SuperLU_DIST structures (also frees user allocated arrays - // A_data, A_col_idxs, and A_row_ptrs) + // Free the SuperLU_DIST structures (also frees user allocated arrays + // A_data, A_col_idxs, and A_row_ptrs) #if defined(USE_SUPERLU_DIST) - if (uopts.ls == "sludist") - { - PStatFree(&A_stat); - dScalePermstructFree(&A_scaleperm); - dLUstructFree(&A_lu); - Destroy_CompRowLoc_Matrix_dist(&A_super); - superlu_gridexit(&grid); - } + if (uopts.ls == "sludist") + { + PStatFree(&A_stat); + dScalePermstructFree(&A_scaleperm); + dLUstructFree(&A_lu); + Destroy_CompRowLoc_Matrix_dist(&A_super); + superlu_gridexit(&grid); + } #endif + } // Free vectors #if defined(USE_HIP) || defined(USE_CUDA) @@ -409,6 +475,26 @@ int main(int argc, char* argv[]) return 0; } +// ----------------------------------------------------------------------------- +// Dominant eigenvalue estimation function +// ----------------------------------------------------------------------------- + +static int dom_eig(sunrealtype t, N_Vector y, N_Vector fn, sunrealtype* lambdaR, + sunrealtype* lambdaI, void* user_data, N_Vector temp1, + N_Vector temp2, N_Vector temp3) +{ + // Access problem data + UserData* udata = (UserData*)user_data; + + // Fill in spectral radius value + *lambdaR = -SUN_RCONST(8.0) * std::max(udata->kx / udata->dx / udata->dx, + udata->ky / udata->dy / udata->dy); + *lambdaI = SUN_RCONST(0.0); + + // return with success + return 0; +} + // ----------------------------------------------------------------------------- // UserOptions Helper functions // ----------------------------------------------------------------------------- @@ -445,6 +531,13 @@ int UserOptions::parse_args(vector& args, bool outproc) args.erase(it, it + 2); } + it = find(args.begin(), args.end(), "--explicitSTS"); + if (it != args.end()) + { + implicit = false; + args.erase(it); + } + it = find(args.begin(), args.end(), "--order"); if (it != args.end()) { @@ -452,6 +545,13 @@ int UserOptions::parse_args(vector& args, bool outproc) args.erase(it, it + 2); } + it = find(args.begin(), args.end(), "--lsrkmethod"); + if (it != args.end()) + { + lsrkmethod = (ARKODE_LSRKMethodType)stoi(*(it + 1)); + args.erase(it, it + 2); + } + it = find(args.begin(), args.end(), "--controller"); if (it != args.end()) { @@ -525,16 +625,22 @@ void UserOptions::help() cout << "Integrator command line options:" << endl; cout << " --rtol : relative tolerance" << endl; cout << " --atol : absolute tolerance" << endl; + cout << " --controller : time step adaptivity controller" << endl; + cout << " --fixedstep : used fixed step size" << endl; + cout << " --explicitSTS : use LSRKStep (instead of ARKStep)" << endl; + cout << endl; + cout << "Implicit (ARKStep) solver command line options:" << endl; cout << " --nonlinear : disable linearly implicit flag" << endl; cout << " --order : method order" << endl; - cout << " --fixedstep : used fixed step size" << endl; - cout << " --controller : time step adaptivity controller" << endl; cout << " --ls : linear solver" << endl; cout << " --lsinfo : output residual history" << endl; cout << " --liniters : max number of iterations" << endl; cout << " --epslin : linear tolerance factor" << endl; cout << " --noprec : disable preconditioner" << endl; cout << " --msbp : max steps between prec setups" << endl; + cout << endl; + cout << "Explicit STS (LSRKStep) solver command line options:" << endl; + cout << " --lsrkmethod : LSRK method choice" << endl; } // Print user options @@ -545,39 +651,57 @@ void UserOptions::print() cout << " --------------------------------- " << endl; cout << " rtol = " << rtol << endl; cout << " atol = " << atol << endl; - cout << " hfixed = " << hfixed << endl; - cout << " order = " << order << endl; cout << " controller = " << controller << endl; - cout << " max steps = " << maxsteps << endl; - cout << " linear RHS = " << linear << endl; + cout << " hfixed = " << hfixed << endl; cout << " --------------------------------- " << endl; - cout << endl; - if (ls == "sludist") + if (implicit) { - cout << " Linear solver options:" << endl; + cout << " ARKStep options:" << endl; cout << " --------------------------------- " << endl; + cout << " order = " << order << endl; + cout << " max steps = " << maxsteps << endl; + cout << " linear RHS = " << linear << endl; + cout << " --------------------------------- " << endl; + cout << endl; + if (ls == "sludist") + { + cout << " Linear solver options:" << endl; + cout << " --------------------------------- " << endl; #if defined(HAVE_HIP) - cout << " LS = SuperLU_DIST (HIP enabled)" << endl; + cout << " LS = SuperLU_DIST (HIP enabled)" << endl; #elif defined(HAVE_CUDA) - cout << " LS = SuperLU_DIST (CUDA enabled)" << endl; + cout << " LS = SuperLU_DIST (CUDA enabled)" << endl; #else - cout << " LS = SuperLU_DIST" << endl; + cout << " LS = SuperLU_DIST" << endl; #endif - cout << " LS info = " << lsinfo << endl; - cout << " msbp = " << msbp << endl; - cout << " --------------------------------- " << endl; + cout << " LS info = " << lsinfo << endl; + cout << " msbp = " << msbp << endl; + cout << " --------------------------------- " << endl; + } + else + { + cout << " Linear solver options:" << endl; + cout << " --------------------------------- " << endl; + cout << " LS = " << ls << endl; + cout << " precond = " << preconditioning << endl; + cout << " LS info = " << lsinfo << endl; + cout << " LS iters = " << liniters << endl; + cout << " msbp = " << msbp << endl; + cout << " epslin = " << epslin << endl; + cout << " --------------------------------- " << endl; + } } else { - cout << " Linear solver options:" << endl; + cout << " LSRKStep options:" << endl; cout << " --------------------------------- " << endl; - cout << " LS = " << ls << endl; - cout << " precond = " << preconditioning << endl; - cout << " LS info = " << lsinfo << endl; - cout << " LS iters = " << liniters << endl; - cout << " msbp = " << msbp << endl; - cout << " epslin = " << epslin << endl; + switch (lsrkmethod) + { + case (ARKODE_LSRK_RKC_2): cout << " method = RKC_2 " << endl; break; + case (ARKODE_LSRK_RKL_2): cout << " method = RKL_2 " << endl; break; + default: cout << " ERROR: illegal lsrkmethod " << endl; + } cout << " --------------------------------- " << endl; } } diff --git a/doc/arkode/guide/source/Constants.rst b/doc/arkode/guide/source/Constants.rst index 22363ef243..106444dec9 100644 --- a/doc/arkode/guide/source/Constants.rst +++ b/doc/arkode/guide/source/Constants.rst @@ -317,6 +317,25 @@ contains the ARKODE output constants. +-----------------------------------------------+------------------------------------------------------------+ | | | +-----------------------------------------------+------------------------------------------------------------+ + | **LSRK method types** | | + +-----------------------------------------------+------------------------------------------------------------+ + | :index:`ARKODE_LSRK_RKC_2` | 2nd order Runge-Kutta-Chebyshev (RKC) method | + | | :c:enumerator:`ARKODE_LSRK_RKC_2` | + +-----------------------------------------------+------------------------------------------------------------+ + | :index:`ARKODE_LSRK_RKL_2` | 2nd order Runge-Kutta-Legendre (RKL) method | + | | :c:enumerator:`ARKODE_LSRK_RKL_2` | + +-----------------------------------------------+------------------------------------------------------------+ + | :index:`ARKODE_LSRK_SSP_S_2` | Optimal 2nd order s-stage SSP RK method | + | | :c:enumerator:`ARKODE_LSRK_SSP_S_2` | + +-----------------------------------------------+------------------------------------------------------------+ + | :index:`ARKODE_LSRK_SSP_S_3` | Optimal 3rd order s-stage SSP RK method | + | | :c:enumerator:`ARKODE_LSRK_SSP_S_3` | + +-----------------------------------------------+------------------------------------------------------------+ + | :index:`ARKODE_LSRK_SSP_10_4` | Optimal 4th order 10-stage SSP RK method | + | | :c:enumerator:`ARKODE_LSRK_SSP_10_4` | + +-----------------------------------------------+------------------------------------------------------------+ + | | | + +-----------------------------------------------+------------------------------------------------------------+ | **MRI method types** | | +-----------------------------------------------+------------------------------------------------------------+ | :index:`MRISTEP_EXPLICIT` | Use an explicit (at the slow time scale) MRI method. | diff --git a/doc/arkode/guide/source/Mathematics.rst b/doc/arkode/guide/source/Mathematics.rst index 9d9835c53d..c19e8a9fea 100644 --- a/doc/arkode/guide/source/Mathematics.rst +++ b/doc/arkode/guide/source/Mathematics.rst @@ -61,13 +61,15 @@ for interpolated solution output. We then discuss the current suite of time-stepping modules supplied with ARKODE, including the ARKStep module for :ref:`additive Runge--Kutta methods `, the ERKStep module that is optimized for :ref:`explicit Runge--Kutta -methods `, and the MRIStep module for :ref:`multirate -infinitesimal step (MIS), multirate infinitesimal GARK (MRI-GARK), and -implicit-explicit MRI-GARK (IMEX-MRI-GARK) methods `. -We then discuss the :ref:`adaptive temporal error controllers -` shared by the time-stepping modules, including -discussion of our choice of norms for measuring errors within various components -of the solver. +methods `, the SPRKStep module that provides +:ref:`symplectic partitioned Runge--Kutta methods ` +for Hamiltonian dynamics, the MRIStep module that provides a variety of +:ref:`multirate infinitesimal methods `, and the +LSRKStep module that supports :ref:`low-storage Runge--Kutta methods +`. We then discuss the :ref:`adaptive temporal +error controllers ` shared by the time-stepping +modules, including discussion of our choice of norms for measuring errors +within various components of the solver. We then discuss the nonlinear and linear solver strategies used by ARKODE for solving implicit algebraic systems that arise in computing each @@ -729,6 +731,61 @@ may supply their own method by defining and attaching a coupling table, see :numref:`ARKODE.Usage.MRIStep.MRIStepCoupling` for more information. +.. _ARKODE.Mathematics.LSRK: + +LSRKStep -- Low-Storage Runge--Kutta methods +============================================ + +The LSRKStep time-stepping module in ARKODE supports a variety of so-called +"low-storage" Runge--Kutta (LSRK) methods, :cite:p:`VSH:04, MBA:14, K:08, FCS:22`. This category includes traditional explicit +fixed-step and low-storage Runge--Kutta methods, adaptive +low-storage Runge--Kutta methods, and others. These are characterized by coefficient tables +that have an exploitable structure, such that their implementation does not require +that all stages be stored simultaneously. At present, this module supports explicit, +adaptive "super-time-stepping" (STS) and "strong-stability-preserving" (SSP) methods. + +The LSRK time-stepping module in ARKODE currently supports IVP +of the form :eq:`ARKODE_IVP_simple_explicit`, i.e., unlike the more general problem form :eq:`ARKODE_IMEX_IVP`, LSRKStep +requires that problems have an identity mass matrix (i.e., :math:`M(t)=I`) +and that the right-hand side function is not split into separate +components. + +LSRKStep currently supports two families of second-order, explicit, and temporally adaptive STS methods: +Runge--Kutta--Chebyshev (RKC), :cite:p:`VSH:04` and Runge--Kutta--Legendre (RKL), :cite:p:`MBA:14`. These methods have the form + +.. math:: + z_0 &= y_n,\\ + z_1 &= z_0 + h \tilde{\mu}_1 f(t_n,z_0),\\ + z_j &= (1-\mu_j-\nu_j)z_0 + \mu_j z_{j-1} + \nu_jz_{j-2} + h \tilde{\gamma}_j f(t_n,z_0) + h \tilde{\mu}_j f(t_n + c_{j-1}h, z_{j-1}) \\ + y_{n+1} &= z_s. + :label: ARKODE_RKC_RKL + +The corresponding coefficients can be found in :cite:p:`VSH:04` and :cite:p:`MBA:14`, respectively. + +LSRK methods of STS type are designed for stiff problems characterized by +having Jacobians with eigenvalues that have large real and small imaginary parts. +While those problems are traditionally treated using implicit methods, STS methods +are explicit. To achieve stability for these stiff problems, STS methods use more stages than +conventional Runge-Kutta methods to extend the stability region along the negative +real axis. The extent of this stability region is proportional to the square of the number +of stages used. + +LSRK methods of the SSP type are designed to preserve the so-called "strong-stability" properties of advection-type equations. +For details, see :cite:p:`K:08`. +The SSPRK methods in ARKODE use the following Shu--Osher representation :cite:p:`SO:88` of explicit Runge--Kutta methods: + +.. math:: + z_1 &= y_n,\\ + z_i &= \sum_{j = 1}^{i-1} \left(\alpha_{i,j}y_j + \beta_{i,j}h f(t_n + c_jh, z_j)\right),\\ + y_{n+1} &= z_s. + :label: ARKODE_SSP + +The coefficients of Shu--Osher representation is not uniquely determined by Butcher table :cite:p:`SR:02`. +In particular, the methods SSP(s,2), SSP(s,3), and SSP(10,4) implemented herein and presented in +:cite:p:`K:08` have "almost" all zero coefficients appearing in :math:`\alpha_{i,i-1}` and +:math:`\beta_{i,i-1}`. This feature facilitates their implementation in a low-storage manner. The +corresponding coefficients and embedding weights can be found in :cite:p:`K:08` and :cite:p:`FCS:22`, respectively. + .. _ARKODE.Mathematics.Error.Norm: Error norms diff --git a/doc/arkode/guide/source/Usage/ARKStep/User_callable.rst b/doc/arkode/guide/source/Usage/ARKStep/User_callable.rst index bef540030c..0795c81651 100644 --- a/doc/arkode/guide/source/Usage/ARKStep/User_callable.rst +++ b/doc/arkode/guide/source/Usage/ARKStep/User_callable.rst @@ -98,7 +98,7 @@ ARKStep tolerance specification functions * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory was ``NULL`` * *ARK_NO_MALLOC* if the ARKStep memory was not allocated by the time-stepping module - * *ARK_ILL_INPUT* if an argument has an illegal value (e.g. a negative tolerance). + * *ARK_ILL_INPUT* if an argument had an illegal value (e.g. a negative tolerance). .. deprecated:: 6.1.0 @@ -121,7 +121,7 @@ ARKStep tolerance specification functions * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory was ``NULL`` * *ARK_NO_MALLOC* if the ARKStep memory was not allocated by the time-stepping module - * *ARK_ILL_INPUT* if an argument has an illegal value (e.g. a negative tolerance). + * *ARK_ILL_INPUT* if an argument had an illegal value (e.g. a negative tolerance). .. deprecated:: 6.1.0 @@ -161,7 +161,7 @@ ARKStep tolerance specification functions * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory was ``NULL`` * *ARK_NO_MALLOC* if the ARKStep memory was not allocated by the time-stepping module - * *ARK_ILL_INPUT* if an argument has an illegal value (e.g. a negative tolerance). + * *ARK_ILL_INPUT* if an argument had an illegal value (e.g. a negative tolerance). .. deprecated:: 6.1.0 @@ -181,7 +181,7 @@ ARKStep tolerance specification functions * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory was ``NULL`` * *ARK_NO_MALLOC* if the ARKStep memory was not allocated by the time-stepping module - * *ARK_ILL_INPUT* if an argument has an illegal value (e.g. a negative tolerance). + * *ARK_ILL_INPUT* if an argument had an illegal value (e.g. a negative tolerance). .. deprecated:: 6.1.0 @@ -536,7 +536,7 @@ Optional inputs for ARKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Does not change the *user_data* pointer or any @@ -572,7 +572,7 @@ Optional inputs for ARKStep * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory or interpolation module are ``NULL`` * *ARK_INTERP_FAIL* if this is called after :c:func:`ARKStepEvolve()` - * *ARK_ILL_INPUT* if an argument has an illegal value or the + * *ARK_ILL_INPUT* if an argument had an illegal value or the interpolation module has already been initialized **Notes:** @@ -622,7 +622,7 @@ Optional inputs for ARKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** This parameter can be ``stdout`` or ``stderr``, although the @@ -651,7 +651,7 @@ Optional inputs for ARKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Pass 0.0 to return ARKStep to the default (adaptive-step) mode. @@ -712,7 +712,7 @@ Optional inputs for ARKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Pass 0.0 to use the default value. @@ -742,7 +742,7 @@ Optional inputs for ARKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** The default value is 10; set *mxhnil* to zero to specify @@ -768,7 +768,7 @@ Optional inputs for ARKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Passing *mxsteps* = 0 results in ARKStep using the @@ -792,7 +792,7 @@ Optional inputs for ARKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Pass *hmax* :math:`\le 0.0` to set the default value of :math:`\infty`. @@ -813,7 +813,7 @@ Optional inputs for ARKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Pass *hmin* :math:`\le 0.0` to set the default value of 0. @@ -835,7 +835,7 @@ Optional inputs for ARKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** The default is that no stop time is imposed. @@ -908,7 +908,7 @@ Optional inputs for ARKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** If specified, the pointer to *user_data* is passed to all @@ -936,7 +936,7 @@ Optional inputs for ARKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** The default value is 7; set *maxnef* :math:`\le 0` @@ -959,7 +959,7 @@ Optional inputs for ARKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Should only be called after the method order and integration @@ -1214,7 +1214,7 @@ Set additive RK tables via their names :c:func:`ARKStepSetTableName()` int **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** For explicit methods, the allowed values are :math:`2 \le` @@ -1243,7 +1243,7 @@ Set additive RK tables via their names :c:func:`ARKStepSetTableName()` int **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** This is automatically deduced when neither of the function @@ -1263,7 +1263,7 @@ Set additive RK tables via their names :c:func:`ARKStepSetTableName()` int **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** This is automatically deduced when the function pointer `fi` @@ -1286,7 +1286,7 @@ Set additive RK tables via their names :c:func:`ARKStepSetTableName()` int **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** This is automatically deduced when the function pointer `fe` @@ -1309,7 +1309,7 @@ Set additive RK tables via their names :c:func:`ARKStepSetTableName()` int **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** For a description of the :c:type:`ARKodeButcherTable` type and related @@ -1366,7 +1366,7 @@ Set additive RK tables via their names :c:func:`ARKStepSetTableName()` int **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** The allowable values for both the *itable* and *etable* arguments @@ -1408,7 +1408,7 @@ Set additive RK tables via their names :c:func:`ARKStepSetTableName()` int **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** The allowable values for both the *itable* and *etable* arguments @@ -1477,7 +1477,7 @@ Optional inputs for time step adaptivity **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** This function should focus on accuracy-based time step @@ -1515,7 +1515,7 @@ Optional inputs for time step adaptivity **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** If custom parameters are supplied, they will be checked @@ -1549,7 +1549,7 @@ Optional inputs for time step adaptivity **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** This should be called prior to calling :c:func:`ARKStepEvolve()`, and can only be @@ -1574,7 +1574,7 @@ Optional inputs for time step adaptivity **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any non-positive parameter will imply a reset to the default @@ -1598,7 +1598,7 @@ Optional inputs for time step adaptivity **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any value below 1.0 will imply a reset to the default value. @@ -1624,7 +1624,7 @@ Optional inputs for time step adaptivity **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any interval *not* containing 1.0 will imply a reset to the default values. @@ -1648,7 +1648,7 @@ Optional inputs for time step adaptivity **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any value outside the interval :math:`(0,1]` will imply a reset to the default value. @@ -1670,7 +1670,7 @@ Optional inputs for time step adaptivity **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any value outside the interval :math:`(0,1]` will imply a reset to the default value. @@ -1693,7 +1693,7 @@ Optional inputs for time step adaptivity **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any value :math:`\le 1.0` will imply a reset to the default value. @@ -1715,7 +1715,7 @@ Optional inputs for time step adaptivity **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any value :math:`\le 1.0` will imply a reset to the default @@ -1740,7 +1740,7 @@ Optional inputs for time step adaptivity **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any value outside the interval :math:`(0,1)` will imply a reset to @@ -1763,7 +1763,7 @@ Optional inputs for time step adaptivity **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any value :math:`\le 0` will imply a reset to the default @@ -1787,7 +1787,7 @@ Optional inputs for time step adaptivity **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any value :math:`\le 0` will imply a reset to the default value. @@ -1811,7 +1811,7 @@ Optional inputs for time step adaptivity **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** This function should return an estimate of the absolute @@ -1845,7 +1845,7 @@ Optional inputs for implicit stage solves **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Tightens the linear solver tolerances and takes only a @@ -1876,7 +1876,7 @@ Optional inputs for implicit stage solves **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** This is the default behavior of ARKStep, so the function @@ -1922,7 +1922,7 @@ Optional inputs for implicit stage solves **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** The default value is 0. If *method* is set to an @@ -2001,7 +2001,7 @@ Optional inputs for implicit stage solves **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value or if the SUNNONLINSOL module is ``NULL`` + * *ARK_ILL_INPUT* if an argument had an illegal value or if the SUNNONLINSOL module is ``NULL`` * *ARK_NLS_OP_ERR* if the SUNNONLINSOL object returned a failure flag **Notes:** @@ -2025,7 +2025,7 @@ Optional inputs for implicit stage solves **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** The default value is 0.1; set *nlscoef* :math:`\le 0` @@ -2047,7 +2047,7 @@ Optional inputs for implicit stage solves **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any non-positive parameter will imply a reset to the default value. @@ -2070,7 +2070,7 @@ Optional inputs for implicit stage solves **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any non-positive parameter will imply a reset to the default value. @@ -2095,7 +2095,7 @@ Optional inputs for implicit stage solves **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** The default value is 10; set *maxncf* :math:`\le 0` @@ -2162,7 +2162,7 @@ Optional inputs for the ARKLS linear solver interface **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any non-positive parameter will imply a reset to the default value. @@ -2321,7 +2321,7 @@ Optional inputs for matrix-based ``SUNLinearSolver`` modules * *ARKLS_SUCCESS* if successful * *ARKLS_MEM_NULL* if the ARKStep memory was ``NULL`` * *ARKLS_MASSMEM_NULL* if the mass matrix solver memory was ``NULL`` - * *ARKLS_ILL_INPUT* if an argument has an illegal value + * *ARKLS_ILL_INPUT* if an argument had an illegal value **Notes:** This routine must be called after the ARKLS mass matrix @@ -2701,7 +2701,7 @@ Rootfinding optional input functions **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** The default behavior is to monitor for both zero-crossing directions. @@ -4225,7 +4225,7 @@ vector. * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory was ``NULL`` * *ARK_MEM_FAIL* if a memory allocation failed - * *ARK_ILL_INPUT* if an argument has an illegal value. + * *ARK_ILL_INPUT* if an argument had an illegal value. **Notes:** All previously set options are retained but may be updated by calling @@ -4257,7 +4257,7 @@ ARKStep reset function * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory was ``NULL`` * *ARK_MEM_FAIL* if a memory allocation failed - * *ARK_ILL_INPUT* if an argument has an illegal value. + * *ARK_ILL_INPUT* if an argument had an illegal value. **Notes:** By default the next call to :c:func:`ARKStepEvolve()` will use the step size @@ -4304,7 +4304,7 @@ ARKStep system resize function * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ARKStep memory was ``NULL`` * *ARK_NO_MALLOC* if *arkode_mem* was not allocated. - * *ARK_ILL_INPUT* if an argument has an illegal value. + * *ARK_ILL_INPUT* if an argument had an illegal value. **Notes:** If an error occurred, :c:func:`ARKStepResize()` also sends an error @@ -4375,7 +4375,7 @@ wrap an ARKStep memory block as an :c:type:`MRIStepInnerStepper`. **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_FAIL* if a memory allocation failed - * *ARK_ILL_INPUT* if an argument has an illegal value. + * *ARK_ILL_INPUT* if an argument had an illegal value. **Example usage:** .. code-block:: C diff --git a/doc/arkode/guide/source/Usage/ERKStep/User_callable.rst b/doc/arkode/guide/source/Usage/ERKStep/User_callable.rst index 6c7b7b9cdf..1f4dc8f140 100644 --- a/doc/arkode/guide/source/Usage/ERKStep/User_callable.rst +++ b/doc/arkode/guide/source/Usage/ERKStep/User_callable.rst @@ -95,7 +95,7 @@ ERKStep tolerance specification functions * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory was ``NULL`` * *ARK_NO_MALLOC* if the ERKStep memory was not allocated by the time-stepping module - * *ARK_ILL_INPUT* if an argument has an illegal value (e.g. a negative tolerance). + * *ARK_ILL_INPUT* if an argument had an illegal value (e.g. a negative tolerance). .. deprecated:: 6.1.0 @@ -119,7 +119,7 @@ ERKStep tolerance specification functions * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory was ``NULL`` * *ARK_NO_MALLOC* if the ERKStep memory was not allocated by the time-stepping module - * *ARK_ILL_INPUT* if an argument has an illegal value (e.g. a negative tolerance). + * *ARK_ILL_INPUT* if an argument had an illegal value (e.g. a negative tolerance). .. deprecated:: 6.1.0 @@ -312,7 +312,7 @@ Optional inputs for ERKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Does not change problem-defining function pointer *f* @@ -350,7 +350,7 @@ Optional inputs for ERKStep * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory or interpolation module are ``NULL`` * *ARK_INTERP_FAIL* if this is called after :c:func:`ERKStepEvolve()` - * *ARK_ILL_INPUT* if an argument has an illegal value or the + * *ARK_ILL_INPUT* if an argument had an illegal value or the interpolation module has already been initialized **Notes:** @@ -402,7 +402,7 @@ Optional inputs for ERKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** This parameter can be ``stdout`` or ``stderr``, although the @@ -431,7 +431,7 @@ Optional inputs for ERKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Pass 0.0 to return ERKStep to the default (adaptive-step) mode. @@ -491,7 +491,7 @@ Optional inputs for ERKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Pass 0.0 to use the default value. @@ -522,7 +522,7 @@ Optional inputs for ERKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** The default value is 10; set *mxhnil* to zero to specify @@ -549,7 +549,7 @@ Optional inputs for ERKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Passing *mxsteps* = 0 results in ERKStep using the @@ -574,7 +574,7 @@ Optional inputs for ERKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Pass *hmax* :math:`\le 0.0` to set the default value of :math:`\infty`. @@ -596,7 +596,7 @@ Optional inputs for ERKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Pass *hmin* :math:`\le 0.0` to set the default value of 0. @@ -619,7 +619,7 @@ Optional inputs for ERKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** The default is that no stop time is imposed. @@ -650,7 +650,7 @@ Optional inputs for ERKStep **Return value:** * *ARK_SUCCESS* if successful - * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` + * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` .. versionadded:: 5.6.0 @@ -695,7 +695,7 @@ Optional inputs for ERKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** If specified, the pointer to *user_data* is passed to all @@ -720,7 +720,7 @@ Optional inputs for ERKStep **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** The default value is 7; set *maxnef* :math:`\le 0` @@ -836,7 +836,7 @@ Optional inputs for IVP method selection **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** The allowed values are :math:`2 \le` *ord* :math:`\le @@ -863,7 +863,7 @@ Optional inputs for IVP method selection **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** @@ -896,7 +896,7 @@ Optional inputs for IVP method selection **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** *etable* should match an existing explicit method from @@ -919,7 +919,7 @@ Optional inputs for IVP method selection **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** *etable* should match an existing explicit method from @@ -982,7 +982,7 @@ the code, is provided in :numref:`ARKODE.Mathematics.Adaptivity`. **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** This function should focus on accuracy-based time step @@ -1019,7 +1019,7 @@ the code, is provided in :numref:`ARKODE.Mathematics.Adaptivity`. **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** If custom parameters are supplied, they will be checked @@ -1051,7 +1051,7 @@ the code, is provided in :numref:`ARKODE.Mathematics.Adaptivity`. **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** This should be called prior to calling :c:func:`ERKStepEvolve()`, and can only be @@ -1076,7 +1076,7 @@ the code, is provided in :numref:`ARKODE.Mathematics.Adaptivity`. **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any non-positive parameter will imply a reset to the default @@ -1101,7 +1101,7 @@ the code, is provided in :numref:`ARKODE.Mathematics.Adaptivity`. **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any value below 1.0 will imply a reset to the default value. @@ -1128,7 +1128,7 @@ the code, is provided in :numref:`ARKODE.Mathematics.Adaptivity`. **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any interval *not* containing 1.0 will imply a reset to the default values. @@ -1151,7 +1151,7 @@ the code, is provided in :numref:`ARKODE.Mathematics.Adaptivity`. **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any value outside the interval :math:`(0,1]` will imply a reset to the default value. @@ -1175,7 +1175,7 @@ the code, is provided in :numref:`ARKODE.Mathematics.Adaptivity`. **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any value :math:`\le 1.0` will imply a reset to the default value. @@ -1198,7 +1198,7 @@ the code, is provided in :numref:`ARKODE.Mathematics.Adaptivity`. **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any value :math:`\le 1.0` will imply a reset to the default @@ -1224,7 +1224,7 @@ the code, is provided in :numref:`ARKODE.Mathematics.Adaptivity`. **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any value :math:`\ge 1.0` or :math:`\le 0.0` will imply a reset to @@ -1248,7 +1248,7 @@ the code, is provided in :numref:`ARKODE.Mathematics.Adaptivity`. **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any non-positive parameter will imply a reset to the default @@ -1273,7 +1273,7 @@ the code, is provided in :numref:`ARKODE.Mathematics.Adaptivity`. **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any non-positive parameter will imply a reset to the default value. @@ -1298,7 +1298,7 @@ the code, is provided in :numref:`ARKODE.Mathematics.Adaptivity`. **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** This function should return an estimate of the absolute @@ -1338,7 +1338,7 @@ Rootfinding optional input functions **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** The default behavior is to monitor for both zero-crossing directions. @@ -1877,7 +1877,7 @@ Main solver optional output functions **Return value:** * *ARK_SUCCESS* if successful - * *ARK_MEM_NULL* if the ARKStep memory was ``NULL`` + * *ARK_MEM_NULL* if the ERKStep memory was ``NULL`` .. versionadded:: 5.3.0 @@ -2016,7 +2016,7 @@ user must call the function :c:func:`ERKStepReInit()`. The new problem must have the same size as the previous one. This routine retains the current settings for all ERKstep module options and performs the same input checking and initializations that are done in -:c:func:`ERKStepCreate`, but it performs no memory allocation as is +:c:func:`ERKStepCreate`, but it performs no memory allocation as it assumes that the existing internal memory is sufficient for the new problem. A call to this re-initialization routine deletes the solution history that was stored internally during the previous @@ -2065,7 +2065,7 @@ vector. * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory was ``NULL`` * *ARK_MEM_FAIL* if a memory allocation failed - * *ARK_ILL_INPUT* if an argument has an illegal value. + * *ARK_ILL_INPUT* if an argument had an illegal value. **Notes:** All previously set options are retained but may be updated by calling @@ -2097,7 +2097,7 @@ ERKStep reset function * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory was ``NULL`` * *ARK_MEM_FAIL* if a memory allocation failed - * *ARK_ILL_INPUT* if an argument has an illegal value. + * *ARK_ILL_INPUT* if an argument had an illegal value. **Notes:** By default the next call to :c:func:`ERKStepEvolve()` will use the step size @@ -2146,7 +2146,7 @@ ERKStep system resize function * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the ERKStep memory was ``NULL`` * *ARK_NO_MALLOC* if *arkode_mem* was not allocated. - * *ARK_ILL_INPUT* if an argument has an illegal value. + * *ARK_ILL_INPUT* if an argument had an illegal value. **Notes:** If an error occurred, :c:func:`ERKStepResize()` also sends an error diff --git a/doc/arkode/guide/source/Usage/General.rst b/doc/arkode/guide/source/Usage/General.rst index c0980dfcee..745e9e9ed9 100644 --- a/doc/arkode/guide/source/Usage/General.rst +++ b/doc/arkode/guide/source/Usage/General.rst @@ -56,7 +56,8 @@ to the SUNDIALS core header file. #include // ARKStep provides explicit, implicit, IMEX additive RK methods. #include // MRIStep provides mutlirate RK methods. #include // SPRKStep provides symplectic partition RK methods. - + #include // LSRKStep provides low-storage RK methods. + Each of these define several types and various constants, include function prototypes, and include the shared ``arkode/arkode.h`` and ``arkode/arkode_ls.h`` header files. No other header files are required to be diff --git a/doc/arkode/guide/source/Usage/LSRKStep/User_callable.rst b/doc/arkode/guide/source/Usage/LSRKStep/User_callable.rst new file mode 100644 index 0000000000..8156ef0341 --- /dev/null +++ b/doc/arkode/guide/source/Usage/LSRKStep/User_callable.rst @@ -0,0 +1,360 @@ +.. ---------------------------------------------------------------- + Programmer(s): Mustafa Aggul @ SMU + ---------------------------------------------------------------- + SUNDIALS Copyright Start + Copyright (c) 2002-2024, Lawrence Livermore National Security + and Southern Methodist University. + All rights reserved. + + See the top-level LICENSE and NOTICE files for details. + + SPDX-License-Identifier: BSD-3-Clause + SUNDIALS Copyright End + ---------------------------------------------------------------- + +.. _ARKODE.Usage.LSRKStep.UserCallable: + +LSRKStep User-callable functions +================================== + +This section describes the LSRKStep-specific functions that may be called +by the user to setup and then solve an IVP using the LSRKStep time-stepping +module. As mentioned in Section :numref:`ARKODE.Usage.UserCallable`, +shared ARKODE-level routines may be used for the large majority of LSRKStep +configuration and use. In this section, we describe only those routines +that are specific to LSRKStep. + +As discussed in the main :ref:`ARKODE user-callable function introduction +`, each of ARKODE's time-stepping modules +clarifies the categories of user-callable functions that it supports. +LSRKStep supports the following categories: + +* temporal adaptivity + + + +.. _ARKODE.Usage.LSRKStep.Initialization: + +LSRKStep initialization functions +--------------------------------- + + +.. c:function:: void* LSRKStepCreateSTS(ARKRhsFn rhs, sunrealtype t0, N_Vector y0, SUNContext sunctx); + + This function allocates and initializes memory for a problem to + be solved using the LSRKStep time-stepping module in ARKODE. + + **Arguments:** + * *rhs* -- the name of the C function (of type :c:func:`ARKRhsFn()`) + defining the right-hand side function. + * *t0* -- the initial value of :math:`t`. + * *y0* -- the initial condition vector :math:`y(t_0)`. + * *sunctx* -- the :c:type:`SUNContext` object (see :numref:`SUNDIALS.SUNContext`) + + **Return value:** + If successful, a pointer to initialized problem memory + of type ``void*``, to be passed to all user-facing LSRKStep routines + listed below. If unsuccessful, a ``NULL`` pointer will be + returned, and an error message will be printed to ``stderr``. + + +.. c:function:: void* LSRKStepCreateSSP(ARKRhsFn rhs, sunrealtype t0, N_Vector y0, SUNContext sunctx); + + This function allocates and initializes memory for a problem to + be solved using the LSRKStep time-stepping module in ARKODE. + + **Arguments:** + * *rhs* -- the name of the C function (of type :c:func:`ARKRhsFn()`) + defining the right-hand side function. + * *t0* -- the initial value of :math:`t`. + * *y0* -- the initial condition vector :math:`y(t_0)`. + * *sunctx* -- the :c:type:`SUNContext` object (see :numref:`SUNDIALS.SUNContext`) + + **Return value:** + If successful, a pointer to initialized problem memory + of type ``void*``, to be passed to all user-facing LSRKStep routines + listed below. If unsuccessful, a ``NULL`` pointer will be + returned, and an error message will be printed to ``stderr``. + + +.. _ARKODE.Usage.LSRKStep.OptionalInputs: + +Optional input functions +------------------------- + + +.. c:function:: int LSRKStepSetMethod(void* arkode_mem, ARKODE_LSRKMethodType method); + + This function selects the LSRK method that should be used. The list of allowable + values for this input is below. + + **Arguments:** + * *arkode_mem* -- pointer to the LSRKStep memory block. + * *method* -- Type of the method. + + **Return value:** + * *ARK_SUCCESS* if successful + * *ARK_ILL_INPUT* if an argument had an illegal value (e.g. typo in the method type). + + .. note:: If this routine is not called, then LSRKStep will use the Runge--Kutta--Chebyshev method by default. + + +Allowable Method Families + +.. c:enum:: ARKODE_LSRKMethodType + + .. c:enumerator:: ARKODE_LSRK_RKC_2 + + Second order Runge--Kutta--Chebyshev method + + .. c:enumerator:: ARKODE_LSRK_RKL_2 + + Second order Runge--Kutta--Legendre method + + .. c:enumerator:: ARKODE_LSRK_SSP_S_2 + + Second order, s-stage SSP(s,2) method + + .. c:enumerator:: ARKODE_LSRK_SSP_S_3 + + Third order, s-stage SSP(s,3) method + + .. c:enumerator:: ARKODE_LSRK_SSP_10_4 + + Fourth order, 10-stage SSP(10,4) method + + +.. c:function:: int LSRKStepSetMethodByName(void* arkode_mem, const char* emethod); + + This function selects the LSRK method by name. The list of allowable values for this input is below. + + **Arguments:** + * *arkode_mem* -- pointer to the LSRKStep memory block. + * *emethod* -- Type of the method in strings. + + **Return value:** + * *ARK_SUCCESS* if successful + * *ARK_ILL_INPUT* if an argument had an illegal value (e.g. typo in the method type). + + .. note:: If one of these set method routines is not called, then LSRKStep will use the + Runge--Kutta--Chebyshev method by default. + + +.. c:function:: int LSRKStepSetDomEigFn(void* arkode_mem, ARKDomEigFn dom_eig); + + Specifies the dominant eigenvalue approximation routine to + be used for determining the number of stages that will be used by either the + RKC or RKL methods. + + **Arguments:** + * *arkode_mem* -- pointer to the LSRKStep memory block. + * *dom_eig* -- name of user-supplied dominant eigenvalue approximation function (of type :c:func:`ARKDomEigFn()`). + + **Return value:** + * *ARK_SUCCESS* if successful + * *ARKLS_MEM_NULL* if ``arkode_mem`` was ``NULL``. + * *ARK_ILL_INPUT* ``dom_eig = NULL`` and LSRKStep does not currently estimate this internally. + + .. note:: This function is currently required when either the RKC or RKL methods are used. + + +.. c:function:: int LSRKStepSetDomEigFrequency(void* arkode_mem, int nsteps); + + Specifies the number of steps after which the dominant eigenvalue information is + considered out-of-date, and should be recomputed. This only applies to RKL and RKC methods. + + **Arguments:** + * *arkode_mem* -- pointer to the LSRKStep memory block. + * *nsteps* -- the dominant eigenvalue re-computation update frequency. A value ``nsteps = 0`` indicates that the dominant eigenvalue will not change throughout the simulation. + + **Return value:** + * *ARK_SUCCESS* if successful + * *ARKLS_MEM_NULL* if ``arkode_mem`` was ``NULL``. + +.. note:: If LSRKStepSetDomEigFrequency routine is not called, then the default ``nsteps`` is set to :math:`25` as recommended in :cite:p:`VSH:04`. + Calling this function with ``nsteps < 0`` resets the default value while ``nsteps = 0`` refers to constant dominant eigenvalue. + + +.. c:function:: int LSRKStepSetMaxNumStages(void* arkode_mem, int stage_max_limit); + + Specifies the maximum number of stages allowed within each time step. This bound only applies to + RKL and RKC methods. + + **Arguments:** + * *arkode_mem* -- pointer to the LSRKStep memory block. + * *stage_max_limit* -- maximum allowed number of stages :math:`(>=2)` and :math:`(<=10000)`. + + **Return value:** + * *ARK_SUCCESS* if successful + * *ARKLS_MEM_NULL* if ``arkode_mem`` was ``NULL``. + +.. note:: If LSRKStepSetMaxNumStages routine is not called, then the default ``stage_max_limit`` is + set to :math:`200`. Calling this function with ``stage_max_limit < 2`` or ``stage_max_limit > 10000`` resets the default value. + + +.. c:function:: int LSRKStepSetDomEigSafetyFactor(void* arkode_mem, sunrealtype dom_eig_safety); + + Specifies a safety factor to use for the result of the dominant eigenvalue estimation function. + This value is used to scale the magnitude of the dominant eigenvalue, in the hope of ensuring + a sufficient number of stages for the method to be stable. This input is only used for RKC + and RKL methods. + + **Arguments:** + * *arkode_mem* -- pointer to the LSRKStep memory block. + * *dom_eig_safety* -- safety factor :math:`(\ge 1)`. + + **Return value:** + * *ARK_SUCCESS* if successful + * *ARKLS_MEM_NULL* if ``arkode_mem`` was ``NULL``. + +.. note:: If LSRKStepSetDomEigSafetyFactor routine is not called, then the default ``dom_eig_safety`` is + set to :math:`1.01`. Calling this function with ``dom_eig_safety < 1`` resets the default value. + + +.. c:function:: int LSRKStepSetSSPStageNum(void* arkode_mem, int num_of_stages); + + Sets the number of stages, ``s`` in ``SSP(s, p)`` methods. This input is only utilized by SSPRK methods. Thus, + this set routine must be called after calling LSRKStepSetMethod with an SSPRK method. + +* :c:enumerator:`ARKODE_LSRK_SSP_S_2` -- ``num_of_stages`` must be greater than or equal to 2 +* :c:enumerator:`ARKODE_LSRK_SSP_S_3` -- ``num_of_stages`` must be a perfect-square greater than or equal to 9 +* :c:enumerator:`ARKODE_LSRK_SSP_10_4` -- ``num_of_stages`` cannot be modified from 10, so this function should not be called. + + **Arguments:** + * *arkode_mem* -- pointer to the LSRKStep memory block. + * *num_of_stages* -- number of stages :math:`(>1)` for ``SSP(s,2)`` and :math:`(n^2 = s \geq 4)` for ``SSP(s,3)``. + + **Return value:** + * *ARK_SUCCESS* if successful + * *ARKLS_MEM_NULL* if ``arkode_mem`` was ``NULL``. + * *ARK_ILL_INPUT* if an argument had an illegal value (e.g. SSP method is not declared) + +.. note:: If LSRKStepSetSSPStageNum routine is not called, then the default ``num_of_stages`` is + set. Calling this function with ``num_of_stages <= 0`` resets the default values: + + * ``num_of_stages = 10`` for :c:enumerator:`ARKODE_LSRK_SSP_S_2` + * ``num_of_stages = 9`` for :c:enumerator:`ARKODE_LSRK_SSP_S_3` + * ``num_of_stages = 10`` for :c:enumerator:`ARKODE_LSRK_SSP_10_4` + +.. _ARKODE.Usage.LSRKStep.OptionalOutputs: + +Optional output functions +------------------------------ + +.. c:function:: int LSRKStepGetNumDomEigUpdates(void* arkode_mem, long int* dom_eig_num_evals); + + Returns the number of dominant eigenvalue evaluations (so far). + + **Arguments:** + * *arkode_mem* -- pointer to the LSRKStep memory block. + * *dom_eig_num_evals* -- number of calls to the user's ``dom_eig`` function. + + **Return value:** + * *ARK_SUCCESS* if successful + * *ARK_MEM_NULL* if the LSRKStep memory was ``NULL`` + + +.. c:function:: int LSRKStepGetMaxNumStages(void* arkode_mem, int* stage_max); + + Returns the max number of stages used in any single step (so far). + + **Arguments:** + * *arkode_mem* -- pointer to the LSRKStep memory block. + * *stage_max* -- max number of stages used. + + **Return value:** + * *ARK_SUCCESS* if successful + * *ARK_MEM_NULL* if the LSRKStep memory was ``NULL`` + +.. _ARKODE.Usage.LSRKStep.Reinitialization: + +LSRKStep re-initialization function +------------------------------------- + +To reinitialize the LSRKStep module for the solution of a new problem, +where a prior call to :c:func:`LSRKStepCreateSTS` or :c:func:`LSRKStepCreateSSP` +has been made, the user must call the function :c:func:`LSRKStepReInitSTS()` +or :c:func:`LSRKStepReInitSSP()`, accordingly. The new problem must have +the same size as the previous one. This routine retains the current settings +for all LSRKstep module options and performs the same input checking and +initializations that are done in :c:func:`LSRKStepCreateSTS` or +:c:func:`LSRKStepCreateSSP`, but it performs no memory allocation as it +assumes that the existing internal memory is sufficient for the new problem. +A call to this re-initialization routine deletes the solution history that +was stored internally during the previous integration, and deletes any +previously-set *tstop* value specified via a call to +:c:func:`ARKodeSetStopTime()`. Following a successful call to +:c:func:`LSRKStepReInitSTS()` or :c:func:`LSRKStepReInitSSP()`, +call :c:func:`ARKodeEvolve()` again for the solution of the new problem. + +One important use of the :c:func:`LSRKStepReInitSTS()` and +:c:func:`LSRKStepReInitSSP()` function is in the treating of jump +discontinuities in the RHS function. Except in cases of fairly small +jumps, it is usually more efficient to stop at each point of discontinuity +and restart the integrator with a readjusted ODE model, using a call to this +routine. To stop when the location of the discontinuity is known, simply +make that location a value of ``tout``. To stop when the location of +the discontinuity is determined by the solution, use the rootfinding feature. +In either case, it is critical that the RHS function *not* incorporate the +discontinuity, but rather have a smooth extension over the discontinuity, +so that the step across it (and subsequent rootfinding, if used) can be done +efficiently. Then use a switch within the RHS function (communicated through +``user_data``) that can be flipped between the stopping of the integration +and the restart, so that the restarted problem uses the new values (which +have jumped). Similar comments apply if there is to be a jump in the +dependent variable vector. + + +.. c:function:: int LSRKStepReInitSTS(void* arkode_mem, ARKRhsFn rhs, sunrealtype t0, N_Vector y0); + + Provides required problem specifications and re-initializes the + LSRKStep time-stepper module. + + **Arguments:** + * *arkode_mem* -- pointer to the LSRKStep memory block. + * *rhs* -- the name of the C function (of type :c:func:`ARKRhsFn()`) + defining the right-hand side function. + * *t0* -- the initial value of :math:`t`. + * *y0* -- the initial condition vector :math:`y(t_0)`. + + **Return value:** + * *ARK_SUCCESS* if successful + * *ARK_MEM_NULL* if the LSRKStep memory was ``NULL`` + * *ARK_MEM_FAIL* if memory allocation failed + * *ARK_NO_MALLOC* if memory allocation failed + * *ARK_CONTROLLER_ERR* if unable to reset error controller object + * *ARK_ILL_INPUT* if an argument had an illegal value. + + .. note:: + All previously set options are retained but may be updated by calling + the appropriate "Set" functions. + + If an error occurred, :c:func:`LSRKStepReInitSTS()` also + sends an error message to the error handler function. + +.. c:function:: int LSRKStepReInitSSP(void* arkode_mem, ARKRhsFn rhs, sunrealtype t0, N_Vector y0); + + Provides required problem specifications and re-initializes the + LSRKStep time-stepper module. + + **Arguments:** + * *arkode_mem* -- pointer to the LSRKStep memory block. + * *rhs* -- the name of the C function (of type :c:func:`ARKRhsFn()`) + defining the right-hand side function. + * *t0* -- the initial value of :math:`t`. + * *y0* -- the initial condition vector :math:`y(t_0)`. + + **Return value:** + * *ARK_SUCCESS* if successful + * *ARK_MEM_NULL* if the LSRKStep memory was ``NULL`` + * *ARK_MEM_FAIL* if memory allocation failed + * *ARK_NO_MALLOC* if memory allocation failed + * *ARK_CONTROLLER_ERR* if unable to reset error controller object + * *ARK_ILL_INPUT* if an argument had an illegal value. + + .. note:: + All previously set options are retained but may be updated by calling + the appropriate "Set" functions. + + If an error occurred, :c:func:`LSRKStepReInitSSP()` also + sends an error message to the error handler function. \ No newline at end of file diff --git a/doc/arkode/guide/source/Usage/LSRKStep/User_supplied.rst b/doc/arkode/guide/source/Usage/LSRKStep/User_supplied.rst new file mode 100644 index 0000000000..c9820d4598 --- /dev/null +++ b/doc/arkode/guide/source/Usage/LSRKStep/User_supplied.rst @@ -0,0 +1,52 @@ +.. ---------------------------------------------------------------- + Programmer(s): Mustafa Aggul @ SMU + ---------------------------------------------------------------- + SUNDIALS Copyright Start + Copyright (c) 2002-2024, Lawrence Livermore National Security + and Southern Methodist University. + All rights reserved. + + See the top-level LICENSE and NOTICE files for details. + + SPDX-License-Identifier: BSD-3-Clause + SUNDIALS Copyright End + ---------------------------------------------------------------- + +.. _LSRKSTEP.Usage.UserSupplied: + +User-supplied functions +============================= + +In addition to the required :c:type:`ARKRhsFn` arguments that define the IVP, +RKL and RKC methods additionally require an :c:type:`ARKDomEigFn`: function to +estimate the dominant eigenvalue. + + + + +.. _LSRKStep.Usage.dom_eig: + +The dominant eigenvalue estimation +---------------------------------- + +When running LSRKStep with either the RKC or RKL methods, the user must supply +a dominant eigenvalue estimation function of type :c:type:`ARKDomEigFn`: + +.. c:type:: int (*ARKDomEigFn)(sunrealtype t, N_Vector y, N_Vector fn, sunrealtype* lambdaR, sunrealtype* lambdaI, void* user_data, N_Vector temp1, N_Vector temp2, N_Vector temp3); + + These functions compute the dominant eigenvalue of the Jacobian of the ODE + right-hand side for a given value of the independent variable :math:`t` and + state vector :math:`y`. + + :param t: the current value of the independent variable. + :param y: the current value of the dependent variable vector. + :param fn: the current value of the vector :math:`f(t,y)`. + :param lambdaR: The real part of the dominant eigenvalue. + :param lambdaI: The imaginary part of the dominant eigenvalue. + :param user_data: the `user_data` pointer that was passed to + :c:func:`ARKodeSetUserData`. + :param tmp*: pointers to memory allocated to + variables of type ``N_Vector`` which can be used by an + ARKDomEigFn as temporary storage or work space. + + :return: An *ARKDomEigFn* should return 0 if successful and any nonzero for a failure. \ No newline at end of file diff --git a/doc/arkode/guide/source/Usage/LSRKStep/index.rst b/doc/arkode/guide/source/Usage/LSRKStep/index.rst new file mode 100644 index 0000000000..ca07bfffc3 --- /dev/null +++ b/doc/arkode/guide/source/Usage/LSRKStep/index.rst @@ -0,0 +1,32 @@ +.. ---------------------------------------------------------------- + Programmer(s): Mustafa Aggul @ SMU + ---------------------------------------------------------------- + SUNDIALS Copyright Start + Copyright (c) 2002-2024, Lawrence Livermore National Security + and Southern Methodist University. + All rights reserved. + + See the top-level LICENSE and NOTICE files for details. + + SPDX-License-Identifier: BSD-3-Clause + SUNDIALS Copyright End + ---------------------------------------------------------------- + +.. _ARKODE.Usage.LSRKStep: + +======================================= +Using the LSRKStep time-stepping module +======================================= + +This section is concerned with the use of the LSRKStep time-stepping +module for the solution of initial value problems (IVPs) in a C or C++ +language setting. Usage of LSRKStep follows that of the rest of ARKODE, +and so in this section we primarily focus on those usage aspects that +are specific to LSRKStep. + +.. toctree:: + :maxdepth: 1 + + User_callable + User_supplied + \ No newline at end of file diff --git a/doc/arkode/guide/source/Usage/MRIStep/User_callable.rst b/doc/arkode/guide/source/Usage/MRIStep/User_callable.rst index a6354a6037..967e8e7c12 100644 --- a/doc/arkode/guide/source/Usage/MRIStep/User_callable.rst +++ b/doc/arkode/guide/source/Usage/MRIStep/User_callable.rst @@ -134,7 +134,7 @@ MRIStep tolerance specification functions * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the MRIStep memory was ``NULL`` * *ARK_NO_MALLOC* if the MRIStep memory was not allocated by the time-stepping module - * *ARK_ILL_INPUT* if an argument has an illegal value (e.g. a negative tolerance). + * *ARK_ILL_INPUT* if an argument had an illegal value (e.g. a negative tolerance). .. deprecated:: 6.1.0 @@ -158,7 +158,7 @@ MRIStep tolerance specification functions * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the MRIStep memory was ``NULL`` * *ARK_NO_MALLOC* if the MRIStep memory was not allocated by the time-stepping module - * *ARK_ILL_INPUT* if an argument has an illegal value (e.g. a negative tolerance). + * *ARK_ILL_INPUT* if an argument had an illegal value (e.g. a negative tolerance). .. deprecated:: 6.1.0 @@ -361,7 +361,7 @@ MRIStep solver function internal time-stepping. (b) The linear solver initialization function (called by the - user after calling :c:func:`ARKStepCreate`) failed to set + user after calling :c:func:`MRIStepCreate`) failed to set the linear solver-specific *lsolve* field in *arkode_mem*. @@ -447,7 +447,7 @@ Optional inputs for MRIStep * *ARK_MEM_NULL* if the MRIStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** This function does not change problem-defining function pointers *fs* and *ff* or the *user_data* pointer. It also does not affect any data @@ -489,7 +489,7 @@ Optional inputs for MRIStep * *ARK_INTERP_FAIL* if this is called after :c:func:`MRIStepEvolve()` - * *ARK_ILL_INPUT* if an argument has an illegal value or the + * *ARK_ILL_INPUT* if an argument had an illegal value or the interpolation module has already been initialized **Notes:** Allowed values are between 0 and 5. @@ -544,7 +544,7 @@ Optional inputs for MRIStep * *ARK_MEM_NULL* if the MRIStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** This parameter can be ``stdout`` or ``stderr``, although the suggested approach is to specify a pointer to a unique file opened @@ -577,7 +577,7 @@ Optional inputs for MRIStep * *ARK_MEM_NULL* if the MRIStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** @@ -608,7 +608,7 @@ Optional inputs for MRIStep * *ARK_MEM_NULL* if the MRIStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** The default value is 10; set *mxhnil* to zero to specify this default. @@ -639,7 +639,7 @@ Optional inputs for MRIStep * *ARK_MEM_NULL* if the MRIStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Passing *mxsteps* = 0 results in MRIStep using the default value (500). @@ -669,7 +669,7 @@ Optional inputs for MRIStep * *ARK_MEM_NULL* if the MRIStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** @@ -701,7 +701,7 @@ Optional inputs for MRIStep **Return value:** * *ARK_SUCCESS* if successful - * *ARK_MEM_NULL* if the ARKStep memory is ``NULL`` + * *ARK_MEM_NULL* if the MRIStep memory is ``NULL`` .. versionadded:: 5.6.0 @@ -751,17 +751,16 @@ Optional inputs for MRIStep * *ARK_MEM_NULL* if the MRIStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** If specified, the pointer to *user_data* is passed to all user-supplied functions called by the outer integrator for which it is an argument; otherwise ``NULL`` is passed. - To attach a user data block to the inner integrator call the appropriate - *SetUserData* function for the inner integrator memory structure (e.g., - :c:func:`ARKStepSetUserData()` if the inner stepper is ARKStep). This pointer - may be the same as or different from the pointer attached to the outer - integrator depending on what is required by the user code. + To attach a user data block to the inner integrator call :c:func:`ARKodeSetUserData` + for the inner integrator memory structure. This pointer may be the same as or + different from the pointer attached to the outer integrator depending on what is + required by the user code. .. deprecated:: 6.1.0 @@ -867,7 +866,7 @@ Optional inputs for IVP method selection * *ARK_MEM_NULL* if the MRIStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** @@ -900,7 +899,7 @@ Optional inputs for implicit stage solves **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the MRIStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Tightens the linear solver tolerances and takes only a single Newton iteration. Calls :c:func:`MRIStepSetDeltaGammaMax()` @@ -930,7 +929,7 @@ Optional inputs for implicit stage solves **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the MRIStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** This is the default behavior of MRIStep, so the function is primarily useful to undo a previous call to @@ -970,7 +969,7 @@ Optional inputs for implicit stage solves **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the MRIStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** The default value is 0. If *method* is set to an undefined value, this default predictor will be used. @@ -996,7 +995,7 @@ Optional inputs for implicit stage solves **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the MRIStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value or if the SUNNONLINSOL module is ``NULL`` + * *ARK_ILL_INPUT* if an argument had an illegal value or if the SUNNONLINSOL module is ``NULL`` * *ARK_NLS_OP_ERR* if the SUNNONLINSOL object returned a failure flag **Notes:** The default value is 3; set *maxcor* :math:`\le 0` @@ -1019,7 +1018,7 @@ Optional inputs for implicit stage solves **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the MRIStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** The default value is 0.1; set *nlscoef* :math:`\le 0` to specify this default. @@ -1041,7 +1040,7 @@ Optional inputs for implicit stage solves **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the MRIStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any non-positive parameter will imply a reset to the default value. @@ -1064,7 +1063,7 @@ Optional inputs for implicit stage solves **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the MRIStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any non-positive parameter will imply a reset to the default value. @@ -1173,7 +1172,7 @@ Optional inputs for the ARKLS linear solver interface **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the MRIStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** Any non-positive parameter will imply a reset to the default value. @@ -1534,7 +1533,7 @@ Rootfinding optional input functions **Return value:** * *ARK_SUCCESS* if successful * *ARK_MEM_NULL* if the MRIStep memory is ``NULL`` - * *ARK_ILL_INPUT* if an argument has an illegal value + * *ARK_ILL_INPUT* if an argument had an illegal value **Notes:** The default behavior is to monitor for both zero-crossing directions. @@ -1990,7 +1989,7 @@ Main solver optional output functions * *ARK_SUCCESS* if successful - * *ARK_MEM_NULL* if the ARKStep memory was ``NULL`` + * *ARK_MEM_NULL* if the MRIStep memory was ``NULL`` .. versionadded:: 5.3.0 @@ -2645,7 +2644,7 @@ To reinitialize the MRIStep module for the solution of a new problem, where a prior call to :c:func:`MRIStepCreate()` has been made, the user must call the function :c:func:`MRIStepReInit()`. The new problem must have the same size as the previous one. This routine -retains the current settings for all ARKstep module options and +retains the current settings for all MRIStep module options and performs the same input checking and initializations that are done in :c:func:`MRIStepCreate()`, but it performs no memory allocation as is assumes that the existing internal memory is sufficient for the new @@ -2708,7 +2707,7 @@ vector. * *ARK_MEM_FAIL* if a memory allocation failed - * *ARK_ILL_INPUT* if an argument has an illegal value. + * *ARK_ILL_INPUT* if an argument had an illegal value. **Notes:** If the inner (fast) stepper also needs to be reinitialized, its @@ -2749,7 +2748,7 @@ MRIStep reset function * *ARK_MEM_FAIL* if a memory allocation failed - * *ARK_ILL_INPUT* if an argument has an illegal value. + * *ARK_ILL_INPUT* if an argument had an illegal value. **Notes:** If the inner (fast) stepper also needs to be reset, its reset function should @@ -2806,7 +2805,7 @@ MRIStep system resize function * *ARK_NO_MALLOC* if *arkode_mem* was not allocated. - * *ARK_ILL_INPUT* if an argument has an illegal value. + * *ARK_ILL_INPUT* if an argument had an illegal value. **Notes:** If an error occurred, :c:func:`MRIStepResize()` also sends an error message to the error handler function. diff --git a/doc/arkode/guide/source/Usage/SPRKStep/User_callable.rst b/doc/arkode/guide/source/Usage/SPRKStep/User_callable.rst index 73bde21f90..e4a88b1d18 100644 --- a/doc/arkode/guide/source/Usage/SPRKStep/User_callable.rst +++ b/doc/arkode/guide/source/Usage/SPRKStep/User_callable.rst @@ -227,7 +227,7 @@ Optional inputs for SPRKStep :retval ARK_SUCCESS: if successful :retval ARK_MEM_NULL: if the SPRKStep memory is ``NULL`` - :retval ARK_ILL_INPUT: if an argument has an illegal value + :retval ARK_ILL_INPUT: if an argument had an illegal value .. note:: @@ -262,7 +262,7 @@ Optional inputs for SPRKStep :retval ARK_SUCCESS: if successful :retval ARK_MEM_NULL: if the SPRKStep memory or interpolation module are ``NULL`` :retval ARK_INTERP_FAIL: if this is called after :c:func:`SPRKStepEvolve()` - :retval ARK_ILL_INPUT: if an argument has an illegal value or the + :retval ARK_ILL_INPUT: if an argument had an illegal value or the interpolation module has already been initialized .. note:: @@ -298,7 +298,7 @@ Optional inputs for SPRKStep :retval ARK_SUCCESS: if successful :retval ARK_MEM_NULL: if the SPRKStep memory is ``NULL`` - :retval ARK_ILL_INPUT: if an argument has an illegal value + :retval ARK_ILL_INPUT: if an argument had an illegal value .. deprecated:: 6.1.0 @@ -321,7 +321,7 @@ Optional inputs for SPRKStep :retval ARK_SUCCESS: if successful :retval ARK_MEM_NULL: if the SPRKStep memory is ``NULL`` - :retval ARK_ILL_INPUT: if an argument has an illegal value + :retval ARK_ILL_INPUT: if an argument had an illegal value .. deprecated:: 6.1.0 @@ -348,7 +348,7 @@ Optional inputs for SPRKStep :retval ARK_SUCCESS: if successful :retval ARK_MEM_NULL: if the SPRKStep memory is ``NULL`` - :retval ARK_ILL_INPUT: if an argument has an illegal value + :retval ARK_ILL_INPUT: if an argument had an illegal value .. deprecated:: 6.1.0 @@ -386,7 +386,7 @@ Optional inputs for SPRKStep :retval ARK_SUCCESS: if successful :retval ARK_MEM_NULL: if the SPRKStep memory is ``NULL`` - :retval ARK_ILL_INPUT: if an argument has an illegal value + :retval ARK_ILL_INPUT: if an argument had an illegal value .. deprecated:: 6.1.0 @@ -430,7 +430,7 @@ Optional inputs for IVP method selection :retval ARK_SUCCESS: if successful :retval ARK_MEM_NULL: if the SPRKStep memory is ``NULL`` - :retval ARK_ILL_INPUT: if an argument has an illegal value + :retval ARK_ILL_INPUT: if an argument had an illegal value .. warning:: @@ -451,7 +451,7 @@ Optional inputs for IVP method selection :retval ARK_SUCCESS: if successful :retval ARK_MEM_NULL: if the SPRKStep memory is ``NULL`` - :retval ARK_ILL_INPUT: if an argument has an illegal value + :retval ARK_ILL_INPUT: if an argument had an illegal value .. note:: @@ -472,7 +472,7 @@ Optional inputs for IVP method selection :retval ARK_SUCCESS: if successful :retval ARK_MEM_NULL: if the SPRKStep memory is ``NULL`` - :retval ARK_ILL_INPUT: if an argument has an illegal value + :retval ARK_ILL_INPUT: if an argument had an illegal value .. warning:: @@ -494,7 +494,7 @@ Optional inputs for IVP method selection :retval ARK_SUCCESS: if successful :retval ARK_MEM_NULL: if the SPRKStep memory is ``NULL`` - :retval ARK_ILL_INPUT: if an argument has an illegal value + :retval ARK_ILL_INPUT: if an argument had an illegal value @@ -522,7 +522,7 @@ Rootfinding optional input functions :retval ARK_SUCCESS: if successful :retval ARK_MEM_NULL: if the SPRKStep memory is ``NULL`` - :retval ARK_ILL_INPUT: if an argument has an illegal value + :retval ARK_ILL_INPUT: if an argument had an illegal value .. deprecated:: 6.1.0 @@ -809,7 +809,7 @@ Main solver optional output functions :param user_data: memory reference to a user data pointer :retval ARK_SUCCESS: if successful - :retval ARK_MEM_NULL: if the ARKStep memory was ``NULL`` + :retval ARK_MEM_NULL: if the SPRKStep memory was ``NULL`` .. deprecated:: 6.1.0 @@ -945,7 +945,7 @@ the RHS function should not incorporate the discontinuity. :retval ARK_SUCCESS: if successful :retval ARK_MEM_NULL: if the SPRKStep memory was ``NULL`` :retval ARK_MEM_FAIL: if a memory allocation failed - :retval ARK_ILL_INPUT: if an argument has an illegal value. + :retval ARK_ILL_INPUT: if an argument had an illegal value. .. _ARKODE.Usage.SPRKStep.Reset: @@ -971,7 +971,7 @@ SPRKStep reset function :retval ARK_SUCCESS: if successful :retval ARK_MEM_NULL: if the SPRKStep memory was ``NULL`` :retval ARK_MEM_FAIL: if a memory allocation failed - :retval ARK_ILL_INPUTL: if an argument has an illegal value. + :retval ARK_ILL_INPUTL: if an argument had an illegal value. .. note:: diff --git a/doc/arkode/guide/source/Usage/index.rst b/doc/arkode/guide/source/Usage/index.rst index 93b63d84b8..88cbfa7d3d 100644 --- a/doc/arkode/guide/source/Usage/index.rst +++ b/doc/arkode/guide/source/Usage/index.rst @@ -37,8 +37,8 @@ ARKODE's time stepping modules, including "relaxation" methods and preconitioners. Following our discussion of these commonalities, we separately discuss the usage details that that are specific to each of ARKODE's time stepping modules: :ref:`ARKStep `, -:ref:`ERKStep `, :ref:`SPRKStep ` -and :ref:`MRIStep `. +:ref:`ERKStep `, :ref:`SPRKStep `, +:ref:`MRIStep `, and :ref:`LSRKStep `. ARKODE also uses various input and output constants; these are defined as needed throughout this chapter, but for convenience the full list is provided @@ -77,5 +77,6 @@ ARKBBDPRE can only be used with NVECTOR_PARALLEL. Preconditioners ARKStep/index.rst ERKStep/index.rst + LSRKStep/index.rst SPRKStep/index.rst MRIStep/index.rst diff --git a/doc/shared/RecentChanges.rst b/doc/shared/RecentChanges.rst index 875d45c9b7..a614cb5827 100644 --- a/doc/shared/RecentChanges.rst +++ b/doc/shared/RecentChanges.rst @@ -8,6 +8,12 @@ only valid for Volta GPUs while the automatically selected value will vary across compilers and compiler versions. As such, users are encouraged to override this value with the architecture for their system. +Added a time-stepping module to ARKODE for low storage Runge--Kutta methods, +:ref:`LSRKStep `. This currently supports five explicit low-storage +methods: the second-order Runge--Kutta--Chebyshev and Runge--Kutta--Legendre methods, +and the second- through fourth-order optimal strong stability preserving Runge--Kutta methods. +All methods include embeddings for temporal adaptivity. + The Trilinos Tpetra NVector interface has been updated to utilize CMake imported targets added in Trilinos 14 to improve support for different Kokkos backends with Trilinos. As such, Trilinos 14 or newer is required and the diff --git a/doc/shared/sundials.bib b/doc/shared/sundials.bib index b000f2c3dc..62106338e6 100644 --- a/doc/shared/sundials.bib +++ b/doc/shared/sundials.bib @@ -2369,3 +2369,66 @@ @article{edwards2014kokkos issn = {0743-7315}, doi = {10.1016/j.jpdc.2014.07.003} } + +@article{VSH:04, +title = {{RKC} time-stepping for advection–diffusion–reaction problems}, +journal = {Journal of Computational Physics}, +volume = {201}, +number = {1}, +pages = {61-79}, +year = {2004}, +issn = {0021-9991}, +doi = {https://doi.org/10.1016/j.jcp.2004.05.002}, +author = {J.G. Verwer and B.P. Sommeijer and W. Hundsdorfer}} + +@article{MBA:14, +title = {A stabilized {Runge--Kutta--Legendre} method for explicit super-time-stepping of parabolic and mixed equations}, +journal = {Journal of Computational Physics}, +volume = {257}, +pages = {594-626}, +year = {2014}, +issn = {0021-9991}, +doi = {https://doi.org/10.1016/j.jcp.2013.08.021}, +author = {Chad D. Meyer and Dinshaw S. Balsara and Tariq D. Aslam}} + +@article{K:08, +title = {Highly Efficient Strong Stability-Preserving {Runge--Kutta} Methods with Low-Storage Implementations}, +journal = {SIAM Journal on Scientific Computing}, +volume = {30}, +number = {4}, +pages = {2113-2136}, +year = {2008}, +doi = {10.1137/07070485X}, +author = {Ketcheson, David I.}} + +@article{FCS:22, +title = {Embedded pairs for optimal explicit strong stability preserving {Runge--Kutta} methods}, +journal = {Journal of Computational and Applied Mathematics}, +volume = {412}, +pages = {114325}, +year = {2022}, +issn = {0377-0427}, +doi = {https://doi.org/10.1016/j.cam.2022.114325}, +author = {Imre Fekete and Sidafa Conde and John N. Shadid}} + +@article{SO:88, +title={Efficient implementation of essentially non-oscillatory shock-capturing schemes}, +journal={Journal of computational physics}, +volume={77}, +number={2}, +pages={439--471}, +year={1988}, +publisher={Elsevier}, +doi={10.1016/0021-9991(88)90177-5}, +author={Shu, Chi-Wang and Osher, Stanley}} + +@article{SR:02, +title={A new class of optimal high-order strong-stability-preserving time discretization methods}, +journal={SIAM Journal on Numerical Analysis}, +volume={40}, +number={2}, +pages={469--491}, +year={2002}, +publisher={SIAM}, +doi={10.1137/S0036142901389025}, +author={Spiteri, Raymond J and Ruuth, Steven J}} \ No newline at end of file diff --git a/examples/arkode/CMakeLists.txt b/examples/arkode/CMakeLists.txt index bfcbe267fd..611ac56cd8 100644 --- a/examples/arkode/CMakeLists.txt +++ b/examples/arkode/CMakeLists.txt @@ -42,6 +42,9 @@ endif() # C++ examples if(EXAMPLES_ENABLE_CXX) add_subdirectory(CXX_serial) + if(BUILD_NVECTOR_MANYVECTOR) + add_subdirectory(CXX_manyvector) + endif() if(ENABLE_MPI AND MPI_CXX_FOUND) add_subdirectory(CXX_parallel) endif() diff --git a/examples/arkode/CXX_manyvector/CMakeLists.txt b/examples/arkode/CXX_manyvector/CMakeLists.txt new file mode 100644 index 0000000000..037f1fd83b --- /dev/null +++ b/examples/arkode/CXX_manyvector/CMakeLists.txt @@ -0,0 +1,84 @@ +# --------------------------------------------------------------- +# Programmer(s): Daniel R. Reynolds @ SMU +# --------------------------------------------------------------- +# SUNDIALS Copyright Start +# Copyright (c) 2002-2024, Lawrence Livermore National Security +# and Southern Methodist University. +# All rights reserved. +# +# See the top-level LICENSE and NOTICE files for details. +# +# SPDX-License-Identifier: BSD-3-Clause +# SUNDIALS Copyright End +# --------------------------------------------------------------- +# CMakeLists.txt file for ARKODE C++ manyvector examples +# --------------------------------------------------------------- + +# Example lists are tuples "name\;args\;type" where the type is 'develop' for +# examples excluded from 'make test' in releases + +# Examples using SUNDIALS linear solvers +set(ARKODE_examples "ark_sod_lsrk.cpp\;\;exclude-single") + +# Header files to install +set(ARKODE_headers) + +# Auxiliary files to install +set(ARKODE_extras plot_sod.py) + +# Add the build and install targets for each example +foreach(example_tuple ${ARKODE_examples}) + + # parse the example tuple + list(GET example_tuple 0 example) + list(GET example_tuple 1 example_args) + list(GET example_tuple 2 example_type) + + # extract the file name without extension + get_filename_component(example_target ${example} NAME_WE) + + # add example executable + if(NOT TARGET ${example_target}) + add_executable(${example_target} ${example}) + + set_target_properties(${example_target} PROPERTIES FOLDER "Examples") + + # directories to include + target_include_directories(${example_target} + PRIVATE ${PROJECT_SOURCE_DIR}/examples/utilities) + + # libraries to link against + target_link_libraries(${example_target} sundials_arkode sundials_nvecserial + sundials_nvecmanyvector ${EXE_EXTRA_LINK_LIBS}) + + endif() + + # check if example args are provided and set the test name + if("${example_args}" STREQUAL "") + set(test_name ${example_target}) + else() + string(REGEX REPLACE " " "_" test_name ${example_target}_${example_args}) + endif() + + # add example to regression tests + sundials_add_test( + ${test_name} ${example_target} + TEST_ARGS ${example_args} + ANSWER_DIR ${CMAKE_CURRENT_SOURCE_DIR} + ANSWER_FILE ${test_name}.out + EXAMPLE_TYPE ${example_type}) + +endforeach(example_tuple ${ARKODE_examples}) + +if(EXAMPLES_INSTALL) + + sundials_install_examples( + arkode ARKODE_examples + CMAKE_TEMPLATE cmakelists_CXX_ex.in + MAKE_TEMPLATE makefile_serial_CXX_ex.in + SUNDIALS_TARGETS nvecserial nvecmanyvector arkode + DESTINATION arkode/CXX_manyvector + EXTRA_FILES ${ARKODE_extras} ${ARKODE_headers} + TEST_INSTALL CXX_manyvector) + +endif() diff --git a/examples/arkode/CXX_manyvector/ark_sod_lsrk.cpp b/examples/arkode/CXX_manyvector/ark_sod_lsrk.cpp new file mode 100644 index 0000000000..1bffa79cb0 --- /dev/null +++ b/examples/arkode/CXX_manyvector/ark_sod_lsrk.cpp @@ -0,0 +1,580 @@ +/* ----------------------------------------------------------------------------- + * Programmer(s): Daniel R. Reynolds @ SMU + * ----------------------------------------------------------------------------- + * SUNDIALS Copyright Start + * Copyright (c) 2002-2024, Lawrence Livermore National Security + * and Southern Methodist University. + * All rights reserved. + * + * See the top-level LICENSE and NOTICE files for details. + * + * SPDX-License-Identifier: BSD-3-Clause + * SUNDIALS Copyright End + * ----------------------------------------------------------------------------- + * This example performs a standard 1D sod shock tube (see section 5.2 of J.A. + * Greenough and W.J. Rider, "A quantitative comparison of numerical methods for + * the compressible Euler equations: fifth-order WENO and piecewise-linear + * Godunov," J. Comput. Phys., 196:259-281, 2004) + * [rho, vx, p] = { [1, 0, 1] if x < 0.5 + * { [0.125, 0, 0.1] if x > 0.5 + * + * The code solves the 1D compressible Euler equations in conserved variables, + * over the domain (t,x) in [0, 0.2] x [0, 1]. + * + * Since the Sod shock tube is specified in terms of primitive variables, we + * convert between primitive and conserved variables for the initial conditions + * and accuracy results. + * + * This problem should be run with homogeneous Neumann boundary conditions. + * + * The system is advanced in time using one of the strong-stability-preserving + * Runge--Kutta methods from LSRKStep, based on the --integrator METHOD input + * value. The following options are available: + * + * SSPRK(s,2) -- specified via METHOD = ARKODE_LSRK_SSP_S_2. The number of + * stages to use defaults to 10. + * + * SSPRK(s,3) -- specified via METHOD = ARKODE_LSRK_SSP_S_3. The number of + * stages to use defaults to 9. + * + * SSPRK(10,4) -- specified via METHOD = ARKODE_LSRK_SSP_10_4. + * + * Both the SSPRK(s,2) and SSPRK(s,3) methods allow specification of a + * non-default number of stages. This may be specified using the --stages S + * input value, where S is an integer. Note: SSPRK(s,2) requires S at least + * 2, and SSPRK(s,3) requires S be a perfect square, with S at least 4. + * + * Alternately, if METHOD corresponds with a valid ARKODE_ERKTableID then + * the system will be advanced using that method in ERKStep. + * + * Several additional command line options are available to change the + * and integrator settings. Use the flag --help for more information. + * ---------------------------------------------------------------------------*/ + +#include "ark_sod_lsrk.hpp" +using namespace std; + +int main(int argc, char* argv[]) +{ + // SUNDIALS context object for this simulation + sundials::Context ctx; + + // ----------------- + // Setup the problem + // ----------------- + + EulerData udata; + ARKODEParameters uopts; + + vector args(argv + 1, argv + argc); + + int flag = ReadInputs(args, udata, uopts, ctx); + if (check_flag(flag, "ReadInputs")) { return 1; } + if (flag > 0) { return 0; } + + flag = PrintSetup(udata, uopts); + if (check_flag(flag, "PrintSetup")) { return 1; } + + // Create state vector and set initial condition + N_Vector vecs[NSPECIES]; + for (int i = 0; i < NSPECIES; i++) + { + vecs[i] = N_VNew_Serial((sunindextype)udata.nx, ctx); // rho (density) + if (check_ptr(vecs[i], "N_VNew_Serial")) { return 1; } + } + N_Vector y = N_VNew_ManyVector(NSPECIES, vecs, ctx); + if (check_ptr(y, "N_VNew_ManyVector")) { return 1; } + + flag = SetIC(y, udata); + if (check_flag(flag, "SetIC")) { return 1; } + + // -------------------- + // Setup the integrator + // -------------------- + void* arkode_mem = nullptr; + + // Determine type (LSRKStep vs ERKStep) + bool lsrk = false; + if (uopts.integrator == "ARKODE_LSRK_SSP_S_2") { lsrk = true; } + if (uopts.integrator == "ARKODE_LSRK_SSP_S_3") { lsrk = true; } + if (uopts.integrator == "ARKODE_LSRK_SSP_10_4") { lsrk = true; } + + if (lsrk) // Setup LSRKStep + { + // ARKODE memory structure + arkode_mem = LSRKStepCreateSSP(frhs, udata.t0, y, ctx); + if (check_ptr(arkode_mem, "LSRKStepCreateSSP")) { return 1; } + + // Select SSPRK method type + flag = LSRKStepSetMethodByName(arkode_mem, uopts.integrator.c_str()); + if (check_flag(flag, "LSRKStepSetMethodByName")) { return 1; } + + // Select number of SSPRK stages + if (uopts.stages > 0) + { + flag = LSRKStepSetSSPStageNum(arkode_mem, uopts.stages); + if (check_flag(flag, "LSRKStepSetSSPStageNum")) { return 1; } + } + } + else + { // Setup ERKStep + + // ARKODE memory structure + arkode_mem = ERKStepCreate(frhs, udata.t0, y, ctx); + if (check_ptr(arkode_mem, "ERKStepCreate")) { return 1; } + + // Select ERK method + flag = ERKStepSetTableName(arkode_mem, uopts.integrator.c_str()); + if (check_flag(flag, "ERKStepSetTableName")) { return 1; } + } + + // Shared setup + + // Specify tolerances + flag = ARKodeSStolerances(arkode_mem, uopts.rtol, uopts.atol); + if (check_flag(flag, "ARKodeSStolerances")) { return 1; } + + // Attach user data + flag = ARKodeSetUserData(arkode_mem, &udata); + if (check_flag(flag, "ARKodeSetUserData")) { return 1; } + + // Set fixed step size or adaptivity method + if (uopts.fixed_h > ZERO) + { + flag = ARKodeSetFixedStep(arkode_mem, uopts.fixed_h); + if (check_flag(flag, "ARKodeSetFixedStep")) { return 1; } + } + + // Set max steps between outputs + flag = ARKodeSetMaxNumSteps(arkode_mem, uopts.maxsteps); + if (check_flag(flag, "ARKodeSetMaxNumSteps")) { return 1; } + + // Set stopping time + flag = ARKodeSetStopTime(arkode_mem, udata.tf); + if (check_flag(flag, "ARKodeSetStopTime")) { return 1; } + + // ---------------------- + // Evolve problem in time + // ---------------------- + + // Initial time, time between outputs, output time + sunrealtype t = ZERO; + sunrealtype dTout = udata.tf / uopts.nout; + sunrealtype tout = dTout; + + // initial output + flag = OpenOutput(udata, uopts); + if (check_flag(flag, "OpenOutput")) { return 1; } + + flag = WriteOutput(t, y, udata, uopts); + if (check_flag(flag, "WriteOutput")) { return 1; } + + // Loop over output times + for (int iout = 0; iout < uopts.nout; iout++) + { + // Evolve + if (uopts.output == 3) + { + // Stop at output time (do not interpolate output) + flag = ARKodeSetStopTime(arkode_mem, tout); + if (check_flag(flag, "ARKodeSetStopTime")) { return 1; } + } + + // Advance in time + flag = ARKodeEvolve(arkode_mem, tout, y, &t, ARK_NORMAL); + if (check_flag(flag, "ARKodeEvolve")) { break; } + + // Output solution + flag = WriteOutput(t, y, udata, uopts); + if (check_flag(flag, "WriteOutput")) { return 1; } + + // Update output time + tout += dTout; + tout = (tout > udata.tf) ? udata.tf : tout; + } + + // Close output + flag = CloseOutput(uopts); + if (check_flag(flag, "CloseOutput")) { return 1; } + + // ------------ + // Output stats + // ------------ + + if (uopts.output) + { + cout << "Final integrator statistics:" << endl; + flag = ARKodePrintAllStats(arkode_mem, stdout, SUN_OUTPUTFORMAT_TABLE); + } + + // -------- + // Clean up + // -------- + + ARKodeFree(&arkode_mem); + for (int i = 0; i < NSPECIES; i++) { N_VDestroy(vecs[i]); } + N_VDestroy(y); + + return 0; +} + +// ----------------------------------------------------------------------------- +// Functions called by the integrator +// ----------------------------------------------------------------------------- + +// ODE RHS function +int frhs(sunrealtype t, N_Vector y, N_Vector f, void* user_data) +{ + // Access problem data + EulerData* udata = (EulerData*)user_data; + + // initialize output to zeros + N_VConst(ZERO, f); + + // Access data arrays + sunrealtype* rho = N_VGetSubvectorArrayPointer_ManyVector(y, 0); + if (check_ptr(rho, "N_VGetSubvectorArrayPointer_ManyVector")) { return -1; } + sunrealtype* mx = N_VGetSubvectorArrayPointer_ManyVector(y, 1); + if (check_ptr(mx, "N_VGetSubvectorArrayPointer_ManyVector")) { return -1; } + sunrealtype* my = N_VGetSubvectorArrayPointer_ManyVector(y, 2); + if (check_ptr(my, "N_VGetSubvectorArrayPointer_ManyVector")) { return -1; } + sunrealtype* mz = N_VGetSubvectorArrayPointer_ManyVector(y, 3); + if (check_ptr(mz, "N_VGetSubvectorArrayPointer_ManyVector")) { return -1; } + sunrealtype* et = N_VGetSubvectorArrayPointer_ManyVector(y, 4); + if (check_ptr(et, "N_VGetSubvectorArrayPointer_ManyVector")) { return -1; } + + sunrealtype* rhodot = N_VGetSubvectorArrayPointer_ManyVector(f, 0); + if (check_ptr(rhodot, "N_VGetSubvectorArrayPointer_ManyVector")) + { + return -1; + } + sunrealtype* mxdot = N_VGetSubvectorArrayPointer_ManyVector(f, 1); + if (check_ptr(mxdot, "N_VGetSubvectorArrayPointer_ManyVector")) { return -1; } + sunrealtype* mydot = N_VGetSubvectorArrayPointer_ManyVector(f, 2); + if (check_ptr(mydot, "N_VGetSubvectorArrayPointer_ManyVector")) { return -1; } + sunrealtype* mzdot = N_VGetSubvectorArrayPointer_ManyVector(f, 3); + if (check_ptr(mzdot, "N_VGetSubvectorArrayPointer_ManyVector")) { return -1; } + sunrealtype* etdot = N_VGetSubvectorArrayPointer_ManyVector(f, 4); + if (check_ptr(etdot, "N_VGetSubvectorArrayPointer_ManyVector")) { return -1; } + + // Set shortcut variables + const long int nx = udata->nx; + const sunrealtype dx = udata->dx; + sunrealtype* flux = udata->flux; + + // compute face-centered fluxes over domain interior: pack 1D x-directional array + // of variable shortcuts, and compute flux at lower x-directional face + for (long int i = 3; i < nx - 2; i++) + { + udata->pack1D(rho, mx, my, mz, et, i); + face_flux(udata->w1d, &(flux[i * NSPECIES]), *udata); + } + + // compute face-centered fluxes at left boundary + for (long int i = 0; i < 3; i++) + { + udata->pack1D_bdry(rho, mx, my, mz, et, i); + face_flux(udata->w1d, &(flux[i * NSPECIES]), *udata); + } + + // compute face-centered fluxes at right boundary + for (long int i = nx - 2; i <= nx; i++) + { + udata->pack1D_bdry(rho, mx, my, mz, et, i); + face_flux(udata->w1d, &(flux[i * NSPECIES]), *udata); + } + + // iterate over subdomain, updating RHS + for (long int i = 0; i < nx; i++) + { + rhodot[i] -= (flux[(i + 1) * NSPECIES + 0] - flux[i * NSPECIES + 0]) / dx; + mxdot[i] -= (flux[(i + 1) * NSPECIES + 1] - flux[i * NSPECIES + 1]) / dx; + mydot[i] -= (flux[(i + 1) * NSPECIES + 2] - flux[i * NSPECIES + 2]) / dx; + mzdot[i] -= (flux[(i + 1) * NSPECIES + 3] - flux[i * NSPECIES + 3]) / dx; + etdot[i] -= (flux[(i + 1) * NSPECIES + 4] - flux[i * NSPECIES + 4]) / dx; + } + + return 0; +} + +// given a 6-point stencil of solution values, +// w(x_{j-3}) w(x_{j-2}) w(x_{j-1}), w(x_j), w(x_{j+1}), w(x_{j+2}), +// compute the face-centered flux (flux) at the center of the stencil, x_{j-1/2}. +// +// This precisely follows the recipe laid out in: +// Chi-Wang Shu (2003) "High-order Finite Difference and Finite Volume WENO +// Schemes and Discontinuous Galerkin Methods for CFD," International Journal of +// Computational Fluid Dynamics, 17:2, 107-118, DOI: 10.1080/1061856031000104851 +// with the only change that since this is 1D, we manually set the y- and +// z-velocities, v and w, to zero. +void face_flux(sunrealtype (&w1d)[STSIZE][NSPECIES], sunrealtype* f_face, + const EulerData& udata) +{ + // local data + int i, j; + sunrealtype rhosqrL, rhosqrR, rhosqrbar, u, v, w, H, qsq, csnd, cinv, gamm, + alpha, beta1, beta2, beta3, w1, w2, w3, f1, f2, f3; + sunrealtype RV[5][5], LV[5][5], p[STSIZE], flux[STSIZE][NSPECIES], + fproj[5][NSPECIES], fs[5][NSPECIES], ff[NSPECIES]; + const sunrealtype bc = + SUN_RCONST(1.083333333333333333333333333333333333333); // 13/12 + const sunrealtype epsilon = SUN_RCONST(1.0e-6); + + // compute pressures over stencil + for (i = 0; i < STSIZE; i++) + { + p[i] = udata.eos(w1d[i][0], w1d[i][1], w1d[i][2], w1d[i][3], w1d[i][4]); + } + + // compute Roe-average state at face: + // wbar = [sqrt(rho), sqrt(rho)*vx, sqrt(rho)*vy, sqrt(rho)*vz, (e+p)/sqrt(rho)] + // [sqrt(rho), mx/sqrt(rho), my/sqrt(rho), mz/sqrt(rho), (e+p)/sqrt(rho)] + // u = wbar_2 / wbar_1 + // v = wbar_3 / wbar_1 + // w = wbar_4 / wbar_1 + // H = wbar_5 / wbar_1 + rhosqrL = sqrt(w1d[2][0]); + rhosqrR = sqrt(w1d[3][0]); + rhosqrbar = HALF * (rhosqrL + rhosqrR); + u = HALF * (w1d[2][1] / rhosqrL + w1d[3][1] / rhosqrR) / rhosqrbar; + v = HALF * (w1d[2][2] / rhosqrL + w1d[3][2] / rhosqrR) / rhosqrbar; + w = HALF * (w1d[2][3] / rhosqrL + w1d[3][3] / rhosqrR) / rhosqrbar; + H = HALF * ((p[2] + w1d[2][4]) / rhosqrL + (p[3] + w1d[3][4]) / rhosqrR) / + rhosqrbar; + + // compute eigenvectors at face (note: eigenvectors for tracers are just identity) + qsq = u * u + v * v + w * w; + gamm = udata.gamma - ONE; + csnd = gamm * (H - HALF * qsq); + cinv = ONE / csnd; + for (i = 0; i < 5; i++) + { + for (j = 0; j < 5; j++) + { + RV[i][j] = ZERO; + LV[i][j] = ZERO; + } + } + + RV[0][0] = ONE; + RV[0][3] = ONE; + RV[0][4] = ONE; + + RV[1][0] = u - csnd; + RV[1][3] = u; + RV[1][4] = u + csnd; + + RV[2][0] = v; + RV[2][1] = ONE; + RV[2][3] = v; + RV[2][4] = v; + + RV[3][0] = w; + RV[3][2] = ONE; + RV[3][3] = w; + RV[3][4] = w; + + RV[4][0] = H - u * csnd; + RV[4][1] = v; + RV[4][2] = w; + RV[4][3] = HALF * qsq; + RV[4][4] = H + u * csnd; + + LV[0][0] = HALF * cinv * (u + HALF * gamm * qsq); + LV[0][1] = -HALF * cinv * (gamm * u + ONE); + LV[0][2] = -HALF * v * gamm * cinv; + LV[0][3] = -HALF * w * gamm * cinv; + LV[0][4] = HALF * gamm * cinv; + + LV[1][0] = -v; + LV[1][2] = ONE; + + LV[2][0] = -w; + LV[2][3] = ONE; + + LV[3][0] = -gamm * cinv * (qsq - H); + LV[3][1] = u * gamm * cinv; + LV[3][2] = v * gamm * cinv; + LV[3][3] = w * gamm * cinv; + LV[3][4] = -gamm * cinv; + + LV[4][0] = -HALF * cinv * (u - HALF * gamm * qsq); + LV[4][1] = -HALF * cinv * (gamm * u - ONE); + LV[4][2] = -HALF * v * gamm * cinv; + LV[4][3] = -HALF * w * gamm * cinv; + LV[4][4] = HALF * gamm * cinv; + + // compute fluxes and max wave speed over stencil + alpha = ZERO; + for (j = 0; j < STSIZE; j++) + { + u = w1d[j][1] / w1d[j][0]; // u = vx = mx/rho + flux[j][0] = w1d[j][1]; // f_rho = rho*u = mx + flux[j][1] = u * w1d[j][1] + p[j]; // f_mx = rho*u*u + p = mx*u + p + flux[j][2] = u * w1d[j][2]; // f_my = rho*v*u = my*u + flux[j][3] = u * w1d[j][3]; // f_mz = rho*w*u = mz*u + flux[j][4] = u * (w1d[j][4] + p[j]); // f_et = u*(et + p) + csnd = sqrt(udata.gamma * p[j] / w1d[j][0]); // csnd = sqrt(gamma*p/rho) + alpha = max(alpha, abs(u) + csnd); + } + + // compute flux from right side of face at x_{i+1/2}: + + // compute right-shifted Lax-Friedrichs flux over left portion of patch + for (j = 0; j < 5; j++) + { + for (i = 0; i < NSPECIES; i++) + { + fs[j][i] = HALF * (flux[j][i] + alpha * w1d[j][i]); + } + } + + // compute projected flux for fluid fields + for (j = 0; j < 5; j++) + { + for (i = 0; i < 5; i++) + { + fproj[j][i] = LV[i][0] * fs[j][0] + LV[i][1] * fs[j][1] + + LV[i][2] * fs[j][2] + LV[i][3] * fs[j][3] + + LV[i][4] * fs[j][4]; + } + } + + // compute WENO signed flux + for (i = 0; i < NSPECIES; i++) + { + // smoothness indicators + beta1 = bc * pow(fproj[2][i] - SUN_RCONST(2.0) * fproj[3][i] + fproj[4][i], + 2) + + FOURTH * pow(SUN_RCONST(3.0) * fproj[2][i] - + SUN_RCONST(4.0) * fproj[3][i] + fproj[4][i], + 2); + beta2 = bc * pow(fproj[1][i] - SUN_RCONST(2.0) * fproj[2][i] + fproj[3][i], + 2) + + FOURTH * pow(fproj[1][i] - fproj[3][i], 2); + beta3 = bc * pow(fproj[0][i] - SUN_RCONST(2.0) * fproj[1][i] + fproj[2][i], + 2) + + FOURTH * pow(fproj[0][i] - SUN_RCONST(4.0) * fproj[1][i] + + SUN_RCONST(3.0) * fproj[2][i], + 2); + // nonlinear weights + w1 = SUN_RCONST(0.3) / ((epsilon + beta1) * (epsilon + beta1)); + w2 = SUN_RCONST(0.6) / ((epsilon + beta2) * (epsilon + beta2)); + w3 = SUN_RCONST(0.1) / ((epsilon + beta3) * (epsilon + beta3)); + // flux stencils + f1 = SUN_RCONST(0.3333333333333333333333333333333333333333) * fproj[2][i] + + SUN_RCONST(0.8333333333333333333333333333333333333333) * fproj[3][i] - + SUN_RCONST(0.1666666666666666666666666666666666666667) * fproj[4][i]; + f2 = -SUN_RCONST(0.1666666666666666666666666666666666666667) * fproj[1][i] + + SUN_RCONST(0.8333333333333333333333333333333333333333) * fproj[2][i] + + SUN_RCONST(0.3333333333333333333333333333333333333333) * fproj[3][i]; + f3 = SUN_RCONST(0.3333333333333333333333333333333333333333) * fproj[0][i] - + SUN_RCONST(1.166666666666666666666666666666666666667) * fproj[1][i] + + SUN_RCONST(1.833333333333333333333333333333333333333) * fproj[2][i]; + // resulting signed flux at face + ff[i] = (f1 * w1 + f2 * w2 + f3 * w3) / (w1 + w2 + w3); + } + + // compute flux from left side of face at x_{i+1/2}: + + // compute left-shifted Lax-Friedrichs flux over right portion of patch + for (j = 0; j < 5; j++) + { + for (i = 0; i < NSPECIES; i++) + { + fs[j][i] = HALF * (flux[j + 1][i] - alpha * w1d[j + 1][i]); + } + } + + // compute projected flux for fluid fields + for (j = 0; j < 5; j++) + { + for (i = 0; i < 5; i++) + { + fproj[j][i] = LV[i][0] * fs[j][0] + LV[i][1] * fs[j][1] + + LV[i][2] * fs[j][2] + LV[i][3] * fs[j][3] + + LV[i][4] * fs[j][4]; + } + } + + // compute WENO signed fluxes + for (i = 0; i < NSPECIES; i++) + { + // smoothness indicators + beta1 = bc * pow(fproj[2][i] - SUN_RCONST(2.0) * fproj[3][i] + fproj[4][i], + 2) + + FOURTH * pow(SUN_RCONST(3.0) * fproj[2][i] - + SUN_RCONST(4.0) * fproj[3][i] + fproj[4][i], + 2); + beta2 = bc * pow(fproj[1][i] - SUN_RCONST(2.0) * fproj[2][i] + fproj[3][i], + 2) + + FOURTH * pow(fproj[1][i] - fproj[3][i], 2); + beta3 = bc * pow(fproj[0][i] - SUN_RCONST(2.0) * fproj[1][i] + fproj[2][i], + 2) + + FOURTH * pow(fproj[0][i] - SUN_RCONST(4.0) * fproj[1][i] + + SUN_RCONST(3.0) * fproj[2][i], + 2); + // nonlinear weights + w1 = SUN_RCONST(0.1) / ((epsilon + beta1) * (epsilon + beta1)); + w2 = SUN_RCONST(0.6) / ((epsilon + beta2) * (epsilon + beta2)); + w3 = SUN_RCONST(0.3) / ((epsilon + beta3) * (epsilon + beta3)); + // flux stencils + f1 = SUN_RCONST(1.833333333333333333333333333333333333333) * fproj[2][i] - + SUN_RCONST(1.166666666666666666666666666666666666667) * fproj[3][i] + + SUN_RCONST(0.3333333333333333333333333333333333333333) * fproj[4][i]; + f2 = SUN_RCONST(0.3333333333333333333333333333333333333333) * fproj[1][i] + + SUN_RCONST(0.8333333333333333333333333333333333333333) * fproj[2][i] - + SUN_RCONST(0.1666666666666666666666666666666666666667) * fproj[3][i]; + f3 = -SUN_RCONST(0.1666666666666666666666666666666666666667) * fproj[0][i] + + SUN_RCONST(0.8333333333333333333333333333333333333333) * fproj[1][i] + + SUN_RCONST(0.3333333333333333333333333333333333333333) * fproj[2][i]; + // resulting signed flux (add to ff) + ff[i] += (f1 * w1 + f2 * w2 + f3 * w3) / (w1 + w2 + w3); + } + + // combine signed fluxes into output, converting back to conserved variables + for (i = 0; i < NSPECIES; i++) + { + f_face[i] = RV[i][0] * ff[0] + RV[i][1] * ff[1] + RV[i][2] * ff[2] + + RV[i][3] * ff[3] + RV[i][4] * ff[4]; + } + return; +} + +// Compute the initial condition +int SetIC(N_Vector y, EulerData& udata) +{ + sunrealtype* rho = N_VGetSubvectorArrayPointer_ManyVector(y, 0); + if (check_ptr(rho, "N_VGetSubvectorArrayPointer_ManyVector")) { return -1; } + sunrealtype* mx = N_VGetSubvectorArrayPointer_ManyVector(y, 1); + if (check_ptr(mx, "N_VGetSubvectorArrayPointer_ManyVector")) { return -1; } + sunrealtype* my = N_VGetSubvectorArrayPointer_ManyVector(y, 2); + if (check_ptr(my, "N_VGetSubvectorArrayPointer_ManyVector")) { return -1; } + sunrealtype* mz = N_VGetSubvectorArrayPointer_ManyVector(y, 3); + if (check_ptr(mz, "N_VGetSubvectorArrayPointer_ManyVector")) { return -1; } + sunrealtype* et = N_VGetSubvectorArrayPointer_ManyVector(y, 4); + if (check_ptr(et, "N_VGetSubvectorArrayPointer_ManyVector")) { return -1; } + + for (long int i = 0; i < udata.nx; i++) + { + sunrealtype xloc = ((sunrealtype)i + HALF) * udata.dx + udata.xl; + if (xloc < HALF) + { + rho[i] = rhoL; + et[i] = udata.eos_inv(rhoL, uL, ZERO, ZERO, pL); + mx[i] = rhoL * uL; + } + else + { + rho[i] = rhoR; + et[i] = udata.eos_inv(rhoR, uR, ZERO, ZERO, pR); + mx[i] = rhoR * uR; + } + my[i] = ZERO; + mz[i] = ZERO; + } + + return 0; +} + +//---- end of file ---- diff --git a/examples/arkode/CXX_manyvector/ark_sod_lsrk.hpp b/examples/arkode/CXX_manyvector/ark_sod_lsrk.hpp new file mode 100644 index 0000000000..0ccd7959df --- /dev/null +++ b/examples/arkode/CXX_manyvector/ark_sod_lsrk.hpp @@ -0,0 +1,538 @@ +/* ----------------------------------------------------------------------------- + * Programmer(s): Daniel R. Reynolds @ SMU + * ----------------------------------------------------------------------------- + * SUNDIALS Copyright Start + * Copyright (c) 2002-2024, Lawrence Livermore National Security + * and Southern Methodist University. + * All rights reserved. + * + * See the top-level LICENSE and NOTICE files for details. + * + * SPDX-License-Identifier: BSD-3-Clause + * SUNDIALS Copyright End + * ----------------------------------------------------------------------------- + * Header file for ARKODE LSRKStep Sod shock tube example, see + * ark_sod_lsrk.cpp for more details. + * ---------------------------------------------------------------------------*/ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// Include desired integrators, vectors, linear solvers, and nonlinear solvers +#include "arkode/arkode_erkstep.h" +#include "arkode/arkode_lsrkstep.h" +#include "nvector/nvector_manyvector.h" +#include "nvector/nvector_serial.h" +#include "sundials/sundials_core.hpp" + +// Macros for problem constants +#define rhoL SUN_RCONST(1.0) +#define rhoR SUN_RCONST(0.125) +#define pL SUN_RCONST(1.0) +#define pR SUN_RCONST(0.1) +#define uL SUN_RCONST(0.0) +#define uR SUN_RCONST(0.0) +#define HALF SUN_RCONST(0.5) +#define ZERO SUN_RCONST(0.0) +#define ONE SUN_RCONST(1.0) +#define TWO SUN_RCONST(2.0) +#define FOURTH SUN_RCONST(0.25) + +#define NSPECIES 5 +#define STSIZE 6 + +#define WIDTH (10 + std::numeric_limits::digits10) + +// ----------------------------------------------------------------------------- +// Problem options +// ----------------------------------------------------------------------------- + +class ARKODEParameters +{ +public: + // Integration method (ARKODE_LSRK_SSP_S_2, ARKODE_LSRK_SSP_S_3, ARKODE_LSRK_SSP_10_4, + // or any valid ARKODE_ERKTableID for ERK methods) + std::string integrator; + + // Method stages (0 => to use the default; ignored if using ARKODE_LSRK_SSP_10_4 or + // an ERK method) + int stages; + + // Relative and absolute tolerances + sunrealtype rtol; + sunrealtype atol; + + // Step size selection (ZERO = adaptive steps) + sunrealtype fixed_h; + + // Maximum number of time steps between outputs + int maxsteps; + + // Output-related information + int output; // 0 = none, 1 = stats, 2 = disk, 3 = disk with tstop + int nout; // number of output times + std::ofstream uout; // output file stream + + // constructor (with default values) + ARKODEParameters() + : integrator("ARKODE_LSRK_SSP_10_4"), + stages(0), + rtol(SUN_RCONST(1.e-4)), + atol(SUN_RCONST(1.e-11)), + fixed_h(ZERO), + maxsteps(10000), + output(1), + nout(10){}; + +}; // end ARKODEParameters + +// ----------------------------------------------------------------------------- +// Problem parameters +// ----------------------------------------------------------------------------- + +// user data class +class EulerData +{ +public: + ///// domain related data ///// + long int nx; // global number of x grid points + sunrealtype t0; // time domain extents + sunrealtype tf; + sunrealtype xl; // spatial domain extents + sunrealtype xr; + sunrealtype dx; // spatial mesh spacing + + ///// problem-defining data ///// + sunrealtype gamma; // ratio of specific heat capacities, cp/cv + + ///// reusable arrays for WENO flux calculations ///// + sunrealtype* flux; + sunrealtype w1d[STSIZE][NSPECIES]; + + ///// class operations ///// + + // constructor + EulerData() + : nx(512), + t0(ZERO), + tf(SUN_RCONST(0.1)), + xl(ZERO), + xr(ONE), + dx(ZERO), + gamma(SUN_RCONST(1.4)), + flux(nullptr){}; + + // manual destructor + void FreeData() + { + delete[] flux; + flux = nullptr; + }; + + // destructor + ~EulerData() { this->FreeData(); }; + + // Utility routine to pack 1-dimensional data for *interior only* data; + // e.g., in the x-direction given a location (i), we return values at + // the 6 nodal values closest to the (i-1/2) face along the x-direction, + // {w(i-3), w(i-2), w(i-1), w(i), w(i+1), w(i+2)}. + inline void pack1D(const sunrealtype* rho, const sunrealtype* mx, + const sunrealtype* my, const sunrealtype* mz, + const sunrealtype* et, const long int& i) + { + for (int l = 0; l < STSIZE; l++) { this->w1d[l][0] = rho[i - 3 + l]; } + for (int l = 0; l < STSIZE; l++) { this->w1d[l][1] = mx[i - 3 + l]; } + for (int l = 0; l < STSIZE; l++) { this->w1d[l][2] = my[i - 3 + l]; } + for (int l = 0; l < STSIZE; l++) { this->w1d[l][3] = mz[i - 3 + l]; } + for (int l = 0; l < STSIZE; l++) { this->w1d[l][4] = et[i - 3 + l]; } + } + + // Utility routine to pack 1-dimensional data for locations near the + // boundary; like the routine above this packs the 6 closest + // entries aligned with, e.g., the (i-1/2) face, but now some entries + // are set to satisfy homogeneous Neumann boundary conditions. + inline void pack1D_bdry(const sunrealtype* rho, const sunrealtype* mx, + const sunrealtype* my, const sunrealtype* mz, + const sunrealtype* et, const long int& i) + { + for (int l = 0; l < 3; l++) + { + this->w1d[l][0] = (i < (3 - l)) ? rho[2 - (i + l)] : rho[i - 3 + l]; + } + for (int l = 0; l < 3; l++) + { + this->w1d[l][1] = (i < (3 - l)) ? mx[2 - (i + l)] : mx[i - 3 + l]; + } + for (int l = 0; l < 3; l++) + { + this->w1d[l][2] = (i < (3 - l)) ? my[2 - (i + l)] : my[i - 3 + l]; + } + for (int l = 0; l < 3; l++) + { + this->w1d[l][3] = (i < (3 - l)) ? mz[2 - (i + l)] : mz[i - 3 + l]; + } + for (int l = 0; l < 3; l++) + { + this->w1d[l][4] = (i < (3 - l)) ? et[2 - (i + l)] : et[i - 3 + l]; + } + for (int l = 0; l < 3; l++) + { + this->w1d[l + 3][0] = (i > (nx - l - 1)) ? rho[i + l - 3] : rho[i + l]; + } + for (int l = 0; l < 3; l++) + { + this->w1d[l + 3][1] = (i > (nx - l - 1)) ? mx[i + l - 3] : mx[i + l]; + } + for (int l = 0; l < 3; l++) + { + this->w1d[l + 3][2] = (i > (nx - l - 1)) ? my[i + l - 3] : my[i + l]; + } + for (int l = 0; l < 3; l++) + { + this->w1d[l + 3][3] = (i > (nx - l - 1)) ? mz[i + l - 3] : mz[i + l]; + } + for (int l = 0; l < 3; l++) + { + this->w1d[l + 3][4] = (i > (nx - l - 1)) ? et[i + l - 3] : et[i + l]; + } + } + + // Equation of state -- compute and return pressure, + // p = (gamma-1)*(e - rho/2*(vx^2+vy^2+vz^2)), or equivalently + // p = (gamma-1)*(e - (mx^2+my^2+mz^2)/(2*rho)) + inline sunrealtype eos(const sunrealtype& rho, const sunrealtype& mx, + const sunrealtype& my, const sunrealtype& mz, + const sunrealtype& et) const + { + return ((gamma - ONE) * (et - (mx * mx + my * my + mz * mz) * HALF / rho)); + } + + // Equation of state inverse -- compute and return energy, + // e_t = p/(gamma-1) + rho/2*(vx^2+vy^2+vz^2), or equivalently + // e_t = p/(gamma-1) + (mx^2+my^2+mz^2)/(2*rho) + inline sunrealtype eos_inv(const sunrealtype& rho, const sunrealtype& mx, + const sunrealtype& my, const sunrealtype& mz, + const sunrealtype& pr) const + { + return (pr / (gamma - ONE) + (mx * mx + my * my + mz * mz) * HALF / rho); + } + +}; // end EulerData; + +// ----------------------------------------------------------------------------- +// Functions provided to the SUNDIALS integrators +// ----------------------------------------------------------------------------- + +// ODE right hand side (RHS) functions +int frhs(sunrealtype t, N_Vector y, N_Vector f, void* user_data); + +// ----------------------------------------------------------------------------- +// Helper functions +// ----------------------------------------------------------------------------- + +// WENO flux calculation helper function +void face_flux(sunrealtype (&w1d)[6][NSPECIES], sunrealtype* f_face, + const EulerData& udata); + +// Compute the initial condition +int SetIC(N_Vector y, EulerData& udata); + +// ----------------------------------------------------------------------------- +// Output and utility functions +// ----------------------------------------------------------------------------- + +// Check function return flag +static int check_flag(int flag, const std::string funcname) +{ + if (flag < 0) + { + std::cerr << "ERROR: " << funcname << " returned " << flag << std::endl; + return 1; + } + return 0; +} + +// Check if a function returned a NULL pointer +static int check_ptr(void* ptr, const std::string funcname) +{ + if (ptr) { return 0; } + std::cerr << "ERROR: " << funcname << " returned NULL" << std::endl; + return 1; +} + +// Print command line options +static void InputHelp() +{ + std::cout << std::endl; + std::cout << "Command line options:" << std::endl; + std::cout << " --integrator : method (ARKODE_LSRK_SSP_S_2, " + "ARKODE_LSRK_SSP_S_3, " + "ARKODE_LSRK_SSP_10_4, or any valid ARKODE_ERKTableID)\n"; + std::cout << " --stages : number of stages (ignored for " + "ARKODE_LSRK_SSP_10_4 and ERK)\n"; + std::cout << " --tf : final time\n"; + std::cout << " --xl : domain lower boundary\n"; + std::cout << " --xr : domain upper boundary\n"; + std::cout << " --gamma : ideal gas constant\n"; + std::cout << " --nx : number of mesh points\n"; + std::cout << " --rtol : relative tolerance\n"; + std::cout << " --atol : absolute tolerance\n"; + std::cout << " --fixed_h : fixed step size\n"; + std::cout << " --maxsteps : max steps between outputs\n"; + std::cout << " --output : output level\n"; + std::cout << " --nout : number of outputs\n"; + std::cout << " --help : print options and exit\n"; +} + +inline void find_arg(std::vector& args, const std::string key, + sunrealtype& dest) +{ + auto it = find(args.begin(), args.end(), key); + if (it != args.end()) + { +#if defined(SUNDIALS_SINGLE_PRECISION) + dest = stof(*(it + 1)); +#elif defined(SUNDIALS_DOUBLE_PRECISION) + dest = stod(*(it + 1)); +#elif defined(SUNDIALS_EXTENDED_PRECISION) + dest = stold(*(it + 1)); +#endif + args.erase(it, it + 2); + } +} + +inline void find_arg(std::vector& args, const std::string key, + long int& dest) +{ + auto it = find(args.begin(), args.end(), key); + if (it != args.end()) + { + dest = stoll(*(it + 1)); + args.erase(it, it + 2); + } +} + +inline void find_arg(std::vector& args, const std::string key, + int& dest) +{ + auto it = find(args.begin(), args.end(), key); + if (it != args.end()) + { + dest = stoi(*(it + 1)); + args.erase(it, it + 2); + } +} + +inline void find_arg(std::vector& args, const std::string key, + std::string& dest) +{ + auto it = find(args.begin(), args.end(), key); + if (it != args.end()) + { + dest = *(it + 1); + args.erase(it, it + 2); + } +} + +inline void find_arg(std::vector& args, const std::string key, + bool& dest, bool store = true) +{ + auto it = find(args.begin(), args.end(), key); + if (it != args.end()) + { + dest = store; + args.erase(it); + } +} + +static int ReadInputs(std::vector& args, EulerData& udata, + ARKODEParameters& uopts, SUNContext ctx) +{ + if (find(args.begin(), args.end(), "--help") != args.end()) + { + InputHelp(); + return 1; + } + + // Problem parameters + find_arg(args, "--gamma", udata.gamma); + find_arg(args, "--tf", udata.tf); + find_arg(args, "--xl", udata.xl); + find_arg(args, "--xr", udata.xr); + find_arg(args, "--nx", udata.nx); + + // Integrator options + find_arg(args, "--integrator", uopts.integrator); + find_arg(args, "--stages", uopts.stages); + find_arg(args, "--rtol", uopts.rtol); + find_arg(args, "--atol", uopts.atol); + find_arg(args, "--fixed_h", uopts.fixed_h); + find_arg(args, "--maxsteps", uopts.maxsteps); + find_arg(args, "--output", uopts.output); + find_arg(args, "--nout", uopts.nout); + + // Recompute mesh spacing and [re]allocate flux array + udata.dx = (udata.xr - udata.xl) / ((sunrealtype)udata.nx); + if (udata.flux) { delete[] udata.flux; } + udata.flux = new sunrealtype[NSPECIES * (udata.nx + 1)]; + + if (uopts.stages < 0) + { + std::cerr << "ERROR: Invalid number of stages" << std::endl; + return -1; + } + + return 0; +} + +// Print user data +static int PrintSetup(EulerData& udata, ARKODEParameters& uopts) +{ + std::cout << std::endl; + std::cout << "Problem parameters and options:" << std::endl; + std::cout << " --------------------------------- " << std::endl; + std::cout << " gamma = " << udata.gamma << std::endl; + std::cout << " --------------------------------- " << std::endl; + std::cout << " tf = " << udata.tf << std::endl; + std::cout << " xl = " << udata.xl << std::endl; + std::cout << " xr = " << udata.xr << std::endl; + std::cout << " nx = " << udata.nx << std::endl; + std::cout << " dx = " << udata.dx << std::endl; + std::cout << " --------------------------------- " << std::endl; + std::cout << " integrator = " << uopts.integrator << std::endl; + if (uopts.stages > 0) + { + std::cout << " stages = " << uopts.stages << std::endl; + } + std::cout << " rtol = " << uopts.rtol << std::endl; + std::cout << " atol = " << uopts.atol << std::endl; + std::cout << " fixed h = " << uopts.fixed_h << std::endl; + std::cout << " --------------------------------- " << std::endl; + std::cout << " output = " << uopts.output << std::endl; + std::cout << " --------------------------------- " << std::endl; + std::cout << std::endl; + + return 0; +} + +// Initialize output +static int OpenOutput(EulerData& udata, ARKODEParameters& uopts) +{ + // Header for status output + if (uopts.output) + { + std::cout << std::scientific; + std::cout << std::setprecision(std::numeric_limits::digits10); + std::cout << " t " + << " ||rho|| " + << " ||mx|| " + << " ||my|| " + << " ||mz|| " + << " ||et||" << std::endl; + std::cout + << " -----------------------------------------------------------------" + "---------" + << std::endl; + } + + // Open output stream and output problem information + if (uopts.output >= 2) + { + // Open output stream + std::stringstream fname; + fname << "sod.out"; + uopts.uout.open(fname.str()); + + uopts.uout << std::scientific; + uopts.uout << std::setprecision(std::numeric_limits::digits10); + uopts.uout << "# title Sod Shock Tube" << std::endl; + uopts.uout << "# nvar 5" << std::endl; + uopts.uout << "# vars rho mx my mz et" << std::endl; + uopts.uout << "# nt " << uopts.nout + 1 << std::endl; + uopts.uout << "# nx " << udata.nx << std::endl; + uopts.uout << "# xl " << udata.xl << std::endl; + uopts.uout << "# xr " << udata.xr << std::endl; + } + + return 0; +} + +// Write output +static int WriteOutput(sunrealtype t, N_Vector y, EulerData& udata, + ARKODEParameters& uopts) +{ + if (uopts.output) + { + // Compute rms norm of the state + N_Vector rho = N_VGetSubvector_ManyVector(y, 0); + N_Vector mx = N_VGetSubvector_ManyVector(y, 1); + N_Vector my = N_VGetSubvector_ManyVector(y, 2); + N_Vector mz = N_VGetSubvector_ManyVector(y, 3); + N_Vector et = N_VGetSubvector_ManyVector(y, 4); + sunrealtype rhorms = sqrt(N_VDotProd(rho, rho) / (sunrealtype)udata.nx); + sunrealtype mxrms = sqrt(N_VDotProd(mx, mx) / (sunrealtype)udata.nx); + sunrealtype myrms = sqrt(N_VDotProd(my, my) / (sunrealtype)udata.nx); + sunrealtype mzrms = sqrt(N_VDotProd(mz, mz) / (sunrealtype)udata.nx); + sunrealtype etrms = sqrt(N_VDotProd(et, et) / (sunrealtype)udata.nx); + std::cout << std::setprecision(2) << " " << t << std::setprecision(5) + << " " << rhorms << " " << mxrms << " " << myrms << " " + << mzrms << " " << etrms << std::endl; + + // Write solution to disk + if (uopts.output >= 2) + { + sunrealtype* rhodata = N_VGetArrayPointer(rho); + if (check_ptr(rhodata, "N_VGetArrayPointer")) { return -1; } + sunrealtype* mxdata = N_VGetArrayPointer(mx); + if (check_ptr(mxdata, "N_VGetArrayPointer")) { return -1; } + sunrealtype* mydata = N_VGetArrayPointer(my); + if (check_ptr(mydata, "N_VGetArrayPointer")) { return -1; } + sunrealtype* mzdata = N_VGetArrayPointer(mz); + if (check_ptr(mzdata, "N_VGetArrayPointer")) { return -1; } + sunrealtype* etdata = N_VGetArrayPointer(et); + if (check_ptr(etdata, "N_VGetArrayPointer")) { return -1; } + + uopts.uout << t; + for (sunindextype i = 0; i < udata.nx; i++) + { + uopts.uout << std::setw(WIDTH) << rhodata[i]; + uopts.uout << std::setw(WIDTH) << mxdata[i]; + uopts.uout << std::setw(WIDTH) << mydata[i]; + uopts.uout << std::setw(WIDTH) << mzdata[i]; + uopts.uout << std::setw(WIDTH) << etdata[i]; + } + uopts.uout << std::endl; + } + } + + return 0; +} + +// Finalize output +static int CloseOutput(ARKODEParameters& uopts) +{ + // Footer for status output + if (uopts.output) + { + std::cout + << " -----------------------------------------------------------------" + "---------" + << std::endl; + std::cout << std::endl; + } + + // Close output streams + if (uopts.output >= 2) { uopts.uout.close(); } + + return 0; +} + +//---- end of file ---- diff --git a/examples/arkode/CXX_manyvector/ark_sod_lsrk.out b/examples/arkode/CXX_manyvector/ark_sod_lsrk.out new file mode 100644 index 0000000000..24abdf255d --- /dev/null +++ b/examples/arkode/CXX_manyvector/ark_sod_lsrk.out @@ -0,0 +1,48 @@ + +Problem parameters and options: + --------------------------------- + gamma = 1.4 + --------------------------------- + tf = 0.1 + xl = 0 + xr = 1 + nx = 512 + dx = 0.00195312 + --------------------------------- + integrator = ARKODE_LSRK_SSP_10_4 + rtol = 0.0001 + atol = 1e-11 + fixed h = 0 + --------------------------------- + output = 1 + --------------------------------- + + t ||rho|| ||mx|| ||my|| ||mz|| ||et|| + -------------------------------------------------------------------------- + 0.00e+00 7.12610e-01 0.00000e+00 0.00000e+00 0.00000e+00 1.77658e+00 + 1.00e-02 7.09460e-01 5.28096e-02 0.00000e+00 0.00000e+00 1.76732e+00 + 2.00e-02 7.06459e-01 7.61804e-02 0.00000e+00 0.00000e+00 1.75863e+00 + 3.00e-02 7.03444e-01 9.39395e-02 0.00000e+00 0.00000e+00 1.74990e+00 + 4.00e-02 7.00418e-01 1.08850e-01 0.00000e+00 0.00000e+00 1.74113e+00 + 5.00e-02 6.97379e-01 1.21955e-01 0.00000e+00 0.00000e+00 1.73232e+00 + 6.00e-02 6.94326e-01 1.33784e-01 0.00000e+00 0.00000e+00 1.72346e+00 + 7.00e-02 6.91261e-01 1.44650e-01 0.00000e+00 0.00000e+00 1.71456e+00 + 8.00e-02 6.88182e-01 1.54756e-01 0.00000e+00 0.00000e+00 1.70561e+00 + 9.00e-02 6.85089e-01 1.64242e-01 0.00000e+00 0.00000e+00 1.69661e+00 + 1.00e-01 6.81982e-01 1.73210e-01 0.00000e+00 0.00000e+00 1.68756e+00 + -------------------------------------------------------------------------- + +Final integrator statistics: +Current time = 0.1 +Steps = 382 +Step attempts = 397 +Stability limited steps = 0 +Accuracy limited steps = 397 +Error test fails = 15 +NLS step fails = 0 +Inequality constraint fails = 0 +Initial step size = 1.122284259749542e-14 +Last step size = 0.0002334416078592081 +Current step size = 0.0002334416078592081 +RHS fn evals = 3957 +Number of stages used = 10 diff --git a/examples/arkode/CXX_manyvector/plot_sod.py b/examples/arkode/CXX_manyvector/plot_sod.py new file mode 100755 index 0000000000..2e87566718 --- /dev/null +++ b/examples/arkode/CXX_manyvector/plot_sod.py @@ -0,0 +1,231 @@ +#!/usr/bin/env python3 +# ------------------------------------------------------------------------------ +# Programmer(s): Daniel R. Reynolds @ SMU +# ------------------------------------------------------------------------------ +# SUNDIALS Copyright Start +# Copyright (c) 2002-2024, Lawrence Livermore National Security +# and Southern Methodist University. +# All rights reserved. +# +# See the top-level LICENSE and NOTICE files for details. +# +# SPDX-License-Identifier: BSD-3-Clause +# SUNDIALS Copyright End +# ------------------------------------------------------------------------------ +# matplotlib-based plotting script for the serial ark_sod_lsrk example +# ------------------------------------------------------------------------------ + +# imports +import sys, os +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.gridspec import GridSpec + +# data file name +datafile = "sod.out" + +# return with an error if the file does not exist +if not os.path.isfile(datafile): + msg = "Error: file " + datafile + " does not exist" + sys.exit(msg) + +# read solution file, storing each line as a string in a list +with open(datafile, "r") as file: + lines = file.readlines() + + # extract header information + title = lines.pop(0) + nvar = int((lines.pop(0).split())[2]) + varnames = lines.pop(0) + nt = int((lines.pop(0).split())[2]) + nx = int((lines.pop(0).split())[2]) + xl = float((lines.pop(0).split())[2]) + xr = float((lines.pop(0).split())[2]) + + # allocate solution data as 2D Python arrays + t = np.zeros((nt), dtype=float) + rho = np.zeros((nt, nx), dtype=float) + mx = np.zeros((nt, nx), dtype=float) + my = np.zeros((nt, nx), dtype=float) + mz = np.zeros((nt, nx), dtype=float) + et = np.zeros((nt, nx), dtype=float) + x = np.linspace(xl, xr, nx) + + # store remaining data into numpy arrays + for it in range(nt): + line = (lines.pop(0)).split() + t[it] = line.pop(0) + for ix in range(nx): + rho[it, ix] = line.pop(0) + mx[it, ix] = line.pop(0) + my[it, ix] = line.pop(0) + mz[it, ix] = line.pop(0) + et[it, ix] = line.pop(0) + + gamma = 1.4 + u = mx / rho + p = (gamma - 1.0) * (et - (mx * mx + my * my + mz * mz) / (2.0 * rho)) + + +# utility routines for computing the analytical solution +def fsecant(p4, p1, p5, rho1, rho5, gamma): + """ + Utility routine for exact_Riemann function below + """ + z = p4 / p5 - 1.0 + c1 = np.sqrt(gamma * p1 / rho1) + c5 = np.sqrt(gamma * p5 / rho5) + gm1 = gamma - 1.0 + gp1 = gamma + 1.0 + g2 = 2.0 * gamma + fact = gm1 / g2 * (c5 / c1) * z / np.sqrt(1.0 + gp1 / g2 * z) + fact = (1.0 - fact) ** (g2 / gm1) + return p1 * fact - p4 + + +def exact_Riemann(t, x, xI): + """ + Exact 1D Riemann problem solver (retrieves domain from EulerData structure), + based on Fortran code at http://cococubed.asu.edu/codes/riemann/exact_riemann.f + + Inputs: (t,x) location for desired solution, + xI location of discontinuity at t=0, + gamma parameter for gas equation of state + Outputs: rho, u, p (density, velocity, and pressure at (t,x)) + """ + + # begin solution + rho1 = 1.0 + p1 = 1.0 + u1 = 0.0 + rho5 = 0.125 + p5 = 0.1 + u5 = 0.0 + + # solve for post-shock pressure by secant method initial guesses + p40 = p1 + p41 = p5 + f0 = fsecant(p40, p1, p5, rho1, rho5, gamma) + itmax = 50 + eps = 1.0e-14 + for iter in range(itmax): + f1 = fsecant(p41, p1, p5, rho1, rho5, gamma) + if f1 == f0: + break + p4 = p41 - (p41 - p40) * f1 / (f1 - f0) + if (np.abs(p4 - p41) / np.abs(p41)) < eps: + break + p40 = p41 + p41 = p4 + f0 = f1 + if iter == itmax - 1: + raise ValueError("exact_Riemann iteration failed to converge") + + # compute post-shock density and velocity + z = p4 / p5 - 1.0 + c5 = np.sqrt(gamma * p5 / rho5) + + gm1 = gamma - 1.0 + gp1 = gamma + 1.0 + + fact = np.sqrt(1.0 + 0.5 * gp1 * z / gamma) + + u4 = c5 * z / (gamma * fact) + rho4 = rho5 * (1.0 + 0.5 * gp1 * z / gamma) / (1.0 + 0.5 * gm1 * z / gamma) + + # shock speed + w = c5 * fact + + # compute values at foot of rarefaction + p3 = p4 + u3 = u4 + rho3 = rho1 * (p3 / p1) ** (1.0 / gamma) + + # compute positions of waves + c1 = np.sqrt(gamma * p1 / rho1) + c3 = np.sqrt(gamma * p3 / rho3) + + xsh = xI + w * t + xcd = xI + u3 * t + xft = xI + (u3 - c3) * t + xhd = xI - c1 * t + + # compute solution as a function of position + if x < xhd: + rho = rho1 + p = p1 + u = u1 + elif x < xft: + u = 2.0 / gp1 * (c1 + (x - xI) / t) + fact = 1.0 - 0.5 * gm1 * u / c1 + rho = rho1 * fact ** (2.0 / gm1) + p = p1 * fact ** (2.0 * gamma / gm1) + elif x < xcd: + rho = rho3 + p = p3 + u = u3 + elif x < xsh: + rho = rho4 + p = p4 + u = u4 + else: + rho = rho5 + p = p5 + u = u5 + + # return with success + return rho, u, p + + +# generate analytical solutions over same mesh and times as loaded from data file +rhotrue = np.zeros((nt, nx), dtype=float) +utrue = np.zeros((nt, nx), dtype=float) +ptrue = np.zeros((nt, nx), dtype=float) +for it in range(nt): + for ix in range(nx): + rhotrue[it, ix], utrue[it, ix], ptrue[it, ix] = exact_Riemann(t[it], x[ix], 0.5) + +# plot defaults: increase default font size, increase plot width, enable LaTeX rendering +plt.rc("font", size=15) +plt.rcParams["figure.figsize"] = [7.2, 4.8] +plt.rcParams["text.usetex"] = True +plt.rcParams["figure.constrained_layout.use"] = True + +# subplots with time snapshots of the density, x-velocity, and pressure +fig = plt.figure(figsize=(10, 5)) +gs = GridSpec(3, 3, figure=fig) +ax00 = fig.add_subplot(gs[0, 0]) # left column +ax10 = fig.add_subplot(gs[1, 0]) +ax20 = fig.add_subplot(gs[2, 0]) +ax01 = fig.add_subplot(gs[0, 1]) # middle column +ax11 = fig.add_subplot(gs[1, 1]) +ax21 = fig.add_subplot(gs[2, 1]) +ax02 = fig.add_subplot(gs[0, 2]) # right column +ax12 = fig.add_subplot(gs[1, 2]) +ax22 = fig.add_subplot(gs[2, 2]) +it = 0 +ax00.plot(x, rho[it, :], "-b", x, rhotrue[it, :], ":k") +ax10.plot(x, u[it, :], "-b", x, utrue[it, :], ":k") +ax20.plot(x, p[it, :], "-b", x, ptrue[it, :], ":k") +ax00.set_title(r"$t =$ " + repr(t[it]).zfill(3)) +ax00.set_ylabel(r"$\rho$") +ax10.set_ylabel(r"$v_x$") +ax20.set_ylabel(r"$p$") +ax20.set_xlabel(r"$x$") +it = nt // 2 +ax01.plot(x, rho[it, :], "-b", x, rhotrue[it, :], ":k") +ax11.plot(x, u[it, :], "-b", x, utrue[it, :], ":k") +ax21.plot(x, p[it, :], "-b", x, ptrue[it, :], ":k") +ax01.set_title(r"$t =$ " + repr(t[it]).zfill(3)) +ax21.set_xlabel(r"$x$") +it = nt - 1 +ax02.plot(x, rho[it, :], "-b", x, rhotrue[it, :], ":k") +ax12.plot(x, u[it, :], "-b", x, utrue[it, :], ":k") +ax22.plot(x, p[it, :], "-b", x, ptrue[it, :], ":k") +ax02.set_title(r"$t =$ " + repr(t[it]).zfill(3)) +ax22.set_xlabel(r"$x$") +plt.savefig("sod_frames.png") + +plt.show() + +##### end of script ##### diff --git a/examples/arkode/CXX_parallel/CMakeLists.txt b/examples/arkode/CXX_parallel/CMakeLists.txt index 05b2a87861..45b01819d5 100644 --- a/examples/arkode/CXX_parallel/CMakeLists.txt +++ b/examples/arkode/CXX_parallel/CMakeLists.txt @@ -121,9 +121,12 @@ set(OTHER_LIBS ${EXE_EXTRA_LINK_LIBS}) # MPI examples # ------------ -set(serial_examples "ark_heat2D_p.cpp\;\;--np 2 2\;1\;4\;develop\;default") +set(parallel_examples + "ark_heat2D_p.cpp\;\;--np 2 2\;1\;4\;develop\;default" + "ark_heat2D_lsrk_p.cpp\;\;--np 2 2\;1\;4\;exclude-single\;default") + set(SUNDIALS_LIBS sundials_arkode sundials_nvecparallel) -build_examples(serial_examples CXX) +build_examples(parallel_examples CXX) # Auxiliary files to install list(APPEND ARKODE_extras plot_heat2D_p.py) @@ -238,7 +241,7 @@ endif() if(EXAMPLES_INSTALL) - set(examples_to_install "${serial_examples}") + set(examples_to_install "${parallel_examples}") set(_sundials_targets arkode nvecparallel) if(examples_cvode) diff --git a/examples/arkode/CXX_parallel/ark_heat2D_lsrk_p.cpp b/examples/arkode/CXX_parallel/ark_heat2D_lsrk_p.cpp new file mode 100644 index 0000000000..c5e7bed661 --- /dev/null +++ b/examples/arkode/CXX_parallel/ark_heat2D_lsrk_p.cpp @@ -0,0 +1,1743 @@ +/* ----------------------------------------------------------------------------- + * Programmer(s): Daniel R. Reynolds @ SMU + * + * (adapted from ark_heat2D_p.cpp, co-authored by Daniel Reynolds and David + * Gardner (LLNL)) + * ----------------------------------------------------------------------------- + * SUNDIALS Copyright Start + * Copyright (c) 2002-2024, Lawrence Livermore National Security + * and Southern Methodist University. + * All rights reserved. + * + * See the top-level LICENSE and NOTICE files for details. + * + * SPDX-License-Identifier: BSD-3-Clause + * SUNDIALS Copyright End + * ----------------------------------------------------------------------------- + * Example problem: + * + * The following test simulates a simple anisotropic 2D heat equation, + * + * u_t = kx u_xx + ky u_yy + b, + * + * for t in [0, 1] and (x,y) in [0, 1]^2, with initial conditions + * + * u(0,x,y) = sin^2(pi x) sin^2(pi y), + * + * stationary boundary conditions + * + * u_t(t,0,y) = u_t(t,1,y) = u_t(t,x,0) = u_t(t,x,1) = 0, + * + * and the heat source + * + * b(t,x,y) = -2 pi sin^2(pi x) sin^2(pi y) sin(pi t) cos(pi t) + * - kx 2 pi^2 (cos^2(pi x) - sin^2(pi x)) sin^2(pi y) cos^2(pi t) + * - ky 2 pi^2 (cos^2(pi y) - sin^2(pi y)) sin^2(pi x) cos^2(pi t). + * + * Under this setup, the problem has the analytical solution + * + * u(t,x,y) = sin^2(pi x) sin^2(pi y) cos^2(pi t). + * + * The spatial derivatives are computed using second-order centered differences, + * with the data distributed over nx * ny points on a uniform spatial grid. The + * problem is advanced in time with the LSRKStep module in ARKODE. + * Several command line options are available to change the problem parameters + * and ARKODE settings. Use the flag --help for more information. + * ---------------------------------------------------------------------------*/ + +#include +#include +#include +#include +#include +#include +#include + +#include "arkode/arkode_lsrkstep.h" // access to LSRKStep +#include "mpi.h" // MPI header file +#include "nvector/nvector_parallel.h" // access to the MPI N_Vector +#include "sunadaptcontroller/sunadaptcontroller_imexgus.h" +#include "sunadaptcontroller/sunadaptcontroller_soderlind.h" + +// Macros for problem constants +#define PI SUN_RCONST(3.141592653589793238462643383279502884197169) +#define ZERO SUN_RCONST(0.0) +#define ONE SUN_RCONST(1.0) +#define TWO SUN_RCONST(2.0) +#define EIGHT SUN_RCONST(8.0) + +// Macro to access (x,y) location in 1D NVector array +#define IDX(x, y, n) ((n) * (y) + (x)) + +using namespace std; + +// ----------------------------------------------------------------------------- +// User data structure +// ----------------------------------------------------------------------------- + +struct UserData +{ + // Diffusion coefficients in the x and y directions + sunrealtype kx; + sunrealtype ky; + + // Enable/disable forcing + bool forcing; + + // Final time + sunrealtype tf; + + // Upper bounds in x and y directions + sunrealtype xu; + sunrealtype yu; + + // Global number of nodes in the x and y directions + sunindextype nx; + sunindextype ny; + + // Global total number of nodes + sunindextype nodes; + + // Mesh spacing in the x and y directions + sunrealtype dx; + sunrealtype dy; + + // Local number of nodes in the x and y directions + sunindextype nx_loc; + sunindextype ny_loc; + + // Overall number of local nodes + sunindextype nodes_loc; + + // Global x and y indices of this subdomain + sunindextype is; // x starting index + sunindextype ie; // x ending index + sunindextype js; // y starting index + sunindextype je; // y ending index + + // MPI variables + MPI_Comm comm_c; // Cartesian communicator in space + + int nprocs_w; // total number of MPI processes in Comm world + int npx; // number of MPI processes in the x-direction + int npy; // number of MPI processes in the y-direction + + int myid_c; // process ID in Cartesian communicator + + // Flags denoting if this process has a neighbor + bool HaveNbrW; + bool HaveNbrE; + bool HaveNbrS; + bool HaveNbrN; + + // Neighbor IDs for exchange + int ipW; + int ipE; + int ipS; + int ipN; + + // Receive buffers for neighbor exchange + sunrealtype* Wrecv; + sunrealtype* Erecv; + sunrealtype* Srecv; + sunrealtype* Nrecv; + + // Receive requests for neighbor exchange + MPI_Request reqRW; + MPI_Request reqRE; + MPI_Request reqRS; + MPI_Request reqRN; + + // Send buffers for neighbor exchange + sunrealtype* Wsend; + sunrealtype* Esend; + sunrealtype* Ssend; + sunrealtype* Nsend; + + // Send requests for neighbor exchange + MPI_Request reqSW; + MPI_Request reqSE; + MPI_Request reqSS; + MPI_Request reqSN; + + // Integrator settings + sunrealtype rtol; // relative tolerance + sunrealtype atol; // absolute tolerance + sunrealtype hfixed; // fixed step size + int controller; // step size adaptivity method + int maxsteps; // max number of steps between outputs + bool diagnostics; // output diagnostics + + // LSRKStep options + ARKODE_LSRKMethodType method; // LSRK method choice + int eigfrequency; // dominant eigenvalue update frequency + int stage_max_limit; // maximum number of stages per step + sunrealtype eigsafety; // dominant eigenvalue safety factor + + // Output variables + int output; // output level + int nout; // number of output times + ofstream uout; // output file stream + ofstream eout; // error file stream + N_Vector e; // error vector + + // Timing variables + bool timing; // print timings + double evolvetime; + double rhstime; + double exchangetime; +}; + +// ----------------------------------------------------------------------------- +// Functions provided to the SUNDIALS integrator +// ----------------------------------------------------------------------------- + +// ODE right hand side function +static int f(sunrealtype t, N_Vector u, N_Vector f, void* user_data); + +// Spectral radius estimation routine +static int eig(sunrealtype t, N_Vector y, N_Vector fn, sunrealtype* lambdaR, + sunrealtype* lambdaI, void* user_data, N_Vector temp1, + N_Vector temp2, N_Vector temp3); + +// ----------------------------------------------------------------------------- +// Helper functions +// ----------------------------------------------------------------------------- + +// Setup the parallel decomposition +static int SetupDecomp(MPI_Comm comm_w, UserData* udata); + +// Perform neighbor exchange +static int PostRecv(UserData* udata); +static int SendData(N_Vector y, UserData* udata); +static int WaitRecv(UserData* udata); + +// ----------------------------------------------------------------------------- +// UserData and input functions +// ----------------------------------------------------------------------------- + +// Set the default values in the UserData structure +static int InitUserData(UserData* udata); + +// Free memory allocated within UserData +static int FreeUserData(UserData* udata); + +// Read the command line inputs and set UserData values +static int ReadInputs(int* argc, char*** argv, UserData* udata, bool outproc); + +// ----------------------------------------------------------------------------- +// Output and utility functions +// ----------------------------------------------------------------------------- + +// Compute the true solution +static int Solution(sunrealtype t, N_Vector u, UserData* udata); + +// Compute the solution error solution +static int SolutionError(sunrealtype t, N_Vector u, N_Vector e, UserData* udata); + +// Print the command line options +static void InputHelp(); + +// Print some UserData information +static int PrintUserData(UserData* udata); + +// Output solution and error +static int OpenOutput(UserData* udata); +static int WriteOutput(sunrealtype t, N_Vector u, UserData* udata); +static int CloseOutput(UserData* udata); + +// Print integration timing +static int OutputTiming(UserData* udata); + +// Check function return values +static int check_flag(void* flagvalue, const string funcname, int opt); + +// ----------------------------------------------------------------------------- +// Main Program +// ----------------------------------------------------------------------------- + +int main(int argc, char* argv[]) +{ + int flag; // reusable error-checking flag + UserData* udata = NULL; // user data structure + N_Vector u = NULL; // vector for storing solution + void* arkode_mem = NULL; // ARKODE memory structure + SUNAdaptController C = NULL; // timestep adaptivity controller + + // Timing variables + double t1 = 0.0; + double t2 = 0.0; + + // MPI variables + MPI_Comm comm_w = MPI_COMM_WORLD; // MPI communicator + int myid; // MPI process ID + + // Initialize MPI + flag = MPI_Init(&argc, &argv); + if (check_flag(&flag, "MPI_Init", 1)) { return 1; } + + flag = MPI_Comm_rank(comm_w, &myid); + if (check_flag(&flag, "MPI_Comm_rank", 1)) { return 1; } + + // Create the SUNDIALS context object for this simulation + SUNContext ctx; + flag = SUNContext_Create(comm_w, &ctx); + if (check_flag(&flag, "SUNContext_Create", 1)) { return 1; } + + // Set output process flag + bool outproc = (myid == 0); + + // ------------------------------------------ + // Setup UserData and parallel decomposition + // ------------------------------------------ + + // Allocate and initialize user data structure with default values. The + // defaults may be overwritten by command line inputs in ReadInputs below. + udata = new UserData; + flag = InitUserData(udata); + if (check_flag(&flag, "InitUserData", 1)) { return 1; } + + // Parse command line inputs + flag = ReadInputs(&argc, &argv, udata, outproc); + if (flag != 0) { return 1; } + + // Setup parallel decomposition + flag = SetupDecomp(comm_w, udata); + if (check_flag(&flag, "SetupDecomp", 1)) { return 1; } + + // Output problem setup/options + if (outproc) + { + flag = PrintUserData(udata); + if (check_flag(&flag, "PrintUserData", 1)) { return 1; } + } + + if (udata->diagnostics) + { + SUNLogger logger = NULL; + + flag = SUNContext_GetLogger(ctx, &logger); + if (check_flag(&flag, "SUNContext_GetLogger", 1)) { return 1; } + + flag = SUNLogger_SetInfoFilename(logger, "diagnostics.txt"); + if (check_flag(&flag, "SUNLogger_SetInfoFilename", 1)) { return 1; } + + flag = SUNLogger_SetDebugFilename(logger, "diagnostics.txt"); + if (check_flag(&flag, "SUNLogger_SetDebugFilename", 1)) { return 1; } + } + + // ------------------------ + // Create parallel vectors + // ------------------------ + + // Create vector for solution + u = N_VNew_Parallel(udata->comm_c, udata->nodes_loc, udata->nodes, ctx); + if (check_flag((void*)u, "N_VNew_Parallel", 0)) { return 1; } + + // Set initial condition + flag = Solution(ZERO, u, udata); + if (check_flag(&flag, "Solution", 1)) { return 1; } + + // Create vector for error + udata->e = N_VClone(u); + if (check_flag((void*)(udata->e), "N_VClone", 0)) { return 1; } + + // -------------- + // Setup ARKODE + // -------------- + + // Create integrator + arkode_mem = LSRKStepCreateSTS(f, ZERO, u, ctx); + if (check_flag((void*)arkode_mem, "LSRKStepCreateSTS", 0)) { return 1; } + + // Specify tolerances + flag = ARKodeSStolerances(arkode_mem, udata->rtol, udata->atol); + if (check_flag(&flag, "ARKodeSStolerances", 1)) { return 1; } + + // Attach user data + flag = ARKodeSetUserData(arkode_mem, (void*)udata); + if (check_flag(&flag, "ARKodeSetUserData", 1)) { return 1; } + + // Select LSRK method + flag = LSRKStepSetMethod(arkode_mem, udata->method); + if (check_flag(&flag, "LSRKStepSetMethod", 1)) { return 1; } + + // Select LSRK spectral radius function and options + flag = LSRKStepSetDomEigFn(arkode_mem, eig); + if (check_flag(&flag, "LSRKStepSetDomEigFn", 1)) { return 1; } + flag = LSRKStepSetDomEigFrequency(arkode_mem, udata->eigfrequency); + if (check_flag(&flag, "LSRKStepSetDomEigFrequency", 1)) { return 1; } + + // Set maximum number of stages per step + flag = LSRKStepSetMaxNumStages(arkode_mem, udata->stage_max_limit); + if (check_flag(&flag, "LSRKStepSetMaxNumStages", 1)) { return 1; } + + // Set spectral radius safety factor + flag = LSRKStepSetDomEigSafetyFactor(arkode_mem, udata->eigsafety); + if (check_flag(&flag, "LSRKStepSetDomEigSafetyFactor", 1)) { return 1; } + + // Set fixed step size or adaptivity method + if (udata->hfixed > ZERO) + { + flag = ARKodeSetFixedStep(arkode_mem, udata->hfixed); + if (check_flag(&flag, "ARKodeSetFixedStep", 1)) { return 1; } + } + else + { + switch (udata->controller) + { + case (ARK_ADAPT_PID): C = SUNAdaptController_PID(ctx); break; + case (ARK_ADAPT_PI): C = SUNAdaptController_PI(ctx); break; + case (ARK_ADAPT_I): C = SUNAdaptController_I(ctx); break; + case (ARK_ADAPT_EXP_GUS): C = SUNAdaptController_ExpGus(ctx); break; + case (ARK_ADAPT_IMP_GUS): C = SUNAdaptController_ImpGus(ctx); break; + case (ARK_ADAPT_IMEX_GUS): C = SUNAdaptController_ImExGus(ctx); break; + } + flag = ARKodeSetAdaptController(arkode_mem, C); + if (check_flag(&flag, "ARKodeSetAdaptController", 1)) { return 1; } + } + + // Set max steps between outputs + flag = ARKodeSetMaxNumSteps(arkode_mem, udata->maxsteps); + if (check_flag(&flag, "ARKodeSetMaxNumSteps", 1)) { return 1; } + + // Set stopping time + flag = ARKodeSetStopTime(arkode_mem, udata->tf); + if (check_flag(&flag, "ARKodeSetStopTime", 1)) { return 1; } + + // ----------------------- + // Loop over output times + // ----------------------- + + sunrealtype t = ZERO; + sunrealtype dTout = udata->tf / udata->nout; + sunrealtype tout = dTout; + + // Initial output + flag = OpenOutput(udata); + if (check_flag(&flag, "OpenOutput", 1)) { return 1; } + + flag = WriteOutput(t, u, udata); + if (check_flag(&flag, "WriteOutput", 1)) { return 1; } + + for (int iout = 0; iout < udata->nout; iout++) + { + // Start timer + t1 = MPI_Wtime(); + + // Evolve in time + flag = ARKodeEvolve(arkode_mem, tout, u, &t, ARK_NORMAL); + if (check_flag(&flag, "ARKodeEvolve", 1)) { break; } + + // Stop timer + t2 = MPI_Wtime(); + + // Update timer + udata->evolvetime += t2 - t1; + + // Output solution and error + flag = WriteOutput(t, u, udata); + if (check_flag(&flag, "WriteOutput", 1)) { return 1; } + + // Update output time + tout += dTout; + tout = (tout > udata->tf) ? udata->tf : tout; + } + + // Close output + flag = CloseOutput(udata); + if (check_flag(&flag, "CloseOutput", 1)) { return 1; } + + // -------------- + // Final outputs + // -------------- + + // Print final integrator stats + if (udata->output > 0 && outproc) + { + cout << "Final integrator statistics:" << endl; + flag = ARKodePrintAllStats(arkode_mem, stdout, SUN_OUTPUTFORMAT_TABLE); + if (check_flag(&flag, "ARKodePrintAllStats", 1)) { return 1; } + } + + if (udata->forcing) + { + // Output final error + flag = SolutionError(t, u, udata->e, udata); + if (check_flag(&flag, "SolutionError", 1)) { return 1; } + + sunrealtype maxerr = N_VMaxNorm(udata->e); + + if (outproc) + { + cout << scientific; + cout << setprecision(numeric_limits::digits10); + cout << " Max error = " << maxerr << endl; + } + } + + // Print timing + if (udata->timing) + { + flag = OutputTiming(udata); + if (check_flag(&flag, "OutputTiming", 1)) { return 1; } + } + + // -------------------- + // Clean up and return + // -------------------- + + ARKodeFree(&arkode_mem); // Free integrator memory + N_VDestroy(u); // Free vectors + FreeUserData(udata); // Free user data + delete udata; + (void)SUNAdaptController_Destroy(C); // Free timestep adaptivity controller + SUNContext_Free(&ctx); // Free context + flag = MPI_Finalize(); // Finalize MPI + return 0; +} + +// ----------------------------------------------------------------------------- +// Setup the parallel decomposition +// ----------------------------------------------------------------------------- + +static int SetupDecomp(MPI_Comm comm_w, UserData* udata) +{ + int flag; + + // Check that this has not been called before + if (udata->Erecv != NULL || udata->Wrecv != NULL || udata->Srecv != NULL || + udata->Nrecv != NULL) + { + cerr << "SetupDecomp error: parallel decomposition already set up" << endl; + return -1; + } + + // Get the number of processes + flag = MPI_Comm_size(comm_w, &(udata->nprocs_w)); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Comm_size = " << flag << endl; + return -1; + } + + // Check the processor grid + if ((udata->npx * udata->npy) != udata->nprocs_w) + { + cerr << "Error: npx * npy != nproc" << endl; + return -1; + } + + // Set up 2D Cartesian communicator + int dims[2]; + dims[0] = udata->npx; + dims[1] = udata->npy; + + int periods[2]; + periods[0] = 0; + periods[1] = 0; + + flag = MPI_Cart_create(comm_w, 2, dims, periods, 0, &(udata->comm_c)); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Cart_create = " << flag << endl; + return -1; + } + + // Get my rank in the new Cartesian communicator + flag = MPI_Comm_rank(udata->comm_c, &(udata->myid_c)); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Comm_rank = " << flag << endl; + return -1; + } + + // Get dimension of the Cartesian communicator and my coordinates + int coords[2]; + flag = MPI_Cart_get(udata->comm_c, 2, dims, periods, coords); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Cart_get = " << flag << endl; + return -1; + } + + // Determine local extents in x-direction + int idx = coords[0]; + sunindextype qx = udata->nx / dims[0]; + sunindextype rx = udata->nx % dims[0]; + + udata->is = qx * idx + (idx < rx ? idx : rx); + udata->ie = udata->is + qx - 1 + (idx < rx ? 1 : 0); + + // Sanity check + if (udata->ie > (udata->nx - 1)) + { + cerr << "Error ie > nx - 1" << endl; + return -1; + } + + // Determine local extents in y-direction + int idy = coords[1]; + sunindextype qy = udata->ny / dims[1]; + sunindextype ry = udata->ny % dims[1]; + + udata->js = qy * idy + (idy < ry ? idy : ry); + udata->je = udata->js + qy - 1 + (idy < ry ? 1 : 0); + + // Sanity check + if (udata->je > (udata->ny - 1)) + { + cerr << "Error je > ny - 1" << endl; + return -1; + } + + // Number of local nodes + udata->nx_loc = (udata->ie) - (udata->is) + 1; + udata->ny_loc = (udata->je) - (udata->js) + 1; + + // Initialize global and local vector lengths + udata->nodes = udata->nx * udata->ny; + udata->nodes_loc = udata->nx_loc * udata->ny_loc; + + // Determine if this proc has neighbors + udata->HaveNbrW = (udata->is != 0); + udata->HaveNbrE = (udata->ie != udata->nx - 1); + udata->HaveNbrS = (udata->js != 0); + udata->HaveNbrN = (udata->je != udata->ny - 1); + + // Allocate exchange buffers if necessary + if (udata->HaveNbrW) + { + udata->Wrecv = new sunrealtype[udata->ny_loc]; + udata->Wsend = new sunrealtype[udata->ny_loc]; + } + if (udata->HaveNbrE) + { + udata->Erecv = new sunrealtype[udata->ny_loc]; + udata->Esend = new sunrealtype[udata->ny_loc]; + } + if (udata->HaveNbrS) + { + udata->Srecv = new sunrealtype[udata->nx_loc]; + udata->Ssend = new sunrealtype[udata->nx_loc]; + } + if (udata->HaveNbrN) + { + udata->Nrecv = new sunrealtype[udata->nx_loc]; + udata->Nsend = new sunrealtype[udata->nx_loc]; + } + + // MPI neighborhood information + int nbcoords[2]; + + // West neighbor + if (udata->HaveNbrW) + { + nbcoords[0] = coords[0] - 1; + nbcoords[1] = coords[1]; + flag = MPI_Cart_rank(udata->comm_c, nbcoords, &(udata->ipW)); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Cart_rank = " << flag << endl; + return -1; + } + } + + // East neighbor + if (udata->HaveNbrE) + { + nbcoords[0] = coords[0] + 1; + nbcoords[1] = coords[1]; + flag = MPI_Cart_rank(udata->comm_c, nbcoords, &(udata->ipE)); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Cart_rank = " << flag << endl; + return -1; + } + } + + // South neighbor + if (udata->HaveNbrS) + { + nbcoords[0] = coords[0]; + nbcoords[1] = coords[1] - 1; + flag = MPI_Cart_rank(udata->comm_c, nbcoords, &(udata->ipS)); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Cart_rank = " << flag << endl; + return -1; + } + } + + // North neighbor + if (udata->HaveNbrN) + { + nbcoords[0] = coords[0]; + nbcoords[1] = coords[1] + 1; + flag = MPI_Cart_rank(udata->comm_c, nbcoords, &(udata->ipN)); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Cart_rank = " << flag << endl; + return -1; + } + } + + // Return success + return 0; +} + +// ----------------------------------------------------------------------------- +// Functions called by the integrator +// ----------------------------------------------------------------------------- + +// f routine to compute the ODE RHS function f(t,y). +static int f(sunrealtype t, N_Vector u, N_Vector f, void* user_data) +{ + int flag; + sunindextype i, j; + + // Start timer + double t1 = MPI_Wtime(); + + // Access problem data + UserData* udata = (UserData*)user_data; + + // Open exchange receives + flag = PostRecv(udata); + if (check_flag(&flag, "PostRecv", 1)) { return -1; } + + // Send exchange data + flag = SendData(u, udata); + if (check_flag(&flag, "SendData", 1)) { return -1; } + + // Shortcuts to local number of nodes + sunindextype nx_loc = udata->nx_loc; + sunindextype ny_loc = udata->ny_loc; + + // Determine iteration range excluding the overall domain boundary + sunindextype istart = (udata->HaveNbrW) ? 0 : 1; + sunindextype iend = (udata->HaveNbrE) ? nx_loc : nx_loc - 1; + sunindextype jstart = (udata->HaveNbrS) ? 0 : 1; + sunindextype jend = (udata->HaveNbrN) ? ny_loc : ny_loc - 1; + + // Constants for computing diffusion term + sunrealtype cx = udata->kx / (udata->dx * udata->dx); + sunrealtype cy = udata->ky / (udata->dy * udata->dy); + sunrealtype cc = -TWO * (cx + cy); + + // Access data arrays + sunrealtype* uarray = N_VGetArrayPointer(u); + if (check_flag((void*)uarray, "N_VGetArrayPointer", 0)) { return -1; } + + sunrealtype* farray = N_VGetArrayPointer(f); + if (check_flag((void*)farray, "N_VGetArrayPointer", 0)) { return -1; } + + // Initialize rhs vector to zero (handles boundary conditions) + N_VConst(ZERO, f); + + // Iterate over subdomain and compute rhs forcing term + if (udata->forcing) + { + sunrealtype x, y; + sunrealtype sin_sqr_x, sin_sqr_y; + sunrealtype cos_sqr_x, cos_sqr_y; + + sunrealtype bx = (udata->kx) * TWO * PI * PI; + sunrealtype by = (udata->ky) * TWO * PI * PI; + + sunrealtype sin_t_cos_t = sin(PI * t) * cos(PI * t); + sunrealtype cos_sqr_t = cos(PI * t) * cos(PI * t); + + for (j = jstart; j < jend; j++) + { + for (i = istart; i < iend; i++) + { + x = (udata->is + i) * udata->dx; + y = (udata->js + j) * udata->dy; + + sin_sqr_x = sin(PI * x) * sin(PI * x); + sin_sqr_y = sin(PI * y) * sin(PI * y); + + cos_sqr_x = cos(PI * x) * cos(PI * x); + cos_sqr_y = cos(PI * y) * cos(PI * y); + + farray[IDX(i, j, nx_loc)] = + -TWO * PI * sin_sqr_x * sin_sqr_y * sin_t_cos_t - + bx * (cos_sqr_x - sin_sqr_x) * sin_sqr_y * cos_sqr_t - + by * (cos_sqr_y - sin_sqr_y) * sin_sqr_x * cos_sqr_t; + } + } + } + + // Iterate over subdomain interior and add rhs diffusion term + for (j = 1; j < ny_loc - 1; j++) + { + for (i = 1; i < nx_loc - 1; i++) + { + farray[IDX(i, j, nx_loc)] += + cc * uarray[IDX(i, j, nx_loc)] + + cx * (uarray[IDX(i - 1, j, nx_loc)] + uarray[IDX(i + 1, j, nx_loc)]) + + cy * (uarray[IDX(i, j - 1, nx_loc)] + uarray[IDX(i, j + 1, nx_loc)]); + } + } + + // Wait for exchange receives + flag = WaitRecv(udata); + if (check_flag(&flag, "WaitRecv", 1)) { return -1; } + + // Iterate over subdomain boundaries and add rhs diffusion term + sunrealtype* Warray = udata->Wrecv; + sunrealtype* Earray = udata->Erecv; + sunrealtype* Sarray = udata->Srecv; + sunrealtype* Narray = udata->Nrecv; + + // West face (updates south-west and north-west corners if necessary) + if (udata->HaveNbrW) + { + i = 0; + if (udata->HaveNbrS) // South-West corner + { + j = 0; + farray[IDX(i, j, nx_loc)] += + cc * uarray[IDX(i, j, nx_loc)] + + cx * (Warray[j] + uarray[IDX(i + 1, j, nx_loc)]) + + cy * (Sarray[i] + uarray[IDX(i, j + 1, nx_loc)]); + } + + for (j = 1; j < ny_loc - 1; j++) + { + farray[IDX(i, j, nx_loc)] += + cc * uarray[IDX(i, j, nx_loc)] + + cx * (Warray[j] + uarray[IDX(i + 1, j, nx_loc)]) + + cy * (uarray[IDX(i, j - 1, nx_loc)] + uarray[IDX(i, j + 1, nx_loc)]); + } + + if (udata->HaveNbrN) // North-West corner + { + j = ny_loc - 1; + farray[IDX(i, j, nx_loc)] += + cc * uarray[IDX(i, j, nx_loc)] + + cx * (Warray[j] + uarray[IDX(i + 1, j, nx_loc)]) + + cy * (uarray[IDX(i, j - 1, nx_loc)] + Narray[i]); + } + } + + // East face (updates south-east and north-east corners if necessary) + if (udata->HaveNbrE) + { + i = nx_loc - 1; + if (udata->HaveNbrS) // South-East corner + { + j = 0; + farray[IDX(i, j, nx_loc)] += + cc * uarray[IDX(i, j, nx_loc)] + + cx * (uarray[IDX(i - 1, j, nx_loc)] + Earray[j]) + + cy * (Sarray[i] + uarray[IDX(i, j + 1, nx_loc)]); + } + + for (j = 1; j < ny_loc - 1; j++) + { + farray[IDX(i, j, nx_loc)] += + cc * uarray[IDX(i, j, nx_loc)] + + cx * (uarray[IDX(i - 1, j, nx_loc)] + Earray[j]) + + cy * (uarray[IDX(i, j - 1, nx_loc)] + uarray[IDX(i, j + 1, nx_loc)]); + } + + if (udata->HaveNbrN) // North-East corner + { + j = ny_loc - 1; + farray[IDX(i, j, nx_loc)] += + cc * uarray[IDX(i, j, nx_loc)] + + cx * (uarray[IDX(i - 1, j, nx_loc)] + Earray[j]) + + cy * (uarray[IDX(i, j - 1, nx_loc)] + Narray[i]); + } + } + + // South face (excludes corners) + if (udata->HaveNbrS) + { + j = 0; + for (i = 1; i < nx_loc - 1; i++) + { + farray[IDX(i, j, nx_loc)] += + cc * uarray[IDX(i, j, nx_loc)] + + cx * (uarray[IDX(i - 1, j, nx_loc)] + uarray[IDX(i + 1, j, nx_loc)]) + + cy * (Sarray[i] + uarray[IDX(i, j + 1, nx_loc)]); + } + } + + // North face (excludes corners) + if (udata->HaveNbrN) + { + j = udata->ny_loc - 1; + for (i = 1; i < nx_loc - 1; i++) + { + farray[IDX(i, j, nx_loc)] += + cc * uarray[IDX(i, j, nx_loc)] + + cx * (uarray[IDX(i - 1, j, nx_loc)] + uarray[IDX(i + 1, j, nx_loc)]) + + cy * (uarray[IDX(i, j - 1, nx_loc)] + Narray[i]); + } + } + + // Stop timer + double t2 = MPI_Wtime(); + + // Update timer + udata->rhstime += t2 - t1; + + // Return success + return 0; +} + +// ----------------------------------------------------------------------------- +// RHS helper functions +// ----------------------------------------------------------------------------- + +// Post exchange receives +static int PostRecv(UserData* udata) +{ + int flag; + + // Start timer + double t1 = MPI_Wtime(); + + // Open Irecv buffers + if (udata->HaveNbrW) + { + flag = MPI_Irecv(udata->Wrecv, (int)udata->ny_loc, MPI_SUNREALTYPE, + udata->ipW, MPI_ANY_TAG, udata->comm_c, &(udata->reqRW)); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Irecv = " << flag << endl; + return -1; + } + } + + if (udata->HaveNbrE) + { + flag = MPI_Irecv(udata->Erecv, (int)udata->ny_loc, MPI_SUNREALTYPE, + udata->ipE, MPI_ANY_TAG, udata->comm_c, &(udata->reqRE)); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Irecv = " << flag << endl; + return -1; + } + } + + if (udata->HaveNbrS) + { + flag = MPI_Irecv(udata->Srecv, (int)udata->nx_loc, MPI_SUNREALTYPE, + udata->ipS, MPI_ANY_TAG, udata->comm_c, &(udata->reqRS)); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Irecv = " << flag << endl; + return -1; + } + } + + if (udata->HaveNbrN) + { + flag = MPI_Irecv(udata->Nrecv, (int)udata->nx_loc, MPI_SUNREALTYPE, + udata->ipN, MPI_ANY_TAG, udata->comm_c, &(udata->reqRN)); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Irecv = " << flag << endl; + return -1; + } + } + + // Stop timer + double t2 = MPI_Wtime(); + + // Update timer + udata->exchangetime += t2 - t1; + + // Return success + return 0; +} + +// Send exchange data +static int SendData(N_Vector y, UserData* udata) +{ + int flag, i; + sunindextype ny_loc = udata->ny_loc; + sunindextype nx_loc = udata->nx_loc; + + // Start timer + double t1 = MPI_Wtime(); + + // Access data array + sunrealtype* Y = N_VGetArrayPointer(y); + if (check_flag((void*)Y, "N_VGetArrayPointer", 0)) { return -1; } + + // Send data + if (udata->HaveNbrW) + { + for (i = 0; i < ny_loc; i++) { udata->Wsend[i] = Y[IDX(0, i, nx_loc)]; } + flag = MPI_Isend(udata->Wsend, (int)udata->ny_loc, MPI_SUNREALTYPE, + udata->ipW, 0, udata->comm_c, &(udata->reqSW)); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Isend = " << flag << endl; + return -1; + } + } + + if (udata->HaveNbrE) + { + for (i = 0; i < ny_loc; i++) + { + udata->Esend[i] = Y[IDX(nx_loc - 1, i, nx_loc)]; + } + flag = MPI_Isend(udata->Esend, (int)udata->ny_loc, MPI_SUNREALTYPE, + udata->ipE, 1, udata->comm_c, &(udata->reqSE)); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Isend = " << flag << endl; + return -1; + } + } + + if (udata->HaveNbrS) + { + for (i = 0; i < nx_loc; i++) { udata->Ssend[i] = Y[IDX(i, 0, nx_loc)]; } + flag = MPI_Isend(udata->Ssend, (int)udata->nx_loc, MPI_SUNREALTYPE, + udata->ipS, 2, udata->comm_c, &(udata->reqSS)); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Isend = " << flag << endl; + return -1; + } + } + + if (udata->HaveNbrN) + { + for (i = 0; i < nx_loc; i++) + { + udata->Nsend[i] = Y[IDX(i, ny_loc - 1, nx_loc)]; + } + flag = MPI_Isend(udata->Nsend, (int)udata->nx_loc, MPI_SUNREALTYPE, + udata->ipN, 3, udata->comm_c, &(udata->reqSN)); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Isend = " << flag << endl; + return -1; + } + } + + // Stop timer + double t2 = MPI_Wtime(); + + // Update timer + udata->exchangetime += t2 - t1; + + // Return success + return 0; +} + +// Wait for exchange data +static int WaitRecv(UserData* udata) +{ + // Local variables + int flag; + MPI_Status stat; + + // Start timer + double t1 = MPI_Wtime(); + + // Wait for messages to finish + if (udata->HaveNbrW) + { + flag = MPI_Wait(&(udata->reqRW), &stat); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Wait = " << flag << endl; + return -1; + } + flag = MPI_Wait(&(udata->reqSW), &stat); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Wait = " << flag << endl; + return -1; + } + } + + if (udata->HaveNbrE) + { + flag = MPI_Wait(&(udata->reqRE), &stat); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Wait = " << flag << endl; + return -1; + } + flag = MPI_Wait(&(udata->reqSE), &stat); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Wait = " << flag << endl; + return -1; + } + } + + if (udata->HaveNbrS) + { + flag = MPI_Wait(&(udata->reqRS), &stat); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Wait = " << flag << endl; + return -1; + } + flag = MPI_Wait(&(udata->reqSS), &stat); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Wait = " << flag << endl; + return -1; + } + } + + if (udata->HaveNbrN) + { + flag = MPI_Wait(&(udata->reqRN), &stat); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Wait = " << flag << endl; + return -1; + } + flag = MPI_Wait(&(udata->reqSN), &stat); + if (flag != MPI_SUCCESS) + { + cerr << "Error in MPI_Wait = " << flag << endl; + return -1; + } + } + + // Stop timer + double t2 = MPI_Wtime(); + + // Update timer + udata->exchangetime += t2 - t1; + + // Return success + return 0; +} + +// Spectral radius estimation routine +static int eig(sunrealtype t, N_Vector y, N_Vector fn, sunrealtype* lambdaR, + sunrealtype* lambdaI, void* user_data, N_Vector temp1, + N_Vector temp2, N_Vector temp3) +{ + // Access problem data + UserData* udata = (UserData*)user_data; + + // Fill in spectral radius value + *lambdaR = -SUN_RCONST(8.0) * max(udata->kx / udata->dx / udata->dx, + udata->ky / udata->dy / udata->dy); + *lambdaI = SUN_RCONST(0.0); + + // return with success + return 0; +} + +// ----------------------------------------------------------------------------- +// UserData and input functions +// ----------------------------------------------------------------------------- + +// Initialize memory allocated within Userdata +static int InitUserData(UserData* udata) +{ + // Diffusion coefficient + udata->kx = SUN_RCONST(10.0); + udata->ky = SUN_RCONST(10.0); + + // Enable forcing + udata->forcing = true; + + // Final time + udata->tf = ONE; + + // Upper bounds in x and y directions + udata->xu = ONE; + udata->yu = ONE; + + // Global number of nodes in the x and y directions + udata->nx = 64; + udata->ny = 64; + udata->nodes = udata->nx * udata->ny; + + // Mesh spacing in the x and y directions + udata->dx = udata->xu / (udata->nx - 1); + udata->dy = udata->yu / (udata->ny - 1); + + // Locals number of nodes in the x and y directions (set in SetupDecomp) + udata->nx_loc = 0; + udata->ny_loc = 0; + udata->nodes_loc = 0; + + // Global indices of this subdomain (set in SetupDecomp) + udata->is = 0; + udata->ie = 0; + udata->js = 0; + udata->je = 0; + + // MPI variables (set in SetupDecomp) + udata->comm_c = MPI_COMM_NULL; + + udata->nprocs_w = 1; + udata->npx = 1; + udata->npy = 1; + + udata->myid_c = 0; + + // Flags denoting neighbors (set in SetupDecomp) + udata->HaveNbrW = true; + udata->HaveNbrE = true; + udata->HaveNbrS = true; + udata->HaveNbrN = true; + + // Exchange receive buffers (allocated in SetupDecomp) + udata->Erecv = NULL; + udata->Wrecv = NULL; + udata->Nrecv = NULL; + udata->Srecv = NULL; + + // Exchange send buffers (allocated in SetupDecomp) + udata->Esend = NULL; + udata->Wsend = NULL; + udata->Nsend = NULL; + udata->Ssend = NULL; + + // Neighbors IDs (set in SetupDecomp) + udata->ipW = -1; + udata->ipE = -1; + udata->ipS = -1; + udata->ipN = -1; + + // Integrator settings + udata->rtol = SUN_RCONST(1.e-5); // relative tolerance + udata->atol = SUN_RCONST(1.e-10); // absolute tolerance + udata->hfixed = ZERO; // using adaptive step sizes + udata->controller = 0; // PID controller + udata->maxsteps = 0; // use default + udata->diagnostics = false; // output diagnostics + + // LSRKStep options + udata->method = ARKODE_LSRK_RKC_2; // RKC + udata->eigfrequency = 25; // update eigenvalue at least every 20 steps + udata->stage_max_limit = 1000; // allow up to 1000 stages/step + udata->eigsafety = SUN_RCONST(1.01); // 1% safety factor + + // Output variables + udata->output = 1; // 0 = no output, 1 = stats output, 2 = output to disk + udata->nout = 20; // Number of output times + udata->e = NULL; + + // Timing variables + udata->timing = false; + udata->evolvetime = 0.0; + udata->rhstime = 0.0; + udata->exchangetime = 0.0; + + // Return success + return 0; +} + +// Free memory allocated within Userdata +static int FreeUserData(UserData* udata) +{ + // Free exchange buffers + if (udata->Wrecv != NULL) { delete[] udata->Wrecv; } + if (udata->Wsend != NULL) { delete[] udata->Wsend; } + if (udata->Erecv != NULL) { delete[] udata->Erecv; } + if (udata->Esend != NULL) { delete[] udata->Esend; } + if (udata->Srecv != NULL) { delete[] udata->Srecv; } + if (udata->Ssend != NULL) { delete[] udata->Ssend; } + if (udata->Nrecv != NULL) { delete[] udata->Nrecv; } + if (udata->Nsend != NULL) { delete[] udata->Nsend; } + + // Free MPI Cartesian communicator + if (udata->comm_c != MPI_COMM_NULL) { MPI_Comm_free(&(udata->comm_c)); } + + // Free error vector + if (udata->e) + { + N_VDestroy(udata->e); + udata->e = NULL; + } + + // Return success + return 0; +} + +// Read command line inputs +static int ReadInputs(int* argc, char*** argv, UserData* udata, bool outproc) +{ + // Check for input args + int arg_idx = 1; + + while (arg_idx < (*argc)) + { + string arg = (*argv)[arg_idx++]; + + // Mesh points + if (arg == "--mesh") + { + udata->nx = stoi((*argv)[arg_idx++]); + udata->ny = stoi((*argv)[arg_idx++]); + } + // MPI processes + else if (arg == "--np") + { + udata->npx = stoi((*argv)[arg_idx++]); + udata->npy = stoi((*argv)[arg_idx++]); + } + // Domain upper bounds + else if (arg == "--domain") + { + udata->xu = stoi((*argv)[arg_idx++]); + udata->yu = stoi((*argv)[arg_idx++]); + } + // Diffusion parameters + else if (arg == "--k") + { + udata->kx = stod((*argv)[arg_idx++]); + udata->ky = stod((*argv)[arg_idx++]); + } + // Disable forcing + else if (arg == "--noforcing") { udata->forcing = false; } + // Temporal domain settings + else if (arg == "--tf") { udata->tf = stod((*argv)[arg_idx++]); } + // Integrator settings + else if (arg == "--rtol") { udata->rtol = stod((*argv)[arg_idx++]); } + else if (arg == "--atol") { udata->atol = stod((*argv)[arg_idx++]); } + else if (arg == "--fixedstep") { udata->hfixed = stod((*argv)[arg_idx++]); } + else if (arg == "--controller") + { + udata->controller = stoi((*argv)[arg_idx++]); + } + else if (arg == "--diagnostics") { udata->diagnostics = true; } + else if (arg == "--method") + { + udata->method = (ARKODE_LSRKMethodType)stoi((*argv)[arg_idx++]); + } + else if (arg == "--eigfrequency") + { + udata->eigfrequency = stoi((*argv)[arg_idx++]); + } + else if (arg == "--stage_max_limit") + { + udata->stage_max_limit = stoi((*argv)[arg_idx++]); + } + else if (arg == "--eigsafety") + { + udata->eigsafety = stod((*argv)[arg_idx++]); + } + // Output settings + else if (arg == "--output") { udata->output = stoi((*argv)[arg_idx++]); } + else if (arg == "--nout") { udata->nout = stoi((*argv)[arg_idx++]); } + else if (arg == "--maxsteps") + { + udata->maxsteps = stoi((*argv)[arg_idx++]); + } + else if (arg == "--timing") { udata->timing = true; } + // Help + else if (arg == "--help") + { + if (outproc) { InputHelp(); } + return -1; + } + // Unknown input + else + { + if (outproc) + { + cerr << "ERROR: Invalid input " << arg << endl; + InputHelp(); + } + return -1; + } + } + + // Recompute total number of nodes + udata->nodes = udata->nx * udata->ny; + + // Recompute x and y mesh spacing + udata->dx = (udata->xu) / (udata->nx - 1); + udata->dy = (udata->yu) / (udata->ny - 1); + + // Return success + return 0; +} + +// ----------------------------------------------------------------------------- +// Output and utility functions +// ----------------------------------------------------------------------------- + +// Compute the exact solution +static int Solution(sunrealtype t, N_Vector u, UserData* udata) +{ + sunrealtype x, y; + sunrealtype cos_sqr_t; + sunrealtype sin_sqr_x, sin_sqr_y; + + // Constants for computing solution + cos_sqr_t = cos(PI * t) * cos(PI * t); + + // Initialize u to zero (handles boundary conditions) + N_VConst(ZERO, u); + + // Iterative over domain interior + sunindextype istart = (udata->HaveNbrW) ? 0 : 1; + sunindextype iend = (udata->HaveNbrE) ? udata->nx_loc : udata->nx_loc - 1; + + sunindextype jstart = (udata->HaveNbrS) ? 0 : 1; + sunindextype jend = (udata->HaveNbrN) ? udata->ny_loc : udata->ny_loc - 1; + + sunrealtype* uarray = N_VGetArrayPointer(u); + if (check_flag((void*)uarray, "N_VGetArrayPointer", 0)) { return -1; } + + for (sunindextype j = jstart; j < jend; j++) + { + for (sunindextype i = istart; i < iend; i++) + { + x = (udata->is + i) * udata->dx; + y = (udata->js + j) * udata->dy; + + sin_sqr_x = sin(PI * x) * sin(PI * x); + sin_sqr_y = sin(PI * y) * sin(PI * y); + + uarray[IDX(i, j, udata->nx_loc)] = sin_sqr_x * sin_sqr_y * cos_sqr_t; + } + } + + return 0; +} + +// Compute the solution error +static int SolutionError(sunrealtype t, N_Vector u, N_Vector e, UserData* udata) +{ + // Compute true solution + int flag = Solution(t, e, udata); + if (flag != 0) { return -1; } + + // Compute absolute error + N_VLinearSum(ONE, u, -ONE, e, e); + N_VAbs(e, e); + + return 0; +} + +// Print command line options +static void InputHelp() +{ + cout << endl; + cout << "Command line options:" << endl; + cout << " --mesh : mesh points in the x and y directions" + << endl; + cout << " --np : number of MPI processes in the x and y " + "directions" + << endl; + cout + << " --domain : domain upper bound in the x and y direction" + << endl; + cout << " --k : diffusion coefficients" << endl; + cout << " --noforcing : disable forcing term" << endl; + cout << " --tf