From ec80fb533f73e3a7b4d547a8fc498c53caf915aa Mon Sep 17 00:00:00 2001 From: maggul Date: Fri, 25 Oct 2024 10:27:31 -0700 Subject: [PATCH 01/28] ARKODE_MRI interface --- CMakeLists.txt | 2 + include/bout/mesh.hxx | 34 +++-- include/bout/monitor.hxx | 14 +- include/bout/physicsmodel.hxx | 35 +++++ include/bout/solver.hxx | 28 ++++ src/bout++.cxx | 62 ++++++--- src/mesh/impls/bout/boutmesh.cxx | 1 - src/mesh/mesh.cxx | 6 +- src/physics/physicsmodel.cxx | 9 +- src/solver/impls/arkode/arkode.cxx | 100 +++++++------- src/solver/solver.cxx | 214 ++++++++++++++++++++++++++++- tests/integrated/CMakeLists.txt | 1 + tests/unit/test_extras.hxx | 1 - 13 files changed, 415 insertions(+), 92 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 5e1924ba1e..0de058a139 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -289,6 +289,8 @@ set(BOUT_SOURCES ./src/solver/impls/adams_bashforth/adams_bashforth.hxx ./src/solver/impls/arkode/arkode.cxx ./src/solver/impls/arkode/arkode.hxx + ./src/solver/impls/arkode/arkode_mri.cxx + ./src/solver/impls/arkode/arkode_mri.hxx ./src/solver/impls/cvode/cvode.cxx ./src/solver/impls/cvode/cvode.hxx ./src/solver/impls/euler/euler.cxx diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index 7a79359ecd..c80716fc12 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -43,33 +43,43 @@ class Mesh; #ifndef BOUT_MESH_H #define BOUT_MESH_H -#include "bout/bout_enum_class.hxx" +#include "mpi.h" + +#include +#include +#include + #include "bout/bout_types.hxx" -#include "bout/coordinates.hxx" // Coordinates class #include "bout/field2d.hxx" #include "bout/field3d.hxx" #include "bout/field_data.hxx" -#include "bout/fieldgroup.hxx" -#include "bout/generic_factory.hxx" -#include "bout/index_derivs_interface.hxx" -#include "bout/mpi_wrapper.hxx" #include "bout/options.hxx" -#include "bout/region.hxx" + +#include "bout/fieldgroup.hxx" + +class BoundaryRegion; +class BoundaryRegionPar; + #include "bout/sys/range.hxx" // RangeIterator + +#include + +#include "bout/coordinates.hxx" // Coordinates class + #include "bout/unused.hxx" -#include "mpi.h" +#include "bout/generic_factory.hxx" +#include + +#include +#include #include #include #include #include #include -class BoundaryRegion; -class BoundaryRegionPar; -class GridDataSource; - class MeshFactory : public Factory { public: static constexpr auto type_name = "Mesh"; diff --git a/include/bout/monitor.hxx b/include/bout/monitor.hxx index 359096e74f..991cbd65c5 100644 --- a/include/bout/monitor.hxx +++ b/include/bout/monitor.hxx @@ -86,10 +86,18 @@ public: /// number of RHS calls int ncalls = 0; - /// number of RHS calls for fast timescale + /// number of RHS calls for explicit timescale int ncalls_e = 0; - /// number of RHS calls for slow timescale + /// number of RHS calls for implicit timescale int ncalls_i = 0; + /// number of RHS calls for slow explicit timescale + int ncalls_se = 0; + /// number of RHS calls for slow implicit timescale + int ncalls_si = 0; + /// number of RHS calls for fast explicit timescale + int ncalls_fe = 0; + /// number of RHS calls for fast implicit timescale + int ncalls_fi = 0; /// wall time spent calculating RHS BoutReal wtime_rhs = 0; @@ -122,7 +130,7 @@ public: /*! * Write job progress to screen */ - void writeProgress(BoutReal simtime, bool output_split); + void writeProgress(BoutReal simtime, bool output_split, bool output_splitmri); }; #endif // BOUT_MONITOR_H diff --git a/include/bout/physicsmodel.hxx b/include/bout/physicsmodel.hxx index 9fa25d8b0f..97c3f09a40 100644 --- a/include/bout/physicsmodel.hxx +++ b/include/bout/physicsmodel.hxx @@ -169,6 +169,12 @@ public: * * Returns a flag: 0 indicates success, non-zero an error flag */ + int runRHS_se(BoutReal time, bool linear = false); + int runRHS_si(BoutReal time, bool linear = false); + int runRHS_fe(BoutReal time, bool linear = false); + int runRHS_fi(BoutReal time, bool linear = false); + int runRHS_s(BoutReal time, bool linear = false); + int runRHS_f(BoutReal time, bool linear = false); int runRHS(BoutReal time, bool linear = false); /*! @@ -176,6 +182,11 @@ public: */ bool splitOperator(); + /*! + * True if this model uses split operators + */ + bool splitOperatorMRI(); + /*! * Run the convective (usually explicit) part of the model * @@ -267,6 +278,24 @@ protected: * which is set to true when the rhs() function can be * linearised. This is used in e.g. linear iterative solves. */ + virtual int rhs_se(BoutReal UNUSED(t)) { return 1; } + virtual int rhs_se(BoutReal t, bool UNUSED(linear)) { return rhs_se(t); } + + virtual int rhs_si(BoutReal UNUSED(t)) { return 1; } + virtual int rhs_si(BoutReal t, bool UNUSED(linear)) { return rhs_si(t); } + + virtual int rhs_fe(BoutReal UNUSED(t)) { return 1; } + virtual int rhs_fe(BoutReal t, bool UNUSED(linear)) { return rhs_fe(t); } + + virtual int rhs_fi(BoutReal UNUSED(t)) { return 1; } + virtual int rhs_fi(BoutReal t, bool UNUSED(linear)) { return rhs_fi(t); } + + virtual int rhs_s(BoutReal UNUSED(t)) { return 1; } + virtual int rhs_s(BoutReal t, bool UNUSED(linear)) { return rhs_s(t); } + + virtual int rhs_f(BoutReal UNUSED(t)) { return 1; } + virtual int rhs_f(BoutReal t, bool UNUSED(linear)) { return rhs_f(t); } + virtual int rhs(BoutReal UNUSED(t)) { return 1; } virtual int rhs(BoutReal t, bool UNUSED(linear)) { return rhs(t); } @@ -309,6 +338,10 @@ protected: /// Specify that this model is split into a convective and diffusive part void setSplitOperator(bool split = true) { splitop = split; } + /// Specify that this model is split into a convective and diffusive part + void setSplitOperatorMRI(bool split = true) { splitopmri = split; } + + /// Specify a preconditioner function void setPrecon(preconfunc pset) { userprecon = pset; } template @@ -391,6 +424,8 @@ private: bool restart_enabled{true}; /// Split operator model? bool splitop{false}; + /// Split operator model? + bool splitopmri{false}; /// Pointer to user-supplied preconditioner function preconfunc userprecon{nullptr}; /// Pointer to user-supplied Jacobian-vector multiply function diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index e35b48c620..29cab9c620 100644 --- a/include/bout/solver.hxx +++ b/include/bout/solver.hxx @@ -87,6 +87,7 @@ constexpr auto SOLVEREULER = "euler"; constexpr auto SOLVERRK3SSP = "rk3ssp"; constexpr auto SOLVERPOWER = "power"; constexpr auto SOLVERARKODE = "arkode"; +constexpr auto SOLVERARKODEMRI = "arkodemri"; constexpr auto SOLVERIMEXBDF2 = "imexbdf2"; constexpr auto SOLVERSNES = "snes"; constexpr auto SOLVERRKGENERIC = "rkgeneric"; @@ -310,9 +311,21 @@ public: /// Same but fur implicit timestep counter - for IMEX int resetRHSCounter_i(); + /// Same but for slow explicit timestep counter - for MRI IMEX + int resetRHSCounter_se(); + /// Same but for slow implicit timestep counter - for MRI IMEX + int resetRHSCounter_si(); + /// Same but for fast explicit timestep counter - for MRI IMEX + int resetRHSCounter_fe(); + /// Same but for fast implicit timestep counter - for MRI IMEX + int resetRHSCounter_fi(); + /// Test if this solver supports split operators (e.g. implicit/explicit) bool splitOperator(); + /// Test if this solver supports split operators (e.g. implicit/explicit) + bool splitOperatorMRI(); + bool canReset{false}; /// Add evolving variables to output (dump) file or restart file @@ -438,6 +451,12 @@ protected: BoutReal simtime{0.0}; /// Run the user's RHS function + int run_rhs_se(BoutReal t, bool linear = false); + int run_rhs_si(BoutReal t, bool linear = false); + int run_rhs_fe(BoutReal t, bool linear = false); + int run_rhs_fi(BoutReal t, bool linear = false); + int run_rhs_s(BoutReal t, bool linear = false); + int run_rhs_f(BoutReal t, bool linear = false); int run_rhs(BoutReal t, bool linear = false); /// Calculate only the convective parts int run_convective(BoutReal t, bool linear = false); @@ -542,6 +561,15 @@ private: int rhs_ncalls_e{0}; /// Number of calls to the implicit (diffusive) RHS function int rhs_ncalls_i{0}; + /// number of RHS calls for slow explicit timescale + int rhs_ncalls_se = 0; + /// number of RHS calls for slow implicit timescale + int rhs_ncalls_si = 0; + /// number of RHS calls for fast explicit timescale + int rhs_ncalls_fe = 0; + /// number of RHS calls for fast implicit timescale + int rhs_ncalls_fi = 0; + /// Default sampling rate at which to call monitors - same as output to screen int default_monitor_period{1}; /// timestep - shouldn't be changed after init is called. diff --git a/src/bout++.cxx b/src/bout++.cxx index 7f23cf5f91..19a9f4df41 100644 --- a/src/bout++.cxx +++ b/src/bout++.cxx @@ -849,8 +849,12 @@ int BoutMonitor::call(Solver* solver, BoutReal t, [[maybe_unused]] int iter, int run_data.ncalls = solver->resetRHSCounter(); run_data.ncalls_e = solver->resetRHSCounter_e(); run_data.ncalls_i = solver->resetRHSCounter_i(); + + run_data.ncalls_se = solver->resetRHSCounter_se(); + run_data.ncalls_si = solver->resetRHSCounter_si(); + run_data.ncalls_fe = solver->resetRHSCounter_fe(); + run_data.ncalls_fi = solver->resetRHSCounter_fi(); - const bool output_split = solver->splitOperator(); run_data.wtime_rhs = Timer::resetTime("rhs"); run_data.wtime_invert = Timer::resetTime("invert"); // Time spent communicating (part of RHS) @@ -872,16 +876,23 @@ int BoutMonitor::call(Solver* solver, BoutReal t, [[maybe_unused]] int iter, int first_time = false; // Print the column header for timing info - if (!output_split) { - output_progress.write(_("Sim Time | RHS evals | Wall Time | Calc Inv Comm " - " I/O SOLVER\n\n")); - } else { + if (solver->splitOperator()) { output_progress.write(_("Sim Time | RHS_e evals | RHS_I evals | Wall Time | " "Calc Inv Comm I/O SOLVER\n\n")); } + else if (solver->splitOperatorMRI()) { + output_progress.write(_("Sim Time | RHS_se evals | RHS_si evals | RHS_fe evals |" + "RHS_fi evals | Wall Time | " + "Calc Inv Comm I/O SOLVER\n\n")); + } + else + { + output_progress.write(_("Sim Time | RHS evals | Wall Time | Calc Inv Comm " + " I/O SOLVER\n\n")); + } } - run_data.writeProgress(simtime, output_split); + run_data.writeProgress(simtime, solver->splitOperator(), solver->splitOperatorMRI()); // This bit only to screen, not log file @@ -1011,6 +1022,10 @@ void RunMetrics::outputVars(Options& output_options) const { output_options["wtime"].assignRepeat(wtime, "t", true, "Output"); output_options["ncalls"].assignRepeat(ncalls, "t", true, "Output"); output_options["ncalls_e"].assignRepeat(ncalls_e, "t", true, "Output"); + output_options["ncalls_se"].assignRepeat(ncalls_se, "t", true, "Output"); + output_options["ncalls_si"].assignRepeat(ncalls_si, "t", true, "Output"); + output_options["ncalls_fe"].assignRepeat(ncalls_fe, "t", true, "Output"); + output_options["ncalls_fi"].assignRepeat(ncalls_fi, "t", true, "Output"); output_options["ncalls_i"].assignRepeat(ncalls_i, "t", true, "Output"); output_options["wtime_rhs"].assignRepeat(wtime_rhs, "t", true, "Output"); output_options["wtime_invert"].assignRepeat(wtime_invert, "t", true, "Output"); @@ -1036,17 +1051,8 @@ void RunMetrics::calculateDerivedMetrics() { wtime_per_rhs_i = wtime / ncalls_i; } -void RunMetrics::writeProgress(BoutReal simtime, bool output_split) { - if (!output_split) { - output_progress.write( - "{:.3e} {:5d} {:.2e} {:5.1f} {:5.1f} {:5.1f} {:5.1f} {:5.1f}\n", - simtime, ncalls, wtime, 100. * (wtime_rhs - wtime_comms - wtime_invert) / wtime, - 100. * wtime_invert / wtime, // Inversions - 100. * wtime_comms / wtime, // Communications - 100. * wtime_io / wtime, // I/O - 100. * (wtime - wtime_io - wtime_rhs) / wtime); // Everything else - - } else { +void RunMetrics::writeProgress(BoutReal simtime, bool output_split, bool output_splitmri) { + if (output_split) { output_progress.write("{:.3e} {:5d} {:5d} {:.2e} {:5.1f} " "{:5.1f} {:5.1f} {:5.1f} {:5.1f}\n", simtime, ncalls_e, ncalls_i, wtime, @@ -1057,4 +1063,26 @@ void RunMetrics::writeProgress(BoutReal simtime, bool output_split) { 100. * (wtime - wtime_io - wtime_rhs) / wtime); // Everything else } + else if (output_splitmri) { + output_progress.write("{:.3e} {:8d} {:8d} {:8d} {:8d} {:.2e} {:5.1f} " + "{:5.1f} {:5.1f} {:5.1f} {:5.1f}\n", + simtime, ncalls_se, ncalls_si, ncalls_fe, ncalls_fi, wtime, + 100. * (wtime_rhs - wtime_comms - wtime_invert) / wtime, + 100. * wtime_invert / wtime, // Inversions + 100. * wtime_comms / wtime, // Communications + 100. * wtime_io / wtime, // I/O + 100. * (wtime - wtime_io - wtime_rhs) + / wtime); // Everything else + } + else + { + output_progress.write( + "{:.3e} {:5d} {:.2e} {:5.1f} {:5.1f} {:5.1f} {:5.1f} {:5.1f}\n", + simtime, ncalls, wtime, 100. * (wtime_rhs - wtime_comms - wtime_invert) / wtime, + 100. * wtime_invert / wtime, // Inversions + 100. * wtime_comms / wtime, // Communications + 100. * wtime_io / wtime, // I/O + 100. * (wtime - wtime_io - wtime_rhs) / wtime); // Everything else + + } } diff --git a/src/mesh/impls/bout/boutmesh.cxx b/src/mesh/impls/bout/boutmesh.cxx index 31319b87d1..8ef651b0cd 100644 --- a/src/mesh/impls/bout/boutmesh.cxx +++ b/src/mesh/impls/bout/boutmesh.cxx @@ -42,7 +42,6 @@ #include #include #include -#include #include #include #include diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx index 6d7a5de512..870f3413cd 100644 --- a/src/mesh/mesh.cxx +++ b/src/mesh/mesh.cxx @@ -1,15 +1,15 @@ -#include #include #include #include -#include #include #include -#include #include #include +#include +#include + #include "impls/bout/boutmesh.hxx" MeshFactory::ReturnType MeshFactory::create(Options* options, diff --git a/src/physics/physicsmodel.cxx b/src/physics/physicsmodel.cxx index 8936982561..85de541410 100644 --- a/src/physics/physicsmodel.cxx +++ b/src/physics/physicsmodel.cxx @@ -111,9 +111,16 @@ void PhysicsModel::initialise(Solver* s) { } } -int PhysicsModel::runRHS(BoutReal time, bool linear) { return rhs(time, linear); } +int PhysicsModel::runRHS_se(BoutReal time, bool linear) { return PhysicsModel::rhs_se(time, linear); } +int PhysicsModel::runRHS_si(BoutReal time, bool linear) { return PhysicsModel::rhs_si(time, linear); } +int PhysicsModel::runRHS_fe(BoutReal time, bool linear) { return PhysicsModel::rhs_fe(time, linear); } +int PhysicsModel::runRHS_fi(BoutReal time, bool linear) { return PhysicsModel::rhs_fi(time, linear); } +int PhysicsModel::runRHS_s(BoutReal time, bool linear) { return PhysicsModel::rhs_s(time, linear); } +int PhysicsModel::runRHS_f(BoutReal time, bool linear) { return PhysicsModel::rhs_f(time, linear); } +int PhysicsModel::runRHS(BoutReal time, bool linear) { return PhysicsModel::rhs(time, linear); } bool PhysicsModel::splitOperator() { return splitop; } +bool PhysicsModel::splitOperatorMRI() { return splitopmri; } int PhysicsModel::runConvective(BoutReal time, bool linear) { return convective(time, linear); diff --git a/src/solver/impls/arkode/arkode.cxx b/src/solver/impls/arkode/arkode.cxx index 440f8f54f1..34a99e2e3c 100644 --- a/src/solver/impls/arkode/arkode.cxx +++ b/src/solver/impls/arkode/arkode.cxx @@ -151,7 +151,7 @@ ArkodeSolver::ArkodeSolver(Options* opts) ArkodeSolver::~ArkodeSolver() { N_VDestroy(uvec); - ARKStepFree(&arkode_mem); + ARKodeFree(&arkode_mem); SUNLinSolFree(sun_solver); SUNNonlinSolFree(nonlinear_solver); @@ -221,13 +221,13 @@ int ArkodeSolver::init() { } break; case Treatment::Explicit: - output_info.write("\tUsing ARKStep Explicit solver \n"); + output_info.write("\tUsing ARKode Explicit solver \n"); if (ARKStepSetExplicit(arkode_mem) != ARK_SUCCESS) { throw BoutException("ARKStepSetExplicit failed\n"); } break; case Treatment::Implicit: - output_info.write("\tUsing ARKStep Implicit solver \n"); + output_info.write("\tUsing ARKode Implicit solver \n"); if (ARKStepSetImplicit(arkode_mem) != ARK_SUCCESS) { throw BoutException("ARKStepSetImplicit failed\n"); } @@ -237,24 +237,24 @@ int ArkodeSolver::init() { } // For callbacks, need pointer to solver object - if (ARKStepSetUserData(arkode_mem, this) != ARK_SUCCESS) { - throw BoutException("ARKStepSetUserData failed\n"); + if (ARKodeSetUserData(arkode_mem, this) != ARK_SUCCESS) { + throw BoutException("ARKodeSetUserData failed\n"); } - if (ARKStepSetLinear(arkode_mem, set_linear) != ARK_SUCCESS) { - throw BoutException("ARKStepSetLinear failed\n"); + if (ARKodeSetLinear(arkode_mem, set_linear) != ARK_SUCCESS) { + throw BoutException("ARKodeSetLinear failed\n"); } if (fixed_step) { // If not given, default to adaptive timestepping const auto fixed_timestep = (*options)["timestep"].withDefault(0.0); - if (ARKStepSetFixedStep(arkode_mem, fixed_timestep) != ARK_SUCCESS) { - throw BoutException("ARKStepSetFixedStep failed\n"); + if (ARKodeSetFixedStep(arkode_mem, fixed_timestep) != ARK_SUCCESS) { + throw BoutException("ARKodeSetFixedStep failed\n"); } } - if (ARKStepSetOrder(arkode_mem, order) != ARK_SUCCESS) { - throw BoutException("ARKStepSetOrder failed\n"); + if (ARKodeSetOrder(arkode_mem, order) != ARK_SUCCESS) { + throw BoutException("ARKodeSetOrder failed\n"); } #if SUNDIALS_TABLE_BY_NAME_SUPPORT @@ -269,8 +269,8 @@ int ArkodeSolver::init() { } #endif - if (ARKStepSetCFLFraction(arkode_mem, cfl_frac) != ARK_SUCCESS) { - throw BoutException("ARKStepSetCFLFraction failed\n"); + if (ARKodeSetCFLFraction(arkode_mem, cfl_frac) != ARK_SUCCESS) { + throw BoutException("ARKodeSetCFLFraction failed\n"); } #if SUNDIALS_CONTROLLER_SUPPORT @@ -297,12 +297,12 @@ int ArkodeSolver::init() { throw BoutException("Invalid adap_method\n"); } - if (ARKStepSetAdaptController(arkode_mem, controller) != ARK_SUCCESS) { - throw BoutException("ARKStepSetAdaptController failed\n"); + if (ARKodeSetAdaptController(arkode_mem, controller) != ARK_SUCCESS) { + throw BoutException("ARKodeSetAdaptController failed\n"); } - if (ARKStepSetAdaptivityAdjustment(arkode_mem, 0) != ARK_SUCCESS) { - throw BoutException("ARKStepSetAdaptivityAdjustment failed\n"); + if (ARKodeSetAdaptivityAdjustment(arkode_mem, 0) != ARK_SUCCESS) { + throw BoutException("ARKodeSetAdaptivityAdjustment failed\n"); } #else int adap_method_int; @@ -330,9 +330,9 @@ int ArkodeSolver::init() { throw BoutException("Invalid adap_method\n"); } - if (ARKStepSetAdaptivityMethod(arkode_mem, adap_method_int, 1, 1, nullptr) + if (ARKodeSetAdaptivityMethod(arkode_mem, adap_method_int, 1, 1, nullptr) != ARK_SUCCESS) { - throw BoutException("ARKStepSetAdaptivityMethod failed\n"); + throw BoutException("ARKodeSetAdaptivityMethod failed\n"); } #endif @@ -365,36 +365,36 @@ int ArkodeSolver::init() { set_abstol_values(N_VGetArrayPointer(abstolvec), f2dtols, f3dtols); - if (ARKStepSVtolerances(arkode_mem, reltol, abstolvec) != ARK_SUCCESS) { - throw BoutException("ARKStepSVtolerances failed\n"); + if (ARKodeSVtolerances(arkode_mem, reltol, abstolvec) != ARK_SUCCESS) { + throw BoutException("ARKodeSVtolerances failed\n"); } N_VDestroy(abstolvec); } else { - if (ARKStepSStolerances(arkode_mem, reltol, abstol) != ARK_SUCCESS) { - throw BoutException("ARKStepSStolerances failed\n"); + if (ARKodeSStolerances(arkode_mem, reltol, abstol) != ARK_SUCCESS) { + throw BoutException("ARKodeSStolerances failed\n"); } } - if (ARKStepSetMaxNumSteps(arkode_mem, mxsteps) != ARK_SUCCESS) { - throw BoutException("ARKStepSetMaxNumSteps failed\n"); + if (ARKodeSetMaxNumSteps(arkode_mem, mxsteps) != ARK_SUCCESS) { + throw BoutException("ARKodeSetMaxNumSteps failed\n"); } if (max_timestep > 0.0) { - if (ARKStepSetMaxStep(arkode_mem, max_timestep) != ARK_SUCCESS) { - throw BoutException("ARKStepSetMaxStep failed\n"); + if (ARKodeSetMaxStep(arkode_mem, max_timestep) != ARK_SUCCESS) { + throw BoutException("ARKodeSetMaxStep failed\n"); } } if (min_timestep > 0.0) { - if (ARKStepSetMinStep(arkode_mem, min_timestep) != ARK_SUCCESS) { - throw BoutException("ARKStepSetMinStep failed\n"); + if (ARKodeSetMinStep(arkode_mem, min_timestep) != ARK_SUCCESS) { + throw BoutException("ARKodeSetMinStep failed\n"); } } if (start_timestep > 0.0) { - if (ARKStepSetInitStep(arkode_mem, start_timestep) != ARK_SUCCESS) { - throw BoutException("ARKStepSetInitStep failed"); + if (ARKodeSetInitStep(arkode_mem, start_timestep) != ARK_SUCCESS) { + throw BoutException("ARKodeSetInitStep failed"); } } @@ -405,8 +405,8 @@ int ArkodeSolver::init() { if (nonlinear_solver == nullptr) { throw BoutException("Creating SUNDIALS fixed point nonlinear solver failed\n"); } - if (ARKStepSetNonlinearSolver(arkode_mem, nonlinear_solver) != ARK_SUCCESS) { - throw BoutException("ARKStepSetNonlinearSolver failed\n"); + if (ARKodeSetNonlinearSolver(arkode_mem, nonlinear_solver) != ARK_SUCCESS) { + throw BoutException("ARKodeSetNonlinearSolver failed\n"); } } else { output.write("\tUsing Newton iteration\n"); @@ -417,8 +417,8 @@ int ArkodeSolver::init() { if (sun_solver == nullptr) { throw BoutException("Creating SUNDIALS linear solver failed\n"); } - if (ARKStepSetLinearSolver(arkode_mem, sun_solver, nullptr) != ARKLS_SUCCESS) { - throw BoutException("ARKStepSetLinearSolver failed\n"); + if (ARKodeSetLinearSolver(arkode_mem, sun_solver, nullptr) != ARKLS_SUCCESS) { + throw BoutException("ARKodeSetLinearSolver failed\n"); } /// Set Preconditioner @@ -426,9 +426,9 @@ int ArkodeSolver::init() { if (hasPreconditioner()) { output.write("\tUsing user-supplied preconditioner\n"); - if (ARKStepSetPreconditioner(arkode_mem, nullptr, arkode_pre) + if (ARKodeSetPreconditioner(arkode_mem, nullptr, arkode_pre) != ARKLS_SUCCESS) { - throw BoutException("ARKStepSetPreconditioner failed\n"); + throw BoutException("ARKodeSetPreconditioner failed\n"); } } else { output.write("\tUsing BBD preconditioner\n"); @@ -480,8 +480,8 @@ int ArkodeSolver::init() { if (use_jacobian and hasJacobian()) { output.write("\tUsing user-supplied Jacobian function\n"); - if (ARKStepSetJacTimes(arkode_mem, nullptr, arkode_jac) != ARKLS_SUCCESS) { - throw BoutException("ARKStepSetJacTimes failed\n"); + if (ARKodeSetJacTimes(arkode_mem, nullptr, arkode_jac) != ARKLS_SUCCESS) { + throw BoutException("ARKodeSetJacTimes failed\n"); } } else { output.write("\tUsing difference quotient approximation for Jacobian\n"); @@ -523,17 +523,17 @@ int ArkodeSolver::run() { // Get additional diagnostics long int temp_long_int, temp_long_int2; - ARKStepGetNumSteps(arkode_mem, &temp_long_int); + ARKodeGetNumSteps(arkode_mem, &temp_long_int); nsteps = int(temp_long_int); ARKStepGetNumRhsEvals(arkode_mem, &temp_long_int, &temp_long_int2); nfe_evals = int(temp_long_int); nfi_evals = int(temp_long_int2); if (treatment == Treatment::ImEx or treatment == Treatment::Implicit) { - ARKStepGetNumNonlinSolvIters(arkode_mem, &temp_long_int); + ARKodeGetNumNonlinSolvIters(arkode_mem, &temp_long_int); nniters = int(temp_long_int); - ARKStepGetNumPrecEvals(arkode_mem, &temp_long_int); + ARKodeGetNumPrecEvals(arkode_mem, &temp_long_int); npevals = int(temp_long_int); - ARKStepGetNumLinIters(arkode_mem, &temp_long_int); + ARKodeGetNumLinIters(arkode_mem, &temp_long_int); nliters = int(temp_long_int); } @@ -571,15 +571,15 @@ BoutReal ArkodeSolver::run(BoutReal tout) { int flag; if (!monitor_timestep) { // Run in normal mode - flag = ARKStepEvolve(arkode_mem, tout, uvec, &simtime, ARK_NORMAL); + flag = ARKodeEvolve(arkode_mem, tout, uvec, &simtime, ARK_NORMAL); } else { // Run in single step mode, to call timestep monitors BoutReal internal_time; - ARKStepGetCurrentTime(arkode_mem, &internal_time); + ARKodeGetCurrentTime(arkode_mem, &internal_time); while (internal_time < tout) { // Run another step const BoutReal last_time = internal_time; - flag = ARKStepEvolve(arkode_mem, tout, uvec, &internal_time, ARK_ONE_STEP); + flag = ARKodeEvolve(arkode_mem, tout, uvec, &internal_time, ARK_ONE_STEP); if (flag != ARK_SUCCESS) { output_error.write("ERROR ARKODE solve failed at t = {:e}, flag = {:d}\n", @@ -591,7 +591,7 @@ BoutReal ArkodeSolver::run(BoutReal tout) { call_timestep_monitors(internal_time, internal_time - last_time); } // Get output at the desired time - flag = ARKStepGetDky(arkode_mem, tout, 0, uvec); + flag = ARKodeGetDky(arkode_mem, tout, 0, uvec); simtime = tout; } @@ -621,7 +621,7 @@ void ArkodeSolver::rhs_e(BoutReal t, BoutReal* udata, BoutReal* dudata) { // Get the current timestep // Note: ARKodeGetCurrentStep updated too late in older versions - ARKStepGetLastStep(arkode_mem, &hcur); + ARKodeGetLastStep(arkode_mem, &hcur); // Call RHS function run_convective(t); @@ -638,7 +638,7 @@ void ArkodeSolver::rhs_i(BoutReal t, BoutReal* udata, BoutReal* dudata) { TRACE("Running RHS: ArkodeSolver::rhs_i({:e})", t); load_vars(udata); - ARKStepGetLastStep(arkode_mem, &hcur); + ARKodeGetLastStep(arkode_mem, &hcur); // Call Implicit RHS function run_diffusive(t); save_derivs(dudata); @@ -651,7 +651,7 @@ void ArkodeSolver::rhs(BoutReal t, BoutReal* udata, BoutReal* dudata) { TRACE("Running RHS: ArkodeSolver::rhs({:e})", t); load_vars(udata); - ARKStepGetLastStep(arkode_mem, &hcur); + ARKodeGetLastStep(arkode_mem, &hcur); // Call Implicit RHS function run_rhs(t); save_derivs(dudata); diff --git a/src/solver/solver.cxx b/src/solver/solver.cxx index 24ef0dd319..2aad36dd5d 100644 --- a/src/solver/solver.cxx +++ b/src/solver/solver.cxx @@ -45,6 +45,7 @@ // Implementations: #include "impls/adams_bashforth/adams_bashforth.hxx" #include "impls/arkode/arkode.hxx" +#include "impls/arkode/arkode_mri.hxx" #include "impls/cvode/cvode.hxx" #include "impls/euler/euler.hxx" #include "impls/ida/ida.hxx" @@ -941,7 +942,32 @@ int Solver::resetRHSCounter_e() { return t; } +int Solver::resetRHSCounter_se() { + int t = rhs_ncalls_se; + rhs_ncalls_se = 0; + return t; +} + +int Solver::resetRHSCounter_si() { + int t = rhs_ncalls_si; + rhs_ncalls_si = 0; + return t; +} + +int Solver::resetRHSCounter_fe() { + int t = rhs_ncalls_fe; + rhs_ncalls_fe = 0; + return t; +} + +int Solver::resetRHSCounter_fi() { + int t = rhs_ncalls_fi; + rhs_ncalls_fi = 0; + return t; +} + bool Solver::splitOperator() { return model->splitOperator(); } +bool Solver::splitOperatorMRI() { return model->splitOperatorMRI(); } ///////////////////////////////////////////////////// @@ -1357,6 +1383,140 @@ Field3D Solver::globalIndex(int localStart) { * Running user-supplied functions **************************************************************************/ +int Solver::run_rhs_se(BoutReal t, bool linear) { + int status; + + Timer timer("rhs"); + + if (first_rhs_call) { + // Ensure that nonlinear terms are calculated on first call + linear = false; + first_rhs_call = false; + } + + pre_rhs(t); + status = model->runRHS_se(t, linear); + post_rhs(t); + + // If using Method of Manufactured Solutions + add_mms_sources(t); + + rhs_ncalls_se++; + return status; +} + +int Solver::run_rhs_si(BoutReal t, bool linear) { + int status; + + Timer timer("rhs"); + + if (first_rhs_call) { + // Ensure that nonlinear terms are calculated on first call + linear = false; + first_rhs_call = false; + } + + pre_rhs(t); + status = model->runRHS_si(t, linear); + post_rhs(t); + + // If using Method of Manufactured Solutions + add_mms_sources(t); + + rhs_ncalls_si++; + return status; +} + +int Solver::run_rhs_fe(BoutReal t, bool linear) { + int status; + + Timer timer("rhs"); + + if (first_rhs_call) { + // Ensure that nonlinear terms are calculated on first call + linear = false; + first_rhs_call = false; + } + + pre_rhs(t); + status = model->runRHS_fe(t, linear); + post_rhs(t); + + // If using Method of Manufactured Solutions + add_mms_sources(t); + + rhs_ncalls_fe++; + return status; +} + +int Solver::run_rhs_fi(BoutReal t, bool linear) { + int status; + + Timer timer("rhs"); + + if (first_rhs_call) { + // Ensure that nonlinear terms are calculated on first call + linear = false; + first_rhs_call = false; + } + + pre_rhs(t); + status = model->runRHS_fi(t, linear); + post_rhs(t); + + // If using Method of Manufactured Solutions + add_mms_sources(t); + + rhs_ncalls_fi++; + return status; +} + +int Solver::run_rhs_s(BoutReal t, bool linear) { + int status; + + Timer timer("rhs"); + + if (first_rhs_call) { + // Ensure that nonlinear terms are calculated on first call + linear = false; + first_rhs_call = false; + } + + pre_rhs(t); + status = model->runRHS_s(t, linear); + post_rhs(t); + + // If using Method of Manufactured Solutions + add_mms_sources(t); + + rhs_ncalls_se++; + rhs_ncalls_si++; + return status; +} + +int Solver::run_rhs_f(BoutReal t, bool linear) { + int status; + + Timer timer("rhs"); + + if (first_rhs_call) { + // Ensure that nonlinear terms are calculated on first call + linear = false; + first_rhs_call = false; + } + + pre_rhs(t); + status = model->runRHS_f(t, linear); + post_rhs(t); + + // If using Method of Manufactured Solutions + add_mms_sources(t); + + rhs_ncalls_fe++; + rhs_ncalls_fi++; + return status; +} + int Solver::run_rhs(BoutReal t, bool linear) { int status; @@ -1368,7 +1528,50 @@ int Solver::run_rhs(BoutReal t, bool linear) { first_rhs_call = false; } - if (model->splitOperator()) { + + // TO DO: Replace true with an appropriate boolean and check for efficiency + if (model->splitOperatorMRI()) { + // Run all four parts + + int nv = getLocalN(); + // Create temporary arrays for system state + Array tmp(nv); + Array tmp2(nv); + Array tmp3(nv); + + save_vars(tmp.begin()); // Copy variables into tmp + + pre_rhs(t); + status = model->runRHS_fe(t, linear); + post_rhs(t); // Check variables, apply boundary conditions + save_derivs(tmp2.begin()); // Save time derivatives + + pre_rhs(t); + status = model->runRHS_fi(t, linear); + post_rhs(t); + save_derivs(tmp3.begin()); // Save time derivatives + for (BoutReal *t3 = tmp3.begin(), *t2 = tmp2.begin(); t3 != tmp3.end(); ++t3, ++t2) { + *t3 += *t2; + } + + pre_rhs(t); + status = model->runRHS_se(t, linear); + post_rhs(t); + save_derivs(tmp2.begin()); // Save time derivatives + for (BoutReal *t3 = tmp3.begin(), *t2 = tmp2.begin(); t3 != tmp3.end(); ++t3, ++t2) { + *t3 += *t2; + } + + pre_rhs(t); + status = model->runRHS_si(t, linear); + post_rhs(t); + save_derivs(tmp2.begin()); // Save time derivatives + for (BoutReal *t3 = tmp3.begin(), *t2 = tmp2.begin(); t3 != tmp3.end(); ++t3, ++t2) { + *t3 += *t2; + } + load_derivs(tmp3.begin()); // Put back time-derivatives + } + else if (model->splitOperator()) { // Run both parts int nv = getLocalN(); @@ -1391,18 +1594,21 @@ int Solver::run_rhs(BoutReal t, bool linear) { *t += *t2; } load_derivs(tmp.begin()); // Put back time-derivatives + rhs_ncalls++; + rhs_ncalls_e++; + rhs_ncalls_i++; } else { pre_rhs(t); status = model->runRHS(t, linear); post_rhs(t); + rhs_ncalls++; + rhs_ncalls_e++; + rhs_ncalls_i++; } // If using Method of Manufactured Solutions add_mms_sources(t); - rhs_ncalls++; - rhs_ncalls_e++; - rhs_ncalls_i++; return status; } diff --git a/tests/integrated/CMakeLists.txt b/tests/integrated/CMakeLists.txt index ef173db7df..a67c6dbeea 100644 --- a/tests/integrated/CMakeLists.txt +++ b/tests/integrated/CMakeLists.txt @@ -21,6 +21,7 @@ add_subdirectory(test-interpolate) add_subdirectory(test-interpolate-z) add_subdirectory(test-invertable-operator) add_subdirectory(test-invpar) +add_subdirectory(test-kpr_mri) add_subdirectory(test-laplace) add_subdirectory(test-laplace-petsc3d) add_subdirectory(test-laplace-hypre3d) diff --git a/tests/unit/test_extras.hxx b/tests/unit/test_extras.hxx index 8db2d0aba4..700b977ac8 100644 --- a/tests/unit/test_extras.hxx +++ b/tests/unit/test_extras.hxx @@ -12,7 +12,6 @@ #include "bout/boutcomm.hxx" #include "bout/coordinates.hxx" #include "bout/field3d.hxx" -#include "bout/griddata.hxx" #include "bout/mesh.hxx" #include "bout/mpi_wrapper.hxx" #include "bout/operatorstencil.hxx" From c14dcd4740e767cb62ff5634d6af367ea4fc97ce Mon Sep 17 00:00:00 2001 From: maggul Date: Fri, 25 Oct 2024 11:17:48 -0700 Subject: [PATCH 02/28] add missing files --- src/solver/impls/arkode/arkode_mri.cxx | 1090 +++++++++++++++++ src/solver/impls/arkode/arkode_mri.hxx | 176 +++ tests/integrated/test-kpr_mri/CMakeLists.txt | 2 + tests/integrated/test-kpr_mri/README.md | 37 + tests/integrated/test-kpr_mri/makefile | 5 + tests/integrated/test-kpr_mri/runtest | 23 + .../integrated/test-kpr_mri/test_kpr_mri.cxx | 183 +++ 7 files changed, 1516 insertions(+) create mode 100644 src/solver/impls/arkode/arkode_mri.cxx create mode 100644 src/solver/impls/arkode/arkode_mri.hxx create mode 100644 tests/integrated/test-kpr_mri/CMakeLists.txt create mode 100644 tests/integrated/test-kpr_mri/README.md create mode 100644 tests/integrated/test-kpr_mri/makefile create mode 100755 tests/integrated/test-kpr_mri/runtest create mode 100644 tests/integrated/test-kpr_mri/test_kpr_mri.cxx diff --git a/src/solver/impls/arkode/arkode_mri.cxx b/src/solver/impls/arkode/arkode_mri.cxx new file mode 100644 index 0000000000..1fb9a79a24 --- /dev/null +++ b/src/solver/impls/arkode/arkode_mri.cxx @@ -0,0 +1,1090 @@ +/************************************************************************** + * Experimental interface to SUNDIALS ARKode MRI solver + * + * NOTE: ARKode is still in beta testing so use with cautious optimism + * + ************************************************************************** + * Copyright 2010-2024 BOUT++ contributors + * + * This file is part of BOUT++. + * + * BOUT++ is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * BOUT++ is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with BOUT++. If not, see . + * + **************************************************************************/ + +#include "bout/build_config.hxx" + +#include "arkode_mri.hxx" + +#if BOUT_HAS_ARKODE + +#include "bout/bout_enum_class.hxx" +#include "bout/boutcomm.hxx" +#include "bout/boutexception.hxx" +#include "bout/field3d.hxx" +#include "bout/mesh.hxx" +#include "bout/msg_stack.hxx" +#include "bout/options.hxx" +#include "bout/output.hxx" +#include "bout/unused.hxx" +#include "bout/utils.hxx" + +#include +#include +#include +#include + +#include +#include + +class Field2D; + +// NOLINTBEGIN(readability-identifier-length) +namespace { +int arkode_rhs_s_explicit(BoutReal t, N_Vector u, N_Vector du, void* user_data); +int arkode_rhs_s_implicit(BoutReal t, N_Vector u, N_Vector du, void* user_data); +int arkode_rhs_f_explicit(BoutReal t, N_Vector u, N_Vector du, void* user_data); +int arkode_rhs_f_implicit(BoutReal t, N_Vector u, N_Vector du, void* user_data); +int arkode_s_rhs(BoutReal t, N_Vector u, N_Vector du, void* user_data); +int arkode_f_rhs(BoutReal t, N_Vector u, N_Vector du, void* user_data); + +int arkode_s_bbd_rhs(sunindextype Nlocal, BoutReal t, N_Vector u, N_Vector du, + void* user_data); +int arkode_f_bbd_rhs(sunindextype Nlocal, BoutReal t, N_Vector u, N_Vector du, + void* user_data); +int arkode_s_pre(BoutReal t, N_Vector yy, N_Vector yp, N_Vector rvec, N_Vector zvec, + BoutReal gamma, BoutReal delta, int lr, void* user_data); +int arkode_f_pre(BoutReal t, N_Vector yy, N_Vector yp, N_Vector rvec, N_Vector zvec, + BoutReal gamma, BoutReal delta, int lr, void* user_data); + +} // namespace +// NOLINTEND(readability-identifier-length) + +ArkodeMRISolver::ArkodeMRISolver(Options* opts) + : Solver(opts), diagnose((*options)["diagnose"] + .doc("Print some additional diagnostics") + .withDefault(false)), + mxsteps((*options)["mxstep"] + .doc("Maximum number of steps to take between outputs") + .withDefault(500)), + treatment((*options)["treatment"] + .doc("Use default capability (imex) or provide a specific treatment: " + "implicit or explicit") + .withDefault(MRI_Treatment::ImEx)), + inner_treatment((*options)["inner_treatment"] + .doc("Use default capability (imex) or provide a specific inner_treatment: " + "implicit or explicit") + .withDefault(MRI_Treatment::ImEx)), + set_linear( + (*options)["set_linear"] + .doc("Use linear implicit solver (only evaluates jacobian inversion once)") + .withDefault(false)), + inner_set_linear( + (*options)["inner_set_linear"] + .doc("Use linear implicit solver (only evaluates jacobian inversion once)") + .withDefault(false)), + fixed_step((*options)["fixed_step"] + .doc("Solve explicit portion in fixed timestep mode. NOTE: This is " + "not recommended except for code comparison") + .withDefault(false)), + order((*options)["order"].doc("Order of internal step").withDefault(4)), + adap_method( + (*options)["adap_method"] + .doc("Set timestep adaptivity function: pid, pi, i, explicit_gustafsson, " + "implicit_gustafsson, imex_gustafsson.") + .withDefault(MRI_AdapMethod::PID)), + abstol((*options)["atol"].doc("Absolute tolerance").withDefault(1.0e-12)), + reltol((*options)["rtol"].doc("Relative tolerance").withDefault(1.0e-5)), + use_vector_abstol((*options)["use_vector_abstol"] + .doc("Use separate absolute tolerance for each field") + .withDefault(false)), + use_precon((*options)["use_precon"] + .doc("Use user-supplied preconditioner function") + .withDefault(false)), + inner_use_precon((*options)["inner_use_precon"] + .doc("Use user-supplied preconditioner function") + .withDefault(false)), + maxl( + (*options)["maxl"].doc("Number of Krylov basis vectors to use").withDefault(0)), + inner_maxl( + (*options)["inner_maxl"].doc("Number of Krylov basis vectors to use").withDefault(0)), + rightprec((*options)["rightprec"] + .doc("Use right preconditioning instead of left preconditioning") + .withDefault(false)), + suncontext(createSUNContext(BoutComm::get())) { + has_constraints = false; // This solver doesn't have constraints + + // Add diagnostics to output + add_int_diagnostic(nsteps, "arkode_nsteps", "Cumulative number of internal steps"); + add_int_diagnostic(inner_nsteps, "arkode_inner_nsteps", "Cumulative number of inner internal steps"); + add_int_diagnostic(nfe_evals, "arkode_nfe_evals", + "No. of calls to fe (explicit portion of the right-hand-side " + "function) function"); + add_int_diagnostic(inner_nfe_evals, "arkode_inner_nfe_evals", + "No. of calls to fe (explicit portion of the inner right-hand-side " + "function) function"); + add_int_diagnostic(nfi_evals, "arkode_nfi_evals", + "No. of calls to fi (implicit portion of the right-hand-side " + "function) function"); + add_int_diagnostic(inner_nfi_evals, "arkode_inner_nfi_evals", + "No. of calls to fi (implicit portion of inner the right-hand-side " + "function) function"); + add_int_diagnostic(nniters, "arkode_nniters", "No. of nonlinear solver iterations"); + add_int_diagnostic(inner_nniters, "arkode_inner_nniters", "No. of inner nonlinear solver iterations"); + add_int_diagnostic(npevals, "arkode_npevals", "No. of preconditioner evaluations"); + add_int_diagnostic(inner_npevals, "arkode_inner_npevals", "No. of inner preconditioner evaluations"); + add_int_diagnostic(nliters, "arkode_nliters", "No. of linear iterations"); + add_int_diagnostic(inner_nliters, "arkode_inner_nliters", "No. of inner linear iterations"); +} + +ArkodeMRISolver::~ArkodeMRISolver() { + N_VDestroy(uvec); + ARKodeFree(&arkode_mem); + ARKodeFree(&inner_arkode_mem); + SUNLinSolFree(sun_solver); + SUNLinSolFree(inner_sun_solver); + SUNNonlinSolFree(nonlinear_solver); + SUNNonlinSolFree(inner_nonlinear_solver); + MRIStepInnerStepper_Free(&inner_stepper); + + SUNAdaptController_Destroy(controller); + SUNAdaptController_Destroy(inner_controller); +} + +/************************************************************************** + * Initialise + **************************************************************************/ + +int ArkodeMRISolver::init() { + TRACE("Initialising ARKODE MRI solver"); + + Solver::init(); + + output.write("Initialising SUNDIALS' ARKODE MRI solver\n"); + + // Calculate number of variables (in generic_solver) + const int local_N = getLocalN(); + + // Get total problem size + int neq; + if (bout::globals::mpi->MPI_Allreduce(&local_N, &neq, 1, MPI_INT, MPI_SUM, + BoutComm::get())) { + throw BoutException("Allreduce localN -> GlobalN failed!\n"); + } + + output.write("\t3d fields = {:d}, 2d fields = {:d} neq={:d}, local_N={:d}\n", n3Dvars(), + n2Dvars(), neq, local_N); + + // Allocate memory + uvec = callWithSUNContext(N_VNew_Parallel, suncontext, BoutComm::get(), local_N, neq); + if (uvec == nullptr) { + throw BoutException("SUNDIALS memory allocation failed\n"); + } + + // Put the variables into uvec + save_vars(N_VGetArrayPointer(uvec)); + + switch (inner_treatment) { + case MRI_Treatment::ImEx: + inner_arkode_mem = callWithSUNContext(ARKStepCreate, suncontext, arkode_rhs_f_explicit, + arkode_rhs_f_implicit, simtime, uvec); + output_info.write("\tUsing ARKode ImEx inner solver \n"); + break; + case MRI_Treatment::Explicit: + inner_arkode_mem = + callWithSUNContext(ARKStepCreate, suncontext, arkode_f_rhs, nullptr, simtime, uvec); + output_info.write("\tUsing ARKode Explicit inner solver \n"); + break; + case MRI_Treatment::Implicit: + inner_arkode_mem = + callWithSUNContext(ARKStepCreate, suncontext, nullptr, arkode_f_rhs, simtime, uvec); + output_info.write("\tUsing ARKode Implicit inner solver \n"); + break; + default: + throw BoutException("Invalid inner_treatment: {}\n", toString(inner_treatment)); + } + if (inner_arkode_mem == nullptr) { + throw BoutException("ARKStepCreate failed\n"); + } + + // For callbacks, need pointer to solver object + if (ARKodeSetUserData(inner_arkode_mem, this) != ARK_SUCCESS) { + throw BoutException("ARKodeSetUserData failed\n"); + } + + if(inner_treatment != MRI_Treatment::Explicit) + if (ARKodeSetLinear(inner_arkode_mem, inner_set_linear) != ARK_SUCCESS) { + throw BoutException("ARKodeSetLinear failed\n"); + } + + if (fixed_step) { + // If not given, default to adaptive timestepping + const BoutReal inner_fixed_timestep = (*options)["inner_timestep"]; + if (ARKodeSetFixedStep(inner_arkode_mem, inner_fixed_timestep) != ARK_SUCCESS) { + throw BoutException("ARKodeSetFixedStep failed\n"); + } + } + + if (ARKodeSetOrder(inner_arkode_mem, order) != ARK_SUCCESS) { + throw BoutException("ARKodeSetOrder failed\n"); + } + + if (ARKStepCreateMRIStepInnerStepper(inner_arkode_mem, &inner_stepper) != ARK_SUCCESS) { + throw BoutException("ARKStepCreateMRIStepInnerStepper failed\n"); + } + + // Initialize the slow integrator. Specify the explicit slow right-hand side + // function in y'=fe(t,y)+fi(t,y)+ff(t,y), the inital time T0, the + // initial dependent variable vector y, and the fast integrator. + + switch (treatment) { + case MRI_Treatment::ImEx: + arkode_mem = callWithSUNContext(MRIStepCreate, suncontext, arkode_rhs_s_explicit, arkode_rhs_s_implicit, + simtime, uvec, inner_stepper); + output_info.write("\tUsing ARKode ImEx solver \n"); + break; + case MRI_Treatment::Explicit: + arkode_mem = callWithSUNContext(MRIStepCreate, suncontext, arkode_s_rhs, nullptr, + simtime, uvec, inner_stepper); + output_info.write("\tUsing ARKode Explicit solver \n"); + break; + case MRI_Treatment::Implicit: + arkode_mem = callWithSUNContext(MRIStepCreate, suncontext, nullptr, arkode_s_rhs, + simtime, uvec, inner_stepper); + output_info.write("\tUsing ARKode Implicit solver \n"); + break; + default: + throw BoutException("Invalid treatment: {}\n", toString(treatment)); + } + if (arkode_mem == nullptr) { + throw BoutException("MRIStepCreate failed\n"); + } + + // For callbacks, need pointer to solver object + if (ARKodeSetUserData(arkode_mem, this) != ARK_SUCCESS) { + throw BoutException("ARKodeSetUserData failed\n"); + } + + if(treatment != MRI_Treatment::Explicit) + if (ARKodeSetLinear(arkode_mem, set_linear) != ARK_SUCCESS) { + throw BoutException("ARKodeSetLinear failed\n"); + } + + if (fixed_step) { + // If not given, default to adaptive timestepping + const BoutReal fixed_timestep = (*options)["timestep"]; + if (ARKodeSetFixedStep(arkode_mem, fixed_timestep) != ARK_SUCCESS) { + throw BoutException("ARKodeSetFixedStep failed\n"); + } + } + + if (ARKodeSetOrder(arkode_mem, order) != ARK_SUCCESS) { + throw BoutException("ARKodeSetOrder failed\n"); + } + + switch (adap_method) { + case MRI_AdapMethod::PID: + controller = SUNAdaptController_PID(suncontext); + inner_controller = SUNAdaptController_PID(suncontext); + break; + case MRI_AdapMethod::PI: + controller = SUNAdaptController_PI(suncontext); + inner_controller = SUNAdaptController_PI(suncontext); + break; + case MRI_AdapMethod::I: + controller = SUNAdaptController_I(suncontext); + inner_controller = SUNAdaptController_I(suncontext); + break; + case MRI_AdapMethod::Explicit_Gustafsson: + controller = SUNAdaptController_ExpGus(suncontext); + inner_controller = SUNAdaptController_ExpGus(suncontext); + break; + case MRI_AdapMethod::Implicit_Gustafsson: + controller = SUNAdaptController_ImpGus(suncontext); + inner_controller = SUNAdaptController_ImpGus(suncontext); + break; + case MRI_AdapMethod::ImEx_Gustafsson: + controller = SUNAdaptController_ImExGus(suncontext); + inner_controller = SUNAdaptController_ImExGus(suncontext); + break; + default: + throw BoutException("Invalid adap_method\n"); + } + + // if (ARKodeSetAdaptController(arkode_mem, controller) != ARK_SUCCESS) { + // throw BoutException("ARKodeSetAdaptController failed\n"); + // } + + // if (ARKodeSetAdaptivityAdjustment(arkode_mem, 0) != ARK_SUCCESS) { + // throw BoutException("ARKodeSetAdaptivityAdjustment failed\n"); + // } + + if (ARKodeSetFixedStep(arkode_mem, 0.001) != ARK_SUCCESS) { + throw BoutException("ARKodeSetAdaptController failed\n"); + } + + if (ARKodeSetAdaptController(inner_arkode_mem, controller) != ARK_SUCCESS) { + throw BoutException("ARKodeSetAdaptController failed\n"); + } + + if (ARKodeSetAdaptivityAdjustment(inner_arkode_mem, 0) != ARK_SUCCESS) { + throw BoutException("ARKodeSetAdaptivityAdjustment failed\n"); + } + + if (use_vector_abstol) { + std::vector f2dtols; + f2dtols.reserve(f2d.size()); + std::transform(begin(f2d), end(f2d), std::back_inserter(f2dtols), + [abstol = abstol](const VarStr& f2) { + auto& f2_options = Options::root()[f2.name]; + const auto wrong_name = f2_options.isSet("abstol"); + if (wrong_name) { + output_warn << "WARNING: Option 'abstol' for field " << f2.name + << " is deprecated. Please use 'atol' instead\n"; + } + const std::string atol_name = wrong_name ? "abstol" : "atol"; + return f2_options[atol_name].withDefault(abstol); + }); + + std::vector f3dtols; + f3dtols.reserve(f3d.size()); + std::transform(begin(f3d), end(f3d), std::back_inserter(f3dtols), + [abstol = abstol](const VarStr& f3) { + return Options::root()[f3.name]["atol"].withDefault(abstol); + }); + + N_Vector abstolvec = N_VClone(uvec); + if (abstolvec == nullptr) { + throw BoutException("SUNDIALS memory allocation (abstol vector) failed\n"); + } + + set_abstol_values(N_VGetArrayPointer(abstolvec), f2dtols, f3dtols); + + if (ARKodeSVtolerances(arkode_mem, reltol, abstolvec) != ARK_SUCCESS) { + throw BoutException("ARKodeSVtolerances failed\n"); + } + if (ARKodeSVtolerances(inner_arkode_mem, reltol, abstolvec) != ARK_SUCCESS) { + throw BoutException("ARKodeSVtolerances failed\n"); + } + + N_VDestroy(abstolvec); + } else { + if (ARKodeSStolerances(arkode_mem, reltol, abstol) != ARK_SUCCESS) { + throw BoutException("ARKodeSStolerances failed\n"); + } + if (ARKodeSStolerances(inner_arkode_mem, reltol, abstol) != ARK_SUCCESS) { + throw BoutException("ARKodeSStolerances failed\n"); + } + } + + if (ARKodeSetMaxNumSteps(arkode_mem, mxsteps) != ARK_SUCCESS) { + throw BoutException("ARKodeSetMaxNumSteps failed\n"); + } + if (ARKodeSetMaxNumSteps(inner_arkode_mem, mxsteps) != ARK_SUCCESS) { + throw BoutException("ARKodeSetMaxNumSteps failed\n"); + } + + if (inner_treatment == MRI_Treatment::ImEx or inner_treatment == MRI_Treatment::Implicit) { + { + output.write("\tUsing Newton iteration for inner solver\n"); + + const auto prectype = + inner_use_precon ? (rightprec ? SUN_PREC_RIGHT : SUN_PREC_LEFT) : SUN_PREC_NONE; + inner_sun_solver = callWithSUNContext(SUNLinSol_SPGMR, suncontext, uvec, prectype, inner_maxl); + if (inner_sun_solver == nullptr) { + throw BoutException("Creating SUNDIALS inner linear solver failed\n"); + } + if (ARKodeSetLinearSolver(inner_arkode_mem, inner_sun_solver, nullptr) != ARKLS_SUCCESS) { + throw BoutException("ARKodeSetLinearSolver failed for inner solver\n"); + } + + /// Set Preconditioner + if (inner_use_precon) { + if (hasPreconditioner()) { // change to inner_hasPreconditioner when it is available + output.write("\tUsing user-supplied preconditioner for inner solver\n"); + + if (ARKodeSetPreconditioner(inner_arkode_mem, nullptr, arkode_f_pre) + != ARKLS_SUCCESS) { + throw BoutException("ARKodeSetPreconditioner failed for inner solver\n"); + } + } else { + output.write("\tUsing BBD preconditioner for inner solver\n"); + + /// Get options + // Compute band_width_default from actually added fields, to allow for multiple + // Mesh objects + // + // Previous implementation was equivalent to: + // int MXSUB = mesh->xend - mesh->xstart + 1; + // int band_width_default = n3Dvars()*(MXSUB+2); + const int band_width_default = std::accumulate( + begin(f3d), end(f3d), 0, [](int acc, const VarStr& fvar) { + Mesh* localmesh = fvar.var->getMesh(); + return acc + localmesh->xend - localmesh->xstart + 3; + }); + + const auto mudq = (*options)["mudq"] + .doc("Upper half-bandwidth to be used in the difference " + "quotient Jacobian approximation") + .withDefault(band_width_default); + const auto mldq = (*options)["mldq"] + .doc("Lower half-bandwidth to be used in the difference " + "quotient Jacobian approximation") + .withDefault(band_width_default); + const auto mukeep = (*options)["mukeep"] + .doc("Upper half-bandwidth of the retained banded " + "approximate Jacobian block") + .withDefault(n3Dvars() + n2Dvars()); + const auto mlkeep = (*options)["mlkeep"] + .doc("Lower half-bandwidth of the retained banded " + "approximate Jacobian block") + .withDefault(n3Dvars() + n2Dvars()); + + if (ARKBBDPrecInit(inner_arkode_mem, local_N, mudq, mldq, mukeep, mlkeep, 0, + arkode_f_bbd_rhs, nullptr) + != ARKLS_SUCCESS) { + throw BoutException("ARKBBDPrecInit failed for inner solver\n"); + } + } + } else { + // Not using preconditioning + output.write("\tNo inner preconditioning\n"); + } + } + + /// Set Jacobian-vector multiplication function + output.write("\tUsing difference quotient approximation for Jacobian in the inner solver\n"); + } + + + if (treatment == MRI_Treatment::ImEx or treatment == MRI_Treatment::Implicit) { + { + output.write("\tUsing Newton iteration\n"); + + const auto prectype = + use_precon ? (rightprec ? SUN_PREC_RIGHT : SUN_PREC_LEFT) : SUN_PREC_NONE; + sun_solver = callWithSUNContext(SUNLinSol_SPGMR, suncontext, uvec, prectype, maxl); + if (sun_solver == nullptr) { + throw BoutException("Creating SUNDIALS linear solver failed\n"); + } + if (ARKodeSetLinearSolver(arkode_mem, sun_solver, nullptr) != ARKLS_SUCCESS) { + throw BoutException("ARKodeSetLinearSolver failed\n"); + } + + /// Set Preconditioner + if (use_precon) { + if (hasPreconditioner()) { + output.write("\tUsing user-supplied preconditioner\n"); + + if (ARKodeSetPreconditioner(arkode_mem, nullptr, arkode_s_pre) + != ARKLS_SUCCESS) { + throw BoutException("ARKodeSetPreconditioner failed\n"); + } + } else { + output.write("\tUsing BBD preconditioner\n"); + + /// Get options + // Compute band_width_default from actually added fields, to allow for multiple + // Mesh objects + // + // Previous implementation was equivalent to: + // int MXSUB = mesh->xend - mesh->xstart + 1; + // int band_width_default = n3Dvars()*(MXSUB+2); + const int band_width_default = std::accumulate( + begin(f3d), end(f3d), 0, [](int acc, const VarStr& fvar) { + Mesh* localmesh = fvar.var->getMesh(); + return acc + localmesh->xend - localmesh->xstart + 3; + }); + + const auto mudq = (*options)["mudq"] + .doc("Upper half-bandwidth to be used in the difference " + "quotient Jacobian approximation") + .withDefault(band_width_default); + const auto mldq = (*options)["mldq"] + .doc("Lower half-bandwidth to be used in the difference " + "quotient Jacobian approximation") + .withDefault(band_width_default); + const auto mukeep = (*options)["mukeep"] + .doc("Upper half-bandwidth of the retained banded " + "approximate Jacobian block") + .withDefault(n3Dvars() + n2Dvars()); + const auto mlkeep = (*options)["mlkeep"] + .doc("Lower half-bandwidth of the retained banded " + "approximate Jacobian block") + .withDefault(n3Dvars() + n2Dvars()); + + if (ARKBBDPrecInit(arkode_mem, local_N, mudq, mldq, mukeep, mlkeep, 0, + arkode_s_bbd_rhs, nullptr) + != ARKLS_SUCCESS) { + throw BoutException("ARKBBDPrecInit failed\n"); + } + } + } else { + // Not using preconditioning + output.write("\tNo preconditioning\n"); + } + } + + /// Set Jacobian-vector multiplication function + output.write("\tUsing difference quotient approximation for Jacobian\n"); + } + + return 0; +} + +/************************************************************************** + * Run - Advance time + **************************************************************************/ + +int ArkodeMRISolver::run() { + TRACE("ArkodeMRISolver::run()"); + + if (!initialised) { + throw BoutException("ArkodeMRISolver not initialised\n"); + } + + for (int i = 0; i < getNumberOutputSteps(); i++) { + + /// Run the solver for one output timestep + simtime = run(simtime + getOutputTimestep()); + + /// Check if the run succeeded + if (simtime < 0.0) { + // Step failed + output.write("Timestep failed. Aborting\n"); + + throw BoutException("ARKode timestep failed\n"); + } + + // Get additional diagnostics + long int temp_long_int, temp_long_int2; + ARKodeGetNumSteps(arkode_mem, &temp_long_int); + nsteps = int(temp_long_int); + MRIStepGetNumRhsEvals(arkode_mem, &temp_long_int, &temp_long_int2); //Change after the release + nfe_evals = int(temp_long_int); + nfi_evals = int(temp_long_int2); + if (treatment == MRI_Treatment::ImEx or treatment == MRI_Treatment::Implicit) { + ARKodeGetNumNonlinSolvIters(arkode_mem, &temp_long_int); + nniters = int(temp_long_int); + ARKodeGetNumPrecEvals(arkode_mem, &temp_long_int); + npevals = int(temp_long_int); + ARKodeGetNumLinIters(arkode_mem, &temp_long_int); + nliters = int(temp_long_int); + } + + ARKodeGetNumSteps(inner_arkode_mem, &temp_long_int); + inner_nsteps = int(temp_long_int); + ARKStepGetNumRhsEvals(inner_arkode_mem, &temp_long_int, &temp_long_int2); //Change after the release + inner_nfe_evals = int(temp_long_int); + inner_nfi_evals = int(temp_long_int2); + if (inner_treatment == MRI_Treatment::ImEx or inner_treatment == MRI_Treatment::Implicit) { + ARKodeGetNumNonlinSolvIters(inner_arkode_mem, &temp_long_int); + inner_nniters = int(temp_long_int); + ARKodeGetNumPrecEvals(inner_arkode_mem, &temp_long_int); + inner_npevals = int(temp_long_int); + ARKodeGetNumLinIters(inner_arkode_mem, &temp_long_int); + inner_nliters = int(temp_long_int); + } + + if (diagnose) { + output.write("\nARKODE: nsteps {:d}, nfe_evals {:d}, nfi_evals {:d}, nniters {:d}, " + "npevals {:d}, nliters {:d}\n", + nsteps, nfe_evals, nfi_evals, nniters, npevals, nliters); + if (treatment == MRI_Treatment::ImEx or treatment == MRI_Treatment::Implicit) { + output.write(" -> Newton iterations per step: {:e}\n", + static_cast(nniters) / static_cast(nsteps)); + output.write(" -> Linear iterations per Newton iteration: {:e}\n", + static_cast(nliters) / static_cast(nniters)); + output.write(" -> Preconditioner evaluations per Newton: {:e}\n", + static_cast(npevals) / static_cast(nniters)); + } + + output.write("\nARKODE Inner: inner_nsteps {:d}, inner_nfe_evals {:d}, inner_nfi_evals {:d}, inner_nniters {:d}, " + "inner_npevals {:d}, inner_nliters {:d}\n", + inner_nsteps, inner_nfe_evals, inner_nfi_evals, inner_nniters, inner_npevals, inner_nliters); + if (inner_treatment == MRI_Treatment::ImEx or inner_treatment == MRI_Treatment::Implicit) { + output.write(" -> Inner Newton iterations per step: {:e}\n", + static_cast(inner_nniters) / static_cast(inner_nsteps)); + output.write(" -> Inner Linear iterations per Newton iteration: {:e}\n", + static_cast(inner_nliters) / static_cast(inner_nniters)); + output.write(" -> Inner Preconditioner evaluations per Newton: {:e}\n", + static_cast(inner_npevals) / static_cast(inner_nniters)); + } + } + + if (call_monitors(simtime, i, getNumberOutputSteps())) { + // User signalled to quit + break; + } + } + + return 0; +} + +BoutReal ArkodeMRISolver::run(BoutReal tout) { + TRACE("Running solver: solver::run({:e})", tout); + + bout::globals::mpi->MPI_Barrier(BoutComm::get()); + + pre_Wtime_s = 0.0; + pre_ncalls_s = 0; + + int flag; + if (!monitor_timestep) { + // Run in normal mode + flag = ARKodeEvolve(arkode_mem, tout, uvec, &simtime, ARK_NORMAL); + } else { + // Run in single step mode, to call timestep monitors + BoutReal internal_time; + ARKodeGetCurrentTime(arkode_mem, &internal_time); + while (internal_time < tout) { + // Run another step + const BoutReal last_time = internal_time; + flag = ARKodeEvolve(arkode_mem, tout, uvec, &internal_time, ARK_ONE_STEP); + + if (flag != ARK_SUCCESS) { + output_error.write("ERROR ARKODE solve failed at t = {:e}, flag = {:d}\n", + internal_time, flag); + return -1.0; + } + + // Call timestep monitor + call_timestep_monitors(internal_time, internal_time - last_time); + } + // Get output at the desired time + flag = ARKodeGetDky(arkode_mem, tout, 0, uvec); + simtime = tout; + } + + // Copy variables + load_vars(N_VGetArrayPointer(uvec)); + // Call rhs function to get extra variables at this time + run_rhs(simtime); + // run_diffusive(simtime); + if (flag != ARK_SUCCESS) { + output_error.write("ERROR ARKODE solve failed at t = {:e}, flag = {:d}\n", simtime, + flag); + return -1.0; + } + + return simtime; +} + +/************************************************************************** + * Explicit RHS function du = F^s_E(t, u) + **************************************************************************/ + +void ArkodeMRISolver::rhs_se(BoutReal t, BoutReal* udata, BoutReal* dudata) { + TRACE("Running RHS: ArkodeMRISolver::rhs_e({:e})", t); + + // Load state from udata + load_vars(udata); + + // Get the current timestep + // Note: ARKodeGetCurrentStep updated too late in older versions + ARKodeGetLastStep(arkode_mem, &hcur); + + // Call RHS function + run_rhs_se(t); + + // Save derivatives to dudata + save_derivs(dudata); +} + +/************************************************************************** + * Implicit RHS function du = F^s_I(t, u) + **************************************************************************/ + +void ArkodeMRISolver::rhs_si(BoutReal t, BoutReal* udata, BoutReal* dudata) { + TRACE("Running RHS: ArkodeMRISolver::rhs_si({:e})", t); + + load_vars(udata); + ARKodeGetLastStep(arkode_mem, &hcur); + // Call Implicit RHS function + run_rhs_si(t); + save_derivs(dudata); +} + +/************************************************************************** + * Explicit RHS function du = F^f_E(t, u) + **************************************************************************/ + +void ArkodeMRISolver::rhs_fe(BoutReal t, BoutReal* udata, BoutReal* dudata) { + TRACE("Running RHS: ArkodeMRISolver::rhs_e({:e})", t); + + // Load state from udata + load_vars(udata); + + // Get the current timestep + // Note: ARKodeGetCurrentStep updated too late in older versions + ARKodeGetLastStep(arkode_mem, &hcur); + + // Call RHS function + run_rhs_fe(t); + + // Save derivatives to dudata + save_derivs(dudata); +} + +/************************************************************************** + * Implicit RHS function du = F^f_I(t, u) + **************************************************************************/ + +void ArkodeMRISolver::rhs_fi(BoutReal t, BoutReal* udata, BoutReal* dudata) { + TRACE("Running RHS: ArkodeMRISolver::rhs_si({:e})", t); + + load_vars(udata); + ARKodeGetLastStep(arkode_mem, &hcur); + // Call Implicit RHS function + run_rhs_fi(t); + save_derivs(dudata); +} + +/************************************************************************** + * Slow RHS function du = F^s(t, u) + **************************************************************************/ + +void ArkodeMRISolver::rhs_s(BoutReal t, BoutReal* udata, BoutReal* dudata) { + TRACE("Running RHS: ArkodeMRISolver::rhs_e({:e})", t); + + // Load state from udata + load_vars(udata); + + // Get the current timestep + // Note: ARKodeGetCurrentStep updated too late in older versions + ARKodeGetLastStep(arkode_mem, &hcur); + + // Call RHS function + // run_rhs_s(t); + run_rhs_s(t); + + // Save derivatives to dudata + save_derivs(dudata); +} + +/************************************************************************** + * Fast RHS function du = F^f(t, u) + **************************************************************************/ + +void ArkodeMRISolver::rhs_f(BoutReal t, BoutReal* udata, BoutReal* dudata) { + TRACE("Running RHS: ArkodeMRISolver::rhs_e({:e})", t); + + // Load state from udata + load_vars(udata); + + // Get the current timestep + // Note: ARKodeGetCurrentStep updated too late in older versions + ARKodeGetLastStep(arkode_mem, &hcur); + + // Call RHS function + // run_rhs_f(t); + run_rhs_f(t); + + // Save derivatives to dudata + save_derivs(dudata); +} + +/************************************************************************** + * Preconditioner functions + **************************************************************************/ + +void ArkodeMRISolver::pre_s(BoutReal t, BoutReal gamma, BoutReal delta, BoutReal* udata, + BoutReal* rvec, BoutReal* zvec) { + TRACE("Running preconditioner: ArkodeMRISolver::pre({:e})", t); + + const BoutReal tstart = bout::globals::mpi->MPI_Wtime(); + + if (!hasPreconditioner()) { + // Identity (but should never happen) + const auto length = N_VGetLocalLength_Parallel(uvec); + std::copy(rvec, rvec + length, zvec); + return; + } + + // Load state from udata (as with res function) + load_vars(udata); + + // Load vector to be inverted into F_vars + load_derivs(rvec); + + runPreconditioner(t, gamma, delta); + + // Save the solution from F_vars + save_derivs(zvec); + + pre_Wtime_s += bout::globals::mpi->MPI_Wtime() - tstart; + pre_ncalls_s++; +} + +void ArkodeMRISolver::pre_f(BoutReal t, BoutReal gamma, BoutReal delta, BoutReal* udata, + BoutReal* rvec, BoutReal* zvec) { + TRACE("Running preconditioner: ArkodeMRISolver::pre({:e})", t); + + const BoutReal tstart = bout::globals::mpi->MPI_Wtime(); + + if (!hasPreconditioner()) { + // Identity (but should never happen) + const auto length = N_VGetLocalLength_Parallel(uvec); + std::copy(rvec, rvec + length, zvec); + return; + } + + // Load state from udata (as with res function) + load_vars(udata); + + // Load vector to be inverted into F_vars + load_derivs(rvec); + + runPreconditioner(t, gamma, delta); + + // Save the solution from F_vars + save_derivs(zvec); + + pre_Wtime_s += bout::globals::mpi->MPI_Wtime() - tstart; + pre_ncalls_s++; +} + +/************************************************************************** + * Jacobian-vector multiplication functions + **************************************************************************/ + +void ArkodeMRISolver::jac_s(BoutReal t, BoutReal* ydata, BoutReal* vdata, BoutReal* Jvdata) { + TRACE("Running Jacobian: ArkodeMRISolver::jac({:e})", t); + + if (not hasJacobian()) { + throw BoutException("No jacobian function supplied!\n"); + } + + // Load state from ydate + load_vars(ydata); + + // Load vector to be multiplied into F_vars + load_derivs(vdata); + + // Call function + runJacobian(t); + + // Save Jv from vars + save_derivs(Jvdata); +} + +void ArkodeMRISolver::jac_f(BoutReal t, BoutReal* ydata, BoutReal* vdata, BoutReal* Jvdata) { + TRACE("Running Jacobian: ArkodeMRISolver::jac({:e})", t); + + if (not hasJacobian()) { + throw BoutException("No jacobian function supplied!\n"); + } + + // Load state from ydate + load_vars(ydata); + + // Load vector to be multiplied into F_vars + load_derivs(vdata); + + // Call function + runJacobian(t); + + // Save Jv from vars + save_derivs(Jvdata); +} + +/************************************************************************** + * ARKODE explicit RHS functions + **************************************************************************/ + +// NOLINTBEGIN(readability-identifier-length) +namespace { +int arkode_rhs_s_explicit(BoutReal t, N_Vector u, N_Vector du, void* user_data) { + + BoutReal* udata = N_VGetArrayPointer(u); + BoutReal* dudata = N_VGetArrayPointer(du); + + auto* s = static_cast(user_data); + + // Calculate RHS function + try { + s->rhs_se(t, udata, dudata); + } catch (BoutRhsFail& error) { + return 1; + } + return 0; +} + +int arkode_rhs_s_implicit(BoutReal t, N_Vector u, N_Vector du, void* user_data) { + + BoutReal* udata = N_VGetArrayPointer(u); + BoutReal* dudata = N_VGetArrayPointer(du); + + auto* s = static_cast(user_data); + + // Calculate RHS function + try { + s->rhs_si(t, udata, dudata); + } catch (BoutRhsFail& error) { + return 1; + } + return 0; +} + +int arkode_rhs_f_explicit(BoutReal t, N_Vector u, N_Vector du, void* user_data) { + + BoutReal* udata = N_VGetArrayPointer(u); + BoutReal* dudata = N_VGetArrayPointer(du); + + auto* s = static_cast(user_data); + + // Calculate RHS function + try { + s->rhs_fe(t, udata, dudata); + } catch (BoutRhsFail& error) { + return 1; + } + return 0; +} + +int arkode_rhs_f_implicit(BoutReal t, N_Vector u, N_Vector du, void* user_data) { + + BoutReal* udata = N_VGetArrayPointer(u); + BoutReal* dudata = N_VGetArrayPointer(du); + + auto* s = static_cast(user_data); + + // Calculate RHS function + try { + s->rhs_fi(t, udata, dudata); + } catch (BoutRhsFail& error) { + return 1; + } + return 0; +} + +int arkode_s_rhs(BoutReal t, N_Vector u, N_Vector du, void* user_data) { + + BoutReal* udata = N_VGetArrayPointer(u); + BoutReal* dudata = N_VGetArrayPointer(du); + + auto* s = static_cast(user_data); + + // Calculate RHS function + try { + s->rhs_s(t, udata, dudata); + } catch (BoutRhsFail& error) { + return 1; + } + return 0; +} + +int arkode_f_rhs(BoutReal t, N_Vector u, N_Vector du, void* user_data) { + + BoutReal* udata = N_VGetArrayPointer(u); + BoutReal* dudata = N_VGetArrayPointer(du); + + auto* s = static_cast(user_data); + + // Calculate RHS function + try { + s->rhs_f(t, udata, dudata); + } catch (BoutRhsFail& error) { + return 1; + } + return 0; +} + +/// RHS function for BBD preconditioner +int arkode_s_bbd_rhs(sunindextype UNUSED(Nlocal), BoutReal t, N_Vector u, N_Vector du, + void* user_data) { + return arkode_rhs_s_implicit(t, u, du, user_data); +} + +int arkode_f_bbd_rhs(sunindextype UNUSED(Nlocal), BoutReal t, N_Vector u, N_Vector du, + void* user_data) { + return arkode_rhs_f_implicit(t, u, du, user_data); +} + +/// Preconditioner function +int arkode_s_pre(BoutReal t, N_Vector yy, N_Vector UNUSED(yp), N_Vector rvec, N_Vector zvec, + BoutReal gamma, BoutReal delta, int UNUSED(lr), void* user_data) { + BoutReal* udata = N_VGetArrayPointer(yy); + BoutReal* rdata = N_VGetArrayPointer(rvec); + BoutReal* zdata = N_VGetArrayPointer(zvec); + + auto* s = static_cast(user_data); + + // Calculate residuals + s->pre_s(t, gamma, delta, udata, rdata, zdata); + + return 0; +} + +int arkode_f_pre(BoutReal t, N_Vector yy, N_Vector UNUSED(yp), N_Vector rvec, N_Vector zvec, + BoutReal gamma, BoutReal delta, int UNUSED(lr), void* user_data) { + BoutReal* udata = N_VGetArrayPointer(yy); + BoutReal* rdata = N_VGetArrayPointer(rvec); + BoutReal* zdata = N_VGetArrayPointer(zvec); + + auto* s = static_cast(user_data); + + // Calculate residuals + s->pre_f(t, gamma, delta, udata, rdata, zdata); + + return 0; +} + +} // namespace +// NOLINTEND(readability-identifier-length) + +/************************************************************************** + * vector abstol functions + **************************************************************************/ + +void ArkodeMRISolver::set_abstol_values(BoutReal* abstolvec_data, + std::vector& f2dtols, + std::vector& f3dtols) { + int p = 0; // Counter for location in abstolvec_data array + + // All boundaries + for (const auto& i2d : bout::globals::mesh->getRegion2D("RGN_BNDRY")) { + loop_abstol_values_op(i2d, abstolvec_data, p, f2dtols, f3dtols, true); + } + // Bulk of points + for (const auto& i2d : bout::globals::mesh->getRegion2D("RGN_NOBNDRY")) { + loop_abstol_values_op(i2d, abstolvec_data, p, f2dtols, f3dtols, false); + } +} + +void ArkodeMRISolver::loop_abstol_values_op(Ind2D UNUSED(i2d), BoutReal* abstolvec_data, + int& p, std::vector& f2dtols, + std::vector& f3dtols, bool bndry) { + // Loop over 2D variables + for (std::vector::size_type i = 0; i < f2dtols.size(); i++) { + if (bndry && !f2d[i].evolve_bndry) { + continue; + } + abstolvec_data[p] = f2dtols[i]; + p++; + } + + for (int jz = 0; jz < bout::globals::mesh->LocalNz; jz++) { + // Loop over 3D variables + for (std::vector::size_type i = 0; i < f3dtols.size(); i++) { + if (bndry && !f3d[i].evolve_bndry) { + continue; + } + abstolvec_data[p] = f3dtols[i]; + p++; + } + } +} + +#endif diff --git a/src/solver/impls/arkode/arkode_mri.hxx b/src/solver/impls/arkode/arkode_mri.hxx new file mode 100644 index 0000000000..80c9c51ac8 --- /dev/null +++ b/src/solver/impls/arkode/arkode_mri.hxx @@ -0,0 +1,176 @@ +/************************************************************************** + * Interface to ARKODE MRI solver + * NOTE: ARKode is currently in beta testing so use with cautious optimism + * + * NOTE: Only one solver can currently be compiled in + * + ************************************************************************** + * Copyright 2010-2024 BOUT++ contributors + * + * Contact: Ben Dudson, dudson2@llnl.gov + * + * This file is part of BOUT++. + * + * BOUT++ is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * BOUT++ is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with BOUT++. If not, see . + * + **************************************************************************/ + +#ifndef BOUT_ARKODE_MRI_SOLVER_H +#define BOUT_ARKODE_MRI_SOLVER_H + +#include "bout/build_config.hxx" +#include "bout/solver.hxx" + +#if not BOUT_HAS_ARKODE + +namespace { +RegisterUnavailableSolver + registerunavailablearkodemri("arkode_mri", "BOUT++ was not configured with ARKODE/SUNDIALS"); +} + +#else + +#include "bout/bout_enum_class.hxx" +#include "bout/bout_types.hxx" +#include "bout/sundials_backports.hxx" + +#include +#include +#include + +#include + +#include + +class ArkodeMRISolver; +class Options; + +namespace { +RegisterSolver registersolverarkodemri("arkode_mri"); +} + +// enum describing treatment of equations +// Note: Capitalized because `explicit` is a C++ reserved keyword +BOUT_ENUM_CLASS(MRI_Treatment, ImEx, Implicit, Explicit); + +// Adaptivity method +BOUT_ENUM_CLASS(MRI_AdapMethod, PID, PI, I, Explicit_Gustafsson, Implicit_Gustafsson, + ImEx_Gustafsson); + +class ArkodeMRISolver : public Solver { +public: + explicit ArkodeMRISolver(Options* opts = nullptr); + ~ArkodeMRISolver(); + + BoutReal getCurrentTimestep() override { return hcur; } + + int init() override; + + int run() override; + BoutReal run(BoutReal tout); + + // These functions used internally (but need to be public) + void rhs_se(BoutReal t, BoutReal* udata, BoutReal* dudata); + void rhs_si(BoutReal t, BoutReal* udata, BoutReal* dudata); + void rhs_fe(BoutReal t, BoutReal* udata, BoutReal* dudata); + void rhs_fi(BoutReal t, BoutReal* udata, BoutReal* dudata); + void rhs_s(BoutReal t, BoutReal* udata, BoutReal* dudata); + void rhs_f(BoutReal t, BoutReal* udata, BoutReal* dudata); + void pre_s(BoutReal t, BoutReal gamma, BoutReal delta, BoutReal* udata, BoutReal* rvec, + BoutReal* zvec); + void pre_f(BoutReal t, BoutReal gamma, BoutReal delta, BoutReal* udata, BoutReal* rvec, + BoutReal* zvec); + void jac_s(BoutReal t, BoutReal* ydata, BoutReal* vdata, BoutReal* Jvdata); + void jac_f(BoutReal t, BoutReal* ydata, BoutReal* vdata, BoutReal* Jvdata); + +private: + BoutReal hcur; //< Current internal timestep + + bool diagnose{false}; //< Output additional diagnostics + + N_Vector uvec{nullptr}; //< Values + void* arkode_mem{nullptr}; //< ARKODE internal memory block + void* inner_arkode_mem{nullptr}; //< ARKODE internal memory block + MRIStepInnerStepper inner_stepper{nullptr}; //< inner stepper + + BoutReal pre_Wtime_s{0.0}; //< Time in preconditioner + BoutReal pre_Wtime_f{0.0}; //< Time in preconditioner + int pre_ncalls_s{0}; //< Number of calls to preconditioner + int pre_ncalls_f{0}; //< Number of calls to preconditioner + + /// Maximum number of steps to take between outputs + int mxsteps; + /// Integrator treatment enum: IMEX, Implicit or Explicit + MRI_Treatment treatment; + MRI_Treatment inner_treatment; + /// Use linear implicit solver (only evaluates jacobian inversion once) + bool set_linear; + bool inner_set_linear; + /// Solve explicit portion in fixed timestep mode. NOTE: This is not recommended except + /// for code comparison + bool fixed_step; + /// Order of the internal step + int order; + /// Timestep adaptivity function + MRI_AdapMethod adap_method; + /// Absolute tolerance + BoutReal abstol; + /// Relative tolerance + BoutReal reltol; + /// Use separate absolute tolerance for each field + bool use_vector_abstol; + /// Maximum timestep (only used if greater than zero) + bool use_precon; + bool inner_use_precon; + /// Number of Krylov basis vectors to use + int maxl; + int inner_maxl; + /// Use right preconditioning instead of left preconditioning + bool rightprec; + + // Diagnostics from ARKODE MRI + int nsteps{0}; + int nfe_evals{0}; + int nfi_evals{0}; + int nniters{0}; + int npevals{0}; + int nliters{0}; + int inner_nsteps{0}; + int inner_nfe_evals{0}; + int inner_nfi_evals{0}; + int inner_nniters{0}; + int inner_npevals{0}; + int inner_nliters{0}; + + void set_abstol_values(BoutReal* abstolvec_data, std::vector& f2dtols, + std::vector& f3dtols); + void loop_abstol_values_op(Ind2D i2d, BoutReal* abstolvec_data, int& p, + std::vector& f2dtols, + std::vector& f3dtols, bool bndry); + + /// SPGMR solver structure + SUNLinearSolver sun_solver{nullptr}; + SUNLinearSolver inner_sun_solver{nullptr}; + /// Solver for implicit stages + SUNNonlinearSolver nonlinear_solver{nullptr}; + SUNNonlinearSolver inner_nonlinear_solver{nullptr}; + /// Timestep controller + SUNAdaptController controller{nullptr}; + SUNAdaptController inner_controller{nullptr}; + /// Context for SUNDIALS memory allocations + sundials::Context suncontext; +}; + +#endif // BOUT_HAS_ARKODE +#endif // BOUT_ARKODE_MRI_SOLVER_H diff --git a/tests/integrated/test-kpr_mri/CMakeLists.txt b/tests/integrated/test-kpr_mri/CMakeLists.txt new file mode 100644 index 0000000000..66a4474bb0 --- /dev/null +++ b/tests/integrated/test-kpr_mri/CMakeLists.txt @@ -0,0 +1,2 @@ +bout_add_integrated_test(test_kpr_mri SOURCES test_kpr_mri.cxx +REQUIRES BOUT_HAS_SUNDIALS) diff --git a/tests/integrated/test-kpr_mri/README.md b/tests/integrated/test-kpr_mri/README.md new file mode 100644 index 0000000000..ffcef3ded9 --- /dev/null +++ b/tests/integrated/test-kpr_mri/README.md @@ -0,0 +1,37 @@ +test-kpr_mri +=========== + + Multirate nonlinear Kvaerno-Prothero-Robinson ODE test problem: + + [f]' = [ G e ] [(-1+f^2-r)/(2f)] + [ r'(t)/(2f) ] + [g] [ e -1 ] [(-2+g^2-s)/(2g)] [ s'(t)/(2*sqrt(2+s(t))) ] + = [ fs(t,f,g) ] + [ ff(t,f,g) ] + + where r(t) = 0.5 * cos(t), s(t) = cos(w * t), 0 < t < 5. + + This problem has analytical solution given by + f(t) = sqrt(1+r(t)), g(t) = sqrt(2+s(t)). + + We use the parameters: + e = 0.5 (fast/slow coupling strength) [default] + G = -100 (stiffness at slow time scale) [default] + w = 100 (time-scale separation factor) [default] + + The stiffness of the slow time scale is essentially determined + by G, for |G| > 50 it is 'stiff' and ideally suited to a + multirate method that is implicit at the slow time scale. + +MRI implementations of the functions are as follows: + +The slow explicit RHS function: + [-0.5 * sin(t)/(2 * f)] + [ 0 ] + +The slow implicit RHS function: + [G e] * [(-1 + f^2 - 0.5 * cos(t))/(2 * f) ] + [0 0] [(-2 + g^2 - cos(w * t))/(2 * g) ] + +The fast implicit RHS function: + [0 0] * [(-1 + f^2 - 0.5 * cos(t))/(2 * f)] + [ 0 ] + [e -1] [(-2 + g^2 - cos(w * t))/(2 * g) ] [-w * sin(w * t)/(2 * sqrt(2 + cos(w*t)))] diff --git a/tests/integrated/test-kpr_mri/makefile b/tests/integrated/test-kpr_mri/makefile new file mode 100644 index 0000000000..84db76e5c3 --- /dev/null +++ b/tests/integrated/test-kpr_mri/makefile @@ -0,0 +1,5 @@ +BOUT_TOP = ../../.. + +SOURCEC = test_kpr_mri.cxx + +include $(BOUT_TOP)/make.config diff --git a/tests/integrated/test-kpr_mri/runtest b/tests/integrated/test-kpr_mri/runtest new file mode 100755 index 0000000000..47f335dc94 --- /dev/null +++ b/tests/integrated/test-kpr_mri/runtest @@ -0,0 +1,23 @@ +#!/usr/bin/env python3 + +# requires: petsc + +from boututils.run_wrapper import shell_safe, launch_safe + +from sys import exit + +nthreads = 1 +nproc = 1 + +print("Making solver test") +shell_safe("make > make.log") + +print("Running solver test") +status, out = launch_safe("./test_kpr_mri", nproc=nproc, mthread=nthreads, pipe=True) +with open("run.log", "w") as f: + f.write(out) + +if status: + print(out) + +exit(status) diff --git a/tests/integrated/test-kpr_mri/test_kpr_mri.cxx b/tests/integrated/test-kpr_mri/test_kpr_mri.cxx new file mode 100644 index 0000000000..c071f25ef2 --- /dev/null +++ b/tests/integrated/test-kpr_mri/test_kpr_mri.cxx @@ -0,0 +1,183 @@ +#include "bout/physicsmodel.hxx" +#include "bout/solver.hxx" + +#include +#include +#include +#include + +// A simple phyics model with a manufactured true solution +// +class TestSolver : public PhysicsModel { +public: + Field3D f, g; + + BoutReal e = 0.5; /* fast/slow coupling strength */ + BoutReal G = -100.0; /* stiffness at slow time scale */ + BoutReal w = 100.0; /* time-scale separation factor */ + + int init(bool UNUSED(restarting)) override { + solver->add(f, "f"); + solver->add(g, "g"); + + f = sqrt(3.0/2.0); + g = sqrt(3.0); + + setSplitOperatorMRI(); + + return 0; + } + + int rhs_se(BoutReal t) override { + /* fill in the slow explicit RHS function: + [-0.5*sin(t)/(2*f)] + [ 0 ] */ + ddt(f) = -0.5*sin(t)/(2.0*f(1,1,0)); + ddt(g) = 0.0; + + return 0; + } + + int rhs_si(BoutReal t) override { + /* fill in the slow implicit RHS function: + [G e]*[(-1+f^2-0.5*cos(t))/(2*f)] + [0 0] [(-2+g^2-cos(w*t))/(2*g) ] */ + BoutReal tmp1 = (-1.0 + f(1,1,0) * f(1,1,0) - 0.5*cos(t)) / (2.0 * f(1,1,0)); + BoutReal tmp2 = (-2.0 + g(1,1,0) * g(1,1,0) - cos(w*t)) / (2.0 * g(1,1,0)); + ddt(f) = G * tmp1 + e * tmp2; + ddt(g) = 0.0; + + return 0; + } + + int rhs_fe(BoutReal UNUSED(t)) override { + + ddt(f) = 0.0; + ddt(g) = 0.0; + + return 0; + } + + int rhs_fi(BoutReal t) override { + /* fill in the fast implicit RHS function: + [0 0]*[(-1+f^2-0.5*cos(t))/(2*f)] + [ 0 ] + [e -1] [(-2+g^2-cos(w*t))/(2*g) ] [-w*sin(w*t)/(2*sqrt(2+cos(w*t)))] */ + BoutReal tmp1 = (-1.0 + f(1,1,0) * f(1,1,0) - 0.5*cos(t)) / (2.0 * f(1,1,0)); + BoutReal tmp2 = (-2.0 + g(1,1,0) * g(1,1,0) - cos(w*t)) / (2.0 * g(1,1,0)); + ddt(f) = 0.0; + ddt(g) = e * tmp1 - tmp2 - w * sin(w*t) / (2.0 * sqrt(2.0 + cos(w * t))); + + return 0; + } + + int rhs_s(BoutReal t) override { + /* fill in the RHS function: + [G e]*[(-1+f^2-0.5*cos(t))/(2*f)] + [-0.5*sin(t)/(2*f)] + [0 0] [(-2+g^2-cos(w*t))/(2*g) ] [ 0 ] */ + BoutReal tmp1 = (-1.0 + f(1,1,0) * f(1,1,0) - 0.5*cos(t)) / (2.0 * f(1,1,0)); + BoutReal tmp2 = (-2.0 + g(1,1,0) * g(1,1,0) - cos(w*t)) / (2.0 * g(1,1,0)); + ddt(f) = G * tmp1 + e * tmp2 - 0.5*sin(t) / (2.0 * f(1,1,0)); + ddt(g) = 0.0; + + return 0; + } + + int rhs_f(BoutReal t) override { + /* fill in the RHS function: + [0 0]*[(-1+f^2-0.5*cos(t))/(2*f)] + [ 0 ] + [e -1] [(-2+g^2-cos(w*t))/(2*g) ] [-w*sin(w*t)/(2*sqrt(2+cos(w*t)))] */ + BoutReal tmp1 = (-1.0 + f(1,1,0) * f(1,1,0) - 0.5*cos(t)) / (2.0 * f(1,1,0)); + BoutReal tmp2 = (-2.0 + g(1,1,0) * g(1,1,0) - cos(w*t)) / (2.0 * g(1,1,0)); + ddt(f) = 0.0; + ddt(g) = e * tmp1 - tmp2 - w * sin(w*t) / (2.0 * sqrt(2.0 + cos(w * t))); + + return 0; + } + + // int rhs(BoutReal t) override { + // /* fill in the RHS function: + // [G e]*[(-1+f^2-0.5*cos(t))/(2*f)] + [-0.5*sin(t)/(2*f) ] + // [e -1] [(-2+g^2-cos(w*t))/(2*g) ] [-w*sin(w*t)/(2*sqrt(2+cos(w*t)))] */ + // BoutReal tmp1 = (-1.0 + f(1,1,0) * f(1,1,0) - 0.5*cos(t)) / (2.0 * f(1,1,0)); + // BoutReal tmp2 = (-2.0 + g(1,1,0) * g(1,1,0) - cos(w*t)) / (2.0 * g(1,1,0)); + + // ddt(f) = G * tmp1 + e * tmp2 - 0.5*sin(t) / (2.0 * f(1,1,0)); + // ddt(g) = e * tmp1 - tmp2 - w * sin(w*t) / (2.0 * sqrt(2.0 + cos(w * t))); + + // return 0; + // } + + bool check_solution(BoutReal atol, BoutReal t) { + // Return true if correct solution + return ((std::abs(sqrt(0.5*cos(t) + 1.0) - f(1,1,0)) < atol) and (std::abs(sqrt(cos(w*t) + 2.0) - g(1,1,0)) < atol)); + } + + BoutReal compute_error(BoutReal t) + { + /* Compute the error with the true solution: + f(t) = sqrt(0.5*cos(t) + 1.0) + g(t) = sqrt(cos(w*t) + 2.0) */ + return sqrt( pow(sqrt(0.5*cos(t) + 1.0) - f(1,1,0), 2.0) + + pow(sqrt(cos(w*t) + 2.0) - g(1,1,0), 2.0)); + } + + // Don't need any restarting, or options to control data paths + int postInit(bool) override { return 0; } +}; + +int main(int argc, char** argv) { + // Absolute tolerance for difference between the actual value and the + // expected value + constexpr BoutReal tolerance = 1.e-5; + + // Our own output to stdout, as main library will only be writing to log files + Output output_test; + + auto& root = Options::root(); + + root["mesh"]["MXG"] = 1; + root["mesh"]["MYG"] = 1; + root["mesh"]["nx"] = 3; + root["mesh"]["ny"] = 1; + root["mesh"]["nz"] = 1; + + root["output"]["enabled"] = false; + root["restart_files"]["enabled"] = false; + + Solver::setArgs(argc, argv); + BoutComm::setArgs(argc, argv); + + bout::globals::mpi = new MpiWrapper(); + + bout::globals::mesh = Mesh::create(); + bout::globals::mesh->load(); + + // Global options + root["nout"] = 100; + root["timestep"] = 0.05; + + // Get specific options section for this solver. Can't just use default + // "solver" section, as we run into problems when solvers use the same + // name for an option with inconsistent defaults + auto options = Options::getRoot()->getSection("arkode_mri"); + auto solver = std::unique_ptr{Solver::create("arkode_mri", options)}; + + TestSolver model{}; + solver->setModel(&model); + + BoutMonitor bout_monitor{}; + solver->addMonitor(&bout_monitor, Solver::BACK); + + solver->solve(); + + BoutReal error = model.compute_error(5.0); + + std::cout << "error = " << error << std::endl; + + if (model.check_solution(tolerance, 5.0)) { + output_test << " PASSED\n"; + return 0; + } + output_test << " FAILED\n"; + return 1; +} From fbba15f7a7807d6ffe6f3f86e69ed65ccba346f6 Mon Sep 17 00:00:00 2001 From: maggul Date: Fri, 25 Oct 2024 11:37:25 -0700 Subject: [PATCH 03/28] added treatment options --- tests/integrated/test-kpr_mri/test_kpr_mri.cxx | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tests/integrated/test-kpr_mri/test_kpr_mri.cxx b/tests/integrated/test-kpr_mri/test_kpr_mri.cxx index c071f25ef2..250f874fe6 100644 --- a/tests/integrated/test-kpr_mri/test_kpr_mri.cxx +++ b/tests/integrated/test-kpr_mri/test_kpr_mri.cxx @@ -24,7 +24,7 @@ class TestSolver : public PhysicsModel { g = sqrt(3.0); setSplitOperatorMRI(); - + return 0; } @@ -155,6 +155,8 @@ int main(int argc, char** argv) { // Global options root["nout"] = 100; root["timestep"] = 0.05; + root["arkode_mri"]["treatment"] = "explicit"; + root["arkode_mri"]["inner_treatment"] = "explicit"; // Get specific options section for this solver. Can't just use default // "solver" section, as we run into problems when solvers use the same From 0861085f912d0b2034d18f0d58f4cac4d962bda5 Mon Sep 17 00:00:00 2001 From: maggul Date: Mon, 18 Nov 2024 09:33:09 -0600 Subject: [PATCH 04/28] remove adaptivity --- src/solver/impls/arkode/arkode_mri.cxx | 55 +------------------------- src/solver/impls/arkode/arkode_mri.hxx | 11 ------ 2 files changed, 1 insertion(+), 65 deletions(-) diff --git a/src/solver/impls/arkode/arkode_mri.cxx b/src/solver/impls/arkode/arkode_mri.cxx index 1fb9a79a24..e092f7be28 100644 --- a/src/solver/impls/arkode/arkode_mri.cxx +++ b/src/solver/impls/arkode/arkode_mri.cxx @@ -99,11 +99,6 @@ ArkodeMRISolver::ArkodeMRISolver(Options* opts) "not recommended except for code comparison") .withDefault(false)), order((*options)["order"].doc("Order of internal step").withDefault(4)), - adap_method( - (*options)["adap_method"] - .doc("Set timestep adaptivity function: pid, pi, i, explicit_gustafsson, " - "implicit_gustafsson, imex_gustafsson.") - .withDefault(MRI_AdapMethod::PID)), abstol((*options)["atol"].doc("Absolute tolerance").withDefault(1.0e-12)), reltol((*options)["rtol"].doc("Relative tolerance").withDefault(1.0e-5)), use_vector_abstol((*options)["use_vector_abstol"] @@ -157,9 +152,6 @@ ArkodeMRISolver::~ArkodeMRISolver() { SUNNonlinSolFree(nonlinear_solver); SUNNonlinSolFree(inner_nonlinear_solver); MRIStepInnerStepper_Free(&inner_stepper); - - SUNAdaptController_Destroy(controller); - SUNAdaptController_Destroy(inner_controller); } /************************************************************************** @@ -293,53 +285,8 @@ int ArkodeMRISolver::init() { throw BoutException("ARKodeSetOrder failed\n"); } - switch (adap_method) { - case MRI_AdapMethod::PID: - controller = SUNAdaptController_PID(suncontext); - inner_controller = SUNAdaptController_PID(suncontext); - break; - case MRI_AdapMethod::PI: - controller = SUNAdaptController_PI(suncontext); - inner_controller = SUNAdaptController_PI(suncontext); - break; - case MRI_AdapMethod::I: - controller = SUNAdaptController_I(suncontext); - inner_controller = SUNAdaptController_I(suncontext); - break; - case MRI_AdapMethod::Explicit_Gustafsson: - controller = SUNAdaptController_ExpGus(suncontext); - inner_controller = SUNAdaptController_ExpGus(suncontext); - break; - case MRI_AdapMethod::Implicit_Gustafsson: - controller = SUNAdaptController_ImpGus(suncontext); - inner_controller = SUNAdaptController_ImpGus(suncontext); - break; - case MRI_AdapMethod::ImEx_Gustafsson: - controller = SUNAdaptController_ImExGus(suncontext); - inner_controller = SUNAdaptController_ImExGus(suncontext); - break; - default: - throw BoutException("Invalid adap_method\n"); - } - - // if (ARKodeSetAdaptController(arkode_mem, controller) != ARK_SUCCESS) { - // throw BoutException("ARKodeSetAdaptController failed\n"); - // } - - // if (ARKodeSetAdaptivityAdjustment(arkode_mem, 0) != ARK_SUCCESS) { - // throw BoutException("ARKodeSetAdaptivityAdjustment failed\n"); - // } - if (ARKodeSetFixedStep(arkode_mem, 0.001) != ARK_SUCCESS) { - throw BoutException("ARKodeSetAdaptController failed\n"); - } - - if (ARKodeSetAdaptController(inner_arkode_mem, controller) != ARK_SUCCESS) { - throw BoutException("ARKodeSetAdaptController failed\n"); - } - - if (ARKodeSetAdaptivityAdjustment(inner_arkode_mem, 0) != ARK_SUCCESS) { - throw BoutException("ARKodeSetAdaptivityAdjustment failed\n"); + throw BoutException("ARKodeSetFixedStep failed\n"); } if (use_vector_abstol) { diff --git a/src/solver/impls/arkode/arkode_mri.hxx b/src/solver/impls/arkode/arkode_mri.hxx index 80c9c51ac8..9cd33c2df8 100644 --- a/src/solver/impls/arkode/arkode_mri.hxx +++ b/src/solver/impls/arkode/arkode_mri.hxx @@ -49,8 +49,6 @@ RegisterUnavailableSolver #include #include -#include - #include class ArkodeMRISolver; @@ -64,10 +62,6 @@ RegisterSolver registersolverarkodemri("arkode_mri"); // Note: Capitalized because `explicit` is a C++ reserved keyword BOUT_ENUM_CLASS(MRI_Treatment, ImEx, Implicit, Explicit); -// Adaptivity method -BOUT_ENUM_CLASS(MRI_AdapMethod, PID, PI, I, Explicit_Gustafsson, Implicit_Gustafsson, - ImEx_Gustafsson); - class ArkodeMRISolver : public Solver { public: explicit ArkodeMRISolver(Options* opts = nullptr); @@ -122,8 +116,6 @@ private: bool fixed_step; /// Order of the internal step int order; - /// Timestep adaptivity function - MRI_AdapMethod adap_method; /// Absolute tolerance BoutReal abstol; /// Relative tolerance @@ -165,9 +157,6 @@ private: /// Solver for implicit stages SUNNonlinearSolver nonlinear_solver{nullptr}; SUNNonlinearSolver inner_nonlinear_solver{nullptr}; - /// Timestep controller - SUNAdaptController controller{nullptr}; - SUNAdaptController inner_controller{nullptr}; /// Context for SUNDIALS memory allocations sundials::Context suncontext; }; From 9795739fb170b93c955af14f80fe60c6ed0610bd Mon Sep 17 00:00:00 2001 From: maggul Date: Mon, 18 Nov 2024 10:14:20 -0600 Subject: [PATCH 05/28] bulk revision --- src/solver/impls/arkode/arkode_mri.cxx | 91 +++++-------------- src/solver/impls/arkode/arkode_mri.hxx | 2 +- tests/integrated/test-kpr_mri/README.md | 6 +- .../integrated/test-kpr_mri/test_kpr_mri.cxx | 41 +++------ 4 files changed, 40 insertions(+), 100 deletions(-) diff --git a/src/solver/impls/arkode/arkode_mri.cxx b/src/solver/impls/arkode/arkode_mri.cxx index e092f7be28..d68748411f 100644 --- a/src/solver/impls/arkode/arkode_mri.cxx +++ b/src/solver/impls/arkode/arkode_mri.cxx @@ -1,7 +1,7 @@ /************************************************************************** - * Experimental interface to SUNDIALS ARKode MRI solver + * Experimental interface to SUNDIALS ARKODE MRI solver * - * NOTE: ARKode is still in beta testing so use with cautious optimism + * NOTE: ARKODE is still in beta testing so use with cautious optimism * ************************************************************************** * Copyright 2010-2024 BOUT++ contributors @@ -62,7 +62,7 @@ int arkode_f_rhs(BoutReal t, N_Vector u, N_Vector du, void* user_data); int arkode_s_bbd_rhs(sunindextype Nlocal, BoutReal t, N_Vector u, N_Vector du, void* user_data); int arkode_f_bbd_rhs(sunindextype Nlocal, BoutReal t, N_Vector u, N_Vector du, - void* user_data); + void* user_data); int arkode_s_pre(BoutReal t, N_Vector yy, N_Vector yp, N_Vector rvec, N_Vector zvec, BoutReal gamma, BoutReal delta, int lr, void* user_data); int arkode_f_pre(BoutReal t, N_Vector yy, N_Vector yp, N_Vector rvec, N_Vector zvec, @@ -191,17 +191,17 @@ int ArkodeMRISolver::init() { case MRI_Treatment::ImEx: inner_arkode_mem = callWithSUNContext(ARKStepCreate, suncontext, arkode_rhs_f_explicit, arkode_rhs_f_implicit, simtime, uvec); - output_info.write("\tUsing ARKode ImEx inner solver \n"); + output_info.write("\tUsing ARKODE ImEx inner solver \n"); break; case MRI_Treatment::Explicit: inner_arkode_mem = callWithSUNContext(ARKStepCreate, suncontext, arkode_f_rhs, nullptr, simtime, uvec); - output_info.write("\tUsing ARKode Explicit inner solver \n"); + output_info.write("\tUsing ARKODE Explicit inner solver \n"); break; case MRI_Treatment::Implicit: inner_arkode_mem = callWithSUNContext(ARKStepCreate, suncontext, nullptr, arkode_f_rhs, simtime, uvec); - output_info.write("\tUsing ARKode Implicit inner solver \n"); + output_info.write("\tUsing ARKODE Implicit inner solver \n"); break; default: throw BoutException("Invalid inner_treatment: {}\n", toString(inner_treatment)); @@ -242,19 +242,19 @@ int ArkodeMRISolver::init() { switch (treatment) { case MRI_Treatment::ImEx: - arkode_mem = callWithSUNContext(MRIStepCreate, suncontext, arkode_rhs_s_explicit, arkode_rhs_s_implicit, + arkode_mem = callWithSUNContext(MRIStepCreate, suncontext, arkode_rhs_s_explicit, arkode_rhs_s_implicit, simtime, uvec, inner_stepper); - output_info.write("\tUsing ARKode ImEx solver \n"); + output_info.write("\tUsing ARKODE-MRI ImEx solver \n"); break; case MRI_Treatment::Explicit: - arkode_mem = callWithSUNContext(MRIStepCreate, suncontext, arkode_s_rhs, nullptr, + arkode_mem = callWithSUNContext(MRIStepCreate, suncontext, arkode_s_rhs, nullptr, simtime, uvec, inner_stepper); - output_info.write("\tUsing ARKode Explicit solver \n"); + output_info.write("\tUsing ARKODE-MRI Explicit solver \n"); break; case MRI_Treatment::Implicit: arkode_mem = callWithSUNContext(MRIStepCreate, suncontext, nullptr, arkode_s_rhs, simtime, uvec, inner_stepper); - output_info.write("\tUsing ARKode Implicit solver \n"); + output_info.write("\tUsing ARKODE-MRI Implicit solver \n"); break; default: throw BoutException("Invalid treatment: {}\n", toString(treatment)); @@ -381,24 +381,24 @@ int ArkodeMRISolver::init() { return acc + localmesh->xend - localmesh->xstart + 3; }); - const auto mudq = (*options)["mudq"] + const auto inner_mudq = (*options)["inner_mudq"] .doc("Upper half-bandwidth to be used in the difference " "quotient Jacobian approximation") .withDefault(band_width_default); - const auto mldq = (*options)["mldq"] + const auto inner_mldq = (*options)["inner_mldq"] .doc("Lower half-bandwidth to be used in the difference " "quotient Jacobian approximation") .withDefault(band_width_default); - const auto mukeep = (*options)["mukeep"] + const auto inner_mukeep = (*options)["inner_mukeep"] .doc("Upper half-bandwidth of the retained banded " "approximate Jacobian block") .withDefault(n3Dvars() + n2Dvars()); - const auto mlkeep = (*options)["mlkeep"] + const auto inner_mlkeep = (*options)["mlkeep"] .doc("Lower half-bandwidth of the retained banded " "approximate Jacobian block") .withDefault(n3Dvars() + n2Dvars()); - if (ARKBBDPrecInit(inner_arkode_mem, local_N, mudq, mldq, mukeep, mlkeep, 0, + if (ARKBBDPrecInit(inner_arkode_mem, local_N, inner_mudq, inner_mldq, inner_mukeep, inner_mlkeep, 0, arkode_f_bbd_rhs, nullptr) != ARKLS_SUCCESS) { throw BoutException("ARKBBDPrecInit failed for inner solver\n"); @@ -511,7 +511,7 @@ int ArkodeMRISolver::run() { // Step failed output.write("Timestep failed. Aborting\n"); - throw BoutException("ARKode timestep failed\n"); + throw BoutException("ARKODE timestep failed\n"); } // Get additional diagnostics @@ -609,8 +609,7 @@ BoutReal ArkodeMRISolver::run(BoutReal tout) { // Call timestep monitor call_timestep_monitors(internal_time, internal_time - last_time); } - // Get output at the desired time - flag = ARKodeGetDky(arkode_mem, tout, 0, uvec); + // Update the current simulation time simtime = tout; } @@ -713,7 +712,6 @@ void ArkodeMRISolver::rhs_s(BoutReal t, BoutReal* udata, BoutReal* dudata) { ARKodeGetLastStep(arkode_mem, &hcur); // Call RHS function - // run_rhs_s(t); run_rhs_s(t); // Save derivatives to dudata @@ -735,7 +733,6 @@ void ArkodeMRISolver::rhs_f(BoutReal t, BoutReal* udata, BoutReal* dudata) { ARKodeGetLastStep(arkode_mem, &hcur); // Call RHS function - // run_rhs_f(t); run_rhs_f(t); // Save derivatives to dudata @@ -748,7 +745,7 @@ void ArkodeMRISolver::rhs_f(BoutReal t, BoutReal* udata, BoutReal* dudata) { void ArkodeMRISolver::pre_s(BoutReal t, BoutReal gamma, BoutReal delta, BoutReal* udata, BoutReal* rvec, BoutReal* zvec) { - TRACE("Running preconditioner: ArkodeMRISolver::pre({:e})", t); + TRACE("Running slow preconditioner: ArkodeMRISolver::pre_s({:e})", t); const BoutReal tstart = bout::globals::mpi->MPI_Wtime(); @@ -776,7 +773,7 @@ void ArkodeMRISolver::pre_s(BoutReal t, BoutReal gamma, BoutReal delta, BoutReal void ArkodeMRISolver::pre_f(BoutReal t, BoutReal gamma, BoutReal delta, BoutReal* udata, BoutReal* rvec, BoutReal* zvec) { - TRACE("Running preconditioner: ArkodeMRISolver::pre({:e})", t); + TRACE("Running fast preconditioner: ArkodeMRISolver::pre_f({:e})", t); const BoutReal tstart = bout::globals::mpi->MPI_Wtime(); @@ -802,50 +799,6 @@ void ArkodeMRISolver::pre_f(BoutReal t, BoutReal gamma, BoutReal delta, BoutReal pre_ncalls_s++; } -/************************************************************************** - * Jacobian-vector multiplication functions - **************************************************************************/ - -void ArkodeMRISolver::jac_s(BoutReal t, BoutReal* ydata, BoutReal* vdata, BoutReal* Jvdata) { - TRACE("Running Jacobian: ArkodeMRISolver::jac({:e})", t); - - if (not hasJacobian()) { - throw BoutException("No jacobian function supplied!\n"); - } - - // Load state from ydate - load_vars(ydata); - - // Load vector to be multiplied into F_vars - load_derivs(vdata); - - // Call function - runJacobian(t); - - // Save Jv from vars - save_derivs(Jvdata); -} - -void ArkodeMRISolver::jac_f(BoutReal t, BoutReal* ydata, BoutReal* vdata, BoutReal* Jvdata) { - TRACE("Running Jacobian: ArkodeMRISolver::jac({:e})", t); - - if (not hasJacobian()) { - throw BoutException("No jacobian function supplied!\n"); - } - - // Load state from ydate - load_vars(ydata); - - // Load vector to be multiplied into F_vars - load_derivs(vdata); - - // Call function - runJacobian(t); - - // Save Jv from vars - save_derivs(Jvdata); -} - /************************************************************************** * ARKODE explicit RHS functions **************************************************************************/ @@ -968,7 +921,7 @@ int arkode_s_pre(BoutReal t, N_Vector yy, N_Vector UNUSED(yp), N_Vector rvec, N_ auto* s = static_cast(user_data); - // Calculate residuals + // Run slow preconditioner s->pre_s(t, gamma, delta, udata, rdata, zdata); return 0; @@ -982,7 +935,7 @@ int arkode_f_pre(BoutReal t, N_Vector yy, N_Vector UNUSED(yp), N_Vector rvec, N_ auto* s = static_cast(user_data); - // Calculate residuals + // Run fast preconditioner s->pre_f(t, gamma, delta, udata, rdata, zdata); return 0; diff --git a/src/solver/impls/arkode/arkode_mri.hxx b/src/solver/impls/arkode/arkode_mri.hxx index 9cd33c2df8..e1bd167db9 100644 --- a/src/solver/impls/arkode/arkode_mri.hxx +++ b/src/solver/impls/arkode/arkode_mri.hxx @@ -1,6 +1,6 @@ /************************************************************************** * Interface to ARKODE MRI solver - * NOTE: ARKode is currently in beta testing so use with cautious optimism + * NOTE: ARKODE is currently in beta testing so use with cautious optimism * * NOTE: Only one solver can currently be compiled in * diff --git a/tests/integrated/test-kpr_mri/README.md b/tests/integrated/test-kpr_mri/README.md index ffcef3ded9..2233782dd8 100644 --- a/tests/integrated/test-kpr_mri/README.md +++ b/tests/integrated/test-kpr_mri/README.md @@ -14,9 +14,9 @@ test-kpr_mri f(t) = sqrt(1+r(t)), g(t) = sqrt(2+s(t)). We use the parameters: - e = 0.5 (fast/slow coupling strength) [default] - G = -100 (stiffness at slow time scale) [default] - w = 100 (time-scale separation factor) [default] + e = 0.5 (fast/slow coupling strength) + G = -100 (stiffness at slow time scale) + w = 100 (time-scale separation factor) The stiffness of the slow time scale is essentially determined by G, for |G| > 50 it is 'stiff' and ideally suited to a diff --git a/tests/integrated/test-kpr_mri/test_kpr_mri.cxx b/tests/integrated/test-kpr_mri/test_kpr_mri.cxx index 250f874fe6..a0e7bcf045 100644 --- a/tests/integrated/test-kpr_mri/test_kpr_mri.cxx +++ b/tests/integrated/test-kpr_mri/test_kpr_mri.cxx @@ -6,7 +6,7 @@ #include #include -// A simple phyics model with a manufactured true solution +// A simple phyics model with an analytical solution // class TestSolver : public PhysicsModel { public: @@ -22,7 +22,7 @@ class TestSolver : public PhysicsModel { f = sqrt(3.0/2.0); g = sqrt(3.0); - + setSplitOperatorMRI(); return 0; @@ -50,16 +50,8 @@ class TestSolver : public PhysicsModel { return 0; } - int rhs_fe(BoutReal UNUSED(t)) override { - - ddt(f) = 0.0; - ddt(g) = 0.0; - - return 0; - } - - int rhs_fi(BoutReal t) override { - /* fill in the fast implicit RHS function: + int rhs_fe(BoutReal t) override { + /* fill in the fast explicit RHS function: [0 0]*[(-1+f^2-0.5*cos(t))/(2*f)] + [ 0 ] [e -1] [(-2+g^2-cos(w*t))/(2*g) ] [-w*sin(w*t)/(2*sqrt(2+cos(w*t)))] */ BoutReal tmp1 = (-1.0 + f(1,1,0) * f(1,1,0) - 0.5*cos(t)) / (2.0 * f(1,1,0)); @@ -70,6 +62,14 @@ class TestSolver : public PhysicsModel { return 0; } + int rhs_fi(BoutReal UNUSED(t)) override { + + ddt(f) = 0.0; + ddt(g) = 0.0; + + return 0; + } + int rhs_s(BoutReal t) override { /* fill in the RHS function: [G e]*[(-1+f^2-0.5*cos(t))/(2*f)] + [-0.5*sin(t)/(2*f)] @@ -94,19 +94,6 @@ class TestSolver : public PhysicsModel { return 0; } - // int rhs(BoutReal t) override { - // /* fill in the RHS function: - // [G e]*[(-1+f^2-0.5*cos(t))/(2*f)] + [-0.5*sin(t)/(2*f) ] - // [e -1] [(-2+g^2-cos(w*t))/(2*g) ] [-w*sin(w*t)/(2*sqrt(2+cos(w*t)))] */ - // BoutReal tmp1 = (-1.0 + f(1,1,0) * f(1,1,0) - 0.5*cos(t)) / (2.0 * f(1,1,0)); - // BoutReal tmp2 = (-2.0 + g(1,1,0) * g(1,1,0) - cos(w*t)) / (2.0 * g(1,1,0)); - - // ddt(f) = G * tmp1 + e * tmp2 - 0.5*sin(t) / (2.0 * f(1,1,0)); - // ddt(g) = e * tmp1 - tmp2 - w * sin(w*t) / (2.0 * sqrt(2.0 + cos(w * t))); - - // return 0; - // } - bool check_solution(BoutReal atol, BoutReal t) { // Return true if correct solution return ((std::abs(sqrt(0.5*cos(t) + 1.0) - f(1,1,0)) < atol) and (std::abs(sqrt(cos(w*t) + 2.0) - g(1,1,0)) < atol)); @@ -117,7 +104,7 @@ class TestSolver : public PhysicsModel { /* Compute the error with the true solution: f(t) = sqrt(0.5*cos(t) + 1.0) g(t) = sqrt(cos(w*t) + 2.0) */ - return sqrt( pow(sqrt(0.5*cos(t) + 1.0) - f(1,1,0), 2.0) + + return sqrt( pow(sqrt(0.5*cos(t) + 1.0) - f(1,1,0), 2.0) + pow(sqrt(cos(w*t) + 2.0) - g(1,1,0), 2.0)); } @@ -163,7 +150,7 @@ int main(int argc, char** argv) { // name for an option with inconsistent defaults auto options = Options::getRoot()->getSection("arkode_mri"); auto solver = std::unique_ptr{Solver::create("arkode_mri", options)}; - + TestSolver model{}; solver->setModel(&model); From 87baf52304d5bf012ca75dad0f370386759bd811 Mon Sep 17 00:00:00 2001 From: maggul Date: Thu, 9 Jan 2025 14:09:31 -0600 Subject: [PATCH 06/28] adaptive step and SUNDIALS release updates --- src/solver/impls/arkode/arkode_mri.cxx | 23 ++++++++++++----------- src/solver/impls/arkode/arkode_mri.hxx | 7 +++++-- 2 files changed, 17 insertions(+), 13 deletions(-) diff --git a/src/solver/impls/arkode/arkode_mri.cxx b/src/solver/impls/arkode/arkode_mri.cxx index d68748411f..ff9d800bcf 100644 --- a/src/solver/impls/arkode/arkode_mri.cxx +++ b/src/solver/impls/arkode/arkode_mri.cxx @@ -95,12 +95,15 @@ ArkodeMRISolver::ArkodeMRISolver(Options* opts) .doc("Use linear implicit solver (only evaluates jacobian inversion once)") .withDefault(false)), fixed_step((*options)["fixed_step"] - .doc("Solve explicit portion in fixed timestep mode. NOTE: This is " + .doc("Solve both fast and slow time scales using fixed time step sizes. NOTE: This is " "not recommended except for code comparison") .withDefault(false)), - order((*options)["order"].doc("Order of internal step").withDefault(4)), + inner_timestep((*options)["inner_timestep"] + .doc("Fixed step size to use for inner solver when running in fixed time step mode.") + .withDefault(1.0e-4)), + order((*options)["order"].doc("Order of internal step").withDefault(3)), abstol((*options)["atol"].doc("Absolute tolerance").withDefault(1.0e-12)), - reltol((*options)["rtol"].doc("Relative tolerance").withDefault(1.0e-5)), + reltol((*options)["rtol"].doc("Relative tolerance").withDefault(1.0e-10)), use_vector_abstol((*options)["use_vector_abstol"] .doc("Use separate absolute tolerance for each field") .withDefault(false)), @@ -232,8 +235,8 @@ int ArkodeMRISolver::init() { throw BoutException("ARKodeSetOrder failed\n"); } - if (ARKStepCreateMRIStepInnerStepper(inner_arkode_mem, &inner_stepper) != ARK_SUCCESS) { - throw BoutException("ARKStepCreateMRIStepInnerStepper failed\n"); + if (ARKodeCreateMRIStepInnerStepper(inner_arkode_mem, &inner_stepper) != ARK_SUCCESS) { + throw BoutException("ARKodeCreateMRIStepInnerStepper failed\n"); } // Initialize the slow integrator. Specify the explicit slow right-hand side @@ -285,10 +288,6 @@ int ArkodeMRISolver::init() { throw BoutException("ARKodeSetOrder failed\n"); } - if (ARKodeSetFixedStep(arkode_mem, 0.001) != ARK_SUCCESS) { - throw BoutException("ARKodeSetFixedStep failed\n"); - } - if (use_vector_abstol) { std::vector f2dtols; f2dtols.reserve(f2d.size()); @@ -518,7 +517,8 @@ int ArkodeMRISolver::run() { long int temp_long_int, temp_long_int2; ARKodeGetNumSteps(arkode_mem, &temp_long_int); nsteps = int(temp_long_int); - MRIStepGetNumRhsEvals(arkode_mem, &temp_long_int, &temp_long_int2); //Change after the release + ARKodeGetNumRhsEvals(arkode_mem, 0, &temp_long_int); + ARKodeGetNumRhsEvals(arkode_mem, 1, &temp_long_int2); nfe_evals = int(temp_long_int); nfi_evals = int(temp_long_int2); if (treatment == MRI_Treatment::ImEx or treatment == MRI_Treatment::Implicit) { @@ -532,7 +532,8 @@ int ArkodeMRISolver::run() { ARKodeGetNumSteps(inner_arkode_mem, &temp_long_int); inner_nsteps = int(temp_long_int); - ARKStepGetNumRhsEvals(inner_arkode_mem, &temp_long_int, &temp_long_int2); //Change after the release + ARKodeGetNumRhsEvals(inner_arkode_mem, 0, &temp_long_int); + ARKodeGetNumRhsEvals(inner_arkode_mem, 1, &temp_long_int2); inner_nfe_evals = int(temp_long_int); inner_nfi_evals = int(temp_long_int2); if (inner_treatment == MRI_Treatment::ImEx or inner_treatment == MRI_Treatment::Implicit) { diff --git a/src/solver/impls/arkode/arkode_mri.hxx b/src/solver/impls/arkode/arkode_mri.hxx index e1bd167db9..2fad1652dd 100644 --- a/src/solver/impls/arkode/arkode_mri.hxx +++ b/src/solver/impls/arkode/arkode_mri.hxx @@ -111,9 +111,12 @@ private: /// Use linear implicit solver (only evaluates jacobian inversion once) bool set_linear; bool inner_set_linear; - /// Solve explicit portion in fixed timestep mode. NOTE: This is not recommended except - /// for code comparison + /// Solve both fast and slow portion in fixed timestep mode. + /// NOTE: This is not recommended except for code comparison bool fixed_step; + /// Fixed step size to use for inner solver when running in fixed + /// time step mode. + BoutReal inner_timestep; /// Order of the internal step int order; /// Absolute tolerance From 9c5d88f0cec589630009571b581c21defbd1d8ca Mon Sep 17 00:00:00 2001 From: maggul Date: Thu, 16 Jan 2025 05:42:54 -0600 Subject: [PATCH 07/28] added SUNDIALS version and removed unnecessary flags --- src/solver/impls/arkode/arkode_mri.hxx | 3 +++ src/solver/solver.cxx | 16 ++-------------- 2 files changed, 5 insertions(+), 14 deletions(-) diff --git a/src/solver/impls/arkode/arkode_mri.hxx b/src/solver/impls/arkode/arkode_mri.hxx index 2fad1652dd..a91aac5a81 100644 --- a/src/solver/impls/arkode/arkode_mri.hxx +++ b/src/solver/impls/arkode/arkode_mri.hxx @@ -32,6 +32,8 @@ #include "bout/build_config.hxx" #include "bout/solver.hxx" +#if ((SUNDIALS_VERSION_MAJOR == 7 && SUNDIALS_VERSION_MINOR >= 2) || SUNDIALS_VERSION_MAJOR > 8) +#else #if not BOUT_HAS_ARKODE namespace { @@ -165,4 +167,5 @@ private: }; #endif // BOUT_HAS_ARKODE +#endif // SUNDIALS_VERSION CHECK #endif // BOUT_ARKODE_MRI_SOLVER_H diff --git a/src/solver/solver.cxx b/src/solver/solver.cxx index 2aad36dd5d..356a04287c 100644 --- a/src/solver/solver.cxx +++ b/src/solver/solver.cxx @@ -2,7 +2,7 @@ * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu * * Contact: Ben Dudson, bd512@york.ac.uk - * + * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -1388,12 +1388,6 @@ int Solver::run_rhs_se(BoutReal t, bool linear) { Timer timer("rhs"); - if (first_rhs_call) { - // Ensure that nonlinear terms are calculated on first call - linear = false; - first_rhs_call = false; - } - pre_rhs(t); status = model->runRHS_se(t, linear); post_rhs(t); @@ -1432,12 +1426,6 @@ int Solver::run_rhs_fe(BoutReal t, bool linear) { Timer timer("rhs"); - if (first_rhs_call) { - // Ensure that nonlinear terms are calculated on first call - linear = false; - first_rhs_call = false; - } - pre_rhs(t); status = model->runRHS_fe(t, linear); post_rhs(t); @@ -1540,7 +1528,7 @@ int Solver::run_rhs(BoutReal t, bool linear) { Array tmp3(nv); save_vars(tmp.begin()); // Copy variables into tmp - + pre_rhs(t); status = model->runRHS_fe(t, linear); post_rhs(t); // Check variables, apply boundary conditions From d23297b96d6d6fd4480b8479330f88cb38900f9a Mon Sep 17 00:00:00 2001 From: maggul Date: Thu, 16 Jan 2025 06:01:15 -0600 Subject: [PATCH 08/28] moved sundials_config.h up before the version check --- src/solver/impls/arkode/arkode_mri.hxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solver/impls/arkode/arkode_mri.hxx b/src/solver/impls/arkode/arkode_mri.hxx index a91aac5a81..f915b0df64 100644 --- a/src/solver/impls/arkode/arkode_mri.hxx +++ b/src/solver/impls/arkode/arkode_mri.hxx @@ -31,9 +31,9 @@ #include "bout/build_config.hxx" #include "bout/solver.hxx" +#include #if ((SUNDIALS_VERSION_MAJOR == 7 && SUNDIALS_VERSION_MINOR >= 2) || SUNDIALS_VERSION_MAJOR > 8) -#else #if not BOUT_HAS_ARKODE namespace { @@ -48,7 +48,6 @@ RegisterUnavailableSolver #include "bout/sundials_backports.hxx" #include -#include #include #include @@ -167,5 +166,6 @@ private: }; #endif // BOUT_HAS_ARKODE +#else #endif // SUNDIALS_VERSION CHECK #endif // BOUT_ARKODE_MRI_SOLVER_H From 340cdc00e1f8447899f7b2cdabdeae1df679c465 Mon Sep 17 00:00:00 2001 From: maggul Date: Thu, 16 Jan 2025 13:09:38 -0600 Subject: [PATCH 09/28] update the SUNDIALS version check --- src/solver/impls/arkode/arkode_mri.hxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/solver/impls/arkode/arkode_mri.hxx b/src/solver/impls/arkode/arkode_mri.hxx index f915b0df64..098548b40c 100644 --- a/src/solver/impls/arkode/arkode_mri.hxx +++ b/src/solver/impls/arkode/arkode_mri.hxx @@ -31,9 +31,9 @@ #include "bout/build_config.hxx" #include "bout/solver.hxx" -#include +#include "bout/sundials_backports.hxx" -#if ((SUNDIALS_VERSION_MAJOR == 7 && SUNDIALS_VERSION_MINOR >= 2) || SUNDIALS_VERSION_MAJOR > 8) +#if SUNDIALS_VERSION_AT_LEAST(7, 2, 0) #if not BOUT_HAS_ARKODE namespace { @@ -45,11 +45,10 @@ RegisterUnavailableSolver #include "bout/bout_enum_class.hxx" #include "bout/bout_types.hxx" -#include "bout/sundials_backports.hxx" #include #include - +#include #include class ArkodeMRISolver; @@ -169,3 +168,4 @@ private: #else #endif // SUNDIALS_VERSION CHECK #endif // BOUT_ARKODE_MRI_SOLVER_H + From 528135c4922187effc53249fe87e29c809c7dd41 Mon Sep 17 00:00:00 2001 From: maggul Date: Fri, 17 Jan 2025 11:29:52 -0600 Subject: [PATCH 10/28] Updated include order --- src/solver/impls/arkode/arkode_mri.hxx | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/solver/impls/arkode/arkode_mri.hxx b/src/solver/impls/arkode/arkode_mri.hxx index 098548b40c..3e5858c7bf 100644 --- a/src/solver/impls/arkode/arkode_mri.hxx +++ b/src/solver/impls/arkode/arkode_mri.hxx @@ -29,11 +29,9 @@ #ifndef BOUT_ARKODE_MRI_SOLVER_H #define BOUT_ARKODE_MRI_SOLVER_H -#include "bout/build_config.hxx" +#include "bout/build_defines.hxx" #include "bout/solver.hxx" -#include "bout/sundials_backports.hxx" -#if SUNDIALS_VERSION_AT_LEAST(7, 2, 0) #if not BOUT_HAS_ARKODE namespace { @@ -45,6 +43,10 @@ RegisterUnavailableSolver #include "bout/bout_enum_class.hxx" #include "bout/bout_types.hxx" +#include "bout/region.hxx" +#include "bout/sundials_backports.hxx" + +#if SUNDIALS_VERSION_AT_LEAST(7, 2, 0) #include #include @@ -164,8 +166,8 @@ private: sundials::Context suncontext; }; -#endif // BOUT_HAS_ARKODE #else #endif // SUNDIALS_VERSION CHECK +#endif // BOUT_HAS_ARKODE #endif // BOUT_ARKODE_MRI_SOLVER_H From 250b1f5093b320c5bc6c119c886200d808974860 Mon Sep 17 00:00:00 2001 From: maggul Date: Fri, 17 Jan 2025 11:50:11 -0600 Subject: [PATCH 11/28] Fixed a function name --- src/solver/impls/arkode/arkode.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solver/impls/arkode/arkode.cxx b/src/solver/impls/arkode/arkode.cxx index cb474878b6..23883cc043 100644 --- a/src/solver/impls/arkode/arkode.cxx +++ b/src/solver/impls/arkode/arkode.cxx @@ -319,9 +319,9 @@ int ArkodeSolver::init() { throw BoutException("Invalid adap_method\n"); } - if (ARKodeSetAdaptivityMethod(arkode_mem, adap_method_int, 1, 1, nullptr) + if (ARKStepSetAdaptivityMethod(arkode_mem, adap_method_int, 1, 1, nullptr) != ARK_SUCCESS) { - throw BoutException("ARKodeSetAdaptivityMethod failed\n"); + throw BoutException("ARKStepSetAdaptivityMethod failed\n"); } #endif From 052be870cc64b222e320f2f7f933de5a08e40364 Mon Sep 17 00:00:00 2001 From: maggul Date: Fri, 17 Jan 2025 12:55:01 -0600 Subject: [PATCH 12/28] update version check --- src/solver/impls/arkode/arkode_mri.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solver/impls/arkode/arkode_mri.hxx b/src/solver/impls/arkode/arkode_mri.hxx index 3e5858c7bf..f18a41173d 100644 --- a/src/solver/impls/arkode/arkode_mri.hxx +++ b/src/solver/impls/arkode/arkode_mri.hxx @@ -46,7 +46,7 @@ RegisterUnavailableSolver #include "bout/region.hxx" #include "bout/sundials_backports.hxx" -#if SUNDIALS_VERSION_AT_LEAST(7, 2, 0) +#if ((SUNDIALS_VERSION_MAJOR == 7 && SUNDIALS_VERSION_MINOR >= 2) || SUNDIALS_VERSION_MAJOR > 8) #include #include From 366aad03ff967a5f1bb51319b30a0964f1e81ba8 Mon Sep 17 00:00:00 2001 From: maggul Date: Fri, 17 Jan 2025 23:48:08 -0600 Subject: [PATCH 13/28] update header files --- src/solver/impls/arkode/arkode_mri.cxx | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/src/solver/impls/arkode/arkode_mri.cxx b/src/solver/impls/arkode/arkode_mri.cxx index ff9d800bcf..66bd20959e 100644 --- a/src/solver/impls/arkode/arkode_mri.cxx +++ b/src/solver/impls/arkode/arkode_mri.cxx @@ -23,32 +23,41 @@ * **************************************************************************/ -#include "bout/build_config.hxx" +#include "bout/build_defines.hxx" + +#if BOUT_HAS_ARKODE #include "arkode_mri.hxx" -#if BOUT_HAS_ARKODE +#if ((SUNDIALS_VERSION_MAJOR == 7 && SUNDIALS_VERSION_MINOR >= 2) || SUNDIALS_VERSION_MAJOR > 8) -#include "bout/bout_enum_class.hxx" +#include "bout/assert.hxx" +#include "bout/bout_types.hxx" #include "bout/boutcomm.hxx" #include "bout/boutexception.hxx" +#include "bout/field2d.hxx" #include "bout/field3d.hxx" +#include "bout/globals.hxx" #include "bout/mesh.hxx" +#include "bout/mpi_wrapper.hxx" #include "bout/msg_stack.hxx" #include "bout/options.hxx" #include "bout/output.hxx" +#include "bout/solver.hxx" +#include "bout/sundials_backports.hxx" #include "bout/unused.hxx" -#include "bout/utils.hxx" +#include #include #include +#include + #include -#include #include +#include #include - -class Field2D; +#include // NOLINTBEGIN(readability-identifier-length) namespace { @@ -988,4 +997,5 @@ void ArkodeMRISolver::loop_abstol_values_op(Ind2D UNUSED(i2d), BoutReal* abstolv } } -#endif +#endif // SUNDIALS_VERSION CHECK +#endif // BOUT_HAS_ARKODE \ No newline at end of file From 92db7e1b05364841728f73eb7e6e952924aabf62 Mon Sep 17 00:00:00 2001 From: maggul Date: Sat, 18 Jan 2025 00:26:46 -0600 Subject: [PATCH 14/28] fixed the arkode_mri solver name --- include/bout/solver.hxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index 33bb7c8e12..ff3b242eaf 100644 --- a/include/bout/solver.hxx +++ b/include/bout/solver.hxx @@ -87,7 +87,7 @@ constexpr auto SOLVEREULER = "euler"; constexpr auto SOLVERRK3SSP = "rk3ssp"; constexpr auto SOLVERPOWER = "power"; constexpr auto SOLVERARKODE = "arkode"; -constexpr auto SOLVERARKODEMRI = "arkodemri"; +constexpr auto SOLVERARKODEMRI = "arkode_mri"; constexpr auto SOLVERIMEXBDF2 = "imexbdf2"; constexpr auto SOLVERSNES = "snes"; constexpr auto SOLVERRKGENERIC = "rkgeneric"; @@ -456,7 +456,7 @@ protected: int run_rhs_fe(BoutReal t, bool linear = false); int run_rhs_fi(BoutReal t, bool linear = false); int run_rhs_s(BoutReal t, bool linear = false); - int run_rhs_f(BoutReal t, bool linear = false); + int run_rhs_f(BoutReal t, bool linear = false); int run_rhs(BoutReal t, bool linear = false); /// Calculate only the convective parts int run_convective(BoutReal t, bool linear = false); From c680123fd996a5a70067f6230a31de16afae4099 Mon Sep 17 00:00:00 2001 From: maggul Date: Mon, 20 Jan 2025 12:40:00 -0600 Subject: [PATCH 15/28] version check update --- src/solver/impls/arkode/arkode_mri.cxx | 2 +- src/solver/impls/arkode/arkode_mri.hxx | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solver/impls/arkode/arkode_mri.cxx b/src/solver/impls/arkode/arkode_mri.cxx index 66bd20959e..e1e5ed14f3 100644 --- a/src/solver/impls/arkode/arkode_mri.cxx +++ b/src/solver/impls/arkode/arkode_mri.cxx @@ -29,7 +29,7 @@ #include "arkode_mri.hxx" -#if ((SUNDIALS_VERSION_MAJOR == 7 && SUNDIALS_VERSION_MINOR >= 2) || SUNDIALS_VERSION_MAJOR > 8) +SUNDIALS_VERSION_AT_LEAST(7, 2, 0) #include "bout/assert.hxx" #include "bout/bout_types.hxx" diff --git a/src/solver/impls/arkode/arkode_mri.hxx b/src/solver/impls/arkode/arkode_mri.hxx index f18a41173d..fd0cac3008 100644 --- a/src/solver/impls/arkode/arkode_mri.hxx +++ b/src/solver/impls/arkode/arkode_mri.hxx @@ -46,7 +46,7 @@ RegisterUnavailableSolver #include "bout/region.hxx" #include "bout/sundials_backports.hxx" -#if ((SUNDIALS_VERSION_MAJOR == 7 && SUNDIALS_VERSION_MINOR >= 2) || SUNDIALS_VERSION_MAJOR > 8) +SUNDIALS_VERSION_AT_LEAST(7, 2, 0) #include #include From 41830beb7e44b3e2e0a37a9fd88bc0c50e9d7fe8 Mon Sep 17 00:00:00 2001 From: maggul Date: Mon, 27 Jan 2025 16:19:27 -0600 Subject: [PATCH 16/28] removed unnecessary endif --- src/solver/impls/arkode/arkode_mri.cxx | 1 - src/solver/impls/arkode/arkode_mri.hxx | 1 - 2 files changed, 2 deletions(-) diff --git a/src/solver/impls/arkode/arkode_mri.cxx b/src/solver/impls/arkode/arkode_mri.cxx index e1e5ed14f3..f197d1065c 100644 --- a/src/solver/impls/arkode/arkode_mri.cxx +++ b/src/solver/impls/arkode/arkode_mri.cxx @@ -997,5 +997,4 @@ void ArkodeMRISolver::loop_abstol_values_op(Ind2D UNUSED(i2d), BoutReal* abstolv } } -#endif // SUNDIALS_VERSION CHECK #endif // BOUT_HAS_ARKODE \ No newline at end of file diff --git a/src/solver/impls/arkode/arkode_mri.hxx b/src/solver/impls/arkode/arkode_mri.hxx index fd0cac3008..c59a3da002 100644 --- a/src/solver/impls/arkode/arkode_mri.hxx +++ b/src/solver/impls/arkode/arkode_mri.hxx @@ -167,7 +167,6 @@ private: }; #else -#endif // SUNDIALS_VERSION CHECK #endif // BOUT_HAS_ARKODE #endif // BOUT_ARKODE_MRI_SOLVER_H From aa7028e4702ba0c8fea1e58f4db37c2d791ea4e4 Mon Sep 17 00:00:00 2001 From: maggul Date: Mon, 27 Jan 2025 16:38:52 -0600 Subject: [PATCH 17/28] sundials version check resolved --- src/solver/impls/arkode/arkode_mri.cxx | 4 +++- src/solver/impls/arkode/arkode_mri.hxx | 3 ++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/src/solver/impls/arkode/arkode_mri.cxx b/src/solver/impls/arkode/arkode_mri.cxx index f197d1065c..bad16362da 100644 --- a/src/solver/impls/arkode/arkode_mri.cxx +++ b/src/solver/impls/arkode/arkode_mri.cxx @@ -29,7 +29,7 @@ #include "arkode_mri.hxx" -SUNDIALS_VERSION_AT_LEAST(7, 2, 0) +#if SUNDIALS_VERSION_AT_LEAST(7, 2, 0) #include "bout/assert.hxx" #include "bout/bout_types.hxx" @@ -997,4 +997,6 @@ void ArkodeMRISolver::loop_abstol_values_op(Ind2D UNUSED(i2d), BoutReal* abstolv } } +#else +#endif // SUNDIALS_VERSION CHECK #endif // BOUT_HAS_ARKODE \ No newline at end of file diff --git a/src/solver/impls/arkode/arkode_mri.hxx b/src/solver/impls/arkode/arkode_mri.hxx index c59a3da002..3e5858c7bf 100644 --- a/src/solver/impls/arkode/arkode_mri.hxx +++ b/src/solver/impls/arkode/arkode_mri.hxx @@ -46,7 +46,7 @@ RegisterUnavailableSolver #include "bout/region.hxx" #include "bout/sundials_backports.hxx" -SUNDIALS_VERSION_AT_LEAST(7, 2, 0) +#if SUNDIALS_VERSION_AT_LEAST(7, 2, 0) #include #include @@ -167,6 +167,7 @@ private: }; #else +#endif // SUNDIALS_VERSION CHECK #endif // BOUT_HAS_ARKODE #endif // BOUT_ARKODE_MRI_SOLVER_H From 52802c1cc50ccd17f322ec4eb5dc86a55c7cf519 Mon Sep 17 00:00:00 2001 From: maggul Date: Fri, 7 Feb 2025 14:13:41 -0600 Subject: [PATCH 18/28] flag and test-kpr updates --- include/bout/solver.hxx | 2 ++ src/solver/impls/arkode/arkode_mri.cxx | 20 ++++++++++--------- src/solver/solver.cxx | 8 ++++---- tests/integrated/test-kpr_mri/CMakeLists.txt | 2 +- .../integrated/test-kpr_mri/test_kpr_mri.cxx | 10 ++++++++++ 5 files changed, 28 insertions(+), 14 deletions(-) diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index ff3b242eaf..d2e63818f1 100644 --- a/include/bout/solver.hxx +++ b/include/bout/solver.hxx @@ -446,6 +446,8 @@ protected: bool initialised{false}; /// If calling user RHS for the first time bool first_rhs_call{true}; + bool first_rhs_s_call{false}; + bool first_rhs_f_call{false}; /// Current simulation time BoutReal simtime{0.0}; diff --git a/src/solver/impls/arkode/arkode_mri.cxx b/src/solver/impls/arkode/arkode_mri.cxx index bad16362da..915b24580b 100644 --- a/src/solver/impls/arkode/arkode_mri.cxx +++ b/src/solver/impls/arkode/arkode_mri.cxx @@ -648,8 +648,8 @@ void ArkodeMRISolver::rhs_se(BoutReal t, BoutReal* udata, BoutReal* dudata) { load_vars(udata); // Get the current timestep - // Note: ARKodeGetCurrentStep updated too late in older versions - ARKodeGetLastStep(arkode_mem, &hcur); + // TO DO: Check to identify which time step is needed current/last, then correct accordingly + ARKodeGetCurrentStep(arkode_mem, &hcur); // Call RHS function run_rhs_se(t); @@ -666,7 +666,8 @@ void ArkodeMRISolver::rhs_si(BoutReal t, BoutReal* udata, BoutReal* dudata) { TRACE("Running RHS: ArkodeMRISolver::rhs_si({:e})", t); load_vars(udata); - ARKodeGetLastStep(arkode_mem, &hcur); + // TO DO: Check to identify which time step is needed current/last, then correct accordingly + ARKodeGetCurrentStep(arkode_mem, &hcur); // Call Implicit RHS function run_rhs_si(t); save_derivs(dudata); @@ -683,8 +684,8 @@ void ArkodeMRISolver::rhs_fe(BoutReal t, BoutReal* udata, BoutReal* dudata) { load_vars(udata); // Get the current timestep - // Note: ARKodeGetCurrentStep updated too late in older versions - ARKodeGetLastStep(arkode_mem, &hcur); + // TO DO: Check to identify which time step is needed current/last, then correct accordingly + ARKodeGetCurrentStep(arkode_mem, &hcur); // Call RHS function run_rhs_fe(t); @@ -701,7 +702,8 @@ void ArkodeMRISolver::rhs_fi(BoutReal t, BoutReal* udata, BoutReal* dudata) { TRACE("Running RHS: ArkodeMRISolver::rhs_si({:e})", t); load_vars(udata); - ARKodeGetLastStep(arkode_mem, &hcur); + // TO DO: Check to identify which time step is needed current/last, then correct accordingly + ARKodeGetCurrentStep(arkode_mem, &hcur); // Call Implicit RHS function run_rhs_fi(t); save_derivs(dudata); @@ -718,8 +720,8 @@ void ArkodeMRISolver::rhs_s(BoutReal t, BoutReal* udata, BoutReal* dudata) { load_vars(udata); // Get the current timestep - // Note: ARKodeGetCurrentStep updated too late in older versions - ARKodeGetLastStep(arkode_mem, &hcur); + // TO DO: Check to identify which time step is needed current/last, then correct accordingly + ARKodeGetCurrentStep(arkode_mem, &hcur); // Call RHS function run_rhs_s(t); @@ -740,7 +742,7 @@ void ArkodeMRISolver::rhs_f(BoutReal t, BoutReal* udata, BoutReal* dudata) { // Get the current timestep // Note: ARKodeGetCurrentStep updated too late in older versions - ARKodeGetLastStep(arkode_mem, &hcur); + ARKodeGetCurrentStep(arkode_mem, &hcur); // Call RHS function run_rhs_f(t); diff --git a/src/solver/solver.cxx b/src/solver/solver.cxx index e9ee7be6cc..2a7d16f319 100644 --- a/src/solver/solver.cxx +++ b/src/solver/solver.cxx @@ -1404,10 +1404,10 @@ int Solver::run_rhs_si(BoutReal t, bool linear) { Timer timer("rhs"); - if (first_rhs_call) { + if (first_rhs_s_call) { // Ensure that nonlinear terms are calculated on first call linear = false; - first_rhs_call = false; + first_rhs_s_call = false; } pre_rhs(t); @@ -1442,10 +1442,10 @@ int Solver::run_rhs_fi(BoutReal t, bool linear) { Timer timer("rhs"); - if (first_rhs_call) { + if (first_rhs_f_call) { // Ensure that nonlinear terms are calculated on first call linear = false; - first_rhs_call = false; + first_rhs_f_call = false; } pre_rhs(t); diff --git a/tests/integrated/test-kpr_mri/CMakeLists.txt b/tests/integrated/test-kpr_mri/CMakeLists.txt index 66a4474bb0..45523803db 100644 --- a/tests/integrated/test-kpr_mri/CMakeLists.txt +++ b/tests/integrated/test-kpr_mri/CMakeLists.txt @@ -1,2 +1,2 @@ -bout_add_integrated_test(test_kpr_mri SOURCES test_kpr_mri.cxx +bout_add_integrated_test(test_kpr_mri SOURCES test_kpr_mri.cxx REQUIRES BOUT_HAS_SUNDIALS) diff --git a/tests/integrated/test-kpr_mri/test_kpr_mri.cxx b/tests/integrated/test-kpr_mri/test_kpr_mri.cxx index a0e7bcf045..fcb590dee0 100644 --- a/tests/integrated/test-kpr_mri/test_kpr_mri.cxx +++ b/tests/integrated/test-kpr_mri/test_kpr_mri.cxx @@ -1,3 +1,5 @@ +#if SUNDIALS_VERSION_AT_LEAST(7, 2, 0) + #include "bout/physicsmodel.hxx" #include "bout/solver.hxx" @@ -170,3 +172,11 @@ int main(int argc, char** argv) { output_test << " FAILED\n"; return 1; } +#else +// ARKODE-MRI option for BOUT++ +// is not available with this configuration +// do nothing for this test +int main(int argc, char** argv) { + return 0; +} +#endif \ No newline at end of file From 63db0f46220ba374033bbfcb62e0be5585fa2d53 Mon Sep 17 00:00:00 2001 From: maggul Date: Fri, 7 Feb 2025 18:40:08 -0600 Subject: [PATCH 19/28] added required header file --- tests/integrated/test-kpr_mri/test_kpr_mri.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/integrated/test-kpr_mri/test_kpr_mri.cxx b/tests/integrated/test-kpr_mri/test_kpr_mri.cxx index fcb590dee0..50b94ee55f 100644 --- a/tests/integrated/test-kpr_mri/test_kpr_mri.cxx +++ b/tests/integrated/test-kpr_mri/test_kpr_mri.cxx @@ -1,3 +1,5 @@ +#include "bout/sundials_backports.hxx" + #if SUNDIALS_VERSION_AT_LEAST(7, 2, 0) #include "bout/physicsmodel.hxx" From 30c16d91411a20a72e5425787f12ef954fdb7005 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 7 Feb 2025 17:15:07 -0800 Subject: [PATCH 20/28] ARKode-MRI: Register unavailable if sundials version is too low Provides an error message if arkode_mri solver is requested but BOUT++ was built with a version of SUNDIALS lower than 7.2.0 --- src/solver/impls/arkode/arkode_mri.hxx | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/src/solver/impls/arkode/arkode_mri.hxx b/src/solver/impls/arkode/arkode_mri.hxx index 3e5858c7bf..d10794a120 100644 --- a/src/solver/impls/arkode/arkode_mri.hxx +++ b/src/solver/impls/arkode/arkode_mri.hxx @@ -5,7 +5,7 @@ * NOTE: Only one solver can currently be compiled in * ************************************************************************** - * Copyright 2010-2024 BOUT++ contributors + * Copyright 2010-2025 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov * @@ -41,13 +41,20 @@ RegisterUnavailableSolver #else -#include "bout/bout_enum_class.hxx" -#include "bout/bout_types.hxx" -#include "bout/region.hxx" #include "bout/sundials_backports.hxx" -#if SUNDIALS_VERSION_AT_LEAST(7, 2, 0) +#if SUNDIALS_VERSION_LESS_THAN(7, 2, 0) + +namespace { +RegisterUnavailableSolver + registerunavailablearkodemri("arkode_mri", "BOUT++ configured with ARKODE/SUNDIALS version less than 7.2.0"); +} + +#else // SUNDIALS_VERSION check +#include "bout/bout_enum_class.hxx" +#include "bout/bout_types.hxx" +#include "bout/region.hxx" #include #include #include @@ -166,7 +173,6 @@ private: sundials::Context suncontext; }; -#else #endif // SUNDIALS_VERSION CHECK #endif // BOUT_HAS_ARKODE #endif // BOUT_ARKODE_MRI_SOLVER_H From 462cfb3833df18e01d6eea45fd3aa9867feaac82 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 7 Feb 2025 17:36:21 -0800 Subject: [PATCH 21/28] Add fast and slow preconditioner callbacks PhysicsModel now has fast and slow preconditioners, and these can be called by the Arkode MRI solver. --- include/bout/physicsmodel.hxx | 42 +++++++++++++++++--------- include/bout/solver.hxx | 5 +++ src/physics/physicsmodel.cxx | 15 ++++++++- src/solver/impls/arkode/arkode_mri.cxx | 21 ++++++------- src/solver/solver.cxx | 27 ++++++++++++----- 5 files changed, 76 insertions(+), 34 deletions(-) diff --git a/include/bout/physicsmodel.hxx b/include/bout/physicsmodel.hxx index 97c3f09a40..21fe1f15c0 100644 --- a/include/bout/physicsmodel.hxx +++ b/include/bout/physicsmodel.hxx @@ -4,16 +4,10 @@ * @brief Base class for Physics Models * * - * - * Changelog: - * - * 2013-08 Ben Dudson - * * Initial version - * ************************************************************************** - * Copyright 2013 B.D.Dudson + * Copyright 2013-2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -179,11 +173,13 @@ public: /*! * True if this model uses split operators + * RHS = convective + diffusive */ bool splitOperator(); /*! - * True if this model uses split operators + * True if this model uses Multi-Rate Integrator (MRI) split operators + * RHS = rhs_se + rhs_si + rhs_fe + rhs_fi */ bool splitOperatorMRI(); @@ -208,7 +204,9 @@ public: /*! * True if a preconditioner has been defined */ - bool hasPrecon(); + bool hasPrecon() const { return (userprecon != nullptr); } + bool hasPreconFast() const { return (userprecon_f != nullptr); } + bool hasPreconSlow() const { return (userprecon_s != nullptr); } /*! * Run the preconditioner. The system state should be in the @@ -219,7 +217,9 @@ public: * */ int runPrecon(BoutReal t, BoutReal gamma, BoutReal delta); - + int runPreconFast(BoutReal t, BoutReal gamma, BoutReal delta); + int runPreconSlow(BoutReal t, BoutReal gamma, BoutReal delta); + /*! * True if a Jacobian function has been defined */ @@ -338,10 +338,9 @@ protected: /// Specify that this model is split into a convective and diffusive part void setSplitOperator(bool split = true) { splitop = split; } - /// Specify that this model is split into a convective and diffusive part + /// Specify that this model is split into fast and slow parts void setSplitOperatorMRI(bool split = true) { splitopmri = split; } - /// Specify a preconditioner function void setPrecon(preconfunc pset) { userprecon = pset; } template @@ -349,6 +348,19 @@ protected: userprecon = static_cast(preconditioner); } + /// Preconditioner for fast implicit RHS + void setPreconFast(preconfunc pset) { userprecon_f = pset; } + template + void setPreconFast(ModelPreconFunc preconditioner) { + userprecon_f = static_cast(preconditioner); + } + /// Preconditioner for slow implicit RHS + void setPreconSlow(preconfunc pset) { userprecon_s = pset; } + template + void setPreconSlow(ModelPreconFunc preconditioner) { + userprecon_s = static_cast(preconditioner); + } + /// Specify a Jacobian-vector multiply function void setJacobian(jacobianfunc jset) { userjacobian = jset; } template @@ -424,10 +436,12 @@ private: bool restart_enabled{true}; /// Split operator model? bool splitop{false}; - /// Split operator model? + /// MPI fast/slow split operator model? bool splitopmri{false}; /// Pointer to user-supplied preconditioner function preconfunc userprecon{nullptr}; + preconfunc userprecon_f{nullptr}; + preconfunc userprecon_s{nullptr}; /// Pointer to user-supplied Jacobian-vector multiply function jacobianfunc userjacobian{nullptr}; /// True if model already initialised diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index d2e63818f1..5b14fc09ed 100644 --- a/include/bout/solver.hxx +++ b/include/bout/solver.hxx @@ -493,8 +493,13 @@ protected: /// Do we have a user preconditioner? bool hasPreconditioner(); + bool hasPreconditionerFast(); + bool hasPreconditionerSlow(); + /// Run the user preconditioner int runPreconditioner(BoutReal time, BoutReal gamma, BoutReal delta); + int runPreconditionerFast(BoutReal time, BoutReal gamma, BoutReal delta); + int runPreconditionerSlow(BoutReal time, BoutReal gamma, BoutReal delta); /// Do we have a user Jacobian? bool hasJacobian(); diff --git a/src/physics/physicsmodel.cxx b/src/physics/physicsmodel.cxx index 9d2ef84497..8eb729faef 100644 --- a/src/physics/physicsmodel.cxx +++ b/src/physics/physicsmodel.cxx @@ -131,7 +131,6 @@ int PhysicsModel::runDiffusive(BoutReal time, bool linear) { return diffusive(time, linear); } -bool PhysicsModel::hasPrecon() { return (userprecon != nullptr); } int PhysicsModel::runPrecon(BoutReal t, BoutReal gamma, BoutReal delta) { if (!userprecon) { @@ -140,6 +139,20 @@ int PhysicsModel::runPrecon(BoutReal t, BoutReal gamma, BoutReal delta) { return (*this.*userprecon)(t, gamma, delta); } +int PhysicsModel::runPreconFast(BoutReal t, BoutReal gamma, BoutReal delta) { + if (!userprecon_f) { + return 1; + } + return (*this.*userprecon_f)(t, gamma, delta); +} + +int PhysicsModel::runPreconSlow(BoutReal t, BoutReal gamma, BoutReal delta) { + if (!userprecon_s) { + return 1; + } + return (*this.*userprecon_s)(t, gamma, delta); +} + bool PhysicsModel::hasJacobian() { return (userjacobian != nullptr); } int PhysicsModel::runJacobian(BoutReal t) { diff --git a/src/solver/impls/arkode/arkode_mri.cxx b/src/solver/impls/arkode/arkode_mri.cxx index 915b24580b..71b3da2343 100644 --- a/src/solver/impls/arkode/arkode_mri.cxx +++ b/src/solver/impls/arkode/arkode_mri.cxx @@ -4,7 +4,7 @@ * NOTE: ARKODE is still in beta testing so use with cautious optimism * ************************************************************************** - * Copyright 2010-2024 BOUT++ contributors + * Copyright 2010-2025 BOUT++ contributors * * This file is part of BOUT++. * @@ -366,15 +366,15 @@ int ArkodeMRISolver::init() { /// Set Preconditioner if (inner_use_precon) { - if (hasPreconditioner()) { // change to inner_hasPreconditioner when it is available - output.write("\tUsing user-supplied preconditioner for inner solver\n"); + if (hasPreconditionerFast()) { + output.write("\tUsing user-supplied preconditioner for inner (fast) solver\n"); if (ARKodeSetPreconditioner(inner_arkode_mem, nullptr, arkode_f_pre) != ARKLS_SUCCESS) { throw BoutException("ARKodeSetPreconditioner failed for inner solver\n"); } } else { - output.write("\tUsing BBD preconditioner for inner solver\n"); + output.write("\tUsing BBD preconditioner for inner (fast) solver\n"); /// Get options // Compute band_width_default from actually added fields, to allow for multiple @@ -422,7 +422,6 @@ int ArkodeMRISolver::init() { output.write("\tUsing difference quotient approximation for Jacobian in the inner solver\n"); } - if (treatment == MRI_Treatment::ImEx or treatment == MRI_Treatment::Implicit) { { output.write("\tUsing Newton iteration\n"); @@ -439,7 +438,7 @@ int ArkodeMRISolver::init() { /// Set Preconditioner if (use_precon) { - if (hasPreconditioner()) { + if (hasPreconditionerSlow()) { output.write("\tUsing user-supplied preconditioner\n"); if (ARKodeSetPreconditioner(arkode_mem, nullptr, arkode_s_pre) @@ -761,7 +760,7 @@ void ArkodeMRISolver::pre_s(BoutReal t, BoutReal gamma, BoutReal delta, BoutReal const BoutReal tstart = bout::globals::mpi->MPI_Wtime(); - if (!hasPreconditioner()) { + if (!hasPreconditionerSlow()) { // Identity (but should never happen) const auto length = N_VGetLocalLength_Parallel(uvec); std::copy(rvec, rvec + length, zvec); @@ -774,7 +773,7 @@ void ArkodeMRISolver::pre_s(BoutReal t, BoutReal gamma, BoutReal delta, BoutReal // Load vector to be inverted into F_vars load_derivs(rvec); - runPreconditioner(t, gamma, delta); + runPreconditionerSlow(t, gamma, delta); // Save the solution from F_vars save_derivs(zvec); @@ -789,7 +788,7 @@ void ArkodeMRISolver::pre_f(BoutReal t, BoutReal gamma, BoutReal delta, BoutReal const BoutReal tstart = bout::globals::mpi->MPI_Wtime(); - if (!hasPreconditioner()) { + if (!hasPreconditionerFast()) { // Identity (but should never happen) const auto length = N_VGetLocalLength_Parallel(uvec); std::copy(rvec, rvec + length, zvec); @@ -802,7 +801,7 @@ void ArkodeMRISolver::pre_f(BoutReal t, BoutReal gamma, BoutReal delta, BoutReal // Load vector to be inverted into F_vars load_derivs(rvec); - runPreconditioner(t, gamma, delta); + runPreconditionerFast(t, gamma, delta); // Save the solution from F_vars save_derivs(zvec); @@ -1001,4 +1000,4 @@ void ArkodeMRISolver::loop_abstol_values_op(Ind2D UNUSED(i2d), BoutReal* abstolv #else #endif // SUNDIALS_VERSION CHECK -#endif // BOUT_HAS_ARKODE \ No newline at end of file +#endif // BOUT_HAS_ARKODE diff --git a/src/solver/solver.cxx b/src/solver/solver.cxx index 3dd8e77894..39fcc493ca 100644 --- a/src/solver/solver.cxx +++ b/src/solver/solver.cxx @@ -1,7 +1,7 @@ /************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -1734,20 +1734,31 @@ void Solver::post_rhs(BoutReal UNUSED(t)) { #endif } -bool Solver::varAdded(const std::string& name) { - return contains(f2d, name) || contains(f3d, name) || contains(v2d, name) - || contains(v3d, name); +bool Solver::hasPreconditioner() { return model->hasPrecon(); } +bool Solver::hasPreconditionerFast() { return model->hasPreconFast(); } +bool Solver::hasPreconditionerSlow() { return model->hasPreconSlow(); } + +int Solver::runPreconditioner(BoutReal time, BoutReal gamma, BoutReal delta) { + return model->runPrecon(time, gamma, delta); } -bool Solver::hasPreconditioner() { return model->hasPrecon(); } +int Solver::runPreconditionerFast(BoutReal time, BoutReal gamma, BoutReal delta) { + return model->runPreconFast(time, gamma, delta); +} -int Solver::runPreconditioner(BoutReal t, BoutReal gamma, BoutReal delta) { - return model->runPrecon(t, gamma, delta); +int Solver::runPreconditionerSlow(BoutReal time, BoutReal gamma, BoutReal delta) { + return model->runPreconSlow(time, gamma, delta); } bool Solver::hasJacobian() { return model->hasJacobian(); } + int Solver::runJacobian(BoutReal time) { return model->runJacobian(time); } +bool Solver::varAdded(const std::string& name) { + return contains(f2d, name) || contains(f3d, name) || contains(v2d, name) + || contains(v3d, name); +} + // Add source terms to time derivatives void Solver::add_mms_sources(BoutReal t) { if (!mms) { From 7731d76b9decb57c8ae75faabbbdfcdffc648273 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 7 Feb 2025 18:04:47 -0800 Subject: [PATCH 22/28] ARKode solver: Use ARKodeGetNumRhsEvals from v7.2.0 Previous use of ARKStepGetNumRhsEvals is deprecated from SUNDIALS 7.2.0, replaced by ARKodeGetNumRhsEvals https://github.com/LLNL/sundials/blob/main/CHANGELOG.md#deprecation-notices --- src/solver/impls/arkode/arkode.cxx | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/solver/impls/arkode/arkode.cxx b/src/solver/impls/arkode/arkode.cxx index 23883cc043..289fa10284 100644 --- a/src/solver/impls/arkode/arkode.cxx +++ b/src/solver/impls/arkode/arkode.cxx @@ -517,7 +517,13 @@ int ArkodeSolver::run() { long int temp_long_int, temp_long_int2; ARKodeGetNumSteps(arkode_mem, &temp_long_int); nsteps = int(temp_long_int); +#if SUNDIALS_VERSION_AT_LEAST(7, 2, 0) + ARKodeGetNumRhsEvals(arkode_mem, 0, &temp_long_int); // Explicit + ARKodeGetNumRhsEvals(arkode_mem, 1, &temp_long_int2); // Implicit +#else + // This function was deprecated in 7.2.0 ARKStepGetNumRhsEvals(arkode_mem, &temp_long_int, &temp_long_int2); +#endif nfe_evals = int(temp_long_int); nfi_evals = int(temp_long_int2); if (treatment == Treatment::ImEx or treatment == Treatment::Implicit) { From 935131dc3e269e2af372427cd3d74c673fdf3df8 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 7 Feb 2025 18:07:10 -0800 Subject: [PATCH 23/28] ARKode-MRI: Only call ARKodeGetCurrentStep if user requests it Rather than calling `ARKodeGetCurrentStep` every rhs() call, only call when the `getCurrentTimestep()` method is called. --- src/solver/impls/arkode/arkode_mri.cxx | 26 ++++++-------------------- src/solver/impls/arkode/arkode_mri.hxx | 5 ++--- 2 files changed, 8 insertions(+), 23 deletions(-) diff --git a/src/solver/impls/arkode/arkode_mri.cxx b/src/solver/impls/arkode/arkode_mri.cxx index 71b3da2343..a612a4ac1d 100644 --- a/src/solver/impls/arkode/arkode_mri.cxx +++ b/src/solver/impls/arkode/arkode_mri.cxx @@ -166,6 +166,12 @@ ArkodeMRISolver::~ArkodeMRISolver() { MRIStepInnerStepper_Free(&inner_stepper); } +BoutReal ArkodeMRISolver::getCurrentTimestep() { + BoutReal hcur; + ARKodeGetCurrentStep(arkode_mem, &hcur); + return hcur; +} + /************************************************************************** * Initialise **************************************************************************/ @@ -646,10 +652,6 @@ void ArkodeMRISolver::rhs_se(BoutReal t, BoutReal* udata, BoutReal* dudata) { // Load state from udata load_vars(udata); - // Get the current timestep - // TO DO: Check to identify which time step is needed current/last, then correct accordingly - ARKodeGetCurrentStep(arkode_mem, &hcur); - // Call RHS function run_rhs_se(t); @@ -665,8 +667,6 @@ void ArkodeMRISolver::rhs_si(BoutReal t, BoutReal* udata, BoutReal* dudata) { TRACE("Running RHS: ArkodeMRISolver::rhs_si({:e})", t); load_vars(udata); - // TO DO: Check to identify which time step is needed current/last, then correct accordingly - ARKodeGetCurrentStep(arkode_mem, &hcur); // Call Implicit RHS function run_rhs_si(t); save_derivs(dudata); @@ -682,10 +682,6 @@ void ArkodeMRISolver::rhs_fe(BoutReal t, BoutReal* udata, BoutReal* dudata) { // Load state from udata load_vars(udata); - // Get the current timestep - // TO DO: Check to identify which time step is needed current/last, then correct accordingly - ARKodeGetCurrentStep(arkode_mem, &hcur); - // Call RHS function run_rhs_fe(t); @@ -701,8 +697,6 @@ void ArkodeMRISolver::rhs_fi(BoutReal t, BoutReal* udata, BoutReal* dudata) { TRACE("Running RHS: ArkodeMRISolver::rhs_si({:e})", t); load_vars(udata); - // TO DO: Check to identify which time step is needed current/last, then correct accordingly - ARKodeGetCurrentStep(arkode_mem, &hcur); // Call Implicit RHS function run_rhs_fi(t); save_derivs(dudata); @@ -718,10 +712,6 @@ void ArkodeMRISolver::rhs_s(BoutReal t, BoutReal* udata, BoutReal* dudata) { // Load state from udata load_vars(udata); - // Get the current timestep - // TO DO: Check to identify which time step is needed current/last, then correct accordingly - ARKodeGetCurrentStep(arkode_mem, &hcur); - // Call RHS function run_rhs_s(t); @@ -739,10 +729,6 @@ void ArkodeMRISolver::rhs_f(BoutReal t, BoutReal* udata, BoutReal* dudata) { // Load state from udata load_vars(udata); - // Get the current timestep - // Note: ARKodeGetCurrentStep updated too late in older versions - ARKodeGetCurrentStep(arkode_mem, &hcur); - // Call RHS function run_rhs_f(t); diff --git a/src/solver/impls/arkode/arkode_mri.hxx b/src/solver/impls/arkode/arkode_mri.hxx index d10794a120..0890c7e06f 100644 --- a/src/solver/impls/arkode/arkode_mri.hxx +++ b/src/solver/impls/arkode/arkode_mri.hxx @@ -76,7 +76,8 @@ public: explicit ArkodeMRISolver(Options* opts = nullptr); ~ArkodeMRISolver(); - BoutReal getCurrentTimestep() override { return hcur; } + // Get the current timestep + BoutReal getCurrentTimestep() override; int init() override; @@ -98,8 +99,6 @@ public: void jac_f(BoutReal t, BoutReal* ydata, BoutReal* vdata, BoutReal* Jvdata); private: - BoutReal hcur; //< Current internal timestep - bool diagnose{false}; //< Output additional diagnostics N_Vector uvec{nullptr}; //< Values From 702cefb044c591317c63f650419677f7bd6b2111 Mon Sep 17 00:00:00 2001 From: bendudson <219233+bendudson@users.noreply.github.com> Date: Tue, 11 Feb 2025 22:06:49 +0000 Subject: [PATCH 24/28] Apply black changes --- tests/integrated/test_suite | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/integrated/test_suite b/tests/integrated/test_suite index 307a8d84b3..77ad7882c4 100755 --- a/tests/integrated/test_suite +++ b/tests/integrated/test_suite @@ -188,7 +188,7 @@ class Test(threading.Thread): self.output += "\n(It is likely that a timeout occured)" else: # ❌ Failed - print("\u274C", end="") # No newline + print("\u274c", end="") # No newline print(" %7.3f s" % (time.time() - self.local.start_time), flush=True) def _cost(self): From 1d1677b3f9ffaad094fe9e2db891070c35e9d835 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Tue, 11 Feb 2025 16:36:49 -0800 Subject: [PATCH 25/28] Increment MRI rhs counts When run_rhs is called (e.g. by CVODE), all rhs functions are called (rhs_se, rhs_si, rhs_fe, rhs_fi), so all counters need to be incremented. --- src/solver/solver.cxx | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/solver/solver.cxx b/src/solver/solver.cxx index 39fcc493ca..4b6fb7b3b2 100644 --- a/src/solver/solver.cxx +++ b/src/solver/solver.cxx @@ -1566,8 +1566,13 @@ int Solver::run_rhs(BoutReal t, bool linear) { *t3 += *t2; } load_derivs(tmp3.begin()); // Put back time-derivatives - } - else if (model->splitOperator()) { + + rhs_ncalls_fe++; + rhs_ncalls_fi++; + rhs_ncalls_se++; + rhs_ncalls_si++; + + } else if (model->splitOperator()) { // Run both parts int nv = getLocalN(); From 3ccd334c4b9b68c57758f3bf066cda5ee50c1e1e Mon Sep 17 00:00:00 2001 From: maggul Date: Mon, 14 Apr 2025 02:06:10 -0500 Subject: [PATCH 26/28] run_rhs_f/s internally --- src/solver/impls/arkode/arkode_mri.cxx | 7 +-- src/solver/impls/arkode/arkode_mri.hxx | 3 - src/solver/solver.cxx | 46 +++++++++++++- .../integrated/test-kpr_mri/test_kpr_mri.cxx | 61 +++++-------------- 4 files changed, 62 insertions(+), 55 deletions(-) diff --git a/src/solver/impls/arkode/arkode_mri.cxx b/src/solver/impls/arkode/arkode_mri.cxx index a612a4ac1d..c4c17c54ec 100644 --- a/src/solver/impls/arkode/arkode_mri.cxx +++ b/src/solver/impls/arkode/arkode_mri.cxx @@ -107,9 +107,6 @@ ArkodeMRISolver::ArkodeMRISolver(Options* opts) .doc("Solve both fast and slow time scales using fixed time step sizes. NOTE: This is " "not recommended except for code comparison") .withDefault(false)), - inner_timestep((*options)["inner_timestep"] - .doc("Fixed step size to use for inner solver when running in fixed time step mode.") - .withDefault(1.0e-4)), order((*options)["order"].doc("Order of internal step").withDefault(3)), abstol((*options)["atol"].doc("Absolute tolerance").withDefault(1.0e-12)), reltol((*options)["rtol"].doc("Relative tolerance").withDefault(1.0e-10)), @@ -240,7 +237,7 @@ int ArkodeMRISolver::init() { if (fixed_step) { // If not given, default to adaptive timestepping - const BoutReal inner_fixed_timestep = (*options)["inner_timestep"]; + const BoutReal inner_fixed_timestep = (*options)["inner_timestep"].withDefault(1.0e-5); if (ARKodeSetFixedStep(inner_arkode_mem, inner_fixed_timestep) != ARK_SUCCESS) { throw BoutException("ARKodeSetFixedStep failed\n"); } @@ -293,7 +290,7 @@ int ArkodeMRISolver::init() { if (fixed_step) { // If not given, default to adaptive timestepping - const BoutReal fixed_timestep = (*options)["timestep"]; + const BoutReal fixed_timestep = (*options)["timestep"].withDefault(1.0e-4); if (ARKodeSetFixedStep(arkode_mem, fixed_timestep) != ARK_SUCCESS) { throw BoutException("ARKodeSetFixedStep failed\n"); } diff --git a/src/solver/impls/arkode/arkode_mri.hxx b/src/solver/impls/arkode/arkode_mri.hxx index 0890c7e06f..15e35d624e 100644 --- a/src/solver/impls/arkode/arkode_mri.hxx +++ b/src/solver/impls/arkode/arkode_mri.hxx @@ -122,9 +122,6 @@ private: /// Solve both fast and slow portion in fixed timestep mode. /// NOTE: This is not recommended except for code comparison bool fixed_step; - /// Fixed step size to use for inner solver when running in fixed - /// time step mode. - BoutReal inner_timestep; /// Order of the internal step int order; /// Absolute tolerance diff --git a/src/solver/solver.cxx b/src/solver/solver.cxx index 4b6fb7b3b2..ef4c4354c9 100644 --- a/src/solver/solver.cxx +++ b/src/solver/solver.cxx @@ -1479,8 +1479,29 @@ int Solver::run_rhs_s(BoutReal t, bool linear) { } pre_rhs(t); - status = model->runRHS_s(t, linear); + + // Run both parts + + int nv = getLocalN(); + // Create two temporary arrays for system state + Array tmp(nv); + Array tmp2(nv); + + save_vars(tmp.begin()); // Copy variables into tmp + pre_rhs(t); + status = model->runRHS_se(t, linear); + post_rhs(t); // Check variables, apply boundary conditions + + load_vars(tmp.begin()); // Reset variables + save_derivs(tmp.begin()); // Save time derivatives + pre_rhs(t); + status = model->runRHS_si(t, linear); post_rhs(t); + save_derivs(tmp2.begin()); // Save time derivatives + for (BoutReal *t = tmp.begin(), *t2 = tmp2.begin(); t != tmp.end(); ++t, ++t2) { + *t += *t2; + } + load_derivs(tmp.begin()); // Put back time-derivatives // If using Method of Manufactured Solutions add_mms_sources(t); @@ -1502,8 +1523,29 @@ int Solver::run_rhs_f(BoutReal t, bool linear) { } pre_rhs(t); - status = model->runRHS_f(t, linear); + + // Run both parts + + int nv = getLocalN(); + // Create two temporary arrays for system state + Array tmp(nv); + Array tmp2(nv); + + save_vars(tmp.begin()); // Copy variables into tmp + pre_rhs(t); + status = model->runRHS_fe(t, linear); + post_rhs(t); // Check variables, apply boundary conditions + + load_vars(tmp.begin()); // Reset variables + save_derivs(tmp.begin()); // Save time derivatives + pre_rhs(t); + status = model->runRHS_fi(t, linear); post_rhs(t); + save_derivs(tmp2.begin()); // Save time derivatives + for (BoutReal *t = tmp.begin(), *t2 = tmp2.begin(); t != tmp.end(); ++t, ++t2) { + *t += *t2; + } + load_derivs(tmp.begin()); // Put back time-derivatives // If using Method of Manufactured Solutions add_mms_sources(t); diff --git a/tests/integrated/test-kpr_mri/test_kpr_mri.cxx b/tests/integrated/test-kpr_mri/test_kpr_mri.cxx index 50b94ee55f..b736f5428b 100644 --- a/tests/integrated/test-kpr_mri/test_kpr_mri.cxx +++ b/tests/integrated/test-kpr_mri/test_kpr_mri.cxx @@ -74,35 +74,6 @@ class TestSolver : public PhysicsModel { return 0; } - int rhs_s(BoutReal t) override { - /* fill in the RHS function: - [G e]*[(-1+f^2-0.5*cos(t))/(2*f)] + [-0.5*sin(t)/(2*f)] - [0 0] [(-2+g^2-cos(w*t))/(2*g) ] [ 0 ] */ - BoutReal tmp1 = (-1.0 + f(1,1,0) * f(1,1,0) - 0.5*cos(t)) / (2.0 * f(1,1,0)); - BoutReal tmp2 = (-2.0 + g(1,1,0) * g(1,1,0) - cos(w*t)) / (2.0 * g(1,1,0)); - ddt(f) = G * tmp1 + e * tmp2 - 0.5*sin(t) / (2.0 * f(1,1,0)); - ddt(g) = 0.0; - - return 0; - } - - int rhs_f(BoutReal t) override { - /* fill in the RHS function: - [0 0]*[(-1+f^2-0.5*cos(t))/(2*f)] + [ 0 ] - [e -1] [(-2+g^2-cos(w*t))/(2*g) ] [-w*sin(w*t)/(2*sqrt(2+cos(w*t)))] */ - BoutReal tmp1 = (-1.0 + f(1,1,0) * f(1,1,0) - 0.5*cos(t)) / (2.0 * f(1,1,0)); - BoutReal tmp2 = (-2.0 + g(1,1,0) * g(1,1,0) - cos(w*t)) / (2.0 * g(1,1,0)); - ddt(f) = 0.0; - ddt(g) = e * tmp1 - tmp2 - w * sin(w*t) / (2.0 * sqrt(2.0 + cos(w * t))); - - return 0; - } - - bool check_solution(BoutReal atol, BoutReal t) { - // Return true if correct solution - return ((std::abs(sqrt(0.5*cos(t) + 1.0) - f(1,1,0)) < atol) and (std::abs(sqrt(cos(w*t) + 2.0) - g(1,1,0)) < atol)); - } - BoutReal compute_error(BoutReal t) { /* Compute the error with the true solution: @@ -117,9 +88,9 @@ class TestSolver : public PhysicsModel { }; int main(int argc, char** argv) { - // Absolute tolerance for difference between the actual value and the - // expected value - constexpr BoutReal tolerance = 1.e-5; + // Relative and Absolute tolerances + constexpr BoutReal rtol = 1.e-10; + constexpr BoutReal atol = 1.e-12; // Our own output to stdout, as main library will only be writing to log files Output output_test; @@ -143,17 +114,22 @@ int main(int argc, char** argv) { bout::globals::mesh = Mesh::create(); bout::globals::mesh->load(); + std::string sunsolver = "arkode_mri"; + int nout = 100; + BoutReal timestep = 0.05; + BoutReal finaltime = nout*timestep; + // Global options - root["nout"] = 100; - root["timestep"] = 0.05; - root["arkode_mri"]["treatment"] = "explicit"; - root["arkode_mri"]["inner_treatment"] = "explicit"; + root["nout"] = nout; + root["timestep"] = timestep; + root[sunsolver]["rtol"] = rtol; + root[sunsolver]["atol"] = atol; // Get specific options section for this solver. Can't just use default // "solver" section, as we run into problems when solvers use the same // name for an option with inconsistent defaults - auto options = Options::getRoot()->getSection("arkode_mri"); - auto solver = std::unique_ptr{Solver::create("arkode_mri", options)}; + auto options = Options::getRoot()->getSection(sunsolver); + auto solver = std::unique_ptr{Solver::create(sunsolver, options)}; TestSolver model{}; solver->setModel(&model); @@ -163,16 +139,11 @@ int main(int argc, char** argv) { solver->solve(); - BoutReal error = model.compute_error(5.0); + BoutReal error = model.compute_error(finaltime); std::cout << "error = " << error << std::endl; - if (model.check_solution(tolerance, 5.0)) { - output_test << " PASSED\n"; - return 0; - } - output_test << " FAILED\n"; - return 1; + return 0; } #else // ARKODE-MRI option for BOUT++ From 5c98b44b4bdaed9a7cdc06d659ab69ef08e0706a Mon Sep 17 00:00:00 2001 From: maggul Date: Mon, 21 Apr 2025 17:16:37 -0500 Subject: [PATCH 27/28] ARKodeSetStopTime before each evolve call --- src/solver/impls/arkode/arkode_mri.cxx | 28 +++++++++++++++++++++++++- src/solver/impls/arkode/arkode_mri.hxx | 5 +++++ 2 files changed, 32 insertions(+), 1 deletion(-) diff --git a/src/solver/impls/arkode/arkode_mri.cxx b/src/solver/impls/arkode/arkode_mri.cxx index c4c17c54ec..eb0124d735 100644 --- a/src/solver/impls/arkode/arkode_mri.cxx +++ b/src/solver/impls/arkode/arkode_mri.cxx @@ -87,6 +87,12 @@ ArkodeMRISolver::ArkodeMRISolver(Options* opts) mxsteps((*options)["mxstep"] .doc("Maximum number of steps to take between outputs") .withDefault(500)), + // mxstepsize((*options)["mxstepsize"] + // .doc("Maximum step size") + // .withDefault(0.1)), + // inner_mxstepsize((*options)["inner_mxstepsize"] + // .doc("Maximum step size") + // .withDefault(0.1)), treatment((*options)["treatment"] .doc("Use default capability (imex) or provide a specific treatment: " "implicit or explicit") @@ -353,6 +359,13 @@ int ArkodeMRISolver::init() { throw BoutException("ARKodeSetMaxNumSteps failed\n"); } + // if (ARKodeSetMaxStep(arkode_mem, mxstepsize) != ARK_SUCCESS) { + // throw BoutException("ARKodeSetMaxStep failed\n"); + // } + // if (ARKodeSetMaxStep(inner_arkode_mem, inner_mxstepsize) != ARK_SUCCESS) { + // throw BoutException("ARKodeSetMaxStep failed\n"); + // } + if (inner_treatment == MRI_Treatment::ImEx or inner_treatment == MRI_Treatment::Implicit) { { output.write("\tUsing Newton iteration for inner solver\n"); @@ -599,7 +612,20 @@ BoutReal ArkodeMRISolver::run(BoutReal tout) { pre_Wtime_s = 0.0; pre_ncalls_s = 0; - int flag; + int flag = ARKodeSetStopTime(arkode_mem, 1.0001*tout); + if (flag != ARK_SUCCESS) { + output_error.write("ERROR ARKodeSetStopTime failed at t = {:e}, flag = {:d}\n", + simtime, flag); + return -1.0; + } + + flag = ARKodeSetStopTime(inner_arkode_mem, 1.0001*tout); + if (flag != ARK_SUCCESS) { + output_error.write("ERROR ARKodeSetStopTime failed at t = {:e}, flag = {:d}\n", + simtime, flag); + return -1.0; + } + if (!monitor_timestep) { // Run in normal mode flag = ARKodeEvolve(arkode_mem, tout, uvec, &simtime, ARK_NORMAL); diff --git a/src/solver/impls/arkode/arkode_mri.hxx b/src/solver/impls/arkode/arkode_mri.hxx index 15e35d624e..ced2590f15 100644 --- a/src/solver/impls/arkode/arkode_mri.hxx +++ b/src/solver/impls/arkode/arkode_mri.hxx @@ -113,6 +113,11 @@ private: /// Maximum number of steps to take between outputs int mxsteps; + + // /// Maximum step sizes + // double mxstepsize; + // double inner_mxstepsize; + /// Integrator treatment enum: IMEX, Implicit or Explicit MRI_Treatment treatment; MRI_Treatment inner_treatment; From 1bd0b0b738154c97944c4b3e74e24016d1b8f3c1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mustafa=20A=C4=9Fg=C3=BCl?= <33010171+maggul@users.noreply.github.com> Date: Tue, 5 Aug 2025 16:47:07 -0500 Subject: [PATCH 28/28] SUNDIALS version update --- cmake/SetupBOUTThirdParty.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/SetupBOUTThirdParty.cmake b/cmake/SetupBOUTThirdParty.cmake index 10942f8aa9..f49131ddf1 100644 --- a/cmake/SetupBOUTThirdParty.cmake +++ b/cmake/SetupBOUTThirdParty.cmake @@ -278,7 +278,7 @@ if (BOUT_USE_SUNDIALS) FetchContent_Declare( sundials GIT_REPOSITORY https://github.com/LLNL/sundials - GIT_TAG v7.2.1 + GIT_TAG v7.4.0 ) # Note: These are settings for building SUNDIALS set(EXAMPLES_ENABLE_C OFF CACHE BOOL "" FORCE)