Use solutionState() in all solution(), solutionOld(), etc

refs idaholab#13964

framework/include/problems/DisplacedProblem.h

@@ -315,8 +315,6 @@ class DisplacedProblem : public SubProblem
 
   const CouplingMatrix * couplingMatrix() const override;
 
-  virtual void needOldSolutionState(const unsigned int state) override;
-
 protected:
   FEProblemBase & _mproblem;
   MooseMesh & _mesh;

framework/include/problems/FEProblemBase.h

@@ -1770,8 +1770,6 @@ class FEProblemBase : public SubProblem, public Restartable
    */
   void uniformRefine();
 
-  virtual void needOldSolutionState(const unsigned int state) override;
-
   using SubProblem::automaticScaling;
   void automaticScaling(bool automatic_scaling) override;
 

framework/include/problems/SubProblem.h

@@ -677,11 +677,6 @@ class SubProblem : public Problem
    */
   void addAlgebraicGhostingFunctor(GhostingFunctor & algebraic_gf, bool to_mesh = true);
 
-  /**
-   * Register with the Systems that a certain solution state is needed (1 = old, 2 = older, etc)
-   */
-  virtual void needOldSolutionState(const unsigned int state) = 0;
-
   /**
    * Automatic scaling setter
    * @param automatic_scaling A boolean representing whether we are performing automatic scaling

framework/include/systems/SystemBase.h

@@ -180,7 +180,7 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface
   /**
    * Initialize the system
    */
-  virtual void init();
+  virtual void init(){};
 
   /**
    * Called only once, just before the solve begins so objects can do some precalculations
@@ -239,11 +239,21 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface
   const NumericVector<Number> & solutionState(const unsigned int state) const;
 
   /**
-   * Get a state of the solution.
-   * State 0 is the current solution, state 1 is the old solution, etc...
+   * 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; }
   virtual const Number & duDotDu() const { return _du_dot_du; }
@@ -829,10 +839,20 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface
   /// caches the dof indices of provided variables in MooseMesh's FaceInfo data structure
   void cacheVarIndicesByFace(const std::vector<VariableName> & vars);
 
+protected:
   /**
-   * Register that a certain solution state is needed (1 = old, 2 = older, etc)
+   * 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.
    */
-  void needOldSolutionState(const unsigned int state);
+  ///@{
+  virtual NumericVector<Number> & solutionInternal() const = 0;
+  virtual NumericVector<Number> & solutionOldInternal() const = 0;
+  virtual NumericVector<Number> & solutionOlderInternal() const = 0;
+  ///@}
 
 protected:
   /**
@@ -903,13 +923,6 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface
   /// Map variable number to its pointer
   std::vector<std::vector<MooseVariableFEBase *>> _numbered_vars;
 
-  /// How many old solutions states are needed (1 = old, 2 = older, etc)
-  unsigned int _old_solution_states_needed;
-  /// 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<Real> *> _saved_solution_states;
-
   /// Storage for MooseVariable objects
   MooseObjectWarehouseBase<MooseVariableBase> _variable_warehouse;
 

framework/src/problems/DisplacedProblem.C

@@ -139,13 +139,6 @@ DisplacedProblem::initAdaptivity()
 {
 }
 
-void
-DisplacedProblem::needOldSolutionState(const unsigned int state)
-{
-  _displaced_nl.needOldSolutionState(state);
-  _displaced_aux.needOldSolutionState(state);
-}
-
 void
 DisplacedProblem::saveOldSolutions()
 {

framework/src/problems/FEProblemBase.C

@@ -4786,10 +4786,6 @@ FEProblemBase::init()
   if (solverParams()._type == Moose::ST_JFNK)
     _nl->turnOffJacobian();
 
-  // Request at least 2 old solution states (old and older)
-  // These get initialized in the init() calls for the systems that follow
-  needOldSolutionState(2);
-
   _nl->init();
   _aux->init();
 
@@ -6768,16 +6764,6 @@ FEProblemBase::uniformRefine()
   meshChangedHelper(/*intermediate_change=*/false);
 }
 
