Multiple time integrators in an input

framework/doc/content/syntax/Executioner/TimeIntegrators/index.md

@@ -0,0 +1,18 @@
+# Multiple Time Integrators
+
+The time integrator system is described in [TimeIntegrator/index.md]. Multiple
+time integrators in a single input are supported using
+`Executioner/TimeIntegrators` syntax as illustrated below:
+
+!listing test/tests/time_integrators/multiple-integrators/test.i block=Executioner
+
+where each individual time integrator has specific `variables` assigned to
+it.
+
+!syntax list /Executioner/TimeIntegrators objects=True actions=False subsystems=False
+
+!syntax list /Executioner/TimeIntegrators objects=False actions=False subsystems=True
+
+!syntax list /Executioner/TimeIntegrators objects=False actions=True subsystems=False
+
+!syntax parameters /Executioner/TimeIntegrators

framework/include/executioners/Executioner.h

@@ -108,7 +108,7 @@ class Executioner : public MooseObject,
    * This is an empty string for non-Transient executioners
    * @return A string of giving the TimeIntegrator name
    */
-  virtual std::string getTimeIntegratorName() const { return std::string(); }
+  virtual std::vector<std::string> getTimeIntegratorNames() const { return {}; }
 
   /**
    * Can be used by subclasses to call parentOutputPositionChanged()

framework/include/executioners/Transient.h

@@ -138,7 +138,7 @@ class Transient : public Executioner
    * Get the name of the time integrator (time integration scheme) used
    * @return string with the time integration scheme name
    */
