Merge pull request #29156 from GiudGiud/PR_scalars_sys

Add scalar systems to linearFV version of SIMPLE

framework/include/actioncomponents/ActionComponent.h

@@ -61,7 +61,7 @@ class ActionComponent : public Action
   virtual void setupComponent() {}
 
   // These routines can help define a component that also defines a Physics
-  virtual void addNonlinearVariables() {}
+  virtual void addSolverVariables() {}
 
   /// Used to add one or more Physics to be active on the component.
   /// We recommend using the PhysicsComponentInterface instead of overriding this directly

framework/include/physics/DiffusionCG.h

@@ -23,7 +23,7 @@ class DiffusionCG : public DiffusionPhysicsBase
   DiffusionCG(const InputParameters & parameters);
 
 private:
-  virtual void addNonlinearVariables() override;
+  virtual void addSolverVariables() override;
   virtual void addFEKernels() override;
   virtual void addFEBCs() override;
 

framework/include/physics/DiffusionFV.h

@@ -23,7 +23,7 @@ class DiffusionFV : public DiffusionPhysicsBase
   DiffusionFV(const InputParameters & parameters);
 
 private:
-  virtual void addNonlinearVariables() override;
+  virtual void addSolverVariables() override;
   virtual void addFVKernels() override;
   virtual void addFVBCs() override;
   virtual void initializePhysicsAdditional() override;

framework/include/physics/PhysicsBase.h

@@ -198,7 +198,7 @@ class PhysicsBase : public Action, public InputParametersChecksUtils<PhysicsBase
 
   /// The default implementation of these routines will do nothing as we do not expect all Physics
   /// to be defining an object of every type
-  virtual void addNonlinearVariables() {}
+  virtual void addSolverVariables() {}
   virtual void addAuxiliaryVariables() {}
   virtual void addInitialConditions() {}
   virtual void addFEKernels() {}

framework/src/actioncomponents/ActionComponent.C

@@ -48,7 +48,7 @@ ActionComponent::act()
   // These combined components will likely still exist in the future, when it makes more
   // sense to include the physics than to split it off into its own block
   else if (_current_task == "add_variable")
-    addNonlinearVariables();
+    addSolverVariables();
   else
     // For a new task that isn't registered to ActionComponent in the framework
     actOnAdditionalTasks();

framework/src/physics/DiffusionCG.C

@@ -213,7 +213,7 @@ DiffusionCG::addFEBCs()
 }
 
 void
-DiffusionCG::addNonlinearVariables()
+DiffusionCG::addSolverVariables()
 {
   // If the variable was added outside the Physics
   if (variableExists(_var_name, /*error_if_aux*/ true))

framework/src/physics/DiffusionFV.C

@@ -175,7 +175,7 @@ DiffusionFV::addFVBCs()
 }
 
 void
-DiffusionFV::addNonlinearVariables()
+DiffusionFV::addSolverVariables()
 {
   if (variableExists(_var_name, true))
     return;

framework/src/physics/PhysicsBase.C

@@ -81,7 +81,7 @@ PhysicsBase::act()
   if (_current_task == "init_physics")
     initializePhysics();
   else if (_current_task == "add_variable")
-    addNonlinearVariables();
+    addSolverVariables();
   else if (_current_task == "add_ic")
     addInitialConditions();
 

modules/heat_transfer/include/physics/HeatConductionCG.h

@@ -22,7 +22,7 @@ class HeatConductionCG : public HeatConductionPhysicsBase
   HeatConductionCG(const InputParameters & parameters);
 
 private:
-  void addNonlinearVariables() override;
+  void addSolverVariables() override;
   void addFEKernels() override;
   void addFEBCs() override;
 };

modules/heat_transfer/include/physics/HeatConductionFV.h

@@ -26,6 +26,6 @@ class HeatConductionFV : public HeatConductionPhysicsBase
 
 private:
   virtual void initializePhysicsAdditional() override;
-  virtual void addNonlinearVariables() override;
+  virtual void addSolverVariables() override;
   virtual void addFVKernels() override;
 };

modules/heat_transfer/src/physics/HeatConductionCG.C

@@ -197,7 +197,7 @@ HeatConductionCG::addFEBCs()
 }
 
 void
-HeatConductionCG::addNonlinearVariables()
+HeatConductionCG::addSolverVariables()
 {
   if (variableExists(_temperature_name, /*error_if_aux=*/true))
     return;

modules/heat_transfer/src/physics/HeatConductionFV.C

@@ -209,7 +209,7 @@ HeatConductionFV::addFVBCs()
 }
 
 void
