diff --git a/Docs/sphinx_documentation/source/TimeIntegration_Chapter.rst b/Docs/sphinx_documentation/source/TimeIntegration_Chapter.rst index 720c312e2d..ef311de043 100644 --- a/Docs/sphinx_documentation/source/TimeIntegration_Chapter.rst +++ b/Docs/sphinx_documentation/source/TimeIntegration_Chapter.rst @@ -16,13 +16,8 @@ A Simple Time Integrator Setup ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ This is best shown with some sample code that sets up a time integrator and -asks it to step forwards by some interval ``dt``. The user needs to supply at -minimum, the right-hand side function using the ``TimeIntegrator::set_rhs()`` -function. By using the ``TimeIntegrator::set_post_update()`` function, a user -can also supply a post update function which is called on state data immediately -before evaluating the right-hand side. This post update function is a good -opportunity to fill boundary conditions for Runge-Kutta stage solution data so that -ghost cells are filled when the right hand side function is called on that solution data. +asks it to step forward by some interval ``dt``. The user needs to supply the +right-hand side function using the ``TimeIntegrator::set_rhs()`` function. .. highlight:: c++ @@ -30,171 +25,26 @@ ghost cells are filled when the right hand side function is called on that solut #include - MultiFab Sborder; // MultiFab containing old-time state data and ghost cells - MultiFab Snew; // MultiFab where we want new-time state data - Geometry geom; // The domain (or level) geometry + MultiFab Sold; // MultiFab containing old-time state data + MultiFab Snew; // MultiFab where we want new-time state data - // [Fill Sborder here] + // [Fill Sold here] // Create a time integrator that will work with // MultiFabs with the same BoxArray, DistributionMapping, // and number of components as the state_data MultiFab. - TimeIntegrator integrator(Sborder); + TimeIntegrator integrator(Sold); - // Create a RHS source function we will integrate - auto source_fun = [&](MultiFab& rhs, const MultiFab& state, const Real time){ - // User function to calculate the rhs MultiFab given the state MultiFab - fill_rhs(rhs, state, time); + // Create a function that fills the state BCs and computes the RHS + auto rhs_fun = [&](MultiFab& rhs, MultiFab& state, const Real time){ + // [Calculate the rhs MultiFab given the state MultiFab] }; - // Create a function to call after updating a state - auto post_update_fun = [&](MultiFab& S_data, const Real time) { - // Call user function to update state MultiFab, e.g. fill BCs - post_update(S_data, time, geom); - }; - - // Attach the right hand side and post-update functions - // to the integrator + // Attach the right hand side function to the integrator integrator.set_rhs(source_fun); - integrator.set_post_update(post_update_fun); - - // integrate forward one step from `time` by `dt` to fill S_new - integrator.advance(Sborder, S_new, time, dt); - -.. _sec:time_int:sundials: - -Using SUNDIALS -^^^^^^^^^^^^^^ - -The AMReX Time Integration interface also supports a SUNDIALS backend that -wraps both the explicit Runge-Kutta (ERK) and multirate (MRI) integration -schemes in the SUNDIALS ARKODE package. To use either of them, the user needs -to compile AMReX with `USE_SUNDIALS=TRUE` and use SUNDIALS v. 6.0 or later. - -There are only minor changes to the code above required to use the SUNDIALS -interface. The first change is that the integration datatype is now a -`Vector` type instead of simply `MultiFab`. The reason for -introducing a `Vector` in this case, is to permit integrating state -data with different spatial centering (e.g. cell centered, face centered, node -centered) concurrently. Shown here is sample code equivalent to the code above, -suitable for the SUNDIALS explicit Runge-Kutta integrator: - -.. highlight:: c++ - -:: - - #include - - Vector Sborder; // MultiFab(s) containing old-time state data and ghost cells - Vector Snew; // MultiFab(s) where we want new-time state data - Geometry geom; // The domain (or level) geometry - - // [Fill Sborder here] - - // Create a time integrator that will work with - // MultiFabs with the same BoxArray, DistributionMapping, - // and number of components as the state_data MultiFab. - TimeIntegrator > integrator(Sborder); - - // Create a RHS source function we will integrate - auto source_fun = [&](Vector& rhs, const Vector& state, const Real time){ - // User function to calculate the rhs MultiFab given the state MultiFab - fill_rhs(rhs, state, time); - }; - - // Create a function to call after updating a state - auto post_update_fun = [&](Vector& S_data, const Real time) { - // Call user function to update state MultiFab, e.g. fill BCs - post_update(S_data, time, geom); - }; - - // Attach the right hand side and post-update functions - // to the integrator - integrator.set_rhs(source_fun); - integrator.set_post_update(post_update_fun); - - // integrate forward one step from `time` by `dt` to fill S_new - integrator.advance(Sborder, S_new, time, dt); - -Afterwards, to select the ERK integrator, one needs only to add the following -two input parameters at runtime: - -:: - - integration.type = SUNDIALS - integration.sundials.strategy = ERK - -If instead one wishes to use the SUNDIALS multirate integrator, then the user -will need to use the following runtime inputs parameters: - -:: - - integration.type = SUNDIALS - integration.sundials.strategy = MRI - -In addition, to set up the multirate problem, the user needs to supply a fast -timescale right-hand-side function in addition to the usual right hand side -function (which is interpreted as the slow timescale right-hand side). The user -will also need to supply the ratio of the slow timestep size to the fast -timestep size, which is an integer corresponding to the number of fast -timesteps the integrator will take per every slow timestep. An example code -snippet would look as follows: - -.. highlight:: c++ - -:: - - #include - - Vector Sborder; // Vector of MultiFab(s) containing old-time state data and ghost cells - Vector Snew; // Vector of MultiFab(s) where we want new-time state data - Geometry geom; // The domain (or level) geometry - - // [Fill Sborder here] - - // Create a time integrator that will work with - // MultiFabs with the same BoxArray, DistributionMapping, - // and number of components as the state_data MultiFab. - TimeIntegrator > integrator(Sborder); - - // Create a slow timescale RHS function we will integrate - auto rhs_fun = [&](Vector& rhs, const Vector& state, const Real time){ - // User function to calculate the rhs MultiFab given the state MultiFab(s) - fill_rhs(rhs, state, time); - }; - - // Create a fast timescale RHS function to integrate - auto rhs_fun_fast = [&](Vector& rhs, - const Vector& stage_data, - const Vector& state, const Real time) { - // User function to calculate the fast-timescale rhs MultiFab given - // the state MultiFab and stage_data which holds the previously - // accessed slow-timescale stage state data. - fill_fast_rhs(rhs, stage_data, state, time); - }; - - // The post update function is called after updating state data or - // immediately before using state data to calculate a fast or slow right hand side. - // (it is a good place to e.g. fill boundary conditions) - auto post_update_fun = [&](Vector& S_data, const Real time) { - // Call user function to update state MultiFab(s), e.g. fill BCs - post_update(S_data, time, geom); - }; - - // Attach the slow and fast right hand side functions to integrator - integrator.set_rhs(rhs_fun); - integrator.set_fast_rhs(rhs_fun_fast); - - // This sets the ratio of slow timestep size to fast timestep size as an integer, - // or equivalently, the number of fast timesteps per slow timestep. - integrator.set_slow_fast_timestep_ratio(2); - - // Attach the post update function to the integrator - integrator.set_post_update(post_update_fun); - - // integrate forward one step from `time` by `dt` to fill S_new - integrator.advance(Sborder, S_new, time, dt); + // integrate forward one step from `time` by `dt` to fill Snew + integrator.advance(Sold, Snew, time, dt); Picking A Time Integration Method ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -205,12 +55,6 @@ a generic explicit Runge-Kutta method. If Runge-Kutta is selected, then the user can choose which of a set of predefined Butcher Tables to use, or can choose to use a custom table and supply it manually. -When AMReX is compiled with SUNDIALS v.6 or later, the user also has an option -to use the SUNDIALS ARKODE integrator as a backend for the AMReX Time Integrator -class. The features of this interface evolve with the needs of our codes, so -they may not yet support all SUNDIALS configurations available. If you find you -need SUNDIALS options we have not implemented, please let us know. - The full set of integrator options are detailed as follows: :: @@ -222,7 +66,7 @@ The full set of integrator options are detailed as follows: ## (without the quotation marks) ## "ForwardEuler" or "0" = Native Forward Euler Integrator ## "RungeKutta" or "1" = Native Explicit Runge Kutta - ## "SUNDIALS" or "2" = SUNDIALS ARKODE Integrator + ## "SUNDIALS" or "2" = SUNDIALS Integrator ## for example: integration.type = RungeKutta @@ -246,35 +90,104 @@ The full set of integrator options are detailed as follows: integration.rk.nodes = 0 integration.rk.tableau = 0.0 - ## *** Parameters Needed For SUNDIALS ARKODE Integrator *** - ## integration.sundials.strategy specifies which ARKODE strategy to use. - ## The available options are (without the quotations): - ## "ERK" = Explicit Runge Kutta - ## "MRI" = Multirate Integrator - ## "MRITEST" = Tests the Multirate Integrator by setting a zero-valued fast RHS function - ## for example: - integration.sundials.strategy = ERK - - ## *** Parameters Specific to SUNDIALS ERK Strategy *** - ## (Requires integration.type=SUNDIALS and integration.sundials.strategy=ERK) - ## integration.sundials.erk.method specifies which explicit Runge Kutta method - ## for SUNDIALS to use. The following options are supported: - ## "SSPRK3" = 3rd order strong stability preserving RK (default) - ## "Trapezoid" = 2nd order trapezoidal rule - ## "ForwardEuler" = 1st order forward euler - ## for example: - integration.sundials.erk.method = SSPRK3 - - ## *** Parameters Specific to SUNDIALS MRI Strategy *** - ## (Requires integration.type=SUNDIALS and integration.sundials.strategy=MRI) - ## integration.sundials.mri.implicit_inner specifies whether or not to use an implicit inner solve - ## integration.sundials.mri.outer_method specifies which outer (slow) method to use - ## integration.sundials.mri.inner_method specifies which inner (fast) method to use - ## The following options are supported for both the inner and outer methods: - ## "KnothWolke3" = 3rd order Knoth-Wolke method (default for outer method) - ## "Trapezoid" = 2nd order trapezoidal rule - ## "ForwardEuler" = 1st order forward euler (default for inner method) - ## for example: - integration.sundials.mri.implicit_inner = false - integration.sundials.mri.outer_method = KnothWolke3 - integration.sundials.mri.inner_method = Trapezoid +.. _sec:time_int:sundials: + +Using SUNDIALS +^^^^^^^^^^^^^^ + +The AMReX Time Integration interface also supports a SUNDIALS backend that +provides explicit, implicit, and implicit-explicit (ImEx) Runge-Kutta methods +as well a multirate (MRI) methods from the ARKODE package in SUNDIALS. +To use SUNDIALS integrators, the user needs to compile AMReX with +``USE_SUNDIALS=TRUE`` and use SUNDIALS v6.0 or later. + +The SUNDIALS interface supports ``MultiFab`` or ``Vector`` data +types. Using a ``Vector`` permits integrating state data with +different spatial centering (e.g. cell centered, face centered, node centered) +concurrently. + +The same code as above can be used with SUNDIALS explicit or implicit +Runge-Kutta methods without any modification. To select a SUNDIALS explicit +Runge-Kutta integrator, one needs only to add the following two input parameters +at runtime: + +:: + + integration.type = SUNDIALS + integration.sundials.type = ERK + +One can select a different method type by changing the value of +``integration.sundials.type`` to one of the following values: + ++------------------------+--------------------------------------------------+ +| Input Option | SUNDIALS Method Type | ++========================+==================================================+ +| ERK | Explicit Runge-Kutta method | ++------------------------+--------------------------------------------------+ +| DIRK | Diagonally Implicit Runge-Kutta method | ++------------------------+--------------------------------------------------+ +| IMEX-RK | Implicit-Explicit Additive Runge-Kutta method | ++------------------------+--------------------------------------------------+ +| EX-MRI | Explicit Multirate Infinitesimal method | ++------------------------+--------------------------------------------------+ +| IM-MRI | Implicit Multirate Infinitesimal method | ++------------------------+--------------------------------------------------+ +| IMEX-MRI | Implicit-Explicit Multirate Infinitesimal method | ++------------------------+--------------------------------------------------+ + +For ImEx methods, the user needs to supply two right-hand side functions, an +implicit and an explicit function, using the function +``TimeIntegrator::set_imex_rhs()``. Similarly for multirate methods, the user +needs to supply slow and fast right-hand side functions using +``TimeIntegrator::set_rhs()`` to set the slow function and +``TimeIntegrator::set_fast_rhs()`` to set the fast function. With multirate +methods, one also needs to select the fast time scale method type using the +input option ``integration.sundials.fast_type`` which maybe set to ``ERK`` or +``DIRK``. + +To select a specific SUNDIALS method use the input option +``integration.sundials.method`` for ERK and DIRK methods as well as the slow +time scale method with an MRI integrator, use ``integration.sundials.method_i`` +and ``integration.sundials.method_e`` to set the implicit and explicit method in +an ImEx method, and ``integration.sundials.fast_method`` to set the ERK or DIRK +method used at the fast time scale with an MRI integrator. These options may be +set to any valid SUNDIALS method name, see the following sections in the +SUNDIALS documentation for more details: + +* `ERK methods `_ +* `DIRK methods `_ +* `ImEx methods `_ +* `MRI methods `_ + +The full set of integrator options are detailed as follows: + +:: + + # INTEGRATION WITH SUNDIALS + + # *** Select the SUNDIALS integrator backend *** + integration.type = SUNDIALS + + # *** Select the SUNDIALS method type *** + # ERK = Explicit Runge-Kutta method + # DIRK = Diagonally Implicit Runge-Kutta method + # IMEX-RK = Implicit-Explicit Additive Runge-Kutta method + # EX-MRI = Explicit Multirate Infatesimal method + # IM-MRI = Implicit Multirate Infatesimal method + # IMEX-MRI = Implicit-Explicit Multirate Infatesimal method + integration.sundials.type = ERK + + # *** Select a specific SUNDIALS ERK method *** + integration.sundials.method = ARKODE_BOGACKI_SHAMPINE_4_2_3 + + # *** Select a specific SUNDIALS ImEx method *** + integration.sundials.method_i = ARKODE_ARK2_DIRK_3_1_2 + integration.sundials.method_e = ARKODE_ARK2_ERK_3_1_2 + + # *** Select a specific SUNDIALS MRI method *** + integration.sundials.method = ARKODE_MIS_KW3 + integration.sundials.fast_method = ARKODE_KNOTH_WOLKE_3_3 + +The features of this interface evolve with the needs of our codes, so they may +not yet support all SUNDIALS configurations available. If you find you need +SUNDIALS options we have not implemented, please let us know. diff --git a/Src/Base/AMReX_FEIntegrator.H b/Src/Base/AMReX_FEIntegrator.H index f8a002ef53..0e89ffb9d5 100644 --- a/Src/Base/AMReX_FEIntegrator.H +++ b/Src/Base/AMReX_FEIntegrator.H @@ -53,9 +53,6 @@ public: // So we initialize S_new by copying the old state. IntegratorOps::Copy(S_new, S_old); - // Call the pre RHS hook - BaseT::pre_rhs_action(S_new, time); - // F = RHS(S, t) T& F = *F_nodes[0]; BaseT::Rhs(F, S_new, time); diff --git a/Src/Base/AMReX_IntegratorBase.H b/Src/Base/AMReX_IntegratorBase.H index d9af8053d7..af29dc82cc 100644 --- a/Src/Base/AMReX_IntegratorBase.H +++ b/Src/Base/AMReX_IntegratorBase.H @@ -165,31 +165,25 @@ protected: /** * \brief Rhs is the right-hand-side function the integrator will use. */ - std::function Rhs; + std::function Rhs; /** * \brief RhsIm is the implicit right-hand-side function an ImEx integrator * will use. */ - std::function RhsIm; + std::function RhsIm; /** * \brief RhsEx is the explicit right-hand-side function an ImEx integrator * will use. */ - std::function RhsEx; + std::function RhsEx; /** * \brief RhsFast is the fast timescale right-hand-side function a multirate * integrator will use. */ - std::function RhsFast; - - /** - * \brief The pre_rhs_action function is called by the integrator on state - * data before using it to evaluate a right-hand side. - */ - std::function pre_rhs_action; + std::function RhsFast; /** * \brief The post_stage_action function is called by the integrator on @@ -283,28 +277,23 @@ public: virtual ~IntegratorBase () = default; - void set_rhs (std::function F) + void set_rhs (std::function F) { Rhs = F; } - void set_imex_rhs (std::function Fi, - std::function Fe) + void set_imex_rhs (std::function Fi, + std::function Fe) { RhsIm = Fi; RhsEx = Fe; } - void set_fast_rhs (std::function F) + void set_fast_rhs (std::function F) { RhsFast = F; } - void set_pre_rhs_action (std::function A) - { - pre_rhs_action = A; - } - void set_post_stage_action (std::function A) { post_stage_action = A; @@ -325,12 +314,6 @@ public: post_fast_step_action = A; } - void set_post_update (std::function A) - { - set_post_stage_action(A); - set_post_step_action(A); - } - amrex::Real get_time_step () { return time_step; diff --git a/Src/Base/AMReX_RKIntegrator.H b/Src/Base/AMReX_RKIntegrator.H index f72890c808..a6efd02853 100644 --- a/Src/Base/AMReX_RKIntegrator.H +++ b/Src/Base/AMReX_RKIntegrator.H @@ -217,9 +217,6 @@ public: BaseT::post_stage_action(S_new, stage_time); } - // Call the update hook for the stage state value - BaseT::pre_rhs_action(S_new, stage_time); - // Fill F[i], the RHS at the current stage // F[i] = RHS(y, t) at y = stage_value, t = stage_time BaseT::Rhs(*F_nodes[i], S_new, stage_time); diff --git a/Src/Base/AMReX_TimeIntegrator.H b/Src/Base/AMReX_TimeIntegrator.H index 1044336153..316ad0ff31 100644 --- a/Src/Base/AMReX_TimeIntegrator.H +++ b/Src/Base/AMReX_TimeIntegrator.H @@ -65,17 +65,14 @@ private: void set_default_functions () { // By default, do nothing in the RHS - set_rhs([](T& /* S_rhs */, const T& /* S_data */, const amrex::Real /* time */){}); - set_imex_rhs([](T& /* S_rhs */, const T& /* S_data */, const amrex::Real /* time */){}, - [](T& /* S_rhs */, const T& /* S_data */, const amrex::Real /* time */){}); - set_fast_rhs([](T& /* S_rhs */, const T& /* S_data */, const amrex::Real /* time */){}); + set_rhs([](T& /* S_rhs */, T& /* S_data */, const amrex::Real /* time */){}); + set_imex_rhs([](T& /* S_rhs */, T& /* S_data */, const amrex::Real /* time */){}, + [](T& /* S_rhs */, T& /* S_data */, const amrex::Real /* time */){}); + set_fast_rhs([](T& /* S_rhs */, T& /* S_data */, const amrex::Real /* time */){}); // In general, the following functions can be used to fill BCs. Which // function to set will depend on the method type and intended use case - // By default, do nothing before calling the RHS - set_pre_rhs_action([](T& /* S_data */, amrex::Real /* time */){}); - // By default, do nothing after a stage or step set_post_stage_action([](T& /* S_data */, const amrex::Real /* time */){}); set_post_step_action([](T& /* S_data */, const amrex::Real /* time */){}); @@ -134,27 +131,22 @@ public: } } - void set_rhs (std::function F) + void set_rhs (std::function F) { integrator_ptr->set_rhs(F); } - void set_imex_rhs (std::function Fi, - std::function Fe) + void set_imex_rhs (std::function Fi, + std::function Fe) { integrator_ptr->set_imex_rhs(Fi, Fe); } - void set_fast_rhs (std::function F) + void set_fast_rhs (std::function F) { integrator_ptr->set_fast_rhs(F); } - void set_pre_rhs_action (std::function A) - { - integrator_ptr->set_pre_rhs_action(A); - } - void set_post_stage_action (std::function A) { integrator_ptr->set_post_stage_action(A); @@ -175,11 +167,6 @@ public: integrator_ptr->set_post_fast_step_action(A); } - void set_post_update (std::function A) - { - integrator_ptr->set_post_update(A); - } - amrex::Real get_time_step () { return integrator_ptr->get_time_step(); diff --git a/Src/Extern/SUNDIALS/AMReX_SundialsIntegrator.H b/Src/Extern/SUNDIALS/AMReX_SundialsIntegrator.H index 30ff30a499..47a028f852 100644 --- a/Src/Extern/SUNDIALS/AMReX_SundialsIntegrator.H +++ b/Src/Extern/SUNDIALS/AMReX_SundialsIntegrator.H @@ -428,7 +428,6 @@ public: T S_rhs; unpack_vector(y_rhs, S_rhs); - BaseT::pre_rhs_action(S_data, rhs_time); BaseT::Rhs(S_rhs, S_data, rhs_time); return 0; @@ -443,7 +442,6 @@ public: T S_rhs; unpack_vector(y_rhs, S_rhs); - BaseT::pre_rhs_action(S_data, rhs_time); BaseT::RhsIm(S_rhs, S_data, rhs_time); return 0; @@ -458,7 +456,6 @@ public: T S_rhs; unpack_vector(y_rhs, S_rhs); - BaseT::pre_rhs_action(S_data, rhs_time); BaseT::RhsEx(S_rhs, S_data, rhs_time); return 0; @@ -473,7 +470,6 @@ public: T S_rhs; unpack_vector(y_rhs, S_rhs); - BaseT::pre_rhs_action(S_data, rhs_time); BaseT::RhsFast(S_rhs, S_data, rhs_time); return 0;