Merge pull request #14658 from cbolisetti/explicit-dynamics

Explicit dynamics

framework/doc/content/source/timeintegrators/CentralDifference.md

@@ -0,0 +1,22 @@
+# CentralDifference
+
+!syntax description /Executioner/TimeIntegrator/CentralDifference
+
+## Description
+
+Central difference integration is a common explicit integration scheme, typically used in structural dynamics. For a variable, $u$, at time $t$ with a time step, $\Delta t$, the central difference approximations for the first and second time derivatives, $\dot{u}$ and $\ddot{u}$, are given as,
+
+\begin{equation}
+\begin{aligned}
+\mathbf{\dot{u}}(t) &=& \frac{\mathbf{u}(t+\Delta t)-\mathbf{u}(t-\Delta t)}{2 \Delta t} \\
+\mathbf{\ddot{u}}(t) &=& \frac{\mathbf{u}(t+\Delta t)-2\mathbf{u}(t)+\mathbf{u}(t-\Delta t)}{\Delta t^2}
+\end{aligned}
+\end{equation}
+
+The central difference time integrator derives from the [`ActuallyExplicitEuler`](ActuallyExplicitEuler.md) class and therefore circumvents the nonlinear solver. It can be used with `consistent`, `lumped` or `lump_preconditioned`, `solve_type` options. Information on these solve options can be found on the [`ActuallyExplicitEuler`](ActuallyExplicitEuler.md) page.
+
+!syntax parameters /Executioner/TimeIntegrator/CentralDifference
+
+!syntax inputs /Executioner/TimeIntegrator/CentralDifference
+
+!syntax children /Executioner/TimeIntegrator/CentralDifference

framework/doc/content/source/timeintegrators/ExplicitTimeIntegrator.md

@@ -65,11 +65,11 @@ is used in the linear system, but the linear solve is preconditioned using the
 lumped mass matrix. This compromise retains the accuracy of the `consistent`
 option while benefiting from some of the speedup offered by the `lumped` option.
 
-The lumped mass matrix preconditioner is applied with a helper class, `LumpedPreconditioner`:
+The lumped mass matrix preconditioner is applied with the class, `LumpedPreconditioner`:
 
-!listing framework/include/timeintegrators/ExplicitTimeIntegrator.h line=class LumpedPreconditioner
+!listing framework/include/timeintegrators/LumpedPreconditioner.h
 
-This helper object simply applies the inverse of the diagonal, lumped
+This object simply applies the inverse of the diagonal, lumped
 mass-matrix as the preconditioner for the linear solve, which is very
 efficient. Note that when this option is applied you shouldn't specify any other
 preconditioners using command-line syntax or they will override this option.

framework/include/interfaces/Coupleable.h

@@ -164,6 +164,16 @@ class Coupleable
   virtual const VariableValue &
   coupledVectorTagValue(const std::string & var_name, TagID tag, unsigned int comp = 0);
 
+  /**
+   * Returns dof value of a coupled variable for a given tag
+   * @param var_name Name of coupled variable
+   * @param tag vector tag ID
+   * @param comp Component number for vector of coupled variables
+   * @return Reference to a DofValue for the coupled variable
+   */
+  virtual const VariableValue &
+  coupledVectorTagDofValue(const std::string & var_name, TagID tag, unsigned int comp = 0);
+
   /**
    * Returns value of a coupled variable for a given tag. This couples the diag vector of matrix
    * @param var_name Name of coupled variable
@@ -514,23 +524,25 @@ class Coupleable
    * 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);
 
   /**
    * 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);

framework/include/kernels/TimeKernel.h

@@ -30,6 +30,14 @@ class TimeKernel : public Kernel
 
   virtual void computeResidual() override;
 
+  /**
+   * Entry point for additional computation for the local residual after the standard calls
+   *
+   * NOTE: This is an advanced interface and in nearly all cases you should just override
+   * computeQpResidual()!
+   */
+  virtual void computeResidualAdditional() {}
+
 protected:
   /// Time derivative of u
   const VariableValue & _u_dot;

framework/include/systems/AuxiliarySystem.h

@@ -43,7 +43,6 @@ class AuxiliarySystem : public SystemBase, public PerfGraphInterface
   AuxiliarySystem(FEProblemBase & subproblem, const std::string & name);
   virtual ~AuxiliarySystem();
 
