diff --git a/framework/include/interfaces/Coupleable.h b/framework/include/interfaces/Coupleable.h index 2def8293bad9..a32737e6c0b0 100644 --- a/framework/include/interfaces/Coupleable.h +++ b/framework/include/interfaces/Coupleable.h @@ -510,19 +510,42 @@ class Coupleable */ virtual const VariableValue & coupledDot(const std::string & var_name, unsigned int comp = 0); + /** + * Residual corresponding to the time derivative of a coupled variable + * Different from time derivative of a coupled variable for explicit solvers + * @param var_name Name of coupled variable + * @param comp Component number for vector of coupled variables + * @return Reference to a VariableValue containing the residual corresponding + * to the time derivative of the coupled variable + */ + virtual const VariableValue & coupledDotResidual(const std::string & var_name, + unsigned int comp = 0); + /** * Second time derivative of a coupled variable * @param var_name Name of coupled variable * @param comp Component number for vector of coupled variables - * @return Reference to a VariableValue containing the time derivative of the coupled variable + * @return Reference to a VariableValue containing the second time derivative of the coupled + * variable */ virtual const VariableValue & coupledDotDot(const std::string & var_name, unsigned int comp = 0); + /** + * Residual corresponding to the second time derivative of a coupled variable + * Different from second time derivative of a coupled variable for explicit solvers + * @param var_name Name of coupled variable + * @param comp Component number for vector of coupled variables + * @ return Reference to a VariableValue containing the residual corresponding + * to the second time derivative of the coupled variable. + */ + virtual const VariableValue & coupledDotDotResidual(const std::string & var_name, + unsigned int comp = 0); + /** * Old time derivative of a coupled variable * @param var_name Name of coupled variable * @param comp Component number for vector of coupled variables - * @return Reference to a VariableValue containing the time derivative of the coupled variable + * @return Reference to a VariableValue containing the old time derivative of the coupled variable */ virtual const VariableValue & coupledDotOld(const std::string & var_name, unsigned int comp = 0); @@ -530,7 +553,8 @@ class Coupleable * Old second time derivative of a coupled variable * @param var_name Name of coupled variable * @param comp Component number for vector of coupled variables - * @return Reference to a VariableValue containing the time derivative of the coupled variable + * @return Reference to a VariableValue containing the old second time derivative of the coupled + * variable */ virtual const VariableValue & coupledDotDotOld(const std::string & var_name, unsigned int comp = 0); @@ -564,7 +588,8 @@ class Coupleable */ virtual const VectorVariableValue & coupledVectorDot(const std::string & var_name, unsigned int comp = 0); - + virtual const VectorVariableValue & coupledVectorDotResidual(const std::string & var_name, + unsigned int comp = 0); /** * Second time derivative of a coupled vector variable * @param var_name Name of coupled vector variable @@ -574,7 +599,8 @@ class Coupleable */ virtual const VectorVariableValue & coupledVectorDotDot(const std::string & var_name, unsigned int comp = 0); - + virtual const VectorVariableValue & coupledVectorDotDotResidual(const std::string & var_name, + unsigned int comp = 0); /** * Old time derivative of a coupled vector variable * @param var_name Name of coupled vector variable @@ -624,6 +650,8 @@ class Coupleable */ virtual const ArrayVariableValue & coupledArrayDot(const std::string & var_name, unsigned int comp = 0); + virtual const ArrayVariableValue & coupledArrayDotResidual(const std::string & var_name, + unsigned int comp = 0); /** * Second time derivative of a coupled array variable @@ -634,6 +662,8 @@ class Coupleable */ virtual const ArrayVariableValue & coupledArrayDotDot(const std::string & var_name, unsigned int comp = 0); + virtual const ArrayVariableValue & coupledArrayDotDotResidual(const std::string & var_name, + unsigned int comp = 0); /** * Old time derivative of a coupled array variable @@ -729,7 +759,8 @@ class Coupleable */ template const T & coupledNodalDot(const std::string & var_name, unsigned int comp = 0); - + template + const T & coupledNodalDotResidual(const std::string & var_name, unsigned int comp = 0); /** * Nodal values of second time derivative of a coupled variable * @param var_name Name of coupled variable @@ -739,7 +770,8 @@ class Coupleable */ virtual const VariableValue & coupledNodalDotDot(const std::string & var_name, unsigned int comp = 0); - + virtual const VariableValue & coupledNodalDotDotResidual(const std::string & var_name, + unsigned int comp = 0); /** * Nodal values of old time derivative of a coupled variable * @param var_name Name of coupled variable diff --git a/framework/include/interfaces/NeighborCoupleable.h b/framework/include/interfaces/NeighborCoupleable.h index 1d2dcb59f54b..fc1c29155eda 100644 --- a/framework/include/interfaces/NeighborCoupleable.h +++ b/framework/include/interfaces/NeighborCoupleable.h @@ -37,6 +37,8 @@ class NeighborCoupleable : public Coupleable unsigned int comp = 0); virtual const VariableValue & coupledNeighborValueDot(const std::string & var_name, unsigned int comp = 0); + virtual const VariableValue & coupledNeighborValueDotResidual(const std::string & var_name, + unsigned int comp = 0); virtual const VariableValue & coupledNeighborValueDotDu(const std::string & var_name, unsigned int comp = 0); virtual const VariableValue & coupledNeighborValueOld(const std::string & var_name, diff --git a/framework/include/interfaces/ScalarCoupleable.h b/framework/include/interfaces/ScalarCoupleable.h index 62459b4de506..9e581aa31bcc 100644 --- a/framework/include/interfaces/ScalarCoupleable.h +++ b/framework/include/interfaces/ScalarCoupleable.h @@ -149,6 +149,17 @@ class ScalarCoupleable */ virtual VariableValue & coupledScalarDot(const std::string & var_name, unsigned int comp = 0); + /** + * Returns the residual corresponding to the time derivative of a scalar coupled variable + * Different from the time derivative of the scalar variable for explicit solvers + * @param var_name Name of coupled variable + * @param comp Component number for vector of coupled variables + * @return Reference to a VariableValue containing the residual corresponding + * to the time derivative for the coupled variable + */ + virtual VariableValue & coupledScalarDotResidual(const std::string & var_name, + unsigned int comp = 0); + /** * Returns the second time derivative of a scalar coupled variable * @param var_name Name of coupled variable @@ -157,6 +168,17 @@ class ScalarCoupleable */ virtual VariableValue & coupledScalarDotDot(const std::string & var_name, unsigned int comp = 0); + /** + * Returns the residual corresponding to the second time derivative of a scalar coupled variable + * Different from the second time derivative of the scalar coupled variable for explicit solvers + * @param var_name Name of coupled variable + * @param comp Component number for vector of coupled variables + * @return Reference to a VariableValue containing the residual corresponding + * to the second time derivative of the coupled variable + */ + virtual VariableValue & coupledScalarDotDotResidual(const std::string & var_name, + unsigned int comp = 0); + /** * Returns the old time derivative of a scalar coupled variable * @param var_name Name of coupled variable diff --git a/framework/include/problems/FEProblemBase.h b/framework/include/problems/FEProblemBase.h index 5fb2dfb61272..403a1ff635b8 100644 --- a/framework/include/problems/FEProblemBase.h +++ b/framework/include/problems/FEProblemBase.h @@ -1731,6 +1731,12 @@ class FEProblemBase : public SubProblem, public Restartable using SubProblem::haveADObjects; void haveADObjects(bool have_ad_objects) override; + /// If solution older than 2 time steps is required, set solution_state to the number of old timesteps required. + virtual void setSolutionState(unsigned int solution_state) { _solution_state = solution_state; } + + /// Get the number of old solution states to be stored + virtual unsigned int getSolutionState() const { return _solution_state; } + // Whether or not we should solve this system bool shouldSolve() const { return _solve; } @@ -2141,6 +2147,9 @@ class FEProblemBase : public SubProblem, public Restartable /// Whether old solution second time derivative needs to be stored bool _u_dotdot_old_requested; + /// Number of old solution states to be stored; + unsigned int _solution_state; + friend class AuxiliarySystem; friend class NonlinearSystemBase; friend class MooseEigenSystem; diff --git a/framework/include/systems/AuxiliarySystem.h b/framework/include/systems/AuxiliarySystem.h index 226e37c513f0..49bee35f09eb 100644 --- a/framework/include/systems/AuxiliarySystem.h +++ b/framework/include/systems/AuxiliarySystem.h @@ -179,6 +179,8 @@ class AuxiliarySystem : public SystemBase, public PerfGraphInterface virtual System & system() override { return _sys; } virtual const System & system() const override { return _sys; } + virtual NumericVector * solutionState(unsigned int i) override; + virtual void setPreviousNewtonSolution(); void setScalarVariableCoupleableTags(ExecFlagType type); @@ -212,8 +214,6 @@ class AuxiliarySystem : public SystemBase, public PerfGraphInterface NumericVector & _serialized_solution; /// Solution vector of the previous nonlinear iterate NumericVector * _solution_previous_nl; - /// Time integrator - std::shared_ptr _time_integrator; /// solution vector for u^dot NumericVector * _u_dot; /// solution vector for u^dotdot @@ -224,6 +224,8 @@ class AuxiliarySystem : public SystemBase, public PerfGraphInterface /// Old solution vector for u^dotdot NumericVector * _u_dotdot_old; + std::vector *> _solution_state; + /// Whether or not a copy of the residual needs to be made bool _need_serialized_solution; diff --git a/framework/include/systems/DisplacedSystem.h b/framework/include/systems/DisplacedSystem.h index 3b25e8045795..df4234d88381 100644 --- a/framework/include/systems/DisplacedSystem.h +++ b/framework/include/systems/DisplacedSystem.h @@ -209,6 +209,11 @@ class DisplacedSystem : public SystemBase return _undisplaced_system.getMatrix(tag); } + virtual NumericVector * solutionState(unsigned int i) override + { + return _undisplaced_system.solutionState(i); + } + virtual TransientExplicitSystem & sys() { return _sys; } virtual System & system() override; diff --git a/framework/include/systems/NonlinearSystemBase.h b/framework/include/systems/NonlinearSystemBase.h index fc035594a7a3..4e667675c8e1 100644 --- a/framework/include/systems/NonlinearSystemBase.h +++ b/framework/include/systems/NonlinearSystemBase.h @@ -617,10 +617,14 @@ class NonlinearSystemBase : public SystemBase, public PerfGraphInterface return _solution_previous_nl; } + virtual NumericVector * solutionState(unsigned int i) override; + virtual void setSolutionUDotOld(const NumericVector & u_dot_old); virtual void setSolutionUDotDotOld(const NumericVector & u_dotdot_old); + virtual void setSolutionState(const std::vector> & solution_state); + virtual void setPreviousNewtonSolution(const NumericVector & soln); virtual TagID timeVectorTag() override { return _Re_time_tag; } @@ -931,4 +935,6 @@ class NonlinearSystemBase : public SystemBase, public PerfGraphInterface /// The required size of the derivative storage array std::size_t _required_derivative_size; #endif + + std::vector *> _solution_state; }; diff --git a/framework/include/systems/SystemBase.h b/framework/include/systems/SystemBase.h index 24a93f46a0e1..5cbedca55524 100644 --- a/framework/include/systems/SystemBase.h +++ b/framework/include/systems/SystemBase.h @@ -215,6 +215,7 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface virtual NumericVector & solution() = 0; virtual NumericVector & solutionOld() = 0; virtual NumericVector & solutionOlder() = 0; + virtual NumericVector * solutionState(unsigned int i) = 0; virtual NumericVector * solutionPreviousNewton() = 0; virtual const NumericVector & solution() const = 0; virtual const NumericVector & solutionOld() const = 0; @@ -861,6 +862,9 @@ class SystemBase : public libMesh::ParallelObject, public ConsoleStreamInterface /// Map variable number to its pointer std::vector> _numbered_vars; + std::vector *> _saved_solution_state; + unsigned int _solution_state_size; + /// Storage for MooseVariable objects MooseObjectWarehouseBase _variable_warehouse; diff --git a/framework/include/variables/MooseVariableData.h b/framework/include/variables/MooseVariableData.h index 2747ffad6dbc..4e97f14b54a7 100644 --- a/framework/include/variables/MooseVariableData.h +++ b/framework/include/variables/MooseVariableData.h @@ -317,8 +317,12 @@ class MooseVariableData const FieldVariableValue & uDot() const; + const FieldVariableValue & uDotResidual() const; + const FieldVariableValue & uDotDot() const; + const FieldVariableValue & uDotDotResidual() const; + const FieldVariableValue & uDotOld() const; const FieldVariableValue & uDotDotOld() const; @@ -341,6 +345,8 @@ class MooseVariableData const MooseArray & nodalValueArray(Moose::SolutionState state) const; const OutputType & nodalValueDot() const; const OutputType & nodalValueDotDot() const; + const OutputType & nodalValueDotResidual() const; + const OutputType & nodalValueDotDotResidual() const; const OutputType & nodalValueDotOld() const; const OutputType & nodalValueDotDotOld() const; const OutputType & nodalValueDuDotDu() const; @@ -429,8 +435,10 @@ class MooseVariableData const DoFValue & dofValuesOlder() const; const DoFValue & dofValuesPreviousNL() const; const DoFValue & dofValuesDot() const; + const DoFValue & dofValuesDotResidual() const; const DoFValue & dofValuesDotOld() const; const DoFValue & dofValuesDotDot() const; + const DoFValue & dofValuesDotDotResidual() const; const DoFValue & dofValuesDotDotOld() const; const MooseArray & dofValuesDuDotDu() const; const MooseArray & dofValuesDuDotDotDu() const; @@ -547,6 +555,10 @@ class MooseVariableData OutputType _nodal_value_dot; /// nodal values of u_dotdot OutputType _nodal_value_dotdot; + /// nodal values of u_dot_residual + OutputType _nodal_value_dot_residual; + /// nodal values of u_dotdot_residual + OutputType _nodal_value_dotdot_residual; /// nodal values of u_dot_old OutputType _nodal_value_dot_old; /// nodal values of u_dotdot_old @@ -572,6 +584,8 @@ class MooseVariableData mutable bool _need_u_dot; mutable bool _need_ad_u_dot; mutable bool _need_u_dotdot; + mutable bool _need_u_dot_residual; + mutable bool _need_u_dotdot_residual; mutable bool _need_u_dot_old; mutable bool _need_u_dotdot_old; mutable bool _need_du_dot_du; @@ -610,6 +624,8 @@ class MooseVariableData mutable bool _need_dof_values_previous_nl; mutable bool _need_dof_values_dot; mutable bool _need_dof_values_dotdot; + mutable bool _need_dof_values_dot_residual; + mutable bool _need_dof_values_dotdot_residual; mutable bool _need_dof_values_dot_old; mutable bool _need_dof_values_dotdot_old; mutable bool _need_dof_du_dot_du; @@ -628,6 +644,10 @@ class MooseVariableData DoFValue _dof_values_dot; /// nodal values of u_dotdot DoFValue _dof_values_dotdot; + /// nodal values of u_dot_residual + DoFValue _dof_values_dot_residual; + /// nodal values of u_dotdot_residual + DoFValue _dof_values_dotdot_residual; /// nodal values of u_dot_old DoFValue _dof_values_dot_old; /// nodal values of u_dotdot_old @@ -680,6 +700,12 @@ class MooseVariableData /// u_dotdot (second time derivative) FieldVariableValue _u_dotdot, _u_dotdot_bak; + /// u_dot_residual (residual corresponding to first time derivative) + FieldVariableValue _u_dot_residual, _u_dot_residual_bak; + + /// u_dotdot_residual (residual corresponding to second time derivative) + FieldVariableValue _u_dotdot_residual, _u_dotdot_residual_bak; + /// u_dot_old (time derivative) FieldVariableValue _u_dot_old, _u_dot_old_bak; diff --git a/framework/include/variables/MooseVariableFE.h b/framework/include/variables/MooseVariableFE.h index b52cbb1ef300..25ca151e9449 100644 --- a/framework/include/variables/MooseVariableFE.h +++ b/framework/include/variables/MooseVariableFE.h @@ -336,6 +336,8 @@ class MooseVariableFE : public MooseVariableField /// element dots const FieldVariableValue & uDot() const { return _element_data->uDot(); } const FieldVariableValue & uDotDot() const { return _element_data->uDotDot(); } + const FieldVariableValue & uDotResidual() const { return _element_data->uDotResidual(); } + const FieldVariableValue & uDotDotResidual() const { return _element_data->uDotDotResidual(); } const FieldVariableValue & uDotOld() const { return _element_data->uDotOld(); } const FieldVariableValue & uDotDotOld() const { return _element_data->uDotDotOld(); } const VariableValue & duDotDu() const { return _element_data->duDotDu(); } @@ -398,10 +400,12 @@ class MooseVariableFE : public MooseVariableField { return _neighbor_data->curlSln(Moose::Current); } + const FieldVariableCurl & curlSlnOldNeighbor() const { return _neighbor_data->curlSln(Moose::Old); } + const FieldVariableCurl & curlSlnOlderNeighbor() const { return _neighbor_data->curlSln(Moose::Older); @@ -410,6 +414,11 @@ class MooseVariableFE : public MooseVariableField /// neighbor dots const FieldVariableValue & uDotNeighbor() const { return _neighbor_data->uDot(); } const FieldVariableValue & uDotDotNeighbor() const { return _neighbor_data->uDotDot(); } + const FieldVariableValue & uDotNeighborResidual() const { return _neighbor_data->uDotResidual(); } + const FieldVariableValue & uDotDotNeighborResidual() const + { + return _neighbor_data->uDotDotResidual(); + } const FieldVariableValue & uDotOldNeighbor() const { return _neighbor_data->uDotOld(); } const FieldVariableValue & uDotDotOldNeighbor() const { return _neighbor_data->uDotDotOld(); } const VariableValue & duDotDuNeighbor() const { return _neighbor_data->duDotDu(); } @@ -500,11 +509,15 @@ class MooseVariableFE : public MooseVariableField const DoFValue & dofValuesOlderNeighbor(); const DoFValue & dofValuesPreviousNLNeighbor(); const DoFValue & dofValuesDot(); + const DoFValue & dofValuesDotResidual(); const DoFValue & dofValuesDotNeighbor(); + const DoFValue & dofValuesDotNeighborResidual(); const DoFValue & dofValuesDotOld(); const DoFValue & dofValuesDotOldNeighbor(); const DoFValue & dofValuesDotDot(); + const DoFValue & dofValuesDotDotResidual(); const DoFValue & dofValuesDotDotNeighbor(); + const DoFValue & dofValuesDotDotNeighborResidual(); const DoFValue & dofValuesDotDotOld(); const DoFValue & dofValuesDotDotOldNeighbor(); const MooseArray & dofValuesDuDotDu(); @@ -573,7 +586,9 @@ class MooseVariableFE : public MooseVariableField const OutputType & nodalValueOlder(); const OutputType & nodalValuePreviousNL(); const OutputType & nodalValueDot(); + const OutputType & nodalValueDotResidual(); const OutputType & nodalValueDotDot(); + const OutputType & nodalValueDotDotResidual(); const OutputType & nodalValueDotOld(); const OutputType & nodalValueDotDotOld(); const OutputType & nodalValueDuDotDu(); @@ -583,7 +598,9 @@ class MooseVariableFE : public MooseVariableField const OutputType & nodalValueOlderNeighbor(); const OutputType & nodalValuePreviousNLNeighbor(); const OutputType & nodalValueDotNeighbor(); + const OutputType & nodalValueDotNeighborResidual(); const OutputType & nodalValueDotDotNeighbor(); + const OutputType & nodalValueDotDotNeighborResidual(); const OutputType & nodalValueDotOldNeighbor(); const OutputType & nodalValueDotDotOldNeighbor(); const OutputType & nodalValueDuDotDuNeighbor(); diff --git a/framework/include/variables/MooseVariableInterface.h b/framework/include/variables/MooseVariableInterface.h index 6cc04eea0a05..47e0ae091aca 100644 --- a/framework/include/variables/MooseVariableInterface.h +++ b/framework/include/variables/MooseVariableInterface.h @@ -82,6 +82,11 @@ class MooseVariableInterface */ virtual const typename OutputTools::VariableValue & dot(); + /** + * The residual correponding to the time derivative of the variable this object is operating on. + */ + virtual const typename OutputTools::VariableValue & dotResidual(); + /** * The second time derivative of the variable this object is operating on. * @@ -89,6 +94,12 @@ class MooseVariableInterface */ virtual const typename OutputTools::VariableValue & dotDot(); + /** + * The residual corresponding to the second time derivative of the variable this object is + * operating on. + */ + virtual const typename OutputTools::VariableValue & dotDotResidual(); + /** * The old time derivative of the variable this object is operating on. * diff --git a/framework/include/variables/MooseVariableScalar.h b/framework/include/variables/MooseVariableScalar.h index c0b1466d1498..fb5c069ef4b0 100644 --- a/framework/include/variables/MooseVariableScalar.h +++ b/framework/include/variables/MooseVariableScalar.h @@ -24,6 +24,9 @@ class MooseVariableScalar; template <> InputParameters validParams(); +class Assembly; +class TimeIntegrator; + /** * Class for scalar variables (they are different). */ @@ -75,6 +78,18 @@ class MooseVariableScalar : public MooseVariableBase "set uDotRequested() to true in FEProblemBase before requesting `u_dot`."); } + VariableValue & uDotResidual() + { + if (_sys.solutionUDot()) + { + _need_u_dot_residual = true; + return _u_dot_residual; + } + else + mooseError("MooseVariableScalar: Time derivative of solution (`u_dot`) is not stored. Please " + "set uDotRequested() to true in FEProblemBase before requesting `u_dot`."); + } + VariableValue & uDotDot() { if (_sys.solutionUDotDot()) @@ -88,6 +103,19 @@ class MooseVariableScalar : public MooseVariableBase "`u_dotdot`."); } + VariableValue & uDotDotResidual() + { + if (_sys.solutionUDotDot()) + { + _need_u_dotdot_residual = true; + return _u_dotdot_residual; + } + else + mooseError("MooseVariableScalar: Second time derivative of solution (`u_dotdot`) is not " + "stored. Please set uDotDotRequested() to true in FEProblemBase before requesting " + "`u_dotdot`."); + } + VariableValue & uDotOld() { if (_sys.solutionUDotOld()) @@ -155,14 +183,18 @@ class MooseVariableScalar : public MooseVariableBase std::vector _need_matrix_tag_u; VariableValue _u_dot; + VariableValue _u_dot_residual; VariableValue _u_dotdot; + VariableValue _u_dotdot_residual; VariableValue _u_dot_old; VariableValue _u_dotdot_old; VariableValue _du_dot_du; VariableValue _du_dotdot_du; bool _need_u_dot; + bool _need_u_dot_residual; bool _need_u_dotdot; + bool _need_u_dotdot_residual; bool _need_u_dot_old; bool _need_u_dotdot_old; bool _need_du_dot_du; @@ -175,6 +207,9 @@ class MooseVariableScalar : public MooseVariableBase /// The scalar solution with derivative information ADVariableValue _dual_u; + /// A pointer to TimeIntegrator. nullptr if _sys is not a NonlinearSystemBase + TimeIntegrator * _time_integrator; + private: /** * Adds derivative information to the scalar variable value arrays diff --git a/framework/src/interfaces/Coupleable.C b/framework/src/interfaces/Coupleable.C index 79a4151b123d..76c81e3a008c 100644 --- a/framework/src/interfaces/Coupleable.C +++ b/framework/src/interfaces/Coupleable.C @@ -658,6 +658,34 @@ Coupleable::coupledDot(const std::string & var_name, unsigned int comp) } } +const VariableValue & +Coupleable::coupledDotResidual(const std::string & var_name, unsigned int comp) +{ + checkVar(var_name); + if (!isCoupled(var_name)) // Return default 0 + return _default_value_zero; + + validateExecutionerType(var_name, "coupledDotResidual"); + MooseVariable * var = getVar(var_name, comp); + if (var == NULL) + mooseError("Call corresponding vector variable method"); + + if (!_coupleable_neighbor) + { + if (_c_nodal) + return var->dofValuesDotResidual(); + else + return var->uDotResidual(); + } + else + { + if (_c_nodal) + return var->dofValuesDotNeighborResidual(); + else + return var->uDotNeighborResidual(); + } +} + const VariableValue & Coupleable::coupledDotDot(const std::string & var_name, unsigned int comp) { @@ -680,6 +708,34 @@ Coupleable::coupledDotDot(const std::string & var_name, unsigned int comp) } } +const VariableValue & +Coupleable::coupledDotDotResidual(const std::string & var_name, unsigned int comp) +{ + checkVar(var_name); + if (!isCoupled(var_name)) // Return default 0 + return _default_value_zero; + + validateExecutionerType(var_name, "coupledDotDotResidual"); + MooseVariable * var = getVar(var_name, comp); + if (var == NULL) + mooseError("Call corresponding vector variable method"); + + if (!_coupleable_neighbor) + { + if (_c_nodal) + return var->dofValuesDotDotResidual(); + else + return var->uDotDotResidual(); + } + else + { + if (_c_nodal) + return var->dofValuesDotDotNeighborResidual(); + else + return var->uDotDotNeighborResidual(); + } +} + const VariableValue & Coupleable::coupledDotOld(const std::string & var_name, unsigned int comp) { @@ -737,6 +793,19 @@ Coupleable::coupledVectorDot(const std::string & var_name, unsigned int comp) return var->uDotNeighbor(); } +const VectorVariableValue & +Coupleable::coupledVectorDotResidual(const std::string & var_name, unsigned int comp) +{ + VectorMooseVariable * var = getVectorVar(var_name, comp); + if (!var) + return _default_vector_value_zero; + checkFuncType(var_name, VarType::Dot, FuncAge::Curr); + + if (!_coupleable_neighbor) + return var->uDotResidual(); + return var->uDotNeighborResidual(); +} + const VectorVariableValue & Coupleable::coupledVectorDotDot(const std::string & var_name, unsigned int comp) { @@ -750,6 +819,19 @@ Coupleable::coupledVectorDotDot(const std::string & var_name, unsigned int comp) return var->uDotDotNeighbor(); } +const VectorVariableValue & +Coupleable::coupledVectorDotDotResidual(const std::string & var_name, unsigned int comp) +{ + VectorMooseVariable * var = getVectorVar(var_name, comp); + if (!var) + return _default_vector_value_zero; + checkFuncType(var_name, VarType::Dot, FuncAge::Curr); + + if (!_coupleable_neighbor) + return var->uDotDotResidual(); + return var->uDotDotNeighborResidual(); +} + const VectorVariableValue & Coupleable::coupledVectorDotOld(const std::string & var_name, unsigned int comp) { @@ -824,6 +906,28 @@ Coupleable::coupledArrayDot(const std::string & var_name, unsigned int comp) } } +const ArrayVariableValue & +Coupleable::coupledArrayDotResidual(const std::string & var_name, unsigned int comp) +{ + ArrayMooseVariable * var = getArrayVar(var_name, comp); + if (!var) + return _default_array_value_zero; + checkFuncType(var_name, VarType::Dot, FuncAge::Curr); + + if (!_coupleable_neighbor) + { + if (_c_nodal) + return var->dofValuesDotResidual(); + return var->uDotResidual(); + } + else + { + if (_c_nodal) + return var->dofValuesDotNeighborResidual(); + return var->uDotNeighborResidual(); + } +} + const ArrayVariableValue & Coupleable::coupledArrayDotDot(const std::string & var_name, unsigned int comp) { @@ -846,6 +950,28 @@ Coupleable::coupledArrayDotDot(const std::string & var_name, unsigned int comp) } } +const ArrayVariableValue & +Coupleable::coupledArrayDotDotResidual(const std::string & var_name, unsigned int comp) +{ + ArrayMooseVariable * var = getArrayVar(var_name, comp); + if (!var) + return _default_array_value_zero; + checkFuncType(var_name, VarType::Dot, FuncAge::Curr); + + if (!_coupleable_neighbor) + { + if (_c_nodal) + return var->dofValuesDotDotResidual(); + return var->uDotDotResidual(); + } + else + { + if (_c_nodal) + return var->dofValuesDotDotNeighborResidual(); + return var->uDotDotNeighborResidual(); + } +} + const ArrayVariableValue & Coupleable::coupledArrayDotOld(const std::string & var_name, unsigned int comp) { @@ -1274,6 +1400,27 @@ Coupleable::coupledNodalDot(const std::string & var_name, unsigned int comp) mooseError("Neighbor version not implemented"); } +template +const T & +Coupleable::coupledNodalDotResidual(const std::string & var_name, unsigned int comp) +{ + checkVar(var_name); + static const T zero = 0; + if (!isCoupled(var_name)) // Return default 0 + return zero; + + validateExecutionerType(var_name, "coupledNodalDotResidual"); + coupledCallback(var_name, false); + MooseVariableFE * var = getVarHelper(var_name, comp); + if (var == NULL) + mooseError("Call corresponding vector variable method"); + + if (!_coupleable_neighbor) + return var->nodalValueDotResidual(); + else + mooseError("Neighbor version not implemented"); +} + const VariableValue & Coupleable::coupledNodalDotDot(const std::string & var_name, unsigned int comp) { @@ -1287,6 +1434,25 @@ Coupleable::coupledNodalDotDot(const std::string & var_name, unsigned int comp) return var->dofValuesDotDotNeighbor(); } +const VariableValue & +Coupleable::coupledNodalDotDotResidual(const std::string & var_name, unsigned int comp) +{ + checkVar(var_name); + if (!isCoupled(var_name)) // Return default 0 + return _default_value_zero; + + validateExecutionerType(var_name, "coupledNodalDotDotResidual"); + coupledCallback(var_name, false); + MooseVariable * var = getVar(var_name, comp); + if (var == NULL) + mooseError("Call corresponding vector variable method"); + + if (!_coupleable_neighbor) + return var->dofValuesDotDotResidual(); + else + return var->dofValuesDotDotNeighborResidual(); +} + const VariableValue & Coupleable::coupledNodalDotOld(const std::string & var_name, unsigned int comp) { @@ -1614,5 +1780,24 @@ Coupleable::coupledNodalValuePreviousNL(const std::string & var unsigned int comp); template const Real & Coupleable::coupledNodalDot(const std::string & var_name, unsigned int comp); +template const Real & Coupleable::coupledNodalDotResidual(const std::string & var_name, + unsigned int comp); template const RealVectorValue & Coupleable::coupledNodalDot(const std::string & var_name, unsigned int comp); + +template const RealVectorValue & +Coupleable::coupledNodalDotResidual(const std::string & var_name, + unsigned int comp); + +template const Real & +Coupleable::adCoupledNodalValueTemplate(const std::string & var_name, + unsigned int comp); +template const RealVectorValue & +Coupleable::adCoupledNodalValueTemplate(const std::string & var_name, + unsigned int comp); +template const DualReal & +Coupleable::adCoupledNodalValueTemplate(const std::string & var_name, + unsigned int comp); +template const libMesh::VectorValue & +Coupleable::adCoupledNodalValueTemplate(const std::string & var_name, + unsigned int comp); diff --git a/framework/src/interfaces/NeighborCoupleable.C b/framework/src/interfaces/NeighborCoupleable.C index 7b5b8598cbd5..cb00be882858 100644 --- a/framework/src/interfaces/NeighborCoupleable.C +++ b/framework/src/interfaces/NeighborCoupleable.C @@ -43,6 +43,16 @@ NeighborCoupleable::coupledNeighborValueDot(const std::string & var_name, unsign return var->uDotNeighbor(); } +const VariableValue & +NeighborCoupleable::coupledNeighborValueDotResidual(const std::string & var_name, unsigned int comp) +{ + MooseVariable * var = getVar(var_name, comp); + if (_neighbor_nodal) + return var->dofValuesDotNeighborResidual(); + else + return var->uDotNeighborResidual(); +} + const VariableValue & NeighborCoupleable::coupledNeighborValueDotDu(const std::string & var_name, unsigned int comp) { diff --git a/framework/src/interfaces/ScalarCoupleable.C b/framework/src/interfaces/ScalarCoupleable.C index 794bac08f407..5c990419d9cc 100644 --- a/framework/src/interfaces/ScalarCoupleable.C +++ b/framework/src/interfaces/ScalarCoupleable.C @@ -236,6 +236,15 @@ ScalarCoupleable::coupledScalarDot(const std::string & var_name, unsigned int co return var->uDot(); } +VariableValue & +ScalarCoupleable::coupledScalarDotResidual(const std::string & var_name, unsigned int comp) +{ + checkVar(var_name); + validateExecutionerType(var_name, "coupledScalarDotResidual"); + MooseVariableScalar * var = getScalarVar(var_name, comp); + return var->uDotResidual(); +} + VariableValue & ScalarCoupleable::coupledScalarDotDot(const std::string & var_name, unsigned int comp) { @@ -245,6 +254,15 @@ ScalarCoupleable::coupledScalarDotDot(const std::string & var_name, unsigned int return var->uDotDot(); } +VariableValue & +ScalarCoupleable::coupledScalarDotDotResidual(const std::string & var_name, unsigned int comp) +{ + checkVar(var_name); + validateExecutionerType(var_name, "coupledScalarDotDotResidual"); + MooseVariableScalar * var = getScalarVar(var_name, comp); + return var->uDotDotResidual(); +} + VariableValue & ScalarCoupleable::coupledScalarDotOld(const std::string & var_name, unsigned int comp) { diff --git a/framework/src/nodalkernels/TimeNodalKernel.C b/framework/src/nodalkernels/TimeNodalKernel.C index 80f071ca8e14..28e5b5147120 100644 --- a/framework/src/nodalkernels/TimeNodalKernel.C +++ b/framework/src/nodalkernels/TimeNodalKernel.C @@ -22,7 +22,7 @@ TimeNodalKernel::validParams() InputParameters params = NodalKernel::validParams(); params.set("vector_tags") = "time"; - params.set("matrix_tags") = "system"; + params.set("matrix_tags") = "system time"; return params; } diff --git a/framework/src/problems/DisplacedProblem.C b/framework/src/problems/DisplacedProblem.C index ec9238c0f9c1..d8454139dbcf 100644 --- a/framework/src/problems/DisplacedProblem.C +++ b/framework/src/problems/DisplacedProblem.C @@ -67,6 +67,7 @@ DisplacedProblem::DisplacedProblem(const InputParameters & parameters) _assembly.emplace_back(libmesh_make_unique(_displaced_nl, i)); _displaced_nl.addTimeIntegrator(_mproblem.getNonlinearSystemBase().getSharedTimeIntegrator()); + _displaced_aux.addTimeIntegrator(_mproblem.getAuxiliarySystem().getSharedTimeIntegrator()); if (!_default_ghosting) _mesh.getMesh().remove_ghosting_functor(_mesh.getMesh().default_ghosting()); diff --git a/framework/src/problems/FEProblemBase.C b/framework/src/problems/FEProblemBase.C index c457614ede46..61860278ef6b 100644 --- a/framework/src/problems/FEProblemBase.C +++ b/framework/src/problems/FEProblemBase.C @@ -336,6 +336,7 @@ FEProblemBase::FEProblemBase(const InputParameters & parameters) _u_dotdot_requested(false), _u_dot_old_requested(false), _u_dotdot_old_requested(false), + _solution_state(0), _has_mortar(false), _num_grid_steps(0), _displaced_neighbor_ref_pts("invert_elem_phys use_undisplaced_ref unset", "unset") @@ -350,6 +351,7 @@ FEProblemBase::FEProblemBase(const InputParameters & parameters) _t_step = 0; _dt = 0; _dt_old = _dt; + _solution_state = 0; unsigned int n_threads = libMesh::n_threads(); @@ -5032,7 +5034,8 @@ FEProblemBase::addTimeIntegrator(const std::string & type, _nl->addTimeIntegrator(type, name, parameters); _has_time_integrator = true; - // add vectors to store u_dot, u_dotdot, udot_old and u_dotdot_old if requested by the time + // add vectors to store u_dot, u_dotdot, udot_old, u_dotdot_old and + // solution vectors older than 2 time steps, if requested by the time // integrator _aux->addDotVectors(); _nl->addDotVectors(); diff --git a/framework/src/systems/AuxiliarySystem.C b/framework/src/systems/AuxiliarySystem.C index bdc763dd393e..7b476c9357b1 100644 --- a/framework/src/systems/AuxiliarySystem.C +++ b/framework/src/systems/AuxiliarySystem.C @@ -43,6 +43,7 @@ AuxiliarySystem::AuxiliarySystem(FEProblemBase & subproblem, const std::string & _u_dotdot(NULL), _u_dot_old(NULL), _u_dotdot_old(NULL), + _solution_state(0), _need_serialized_solution(false), _aux_scalar_storage(_app.getExecuteOnEnum()), _nodal_aux_storage(_app.getExecuteOnEnum()), @@ -88,6 +89,14 @@ AuxiliarySystem::addDotVectors() _u_dot_old = &addVector("u_dot_old", true, GHOSTED); if (_fe_problem.uDotDotOldRequested()) _u_dotdot_old = &addVector("u_dotdot_old", true, GHOSTED); + + if (_fe_problem.getSolutionState() > 3) + { + _solution_state_size = _fe_problem.getSolutionState(); + _solution_state.resize(_fe_problem.getSolutionState()); + for (unsigned int i = 0; i < _fe_problem.getSolutionState(); ++i) + _solution_state[i] = &addVector("solution_state_" + std::to_string(i), true, GHOSTED); + } } void @@ -242,7 +251,8 @@ AuxiliarySystem::addTimeIntegrator(const std::string & type, InputParameters & parameters) { parameters.set("_sys") = this; - _time_integrator = _factory.create(type, name, parameters); + std::shared_ptr ti = _factory.create(type, name, parameters); + _time_integrator = ti; } void @@ -637,6 +647,20 @@ AuxiliarySystem::setPreviousNewtonSolution() compute(EXEC_LINEAR); } +NumericVector * +AuxiliarySystem::solutionState(unsigned int i) +{ + if (_fe_problem.getSolutionState() > 3 && i < _fe_problem.getSolutionState()) + return _solution_state[i]; + else if (_fe_problem.getSolutionState() <= 3) + mooseError("AuxiliarySystem: If only current, old or older solution is required, please use " + "solution(), solutionOld() or solutionOlder(), respectively."); + else + mooseError("AuxiliarySystem: Requested solutionState is not saved, please make sure solution " + "state is less than" + + std::to_string(_fe_problem.getSolutionState())); +} + template void AuxiliarySystem::computeElementalVarsHelper( diff --git a/framework/src/systems/NonlinearSystemBase.C b/framework/src/systems/NonlinearSystemBase.C index b963fa265e77..26f91cfcd690 100644 --- a/framework/src/systems/NonlinearSystemBase.C +++ b/framework/src/systems/NonlinearSystemBase.C @@ -173,8 +173,9 @@ NonlinearSystemBase::NonlinearSystemBase(FEProblemBase & fe_problem, _resid_vs_jac_scaling_param(0) #ifndef MOOSE_SPARSE_AD , - _required_derivative_size(0) + _required_derivative_size(0), #endif + _solution_state(0) { getResidualNonTimeVector(); // Don't need to add the matrix - it already exists (for now) @@ -226,6 +227,14 @@ NonlinearSystemBase::addDotVectors() _u_dotdot = &addVector("u_dotdot", true, GHOSTED); if (_fe_problem.uDotDotOldRequested()) _u_dotdot_old = &addVector("u_dotdot_old", true, GHOSTED); + + if (_fe_problem.getSolutionState() > 3) + { + _solution_state_size = _fe_problem.getSolutionState(); + _solution_state.resize(_fe_problem.getSolutionState()); + for (unsigned int i = 0; i < _fe_problem.getSolutionState(); ++i) + _solution_state[i] = &addVector("solution_state_" + std::to_string(i), true, GHOSTED); + } } void @@ -3005,6 +3014,17 @@ NonlinearSystemBase::setSolutionUDotDotOld(const NumericVector & u_dotdo *_u_dotdot_old = u_dotdot_old; } +void +NonlinearSystemBase::setSolutionState(const std::vector> & solution_state) +{ + if (solution_state.size() != _fe_problem.getSolutionState()) + mooseError("NonlinearSystemBase: Size of solution state should match the SolutionState " + "provided by FeProblem."); + + for (unsigned int i = 0; i < _fe_problem.getSolutionState(); ++i) + *_solution_state[i] = solution_state[i]; +} + NumericVector & NonlinearSystemBase::serializedSolution() { @@ -3216,3 +3236,17 @@ NonlinearSystemBase::computeScaling() mooseWarning("The NonlinearSystemBase derived class that is being used does not currently " "support automatic scaling."); } + +NumericVector * +NonlinearSystemBase::solutionState(unsigned int i) +{ + if (_fe_problem.getSolutionState() > 3 && i < _fe_problem.getSolutionState()) + return _solution_state[i]; + else if (_fe_problem.getSolutionState() <= 3) + mooseError("NonlinearSystemBase: If only current, old or older solution is required, please " + "use solution(), solutionOld() or solutionOlder(), respectively."); + else + mooseError("NonlinearSystemBase: Requested solutionState is not saved, please make sure " + "solution state is less than" + + std::to_string(_fe_problem.getSolutionState())); +} diff --git a/framework/src/systems/SystemBase.C b/framework/src/systems/SystemBase.C index 45b4029aff70..fef3c202259c 100644 --- a/framework/src/systems/SystemBase.C +++ b/framework/src/systems/SystemBase.C @@ -104,6 +104,8 @@ SystemBase::SystemBase(SubProblem & subproblem, _max_var_n_dofs_per_elem(0), _max_var_n_dofs_per_node(0), _time_integrator(nullptr), + _saved_solution_state(0), + _solution_state_size(0), _computing_scaling_jacobian(false), _computing_scaling_residual(false), _automatic_scaling(false), @@ -525,6 +527,20 @@ SystemBase::saveOldSolutions() if (solutionUDotDotOld()) *_saved_dotdot_old = *solutionUDotDotOld(); + + if (_solution_state_size > 3) + { + if (_saved_solution_state.size() == 0) + { + _saved_solution_state.resize(_solution_state_size); + for (unsigned int i = 0; i < _solution_state_size; ++i) + _saved_solution_state[i] = + &addVector("save_solution_state_" + std::to_string(i), false, PARALLEL); + } + + for (unsigned int i = 0; i < _solution_state_size; ++i) + *(_saved_solution_state[i]) = *solutionState(i); + } } /** @@ -557,6 +573,15 @@ SystemBase::restoreOldSolutions() removeVector("save_solution_dotdot_old"); _saved_dotdot_old = NULL; } + if (_saved_solution_state.size() != 0 && _solution_state_size > 3) + { + for (unsigned int i = 0; i < _solution_state_size; ++i) + { + *solutionState(i) = *(_saved_solution_state[i]); + removeVector("save_solution_state" + std::to_string(i)); + } + _saved_solution_state.clear(); + } } NumericVector & @@ -1118,6 +1143,10 @@ SystemBase::copySolutionsBackwards() *solutionUDotDotOld() = *solutionUDotDot(); if (solutionPreviousNewton()) *solutionPreviousNewton() = *currentSolution(); + + if (_solution_state_size > 3) + for (unsigned int i = 0; i < _solution_state_size; ++i) + *solutionState(i) = *currentSolution(); } /** @@ -1134,6 +1163,15 @@ SystemBase::copyOldSolutions() *solutionUDotDotOld() = *solutionUDotDot(); if (solutionPreviousNewton()) *solutionPreviousNewton() = *currentSolution(); + + if (_solution_state_size > 3) + { + for (unsigned int i = _solution_state_size - 1; i > 1; --i) + *solutionState(i) = *solutionState(i - 1); + + *solutionState(1) = *currentSolution(); + *solutionState(0) = *currentSolution(); + } } /** @@ -1150,6 +1188,10 @@ SystemBase::restoreSolutions() *solutionUDotDot() = *solutionUDotDotOld(); if (solutionPreviousNewton()) *solutionPreviousNewton() = solutionOld(); + + if (_solution_state_size > 3) + *solutionState(0) = *solutionState(1); + system().update(); } diff --git a/framework/src/variables/MooseVariableData.C b/framework/src/variables/MooseVariableData.C index 4db6b92e009b..542133bb04d9 100644 --- a/framework/src/variables/MooseVariableData.C +++ b/framework/src/variables/MooseVariableData.C @@ -48,6 +48,8 @@ MooseVariableData::MooseVariableData(const MooseVariableFE::MooseVariableData(const MooseVariableFE::uDot() const "uDotRequested() to true in FEProblemBase before requesting `u_dot`."); } +template +const typename MooseVariableData::FieldVariableValue & +MooseVariableData::uDotResidual() const +{ + if (_sys.solutionUDot()) + { + _need_u_dot_residual = true; + return _u_dot_residual; + } + else + mooseError("MooseVariableFE: Time derivative of solution (`u_dot`) is not stored. Please set " + "uDotRequested() to true in FEProblemBase before requesting `u_dot`."); +} + template const typename MooseVariableData::FieldVariableValue & MooseVariableData::uDotDot() const @@ -272,6 +290,21 @@ MooseVariableData::uDotDot() const "`u_dotdot`."); } +template +const typename MooseVariableData::FieldVariableValue & +MooseVariableData::uDotDotResidual() const +{ + if (_sys.solutionUDotDot()) + { + _need_u_dotdot_residual = true; + return _u_dotdot_residual; + } + else + mooseError("MooseVariableFE: Second time derivative of solution (`u_dotdot`) is not stored. " + "Please set uDotDotRequested() to true in FEProblemBase before requesting " + "`u_dotdot`."); +} + template const typename MooseVariableData::FieldVariableValue & MooseVariableData::uDotOld() const @@ -517,6 +550,12 @@ MooseVariableData::computeValues() if (_need_u_dotdot) _u_dotdot.resize(nqp); + if (_need_u_dot_residual) + _u_dot_residual.resize(nqp); + + if (_need_u_dotdot_residual) + _u_dotdot_residual.resize(nqp); + if (_need_u_dot_old) _u_dot_old.resize(nqp); @@ -593,6 +632,12 @@ MooseVariableData::computeValues() if (_need_u_dotdot) _u_dotdot[i] = 0; + if (_need_u_dot_residual) + _u_dot_residual[i] = 0; + + if (_need_u_dotdot_residual) + _u_dotdot_residual[i] = 0; + if (_need_u_dot_old) _u_dot_old[i] = 0; @@ -669,6 +714,12 @@ MooseVariableData::computeValues() if (_need_u_dotdot) _u_dotdot[qp] += phi_local * _dof_values_dotdot[i]; + if (_need_u_dot_residual) + _u_dot_residual[qp] += phi_local * _dof_values_dot_residual[i]; + + if (_need_u_dotdot_residual) + _u_dotdot_residual[qp] += phi_local * _dof_values_dotdot_residual[i]; + if (_need_u_dot_old) _u_dot_old[qp] += phi_local * _dof_values_dot_old[i]; @@ -1095,6 +1146,12 @@ MooseVariableData::computeMonomialValues() if (_need_u_dotdot) _u_dotdot.resize(nqp); + if (_need_u_dot_residual) + _u_dot_residual.resize(nqp); + + if (_need_u_dotdot_residual) + _u_dotdot_residual.resize(nqp); + if (_need_u_dot_old) _u_dot_old.resize(nqp); @@ -1142,6 +1199,10 @@ MooseVariableData::computeMonomialValues() _dof_values_dot.resize(1); if (_need_dof_values_dotdot) _dof_values_dotdot.resize(1); + if (_need_dof_values_dot_residual) + _dof_values_dot_residual.resize(1); + if (_need_dof_values_dotdot_residual) + _dof_values_dotdot_residual.resize(1); if (_need_dof_values_dot_old) _dof_values_dot_old.resize(1); if (_need_dof_values_dotdot_old) @@ -1155,6 +1216,8 @@ MooseVariableData::computeMonomialValues() Real soln_previous_nl = 0; Real u_dot = 0; Real u_dotdot = 0; + Real u_dot_residual = 0; + Real u_dotdot_residual = 0; Real u_dot_old = 0; Real u_dotdot_old = 0; const Real & du_dot_du = _sys.duDotDu(); @@ -1185,9 +1248,15 @@ MooseVariableData::computeMonomialValues() _dof_values_older[0] = soln_older; if (_sys.solutionUDot()) + { u_dot = (*_sys.solutionUDot())(idx); + u_dot_residual = (_time_integrator->uDotResidual())(idx); + } if (_sys.solutionUDotDot()) + { u_dotdot = (*_sys.solutionUDotDot())(idx); + u_dotdot_residual = (_time_integrator->uDotDotResidual())(idx); + } if (_sys.solutionUDotOld()) u_dot_old = (*_sys.solutionUDotOld())(idx); if (_sys.solutionUDotDotOld()) @@ -1198,6 +1267,12 @@ MooseVariableData::computeMonomialValues() if (_need_dof_values_dotdot) _dof_values_dotdot[0] = u_dotdot; + + if (_need_dof_values_dot_residual) + _dof_values_dot_residual[0] = u_dot_residual; + + if (_need_dof_values_dotdot_residual) + _dof_values_dotdot_residual[0] = u_dotdot_residual; } auto phi = (*_current_phi)[0][0]; @@ -1215,6 +1290,12 @@ MooseVariableData::computeMonomialValues() if (_need_u_dotdot) _u_dotdot[0] = phi * u_dotdot; + if (_need_u_dot_residual) + _u_dot_residual[0] = phi * u_dot_residual; + + if (_need_u_dotdot_residual) + _u_dotdot_residual[0] = phi * u_dotdot_residual; + if (_need_u_dot_old) _u_dot_old[0] = phi * u_dot_old; @@ -1249,6 +1330,12 @@ MooseVariableData::computeMonomialValues() if (_need_u_dotdot) _u_dotdot[qp] = _u_dotdot[0]; + if (_need_u_dot_residual) + _u_dot_residual[qp] = _u_dot_residual[0]; + + if (_need_u_dotdot_residual) + _u_dotdot_residual[qp] = _u_dotdot_residual[0]; + if (_need_u_dot_old) _u_dot_old[qp] = _u_dot_old[0]; @@ -1779,6 +1866,35 @@ MooseVariableData::dofValuesDotDot() const "`u_dotdot`."); } +template +const typename MooseVariableData::DoFValue & +MooseVariableData::dofValuesDotResidual() const +{ + if (_sys.solutionUDot()) + { + _need_dof_values_dot_residual = true; + return _dof_values_dot_residual; + } + else + mooseError("MooseVariableFE: Time derivative of solution (`u_dot`) is not stored. Please set " + "uDotRequested() to true in FEProblemBase before requesting `u_dot`."); +} + +template +const typename MooseVariableData::DoFValue & +MooseVariableData::dofValuesDotDotResidual() const +{ + if (_sys.solutionUDotDot()) + { + _need_dof_values_dotdot_residual = true; + return _dof_values_dotdot_residual; + } + else + mooseError("MooseVariableFE: Second time derivative of solution (`u_dotdot`) is not stored. " + "Please set uDotDotRequested() to true in FEProblemBase before requesting " + "`u_dotdot`."); +} + template const typename MooseVariableData::DoFValue & MooseVariableData::dofValuesDotOld() const @@ -2039,6 +2155,49 @@ MooseVariableData::nodalValueDotDot() const "' is not nodal."); } +template +const OutputType & +MooseVariableData::nodalValueDotResidual() const +{ + if (isNodal()) + { + if (_sys.solutionUDot()) + { + _need_dof_values_dot_residual = true; + return _nodal_value_dot_residual; + } + else + mooseError("MooseVariableFE: Time derivative of solution (`u_dot`) is not stored. Please set " + "uDotRequested() to true in FEProblemBase before requesting `u_dot`."); + } + else + mooseError("Nodal values can be requested only on nodal variables, variable '", + _var.name(), + "' is not nodal."); +} + +template +const OutputType & +MooseVariableData::nodalValueDotDotResidual() const +{ + if (isNodal()) + { + if (_sys.solutionUDotDot()) + { + _need_dof_values_dotdot_residual = true; + return _nodal_value_dotdot_residual; + } + else + mooseError("MooseVariableFE: Second time derivative of solution (`u_dotdot`) is not stored. " + "Please set uDotDotRequested() to true in FEProblemBase before requesting " + "`u_dotdot`."); + } + else + mooseError("Nodal values can be requested only on nodal variables, variable '", + _var.name(), + "' is not nodal."); +} + template const OutputType & MooseVariableData::nodalValueDotOld() const @@ -2180,6 +2339,18 @@ MooseVariableData::fetchDoFValues() _dof_values_dotdot.resize(n); _sys.solutionUDotDot()->get(_dof_indices, &_dof_values_dotdot[0]); } + if (_need_u_dot_residual || _need_dof_values_dot_residual) + { + libmesh_assert(_sys.solutionUDot()); + _dof_values_dot_residual.resize(n); + _time_integrator->uDotResidual().get(_dof_indices, &_dof_values_dot_residual[0]); + } + if (_need_u_dotdot_residual || _need_dof_values_dotdot_residual) + { + libmesh_assert(_sys.solutionUDotDot()); + _dof_values_dotdot_residual.resize(n); + _time_integrator->uDotDotResidual().get(_dof_indices, &_dof_values_dotdot_residual[0]); + } if (_need_u_dot_old || _need_dof_values_dot_old) { libmesh_assert(_sys.solutionUDotOld()); @@ -2379,6 +2550,8 @@ MooseVariableData::zeroSizeDofValues() _dof_values_older.resize(0); _dof_values_dot.resize(0); _dof_values_dotdot.resize(0); + _dof_values_dot_residual.resize(0); + _dof_values_dotdot_residual.resize(0); _dof_values_dot_old.resize(0); _dof_values_dotdot_old.resize(0); _dof_du_dot_du.resize(0); @@ -2445,6 +2618,10 @@ MooseVariableData::assignNodalValue() _nodal_value_dot = _dof_values_dot[0]; if (_need_dof_values_dotdot) _nodal_value_dotdot = _dof_values_dotdot[0]; + if (_need_dof_values_dot_residual) + _nodal_value_dot_residual = _dof_values_dot_residual[0]; + if (_need_dof_values_dotdot_residual) + _nodal_value_dotdot_residual = _dof_values_dotdot_residual[0]; if (_need_dof_values_dot_old) _nodal_value_dot_old = _dof_values_dot_old[0]; if (_need_dof_values_dotdot_old) @@ -2480,6 +2657,12 @@ MooseVariableData::assignNodalValue() if (_need_dof_values_dotdot) for (decltype(n) i = 0; i < n; ++i) _nodal_value_dotdot(i) = _dof_values_dotdot[i]; + if (_need_dof_values_dot_residual) + for (decltype(n) i = 0; i < n; ++i) + _nodal_value_dot_residual(i) = _dof_values_dot_residual[i]; + if (_need_dof_values_dotdot_residual) + for (decltype(n) i = 0; i < n; ++i) + _nodal_value_dotdot_residual(i) = _dof_values_dotdot_residual[i]; if (_need_dof_values_dot_old) for (decltype(n) i = 0; i < n; ++i) _nodal_value_dot_old(i) = _dof_values_dot_old[i]; diff --git a/framework/src/variables/MooseVariableFE.C b/framework/src/variables/MooseVariableFE.C index c48f77a221d9..9e79f535afca 100644 --- a/framework/src/variables/MooseVariableFE.C +++ b/framework/src/variables/MooseVariableFE.C @@ -296,6 +296,13 @@ MooseVariableFE::dofValuesDot() return _element_data->dofValuesDot(); } +template +const typename MooseVariableFE::DoFValue & +MooseVariableFE::dofValuesDotResidual() +{ + return _element_data->dofValuesDotResidual(); +} + template const typename MooseVariableFE::DoFValue & MooseVariableFE::dofValuesDotDot() @@ -303,6 +310,13 @@ MooseVariableFE::dofValuesDotDot() return _element_data->dofValuesDotDot(); } +template +const typename MooseVariableFE::DoFValue & +MooseVariableFE::dofValuesDotDotResidual() +{ + return _element_data->dofValuesDotDotResidual(); +} + template const typename MooseVariableFE::DoFValue & MooseVariableFE::dofValuesDotOld() @@ -324,6 +338,13 @@ MooseVariableFE::dofValuesDotNeighbor() return _neighbor_data->dofValuesDot(); } +template +const typename MooseVariableFE::DoFValue & +MooseVariableFE::dofValuesDotNeighborResidual() +{ + return _neighbor_data->dofValuesDotResidual(); +} + template const typename MooseVariableFE::DoFValue & MooseVariableFE::dofValuesDotDotNeighbor() @@ -331,6 +352,13 @@ MooseVariableFE::dofValuesDotDotNeighbor() return _neighbor_data->dofValuesDotDot(); } +template +const typename MooseVariableFE::DoFValue & +MooseVariableFE::dofValuesDotDotNeighborResidual() +{ + return _neighbor_data->dofValuesDotDotResidual(); +} + template const typename MooseVariableFE::DoFValue & MooseVariableFE::dofValuesDotOldNeighbor() @@ -634,6 +662,20 @@ MooseVariableFE::nodalValueDotDot() return _element_data->nodalValueDotDot(); } +template +const OutputType & +MooseVariableFE::nodalValueDotResidual() +{ + return _element_data->nodalValueDotResidual(); +} + +template +const OutputType & +MooseVariableFE::nodalValueDotDotResidual() +{ + return _element_data->nodalValueDotDotResidual(); +} + template const OutputType & MooseVariableFE::nodalValueDotOld() diff --git a/framework/src/variables/MooseVariableInterface.C b/framework/src/variables/MooseVariableInterface.C index 5d3262085266..469aa2cc31b3 100644 --- a/framework/src/variables/MooseVariableInterface.C +++ b/framework/src/variables/MooseVariableInterface.C @@ -155,6 +155,16 @@ MooseVariableInterface::dot() return _variable->uDot(); } +template +const typename OutputTools::VariableValue & +MooseVariableInterface::dotResidual() +{ + if (_nodal) + return _variable->dofValuesDotResidual(); + else + return _variable->uDotResidual(); +} + template const typename OutputTools::VariableValue & MooseVariableInterface::dotDot() @@ -165,6 +175,16 @@ MooseVariableInterface::dotDot() return _variable->uDotDot(); } +template +const typename OutputTools::VariableValue & +MooseVariableInterface::dotDotResidual() +{ + if (_nodal) + return _variable->dofValuesDotDotResidual(); + else + return _variable->uDotDotResidual(); +} + template const typename OutputTools::VariableValue & MooseVariableInterface::dotOld() @@ -195,6 +215,16 @@ MooseVariableInterface::dot() return _variable->uDot(); } +template <> +const VectorVariableValue & +MooseVariableInterface::dotResidual() +{ + if (_nodal) + mooseError("Dofs are scalars while vector variables have vector values. Mismatch"); + else + return _variable->uDotResidual(); +} + template <> const VectorVariableValue & MooseVariableInterface::dotDot() @@ -205,6 +235,16 @@ MooseVariableInterface::dotDot() return _variable->uDotDot(); } +template <> +const VectorVariableValue & +MooseVariableInterface::dotDotResidual() +{ + if (_nodal) + mooseError("Dofs are scalars while vector variables have vector values. Mismatch"); + else + return _variable->uDotDotResidual(); +} + template <> const VectorVariableValue & MooseVariableInterface::dotOld() diff --git a/framework/src/variables/MooseVariableScalar.C b/framework/src/variables/MooseVariableScalar.C index 0deb30e8d444..3faff28af114 100644 --- a/framework/src/variables/MooseVariableScalar.C +++ b/framework/src/variables/MooseVariableScalar.C @@ -12,6 +12,7 @@ #include "SystemBase.h" #include "Assembly.h" #include "SystemBase.h" +#include "TimeIntegrator.h" // libMesh #include "libmesh/numeric_vector.h" @@ -35,13 +36,16 @@ MooseVariableScalar::validParams() MooseVariableScalar::MooseVariableScalar(const InputParameters & parameters) : MooseVariableBase(parameters), _need_u_dot(false), + _need_u_dot_residual(false), _need_u_dotdot(false), + _need_u_dotdot_residual(false), _need_u_dot_old(false), _need_u_dotdot_old(false), _need_du_dot_du(false), _need_du_dotdot_du(false), _need_dual(false), - _need_dual_u(false) + _need_dual_u(false), + _time_integrator(nullptr) { auto num_vector_tags = _sys.subproblem().numVectorTags(); @@ -52,6 +56,8 @@ MooseVariableScalar::MooseVariableScalar(const InputParameters & parameters) _matrix_tag_u.resize(num_matrix_tags); _need_matrix_tag_u.resize(num_matrix_tags); + + _time_integrator = _sys.getTimeIntegrator(); } MooseVariableScalar::~MooseVariableScalar() @@ -61,7 +67,9 @@ MooseVariableScalar::~MooseVariableScalar() _u_older.release(); _u_dot.release(); + _u_dot_residual.release(); _u_dotdot.release(); + _u_dotdot_residual.release(); _u_dot_old.release(); _u_dotdot_old.release(); _du_dot_du.release(); @@ -125,9 +133,15 @@ MooseVariableScalar::reinit(bool reinit_for_derivative_reordering /* = false*/) if (_need_u_dot) _u_dot.resize(n); + if (_need_u_dot_residual) + _u_dot_residual.resize(n); + if (_need_u_dotdot) _u_dotdot.resize(n); + if (_need_u_dotdot_residual) + _u_dotdot_residual.resize(n); + if (_need_u_dot_old) _u_dot_old.resize(n); @@ -177,9 +191,15 @@ MooseVariableScalar::reinit(bool reinit_for_derivative_reordering /* = false*/) if (_need_u_dot) (*u_dot).get(_dof_indices, &_u_dot[0]); + if (_need_u_dot_residual) + (_time_integrator->uDotResidual()).get(_dof_indices, &_u_dot_residual[0]); + if (_need_u_dotdot) (*u_dotdot).get(_dof_indices, &_u_dotdot[0]); + if (_need_u_dotdot_residual) + (_time_integrator->uDotDotResidual()).get(_dof_indices, &_u_dotdot_residual[0]); + if (_need_u_dot_old) (*u_dot_old).get(_dof_indices, &_u_dot_old[0]); @@ -220,9 +240,15 @@ MooseVariableScalar::reinit(bool reinit_for_derivative_reordering /* = false*/) if (_need_u_dot) (*u_dot).get(one_dof_index, &_u_dot[i]); + if (_need_u_dot_residual) + (_time_integrator->uDotResidual()).get(one_dof_index, &_u_dot_residual[i]); + if (_need_u_dotdot) (*u_dotdot).get(one_dof_index, &_u_dotdot[i]); + if (_need_u_dotdot_residual) + (_time_integrator->uDotDotResidual()).get(one_dof_index, &_u_dotdot_residual[i]); + if (_need_u_dot_old) (*u_dot_old).get(one_dof_index, &_u_dot_old[i]); @@ -250,9 +276,15 @@ MooseVariableScalar::reinit(bool reinit_for_derivative_reordering /* = false*/) if (_need_u_dot) _u_dot.resize(i); + if (_need_u_dot_residual) + _u_dot_residual.resize(i); + if (_need_u_dotdot) _u_dotdot.resize(i); + if (_need_u_dotdot_residual) + _u_dotdot_residual.resize(i); + if (_need_u_dot_old) _u_dot_old.resize(i); @@ -277,9 +309,15 @@ MooseVariableScalar::reinit(bool reinit_for_derivative_reordering /* = false*/) if (_need_u_dot) _u_dot[i] = std::numeric_limits::quiet_NaN(); + if (_need_u_dot_residual) + _u_dot_residual[i] = std::numeric_limits::quiet_NaN(); + if (_need_u_dotdot) _u_dotdot[i] = std::numeric_limits::quiet_NaN(); + if (_need_u_dotdot_residual) + _u_dotdot_residual[i] = std::numeric_limits::quiet_NaN(); + if (_need_u_dot_old) _u_dot_old[i] = std::numeric_limits::quiet_NaN();