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 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/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/systems/AuxiliarySystem.h b/framework/include/systems/AuxiliarySystem.h index 60e456054b0a..660d3106e04e 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..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 Number & 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 override { return _undisplaced_system.duDotDu(); } + 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(); } @@ -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: @@ -297,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/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..b3be2f184dd6 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 & duDotDus() { return _du_dot_du; } virtual Number & duDotDotDu() { return _du_dotdot_du; } - virtual const Number & duDotDu() const { return _du_dot_du; } + 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; } @@ -878,18 +878,9 @@ 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(); } - - std::shared_ptr getSharedTimeIntegrator() { return _time_integrator; } + void addTimeIntegrator(const std::string & type, + const std::string & name, + InputParameters & parameters); /// 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 +941,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 * queryTimeIntegrator(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. @@ -990,7 +1003,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) @@ -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/BDF2.h b/framework/include/timeintegrators/BDF2.h index 50ffc9928a6e..5893756dbaae 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: /** @@ -38,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 8b9abfbe0db1..90e403707191 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: /** @@ -42,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/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/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..5176f9130d90 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: /** @@ -37,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/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/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/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/NewmarkBeta.h b/framework/include/timeintegrators/NewmarkBeta.h index 03ac4af36d54..aabfc15376b9 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: /** @@ -41,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 c7aade25760a..21c2896d322e 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() {} @@ -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; @@ -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 = 0; + protected: /** * Gets the number of nonlinear iterations in the most recent solve. @@ -164,6 +179,23 @@ 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(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; NonlinearSystemBase & _nl; @@ -177,11 +209,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; // @@ -200,4 +236,17 @@ 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; + + /// A vector that is used for creating 'from' subvectors in the \p copyVector() method + std::unique_ptr> _from_subvector; }; 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/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/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 d5134e07e5fc..ab2324027043 100644 --- a/framework/src/problems/FEProblemBase.C +++ b/framework/src/problems/FEProblemBase.C @@ -1149,12 +1149,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(); } } @@ -6501,10 +6501,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/restart/DataIO.C b/framework/src/restart/DataIO.C index a041da188f5a..c5dd8abdbdc4 100644 --- a/framework/src/restart/DataIO.C +++ b/framework/src/restart/DataIO.C @@ -326,7 +326,14 @@ dataStore(std::ostream & stream, std::unique_ptr> & v, void * context) { - mooseAssert(v, "Null vector"); + // 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) + return; + mooseAssert(context, "Needs a context of the communicator"); const auto & comm = *static_cast(context); mooseAssert(&comm == &v->comm(), "Inconsistent communicator"); @@ -669,6 +676,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/systems/AuxiliarySystem.C b/framework/src/systems/AuxiliarySystem.C index 55244e0050d7..42e246bc6003 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..b42544c57568 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" @@ -76,7 +77,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) @@ -797,6 +797,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 @@ -1579,6 +1580,55 @@ 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::queryTimeIntegrator(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 = queryTimeIntegrator(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; +} + +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/AStableDirk4.C b/framework/src/timeintegrators/AStableDirk4.C index 5e6cd0af00cd..08bbbe5e9ac3 100644 --- a/framework/src/timeintegrators/AStableDirk4.C +++ b/framework/src/timeintegrators/AStableDirk4.C @@ -92,7 +92,7 @@ AStableDirk4::computeTimeDerivatives() u_dot = *_solution; computeTimeDerivativeHelper(u_dot, _solution_old); u_dot.close(); - _du_dot_du = 1. / _dt; + computeDuDotDu(); } void diff --git a/framework/src/timeintegrators/ActuallyExplicitEuler.C b/framework/src/timeintegrators/ActuallyExplicitEuler.C index 0d3f332a16e4..a02ef4f262f0 100644 --- a/framework/src/timeintegrators/ActuallyExplicitEuler.C +++ b/framework/src/timeintegrators/ActuallyExplicitEuler.C @@ -49,8 +49,7 @@ ActuallyExplicitEuler::computeTimeDerivatives() u_dot = *_solution; computeTimeDerivativeHelper(u_dot, _solution_old); u_dot.close(); - - _du_dot_du = 1.0 / _dt; + computeDuDotDu(); } void diff --git a/framework/src/timeintegrators/BDF2.C b/framework/src/timeintegrators/BDF2.C index 416961d7f432..87f5123df0a1 100644 --- a/framework/src/timeintegrators/BDF2.C +++ b/framework/src/timeintegrators/BDF2.C @@ -50,17 +50,12 @@ BDF2::computeTimeDerivatives() NumericVector & u_dot = *_sys.solutionUDot(); if (_t_step == 1) - { u_dot = *_solution; - _du_dot_du = 1. / _dt; - } else - { u_dot.zero(); - _du_dot_du = _weight[0] / _dt; - } computeTimeDerivativeHelper(u_dot, *_solution, _solution_old, _solution_older); u_dot.close(); + computeDuDotDu(); } void @@ -81,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 a3cfb30eff89..1e77b4957930 100644 --- a/framework/src/timeintegrators/CentralDifference.C +++ b/framework/src/timeintegrators/CentralDifference.C @@ -87,7 +87,7 @@ CentralDifference::computeTimeDerivatives() u_dot.close(); // used for Jacobian calculations - _du_dot_du = 1.0 / (2 * _dt); + computeDuDotDu(); _du_dotdot_du = 1.0 / (_dt * _dt); // Computing udotdot "factor" @@ -114,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 3126bf136bdd..78cc9a5efa81 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 @@ -52,6 +67,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`."); @@ -59,7 +76,7 @@ CrankNicolson::init() // time derivative is assumed to be zero on initial NumericVector & u_dot = *_sys.solutionUDot(); u_dot.zero(); - _du_dot_du = 0; + computeDuDotDu(); // compute residual for the initial time step // Note: we can not directly pass _residual_old in computeResidualTag because @@ -70,7 +87,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 +111,37 @@ 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); +} + +Real +CrankNicolson::duDotDuCoeff() const +{ + return 2; } diff --git a/framework/src/timeintegrators/ExplicitEuler.C b/framework/src/timeintegrators/ExplicitEuler.C index 1639be7cc1ad..26f0cd76d98e 100644 --- a/framework/src/timeintegrators/ExplicitEuler.C +++ b/framework/src/timeintegrators/ExplicitEuler.C @@ -43,8 +43,7 @@ ExplicitEuler::computeTimeDerivatives() u_dot = *_solution; computeTimeDerivativeHelper(u_dot, _solution_old); u_dot.close(); - - _du_dot_du = 1.0 / _dt; + computeDuDotDu(); } void diff --git a/framework/src/timeintegrators/ExplicitRK2.C b/framework/src/timeintegrators/ExplicitRK2.C index 8f8d95889463..d1ae060345a8 100644 --- a/framework/src/timeintegrators/ExplicitRK2.C +++ b/framework/src/timeintegrators/ExplicitRK2.C @@ -53,9 +53,7 @@ ExplicitRK2::computeTimeDerivatives() NumericVector & u_dot = *_sys.solutionUDot(); u_dot = *_solution; computeTimeDerivativeHelper(u_dot, _solution_old, _solution_older); - - _du_dot_du = 1. / _dt; - u_dot.close(); + computeDuDotDu(); } void diff --git a/framework/src/timeintegrators/ExplicitSSPRungeKutta.C b/framework/src/timeintegrators/ExplicitSSPRungeKutta.C index 9e1b54196188..1dc48694012b 100644 --- a/framework/src/timeintegrators/ExplicitSSPRungeKutta.C +++ b/framework/src/timeintegrators/ExplicitSSPRungeKutta.C @@ -69,7 +69,7 @@ void ExplicitSSPRungeKutta::computeTimeDerivatives() { // Only the Jacobian needs to be computed, since the mass matrix needs it - _du_dot_du = 1.0 / (_b[_stage] * _dt); + computeDuDotDu(); } void @@ -207,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 1a625c6d904b..01d5a747b924 100644 --- a/framework/src/timeintegrators/ExplicitTVDRK2.C +++ b/framework/src/timeintegrators/ExplicitTVDRK2.C @@ -56,7 +56,7 @@ ExplicitTVDRK2::computeTimeDerivatives() u_dot = *_solution; computeTimeDerivativeHelper(u_dot, _solution_old, _solution_older); - _du_dot_du = 1. / _dt; + computeDuDotDu(); u_dot.close(); } diff --git a/framework/src/timeintegrators/ImplicitEuler.C b/framework/src/timeintegrators/ImplicitEuler.C index 0cc7df9a456b..5a2a6a447781 100644 --- a/framework/src/timeintegrators/ImplicitEuler.C +++ b/framework/src/timeintegrators/ImplicitEuler.C @@ -32,11 +32,24 @@ 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; + computeDuDotDu(); } void @@ -50,7 +63,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/ImplicitMidpoint.C b/framework/src/timeintegrators/ImplicitMidpoint.C index 7d7daa654ebe..f72077c24538 100644 --- a/framework/src/timeintegrators/ImplicitMidpoint.C +++ b/framework/src/timeintegrators/ImplicitMidpoint.C @@ -45,7 +45,7 @@ ImplicitMidpoint::computeTimeDerivatives() u_dot = *_solution; computeTimeDerivativeHelper(u_dot, _solution_old); u_dot.close(); - _du_dot_du = 1. / _dt; + computeDuDotDu(); } void diff --git a/framework/src/timeintegrators/LStableDirk2.C b/framework/src/timeintegrators/LStableDirk2.C index 5904203e5ea2..06cc76457950 100644 --- a/framework/src/timeintegrators/LStableDirk2.C +++ b/framework/src/timeintegrators/LStableDirk2.C @@ -48,7 +48,7 @@ LStableDirk2::computeTimeDerivatives() u_dot = *_solution; computeTimeDerivativeHelper(u_dot, _solution_old); u_dot.close(); - _du_dot_du = 1. / _dt; + computeDuDotDu(); } void diff --git a/framework/src/timeintegrators/LStableDirk3.C b/framework/src/timeintegrators/LStableDirk3.C index 9fec25a0deca..006c5e1748d4 100644 --- a/framework/src/timeintegrators/LStableDirk3.C +++ b/framework/src/timeintegrators/LStableDirk3.C @@ -67,7 +67,7 @@ LStableDirk3::computeTimeDerivatives() u_dot = *_solution; computeTimeDerivativeHelper(u_dot, _solution_old); u_dot.close(); - _du_dot_du = 1. / _dt; + computeDuDotDu(); } void diff --git a/framework/src/timeintegrators/LStableDirk4.C b/framework/src/timeintegrators/LStableDirk4.C index f144ced1ee8b..df5a8f88f788 100644 --- a/framework/src/timeintegrators/LStableDirk4.C +++ b/framework/src/timeintegrators/LStableDirk4.C @@ -62,7 +62,7 @@ LStableDirk4::computeTimeDerivatives() u_dot = *_solution; computeTimeDerivativeHelper(u_dot, _solution_old); u_dot.close(); - _du_dot_du = 1. / _dt; + computeDuDotDu(); } void diff --git a/framework/src/timeintegrators/NewmarkBeta.C b/framework/src/timeintegrators/NewmarkBeta.C index 4fa100cbdef1..4ff7475fe254 100644 --- a/framework/src/timeintegrators/NewmarkBeta.C +++ b/framework/src/timeintegrators/NewmarkBeta.C @@ -92,7 +92,7 @@ NewmarkBeta::computeTimeDerivatives() // used for Jacobian calculations _du_dotdot_du = 1.0 / _beta / _dt / _dt; - _du_dot_du = _gamma / _beta / _dt; + computeDuDotDu(); } void @@ -118,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 ae095b38324f..ebf2b63f1a99 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; } @@ -33,9 +35,13 @@ 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>>( + "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()), @@ -43,16 +49,68 @@ 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(declareRestartableData( + "var_restriction", !getParam>("variables").empty())), + _local_indices(declareRestartableData>("local_indices")), + _vars(declareRestartableData>("vars")), + _from_subvector(NumericVector::build(this->comm())) { _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)); + } + + _solution_sub = NumericVector::build(_solution->comm()); + _solution_old_sub = NumericVector::build(_solution_old.comm()); +} + void TimeIntegrator::solve() { - _nl.system().solve(); + mooseError("Calling TimeIntegrator::solve() is no longer supported"); +} +void +TimeIntegrator::setNumIterationsLastSolve() +{ _n_nonlinear_iterations = getNumNonlinearIterationsLastSolve(); _n_linear_iterations = getNumLinearIterationsLastSolve(); } @@ -66,8 +124,40 @@ 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(const NumericVector & from, NumericVector & to) +{ + if (!_var_restriction) + to = from; + else + { + auto to_sub = to.get_subvector(_local_indices); + from.create_subvector(*_from_subvector, _local_indices, false); + *to_sub = *_from_subvector; + to.restore_subvector(std::move(to_sub), _local_indices); + } +} - return nonlinear_solver.get_total_linear_iterations(); +bool +TimeIntegrator::integratesVar(const unsigned int var_num) const +{ + if (!_var_restriction) + return true; + + 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/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..28bf58d72e7a 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -77,7 +77,7 @@ MooseVariableData::MooseVariableData(const MooseVariableField::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/MooseVariableDataFV.C b/framework/src/variables/MooseVariableDataFV.C index 6dd13a0c526d..e0840ad3b0ae 100644 --- a/framework/src/variables/MooseVariableDataFV.C +++ b/framework/src/variables/MooseVariableDataFV.C @@ -61,7 +61,7 @@ MooseVariableDataFV::MooseVariableDataFV(const MooseVariableFV(&_sys) ? true : false), _qrule(nullptr) diff --git a/framework/src/variables/MooseVariableDataLinearFV.C b/framework/src/variables/MooseVariableDataLinearFV.C index 5cc60e6c9efb..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.getTimeIntegrator()), + _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 da3685e09254..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.getTimeIntegrator()) + _time_integrator(_sys.queryTimeIntegrator(_var_num)) { } 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(); 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/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."); 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/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; 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..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).getTimeIntegrator(); + 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) || 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.' + [] +[] 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);