ARKODE_MRI interface

CMakeLists.txt

@@ -290,6 +290,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

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)

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 <bout/deriv_store.hxx>
+#include <bout/index_derivs_interface.hxx>
+#include <bout/mpi_wrapper.hxx>
+
 #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 <bout/griddata.hxx>
+
+#include "bout/coordinates.hxx" // Coordinates class
+
 #include "bout/unused.hxx"
 
-#include "mpi.h"
+#include "bout/generic_factory.hxx"
+#include <bout/region.hxx>
+
+#include <bout/bout_enum_class.hxx>
 
+#include <list>
 #include <map>
 #include <memory>
 #include <optional>
 #include <set>
 #include <string>
 
-class BoundaryRegion;
-class BoundaryRegionPar;
-class GridDataSource;
-
 class MeshFactory : public Factory<Mesh, MeshFactory, GridDataSource*, Options*> {
 public:
   static constexpr auto type_name = "Mesh";

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

include/bout/physicsmodel.hxx

@@ -4,16 +4,10 @@
  * @brief Base class for Physics Models
  * 
  * 
- *
- * Changelog:
- * 
- * 2013-08 Ben Dudson <benjamin.dudson@york.ac.uk>
- *    * 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++.
  *
@@ -169,13 +163,26 @@ 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);
 
   /*!
    * True if this model uses split operators
+   * RHS = convective + diffusive
    */
   bool splitOperator();
 
+  /*!
+   * True if this model uses Multi-Rate Integrator (MRI) split operators
+   * RHS =  rhs_se + rhs_si + rhs_fe + rhs_fi
+   */
+  bool splitOperatorMRI();
+
   /*!
    * Run the convective (usually explicit) part of the model
    *
@@ -197,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 
@@ -208,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
    */
@@ -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,13 +338,29 @@ 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 fast and slow parts
+  void setSplitOperatorMRI(bool split = true) { splitopmri = split; }
+
   /// Specify a preconditioner function
   void setPrecon(preconfunc pset) { userprecon = pset; }
   template <class Model>
   void setPrecon(ModelPreconFunc<Model> preconditioner) {
     userprecon = static_cast<preconfunc>(preconditioner);
   }
 
+  /// Preconditioner for fast implicit RHS
+  void setPreconFast(preconfunc pset) { userprecon_f = pset; }
+  template <class Model>
+  void setPreconFast(ModelPreconFunc<Model> preconditioner) {
+    userprecon_f = static_cast<preconfunc>(preconditioner);
+  }
+  /// Preconditioner for slow implicit RHS
+  void setPreconSlow(preconfunc pset) { userprecon_s = pset; }
+  template <class Model>
+  void setPreconSlow(ModelPreconFunc<Model> preconditioner) {
+    userprecon_s = static_cast<preconfunc>(preconditioner);
+  }
+
   /// Specify a Jacobian-vector multiply function
   void setJacobian(jacobianfunc jset) { userjacobian = jset; }
   template <class Model>
@@ -391,8 +436,12 @@ private:
   bool restart_enabled{true};
   /// Split operator model?
   bool splitop{false};
+  /// 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

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 = "arkode_mri";
 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
@@ -433,11 +446,19 @@ 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};
 
   /// 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);
@@ -472,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();
@@ -542,6 +568,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.