-  virtual std::string getTimeIntegratorName() const override;
+  virtual std::vector<std::string> getTimeIntegratorNames() const override;
 
   /**
    * Get the time scheme used

framework/include/outputs/JsonIO.h

@@ -26,6 +26,8 @@ template <typename T>
 class DenseVector;
 template <typename T>
 class DenseMatrix;
+template <typename T>
+class NumericVector;
 }
 
 // Overloads for to_json, which _must_ be overloaded in the namespace
@@ -38,6 +40,7 @@ namespace libMesh
 void to_json(nlohmann::json & json, const Point & p);
 void to_json(nlohmann::json & json, const DenseVector<Real> & vector);
 void to_json(nlohmann::json & json, const DenseMatrix<Real> & matrix);
+void to_json(nlohmann::json & json, const std::unique_ptr<NumericVector<Number>> & vector);
 }
 
 namespace nlohmann

framework/include/restart/DataIO.h

@@ -35,6 +35,7 @@
 #include <iostream>
 #include <map>
 #include <unordered_map>
+#include <unordered_set>
 #include <memory>
 #include <optional>
 
@@ -341,6 +342,18 @@ dataStore(std::ostream & stream, std::unordered_map<T, U> & m, void * context)
   }
 }
 
+template <typename T>
+inline void
+dataStore(std::ostream & stream, std::unordered_set<T> & s, void * context)
+{
+  // First store the size of the set
+  std::size_t size = s.size();
+  dataStore(stream, size, nullptr);
+
+  for (auto & element : s)
+    dataStore(stream, element, context);
+}
+
 template <typename T>
 inline void
 dataStore(std::ostream & stream, std::optional<T> & m, void * context)
@@ -406,6 +419,7 @@ template <>
 void dataStore(std::ostream & stream, RealEigenMatrix & v, void * context);
 template <>
 void dataStore(std::ostream & stream, libMesh::Parameters & p, void * context);
+
 template <>
 /**
  * Stores an owned numeric vector.
@@ -666,6 +680,25 @@ dataLoad(std::istream & stream, std::unordered_map<T, U> & m, void * context)
   }
 }
 
+template <typename T>
+inline void
+dataLoad(std::istream & stream, std::unordered_set<T> & s, void * context)
+{
+  s.clear();
+
+  // First read the size of the set
+  std::size_t size = 0;
+  dataLoad(stream, size, nullptr);
+  s.reserve(size);
+
+  for (std::size_t i = 0; i < size; i++)
+  {
+    T element;
+    dataLoad(stream, element, context);
+    s.insert(element);
+  }
+}
+
 template <typename T>
 inline void
 dataLoad(std::istream & stream, std::optional<T> & m, void * context)

framework/include/systems/AuxiliarySystem.h

@@ -55,16 +55,6 @@ class AuxiliarySystem : public SystemBase, public PerfGraphInterface
   virtual void addVariable(const std::string & var_type,
                            const std::string & name,
                            InputParameters & parameters) override;
-  /**
-   * Add a time integrator
-   * @param type Type of the integrator
-   * @param name The name of the integrator
-   * @param parameters Integrator params
-   */
-  void addTimeIntegrator(const std::string & type,
-                         const std::string & name,
-                         InputParameters & parameters) override;
-  using SystemBase::addTimeIntegrator;
 
   /**
    * Adds an auxiliary kernel

framework/include/systems/DisplacedSystem.h

@@ -149,9 +149,9 @@ class DisplacedSystem : public SystemBase
     return _undisplaced_system.solutionUDotDotOld();
   }
 
-  virtual Number & duDotDu() override { return _undisplaced_system.duDotDu(); }
+  virtual std::vector<Number> & duDotDus() override { return _undisplaced_system.duDotDus(); }
   virtual Number & duDotDotDu() override { return _undisplaced_system.duDotDotDu(); }
-  virtual const Number & duDotDu() const override { return _undisplaced_system.duDotDu(); }
+  virtual const Number & duDotDu(unsigned int var_num = 0) const override;
   virtual const Number & duDotDotDu() const override { return _undisplaced_system.duDotDotDu(); }
 
   virtual void addDotVectors() override { _undisplaced_system.addDotVectors(); }
@@ -249,9 +249,6 @@ class DisplacedSystem : public SystemBase
   virtual System & system() override;
   virtual const System & system() const override;
 
-  using SystemBase::addTimeIntegrator;
-  void addTimeIntegrator(std::shared_ptr<TimeIntegrator> ti) override;
-
   virtual void compute(ExecFlagType) override {}
 
 protected:
@@ -297,3 +294,9 @@ DisplacedSystem::hasSolutionState(const unsigned int state,
 {
   return _undisplaced_system.hasSolutionState(state, iteration_type);
 }
+
+inline const Number &
+DisplacedSystem::duDotDu(const unsigned int var_num) const
+{
+  return _undisplaced_system.duDotDu(var_num);
+}

framework/include/systems/LinearSystem.h

@@ -40,11 +40,6 @@ class LinearSystem : public SolverSystem, public PerfGraphInterface
 
   virtual void solve() override;
 
-  virtual void addTimeIntegrator(const std::string & type,
-                                 const std::string & name,
-                                 InputParameters & parameters) override;
-  using SystemBase::addTimeIntegrator;
-
   virtual void initialSetup() override;
 
   // Overriding these to make sure the linear systems don't do anything during

framework/include/systems/NonlinearSystemBase.h

@@ -115,17 +115,6 @@ class NonlinearSystemBase : public SolverSystem, public PerfGraphInterface
    */
   virtual bool converged() = 0;
 
-  /**
-   * Add a time integrator
-   * @param type Type of the integrator
-   * @param name The name of the integrator
-   * @param parameters Integrator params
-   */
-  void addTimeIntegrator(const std::string & type,
-                         const std::string & name,
-                         InputParameters & parameters) override;
-  using SystemBase::addTimeIntegrator;
-
   /**
    * Adds a kernel
    * @param kernel_name The type of the kernel

framework/include/systems/SystemBase.h

@@ -244,9 +244,9 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface
    */
   virtual void addDotVectors();
 
-  virtual Number & duDotDu() { return _du_dot_du; }
+  virtual std::vector<Number> & duDotDus() { return _du_dot_du; }
   virtual Number & duDotDotDu() { return _du_dotdot_du; }
-  virtual const Number & duDotDu() const { return _du_dot_du; }
+  virtual const Number & duDotDu(unsigned int var_num = 0) const;
   virtual const Number & duDotDotDu() const { return _du_dotdot_du; }
 
   virtual NumericVector<Number> * solutionUDot() { return _u_dot; }