-HeatConductionFV::addNonlinearVariables()
+HeatConductionFV::addSolverVariables()
 {
   if (variableExists(_temperature_name, /*error_if_aux=*/true))
     return;

modules/navier_stokes/doc/content/source/executioners/SIMPLE.md

@@ -146,6 +146,13 @@ As a last step, we add the SIMPLE executioner:
 
 !listing modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/2d-velocity-pressure.i block=Executioner
 
+## Passive scalar advection
+
+The `SIMPLE` executioner can be used to solve coupled problems involving both flow and passive scalar advection.
+Advected passive scalars do not affect the flow distribution, and therefore can be solved after the velocity and
+pressure fields have been computed using the `SIMPLE` algorithm.
+Several systems may be used, for each passive scalar.
+
 !syntax parameters /Executioner/SIMPLE
 
 !syntax inputs /Executioner/SIMPLE

modules/navier_stokes/doc/content/source/linearfvkernels/LinearFVScalarAdvection.md

@@ -0,0 +1,23 @@
+# LinearFVScalarAdvection
+
+This kernel adds the contributions of the scalar advection term to the matrix and right hand side of the scalar equation system for the finite volume SIMPLE segregated solver [SIMPLE.md].
+
+This term is described by $\nabla \cdot \left(\vec{u} C_i \right)$ present in the scalar equation conservation for an incompressible/weakly-compressible formulation.
+
+For FV, the integral of the advection term of scalar $C_i$ over a cell can be expressed as:
+
+\begin{equation}
+\int\limits_{V_C} \nabla \cdot \left(\vec{u} C_i \right) dV \approx \sum\limits_f ( \vec{u}\cdot \vec{n})_{RC} C_if |S_f| \,
+\end{equation}
+
+where $C_{if}$ is the face value of the scalar concentration. An interpolation scheme (e.g. upwind) can be used to compute the face value. This kernel adds the face contribution for each face $f$ to the right hand side and matrix.
+
+The volumetric face flux $(\vec{u}\cdot \vec{n})_{RC}$ is provided by the [RhieChowMassFlux.md] object which uses pressure
+gradients and the discrete momentum equation to compute face velocities and mass fluxes.
+For more information on the expression that is used, see [SIMPLE.md].
+
+!syntax parameters /LinearFVKernels/LinearFVScalarAdvection
+
+!syntax inputs /LinearFVKernels/LinearFVScalarAdvection
+
+!syntax children /LinearFVKernels/LinearFVScalarAdvection

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

@@ -0,0 +1,43 @@
+# Navier Stokes Linear Scalar Transport / WCNSLinearFVScalarTransportPhysics
+
+!syntax description /Physics/NavierStokes/ScalarTransportSegregated/WCNSLinearFVScalarTransportPhysics
+
+## Equation
+
+This [Physics](Physics/index.md) object creates the kernels and boundary conditions to solve the advection-diffusion-reaction
+equation for several scalar quantities advected by the flow.
+
+!equation
+\dfrac{\partial \phi_i}{\partial t} + \nabla \cdot (\phi_i \mathbf{v}) - \nabla \cdot (k_i \nabla \phi_i) - Q_i - \lambda_i \phi_i = 0
+
+where:
+
+- $\phi_i$ is the i-th scalar quantity
+- \mathbf{v} is the advecting velocity
+- $k_i$ the i-th scalar diffusivity
+- $Q_i$ is the i-th scalar source
+- $\lambda_i$ is a reaction coefficient. It should be negative for a loss term
+
+The kernels created are:
+
+- [LinearFVScalarAdvection.md] for the scalar advection term
+- [LinearFVDiffusion.md] for the scalar diffusion term
+- [LinearFVSource.md] for the source terms
+
+Reaction terms can be expressed as source terms by using a negative coefficient in the
+[!param](/Physics/NavierStokes/ScalarTransportSegregated/WCNSLinearFVScalarTransportPhysics/passive_scalar_coupled_source_coeff)
+parameter, and the scalar variable as one of the
+[!param](/Physics/NavierStokes/ScalarTransportSegregated/WCNSLinearFVScalarTransportPhysics/passive_scalar_coupled_source)(s).
+
+## Coupling with other Physics
+
+Scalar advection equations can be solved concurrently with the flow equations by combining the [WCNSLinearFVFlowPhysics.md] with the `WCNSLinearFVScalarTransportPhysics`.
+The following input performs this coupling for incompressible flow in a 2D channel.
+
+!listing test/tests/finite_volume/ins/channel-flow/linear-segregated/2d-scalar/channel-physics.i block=Physics
+
+!syntax parameters /Physics/NavierStokes/ScalarTransportSegregated/WCNSLinearFVScalarTransportPhysics
+
+!syntax inputs /Physics/NavierStokes/ScalarTransportSegregated/WCNSLinearFVScalarTransportPhysics
+
+!syntax children /Physics/NavierStokes/ScalarTransportSegregated/WCNSLinearFVScalarTransportPhysics

modules/navier_stokes/doc/content/syntax/Physics/NavierStokes/ScalarTransportSegregated/index.md

@@ -0,0 +1,11 @@
+# Navier Stokes Scalar Transport Linear Physics
+
+This syntax was created for the [WCNSLinearFVScalarTransportPhysics.md] `Physics`.
+The additional nesting is intended to allow the definition of multiple instances of this `Physics`,
+with different block restrictions.
+
+!syntax list /Physics/NavierStokes/ScalarTransportSegregated objects=True actions=False subsystems=False
+
+!syntax list /Physics/NavierStokes/ScalarTransportSegregated objects=False actions=False subsystems=True
+
+!syntax list /Physics/NavierStokes/ScalarTransportSegregated objects=False actions=True subsystems=False

modules/navier_stokes/include/executioners/SIMPLESolve.h

@@ -68,6 +68,9 @@ class SIMPLESolve : public SIMPLESolveBase
   /// Pointer to the nonlinear system corresponding to the fluid energy equation
   LinearSystem * _energy_system;
 
+  /// Pointer(s) to the system(s) corresponding to the passive scalar equation(s)
+  std::vector<LinearSystem *> _passive_scalar_systems;
+
   /// Pointer to the segregated RhieChow interpolation object
   RhieChowMassFlux * _rc_uo;
 };

