Framework changes for CentralDifference timeintegrator. Refs idaholab…

…#13964.

framework/include/interfaces/Coupleable.h

@@ -510,27 +510,51 @@ class Coupleable
    */
   virtual const VariableValue & coupledDot(const std::string & var_name, unsigned int comp = 0);
 
+  /**
+   * Residual corresponding to the time derivative of a coupled variable
+   * Different from time derivative of a coupled variable for explicit solvers
+   * @param var_name Name of coupled variable
+   * @param comp Component number for vector of coupled variables
+   * @return Reference to a VariableValue containing the residual corresponding
+   * to the time derivative of the coupled variable
+   */
+  virtual const VariableValue & coupledDotResidual(const std::string & var_name,
+                                                   unsigned int comp = 0);
+
   /**
    * Second time derivative of a coupled variable
    * @param var_name Name of coupled variable
    * @param comp Component number for vector of coupled variables
-   * @return Reference to a VariableValue containing the time derivative of the coupled variable
+   * @return Reference to a VariableValue containing the second time derivative of the coupled
+   * variable
    */
   virtual const VariableValue & coupledDotDot(const std::string & var_name, unsigned int comp = 0);
 
+  /**
+   * Residual corresponding to the second time derivative of a coupled variable
+   * Different from second time derivative of a coupled variable for explicit solvers
+   * @param var_name Name of coupled variable
+   * @param comp Component number for vector of coupled variables
+   * @ return Reference to a VariableValue containing the residual corresponding
+   * to the second time derivative of the coupled variable.
+   */
+  virtual const VariableValue & coupledDotDotResidual(const std::string & var_name,
+                                                      unsigned int comp = 0);
+
   /**
    * Old time derivative of a coupled variable
    * @param var_name Name of coupled variable
    * @param comp Component number for vector of coupled variables
-   * @return Reference to a VariableValue containing the time derivative of the coupled variable
+   * @return Reference to a VariableValue containing the old time derivative of the coupled variable
    */
   virtual const VariableValue & coupledDotOld(const std::string & var_name, unsigned int comp = 0);
 
   /**
    * Old second time derivative of a coupled variable
    * @param var_name Name of coupled variable
    * @param comp Component number for vector of coupled variables
-   * @return Reference to a VariableValue containing the time derivative of the coupled variable
+   * @return Reference to a VariableValue containing the old second time derivative of the coupled
+   * variable
    */
   virtual const VariableValue & coupledDotDotOld(const std::string & var_name,
                                                  unsigned int comp = 0);
@@ -564,7 +588,8 @@ class Coupleable
    */
   virtual const VectorVariableValue & coupledVectorDot(const std::string & var_name,
                                                        unsigned int comp = 0);
-
+  virtual const VectorVariableValue & coupledVectorDotResidual(const std::string & var_name,
+                                                               unsigned int comp = 0);
   /**
    * Second time derivative of a coupled vector variable
    * @param var_name Name of coupled vector variable
@@ -574,7 +599,8 @@ class Coupleable
    */
   virtual const VectorVariableValue & coupledVectorDotDot(const std::string & var_name,
                                                           unsigned int comp = 0);
-
+  virtual const VectorVariableValue & coupledVectorDotDotResidual(const std::string & var_name,
+                                                                  unsigned int comp = 0);
   /**
    * Old time derivative of a coupled vector variable
    * @param var_name Name of coupled vector variable
@@ -624,6 +650,8 @@ class Coupleable
    */
   virtual const ArrayVariableValue & coupledArrayDot(const std::string & var_name,
                                                      unsigned int comp = 0);
+  virtual const ArrayVariableValue & coupledArrayDotResidual(const std::string & var_name,
+                                                             unsigned int comp = 0);
 
   /**
    * Second time derivative of a coupled array variable
@@ -634,6 +662,8 @@ class Coupleable
    */
   virtual const ArrayVariableValue & coupledArrayDotDot(const std::string & var_name,
                                                         unsigned int comp = 0);