@@ -878,18 +878,9 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface
    */
   virtual void copySolutionsBackwards();
 
-  virtual void addTimeIntegrator(const std::string & /*type*/,
-                                 const std::string & /*name*/,
-                                 InputParameters & /*parameters*/)
-  {
-  }
-
-  virtual void addTimeIntegrator(std::shared_ptr<TimeIntegrator> /*ti*/) {}
-
-  TimeIntegrator * getTimeIntegrator() { return _time_integrator.get(); }
-  const TimeIntegrator * getTimeIntegrator() const { return _time_integrator.get(); }
-
-  std::shared_ptr<TimeIntegrator> getSharedTimeIntegrator() { return _time_integrator; }
+  void addTimeIntegrator(const std::string & type,
+                         const std::string & name,
+                         InputParameters & parameters);
 
   /// Whether or not there are variables to be restarted from an Exodus mesh file
   bool hasVarCopy() const { return _var_to_copy.size() > 0; }
@@ -950,6 +941,28 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface
    */
   virtual void compute(ExecFlagType type) = 0;
 
+  /**
+   * Copy time integrators from another system
+   */
+  void copyTimeIntegrators(const SystemBase & other_sys);
+
+  /**
+   * Retrieve the time integrator that integrates the given variable's equation
+   */
+  const TimeIntegrator & getTimeIntegrator(const unsigned int var_num) const;
+
+  /**
+   * Retrieve the time integrator that integrates the given variable's equation. If no suitable time
+   * integrator is found (this could happen for instance if we're solving a non-transient problem),
+   * then a nullptr will be returned
+   */
+  const TimeIntegrator * queryTimeIntegrator(const unsigned int var_num) const;
+
+  /**
+   * @returns All the time integrators owned by this system
+   */
+  const std::vector<std::shared_ptr<TimeIntegrator>> & getTimeIntegrators();
+
 protected:
   /**
    * Internal getter for solution owned by libMesh.
@@ -990,7 +1003,9 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface
   /// old solution vector for u^dotdot
   NumericVector<Number> * _u_dotdot_old;
 
-  Real _du_dot_du;
+  /// Derivative of time derivative of u with respect to uj. This depends on the time integration
+  /// scheme
+  std::vector<Real> _du_dot_du;
   Real _du_dotdot_du;
 
   /// Tagged vectors (pointer)
@@ -1020,7 +1035,7 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface
   size_t _max_var_n_dofs_per_node;
 
   /// Time integrator
-  std::shared_ptr<TimeIntegrator> _time_integrator;
+  std::vector<std::shared_ptr<TimeIntegrator>> _time_integrators;
 
   /// Map variable number to its pointer
   std::vector<std::vector<MooseVariableFieldBase *>> _numbered_vars;

framework/include/timeintegrators/AStableDirk4.h

@@ -62,6 +62,7 @@ class AStableDirk4 : public TimeIntegrator
                                 ADReal & ad_u_dotdot) const override;
   virtual void solve() override;
   virtual void postResidual(NumericVector<Number> & residual) override;
+  virtual bool overridesSolve() const override { return true; }
 
 protected:
   /**

framework/include/timeintegrators/ActuallyExplicitEuler.h

@@ -29,6 +29,7 @@ class ActuallyExplicitEuler : public ExplicitTimeIntegrator
                                 ADReal & ad_u_dotdot) const override;
   virtual void solve() override;
   virtual void postResidual(NumericVector<Number> & residual) override;
+  virtual bool overridesSolve() const override { return true; }
 
 protected:
   /**

framework/include/timeintegrators/BDF2.h

@@ -29,6 +29,7 @@ class BDF2 : public TimeIntegrator
                                 const dof_id_type & dof,
                                 ADReal & ad_u_dotdot) const override;
   virtual void postResidual(NumericVector<Number> & residual) override;
+  virtual bool overridesSolve() const override { return false; }
 
 protected:
   /**
@@ -38,6 +39,8 @@ class BDF2 : public TimeIntegrator
   void
   computeTimeDerivativeHelper(T & u_dot, const T2 & u, const T3 & u_old, const T4 & u_older) const;
 
+  virtual Real duDotDuCoeff() const override;
+
   std::vector<Real> & _weight;
 
   /// The older solution

framework/include/timeintegrators/CentralDifference.h

@@ -31,6 +31,8 @@ class CentralDifference : public ActuallyExplicitEuler
                                 ADReal & ad_u_dotdot) const override;
 
 protected:
+  virtual Real duDotDuCoeff() const override;
+
   /// solution vector for \f$ {du^dotdot}\over{du} \f$
   Real & _du_dotdot_du;
 

framework/include/timeintegrators/CrankNicolson.h

@@ -34,6 +34,7 @@ class CrankNicolson : public TimeIntegrator
                                 ADReal & ad_u_dotdot) const override;
   virtual void postResidual(NumericVector<Number> & residual) override;
   virtual void postStep() override;
+  virtual bool overridesSolve() const override { return false; }
 
 protected:
   /**
@@ -42,6 +43,8 @@ class CrankNicolson : public TimeIntegrator
   template <typename T, typename T2>
   void computeTimeDerivativeHelper(T & u_dot, const T2 & u_old) const;
 
+  virtual Real duDotDuCoeff() const override;
+
   NumericVector<Number> & _residual_old;
 };
 

framework/include/timeintegrators/ExplicitEuler.h

@@ -28,6 +28,7 @@ class ExplicitEuler : public TimeIntegrator
                                 const dof_id_type & dof,
                                 ADReal & ad_u_dotdot) const override;
   virtual void postResidual(NumericVector<Number> & residual) override;
+  virtual bool overridesSolve() const override { return false; }
 
 protected:
   /**

framework/include/timeintegrators/ExplicitRK2.h

@@ -67,6 +67,7 @@ class ExplicitRK2 : public TimeIntegrator
                                 ADReal & ad_u_dotdot) const override;
   virtual void solve() override;
   virtual void postResidual(NumericVector<Number> & residual) override;
+  virtual bool overridesSolve() const override { return true; }
 
 protected:
   /**

framework/include/timeintegrators/ExplicitSSPRungeKutta.h

@@ -28,6 +28,7 @@ class ExplicitSSPRungeKutta : public ExplicitTimeIntegrator
   virtual void solve() override;
   virtual void postResidual(NumericVector<Number> & residual) override;
   virtual int order() override { return _order; }
+  virtual bool overridesSolve() const override { return true; }
 
 protected:
   /**
@@ -37,6 +38,8 @@ class ExplicitSSPRungeKutta : public ExplicitTimeIntegrator
    */
   bool solveStage();
 
