diff --git a/framework/include/actioncomponents/ActionComponent.h b/framework/include/actioncomponents/ActionComponent.h index 220192805a32..54681354752b 100644 --- a/framework/include/actioncomponents/ActionComponent.h +++ b/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 diff --git a/framework/include/physics/DiffusionCG.h b/framework/include/physics/DiffusionCG.h index 1a0e9a61f242..fa8021fe5ee2 100644 --- a/framework/include/physics/DiffusionCG.h +++ b/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; diff --git a/framework/include/physics/DiffusionFV.h b/framework/include/physics/DiffusionFV.h index 52eb21920388..8e0dbbff8d50 100644 --- a/framework/include/physics/DiffusionFV.h +++ b/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; diff --git a/framework/include/physics/PhysicsBase.h b/framework/include/physics/PhysicsBase.h index 93213f0bd617..27cd76b9e2c0 100644 --- a/framework/include/physics/PhysicsBase.h +++ b/framework/include/physics/PhysicsBase.h @@ -198,7 +198,7 @@ class PhysicsBase : public Action, public InputParametersChecksUtils _passive_scalar_systems; + /// Pointer to the segregated RhieChow interpolation object RhieChowMassFlux * _rc_uo; }; diff --git a/modules/navier_stokes/include/executioners/SIMPLESolveBase.h b/modules/navier_stokes/include/executioners/SIMPLESolveBase.h index 57437afbe3ef..3a1cc7916353 100644 --- a/modules/navier_stokes/include/executioners/SIMPLESolveBase.h +++ b/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 & _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 _passive_scalar_system_numbers; + + /// The user-defined relaxation parameter(s) for the passive scalar equation(s) + const std::vector _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 _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; diff --git a/modules/navier_stokes/include/executioners/SIMPLESolveNonlinearAssembly.h b/modules/navier_stokes/include/executioners/SIMPLESolveNonlinearAssembly.h index f02f3115f9a5..bb70e89ea3ba 100644 --- a/modules/navier_stokes/include/executioners/SIMPLESolveNonlinearAssembly.h +++ b/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 & _passive_scalar_system_names; - - // The number(s) of the system(s) corresponding to the passive scalar equation(s) - std::vector _passive_scalar_system_numbers; + // ********************* Passive Scalar Eq. Variables *********************** // /// Pointer(s) to the system(s) corresponding to the passive scalar equation(s) std::vector _passive_scalar_systems; - /// The user-defined relaxation parameter(s) for the passive scalar equation(s) - const std::vector _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 _passive_scalar_absolute_tolerance; - /// The user-defined absolute tolerance for determining the convergence in turbulence equations const std::vector _turbulence_absolute_tolerance; diff --git a/modules/navier_stokes/include/linearfvkernels/LinearFVScalarAdvection.h b/modules/navier_stokes/include/linearfvkernels/LinearFVScalarAdvection.h new file mode 100644 index 000000000000..3dffd7e69446 --- /dev/null +++ b/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 _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; +}; diff --git a/modules/navier_stokes/include/physics/PNSFVSolidHeatTransferPhysics.h b/modules/navier_stokes/include/physics/PNSFVSolidHeatTransferPhysics.h index 7d0150ae88d6..3bf4a9faaf89 100644 --- a/modules/navier_stokes/include/physics/PNSFVSolidHeatTransferPhysics.h +++ b/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; diff --git a/modules/navier_stokes/include/physics/WCNSFVCoupledAdvectionPhysicsHelper.h b/modules/navier_stokes/include/physics/WCNSFVCoupledAdvectionPhysicsHelper.h index 07b652a483d3..18936fec8254 100644 --- a/modules/navier_stokes/include/physics/WCNSFVCoupledAdvectionPhysicsHelper.h +++ b/modules/navier_stokes/include/physics/WCNSFVCoupledAdvectionPhysicsHelper.h @@ -12,7 +12,7 @@ #include "PhysicsBase.h" class NavierStokesPhysicsBase; -class WCNSFVFlowPhysics; +class WCNSFVFlowPhysicsBase; class WCNSFVTurbulencePhysics; /** @@ -26,7 +26,7 @@ class WCNSFVCoupledAdvectionPhysicsHelper WCNSFVCoupledAdvectionPhysicsHelper(const NavierStokesPhysicsBase * derived_physics); - const WCNSFVFlowPhysics * getCoupledFlowPhysics() const; + const WCNSFVFlowPhysicsBase * getCoupledFlowPhysics() const; const WCNSFVTurbulencePhysics * getCoupledTurbulencePhysics() const; /// Return the porosity functor name. @@ -40,7 +40,7 @@ class WCNSFVCoupledAdvectionPhysicsHelper /// The Physics class using this helper const NavierStokesPhysicsBase * _advection_physics; /// Flow physics - const WCNSFVFlowPhysics * _flow_equations_physics; + const WCNSFVFlowPhysicsBase * _flow_equations_physics; /// Turbulence const WCNSFVTurbulencePhysics * _turbulence_physics; diff --git a/modules/navier_stokes/include/physics/WCNSFVFlowPhysics.h b/modules/navier_stokes/include/physics/WCNSFVFlowPhysics.h index 4230e424ef5d..f693e95167b1 100644 --- a/modules/navier_stokes/include/physics/WCNSFVFlowPhysics.h +++ b/modules/navier_stokes/include/physics/WCNSFVFlowPhysics.h @@ -23,14 +23,14 @@ class WCNSFVFlowPhysics final : public WCNSFVFlowPhysicsBase WCNSFVFlowPhysics(const InputParameters & parameters); /// Get the name of the linear friction coefficient. Returns an empty string if no friction. - MooseFunctorName getLinearFrictionCoefName() const; + virtual MooseFunctorName getLinearFrictionCoefName() const override; /// Return the name of the Rhie Chow user object UserObjectName rhieChowUOName() const override; /// Return the number of algebraic ghosting layers needed unsigned short getNumberAlgebraicGhostingLayersNeeded() const override; private: - virtual void addNonlinearVariables() override; + virtual void addSolverVariables() override; virtual void addFVKernels() override; virtual void addUserObjects() override; virtual void addCorrectors() override; diff --git a/modules/navier_stokes/include/physics/WCNSFVFlowPhysicsBase.h b/modules/navier_stokes/include/physics/WCNSFVFlowPhysicsBase.h index a396f7bb459e..24a677701c8a 100644 --- a/modules/navier_stokes/include/physics/WCNSFVFlowPhysicsBase.h +++ b/modules/navier_stokes/include/physics/WCNSFVFlowPhysicsBase.h @@ -70,6 +70,8 @@ class WCNSFVFlowPhysicsBase : public NavierStokesPhysicsBase const std::vector & getFluxInletDirections() const { return _flux_inlet_directions; } /// Get the inlet flux postprocessor if using a flux inlet const std::vector & getFluxInletPPs() const { return _flux_inlet_pps; } + /// Get the name of the linear friction coefficient. Returns an empty string if no friction. + virtual MooseFunctorName getLinearFrictionCoefName() const = 0; /// Return the name of the Rhie Chow user object virtual UserObjectName rhieChowUOName() const = 0; /// Return the number of algebraic ghosting layers needed @@ -78,7 +80,7 @@ class WCNSFVFlowPhysicsBase : public NavierStokesPhysicsBase protected: virtual void initializePhysicsAdditional() override; virtual void actOnAdditionalTasks() override; - virtual void addNonlinearVariables() override = 0; + virtual void addSolverVariables() override = 0; virtual void addInitialConditions() override; virtual void addFVKernels() override = 0; virtual void addFVBCs() override; diff --git a/modules/navier_stokes/include/physics/WCNSFVFluidHeatTransferPhysics.h b/modules/navier_stokes/include/physics/WCNSFVFluidHeatTransferPhysics.h index abfe625caaaf..9a452fbd9961 100644 --- a/modules/navier_stokes/include/physics/WCNSFVFluidHeatTransferPhysics.h +++ b/modules/navier_stokes/include/physics/WCNSFVFluidHeatTransferPhysics.h @@ -52,7 +52,7 @@ class WCNSFVFluidHeatTransferPhysics final : public NavierStokesPhysicsBase, protected: private: void actOnAdditionalTasks() override; - void addNonlinearVariables() override; + void addSolverVariables() override; void addInitialConditions() override; void addFVKernels() override; void addFVBCs() override; diff --git a/modules/navier_stokes/include/physics/WCNSFVScalarTransportPhysics.h b/modules/navier_stokes/include/physics/WCNSFVScalarTransportPhysics.h index 74066dc3e0fa..d5b178b77195 100644 --- a/modules/navier_stokes/include/physics/WCNSFVScalarTransportPhysics.h +++ b/modules/navier_stokes/include/physics/WCNSFVScalarTransportPhysics.h @@ -9,8 +9,7 @@ #pragma once -#include "NavierStokesPhysicsBase.h" -#include "WCNSFVCoupledAdvectionPhysicsHelper.h" +#include "WCNSFVScalarTransportPhysicsBase.h" #define registerWCNSFVScalarTransportBaseTasks(app_name, derived_name) \ registerMooseAction(app_name, derived_name, "add_variable"); \ @@ -20,44 +19,17 @@ /** * Creates all the objects needed to solve the Navier Stokes scalar transport equations + * using the nonlinear finite volume weakly-compressible discretization (WCNSFV) */ -class WCNSFVScalarTransportPhysics : public NavierStokesPhysicsBase, - public WCNSFVCoupledAdvectionPhysicsHelper +class WCNSFVScalarTransportPhysics : public WCNSFVScalarTransportPhysicsBase { public: static InputParameters validParams(); WCNSFVScalarTransportPhysics(const InputParameters & parameters); - /// Get the names of the advected scalar quantity variables - const std::vector & getAdvectedScalarNames() const - { - return _passive_scalar_names; - } - - /// Whether the physics is actually creating the scalar advection equations - bool hasScalarEquations() const { return _has_scalar_equation; } - -protected: - virtual void addFVKernels() override; - virtual void addFVBCs() override; - - /// Names of the passive scalar variables - std::vector _passive_scalar_names; - /// A boolean to help compatibility with the old Modules/NavierStokesFV syntax - /// or to deliberately skip adding the equations (for example for mixtures with a stationary phase) - const bool _has_scalar_equation; - - /// Passive scalar inlet boundary types - MultiMooseEnum _passive_scalar_inlet_types; - /// Functors describing the inlet boundary values. See passive_scalar_inlet_types for what the functors actually represent - std::vector> _passive_scalar_inlet_functors; - private: - void addNonlinearVariables() override; - void addInitialConditions() override; - - unsigned short getNumberAlgebraicGhostingLayersNeeded() const override; + virtual void addSolverVariables() override; /** * Functions adding kernels for the incompressible / weakly-compressible scalar transport @@ -65,22 +37,15 @@ class WCNSFVScalarTransportPhysics : public NavierStokesPhysicsBase, * If the material properties are not constant, some of these can be used for * weakly-compressible simulations as well. */ - void addScalarTimeKernels(); - void addScalarDiffusionKernels(); - void addScalarAdvectionKernels(); - virtual void setSlipVelocityParams(InputParameters & /* params */) const {} + virtual void addScalarTimeKernels() override; + virtual void addScalarDiffusionKernels() override; + virtual void addScalarAdvectionKernels() override; /// Equivalent of NSFVAction addScalarCoupledSourceKernels - void addScalarSourceKernels(); + virtual void addScalarSourceKernels() override; /// Functions adding boundary conditions for the incompressible simulation. /// These are used for weakly-compressible simulations as well. - void addScalarInletBC(); - void addScalarWallBC(); - - /// Functors for the passive scalar sources. Indexing is scalar variable index - std::vector _passive_scalar_sources; - /// Functors for the passive scalar (coupled) sources. Inner indexing is scalar variable index - std::vector> _passive_scalar_coupled_sources; - /// Coefficients multiplying for the passive scalar sources. Inner indexing is scalar variable index - std::vector> _passive_scalar_sources_coef; + virtual void addScalarInletBC() override; + virtual void addScalarWallBC() override{}; + virtual void addScalarOutletBC() override; }; diff --git a/modules/navier_stokes/include/physics/WCNSFVScalarTransportPhysicsBase.h b/modules/navier_stokes/include/physics/WCNSFVScalarTransportPhysicsBase.h new file mode 100644 index 000000000000..a08fd4f42693 --- /dev/null +++ b/modules/navier_stokes/include/physics/WCNSFVScalarTransportPhysicsBase.h @@ -0,0 +1,86 @@ +//* 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 "NavierStokesPhysicsBase.h" +#include "WCNSFVCoupledAdvectionPhysicsHelper.h" + +#define registerWCNSFVScalarTransportBaseTasks(app_name, derived_name) \ + registerMooseAction(app_name, derived_name, "add_variable"); \ + registerMooseAction(app_name, derived_name, "add_ic"); \ + registerMooseAction(app_name, derived_name, "add_fv_kernel"); \ + registerMooseAction(app_name, derived_name, "add_fv_bc") + +/** + * Creates all the objects needed to solve the Navier Stokes scalar transport equations + */ +class WCNSFVScalarTransportPhysicsBase : public NavierStokesPhysicsBase, + public WCNSFVCoupledAdvectionPhysicsHelper +{ +public: + static InputParameters validParams(); + + WCNSFVScalarTransportPhysicsBase(const InputParameters & parameters); + + /// Get the names of the advected scalar quantity variables + const std::vector & getAdvectedScalarNames() const + { + return _passive_scalar_names; + } + + /// Whether the physics is actually creating the scalar advection equations + bool hasScalarEquations() const { return _has_scalar_equation; } + +protected: + virtual void addFVKernels() override; + virtual void addFVBCs() override; + virtual void setSlipVelocityParams(InputParameters & /* params */) const {} + + /// Names of the passive scalar variables + std::vector _passive_scalar_names; + /// A boolean to help compatibility with the old Modules/NavierStokesFV syntax + /// or to deliberately skip adding the equations (for example for mixtures with a stationary phase) + const bool _has_scalar_equation; + + /// Passive scalar inlet boundary types + MultiMooseEnum _passive_scalar_inlet_types; + /// Functors describing the inlet boundary values. See passive_scalar_inlet_types for what the functors actually represent + std::vector> _passive_scalar_inlet_functors; + + /// Functors for the passive scalar sources. Indexing is scalar variable index + std::vector _passive_scalar_sources; + /// Functors for the passive scalar (coupled) sources. Inner indexing is scalar variable index + std::vector> _passive_scalar_coupled_sources; + /// Coefficients multiplying for the passive scalar sources. Inner indexing is scalar variable index + std::vector> _passive_scalar_sources_coef; + +private: + virtual void addInitialConditions() override; + + /** + * Functions adding kernels for the incompressible / weakly-compressible scalar transport + * equation + * If the material properties are not constant, some of these can be used for + * weakly-compressible simulations as well. + */ + virtual void addScalarTimeKernels() = 0; + virtual void addScalarDiffusionKernels() = 0; + virtual void addScalarAdvectionKernels() = 0; + /// Equivalent of NSFVAction addScalarCoupledSourceKernels + virtual void addScalarSourceKernels() = 0; + + /// Functions adding boundary conditions for the incompressible simulation. + /// These are used for weakly-compressible simulations as well. + virtual void addScalarInletBC() = 0; + virtual void addScalarWallBC() = 0; + virtual void addScalarOutletBC() = 0; + + virtual unsigned short getNumberAlgebraicGhostingLayersNeeded() const override; +}; diff --git a/modules/navier_stokes/include/physics/WCNSFVTurbulencePhysics.h b/modules/navier_stokes/include/physics/WCNSFVTurbulencePhysics.h index b89bfd09321e..fb5fef49ac4b 100644 --- a/modules/navier_stokes/include/physics/WCNSFVTurbulencePhysics.h +++ b/modules/navier_stokes/include/physics/WCNSFVTurbulencePhysics.h @@ -49,7 +49,7 @@ class WCNSFVTurbulencePhysics final : public NavierStokesPhysicsBase, /// to add the relevant terms (turbulent diffusion notably) void retrieveCoupledPhysics(); - virtual void addNonlinearVariables() override; + virtual void addSolverVariables() override; virtual void addAuxiliaryVariables() override; virtual void addFVKernels() override; virtual void addFVBCs() override; diff --git a/modules/navier_stokes/include/physics/WCNSLinearFVFlowPhysics.h b/modules/navier_stokes/include/physics/WCNSLinearFVFlowPhysics.h index 5c6ac8cea234..d44c79ffedd7 100644 --- a/modules/navier_stokes/include/physics/WCNSLinearFVFlowPhysics.h +++ b/modules/navier_stokes/include/physics/WCNSLinearFVFlowPhysics.h @@ -32,7 +32,7 @@ class WCNSLinearFVFlowPhysics final : public WCNSFVFlowPhysicsBase virtual void initializePhysicsAdditional() override; private: - virtual void addNonlinearVariables() override; + virtual void addSolverVariables() override; virtual void addFVKernels() override; virtual void addUserObjects() override; @@ -57,6 +57,10 @@ class WCNSLinearFVFlowPhysics final : public WCNSFVFlowPhysicsBase virtual void addRhieChowUserObjects() override; + virtual MooseFunctorName getLinearFrictionCoefName() const override + { + mooseError("Not implemented"); + }; UserObjectName rhieChowUOName() const override; unsigned short getNumberAlgebraicGhostingLayersNeeded() const override; diff --git a/modules/navier_stokes/include/physics/WCNSLinearFVScalarTransportPhysics.h b/modules/navier_stokes/include/physics/WCNSLinearFVScalarTransportPhysics.h new file mode 100644 index 000000000000..b02447a91808 --- /dev/null +++ b/modules/navier_stokes/include/physics/WCNSLinearFVScalarTransportPhysics.h @@ -0,0 +1,44 @@ +//* 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 "WCNSFVScalarTransportPhysicsBase.h" + +/** + * Creates all the objects needed to solve the Navier Stokes scalar transport equations + * using the linear finite volume weakly-compressible discretization (WCNSFV) + */ +class WCNSLinearFVScalarTransportPhysics : public WCNSFVScalarTransportPhysicsBase +{ +public: + static InputParameters validParams(); + + WCNSLinearFVScalarTransportPhysics(const InputParameters & parameters); + +private: + virtual void addSolverVariables() override; + + /** + * Functions adding kernels for the incompressible / weakly-compressible scalar transport + * equation + * If the material properties are not constant, some of these can be used for + * weakly-compressible simulations as well. + */ + virtual void addScalarTimeKernels() override; + virtual void addScalarDiffusionKernels() override; + virtual void addScalarAdvectionKernels() override; + virtual void addScalarSourceKernels() override; + + /// Functions adding boundary conditions for the incompressible simulation. + /// These are used for weakly-compressible simulations as well. + virtual void addScalarInletBC() override; + virtual void addScalarWallBC() override{}; + virtual void addScalarOutletBC() override; +}; diff --git a/modules/navier_stokes/include/userobjects/RhieChowMassFlux.h b/modules/navier_stokes/include/userobjects/RhieChowMassFlux.h index 973e93db2e01..b88d089a3a02 100644 --- a/modules/navier_stokes/include/userobjects/RhieChowMassFlux.h +++ b/modules/navier_stokes/include/userobjects/RhieChowMassFlux.h @@ -44,8 +44,10 @@ class RhieChowMassFlux : public GeneralUserObject, static InputParameters validParams(); RhieChowMassFlux(const InputParameters & params); - /// Get the face velocity (used in advection terms) + /// Get the face velocity times density (used in advection terms) Real getMassFlux(const FaceInfo & fi) const; + /// Get the volumetric face flux (used in advection terms) + Real getVolumetricFaceFlux(const FaceInfo & fi) const; /// Initialize the container for face velocities void initFaceMassFlux(); diff --git a/modules/navier_stokes/src/base/NavierStokesApp.C b/modules/navier_stokes/src/base/NavierStokesApp.C index f7388ba05b5a..1b6c7ad20d00 100644 --- a/modules/navier_stokes/src/base/NavierStokesApp.C +++ b/modules/navier_stokes/src/base/NavierStokesApp.C @@ -54,6 +54,8 @@ associateSyntaxInner(Syntax & syntax, ActionFactory & /*action_factory*/) registerSyntax("WCNSLinearFVFlowPhysics", "Physics/NavierStokes/FlowSegregated/*"); registerSyntax("WCNSFVFluidHeatTransferPhysics", "Physics/NavierStokes/FluidHeatTransfer/*"); registerSyntax("WCNSFVScalarTransportPhysics", "Physics/NavierStokes/ScalarTransport/*"); + registerSyntax("WCNSLinearFVScalarTransportPhysics", + "Physics/NavierStokes/ScalarTransportSegregated/*"); registerSyntax("WCNSFVTurbulencePhysics", "Physics/NavierStokes/Turbulence/*"); registerSyntax("PNSFVSolidHeatTransferPhysics", "Physics/NavierStokes/SolidHeatTransfer/*"); registerSyntax("WCNSFVTwoPhaseMixturePhysics", "Physics/NavierStokes/TwoPhaseMixture/*"); diff --git a/modules/navier_stokes/src/executioners/SIMPLESolve.C b/modules/navier_stokes/src/executioners/SIMPLESolve.C index 9508144bb306..75e87778d559 100644 --- a/modules/navier_stokes/src/executioners/SIMPLESolve.C +++ b/modules/navier_stokes/src/executioners/SIMPLESolve.C @@ -34,6 +34,15 @@ SIMPLESolve::SIMPLESolve(Executioner & ex) _momentum_system_numbers.push_back(_problem.linearSysNum(_momentum_system_names[system_i])); _momentum_systems.push_back(&_problem.getLinearSystem(_momentum_system_numbers[system_i])); } + // and for the passive scalar equations + if (_has_passive_scalar_systems) + for (auto system_i : index_range(_passive_scalar_system_names)) + { + _passive_scalar_system_numbers.push_back( + _problem.linearSysNum(_passive_scalar_system_names[system_i])); + _passive_scalar_systems.push_back( + &_problem.getLinearSystem(_passive_scalar_system_numbers[system_i])); + } } void @@ -365,6 +374,48 @@ SIMPLESolve::solve() converged = NS::FV::converged(ns_residuals, ns_abs_tols); } + // If we have passive scalar equations, solve them here. We assume the material properties in the + // Navier-Stokes equations do not depend on passive scalars, as they are passive, therefore we + // solve outside of the velocity-pressure loop + if (_has_passive_scalar_systems && converged) + { + // The reason why we need more than one iteration is due to the matrix relaxation + // which can be used to stabilize the equations + bool passive_scalar_converged = false; + iteration_counter = 0; + while (iteration_counter < _num_iterations && !passive_scalar_converged) + { + std::vector> scalar_residuals( + _passive_scalar_system_names.size(), std::make_pair(0, 1.0)); + std::vector scalar_abs_tols; + for (const auto scalar_tol : _passive_scalar_absolute_tolerance) + scalar_abs_tols.push_back(scalar_tol); + + // We set the preconditioner/controllable parameters through petsc options. Linear + // tolerances will be overridden within the solver. + Moose::PetscSupport::petscSetOptions(_passive_scalar_petsc_options, solver_params); + for (const auto i : index_range(_passive_scalar_system_names)) + scalar_residuals[i] = solveAdvectedSystem(_passive_scalar_system_numbers[i], + *_passive_scalar_systems[i], + _passive_scalar_equation_relaxation[i], + _passive_scalar_linear_control, + _passive_scalar_l_abs_tol); + + if (_has_passive_scalar_systems) + _console << " Passive scalar equations:\n"; + for (const auto i : index_range(_passive_scalar_system_names)) + _console << COLOR_GREEN << " " << _passive_scalar_system_names[i] << ": " + << scalar_residuals[i].second << COLOR_DEFAULT + << " Linear its: " << scalar_residuals[i].first << std::endl; + + passive_scalar_converged = NS::FV::converged(scalar_residuals, scalar_abs_tols); + iteration_counter++; + } + + // Both flow and scalars must converge + converged = passive_scalar_converged && converged; + } + converged = _continue_on_max_its ? true : converged; return converged; diff --git a/modules/navier_stokes/src/executioners/SIMPLESolveBase.C b/modules/navier_stokes/src/executioners/SIMPLESolveBase.C index 393fd6f39555..9006ef8a7a1c 100644 --- a/modules/navier_stokes/src/executioners/SIMPLESolveBase.C +++ b/modules/navier_stokes/src/executioners/SIMPLESolveBase.C @@ -25,6 +25,8 @@ SIMPLESolveBase::validParams() params.addRequiredParam("pressure_system", "The solver system for the pressure equation."); params.addParam("energy_system", "The solver system for the energy equation."); + params.addParam>( + "passive_scalar_systems", {}, "The solver system for each scalar advection equation."); /* * Parameters to control the solution of the momentum equation @@ -183,6 +185,59 @@ SIMPLESolveBase::validParams() "0>("passive_scalar_equation_relaxation", + std::vector(), + "The relaxation which should be used for the passive scalar " + "equations. (=1 for no relaxation, " + "diagonal dominance will still be enforced)"); + + params.addParam("passive_scalar_petsc_options", + Moose::PetscSupport::getCommonPetscFlags(), + "Singleton PETSc options for the passive scalar equation(s)"); + params.addParam( + "passive_scalar_petsc_options_iname", + Moose::PetscSupport::getCommonPetscKeys(), + "Names of PETSc name/value pairs for the passive scalar equation(s)"); + params.addParam>( + "passive_scalar_petsc_options_value", + "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\" for the " + "passive scalar equation(s)"); + params.addParam>( + "passive_scalar_absolute_tolerance", + std::vector(), + "The absolute tolerance(s) on the normalized residual(s) of the passive scalar equation(s)."); + params.addRangeCheckedParam("passive_scalar_l_tol", + 1e-5, + "0.0<=passive_scalar_l_tol & passive_scalar_l_tol<1.0", + "The relative tolerance on the normalized residual in the " + "linear solver of the passive scalar equation(s)."); + params.addRangeCheckedParam("passive_scalar_l_abs_tol", + 1e-10, + "0.0( + "passive_scalar_l_max_its", + 10000, + "The maximum allowed iterations in the linear solver of the turbulence equation."); + + params.addParamNamesToGroup( + "passive_scalar_systems passive_scalar_equation_relaxation passive_scalar_petsc_options " + "passive_scalar_petsc_options_iname " + "passive_scalar_petsc_options_value passive_scalar_petsc_options_value " + "passive_scalar_absolute_tolerance " + "passive_scalar_l_tol passive_scalar_l_abs_tol passive_scalar_l_max_its", + "passive_scalar Equation"); + /* * SIMPLE iteration control */ @@ -214,10 +269,17 @@ SIMPLESolveBase::SIMPLESolveBase(Executioner & ex) _has_energy_system(isParamValid("energy_system")), _energy_equation_relaxation(getParam("energy_equation_relaxation")), _energy_l_abs_tol(getParam("energy_l_abs_tol")), - _num_iterations(getParam("num_iterations")), + _passive_scalar_system_names(getParam>("passive_scalar_systems")), + _has_passive_scalar_systems(!_passive_scalar_system_names.empty()), + _passive_scalar_equation_relaxation( + getParam>("passive_scalar_equation_relaxation")), + _passive_scalar_l_abs_tol(getParam("passive_scalar_l_abs_tol")), _momentum_absolute_tolerance(getParam("momentum_absolute_tolerance")), _pressure_absolute_tolerance(getParam("pressure_absolute_tolerance")), _energy_absolute_tolerance(getParam("energy_absolute_tolerance")), + _passive_scalar_absolute_tolerance( + getParam>("passive_scalar_absolute_tolerance")), + _num_iterations(getParam("num_iterations")), _continue_on_max_its(getParam("continue_on_max_its")), _print_fields(getParam("print_fields")) { @@ -274,6 +336,50 @@ SIMPLESolveBase::SIMPLESolveBase(Executioner & ex) "energy_absolute_tolerance", "energy_equation_relaxation"}, false); + + // We check for input errors with regards to the passive scalar equations. At the same time, we + // set up the corresponding system numbers + if (_has_passive_scalar_systems) + { + if (_passive_scalar_system_names.size() != _passive_scalar_equation_relaxation.size()) + paramError("passive_scalar_equation_relaxation", + "The number of equation relaxation parameters does not match the number of " + "passive scalar equations!"); + if (_passive_scalar_system_names.size() != _passive_scalar_absolute_tolerance.size()) + paramError("passive_scalar_absolute_tolerance", + "The number of absolute tolerances does not match the number of " + "passive scalar equations!"); + } + if (_has_passive_scalar_systems) + { + const auto & passive_scalar_petsc_options = + getParam("passive_scalar_petsc_options"); + const auto & passive_scalar_petsc_pair_options = getParam( + "passive_scalar_petsc_options_iname", "passive_scalar_petsc_options_value"); + Moose::PetscSupport::processPetscFlags(passive_scalar_petsc_options, + _passive_scalar_petsc_options); + Moose::PetscSupport::processPetscPairs(passive_scalar_petsc_pair_options, + _problem.mesh().dimension(), + _passive_scalar_petsc_options); + + _passive_scalar_linear_control.real_valued_data["rel_tol"] = + getParam("passive_scalar_l_tol"); + _passive_scalar_linear_control.real_valued_data["abs_tol"] = + getParam("passive_scalar_l_abs_tol"); + _passive_scalar_linear_control.int_valued_data["max_its"] = + getParam("passive_scalar_l_max_its"); + } + else + checkDependentParameterError("passive_scalar_systems", + {"passive_scalar_petsc_options", + "passive_scalar_petsc_options_iname", + "passive_scalar_petsc_options_value", + "passive_scalar_l_tol", + "passive_scalar_l_abs_tol", + "passive_scalar_l_max_its", + "passive_scalar_equation_relaxation", + "passive_scalar_absolute_tolerance"}, + false); } void diff --git a/modules/navier_stokes/src/executioners/SIMPLESolveNonlinearAssembly.C b/modules/navier_stokes/src/executioners/SIMPLESolveNonlinearAssembly.C index c7cb52af6b71..e049196173c0 100644 --- a/modules/navier_stokes/src/executioners/SIMPLESolveNonlinearAssembly.C +++ b/modules/navier_stokes/src/executioners/SIMPLESolveNonlinearAssembly.C @@ -24,23 +24,14 @@ SIMPLESolveNonlinearAssembly::validParams() /* * The names of the different systems in the segregated solver */ - params.addParam("energy_system", "The solver system for the energy equation."); params.addParam("solid_energy_system", "The solver system for the solid energy equation."); - params.addParam>( - "passive_scalar_systems", {}, "The solver system(s) for the passive scalar equation(s)."); params.addParam>( "turbulence_systems", {}, "The solver system(s) for the turbulence equation(s)."); /* * Relaxation parameters for the different system */ - - params.addParam>("passive_scalar_equation_relaxation", - std::vector(), - "The relaxation which should be used for the passive scalar " - "equations. (=1 for no relaxation, " - "diagonal dominance will still be enforced)"); params.addParam>( "turbulence_equation_relaxation", std::vector(), @@ -54,7 +45,7 @@ SIMPLESolveNonlinearAssembly::validParams() "The lower limit imposed on turbulent quantities. The recommended value for robustness " "is 1e-8."); - params.addParamNamesToGroup("energy_equation_relaxation passive_scalar_equation_relaxation " + params.addParamNamesToGroup("energy_equation_relaxation " "turbulence_equation_relaxation", "Relaxation"); @@ -72,18 +63,6 @@ SIMPLESolveNonlinearAssembly::validParams() "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\" for the " "solid energy equation"); - params.addParam("passive_scalar_petsc_options", - Moose::PetscSupport::getCommonPetscFlags(), - "Singleton PETSc options for the passive scalar equation(s)"); - params.addParam( - "passive_scalar_petsc_options_iname", - Moose::PetscSupport::getCommonPetscKeys(), - "Names of PETSc name/value pairs for the passive scalar equation(s)"); - params.addParam>( - "passive_scalar_petsc_options_value", - "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\" for the " - "passive scalar equation(s)"); - params.addParam("turbulence_petsc_options", Moose::PetscSupport::getCommonPetscFlags(), "Singleton PETSc options for the turbulence equation(s)"); @@ -97,8 +76,7 @@ SIMPLESolveNonlinearAssembly::validParams() params.addParamNamesToGroup( "solid_energy_petsc_options solid_energy_petsc_options_iname " - "solid_energy_petsc_options_value passive_scalar_petsc_options " - "passive_scalar_petsc_options_iname passive_scalar_petsc_options_value " + "solid_energy_petsc_options_value " "turbulence_petsc_options turbulence_petsc_options_iname turbulence_petsc_options_value", "PETSc Control"); @@ -110,17 +88,12 @@ SIMPLESolveNonlinearAssembly::validParams() 1e-5, "0.0>( - "passive_scalar_absolute_tolerance", - std::vector(), - "The absolute tolerance(s) on the normalized residual(s) of the passive scalar equation(s)."); params.addParam>( "turbulence_absolute_tolerance", std::vector(), "The absolute tolerance(s) on the normalized residual(s) of the turbulence equation(s)."); - params.addParamNamesToGroup("solid_energy_absolute_tolerance passive_scalar_absolute_tolerance " - "turbulence_absolute_tolerance", + params.addParamNamesToGroup("solid_energy_absolute_tolerance turbulence_absolute_tolerance", "Iteration Control"); /* * Linear iteration tolerances for the different equations @@ -140,20 +113,6 @@ SIMPLESolveNonlinearAssembly::validParams() 10000, "0("passive_scalar_l_tol", - 1e-5, - "0.0<=passive_scalar_l_tol & passive_scalar_l_tol<1.0", - "The relative tolerance on the normalized residual in the " - "linear solver of the passive scalar equation(s)."); - params.addRangeCheckedParam("passive_scalar_l_abs_tol", - 1e-10, - "0.0( - "passive_scalar_l_max_its", - 10000, - "The maximum allowed iterations in the linear solver of the turbulence equation."); params.addRangeCheckedParam("turbulence_l_tol", 1e-5, "0.0<=turbulence_l_tol & turbulence_l_tol<1.0", @@ -169,10 +128,10 @@ SIMPLESolveNonlinearAssembly::validParams() 10000, "The maximum allowed iterations in the linear solver of the turbulence equation(s)."); - params.addParamNamesToGroup("solid_energy_l_tol solid_energy_l_abs_tol solid_energy_l_max_its " - "passive_scalar_l_abs_tol passive_scalar_l_max_its turbulence_l_tol " - "turbulence_l_abs_tol turbulence_l_max_its", - "Linear Iteration Control"); + params.addParamNamesToGroup( + "solid_energy_l_tol solid_energy_l_abs_tol solid_energy_l_max_its turbulence_l_tol " + "turbulence_l_abs_tol turbulence_l_max_its", + "Linear Iteration Control"); return params; } @@ -182,8 +141,6 @@ SIMPLESolveNonlinearAssembly::SIMPLESolveNonlinearAssembly(Executioner & ex) _pressure_sys_number(_problem.nlSysNum(getParam("pressure_system"))), _pressure_system(_problem.getNonlinearSystemBase(_pressure_sys_number)), _has_solid_energy_system(_has_energy_system && isParamValid("solid_energy_system")), - _has_passive_scalar_systems( - !getParam>("passive_scalar_systems").empty()), _has_turbulence_systems(!getParam>("turbulence_systems").empty()), _energy_sys_number(_has_energy_system ? _problem.nlSysNum(getParam("energy_system")) @@ -198,17 +155,11 @@ SIMPLESolveNonlinearAssembly::SIMPLESolveNonlinearAssembly(Executioner & ex) ? &_problem.getNonlinearSystemBase(_solid_energy_sys_number) : nullptr), _solid_energy_l_abs_tol(getParam("solid_energy_l_abs_tol")), - _passive_scalar_system_names(getParam>("passive_scalar_systems")), - _passive_scalar_equation_relaxation( - getParam>("passive_scalar_equation_relaxation")), - _passive_scalar_l_abs_tol(getParam("passive_scalar_l_abs_tol")), _turbulence_system_names(getParam>("turbulence_systems")), _turbulence_equation_relaxation(getParam>("turbulence_equation_relaxation")), _turbulence_field_min_limit(getParam>("turbulence_field_min_limit")), _turbulence_l_abs_tol(getParam("turbulence_l_abs_tol")), _solid_energy_absolute_tolerance(getParam("solid_energy_absolute_tolerance")), - _passive_scalar_absolute_tolerance( - getParam>("passive_scalar_absolute_tolerance")), _turbulence_absolute_tolerance(getParam>("turbulence_absolute_tolerance")), _pressure_tag_name(getParam("pressure_gradient_tag")), _pressure_tag_id(_problem.addVectorTag(_pressure_tag_name)) @@ -240,20 +191,6 @@ SIMPLESolveNonlinearAssembly::SIMPLESolveNonlinearAssembly(Executioner & ex) &_problem.getNonlinearSystemBase(_turbulence_system_numbers[system_i])); } - // We check for input errors with regards to the passive scalar equations. At the same time, we - // set up the corresponding system numbers - if (_has_passive_scalar_systems) - { - if (_passive_scalar_system_names.size() != _passive_scalar_equation_relaxation.size()) - paramError("passive_scalar_equation_relaxation", - "The number of equation relaxation parameters does not match the number of " - "passive scalar equations!"); - if (_passive_scalar_system_names.size() != _passive_scalar_absolute_tolerance.size()) - paramError("passive_scalar_absolute_tolerance", - "The number of absolute tolerances does not match the number of " - "passive scalar equations!"); - } - // We check for input errors with regards to the turbulence equations. At the same time, we // set up the corresponding system numbers if (_has_turbulence_systems) @@ -314,37 +251,6 @@ SIMPLESolveNonlinearAssembly::SIMPLESolveNonlinearAssembly(Executioner & ex) false); } - if (_has_passive_scalar_systems) - { - const auto & passive_scalar_petsc_options = - getParam("passive_scalar_petsc_options"); - const auto & passive_scalar_petsc_pair_options = getParam( - "passive_scalar_petsc_options_iname", "passive_scalar_petsc_options_value"); - Moose::PetscSupport::processPetscFlags(passive_scalar_petsc_options, - _passive_scalar_petsc_options); - Moose::PetscSupport::processPetscPairs(passive_scalar_petsc_pair_options, - _problem.mesh().dimension(), - _passive_scalar_petsc_options); - - _passive_scalar_linear_control.real_valued_data["rel_tol"] = - getParam("passive_scalar_l_tol"); - _passive_scalar_linear_control.real_valued_data["abs_tol"] = - getParam("passive_scalar_l_abs_tol"); - _passive_scalar_linear_control.int_valued_data["max_its"] = - getParam("passive_scalar_l_max_its"); - } - else - checkDependentParameterError("passive_scalar_system", - {"passive_scalar_petsc_options", - "passive_scalar_petsc_options_iname", - "passive_scalar_petsc_options_value", - "passive_scalar_l_tol", - "passive_scalar_l_abs_tol", - "passive_scalar_l_max_its", - "passive_scalar_equation_relaxation", - "passive_scalar_absolute_tolerance"}, - false); - if (_has_turbulence_systems) { const auto & turbulence_petsc_options = getParam("turbulence_petsc_options"); diff --git a/modules/navier_stokes/src/linearfvkernels/LinearFVScalarAdvection.C b/modules/navier_stokes/src/linearfvkernels/LinearFVScalarAdvection.C new file mode 100644 index 000000000000..175c5057934b --- /dev/null +++ b/modules/navier_stokes/src/linearfvkernels/LinearFVScalarAdvection.C @@ -0,0 +1,103 @@ +//* 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 + +#include "LinearFVScalarAdvection.h" +#include "MooseLinearVariableFV.h" +#include "NSFVUtils.h" +#include "NS.h" + +registerMooseObject("NavierStokesApp", LinearFVScalarAdvection); + +InputParameters +LinearFVScalarAdvection::validParams() +{ + InputParameters params = LinearFVFluxKernel::validParams(); + params.addClassDescription("Represents the matrix and right hand side contributions of an " + "advection term for a passive scalar."); + params.addRequiredParam( + "rhie_chow_user_object", + "The rhie-chow user-object which is used to determine the face velocity."); + params += Moose::FV::advectedInterpolationParameter(); + return params; +} + +LinearFVScalarAdvection::LinearFVScalarAdvection(const InputParameters & params) + : LinearFVFluxKernel(params), + _mass_flux_provider(getUserObject("rhie_chow_user_object")), + _advected_interp_coeffs(std::make_pair(0, 0)), + _volumetric_face_flux(0.0) +{ + Moose::FV::setInterpolationMethod(*this, _advected_interp_method, "advected_interp_method"); +} + +Real +LinearFVScalarAdvection::computeElemMatrixContribution() +{ + return _advected_interp_coeffs.first * _volumetric_face_flux * _current_face_area; +} + +Real +LinearFVScalarAdvection::computeNeighborMatrixContribution() +{ + return _advected_interp_coeffs.second * _volumetric_face_flux * _current_face_area; +} + +Real +LinearFVScalarAdvection::computeElemRightHandSideContribution() +{ + return 0.0; +} + +Real +LinearFVScalarAdvection::computeNeighborRightHandSideContribution() +{ + return 0.0; +} + +Real +LinearFVScalarAdvection::computeBoundaryMatrixContribution(const LinearFVBoundaryCondition & bc) +{ + const auto * const adv_bc = static_cast(&bc); + mooseAssert(adv_bc, "This should be a valid BC!"); + + const auto boundary_value_matrix_contrib = adv_bc->computeBoundaryValueMatrixContribution(); + + // We support internal boundaries too so we have to make sure the normal points always outward + const auto factor = (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM) ? 1.0 : -1.0; + + return boundary_value_matrix_contrib * factor * _volumetric_face_flux * _current_face_area; +} + +Real +LinearFVScalarAdvection::computeBoundaryRHSContribution(const LinearFVBoundaryCondition & bc) +{ + const auto * const adv_bc = static_cast(&bc); + mooseAssert(adv_bc, "This should be a valid BC!"); + + // We support internal boundaries too so we have to make sure the normal points always outward + const auto factor = (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM ? 1.0 : -1.0); + + const auto boundary_value_rhs_contrib = adv_bc->computeBoundaryValueRHSContribution(); + return -boundary_value_rhs_contrib * factor * _volumetric_face_flux * _current_face_area; +} + +void +LinearFVScalarAdvection::setupFaceData(const FaceInfo * face_info) +{ + LinearFVFluxKernel::setupFaceData(face_info); + + // Caching the velocity on the face which will be reused in the advection term's matrix and right + // hand side contributions + _volumetric_face_flux = _mass_flux_provider.getVolumetricFaceFlux(*face_info); + + // Caching the interpolation coefficients so they will be reused for the matrix and right hand + // side terms + _advected_interp_coeffs = + interpCoeffs(_advected_interp_method, *_current_face_info, true, _volumetric_face_flux); +} diff --git a/modules/navier_stokes/src/physics/PNSFVSolidHeatTransferPhysics.C b/modules/navier_stokes/src/physics/PNSFVSolidHeatTransferPhysics.C index b7b41e20118c..e8ab4a785299 100644 --- a/modules/navier_stokes/src/physics/PNSFVSolidHeatTransferPhysics.C +++ b/modules/navier_stokes/src/physics/PNSFVSolidHeatTransferPhysics.C @@ -154,7 +154,7 @@ PNSFVSolidHeatTransferPhysics::PNSFVSolidHeatTransferPhysics(const InputParamete } void -PNSFVSolidHeatTransferPhysics::addNonlinearVariables() +PNSFVSolidHeatTransferPhysics::addSolverVariables() { // Dont add if the user already defined the variable if (variableExists(_solid_temperature_name, diff --git a/modules/navier_stokes/src/physics/WCNSFVCoupledAdvectionPhysicsHelper.C b/modules/navier_stokes/src/physics/WCNSFVCoupledAdvectionPhysicsHelper.C index 6f7ccca4cd62..de095781f1fa 100644 --- a/modules/navier_stokes/src/physics/WCNSFVCoupledAdvectionPhysicsHelper.C +++ b/modules/navier_stokes/src/physics/WCNSFVCoupledAdvectionPhysicsHelper.C @@ -10,7 +10,7 @@ #include "WCNSFVCoupledAdvectionPhysicsHelper.h" #include "INSFVRhieChowInterpolator.h" #include "RelationshipManager.h" -#include "WCNSFVFlowPhysics.h" +#include "WCNSFVFlowPhysicsBase.h" #include "WCNSFVTurbulencePhysics.h" InputParameters @@ -49,17 +49,18 @@ WCNSFVCoupledAdvectionPhysicsHelper::getPorosityFunctorName(bool smoothed) const return _flow_equations_physics->getPorosityFunctorName(smoothed); } -const WCNSFVFlowPhysics * +const WCNSFVFlowPhysicsBase * WCNSFVCoupledAdvectionPhysicsHelper::getCoupledFlowPhysics() const { // User passed it, just use that if (_advection_physics->isParamValid("coupled_flow_physics")) - return _advection_physics->getCoupledPhysics( + return _advection_physics->getCoupledPhysics( _advection_physics->getParam("coupled_flow_physics")); // Look for any physics of the right type, and check the block restriction else { - const auto all_flow_physics = _advection_physics->getCoupledPhysics(); + const auto all_flow_physics = + _advection_physics->getCoupledPhysics(); for (const auto physics : all_flow_physics) if (_advection_physics->checkBlockRestrictionIdentical( physics->name(), physics->blocks(), /*error_if_not_identical=*/false)) @@ -67,9 +68,9 @@ WCNSFVCoupledAdvectionPhysicsHelper::getCoupledFlowPhysics() const return physics; } } - mooseError( - "No coupled flow Physics found of type 'WCNSFVFlowPhysics'. Use the 'coupled_flow_physics' " - "parameter to give the name of the desired WCNSFVFlowPhysics to couple with"); + mooseError("No coupled flow Physics found of type derived from 'WCNSFVFlowPhysicsBase'. Use the " + "'coupled_flow_physics' parameter to give the name of the desired " + "WCNSFVFlowPhysicsBase-derived Physics to couple with"); } const WCNSFVTurbulencePhysics * diff --git a/modules/navier_stokes/src/physics/WCNSFVFlowPhysics.C b/modules/navier_stokes/src/physics/WCNSFVFlowPhysics.C index d10855f61d1c..307147cb578e 100644 --- a/modules/navier_stokes/src/physics/WCNSFVFlowPhysics.C +++ b/modules/navier_stokes/src/physics/WCNSFVFlowPhysics.C @@ -145,7 +145,7 @@ WCNSFVFlowPhysics::WCNSFVFlowPhysics(const InputParameters & parameters) } void -WCNSFVFlowPhysics::addNonlinearVariables() +WCNSFVFlowPhysics::addSolverVariables() { if (!_has_flow_equations) return; diff --git a/modules/navier_stokes/src/physics/WCNSFVFluidHeatTransferPhysics.C b/modules/navier_stokes/src/physics/WCNSFVFluidHeatTransferPhysics.C index 0ccaee22b1ea..c292fb63de9f 100644 --- a/modules/navier_stokes/src/physics/WCNSFVFluidHeatTransferPhysics.C +++ b/modules/navier_stokes/src/physics/WCNSFVFluidHeatTransferPhysics.C @@ -115,7 +115,7 @@ WCNSFVFluidHeatTransferPhysics::WCNSFVFluidHeatTransferPhysics(const InputParame } void -WCNSFVFluidHeatTransferPhysics::addNonlinearVariables() +WCNSFVFluidHeatTransferPhysics::addSolverVariables() { // For compatibility with Modules/NavierStokesFV syntax if (!_has_energy_equation) diff --git a/modules/navier_stokes/src/physics/WCNSFVScalarTransportPhysics.C b/modules/navier_stokes/src/physics/WCNSFVScalarTransportPhysics.C index faa573c393f9..7709d18a7baa 100644 --- a/modules/navier_stokes/src/physics/WCNSFVScalarTransportPhysics.C +++ b/modules/navier_stokes/src/physics/WCNSFVScalarTransportPhysics.C @@ -8,8 +8,7 @@ //* https://www.gnu.org/licenses/lgpl-2.1.html #include "WCNSFVScalarTransportPhysics.h" -#include "WCNSFVCoupledAdvectionPhysicsHelper.h" -#include "WCNSFVFlowPhysics.h" +#include "WCNSFVFlowPhysicsBase.h" #include "NSFVBase.h" #include "NS.h" @@ -19,107 +18,21 @@ registerWCNSFVScalarTransportBaseTasks("NavierStokesApp", WCNSFVScalarTransportP InputParameters WCNSFVScalarTransportPhysics::validParams() { - InputParameters params = NavierStokesPhysicsBase::validParams(); - params += WCNSFVCoupledAdvectionPhysicsHelper::validParams(); - params.addClassDescription( - "Define the Navier Stokes weakly-compressible scalar field transport equation(s)"); - - params += NSFVBase::commonScalarFieldAdvectionParams(); - - // TODO Remove the parameter once NavierStokesFV syntax has been removed - params.addParam( - "add_scalar_equation", - "Whether to add the scalar transport equation. This parameter is not necessary if " - "using the Physics syntax"); - - // These parameters are not shared because the NSFVPhysics use functors - params.addParam>>( - "passive_scalar_inlet_function", - std::vector>(), - "Functors for inlet boundaries in the passive scalar equations."); - - // New functor boundary conditions - params.deprecateParam( - "passive_scalar_inlet_function", "passive_scalar_inlet_functors", "01/01/2025"); - - // No need for the duplication - params.addParam>("passive_scalar_source", "Passive scalar sources"); - - // Spatial finite volume discretization scheme - params.transferParam(NSFVBase::validParams(), - "passive_scalar_advection_interpolation"); + InputParameters params = WCNSFVScalarTransportPhysicsBase::validParams(); + params.addClassDescription("Define the Navier Stokes weakly-compressible scalar field transport " + "equation(s) using the nonlinear finite volume discretization"); params.transferParam(NSFVBase::validParams(), "passive_scalar_face_interpolation"); - params.transferParam(NSFVBase::validParams(), "passive_scalar_two_term_bc_expansion"); - - // Nonlinear equation solver scaling - params.addRangeCheckedParam>( - "passive_scalar_scaling", - "passive_scalar_scaling > 0.0", - "The scaling factor for the passive scalar field variables."); - - // Parameter groups - params.addParamNamesToGroup("passive_scalar_names initial_scalar_variables", "Variable"); - params.addParamNamesToGroup( - "passive_scalar_advection_interpolation passive_scalar_face_interpolation " - "passive_scalar_two_term_bc_expansion passive_scalar_scaling", - "Numerical scheme"); - params.addParamNamesToGroup("passive_scalar_inlet_types passive_scalar_inlet_functors", - "Inlet boundary"); - + params.addParamNamesToGroup("passive_scalar_face_interpolation", "Numerical scheme"); return params; } WCNSFVScalarTransportPhysics::WCNSFVScalarTransportPhysics(const InputParameters & parameters) - : NavierStokesPhysicsBase(parameters), - WCNSFVCoupledAdvectionPhysicsHelper(this), - _passive_scalar_names(getParam>("passive_scalar_names")), - _has_scalar_equation(isParamValid("add_scalar_equation") ? getParam("add_scalar_equation") - : !usingNavierStokesFVSyntax()), - _passive_scalar_inlet_types(getParam("passive_scalar_inlet_types")), - _passive_scalar_inlet_functors( - getParam>>("passive_scalar_inlet_functors")), - _passive_scalar_sources(getParam>("passive_scalar_source")), - _passive_scalar_coupled_sources( - getParam>>("passive_scalar_coupled_source")), - _passive_scalar_sources_coef( - getParam>>("passive_scalar_coupled_source_coeff")) + : WCNSFVScalarTransportPhysicsBase(parameters) { - if (_has_scalar_equation) - for (const auto & scalar_name : _passive_scalar_names) - saveSolverVariableName(scalar_name); - - // For compatibility with Modules/NavierStokesFV syntax - if (!_has_scalar_equation) - return; - - // These parameters must be passed for every passive scalar at a time - checkVectorParamsSameLengthIfSet( - "passive_scalar_names", "passive_scalar_diffusivity", true); - checkVectorParamsSameLengthIfSet>( - "passive_scalar_names", "passive_scalar_coupled_source", true); - checkVectorParamsSameLengthIfSet( - "passive_scalar_names", "passive_scalar_source", true); - checkVectorParamsSameLengthIfSet( - "passive_scalar_names", "passive_scalar_scaling", true); - checkVectorParamsSameLengthIfSet( - "passive_scalar_names", "initial_scalar_variables", true); - checkVectorParamsSameLengthIfSet>( - "passive_scalar_names", "passive_scalar_inlet_functors", true); - if (_passive_scalar_inlet_functors.size()) - checkTwoDVectorParamMultiMooseEnumSameLength( - "passive_scalar_inlet_functors", "passive_scalar_inlet_types", false); - - if (_passive_scalar_sources_coef.size()) - checkTwoDVectorParamsSameLength("passive_scalar_coupled_source", - "passive_scalar_coupled_source_coeff"); - - if (_porous_medium_treatment) - _flow_equations_physics->paramError("porous_medium_treatment", - "Porous media scalar advection is currently unimplemented"); } void -WCNSFVScalarTransportPhysics::addNonlinearVariables() +WCNSFVScalarTransportPhysics::addSolverVariables() { // For compatibility with Modules/NavierStokesFV syntax if (!_has_scalar_equation) @@ -152,22 +65,6 @@ WCNSFVScalarTransportPhysics::addNonlinearVariables() } } -void -WCNSFVScalarTransportPhysics::addFVKernels() -{ - // For compatibility with Modules/NavierStokesFV syntax - if (!_has_scalar_equation) - return; - - if (isTransient()) - addScalarTimeKernels(); - - addScalarAdvectionKernels(); - addScalarDiffusionKernels(); - if (_passive_scalar_sources.size() || _passive_scalar_coupled_sources.size()) - addScalarSourceKernels(); -} - void WCNSFVScalarTransportPhysics::addScalarTimeKernels() { @@ -258,19 +155,6 @@ WCNSFVScalarTransportPhysics::addScalarSourceKernels() } } -void -WCNSFVScalarTransportPhysics::addFVBCs() -{ - // For compatibility with Modules/NavierStokesFV syntax - if (!_has_scalar_equation) - return; - - addScalarInletBC(); - // There is typically no wall flux of passive scalars, similarly we rarely know - // their concentrations at the outlet at the beginning of the simulation - // TODO: we will know the outlet values in case of flow reversal. Implement scalar outlet -} - void WCNSFVScalarTransportPhysics::addScalarInletBC() { @@ -361,48 +245,8 @@ WCNSFVScalarTransportPhysics::addScalarInletBC() } void -WCNSFVScalarTransportPhysics::addInitialConditions() +WCNSFVScalarTransportPhysics::addScalarOutletBC() { - // For compatibility with Modules/NavierStokesFV syntax - if (!_has_scalar_equation) - return; - if (!_define_variables && parameters().isParamSetByUser("initial_scalar_variables")) - paramError("initial_scalar_variables", - "Scalar variables are defined externally of NavierStokesFV, so should their inital " - "conditions"); - // do not set initial conditions if we load from file - if (getParam("initialize_variables_from_mesh_file")) - return; - // do not set initial conditions if we are not defining variables - if (!_define_variables) - return; - - InputParameters params = getFactory().getValidParams("FunctionIC"); - assignBlocks(params, _blocks); - - // There are no default initial conditions for passive scalar variables, we however - // must obey the user-defined initial conditions, even if we are restarting - if (parameters().isParamSetByUser("initial_scalar_variables")) - { - for (unsigned int name_i = 0; name_i < _passive_scalar_names.size(); ++name_i) - { - params.set("variable") = _passive_scalar_names[name_i]; - params.set("function") = - getParam>("initial_scalar_variables")[name_i]; - - getProblem().addInitialCondition("FunctionIC", _passive_scalar_names[name_i] + "_ic", params); - } - } -} - -unsigned short -WCNSFVScalarTransportPhysics::getNumberAlgebraicGhostingLayersNeeded() const -{ - unsigned short necessary_layers = getParam("ghost_layers"); - necessary_layers = - std::max(necessary_layers, _flow_equations_physics->getNumberAlgebraicGhostingLayersNeeded()); - if (getParam("passive_scalar_face_interpolation") == "skewness-corrected") - necessary_layers = std::max(necessary_layers, (unsigned short)3); - - return necessary_layers; + // Advection outlet is naturally handled by the advection flux kernel + return; } diff --git a/modules/navier_stokes/src/physics/WCNSFVScalarTransportPhysicsBase.C b/modules/navier_stokes/src/physics/WCNSFVScalarTransportPhysicsBase.C new file mode 100644 index 000000000000..2926a8b58072 --- /dev/null +++ b/modules/navier_stokes/src/physics/WCNSFVScalarTransportPhysicsBase.C @@ -0,0 +1,193 @@ +//* 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 + +#include "WCNSFVScalarTransportPhysicsBase.h" +#include "WCNSFVCoupledAdvectionPhysicsHelper.h" +#include "WCNSFVFlowPhysics.h" +#include "NSFVBase.h" +#include "NS.h" + +InputParameters +WCNSFVScalarTransportPhysicsBase::validParams() +{ + InputParameters params = NavierStokesPhysicsBase::validParams(); + params += WCNSFVCoupledAdvectionPhysicsHelper::validParams(); + params.addClassDescription( + "Define the Navier Stokes weakly-compressible scalar field transport equation(s)"); + + params += NSFVBase::commonScalarFieldAdvectionParams(); + + // TODO Remove the parameter once NavierStokesFV syntax has been removed + params.addParam( + "add_scalar_equation", + "Whether to add the scalar transport equation. This parameter is not necessary if " + "using the Physics syntax"); + + // These parameters are not shared because the NSFVPhysics use functors + params.addParam>>( + "passive_scalar_inlet_function", + std::vector>(), + "Functors for inlet boundaries in the passive scalar equations."); + + // New functor boundary conditions + params.deprecateParam( + "passive_scalar_inlet_function", "passive_scalar_inlet_functors", "01/01/2025"); + + // No need for the duplication + params.addParam>("passive_scalar_source", "Passive scalar sources"); + + // Spatial finite volume discretization scheme + params.transferParam(NSFVBase::validParams(), + "passive_scalar_advection_interpolation"); + params.transferParam(NSFVBase::validParams(), "passive_scalar_two_term_bc_expansion"); + + // Nonlinear equation solver scaling + params.addRangeCheckedParam>( + "passive_scalar_scaling", + "passive_scalar_scaling > 0.0", + "The scaling factor for the passive scalar field variables."); + + // Parameter groups + params.addParamNamesToGroup("passive_scalar_names initial_scalar_variables", "Variable"); + params.addParamNamesToGroup("passive_scalar_advection_interpolation passive_scalar_scaling " + "passive_scalar_two_term_bc_expansion", + "Numerical scheme"); + params.addParamNamesToGroup("passive_scalar_inlet_types passive_scalar_inlet_functors", + "Inlet boundary"); + + return params; +} + +WCNSFVScalarTransportPhysicsBase::WCNSFVScalarTransportPhysicsBase( + const InputParameters & parameters) + : NavierStokesPhysicsBase(parameters), + WCNSFVCoupledAdvectionPhysicsHelper(this), + _passive_scalar_names(getParam>("passive_scalar_names")), + _has_scalar_equation(isParamValid("add_scalar_equation") ? getParam("add_scalar_equation") + : !usingNavierStokesFVSyntax()), + _passive_scalar_inlet_types(getParam("passive_scalar_inlet_types")), + _passive_scalar_inlet_functors( + getParam>>("passive_scalar_inlet_functors")), + _passive_scalar_sources(getParam>("passive_scalar_source")), + _passive_scalar_coupled_sources( + getParam>>("passive_scalar_coupled_source")), + _passive_scalar_sources_coef( + getParam>>("passive_scalar_coupled_source_coeff")) +{ + if (_has_scalar_equation) + for (const auto & scalar_name : _passive_scalar_names) + saveSolverVariableName(scalar_name); + + // For compatibility with Modules/NavierStokesFV syntax + if (!_has_scalar_equation) + return; + + // These parameters must be passed for every passive scalar at a time + checkVectorParamsSameLengthIfSet( + "passive_scalar_names", "passive_scalar_diffusivity", true); + checkVectorParamsSameLengthIfSet>( + "passive_scalar_names", "passive_scalar_coupled_source", true); + checkVectorParamsSameLengthIfSet( + "passive_scalar_names", "passive_scalar_source", true); + checkVectorParamsSameLengthIfSet( + "passive_scalar_names", "passive_scalar_scaling", true); + checkVectorParamsSameLengthIfSet( + "passive_scalar_names", "initial_scalar_variables", true); + checkVectorParamsSameLengthIfSet>( + "passive_scalar_names", "passive_scalar_inlet_functors", true); + if (_passive_scalar_inlet_functors.size()) + checkTwoDVectorParamMultiMooseEnumSameLength( + "passive_scalar_inlet_functors", "passive_scalar_inlet_types", false); + + if (_passive_scalar_sources_coef.size()) + checkTwoDVectorParamsSameLength("passive_scalar_coupled_source", + "passive_scalar_coupled_source_coeff"); + + if (_porous_medium_treatment) + _flow_equations_physics->paramError("porous_medium_treatment", + "Porous media scalar advection is currently unimplemented"); +} + +void +WCNSFVScalarTransportPhysicsBase::addFVKernels() +{ + // For compatibility with Modules/NavierStokesFV syntax + if (!_has_scalar_equation) + return; + + if (isTransient()) + addScalarTimeKernels(); + + addScalarAdvectionKernels(); + addScalarDiffusionKernels(); + if (_passive_scalar_sources.size() || _passive_scalar_coupled_sources.size()) + addScalarSourceKernels(); +} + +void +WCNSFVScalarTransportPhysicsBase::addFVBCs() +{ + // For compatibility with Modules/NavierStokesFV syntax + if (!_has_scalar_equation) + return; + + addScalarInletBC(); + // There is typically no wall flux of passive scalars, similarly we rarely know + // their concentrations at the outlet at the beginning of the simulation + // TODO: we will know the outlet values in case of flow reversal. Implement scalar outlet + addScalarWallBC(); + addScalarOutletBC(); +} + +void +WCNSFVScalarTransportPhysicsBase::addInitialConditions() +{ + // For compatibility with Modules/NavierStokesFV syntax + if (!_has_scalar_equation) + return; + if (!_define_variables && parameters().isParamSetByUser("initial_scalar_variables")) + paramError("initial_scalar_variables", + "Scalar variables are defined externally of NavierStokesFV, so should their inital " + "conditions"); + // do not set initial conditions if we load from file + if (getParam("initialize_variables_from_mesh_file")) + return; + // do not set initial conditions if we are not defining variables + if (!_define_variables) + return; + + InputParameters params = getFactory().getValidParams("FunctionIC"); + assignBlocks(params, _blocks); + + // There are no default initial conditions for passive scalar variables, we however + // must obey the user-defined initial conditions, even if we are restarting + if (parameters().isParamSetByUser("initial_scalar_variables")) + { + for (unsigned int name_i = 0; name_i < _passive_scalar_names.size(); ++name_i) + { + params.set("variable") = _passive_scalar_names[name_i]; + params.set("function") = + getParam>("initial_scalar_variables")[name_i]; + + getProblem().addInitialCondition("FunctionIC", _passive_scalar_names[name_i] + "_ic", params); + } + } +} + +unsigned short +WCNSFVScalarTransportPhysicsBase::getNumberAlgebraicGhostingLayersNeeded() const +{ + unsigned short necessary_layers = getParam("ghost_layers"); + necessary_layers = + std::max(necessary_layers, _flow_equations_physics->getNumberAlgebraicGhostingLayersNeeded()); + if (getParam("passive_scalar_advection_interpolation") == "skewness-corrected") + necessary_layers = std::max(necessary_layers, (unsigned short)3); + + return necessary_layers; +} diff --git a/modules/navier_stokes/src/physics/WCNSFVTurbulencePhysics.C b/modules/navier_stokes/src/physics/WCNSFVTurbulencePhysics.C index 5bf5a7a7a1ba..15898c2aa58d 100644 --- a/modules/navier_stokes/src/physics/WCNSFVTurbulencePhysics.C +++ b/modules/navier_stokes/src/physics/WCNSFVTurbulencePhysics.C @@ -330,7 +330,7 @@ WCNSFVTurbulencePhysics::retrieveCoupledPhysics() } void -WCNSFVTurbulencePhysics::addNonlinearVariables() +WCNSFVTurbulencePhysics::addSolverVariables() { if (_turbulence_model == "mixing-length" || _turbulence_model == "none") return; diff --git a/modules/navier_stokes/src/physics/WCNSLinearFVFlowPhysics.C b/modules/navier_stokes/src/physics/WCNSLinearFVFlowPhysics.C index 555d583905fe..0a3046682feb 100644 --- a/modules/navier_stokes/src/physics/WCNSLinearFVFlowPhysics.C +++ b/modules/navier_stokes/src/physics/WCNSLinearFVFlowPhysics.C @@ -83,7 +83,7 @@ WCNSLinearFVFlowPhysics::initializePhysicsAdditional() } void -WCNSLinearFVFlowPhysics::addNonlinearVariables() +WCNSLinearFVFlowPhysics::addSolverVariables() { if (!_has_flow_equations) return; @@ -279,8 +279,8 @@ WCNSLinearFVFlowPhysics::addINSInletBC() { // Check the size of the BC parameters unsigned int num_velocity_functor_inlets = 0; - for (const auto & [bdy, momentum_outlet_type] : _momentum_inlet_types) - if (momentum_outlet_type == "fixed-velocity" || momentum_outlet_type == "fixed-pressure") + for (const auto & [bdy, momentum_inlet_type] : _momentum_inlet_types) + if (momentum_inlet_type == "fixed-velocity" || momentum_inlet_type == "fixed-pressure") num_velocity_functor_inlets++; if (num_velocity_functor_inlets != _momentum_inlet_functors.size()) diff --git a/modules/navier_stokes/src/physics/WCNSLinearFVScalarTransportPhysics.C b/modules/navier_stokes/src/physics/WCNSLinearFVScalarTransportPhysics.C new file mode 100644 index 000000000000..7fd0e90fc64e --- /dev/null +++ b/modules/navier_stokes/src/physics/WCNSLinearFVScalarTransportPhysics.C @@ -0,0 +1,228 @@ +//* 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 + +#include "WCNSLinearFVScalarTransportPhysics.h" +#include "WCNSFVFlowPhysicsBase.h" +#include "NS.h" + +registerNavierStokesPhysicsBaseTasks("NavierStokesApp", WCNSLinearFVScalarTransportPhysics); +registerWCNSFVScalarTransportBaseTasks("NavierStokesApp", WCNSLinearFVScalarTransportPhysics); + +InputParameters +WCNSLinearFVScalarTransportPhysics::validParams() +{ + InputParameters params = WCNSFVScalarTransportPhysicsBase::validParams(); + params.addClassDescription("Define the Navier Stokes weakly-compressible scalar field transport " + "equation(s) using the linear finite volume discretization"); + params.addParam("use_nonorthogonal_correction", + true, + "If the nonorthogonal correction should be used when computing the normal " + "gradient, notably in the diffusion term."); + return params; +} + +WCNSLinearFVScalarTransportPhysics::WCNSLinearFVScalarTransportPhysics( + const InputParameters & parameters) + : WCNSFVScalarTransportPhysicsBase(parameters) +{ +} + +void +WCNSLinearFVScalarTransportPhysics::addSolverVariables() +{ + // For compatibility with Modules/NavierStokesFV syntax + if (!_has_scalar_equation) + return; + + auto params = getFactory().getValidParams("MooseLinearVariableFVReal"); + assignBlocks(params, _blocks); + + for (const auto name_i : index_range(_passive_scalar_names)) + { + // Dont add if the user already defined the variable + if (variableExists(_passive_scalar_names[name_i], /*error_if_aux=*/true)) + { + checkBlockRestrictionIdentical( + _passive_scalar_names[name_i], + getProblem().getVariable(0, _passive_scalar_names[name_i]).blocks()); + continue; + } + + params.set("solver_sys") = getSolverSystem(name_i); + if (isParamValid("passive_scalar_scaling")) + params.set>("scaling") = { + getParam>("passive_scalar_scaling")[name_i]}; + + getProblem().addVariable("MooseLinearVariableFVReal", _passive_scalar_names[name_i], params); + } +} + +void +WCNSLinearFVScalarTransportPhysics::addScalarTimeKernels() +{ + paramError("transient", "Transient simulations are not supported at this time"); +} + +void +WCNSLinearFVScalarTransportPhysics::addScalarAdvectionKernels() +{ + const std::string kernel_type = "LinearFVScalarAdvection"; + InputParameters params = getFactory().getValidParams(kernel_type); + + assignBlocks(params, _blocks); + params.set("rhie_chow_user_object") = _flow_equations_physics->rhieChowUOName(); + params.set("advected_interp_method") = + getParam("passive_scalar_advection_interpolation"); + setSlipVelocityParams(params); + + for (const auto & vname : _passive_scalar_names) + { + params.set("variable") = vname; + getProblem().addLinearFVKernel(kernel_type, prefix() + "ins_" + vname + "_advection", params); + } +} + +void +WCNSLinearFVScalarTransportPhysics::addScalarDiffusionKernels() +{ + // Direct specification of diffusion term + const auto passive_scalar_diffusivities = + getParam>("passive_scalar_diffusivity"); + + if (passive_scalar_diffusivities.size()) + { + const std::string kernel_type = "LinearFVDiffusion"; + InputParameters params = getFactory().getValidParams(kernel_type); + assignBlocks(params, _blocks); + params.set("use_nonorthogonal_correction") = + getParam("use_nonorthogonal_correction"); + for (const auto name_i : index_range(_passive_scalar_names)) + { + params.set("variable") = _passive_scalar_names[name_i]; + params.set("diffusion_coeff") = passive_scalar_diffusivities[name_i]; + getProblem().addLinearFVKernel( + kernel_type, prefix() + "ins_" + _passive_scalar_names[name_i] + "_diffusion", params); + } + } +} + +void +WCNSLinearFVScalarTransportPhysics::addScalarSourceKernels() +{ + const std::string kernel_type = "LinearFVSource"; + InputParameters params = getFactory().getValidParams(kernel_type); + assignBlocks(params, _blocks); + + for (const auto scalar_i : index_range(_passive_scalar_names)) + { + params.set("variable") = _passive_scalar_names[scalar_i]; + + if (_passive_scalar_sources.size()) + { + // Added for backward compatibility with former Modules/NavierStokesFV syntax + params.set("source_density") = _passive_scalar_sources[scalar_i]; + getProblem().addFVKernel( + kernel_type, prefix() + "ins_" + _passive_scalar_names[scalar_i] + "_source", params); + } + + if (_passive_scalar_coupled_sources.size()) + for (const auto i : index_range(_passive_scalar_coupled_sources[scalar_i])) + { + params.set("source_density") = + _passive_scalar_coupled_sources[scalar_i][i]; + if (_passive_scalar_sources_coef.size()) + params.set("scaling_factor") = _passive_scalar_sources_coef[scalar_i][i]; + + getProblem().addLinearFVKernel(kernel_type, + prefix() + "ins_" + _passive_scalar_names[scalar_i] + + "_coupled_source_" + std::to_string(i), + params); + } + } +} + +void +WCNSLinearFVScalarTransportPhysics::addScalarInletBC() +{ + const auto & inlet_boundaries = _flow_equations_physics->getInletBoundaries(); + if (inlet_boundaries.empty()) + return; + + // Boundary checks + // TODO: once we have vectors of MooseEnum, we could use the same templated check for types and + // functors + if (inlet_boundaries.size() * _passive_scalar_names.size() != _passive_scalar_inlet_types.size()) + paramError( + "passive_scalar_inlet_types", + "The number of scalar inlet types (" + std::to_string(_passive_scalar_inlet_types.size()) + + ") is not equal to the number of inlet boundaries (" + + std::to_string(inlet_boundaries.size()) + ") times the number of passive scalars (" + + std::to_string(_passive_scalar_names.size()) + ")"); + if (_passive_scalar_names.size() != _passive_scalar_inlet_functors.size()) + paramError("passive_scalar_inlet_functors", + "The number of groups of inlet functors (" + + std::to_string(_passive_scalar_inlet_functors.size()) + + ") is not equal to the number of passive scalars (" + + std::to_string(_passive_scalar_names.size()) + ")"); + + for (const auto name_i : index_range(_passive_scalar_names)) + { + if (inlet_boundaries.size() != _passive_scalar_inlet_functors[name_i].size()) + paramError("passive_scalar_inlet_functors", + "The number of inlet boundary functors for scalar '" + + _passive_scalar_names[name_i] + + "' does not match the number of inlet boundaries (" + + std::to_string(_passive_scalar_inlet_functors[name_i].size()) + ")"); + + unsigned int num_inlets = inlet_boundaries.size(); + for (unsigned int bc_ind = 0; bc_ind < num_inlets; ++bc_ind) + { + if (_passive_scalar_inlet_types[name_i * num_inlets + bc_ind] == "fixed-value") + { + const std::string bc_type = "LinearFVAdvectionDiffusionFunctorDirichletBC"; + InputParameters params = getFactory().getValidParams(bc_type); + params.set("variable") = _passive_scalar_names[name_i]; + params.set("functor") = _passive_scalar_inlet_functors[name_i][bc_ind]; + params.set>("boundary") = {inlet_boundaries[bc_ind]}; + + getProblem().addLinearFVBC( + bc_type, _passive_scalar_names[name_i] + "_" + inlet_boundaries[bc_ind], params); + } + else if (_passive_scalar_inlet_types[name_i * num_inlets + bc_ind] == "flux-mass" || + _passive_scalar_inlet_types[name_i * num_inlets + bc_ind] == "flux-velocity") + { + mooseError("Flux boundary conditions not supported at this time using the linear finite " + "volume discretization"); + } + } + } +} + +void +WCNSLinearFVScalarTransportPhysics::addScalarOutletBC() +{ + const auto & outlet_boundaries = _flow_equations_physics->getOutletBoundaries(); + if (outlet_boundaries.empty()) + return; + + for (const auto & outlet_bdy : outlet_boundaries) + { + const std::string bc_type = "LinearFVAdvectionDiffusionOutflowBC"; + InputParameters params = getFactory().getValidParams(bc_type); + params.set>("boundary") = {outlet_bdy}; + params.set("use_two_term_expansion") = + getParam("passive_scalar_two_term_bc_expansion"); + + for (const auto name_i : index_range(_passive_scalar_names)) + { + params.set("variable") = _passive_scalar_names[name_i]; + getProblem().addLinearFVBC(bc_type, _passive_scalar_names[name_i] + "_" + outlet_bdy, params); + } + } +} diff --git a/modules/navier_stokes/src/userobjects/RhieChowMassFlux.C b/modules/navier_stokes/src/userobjects/RhieChowMassFlux.C index 4e33d9722d40..a0cc5e3a33cd 100644 --- a/modules/navier_stokes/src/userobjects/RhieChowMassFlux.C +++ b/modules/navier_stokes/src/userobjects/RhieChowMassFlux.C @@ -234,6 +234,19 @@ RhieChowMassFlux::getMassFlux(const FaceInfo & fi) const return _face_mass_flux.evaluate(&fi); } +Real +RhieChowMassFlux::getVolumetricFaceFlux(const FaceInfo & fi) const +{ + const Moose::FaceArg face_arg{&fi, + /*limiter_type=*/Moose::FV::LimiterType::CentralDifference, + /*elem_is_upwind=*/true, + /*correct_skewness=*/false, + &fi.elem(), + /*state_limiter*/ nullptr}; + const Real face_rho = _rho(face_arg, Moose::currentState()); + return libmesh_map_find(_face_mass_flux, fi.id()) / face_rho; +} + void RhieChowMassFlux::computeFaceMassFlux() { @@ -293,7 +306,7 @@ RhieChowMassFlux::computeFaceMassFlux() const auto rhs_contribution = _p_diffusion_kernel->computeBoundaryRHSContribution(*bc_pointer); - // On the boundary, only the element side has a constibution + // On the boundary, only the element side has a contribution p_grad_flux = (p_elem_value * matrix_contribution - rhs_contribution); } // Compute the new face flux diff --git a/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d-scalar/channel-physics.i b/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d-scalar/channel-physics.i new file mode 100644 index 000000000000..11f7f2117f3e --- /dev/null +++ b/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d-scalar/channel-physics.i @@ -0,0 +1,104 @@ +mu = 2.6 +rho = 1.0 +advected_interp_method = 'upwind' +k1 = 0.1 +k2 = 0.2 + +[Mesh] + [mesh] + type = CartesianMeshGenerator + dim = 2 + dx = '0.25 0.25' + dy = '0.2' + ix = '5 5' + iy = '5' + subdomain_id = '0 1' + [] +[] + +[Problem] + linear_sys_names = 'u_system v_system pressure_system s1_system s2_system' +[] + + +[Physics] + [NavierStokes] + [FlowSegregated/flow] + velocity_variable = 'vel_x vel_y' + pressure_variable = 'pressure' + + initial_velocity = '0.5 0 0' + initial_pressure = '0.2' + + density = ${rho} + dynamic_viscosity = ${mu} + + # use inlet for moving wall to match the reference input + # we could also use a noslip BC with a velocity wall functor + inlet_boundaries = 'left' + momentum_inlet_types = 'fixed-velocity' + momentum_inlet_functors = '1.1 0' + + wall_boundaries = 'top bottom' + momentum_wall_types = 'noslip noslip' + + outlet_boundaries = 'right' + momentum_outlet_types = 'fixed-pressure' + pressure_functors = '1.4' + + orthogonality_correction = false + pressure_two_term_bc_expansion = false + momentum_advection_interpolation = ${advected_interp_method} + [] + [ScalarTransportSegregated/scalar] + passive_scalar_names = 'scalar1 scalar2' + system_names = 's1_system s2_system' + initial_scalar_variables = '1.1 3' + + passive_scalar_diffusivity = '${k1} ${k2}' + + passive_scalar_inlet_types = 'fixed-value fixed-value' + passive_scalar_inlet_function = '1; 2' + + passive_scalar_advection_interpolation = ${advected_interp_method} + passive_scalar_two_term_bc_expansion = false + use_nonorthogonal_correction = false + [] + [] +[] + +[Executioner] + type = SIMPLE + momentum_l_abs_tol = 1e-13 + pressure_l_abs_tol = 1e-13 + passive_scalar_l_abs_tol = 1e-13 + momentum_l_tol = 0 + pressure_l_tol = 0 + passive_scalar_l_tol = 0 + rhie_chow_user_object = 'ins_rhie_chow_interpolator' + momentum_systems = 'u_system v_system' + pressure_system = 'pressure_system' + passive_scalar_systems = 's1_system s2_system' + momentum_equation_relaxation = 0.8 + passive_scalar_equation_relaxation = '0.9 0.9' + pressure_variable_relaxation = 0.3 + # We need to converge the problem to show conservation + num_iterations = 200 + pressure_absolute_tolerance = 1e-10 + momentum_absolute_tolerance = 1e-10 + passive_scalar_absolute_tolerance = '1e-10 1e-10' + momentum_petsc_options_iname = '-pc_type -pc_hypre_type' + momentum_petsc_options_value = 'hypre boomeramg' + pressure_petsc_options_iname = '-pc_type -pc_hypre_type' + pressure_petsc_options_value = 'hypre boomeramg' + passive_scalar_petsc_options_iname = '-pc_type -pc_hypre_type' + passive_scalar_petsc_options_value = 'hypre boomeramg' + print_fields = false + continue_on_max_its = true +[] + +[Outputs] + exodus = true + execute_on = timestep_end + hide = 'pressure vel_x vel_y' +[] diff --git a/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d-scalar/channel.i b/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d-scalar/channel.i new file mode 100644 index 000000000000..b5537906af8a --- /dev/null +++ b/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d-scalar/channel.i @@ -0,0 +1,241 @@ +mu = 2.6 +rho = 1.0 +advected_interp_method = 'upwind' +k1 = 0.1 +k2 = 0.2 + +[Mesh] + [mesh] + type = CartesianMeshGenerator + dim = 2 + dx = '0.25 0.25' + dy = '0.2' + ix = '5 5' + iy = '5' + subdomain_id = '0 1' + [] +[] + +[Problem] + linear_sys_names = 'u_system v_system pressure_system s1_system s2_system' + previous_nl_solution_required = true +[] + +[UserObjects] + [rc] + type = RhieChowMassFlux + u = vel_x + v = vel_y + pressure = pressure + rho = ${rho} + p_diffusion_kernel = p_diffusion + [] +[] + +[Variables] + [vel_x] + type = MooseLinearVariableFVReal + initial_condition = 0.5 + solver_sys = u_system + [] + [vel_y] + type = MooseLinearVariableFVReal + solver_sys = v_system + initial_condition = 0.0 + [] + [pressure] + type = MooseLinearVariableFVReal + solver_sys = pressure_system + initial_condition = 0.2 + [] + [scalar1] + type = MooseLinearVariableFVReal + solver_sys = s1_system + initial_condition = 1.1 + [] + [scalar2] + type = MooseLinearVariableFVReal + solver_sys = s2_system + initial_condition = 3 + [] +[] + +[LinearFVKernels] + [u_advection_stress] + type = LinearWCNSFVMomentumFlux + variable = vel_x + advected_interp_method = ${advected_interp_method} + mu = ${mu} + u = vel_x + v = vel_y + momentum_component = 'x' + rhie_chow_user_object = 'rc' + use_nonorthogonal_correction = false + [] + [v_advection_stress] + type = LinearWCNSFVMomentumFlux + variable = vel_y + advected_interp_method = ${advected_interp_method} + mu = ${mu} + u = vel_x + v = vel_y + momentum_component = 'y' + rhie_chow_user_object = 'rc' + use_nonorthogonal_correction = false + [] + [u_pressure] + type = LinearFVMomentumPressure + variable = vel_x + pressure = pressure + momentum_component = 'x' + [] + [v_pressure] + type = LinearFVMomentumPressure + variable = vel_y + pressure = pressure + momentum_component = 'y' + [] + + [p_diffusion] + type = LinearFVAnisotropicDiffusion + variable = pressure + diffusion_tensor = Ainv + use_nonorthogonal_correction = false + [] + [HbyA_divergence] + type = LinearFVDivergence + variable = pressure + face_flux = HbyA + force_boundary_execution = true + [] + + [s1_advection] + type = LinearFVScalarAdvection + variable = scalar1 + advected_interp_method = ${advected_interp_method} + rhie_chow_user_object = 'rc' + [] + [s1_diffusion] + type = LinearFVDiffusion + variable = scalar1 + diffusion_coeff = ${k1} + use_nonorthogonal_correction = false + [] + [s2_diffusion] + type = LinearFVScalarAdvection + variable = scalar2 + advected_interp_method = ${advected_interp_method} + rhie_chow_user_object = 'rc' + [] + [s2_conduction] + type = LinearFVDiffusion + variable = scalar2 + diffusion_coeff = ${k2} + use_nonorthogonal_correction = false + [] +[] + +[LinearFVBCs] + [inlet-u] + type = LinearFVAdvectionDiffusionFunctorDirichletBC + boundary = 'left' + variable = vel_x + functor = '1.1' + [] + [inlet-v] + type = LinearFVAdvectionDiffusionFunctorDirichletBC + boundary = 'left' + variable = vel_y + functor = '0.0' + [] + [walls-u] + type = LinearFVAdvectionDiffusionFunctorDirichletBC + boundary = 'top bottom' + variable = vel_x + functor = 0.0 + [] + [walls-v] + type = LinearFVAdvectionDiffusionFunctorDirichletBC + boundary = 'top bottom' + variable = vel_y + functor = 0.0 + [] + [outlet_p] + type = LinearFVAdvectionDiffusionFunctorDirichletBC + boundary = 'right' + variable = pressure + functor = 1.4 + [] + [outlet_u] + type = LinearFVAdvectionDiffusionOutflowBC + variable = vel_x + use_two_term_expansion = false + boundary = right + [] + [outlet_v] + type = LinearFVAdvectionDiffusionOutflowBC + variable = vel_y + use_two_term_expansion = false + boundary = right + [] + [inlet_wall_scalar1] + type = LinearFVAdvectionDiffusionFunctorDirichletBC + variable = scalar1 + functor = 1 + boundary = 'left' + [] + [outlet_scalar1] + type = LinearFVAdvectionDiffusionOutflowBC + variable = scalar1 + use_two_term_expansion = false + boundary = right + [] + [inlet_wall_scalar2] + type = LinearFVAdvectionDiffusionFunctorDirichletBC + variable = scalar2 + functor = 2 + boundary = 'left' + [] + [outlet_scalar2] + type = LinearFVAdvectionDiffusionOutflowBC + variable = scalar2 + use_two_term_expansion = false + boundary = right + [] +[] + +[Executioner] + type = SIMPLE + momentum_l_abs_tol = 1e-13 + pressure_l_abs_tol = 1e-13 + passive_scalar_l_abs_tol = 1e-13 + momentum_l_tol = 0 + pressure_l_tol = 0 + passive_scalar_l_tol = 0 + rhie_chow_user_object = 'rc' + momentum_systems = 'u_system v_system' + pressure_system = 'pressure_system' + passive_scalar_systems = 's1_system s2_system' + momentum_equation_relaxation = 0.8 + passive_scalar_equation_relaxation = '0.9 0.9' + pressure_variable_relaxation = 0.3 + # We need to converge the problem to show conservation + num_iterations = 200 + pressure_absolute_tolerance = 1e-10 + momentum_absolute_tolerance = 1e-10 + passive_scalar_absolute_tolerance = '1e-10 1e-10' + momentum_petsc_options_iname = '-pc_type -pc_hypre_type' + momentum_petsc_options_value = 'hypre boomeramg' + pressure_petsc_options_iname = '-pc_type -pc_hypre_type' + pressure_petsc_options_value = 'hypre boomeramg' + passive_scalar_petsc_options_iname = '-pc_type -pc_hypre_type' + passive_scalar_petsc_options_value = 'hypre boomeramg' + print_fields = false + continue_on_max_its = true +[] + +[Outputs] + exodus = true + execute_on = timestep_end + hide = 'pressure vel_x vel_y' +[] diff --git a/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d-scalar/gold/channel-physics_out.e b/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d-scalar/gold/channel-physics_out.e new file mode 120000 index 000000000000..ffa5923ca413 --- /dev/null +++ b/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d-scalar/gold/channel-physics_out.e @@ -0,0 +1 @@ +channel_out.e \ No newline at end of file diff --git a/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d-scalar/gold/channel_out.e b/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d-scalar/gold/channel_out.e new file mode 100644 index 000000000000..155370847a0d Binary files /dev/null and b/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d-scalar/gold/channel_out.e differ diff --git a/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d-scalar/tests b/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d-scalar/tests new file mode 100644 index 000000000000..aeef8237c5b4 --- /dev/null +++ b/modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d-scalar/tests @@ -0,0 +1,23 @@ +[Tests] + issues = '#28819' + design = 'SIMPLE.md LinearFVScalarAdvection.md' + [channel] + requirement = "The system shall be able to solve a scalar advection problem," + [simple] + type = 'Exodiff' + input = channel.i + exodiff = 'channel_out.e' + recover = false # we don't support recovery for SIMPLE yet + max_threads = 1 # see libmesh issue #3808 + detail = "using the kernel-based syntax," + [] + [physics] + type = 'Exodiff' + input = channel-physics.i + exodiff = 'channel-physics_out.e' + recover = false # we don't support recovery for SIMPLE yet + max_threads = 1 # see libmesh issue #3808 + detail = "using the shorthand Physics syntax." + [] + [] +[] diff --git a/modules/scalar_transport/include/physics/MultiSpeciesDiffusionCG.h b/modules/scalar_transport/include/physics/MultiSpeciesDiffusionCG.h index 8464c3ecf7fd..0d7ee2da158e 100644 --- a/modules/scalar_transport/include/physics/MultiSpeciesDiffusionCG.h +++ b/modules/scalar_transport/include/physics/MultiSpeciesDiffusionCG.h @@ -23,7 +23,7 @@ class MultiSpeciesDiffusionCG : public MultiSpeciesDiffusionPhysicsBase MultiSpeciesDiffusionCG(const InputParameters & parameters); protected: - virtual void addNonlinearVariables() override; + virtual void addSolverVariables() override; virtual void addFEKernels() override; virtual void addFEBCs() override; }; diff --git a/modules/scalar_transport/src/physics/MultiSpeciesDiffusionCG.C b/modules/scalar_transport/src/physics/MultiSpeciesDiffusionCG.C index 614aec8871dd..a660d7ddeb63 100644 --- a/modules/scalar_transport/src/physics/MultiSpeciesDiffusionCG.C +++ b/modules/scalar_transport/src/physics/MultiSpeciesDiffusionCG.C @@ -228,7 +228,7 @@ MultiSpeciesDiffusionCG::addFEBCs() } void -MultiSpeciesDiffusionCG::addNonlinearVariables() +MultiSpeciesDiffusionCG::addSolverVariables() { for (const auto & var_name : _species_names) {