Merge pull request #27614 from GiudGiud/PR_turb_keps

Add physics for fully coupled k-epsilon

framework/include/problems/DumpObjectsProblem.h

@@ -36,6 +36,8 @@ class DumpObjectsProblem : public FEProblemBase
 
   /// output input blocks for a given action path
   void dumpGeneratedSyntax(const std::string path);
+  /// output input blocks for all paths
+  void dumpAllGeneratedSyntax() const;
 
   /// output data in solve
   virtual void solve(unsigned int nl_sys_num) override;

framework/include/problems/FEProblemBase.h

@@ -606,6 +606,9 @@ class FEProblemBase : public SubProblem, public Restartable
    */
   virtual void saveOldSolutions();
 
+  /// Set whether you want the non linear solution to be saved
+  void setSavePreviousNLSolution(bool set) { _previous_nl_solution_required = set; }
+
   /**
    * Restore old solutions from the backup vectors and deallocate them.
    */
@@ -2681,6 +2684,9 @@ class FEProblemBase : public SubProblem, public Restartable
   /// Indicates that we need to compute variable values for previous Newton iteration
   bool _needs_old_newton_iter;
 
+  /// Indicates we need to save the previous NL iteration variable values
+  bool _previous_nl_solution_required;
+
   /// Indicates if nonlocal coupling is required/exists
   bool _has_nonlocal_coupling;
   bool _calculate_jacobian_in_uo;

framework/include/utils/InputParametersChecksUtils.h

@@ -132,6 +132,13 @@ class InputParametersChecksUtils
   void errorDependentParameter(const std::string & param1,
                                const std::string & value_not_set,
                                const std::vector<std::string> & dependent_params) const;
+  /// Error messages for parameters that should depend on another parameter but with a different error message
+  /// @param param1 the parameter has not been set to the desired value (for logging purposes)
+  /// @param value_set the value it has been set to and which is not appropriate (for logging purposes)
+  /// @param dependent_params all the parameters that should not have been since 'param1' was not set to 'value_not_set'
+  void errorInconsistentDependentParameter(const std::string & param1,
+                                           const std::string & value_set,
+                                           const std::vector<std::string> & dependent_params) const;
 
 private:
   // Convenience routines so that defining new checks feels very similar to coding checks in
@@ -554,6 +561,21 @@ InputParametersChecksUtils<C>::errorDependentParameter(
                             "' has not been set to '" + value_not_set + "'");
 }
 
+template <typename C>
+void
+InputParametersChecksUtils<C>::errorInconsistentDependentParameter(
+    const std::string & param1,
+    const std::string & value_set,
+    const std::vector<std::string> & dependent_params) const
+{
+  for (const auto & dependent_param : dependent_params)
+    if (forwardIsParamSetByUser(dependent_param))
+      forwardParamError(dependent_param,
+                        "Parameter '" + dependent_param +
+                            "' should not be set by the user if parameter '" + param1 +
+                            "' has been set to '" + value_set + "'");
+}
+
 template <typename C>
 void
 InputParametersChecksUtils<C>::checkParamsBothSetOrNotSet(const std::string & param1,

framework/include/utils/MooseUtils.h

@@ -1209,6 +1209,21 @@ isDigits(const std::string & str)
 {
   return std::all_of(str.begin(), str.end(), [](unsigned char c) { return std::isdigit(c); });
 }
+
+/**
+ * Courtesy https://stackoverflow.com/a/57163016 and
+ * https://stackoverflow.com/questions/447206/c-isfloat-function
+ * @return Whether the string is convertible to a float
+ */
+inline bool
+isFloat(const std::string & str)
+{
+  if (str.empty())
+    return false;
+  char * ptr;
+  strtof(str.c_str(), &ptr);
+  return (*ptr) == '\0';
+}
 } // MooseUtils namespace
 
 namespace Moose

framework/src/problems/DumpObjectsProblem.C

@@ -23,8 +23,8 @@ DumpObjectsProblem::validParams()
   params.addClassDescription("Single purpose problem object that does not run the given input but "
                              "allows deconstructing actions into their series of underlying Moose "
                              "objects and variables.");
