diff --git a/CMakeLists.txt b/CMakeLists.txt
index c45fca3b72..fa818f7cc0 100644
--- a/CMakeLists.txt
+++ b/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
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)
diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx
index a1c88a2634..7d963a291a 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 <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";
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..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 <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,6 +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 fast and slow parts
+  void setSplitOperatorMRI(bool split = true) { splitopmri = split; }
+
   /// Specify a preconditioner function
   void setPrecon(preconfunc pset) { userprecon = pset; }
   template <class Model>
@@ -316,6 +348,19 @@ protected:
     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
diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx
index 446acefecb..5b14fc09ed 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 = "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.
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 <bout/dcomplex.hxx>
 #include <bout/derivs.hxx>
 #include <bout/fft.hxx>
-#include <bout/griddata.hxx>
 #include <bout/msg_stack.hxx>
 #include <bout/options.hxx>
 #include <bout/output.hxx>
diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx
index cb6cb8410c..720c3ff963 100644
--- a/src/mesh/mesh.cxx
+++ b/src/mesh/mesh.cxx
@@ -1,15 +1,15 @@
-#include <bout/boutcomm.hxx>
 #include <bout/coordinates.hxx>
 #include <bout/derivs.hxx>
 #include <bout/globals.hxx>
-#include <bout/griddata.hxx>
 #include <bout/mesh.hxx>
 #include <bout/msg_stack.hxx>
-#include <bout/output.hxx>
 #include <bout/utils.hxx>
 
 #include <cmath>
 
