Add an option to use a constant DT for RC coefficients time derivative component

modules/navier_stokes/include/fvkernels/INSFVTimeKernel.h

@@ -46,6 +46,16 @@ class INSFVTimeKernel : public FVFunctorTimeKernel, public INSFVMomentumResidual
   /// Whether to contribute to RC coefficients
   const bool _contribute_to_rc_coeffs;
 
+  /// Whether to contribute to RC coefficients with a fixed time step
+  const bool _use_fixed_dt_rc_contrib;
+
+  /// Time at which to start using the fixed time step for RC coefficients
+  const Real _fixed_dt_rc_contrib_start;
+
+  /// Fixed time step to use for RC coefficients
+  const Real & _fixed_dt_rc_contrib;
+
 private:
   using FVFunctorTimeKernel::_current_elem;
+  const Real _zero = 0;
 };

modules/navier_stokes/src/fvkernels/INSFVMomentumTimeDerivative.C

@@ -10,6 +10,7 @@
 #include "INSFVMomentumTimeDerivative.h"
 #include "SystemBase.h"
 #include "NS.h"
+#include "ImplicitEuler.h"
 
 registerMooseObject("NavierStokesApp", INSFVMomentumTimeDerivative);
 
@@ -26,6 +27,11 @@ INSFVMomentumTimeDerivative::validParams()
 INSFVMomentumTimeDerivative::INSFVMomentumTimeDerivative(const InputParameters & params)
   : INSFVTimeKernel(params), _rho(getFunctor<ADReal>(NS::density))
 {
+  // Check time integrator
+  if (_use_fixed_dt_rc_contrib)
+    if (!dynamic_cast<const ImplicitEuler *>(
+            _fe_problem.getNonlinearSystemBase(_sys.number()).queryTimeIntegrator(_var.number())))
+      paramError("fixed_dt_contribution_to_RC", "Only ImplicitEuler time integrator is supported");
 }
 
 void
@@ -38,7 +44,13 @@ INSFVMomentumTimeDerivative::gatherRCData(const Elem & elem)
   const Real a = residual.derivatives()[dof_number];
 
   if (_contribute_to_rc_coeffs)
-    _rc_uo.addToA(&elem, _index, a);
+  {
+    if (_use_fixed_dt_rc_contrib && _fe_problem.time() > _fixed_dt_rc_contrib_start)
+      // Change the time derivative contribution to only depend on the fixed dt
+      _rc_uo.addToA(&elem, _index, a * _dt / _fixed_dt_rc_contrib);
+    else
+      _rc_uo.addToA(&elem, _index, a);
+  }
 
   addResidualAndJacobian(residual, dof_number);
 }

modules/navier_stokes/src/fvkernels/INSFVTimeKernel.C

@@ -18,19 +18,35 @@ INSFVTimeKernel::validParams()
 {
   auto params = FVFunctorTimeKernel::validParams();
   params += INSFVMomentumResidualObject::validParams();
+
+  // These options can be set to facilitate restarts with RC coefficients that do not depend on
+  // the time steps. By default, they are not active
   params.addParam<bool>(
       "contribute_to_rc",
       true,
       "Whether the time derivative term should contribute to Rhie-Chow coefficients");
-  params.addParamNamesToGroup("contribute_to_rc", "Advanced");
+  params.addParam<Real>("start_fixed_dt_contribution_to_RC",
+                        0,
+                        "Time at which to start a fixed dt contribution to the RC coefficients");
+  params.addParam<Real>("fixed_dt_contribution_to_RC",
+                        "Fixed dt for the contribution to the Rhie Chow coefficients");
+  params.declareControllable("fixed_dt_contribution_to_RC");
+  params.addParamNamesToGroup("contribute_to_rc start_fixed_dt_contribution_to_RC", "Advanced");
   return params;
 }
 
 INSFVTimeKernel::INSFVTimeKernel(const InputParameters & params)
   : FVFunctorTimeKernel(params),
     INSFVMomentumResidualObject(*this),
-    _contribute_to_rc_coeffs(getParam<bool>("contribute_to_rc"))
+    _contribute_to_rc_coeffs(getParam<bool>("contribute_to_rc")),
+    _use_fixed_dt_rc_contrib(isParamValid("fixed_dt_contribution_to_RC")),
+    _fixed_dt_rc_contrib_start(getParam<Real>("start_fixed_dt_contribution_to_RC")),
+    _fixed_dt_rc_contrib(_use_fixed_dt_rc_contrib ? getParam<Real>("fixed_dt_contribution_to_RC")
+                                                  : _zero)
 {
+  if (!_use_fixed_dt_rc_contrib && isParamSetByUser("start_fixed_dt_contribution_to_RC"))
+    paramError("start_fixed_dt_contribution_to_RC",
+               "Should not be set unless 'start_fixed_dt_contribution_to_RC' is set");
 }
 
 void

modules/navier_stokes/src/physics/WCNSFVFlowPhysics.C

