From cfb18bd51f16298fb950c6e607fee1cd16480f5a Mon Sep 17 00:00:00 2001 From: Logan Harbour Date: Mon, 18 May 2020 17:30:29 -0600 Subject: [PATCH] Use solutionState() in all solution(), solutionOld(), etc refs #13964 --- framework/include/problems/DisplacedProblem.h | 2 - framework/include/problems/FEProblemBase.h | 2 - framework/include/problems/SubProblem.h | 13 ++-- framework/include/systems/AuxiliarySystem.h | 26 ++++--- framework/include/systems/DisplacedSystem.h | 22 +++--- .../systems/DumpObjectsNonlinearSystem.h | 13 ++-- .../include/systems/NonlinearEigenSystem.h | 32 +++++---- framework/include/systems/NonlinearSystem.h | 24 ++++--- .../include/systems/NonlinearSystemBase.h | 6 +- framework/include/systems/SystemBase.h | 67 ++++++++++++++----- .../timeintegrators/CentralDifference.h | 2 + framework/src/problems/DisplacedProblem.C | 7 -- framework/src/problems/FEProblemBase.C | 14 ---- framework/src/problems/SubProblem.C | 3 +- framework/src/systems/AuxiliarySystem.C | 5 ++ framework/src/systems/DisplacedSystem.C | 13 ++-- .../src/systems/DumpObjectsNonlinearSystem.C | 1 + framework/src/systems/NonlinearEigenSystem.C | 5 ++ framework/src/systems/NonlinearSystem.C | 5 ++ framework/src/systems/NonlinearSystemBase.C | 3 +- framework/src/systems/SystemBase.C | 51 ++++++++++---- .../src/timeintegrators/CentralDifference.C | 13 ++-- 22 files changed, 201 insertions(+), 128 deletions(-) diff --git a/framework/include/problems/DisplacedProblem.h b/framework/include/problems/DisplacedProblem.h index a12b67a1805a..98f94ee9e434 100644 --- a/framework/include/problems/DisplacedProblem.h +++ b/framework/include/problems/DisplacedProblem.h @@ -308,8 +308,6 @@ class DisplacedProblem : public SubProblem const CouplingMatrix * couplingMatrix() const override; - virtual void needOldSolutionState(const unsigned int state) override; - protected: FEProblemBase & _mproblem; MooseMesh & _mesh; diff --git a/framework/include/problems/FEProblemBase.h b/framework/include/problems/FEProblemBase.h index 553c34f5a972..a1978f75e949 100644 --- a/framework/include/problems/FEProblemBase.h +++ b/framework/include/problems/FEProblemBase.h @@ -1761,8 +1761,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; diff --git a/framework/include/problems/SubProblem.h b/framework/include/problems/SubProblem.h index fc0e40b35704..ce214f583801 100644 --- a/framework/include/problems/SubProblem.h +++ b/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 @@ -694,6 +689,11 @@ class SubProblem : public Problem */ bool automaticScaling() const; + /** + * The number of default solution states each system associated with this problem should create + */ + unsigned int defaultSolutionStates() const { return _default_solution_states; } + protected: /** * Helper function called by getVariable that handles the logic for @@ -811,6 +811,9 @@ class SubProblem : public Problem bool _have_fv = false; + /// The number of default solution states that should be created + const unsigned int _default_solution_states; + ///@{ Helper functions for checking MaterialProperties std::string restrictionSubdomainCheckName(SubdomainID check_id); std::string restrictionBoundaryCheckName(BoundaryID check_id); diff --git a/framework/include/systems/AuxiliarySystem.h b/framework/include/systems/AuxiliarySystem.h index 96fe89f37944..c306cc94e68f 100644 --- a/framework/include/systems/AuxiliarySystem.h +++ b/framework/include/systems/AuxiliarySystem.h @@ -157,17 +157,7 @@ class AuxiliarySystem : public SystemBase, public PerfGraphInterface */ bool needMaterialOnSide(BoundaryID bnd_id); - NumericVector & solution() override { return *_sys.solution; } - NumericVector & solutionOld() override { return *_sys.old_local_solution; } - NumericVector & solutionOlder() override { return *_sys.older_local_solution; } NumericVector * solutionPreviousNewton() override { return _solution_previous_nl; } - - const NumericVector & solution() const override { return *_sys.solution; } - const NumericVector & solutionOld() const override { return *_sys.old_local_solution; } - const NumericVector & solutionOlder() const override - { - return *_sys.older_local_solution; - } const NumericVector * solutionPreviousNewton() const override { return _solution_previous_nl; @@ -184,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); @@ -261,4 +251,18 @@ class AuxiliarySystem : public SystemBase, public PerfGraphInterface friend class ComputeNodalKernelBcsThread; friend class ComputeNodalKernelJacobiansThread; friend class ComputeNodalKernelBCJacobiansThread; + + NumericVector & solutionInternal() override { return *_sys.solution; } + NumericVector & solutionOldInternal() override { return *_sys.old_local_solution; } + NumericVector & solutionOlderInternal() override { return *_sys.older_local_solution; } + + const NumericVector & solutionInternal() const override { return *_sys.solution; } + const NumericVector & solutionOldInternal() const override + { + return *_sys.old_local_solution; + } + const NumericVector & solutionOlderInternal() const override + { + return *_sys.older_local_solution; + } }; diff --git a/framework/include/systems/DisplacedSystem.h b/framework/include/systems/DisplacedSystem.h index e188624253c6..8e99729bde5a 100644 --- a/framework/include/systems/DisplacedSystem.h +++ b/framework/include/systems/DisplacedSystem.h @@ -89,15 +89,8 @@ class DisplacedSystem : public SystemBase return _undisplaced_system.currentSolution(); } - NumericVector & solution() override { return _undisplaced_system.solution(); } - NumericVector & solutionOld() override; - NumericVector & solutionOlder() override; - NumericVector * solutionPreviousNewton() override { return NULL; } - - const NumericVector & solution() const override { return _undisplaced_system.solution(); } - const NumericVector & solutionOld() const override; - const NumericVector & solutionOlder() const override; - const NumericVector * solutionPreviousNewton() const override { return NULL; } + NumericVector * solutionPreviousNewton() override { return nullptr; } + const NumericVector * solutionPreviousNewton() const override { return nullptr; } NumericVector * solutionUDot() override { return _undisplaced_system.solutionUDot(); } NumericVector * solutionUDotDot() override @@ -216,6 +209,17 @@ class DisplacedSystem : public SystemBase void addTimeIntegrator(std::shared_ptr ti) override; protected: + NumericVector & solutionInternal() override { return _undisplaced_system.solution(); } + NumericVector & solutionOldInternal() override; + NumericVector & solutionOlderInternal() override; + + const NumericVector & solutionInternal() const override + { + return _undisplaced_system.solution(); + } + const NumericVector & solutionOldInternal() const override; + const NumericVector & solutionOlderInternal() const override; + SystemBase & _undisplaced_system; TransientExplicitSystem & _sys; }; diff --git a/framework/include/systems/DumpObjectsNonlinearSystem.h b/framework/include/systems/DumpObjectsNonlinearSystem.h index f180ee3880a0..2b408a27eda2 100644 --- a/framework/include/systems/DumpObjectsNonlinearSystem.h +++ b/framework/include/systems/DumpObjectsNonlinearSystem.h @@ -32,16 +32,15 @@ class DumpObjectsNonlinearSystem : public NonlinearSystemBase virtual bool converged() override { return true; } virtual NumericVector & RHS() override { return *_dummy; } - NumericVector & solutionOld() override { return *_dummy; } - NumericVector & solutionOlder() override { return *_dummy; } - - const NumericVector & solutionOld() const override { return *_dummy; } - const NumericVector & solutionOlder() const override { return *_dummy; } - virtual unsigned int getCurrentNonlinearIterationNumber() override { return 0; } virtual void setupFiniteDifferencedPreconditioner() override {} protected: + NumericVector & solutionOldInternal() override { return *_dummy; } + NumericVector & solutionOlderInternal() override { return *_dummy; } + + const NumericVector & solutionOldInternal() const override { return *_dummy; } + const NumericVector & solutionOlderInternal() const override { return *_dummy; } + NumericVector * _dummy; }; - diff --git a/framework/include/systems/NonlinearEigenSystem.h b/framework/include/systems/NonlinearEigenSystem.h index 8a8929a2841b..4ac8852705a8 100644 --- a/framework/include/systems/NonlinearEigenSystem.h +++ b/framework/include/systems/NonlinearEigenSystem.h @@ -84,20 +84,6 @@ class NonlinearEigenSystem : public NonlinearSystemBase virtual NonlinearSolver * nonlinearSolver() override; - NumericVector & solutionOld() override { return *_transient_sys.old_local_solution; } - - NumericVector & solutionOlder() override { return *_transient_sys.older_local_solution; } - - const NumericVector & solutionOld() const override - { - return *_transient_sys.old_local_solution; - } - - const NumericVector & solutionOlder() const override - { - return *_transient_sys.older_local_solution; - } - virtual TransientEigenSystem & sys() { return _transient_sys; } /** @@ -147,6 +133,24 @@ class NonlinearEigenSystem : public NonlinearSystemBase TagID nonEigenMatrixTag() { return _A_tag; } protected: + NumericVector & solutionOldInternal() override + { + return *_transient_sys.old_local_solution; + } + const NumericVector & solutionOldInternal() const override + { + return *_transient_sys.old_local_solution; + } + + NumericVector & solutionOlderInternal() override + { + return *_transient_sys.older_local_solution; + } + const NumericVector & solutionOlderInternal() const override + { + return *_transient_sys.older_local_solution; + } + TransientEigenSystem & _transient_sys; EigenProblem & _eigen_problem; std::vector> _eigen_values; diff --git a/framework/include/systems/NonlinearSystem.h b/framework/include/systems/NonlinearSystem.h index 8d8e39d4bf98..9ff1cbf196dd 100644 --- a/framework/include/systems/NonlinearSystem.h +++ b/framework/include/systems/NonlinearSystem.h @@ -64,23 +64,29 @@ class NonlinearSystem : public NonlinearSystemBase return _transient_sys.nonlinear_solver.get(); } - NumericVector & solutionOld() override { return *_transient_sys.old_local_solution; } - const NumericVector & solutionOld() const override + virtual TransientNonlinearImplicitSystem & sys() { return _transient_sys; } + + void computeScaling() override; + +protected: + NumericVector & solutionOldInternal() override + { + return *_transient_sys.old_local_solution; + } + const NumericVector & solutionOldInternal() const override { return *_transient_sys.old_local_solution; } - NumericVector & solutionOlder() override { return *_transient_sys.older_local_solution; } - const NumericVector & solutionOlder() const override + NumericVector & solutionOlderInternal() override + { + return *_transient_sys.older_local_solution; + } + const NumericVector & 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; diff --git a/framework/include/systems/NonlinearSystemBase.h b/framework/include/systems/NonlinearSystemBase.h index e25d4f53588e..f3a373dd16d8 100644 --- a/framework/include/systems/NonlinearSystemBase.h +++ b/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 & solution() override { return *_sys.solution; } - const NumericVector & solution() const override { return *_sys.solution; } - virtual System & system() override { return _sys; } virtual const System & system() const override { return _sys; } @@ -719,6 +716,9 @@ class NonlinearSystemBase : public SystemBase, public PerfGraphInterface void mortarConstraints(bool displaced); protected: + NumericVector & solutionInternal() override { return *_sys.solution; } + const NumericVector & solutionInternal() const override { return *_sys.solution; } + /// solution vector from nonlinear solver const NumericVector * _current_solution; /// ghosted form of the residual diff --git a/framework/include/systems/SystemBase.h b/framework/include/systems/SystemBase.h index 4f1b536af8e0..5ebb680fa144 100644 --- a/framework/include/systems/SystemBase.h +++ b/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 @@ -212,20 +212,26 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface */ virtual const NumericVector * const & currentSolution() const = 0; - virtual NumericVector & solution() = 0; - virtual NumericVector & solutionOld() = 0; - virtual NumericVector & solutionOlder() = 0; - virtual NumericVector * solutionPreviousNewton() = 0; - virtual const NumericVector & solution() const = 0; - virtual const NumericVector & solutionOld() const = 0; - virtual const NumericVector & solutionOlder() const = 0; + NumericVector & solution() { return solutionState(0); } + NumericVector & solutionOld() { return solutionState(1); } + NumericVector & solutionOlder() { return solutionState(2); } + const NumericVector & solution() const { return solutionState(0); } + const NumericVector & solutionOld() const { return solutionState(1); } + const NumericVector & solutionOlder() const { return solutionState(2); } + virtual const NumericVector * solutionPreviousNewton() const = 0; + virtual NumericVector * solutionPreviousNewton() = 0; /** * Get a state of the solution. * State 0 is the current solution, state 1 is the old solution, etc... */ + ///@{ NumericVector & solutionState(const unsigned int state); + const NumericVector & solutionState(const unsigned int state) const; + ///@} + + NumericVector & addSolutionState(const unsigned int state); virtual Number & duDotDu() { return _du_dot_du; } virtual Number & duDotDotDu() { return _du_dotdot_du; } @@ -812,12 +818,36 @@ 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 & vars); +protected: + /** + * Internal getters for the libMesh solution, used only in filling _solution_states + */ + ///@{ + virtual NumericVector & solutionInternal() = 0; + virtual const NumericVector & solutionInternal() const = 0; + ///@} + + /** + * Internal getters for the old libMesh solution, used only in filling _solution_states + */ + ///@{ + virtual NumericVector & solutionOldInternal() = 0; + virtual const NumericVector & solutionOldInternal() const = 0; + ///@} + + /** + * Internal getters for the older libMesh solution, used only in filling _solution_states + */ + ///@{ + virtual NumericVector & solutionOlderInternal() = 0; + virtual const NumericVector & solutionOlderInternal() const = 0; + ///@} + /** * Register that a certain solution state is needed (1 = old, 2 = older, etc) */ void needOldSolutionState(const unsigned int state); - -protected: + SubProblem & _subproblem; MooseApp & _app; @@ -872,13 +902,6 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface /// Map variable number to its pointer std::vector> _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 *> _solution_states; - /// The saved solution states (0 = current, 1 = old, 2 = older, etc) - std::vector *> _saved_solution_states; - /// Storage for MooseVariable objects MooseObjectWarehouseBase _variable_warehouse; @@ -894,11 +917,21 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface /// True if printing out additional information bool _verbose; + /// Whether or not init() has been called + bool _initialized; + private: /* * Sets up the solution states - stores older copies of the solutions in _solution_states) */ void setupSolutionStates(); + + /// 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 *> _solution_states; + /// The saved solution states (0 = current, 1 = old, 2 = older, etc) + std::vector *> _saved_solution_states; }; #define PARALLEL_TRY diff --git a/framework/include/timeintegrators/CentralDifference.h b/framework/include/timeintegrators/CentralDifference.h index 452811b78cc5..fda6150cc694 100644 --- a/framework/include/timeintegrators/CentralDifference.h +++ b/framework/include/timeintegrators/CentralDifference.h @@ -36,6 +36,8 @@ class CentralDifference : public ActuallyExplicitEuler /// solution vector for \f$ {du^dotdot}\over{du} \f$ Real & _du_dotdot_du; + const NumericVector & _solution_old_old_old; + /** * Helper function that actually does the math for computing the time derivative */ diff --git a/framework/src/problems/DisplacedProblem.C b/framework/src/problems/DisplacedProblem.C index bb2500c5415a..d8454139dbcf 100644 --- a/framework/src/problems/DisplacedProblem.C +++ b/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() { diff --git a/framework/src/problems/FEProblemBase.C b/framework/src/problems/FEProblemBase.C index 8dbe5f5cbfe8..4a45b6c4f3ec 100644 --- a/framework/src/problems/FEProblemBase.C +++ b/framework/src/problems/FEProblemBase.C @@ -4748,10 +4748,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(); @@ -6730,16 +6726,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) { diff --git a/framework/src/problems/SubProblem.C b/framework/src/problems/SubProblem.C index 17a650dd4ed7..d478068ae0eb 100644 --- a/framework/src/problems/SubProblem.C +++ b/framework/src/problems/SubProblem.C @@ -52,7 +52,8 @@ SubProblem::SubProblem(const InputParameters & parameters) _safe_access_tagged_matrices(false), _safe_access_tagged_vectors(false), _have_ad_objects(false), - _typed_vector_tags(2) + _typed_vector_tags(2), + _default_solution_states(2) { unsigned int n_threads = libMesh::n_threads(); _active_elemental_moose_variables.resize(n_threads); diff --git a/framework/src/systems/AuxiliarySystem.C b/framework/src/systems/AuxiliarySystem.C index e3ec9d894f90..55f47b67a97a 100644 --- a/framework/src/systems/AuxiliarySystem.C +++ b/framework/src/systems/AuxiliarySystem.C @@ -68,6 +68,11 @@ AuxiliarySystem::AuxiliarySystem(FEProblemBase & subproblem, const std::string & dof_map.remove_algebraic_ghosting_functor(dof_map.default_algebraic_ghosting()); dof_map.set_implicit_neighbor_dofs(false); } + + // By default, we have always had current, old, and older solutions available. With the transition + // to storing these vectors for generalized access in _solution_states in each system, we need to + // setup these default solution states for the objects that request them in their constructors + needOldSolutionState(subproblem.defaultSolutionStates()); } AuxiliarySystem::~AuxiliarySystem() { delete &_serialized_solution; } diff --git a/framework/src/systems/DisplacedSystem.C b/framework/src/systems/DisplacedSystem.C index 9025af6cb3b7..e734c60cd95f 100644 --- a/framework/src/systems/DisplacedSystem.C +++ b/framework/src/systems/DisplacedSystem.C @@ -29,6 +29,11 @@ DisplacedSystem::DisplacedSystem(DisplacedProblem & problem, dof_map.remove_algebraic_ghosting_functor(dof_map.default_algebraic_ghosting()); dof_map.set_implicit_neighbor_dofs(false); } + + // By default, we have always had current, old, and older solutions available. With the transition + // to storing these vectors for generalized access in _solution_states in each system, we need to + // setup these default solution states for the objects that request them in their constructors + needOldSolutionState(problem.defaultSolutionStates()); } DisplacedSystem::~DisplacedSystem() {} @@ -49,25 +54,25 @@ DisplacedSystem::addTimeIntegrator(std::shared_ptr ti) } NumericVector & -DisplacedSystem::solutionOld() +DisplacedSystem::solutionOldInternal() { return *_sys.old_local_solution; } NumericVector & -DisplacedSystem::solutionOlder() +DisplacedSystem::solutionOlderInternal() { return *_sys.older_local_solution; } const NumericVector & -DisplacedSystem::solutionOld() const +DisplacedSystem::solutionOldInternal() const { return *_sys.old_local_solution; } const NumericVector & -DisplacedSystem::solutionOlder() const +DisplacedSystem::solutionOlderInternal() const { return *_sys.older_local_solution; } diff --git a/framework/src/systems/DumpObjectsNonlinearSystem.C b/framework/src/systems/DumpObjectsNonlinearSystem.C index 365e48e7ced2..cc27f74ec853 100644 --- a/framework/src/systems/DumpObjectsNonlinearSystem.C +++ b/framework/src/systems/DumpObjectsNonlinearSystem.C @@ -16,4 +16,5 @@ DumpObjectsNonlinearSystem::DumpObjectsNonlinearSystem(FEProblemBase & problem, problem, problem.es().add_system(name), name), _dummy(NULL) { + needOldSolutionState(problem.defaultSolutionStates()); } diff --git a/framework/src/systems/NonlinearEigenSystem.C b/framework/src/systems/NonlinearEigenSystem.C index 1eb361dcf72b..10a49de00473 100644 --- a/framework/src/systems/NonlinearEigenSystem.C +++ b/framework/src/systems/NonlinearEigenSystem.C @@ -142,6 +142,11 @@ NonlinearEigenSystem::NonlinearEigenSystem(EigenProblem & eigen_problem, const s _A_tag = eigen_problem.addMatrixTag("A_tag"); _B_tag = eigen_problem.addMatrixTag("Eigen"); + + // By default, we have always had current, old, and older solutions available. With the transition + // to storing these vectors for generalized access in _solution_states in each system, we need to + // setup these default solution states for the objects that request them in their constructors + needOldSolutionState(eigen_problem.defaultSolutionStates()); } void diff --git a/framework/src/systems/NonlinearSystem.C b/framework/src/systems/NonlinearSystem.C index 7d8d8228d703..9b2dbf8f8e5a 100644 --- a/framework/src/systems/NonlinearSystem.C +++ b/framework/src/systems/NonlinearSystem.C @@ -115,6 +115,11 @@ NonlinearSystem::NonlinearSystem(FEProblemBase & fe_problem, const std::string & petsc_solver->use_default_monitor(false); } #endif + + // By default, we have always had current, old, and older solutions available. With the transition + // to storing these vectors for generalized access in _solution_states in each system, we need to + // setup these default solution states for the objects that request them in their constructors + needOldSolutionState(fe_problem.defaultSolutionStates()); } NonlinearSystem::~NonlinearSystem() {} diff --git a/framework/src/systems/NonlinearSystemBase.C b/framework/src/systems/NonlinearSystemBase.C index df3e0dff44f3..000445ecb406 100644 --- a/framework/src/systems/NonlinearSystemBase.C +++ b/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) diff --git a/framework/src/systems/SystemBase.C b/framework/src/systems/SystemBase.C index e2e0af73a20a..af629c72a653 100644 --- a/framework/src/systems/SystemBase.C +++ b/framework/src/systems/SystemBase.C @@ -104,37 +104,33 @@ 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), - _verbose(false) + _verbose(false), + _old_solution_states_needed(0) { } -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(); + // 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 (_old_solution_states_needed > 0) - _solution_states[1] = &solutionOld(); + _solution_states[1] = &solutionOldInternal(); if (_old_solution_states_needed > 1) - _solution_states[2] = &solutionOlder(); + _solution_states[2] = &solutionOlderInternal(); // 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); + if (!_solution_states[i]) + _solution_states[i] = &addVector("solution_state_" + std::to_string(i), true, GHOSTED); } MooseVariableFEBase & @@ -1205,7 +1201,10 @@ void SystemBase::needOldSolutionState(const unsigned int state) { if (_old_solution_states_needed < state) + { _old_solution_states_needed = state; + setupSolutionStates(); + } } NumericVector & @@ -1218,9 +1217,33 @@ SystemBase::solutionState(const unsigned int state) _old_solution_states_needed, " is available.\nUse the needOldSolutionState() method in the Problem to request " "older solution states."); + + mooseAssert(_solution_states.size() > state, "_solution_state size mismatch"); + return *_solution_states[state]; +} + +const NumericVector & +SystemBase::solutionState(const unsigned int state) const +{ + 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."); + + mooseAssert(_solution_states.size() > state, "_solution_state size mismatch"); return *_solution_states[state]; } +NumericVector & +SystemBase::addSolutionState(const unsigned int state) +{ + needOldSolutionState(state); + return solutionState(state); +} + void SystemBase::applyScalingFactors(const std::vector & inverse_scaling_factors) { diff --git a/framework/src/timeintegrators/CentralDifference.C b/framework/src/timeintegrators/CentralDifference.C index 98c4d072d846..87db6d2b7f5c 100644 --- a/framework/src/timeintegrators/CentralDifference.C +++ b/framework/src/timeintegrators/CentralDifference.C @@ -30,7 +30,9 @@ validParams() } CentralDifference::CentralDifference(const InputParameters & parameters) - : ActuallyExplicitEuler(parameters), _du_dotdot_du(_sys.duDotDotDu()) + : ActuallyExplicitEuler(parameters), + _du_dotdot_du(_sys.duDotDotDu()), + _solution_old_old_old(_sys.addSolutionState(3)) { _is_explicit = true; if (_solve_type == LUMPED) @@ -39,19 +41,16 @@ CentralDifference::CentralDifference(const InputParameters & parameters) _fe_problem.setUDotOldRequested(true); _fe_problem.setUDotDotRequested(true); _fe_problem.setUDotDotOldRequested(true); - _fe_problem.needOldSolutionState(3); } void CentralDifference::computeADTimeDerivatives(DualReal & ad_u_dot, const dof_id_type & dof) const { - const auto & u_old_old_old = (_sys.solutionState(3))(dof); - auto u_dotdot = ad_u_dot; // TODO: Ask Alex if this has to be ad_u_dotdot. Currently, it will // change u_dotdot when called. computeTimeDerivativeHelper( - ad_u_dot, u_dotdot, _solution_old(dof), _solution_older(dof), u_old_old_old); + ad_u_dot, u_dotdot, _solution_old(dof), _solution_older(dof), _solution_old_old_old(dof)); } void @@ -78,10 +77,10 @@ CentralDifference::computeTimeDerivatives() // Declaring u_dot and u_dotdot auto & u_dot = *_sys.solutionUDot(); auto & u_dotdot = *_sys.solutionUDotDot(); - const auto & u_old_old_old = _sys.solutionState(3); // Computing derivatives - computeTimeDerivativeHelper(u_dot, u_dotdot, _solution_old, _solution_older, u_old_old_old); + computeTimeDerivativeHelper( + u_dot, u_dotdot, _solution_old, _solution_older, _solution_old_old_old); // make sure _u_dotdot and _u_dot are in good state u_dotdot.close();