+  virtual const ArrayVariableValue & coupledArrayDotDotResidual(const std::string & var_name,
+                                                                unsigned int comp = 0);
 
   /**
    * Old time derivative of a coupled array variable
@@ -729,7 +759,8 @@ class Coupleable
    */
   template <typename T>
   const T & coupledNodalDot(const std::string & var_name, unsigned int comp = 0);
-
+  template <typename T>
+  const T & coupledNodalDotResidual(const std::string & var_name, unsigned int comp = 0);
   /**
    * Nodal values of second time derivative of a coupled variable
    * @param var_name Name of coupled variable
@@ -739,7 +770,8 @@ class Coupleable
    */
   virtual const VariableValue & coupledNodalDotDot(const std::string & var_name,
                                                    unsigned int comp = 0);
-
+  virtual const VariableValue & coupledNodalDotDotResidual(const std::string & var_name,
+                                                           unsigned int comp = 0);
   /**
    * Nodal values of old time derivative of a coupled variable
    * @param var_name Name of coupled variable

framework/include/interfaces/NeighborCoupleable.h

@@ -37,6 +37,8 @@ class NeighborCoupleable : public Coupleable
                                                      unsigned int comp = 0);
   virtual const VariableValue & coupledNeighborValueDot(const std::string & var_name,
                                                         unsigned int comp = 0);
+  virtual const VariableValue & coupledNeighborValueDotResidual(const std::string & var_name,
+                                                                unsigned int comp = 0);
   virtual const VariableValue & coupledNeighborValueDotDu(const std::string & var_name,
                                                           unsigned int comp = 0);
   virtual const VariableValue & coupledNeighborValueOld(const std::string & var_name,

framework/include/interfaces/ScalarCoupleable.h

@@ -149,6 +149,17 @@ class ScalarCoupleable
    */
   virtual VariableValue & coupledScalarDot(const std::string & var_name, unsigned int comp = 0);
 
+  /**
+   * Returns the residual corresponding to the time derivative of a scalar coupled variable
+   * Different from the time derivative of the scalar variable for explicit solvers
+   * @param var_name Name of coupled variable
+   * @param comp Component number for vector of coupled variables
+   * @return Reference to a VariableValue containing the residual corresponding
+   * to the time derivative for the coupled variable
+   */
+  virtual VariableValue & coupledScalarDotResidual(const std::string & var_name,
+                                                   unsigned int comp = 0);
+
   /**
    * Returns the second time derivative of a scalar coupled variable
    * @param var_name Name of coupled variable
@@ -157,6 +168,17 @@ class ScalarCoupleable
    */
   virtual VariableValue & coupledScalarDotDot(const std::string & var_name, unsigned int comp = 0);
 