+  virtual Real duDotDuCoeff() const override;
+
   /// Order of time integration
   const MooseEnum & _order;
 

framework/include/timeintegrators/ExplicitTVDRK2.h

@@ -58,6 +58,7 @@ class ExplicitTVDRK2 : public TimeIntegrator
                                 ADReal & ad_u_dotdot) const override;
   virtual void solve() override;
   virtual void postResidual(NumericVector<Number> & residual) override;
+  virtual bool overridesSolve() const override { return true; }
 
 protected:
   /**

framework/include/timeintegrators/ImplicitEuler.h

@@ -28,6 +28,7 @@ class ImplicitEuler : public TimeIntegrator
                                 const dof_id_type & dof,
                                 ADReal & ad_u_dotdot) const override;
   virtual void postResidual(NumericVector<Number> & residual) override;
+  virtual bool overridesSolve() const override { return false; }
 
 protected:
   /**

framework/include/timeintegrators/ImplicitMidpoint.h

@@ -53,6 +53,7 @@ class ImplicitMidpoint : public TimeIntegrator
                                 ADReal & ad_u_dotdot) const override;
   virtual void solve() override;
   virtual void postResidual(NumericVector<Number> & residual) override;
+  virtual bool overridesSolve() const override { return true; }
 
 protected:
   /**