-  virtual void init() override;
   virtual void addExtraVectors() override;
 
   virtual void initialSetup();
@@ -158,17 +157,7 @@ class AuxiliarySystem : public SystemBase, public PerfGraphInterface
    */
   bool needMaterialOnSide(BoundaryID bnd_id);
 
-  NumericVector<Number> & solution() override { return *_sys.solution; }
-  NumericVector<Number> & solutionOld() override { return *_sys.old_local_solution; }
-  NumericVector<Number> & solutionOlder() override { return *_sys.older_local_solution; }
   NumericVector<Number> * solutionPreviousNewton() override { return _solution_previous_nl; }
-
-  const NumericVector<Number> & solution() const override { return *_sys.solution; }
-  const NumericVector<Number> & solutionOld() const override { return *_sys.old_local_solution; }
-  const NumericVector<Number> & solutionOlder() const override
-  {
-    return *_sys.older_local_solution;
-  }
   const NumericVector<Number> * solutionPreviousNewton() const override
   {
     return _solution_previous_nl;
@@ -185,7 +174,7 @@ class AuxiliarySystem : public SystemBase, public PerfGraphInterface
 
   void clearScalarVariableCoupleableTags();
 
-  // protected:
+protected:
   void computeScalarVars(ExecFlagType type);
   void computeNodalVars(ExecFlagType type);
   void computeNodalVecVars(ExecFlagType type);
@@ -212,8 +201,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 +211,9 @@ class AuxiliarySystem : public SystemBase, public PerfGraphInterface
   /// Old solution vector for u^dotdot
   NumericVector<Number> * _u_dotdot_old;
 
+  /// The current states of the solution (0 = current, 1 = old, etc)
+  std::vector<NumericVector<Number> *> _solution_state;
+
   /// Whether or not a copy of the residual needs to be made
   bool _need_serialized_solution;
 
@@ -261,4 +251,11 @@ class AuxiliarySystem : public SystemBase, public PerfGraphInterface
   friend class ComputeNodalKernelBcsThread;
   friend class ComputeNodalKernelJacobiansThread;
   friend class ComputeNodalKernelBCJacobiansThread;
+
+  NumericVector<Number> & solutionInternal() const override { return *_sys.solution; }
+  NumericVector<Number> & solutionOldInternal() const override { return *_sys.old_local_solution; }
+  NumericVector<Number> & solutionOlderInternal() const override
+  {
+    return *_sys.older_local_solution;
+  }
 };

framework/include/systems/DisplacedSystem.h

@@ -30,8 +30,6 @@ class DisplacedSystem : public SystemBase
                   Moose::VarKindType var_kind);
   virtual ~DisplacedSystem();
 
-  virtual void init() override;
-
   virtual NumericVector<Number> & getVector(TagID tag_id) override
   {
     return _undisplaced_system.getVector(tag_id);
@@ -91,15 +89,8 @@ class DisplacedSystem : public SystemBase
     return _undisplaced_system.currentSolution();
   }
 
-  NumericVector<Number> & solution() override { return _undisplaced_system.solution(); }
-  NumericVector<Number> & solutionOld() override;
-  NumericVector<Number> & solutionOlder() override;
-  NumericVector<Number> * solutionPreviousNewton() override { return NULL; }
-
-  const NumericVector<Number> & solution() const override { return _undisplaced_system.solution(); }
-  const NumericVector<Number> & solutionOld() const override;
-  const NumericVector<Number> & solutionOlder() const override;
-  const NumericVector<Number> * solutionPreviousNewton() const override { return NULL; }
+  NumericVector<Number> * solutionPreviousNewton() override { return nullptr; }
+  const NumericVector<Number> * solutionPreviousNewton() const override { return nullptr; }
 
   NumericVector<Number> * solutionUDot() override { return _undisplaced_system.solutionUDot(); }
   NumericVector<Number> * solutionUDotDot() override
@@ -218,7 +209,13 @@ class DisplacedSystem : public SystemBase
   void addTimeIntegrator(std::shared_ptr<TimeIntegrator> ti) override;
 
 protected:
+  NumericVector<Number> & solutionInternal() const override
+  {
+    return _undisplaced_system.solution();
+  }
+  NumericVector<Number> & solutionOldInternal() const override;
+  NumericVector<Number> & solutionOlderInternal() const override;
+
   SystemBase & _undisplaced_system;
   TransientExplicitSystem & _sys;
 };
-

framework/include/systems/DumpObjectsNonlinearSystem.h

@@ -32,16 +32,12 @@ class DumpObjectsNonlinearSystem : public NonlinearSystemBase
   virtual bool converged() override { return true; }
   virtual NumericVector<Number> & RHS() override { return *_dummy; }
 
-  NumericVector<Number> & solutionOld() override { return *_dummy; }
-  NumericVector<Number> & solutionOlder() override { return *_dummy; }
-
-  const NumericVector<Number> & solutionOld() const override { return *_dummy; }
-  const NumericVector<Number> & solutionOlder() const override { return *_dummy; }
-
   virtual unsigned int getCurrentNonlinearIterationNumber() override { return 0; }
   virtual void setupFiniteDifferencedPreconditioner() override {}
 
 protected:
+  NumericVector<Number> & solutionOldInternal() const override { return *_dummy; }
+  NumericVector<Number> & solutionOlderInternal() const override { return *_dummy; }
+
   NumericVector<Number> * _dummy;
 };
-

framework/include/systems/NonlinearEigenSystem.h

@@ -84,20 +84,6 @@ class NonlinearEigenSystem : public NonlinearSystemBase
 
   virtual NonlinearSolver<Number> * nonlinearSolver() override;
 
-  NumericVector<Number> & solutionOld() override { return *_transient_sys.old_local_solution; }
-
-  NumericVector<Number> & solutionOlder() override { return *_transient_sys.older_local_solution; }
-
-  const NumericVector<Number> & solutionOld() const override
-  {
-    return *_transient_sys.old_local_solution;
-  }
-
-  const NumericVector<Number> & solutionOlder() const override
-  {
-    return *_transient_sys.older_local_solution;
-  }
-
   virtual TransientEigenSystem & sys() { return _transient_sys; }
 
   /**
@@ -147,6 +133,16 @@ class NonlinearEigenSystem : public NonlinearSystemBase
   TagID nonEigenMatrixTag() { return _A_tag; }
 
 protected:
+  NumericVector<Number> & solutionOldInternal() const override
+  {
+    return *_transient_sys.old_local_solution;
+  }
+
+  NumericVector<Number> & solutionOlderInternal() const override
+  {
+    return *_transient_sys.older_local_solution;
+  }
+
   TransientEigenSystem & _transient_sys;
   EigenProblem & _eigen_problem;
   std::vector<std::pair<Real, Real>> _eigen_values;

framework/include/systems/NonlinearSystem.h

@@ -64,23 +64,20 @@ class NonlinearSystem : public NonlinearSystemBase
     return _transient_sys.nonlinear_solver.get();
   }
 
-  NumericVector<Number> & solutionOld() override { return *_transient_sys.old_local_solution; }
-  const NumericVector<Number> & solutionOld() const override
+  virtual TransientNonlinearImplicitSystem & sys() { return _transient_sys; }
+
+  void computeScaling() override;
+
+protected:
+  NumericVector<Number> & solutionOldInternal() const override
   {
     return *_transient_sys.old_local_solution;
   }
-
-  NumericVector<Number> & solutionOlder() override { return *_transient_sys.older_local_solution; }
-  const NumericVector<Number> & solutionOlder() const override
+  NumericVector<Number> & solutionOlderInternal() const override
   {
     return *_transient_sys.older_local_solution;
   }
 
-  virtual TransientNonlinearImplicitSystem & sys() { return _transient_sys; }
-
-  void computeScaling() override;
-
-protected:
   TransientNonlinearImplicitSystem & _transient_sys;
   ComputeResidualFunctor _nl_residual_functor;
   ComputeFDResidualFunctor _fd_residual_functor;

framework/include/systems/NonlinearSystemBase.h

@@ -605,9 +605,6 @@ class NonlinearSystemBase : public SystemBase, public PerfGraphInterface
    */
   bool hasDiagSaveIn() const { return _has_diag_save_in || _has_nodalbc_diag_save_in; }
 
