Merge pull request #28736 from GiudGiud/PR_patch_physics_restart

framework/doc/content/syntax/Physics/index.md

@@ -87,3 +87,29 @@ when you are adding a parameter you need to think about:
 - is this more tied to the discretization of the equation? If so then it likely belong in the derived, user-instantiated,
   `XYZPhysics(CG/DG/HDG/FV/LinearFV)` class.
 
+#### Rules for implementation of Physics with regards to restarting variables or using initial conditions
+
+It is often convenient to define initial conditions in the `Physics`, and also to be able to
+restart the variables defined by the `Physics` automatically with minimal user effort. User-defined initial conditions
+are convenient to keep the input syntax compact, and default initial conditions are useful to avoid
+non-physical initial states. However, all these objectives conflict when the user defines parameters for initialization in
+a restarted simulation. To make things simple, developers of `Physics` should follow these rules, which we developed based on user
+feedback.
+
+- if the `initialize_variables_from_mesh_file` parameter is set to true, then:
+  - skip adding initial conditions
+  - error if an initial condition parameter is passed by the user to the `Physics`
+- if the `Physics` is set to use (define kernels for) variables that are defined outside the `Physics`, then:
+  - skip adding initial conditions
+  - error if an initial condition parameter is passed by the user to the `Physics`
+- else, if the user specifies initial conditions for variables in the `Physics`
+  - always obey these parameters and add the initial conditions, even if the simulation is restarting
+  - as a sanity check, the [FEProblemBase.md] will error during restarts, unless [!param](/Problem/FEProblem/allow_initial_conditions_with_restart) is set to true
+- else, if the user does not specify initial conditions in the `Physics`, but the `Physics` does define default values for the initial conditions
+  - if the simulation is restarting (from [Checkpoint.md] notably), skip adding the default initial conditions
+  - (redundant due to the first rule) if the `initialize_variables_from_mesh_file` parameter is set to true, skip adding the default initial conditions
+  - (redundant due to the second rule) if the `Physics` is set to use (define kernels for) variables that are defined outside the `Physics`, skip adding the default initial conditions
+
+!alert note
+For `initialize_variables_from_mesh_file` to work correctly, you must use the `saveNonlinearVariable()` and `saveAuxiliaryVariable()` `Physics` routines
+in the constructor of your `Physics` on any variable that you desire to be restarted.

framework/include/problems/FEProblemBase.h

@@ -1881,6 +1881,26 @@ class FEProblemBase : public SubProblem, public Restartable
     _error_on_jacobian_nonzero_reallocation = state;
   }
 
+  /**
+   * Will return True if the executioner in use requires preserving the sparsity pattern of the
+   * matrices being formed during the solve. This is usually the Jacobian.
+   */
+  bool preserveMatrixSparsityPattern() const { return _preserve_matrix_sparsity_pattern; };
+
+  /// Set whether the sparsity pattern of the matrices being formed during the solve (usually the Jacobian)
+  /// should be preserved. This global setting can be retrieved by kernels, notably those using AD, to decide
+  /// whether to take additional care to preserve the sparsity pattern
+  void setPreserveMatrixSparsityPattern(bool preserve);
+
+  /**
+   * Will return true if zeros in the Jacobian are to be dropped from the sparsity pattern.
+   * Note that this can make preserving the matrix sparsity pattern impossible.
+   */
+  bool ignoreZerosInJacobian() const { return _ignore_zeros_in_jacobian; }
+
+  /// Set whether the zeros in the Jacobian should be dropped from the sparsity pattern
+  void setIgnoreZerosInJacobian(bool state) { _ignore_zeros_in_jacobian = state; }
+
   /**
    * Whether or not to accept the solution based on its invalidity.
    *
@@ -1903,10 +1923,6 @@ class FEProblemBase : public SubProblem, public Restartable
    */
   bool immediatelyPrintInvalidSolution() const { return _immediately_print_invalid_solution; }
 
-  bool ignoreZerosInJacobian() const { return _ignore_zeros_in_jacobian; }
-
-  void setIgnoreZerosInJacobian(bool state) { _ignore_zeros_in_jacobian = state; }
-
   /// Returns whether or not this Problem has a TimeIntegrator
   bool hasTimeIntegrator() const { return _has_time_integrator; }
 
@@ -2847,8 +2863,14 @@ class FEProblemBase : public SubProblem, public Restartable
    */
   virtual void resetState();
 
+  // Parameters handling Jacobian sparsity pattern behavior
+  /// Whether to error when the Jacobian is re-allocated, usually because the sparsity pattern changed
   bool _error_on_jacobian_nonzero_reallocation;