+#include <bout/boutcomm.hxx>
+#include <bout/output.hxx>
+
 #include "impls/bout/boutmesh.hxx"
 
 MeshFactory::ReturnType MeshFactory::create(Options* options,
diff --git a/src/physics/physicsmodel.cxx b/src/physics/physicsmodel.cxx
index f4be3bde2a..8eb729faef 100644
--- a/src/physics/physicsmodel.cxx
+++ b/src/physics/physicsmodel.cxx
@@ -112,9 +112,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);
@@ -124,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) {
@@ -133,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.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) {
diff --git a/src/solver/impls/arkode/arkode_mri.cxx b/src/solver/impls/arkode/arkode_mri.cxx
new file mode 100644
index 0000000000..eb0124d735
--- /dev/null
+++ b/src/solver/impls/arkode/arkode_mri.cxx
@@ -0,0 +1,1012 @@
+/**************************************************************************
+ * Experimental interface to SUNDIALS ARKODE MRI solver
+ *
+ * NOTE: ARKODE is still in beta testing so use with cautious optimism
+ *
+ **************************************************************************
+ * Copyright 2010-2025 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 <http://www.gnu.org/licenses/>.
+ *
+ **************************************************************************/
+
+#include "bout/build_defines.hxx"
+
+#if BOUT_HAS_ARKODE
+
+#include "arkode_mri.hxx"
+
+#if SUNDIALS_VERSION_AT_LEAST(7, 2, 0)
+
+#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 <arkode/arkode.h>
+#include <arkode/arkode_arkstep.h>
+#include <arkode/arkode_bbdpre.h>
+#include <arkode/arkode_ls.h>
+
+#include <sundials/sundials_math.h>
+
+#include <algorithm>
+#include <iterator>
+#include <numeric>
+#include <vector>
+
+// 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)),
+      // 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")
+                    .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 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(3)),
+      abstol((*options)["atol"].doc("Absolute tolerance").withDefault(1.0e-12)),
+      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)),
+      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);
+}
+
+BoutReal ArkodeMRISolver::getCurrentTimestep() {
+  BoutReal hcur;
+  ARKodeGetCurrentStep(arkode_mem, &hcur);
+  return hcur;
+}
+
+/**************************************************************************
+ * 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"].withDefault(1.0e-5);
+    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 (ARKodeCreateMRIStepInnerStepper(inner_arkode_mem, &inner_stepper) != ARK_SUCCESS) {
+    throw BoutException("ARKodeCreateMRIStepInnerStepper 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-MRI 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-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-MRI 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"].withDefault(1.0e-4);
+    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");
+  }
+
+  if (use_vector_abstol) {
+    std::vector<BoutReal> f2dtols;
+    f2dtols.reserve(f2d.size());
+    std::transform(begin(f2d), end(f2d), std::back_inserter(f2dtols),
+                   [abstol = abstol](const VarStr<Field2D>& 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<BoutReal> f3dtols;
+    f3dtols.reserve(f3d.size());
+    std::transform(begin(f3d), end(f3d), std::back_inserter(f3dtols),
+                   [abstol = abstol](const VarStr<Field3D>& 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 (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");
+
+      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 (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 (fast) 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<Field3D>& fvar) {
+                Mesh* localmesh = fvar.var->getMesh();
+                return acc + localmesh->xend - localmesh->xstart + 3;
+              });
+
+          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 inner_mldq = (*options)["inner_mldq"]
+                                .doc("Lower half-bandwidth to be used in the difference "
+                                     "quotient Jacobian approximation")
+                                .withDefault(band_width_default);
+          const auto inner_mukeep = (*options)["inner_mukeep"]
+                                  .doc("Upper half-bandwidth of the retained banded "
+                                       "approximate Jacobian block")
+                                  .withDefault(n3Dvars() + n2Dvars());
+          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, inner_mudq, inner_mldq, inner_mukeep, inner_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 (hasPreconditionerSlow()) {
+          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<Field3D>& 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);
+    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) {
+      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);
+    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) {
+      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<BoutReal>(nniters) / static_cast<BoutReal>(nsteps));
+        output.write("    -> Linear iterations per Newton iteration: {:e}\n",
+                     static_cast<BoutReal>(nliters) / static_cast<BoutReal>(nniters));
+        output.write("    -> Preconditioner evaluations per Newton: {:e}\n",
+                     static_cast<BoutReal>(npevals) / static_cast<BoutReal>(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<BoutReal>(inner_nniters) / static_cast<BoutReal>(inner_nsteps));
+        output.write("    -> Inner Linear iterations per Newton iteration: {:e}\n",
+                     static_cast<BoutReal>(inner_nliters) / static_cast<BoutReal>(inner_nniters));
+        output.write("    -> Inner Preconditioner evaluations per Newton: {:e}\n",
+                     static_cast<BoutReal>(inner_npevals) / static_cast<BoutReal>(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 = 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);
+  } 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);
+    }
+    // Update the current simulation time
+    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);
+
+  // 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);
+  // 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);
+
+  // 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);
+  // 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);
+
+  // Call RHS function
+  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);
+
+  // Call RHS function
+  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 slow preconditioner: ArkodeMRISolver::pre_s({:e})", t);
+
+  const BoutReal tstart = bout::globals::mpi->MPI_Wtime();
+
+  if (!hasPreconditionerSlow()) {
+    // 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);
+
+  runPreconditionerSlow(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 fast preconditioner: ArkodeMRISolver::pre_f({:e})", t);
+
+  const BoutReal tstart = bout::globals::mpi->MPI_Wtime();
+
+  if (!hasPreconditionerFast()) {
+    // 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);
+
+  runPreconditionerFast(t, gamma, delta);
+
+  // Save the solution from F_vars
+  save_derivs(zvec);
+
+  pre_Wtime_s += bout::globals::mpi->MPI_Wtime() - tstart;
+  pre_ncalls_s++;
+}
+
+/**************************************************************************
+ * 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<ArkodeMRISolver*>(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<ArkodeMRISolver*>(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<ArkodeMRISolver*>(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<ArkodeMRISolver*>(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<ArkodeMRISolver*>(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<ArkodeMRISolver*>(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<ArkodeMRISolver*>(user_data);
+
+  // Run slow preconditioner
+  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<ArkodeMRISolver*>(user_data);
+
+  // Run fast preconditioner
+  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<BoutReal>& f2dtols,
+                                     std::vector<BoutReal>& 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<BoutReal>& f2dtols,
+                                         std::vector<BoutReal>& f3dtols, bool bndry) {
+  // Loop over 2D variables
+  for (std::vector<BoutReal>::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<BoutReal>::size_type i = 0; i < f3dtols.size(); i++) {
+      if (bndry && !f3d[i].evolve_bndry) {
+        continue;
+      }
+      abstolvec_data[p] = f3dtols[i];
+      p++;
+    }
+  }
+}
+
+#else
+#endif // SUNDIALS_VERSION CHECK
+#endif // BOUT_HAS_ARKODE
diff --git a/src/solver/impls/arkode/arkode_mri.hxx b/src/solver/impls/arkode/arkode_mri.hxx
new file mode 100644
index 0000000000..ced2590f15
--- /dev/null
+++ b/src/solver/impls/arkode/arkode_mri.hxx
@@ -0,0 +1,180 @@
+/**************************************************************************
+ * 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-2025 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 <http://www.gnu.org/licenses/>.
+ *
+ **************************************************************************/
+
+#ifndef BOUT_ARKODE_MRI_SOLVER_H
+#define BOUT_ARKODE_MRI_SOLVER_H
+
+#include "bout/build_defines.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/sundials_backports.hxx"
+
+#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 <nvector/nvector_parallel.h>
+#include <arkode/arkode_mristep.h>
+#include <sundials/sundials_config.h>
+#include <vector>
+
+class ArkodeMRISolver;
+class Options;
+
+namespace {
+RegisterSolver<ArkodeMRISolver> 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);
+
+class ArkodeMRISolver : public Solver {
+public:
+  explicit ArkodeMRISolver(Options* opts = nullptr);
+  ~ArkodeMRISolver();
+
+  // Get the current timestep
+  BoutReal getCurrentTimestep() override;
+
+  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:
+  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;
+
+  // /// Maximum step sizes
+  // double mxstepsize;
+  // double inner_mxstepsize;
+
+  /// 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 both fast and slow portion in fixed timestep mode.
+  /// NOTE: This is not recommended except for code comparison
+  bool fixed_step;
+  /// Order of the internal step
+  int order;
+  /// 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<BoutReal>& f2dtols,
+                         std::vector<BoutReal>& f3dtols);
+  void loop_abstol_values_op(Ind2D i2d, BoutReal* abstolvec_data, int& p,
+                             std::vector<BoutReal>& f2dtols,
+                             std::vector<BoutReal>& 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};
+  /// Context for SUNDIALS memory allocations
+  sundials::Context suncontext;
+};
+
+#endif // SUNDIALS_VERSION CHECK
+#endif // BOUT_HAS_ARKODE
+#endif // BOUT_ARKODE_MRI_SOLVER_H
+
diff --git a/src/solver/solver.cxx b/src/solver/solver.cxx
index a3f6a874f6..ef4c4354c9 100644
--- a/src/solver/solver.cxx
+++ b/src/solver/solver.cxx
@@ -1,8 +1,8 @@
 /**************************************************************************
- * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu
+ * Copyright 2010 - 2025 BOUT++ contributors
+ *
+ * Contact: Ben Dudson, dudson2@llnl.gov
  *
- * Contact: Ben Dudson, bd512@york.ac.uk
- * 
  * This file is part of BOUT++.
  *
  * BOUT++ is free software: you can redistribute it and/or modify
@@ -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"
@@ -949,7 +950,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(); }
 
 /////////////////////////////////////////////////////
 
@@ -1365,6 +1391,170 @@ Field3D Solver::globalIndex(int localStart) {
  * Running user-supplied functions
  **************************************************************************/
 