modules/navier_stokes/include/executioners/SIMPLESolveBase.h

@@ -136,10 +136,31 @@ class SIMPLESolveBase : public SolveObject, public UserObjectInterface
   /// it needs to be scaled with a representative flux.
   const Real _energy_l_abs_tol;
 
-  // ************************ Iteration control **************************** //
+  // ************************ Passive Scalar Variables ************************ //
 
-  /// The maximum number of momentum-pressure iterations
-  const unsigned int _num_iterations;
+  /// The names of the passive scalar systems
+  const std::vector<SolverSystemName> & _passive_scalar_system_names;
+
+  /// Boolean for easy check if a passive scalar systems shall be solved or not
+  const bool _has_passive_scalar_systems;
+
+  // The number(s) of the system(s) corresponding to the passive scalar equation(s)
+  std::vector<unsigned int> _passive_scalar_system_numbers;
+
+  /// The user-defined relaxation parameter(s) for the passive scalar equation(s)
+  const std::vector<Real> _passive_scalar_equation_relaxation;
+
+  /// Options which hold the petsc settings for the passive scalar equation(s)
+  Moose::PetscSupport::PetscOptions _passive_scalar_petsc_options;
+
+  /// Options for the linear solver of the passive scalar equation(s)
+  SIMPLESolverConfiguration _passive_scalar_linear_control;
+
+  /// Absolute linear tolerance for the passive scalar equation(s). We need to store this, because
+  /// it needs to be scaled with a representative flux.
+  const Real _passive_scalar_l_abs_tol;
+
+  // ************************ Iteration control **************************** //
 
   /// The user-defined absolute tolerance for determining the convergence in momentum
   const Real _momentum_absolute_tolerance;
@@ -150,6 +171,12 @@ class SIMPLESolveBase : public SolveObject, public UserObjectInterface
   /// The user-defined absolute tolerance for determining the convergence in energy
   const Real _energy_absolute_tolerance;
 
+  /// The user-defined absolute tolerance for determining the convergence in passive scalars
+  const std::vector<Real> _passive_scalar_absolute_tolerance;
+
+  /// The maximum number of momentum-pressure iterations
+  const unsigned int _num_iterations;
+
   /// If solve should continue if maximum number of iterations is hit
   const bool _continue_on_max_its;
 

modules/navier_stokes/include/executioners/SIMPLESolveNonlinearAssembly.h

@@ -78,9 +78,6 @@ class SIMPLESolveNonlinearAssembly : public SIMPLESolveBase
   /// Boolean for easy check if a solid energy system shall be solved or not
   const bool _has_solid_energy_system;
 
