Apply suggestions from code review

- Remove commented method - `getPossiblyNullTimeIntegrator` -> `queryTimeIntegrator` - `const` addition to `from` vector in `copyVector` - Add method name to `mooseError` in `TimeIntegrator::solve` - Add comment regarding handling of null vectors - Create computeDuDotDu and duDotDuCoeff methods Co-authored-by: Logan Harbour <[email protected]>
idaholab · Oct 28, 2024 · 69c996e · 69c996e
1 parent af5b6d7
commit 69c996e
Show file tree

Hide file tree

Showing 29 changed files with 98 additions and 72 deletions.
diff --git a/framework/include/systems/SystemBase.h b/framework/include/systems/SystemBase.h
@@ -882,8 +882,6 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface
                          const std::string & name,
                          InputParameters & parameters);
 
-  // virtual void addTimeIntegrator(std::shared_ptr<TimeIntegrator> /*ti*/) {}
-
   /// Whether or not there are variables to be restarted from an Exodus mesh file
   bool hasVarCopy() const { return _var_to_copy.size() > 0; }
 
@@ -958,7 +956,7 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface
    * integrator is found (this could happen for instance if we're solving a non-transient problem),
    * then a nullptr will be returned
    */
-  const TimeIntegrator * getPossiblyNullTimeIntegrator(const unsigned int var_num) const;
+  const TimeIntegrator * queryTimeIntegrator(const unsigned int var_num) const;
 
   /**
    * @returns All the time integrators owned by this system

diff --git a/framework/include/timeintegrators/BDF2.h b/framework/include/timeintegrators/BDF2.h
@@ -39,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

diff --git a/framework/include/timeintegrators/CentralDifference.h b/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;
 

diff --git a/framework/include/timeintegrators/CrankNicolson.h b/framework/include/timeintegrators/CrankNicolson.h
@@ -43,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;
 };
 

diff --git a/framework/include/timeintegrators/ExplicitSSPRungeKutta.h b/framework/include/timeintegrators/ExplicitSSPRungeKutta.h
@@ -38,6 +38,8 @@ class ExplicitSSPRungeKutta : public ExplicitTimeIntegrator
    */
   bool solveStage();
 
+  virtual Real duDotDuCoeff() const override;
+
   /// Order of time integration
   const MooseEnum & _order;
 

diff --git a/framework/include/timeintegrators/NewmarkBeta.h b/framework/include/timeintegrators/NewmarkBeta.h
@@ -42,6 +42,8 @@ class NewmarkBeta : public TimeIntegrator
                                    T4 & u_dotdot,
                                    const T5 & u_dotdot_old) const;
 
+  virtual Real duDotDuCoeff() const override;
+
   /// Newmark time integration parameter-beta
   Real _beta;
 

diff --git a/framework/include/timeintegrators/TimeIntegrator.h b/framework/include/timeintegrators/TimeIntegrator.h
@@ -183,7 +183,18 @@ class TimeIntegrator : public MooseObject, public Restartable
    * Copy from one vector into another. If the time integrator has been restricted to a subset of
    * variables, then this will selectively copy their dofs
    */
-  void copyVector(NumericVector<Number> & from, NumericVector<Number> & to);
+  void copyVector(const NumericVector<Number> & from, NumericVector<Number> & to);
+
+  /**
+   * @returns The \p _du_dot_du multiplicative coefficient, e.g. if \p _du_dot_du is equivalent to
+   * 2/dt, then this method returns 2
+   */
+  virtual Real duDotDuCoeff() const { return 1; }
+
+  /**
+   * Compute \p _du_dot_du
+   */
+  void computeDuDotDu();
 
   FEProblemBase & _fe_problem;
   SystemBase & _sys;
@@ -235,4 +246,7 @@ class TimeIntegrator : public MooseObject, public Restartable
 
   /// The variables that this time integrator integrates
   std::unordered_set<unsigned int> & _vars;