+int Solver::run_rhs_se(BoutReal t, bool linear) {
+  int status;
+
+  Timer timer("rhs");
+
+  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_s_call) {
+    // Ensure that nonlinear terms are calculated on first call
+    linear = false;
+    first_rhs_s_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");
+
+  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_f_call) {
+    // Ensure that nonlinear terms are calculated on first call
+    linear = false;
+    first_rhs_f_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);
+
+  // Run both parts
+
+  int nv = getLocalN();
+  // Create two temporary arrays for system state
+  Array<BoutReal> tmp(nv);
+  Array<BoutReal> 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);
+
+  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);
+
+  // Run both parts
+
+  int nv = getLocalN();
+  // Create two temporary arrays for system state
+  Array<BoutReal> tmp(nv);
+  Array<BoutReal> 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);
+
+  rhs_ncalls_fe++;
+  rhs_ncalls_fi++;
+  return status;
+}
+
 int Solver::run_rhs(BoutReal t, bool linear) {
   int status;
 
@@ -1376,7 +1566,55 @@ 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<BoutReal> tmp(nv);
+    Array<BoutReal> tmp2(nv);
+    Array<BoutReal> 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
+
+    rhs_ncalls_fe++;
+    rhs_ncalls_fi++;
+    rhs_ncalls_se++;
+    rhs_ncalls_si++;
+
+  } else if (model->splitOperator()) {
     // Run both parts
 
     int nv = getLocalN();
@@ -1399,18 +1637,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;
 }
 
@@ -1540,20 +1781,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) {
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/integrated/test-kpr_mri/CMakeLists.txt b/tests/integrated/test-kpr_mri/CMakeLists.txt
new file mode 100644
index 0000000000..45523803db
--- /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..2233782dd8
--- /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)
+   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
+ 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..b736f5428b
--- /dev/null
+++ b/tests/integrated/test-kpr_mri/test_kpr_mri.cxx
@@ -0,0 +1,155 @@
+#include "bout/sundials_backports.hxx"
+
+#if SUNDIALS_VERSION_AT_LEAST(7, 2, 0)
+
+#include "bout/physicsmodel.hxx"
+#include "bout/solver.hxx"
+
+#include <cmath>
+#include <memory>
+#include <string>
+#include <vector>
+
+// A simple phyics model with an analytical 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 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));
+    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_fi(BoutReal UNUSED(t)) override {
+
+    ddt(f) = 0.0;
+    ddt(g) = 0.0;
+
+    return 0;
+  }
+
+  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) {
+  // 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;
+
+  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();
+
+  std::string sunsolver = "arkode_mri";
+  int nout = 100;
+  BoutReal timestep = 0.05;
+  BoutReal finaltime = nout*timestep;
+
+  // Global options
+  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(sunsolver);
+  auto solver = std::unique_ptr<Solver>{Solver::create(sunsolver, options)};
+
+  TestSolver model{};
+  solver->setModel(&model);
+
+  BoutMonitor bout_monitor{};
+  solver->addMonitor(&bout_monitor, Solver::BACK);
+
+  solver->solve();
+
+  BoutReal error = model.compute_error(finaltime);
+
+  std::cout << "error = " << error << std::endl;
+
+  return 0;
+}
+#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
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):