+  /// Whether to ignore zeros in the Jacobian, thereby leading to a reduced sparsity pattern
   bool _ignore_zeros_in_jacobian;
+  /// Whether to preserve the system matrix / Jacobian sparsity pattern, using 0-valued entries usually
+  bool _preserve_matrix_sparsity_pattern;
+
   const bool _force_restart;
   const bool _allow_ics_during_restart;
   const bool _skip_nl_system_check;

framework/src/problems/FEProblemBase.C

@@ -464,6 +464,7 @@ FEProblemBase::FEProblemBase(const InputParameters & parameters)
             ? getParam<bool>("error_on_jacobian_nonzero_reallocation")
             : _app.errorOnJacobianNonzeroReallocation()),
     _ignore_zeros_in_jacobian(getParam<bool>("ignore_zeros_in_jacobian")),
+    _preserve_matrix_sparsity_pattern(true),
     _force_restart(getParam<bool>("force_restart")),
     _allow_ics_during_restart(getParam<bool>("allow_initial_conditions_with_restart")),
     _skip_nl_system_check(getParam<bool>("skip_nl_system_check")),
@@ -3594,6 +3595,17 @@ FEProblemBase::getMaterialData(Moose::MaterialDataType type, const THREAD_ID tid
   mooseError("FEProblemBase::getMaterialData(): Invalid MaterialDataType ", type);
 }
 