+
+  /// A vector that is used for creating 'from' subvectors in the \p copyVector() method
+  std::unique_ptr<NumericVector<Number>> _from_subvector;
 };
diff --git a/framework/src/restart/DataIO.C b/framework/src/restart/DataIO.C
@@ -326,6 +326,9 @@ dataStore(std::ostream & stream,
           std::unique_ptr<libMesh::NumericVector<Number>> & v,
           void * context)
 {
+  // Classes may declare unique pointers to vectors as restartable data and never actually create
+  // vector instances. This happens for example in the `TimeIntegrator` class where subvector
+  // instances are only created if multiple time integrators are present
   bool have_vector = v.get();
   dataStore(stream, have_vector, context);
   if (!have_vector)

diff --git a/framework/src/systems/SystemBase.C b/framework/src/systems/SystemBase.C
@@ -1596,7 +1596,7 @@ SystemBase::copyTimeIntegrators(const SystemBase & other_sys)
 }
 
 const TimeIntegrator *
-SystemBase::getPossiblyNullTimeIntegrator(const unsigned int var_num) const
+SystemBase::queryTimeIntegrator(const unsigned int var_num) const
 {
   for (auto & ti : _time_integrators)
     if (ti->integratesVar(var_num))
@@ -1608,7 +1608,7 @@ SystemBase::getPossiblyNullTimeIntegrator(const unsigned int var_num) const
 const TimeIntegrator &
 SystemBase::getTimeIntegrator(const unsigned int var_num) const
 {
-  const auto * const ti = getPossiblyNullTimeIntegrator(var_num);
+  const auto * const ti = queryTimeIntegrator(var_num);
 
   if (ti)
     return *ti;

diff --git a/framework/src/timeintegrators/AStableDirk4.C b/framework/src/timeintegrators/AStableDirk4.C
@@ -92,9 +92,7 @@ AStableDirk4::computeTimeDerivatives()
   u_dot = *_solution;
   computeTimeDerivativeHelper(u_dot, _solution_old);
   u_dot.close();
-  for (const auto i : index_range(_du_dot_du))
-    if (integratesVar(i))
-      _du_dot_du[i] = 1. / _dt;
+  computeDuDotDu();
 }
 
 void

diff --git a/framework/src/timeintegrators/ActuallyExplicitEuler.C b/framework/src/timeintegrators/ActuallyExplicitEuler.C
@@ -49,10 +49,7 @@ ActuallyExplicitEuler::computeTimeDerivatives()
   u_dot = *_solution;
   computeTimeDerivativeHelper(u_dot, _solution_old);
   u_dot.close();
-
-  for (const auto i : index_range(_du_dot_du))
-    if (integratesVar(i))
-      _du_dot_du[i] = 1.0 / _dt;
+  computeDuDotDu();
 }
 
 void

diff --git a/framework/src/timeintegrators/BDF2.C b/framework/src/timeintegrators/BDF2.C
@@ -50,21 +50,12 @@ BDF2::computeTimeDerivatives()
 
   NumericVector<Number> & u_dot = *_sys.solutionUDot();
   if (_t_step == 1)
-  {
     u_dot = *_solution;
-    for (const auto i : index_range(_du_dot_du))
-      if (integratesVar(i))
-        _du_dot_du[i] = 1. / _dt;
-  }
   else
-  {
     u_dot.zero();
-    for (const auto i : index_range(_du_dot_du))
-      if (integratesVar(i))
-        _du_dot_du[i] = _weight[0] / _dt;
-  }
   computeTimeDerivativeHelper(u_dot, *_solution, _solution_old, _solution_older);
   u_dot.close();
+  computeDuDotDu();
 }
 
 void
@@ -85,3 +76,12 @@ BDF2::postResidual(NumericVector<Number> & residual)
   residual += _Re_non_time;
   residual.close();
 }