-  params.addRequiredParam<std::string>(
-      "dump_path", "Syntax path of the action of which to dump the generated syntax");
+  params.addParam<std::string>(
+      "dump_path", "all", "Syntax path of the action of which to dump the generated syntax");
   params.addParam<bool>(
       "include_all_user_specified_params",
       true,
@@ -144,7 +144,11 @@ DumpObjectsProblem::dumpVariableHelper(const std::string & system,
 void
 DumpObjectsProblem::solve(unsigned int)
 {
-  dumpGeneratedSyntax(getParam<std::string>("dump_path"));
+  const auto path = getParam<std::string>("dump_path");
+  if (path != "all")
+    dumpGeneratedSyntax(path);
+  else
+    dumpAllGeneratedSyntax();
 }
 
 void
@@ -161,6 +165,17 @@ DumpObjectsProblem::dumpGeneratedSyntax(const std::string path)
   Moose::out << std::flush;
 }
 
+void
+DumpObjectsProblem::dumpAllGeneratedSyntax() const
+{
+  Moose::out << "**START DUMP DATA**\n";
+  for (const auto & path : _generated_syntax)
+    for (const auto & system_pair : path.second)
+      Moose::out << '[' << system_pair.first << "]\n" << system_pair.second << "[]\n\n";
+  Moose::out << "**END DUMP DATA**\n";
+  Moose::out << std::flush;
+}
+
 std::map<std::string, std::string>
 DumpObjectsProblem::stringifyParameters(const InputParameters & parameters)
 {

framework/src/problems/FEProblemBase.C

@@ -431,6 +431,7 @@ FEProblemBase::FEProblemBase(const InputParameters & parameters)
     _const_jacobian(false),
     _has_jacobian(false),
     _needs_old_newton_iter(false),
+    _previous_nl_solution_required(getParam<bool>("previous_nl_solution_required")),
     _has_nonlocal_coupling(false),
     _calculate_jacobian_in_uo(false),
     _kernel_coverage_check(
@@ -627,7 +628,7 @@ FEProblemBase::createTagSolutions()
     _aux->addVector(tag, false, GHOSTED);
   }
 
-  if (getParam<bool>("previous_nl_solution_required"))
+  if (_previous_nl_solution_required)
   {
     // We'll populate the zeroth state of the nonlinear iterations with the current solution for
     // ease of use in doing things like copying solutions backwards. We're just storing pointers in

framework/src/utils/FunctionParserUtils.C

@@ -147,7 +147,9 @@ FunctionParserUtils<is_ad>::addFParserConstants(
   // check constant vectors
   unsigned int nconst = constant_expressions.size();
   if (nconst != constant_names.size())
-    mooseError("The parameter vectors constant_names and constant_values must have equal length.");
+    mooseError("The parameter vectors constant_names (size " +
+               std::to_string(constant_names.size()) + ") and constant_values (size " +
+               std::to_string(nconst) + ") must have equal length.");
 
   // previously evaluated constant_expressions may be used in following constant_expressions
   std::vector<Real> constant_values(nconst);

modules/navier_stokes/doc/content/source/physics/WCNSFVTurbulencePhysics.md

@@ -19,16 +19,47 @@ The following kernels are then added:
 !alert note
 These kernels are only added if each of these equations are being defined using their respective `Physics`.
 
+## K-Epsilon turbulence model
+
+The `WCNSFVTurbulencePhysics` can be set to create the k-epsilon two-equation model.
+
+The turbulent viscosity is then computed with:
+
+- a [kEpsilonViscosityAux.md] if [!param](/Physics/NavierStokes/Turbulence/WCNSFVTurbulencePhysics/mu_t_as_aux_variable) is set to true
+- a [INSFVkEpsilonViscosityFunctorMaterial.md] otherwise
+
+
+The k equation is created with:
+
+- a [FVFunctorTimeKernel.md] for the time derivative if simulating a transient
+- a [INSFVTurbulentAdvection.md] for the turbulent kinetic energy advection term
+- a [INSFVTurbulentDiffusion.md] for the turbulent kinetic energy diffusion term
+- a [INSFVTKESourceSink.md] for the turbulent kinetic energy source and dissipation (sink) terms
+
+
+The epsilon equation is created with:
+
+- a [FVFunctorTimeKernel.md] for the time derivative if simulating a transient
+- a [INSFVTurbulentAdvection.md] for the turbulent kinetic energy dissipation rate advection term
+- a [INSFVTurbulentDiffusion.md] for the turbulent kinetic energy dissipation rate diffusion term
+- a [INSFVTKEDSourceSink.md] for the turbulent kinetic energy dissipation rate source and removal (sink) terms
+
+
+The boundary conditions are not set in this object for the `TKE` and `TKED` variables, as they
+are computed by the wall-functions in the relevant kernels. A boundary condition is set for the turbulent
+viscosity when using an auxiliary variable, with a [INSFVTurbulentViscosityWallFunction.md].
+
+
 ## Coupling with other Physics
 
 A turbulence model can be added to a heat advection solve by using both a `WCNSFVTurbulencePhysics` and a [WCNSFVFluidHeatTransferPhysics.md].
-The following input performs this coupling for weakly compressible flow in a 2D channel.
+The following input performs this coupling for weakly compressible flow for the mixing length turbulence model in a 2D channel.
 No system parameters are passed, so the equations are solved in a fully coupled manner in the same [nonlinear system](systems/NonlinearSystem.md).
 
 !listing test/tests/finite_volume/wcns/channel-flow/2d-transient-physics.i block=Physics
 
 A turbulence model can be added to a scalar advection solve by using both a `WCNSFVTurbulencePhysics` and a [WCNSFVScalarTransportPhysics.md].
-The following input performs this coupling for incompressible flow in a 2D channel.
+The following input performs this coupling for incompressible flow for the mixing length turbulence model in a 2D channel.
 No system parameters are passed, so the equations are solved in a fully coupled manner in the same [nonlinear system](systems/NonlinearSystem.md).
 
 !listing test/tests/finite_volume/ins/channel-flow/2d-mixing-length-physics.i block=Physics

modules/navier_stokes/include/auxkernels/kEpsilonViscosityAux.h

@@ -68,6 +68,9 @@ class kEpsilonViscosityAux : public AuxKernel
   /// Method used to limit the k-e time scale
   const MooseEnum _scale_limiter;
 
+  /// Whether we are using a newton solve
+  const bool _newton_solve;
+
   // -- Parameters of the wall function method
 
   /// Maximum number of iterations to find the friction velocity

modules/navier_stokes/include/base/NS.h

@@ -113,20 +113,25 @@ static const std::string drhos_dTs = "drhos_dTs";
 static const std::string dks_dTs = "dks_dTs";
 static const std::string kappa = "kappa";
 static const std::string kappa_s = "kappa_s";
-static const std::string mu_eff = "mu_eff";
 static const std::string rho_s = "rho_s";
 static const std::string cp_s = "cp_s";
 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
 static const std::string mu_t = "mu_t";
+// Effective viscosity = sum of viscosities
+static const std::string mu_eff = "mu_eff";
 static const std::string k = "k";
+// Turbulence 'conductivity'
+static const std::string k_t = "k_t";
 static const std::string thermal_diffusivity = "thermal_diffusivity";
 static const std::string alpha = "alpha";
 static const std::string alpha_wall = "alpha_wall";
 static const std::string solid = "solid";
 static const std::string Prandtl = "Pr";
+static const std::string turbulent_Prandtl = "Pr_t";
 static const std::string Reynolds = "Re";
 static const std::string Reynolds_hydraulic = "Re_h";
 static const std::string Reynolds_interstitial = "Re_i";
@@ -161,10 +166,10 @@ static const std::string Z = "Z";
 static const std::string K = "K";
 static const std::string mass_flux = "mass_flux";
 
-// Turbuelnce
+// Turbulence
 
 // Turbulence variables
-static const std::string TKE = "k";
+static const std::string TKE = "tke";
 static const std::string TKED = "epsilon";
 
 /**
@@ -183,6 +188,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 y_plus
+static constexpr Real min_y_plus = 1e-10;
 }
 
 namespace NS_DEFAULT_VALUES

modules/navier_stokes/include/fvbcs/INSFVTKEDWallFunctionBC.h

@@ -49,4 +49,7 @@ class INSFVTKEDWallFunctionBC : public FVDirichletBCBase
 
   /// C_mu turbulent coefficient
   const Moose::Functor<ADReal> & _C_mu;
+
+  /// Whether we are using a newton solve
+  const bool _newton_solve;
 };

modules/navier_stokes/include/fvbcs/INSFVTurbulentTemperatureWallFunction.h

@@ -52,4 +52,6 @@ class INSFVTurbulentTemperatureWallFunction : public FVFluxBC
   const Real _C_mu;
   /// Method used for wall treatment
   const MooseEnum _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;
 };

modules/navier_stokes/include/fvbcs/INSFVTurbulentViscosityWallFunction.h

@@ -52,6 +52,9 @@ class INSFVTurbulentViscosityWallFunction : public FVDirichletBCBase
   /// Method used for wall treatment
   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;
+
   // Mu_t evaluated at y+=30 for blending purposes
   const Real _mut_30 =
       (NS::von_karman_constant * 30.0 / std::log(NS::E_turb_constant * 30.0) - 1.0);

modules/navier_stokes/include/fvkernels/INSFVMomentumDiffusion.h

@@ -59,4 +59,7 @@ class INSFVMomentumDiffusion : public INSFVFluxKernel, public SolutionInvalidInt
 
   /// dimension
   const unsigned int _dim;
+
+  /// For Newton solves we want to add extra zero-valued terms to avoid sparsity pattern changes
+  const bool _newton_solve;
 };

modules/navier_stokes/include/fvkernels/INSFVTKEDSourceSink.h

@@ -88,4 +88,7 @@ class INSFVTKEDSourceSink : public FVElementalKernel
   std::map<const Elem *, std::vector<Real>> _dist;
   std::map<const Elem *, std::vector<const FaceInfo *>> _face_infos;
   ///@}
+
+  /// Whether a nonlinear Newton-like solver is being used (as opposed to a linearized solver)
+  const bool _newton_solve;
 };

modules/navier_stokes/include/fvkernels/INSFVTKESourceSink.h

@@ -66,6 +66,9 @@ class INSFVTKESourceSink : public FVElementalKernel
   // Production Limiter Constant
   const Real _C_pl;
 
+  /// 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;
+
   ///@{
   /// Maps for wall treatement
   std::map<const Elem *, bool> _wall_bounded;

modules/navier_stokes/include/physics/WCNSFVCoupledAdvectionPhysicsHelper.h

@@ -13,19 +13,21 @@
 
 class NavierStokesPhysicsBase;
 class WCNSFVFlowPhysics;
+class WCNSFVTurbulencePhysics;
 
 /**
- * Helper class to interact with a WCNSFVFlowPhysics
+ * Helper class to interact with a flow and turbulence physics for a Physics that solves
+ * an advection problem using incompressible or weakly-compressible finite volume
  */
 class WCNSFVCoupledAdvectionPhysicsHelper
 {
 public:
   static InputParameters validParams();
 
-  WCNSFVCoupledAdvectionPhysicsHelper(const InputParameters & parameters,
-                                      const NavierStokesPhysicsBase * derived_physics);
+  WCNSFVCoupledAdvectionPhysicsHelper(const NavierStokesPhysicsBase * derived_physics);
 
-  const WCNSFVFlowPhysics * getCoupledFlowPhysics(const InputParameters & parameters) const;
+  const WCNSFVFlowPhysics * getCoupledFlowPhysics() const;
+  const WCNSFVTurbulencePhysics * getCoupledTurbulencePhysics() const;
 
   /// Return the porosity functor name.
   /// It is important to forward to the Physics so we do not get the smoothing status wrong
@@ -39,6 +41,8 @@ class WCNSFVCoupledAdvectionPhysicsHelper
   const NavierStokesPhysicsBase * _advection_physics;
   /// Flow physics
   const WCNSFVFlowPhysics * _flow_equations_physics;
+  /// Turbulence
+  const WCNSFVTurbulencePhysics * _turbulence_physics;
 
   /// Compressibility type, can be compressible, incompressible or weakly-compressible
   const MooseEnum _compressibility;

modules/navier_stokes/include/physics/WCNSFVFlowPhysics.h

@@ -188,6 +188,8 @@ class WCNSFVFlowPhysics final : public NavierStokesPhysicsBase
   std::map<BoundaryName, MooseEnum> _momentum_outlet_types;
   /// Momentum wall boundary types
   MultiMooseEnum _momentum_wall_types;
+  /// Momentum wall functors
+  std::vector<MooseFunctorName> _momentum_wall_functors;
 
   /// Postprocessors describing the momentum inlet for each boundary. Indexing based on the number of flux boundaries
   std::vector<PostprocessorName> _flux_inlet_pps;

modules/navier_stokes/include/physics/WCNSFVFluidHeatTransferPhysics.h

@@ -30,6 +30,11 @@ class WCNSFVFluidHeatTransferPhysics final : public NavierStokesPhysicsBase,
   /// Get the name of the specific heat material property
   const MooseFunctorName & getSpecificHeatName() const { return _specific_heat_name; }
   MooseFunctorName getSpecificEnthalpyName() const { return NS::specific_enthalpy; }
+  const std::vector<MooseFunctorName> & getThermalConductivityName() const
+  {
+    return _thermal_conductivity_name;
+  }
+
   /// Get the ambient convection parameters for parameter checking
   const std::vector<std::vector<SubdomainName>> & getAmbientConvectionBlocks() const
   {
@@ -46,6 +51,7 @@ class WCNSFVFluidHeatTransferPhysics final : public NavierStokesPhysicsBase,
 
 protected:
 private:
+  void actOnAdditionalTasks() override;
   void addNonlinearVariables() override;
   void addInitialConditions() override;
   void addFVKernels() override;
@@ -65,7 +71,6 @@ class WCNSFVFluidHeatTransferPhysics final : public NavierStokesPhysicsBase,
   void addINSEnergyAdvectionKernels();
   void addINSEnergyAmbientConvection();
   void addINSEnergyExternalHeatSource();
-  void addWCNSEnergyMixingLengthKernels();
 
   /// Functions adding boundary conditions for the incompressible simulation.
   /// These are used for weakly-compressible simulations as well.