+  /**
+   * Returns the residual corresponding to the second time derivative of a scalar coupled variable
+   * Different from the second time derivative of the scalar coupled variable for explicit solvers
+   * @param var_name Name of coupled variable
+   * @param comp Component number for vector of coupled variables
+   * @return Reference to a VariableValue containing the residual corresponding
+   * to the second time derivative of the coupled variable
+   */
+  virtual VariableValue & coupledScalarDotDotResidual(const std::string & var_name,
+                                                      unsigned int comp = 0);
+
   /**
    * Returns the old time derivative of a scalar coupled variable
    * @param var_name Name of coupled variable

framework/include/problems/FEProblemBase.h

@@ -1731,6 +1731,12 @@ class FEProblemBase : public SubProblem, public Restartable
   using SubProblem::haveADObjects;
   void haveADObjects(bool have_ad_objects) override;
 
+  /// If solution older than 2 time steps is required, set solution_state to the number of old timesteps required.
+  virtual void setSolutionState(unsigned int solution_state) { _solution_state = solution_state; }
+
+  /// Get the number of old solution states to be stored
+  virtual unsigned int getSolutionState() const { return _solution_state; }
+
   // Whether or not we should solve this system
   bool shouldSolve() const { return _solve; }
 
@@ -2141,6 +2147,9 @@ class FEProblemBase : public SubProblem, public Restartable
   /// Whether old solution second time derivative needs to be stored
   bool _u_dotdot_old_requested;
 
+  /// Number of old solution states to be stored;
+  unsigned int _solution_state;
+
   friend class AuxiliarySystem;
   friend class NonlinearSystemBase;
   friend class MooseEigenSystem;

framework/include/systems/AuxiliarySystem.h

@@ -179,6 +179,8 @@ class AuxiliarySystem : public SystemBase, public PerfGraphInterface
   virtual System & system() override { return _sys; }
   virtual const System & system() const override { return _sys; }
 
+  virtual NumericVector<Number> * solutionState(unsigned int i) override;
+
   virtual void setPreviousNewtonSolution();
 
   void setScalarVariableCoupleableTags(ExecFlagType type);
@@ -212,8 +214,6 @@ class AuxiliarySystem : public SystemBase, public PerfGraphInterface
   NumericVector<Number> & _serialized_solution;
   /// Solution vector of the previous nonlinear iterate
   NumericVector<Number> * _solution_previous_nl;
-  /// Time integrator
-  std::shared_ptr<TimeIntegrator> _time_integrator;
   /// solution vector for u^dot
   NumericVector<Number> * _u_dot;
   /// solution vector for u^dotdot
@@ -224,6 +224,8 @@ class AuxiliarySystem : public SystemBase, public PerfGraphInterface
   /// Old solution vector for u^dotdot
   NumericVector<Number> * _u_dotdot_old;
 
+  std::vector<NumericVector<Number> *> _solution_state;
+
   /// Whether or not a copy of the residual needs to be made
   bool _need_serialized_solution;
 

framework/include/systems/DisplacedSystem.h

@@ -209,6 +209,11 @@ class DisplacedSystem : public SystemBase
     return _undisplaced_system.getMatrix(tag);
   }
 
+  virtual NumericVector<Number> * solutionState(unsigned int i) override
+  {
+    return _undisplaced_system.solutionState(i);
+  }
+
   virtual TransientExplicitSystem & sys() { return _sys; }
 
   virtual System & system() override;

framework/include/systems/NonlinearSystemBase.h

@@ -617,10 +617,14 @@ class NonlinearSystemBase : public SystemBase, public PerfGraphInterface
     return _solution_previous_nl;
   }
 
+  virtual NumericVector<Number> * solutionState(unsigned int i) override;
+
   virtual void setSolutionUDotOld(const NumericVector<Number> & u_dot_old);
 
   virtual void setSolutionUDotDotOld(const NumericVector<Number> & u_dotdot_old);
 
+  virtual void setSolutionState(const std::vector<NumericVector<Number>> & solution_state);
+
   virtual void setPreviousNewtonSolution(const NumericVector<Number> & soln);
 
   virtual TagID timeVectorTag() override { return _Re_time_tag; }
@@ -931,4 +935,6 @@ class NonlinearSystemBase : public SystemBase, public PerfGraphInterface
   /// The required size of the derivative storage array
   std::size_t _required_derivative_size;
 #endif
+
+  std::vector<NumericVector<Number> *> _solution_state;
 };

framework/include/systems/SystemBase.h

@@ -215,6 +215,7 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface
   virtual NumericVector<Number> & solution() = 0;
   virtual NumericVector<Number> & solutionOld() = 0;
   virtual NumericVector<Number> & solutionOlder() = 0;
+  virtual NumericVector<Number> * solutionState(unsigned int i) = 0;
   virtual NumericVector<Number> * solutionPreviousNewton() = 0;
   virtual const NumericVector<Number> & solution() const = 0;
   virtual const NumericVector<Number> & solutionOld() const = 0;
@@ -861,6 +862,9 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface
   /// Map variable number to its pointer
   std::vector<std::vector<MooseVariableFEBase *>> _numbered_vars;
 
+  std::vector<NumericVector<Real> *> _saved_solution_state;
+  unsigned int _solution_state_size;
+
   /// Storage for MooseVariable objects
   MooseObjectWarehouseBase<MooseVariableBase> _variable_warehouse;
 

framework/include/variables/MooseVariableData.h

@@ -317,8 +317,12 @@ class MooseVariableData
 
   const FieldVariableValue & uDot() const;
 
+  const FieldVariableValue & uDotResidual() const;
+
   const FieldVariableValue & uDotDot() const;
 
+  const FieldVariableValue & uDotDotResidual() const;
+
   const FieldVariableValue & uDotOld() const;
 
   const FieldVariableValue & uDotDotOld() const;
@@ -341,6 +345,8 @@ class MooseVariableData
   const MooseArray<OutputType> & nodalValueArray(Moose::SolutionState state) const;
   const OutputType & nodalValueDot() const;
   const OutputType & nodalValueDotDot() const;
+  const OutputType & nodalValueDotResidual() const;
+  const OutputType & nodalValueDotDotResidual() const;
   const OutputType & nodalValueDotOld() const;
   const OutputType & nodalValueDotDotOld() const;
   const OutputType & nodalValueDuDotDu() const;
@@ -429,8 +435,10 @@ class MooseVariableData
   const DoFValue & dofValuesOlder() const;
   const DoFValue & dofValuesPreviousNL() const;
   const DoFValue & dofValuesDot() const;
+  const DoFValue & dofValuesDotResidual() const;
   const DoFValue & dofValuesDotOld() const;
   const DoFValue & dofValuesDotDot() const;
+  const DoFValue & dofValuesDotDotResidual() const;
   const DoFValue & dofValuesDotDotOld() const;
   const MooseArray<Number> & dofValuesDuDotDu() const;
   const MooseArray<Number> & dofValuesDuDotDotDu() const;
@@ -547,6 +555,10 @@ class MooseVariableData
   OutputType _nodal_value_dot;
   /// nodal values of u_dotdot
   OutputType _nodal_value_dotdot;
+  /// nodal values of u_dot_residual
+  OutputType _nodal_value_dot_residual;
+  /// nodal values of u_dotdot_residual
+  OutputType _nodal_value_dotdot_residual;
   /// nodal values of u_dot_old
   OutputType _nodal_value_dot_old;
   /// nodal values of u_dotdot_old
@@ -572,6 +584,8 @@ class MooseVariableData
   mutable bool _need_u_dot;
   mutable bool _need_ad_u_dot;
   mutable bool _need_u_dotdot;
+  mutable bool _need_u_dot_residual;
+  mutable bool _need_u_dotdot_residual;
   mutable bool _need_u_dot_old;
   mutable bool _need_u_dotdot_old;
   mutable bool _need_du_dot_du;
@@ -610,6 +624,8 @@ class MooseVariableData
   mutable bool _need_dof_values_previous_nl;
   mutable bool _need_dof_values_dot;
   mutable bool _need_dof_values_dotdot;
+  mutable bool _need_dof_values_dot_residual;
+  mutable bool _need_dof_values_dotdot_residual;
   mutable bool _need_dof_values_dot_old;
   mutable bool _need_dof_values_dotdot_old;
   mutable bool _need_dof_du_dot_du;
@@ -628,6 +644,10 @@ class MooseVariableData
   DoFValue _dof_values_dot;
   /// nodal values of u_dotdot
   DoFValue _dof_values_dotdot;
+  /// nodal values of u_dot_residual
+  DoFValue _dof_values_dot_residual;
+  /// nodal values of u_dotdot_residual
+  DoFValue _dof_values_dotdot_residual;
   /// nodal values of u_dot_old
   DoFValue _dof_values_dot_old;
   /// nodal values of u_dotdot_old
@@ -680,6 +700,12 @@ class MooseVariableData
   /// u_dotdot (second time derivative)
   FieldVariableValue _u_dotdot, _u_dotdot_bak;
 
+  /// u_dot_residual (residual corresponding to first time derivative)
+  FieldVariableValue _u_dot_residual, _u_dot_residual_bak;
+
+  /// u_dotdot_residual (residual corresponding to second time derivative)
+  FieldVariableValue _u_dotdot_residual, _u_dotdot_residual_bak;
+
   /// u_dot_old (time derivative)
   FieldVariableValue _u_dot_old, _u_dot_old_bak;