-  /// Boolean for easy check if a passive scalar systems shall be solved or not
-  const bool _has_passive_scalar_systems;
-
   /// Boolean for easy check if turbulence systems shall be solved or not
   const bool _has_turbulence_systems;
 
@@ -110,30 +107,11 @@ class SIMPLESolveNonlinearAssembly : public SIMPLESolveBase
   /// it needs to be scaled with a representative flux.
   const Real _solid_energy_l_abs_tol;
 
-  // ******************* Passive scalar Eq Variables *********************** //
-
-  /// The names of the passive scalar systems
-  const std::vector<SolverSystemName> & _passive_scalar_system_names;
-
-  // The number(s) of the system(s) corresponding to the passive scalar equation(s)
-  std::vector<unsigned int> _passive_scalar_system_numbers;
+  // ********************* Passive Scalar Eq. Variables *********************** //
 
   /// Pointer(s) to the system(s) corresponding to the passive scalar equation(s)
   std::vector<NonlinearSystemBase *> _passive_scalar_systems;
 
-  /// The user-defined relaxation parameter(s) for the passive scalar equation(s)
-  const std::vector<Real> _passive_scalar_equation_relaxation;
-
-  /// Options which hold the petsc settings for the passive scalar equation(s)
-  Moose::PetscSupport::PetscOptions _passive_scalar_petsc_options;
-
-  /// Options for the linear solver of the passive scalar equation(s)
-  SIMPLESolverConfiguration _passive_scalar_linear_control;
-
-  /// Absolute linear tolerance for the passive scalar equation(s). We need to store this, because
-  /// it needs to be scaled with a representative flux.
-  const Real _passive_scalar_l_abs_tol;
-
   // ********************* Turbulence Eq Variables ************************* //
 
   /// The names of the turbulence scalar systems
@@ -166,9 +144,6 @@ class SIMPLESolveNonlinearAssembly : public SIMPLESolveBase
   /// The user-defined absolute tolerance for determining the convergence in solid energy
   const Real _solid_energy_absolute_tolerance;
 
-  /// The user-defined absolute tolerance for determining the convergence in passive scalars
-  const std::vector<Real> _passive_scalar_absolute_tolerance;
-
   /// The user-defined absolute tolerance for determining the convergence in turbulence equations
   const std::vector<Real> _turbulence_absolute_tolerance;
 

modules/navier_stokes/include/linearfvkernels/LinearFVScalarAdvection.h

@@ -0,0 +1,54 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "LinearFVFluxKernel.h"
+#include "RhieChowMassFlux.h"
+#include "LinearFVAdvectionDiffusionBC.h"
+
+/**
+ * An advection kernel that implements the advection term for the passive scalar transport equation.
+ */
+class LinearFVScalarAdvection : public LinearFVFluxKernel
+{
+public:
+  static InputParameters validParams();
+  LinearFVScalarAdvection(const InputParameters & params);
+
+  virtual Real computeElemMatrixContribution() override;
+
+  virtual Real computeNeighborMatrixContribution() override;
+
+  virtual Real computeElemRightHandSideContribution() override;
+
+  virtual Real computeNeighborRightHandSideContribution() override;
+
+  virtual Real computeBoundaryMatrixContribution(const LinearFVBoundaryCondition & bc) override;
+
+  virtual Real computeBoundaryRHSContribution(const LinearFVBoundaryCondition & bc) override;
+
+  virtual void setupFaceData(const FaceInfo * face_info) override;
+
+protected:
+  /// The Rhie-Chow user object that provides us with the face velocity
+  const RhieChowMassFlux & _mass_flux_provider;
+
+private:
+  /// Container for the current advected interpolation coefficients on the face to make sure
+  /// we don't compute it multiple times for different terms.
+  std::pair<Real, Real> _advected_interp_coeffs;
+
+  /// Container for the velocity on the face which will be reused in the advection term's
+  /// matrix and right hand side contribution
+  Real _volumetric_face_flux;
+
+  /// The interpolation method to use for the advected quantity
+  Moose::FV::InterpMethod _advected_interp_method;
+};

modules/navier_stokes/include/physics/PNSFVSolidHeatTransferPhysics.h

@@ -24,7 +24,7 @@ class PNSFVSolidHeatTransferPhysics final : public HeatConductionFV
 
 protected:
 private:
-  virtual void addNonlinearVariables() override;
+  virtual void addSolverVariables() override;
   virtual void addFVKernels() override;
   virtual void addMaterials() override;