+
+Real
+BDF2::duDotDuCoeff() const
+{
+  if (_t_step == 1)
+    return 1;
+  else
+    return _weight[0];
+}
diff --git a/framework/src/timeintegrators/CentralDifference.C b/framework/src/timeintegrators/CentralDifference.C
@@ -87,9 +87,7 @@ CentralDifference::computeTimeDerivatives()
   u_dot.close();
 
   // used for Jacobian calculations
-  for (const auto i : index_range(_du_dot_du))
-    if (integratesVar(i))
-      _du_dot_du[i] = 1.0 / (2 * _dt);
+  computeDuDotDu();
   _du_dotdot_du = 1.0 / (_dt * _dt);
 
   // Computing udotdot "factor"
@@ -116,3 +114,9 @@ CentralDifference::computeTimeDerivatives()
     u_dot_factor.close();
   }
 }
+
+Real
+CentralDifference::duDotDuCoeff() const
+{
+  return Real(1) / Real(2);
+}
diff --git a/framework/src/timeintegrators/CrankNicolson.C b/framework/src/timeintegrators/CrankNicolson.C
@@ -76,9 +76,7 @@ CrankNicolson::init()
   // time derivative is assumed to be zero on initial
   NumericVector<Number> & u_dot = *_sys.solutionUDot();
   u_dot.zero();
-  for (const auto i : index_range(_du_dot_du))
-    if (integratesVar(i))
-      _du_dot_du[i] = 0;
+  computeDuDotDu();
 
   // compute residual for the initial time step
   // Note: we can not directly pass _residual_old in computeResidualTag because
@@ -141,3 +139,9 @@ CrankNicolson::postStep()
 {
   copyVector(_Re_non_time, _residual_old);
 }
+
+Real
+CrankNicolson::duDotDuCoeff() const
+{
+  return 2;
+}
diff --git a/framework/src/timeintegrators/ExplicitEuler.C b/framework/src/timeintegrators/ExplicitEuler.C
@@ -43,10 +43,7 @@ ExplicitEuler::computeTimeDerivatives()
   u_dot = *_solution;
   computeTimeDerivativeHelper(u_dot, _solution_old);
   u_dot.close();
-
-  for (const auto i : index_range(_du_dot_du))
-    if (integratesVar(i))
-      _du_dot_du[i] = 1.0 / _dt;
+  computeDuDotDu();
 }
 
 void

diff --git a/framework/src/timeintegrators/ExplicitRK2.C b/framework/src/timeintegrators/ExplicitRK2.C
@@ -53,11 +53,7 @@ ExplicitRK2::computeTimeDerivatives()
   NumericVector<Number> & u_dot = *_sys.solutionUDot();
   u_dot = *_solution;
   computeTimeDerivativeHelper(u_dot, _solution_old, _solution_older);
-
-  for (const auto i : index_range(_du_dot_du))
-    if (integratesVar(i))
-      _du_dot_du[i] = 1. / _dt;
-  u_dot.close();
+  computeDuDotDu();
 }
 
 void

diff --git a/framework/src/timeintegrators/ExplicitSSPRungeKutta.C b/framework/src/timeintegrators/ExplicitSSPRungeKutta.C
@@ -69,9 +69,7 @@ void
 ExplicitSSPRungeKutta::computeTimeDerivatives()
 {
   // Only the Jacobian needs to be computed, since the mass matrix needs it
-  for (const auto i : index_range(_du_dot_du))
-    if (integratesVar(i))
-      _du_dot_du[i] = 1.0 / (_b[_stage] * _dt);
+  computeDuDotDu();
 }
 
 void
@@ -209,3 +207,9 @@ ExplicitSSPRungeKutta::postResidual(NumericVector<Number> & residual)
   // Set time at which to evaluate nodal BCs
   _fe_problem.time() = _current_time;
 }