-void
-FEProblemBase::needOldSolutionState(const unsigned int state)
-{
-  _nl->needOldSolutionState(state);
-  _aux->needOldSolutionState(state);
-
-  if (_displaced_problem)
-    _displaced_problem->needOldSolutionState(state);
-}
-
 void
 FEProblemBase::automaticScaling(bool automatic_scaling)
 {

framework/src/systems/NonlinearSystemBase.C

@@ -173,9 +173,8 @@ NonlinearSystemBase::NonlinearSystemBase(FEProblemBase & fe_problem,
     _resid_vs_jac_scaling_param(0)
 #ifndef MOOSE_SPARSE_AD
     ,
-    _required_derivative_size(0),
+    _required_derivative_size(0)
 #endif
-    _solution_state(0)
 {
   getResidualNonTimeVector();
   // Don't need to add the matrix - it already exists (for now)

framework/src/systems/SystemBase.C

@@ -104,7 +104,6 @@ SystemBase::SystemBase(SubProblem & subproblem,
     _max_var_n_dofs_per_elem(0),
     _max_var_n_dofs_per_node(0),
     _time_integrator(nullptr),
-    _old_solution_states_needed(0),
     _computing_scaling_jacobian(false),
     _computing_scaling_residual(false),
     _automatic_scaling(false),
@@ -113,31 +112,6 @@ SystemBase::SystemBase(SubProblem & subproblem,
 {
 }
 
-void
-SystemBase::init()
-{
-  // At this point, we know how many solution states we need
-  setupSolutionStates();
-}
-
-void
-SystemBase::setupSolutionStates()
-{
-  // We'll always have 1 state for the current solution
-  _solution_states.resize(_old_solution_states_needed + 1);
-
-  // The first three states (now, old, older) will point to the solutions in the libMesh system
-  _solution_states[0] = &solution();
-  if (_old_solution_states_needed > 0)
-    _solution_states[1] = &solutionOld();
-  if (_old_solution_states_needed > 1)
-    _solution_states[2] = &solutionOlder();
-
-  // Create anything that is past older (state of 3+)
-  for (unsigned int i = 3; i <= _old_solution_states_needed; ++i)
-    _solution_states[i] = &addVector("solution_state_" + std::to_string(i), true, GHOSTED);
-}
-
 MooseVariableFEBase &
 SystemBase::getVariable(THREAD_ID tid, const std::string & var_name)
 {
@@ -1260,20 +1234,45 @@ SystemBase::solutionState(const unsigned int state)
 void
 SystemBase::needOldSolutionState(const unsigned int state)
 {
-  if (_old_solution_states_needed < state)
-    _old_solution_states_needed = state;
+  mooseAssert(
+      !_solution_states.empty(),
+      "No solution states available: make sure to init the default states in system constructors");
+
+  if (state >= _solution_states.size())
+    mooseError("Solution state ",
+               state,
+               " was requested in ",
+               name(),
+               " but only up to state ",
+               _solution_states.size() - 1,
+               " is available.");
+
+  return *_solution_states[state];
 }
 
 NumericVector<Number> &
 SystemBase::solutionState(const unsigned int state)
 {
-  if (state > _old_solution_states_needed)
-    mooseError("Solution state ",
-               state,
-               " was requested, but only up to state ",
-               _old_solution_states_needed,
-               " is available.\nUse the needOldSolutionState() method in the Problem to request "
-               "older solution states.");
+  // Create up to the state requested if unavailable
+  if (state >= _solution_states.size())
+  {
+    _solution_states.resize(state + 1);
+
+    // The first three states (now, old, older) will point to the solutions in the libMesh system,
+    // which is why we are using the "internal" calls to these vectors. _solution_states will then
+    // be the forward facing access to these vectors
+    _solution_states[0] = &solutionInternal();
+    if (state > 0)
+      _solution_states[1] = &solutionOldInternal();
+    if (state > 1)
+      _solution_states[2] = &solutionOlderInternal();
+
+    // Create anything that is past older (state of 3+)
+    for (unsigned int i = 3; i <= state; ++i)
+      if (!_solution_states[i])
+        _solution_states[i] = &addVector("solution_state_" + std::to_string(i), true, GHOSTED);
+  }
+
   return *_solution_states[state];
 }