From 328936bc959dbb385510619afa5fec9b8c8fba4f Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Thu, 9 May 2024 06:20:20 -0700 Subject: [PATCH 01/16] Add test --- .../multiple-integrators/gold/test_out.csv | 5 + .../multiple-integrators/test.i | 102 ++++++++++++++++++ .../multiple-integrators/test.py | 37 +++++++ .../multiple-integrators/tests | 26 +++++ 4 files changed, 170 insertions(+) create mode 100644 test/tests/time_integrators/multiple-integrators/gold/test_out.csv create mode 100644 test/tests/time_integrators/multiple-integrators/test.i create mode 100755 test/tests/time_integrators/multiple-integrators/test.py create mode 100644 test/tests/time_integrators/multiple-integrators/tests diff --git a/test/tests/time_integrators/multiple-integrators/gold/test_out.csv b/test/tests/time_integrators/multiple-integrators/gold/test_out.csv new file mode 100644 index 000000000000..0d82e2945d94 --- /dev/null +++ b/test/tests/time_integrators/multiple-integrators/gold/test_out.csv @@ -0,0 +1,5 @@ +time,L2u,L2v +0,0,0 +1,0.0097190858385144,0.02031363438757 +2,0.0017026874015291,0.051701555106182 +3,0.0083496354351904,0.08359162599358 diff --git a/test/tests/time_integrators/multiple-integrators/test.i b/test/tests/time_integrators/multiple-integrators/test.i new file mode 100644 index 000000000000..ce63ecc158b9 --- /dev/null +++ b/test/tests/time_integrators/multiple-integrators/test.i @@ -0,0 +1,102 @@ +[Mesh] + type = GeneratedMesh + dim = 2 + nx = 10 + ny = 10 +[] + +[Variables] + [u][] + [v][] +[] + +[Kernels] + [timeu] + type = TimeDerivative + variable = u + [] + [timev] + type = TimeDerivative + variable = v + [] + [diffu] + type = Diffusion + variable = u + [] + [diffv] + type = Diffusion + variable = v + [] + [forceu] + type = BodyForce + variable = u + function = force + [] + [forcev] + type = BodyForce + variable = v + function = force + [] +[] + +[Functions] + [exact] + type = ParsedFunction + expression = 't^3*x*y' + [] + [force] + type = ParsedFunction + expression = '3*x*y*t^2' + [] +[] + +[BCs] + [allu] + type = FunctionDirichletBC + variable = u + function = exact + boundary = 'left right top bottom' + [] + [allv] + type = FunctionDirichletBC + variable = v + function = exact + boundary = 'left right top bottom' + [] +[] + +[Executioner] + type = Transient + solve_type = 'NEWTON' + petsc_options_iname = '-pc_type' + petsc_options_value = 'hypre' + dt = 1 + end_time = 3 + [TimeIntegrators] + [cn] + type = CrankNicolson + variables = 'u' + [] + [ie] + type = ImplicitEuler + variables = 'v' + [] + [] +[] + +[Postprocessors] + [L2u] + type = ElementL2Error + function = exact + variable = u + [] + [L2v] + type = ElementL2Error + function = exact + variable = v + [] +[] + +[Outputs] + csv = true +[] diff --git a/test/tests/time_integrators/multiple-integrators/test.py b/test/tests/time_integrators/multiple-integrators/test.py new file mode 100755 index 000000000000..bcac6a9eedfa --- /dev/null +++ b/test/tests/time_integrators/multiple-integrators/test.py @@ -0,0 +1,37 @@ +import mms +import unittest +from mooseutils import fuzzyAbsoluteEqual + +class TestMultipleTimeIntegrators(unittest.TestCase): + def test(self): + df1 = mms.run_temporal('test.i', 5, y_pp=['L2u', 'L2v']) + fig = mms.ConvergencePlot(xlabel=r'$\Delta$t', ylabel='$L_2$ Error') + fig.plot(df1, + label=['Crank Nicolson', 'Implicit Euler'], + marker='o', + markersize=8, + slope_precision=1, + num_fitted_points=3) + fig.save('mms_temporal.png') + for label,value in fig.label_to_slope.items(): + if label == 'Implicit Euler': + self.assertTrue(fuzzyAbsoluteEqual(value, 1., .1)) + else: + self.assertTrue(fuzzyAbsoluteEqual(value, 2., .1)) + +class TestMultipleTimeIntegratorsParallel(unittest.TestCase): + def test(self): + df1 = mms.run_temporal('test.i', 5, y_pp=['L2u', 'L2v'], mpi=2) + fig = mms.ConvergencePlot(xlabel=r'$\Delta$t', ylabel='$L_2$ Error') + fig.plot(df1, + label=['Crank Nicolson', 'Implicit Euler'], + marker='o', + markersize=8, + slope_precision=1, + num_fitted_points=3) + fig.save('mms_temporal.png') + for label,value in fig.label_to_slope.items(): + if label == 'Implicit Euler': + self.assertTrue(fuzzyAbsoluteEqual(value, 1., .1)) + else: + self.assertTrue(fuzzyAbsoluteEqual(value, 2., .1)) diff --git a/test/tests/time_integrators/multiple-integrators/tests b/test/tests/time_integrators/multiple-integrators/tests new file mode 100644 index 000000000000..5b21826a9be3 --- /dev/null +++ b/test/tests/time_integrators/multiple-integrators/tests @@ -0,0 +1,26 @@ +[Tests] + issues = '#19228' + design = 'ImplicitEuler.md CrankNicolson.md' + [test] + requirement = 'The system shall support having multiple time integrators in the same input file and yield the same results across thread counts, process counts, and mesh modes.' + type = CSVDiff + input = test.i + csvdiff = test_out.csv + [] + [test_accuracy] + requirement = 'The system shall support having multiple time integrators in the same input file and yield the expected temporal convergence rates for each variable.' + type = PythonUnitTest + input = test.py + test_case = TestMultipleTimeIntegrators + required_python_packages = 'pandas matplotlib' + installation_type = in_tree # see #26480 + prereq = test + [] + [jac] + type = PetscJacobianTester + input = test.i + cli_args = 'Outputs/csv=false' + run_sim = True + requirement = 'The system shall have an exact Jacobian when running with multiple time integrators.' + [] +[] From f399aa9ec14d2a7e1cfbe17889b256765a671dc5 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Tue, 21 May 2024 15:11:10 -0700 Subject: [PATCH 02/16] Plumbing for multiple time integrators Closes #19228 --- framework/include/executioners/Executioner.h | 2 +- framework/include/executioners/Transient.h | 2 +- framework/include/systems/AuxiliarySystem.h | 10 -- framework/include/systems/DisplacedSystem.h | 3 - framework/include/systems/LinearSystem.h | 5 - .../include/systems/NonlinearSystemBase.h | 11 --- framework/include/systems/SystemBase.h | 39 +++++--- .../include/timeintegrators/AStableDirk4.h | 1 + .../timeintegrators/ActuallyExplicitEuler.h | 1 + .../include/timeintegrators/ExplicitRK2.h | 1 + .../timeintegrators/ExplicitSSPRungeKutta.h | 1 + .../include/timeintegrators/ExplicitTVDRK2.h | 1 + .../timeintegrators/ImplicitMidpoint.h | 1 + .../include/timeintegrators/LStableDirk2.h | 1 + .../include/timeintegrators/LStableDirk3.h | 1 + .../include/timeintegrators/LStableDirk4.h | 1 + .../include/timeintegrators/TimeIntegrator.h | 33 ++++++- .../src/actions/SetupTimeIntegratorAction.C | 11 ++- framework/src/base/Moose.C | 8 +- framework/src/executioners/Transient.C | 24 +++-- framework/src/outputs/ConsoleUtils.C | 7 +- framework/src/problems/DisplacedProblem.C | 6 +- framework/src/problems/FEProblemBase.C | 10 +- framework/src/systems/AuxiliarySystem.C | 25 ++--- framework/src/systems/DisplacedSystem.C | 6 -- framework/src/systems/LinearSystem.C | 8 -- framework/src/systems/NonlinearEigenSystem.C | 23 +++-- framework/src/systems/NonlinearSystem.C | 30 ++++-- framework/src/systems/NonlinearSystemBase.C | 29 +++--- framework/src/systems/SolverSystem.C | 13 +-- framework/src/systems/SystemBase.C | 44 ++++++++- framework/src/timeintegrators/CrankNicolson.C | 32 +++++-- framework/src/timeintegrators/ImplicitEuler.C | 20 +++- .../src/timeintegrators/TimeIntegrator.C | 91 +++++++++++++++++-- .../src/timesteppers/AB2PredictorCorrector.C | 6 +- framework/src/variables/MooseVariableData.C | 2 +- framework/src/variables/MooseVariableDataFV.C | 2 +- .../src/variables/MooseVariableDataLinearFV.C | 2 +- framework/src/variables/MooseVariableField.C | 2 +- .../src/bcs/INSADDisplaceBoundaryBC.C | 2 +- 40 files changed, 357 insertions(+), 160 deletions(-) diff --git a/framework/include/executioners/Executioner.h b/framework/include/executioners/Executioner.h index c15238037f04..cf04fc526c67 100644 --- a/framework/include/executioners/Executioner.h +++ b/framework/include/executioners/Executioner.h @@ -108,7 +108,7 @@ class Executioner : public MooseObject, * This is an empty string for non-Transient executioners * @return A string of giving the TimeIntegrator name */ - virtual std::string getTimeIntegratorName() const { return std::string(); } + virtual std::vector getTimeIntegratorNames() const { return {}; } /** * Can be used by subclasses to call parentOutputPositionChanged() diff --git a/framework/include/executioners/Transient.h b/framework/include/executioners/Transient.h index 14d03f2ab419..98911afaa3f0 100644 --- a/framework/include/executioners/Transient.h +++ b/framework/include/executioners/Transient.h @@ -138,7 +138,7 @@ class Transient : public Executioner * Get the name of the time integrator (time integration scheme) used * @return string with the time integration scheme name */ - virtual std::string getTimeIntegratorName() const override; + virtual std::vector getTimeIntegratorNames() const override; /** * Get the time scheme used diff --git a/framework/include/systems/AuxiliarySystem.h b/framework/include/systems/AuxiliarySystem.h index dae30d4f2a11..38c693c44d34 100644 --- a/framework/include/systems/AuxiliarySystem.h +++ b/framework/include/systems/AuxiliarySystem.h @@ -55,16 +55,6 @@ class AuxiliarySystem : public SystemBase, public PerfGraphInterface virtual void addVariable(const std::string & var_type, const std::string & name, InputParameters & parameters) override; - /** - * Add a time integrator - * @param type Type of the integrator - * @param name The name of the integrator - * @param parameters Integrator params - */ - void addTimeIntegrator(const std::string & type, - const std::string & name, - InputParameters & parameters) override; - using SystemBase::addTimeIntegrator; /** * Adds an auxiliary kernel diff --git a/framework/include/systems/DisplacedSystem.h b/framework/include/systems/DisplacedSystem.h index 32b62ad16430..436e97404a28 100644 --- a/framework/include/systems/DisplacedSystem.h +++ b/framework/include/systems/DisplacedSystem.h @@ -249,9 +249,6 @@ class DisplacedSystem : public SystemBase virtual System & system() override; virtual const System & system() const override; - using SystemBase::addTimeIntegrator; - void addTimeIntegrator(std::shared_ptr ti) override; - virtual void compute(ExecFlagType) override {} protected: diff --git a/framework/include/systems/LinearSystem.h b/framework/include/systems/LinearSystem.h index 3fe3f4fd86e8..5c8fa562c35c 100644 --- a/framework/include/systems/LinearSystem.h +++ b/framework/include/systems/LinearSystem.h @@ -40,11 +40,6 @@ class LinearSystem : public SolverSystem, public PerfGraphInterface virtual void solve() override; - virtual void addTimeIntegrator(const std::string & type, - const std::string & name, - InputParameters & parameters) override; - using SystemBase::addTimeIntegrator; - virtual void initialSetup() override; // Overriding these to make sure the linear systems don't do anything during diff --git a/framework/include/systems/NonlinearSystemBase.h b/framework/include/systems/NonlinearSystemBase.h index 43b901581532..fb705d00c9f1 100644 --- a/framework/include/systems/NonlinearSystemBase.h +++ b/framework/include/systems/NonlinearSystemBase.h @@ -115,17 +115,6 @@ class NonlinearSystemBase : public SolverSystem, public PerfGraphInterface */ virtual bool converged() = 0; - /** - * Add a time integrator - * @param type Type of the integrator - * @param name The name of the integrator - * @param parameters Integrator params - */ - void addTimeIntegrator(const std::string & type, - const std::string & name, - InputParameters & parameters) override; - using SystemBase::addTimeIntegrator; - /** * Adds a kernel * @param kernel_name The type of the kernel diff --git a/framework/include/systems/SystemBase.h b/framework/include/systems/SystemBase.h index 985a61acb915..919dc6349184 100644 --- a/framework/include/systems/SystemBase.h +++ b/framework/include/systems/SystemBase.h @@ -878,18 +878,11 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface */ virtual void copySolutionsBackwards(); - virtual void addTimeIntegrator(const std::string & /*type*/, - const std::string & /*name*/, - InputParameters & /*parameters*/) - { - } - - virtual void addTimeIntegrator(std::shared_ptr /*ti*/) {} - - TimeIntegrator * getTimeIntegrator() { return _time_integrator.get(); } - const TimeIntegrator * getTimeIntegrator() const { return _time_integrator.get(); } + void addTimeIntegrator(const std::string & type, + const std::string & name, + InputParameters & parameters); - std::shared_ptr getSharedTimeIntegrator() { return _time_integrator; } + // virtual void addTimeIntegrator(std::shared_ptr /*ti*/) {} /// Whether or not there are variables to be restarted from an Exodus mesh file bool hasVarCopy() const { return _var_to_copy.size() > 0; } @@ -950,6 +943,28 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface */ virtual void compute(ExecFlagType type) = 0; + /** + * Copy time integrators from another system + */ + void copyTimeIntegrators(const SystemBase & other_sys); + + /** + * Retrieve the time integrator that integrates the given variable's equation + */ + const TimeIntegrator & getTimeIntegrator(const unsigned int var_num) const; + + /** + * Retrieve the time integrator that integrates the given variable's equation. If no suitable time + * 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; + + /** + * @returns All the time integrators owned by this system + */ + const std::vector> & getTimeIntegrators(); + protected: /** * Internal getter for solution owned by libMesh. @@ -1020,7 +1035,7 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface size_t _max_var_n_dofs_per_node; /// Time integrator - std::shared_ptr _time_integrator; + std::vector> _time_integrators; /// Map variable number to its pointer std::vector> _numbered_vars; diff --git a/framework/include/timeintegrators/AStableDirk4.h b/framework/include/timeintegrators/AStableDirk4.h index 29ba0474182a..4fb8ba69216c 100644 --- a/framework/include/timeintegrators/AStableDirk4.h +++ b/framework/include/timeintegrators/AStableDirk4.h @@ -62,6 +62,7 @@ class AStableDirk4 : public TimeIntegrator ADReal & ad_u_dotdot) const override; virtual void solve() override; virtual void postResidual(NumericVector & residual) override; + virtual bool overridesSolve() const override { return true; } protected: /** diff --git a/framework/include/timeintegrators/ActuallyExplicitEuler.h b/framework/include/timeintegrators/ActuallyExplicitEuler.h index 6620aaa9f1d7..03182e6312b2 100644 --- a/framework/include/timeintegrators/ActuallyExplicitEuler.h +++ b/framework/include/timeintegrators/ActuallyExplicitEuler.h @@ -29,6 +29,7 @@ class ActuallyExplicitEuler : public ExplicitTimeIntegrator ADReal & ad_u_dotdot) const override; virtual void solve() override; virtual void postResidual(NumericVector & residual) override; + virtual bool overridesSolve() const override { return true; } protected: /** diff --git a/framework/include/timeintegrators/ExplicitRK2.h b/framework/include/timeintegrators/ExplicitRK2.h index 490e916d802f..d4ed2aaed04c 100644 --- a/framework/include/timeintegrators/ExplicitRK2.h +++ b/framework/include/timeintegrators/ExplicitRK2.h @@ -67,6 +67,7 @@ class ExplicitRK2 : public TimeIntegrator ADReal & ad_u_dotdot) const override; virtual void solve() override; virtual void postResidual(NumericVector & residual) override; + virtual bool overridesSolve() const override { return true; } protected: /** diff --git a/framework/include/timeintegrators/ExplicitSSPRungeKutta.h b/framework/include/timeintegrators/ExplicitSSPRungeKutta.h index 8c7aa888bc12..9423eb76acb3 100644 --- a/framework/include/timeintegrators/ExplicitSSPRungeKutta.h +++ b/framework/include/timeintegrators/ExplicitSSPRungeKutta.h @@ -28,6 +28,7 @@ class ExplicitSSPRungeKutta : public ExplicitTimeIntegrator virtual void solve() override; virtual void postResidual(NumericVector & residual) override; virtual int order() override { return _order; } + virtual bool overridesSolve() const override { return true; } protected: /** diff --git a/framework/include/timeintegrators/ExplicitTVDRK2.h b/framework/include/timeintegrators/ExplicitTVDRK2.h index 98caff9453d0..b8cd92be8e2a 100644 --- a/framework/include/timeintegrators/ExplicitTVDRK2.h +++ b/framework/include/timeintegrators/ExplicitTVDRK2.h @@ -58,6 +58,7 @@ class ExplicitTVDRK2 : public TimeIntegrator ADReal & ad_u_dotdot) const override; virtual void solve() override; virtual void postResidual(NumericVector & residual) override; + virtual bool overridesSolve() const override { return true; } protected: /** diff --git a/framework/include/timeintegrators/ImplicitMidpoint.h b/framework/include/timeintegrators/ImplicitMidpoint.h index 017f6a0daa85..bdba146909c3 100644 --- a/framework/include/timeintegrators/ImplicitMidpoint.h +++ b/framework/include/timeintegrators/ImplicitMidpoint.h @@ -53,6 +53,7 @@ class ImplicitMidpoint : public TimeIntegrator ADReal & ad_u_dotdot) const override; virtual void solve() override; virtual void postResidual(NumericVector & residual) override; + virtual bool overridesSolve() const override { return true; } protected: /** diff --git a/framework/include/timeintegrators/LStableDirk2.h b/framework/include/timeintegrators/LStableDirk2.h index 7b93014f0fa7..cc268a50f0bc 100644 --- a/framework/include/timeintegrators/LStableDirk2.h +++ b/framework/include/timeintegrators/LStableDirk2.h @@ -49,6 +49,7 @@ class LStableDirk2 : public TimeIntegrator ADReal & ad_u_dotdot) const override; virtual void solve() override; virtual void postResidual(NumericVector & residual) override; + virtual bool overridesSolve() const override { return true; } protected: /** diff --git a/framework/include/timeintegrators/LStableDirk3.h b/framework/include/timeintegrators/LStableDirk3.h index f437314019e4..ccf2ef76e512 100644 --- a/framework/include/timeintegrators/LStableDirk3.h +++ b/framework/include/timeintegrators/LStableDirk3.h @@ -51,6 +51,7 @@ class LStableDirk3 : public TimeIntegrator ADReal & ad_u_dotdot) const override; virtual void solve() override; virtual void postResidual(NumericVector & residual) override; + virtual bool overridesSolve() const override { return true; } protected: /** diff --git a/framework/include/timeintegrators/LStableDirk4.h b/framework/include/timeintegrators/LStableDirk4.h index 64aa1cb71e00..119cdf68cbd9 100644 --- a/framework/include/timeintegrators/LStableDirk4.h +++ b/framework/include/timeintegrators/LStableDirk4.h @@ -62,6 +62,7 @@ class LStableDirk4 : public TimeIntegrator ADReal & ad_u_dotdot) const override; virtual void solve() override; virtual void postResidual(NumericVector & residual) override; + virtual bool overridesSolve() const override { return true; } protected: /** diff --git a/framework/include/timeintegrators/TimeIntegrator.h b/framework/include/timeintegrators/TimeIntegrator.h index c7aade25760a..99bb604ac1d3 100644 --- a/framework/include/timeintegrators/TimeIntegrator.h +++ b/framework/include/timeintegrators/TimeIntegrator.h @@ -60,7 +60,7 @@ class TimeIntegrator : public MooseObject, public Restartable * Called _only_ before the very first timestep (t_step = 0) * Never called again (not even during recover/restart) */ - virtual void init() {} + virtual void init(); virtual void preSolve() {} virtual void preStep() {} @@ -153,6 +153,21 @@ class TimeIntegrator : public MooseObject, public Restartable */ TagID uDotDotFactorTag() const { return _u_dotdot_factor_tag; } + /** + * @returns whether this integrator integrates the given variable + */ + bool integratesVar(const unsigned int var_num) const; + + /** + * Record the linear and nonlinear iterations from a just finished solve + */ + void setNumIterationsLastSolve(); + + /** + * @returns Whether the virtual solve method is overridden + */ + virtual bool overridesSolve() const { return false; } + protected: /** * Gets the number of nonlinear iterations in the most recent solve. @@ -164,6 +179,12 @@ class TimeIntegrator : public MooseObject, public Restartable */ unsigned int getNumLinearIterationsLastSolve() const; + /** + * 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 & from, NumericVector & to); + FEProblemBase & _fe_problem; SystemBase & _sys; NonlinearSystemBase & _nl; @@ -200,4 +221,14 @@ class TimeIntegrator : public MooseObject, public Restartable const TagID _u_dot_factor_tag; /// The vector tag for the nodal multiplication factor for the residual calculation of the udotdot term const TagID _u_dotdot_factor_tag; + + /// Whether the user has requested that the time integrator be applied to a subset of variables + bool _var_restriction; + + /// The local degree of freedom indices this time integrator is being applied to. If this + /// container is empty then the time integrator is applied to all indices + std::vector _local_indices; + + /// The variables that this time integrator integrates + std::unordered_set _vars; }; diff --git a/framework/src/actions/SetupTimeIntegratorAction.C b/framework/src/actions/SetupTimeIntegratorAction.C index b8fc67d16762..dc5408b7d106 100644 --- a/framework/src/actions/SetupTimeIntegratorAction.C +++ b/framework/src/actions/SetupTimeIntegratorAction.C @@ -12,6 +12,7 @@ #include "Factory.h" registerMooseAction("MooseApp", SetupTimeIntegratorAction, "setup_time_integrator"); +registerMooseAction("MooseApp", SetupTimeIntegratorAction, "setup_time_integrators"); InputParameters SetupTimeIntegratorAction::validParams() @@ -29,5 +30,13 @@ SetupTimeIntegratorAction::SetupTimeIntegratorAction(const InputParameters & par void SetupTimeIntegratorAction::act() { - _problem->addTimeIntegrator(_type, _name, _moose_object_pars); + std::string name; + // Task: setup_time_integrator corresponding to [TimeIntegrator] block + if (_current_task == "setup_time_integrator") + name = _type; + // Task: setup_time_integrators corresponding to [TimeIntegrators] block + else + name = _name; + + _problem->addTimeIntegrator(_type, name, _moose_object_pars); } diff --git a/framework/src/base/Moose.C b/framework/src/base/Moose.C index 502f7de0798a..9343de8feb89 100644 --- a/framework/src/base/Moose.C +++ b/framework/src/base/Moose.C @@ -170,6 +170,7 @@ addActionTypes(Syntax & syntax) registerMooseObjectTask("add_time_steppers", TimeStepper, false); registerMooseObjectTask("add_time_stepper", TimeStepper, false); registerTask ("compose_time_stepper", true); + registerMooseObjectTask("setup_time_integrators", TimeIntegrator, false); registerMooseObjectTask("setup_time_integrator", TimeIntegrator, false); registerMooseObjectTask("add_preconditioning", MoosePreconditioner, false); @@ -327,7 +328,7 @@ addActionTypes(Syntax & syntax) "(init_component_physics)" // components must add their blocks to physics before init_physics "(init_physics)" "(setup_postprocessor_data)" - "(setup_time_integrator)" + "(setup_time_integrator, setup_time_integrators)" "(setup_executioner)" "(setup_executioner_complete)" "(setup_component)" // no particular reason for that placement @@ -567,7 +568,10 @@ associateSyntaxInner(Syntax & syntax, ActionFactory & /*action_factory*/) registerSyntaxTask("AddTimeStepperAction", "Executioner/TimeStepper", "add_time_stepper"); registerSyntaxTask( "ComposeTimeStepperAction", "Executioner/TimeSteppers", "compose_time_stepper"); - registerSyntax("SetupTimeIntegratorAction", "Executioner/TimeIntegrator"); + registerSyntaxTask( + "SetupTimeIntegratorAction", "Executioner/TimeIntegrators/*", "setup_time_integrators"); + registerSyntaxTask( + "SetupTimeIntegratorAction", "Executioner/TimeIntegrator", "setup_time_integrator"); syntax.registerSyntaxType("Executors/*", "ExecutorName"); registerSyntax("SetupQuadratureAction", "Executioner/Quadrature"); diff --git a/framework/src/executioners/Transient.C b/framework/src/executioners/Transient.C index 7081a8ea93b4..3bdc937ab6c8 100644 --- a/framework/src/executioners/Transient.C +++ b/framework/src/executioners/Transient.C @@ -261,7 +261,9 @@ Transient::preExecute() "input file or report an error.\n" "2. If you are developing a new time stepper, make sure that initial time step " "size in your code is computed correctly."); - _nl.getTimeIntegrator()->init(); + const auto & tis = _nl.getTimeIntegrators(); + for (auto & ti : tis) + ti->init(); } } @@ -456,7 +458,9 @@ Transient::endStep(Real input_time) _time = _time_old; else { - _nl.getTimeIntegrator()->postStep(); + const auto & tis = _nl.getTimeIntegrators(); + for (auto & ti : tis) + ti->postStep(); // Compute the Error Indicators and Markers _problem.computeIndicators(); @@ -703,14 +707,18 @@ Transient::getTimeStepperName() const return std::string(); } -std::string -Transient::getTimeIntegratorName() const +std::vector +Transient::getTimeIntegratorNames() const { - const auto * ti = _nl.getTimeIntegrator(); - if (ti) - return ti->type(); - else + const auto & tis = _nl.getTimeIntegrators(); + if (tis.empty()) mooseError("Time integrator has not been built yet so we can't retrieve its name"); + + std::vector ret; + for (const auto & ti : tis) + ret.push_back(ti->type()); + + return ret; } Real diff --git a/framework/src/outputs/ConsoleUtils.C b/framework/src/outputs/ConsoleUtils.C index a05fe72f52cb..7828d8e67a8f 100644 --- a/framework/src/outputs/ConsoleUtils.C +++ b/framework/src/outputs/ConsoleUtils.C @@ -336,9 +336,10 @@ outputExecutionInformation(const MooseApp & app, FEProblemBase & problem) std::string time_stepper = exec->getTimeStepperName(); if (time_stepper != "") oss << std::setw(console_field_width) << " TimeStepper: " << time_stepper << '\n'; - std::string time_integrator = exec->getTimeIntegratorName(); - if (time_integrator != "") - oss << std::setw(console_field_width) << " TimeIntegrator: " << time_integrator << '\n'; + const auto time_integrator_names = exec->getTimeIntegratorNames(); + if (!time_integrator_names.empty()) + oss << std::setw(console_field_width) + << " TimeIntegrator(s): " << MooseUtils::join(time_integrator_names, ", ") << '\n'; oss << std::setw(console_field_width) << " Solver Mode: " << problem.solverTypeString() << '\n'; diff --git a/framework/src/problems/DisplacedProblem.C b/framework/src/problems/DisplacedProblem.C index 58752471682b..900f690defe9 100644 --- a/framework/src/problems/DisplacedProblem.C +++ b/framework/src/problems/DisplacedProblem.C @@ -194,9 +194,9 @@ void DisplacedProblem::addTimeIntegrator() { for (const auto nl_sys_num : make_range(_mproblem.numNonlinearSystems())) - _displaced_solver_systems[nl_sys_num]->addTimeIntegrator( - _mproblem.getNonlinearSystemBase(nl_sys_num).getSharedTimeIntegrator()); - _displaced_aux->addTimeIntegrator(_mproblem.getAuxiliarySystem().getSharedTimeIntegrator()); + _displaced_solver_systems[nl_sys_num]->copyTimeIntegrators( + _mproblem.getNonlinearSystemBase(nl_sys_num)); + _displaced_aux->copyTimeIntegrators(_mproblem.getAuxiliarySystem()); } void diff --git a/framework/src/problems/FEProblemBase.C b/framework/src/problems/FEProblemBase.C index ae099d36bda0..c80a3d1c4a69 100644 --- a/framework/src/problems/FEProblemBase.C +++ b/framework/src/problems/FEProblemBase.C @@ -1148,12 +1148,12 @@ FEProblemBase::initialSetup() for (auto & sys : _solver_systems) { - auto ti = sys->getTimeIntegrator(); + const auto & tis = sys->getTimeIntegrators(); - if (ti) { TIME_SECTION("timeIntegratorInitialSetup", 5, "Initializing Time Integrator"); - ti->initialSetup(); + for (auto & ti : tis) + ti->initialSetup(); } } @@ -6489,10 +6489,10 @@ FEProblemBase::addTimeIntegrator(const std::string & type, { nl->addDotVectors(); - auto tag_udot = nl->getTimeIntegrator()->uDotFactorTag(); + auto tag_udot = nl->getTimeIntegrators()[0]->uDotFactorTag(); if (!nl->hasVector(tag_udot)) nl->associateVectorToTag(*nl->solutionUDot(), tag_udot); - auto tag_udotdot = nl->getTimeIntegrator()->uDotDotFactorTag(); + auto tag_udotdot = nl->getTimeIntegrators()[0]->uDotDotFactorTag(); if (!nl->hasVector(tag_udotdot) && uDotDotRequested()) nl->associateVectorToTag(*nl->solutionUDotDot(), tag_udotdot); } diff --git a/framework/src/systems/AuxiliarySystem.C b/framework/src/systems/AuxiliarySystem.C index 1c8c8b0a136d..b2724842b707 100644 --- a/framework/src/systems/AuxiliarySystem.C +++ b/framework/src/systems/AuxiliarySystem.C @@ -252,16 +252,6 @@ AuxiliarySystem::addVariable(const std::string & var_type, } } -void -AuxiliarySystem::addTimeIntegrator(const std::string & type, - const std::string & name, - InputParameters & parameters) -{ - parameters.set("_sys") = this; - std::shared_ptr ti = _factory.create(type, name, parameters); - _time_integrator = ti; -} - void AuxiliarySystem::addKernel(const std::string & kernel_name, const std::string & name, @@ -380,8 +370,9 @@ void AuxiliarySystem::compute(ExecFlagType type) { // avoid division by dt which might be zero. - if (_fe_problem.dt() > 0. && _time_integrator) - _time_integrator->preStep(); + if (_fe_problem.dt() > 0.) + for (auto & ti : _time_integrators) + ti->preStep(); // We need to compute time derivatives every time each kind of the variables is finished, because: // @@ -395,8 +386,9 @@ AuxiliarySystem::compute(ExecFlagType type) { computeScalarVars(type); // compute time derivatives of scalar aux variables _after_ the values were updated - if (_fe_problem.dt() > 0. && _time_integrator) - _time_integrator->computeTimeDerivatives(); + if (_fe_problem.dt() > 0.) + for (auto & ti : _time_integrators) + ti->computeTimeDerivatives(); } if (_vars[0].fieldVariables().size() > 0) @@ -410,8 +402,9 @@ AuxiliarySystem::compute(ExecFlagType type) computeElementalVars(type); // compute time derivatives of nodal aux variables _after_ the values were updated - if (_fe_problem.dt() > 0. && _time_integrator) - _time_integrator->computeTimeDerivatives(); + if (_fe_problem.dt() > 0.) + for (auto & ti : _time_integrators) + ti->computeTimeDerivatives(); } if (_serialized_solution.get()) diff --git a/framework/src/systems/DisplacedSystem.C b/framework/src/systems/DisplacedSystem.C index a3a5cb494e54..91e55e642d8e 100644 --- a/framework/src/systems/DisplacedSystem.C +++ b/framework/src/systems/DisplacedSystem.C @@ -52,12 +52,6 @@ DisplacedSystem::getVector(const std::string & name) const return _undisplaced_system.getVector(name); } -void -DisplacedSystem::addTimeIntegrator(std::shared_ptr ti) -{ - _time_integrator = ti; -} - System & DisplacedSystem::system() { diff --git a/framework/src/systems/LinearSystem.C b/framework/src/systems/LinearSystem.C index 56ca1259bd13..cbfb5112dd9b 100644 --- a/framework/src/systems/LinearSystem.C +++ b/framework/src/systems/LinearSystem.C @@ -92,14 +92,6 @@ LinearSystem::LinearSystem(FEProblemBase & fe_problem, const std::string & name) LinearSystem::~LinearSystem() = default; -void -LinearSystem::addTimeIntegrator(const std::string & /*type*/, - const std::string & /*name*/, - InputParameters & /*parameters*/) -{ - mooseError("LinearSystem does not support time integrators yet!"); -} - void LinearSystem::initialSetup() { diff --git a/framework/src/systems/NonlinearEigenSystem.C b/framework/src/systems/NonlinearEigenSystem.C index 78bd527763f9..ec5cd6589a8b 100644 --- a/framework/src/systems/NonlinearEigenSystem.C +++ b/framework/src/systems/NonlinearEigenSystem.C @@ -254,16 +254,25 @@ NonlinearEigenSystem::solve() _eigen_sys.set_initial_space(solution()); } - // Solve the transient problem if we have a time integrator; the - // steady problem if not. - if (_time_integrator) - { - _time_integrator->solve(); - _time_integrator->postSolve(); - } + const bool time_integrator_solve = std::any_of(_time_integrators.begin(), + _time_integrators.end(), + [](auto & ti) { return ti->overridesSolve(); }); + if (time_integrator_solve) + mooseAssert(_time_integrators.size() == 1, + "If solve is overridden, then there must be only one time integrator"); + + if (time_integrator_solve) + _time_integrators.front()->solve(); else system().solve(); + for (auto & ti : _time_integrators) + { + if (!ti->overridesSolve()) + ti->setNumIterationsLastSolve(); + ti->postSolve(); + } + // store eigenvalues unsigned int n_converged_eigenvalues = getNumConvergedEigenvalues(); diff --git a/framework/src/systems/NonlinearSystem.C b/framework/src/systems/NonlinearSystem.C index 5a018c24b87a..23bd3b70dd8a 100644 --- a/framework/src/systems/NonlinearSystem.C +++ b/framework/src/systems/NonlinearSystem.C @@ -163,7 +163,7 @@ NonlinearSystem::solve() _nl_implicit_sys.nonlinear_solver->postcheck = Moose::compute_postcheck; // reset solution invalid counter for the time step - if (_time_integrator) + if (!_time_integrators.empty()) _app.solutionInvalidity().resetSolutionInvalidTimeStep(); if (shouldEvaluatePreSMOResidual()) @@ -184,16 +184,32 @@ NonlinearSystem::solve() potentiallySetupFiniteDifferencing(); - if (_time_integrator) + const bool time_integrator_solve = std::any_of(_time_integrators.begin(), + _time_integrators.end(), + [](auto & ti) { return ti->overridesSolve(); }); + if (time_integrator_solve) + mooseAssert(_time_integrators.size() == 1, + "If solve is overridden, then there must be only one time integrator"); + + if (time_integrator_solve) + _time_integrators.front()->solve(); + else + system().solve(); + + for (auto & ti : _time_integrators) { - _time_integrator->solve(); - _time_integrator->postSolve(); - _n_iters = _time_integrator->getNumNonlinearIterations(); - _n_linear_iters = _time_integrator->getNumLinearIterations(); + if (!ti->overridesSolve()) + ti->setNumIterationsLastSolve(); + ti->postSolve(); + } + + if (!_time_integrators.empty()) + { + _n_iters = _time_integrators.front()->getNumNonlinearIterations(); + _n_linear_iters = _time_integrators.front()->getNumLinearIterations(); } else { - system().solve(); _n_iters = _nl_implicit_sys.n_nonlinear_iterations(); _n_linear_iters = _nl_implicit_sys.nonlinear_solver->get_total_linear_iterations(); } diff --git a/framework/src/systems/NonlinearSystemBase.C b/framework/src/systems/NonlinearSystemBase.C index 1cb35ce8f827..f380e541d37b 100644 --- a/framework/src/systems/NonlinearSystemBase.C +++ b/framework/src/systems/NonlinearSystemBase.C @@ -439,17 +439,6 @@ NonlinearSystemBase::setupFieldDecomposition() top_split->setup(*this); } -void -NonlinearSystemBase::addTimeIntegrator(const std::string & type, - const std::string & name, - InputParameters & parameters) -{ - parameters.set("_sys") = this; - - std::shared_ptr ti = _factory.create(type, name, parameters); - _time_integrator = ti; -} - void NonlinearSystemBase::addKernel(const std::string & kernel_name, const std::string & name, @@ -840,8 +829,11 @@ NonlinearSystemBase::computeResidualTags(const std::set & tags) if (required_residual) { auto & residual = getVector(residualVectorTag()); - if (_time_integrator) - _time_integrator->postResidual(residual); + if (!_time_integrators.empty()) + { + for (auto & ti : _time_integrators) + ti->postResidual(residual); + } else residual += *_Re_non_time; residual.close(); @@ -898,8 +890,11 @@ NonlinearSystemBase::computeResidualAndJacobianTags(const std::set & vect if (required_residual) { auto & residual = getVector(residualVectorTag()); - if (_time_integrator) - _time_integrator->postResidual(residual); + if (!_time_integrators.empty()) + { + for (auto & ti : _time_integrators) + ti->postResidual(residual); + } else residual += *_Re_non_time; residual.close(); @@ -920,8 +915,8 @@ NonlinearSystemBase::computeResidualAndJacobianTags(const std::set & vect void NonlinearSystemBase::onTimestepBegin() { - if (_time_integrator) - _time_integrator->preSolve(); + for (auto & ti : _time_integrators) + ti->preSolve(); if (_predictor.get()) _predictor->timestepSetup(); } diff --git a/framework/src/systems/SolverSystem.C b/framework/src/systems/SolverSystem.C index 376c7147b080..a67435f0e586 100644 --- a/framework/src/systems/SolverSystem.C +++ b/framework/src/systems/SolverSystem.C @@ -148,10 +148,11 @@ SolverSystem::compute(const ExecFlagType type) compute_tds = true; } - if (compute_tds && _fe_problem.dt() > 0. && _time_integrator) - { - // avoid division by dt which might be zero. - _time_integrator->preStep(); - _time_integrator->computeTimeDerivatives(); - } + if (compute_tds && _fe_problem.dt() > 0.) + for (auto & ti : _time_integrators) + { + // avoid division by dt which might be zero. + ti->preStep(); + ti->computeTimeDerivatives(); + } } diff --git a/framework/src/systems/SystemBase.C b/framework/src/systems/SystemBase.C index e51aa4a9bbdd..836f02b66de7 100644 --- a/framework/src/systems/SystemBase.C +++ b/framework/src/systems/SystemBase.C @@ -76,7 +76,6 @@ SystemBase::SystemBase(SubProblem & subproblem, _var_kind(var_kind), _max_var_n_dofs_per_elem(0), _max_var_n_dofs_per_node(0), - _time_integrator(nullptr), _automatic_scaling(false), _verbose(false), _solution_states_initialized(false) @@ -1579,6 +1578,49 @@ SystemBase::serializedSolution() return *_serialized_solution; } +void +SystemBase::addTimeIntegrator(const std::string & type, + const std::string & name, + InputParameters & parameters) +{ + parameters.set("_sys") = this; + _time_integrators.push_back(_factory.create(type, name, parameters)); +} + +void +SystemBase::copyTimeIntegrators(const SystemBase & other_sys) +{ + _time_integrators = other_sys._time_integrators; +} + +const TimeIntegrator * +SystemBase::getPossiblyNullTimeIntegrator(const unsigned int var_num) const +{ + for (auto & ti : _time_integrators) + if (ti->integratesVar(var_num)) + return ti.get(); + + return nullptr; +} + +const TimeIntegrator & +SystemBase::getTimeIntegrator(const unsigned int var_num) const +{ + const auto * const ti = getPossiblyNullTimeIntegrator(var_num); + + if (ti) + return *ti; + else + mooseError("No time integrator found that integrates variable number ", + std::to_string(var_num)); +} + +const std::vector> & +SystemBase::getTimeIntegrators() +{ + return _time_integrators; +} + template MooseVariableFE & SystemBase::getFieldVariable(THREAD_ID tid, const std::string & var_name); diff --git a/framework/src/timeintegrators/CrankNicolson.C b/framework/src/timeintegrators/CrankNicolson.C index 3126bf136bdd..be4f3c7ba631 100644 --- a/framework/src/timeintegrators/CrankNicolson.C +++ b/framework/src/timeintegrators/CrankNicolson.C @@ -52,6 +52,8 @@ CrankNicolson::computeADTimeDerivatives(ADReal & ad_u_dot, void CrankNicolson::init() { + TimeIntegrator::init(); + if (!_sys.solutionUDot()) mooseError("CrankNicolson: Time derivative of solution (`u_dot`) is not stored. Please set " "uDotRequested() to true in FEProblemBase befor requesting `u_dot`."); @@ -70,7 +72,8 @@ CrankNicolson::init() _fe_problem.setCurrentResidualVectorTags({_nl.nonTimeVectorTag()}); _nl.computeResidualTag(_nl.RHS(), _nl.nonTimeVectorTag()); _fe_problem.clearCurrentResidualVectorTags(); - _residual_old = _nl.RHS(); + + copyVector(_nl.RHS(), _residual_old); } void @@ -93,14 +96,31 @@ CrankNicolson::postResidual(NumericVector & residual) _Re_non_time.close(); if (!inputs_closed[2]) _residual_old.close(); - residual += _Re_time; - residual += _Re_non_time; - residual += _residual_old; + + if (!_var_restriction) + { + residual += _Re_time; + residual += _Re_non_time; + residual += _residual_old; + } + else + { + auto residual_sub = residual.get_subvector(_local_indices); + auto re_time_sub = _Re_time.get_subvector(_local_indices); + auto re_non_time_sub = _Re_non_time.get_subvector(_local_indices); + auto residual_old_sub = _residual_old.get_subvector(_local_indices); + *residual_sub += *re_time_sub; + *residual_sub += *re_non_time_sub; + *residual_sub += *residual_old_sub; + residual.restore_subvector(std::move(*residual_sub), _local_indices); + _Re_time.restore_subvector(std::move(*re_time_sub), _local_indices); + _Re_non_time.restore_subvector(std::move(*re_non_time_sub), _local_indices); + _residual_old.restore_subvector(std::move(*residual_old_sub), _local_indices); + } } void CrankNicolson::postStep() { - // shift the residual in time - _residual_old = _Re_non_time; + copyVector(_Re_non_time, _residual_old); } diff --git a/framework/src/timeintegrators/ImplicitEuler.C b/framework/src/timeintegrators/ImplicitEuler.C index 0cc7df9a456b..ba3aba153729 100644 --- a/framework/src/timeintegrators/ImplicitEuler.C +++ b/framework/src/timeintegrators/ImplicitEuler.C @@ -50,7 +50,21 @@ ImplicitEuler::computeADTimeDerivatives(ADReal & ad_u_dot, void ImplicitEuler::postResidual(NumericVector & residual) { - residual += _Re_time; - residual += _Re_non_time; - residual.close(); + if (!_var_restriction) + { + residual += _Re_time; + residual += _Re_non_time; + residual.close(); + } + else + { + auto residual_sub = residual.get_subvector(_local_indices); + auto re_time_sub = _Re_time.get_subvector(_local_indices); + auto re_non_time_sub = _Re_non_time.get_subvector(_local_indices); + *residual_sub += *re_time_sub; + *residual_sub += *re_non_time_sub; + residual.restore_subvector(std::move(*residual_sub), _local_indices); + _Re_time.restore_subvector(std::move(*re_time_sub), _local_indices); + _Re_non_time.restore_subvector(std::move(*re_non_time_sub), _local_indices); + } } diff --git a/framework/src/timeintegrators/TimeIntegrator.C b/framework/src/timeintegrators/TimeIntegrator.C index ae095b38324f..311e94eca937 100644 --- a/framework/src/timeintegrators/TimeIntegrator.C +++ b/framework/src/timeintegrators/TimeIntegrator.C @@ -7,18 +7,20 @@ //* Licensed under LGPL 2.1, please see LICENSE for details //* https://www.gnu.org/licenses/lgpl-2.1.html -#include "libmesh/nonlinear_implicit_system.h" -#include "libmesh/petsc_nonlinear_solver.h" - #include "TimeIntegrator.h" #include "FEProblem.h" -#include "SystemBase.h" -#include "NonlinearSystem.h" +#include "NonlinearSystemBase.h" + +#include "libmesh/nonlinear_implicit_system.h" +#include "libmesh/nonlinear_solver.h" +#include "libmesh/dof_map.h" InputParameters TimeIntegrator::validParams() { InputParameters params = MooseObject::validParams(); + params.addParam>( + "variables", {}, "A subset of the variables that this time integrator should be applied to"); params.registerBase("TimeIntegrator"); return params; } @@ -43,16 +45,61 @@ TimeIntegrator::TimeIntegrator(const InputParameters & parameters) _n_linear_iterations(0), _is_lumped(false), _u_dot_factor_tag(_fe_problem.addVectorTag("u_dot_factor", Moose::VECTOR_TAG_SOLUTION)), - _u_dotdot_factor_tag(_fe_problem.addVectorTag("u_dotdot_factor", Moose::VECTOR_TAG_SOLUTION)) + _u_dotdot_factor_tag(_fe_problem.addVectorTag("u_dotdot_factor", Moose::VECTOR_TAG_SOLUTION)), + _var_restriction(!getParam>("variables").empty()) { _fe_problem.setUDotRequested(true); } +void +TimeIntegrator::init() +{ + if (!_var_restriction) + return; + + const auto & var_names = getParam>("variables"); + std::vector var_num_vec; + auto & lm_sys = _sys.system(); + lm_sys.get_all_variable_numbers(var_num_vec); + std::unordered_set var_nums(var_num_vec.begin(), var_num_vec.end()); + for (const auto & var_name : var_names) + if (lm_sys.has_variable(var_name)) + { + const auto var_num = lm_sys.variable_number(var_name); + _vars.insert(var_num); + var_nums.erase(var_num); + } + + // If var_nums is empty then that means the user has specified all the variables in this system + if (var_nums.empty()) + { + _var_restriction = false; + return; + } + + std::vector var_dof_indices, work_vec; + for (const auto var_num : _vars) + { + work_vec = _local_indices; + _local_indices.clear(); + lm_sys.get_dof_map().local_variable_indices(var_dof_indices, lm_sys.get_mesh(), var_num); + std::merge(work_vec.begin(), + work_vec.end(), + var_dof_indices.begin(), + var_dof_indices.end(), + std::back_inserter(_local_indices)); + } +} + void TimeIntegrator::solve() { - _nl.system().solve(); + mooseError("Calling this method is no longer supported"); +} +void +TimeIntegrator::setNumIterationsLastSolve() +{ _n_nonlinear_iterations = getNumNonlinearIterationsLastSolve(); _n_linear_iterations = getNumLinearIterationsLastSolve(); } @@ -66,8 +113,32 @@ TimeIntegrator::getNumNonlinearIterationsLastSolve() const unsigned int TimeIntegrator::getNumLinearIterationsLastSolve() const { - NonlinearSolver & nonlinear_solver = - static_cast &>(*_nonlinear_implicit_system->nonlinear_solver); + auto & nonlinear_solver = _nonlinear_implicit_system->nonlinear_solver; + libmesh_assert(nonlinear_solver); + + return nonlinear_solver->get_total_linear_iterations(); +} + +void +TimeIntegrator::copyVector(NumericVector & from, NumericVector & to) +{ + if (!_var_restriction) + to = from; + else + { + auto to_sub = to.get_subvector(_local_indices); + auto from_sub = from.get_subvector(_local_indices); + *to_sub = *from_sub; + to.restore_subvector(std::move(*to_sub), _local_indices); + from.restore_subvector(std::move(*from_sub), _local_indices); + } +} + +bool +TimeIntegrator::integratesVar(const unsigned int var_num) const +{ + if (!_var_restriction) + return true; - return nonlinear_solver.get_total_linear_iterations(); + return _vars.count(var_num); } diff --git a/framework/src/timesteppers/AB2PredictorCorrector.C b/framework/src/timesteppers/AB2PredictorCorrector.C index f30e54adc834..5a3bff95145f 100644 --- a/framework/src/timesteppers/AB2PredictorCorrector.C +++ b/framework/src/timesteppers/AB2PredictorCorrector.C @@ -166,8 +166,10 @@ Real AB2PredictorCorrector::estimateTimeError(NumericVector & solution) { _pred1 = _fe_problem.getNonlinearSystemBase(/*nl_sys=*/0).getPredictor()->solutionPredictor(); - TimeIntegrator * ti = _fe_problem.getNonlinearSystemBase(/*nl_sys=*/0).getTimeIntegrator(); - auto scheme = Moose::stringToEnum(ti->type()); + const auto & ti = + _fe_problem.getNonlinearSystemBase(/*nl_sys=*/0).getTimeIntegrator(/*var_num=*/0); + + auto scheme = Moose::stringToEnum(ti.type()); Real dt_old = _my_dt_old; if (dt_old == 0) dt_old = _dt; diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index 181c9b46d503..82fb70aeb951 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -77,7 +77,7 @@ MooseVariableData::MooseVariableData(const MooseVariableField::MooseVariableDataFV(const MooseVariableFV(&_sys) ? true : false), _qrule(nullptr) diff --git a/framework/src/variables/MooseVariableDataLinearFV.C b/framework/src/variables/MooseVariableDataLinearFV.C index 5cc60e6c9efb..94eac24b1361 100644 --- a/framework/src/variables/MooseVariableDataLinearFV.C +++ b/framework/src/variables/MooseVariableDataLinearFV.C @@ -41,7 +41,7 @@ MooseVariableDataLinearFV::MooseVariableDataLinearFV( _var_num(_var.number()), _assembly(_subproblem.assembly(_tid, var.kind() == Moose::VAR_SOLVER ? sys.number() : 0)), _element_type(element_type), - _time_integrator(_sys.getTimeIntegrator()), + _time_integrator(_sys.getPossiblyNullTimeIntegrator(_var_num)), _elem(elem), _displaced(dynamic_cast(&_sys) ? true : false), _qrule(nullptr) diff --git a/framework/src/variables/MooseVariableField.C b/framework/src/variables/MooseVariableField.C index da3685e09254..784c128521a6 100644 --- a/framework/src/variables/MooseVariableField.C +++ b/framework/src/variables/MooseVariableField.C @@ -27,7 +27,7 @@ MooseVariableField::MooseVariableField(const InputParameters & param : MooseVariableFieldBase(parameters), Moose::FunctorBase::type>(name()), MeshChangedInterface(parameters), - _time_integrator(_sys.getTimeIntegrator()) + _time_integrator(_sys.getPossiblyNullTimeIntegrator(_var_num)) { } diff --git a/modules/navier_stokes/src/bcs/INSADDisplaceBoundaryBC.C b/modules/navier_stokes/src/bcs/INSADDisplaceBoundaryBC.C index 9b0211fd30f9..e06b695f8ecb 100644 --- a/modules/navier_stokes/src/bcs/INSADDisplaceBoundaryBC.C +++ b/modules/navier_stokes/src/bcs/INSADDisplaceBoundaryBC.C @@ -36,7 +36,7 @@ INSADDisplaceBoundaryBC::INSADDisplaceBoundaryBC(const InputParameters & paramet _component(getParam("component")), _sub_id(_mesh.getSubdomainID(getParam("associated_subdomain"))) { - if (!dynamic_cast(_sys.getTimeIntegrator())) + if (!dynamic_cast(&_sys.getTimeIntegrator(_var.number()))) mooseError("This boundary condition hard-codes a displacement update with the form of an " "implicit Euler discretization. Consequently please use the default time " "integrator, ImplicitEuler."); From 2bc73a685c2563efc49bda18cf298b6d20deac8f Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Tue, 21 May 2024 15:12:20 -0700 Subject: [PATCH 03/16] Don't multiply time derivatives by 2 for Crank-Nicolson --- framework/include/timeintegrators/CrankNicolson.h | 2 +- framework/include/timeintegrators/TimeIntegrator.h | 7 +++++++ framework/src/systems/NonlinearSystemBase.C | 4 ++++ framework/src/timeintegrators/CrankNicolson.C | 14 +++++++++----- 4 files changed, 21 insertions(+), 6 deletions(-) diff --git a/framework/include/timeintegrators/CrankNicolson.h b/framework/include/timeintegrators/CrankNicolson.h index 8b9abfbe0db1..1960e9589777 100644 --- a/framework/include/timeintegrators/CrankNicolson.h +++ b/framework/include/timeintegrators/CrankNicolson.h @@ -50,5 +50,5 @@ void CrankNicolson::computeTimeDerivativeHelper(T & u_dot, const T2 & u_old) const { u_dot -= u_old; - u_dot *= 2. / _dt; + u_dot *= 1. / _dt; } diff --git a/framework/include/timeintegrators/TimeIntegrator.h b/framework/include/timeintegrators/TimeIntegrator.h index 99bb604ac1d3..43b6f1692436 100644 --- a/framework/include/timeintegrators/TimeIntegrator.h +++ b/framework/include/timeintegrators/TimeIntegrator.h @@ -79,6 +79,13 @@ class TimeIntegrator : public MooseObject, public Restartable */ virtual void postResidual(NumericVector & /*residual*/) {} + /** + * Callback to modify the Jacobian. This is useful in time integration schemes which do summations + * of residuals with non-unity coefficients (e.g. Crank-Nicolson in which we add half of the + * non-time residual) + */ + virtual void postJacobian(SparseMatrix & /*jacobian*/) {} + /** * Callback to the TimeIntegrator called immediately after * TimeIntegrator::solve() (so the name does make sense!). See diff --git a/framework/src/systems/NonlinearSystemBase.C b/framework/src/systems/NonlinearSystemBase.C index f380e541d37b..7b0a5f620f59 100644 --- a/framework/src/systems/NonlinearSystemBase.C +++ b/framework/src/systems/NonlinearSystemBase.C @@ -890,10 +890,14 @@ NonlinearSystemBase::computeResidualAndJacobianTags(const std::set & vect if (required_residual) { auto & residual = getVector(residualVectorTag()); + auto & jac = getMatrix(systemMatrixTag()); if (!_time_integrators.empty()) { for (auto & ti : _time_integrators) + { ti->postResidual(residual); + ti->postJacobian(jac); + } } else residual += *_Re_non_time; diff --git a/framework/src/timeintegrators/CrankNicolson.C b/framework/src/timeintegrators/CrankNicolson.C index be4f3c7ba631..b8d14fc5b983 100644 --- a/framework/src/timeintegrators/CrankNicolson.C +++ b/framework/src/timeintegrators/CrankNicolson.C @@ -38,7 +38,7 @@ CrankNicolson::computeTimeDerivatives() computeTimeDerivativeHelper(u_dot, _solution_old); u_dot.close(); - _du_dot_du = 2. / _dt; + _du_dot_du = 1. / _dt; } void @@ -73,6 +73,10 @@ CrankNicolson::init() _nl.computeResidualTag(_nl.RHS(), _nl.nonTimeVectorTag()); _fe_problem.clearCurrentResidualVectorTags(); + // We multiplied the non-time residual by 0.5 in the postResidual method. We need to restore its + // value here + _nl.RHS() *= 2; + copyVector(_nl.RHS(), _residual_old); } @@ -100,8 +104,8 @@ CrankNicolson::postResidual(NumericVector & residual) if (!_var_restriction) { residual += _Re_time; - residual += _Re_non_time; - residual += _residual_old; + residual.add(0.5, _Re_non_time); + residual.add(0.5, _residual_old); } else { @@ -110,8 +114,8 @@ CrankNicolson::postResidual(NumericVector & residual) auto re_non_time_sub = _Re_non_time.get_subvector(_local_indices); auto residual_old_sub = _residual_old.get_subvector(_local_indices); *residual_sub += *re_time_sub; - *residual_sub += *re_non_time_sub; - *residual_sub += *residual_old_sub; + residual_sub->add(0.5, *re_non_time_sub); + residual_sub->add(0.5, *residual_old_sub); residual.restore_subvector(std::move(*residual_sub), _local_indices); _Re_time.restore_subvector(std::move(*re_time_sub), _local_indices); _Re_non_time.restore_subvector(std::move(*re_non_time_sub), _local_indices); From d80a8c494a123387a6f6b9ef058da4e54316f2ac Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Tue, 21 May 2024 15:46:15 -0700 Subject: [PATCH 04/16] Revert "Don't multiply time derivatives by 2 for Crank-Nicolson" This reverts commit 572f7f6d312dcb171e280099d3d837802f42a007. --- framework/include/timeintegrators/CrankNicolson.h | 2 +- framework/include/timeintegrators/TimeIntegrator.h | 7 ------- framework/src/systems/NonlinearSystemBase.C | 4 ---- framework/src/timeintegrators/CrankNicolson.C | 14 +++++--------- 4 files changed, 6 insertions(+), 21 deletions(-) diff --git a/framework/include/timeintegrators/CrankNicolson.h b/framework/include/timeintegrators/CrankNicolson.h index 1960e9589777..8b9abfbe0db1 100644 --- a/framework/include/timeintegrators/CrankNicolson.h +++ b/framework/include/timeintegrators/CrankNicolson.h @@ -50,5 +50,5 @@ void CrankNicolson::computeTimeDerivativeHelper(T & u_dot, const T2 & u_old) const { u_dot -= u_old; - u_dot *= 1. / _dt; + u_dot *= 2. / _dt; } diff --git a/framework/include/timeintegrators/TimeIntegrator.h b/framework/include/timeintegrators/TimeIntegrator.h index 43b6f1692436..99bb604ac1d3 100644 --- a/framework/include/timeintegrators/TimeIntegrator.h +++ b/framework/include/timeintegrators/TimeIntegrator.h @@ -79,13 +79,6 @@ class TimeIntegrator : public MooseObject, public Restartable */ virtual void postResidual(NumericVector & /*residual*/) {} - /** - * Callback to modify the Jacobian. This is useful in time integration schemes which do summations - * of residuals with non-unity coefficients (e.g. Crank-Nicolson in which we add half of the - * non-time residual) - */ - virtual void postJacobian(SparseMatrix & /*jacobian*/) {} - /** * Callback to the TimeIntegrator called immediately after * TimeIntegrator::solve() (so the name does make sense!). See diff --git a/framework/src/systems/NonlinearSystemBase.C b/framework/src/systems/NonlinearSystemBase.C index 7b0a5f620f59..f380e541d37b 100644 --- a/framework/src/systems/NonlinearSystemBase.C +++ b/framework/src/systems/NonlinearSystemBase.C @@ -890,14 +890,10 @@ NonlinearSystemBase::computeResidualAndJacobianTags(const std::set & vect if (required_residual) { auto & residual = getVector(residualVectorTag()); - auto & jac = getMatrix(systemMatrixTag()); if (!_time_integrators.empty()) { for (auto & ti : _time_integrators) - { ti->postResidual(residual); - ti->postJacobian(jac); - } } else residual += *_Re_non_time; diff --git a/framework/src/timeintegrators/CrankNicolson.C b/framework/src/timeintegrators/CrankNicolson.C index b8d14fc5b983..be4f3c7ba631 100644 --- a/framework/src/timeintegrators/CrankNicolson.C +++ b/framework/src/timeintegrators/CrankNicolson.C @@ -38,7 +38,7 @@ CrankNicolson::computeTimeDerivatives() computeTimeDerivativeHelper(u_dot, _solution_old); u_dot.close(); - _du_dot_du = 1. / _dt; + _du_dot_du = 2. / _dt; } void @@ -73,10 +73,6 @@ CrankNicolson::init() _nl.computeResidualTag(_nl.RHS(), _nl.nonTimeVectorTag()); _fe_problem.clearCurrentResidualVectorTags(); - // We multiplied the non-time residual by 0.5 in the postResidual method. We need to restore its - // value here - _nl.RHS() *= 2; - copyVector(_nl.RHS(), _residual_old); } @@ -104,8 +100,8 @@ CrankNicolson::postResidual(NumericVector & residual) if (!_var_restriction) { residual += _Re_time; - residual.add(0.5, _Re_non_time); - residual.add(0.5, _residual_old); + residual += _Re_non_time; + residual += _residual_old; } else { @@ -114,8 +110,8 @@ CrankNicolson::postResidual(NumericVector & residual) auto re_non_time_sub = _Re_non_time.get_subvector(_local_indices); auto residual_old_sub = _residual_old.get_subvector(_local_indices); *residual_sub += *re_time_sub; - residual_sub->add(0.5, *re_non_time_sub); - residual_sub->add(0.5, *residual_old_sub); + *residual_sub += *re_non_time_sub; + *residual_sub += *residual_old_sub; residual.restore_subvector(std::move(*residual_sub), _local_indices); _Re_time.restore_subvector(std::move(*re_time_sub), _local_indices); _Re_non_time.restore_subvector(std::move(*re_non_time_sub), _local_indices); From 0ac761cfc4b89add49d89b67bda1380c6710acd0 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Wed, 22 May 2024 06:46:37 -0700 Subject: [PATCH 05/16] sub-vector time derivative calculation --- framework/include/systems/DisplacedSystem.h | 10 +++++-- framework/include/systems/SystemBase.h | 8 +++--- .../include/timeintegrators/TimeIntegrator.h | 8 ++++-- framework/src/systems/SystemBase.C | 1 + framework/src/timeintegrators/AStableDirk4.C | 4 ++- .../timeintegrators/ActuallyExplicitEuler.C | 4 ++- framework/src/timeintegrators/BDF2.C | 8 ++++-- .../src/timeintegrators/CentralDifference.C | 4 ++- framework/src/timeintegrators/CrankNicolson.C | 27 +++++++++++++++---- framework/src/timeintegrators/ExplicitEuler.C | 4 ++- framework/src/timeintegrators/ExplicitRK2.C | 4 ++- .../timeintegrators/ExplicitSSPRungeKutta.C | 4 ++- .../src/timeintegrators/ExplicitTVDRK2.C | 4 ++- framework/src/timeintegrators/ImplicitEuler.C | 23 +++++++++++++--- .../src/timeintegrators/ImplicitMidpoint.C | 4 ++- framework/src/timeintegrators/LStableDirk2.C | 4 ++- framework/src/timeintegrators/LStableDirk3.C | 4 ++- framework/src/timeintegrators/LStableDirk4.C | 4 ++- framework/src/timeintegrators/NewmarkBeta.C | 4 ++- .../src/timeintegrators/TimeIntegrator.C | 3 +++ framework/src/variables/MooseVariableData.C | 2 +- .../src/variables/MooseVariableDataBase.C | 4 +-- framework/src/variables/MooseVariableScalar.C | 2 +- 23 files changed, 110 insertions(+), 34 deletions(-) diff --git a/framework/include/systems/DisplacedSystem.h b/framework/include/systems/DisplacedSystem.h index 436e97404a28..00ffcc6a5c0c 100644 --- a/framework/include/systems/DisplacedSystem.h +++ b/framework/include/systems/DisplacedSystem.h @@ -149,9 +149,9 @@ class DisplacedSystem : public SystemBase return _undisplaced_system.solutionUDotDotOld(); } - virtual Number & duDotDu() override { return _undisplaced_system.duDotDu(); } + virtual std::vector & duDotDu() override { return _undisplaced_system.duDotDu(); } virtual Number & duDotDotDu() override { return _undisplaced_system.duDotDotDu(); } - virtual const Number & duDotDu() const override { return _undisplaced_system.duDotDu(); } + virtual const Number & duDotDu(const unsigned int var_num) const override; virtual const Number & duDotDotDu() const override { return _undisplaced_system.duDotDotDu(); } virtual void addDotVectors() override { _undisplaced_system.addDotVectors(); } @@ -294,3 +294,9 @@ DisplacedSystem::hasSolutionState(const unsigned int state, { return _undisplaced_system.hasSolutionState(state, iteration_type); } + +inline const Number & +DisplacedSystem::duDotDu(const unsigned int var_num) const +{ + return _undisplaced_system.duDotDu(var_num); +} diff --git a/framework/include/systems/SystemBase.h b/framework/include/systems/SystemBase.h index 919dc6349184..b960b8a750e9 100644 --- a/framework/include/systems/SystemBase.h +++ b/framework/include/systems/SystemBase.h @@ -244,9 +244,9 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface */ virtual void addDotVectors(); - virtual Number & duDotDu() { return _du_dot_du; } + virtual std::vector & duDotDu() { return _du_dot_du; } virtual Number & duDotDotDu() { return _du_dotdot_du; } - virtual const Number & duDotDu() const { return _du_dot_du; } + virtual const Number & duDotDu(const unsigned int var_num) const { return _du_dot_du[var_num]; } virtual const Number & duDotDotDu() const { return _du_dotdot_du; } virtual NumericVector * solutionUDot() { return _u_dot; } @@ -1005,7 +1005,9 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface /// old solution vector for u^dotdot NumericVector * _u_dotdot_old; - Real _du_dot_du; + /// Derivative of time derivative of u with respect to uj. This depends on the time integration + /// scheme + std::vector _du_dot_du; Real _du_dotdot_du; /// Tagged vectors (pointer) diff --git a/framework/include/timeintegrators/TimeIntegrator.h b/framework/include/timeintegrators/TimeIntegrator.h index 99bb604ac1d3..b6d7bfd3ab9a 100644 --- a/framework/include/timeintegrators/TimeIntegrator.h +++ b/framework/include/timeintegrators/TimeIntegrator.h @@ -198,11 +198,15 @@ class TimeIntegrator : public MooseObject, public Restartable /// residual vector for non-time contributions NumericVector & _Re_non_time; - /// Derivative of time derivative with respect to current solution: \f$ {du^dot}\over{du} \f$ - Real & _du_dot_du; + /// Derivative of time derivative with respect to current solution: \f$ {du^dot}\over{du} \f$ for + /// the different variables. We will only modify the elements in this vector corresponding to the + /// variables that we integrate + std::vector & _du_dot_du; /// solution vectors const NumericVector * const & _solution; const NumericVector & _solution_old; + std::unique_ptr> _solution_sub; + std::unique_ptr> _solution_old_sub; // int & _t_step; // diff --git a/framework/src/systems/SystemBase.C b/framework/src/systems/SystemBase.C index 836f02b66de7..72494d1e46d4 100644 --- a/framework/src/systems/SystemBase.C +++ b/framework/src/systems/SystemBase.C @@ -796,6 +796,7 @@ SystemBase::addVariable(const std::string & var_type, // getMaxVariableNumber is an API method used in Rattlesnake if (var_num > _max_var_number) _max_var_number = var_num; + _du_dot_du.resize(var_num + 1); } bool diff --git a/framework/src/timeintegrators/AStableDirk4.C b/framework/src/timeintegrators/AStableDirk4.C index 5e6cd0af00cd..4a23f5580092 100644 --- a/framework/src/timeintegrators/AStableDirk4.C +++ b/framework/src/timeintegrators/AStableDirk4.C @@ -92,7 +92,9 @@ AStableDirk4::computeTimeDerivatives() u_dot = *_solution; computeTimeDerivativeHelper(u_dot, _solution_old); u_dot.close(); - _du_dot_du = 1. / _dt; + for (const auto i : index_range(_du_dot_du)) + if (integratesVar(i)) + _du_dot_du[i] = 1. / _dt; } void diff --git a/framework/src/timeintegrators/ActuallyExplicitEuler.C b/framework/src/timeintegrators/ActuallyExplicitEuler.C index 0d3f332a16e4..5b4a154da796 100644 --- a/framework/src/timeintegrators/ActuallyExplicitEuler.C +++ b/framework/src/timeintegrators/ActuallyExplicitEuler.C @@ -50,7 +50,9 @@ ActuallyExplicitEuler::computeTimeDerivatives() computeTimeDerivativeHelper(u_dot, _solution_old); u_dot.close(); - _du_dot_du = 1.0 / _dt; + for (const auto i : index_range(_du_dot_du)) + if (integratesVar(i)) + _du_dot_du[i] = 1.0 / _dt; } void diff --git a/framework/src/timeintegrators/BDF2.C b/framework/src/timeintegrators/BDF2.C index 416961d7f432..43a3f6bb15ab 100644 --- a/framework/src/timeintegrators/BDF2.C +++ b/framework/src/timeintegrators/BDF2.C @@ -52,12 +52,16 @@ BDF2::computeTimeDerivatives() if (_t_step == 1) { u_dot = *_solution; - _du_dot_du = 1. / _dt; + for (const auto i : index_range(_du_dot_du)) + if (integratesVar(i)) + _du_dot_du[i] = 1. / _dt; } else { u_dot.zero(); - _du_dot_du = _weight[0] / _dt; + 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(); diff --git a/framework/src/timeintegrators/CentralDifference.C b/framework/src/timeintegrators/CentralDifference.C index a3cfb30eff89..869df4988321 100644 --- a/framework/src/timeintegrators/CentralDifference.C +++ b/framework/src/timeintegrators/CentralDifference.C @@ -87,7 +87,9 @@ CentralDifference::computeTimeDerivatives() u_dot.close(); // used for Jacobian calculations - _du_dot_du = 1.0 / (2 * _dt); + for (const auto i : index_range(_du_dot_du)) + if (integratesVar(i)) + _du_dot_du[i] = 1.0 / (2 * _dt); _du_dotdot_du = 1.0 / (_dt * _dt); // Computing udotdot "factor" diff --git a/framework/src/timeintegrators/CrankNicolson.C b/framework/src/timeintegrators/CrankNicolson.C index be4f3c7ba631..4fcf5d149bbb 100644 --- a/framework/src/timeintegrators/CrankNicolson.C +++ b/framework/src/timeintegrators/CrankNicolson.C @@ -34,11 +34,26 @@ CrankNicolson::computeTimeDerivatives() "uDotRequested() to true in FEProblemBase befor requesting `u_dot`."); NumericVector & u_dot = *_sys.solutionUDot(); - u_dot = *_solution; - computeTimeDerivativeHelper(u_dot, _solution_old); - u_dot.close(); + if (!_var_restriction) + { + u_dot = *_solution; + computeTimeDerivativeHelper(u_dot, _solution_old); + } + else + { + auto u_dot_sub = u_dot.get_subvector(_local_indices); + _solution->create_subvector(*_solution_sub, _local_indices, false); + _solution_old.create_subvector(*_solution_old_sub, _local_indices, false); + *u_dot_sub = *_solution_sub; + computeTimeDerivativeHelper(*u_dot_sub, *_solution_old_sub); + u_dot.restore_subvector(std::move(*u_dot_sub), _local_indices); + // Scatter info needed for ghosts + u_dot.close(); + } - _du_dot_du = 2. / _dt; + for (const auto i : index_range(_du_dot_du)) + if (integratesVar(i)) + _du_dot_du[i] = 2. / _dt; } void @@ -61,7 +76,9 @@ CrankNicolson::init() // time derivative is assumed to be zero on initial NumericVector & u_dot = *_sys.solutionUDot(); u_dot.zero(); - _du_dot_du = 0; + for (const auto i : index_range(_du_dot_du)) + if (integratesVar(i)) + _du_dot_du[i] = 0; // compute residual for the initial time step // Note: we can not directly pass _residual_old in computeResidualTag because diff --git a/framework/src/timeintegrators/ExplicitEuler.C b/framework/src/timeintegrators/ExplicitEuler.C index 1639be7cc1ad..bf56516eae59 100644 --- a/framework/src/timeintegrators/ExplicitEuler.C +++ b/framework/src/timeintegrators/ExplicitEuler.C @@ -44,7 +44,9 @@ ExplicitEuler::computeTimeDerivatives() computeTimeDerivativeHelper(u_dot, _solution_old); u_dot.close(); - _du_dot_du = 1.0 / _dt; + for (const auto i : index_range(_du_dot_du)) + if (integratesVar(i)) + _du_dot_du[i] = 1.0 / _dt; } void diff --git a/framework/src/timeintegrators/ExplicitRK2.C b/framework/src/timeintegrators/ExplicitRK2.C index 8f8d95889463..b99af6045e1c 100644 --- a/framework/src/timeintegrators/ExplicitRK2.C +++ b/framework/src/timeintegrators/ExplicitRK2.C @@ -54,7 +54,9 @@ ExplicitRK2::computeTimeDerivatives() u_dot = *_solution; computeTimeDerivativeHelper(u_dot, _solution_old, _solution_older); - _du_dot_du = 1. / _dt; + for (const auto i : index_range(_du_dot_du)) + if (integratesVar(i)) + _du_dot_du[i] = 1. / _dt; u_dot.close(); } diff --git a/framework/src/timeintegrators/ExplicitSSPRungeKutta.C b/framework/src/timeintegrators/ExplicitSSPRungeKutta.C index 9e1b54196188..f238db646c4b 100644 --- a/framework/src/timeintegrators/ExplicitSSPRungeKutta.C +++ b/framework/src/timeintegrators/ExplicitSSPRungeKutta.C @@ -69,7 +69,9 @@ void ExplicitSSPRungeKutta::computeTimeDerivatives() { // Only the Jacobian needs to be computed, since the mass matrix needs it - _du_dot_du = 1.0 / (_b[_stage] * _dt); + for (const auto i : index_range(_du_dot_du)) + if (integratesVar(i)) + _du_dot_du[i] = 1.0 / (_b[_stage] * _dt); } void diff --git a/framework/src/timeintegrators/ExplicitTVDRK2.C b/framework/src/timeintegrators/ExplicitTVDRK2.C index 1a625c6d904b..aed081626516 100644 --- a/framework/src/timeintegrators/ExplicitTVDRK2.C +++ b/framework/src/timeintegrators/ExplicitTVDRK2.C @@ -56,7 +56,9 @@ ExplicitTVDRK2::computeTimeDerivatives() u_dot = *_solution; computeTimeDerivativeHelper(u_dot, _solution_old, _solution_older); - _du_dot_du = 1. / _dt; + for (const auto i : index_range(_du_dot_du)) + if (integratesVar(i)) + _du_dot_du[i] = 1. / _dt; u_dot.close(); } diff --git a/framework/src/timeintegrators/ImplicitEuler.C b/framework/src/timeintegrators/ImplicitEuler.C index ba3aba153729..7cf2c2af486a 100644 --- a/framework/src/timeintegrators/ImplicitEuler.C +++ b/framework/src/timeintegrators/ImplicitEuler.C @@ -32,11 +32,26 @@ ImplicitEuler::computeTimeDerivatives() "uDotRequested() to true in FEProblemBase befor requesting `u_dot`."); NumericVector & u_dot = *_sys.solutionUDot(); - u_dot = *_solution; - computeTimeDerivativeHelper(u_dot, _solution_old); - u_dot.close(); + if (!_var_restriction) + { + u_dot = *_solution; + computeTimeDerivativeHelper(u_dot, _solution_old); + } + else + { + auto u_dot_sub = u_dot.get_subvector(_local_indices); + _solution->create_subvector(*_solution_sub, _local_indices, false); + _solution_old.create_subvector(*_solution_old_sub, _local_indices, false); + *u_dot_sub = *_solution_sub; + computeTimeDerivativeHelper(*u_dot_sub, *_solution_old_sub); + u_dot.restore_subvector(std::move(*u_dot_sub), _local_indices); + // Scatter info needed for ghosts + u_dot.close(); + } - _du_dot_du = 1.0 / _dt; + for (const auto i : index_range(_du_dot_du)) + if (integratesVar(i)) + _du_dot_du[i] = 1.0 / _dt; } void diff --git a/framework/src/timeintegrators/ImplicitMidpoint.C b/framework/src/timeintegrators/ImplicitMidpoint.C index 7d7daa654ebe..8bb834dec8de 100644 --- a/framework/src/timeintegrators/ImplicitMidpoint.C +++ b/framework/src/timeintegrators/ImplicitMidpoint.C @@ -45,7 +45,9 @@ ImplicitMidpoint::computeTimeDerivatives() u_dot = *_solution; computeTimeDerivativeHelper(u_dot, _solution_old); u_dot.close(); - _du_dot_du = 1. / _dt; + for (const auto i : index_range(_du_dot_du)) + if (integratesVar(i)) + _du_dot_du[i] = 1. / _dt; } void diff --git a/framework/src/timeintegrators/LStableDirk2.C b/framework/src/timeintegrators/LStableDirk2.C index 5904203e5ea2..cfd22dc8d25b 100644 --- a/framework/src/timeintegrators/LStableDirk2.C +++ b/framework/src/timeintegrators/LStableDirk2.C @@ -48,7 +48,9 @@ LStableDirk2::computeTimeDerivatives() u_dot = *_solution; computeTimeDerivativeHelper(u_dot, _solution_old); u_dot.close(); - _du_dot_du = 1. / _dt; + for (const auto i : index_range(_du_dot_du)) + if (integratesVar(i)) + _du_dot_du[i] = 1. / _dt; } void diff --git a/framework/src/timeintegrators/LStableDirk3.C b/framework/src/timeintegrators/LStableDirk3.C index 9fec25a0deca..0412e0ae680b 100644 --- a/framework/src/timeintegrators/LStableDirk3.C +++ b/framework/src/timeintegrators/LStableDirk3.C @@ -67,7 +67,9 @@ LStableDirk3::computeTimeDerivatives() u_dot = *_solution; computeTimeDerivativeHelper(u_dot, _solution_old); u_dot.close(); - _du_dot_du = 1. / _dt; + for (const auto i : index_range(_du_dot_du)) + if (integratesVar(i)) + _du_dot_du[i] = 1. / _dt; } void diff --git a/framework/src/timeintegrators/LStableDirk4.C b/framework/src/timeintegrators/LStableDirk4.C index f144ced1ee8b..ea9af65567da 100644 --- a/framework/src/timeintegrators/LStableDirk4.C +++ b/framework/src/timeintegrators/LStableDirk4.C @@ -62,7 +62,9 @@ LStableDirk4::computeTimeDerivatives() u_dot = *_solution; computeTimeDerivativeHelper(u_dot, _solution_old); u_dot.close(); - _du_dot_du = 1. / _dt; + for (const auto i : index_range(_du_dot_du)) + if (integratesVar(i)) + _du_dot_du[i] = 1. / _dt; } void diff --git a/framework/src/timeintegrators/NewmarkBeta.C b/framework/src/timeintegrators/NewmarkBeta.C index 4fa100cbdef1..0fc03741a4f5 100644 --- a/framework/src/timeintegrators/NewmarkBeta.C +++ b/framework/src/timeintegrators/NewmarkBeta.C @@ -92,7 +92,9 @@ NewmarkBeta::computeTimeDerivatives() // used for Jacobian calculations _du_dotdot_du = 1.0 / _beta / _dt / _dt; - _du_dot_du = _gamma / _beta / _dt; + for (const auto i : index_range(_du_dot_du)) + if (integratesVar(i)) + _du_dot_du[i] = _gamma / _beta / _dt; } void diff --git a/framework/src/timeintegrators/TimeIntegrator.C b/framework/src/timeintegrators/TimeIntegrator.C index 311e94eca937..b8451b816e64 100644 --- a/framework/src/timeintegrators/TimeIntegrator.C +++ b/framework/src/timeintegrators/TimeIntegrator.C @@ -89,6 +89,9 @@ TimeIntegrator::init() var_dof_indices.end(), std::back_inserter(_local_indices)); } + + _solution_sub = NumericVector::build(_solution->comm()); + _solution_old_sub = NumericVector::build(_solution_old.comm()); } void diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index 82fb70aeb951..14adab99773b 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -1022,7 +1022,7 @@ MooseVariableData::computeMonomialValues() Real u_dotdot = 0; Real u_dot_old = 0; Real u_dotdot_old = 0; - const Real & du_dot_du = _sys.duDotDu(); + const Real & du_dot_du = _sys.duDotDu(_var_num); const Real & du_dotdot_du = _sys.duDotDotDu(); if (is_transient) diff --git a/framework/src/variables/MooseVariableDataBase.C b/framework/src/variables/MooseVariableDataBase.C index 34917c824b06..da3c3fc96029 100644 --- a/framework/src/variables/MooseVariableDataBase.C +++ b/framework/src/variables/MooseVariableDataBase.C @@ -600,7 +600,7 @@ MooseVariableDataBase::fetchDoFValues() { _dof_du_dot_du.resize(n); for (decltype(n) i = 0; i < n; ++i) - _dof_du_dot_du[i] = _sys.duDotDu(); + _dof_du_dot_du[i] = _sys.duDotDu(_var.number()); } if (_need_du_dotdot_du || _need_dof_du_dotdot_du) { @@ -708,7 +708,7 @@ MooseVariableDataBase::fetchDoFValues() { _dof_du_dot_du.resize(n); for (decltype(n) i = 0; i < n; ++i) - _dof_du_dot_du[i] = _sys.duDotDu(); + _dof_du_dot_du[i] = _sys.duDotDu(_var.number()); } if (_need_du_dotdot_du || _need_dof_du_dotdot_du) { diff --git a/framework/src/variables/MooseVariableScalar.C b/framework/src/variables/MooseVariableScalar.C index 9779c8816f54..f8a9616a9916 100644 --- a/framework/src/variables/MooseVariableScalar.C +++ b/framework/src/variables/MooseVariableScalar.C @@ -100,7 +100,7 @@ MooseVariableScalar::reinit(bool reinit_for_derivative_reordering /* = false*/) const NumericVector * u_dotdot = sys.solutionUDotDot(); const NumericVector * u_dot_old = sys.solutionUDotOld(); const NumericVector * u_dotdot_old = sys.solutionUDotDotOld(); - const Real & du_dot_du = sys.duDotDu(); + const Real & du_dot_du = sys.duDotDu(number()); const Real & du_dotdot_du = sys.duDotDotDu(); auto safe_access_tagged_vectors = sys.subproblem().safeAccessTaggedVectors(); auto safe_access_tagged_matrices = sys.subproblem().safeAccessTaggedMatrices(); From 64bdbde7e35497a7d5539aaea5255dafd2a08649 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Tue, 24 Sep 2024 14:53:54 -0700 Subject: [PATCH 06/16] Update to new restore_subvector API --- framework/src/timeintegrators/CrankNicolson.C | 10 +++++----- framework/src/timeintegrators/ImplicitEuler.C | 8 ++++---- framework/src/timeintegrators/TimeIntegrator.C | 4 ++-- 3 files changed, 11 insertions(+), 11 deletions(-) diff --git a/framework/src/timeintegrators/CrankNicolson.C b/framework/src/timeintegrators/CrankNicolson.C index 4fcf5d149bbb..f5dbbc9c9b43 100644 --- a/framework/src/timeintegrators/CrankNicolson.C +++ b/framework/src/timeintegrators/CrankNicolson.C @@ -46,7 +46,7 @@ CrankNicolson::computeTimeDerivatives() _solution_old.create_subvector(*_solution_old_sub, _local_indices, false); *u_dot_sub = *_solution_sub; computeTimeDerivativeHelper(*u_dot_sub, *_solution_old_sub); - u_dot.restore_subvector(std::move(*u_dot_sub), _local_indices); + u_dot.restore_subvector(std::move(u_dot_sub), _local_indices); // Scatter info needed for ghosts u_dot.close(); } @@ -129,10 +129,10 @@ CrankNicolson::postResidual(NumericVector & residual) *residual_sub += *re_time_sub; *residual_sub += *re_non_time_sub; *residual_sub += *residual_old_sub; - residual.restore_subvector(std::move(*residual_sub), _local_indices); - _Re_time.restore_subvector(std::move(*re_time_sub), _local_indices); - _Re_non_time.restore_subvector(std::move(*re_non_time_sub), _local_indices); - _residual_old.restore_subvector(std::move(*residual_old_sub), _local_indices); + residual.restore_subvector(std::move(residual_sub), _local_indices); + _Re_time.restore_subvector(std::move(re_time_sub), _local_indices); + _Re_non_time.restore_subvector(std::move(re_non_time_sub), _local_indices); + _residual_old.restore_subvector(std::move(residual_old_sub), _local_indices); } } diff --git a/framework/src/timeintegrators/ImplicitEuler.C b/framework/src/timeintegrators/ImplicitEuler.C index 7cf2c2af486a..899865ccf8dc 100644 --- a/framework/src/timeintegrators/ImplicitEuler.C +++ b/framework/src/timeintegrators/ImplicitEuler.C @@ -44,7 +44,7 @@ ImplicitEuler::computeTimeDerivatives() _solution_old.create_subvector(*_solution_old_sub, _local_indices, false); *u_dot_sub = *_solution_sub; computeTimeDerivativeHelper(*u_dot_sub, *_solution_old_sub); - u_dot.restore_subvector(std::move(*u_dot_sub), _local_indices); + u_dot.restore_subvector(std::move(u_dot_sub), _local_indices); // Scatter info needed for ghosts u_dot.close(); } @@ -78,8 +78,8 @@ ImplicitEuler::postResidual(NumericVector & residual) auto re_non_time_sub = _Re_non_time.get_subvector(_local_indices); *residual_sub += *re_time_sub; *residual_sub += *re_non_time_sub; - residual.restore_subvector(std::move(*residual_sub), _local_indices); - _Re_time.restore_subvector(std::move(*re_time_sub), _local_indices); - _Re_non_time.restore_subvector(std::move(*re_non_time_sub), _local_indices); + residual.restore_subvector(std::move(residual_sub), _local_indices); + _Re_time.restore_subvector(std::move(re_time_sub), _local_indices); + _Re_non_time.restore_subvector(std::move(re_non_time_sub), _local_indices); } } diff --git a/framework/src/timeintegrators/TimeIntegrator.C b/framework/src/timeintegrators/TimeIntegrator.C index b8451b816e64..0725d32f8591 100644 --- a/framework/src/timeintegrators/TimeIntegrator.C +++ b/framework/src/timeintegrators/TimeIntegrator.C @@ -132,8 +132,8 @@ TimeIntegrator::copyVector(NumericVector & from, NumericVector & auto to_sub = to.get_subvector(_local_indices); auto from_sub = from.get_subvector(_local_indices); *to_sub = *from_sub; - to.restore_subvector(std::move(*to_sub), _local_indices); - from.restore_subvector(std::move(*from_sub), _local_indices); + to.restore_subvector(std::move(to_sub), _local_indices); + from.restore_subvector(std::move(from_sub), _local_indices); } } From efe12486f93eec1f8278e59016d7a535a2a19d19 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Wed, 9 Oct 2024 14:08:51 -0700 Subject: [PATCH 07/16] Add extra space per Brandon's suggestion --- unit/src/MooseServerTest.C | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/unit/src/MooseServerTest.C b/unit/src/MooseServerTest.C index 683cc26b4cbd..f127beaa4540 100644 --- a/unit/src/MooseServerTest.C +++ b/unit/src/MooseServerTest.C @@ -1457,9 +1457,9 @@ label: v text: v desc: from /Variables/* pos: [16.15]-[16.15] kind: 18 format: r request_char = 6; expect_count = 3; expect_items = R"INPUT( -label: TimeIntegrator text: TimeIntegrator]\n type = \n $0\n[] desc: application named... pos: [21.3]-[21.6] kind: 22 format: snippet -label: TimeStepper text: TimeStepper]\n type = \n $0\n[] desc: application named... pos: [21.3]-[21.6] kind: 22 format: snippet -label: TimeSteppers text: TimeSteppers]\n $0\n[] desc: application named... pos: [21.3]-[21.6] kind: 22 format: snippet +label: TimeIntegrator text: TimeIntegrator]\n type = \n $0\n[] desc: application named... pos: [21.3]-[21.6] kind: 22 format: snippet +label: TimeStepper text: TimeStepper]\n type = \n $0\n[] desc: application named... pos: [21.3]-[21.6] kind: 22 format: snippet +label: TimeSteppers text: TimeSteppers]\n $0\n[] desc: application named... pos: [21.3]-[21.6] kind: 22 format: snippet )INPUT"; check_completions(request_id, doc_uri, request_line, request_char, expect_count, expect_items); From 8722f8b68f16390a2ad47e40d5d94b87bda7cb22 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Wed, 9 Oct 2024 17:14:40 -0700 Subject: [PATCH 08/16] Adapt modules to new time integrator APIs --- modules/fsi/src/kernels/AcousticInertia.C | 8 ++++---- modules/solid_mechanics/include/kernels/InertialForce.h | 2 +- .../include/nodalkernels/NodalRotationalInertia.h | 2 +- .../include/nodalkernels/NodalTranslationalInertia.h | 2 +- modules/solid_mechanics/src/kernels/InertialForce.C | 2 +- .../src/nodalkernels/NodalRotationalInertia.C | 2 +- .../src/nodalkernels/NodalTranslationalInertia.C | 2 +- modules/thermal_hydraulics/src/base/Simulation.C | 2 +- 8 files changed, 11 insertions(+), 11 deletions(-) diff --git a/modules/fsi/src/kernels/AcousticInertia.C b/modules/fsi/src/kernels/AcousticInertia.C index da3a407fb522..556557220cb7 100644 --- a/modules/fsi/src/kernels/AcousticInertia.C +++ b/modules/fsi/src/kernels/AcousticInertia.C @@ -31,11 +31,11 @@ AcousticInertia::AcousticInertia(const InputParameters & parameters) _u_dot_old(_var.uDotOld()), _du_dot_du(_var.duDotDu()), _du_dotdot_du(_var.duDotDotDu()), - _u_dot_factor(_var.vectorTagValue(_sys.getTimeIntegrator()->uDotFactorTag())), - _u_dotdot_factor(_var.vectorTagValue(_sys.getTimeIntegrator()->uDotDotFactorTag())) + _u_dot_factor(_var.vectorTagValue(_sys.getTimeIntegrator(_var.number()).uDotFactorTag())), + _u_dotdot_factor(_var.vectorTagValue(_sys.getTimeIntegrator(_var.number()).uDotDotFactorTag())) { - addFEVariableCoupleableVectorTag(_sys.getTimeIntegrator()->uDotFactorTag()); - addFEVariableCoupleableVectorTag(_sys.getTimeIntegrator()->uDotDotFactorTag()); + addFEVariableCoupleableVectorTag(_sys.getTimeIntegrator(_var.number()).uDotFactorTag()); + addFEVariableCoupleableVectorTag(_sys.getTimeIntegrator(_var.number()).uDotDotFactorTag()); } Real diff --git a/modules/solid_mechanics/include/kernels/InertialForce.h b/modules/solid_mechanics/include/kernels/InertialForce.h index 36e8fc37dfda..3dadd0d874c0 100644 --- a/modules/solid_mechanics/include/kernels/InertialForce.h +++ b/modules/solid_mechanics/include/kernels/InertialForce.h @@ -60,7 +60,7 @@ class InertialForceTempl : public InertialForceParent const VariableValue * _du_dotdot_du; /// The TimeIntegrator - TimeIntegrator & _time_integrator; + const TimeIntegrator & _time_integrator; using InertialForceParent::_dt; using InertialForceParent::_i; diff --git a/modules/solid_mechanics/include/nodalkernels/NodalRotationalInertia.h b/modules/solid_mechanics/include/nodalkernels/NodalRotationalInertia.h index 0b5b54013d31..3fc44860dbda 100644 --- a/modules/solid_mechanics/include/nodalkernels/NodalRotationalInertia.h +++ b/modules/solid_mechanics/include/nodalkernels/NodalRotationalInertia.h @@ -105,5 +105,5 @@ class NodalRotationalInertia : public TimeNodalKernel const VariableValue * _du_dotdot_du; /// The TimeIntegrator - TimeIntegrator & _time_integrator; + const TimeIntegrator & _time_integrator; }; diff --git a/modules/solid_mechanics/include/nodalkernels/NodalTranslationalInertia.h b/modules/solid_mechanics/include/nodalkernels/NodalTranslationalInertia.h index c2801c232cee..25040f0fdcec 100644 --- a/modules/solid_mechanics/include/nodalkernels/NodalTranslationalInertia.h +++ b/modules/solid_mechanics/include/nodalkernels/NodalTranslationalInertia.h @@ -75,5 +75,5 @@ class NodalTranslationalInertia : public TimeNodalKernel const VariableValue * _du_dotdot_du; /// The TimeIntegrator - TimeIntegrator & _time_integrator; + const TimeIntegrator & _time_integrator; }; diff --git a/modules/solid_mechanics/src/kernels/InertialForce.C b/modules/solid_mechanics/src/kernels/InertialForce.C index 1866cf0c8519..2656f1963c73 100644 --- a/modules/solid_mechanics/src/kernels/InertialForce.C +++ b/modules/solid_mechanics/src/kernels/InertialForce.C @@ -61,7 +61,7 @@ InertialForceTempl::InertialForceTempl(const InputParameters & parameters _eta(this->template getGenericMaterialProperty("eta")), _density_scaling(this->template getMaterialProperty("density_scaling")), _alpha(this->template getParam("alpha")), - _time_integrator(*_sys.getTimeIntegrator()) + _time_integrator(_sys.getTimeIntegrator(this->_var.number())) { if (_has_beta && _has_gamma && _has_velocity && _has_acceleration) { diff --git a/modules/solid_mechanics/src/nodalkernels/NodalRotationalInertia.C b/modules/solid_mechanics/src/nodalkernels/NodalRotationalInertia.C index 5e9ec4308d41..4951b3b9d030 100644 --- a/modules/solid_mechanics/src/nodalkernels/NodalRotationalInertia.C +++ b/modules/solid_mechanics/src/nodalkernels/NodalRotationalInertia.C @@ -85,7 +85,7 @@ NodalRotationalInertia::NodalRotationalInertia(const InputParameters & parameter _rot_dot_residual(_nrot), _rot_vel_old_value(_nrot), _rot_dotdot_residual(_nrot), - _time_integrator(*_sys.getTimeIntegrator()) + _time_integrator(_sys.getTimeIntegrator(_var.number())) { if (_has_beta && _has_gamma && _has_rot_velocities && _has_rot_accelerations) { diff --git a/modules/solid_mechanics/src/nodalkernels/NodalTranslationalInertia.C b/modules/solid_mechanics/src/nodalkernels/NodalTranslationalInertia.C index 6700e9bfa1be..52a472c08bf5 100644 --- a/modules/solid_mechanics/src/nodalkernels/NodalTranslationalInertia.C +++ b/modules/solid_mechanics/src/nodalkernels/NodalTranslationalInertia.C @@ -58,7 +58,7 @@ NodalTranslationalInertia::NodalTranslationalInertia(const InputParameters & par _gamma(_has_gamma ? getParam("gamma") : 0.1), _eta(getParam("eta")), _alpha(getParam("alpha")), - _time_integrator(*_sys.getTimeIntegrator()) + _time_integrator(_sys.getTimeIntegrator(_var.number())) { if (_has_beta && _has_gamma && _has_velocity && _has_acceleration) { diff --git a/modules/thermal_hydraulics/src/base/Simulation.C b/modules/thermal_hydraulics/src/base/Simulation.C index 2e3822e9d880..f3625ce87ad4 100644 --- a/modules/thermal_hydraulics/src/base/Simulation.C +++ b/modules/thermal_hydraulics/src/base/Simulation.C @@ -808,7 +808,7 @@ Simulation::couplingMatrixIntegrityCheck() const return; const TimeIntegrator * ti = - _fe_problem.getNonlinearSystemBase(/*nl_sys_num=*/0).getTimeIntegrator(); + _fe_problem.getNonlinearSystemBase(/*nl_sys_num=*/0).getTimeIntegrators()[0].get(); // Yes, this is horrible. Don't ask why... if ((dynamic_cast(ti) != nullptr) || (dynamic_cast(ti) != nullptr) || From fcb8b22b04b3653d7131c7934f5a8dd2fbea96f1 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Thu, 24 Oct 2024 18:25:44 -0700 Subject: [PATCH 09/16] Make subvectors restartable so recover works --- framework/include/outputs/JsonIO.h | 3 ++ framework/include/restart/DataIO.h | 33 +++++++++++++++++++ .../include/timeintegrators/TimeIntegrator.h | 10 +++--- framework/src/outputs/JsonIO.C | 11 +++++++ framework/src/restart/DataIO.C | 12 ++++++- .../src/timeintegrators/TimeIntegrator.C | 9 ++++- 6 files changed, 71 insertions(+), 7 deletions(-) diff --git a/framework/include/outputs/JsonIO.h b/framework/include/outputs/JsonIO.h index 930005e2d3c5..08802c06f90a 100644 --- a/framework/include/outputs/JsonIO.h +++ b/framework/include/outputs/JsonIO.h @@ -26,6 +26,8 @@ template class DenseVector; template class DenseMatrix; +template +class NumericVector; } // Overloads for to_json, which _must_ be overloaded in the namespace @@ -38,6 +40,7 @@ namespace libMesh void to_json(nlohmann::json & json, const Point & p); void to_json(nlohmann::json & json, const DenseVector & vector); void to_json(nlohmann::json & json, const DenseMatrix & matrix); +void to_json(nlohmann::json & json, const std::unique_ptr> & vector); } namespace nlohmann diff --git a/framework/include/restart/DataIO.h b/framework/include/restart/DataIO.h index d3ae60059669..b385b72b94c2 100644 --- a/framework/include/restart/DataIO.h +++ b/framework/include/restart/DataIO.h @@ -35,6 +35,7 @@ #include #include #include +#include #include #include @@ -341,6 +342,18 @@ dataStore(std::ostream & stream, std::unordered_map & m, void * context) } } +template +inline void +dataStore(std::ostream & stream, std::unordered_set & s, void * context) +{ + // First store the size of the set + std::size_t size = s.size(); + dataStore(stream, size, nullptr); + + for (auto & element : s) + dataStore(stream, element, context); +} + template inline void dataStore(std::ostream & stream, std::optional & m, void * context) @@ -406,6 +419,7 @@ template <> void dataStore(std::ostream & stream, RealEigenMatrix & v, void * context); template <> void dataStore(std::ostream & stream, libMesh::Parameters & p, void * context); + template <> /** * Stores an owned numeric vector. @@ -666,6 +680,25 @@ dataLoad(std::istream & stream, std::unordered_map & m, void * context) } } +template +inline void +dataLoad(std::istream & stream, std::unordered_set & s, void * context) +{ + s.clear(); + + // First read the size of the set + std::size_t size = 0; + dataLoad(stream, size, nullptr); + s.reserve(size); + + for (std::size_t i = 0; i < size; i++) + { + T element; + dataLoad(stream, element, context); + s.insert(element); + } +} + template inline void dataLoad(std::istream & stream, std::optional & m, void * context) diff --git a/framework/include/timeintegrators/TimeIntegrator.h b/framework/include/timeintegrators/TimeIntegrator.h index b6d7bfd3ab9a..4e5f82c92c59 100644 --- a/framework/include/timeintegrators/TimeIntegrator.h +++ b/framework/include/timeintegrators/TimeIntegrator.h @@ -205,8 +205,8 @@ class TimeIntegrator : public MooseObject, public Restartable /// solution vectors const NumericVector * const & _solution; const NumericVector & _solution_old; - std::unique_ptr> _solution_sub; - std::unique_ptr> _solution_old_sub; + std::unique_ptr> & _solution_sub; + std::unique_ptr> & _solution_old_sub; // int & _t_step; // @@ -227,12 +227,12 @@ class TimeIntegrator : public MooseObject, public Restartable const TagID _u_dotdot_factor_tag; /// Whether the user has requested that the time integrator be applied to a subset of variables - bool _var_restriction; + bool & _var_restriction; /// The local degree of freedom indices this time integrator is being applied to. If this /// container is empty then the time integrator is applied to all indices - std::vector _local_indices; + std::vector & _local_indices; /// The variables that this time integrator integrates - std::unordered_set _vars; + std::unordered_set & _vars; }; diff --git a/framework/src/outputs/JsonIO.C b/framework/src/outputs/JsonIO.C index e80c6e0fde1e..7d9870d6dc6f 100644 --- a/framework/src/outputs/JsonIO.C +++ b/framework/src/outputs/JsonIO.C @@ -69,4 +69,15 @@ to_json(nlohmann::json & json, const DenseMatrix & matrix) values[i][j] = matrix(i, j); nlohmann::to_json(json, values); } + +void +to_json(nlohmann::json & json, const std::unique_ptr> & vector) +{ + mooseAssert(vector, "vector must be non-null"); + std::vector local_values; + local_values.reserve(vector->local_size()); + for (const auto i : make_range(vector->first_local_index(), vector->last_local_index())) + local_values.push_back((*vector)(i)); + nlohmann::to_json(json, local_values); +} } diff --git a/framework/src/restart/DataIO.C b/framework/src/restart/DataIO.C index a041da188f5a..b47646c7d788 100644 --- a/framework/src/restart/DataIO.C +++ b/framework/src/restart/DataIO.C @@ -326,7 +326,11 @@ dataStore(std::ostream & stream, std::unique_ptr> & v, void * context) { - mooseAssert(v, "Null vector"); + bool have_vector = v.get(); + dataStore(stream, have_vector, context); + if (!have_vector) + return; + mooseAssert(context, "Needs a context of the communicator"); const auto & comm = *static_cast(context); mooseAssert(&comm == &v->comm(), "Inconsistent communicator"); @@ -669,6 +673,12 @@ template <> void dataLoad(std::istream & stream, std::unique_ptr> & v, void * context) { + bool have_vector; + dataLoad(stream, have_vector, context); + + if (!have_vector) + return; + mooseAssert(context, "Needs a context of the communicator"); const auto & comm = *static_cast(context); if (v) diff --git a/framework/src/timeintegrators/TimeIntegrator.C b/framework/src/timeintegrators/TimeIntegrator.C index 0725d32f8591..4db01fe6e335 100644 --- a/framework/src/timeintegrators/TimeIntegrator.C +++ b/framework/src/timeintegrators/TimeIntegrator.C @@ -38,6 +38,10 @@ TimeIntegrator::TimeIntegrator(const InputParameters & parameters) _du_dot_du(_sys.duDotDu()), _solution(_sys.currentSolution()), _solution_old(_sys.solutionState(1)), + _solution_sub(declareRestartableDataWithContext>>( + "solution_sub", &const_cast(this->comm()))), + _solution_old_sub(declareRestartableDataWithContext>>( + "solution_old_sub", &const_cast(this->comm()))), _t_step(_fe_problem.timeStep()), _dt(_fe_problem.dt()), _dt_old(_fe_problem.dtOld()), @@ -46,7 +50,10 @@ TimeIntegrator::TimeIntegrator(const InputParameters & parameters) _is_lumped(false), _u_dot_factor_tag(_fe_problem.addVectorTag("u_dot_factor", Moose::VECTOR_TAG_SOLUTION)), _u_dotdot_factor_tag(_fe_problem.addVectorTag("u_dotdot_factor", Moose::VECTOR_TAG_SOLUTION)), - _var_restriction(!getParam>("variables").empty()) + _var_restriction(declareRestartableData( + "var_restriction", !getParam>("variables").empty())), + _local_indices(declareRestartableData>("local_indices")), + _vars(declareRestartableData>("vars")) { _fe_problem.setUDotRequested(true); } From 0d42adaa7b362d6ae1e0fff6ae738d7e0c141983 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Thu, 24 Oct 2024 21:44:25 -0700 Subject: [PATCH 10/16] Denote that DirectCentralDifference overridesSolve --- .../include/timeintegrators/DirectCentralDifference.h | 1 + 1 file changed, 1 insertion(+) diff --git a/modules/solid_mechanics/include/timeintegrators/DirectCentralDifference.h b/modules/solid_mechanics/include/timeintegrators/DirectCentralDifference.h index 7dd344db4a18..2d61d243fbb7 100644 --- a/modules/solid_mechanics/include/timeintegrators/DirectCentralDifference.h +++ b/modules/solid_mechanics/include/timeintegrators/DirectCentralDifference.h @@ -34,6 +34,7 @@ class DirectCentralDifference : public ExplicitTimeIntegrator virtual void solve() override; virtual void postResidual(NumericVector & residual) override; + virtual bool overridesSolve() const override { return true; } virtual bool performExplicitSolve(SparseMatrix & mass_matrix) override; From a7438aca540d1a430931c23d22b1d8f4fd8b7878 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Thu, 24 Oct 2024 21:48:32 -0700 Subject: [PATCH 11/16] Make overridesSolve pure virtual This forces downstream developers to remember to implement it --- framework/include/timeintegrators/BDF2.h | 1 + framework/include/timeintegrators/CrankNicolson.h | 1 + framework/include/timeintegrators/ExplicitEuler.h | 1 + framework/include/timeintegrators/ImplicitEuler.h | 1 + framework/include/timeintegrators/NewmarkBeta.h | 1 + framework/include/timeintegrators/TimeIntegrator.h | 2 +- 6 files changed, 6 insertions(+), 1 deletion(-) diff --git a/framework/include/timeintegrators/BDF2.h b/framework/include/timeintegrators/BDF2.h index 50ffc9928a6e..e5956adbb89f 100644 --- a/framework/include/timeintegrators/BDF2.h +++ b/framework/include/timeintegrators/BDF2.h @@ -29,6 +29,7 @@ class BDF2 : public TimeIntegrator const dof_id_type & dof, ADReal & ad_u_dotdot) const override; virtual void postResidual(NumericVector & residual) override; + virtual bool overridesSolve() const override { return false; } protected: /** diff --git a/framework/include/timeintegrators/CrankNicolson.h b/framework/include/timeintegrators/CrankNicolson.h index 8b9abfbe0db1..a79b17f7e6ea 100644 --- a/framework/include/timeintegrators/CrankNicolson.h +++ b/framework/include/timeintegrators/CrankNicolson.h @@ -34,6 +34,7 @@ class CrankNicolson : public TimeIntegrator ADReal & ad_u_dotdot) const override; virtual void postResidual(NumericVector & residual) override; virtual void postStep() override; + virtual bool overridesSolve() const override { return false; } protected: /** diff --git a/framework/include/timeintegrators/ExplicitEuler.h b/framework/include/timeintegrators/ExplicitEuler.h index 2a3dabc41dd2..d37ceef7029e 100644 --- a/framework/include/timeintegrators/ExplicitEuler.h +++ b/framework/include/timeintegrators/ExplicitEuler.h @@ -28,6 +28,7 @@ class ExplicitEuler : public TimeIntegrator const dof_id_type & dof, ADReal & ad_u_dotdot) const override; virtual void postResidual(NumericVector & residual) override; + virtual bool overridesSolve() const override { return false; } protected: /** diff --git a/framework/include/timeintegrators/ImplicitEuler.h b/framework/include/timeintegrators/ImplicitEuler.h index 6c8d7ee9bd7f..b71317e5db1b 100644 --- a/framework/include/timeintegrators/ImplicitEuler.h +++ b/framework/include/timeintegrators/ImplicitEuler.h @@ -28,6 +28,7 @@ class ImplicitEuler : public TimeIntegrator const dof_id_type & dof, ADReal & ad_u_dotdot) const override; virtual void postResidual(NumericVector & residual) override; + virtual bool overridesSolve() const override { return false; } protected: /** diff --git a/framework/include/timeintegrators/NewmarkBeta.h b/framework/include/timeintegrators/NewmarkBeta.h index 03ac4af36d54..e2bec22e5d55 100644 --- a/framework/include/timeintegrators/NewmarkBeta.h +++ b/framework/include/timeintegrators/NewmarkBeta.h @@ -29,6 +29,7 @@ class NewmarkBeta : public TimeIntegrator const dof_id_type & dof, ADReal & ad_u_dotdot) const override; virtual void postResidual(NumericVector & residual) override; + virtual bool overridesSolve() const override { return false; } protected: /** diff --git a/framework/include/timeintegrators/TimeIntegrator.h b/framework/include/timeintegrators/TimeIntegrator.h index 4e5f82c92c59..034c8a9bad78 100644 --- a/framework/include/timeintegrators/TimeIntegrator.h +++ b/framework/include/timeintegrators/TimeIntegrator.h @@ -166,7 +166,7 @@ class TimeIntegrator : public MooseObject, public Restartable /** * @returns Whether the virtual solve method is overridden */ - virtual bool overridesSolve() const { return false; } + virtual bool overridesSolve() const = 0; protected: /** From 136a36529640abc0c98382f0c77d30435d063c58 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Thu, 24 Oct 2024 21:57:02 -0700 Subject: [PATCH 12/16] Fix time integrator retrieval in Simulation.C --- modules/thermal_hydraulics/src/base/Simulation.C | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/modules/thermal_hydraulics/src/base/Simulation.C b/modules/thermal_hydraulics/src/base/Simulation.C index f3625ce87ad4..0fc2c001c706 100644 --- a/modules/thermal_hydraulics/src/base/Simulation.C +++ b/modules/thermal_hydraulics/src/base/Simulation.C @@ -807,8 +807,11 @@ Simulation::couplingMatrixIntegrityCheck() const if (!_fe_problem.shouldSolve()) return; - const TimeIntegrator * ti = - _fe_problem.getNonlinearSystemBase(/*nl_sys_num=*/0).getTimeIntegrators()[0].get(); + const TimeIntegrator * ti = nullptr; + const auto & time_integrators = + _fe_problem.getNonlinearSystemBase(/*nl_sys_num=*/0).getTimeIntegrators(); + if (!time_integrators.empty()) + ti = time_integrators.front().get(); // Yes, this is horrible. Don't ask why... if ((dynamic_cast(ti) != nullptr) || (dynamic_cast(ti) != nullptr) || From 236d1fd2f8ca47d576ed22863056a347e8d345c7 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Fri, 25 Oct 2024 15:51:32 -0700 Subject: [PATCH 13/16] Address unity build failure --- framework/src/systems/SystemBase.C | 1 + 1 file changed, 1 insertion(+) diff --git a/framework/src/systems/SystemBase.C b/framework/src/systems/SystemBase.C index 72494d1e46d4..6ddb2c29a0a6 100644 --- a/framework/src/systems/SystemBase.C +++ b/framework/src/systems/SystemBase.C @@ -26,6 +26,7 @@ #include "MooseUtils.h" #include "FVBoundaryCondition.h" #include "FEProblemBase.h" +#include "TimeIntegrator.h" #include "libmesh/dof_map.h" #include "libmesh/string_to_enum.h" From 8bb4b776fc6fb43297398eeff626d4284e0af9c6 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Sat, 26 Oct 2024 11:23:44 -0700 Subject: [PATCH 14/16] Add TimeIntegrators index page --- .../Executioner/TimeIntegrators/index.md | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 framework/doc/content/syntax/Executioner/TimeIntegrators/index.md diff --git a/framework/doc/content/syntax/Executioner/TimeIntegrators/index.md b/framework/doc/content/syntax/Executioner/TimeIntegrators/index.md new file mode 100644 index 000000000000..73e75b6bf13e --- /dev/null +++ b/framework/doc/content/syntax/Executioner/TimeIntegrators/index.md @@ -0,0 +1,18 @@ +# Multiple Time Integrators + +The time integrator system is described in [TimeIntegrator/index.md]. Multiple +time integrators in a single input are supported using +`Executioner/TimeIntegrators` syntax as illustrated below: + +!listing test/tests/time_integrators/multiple-integrators/test.i block=Executioner + +where each individual time integrator has specific `variables` assigned to +it. + +!syntax list /Executioner/TimeIntegrators objects=True actions=False subsystems=False + +!syntax list /Executioner/TimeIntegrators objects=False actions=False subsystems=True + +!syntax list /Executioner/TimeIntegrators objects=False actions=True subsystems=False + +!syntax parameters /Executioner/TimeIntegrators From af5b6d730297233f9a8a79ad4938c8e724871d14 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Sun, 27 Oct 2024 12:35:04 -0700 Subject: [PATCH 15/16] Make changes backwards compatible with Griffin --- framework/include/systems/DisplacedSystem.h | 4 ++-- framework/include/systems/SystemBase.h | 4 ++-- framework/include/timeintegrators/TimeIntegrator.h | 2 +- framework/src/systems/SystemBase.C | 6 ++++++ framework/src/timeintegrators/TimeIntegrator.C | 2 +- 5 files changed, 12 insertions(+), 6 deletions(-) diff --git a/framework/include/systems/DisplacedSystem.h b/framework/include/systems/DisplacedSystem.h index 00ffcc6a5c0c..9a70b8c63dd6 100644 --- a/framework/include/systems/DisplacedSystem.h +++ b/framework/include/systems/DisplacedSystem.h @@ -149,9 +149,9 @@ class DisplacedSystem : public SystemBase return _undisplaced_system.solutionUDotDotOld(); } - virtual std::vector & duDotDu() override { return _undisplaced_system.duDotDu(); } + virtual std::vector & duDotDus() override { return _undisplaced_system.duDotDus(); } virtual Number & duDotDotDu() override { return _undisplaced_system.duDotDotDu(); } - virtual const Number & duDotDu(const unsigned int var_num) const override; + virtual const Number & duDotDu(unsigned int var_num = 0) const override; virtual const Number & duDotDotDu() const override { return _undisplaced_system.duDotDotDu(); } virtual void addDotVectors() override { _undisplaced_system.addDotVectors(); } diff --git a/framework/include/systems/SystemBase.h b/framework/include/systems/SystemBase.h index b960b8a750e9..de914eb7c288 100644 --- a/framework/include/systems/SystemBase.h +++ b/framework/include/systems/SystemBase.h @@ -244,9 +244,9 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface */ virtual void addDotVectors(); - virtual std::vector & duDotDu() { return _du_dot_du; } + virtual std::vector & duDotDus() { return _du_dot_du; } virtual Number & duDotDotDu() { return _du_dotdot_du; } - virtual const Number & duDotDu(const unsigned int var_num) const { return _du_dot_du[var_num]; } + virtual const Number & duDotDu(unsigned int var_num = 0) const; virtual const Number & duDotDotDu() const { return _du_dotdot_du; } virtual NumericVector * solutionUDot() { return _u_dot; } diff --git a/framework/include/timeintegrators/TimeIntegrator.h b/framework/include/timeintegrators/TimeIntegrator.h index 034c8a9bad78..af14a4d2e60d 100644 --- a/framework/include/timeintegrators/TimeIntegrator.h +++ b/framework/include/timeintegrators/TimeIntegrator.h @@ -101,7 +101,7 @@ class TimeIntegrator : public MooseObject, public Restartable * This function is responsible for the following: * - computing the time derivative; a reference is retrieved from _sys.solutionUDot(). * - computing the time derivative Jacobian, stored in _du_dot_du, which is a - * reference to _sys.duDotDu(). + * reference to _sys.duDotDus(). */ virtual void computeTimeDerivatives() = 0; diff --git a/framework/src/systems/SystemBase.C b/framework/src/systems/SystemBase.C index 6ddb2c29a0a6..845ee6525aed 100644 --- a/framework/src/systems/SystemBase.C +++ b/framework/src/systems/SystemBase.C @@ -1623,6 +1623,12 @@ SystemBase::getTimeIntegrators() return _time_integrators; } +const Number & +SystemBase::duDotDu(const unsigned int var_num) const +{ + return _du_dot_du[var_num]; +} + template MooseVariableFE & SystemBase::getFieldVariable(THREAD_ID tid, const std::string & var_name); diff --git a/framework/src/timeintegrators/TimeIntegrator.C b/framework/src/timeintegrators/TimeIntegrator.C index 4db01fe6e335..3ea40728227e 100644 --- a/framework/src/timeintegrators/TimeIntegrator.C +++ b/framework/src/timeintegrators/TimeIntegrator.C @@ -35,7 +35,7 @@ TimeIntegrator::TimeIntegrator(const InputParameters & parameters) _nonlinear_implicit_system(dynamic_cast(&_sys.system())), _Re_time(_nl.getResidualTimeVector()), _Re_non_time(_nl.getResidualNonTimeVector()), - _du_dot_du(_sys.duDotDu()), + _du_dot_du(_sys.duDotDus()), _solution(_sys.currentSolution()), _solution_old(_sys.solutionState(1)), _solution_sub(declareRestartableDataWithContext>>( From 69c996ef13d6fdea3d2e64ec88130810c2610887 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Mon, 28 Oct 2024 11:23:52 -0700 Subject: [PATCH 16/16] 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 --- framework/include/systems/SystemBase.h | 4 +--- framework/include/timeintegrators/BDF2.h | 2 ++ .../timeintegrators/CentralDifference.h | 2 ++ .../include/timeintegrators/CrankNicolson.h | 2 ++ .../timeintegrators/ExplicitSSPRungeKutta.h | 2 ++ .../include/timeintegrators/NewmarkBeta.h | 2 ++ .../include/timeintegrators/TimeIntegrator.h | 16 +++++++++++++- framework/src/restart/DataIO.C | 3 +++ framework/src/systems/SystemBase.C | 4 ++-- framework/src/timeintegrators/AStableDirk4.C | 4 +--- .../timeintegrators/ActuallyExplicitEuler.C | 5 +---- framework/src/timeintegrators/BDF2.C | 20 +++++++++--------- .../src/timeintegrators/CentralDifference.C | 10 ++++++--- framework/src/timeintegrators/CrankNicolson.C | 10 ++++++--- framework/src/timeintegrators/ExplicitEuler.C | 5 +---- framework/src/timeintegrators/ExplicitRK2.C | 6 +----- .../timeintegrators/ExplicitSSPRungeKutta.C | 10 ++++++--- .../src/timeintegrators/ExplicitTVDRK2.C | 4 +--- framework/src/timeintegrators/ImplicitEuler.C | 4 +--- .../src/timeintegrators/ImplicitMidpoint.C | 4 +--- framework/src/timeintegrators/LStableDirk2.C | 4 +--- framework/src/timeintegrators/LStableDirk3.C | 4 +--- framework/src/timeintegrators/LStableDirk4.C | 4 +--- framework/src/timeintegrators/NewmarkBeta.C | 10 ++++++--- .../src/timeintegrators/TimeIntegrator.C | 21 +++++++++++++------ framework/src/variables/MooseVariableData.C | 2 +- framework/src/variables/MooseVariableDataFV.C | 2 +- .../src/variables/MooseVariableDataLinearFV.C | 2 +- framework/src/variables/MooseVariableField.C | 2 +- 29 files changed, 98 insertions(+), 72 deletions(-) diff --git a/framework/include/systems/SystemBase.h b/framework/include/systems/SystemBase.h index de914eb7c288..b3be2f184dd6 100644 --- 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 /*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 index e5956adbb89f..5893756dbaae 100644 --- 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 & _weight; /// The older solution diff --git a/framework/include/timeintegrators/CentralDifference.h b/framework/include/timeintegrators/CentralDifference.h index 175c600c804c..f5f5d40e4af0 100644 --- 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 index a79b17f7e6ea..90e403707191 100644 --- a/framework/include/timeintegrators/CrankNicolson.h +++ b/framework/include/timeintegrators/CrankNicolson.h @@ -43,6 +43,8 @@ class CrankNicolson : public TimeIntegrator template void computeTimeDerivativeHelper(T & u_dot, const T2 & u_old) const; + virtual Real duDotDuCoeff() const override; + NumericVector & _residual_old; }; diff --git a/framework/include/timeintegrators/ExplicitSSPRungeKutta.h b/framework/include/timeintegrators/ExplicitSSPRungeKutta.h index 9423eb76acb3..5176f9130d90 100644 --- 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 index e2bec22e5d55..aabfc15376b9 100644 --- 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 index af14a4d2e60d..21c2896d322e 100644 --- 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 & from, NumericVector & to); + void copyVector(const NumericVector & from, NumericVector & 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 & _vars; + + /// A vector that is used for creating 'from' subvectors in the \p copyVector() method + std::unique_ptr> _from_subvector; }; diff --git a/framework/src/restart/DataIO.C b/framework/src/restart/DataIO.C index b47646c7d788..c5dd8abdbdc4 100644 --- a/framework/src/restart/DataIO.C +++ b/framework/src/restart/DataIO.C @@ -326,6 +326,9 @@ dataStore(std::ostream & stream, std::unique_ptr> & 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 index 845ee6525aed..b42544c57568 100644 --- 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 index 4a23f5580092..08bbbe5e9ac3 100644 --- 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 index 5b4a154da796..a02ef4f262f0 100644 --- 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 index 43a3f6bb15ab..87f5123df0a1 100644 --- a/framework/src/timeintegrators/BDF2.C +++ b/framework/src/timeintegrators/BDF2.C @@ -50,21 +50,12 @@ BDF2::computeTimeDerivatives() NumericVector & 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 & 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 index 869df4988321..1e77b4957930 100644 --- 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 index f5dbbc9c9b43..78cc9a5efa81 100644 --- 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 & 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 index bf56516eae59..26f0cd76d98e 100644 --- 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 index b99af6045e1c..d1ae060345a8 100644 --- a/framework/src/timeintegrators/ExplicitRK2.C +++ b/framework/src/timeintegrators/ExplicitRK2.C @@ -53,11 +53,7 @@ ExplicitRK2::computeTimeDerivatives() NumericVector & 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 index f238db646c4b..1dc48694012b 100644 --- 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 & 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 index aed081626516..01d5a747b924 100644 --- 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 index 899865ccf8dc..5a2a6a447781 100644 --- 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 index 8bb834dec8de..f72077c24538 100644 --- 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 index cfd22dc8d25b..06cc76457950 100644 --- 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 index 0412e0ae680b..006c5e1748d4 100644 --- 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 index ea9af65567da..df5a8f88f788 100644 --- 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 index 0fc03741a4f5..4ff7475fe254 100644 --- 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 & residual) residual += _Re_non_time; residual.close(); } + +Real +NewmarkBeta::duDotDuCoeff() const +{ + return _gamma / _beta; +} diff --git a/framework/src/timeintegrators/TimeIntegrator.C b/framework/src/timeintegrators/TimeIntegrator.C index 3ea40728227e..ebf2b63f1a99 100644 --- a/framework/src/timeintegrators/TimeIntegrator.C +++ b/framework/src/timeintegrators/TimeIntegrator.C @@ -53,7 +53,8 @@ TimeIntegrator::TimeIntegrator(const InputParameters & parameters) _var_restriction(declareRestartableData( "var_restriction", !getParam>("variables").empty())), _local_indices(declareRestartableData>("local_indices")), - _vars(declareRestartableData>("vars")) + _vars(declareRestartableData>("vars")), + _from_subvector(NumericVector::build(this->comm())) { _fe_problem.setUDotRequested(true); } @@ -104,7 +105,7 @@ TimeIntegrator::init() void TimeIntegrator::solve() { - mooseError("Calling this method is no longer supported"); + mooseError("Calling TimeIntegrator::solve() is no longer supported"); } void @@ -130,17 +131,16 @@ TimeIntegrator::getNumLinearIterationsLastSolve() const } void -TimeIntegrator::copyVector(NumericVector & from, NumericVector & to) +TimeIntegrator::copyVector(const NumericVector & from, NumericVector & to) { if (!_var_restriction) to = from; else { auto to_sub = to.get_subvector(_local_indices); - auto from_sub = from.get_subvector(_local_indices); - *to_sub = *from_sub; + from.create_subvector(*_from_subvector, _local_indices, false); + *to_sub = *_from_subvector; to.restore_subvector(std::move(to_sub), _local_indices); - from.restore_subvector(std::move(from_sub), _local_indices); } } @@ -152,3 +152,12 @@ TimeIntegrator::integratesVar(const unsigned int var_num) const return _vars.count(var_num); } + +void +TimeIntegrator::computeDuDotDu() +{ + const auto coeff = duDotDuCoeff(); + for (const auto i : index_range(_du_dot_du)) + if (integratesVar(i)) + _du_dot_du[i] = coeff / _dt; +} diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index 14adab99773b..28bf58d72e7a 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -77,7 +77,7 @@ MooseVariableData::MooseVariableData(const MooseVariableField::MooseVariableDataFV(const MooseVariableFV(&_sys) ? true : false), _qrule(nullptr) diff --git a/framework/src/variables/MooseVariableDataLinearFV.C b/framework/src/variables/MooseVariableDataLinearFV.C index 94eac24b1361..0ab886a2e23e 100644 --- a/framework/src/variables/MooseVariableDataLinearFV.C +++ b/framework/src/variables/MooseVariableDataLinearFV.C @@ -41,7 +41,7 @@ MooseVariableDataLinearFV::MooseVariableDataLinearFV( _var_num(_var.number()), _assembly(_subproblem.assembly(_tid, var.kind() == Moose::VAR_SOLVER ? sys.number() : 0)), _element_type(element_type), - _time_integrator(_sys.getPossiblyNullTimeIntegrator(_var_num)), + _time_integrator(_sys.queryTimeIntegrator(_var_num)), _elem(elem), _displaced(dynamic_cast(&_sys) ? true : false), _qrule(nullptr) diff --git a/framework/src/variables/MooseVariableField.C b/framework/src/variables/MooseVariableField.C index 784c128521a6..979181084650 100644 --- a/framework/src/variables/MooseVariableField.C +++ b/framework/src/variables/MooseVariableField.C @@ -27,7 +27,7 @@ MooseVariableField::MooseVariableField(const InputParameters & param : MooseVariableFieldBase(parameters), Moose::FunctorBase::type>(name()), MeshChangedInterface(parameters), - _time_integrator(_sys.getPossiblyNullTimeIntegrator(_var_num)) + _time_integrator(_sys.queryTimeIntegrator(_var_num)) { }