@@ -34,8 +34,17 @@ WCNSFVFlowPhysics::validParams()
       true,
       "Whether the time derivative term should contribute to the Rhie Chow coefficients. This adds "
       "stabilization, but makes the solution dependent on the time step size");
-  params.addParamNamesToGroup("time_derivative_contributes_to_RC_coefficients characteristic_speed",
-                              "Numerical scheme");
+  params.addParam<Real>(
+      "time_derivative_rc_coef_fixed_dt",
+      "Fixed dt used to make Rhie Chow coefficients independent of the time step before restarts");
+  params.addParam<Real>("time_derivative_rc_coef_fixed_dt_start",
+                        "Time after which to start leveraging the fixed dt to make Rhie Chow "
+                        "coefficients independent of the time step before restarts");
+  params.addParamNamesToGroup("characteristic_speed", "Numerical scheme");
+  params.addParamNamesToGroup("time_derivative_contributes_to_RC_coefficients "
+                              "time_derivative_rc_coef_fixed_dt "
+                              "time_derivative_rc_coef_fixed_dt_start",
+                              "Restart");
 
   // Used for flow mixtures, where one phase is solid / not moving under the action of gravity
   params.addParam<MooseFunctorName>(
@@ -375,6 +384,12 @@ WCNSFVFlowPhysics::addINSMomentumTimeKernels()
   params.set<UserObjectName>("rhie_chow_user_object") = rhieChowUOName();
   params.set<bool>("contribute_to_rc") =
       getParam<bool>("time_derivative_contributes_to_RC_coefficients");
+  if (isParamValid("time_derivative_rc_coef_fixed_dt"))
+    params.set<Real>("fixed_dt_contribution_to_RC") =
+        getParam<Real>("time_derivative_rc_coef_fixed_dt");
+  if (isParamValid("time_derivative_rc_coef_fixed_dt_start"))
+    params.set<Real>("start_fixed_dt_contribution_to_RC") =
+        getParam<Real>("time_derivative_rc_coef_fixed_dt_start");
 
   for (const auto d : make_range(dimension()))
   {

modules/navier_stokes/test/tests/finite_volume/physics/restart/2d_channel_init.i

@@ -0,0 +1,130 @@
+### Thermophysical Properties ###
+rho = 1
+mu = 100
+
+### Simulation parameters
+inlet_velocity = 1
+side_length = 1
+restart_initial_dt = 1
+[Mesh]
+  active = 'gen'
+  [gen]
+    type = GeneratedMeshGenerator
+    dim = 2
+    xmin = 0
+    xmax = 10
+    ymin = 0
+    ymax = ${side_length}
+    nx = 100
+    ny = 20
+  []
+  [fmg_restart]
+    type = FileMeshGenerator
+    file = 2d_channel_init_exodus.e
+    use_for_exodus_restart = true
+  []
+[]
+
+[Physics]
+  [NavierStokes]
+    [Flow]
+      [all_flow]
+        compressibility = 'incompressible'
+
+        density = ${rho}
+        dynamic_viscosity = ${mu}
+
+        inlet_boundaries = 'left'
+        momentum_inlet_types = 'fixed-velocity'
+        momentum_inlet_functors = '${inlet_velocity} 0'
+
+        wall_boundaries = 'top bottom'
+        momentum_wall_types = 'noslip noslip'
+
+        outlet_boundaries = 'right'
+        momentum_outlet_types = 'fixed-pressure'
+        pressure_functors = '0'
+
+        # Make sure restart is perfect by freezing the time step dependence of
+        # the contribution to the RC coefficient by the time derivative kernel
+        time_derivative_rc_coef_fixed_dt = ${restart_initial_dt} # (s)
+        time_derivative_rc_coef_fixed_dt_start = 10
+
+        mass_advection_interpolation = 'upwind'
+        momentum_advection_interpolation = 'upwind'
+      []
+    []
+  []
+[]
+
+[Functions]
+  [grow_dt]
+    type = PiecewiseLinear
+    x = '0 10'
+    y = '0.1 10'
+  []
+  [constant]
+    type = PiecewiseConstant
+    x = '0 10'
+    y = '${restart_initial_dt} ${restart_initial_dt}'
+  []
+[]
+
+[Executioner]
+  type = Transient
+  solve_type = 'NEWTON'
+  petsc_options_iname = '-pc_type -pc_factor_shift_type'
+  petsc_options_value = 'lu NONZERO'
+  nl_abs_tol = 1e-10
+
+  line_search = 'none'
+  [TimeStepper]
+    type = FunctionDT
+    function = 'grow_dt'
+  []
+  # Enough to obtain steady state (no nonlinear iteration needed)
+  end_time = 40
+[]
+
+[Outputs]
+  # Used to set up a restart from exodus file
+  [exodus]
+    type = Exodus
+    execute_on = TIMESTEP_END
+  []
+  # Used to check results
+  csv = true
+[]
+
+[Postprocessors]
+  [min_vel_x]
+    type = ElementExtremeValue
+    variable = 'vel_x'
+    value_type = 'min'
+  []
+  [max_vel_x]
+    type = ElementExtremeValue
+    variable = 'vel_x'
+    value_type = 'max'
+  []
+  [min_vel_y]
+    type = ElementExtremeValue
+    variable = 'vel_y'
+    value_type = 'min'
+  []
+  [max_vel_y]
+    type = ElementExtremeValue
+    variable = 'vel_y'
+    value_type = 'max'
+  []
+  [min_pressure]
+    type = ElementExtremeValue
+    variable = 'pressure'
+    value_type = 'min'
+  []
+  [max_pressure]
+    type = ElementExtremeValue
+    variable = 'pressure'
+    value_type = 'max'
+  []
+[]

modules/navier_stokes/test/tests/finite_volume/physics/restart/gold/2d_channel_init_out.csv

@@ -0,0 +1,12 @@
+time,max_pressure,max_vel_x,max_vel_y,min_pressure,min_vel_x,min_vel_y
+0,0,0,0,0,0,0
+0.1,14212.706788194,1.4919498266697,0.34043711118634,60.299188126278,0.14944554792803,-0.34043711118634
+0.299,14093.89890486,1.492583174832,0.34054313371413,59.70153045387,0.14920089724169,-0.34054313371413
+0.69501,14093.890023708,1.4925835978148,0.34054306059652,59.701492548005,0.14920080194696,-0.34054306059652
+1.4830699,14093.889271499,1.4925835979685,0.34054300989256,59.701492537316,0.14920080190933,-0.34054300989256
+3.051309101,14093.888894536,1.4925835979762,0.34054298441005,59.701492537314,0.14920080190394,-0.34054298441005
+6.17210511099,14093.888705106,1.4925835979801,0.34054297160468,59.701492537314,0.14920080190123,-0.34054297160468
+10,14093.88866976,1.4925835979809,0.34054296921534,59.701492537314,0.14920080190073,-0.34054296921534
+20,14093.889110909,1.4925835979718,0.34054299903646,59.701492537314,0.14920080190704,-0.34054299903646
+30,14093.889110909,1.4925835979718,0.34054299903646,59.701492537314,0.14920080190704,-0.34054299903646
+40,14093.889110909,1.4925835979718,0.34054299903646,59.701492537314,0.14920080190704,-0.34054299903646

modules/navier_stokes/test/tests/finite_volume/physics/restart/gold/transient_from_restart.csv

@@ -0,0 +1,12 @@
+time,max_pressure,max_vel_x,max_vel_y,min_pressure,min_vel_x,min_vel_y
+0,0,0,0,0,0,0
+1,14093.889110909,1.4925835979718,0.34054299903646,59.701492537314,0.14920080190704,-0.34054299903646
+2,14093.889110909,1.4925835979718,0.34054299903646,59.701492537314,0.14920080190704,-0.34054299903646
+3,14093.889110909,1.4925835979718,0.34054299903646,59.701492537314,0.14920080190704,-0.34054299903646
+4,14093.889110909,1.4925835979718,0.34054299903646,59.701492537314,0.14920080190704,-0.34054299903646
+5,14093.889110909,1.4925835979718,0.34054299903646,59.701492537314,0.14920080190704,-0.34054299903646
+6,14093.889110909,1.4925835979718,0.34054299903646,59.701492537314,0.14920080190704,-0.34054299903646
+7,14093.889110909,1.4925835979718,0.34054299903646,59.701492537314,0.14920080190704,-0.34054299903646
+8,14093.889110909,1.4925835979718,0.34054299903646,59.701492537314,0.14920080190704,-0.34054299903646
+9,14093.889110909,1.4925835979718,0.34054299903646,59.701492537314,0.14920080190704,-0.34054299903646
+10,14093.889110909,1.4925835979718,0.34054299903646,59.701492537314,0.14920080190704,-0.34054299903646

modules/navier_stokes/test/tests/finite_volume/physics/restart/tests

@@ -64,4 +64,24 @@
       allow_unused = true
     []
   []
+  [restart_different_dt]
+    requirement = 'The system shall be able to'
+    [init]
+      type = CSVDiff
+      input = '2d_channel_init.i'
+      csvdiff = 2d_channel_init_out.csv'
+      detail = 'initialize a solution for a restart simulation using a fixed time step for the contribution to Rhie Chow coefficients,'
+    []
+    [flat_restart_different_dt]
+      type = CSVDiff
+      input = '2d_channel_init.i'
+      csvdiff = 'transient_from_restart.csv'
+      cli_args = "Physics/NavierStokes/Flow/all_flow/time_derivative_rc_coef_fixed_dt_start=1e8
+                  Mesh/active=fmg_restart
+                  Physics/NavierStokes/Flow/all_flow/initialize_variables_from_mesh_file=true
+                  Outputs/file_base=transient_from_restart
+                  Executioner/TimeStepper/function=constant Executioner/end_time=10"
+      detail = 'and restarting from such solution with a perfect, zero-residual restart.'
+    []
+  []
 []