+
+Real
+ExplicitSSPRungeKutta::duDotDuCoeff() const
+{
+  return Real(1) / _b[_stage];
+}
diff --git a/framework/src/timeintegrators/ExplicitTVDRK2.C b/framework/src/timeintegrators/ExplicitTVDRK2.C
@@ -56,9 +56,7 @@ ExplicitTVDRK2::computeTimeDerivatives()
   u_dot = *_solution;
   computeTimeDerivativeHelper(u_dot, _solution_old, _solution_older);
 
-  for (const auto i : index_range(_du_dot_du))
-    if (integratesVar(i))
-      _du_dot_du[i] = 1. / _dt;
+  computeDuDotDu();
   u_dot.close();
 }
 

diff --git a/framework/src/timeintegrators/ImplicitEuler.C b/framework/src/timeintegrators/ImplicitEuler.C
@@ -49,9 +49,7 @@ ImplicitEuler::computeTimeDerivatives()
     u_dot.close();
   }
 
-  for (const auto i : index_range(_du_dot_du))
-    if (integratesVar(i))
-      _du_dot_du[i] = 1.0 / _dt;
+  computeDuDotDu();
 }
 
 void

diff --git a/framework/src/timeintegrators/ImplicitMidpoint.C b/framework/src/timeintegrators/ImplicitMidpoint.C
@@ -45,9 +45,7 @@ ImplicitMidpoint::computeTimeDerivatives()
   u_dot = *_solution;
   computeTimeDerivativeHelper(u_dot, _solution_old);
   u_dot.close();
-  for (const auto i : index_range(_du_dot_du))
-    if (integratesVar(i))
-      _du_dot_du[i] = 1. / _dt;
+  computeDuDotDu();
 }
 
 void

diff --git a/framework/src/timeintegrators/LStableDirk2.C b/framework/src/timeintegrators/LStableDirk2.C
@@ -48,9 +48,7 @@ LStableDirk2::computeTimeDerivatives()
   u_dot = *_solution;
   computeTimeDerivativeHelper(u_dot, _solution_old);
   u_dot.close();
-  for (const auto i : index_range(_du_dot_du))
-    if (integratesVar(i))
-      _du_dot_du[i] = 1. / _dt;
+  computeDuDotDu();
 }
 
 void

diff --git a/framework/src/timeintegrators/LStableDirk3.C b/framework/src/timeintegrators/LStableDirk3.C
@@ -67,9 +67,7 @@ LStableDirk3::computeTimeDerivatives()
   u_dot = *_solution;
   computeTimeDerivativeHelper(u_dot, _solution_old);
   u_dot.close();
-  for (const auto i : index_range(_du_dot_du))
-    if (integratesVar(i))
-      _du_dot_du[i] = 1. / _dt;
+  computeDuDotDu();
 }
 
 void

diff --git a/framework/src/timeintegrators/LStableDirk4.C b/framework/src/timeintegrators/LStableDirk4.C
@@ -62,9 +62,7 @@ LStableDirk4::computeTimeDerivatives()
   u_dot = *_solution;
   computeTimeDerivativeHelper(u_dot, _solution_old);
   u_dot.close();
-  for (const auto i : index_range(_du_dot_du))
-    if (integratesVar(i))
-      _du_dot_du[i] = 1. / _dt;
+  computeDuDotDu();
 }
 
 void

diff --git a/framework/src/timeintegrators/NewmarkBeta.C b/framework/src/timeintegrators/NewmarkBeta.C
@@ -92,9 +92,7 @@ NewmarkBeta::computeTimeDerivatives()
 
   // used for Jacobian calculations
   _du_dotdot_du = 1.0 / _beta / _dt / _dt;
-  for (const auto i : index_range(_du_dot_du))
-    if (integratesVar(i))
-      _du_dot_du[i] = _gamma / _beta / _dt;
+  computeDuDotDu();
 }
 
 void
@@ -120,3 +118,9 @@ NewmarkBeta::postResidual(NumericVector<Number> & residual)
   residual += _Re_non_time;
   residual.close();
 }
+
+Real
+NewmarkBeta::duDotDuCoeff() const
+{
+  return _gamma / _beta;
+}