+void
+FEProblemBase::setPreserveMatrixSparsityPattern(bool preserve)
+{
+  if (_ignore_zeros_in_jacobian && preserve)
+    paramWarning(
+        "ignore_zeros_in_jacobian",
+        "We likely cannot preserve the sparsity pattern if ignoring zeros in the Jacobian, which "
+        "leads to removing those entries from the Jacobian sparsity pattern");
+  _preserve_matrix_sparsity_pattern = preserve;
+}
+
 bool
 FEProblemBase::acceptInvalidSolution() const
 {

modules/heat_transfer/src/physics/HeatConductionPhysicsBase.C

@@ -67,6 +67,9 @@ HeatConductionPhysicsBase::validParams()
 HeatConductionPhysicsBase::HeatConductionPhysicsBase(const InputParameters & parameters)
   : PhysicsBase(parameters), _temperature_name(getParam<VariableName>("temperature_name"))
 {
+  // Save variables (for initialization from file for example)
+  saveSolverVariableName(_temperature_name);
+
   // Parameter checking
   checkVectorParamsSameLength<BoundaryName, MooseFunctorName>("heat_flux_boundaries",
                                                               "boundary_heat_fluxes");
@@ -83,6 +86,16 @@ HeatConductionPhysicsBase::HeatConductionPhysicsBase(const InputParameters & par
 void
 HeatConductionPhysicsBase::addInitialConditions()
 {
+  // error on inconsistent user selections
+  if (getParam<bool>("initialize_variables_from_mesh_file") &&
+      isParamSetByUser("initial_temperature"))
+    paramError("initial_temperature",
+               "Initial temperature should not be set if the variables should be initialized from "
+               "the mesh file");
+  // do not set initial conditions if we load from file
+  if (getParam<bool>("initialize_variables_from_mesh_file"))
+    return;
+
   // Always obey the user, but dont set a hidden default when restarting
   if (!_app.isRestarting() || parameters().isParamSetByUser("initial_temperature"))
   {

modules/heat_transfer/test/tests/physics/restart/gold/default.csv

@@ -0,0 +1,2 @@
+time,max_T,min_T
+0,300,300

modules/heat_transfer/test/tests/physics/restart/gold/from_file.csv

@@ -0,0 +1,2 @@
+time,max_T,min_T
+0,100,100

modules/heat_transfer/test/tests/physics/restart/gold/restart_user_ics.csv

@@ -0,0 +1,2 @@
+time,max_T,min_T
+1,100,100

modules/heat_transfer/test/tests/physics/restart/gold/user_ics.csv

@@ -0,0 +1,2 @@
+time,max_T,min_T
+0,100,100

modules/heat_transfer/test/tests/physics/restart/test_fv.i

@@ -0,0 +1,86 @@
+[Mesh]
+  active = 'cmg'
+  [cmg]
+    type = CartesianMeshGenerator
+    dim = 2
+    dx = 10
+    dy = 10
+  []
+  [fmg_restart]
+    type = FileMeshGenerator
+    file = user_ics.e
+    use_for_exodus_restart = true
+  []
+[]
+
+[Debug]
+  show_actions=true
+[]
+
+[Physics]
+  [HeatConduction]
+    [FiniteVolume]
+      [h1]
+        temperature_name = 'T'
+
+        # Thermal properties
+        thermal_conductivity_functor = 'k0'
+        specific_heat = 5
+        density = 10
+
+        # Boundary conditions
+        heat_flux_boundaries = 'left right'
+        boundary_heat_fluxes = '0 500'
+        insulated_boundaries = 'top'
+        fixed_temperature_boundaries = 'bottom'
+        boundary_temperatures = '300'
+      []
+    []
+  []
+[]
+
+[Executioner]
+  type = Transient
+  num_steps = 1
+  verbose = true
+[]
+
+[Problem]
+  solve = false
+[]
+
+[FunctorMaterials]
+  [mat_k]
+    type = ADGenericFunctorMaterial
+    prop_names = 'k0'
+    prop_values = '1'
+  []
+[]
+
+[Outputs]
+  # Used to set up a restart from checkpoint
+  checkpoint = true
+  # Used to set up a restart from exodus file
+  [exodus]
+    type = Exodus
+    execute_on = TIMESTEP_END
+  []
+  # Used to check results
+  csv = true
+  execute_on = INITIAL
+[]
+
+[Postprocessors]
+  [min_T]
+    type = ElementExtremeValue
+    variable = 'T'
+    value_type = 'min'
+    execute_on = 'INITIAL'
+  []
+  [max_T]
+    type = ElementExtremeValue
+    variable = 'T'
+    value_type = 'max'
+    execute_on = 'INITIAL'
+  []
+[]

modules/heat_transfer/test/tests/physics/restart/tests

@@ -0,0 +1,45 @@
+[Tests]
+  issues = '#28730'
+  # Note: ICs are handled by the parent class HeatConductionPhysics
+  design = 'HeatConductionFV.md HeatConductionCG.md'
+  [restart]
+    requirement = 'The system shall be able to restart the temperature variable in the shorthand Physics-syntax'
+    [default]
+      type = CSVDiff
+      input = test_fv.i
+      csvdiff = 'default.csv'
+      cli_args = 'Outputs/file_base=default'
+      detail = 'using the default initial condition,'
+    []
+    [user_ics]
+      type = CSVDiff
+      input = test_fv.i
+      csvdiff = 'user_ics.csv'
+      cli_args = 'Physics/HeatConduction/FiniteVolume/h1/initial_temperature=100 Outputs/file_base=user_ics'
+      detail = 'with a user-defined initial condition,'
+    []
+    [restart_with_user_ics]
+      type = CSVDiff
+      input = test_fv.i
+      csvdiff = 'restart_user_ics.csv'
+      cli_args = "Physics/HeatConduction/FiniteVolume/h1/initial_temperature=100
+                  Problem/restart_file_base=default_cp/LATEST Problem/allow_initial_conditions_with_restart=true
+                  Outputs/file_base=restart_user_ics"
+      detail = 'when performing a regular checkpoint restart, but still obeying the user-defined initial condition,'
+    []
+    [restart_from_file]
+      type = CSVDiff
+      input = test_fv.i
+      csvdiff = 'from_file.csv'
+      cli_args = "Mesh/active='fmg_restart' Physics/HeatConduction/FiniteVolume/h1/initialize_variables_from_mesh_file=true Outputs/file_base=from_file"
+      detail = 'when performing manual restart from a mesh file, ignoring the default initial condition.'
+    []
+  []
+  [error]
+    type = RunException
+    input = test_fv.i
+    cli_args = 'Physics/HeatConduction/FiniteVolume/h1/initial_temperature=100 Physics/HeatConduction/FiniteVolume/h1/initialize_variables_from_mesh_file=true'
+    expect_err = 'Initial temperature should not be set if the variables should be initialized from the mesh file'
+    requirement = 'The system shall error if the user specifies initial conditions while also requesting variables be loaded from a mesh file.'
+  []
+[]

modules/navier_stokes/include/base/NS.h

@@ -119,8 +119,10 @@ static const std::string k_s = "k_s";
 static const std::string cp = "cp";
 static const std::string cv = "cv";
 static const std::string mu = "mu";
-// Turbulent viscosity
+// Turbulent dynamic viscosity
 static const std::string mu_t = "mu_t";
+// Turbulent dynamic scalar viscosity
+static const std::string mu_t_passive_scalar = "mu_t_passive_scalar";
 // Effective viscosity = sum of viscosities
 static const std::string mu_eff = "mu_eff";
 static const std::string k = "k";
@@ -188,6 +190,8 @@ static constexpr Real von_karman_constant = 0.4187;
 static constexpr Real E_turb_constant = 9.793;
 // Lower limit for mu_t
 static constexpr Real mu_t_low_limit = 1.0e-8;
+// Lower limit for epsilon in the k-epsilon
+static constexpr Real epsilon_low_limit = 1.0e-8;
 // Lower limit for y_plus
 static constexpr Real min_y_plus = 1e-10;
 }

modules/navier_stokes/include/functormaterials/INSFVkEpsilonViscosityFunctorMaterial.h

@@ -30,4 +30,6 @@ class INSFVkEpsilonViscosityFunctorMaterial : public FunctorMaterial
   const Moose::Functor<ADReal> & _epsilon;
   /// The C_mu
   const Moose::Functor<ADReal> & _C_mu;
+  /// Whether to preserve the sparsity pattern between iterations (needed for Newton solvers)
+  const bool _preserve_sparsity_pattern;
 };

modules/navier_stokes/include/fvbcs/INSFVTurbulentViscosityWallFunction.h

@@ -53,7 +53,7 @@ class INSFVTurbulentViscosityWallFunction : public FVDirichletBCBase
   NS::WallTreatmentEnum _wall_treatment;
 
   /// For Newton solves we want to add extra zero-valued terms regardless of y-plus to avoid sparsity pattern changes as y-plus changes near the walls
-  const bool _newton_solve;
+  const bool _preserve_sparsity_pattern;
 
   // Mu_t evaluated at y+=30 for blending purposes
   const Real _mut_30 =

modules/navier_stokes/include/fvkernels/INSFVTurbulentAdvection.h

@@ -37,6 +37,9 @@ class INSFVTurbulentAdvection : public INSFVAdvectionKernel
   /// Wall boundaries
   const std::vector<BoundaryName> & _wall_boundary_names;
 
-  /// Maps for wall treatement
+  /// Maps for wall treatment
   std::map<const Elem *, bool> _wall_bounded;
+
+  /// Whether to remove the derivative of this term wrt to velocity
+  const bool _neglect_advection_derivatives;
 };

modules/navier_stokes/include/fvkernels/INSFVTurbulentDiffusion.h

@@ -41,6 +41,9 @@ class INSFVTurbulentDiffusion : public FVDiffusion
   /// Wall boundaries
   const std::vector<BoundaryName> & _wall_boundary_names;
 
-  /// Maps for wall treatement
+  /// Maps for wall treatment
   std::map<const Elem *, bool> _wall_bounded;
+
+  /// Whether a Newton's method is being used (and we need to preserve the sparsity pattern in edge cases)
+  const bool _preserve_sparsity_pattern;
 };

modules/navier_stokes/src/auxkernels/kEpsilonViscosityAux.C

@@ -196,9 +196,10 @@ kEpsilonViscosityAux::computeValue()
     {
       time_scale = _k(current_argument, state) / _epsilon(current_argument, state);
     }
+    // For newton solvers, epsilon might not be bounded
     if (_newton_solve)
-      if (MooseUtils::absoluteFuzzyEqual(_epsilon(current_argument, state), 0))
-        time_scale = 1;
+      time_scale = _k(current_argument, state) /
+                   std::max(NS::epsilon_low_limit, _epsilon(current_argument, state));
 
     const ADReal mu_t_nl =
         _rho(current_argument, state) * _C_mu * _k(current_argument, state) * time_scale;

modules/navier_stokes/src/base/NSFVBase.C

@@ -328,10 +328,11 @@ NSFVBase::commonTurbulenceParams()
                                     1,
                                     "turbulent_prandtl > 0",
                                     "Turbulent Prandtl number for energy turbulent diffusion");
-  params.addParam<std::vector<Real>>("passive_scalar_schmidt_number",
-                                     std::vector<Real>(),
-                                     "Schmidt numbers used for the passive scalar fields.");
-
+  params.addParam<std::vector<Real>>(
+      "passive_scalar_schmidt_number",
+      std::vector<Real>(),
+      "Turbulent Schmidt numbers used for the passive scalar fields.");
+  params.deprecateParam("passive_scalar_schmidt_number", "Sc_t", "01/01/2025");
   params.addParamNamesToGroup("mixing_length_walls mixing_length_aux_execute_on von_karman_const "
                               "von_karman_const_0 mixing_length_delta",
                               "Mixing length model");

modules/navier_stokes/src/executioners/SegregatedSolverBase.C

@@ -571,6 +571,9 @@ SegregatedSolverBase::SegregatedSolverBase(const InputParameters & parameters)
                                  false);
 
   _time = 0;
+  // This class of segregated solvers for the Navier Stokes equations does not require preserving
+  // entries in a Jacobian
+  _problem.setPreserveMatrixSparsityPattern(false);
 }
 
 void