-  NumericVector<Number> & solution() override { return *_sys.solution; }
-  const NumericVector<Number> & solution() const override { return *_sys.solution; }
-
   virtual System & system() override { return _sys; }
   virtual const System & system() const override { return _sys; }
 
@@ -719,6 +716,8 @@ class NonlinearSystemBase : public SystemBase, public PerfGraphInterface
   void mortarConstraints(bool displaced);
 
 protected:
+  NumericVector<Number> & solutionInternal() const override { return *_sys.solution; }
+
   /// solution vector from nonlinear solver
   const NumericVector<Number> * _current_solution;
   /// ghosted form of the residual
@@ -931,4 +930,7 @@ class NonlinearSystemBase : public SystemBase, public PerfGraphInterface
   /// The required size of the derivative storage array
   std::size_t _required_derivative_size;
 #endif
+
+  /// The current states of the solution (0 = current, 1 = old, etc)
+  std::vector<NumericVector<Number> *> _solution_state;
 };

framework/include/systems/SystemBase.h

@@ -212,14 +212,31 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface
    */
   virtual const NumericVector<Number> * const & currentSolution() const = 0;
 
-  virtual NumericVector<Number> & solution() = 0;
-  virtual NumericVector<Number> & solutionOld() = 0;
-  virtual NumericVector<Number> & solutionOlder() = 0;
-  virtual NumericVector<Number> * solutionPreviousNewton() = 0;
-  virtual const NumericVector<Number> & solution() const = 0;
-  virtual const NumericVector<Number> & solutionOld() const = 0;
-  virtual const NumericVector<Number> & solutionOlder() const = 0;
+  NumericVector<Number> & solution() { return solutionState(0); }
+  NumericVector<Number> & solutionOld() { return solutionState(1); }
+  NumericVector<Number> & solutionOlder() { return solutionState(2); }
+  const NumericVector<Number> & solution() const { return solutionState(0); }
+  const NumericVector<Number> & solutionOld() const { return solutionState(1); }
+  const NumericVector<Number> & solutionOlder() const { return solutionState(2); }
+
   virtual const NumericVector<Number> * solutionPreviousNewton() const = 0;
+  virtual NumericVector<Number> * solutionPreviousNewton() = 0;
+
+  /**
+   * Get a state of the solution (0 = current, 1 = old, 2 = older, etc).
+   *
+   * If the state does not exist, it will be initialized in addition to any newer
+   * states before it that have not been initialized.
+   */
+  NumericVector<Number> & solutionState(const unsigned int state);
+
+  /**
+   * Get a state of the solution (0 = current, 1 = old, 2 = older, etc).
+   *
+   * By default, up to state _default_solution_states is added. Any older states must be
+   * added using the non-const solutionState().
+   */
+  const NumericVector<Number> & solutionState(const unsigned int state) const;
 
   virtual Number & duDotDu() { return _du_dot_du; }
   virtual Number & duDotDotDu() { return _du_dotdot_du; }
@@ -807,6 +824,20 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface
   void cacheVarIndicesByFace(const std::vector<VariableName> & vars);
 
 protected:
+  /**
+   * Internal getters for the states of the solution as owned by libMesh.
+   *
+   * For the first three states (0 = current, 1 = old, 2 = older), we point directly to the
+   * solutions in libMesh (which is why these virtuals are needed). This allows us to store a more
+   * generalized set of solution states in _solution_states that also enables the addition of older
+   * states if we need them.
+   */
+  ///@{
+  virtual NumericVector<Number> & solutionInternal() const = 0;
+  virtual NumericVector<Number> & solutionOldInternal() const = 0;
+  virtual NumericVector<Number> & solutionOlderInternal() const = 0;
+  ///@}
+
   SubProblem & _subproblem;
 
   MooseApp & _app;
@@ -875,6 +906,15 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface
 
   /// True if printing out additional information
   bool _verbose;
+
+  /// The number of default solution states to store
+  unsigned int _default_solution_states;
+
+private:
+  /// The solution states (0 = current, 1 = old, 2 = older, etc)
+  std::vector<NumericVector<Number> *> _solution_states;
+  /// The saved solution states (0 = current, 1 = old, 2 = older, etc)
+  std::vector<NumericVector<Number> *> _saved_solution_states;
 };
 
 #define PARALLEL_TRY