diff --git a/framework/doc/content/bib/moose.bib b/framework/doc/content/bib/moose.bib index 40183a629b63..6d947f08e73d 100644 --- a/framework/doc/content/bib/moose.bib +++ b/framework/doc/content/bib/moose.bib @@ -621,3 +621,22 @@ @article{rao2018stopping url = {https://www.sciencedirect.com/science/article/pii/S0021999117306939}, author = {Kaustubh Rao and Paul Malan and J. Blair Perot} } + +@inproceedings{venkatakrishnan1993, + title={On the accuracy of limiters and convergence to steady state solutions}, + author={Venkatakrishnan, Venkat}, + booktitle={31st Aerospace Sciences Meeting}, + pages={880}, + year={1993} +} + +@article{harten1997, + title={High resolution schemes for hyperbolic conservation laws}, + author={Harten, Ami}, + journal={Journal of computational physics}, + volume={135}, + number={2}, + pages={260--278}, + year={1997}, + publisher={Elsevier} +} diff --git a/framework/doc/content/syntax/Limiters/index.md b/framework/doc/content/syntax/Limiters/index.md index a08a469d9bce..3e631e869d04 100644 --- a/framework/doc/content/syntax/Limiters/index.md +++ b/framework/doc/content/syntax/Limiters/index.md @@ -20,7 +20,10 @@ and $u_j^n = u(x_j,t^n)$. A numerical method is TVD if TV(u^{n+1}) \leq TV(u^n) \end{equation} -## Limiting Process +Different formulations are used for compressible and incompressible/weakly-compressible +flow. + +## Limiting Process for Compressible Flow Borrowing notation from [!citep](greenshields2010implementation), we will now discuss computation of limited quantities, represented by $\bm{\Psi}_{f\pm}$ where @@ -66,7 +69,7 @@ non-skewed mesh the definition in [eq:weighting] and [!citep](greenshields2010implementation) are the same. The flux limiter function $\beta(r_{\pm})$ takes different forms as shown in -[limiter_summary]. $r_{\pm}$ is computed as follows +[limiter_summary_compressible]. $r_{\pm}$ is computed as follows \begin{equation} r_{\pm} = 2 \frac{\bm{d}_{\pm}\cdot\left(\nabla @@ -81,14 +84,95 @@ The following limiters are available in MOOSE. We have noted the convergence orders of each (when considering that the solution is smooth), whether they are TVD, and what the functional form of the flux limiting function $\beta(r)$ is. -!table id=limiter_summary caption=Limiter Overview -|Limiter class name | Convergence Order | TVD | $\beta(r)$ | -|------------------ | ----------------- | --- | --- | -| `VanLeer` | 2 | Yes | $\frac{r +\text{abs}(r)}{1 + \text{abs}(r)}$ | -| `Upwind` | 1 | Yes | 0 | -| `CentralDifference` | 2 | No | 1 | -| `MinMod` | 2 | Yes | $\text{max}(0, \text{min}(1, r))$ | -| `SOU` | 2 | No | $r$ | -| `QUICK` | 2 | No | $\frac{3+r}{4}$ | +!table id=limiter_summary_compressible caption=Compressible Limiter Overview +| Limiter class name | Convergence Order | TVD | $\beta(r)$ | +| ------------------- | ----------------- | --- | -------------------------------------------- | +| `VanLeer` | 2 | Yes | $\frac{r +\text{abs}(r)}{1 + \text{abs}(r)}$ | +| `Upwind` | 1 | Yes | 0 | +| `CentralDifference` | 2 | No | 1 | +| `MinMod` | 2 | Yes | $\text{max}(0, \text{min}(1, r))$ | +| `SOU` | 2 | No | $r$ | +| `QUICK` | 2 | No | $\frac{3+r}{4}$ | + +## Limiting Process for Incompressible and Weakly-Compressible flow + +A full second-order upwind reconstruction is used for incompressible and weakly-compressible solvers. In this reconstruction, the limited quantity at the face is expressed as follows: + +\begin{equation} +\bm{\Psi}_f = \bm{\Psi}_C + \beta(r) ((\nabla \bm{\Psi})_C \cdot \bm{d}_{fC}) +\end{equation} + +where: + +- $\bm{\Psi}_f$ is the value of the variable at the face +- $\bm{\Psi}_C$ is the value of the variable at the cell +- $(\nabla \bm{\Psi})_C$ is the value of the gradient at the cell, which is computed with second-order methods (Green-Gauss without skewness correction and Least-Squares for skewness corrected) +- $\bm{d}_{fC}$ is the distance vector from the face to the cell used in the interpolation +- $\beta(r)$ is the limiting function + +Two kinds of limiters are supported: slope-limited and face-value limited. These limiters are defined below. + +For slope-limiting, the approximate gradient ratio (or flux limiting ratio) $r$ is defined as follows: + +\begin{equation} +r = 2 \frac{\bm{d}_{NC} \cdot (\nabla \bm{\Psi})_C}{\bm{d}_{NC} \cdot (\nabla \bm{\Psi})_f} - 1 +\end{equation} + +where: + +- $\bm{d}_{NC}$ is the vector between the neighbor and current cell adjacent to the face +- $(\nabla \bm{\Psi})_f$ is the gradient of the variable at the face, which is computed by linear interpolation of second-order gradients at the adjacent cells to the face + +For face-value limiting, the limiting function is defined as follows: + +\begin{equation} +r = +\begin{cases} + \frac{|\Delta_f|}{\Delta_{\text{max}}} & \text{if } \Delta_f > 0 \\ + \frac{|\Delta_f|}{\Delta_{\text{min}}} & \text{if } \Delta_f \leq 0 +\end{cases} +\end{equation} + +where: + +- $\Delta_f = (\nabla \bm{\Psi})_C \cdot \bm{d}_{fC}$ is the increment at the face +- $\Delta_{\text{max}} = \bm{\Psi}_{\text{max}} - \bm{\Psi}_C$ is the maximum increment +- $\Delta_{\text{min}} = \bm{\Psi}_{\text{min}} - \bm{\Psi}_C$ is the minimum increment + +The maximum and minimum variable values, $\Delta_{\text{max}}$ and $\Delta_{\text{min}}$, respectively, are computed with a two-cell stencil. In this method, the maximum value is determined as the maximum cell value of the two faces adjacent to the face and their neighbors, respectively. Similarly, the minimum value is computed as the minimum cell value for these cells. + +Each of the limiters implemented along with the implementation reference, limiting type, whether they are TVD, and the functional form of the flux limiting function $\beta(r)$ is shown in [limiter_summary_incompressible]. + +!table id=limiter_summary_incompressible caption=Incompressible/Weakly-Compressible Limiter Overview +| Limiter class name | Limiting Type | TVD | $\beta(r)$ | +| ----------------------------------------------- | ------------- | --- | ----------------------------------------------------------------- | +| `VanLeer` [!citep](harten1997) | Slope | Yes | $\frac{r +\text{abs}(r)}{1 + \text{abs}(r)}$ | +| `MinMod` [!citep](harten1997) | Slope | Yes | $\text{max}(0, \text{min}(1, r))$ | +| `QUICK` [!citep](harten1997) | Slope | Yes | $\text{min}(1,\text{max}(\text{min}(\frac{1 + 3r}{4}, 2r, 2),0))$ | +| `SOU` [!citep](harten1997) | Face-Value | No | $\text{min}(1,1/r)$ | +| `Venkatakrishnan` [!citep](venkatakrishnan1993) | Face-Value | No | $\frac{2r+1}{r(2r+1)+1}$ | + +To illustrate the performance of the limiters, a dispersion analysis is developedand presented in [dispersion]. +This consists of the advection of a passive scalar in a Cartesian mesh at 45 degrees. +The exact solution, without numerical diffusion, is a straight line at 45 degrees +dividing the regions with a scalar concentration of 1 and 0. + +!alert note +In general, we recomment using `VanLeer` and `MinMod` limiters for most of the +applications considering that they provide truly bounded solutions. + +!media framework/finite_volume/dispersion.png + style=display: block;margin-left:auto;margin-right:auto;width:40%; + id=dispersion + caption=Dispersion problem, advection in a Cartesian mesh at 45 degrees. + +The results and performance of each of the limiters are shown in [dispersion_line]. +This image provides an idea of the limiting action and results that +can be expected for each of the limiters. + +!media framework/finite_volume/dispersion_line.png + style=display: block;margin-left:auto;margin-right:auto;width:40%; + id=dispersion_line + caption=Performance of each of the limiters in a line perpendicular to the advection front. !bibtex bibliography diff --git a/framework/include/base/MooseFunctorArguments.h b/framework/include/base/MooseFunctorArguments.h index c87371ed3c62..74174e8bfdaa 100644 --- a/framework/include/base/MooseFunctorArguments.h +++ b/framework/include/base/MooseFunctorArguments.h @@ -75,6 +75,94 @@ struct ElemPointArg ElemArg makeElem() const { return {elem, correct_skewness}; } }; +/** + * Argument for requesting functor evaluation at quadrature point locations on an element side. + * Data in the argument: + * - The element + * - The element side on which the quadrature points are located + * - The quadrature point index, e.g. if there are \p n quadrature points, we are requesting the\n + * evaluation of the i-th point + * - The quadrature rule that can be used to initialize the functor on the given element and side + */ +struct ElemSideQpArg +{ + /// The element + const libMesh::Elem * elem; + + /// The local side index + unsigned int side; + + /// The quadrature point index + unsigned int qp; + + /// The qudrature rule + const QBase * qrule; + + /// The physical location of the quadrature point + Point point; + + /** + * @returns The conceptual physical location of this data structure + */ + Point getPoint() const { return point; } +}; + +/** + * State argument for evaluating functors. The iteration type indicates whether you want to evaluate + * a functor based on some iterate state of a transient calculation, nonlinear solve, etc. The state + * indicates which iterate of the iterate type we want to evaluate on. A state of 0 indicates + * "current", e.g. the current time or the current nonlinear iteration (which should actually be + * equivalent); a state of 1 indicates the most-recent "old" time or the most recent previous + * nonlinear iteration, etc. + */ + +struct StateArg +{ + /** + * Prevent implicit conversions from boolean to avoid users accidentally constructing a time + * argument when they meant to construct a skewness argument, etc. + */ + StateArg(bool) = delete; + + StateArg(unsigned int state_in) : state(state_in), iteration_type(SolutionIterationType::Time) {} + + StateArg(unsigned int state_in, SolutionIterationType iteration_type_in) + : state(state_in), iteration_type(iteration_type_in) + { + } + + /// The state. Zero represents the most recent state, so for any kind of iteration type, a zero + /// state represents the current state, e.g. current solution + /// One may represent the 'old' value (one before, in the iteration_type specified), and two an 'older' or two steps away state + unsigned int state; + + /// The solution iteration type, e.g. time or nonlinear + SolutionIterationType iteration_type; + +private: + StateArg() : state(0), iteration_type(SolutionIterationType::Time) {} + + friend StateArg currentState(); +}; + +inline StateArg +currentState() +{ + return {}; +} + +inline StateArg +oldState() +{ + return {(unsigned int)1}; +} + +inline StateArg +previousNonlinearState() +{ + return {(unsigned int)1, SolutionIterationType::Nonlinear}; +} + /** * A structure defining a "face" evaluation calling argument for Moose functors */ @@ -104,6 +192,9 @@ struct FaceArg /// on one side of the face. const Elem * face_side; + /// A member that can be used to define the instance in which the limiters are executed + const Moose::StateArg * state_limiter; + /** * @returns The conceptual physical location of this data structure */ @@ -169,92 +260,4 @@ struct ElemQpArg */ Point getPoint() const { return point; } }; - -/** - * Argument for requesting functor evaluation at quadrature point locations on an element side. - * Data in the argument: - * - The element - * - The element side on which the quadrature points are located - * - The quadrature point index, e.g. if there are \p n quadrature points, we are requesting the\n - * evaluation of the i-th point - * - The quadrature rule that can be used to initialize the functor on the given element and side - */ -struct ElemSideQpArg -{ - /// The element - const libMesh::Elem * elem; - - /// The local side index - unsigned int side; - - /// The quadrature point index - unsigned int qp; - - /// The qudrature rule - const QBase * qrule; - - /// The physical location of the quadrature point - Point point; - - /** - * @returns The conceptual physical location of this data structure - */ - Point getPoint() const { return point; } -}; - -/** - * State argument for evaluating functors. The iteration type indicates whether you want to evaluate - * a functor based on some iterate state of a transient calculation, nonlinear solve, etc. The state - * indicates which iterate of the iterate type we want to evaluate on. A state of 0 indicates - * "current", e.g. the current time or the current nonlinear iteration (which should actually be - * equivalent); a state of 1 indicates the most-recent "old" time or the most recent previous - * nonlinear iteration, etc. - */ - -struct StateArg -{ - /** - * Prevent implicit conversions from boolean to avoid users accidentally constructing a time - * argument when they meant to construct a skewness argument, etc. - */ - StateArg(bool) = delete; - - StateArg(unsigned int state_in) : state(state_in), iteration_type(SolutionIterationType::Time) {} - - StateArg(unsigned int state_in, SolutionIterationType iteration_type_in) - : state(state_in), iteration_type(iteration_type_in) - { - } - - /// The state. Zero represents the most recent state, so for any kind of iteration type, a zero - /// state represents the current state, e.g. current solution - /// One may represent the 'old' value (one before, in the iteration_type specified), and two an 'older' or two steps away state - unsigned int state; - - /// The solution iteration type, e.g. time or nonlinear - SolutionIterationType iteration_type; - -private: - StateArg() : state(0), iteration_type(SolutionIterationType::Time) {} - - friend StateArg currentState(); -}; - -inline StateArg -currentState() -{ - return {}; -} - -inline StateArg -oldState() -{ - return {(unsigned int)1}; -} - -inline StateArg -previousNonlinearState() -{ - return {(unsigned int)1, SolutionIterationType::Nonlinear}; -} } diff --git a/framework/include/fvbcs/FVBoundaryCondition.h b/framework/include/fvbcs/FVBoundaryCondition.h index 8875eba4e1fe..276d06ecefe6 100644 --- a/framework/include/fvbcs/FVBoundaryCondition.h +++ b/framework/include/fvbcs/FVBoundaryCondition.h @@ -91,7 +91,8 @@ class FVBoundaryCondition : public MooseObject, Moose::FaceArg singleSidedFaceArg( const FaceInfo * fi = nullptr, Moose::FV::LimiterType limiter_type = Moose::FV::LimiterType::CentralDifference, - bool correct_skewness = false) const; + bool correct_skewness = false, + const Moose::StateArg * state_limiter = nullptr) const; MooseVariableFV & _var; diff --git a/framework/include/fviks/FVInterfaceKernel.h b/framework/include/fviks/FVInterfaceKernel.h index ff6ad444a88b..670d9f26cf47 100644 --- a/framework/include/fviks/FVInterfaceKernel.h +++ b/framework/include/fviks/FVInterfaceKernel.h @@ -160,7 +160,8 @@ class FVInterfaceKernel : public MooseObject, const MooseVariableFV & variable, const FaceInfo * fi = nullptr, Moose::FV::LimiterType limiter_type = Moose::FV::LimiterType::CentralDifference, - bool correct_skewness = false) const; + bool correct_skewness = false, + const Moose::StateArg * state_limiter = nullptr) const; /// To be consistent with FE interfaces we introduce this quadrature point member. However, for FV /// calculations there should every only be one qudrature point and it should be located at the diff --git a/framework/include/fvkernels/FVFluxKernel.h b/framework/include/fvkernels/FVFluxKernel.h index 1812455b13f8..9fa815b1add0 100644 --- a/framework/include/fvkernels/FVFluxKernel.h +++ b/framework/include/fvkernels/FVFluxKernel.h @@ -118,7 +118,8 @@ class FVFluxKernel : public FVKernel, Moose::FaceArg singleSidedFaceArg( const FaceInfo * fi = nullptr, Moose::FV::LimiterType limiter_type = Moose::FV::LimiterType::CentralDifference, - bool correct_skewness = false) const; + bool correct_skewness = false, + const Moose::StateArg * state_limiter = nullptr) const; /** * Returns whether to avoid execution on a boundary diff --git a/framework/include/interfaces/FaceArgInterface.h b/framework/include/interfaces/FaceArgInterface.h index 82f9b9b51c4b..f994ad862286 100644 --- a/framework/include/interfaces/FaceArgInterface.h +++ b/framework/include/interfaces/FaceArgInterface.h @@ -44,7 +44,8 @@ class FaceArgProducerInterface : public FaceArgInterface Moose::FaceArg makeFace(const FaceInfo & fi, const Moose::FV::LimiterType limiter_type, const bool elem_is_upwind, - const bool correct_skewness = false) const; + const bool correct_skewness = false, + const Moose::StateArg * state_limiter = nullptr) const; /** * Make a functor face argument with a central differencing limiter, e.g. compose a face diff --git a/framework/include/limiters/CentralDifferenceLimiter.h b/framework/include/limiters/CentralDifferenceLimiter.h index 9b1c878bd58e..abfa9a66e34b 100644 --- a/framework/include/limiters/CentralDifferenceLimiter.h +++ b/framework/include/limiters/CentralDifferenceLimiter.h @@ -25,8 +25,15 @@ template class CentralDifferenceLimiter : public Limiter { public: - T - limit(const T &, const T &, const VectorValue *, const RealVectorValue &) const override final + T limit(const T &, + const T &, + const VectorValue *, + const VectorValue *, + const RealVectorValue &, + const Real &, + const Real &, + const FaceInfo *, + const bool &) const override final { return 1; } diff --git a/framework/include/limiters/Limiter.h b/framework/include/limiters/Limiter.h index 5b2f35ae8cf0..ae1b3450e796 100644 --- a/framework/include/limiters/Limiter.h +++ b/framework/include/limiters/Limiter.h @@ -12,6 +12,7 @@ #include "ADReal.h" #include "MooseTypes.h" #include "HasMembers.h" +#include "FaceInfo.h" #include class MooseEnum; @@ -29,7 +30,8 @@ enum class LimiterType : int CentralDifference, MinMod, SOU, - QUICK + QUICK, + Venkatakrishnan }; extern const MooseEnum moose_limiter_type; @@ -61,23 +63,147 @@ class Limiter { public: /** - * Defines the slope limiter function $\beta(r_f)$ where $r_f$ represents the ratio of upstream to - * downstream gradients + * This method computes the flux limiting ratio based on the provided scalar values, + * gradient vectors, direction vector, maximum and minimum allowable values, and + * geometric information from the face and cell centroids. It must be overridden + * by any derived class implementing a specific limiting strategy. + * + * @tparam T The data type of the field values and the return type. + * @param phi_upwind The field value at the upwind location. + * @param phi_downwind The field value at the downwind location. + * @param grad_phi_upwind Pointer to the gradient of the field at the upwind location. + * @param grad_phi_downwind Pointer to the gradient of the field at the downwind location. + * @param dCD A constant direction vector representing the direction of the cell face. + * @param max_value The maximum allowable value. + * @param min_value The minimum allowable value. + * @param fi FaceInfo object containing geometric details such as face centroid and cell + * centroids. + * @param fi_elem_is_upwind Boolean indicating if the face info element is upwind. + * @return The computed flux limiting ratio. + * + * This pure virtual function is intended to be defined in derived classes, which will implement + * the specific limiting algorithm. Derived classes will provide the actual logic for computing + * the flux limiting ratio, ensuring that the result adheres to the constraints and properties + * required by the specific limiting method. */ virtual T limit(const T & phi_upwind, const T & phi_downwind, const VectorValue * grad_phi_upwind, - const RealVectorValue & dCD) const = 0; + const VectorValue * grad_phi_downwind, + const RealVectorValue & dCD, + const Real & max_value, + const Real & min_value, + const FaceInfo * fi, + const bool & fi_elem_is_upwind) const = 0; virtual bool constant() const = 0; virtual InterpMethod interpMethod() const = 0; + + /** + * @tparam T The data type of the field values and the return type. + * @param phi_upwind The field value at the upwind location. + * @param phi_downwind The field value at the downwind location. + * @param grad_phi_upwind Pointer to the gradient of the field value at the upwind location. + * @param dCD A constant direction vector representing the direction of the cell face. + * @return The computed limited value, ensuring it is within the range [0, 2]. + */ T operator()(const T & phi_upwind, const T & phi_downwind, const VectorValue * grad_phi_upwind, const RealVectorValue & dCD) const { - return std::max(T(0), std::min(T(2), limit(phi_upwind, phi_downwind, grad_phi_upwind, dCD))); + return std::max(T(0), + std::min(T(2), + limit(phi_upwind, + phi_downwind, + grad_phi_upwind, + nullptr, + dCD, + T(0), + T(0), + nullptr, + false))); + } + + /** + * @tparam T The data type of the field values and the return type. + * @param phi_upwind The field value at the upwind location. + * @param phi_downwind The field value at the downwind location. + * @param grad_phi_upwind Pointer to the gradient at the upwind location. + * @param grad_phi_downwind Pointer to the gradient at the downwind location. + * @param dCD A constant direction vector representing the direction of the cell face. + * @param max_value The maximum allowable value. + * @param min_value The minimum allowable value. + * @param fi FaceInfo object containing geometric details such as face centroid and cell + * centroids. + * @param fi_elem_is_upwind Boolean indicating if the face info element is upwind. + * @return The result of the `limit` function applied to the provided parameters. + */ + T operator()(const T & phi_upwind, + const T & phi_downwind, + const VectorValue * grad_phi_upwind, + const VectorValue * grad_phi_downwind, + const RealVectorValue & dCD, + const Real & max_value, + const Real & min_value, + const FaceInfo * fi, + const bool & fi_elem_is_upwind) const + { + return limit(phi_upwind, + phi_downwind, + grad_phi_upwind, + grad_phi_downwind, + dCD, + max_value, + min_value, + fi, + fi_elem_is_upwind); } + /** + * @tparam T The data type of the gradient values and the return type. + * @param grad_phi_upwind Pointer to the gradient at the upwind location. + * @param grad_phi_downwind Pointer to the gradient at the downwind location. + * @param dCD A constant direction vector representing the direction of the cell face. + * @return The computed flux limiting ratio. + */ + T rf_grad(const VectorValue * grad_phi_upwind, + const VectorValue * grad_phi_downwind, + const RealVectorValue & dCD) const + { + const auto grad_elem = (*grad_phi_upwind) * dCD; + const auto grad_face = (*grad_phi_downwind) * dCD; + const auto grad_ratio = grad_elem / (grad_face + 1e-10); + return std::max(2.0 * grad_ratio - 1.0, 0.0); + }; + + /** + * @tparam T The data type of the field values and the return type. + * @param phi_upwind The field value at the upwind location. + * @param grad_phi_upwind Pointer to the gradient at the upwind location. + * @param max_value The maximum allowable value. + * @param min_value The minimum allowable value. + * @param fi FaceInfo object containing geometric details such as face centroid and cell + * centroids. + * @param fi_elem_is_upwind Boolean indicating if the face info element is upwind. + * @return The computed flux limiting ratio. + */ + T rf_minmax(const T & phi_upwind, + const VectorValue * grad_phi_upwind, + const Real & max_value, + const Real & min_value, + const FaceInfo * fi, + const bool & fi_elem_is_upwind) const + { + const auto face_centroid = fi->faceCentroid(); + const auto cell_centroid = fi_elem_is_upwind ? fi->elemCentroid() : fi->neighborCentroid(); + + const auto delta_face = (*grad_phi_upwind) * (face_centroid - cell_centroid); + const auto delta_max = max_value - phi_upwind + 1e-10; + const auto delta_min = min_value - phi_upwind + 1e-10; + + return delta_face >= 0 ? delta_face / delta_max : delta_face / delta_min; + }; + Limiter() = default; virtual ~Limiter() = default; diff --git a/framework/include/limiters/MinModLimiter.h b/framework/include/limiters/MinModLimiter.h index dfffa029c2a3..a70f3a6814de 100644 --- a/framework/include/limiters/MinModLimiter.h +++ b/framework/include/limiters/MinModLimiter.h @@ -17,25 +17,54 @@ namespace Moose namespace FV { /** - * Implements the Min-Mod limiter, defined by - * $\beta(r_f) = \text{max}(0, \text{min}(1, r_f))$ + * The Min-Mod limiter function $\beta(r_f)$ is defined as: + * + * \f[ + * \beta(r_f) = \text{max}(0, \text{min}(1, r_f)) + * \f] + * + * where \( r_f \) is the ratio of the upwind to downwind gradients. + * + * @tparam T The data type of the scalar values and the return type. */ template class MinModLimiter : public Limiter { public: + /** + * This method overrides the pure virtual `limit` method in the base `Limiter` class. + * It calculates the flux limiting ratio based on the Min-Mod limiter formula. + * + * @param grad_phi_upwind Pointer to the gradient vector at the upwind location. + * @param grad_phi_downwind Pointer to the gradient vector at the downwind location. + * @param dCD A constant direction vector representing the direction of the cell face. + * @return The computed flux limiting ratio. + */ T limit(const T & phi_upwind, const T & phi_downwind, const VectorValue * grad_phi_upwind, - const RealVectorValue & dCD) const override final + const VectorValue * grad_phi_downwind, + const RealVectorValue & dCD, + const Real & /* max_value */, + const Real & /* min_value */, + const FaceInfo * /* fi */, + const bool & /* fi_elem_is_upwind */) const override final { mooseAssert(grad_phi_upwind, "min-mod limiter requires a gradient"); - const auto r_f = Moose::FV::rF(phi_upwind, phi_downwind, *grad_phi_upwind, dCD); - // Dummy addition to avoid new nonzeros + // Compute gradient ratio coefficient + T r_f; + if (grad_phi_downwind) // compute full slope-reconstruction limiter + r_f = this->rf_grad(grad_phi_upwind, grad_phi_downwind, dCD); + else // compute upwind slope reconstruction limiter + r_f = Moose::FV::rF(phi_upwind, phi_downwind, *grad_phi_upwind, dCD); + + // Return limiter value return 0 * r_f + std::max(T(0), std::min(T(1), r_f)); } + bool constant() const override final { return false; } + InterpMethod interpMethod() const override final { return InterpMethod::MinMod; } MinModLimiter() = default; diff --git a/framework/include/limiters/QUICKLimiter.h b/framework/include/limiters/QUICKLimiter.h index 71171929dcd5..38355ddda121 100644 --- a/framework/include/limiters/QUICKLimiter.h +++ b/framework/include/limiters/QUICKLimiter.h @@ -17,24 +17,74 @@ namespace Moose namespace FV { /** - * Implements a limiter which reproduces the QUICK scheme, defined by - * $\beta(r_f) = \frac{3+r_f}{4}$ + * The QUICK (Quadratic Upstream Interpolation for Convective Kinematics) limiter + * function is derived from the following equations: + * + * 1. Calculation of the gradient ratio coefficient \( r_f \): + * \f[ + * r_f = \begin{cases} + * \text{if grad\_phi\_downwind is not null, use } \text{rf\_grad}(\nabla \phi_{\text{upwind}}, + * \nabla \phi_{\text{downwind}}, \mathbf{dCD}) \\ \text{otherwise, use } + * \text{rF}(\phi_{\text{upwind}}, \phi_{\text{downwind}}, \nabla \phi_{\text{upwind}}, + * \mathbf{dCD}) \end{cases} \f] + * + * 2. QUICK limiter formula ensuring TVD compliance: + * \f[ + * \beta(r_f) = \min\left(\beta, \max\left(\min\left(\min\left(\frac{1 + 3.0 r_f}{4.0}, 2.0 + * r_f\right), 2.0\right), 0.0\right)\right) \f] where \(\beta = 1.0\). + * + * @tparam T The data type of the scalar values and the return type. */ template class QUICKLimiter : public Limiter { public: + /** + * This method overrides the pure virtual `limit` method in the base `Limiter` class. + * It calculates the flux limiting ratio based on the QUICK limiter formula. + * + * @param phi_upwind Scalar value at the upwind location. + * @param phi_downwind Scalar value at the downwind location. + * @param grad_phi_upwind Pointer to the gradient vector at the upwind location. + * @param grad_phi_downwind Pointer to the gradient vector at the downwind location. + * @param dCD A constant direction vector representing the direction of the cell face. + * @return The computed flux limiting ratio. + */ T limit(const T & phi_upwind, const T & phi_downwind, const VectorValue * grad_phi_upwind, - const RealVectorValue & dCD) const override final + const VectorValue * grad_phi_downwind, + const RealVectorValue & dCD, + const Real & /* max_value */, + const Real & /* min_value */, + const FaceInfo * /* fi */, + const bool & /* fi_elem_is_upwind */) const override final { mooseAssert(grad_phi_upwind, "QUICK limiter requires a gradient"); - const auto r_f = Moose::FV::rF(phi_upwind, phi_downwind, *grad_phi_upwind, dCD); - return (3. + r_f) / 4.; + // Compute gradient ratio coefficient + T limiter; + if (grad_phi_downwind) // compute full slope-reconstruction limiter for weakly-compressible flow + { + const auto & r_f = this->rf_grad(grad_phi_upwind, grad_phi_downwind, dCD); + const auto & beta = T(1.0); // Ensures TVD compliance + limiter = + 0 * r_f + + std::min(beta, + std::max(std::min(std::min((1 + 3.0 * r_f) / 4.0, 2.0 * r_f), T(2.0)), T(0.0))); + } + else // compute upwind slope reconstruction limiter for compressible flow + { + const auto & r_f = Moose::FV::rF(phi_upwind, phi_downwind, *grad_phi_upwind, dCD); + limiter = (3.0 * r_f) / 4.0; + } + + // Return limiter value + return limiter; } + bool constant() const override final { return false; } + InterpMethod interpMethod() const override final { return InterpMethod::QUICK; } QUICKLimiter() = default; diff --git a/framework/include/limiters/SOULimiter.h b/framework/include/limiters/SOULimiter.h index aaaf08131cca..9192dd354b08 100644 --- a/framework/include/limiters/SOULimiter.h +++ b/framework/include/limiters/SOULimiter.h @@ -17,23 +17,59 @@ namespace Moose namespace FV { /** - * Implements a limiter which reproduces the second-order-upwind scheme, defined by - * $\beta(r_f) = r_f$ + * The SOU limiter is used for reproducing the second-order-upwind scheme. The limiter function + * $\beta(delta_max)$ is defined as: + * + * \f[ + * \beta(delta_max) = min(1, delta_face/delta_max) + * \f] + * + * where: + * \( delta_max \) is the maximum variation admited for the computatioanl stencil + * \( delta_max \) is the variation at the face computed with SOU + * + * @tparam T The data type of the scalar values and the return type. */ template class SOULimiter : public Limiter { public: - T limit(const T & phi_upwind, - const T & phi_downwind, + /** + * This method overrides the pure virtual `limit` method in the base `Limiter` class. + * It calculates the flux limiting ratio based on the SOU limiter formula. + * + * @param grad_phi_upwind Pointer to the gradient vector at the upwind location. + * @param grad_phi_downwind Pointer to the gradient vector at the downwind location. + * @param dCD A constant direction vector representing the direction of the cell face. + * @param max_value The maximum value in the current element. + * @param min_value The minimum value in the current element. + * @param fi Pointer to the face information structure. + * @param fi_elem_is_upwind Boolean flag indicating if the current element is upwind. + * @return The computed flux limiting ratio. + */ + T limit(const T & /* phi_upwind */, + const T & /* phi_downwind */, const VectorValue * grad_phi_upwind, - const RealVectorValue & dCD) const override final + const VectorValue * /* grad_phi_downwind */, + const RealVectorValue & dCD, + const Real & max_value, + const Real & min_value, + const FaceInfo * /* fi */, + const bool & /*fi_elem_is_upwind */) const override final { mooseAssert(grad_phi_upwind, "SOU limiter requires a gradient"); - const auto r_f = Moose::FV::rF(phi_upwind, phi_downwind, *grad_phi_upwind, dCD); - return r_f; + + T delta_face = std::abs((*grad_phi_upwind) * dCD); + T delta_max = std::abs(max_value - min_value); + + if (delta_face > 1e-10) + return std::min(1.0, delta_max / delta_face); + else + return T(1.0); } + bool constant() const override final { return false; } + InterpMethod interpMethod() const override final { return InterpMethod::SOU; } SOULimiter() = default; diff --git a/framework/include/limiters/UpwindLimiter.h b/framework/include/limiters/UpwindLimiter.h index 62e9148a453f..43734b7d6928 100644 --- a/framework/include/limiters/UpwindLimiter.h +++ b/framework/include/limiters/UpwindLimiter.h @@ -24,8 +24,15 @@ template class UpwindLimiter : public Limiter { public: - T - limit(const T &, const T &, const VectorValue *, const RealVectorValue &) const override final + T limit(const T &, + const T &, + const VectorValue *, + const VectorValue *, + const RealVectorValue &, + const Real &, + const Real &, + const FaceInfo *, + const bool &) const override final { return 0; } diff --git a/framework/include/limiters/VanLeerLimiter.h b/framework/include/limiters/VanLeerLimiter.h index f07e8e5a072e..8ed61b409fb0 100644 --- a/framework/include/limiters/VanLeerLimiter.h +++ b/framework/include/limiters/VanLeerLimiter.h @@ -17,23 +17,54 @@ namespace Moose namespace FV { /** - * Implements the Van Leer limiter, defined by - * $\beta(r_f) = \frac{r_f + \text{abs}(r_f)}{1 + \text{abs}(r_f)}$ + * The Van Leer limiter limiter function $\beta(r_f)$ is defined as: + * + * \f[ + * \beta(r_f) = \frac{r_f + \text{abs}(r_f)}{1 + \text{abs}(r_f)} + * \f] + * + * where \( r_f \) is the ratio of the upwind to downwind gradients. + * + * @tparam T The data type of the scalar values and the return type. */ template class VanLeerLimiter : public Limiter { public: + /** + * This method overrides the pure virtual `limit` method in the base `Limiter` class. + * It calculates the flux limiting ratio based on the Van Leer limiter formula. + * + * @param grad_phi_upwind Pointer to the gradient vector at the upwind location. + * @param grad_phi_downwind Pointer to the gradient vector at the downwind location. + * @param dCD A constant direction vector representing the direction of the cell face. + * @return The computed flux limiting ratio. + */ T limit(const T & phi_upwind, const T & phi_downwind, const VectorValue * grad_phi_upwind, - const RealVectorValue & dCD) const override final + const VectorValue * grad_phi_downwind, + const RealVectorValue & dCD, + const Real & /* max_value */, + const Real & /* min_value */, + const FaceInfo * /* fi */, + const bool & /* fi_elem_is_upwind */) const override final { - mooseAssert(grad_phi_upwind, "Van Leer limiter requires a gradient"); - const auto r_f = Moose::FV::rF(phi_upwind, phi_downwind, *grad_phi_upwind, dCD); - return (r_f + std::abs(r_f)) / (1. + std::abs(r_f)); + mooseAssert(grad_phi_upwind, "Van Leer limiter requires the upwind gradient"); + + // Compute gradient ratio coefficient + T r_f; + if (grad_phi_downwind) // compute full slope-reconstruction limiter + r_f = this->rf_grad(grad_phi_upwind, grad_phi_downwind, dCD); + else // compute upwind slope reconstruction limiter + r_f = Moose::FV::rF(phi_upwind, phi_downwind, *grad_phi_upwind, dCD); + + // Return limiter value + return (r_f + std::abs(r_f)) / (1.0 + std::abs(r_f)); } + bool constant() const override final { return false; } + InterpMethod interpMethod() const override final { return InterpMethod::VanLeer; } VanLeerLimiter() = default; diff --git a/framework/include/limiters/VenkatakrishnanLimiter.h b/framework/include/limiters/VenkatakrishnanLimiter.h new file mode 100644 index 000000000000..7688bd8ae85c --- /dev/null +++ b/framework/include/limiters/VenkatakrishnanLimiter.h @@ -0,0 +1,93 @@ +//* 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 "Limiter.h" +#include "MathFVUtils.h" + +namespace Moose +{ +namespace FV +{ +/** + * The Venkatakrishnan limiter is derived from the following equations: + * + * 1. Calculation of the face delta: + * \f[ + * \Delta_{\text{face}} = \nabla \phi_{\text{upwind}} \cdot (\mathbf{x}_{\text{face}} - + * \mathbf{x}_{\text{cell}}) \f] + * + * 2. Calculation of the deltas for maximum and minimum values relative to the upwind value with a + * small perturbation to avoid division by zero: \f[ \Delta_{\text{max}} = \phi_{\text{max}} - + * \phi_{\text{upwind}} + 1e-10 \f] \f[ \Delta_{\text{min}} = \phi_{\text{min}} - + * \phi_{\text{upwind}} + 1e-10 \f] + * + * 3. Calculation of the gradient ratio coefficient \( r_f \): + * \f[ + * r_f = \begin{cases} + * \frac{\Delta_{\text{face}}}{\Delta_{\text{max}}} & \text{if } \Delta_{\text{face}} \geq 0 \\ + * \frac{\Delta_{\text{face}}}{\Delta_{\text{min}}} & \text{otherwise} + * \end{cases} + * \f] + * + * 4. Venkatakrishnan limiter formula (Venkatakrishnan, 1993): + * \f[ + * \beta(r_f) = \frac{2 r_f + 1.0}{r_f (2 r_f + 1.0) + 1.0} + * \f] + * + * @tparam T The data type of the scalar values and the return type. + */ +template +class VenkatakrishnanLimiter : public Limiter +{ +public: + /** + * This method overrides the pure virtual `limit` method in the base `Limiter` class. + * It calculates the flux limiting ratio based on the Venkatakrishnan limiter formula. + * + * @param phi_upwind Scalar value at the upwind location. + * @param grad_phi_upwind Pointer to the gradient vector at the upwind location. + * @param max_value The maximum value in the current element. + * @param min_value The minimum value in the current element. + * @param fi Pointer to the face information structure. + * @param fi_elem_is_upwind Boolean flag indicating if the current element is upwind. + * @return The computed flux limiting ratio. + */ + T limit(const T & phi_upwind, + const T & /* phi_downwind */, + const VectorValue * grad_phi_upwind, + const VectorValue * /* grad_phi_downwind */, + const RealVectorValue & /* dCD */, + const Real & max_value, + const Real & min_value, + const FaceInfo * fi, + const bool & fi_elem_is_upwind) const override final + { + const auto face_centroid = fi->faceCentroid(); + const auto cell_centroid = fi_elem_is_upwind ? fi->elemCentroid() : fi->neighborCentroid(); + + const auto delta_face = (*grad_phi_upwind) * (face_centroid - cell_centroid); + const auto delta_max = std::abs(max_value - phi_upwind) + 1e-10; + const auto delta_min = std::abs(min_value - phi_upwind) + 1e-10; + + const auto rf = + delta_face >= 0 ? std::abs(delta_face) / delta_max : std::abs(delta_face) / delta_min; + + return (2 * rf + 1.0) / (rf * (2 * rf + 1.0) + 1.0); + } + + bool constant() const override final { return false; } + + InterpMethod interpMethod() const override final { return InterpMethod::Venkatakrishnan; } + + VenkatakrishnanLimiter() = default; +}; +} +} diff --git a/framework/include/utils/GreenGaussGradient.h b/framework/include/utils/GreenGaussGradient.h index c8def8cc405c..704fbee3a29f 100644 --- a/framework/include/utils/GreenGaussGradient.h +++ b/framework/include/utils/GreenGaussGradient.h @@ -142,7 +142,8 @@ greenGaussGradient(const ElemArg & elem_arg, Moose::FV::LimiterType::CentralDifference, true, elem_arg.correct_skewness, - elem_arg.elem}, + elem_arg.elem, + nullptr}, state_arg); if (!volume_set) @@ -370,7 +371,8 @@ greenGaussGradient(const ElemArg & elem_arg, Moose::FV::LimiterType::CentralDifference, true, elem_arg.correct_skewness, - elem_arg.elem}, + elem_arg.elem, + nullptr}, state_arg); } diff --git a/framework/include/utils/MathFVUtils.h b/framework/include/utils/MathFVUtils.h index 7ac85fbd5ca0..3e07407219d9 100644 --- a/framework/include/utils/MathFVUtils.h +++ b/framework/include/utils/MathFVUtils.h @@ -47,7 +47,8 @@ enum class InterpMethod VanLeer, MinMod, SOU, - QUICK + QUICK, + Venkatakrishnan }; enum class LinearFVComputationMode @@ -472,12 +473,22 @@ interpCoeffs(const Limiter & limiter, const T & phi_upwind, const T & phi_downwind, const VectorValue * const grad_phi_upwind, + const VectorValue * const grad_phi_face, + const Real & max_value, + const Real & min_value, const FaceInfo & fi, const bool fi_elem_is_upwind) { // Using beta, w_f, g nomenclature from Greenshields - const auto beta = limiter( - phi_upwind, phi_downwind, grad_phi_upwind, fi_elem_is_upwind ? fi.dCN() : Point(-fi.dCN())); + const auto beta = limiter(phi_upwind, + phi_downwind, + grad_phi_upwind, + grad_phi_face, + fi_elem_is_upwind ? fi.dCN() : Point(-fi.dCN()), + max_value, + min_value, + &fi, + fi_elem_is_upwind); const auto w_f = fi_elem_is_upwind ? fi.gC() : (1. - fi.gC()); @@ -500,7 +511,16 @@ interpolate(const Limiter & limiter, const FaceInfo & fi, const bool fi_elem_is_upwind) { - auto pr = interpCoeffs(limiter, phi_upwind, phi_downwind, grad_phi_upwind, fi, fi_elem_is_upwind); + auto pr = + interpCoeffs(limiter, + phi_upwind, + phi_downwind, + grad_phi_upwind, + /*grad_phi_face*/ static_cast *>(nullptr), + /* max_value */ std::numeric_limits::max(), + /* min_value */ std::numeric_limits::min(), + fi, + fi_elem_is_upwind); return pr.first * phi_upwind + pr.second * phi_downwind; } @@ -538,11 +558,225 @@ interpolate(const Limiter & limiter, } /** - * Interpolates with a limiter and a face argument - * @return a pair of pairs. The first pair corresponds to the interpolation coefficients with the - * first corresponding to the face information element and the second corresponding to the face - * information neighbor. This pair should sum to unity. The second pair corresponds to the face - * information functor element value and neighbor + * This function performs a full limited interpolation of a variable, taking into account + * the values and gradients at both upwind and downwind locations, as well as geometric + * information and limits. It applies the specified limiter to ensure the interpolation + * adheres to the constraints imposed by the limiter. + * + * @tparam T The data type of the scalar values and the return type. + * @param limiter The limiter object used to compute the flux limiting ratio. + * @param phi_upwind The field value at the upwind location. + * @param phi_downwind The field value at the downwind location. + * @param grad_phi_upwind Pointer to the gradient at the upwind location. + * @param grad_phi_face Pointer to the gradient at the face location. + * @param fi FaceInfo object containing geometric details such as face centroid and cell centroids. + * @param fi_elem_is_upwind Boolean indicating if the face info element is upwind. + * @param max_value The maximum allowable value. + * @param min_value The minimum allowable value. + * @return The computed limited interpolation value. + */ +template +T +fullLimitedInterpolation(const Limiter & limiter, + const T & phi_upwind, + const T & phi_downwind, + const VectorValue * const grad_phi_upwind, + const VectorValue * const grad_phi_face, + const Real & max_value, + const Real & min_value, + const FaceArg & face) +{ + // Compute the direction vector based on whether the current element is upwind + const auto dCD = face.elem_is_upwind ? face.fi->dCN() : Point(-face.fi->dCN()); + + // Compute the flux limiting ratio using the specified limiter + const auto beta = limiter(phi_upwind, + phi_downwind, + grad_phi_upwind, + grad_phi_face, + dCD, + max_value, + min_value, + face.fi, + face.elem_is_upwind); + + // Determine the face centroid and the appropriate cell centroid + const auto & face_centroid = face.fi->faceCentroid(); + const auto & cell_centroid = + face.elem_is_upwind ? face.fi->elemCentroid() : face.fi->neighborCentroid(); + + // Compute the delta value at the face + T delta_face; + if (face.limiter_type == LimiterType::Venkatakrishnan || face.limiter_type == LimiterType::SOU) + delta_face = (*grad_phi_upwind) * (face_centroid - cell_centroid); + else + delta_face = (*grad_phi_face) * (face_centroid - cell_centroid); + + // Return the interpolated value + return phi_upwind + beta * delta_face; +} + +/** + * This function calculates the minimum and maximum values within a two-cell stencil. The stencil + * includes the immediate neighboring elements of the face's associated element and the neighboring + * elements of those neighbors. It evaluates the values using a provided functor and accounts for + * the valid (non-null) neighbors. + * + * @tparam T The data type for the values being computed. This is typically a scalar type. + * @tparam FEK Enumeration type FunctorEvaluationKind with a default value of + * FunctorEvaluationKind::Value. This determines the kind of evaluation the functor will perform. + * @tparam Enable A type trait used for SFINAE (Substitution Failure Is Not An Error) to ensure that + * T is a valid scalar type as determined by ScalarTraits. + * + * @param functor An object of a functor class derived from FunctorBase. This object provides the + * genericEvaluate method to compute the value at a given element and time. + * @param face An argument representing the face in the computational domain. The face provides + * access to neighboring elements via neighbor_ptr_range(). + * @param time An argument representing the state or time at which the evaluation is performed. + * + * @return std::pair A pair containing the minimum and maximum values computed across the + * two-cell stencil. The first element is the maximum value, and the second element is the minimum + * value. + */ +template ::value>::type> +std::pair +computeMinMaxValue(const FunctorBase & functor, const FaceArg & face, const StateArg & time) +{ + // Initialize max_value to 0 and min_value to a large value + Real max_value(std::numeric_limits::min()), min_value(std::numeric_limits::max()); + + // Iterate over the direct neighbors of the element associated with the face + for (const auto neighbor : (*face.fi).elem().neighbor_ptr_range()) + { + // If not a valid neighbor, skip to the next one + if (neighbor == nullptr) + continue; + + // Evaluate the functor at the neighbor and update max_value and min_value + const ElemArg local_cell_arg = {neighbor, false}; + const auto value = + MetaPhysicL::raw_value(functor.template genericEvaluate(local_cell_arg, time)); + max_value = std::max(max_value, value); + min_value = std::min(min_value, value); + } + + // Iterate over the neighbors of the neighbor + for (const auto neighbor : (*face.fi).neighbor().neighbor_ptr_range()) + { + // If not a valid neighbor, skip to the next one + if (neighbor == nullptr) + continue; + + // Evaluate the functor at the neighbor and update max_value and min_value + const ElemArg local_cell_arg = {neighbor, false}; + const auto value = + MetaPhysicL::raw_value(functor.template genericEvaluate(local_cell_arg, time)); + max_value = std::max(max_value, value); + min_value = std::min(min_value, value); + } + + // Return a pair containing the computed maximum and minimum values + return std::make_pair(max_value, min_value); +} + +/** + * This function calculates the minimum and maximum values of a specified component within a + * two-cell stencil. The stencil includes the immediate neighboring elements of the face's + * associated element and the neighboring elements of those neighbors. It evaluates the values using + * a provided functor and accounts for the valid (non-null) neighbors. + * + * @tparam T The data type for the values being computed. This is typically a scalar type. + * + * @param functor An object of a functor class derived from FunctorBase>. This object + * provides the operator() method to compute the value at a given element and time. + * @param face An argument representing the face in the computational domain. The face provides + * access to neighboring elements via neighbor_ptr_range(). + * @param time An argument representing the state or time at which the evaluation is performed. + * @param component An unsigned integer representing the specific component of the vector to be + * evaluated. + * + * @return std::pair A pair containing the minimum and maximum values of the specified + * component computed across the two-cell stencil. The first element is the maximum value, and the + * second element is the minimum value. + * + * Usage: + * This function is typically used in the finite volume methods for min-max computations over + * stencils (neighborhoods). It helps compute the limiting for limited second order upwind at the + * faces. + */ +template +std::pair +computeMinMaxValue(const FunctorBase> & functor, + const FaceArg & face, + const StateArg & time, + const unsigned int & component) +{ + // Initialize max_value to 0 and min_value to a large value + Real max_value(std::numeric_limits::min()), min_value(std::numeric_limits::max()); + + // Iterate over the direct neighbors of the element associated with the face + for (const auto neighbor : (*face.fi).elem().neighbor_ptr_range()) + { + // If not a valid neighbor, skip to the next one + if (neighbor == nullptr) + continue; + + // Evaluate the functor at the neighbor for the specified component and update max_value and + // min_value + const ElemArg local_cell_arg = {neighbor, false}; + const auto value = MetaPhysicL::raw_value(functor(local_cell_arg, time)(component)); + max_value = std::max(max_value, value); + min_value = std::min(min_value, value); + } + + // Iterate over the neighbors of the neighbor associated with the face + for (const auto neighbor : (*face.fi).neighbor().neighbor_ptr_range()) + { + // If not a valid neighbor, skip to the next one + if (neighbor == nullptr) + continue; + + // Evaluate the functor at the neighbor for the specified component and update max_value and + // min_value + const ElemArg local_cell_arg = {neighbor, false}; + const auto value = MetaPhysicL::raw_value(functor(local_cell_arg, time)(component)); + max_value = std::max(max_value, value); + min_value = std::min(min_value, value); + } + + // Return a pair containing the computed maximum and minimum values + return std::make_pair(max_value, min_value); +} + +/** + * This function interpolates values using a specified limiter and face argument. It evaluates the + * values at upwind and downwind locations and computes interpolation coefficients and advected + * values. + * + * @tparam T The data type for the values being interpolated. This is typically a scalar type. + * @tparam FEK Enumeration type FunctorEvaluationKind with a default value of + * FunctorEvaluationKind::Value. This determines the kind of evaluation the functor will perform. + * @tparam Enable A type trait used for SFINAE (Substitution Failure Is Not An Error) to ensure that + * T is a valid scalar type as determined by ScalarTraits. + * + * @param functor An object of a functor class derived from FunctorBase. This object provides the + * genericEvaluate method to compute the value at a given element and time. + * @param face An argument representing the face in the computational domain. The face provides + * access to neighboring elements and limiter type. + * @param time An argument representing the state or time at which the evaluation is performed. + * + * @return std::pair, std::pair> A pair of pairs. + * - The first pair corresponds to the interpolation coefficients, with the + * first value corresponding to the face information element and the second to the face information + * neighbor. This pair should sum to unity. + * - The second pair corresponds to the face information functor element + * value and neighbor value. + * + * Usage: + * This function is used for interpolating values at faces in a finite volume method, ensuring that + * the interpolation adheres to the constraints imposed by the limiter. */ template , std::pair> interpCoeffsAndAdvected(const FunctorBase & functor, const FaceArg & face, const StateArg & time) { + // Ensure only supported FunctorEvaluationKinds are used static_assert((FEK == FunctorEvaluationKind::Value) || (FEK == FunctorEvaluationKind::Dot), "Only Value and Dot evaluations currently supported"); + // Determine the gradient evaluation kind constexpr FunctorEvaluationKind GFEK = FunctorGradientEvaluationKind::value; typedef typename FunctorBase::GradientType GradientType; static const GradientType zero(0); mooseAssert(face.fi, "this must be non-null"); + + // Construct the limiter based on the face limiter type const auto limiter = Limiter::value_type>::build(face.limiter_type); + // Determine upwind and downwind arguments based on the face element const auto upwind_arg = face.elem_is_upwind ? face.makeElem() : face.makeNeighbor(); const auto downwind_arg = face.elem_is_upwind ? face.makeNeighbor() : face.makeElem(); + + // Evaluate the functor at the upwind and downwind locations auto phi_upwind = functor.template genericEvaluate(upwind_arg, time); auto phi_downwind = functor.template genericEvaluate(downwind_arg, time); + // Initialize the interpolation coefficients std::pair interp_coeffs; + + // Compute interpolation coefficients for Upwind or CentralDifference limiters if (face.limiter_type == LimiterType::Upwind || face.limiter_type == LimiterType::CentralDifference) - interp_coeffs = - interpCoeffs(*limiter, phi_upwind, phi_downwind, &zero, *face.fi, face.elem_is_upwind); + interp_coeffs = interpCoeffs(*limiter, + phi_upwind, + phi_downwind, + &zero, + &zero, + std::numeric_limits::max(), + std::numeric_limits::min(), + *face.fi, + face.elem_is_upwind); else { - const auto grad_phi_upwind = functor.template genericEvaluate(upwind_arg, time); - interp_coeffs = interpCoeffs( - *limiter, phi_upwind, phi_downwind, &grad_phi_upwind, *face.fi, face.elem_is_upwind); + // Determine the time argument for the limiter + auto * time_arg = face.state_limiter; + if (!time_arg) + { + static Moose::StateArg temp_state(0, Moose::SolutionIterationType::Time); + time_arg = &temp_state; + } + + Real max_value(std::numeric_limits::min()), min_value(std::numeric_limits::max()); + + // Compute min-max values for min-max limiters + if (face.limiter_type == LimiterType::Venkatakrishnan || face.limiter_type == LimiterType::SOU) + std::tie(max_value, min_value) = computeMinMaxValue(functor, face, *time_arg); + + // Evaluate the gradient of the functor at the upwind and downwind locations + const auto grad_phi_upwind = functor.template genericEvaluate(upwind_arg, *time_arg); + const auto grad_phi_face = functor.template genericEvaluate(face, *time_arg); + + // Compute the interpolation coefficients using the specified limiter + interp_coeffs = interpCoeffs(*limiter, + phi_upwind, + phi_downwind, + &grad_phi_upwind, + &grad_phi_face, + max_value, + min_value, + *face.fi, + face.elem_is_upwind); } + // Return the interpolation coefficients and advected values if (face.elem_is_upwind) return std::make_pair(std::move(interp_coeffs), std::make_pair(std::move(phi_upwind), std::move(phi_downwind))); @@ -586,12 +863,36 @@ interpCoeffsAndAdvected(const FunctorBase & functor, const FaceArg & face, co std::make_pair(std::move(phi_downwind), std::move(phi_upwind))); } +/** + * This function interpolates values at faces in a computational grid using a specified functor, + * face argument, and evaluation kind. It handles different limiter types and performs + * interpolation accordingly. + * + * @tparam T The data type for the values being interpolated. This is typically a scalar type. + * @tparam FEK Enumeration type FunctorEvaluationKind with a default value of + * FunctorEvaluationKind::Value. This determines the kind of evaluation the functor will perform. + * @tparam Enable A type trait used for SFINAE (Substitution Failure Is Not An Error) to ensure that + * T is a valid scalar type as determined by ScalarTraits. + * + * @param functor An object of a functor class derived from FunctorBase. This object provides the + * genericEvaluate method to compute the value at a given element and time. + * @param face An argument representing the face in the computational domain. The face provides + * access to neighboring elements and limiter type. + * @param time An argument representing the state or time at which the evaluation is performed. + * + * @return T The interpolated value at the face. + * + * Usage: + * This function is used for interpolating values at faces in a finite volume method, ensuring that + * the interpolation adheres to the constraints imposed by the limiter. + */ template ::value>::type> T interpolate(const FunctorBase & functor, const FaceArg & face, const StateArg & time) { + // Ensure only supported FunctorEvaluationKinds are used static_assert((FEK == FunctorEvaluationKind::Value) || (FEK == FunctorEvaluationKind::Dot), "Only Value and Dot evaluations currently supported"); @@ -600,29 +901,103 @@ interpolate(const FunctorBase & functor, const FaceArg & face, const StateArg if (face.limiter_type == LimiterType::CentralDifference) return linearInterpolation(functor, face, time); - const auto [interp_coeffs, advected] = interpCoeffsAndAdvected(functor, face, time); - return interp_coeffs.first * advected.first + interp_coeffs.second * advected.second; + if (face.limiter_type == LimiterType::Upwind || + face.limiter_type == LimiterType::CentralDifference) + { + // Compute interpolation coefficients and advected values + const auto [interp_coeffs, advected] = interpCoeffsAndAdvected(functor, face, time); + // Return the interpolated value + return interp_coeffs.first * advected.first + interp_coeffs.second * advected.second; + } + else + { + // Determine the gradient evaluation kind + constexpr FunctorEvaluationKind GFEK = FunctorGradientEvaluationKind::value; + typedef typename FunctorBase::GradientType GradientType; + static const GradientType zero(0); + const auto & limiter = + Limiter::value_type>::build(face.limiter_type); + + // Determine upwind and downwind arguments based on the face element + const auto & upwind_arg = face.elem_is_upwind ? face.makeElem() : face.makeNeighbor(); + const auto & downwind_arg = face.elem_is_upwind ? face.makeNeighbor() : face.makeElem(); + const auto & phi_upwind = functor.template genericEvaluate(upwind_arg, time); + const auto & phi_downwind = functor.template genericEvaluate(downwind_arg, time); + + // Determine the time argument for the limiter + auto * time_arg = face.state_limiter; + if (!time_arg) + { + static Moose::StateArg temp_state(0, Moose::SolutionIterationType::Time); + time_arg = &temp_state; + } + + // Initialize min-max values + Real max_value(std::numeric_limits::min()), min_value(std::numeric_limits::max()); + if (face.limiter_type == LimiterType::Venkatakrishnan || face.limiter_type == LimiterType::SOU) + std::tie(max_value, min_value) = computeMinMaxValue(functor, face, *time_arg); + + // Evaluate the gradient of the functor at the upwind and downwind locations + const auto & grad_phi_upwind = functor.template genericEvaluate(upwind_arg, *time_arg); + const auto & grad_phi_face = functor.template genericEvaluate(face, *time_arg); + + // Perform full limited interpolation and return the interpolated value + return fullLimitedInterpolation(*limiter, + phi_upwind, + phi_downwind, + &grad_phi_upwind, + &grad_phi_face, + max_value, + min_value, + face); + } } +/** + * This function interpolates vector values at faces in a computational grid using a specified + * functor, face argument, and limiter type. It handles different limiter types and performs + * interpolation accordingly. + * + * @tparam T The data type for the vector values being interpolated. This is typically a scalar + * type. + * + * @param functor An object of a functor class derived from FunctorBase>. This object + * provides the operator() method to compute the value at a given element and time. + * @param face An argument representing the face in the computational domain. The face provides + * access to neighboring elements and limiter type. + * @param time An argument representing the state or time at which the evaluation is performed. + * + * @return VectorValue The interpolated vector value at the face. + * + * Usage: + * This function is used for interpolating vector values at faces in a finite volume method, + * ensuring that the interpolation adheres to the constraints imposed by the limiter. + */ template VectorValue interpolate(const FunctorBase> & functor, const FaceArg & face, const StateArg & time) { + // Define a zero gradient vector for initialization static const VectorValue grad_zero(0); mooseAssert(face.fi, "this must be non-null"); + + // Construct the limiter based on the face limiter type const auto limiter = Limiter::value_type>::build(face.limiter_type); + // Determine upwind and downwind arguments based on the face element const auto upwind_arg = face.elem_is_upwind ? face.makeElem() : face.makeNeighbor(); const auto downwind_arg = face.elem_is_upwind ? face.makeNeighbor() : face.makeElem(); auto phi_upwind = functor(upwind_arg, time); auto phi_downwind = functor(downwind_arg, time); + // Initialize the return vector value VectorValue ret; T coeff_upwind, coeff_downwind; + // Compute interpolation coefficients and advected values for Upwind or CentralDifference limiters if (face.limiter_type == LimiterType::Upwind || face.limiter_type == LimiterType::CentralDifference) { @@ -633,6 +1008,9 @@ interpolate(const FunctorBase> & functor, component_upwind, component_downwind, &grad_zero, + &grad_zero, + std::numeric_limits::max(), + std::numeric_limits::min(), *face.fi, face.elem_is_upwind); ret(i) = coeff_upwind * component_upwind + coeff_downwind * component_downwind; @@ -640,20 +1018,70 @@ interpolate(const FunctorBase> & functor, } else { - const auto grad_phi_upwind = functor.gradient(upwind_arg, time); + // Determine the time argument for the limiter + auto * time_arg = face.state_limiter; + if (!time_arg) + { + static Moose::StateArg temp_state(0, Moose::SolutionIterationType::Time); + time_arg = &temp_state; + } + + // Evaluate the gradient of the functor at the upwind and downwind locations + const auto & grad_phi_upwind = functor.gradient(upwind_arg, *time_arg); + const auto & grad_phi_downwind = functor.gradient(downwind_arg, *time_arg); + + const auto coeffs = interpCoeffs(InterpMethod::Average, *face.fi, face.elem_is_upwind); + + // Compute interpolation coefficients and advected values for each component for (unsigned int i = 0; i < LIBMESH_DIM; ++i) { const auto &component_upwind = phi_upwind(i), component_downwind = phi_downwind(i); - const VectorValue grad = grad_phi_upwind.row(i); - std::tie(coeff_upwind, coeff_downwind) = interpCoeffs( - *limiter, component_upwind, component_downwind, &grad, *face.fi, face.elem_is_upwind); - ret(i) = coeff_upwind * component_upwind + coeff_downwind * component_downwind; + const VectorValue &grad_upwind = grad_phi_upwind.row(i), + grad_face = coeffs.first * grad_phi_upwind.row(i) + + coeffs.second * grad_phi_downwind.row(i); + + // Initialize min-max values + Real max_value(std::numeric_limits::min()), min_value(std::numeric_limits::max()); + if (face.limiter_type == LimiterType::Venkatakrishnan || + face.limiter_type == LimiterType::SOU) + std::tie(max_value, min_value) = computeMinMaxValue(functor, face, *time_arg, i); + + // Perform full limited interpolation for the component + ret(i) = fullLimitedInterpolation(*limiter, + component_upwind, + component_downwind, + &grad_upwind, + &grad_face, + max_value, + min_value, + face); } } + // Return the interpolated vector value return ret; } +/** + * This function interpolates container values at faces in a computational grid using a specified + * functor, face argument, and limiter type. It handles different limiter types and performs + * interpolation accordingly. + * + * @tparam T The data type for the container values being interpolated. This is typically a + * container type such as a vector or array. + * + * @param functor An object of a functor class derived from FunctorBase. This object provides the + * operator() method to compute the value at a given element and time. + * @param face An argument representing the face in the computational domain. The face provides + * access to neighboring elements and limiter type. + * @param time An argument representing the state or time at which the evaluation is performed. + * + * @return T The interpolated container value at the face. + * + * Usage: + * This function is used for interpolating container values at faces in a finite volume method, + * ensuring that the interpolation adheres to the constraints imposed by the limiter. + */ template T containerInterpolate(const FunctorBase & functor, const FaceArg & face, const StateArg & time) @@ -684,6 +1112,9 @@ containerInterpolate(const FunctorBase & functor, const FaceArg & face, const component_upwind, component_downwind, example_gradient, + example_gradient, + std::numeric_limits::max(), + std::numeric_limits::min(), *face.fi, face.elem_is_upwind); ret[i] = coeff_upwind * component_upwind + coeff_downwind * component_downwind; @@ -696,8 +1127,15 @@ containerInterpolate(const FunctorBase & functor, const FaceArg & face, const { const auto &component_upwind = phi_upwind[i], component_downwind = phi_downwind[i]; const auto & grad = grad_phi_upwind[i]; - std::tie(coeff_upwind, coeff_downwind) = interpCoeffs( - *limiter, component_upwind, component_downwind, &grad, *face.fi, face.elem_is_upwind); + std::tie(coeff_upwind, coeff_downwind) = interpCoeffs(*limiter, + component_upwind, + component_downwind, + &grad, + example_gradient, + std::numeric_limits::max(), + std::numeric_limits::min(), + *face.fi, + face.elem_is_upwind); ret[i] = coeff_upwind * component_upwind + coeff_downwind * component_downwind; } } diff --git a/framework/include/variables/MooseVariableFV.h b/framework/include/variables/MooseVariableFV.h index 6084c33d6482..34badfa92c81 100644 --- a/framework/include/variables/MooseVariableFV.h +++ b/framework/include/variables/MooseVariableFV.h @@ -651,6 +651,9 @@ class MooseVariableFV : public MooseVariableField const Elem * elem_side_to_extrapolate_from, const StateArg & state) const; + /// Function to get wether two term boundary expansion is used for the variable + const bool & getTwoTermBoundaryExpansion() const { return _two_term_boundary_expansion; } + protected: /** * clear finite volume caches diff --git a/framework/src/functormaterials/FunctorSmoother.C b/framework/src/functormaterials/FunctorSmoother.C index e457aa729a31..63dd7059bc66 100644 --- a/framework/src/functormaterials/FunctorSmoother.C +++ b/framework/src/functormaterials/FunctorSmoother.C @@ -102,7 +102,7 @@ FunctorSmootherTempl::FunctorSmootherTempl(const InputParameters & parameters _mesh.faceInfo(r_elem->neighbor_ptr(side_index), r_elem->neighbor_ptr(side_index)->which_neighbor_am_i(r_elem)); Moose::FaceArg face_arg{ - fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr}; + fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr}; if (face_arg.fi) { average += functor_in(face_arg, t); diff --git a/framework/src/fvbcs/FVBoundaryCondition.C b/framework/src/fvbcs/FVBoundaryCondition.C index 2a54c9823604..5d992c748a52 100644 --- a/framework/src/fvbcs/FVBoundaryCondition.C +++ b/framework/src/fvbcs/FVBoundaryCondition.C @@ -85,12 +85,13 @@ FVBoundaryCondition::FVBoundaryCondition(const InputParameters & parameters) Moose::FaceArg FVBoundaryCondition::singleSidedFaceArg(const FaceInfo * fi, const Moose::FV::LimiterType limiter_type, - const bool correct_skewness) const + const bool correct_skewness, + const Moose::StateArg * state_limiter) const { if (!fi) fi = _face_info; - return makeFace(*fi, limiter_type, true, correct_skewness); + return makeFace(*fi, limiter_type, true, correct_skewness, state_limiter); } bool diff --git a/framework/src/fvbcs/FVConstantScalarOutflowBC.C b/framework/src/fvbcs/FVConstantScalarOutflowBC.C index 07b2c66344cd..da240085b891 100644 --- a/framework/src/fvbcs/FVConstantScalarOutflowBC.C +++ b/framework/src/fvbcs/FVConstantScalarOutflowBC.C @@ -33,7 +33,10 @@ FVConstantScalarOutflowBC::computeQpResidual() "This boundary condition is for outflow but the flow is in the opposite direction of " "the boundary normal"); + const auto boundary_face = singleSidedFaceArg(); + const auto state = determineState(); + // This will either be second or first order accurate depending on whether the user has asked // for a two term expansion in their input file - return _normal * _velocity * _var.getBoundaryFaceValue(*_face_info, determineState()); + return _normal * _velocity * _var(boundary_face, state); } diff --git a/framework/src/fvbcs/FVFunctorDirichletBC.C b/framework/src/fvbcs/FVFunctorDirichletBC.C index 103f8854a051..3f75382921fa 100644 --- a/framework/src/fvbcs/FVFunctorDirichletBC.C +++ b/framework/src/fvbcs/FVFunctorDirichletBC.C @@ -59,11 +59,13 @@ FVFunctorDirichletBCTempl::boundaryValue(const FaceInfo & fi, if (!_use_other_side) return _functor(sfa, state); else if (fi.elemPtr() == sfa.face_side) - return _functor({&fi, Moose::FV::LimiterType::CentralDifference, true, false, fi.neighborPtr()}, - state); + return _functor( + {&fi, Moose::FV::LimiterType::CentralDifference, true, false, fi.neighborPtr(), nullptr}, + state); else - return _functor({&fi, Moose::FV::LimiterType::CentralDifference, true, false, fi.elemPtr()}, - state); + return _functor( + {&fi, Moose::FV::LimiterType::CentralDifference, true, false, fi.elemPtr(), nullptr}, + state); } template class FVFunctorDirichletBCTempl; diff --git a/framework/src/fviks/FVInterfaceKernel.C b/framework/src/fviks/FVInterfaceKernel.C index 11d831fe4b18..72212f45388b 100644 --- a/framework/src/fviks/FVInterfaceKernel.C +++ b/framework/src/fviks/FVInterfaceKernel.C @@ -241,7 +241,8 @@ Moose::FaceArg FVInterfaceKernel::singleSidedFaceArg(const MooseVariableFV & variable, const FaceInfo * fi, const Moose::FV::LimiterType limiter_type, - const bool correct_skewness) const + const bool correct_skewness, + const Moose::StateArg * state_limiter) const { if (!fi) fi = _face_info; @@ -255,7 +256,7 @@ FVInterfaceKernel::singleSidedFaceArg(const MooseVariableFV & variable, ? nullptr : (defined_on_elem_side ? fi->elemPtr() : fi->neighborPtr()); - return {fi, limiter_type, true, correct_skewness, elem}; + return {fi, limiter_type, true, correct_skewness, elem, state_limiter}; } bool diff --git a/framework/src/fvkernels/FVAdvection.C b/framework/src/fvkernels/FVAdvection.C index a731a791c274..a2bfc4653ec9 100644 --- a/framework/src/fvkernels/FVAdvection.C +++ b/framework/src/fvkernels/FVAdvection.C @@ -8,6 +8,7 @@ //* https://www.gnu.org/licenses/lgpl-2.1.html #include "FVAdvection.h" +#include "Steady.h" registerADMooseObject("MooseApp", FVAdvection); @@ -37,15 +38,33 @@ FVAdvection::FVAdvection(const InputParameters & params) getCheckedPointerParam("_fe_problem_base") ->setErrorOnJacobianNonzeroReallocation(false); } + + if (dynamic_cast(_app.getExecutioner())) + { + const MooseEnum not_available_with_steady("sou min_mod vanLeer quick venkatakrishnan"); + const std::string chosen_scheme = + static_cast(getParam("advected_interp_method")); + if (not_available_with_steady.find(chosen_scheme) != not_available_with_steady.items().end()) + paramError("advected_interp_method", + "The given advected interpolation cannot be used with steady-state runs!"); + } } ADReal FVAdvection::computeQpResidual() { + const auto state = determineState(); + const auto & limiter_time = _subproblem.isTransient() + ? Moose::StateArg(1, Moose::SolutionIterationType::Time) + : Moose::StateArg(1, Moose::SolutionIterationType::Nonlinear); + const bool elem_is_upwind = _velocity * _normal >= 0; - const auto face = - makeFace(*_face_info, Moose::FV::limiterType(_advected_interp_method), elem_is_upwind); - ADReal u_interface = _var(face, determineState()); + const auto face = makeFace(*_face_info, + Moose::FV::limiterType(_advected_interp_method), + elem_is_upwind, + false, + &limiter_time); + ADReal u_interface = _var(face, state); return _normal * _velocity * u_interface; } diff --git a/framework/src/fvkernels/FVFluxKernel.C b/framework/src/fvkernels/FVFluxKernel.C index ac426dd2cfad..0d3dabe8b345 100644 --- a/framework/src/fvkernels/FVFluxKernel.C +++ b/framework/src/fvkernels/FVFluxKernel.C @@ -249,12 +249,13 @@ FVFluxKernel::neighborArg(const bool correct_skewness) const Moose::FaceArg FVFluxKernel::singleSidedFaceArg(const FaceInfo * fi, const Moose::FV::LimiterType limiter_type, - const bool correct_skewness) const + const bool correct_skewness, + const Moose::StateArg * state_limiter) const { if (!fi) fi = _face_info; - return makeFace(*fi, limiter_type, true, correct_skewness); + return makeFace(*fi, limiter_type, true, correct_skewness, state_limiter); } bool diff --git a/framework/src/interfaces/FaceArgInterface.C b/framework/src/interfaces/FaceArgInterface.C index 6e75706117f2..110e371b1875 100644 --- a/framework/src/interfaces/FaceArgInterface.C +++ b/framework/src/interfaces/FaceArgInterface.C @@ -13,7 +13,8 @@ Moose::FaceArg FaceArgProducerInterface::makeFace(const FaceInfo & fi, const Moose::FV::LimiterType limiter_type, const bool elem_is_upwind, - const bool correct_skewness) const + const bool correct_skewness, + const Moose::StateArg * state_limiter) const { const bool defined_on_elem_side = hasFaceSide(fi, true); const bool defined_on_neighbor_side = hasFaceSide(fi, false); @@ -24,5 +25,5 @@ FaceArgProducerInterface::makeFace(const FaceInfo & fi, if (!defined_on_elem_side && !defined_on_neighbor_side) mooseError("No definition on either side"); - return {&fi, limiter_type, elem_is_upwind, correct_skewness, elem}; + return {&fi, limiter_type, elem_is_upwind, correct_skewness, elem, state_limiter}; } diff --git a/framework/src/limiters/Limiter.C b/framework/src/limiters/Limiter.C index a96b192e9cd9..dc5adb679b1a 100644 --- a/framework/src/limiters/Limiter.C +++ b/framework/src/limiters/Limiter.C @@ -14,6 +14,7 @@ #include "MinModLimiter.h" #include "SOULimiter.h" #include "QUICKLimiter.h" +#include "VenkatakrishnanLimiter.h" #include "MooseError.h" #include "MathFVUtils.h" @@ -23,8 +24,8 @@ namespace Moose { namespace FV { -const MooseEnum - moose_limiter_type("vanLeer=0 upwind=1 central_difference=2 min_mod=3 sou=4 quick=5", "upwind"); +const MooseEnum moose_limiter_type( + "vanLeer=0 upwind=1 central_difference=2 min_mod=3 sou=4 quick=5 venkatakrishnan=6", "upwind"); template std::unique_ptr> @@ -50,6 +51,9 @@ Limiter::build(const LimiterType limiter) case LimiterType::QUICK: return std::make_unique>(); + case LimiterType::Venkatakrishnan: + return std::make_unique>(); + default: mooseError("Unrecognized limiter type ", unsigned(limiter)); } @@ -79,6 +83,9 @@ limiterType(const InterpMethod interp_method) case InterpMethod::QUICK: return LimiterType::QUICK; + case InterpMethod::Venkatakrishnan: + return LimiterType::Venkatakrishnan; + default: mooseError("Unrecognized interpolation method type."); } diff --git a/framework/src/postprocessors/SideAdvectiveFluxIntegral.C b/framework/src/postprocessors/SideAdvectiveFluxIntegral.C index 393e999f01a1..3751638a9e39 100644 --- a/framework/src/postprocessors/SideAdvectiveFluxIntegral.C +++ b/framework/src/postprocessors/SideAdvectiveFluxIntegral.C @@ -91,19 +91,21 @@ SideAdvectiveFluxIntegralTempl::computeFaceInfoIntegral(const FaceInfo * // Get face value for velocity const auto vel_x = (_vel_x)(Moose::FaceArg( - {fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr}), + {fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr}), state); const auto vel_y = _vel_y - ? ((*_vel_y)(Moose::FaceArg( - {fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr}), - state)) + ? ((*_vel_y)( + Moose::FaceArg( + {fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr}), + state)) : 0; const auto vel_z = _vel_z - ? ((*_vel_z)(Moose::FaceArg( - {fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr}), - state)) + ? ((*_vel_z)( + Moose::FaceArg( + {fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr}), + state)) : 0; auto fi_normal = _current_elem == fi->elemPtr() ? fi->normal() : Point(-fi->normal()); @@ -111,7 +113,7 @@ SideAdvectiveFluxIntegralTempl::computeFaceInfoIntegral(const FaceInfo * const auto adv_quant_face = (*_adv_quant)( Moose::FaceArg( - {fi, Moose::FV::LimiterType::CentralDifference, elem_is_upwind, false, nullptr}), + {fi, Moose::FV::LimiterType::CentralDifference, elem_is_upwind, false, nullptr, nullptr}), state); return fi_normal * adv_quant_face * RealVectorValue(vel_x, vel_y, vel_z); diff --git a/framework/src/utils/MathFVUtils.C b/framework/src/utils/MathFVUtils.C index ad16b2f7dd89..ba6becde8585 100644 --- a/framework/src/utils/MathFVUtils.C +++ b/framework/src/utils/MathFVUtils.C @@ -60,19 +60,20 @@ onBoundary(const std::set & subs, const FaceInfo & fi) MooseEnum interpolationMethods() { - return MooseEnum("average upwind sou min_mod vanLeer quick skewness-corrected", "upwind"); + return MooseEnum("average upwind sou min_mod vanLeer quick venkatakrishnan skewness-corrected", + "upwind"); } InputParameters advectedInterpolationParameter() { auto params = emptyInputParameters(); - params.addParam( - "advected_interp_method", - interpolationMethods(), - "The interpolation to use for the advected quantity. Options are " - "'upwind', 'average', 'sou' (for second-order upwind), 'min_mod', 'vanLeer', 'quick', and " - "'skewness-corrected' with the default being 'upwind'."); + params.addParam("advected_interp_method", + interpolationMethods(), + "The interpolation to use for the advected quantity. Options are " + "'upwind', 'average', 'sou' (for second-order upwind), 'min_mod', " + "'vanLeer', 'quick', 'venkatakrishnan', and " + "'skewness-corrected' with the default being 'upwind'."); return params; } @@ -97,6 +98,8 @@ selectInterpolationMethod(const std::string & interp_method) return InterpMethod::SOU; else if (interp_method == "quick") return InterpMethod::QUICK; + else if (interp_method == "venkatakrishnan") + return InterpMethod::Venkatakrishnan; else mooseError("Interpolation method ", interp_method, @@ -114,7 +117,8 @@ setInterpolationMethod(const MooseObject & obj, interp_method = selectInterpolationMethod(interp_method_in); if (interp_method == InterpMethod::SOU || interp_method == InterpMethod::MinMod || - interp_method == InterpMethod::VanLeer || interp_method == InterpMethod::QUICK) + interp_method == InterpMethod::VanLeer || interp_method == InterpMethod::QUICK || + interp_method == InterpMethod::Venkatakrishnan) need_more_ghosting = true; return need_more_ghosting; diff --git a/modules/heat_transfer/src/fvbcs/FVFunctorConvectiveHeatFluxBC.C b/modules/heat_transfer/src/fvbcs/FVFunctorConvectiveHeatFluxBC.C index 59381376eb72..ad0daf3060ec 100644 --- a/modules/heat_transfer/src/fvbcs/FVFunctorConvectiveHeatFluxBC.C +++ b/modules/heat_transfer/src/fvbcs/FVFunctorConvectiveHeatFluxBC.C @@ -43,7 +43,7 @@ FVFunctorConvectiveHeatFluxBC::computeQpResidual() // Allow the functors to pick their side evaluation since either T_bulk or T_solid is likely not // defined on this boundary condition's side const Moose::FaceArg face{ - _face_info, Moose::FV::LimiterType::CentralDifference, true, false, nullptr}; + _face_info, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr}; const auto flux = _htc(face, determineState()) * (_T_bulk(face, determineState()) - _T_solid(face, determineState())); if (_is_solid) diff --git a/modules/navier_stokes/include/utils/Reconstructions.h b/modules/navier_stokes/include/utils/Reconstructions.h index 08b49c3c0c99..6a85ef224636 100644 --- a/modules/navier_stokes/include/utils/Reconstructions.h +++ b/modules/navier_stokes/include/utils/Reconstructions.h @@ -61,7 +61,8 @@ interpolateReconstruct(CellCenteredMapFunctor & output_functor, false, input_functor.hasFaceSide(*face, true) ? (input_functor.hasFaceSide(*face, false) ? nullptr : face->elemPtr()) - : face->neighborPtr()}; + : face->neighborPtr(), + nullptr}; auto face_value = input_functor(face_arg, time); std::pair * neighbor_pair = nullptr; if (face->neighborPtr() && face->neighborPtr() != libMesh::remote_elem) diff --git a/modules/navier_stokes/src/fvkernels/INSFVAdvectionKernel.C b/modules/navier_stokes/src/fvkernels/INSFVAdvectionKernel.C index 3d3fc9e45ae4..a1cfc3f7f1b8 100644 --- a/modules/navier_stokes/src/fvkernels/INSFVAdvectionKernel.C +++ b/modules/navier_stokes/src/fvkernels/INSFVAdvectionKernel.C @@ -13,6 +13,8 @@ #include "RelationshipManager.h" #include "NSFVUtils.h" #include "FVBoundaryScalarLagrangeMultiplierConstraint.h" +#include "Limiter.h" +#include "Steady.h" InputParameters INSFVAdvectionKernel::validParams() @@ -55,6 +57,28 @@ INSFVAdvectionKernel::INSFVAdvectionKernel(const InputParameters & params) }; param_check("force_boundary_execution"); + + if (_var.getTwoTermBoundaryExpansion() && + !(_advected_interp_method == Moose::FV::InterpMethod::Upwind || + _advected_interp_method == Moose::FV::InterpMethod::Average || + _advected_interp_method == Moose::FV::InterpMethod::HarmonicAverage || + _advected_interp_method == Moose::FV::InterpMethod::SkewCorrectedAverage)) + mooseWarning( + "Second order upwind limiting is not supported when `two_term_boundary_expansion " + "= true` for the limited variable. Use at your own risk or please consider " + "setting `two_term_boundary_expansion = false` in the advected variable parameters or " + "changing your " + "'advected_interp_method' of the kernel to first order methods (`upwind`, `average`)"); + + if (dynamic_cast(_app.getExecutioner())) + { + const MooseEnum not_available_with_steady("sou min_mod vanLeer quick venkatakrishnan"); + const std::string chosen_scheme = + static_cast(getParam("advected_interp_method")); + if (not_available_with_steady.find(chosen_scheme) != not_available_with_steady.items().end()) + paramError("advected_interp_method", + "The given advected interpolation cannot be used with steady-state runs!"); + } } void diff --git a/modules/navier_stokes/src/fvkernels/INSFVEnergyAdvection.C b/modules/navier_stokes/src/fvkernels/INSFVEnergyAdvection.C index 97ead95634a9..d311363004e0 100644 --- a/modules/navier_stokes/src/fvkernels/INSFVEnergyAdvection.C +++ b/modules/navier_stokes/src/fvkernels/INSFVEnergyAdvection.C @@ -36,9 +36,14 @@ ADReal INSFVEnergyAdvection::computeQpResidual() { const auto v = velocity(); + const auto & limiter_time = _subproblem.isTransient() + ? Moose::StateArg(1, Moose::SolutionIterationType::Time) + : Moose::StateArg(1, Moose::SolutionIterationType::Nonlinear); const auto adv_quant_face = _adv_quant(makeFace(*_face_info, limiterType(_advected_interp_method), - MetaPhysicL::raw_value(v) * _normal > 0), + MetaPhysicL::raw_value(v) * _normal > 0, + false, + &limiter_time), determineState()); return _normal * v * adv_quant_face; } diff --git a/modules/navier_stokes/src/fvkernels/INSFVMomentumAdvection.C b/modules/navier_stokes/src/fvkernels/INSFVMomentumAdvection.C index 5fc190c0c5ef..5f1628c67f09 100644 --- a/modules/navier_stokes/src/fvkernels/INSFVMomentumAdvection.C +++ b/modules/navier_stokes/src/fvkernels/INSFVMomentumAdvection.C @@ -123,8 +123,15 @@ INSFVMomentumAdvection::computeResidualsAndAData(const FaceInfo & fi) else // we are an internal fluid flow face { const bool elem_is_upwind = MetaPhysicL::raw_value(v_face) * _normal > 0; - const Moose::FaceArg advected_face_arg{ - &fi, limiterType(_advected_interp_method), elem_is_upwind, correct_skewness, nullptr}; + const auto & limiter_time = _subproblem.isTransient() + ? Moose::StateArg(1, Moose::SolutionIterationType::Time) + : Moose::StateArg(1, Moose::SolutionIterationType::Nonlinear); + const Moose::FaceArg advected_face_arg{&fi, + limiterType(_advected_interp_method), + elem_is_upwind, + correct_skewness, + nullptr, + &limiter_time}; if (const auto [is_jump, eps_elem_face, eps_neighbor_face] = NS::isPorosityJumpFace(epsilon(), fi, state); is_jump) diff --git a/modules/navier_stokes/src/fvkernels/INSFVScalarFieldAdvection.C b/modules/navier_stokes/src/fvkernels/INSFVScalarFieldAdvection.C index ae30650de2e7..4072387362b2 100644 --- a/modules/navier_stokes/src/fvkernels/INSFVScalarFieldAdvection.C +++ b/modules/navier_stokes/src/fvkernels/INSFVScalarFieldAdvection.C @@ -46,16 +46,24 @@ ADReal INSFVScalarFieldAdvection::computeQpResidual() { const auto state = determineState(); + const auto & limiter_time = _subproblem.isTransient() + ? Moose::StateArg(1, Moose::SolutionIterationType::Time) + : Moose::StateArg(1, Moose::SolutionIterationType::Nonlinear); ADRealVectorValue advection_velocity; if (_add_slip_model) { + Moose::FaceArg face_arg; if (onBoundary(*_face_info)) face_arg = singleSidedFaceArg(); else - face_arg = Moose::FaceArg{ - _face_info, Moose::FV::LimiterType::CentralDifference, true, false, nullptr}; + face_arg = Moose::FaceArg{_face_info, + Moose::FV::LimiterType::CentralDifference, + true, + false, + nullptr, + &limiter_time}; ADRealVectorValue velocity_slip_vel_vec; if (_dim >= 1) @@ -71,7 +79,9 @@ INSFVScalarFieldAdvection::computeQpResidual() advection_velocity += v; const auto var_face = _var(makeFace(*_face_info, limiterType(_advected_interp_method), - MetaPhysicL::raw_value(v) * _normal > 0), + MetaPhysicL::raw_value(v) * _normal > 0, + false, + &limiter_time), state); return _normal * advection_velocity * var_face; diff --git a/modules/navier_stokes/src/fvkernels/INSFVTKEDSourceSink.C b/modules/navier_stokes/src/fvkernels/INSFVTKEDSourceSink.C index ea2a076418e5..31a4a14fd6cd 100644 --- a/modules/navier_stokes/src/fvkernels/INSFVTKEDSourceSink.C +++ b/modules/navier_stokes/src/fvkernels/INSFVTKEDSourceSink.C @@ -149,7 +149,7 @@ INSFVTKEDSourceSink::computeQpResidual() const bool defined_on_elem_side = _var.hasFaceSide(*fi, true); const Elem * const loc_elem = defined_on_elem_side ? &fi->elem() : fi->neighborPtr(); const Moose::FaceArg facearg = { - fi, Moose::FV::LimiterType::CentralDifference, false, false, loc_elem}; + fi, Moose::FV::LimiterType::CentralDifference, false, false, loc_elem, nullptr}; destruction += 2.0 * TKE_old * _mu(facearg, state) / rho / Utility::pow<2>(distance_vec[i]) / tot_weight; } diff --git a/modules/navier_stokes/src/fvkernels/INSFVTKESourceSink.C b/modules/navier_stokes/src/fvkernels/INSFVTKESourceSink.C index b5902edcd259..354ba1775921 100644 --- a/modules/navier_stokes/src/fvkernels/INSFVTKESourceSink.C +++ b/modules/navier_stokes/src/fvkernels/INSFVTKESourceSink.C @@ -163,7 +163,7 @@ INSFVTKESourceSink::computeQpResidual() const bool defined_on_elem_side = _var.hasFaceSide(*fi, true); const Elem * const loc_elem = defined_on_elem_side ? &fi->elem() : fi->neighborPtr(); const Moose::FaceArg facearg = { - fi, Moose::FV::LimiterType::CentralDifference, false, false, loc_elem}; + fi, Moose::FV::LimiterType::CentralDifference, false, false, loc_elem, nullptr}; const ADReal wall_mut = _mu_t(facearg, state); const ADReal wall_mu = _mu(facearg, state); diff --git a/modules/navier_stokes/src/fvkernels/WCNSFV2PMomentumAdvectionSlip.C b/modules/navier_stokes/src/fvkernels/WCNSFV2PMomentumAdvectionSlip.C index 21c20a90ca7b..7dab106d0fe3 100644 --- a/modules/navier_stokes/src/fvkernels/WCNSFV2PMomentumAdvectionSlip.C +++ b/modules/navier_stokes/src/fvkernels/WCNSFV2PMomentumAdvectionSlip.C @@ -73,7 +73,8 @@ WCNSFV2PMomentumAdvectionSlip::computeResidualsAndAData(const FaceInfo & fi) if (on_boundary) face_arg = singleSidedFaceArg(); else - face_arg = Moose::FaceArg{&fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr}; + face_arg = Moose::FaceArg{ + &fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr}; ADRealVectorValue u_slip_vel_vec; if (_dim == 1) @@ -118,8 +119,12 @@ WCNSFV2PMomentumAdvectionSlip::computeResidualsAndAData(const FaceInfo & fi) else // we are an internal fluid flow face { const bool elem_is_upwind = MetaPhysicL::raw_value(u_slip_vel_vec) * _normal > 0; - const Moose::FaceArg advected_face_arg{ - &fi, limiterType(_advected_interp_method), elem_is_upwind, correct_skewness, nullptr}; + const Moose::FaceArg advected_face_arg{&fi, + limiterType(_advected_interp_method), + elem_is_upwind, + correct_skewness, + nullptr, + nullptr}; if (const auto [is_jump, eps_elem_face, eps_neighbor_face] = NS::isPorosityJumpFace(epsilon(), fi, state); is_jump) diff --git a/modules/navier_stokes/src/postprocessors/IntegralDirectedSurfaceForce.C b/modules/navier_stokes/src/postprocessors/IntegralDirectedSurfaceForce.C index 9969dea156a6..604f7ad61442 100644 --- a/modules/navier_stokes/src/postprocessors/IntegralDirectedSurfaceForce.C +++ b/modules/navier_stokes/src/postprocessors/IntegralDirectedSurfaceForce.C @@ -74,8 +74,8 @@ IntegralDirectedSurfaceForce::computeFaceInfoIntegral(const FaceInfo * fi) mooseAssert(fi, "We should have a face info in " + name()); const auto state = determineState(); - const auto face_arg = - Moose::FaceArg({fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr}); + const auto face_arg = Moose::FaceArg( + {fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr}); const auto elem_arg = Moose::ElemArg({fi->elemPtr(), false}); RealTensorValue pressure_term; diff --git a/modules/navier_stokes/src/postprocessors/MassFluxWeightedFlowRate.C b/modules/navier_stokes/src/postprocessors/MassFluxWeightedFlowRate.C index 3d200dbcc3c4..dc614f188e9f 100644 --- a/modules/navier_stokes/src/postprocessors/MassFluxWeightedFlowRate.C +++ b/modules/navier_stokes/src/postprocessors/MassFluxWeightedFlowRate.C @@ -70,7 +70,8 @@ MassFluxWeightedFlowRate::computeFaceInfoIntegral([[maybe_unused]] const FaceInf Moose::FV::limiterType(_advected_interp_method), MetaPhysicL::raw_value(vel) * fi->normal() > 0, correct_skewness, - _adv_quant->hasFaceSide(*fi, true) ? fi->elemPtr() : fi->neighborPtr()}); + _adv_quant->hasFaceSide(*fi, true) ? fi->elemPtr() : fi->neighborPtr(), + nullptr}); auto dens = _density(face_arg, state); const auto adv_quant_face = MetaPhysicL::raw_value(dens * (*_adv_quant)(face_arg, state)); _mdot += fi->faceArea() * fi->faceCoord() * MetaPhysicL::raw_value(dens) * fi->normal() * vel; diff --git a/modules/navier_stokes/src/postprocessors/PressureDrop.C b/modules/navier_stokes/src/postprocessors/PressureDrop.C index 977c2abb30c1..6b2688ce0c0a 100644 --- a/modules/navier_stokes/src/postprocessors/PressureDrop.C +++ b/modules/navier_stokes/src/postprocessors/PressureDrop.C @@ -238,14 +238,15 @@ PressureDrop::computeFaceInfoWeightedPressureIntegral(const FaceInfo * const fi) Moose::FV::limiterType(_weight_interp_method), MetaPhysicL::raw_value((*_weighting_functor)(elem_arg, state)) * fi->normal() > 0, correct_skewness, - elem}); + elem, + nullptr}); const auto face_weighting = MetaPhysicL::raw_value((*_weighting_functor)(ssf, state)); return fi->normal() * face_weighting * _pressure(ssf, state); } else { const auto ssf = Moose::FaceArg( - {fi, Moose::FV::limiterType(_weight_interp_method), true, correct_skewness, elem}); + {fi, Moose::FV::limiterType(_weight_interp_method), true, correct_skewness, elem, nullptr}); return _pressure(ssf, state); } } @@ -268,7 +269,8 @@ PressureDrop::computeFaceInfoWeightIntegral(const FaceInfo * fi) const Moose::FV::limiterType(_weight_interp_method), MetaPhysicL::raw_value((*_weighting_functor)(elem_arg, state)) * fi->normal() > 0, correct_skewness, - _weighting_functor->hasFaceSide(*fi, true) ? fi->elemPtr() : fi->neighborPtr()}); + _weighting_functor->hasFaceSide(*fi, true) ? fi->elemPtr() : fi->neighborPtr(), + nullptr}); return fi->normal() * MetaPhysicL::raw_value((*_weighting_functor)(ssf, state)); } else diff --git a/modules/navier_stokes/src/postprocessors/VolumetricFlowRate.C b/modules/navier_stokes/src/postprocessors/VolumetricFlowRate.C index 7c9220f8809d..540146720b07 100644 --- a/modules/navier_stokes/src/postprocessors/VolumetricFlowRate.C +++ b/modules/navier_stokes/src/postprocessors/VolumetricFlowRate.C @@ -147,7 +147,8 @@ VolumetricFlowRate::computeFaceInfoIntegral(const FaceInfo * fi) Moose::FV::limiterType(_advected_interp_method), MetaPhysicL::raw_value(vel) * fi->normal() > 0, correct_skewness, - elem}), + elem, + nullptr}), state)); return fi->normal() * adv_quant_face * vel; } diff --git a/modules/navier_stokes/src/userobjects/INSFVRhieChowInterpolator.C b/modules/navier_stokes/src/userobjects/INSFVRhieChowInterpolator.C index 619e6335cd75..fa5b49240454 100644 --- a/modules/navier_stokes/src/userobjects/INSFVRhieChowInterpolator.C +++ b/modules/navier_stokes/src/userobjects/INSFVRhieChowInterpolator.C @@ -544,8 +544,12 @@ INSFVRhieChowInterpolator::getVelocity(const Moose::FV::InterpMethod m, if (Moose::FV::onBoundary(*this, fi)) { const Elem * const boundary_elem = hasBlocks(elem->subdomain_id()) ? elem : neighbor; - const Moose::FaceArg boundary_face{ - &fi, Moose::FV::LimiterType::CentralDifference, true, correct_skewness, boundary_elem}; + const Moose::FaceArg boundary_face{&fi, + Moose::FV::LimiterType::CentralDifference, + true, + correct_skewness, + boundary_elem, + nullptr}; auto velocity = vel(boundary_face, time); incorporate_mesh_velocity(boundary_face, velocity); @@ -559,7 +563,7 @@ INSFVRhieChowInterpolator::getVelocity(const Moose::FV::InterpMethod m, VectorValue velocity; Moose::FaceArg face{ - &fi, Moose::FV::LimiterType::CentralDifference, true, correct_skewness, nullptr}; + &fi, Moose::FV::LimiterType::CentralDifference, true, correct_skewness, nullptr, nullptr}; // Create the average face velocity (not corrected using RhieChow yet) velocity(0) = (*u)(face, time); if (v) @@ -625,7 +629,7 @@ INSFVRhieChowInterpolator::getVelocity(const Moose::FV::InterpMethod m, // Compute the corrected interpolated face value Moose::FaceArg face{ - &fi, Moose::FV::LimiterType::CentralDifference, true, correct_skewness, nullptr}; + &fi, Moose::FV::LimiterType::CentralDifference, true, correct_skewness, nullptr, nullptr}; interp_vf = 0.0; for (const auto i : make_range(_volumetric_force.size())) @@ -654,7 +658,7 @@ INSFVRhieChowInterpolator::getVelocity(const Moose::FV::InterpMethod m, elem_has_fi ? side : loc_neighbor->which_neighbor_am_i(elem)); Moose::FaceArg loc_face{ - fi_loc, Moose::FV::LimiterType::CentralDifference, true, correct_skewness, elem}; + fi_loc, Moose::FV::LimiterType::CentralDifference, true, correct_skewness, elem, nullptr}; MooseMeshUtils::coordTransformFactor( elem->vertex_average(), coord_multiplier, coord_type, rz_radial_coord); @@ -689,7 +693,7 @@ INSFVRhieChowInterpolator::getVelocity(const Moose::FV::InterpMethod m, elem_has_fi ? side : loc_elem->which_neighbor_am_i(neighbor)); Moose::FaceArg loc_face{ - fi_loc, Moose::FV::LimiterType::CentralDifference, true, correct_skewness, elem}; + fi_loc, Moose::FV::LimiterType::CentralDifference, true, correct_skewness, elem, nullptr}; MooseMeshUtils::coordTransformFactor( neighbor->vertex_average(), coord_multiplier, coord_type, rz_radial_coord); @@ -730,7 +734,7 @@ INSFVRhieChowInterpolator::getVelocity(const Moose::FV::InterpMethod m, { Real value = 0.0; Moose::FaceArg loc_face{ - &fi, Moose::FV::LimiterType::CentralDifference, true, correct_skewness, nullptr}; + &fi, Moose::FV::LimiterType::CentralDifference, true, correct_skewness, nullptr, nullptr}; for (const auto i : make_range(_volumetric_force.size())) value += (*_volumetric_force[i])(loc_face, time) * (face_normal * fi.normal()); diff --git a/modules/navier_stokes/src/userobjects/INSFVRhieChowInterpolatorSegregated.C b/modules/navier_stokes/src/userobjects/INSFVRhieChowInterpolatorSegregated.C index 3b7da0b91cf1..3fe5dc55962f 100644 --- a/modules/navier_stokes/src/userobjects/INSFVRhieChowInterpolatorSegregated.C +++ b/modules/navier_stokes/src/userobjects/INSFVRhieChowInterpolatorSegregated.C @@ -136,7 +136,7 @@ INSFVRhieChowInterpolatorSegregated::initFaceVelocities() if (_u->isInternalFace(*fi)) { const Moose::FaceArg face{ - fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr}; + fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr}; _face_velocity[fi->id()] = MetaPhysicL::raw_value((*_vel)(face, Moose::currentState())); } @@ -147,7 +147,7 @@ INSFVRhieChowInterpolatorSegregated::initFaceVelocities() hasBlocks(fi->elemPtr()->subdomain_id()) ? fi->elemPtr() : fi->neighborPtr(); const Moose::FaceArg boundary_face{ - fi, Moose::FV::LimiterType::CentralDifference, true, false, boundary_elem}; + fi, Moose::FV::LimiterType::CentralDifference, true, false, boundary_elem, nullptr}; _face_velocity[fi->id()] = MetaPhysicL::raw_value((*_vel)(boundary_face, Moose::currentState())); @@ -185,7 +185,7 @@ INSFVRhieChowInterpolatorSegregated::computeFaceVelocity() if (_u->isInternalFace(*fi)) { const Moose::FaceArg face{ - fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr}; + fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr}; RealVectorValue Ainv; RealVectorValue HbyA = MetaPhysicL::raw_value(_HbyA(face, time_arg)); @@ -207,7 +207,7 @@ INSFVRhieChowInterpolatorSegregated::computeFaceVelocity() const Elem * const boundary_elem = hasBlocks(fi->elemPtr()->subdomain_id()) ? fi->elemPtr() : fi->neighborPtr(); const Moose::FaceArg boundary_face{ - fi, Moose::FV::LimiterType::CentralDifference, true, false, boundary_elem}; + fi, Moose::FV::LimiterType::CentralDifference, true, false, boundary_elem, nullptr}; // If we have a dirichlet boundary conditions, this sill give us the exact value of the // velocity on the face as expected (see populateHbyA()) @@ -310,7 +310,7 @@ INSFVRhieChowInterpolatorSegregated::populateHbyA( hasBlocks(fi->elemPtr()->subdomain_id()) ? fi->elemPtr() : fi->neighborPtr(); const Moose::FaceArg boundary_face{ - fi, Moose::FV::LimiterType::CentralDifference, true, false, boundary_elem}; + fi, Moose::FV::LimiterType::CentralDifference, true, false, boundary_elem, nullptr}; if (_u->isDirichletBoundaryFace(*fi, boundary_elem, Moose::currentState())) _HbyA[fi->id()] = -MetaPhysicL::raw_value((*_vel)(boundary_face, Moose::currentState())); diff --git a/modules/navier_stokes/src/userobjects/RhieChowMassFlux.C b/modules/navier_stokes/src/userobjects/RhieChowMassFlux.C index 288788806778..4e33d9722d40 100644 --- a/modules/navier_stokes/src/userobjects/RhieChowMassFlux.C +++ b/modules/navier_stokes/src/userobjects/RhieChowMassFlux.C @@ -216,7 +216,7 @@ RhieChowMassFlux::initFaceMassFlux() hasBlocks(fi->elemPtr()->subdomain_id()) ? fi->elemPtr() : fi->neighborPtr(); const Moose::FaceArg boundary_face{ - fi, Moose::FV::LimiterType::CentralDifference, true, false, boundary_elem}; + fi, Moose::FV::LimiterType::CentralDifference, true, false, boundary_elem, nullptr}; const Real face_rho = _rho(boundary_face, time_arg); for (const auto dim_i : index_range(_vel)) @@ -405,7 +405,7 @@ RhieChowMassFlux::populateCouplingFunctors( if (_vel[0]->isDirichletBoundaryFace(*fi)) { const Moose::FaceArg boundary_face{ - fi, Moose::FV::LimiterType::CentralDifference, true, false, elem_info.elem()}; + fi, Moose::FV::LimiterType::CentralDifference, true, false, elem_info.elem(), nullptr}; face_rho = _rho(boundary_face, Moose::currentState()); for (const auto dim_i : make_range(_dim)) diff --git a/modules/navier_stokes/src/utils/NSFVUtils.C b/modules/navier_stokes/src/utils/NSFVUtils.C index 747b2e08ad90..2db25a74299e 100644 --- a/modules/navier_stokes/src/utils/NSFVUtils.C +++ b/modules/navier_stokes/src/utils/NSFVUtils.C @@ -68,9 +68,9 @@ isPorosityJumpFace(const Moose::Functor & porosity, "Porosity should have blocks on both elem and neighbor"); const Moose::FaceArg face_elem{ - &fi, Moose::FV::LimiterType::CentralDifference, true, false, fi.elemPtr()}; + &fi, Moose::FV::LimiterType::CentralDifference, true, false, fi.elemPtr(), nullptr}; const Moose::FaceArg face_neighbor{ - &fi, Moose::FV::LimiterType::CentralDifference, true, false, fi.neighborPtr()}; + &fi, Moose::FV::LimiterType::CentralDifference, true, false, fi.neighborPtr(), nullptr}; const auto eps_elem = porosity(face_elem, time), eps_neighbor = porosity(face_neighbor, time); return {!MooseUtils::relativeFuzzyEqual(eps_elem, eps_neighbor), eps_elem, eps_neighbor}; } diff --git a/modules/navier_stokes/src/variables/BernoulliPressureVariable.C b/modules/navier_stokes/src/variables/BernoulliPressureVariable.C index 1a2b2ee75c24..c605dda566da 100644 --- a/modules/navier_stokes/src/variables/BernoulliPressureVariable.C +++ b/modules/navier_stokes/src/variables/BernoulliPressureVariable.C @@ -68,7 +68,8 @@ BernoulliPressureVariable::elemIsUpwind(const Elem & elem, const FaceInfo & fi, const Moose::StateArg & time) const { - const Moose::FaceArg face{&fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr}; + const Moose::FaceArg face{ + &fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr}; const VectorValue vel_face{(*_u)(face, time), (*_v)(face, time), (*_w)(face, time)}; const bool fi_elem_is_upwind = vel_face * fi.normal() > 0; @@ -131,9 +132,9 @@ BernoulliPressureVariable::getDirichletBoundaryFaceValue(const FaceInfo & fi, #endif const Moose::FaceArg face_elem{ - &fi, Moose::FV::LimiterType::CentralDifference, true, false, &fi.elem()}; + &fi, Moose::FV::LimiterType::CentralDifference, true, false, &fi.elem(), nullptr}; const Moose::FaceArg face_neighbor{ - &fi, Moose::FV::LimiterType::CentralDifference, true, false, fi.neighborPtr()}; + &fi, Moose::FV::LimiterType::CentralDifference, true, false, fi.neighborPtr(), nullptr}; const auto [elem_is_upwind, vel_face] = elemIsUpwind(*elem, fi, time); const auto & vel_elem = vel_face; diff --git a/modules/navier_stokes/src/variables/INSFVVelocityVariable.C b/modules/navier_stokes/src/variables/INSFVVelocityVariable.C index f2f17a353ea6..aa214a7af666 100644 --- a/modules/navier_stokes/src/variables/INSFVVelocityVariable.C +++ b/modules/navier_stokes/src/variables/INSFVVelocityVariable.C @@ -273,6 +273,7 @@ INSFVVelocityVariable::adGradSln(const Elem * const elem, Moose::FV::LimiterType::CentralDifference, true, correct_skewness, + nullptr, nullptr}), time); else diff --git a/modules/navier_stokes/test/tests/finite_volume/ins/advection-schemes/gold/min_mod.e b/modules/navier_stokes/test/tests/finite_volume/ins/advection-schemes/gold/min_mod.e deleted file mode 100644 index 49f615ca20a8..000000000000 Binary files a/modules/navier_stokes/test/tests/finite_volume/ins/advection-schemes/gold/min_mod.e and /dev/null differ diff --git a/modules/navier_stokes/test/tests/finite_volume/ins/advection-schemes/gold/quick.e b/modules/navier_stokes/test/tests/finite_volume/ins/advection-schemes/gold/quick.e deleted file mode 100644 index dfecc7a513ad..000000000000 Binary files a/modules/navier_stokes/test/tests/finite_volume/ins/advection-schemes/gold/quick.e and /dev/null differ diff --git a/modules/navier_stokes/test/tests/finite_volume/ins/advection-schemes/gold/sou.e b/modules/navier_stokes/test/tests/finite_volume/ins/advection-schemes/gold/sou.e deleted file mode 100644 index 9009f9b5e75d..000000000000 Binary files a/modules/navier_stokes/test/tests/finite_volume/ins/advection-schemes/gold/sou.e and /dev/null differ diff --git a/modules/navier_stokes/test/tests/finite_volume/ins/advection-schemes/gold/vanLeer.e b/modules/navier_stokes/test/tests/finite_volume/ins/advection-schemes/gold/vanLeer.e deleted file mode 100644 index 40370c331108..000000000000 Binary files a/modules/navier_stokes/test/tests/finite_volume/ins/advection-schemes/gold/vanLeer.e and /dev/null differ diff --git a/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/min_mod.e b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/min_mod.e new file mode 100644 index 000000000000..0c5570f3479f Binary files /dev/null and b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/min_mod.e differ diff --git a/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/min_mod_run.e b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/min_mod_run.e new file mode 100644 index 000000000000..7b045c640e86 Binary files /dev/null and b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/min_mod_run.e differ diff --git a/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/quick.e b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/quick.e new file mode 100644 index 000000000000..b28b836b8878 Binary files /dev/null and b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/quick.e differ diff --git a/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/quick_run.e b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/quick_run.e new file mode 100644 index 000000000000..8a6325e7aa6c Binary files /dev/null and b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/quick_run.e differ diff --git a/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/upwind.e b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/upwind.e new file mode 100644 index 000000000000..c70206540233 Binary files /dev/null and b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/upwind.e differ diff --git a/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/upwind_run.e b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/upwind_run.e new file mode 100644 index 000000000000..935d13262b6f Binary files /dev/null and b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/upwind_run.e differ diff --git a/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/vanLeer.e b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/vanLeer.e new file mode 100644 index 000000000000..7a004e165d69 Binary files /dev/null and b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/vanLeer.e differ diff --git a/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/vanLeer_run.e b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/vanLeer_run.e new file mode 100644 index 000000000000..99280fee2f7f Binary files /dev/null and b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/gold/vanLeer_run.e differ diff --git a/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/lid-driven-segregated.i b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/lid-driven-segregated.i new file mode 100644 index 000000000000..fa95501f020d --- /dev/null +++ b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/lid-driven-segregated.i @@ -0,0 +1,169 @@ +mu = .001 +rho = 1 + +pressure_tag = "pressure_grad" + +[GlobalParams] + rhie_chow_user_object = 'rc' + advected_interp_method = 'min_mod' #average upwind sou min_mod vanLeer quick venkatakrishnan skewness-corrected + velocity_interp_method = 'rc' +[] + +[Mesh] + [gen] + type = GeneratedMeshGenerator + dim = 2 + xmin = 0 + xmax = 1 + ymin = 0 + ymax = 1 + nx = 25 + ny = 25 + [] +[] + +[Problem] + nl_sys_names = 'u_system v_system pressure_system' + previous_nl_solution_required = true +[] + +[UserObjects] + [rc] + type = INSFVRhieChowInterpolatorSegregated + u = vel_x + v = vel_y + pressure = pressure + [] +[] + +[Variables] + [vel_x] + type = INSFVVelocityVariable + initial_condition = 0.0 + solver_sys = u_system + two_term_boundary_expansion = false + [] + [vel_y] + type = INSFVVelocityVariable + initial_condition = 0.0 + solver_sys = v_system + two_term_boundary_expansion = false + [] + [pressure] + type = INSFVPressureVariable + solver_sys = pressure_system + initial_condition = 0.2 + two_term_boundary_expansion = false + [] +[] + +[FVKernels] + [u_advection] + type = INSFVMomentumAdvection + variable = vel_x + rho = ${rho} + momentum_component = 'x' + [] + [u_viscosity] + type = INSFVMomentumDiffusion + variable = vel_x + mu = ${mu} + momentum_component = 'x' + [] + [u_pressure] + type = INSFVMomentumPressure + variable = vel_x + momentum_component = 'x' + pressure = pressure + extra_vector_tags = ${pressure_tag} + [] + [v_advection] + type = INSFVMomentumAdvection + variable = vel_y + rho = ${rho} + momentum_component = 'y' + [] + [v_viscosity] + type = INSFVMomentumDiffusion + variable = vel_y + mu = ${mu} + momentum_component = 'y' + [] + [v_pressure] + type = INSFVMomentumPressure + variable = vel_y + momentum_component = 'y' + pressure = pressure + extra_vector_tags = ${pressure_tag} + [] + [p_diffusion] + type = FVAnisotropicDiffusion + variable = pressure + coeff = "Ainv" + coeff_interp_method = 'average' + [] + [p_source] + type = FVDivergence + variable = pressure + vector_field = "HbyA" + force_boundary_execution = true + [] +[] + +[FVBCs] + [top_x] + type = INSFVNoSlipWallBC + variable = vel_x + boundary = 'top' + function = 1 + [] + [no_slip_x] + type = INSFVNoSlipWallBC + variable = vel_x + boundary = 'left right bottom' + function = 0 + [] + [no_slip_y] + type = INSFVNoSlipWallBC + variable = vel_y + boundary = 'left right top bottom' + function = 0 + [] +[] + +[Executioner] + type = SIMPLENonlinearAssembly + rhie_chow_user_object = 'rc' + momentum_systems = 'u_system v_system' + pressure_system = 'pressure_system' + pressure_gradient_tag = ${pressure_tag} + momentum_equation_relaxation = 0.5 + pressure_variable_relaxation = 0.2 + num_iterations = 1000 + pressure_absolute_tolerance = 1e-13 + momentum_absolute_tolerance = 1e-13 + 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' + + momentum_l_abs_tol = 1e-14 + pressure_l_abs_tol = 1e-14 + momentum_l_max_its = 30 + pressure_l_max_its = 30 + momentum_l_tol = 0.0 + pressure_l_tol = 0.0 + print_fields = false + + pin_pressure = true + pressure_pin_value = 0.0 + pressure_pin_point = '0.01 0.099 0.0' +[] + +[Outputs] + exodus = true + csv = false + perf_graph = false + print_nonlinear_residuals = false + print_linear_residuals = true +[] diff --git a/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/tests b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/tests new file mode 100644 index 000000000000..ecc297a13ab9 --- /dev/null +++ b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven-segregated/tests @@ -0,0 +1,89 @@ +[Tests] + design = 'INSFVMomentumAdvection.md' + issues = '#28891' + [advection_limiting_schemes] + requirement = 'The system shall be able to perform a variety of limiting schemes when solving fluid flow equations in the segreagated solver. These schemes include' + [upwind] + type = 'Exodiff' + input = 'lid-driven-segregated.i' + exodiff = 'upwind.e' + detail = 'and reach converged results with upwind advection scheme.' + cli_args = 'GlobalParams/advected_interp_method=upwind Outputs/file_base=upwind' + heavy = true + recover = false # we don't support recovery for SIMPLE yet + max_threads = 1 # see libmesh issue #3808 + [] + [upwind_run] + type = 'Exodiff' + input = 'lid-driven-segregated.i' + exodiff = 'upwind_run.e' + detail = 'and pass debugging checks with segregated solvers with upwind advection scheme.' + cli_args = 'GlobalParams/advected_interp_method=upwind Executioner/num_iterations=10 Outputs/file_base=upwind_run' + recover = false # we don't support recovery for SIMPLE yet + max_threads = 1 # see libmesh issue #3808 + [] + [vanLeer] + type = 'Exodiff' + input = 'lid-driven-segregated.i' + exodiff = 'vanLeer.e' + detail = 'and reach converged results with van Leer limiter.' + cli_args = 'GlobalParams/advected_interp_method=vanLeer Outputs/file_base=vanLeer Outputs/file_base=vanLeer Mesh/gen/nx=10 Mesh/gen/ny=10' + heavy = true + recover = false # we don't support recovery for SIMPLE yet + max_threads = 1 # see libmesh issue #3808 + [] + [vanLeer_run] + type = 'Exodiff' + input = 'lid-driven-segregated.i' + exodiff = 'vanLeer_run.e' + detail = 'and pass debugging checks with segregated solvers with van Leer limiter.' + cli_args = 'GlobalParams/advected_interp_method=vanLeer Executioner/num_iterations=10 Outputs/file_base=vanLeer_run Mesh/gen/nx=10 Mesh/gen/ny=10' + recover = false # we don't support recovery for SIMPLE yet + max_threads = 1 # see libmesh issue #3808 + [] + [min_mod] + type = 'Exodiff' + input = 'lid-driven-segregated.i' + exodiff = 'min_mod.e' + detail = 'and reach converged results with min-mod limiter.' + cli_args = 'GlobalParams/advected_interp_method=min_mod Outputs/file_base=min_mod' + heavy = true + recover = false # we don't support recovery for SIMPLE yet + max_threads = 1 # see libmesh issue #3808 + [] + [min_mod_run] + type = 'Exodiff' + input = 'lid-driven-segregated.i' + exodiff = 'min_mod_run.e' + detail = 'and pass debugging checks with segregated solvers with min-mod limiter.' + abs_zero = 1e-5 # limiter sensitive to roundoff + rel_err = 1e-4 # limiter sensitive to roundoff + cli_args = 'GlobalParams/advected_interp_method=min_mod Executioner/num_iterations=10 Outputs/file_base=min_mod_run' + recover = false # we don't support recovery for SIMPLE yet + max_threads = 1 # see libmesh issue #3808 + max_parallel = 4 # limiter sensitive to roundoff + [] + [quick] + type = 'Exodiff' + input = 'lid-driven-segregated.i' + exodiff = 'quick.e' + detail = 'and reach converged results with QUICK limiter.' + cli_args = 'GlobalParams/advected_interp_method=quick Outputs/file_base=quick' + heavy = true + recover = false # we don't support recovery for SIMPLE yet + max_threads = 1 # see libmesh issue #3808 + [] + [quick_run] + type = 'Exodiff' + input = 'lid-driven-segregated.i' + exodiff = 'quick_run.e' + detail = 'and pass debugging checks with segregated solvers with QUICK limiter.' + abs_zero = 1e-5 # limiter sensitive to roundoff + rel_err = 1e-4 # limiter sensitive to roundoff + cli_args = 'GlobalParams/advected_interp_method=quick Executioner/num_iterations=10 Outputs/file_base=quick_run' + recover = false # we don't support recovery for SIMPLE yet + max_threads = 1 # see libmesh issue #3808 + max_parallel = 4 # limiter sensitive to roundoff + [] + [] +[] diff --git a/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/gold/min_mod.e b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/gold/min_mod.e new file mode 100644 index 000000000000..530874a05e8e Binary files /dev/null and b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/gold/min_mod.e differ diff --git a/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/gold/quick.e b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/gold/quick.e new file mode 100644 index 000000000000..7ec8a0b259e6 Binary files /dev/null and b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/gold/quick.e differ diff --git a/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/gold/sou.e b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/gold/sou.e new file mode 100644 index 000000000000..ece57e16f3e1 Binary files /dev/null and b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/gold/sou.e differ diff --git a/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/gold/upwind.e b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/gold/upwind.e new file mode 100644 index 000000000000..c3f5fac22fe7 Binary files /dev/null and b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/gold/upwind.e differ diff --git a/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/gold/vanLeer.e b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/gold/vanLeer.e new file mode 100644 index 000000000000..5a1876f2d1f4 Binary files /dev/null and b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/gold/vanLeer.e differ diff --git a/modules/navier_stokes/test/tests/finite_volume/ins/advection-schemes/test.i b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/test.i similarity index 84% rename from modules/navier_stokes/test/tests/finite_volume/ins/advection-schemes/test.i rename to modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/test.i index dcac9d28c458..b71871c695d0 100644 --- a/modules/navier_stokes/test/tests/finite_volume/ins/advection-schemes/test.i +++ b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/test.i @@ -19,12 +19,15 @@ rho=1 [Variables] [u] type = INSFVVelocityVariable + two_term_boundary_expansion = false [] [v] type = INSFVVelocityVariable + two_term_boundary_expansion = false [] [pressure] type = INSFVPressureVariable + two_term_boundary_expansion = false [] [lambda] family = SCALAR @@ -128,13 +131,21 @@ rho=1 [] [Executioner] - type = Steady + type = Transient solve_type = 'NEWTON' petsc_options_iname = '-pc_type -pc_factor_shift_type' - petsc_options_value = 'lu NONZERO' - nl_rel_tol = 1e-12 + petsc_options_value = 'lu NONZERO' + dt = 0.1 + end_time = 5.0 + steady_state_detection = true + steady_state_tolerance = 1e-12 + nl_abs_tol = 1e-12 [] [Outputs] - exodus = true + [out] + type = Exodus + execute_on = 'final' + hide = 'lambda' + [] [] diff --git a/modules/navier_stokes/test/tests/finite_volume/ins/advection-schemes/tests b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/tests similarity index 62% rename from modules/navier_stokes/test/tests/finite_volume/ins/advection-schemes/tests rename to modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/tests index 00b1b96cc8e5..f814763a1f6a 100644 --- a/modules/navier_stokes/test/tests/finite_volume/ins/advection-schemes/tests +++ b/modules/navier_stokes/test/tests/finite_volume/limiters/lid-driven/tests @@ -1,14 +1,23 @@ [Tests] design = 'INSFVMomentumAdvection.md' - issues = '#20493' - [limiting_schemes] + issues = '#20493 #28891' + [advection_limiting_schemes] requirement = 'The system shall be able to perform a variety of limiting schemes when solving fluid flow equations. These schemes include' + [upwind] + type = Exodiff + input = 'test.i' + exodiff = upwind.e + cli_args = 'GlobalParams/advected_interp_method=upwind Outputs/file_base=upwind' + detail = 'first-order upwind' + recover = false # only using final for output in a transient + [] [sou] type = Exodiff input = 'test.i' exodiff = sou.e cli_args = 'GlobalParams/advected_interp_method=sou Outputs/file_base=sou' detail = 'second-order upwind' + recover = false # only using final for output in a transient [] [vanLeer] type = Exodiff @@ -16,6 +25,7 @@ exodiff = vanLeer.e cli_args = 'GlobalParams/advected_interp_method=vanLeer Outputs/file_base=vanLeer' detail = 'van Leer' + recover = false # only using final for output in a transient [] [min_mod] type = Exodiff @@ -23,6 +33,7 @@ exodiff = min_mod.e cli_args = 'GlobalParams/advected_interp_method=min_mod Outputs/file_base=min_mod' detail = 'min-mod' + recover = false # only using final for output in a transient [] [quick] type = Exodiff @@ -30,6 +41,7 @@ exodiff = quick.e cli_args = 'GlobalParams/advected_interp_method=quick Outputs/file_base=quick' detail = 'QUICK' + recover = false # only using final for output in a transient [] [] [] diff --git a/modules/navier_stokes/unit/src/TestFaceCenteredMapFunctor.C b/modules/navier_stokes/unit/src/TestFaceCenteredMapFunctor.C index d74a47274750..7c7def4f9744 100644 --- a/modules/navier_stokes/unit/src/TestFaceCenteredMapFunctor.C +++ b/modules/navier_stokes/unit/src/TestFaceCenteredMapFunctor.C @@ -85,8 +85,8 @@ TEST(FaceCenteredMapFunctorTest, testArgs) for (auto & fi : all_fi) { const auto & face_center = fi.faceCentroid(); - const auto face_arg = - Moose::FaceArg{&fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr}; + const auto face_arg = Moose::FaceArg{ + &fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr}; const auto result = u(face_arg, Moose::currentState()); @@ -128,8 +128,8 @@ TEST(FaceCenteredMapFunctorTest, testArgs) // Arguments for the simple error checks, we use the first face and the corresponding // owner element QGauss qrule(1, CONSTANT); - const auto face_arg = - Moose::FaceArg{&all_fi[0], Moose::FV::LimiterType::CentralDifference, true, false, nullptr}; + const auto face_arg = Moose::FaceArg{ + &all_fi[0], Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr}; const auto elem_arg = ElemArg{all_fi[0].elemPtr(), false}; const auto elem_qp_arg = ElemQpArg({all_fi[0].elemPtr(), 0, &qrule, Point(0)}); const auto elem_side_qp_arg = ElemSideQpArg({all_fi[0].elemPtr(), 0, 0, &qrule, Point(0)}); @@ -148,7 +148,7 @@ TEST(FaceCenteredMapFunctorTest, testArgs) try { unrestricted_error_test( - FaceArg{&all_fi[2], LimiterType::CentralDifference, true, false, nullptr}, + FaceArg{&all_fi[2], LimiterType::CentralDifference, true, false, nullptr, nullptr}, Moose::currentState()); EXPECT_TRUE(false); } @@ -162,8 +162,9 @@ TEST(FaceCenteredMapFunctorTest, testArgs) restricted_error_test(*mesh, {1}, "is_restricted"); try { - restricted_error_test(FaceArg{&all_fi[2], LimiterType::CentralDifference, true, false, nullptr}, - Moose::currentState()); + restricted_error_test( + FaceArg{&all_fi[2], LimiterType::CentralDifference, true, false, nullptr, nullptr}, + Moose::currentState()); EXPECT_TRUE(false); } catch (std::runtime_error & e) diff --git a/modules/navier_stokes/unit/src/TestReconstruction.C b/modules/navier_stokes/unit/src/TestReconstruction.C index 9c9d385980d9..65c381f07d80 100644 --- a/modules/navier_stokes/unit/src/TestReconstruction.C +++ b/modules/navier_stokes/unit/src/TestReconstruction.C @@ -115,8 +115,8 @@ testReconstruction(const Moose::CoordinateSystemType coord_type) { auto moukalled_reconstruct = [&fi](auto & functor, auto & container) { - auto face = - Moose::FaceArg({&fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr}); + auto face = Moose::FaceArg( + {&fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr}); const RealVectorValue uf(functor(face, Moose::currentState())); const Point surface_vector = fi.normal() * fi.faceArea(); auto product = (uf * fi.dCN()) * surface_vector; diff --git a/test/tests/fvkernels/dispersion-test/cartesian_advection.i b/test/tests/fvkernels/dispersion-test/cartesian_advection.i new file mode 100644 index 000000000000..24382b94e1bd --- /dev/null +++ b/test/tests/fvkernels/dispersion-test/cartesian_advection.i @@ -0,0 +1,76 @@ +[GlobalParams] + advected_interp_method = 'min_mod' #average upwind sou min_mod vanLeer quick venkatakrishnan skewness-corrected +[] + +[Mesh] + [gen] + type = GeneratedMeshGenerator + dim = 2 + xmin = 0 + xmax = 1 + ymin = 0 + ymax = 1 + nx = 11 + ny = 11 + [] +[] + +[Variables] + [scalar] + type = MooseVariableFVReal + two_term_boundary_expansion = false + [] +[] + +[FVKernels] + [time_derivative] + type = FVTimeKernel + variable = scalar + [] + [scalar_advection] + type = FVAdvection + variable = scalar + velocity = '1 1 0' + [] +[] + +[FVBCs] + [inflow_1] + type = FVDirichletBC + boundary = 'left' + value = '1' + variable = scalar + [] + [inflow_0] + type = FVDirichletBC + boundary = 'bottom' + value = '0' + variable = scalar + [] + [outflow] + type = FVConstantScalarOutflowBC + variable = scalar + velocity = '1 1 0' + boundary = 'right top' + [] +[] + +[Executioner] + type = Transient + solve_type = 'NEWTON' + petsc_options_iname = '-pc_type -pc_factor_shift_type' + petsc_options_value = 'lu NONZERO' + dt = 0.1 + end_time = 10.0 + steady_state_detection = true + steady_state_tolerance = 1e-12 + nl_abs_tol = 1e-12 + line_search = none +[] + +[Outputs] + [out] + type = Exodus + execute_on = 'final' + [] +[] diff --git a/test/tests/fvkernels/dispersion-test/gold/bottom_left_min_mod.e b/test/tests/fvkernels/dispersion-test/gold/bottom_left_min_mod.e new file mode 100644 index 000000000000..029c5197e5f6 Binary files /dev/null and b/test/tests/fvkernels/dispersion-test/gold/bottom_left_min_mod.e differ diff --git a/test/tests/fvkernels/dispersion-test/gold/bottom_left_quick.e b/test/tests/fvkernels/dispersion-test/gold/bottom_left_quick.e new file mode 100644 index 000000000000..222cb1b385d5 Binary files /dev/null and b/test/tests/fvkernels/dispersion-test/gold/bottom_left_quick.e differ diff --git a/test/tests/fvkernels/dispersion-test/gold/bottom_left_sou.e b/test/tests/fvkernels/dispersion-test/gold/bottom_left_sou.e new file mode 100644 index 000000000000..87477ee54f97 Binary files /dev/null and b/test/tests/fvkernels/dispersion-test/gold/bottom_left_sou.e differ diff --git a/test/tests/fvkernels/dispersion-test/gold/bottom_left_vanLeer.e b/test/tests/fvkernels/dispersion-test/gold/bottom_left_vanLeer.e new file mode 100644 index 000000000000..30d55a606056 Binary files /dev/null and b/test/tests/fvkernels/dispersion-test/gold/bottom_left_vanLeer.e differ diff --git a/test/tests/fvkernels/dispersion-test/gold/bottom_left_venkatakrishnan.e b/test/tests/fvkernels/dispersion-test/gold/bottom_left_venkatakrishnan.e new file mode 100644 index 000000000000..d9e50a1b3aea Binary files /dev/null and b/test/tests/fvkernels/dispersion-test/gold/bottom_left_venkatakrishnan.e differ diff --git a/test/tests/fvkernels/dispersion-test/gold/top_right_min_mod.e b/test/tests/fvkernels/dispersion-test/gold/top_right_min_mod.e new file mode 100644 index 000000000000..46157fc4aff4 Binary files /dev/null and b/test/tests/fvkernels/dispersion-test/gold/top_right_min_mod.e differ diff --git a/test/tests/fvkernels/dispersion-test/gold/top_right_quick.e b/test/tests/fvkernels/dispersion-test/gold/top_right_quick.e new file mode 100644 index 000000000000..369b47b35722 Binary files /dev/null and b/test/tests/fvkernels/dispersion-test/gold/top_right_quick.e differ diff --git a/test/tests/fvkernels/dispersion-test/gold/top_right_sou.e b/test/tests/fvkernels/dispersion-test/gold/top_right_sou.e new file mode 100644 index 000000000000..b0cb59df596f Binary files /dev/null and b/test/tests/fvkernels/dispersion-test/gold/top_right_sou.e differ diff --git a/test/tests/fvkernels/dispersion-test/gold/top_right_vanLeer.e b/test/tests/fvkernels/dispersion-test/gold/top_right_vanLeer.e new file mode 100644 index 000000000000..0cdcd42a1d33 Binary files /dev/null and b/test/tests/fvkernels/dispersion-test/gold/top_right_vanLeer.e differ diff --git a/test/tests/fvkernels/dispersion-test/gold/top_right_venkatakrishnan.e b/test/tests/fvkernels/dispersion-test/gold/top_right_venkatakrishnan.e new file mode 100644 index 000000000000..b93a3270635d Binary files /dev/null and b/test/tests/fvkernels/dispersion-test/gold/top_right_venkatakrishnan.e differ diff --git a/test/tests/fvkernels/dispersion-test/tests b/test/tests/fvkernels/dispersion-test/tests new file mode 100644 index 000000000000..723515b0e7d5 --- /dev/null +++ b/test/tests/fvkernels/dispersion-test/tests @@ -0,0 +1,92 @@ +[Tests] + design = 'FVAdvection.md' + issues = '#28891' + [bottom_left_limited_scalar_advection] + requirement = 'The system shall be able to perform a variety of limiting schemes when solving scalar transport equations in cartesian meshes with bottom-left advection. These schemes include' + [sou] + type = Exodiff + input = 'cartesian_advection.i' + exodiff = bottom_left_sou.e + cli_args = 'GlobalParams/advected_interp_method=sou Outputs/file_base=bottom_left_sou' + detail = 'second-order upwind' + recover = false # only using final for output in a transient + [] + [vanLeer] + type = Exodiff + input = 'cartesian_advection.i' + exodiff = bottom_left_vanLeer.e + cli_args = 'GlobalParams/advected_interp_method=vanLeer Outputs/file_base=bottom_left_vanLeer' + detail = 'van Leer' + recover = false # only using final for output in a transient + [] + [min_mod] + type = Exodiff + input = 'cartesian_advection.i' + exodiff = bottom_left_min_mod.e + cli_args = 'GlobalParams/advected_interp_method=min_mod Outputs/file_base=bottom_left_min_mod' + detail = 'min-mod' + recover = false # only using final for output in a transient + [] + [quick] + type = Exodiff + input = 'cartesian_advection.i' + exodiff = bottom_left_quick.e + cli_args = 'GlobalParams/advected_interp_method=quick Outputs/file_base=bottom_left_quick' + detail = 'QUICK' + recover = false # only using final for output in a transient + [] + [venkatakrishnan] + type = Exodiff + input = 'cartesian_advection.i' + exodiff = bottom_left_venkatakrishnan.e + cli_args = 'GlobalParams/advected_interp_method=venkatakrishnan Outputs/file_base=bottom_left_venkatakrishnan' + abs_zero = 1e-5 + detail = 'venkatakrishnan' + recover = false # only using final for output in a transient + [] + [] + [top_right_limited_scalar_advection] + requirement = 'The system shall be able to perform a variety of limiting schemes when solving scalar transport equations in cartesian meshes with top-right advection. These schemes include' + [sou] + type = Exodiff + input = 'cartesian_advection.i' + exodiff = top_right_sou.e + cli_args = "GlobalParams/advected_interp_method=sou Outputs/file_base=top_right_sou FVKernels/scalar_advection/velocity='-1 -1 0' FVBCs/inflow_1/boundary='top' FVBCs/inflow_0/boundary='right' FVBCs/outflow/boundary='left bottom' FVBCs/outflow/velocity='-1 -1 0'" + detail = 'second-order upwind' + recover = false # only using final for output in a transient + [] + [vanLeer] + type = Exodiff + input = 'cartesian_advection.i' + exodiff = top_right_vanLeer.e + cli_args = "GlobalParams/advected_interp_method=vanLeer Outputs/file_base=top_right_vanLeer FVKernels/scalar_advection/velocity='-1 -1 0' FVBCs/inflow_1/boundary='top' FVBCs/inflow_0/boundary='right' FVBCs/outflow/boundary='left bottom' FVBCs/outflow/velocity='-1 -1 0'" + detail = 'van Leer' + recover = false # only using final for output in a transient + [] + [min_mod] + type = Exodiff + input = 'cartesian_advection.i' + exodiff = top_right_min_mod.e + cli_args = "GlobalParams/advected_interp_method=min_mod Outputs/file_base=top_right_min_mod FVKernels/scalar_advection/velocity='-1 -1 0' FVBCs/inflow_1/boundary='top' FVBCs/inflow_0/boundary='right' FVBCs/outflow/boundary='left bottom' FVBCs/outflow/velocity='-1 -1 0'" + detail = 'min-mod' + recover = false # only using final for output in a transient + [] + [quick] + type = Exodiff + input = 'cartesian_advection.i' + exodiff = top_right_quick.e + cli_args = "GlobalParams/advected_interp_method=quick Outputs/file_base=top_right_quick FVKernels/scalar_advection/velocity='-1 -1 0' FVBCs/inflow_1/boundary='top' FVBCs/inflow_0/boundary='right' FVBCs/outflow/boundary='left bottom' FVBCs/outflow/velocity='-1 -1 0'" + detail = 'QUICK' + recover = false # only using final for output in a transient + [] + [venkatakrishnan] + type = Exodiff + input = 'cartesian_advection.i' + exodiff = top_right_venkatakrishnan.e + cli_args = "GlobalParams/advected_interp_method=venkatakrishnan Outputs/file_base=top_right_venkatakrishnan FVKernels/scalar_advection/velocity='-1 -1 0' FVBCs/inflow_1/boundary='top' FVBCs/inflow_0/boundary='right' FVBCs/outflow/boundary='left bottom' FVBCs/outflow/velocity='-1 -1 0'" + abs_zero = 1e-5 + detail = 'venkatakrishnan' + recover = false # only using final for output in a transient + [] + [] +[] diff --git a/test/tests/fvkernels/mms/advective-outflow/advection-outflow.i b/test/tests/fvkernels/mms/advective-outflow/advection-outflow.i index 07d9208e5139..378f33da74ba 100644 --- a/test/tests/fvkernels/mms/advective-outflow/advection-outflow.i +++ b/test/tests/fvkernels/mms/advective-outflow/advection-outflow.i @@ -93,14 +93,20 @@ a=1.1 [] [Executioner] - type = Steady + type = Transient + dt = 1 + end_time = 10 + steady_state_tolerance = 1e-10 solve_type = 'NEWTON' petsc_options_iname = '-pc_type -pc_factor_shift_type -pc_factor_mat_solver_type' petsc_options_value = 'lu NONZERO mumps' [] [Outputs] - exodus = true + [out] + type = Exodus + execute_on = 'final' + [] csv = true [] diff --git a/test/tests/fvkernels/mms/advective-outflow/test.py b/test/tests/fvkernels/mms/advective-outflow/test.py index 33eb5133f465..535d77e90140 100644 --- a/test/tests/fvkernels/mms/advective-outflow/test.py +++ b/test/tests/fvkernels/mms/advective-outflow/test.py @@ -1,6 +1,8 @@ import mms import unittest -from mooseutils import fuzzyAbsoluteEqual + +def bottom_bound(value, order, tolerance): + return value > (order - tolerance) class TestOutflow(unittest.TestCase): def test(self): @@ -15,9 +17,9 @@ def test(self): fig.save('outflow.png') for label,value in fig.label_to_slope.items(): if label == 'L2u': - self.assertTrue(fuzzyAbsoluteEqual(value, 1., .05)) + self.assertTrue(bottom_bound(value, 1., .05)) else: - self.assertTrue(fuzzyAbsoluteEqual(value, 2., .05)) + self.assertTrue(bottom_bound(value, 2., .05)) class TestOutflowMinMod(unittest.TestCase): def test(self): @@ -39,9 +41,9 @@ def test(self): fig.save('outflow-min-mod.png') for label,value in fig.label_to_slope.items(): if label == 'L2u': - self.assertTrue(fuzzyAbsoluteEqual(value, 1.5, .05)) + self.assertTrue(bottom_bound(value, 1.5, .05)) else: - self.assertTrue(fuzzyAbsoluteEqual(value, 2., .05)) + self.assertTrue(bottom_bound(value, 2., .05)) class TestExtrapolation(unittest.TestCase): def test(self): @@ -56,9 +58,9 @@ def test(self): fig.save('extrapolation.png') for label,value in fig.label_to_slope.items(): if label == 'L2u': - self.assertTrue(fuzzyAbsoluteEqual(value, 1., .05)) + self.assertTrue(bottom_bound(value, 1., .05)) else: - self.assertTrue(fuzzyAbsoluteEqual(value, 2., .05)) + self.assertTrue(bottom_bound(value, 2., .05)) class UpwindLimiter(unittest.TestCase): def test(self): @@ -72,7 +74,7 @@ def test(self): slope_precision=1) fig.save('upwind-limiter.png') for label,value in fig.label_to_slope.items(): - self.assertTrue(fuzzyAbsoluteEqual(value, 1., .05)) + self.assertTrue(bottom_bound(value, 1., .05)) class CentralDifferenceLimiter(unittest.TestCase): def test(self): @@ -86,7 +88,7 @@ def test(self): slope_precision=1) fig.save('cd-limiter.png') for label,value in fig.label_to_slope.items(): - self.assertTrue(fuzzyAbsoluteEqual(value, 2., .05)) + self.assertTrue(bottom_bound(value, 2., .05)) class VanLeerLimiter(unittest.TestCase): def test(self): @@ -100,7 +102,7 @@ def test(self): slope_precision=1) fig.save('vanLeer-limiter.png') for label,value in fig.label_to_slope.items(): - self.assertTrue(fuzzyAbsoluteEqual(value, 2., .05)) + self.assertTrue(bottom_bound(value, 2., .05)) class MinModLimiter(unittest.TestCase): def test(self): @@ -114,7 +116,7 @@ def test(self): slope_precision=1) fig.save('min-mod-limiter.png') for label,value in fig.label_to_slope.items(): - self.assertTrue(fuzzyAbsoluteEqual(value, 2., .05)) + self.assertTrue(bottom_bound(value, 2., .05)) class SOULimiter(unittest.TestCase): def test(self): @@ -128,7 +130,7 @@ def test(self): slope_precision=1) fig.save('sou-limiter.png') for label,value in fig.label_to_slope.items(): - self.assertTrue(fuzzyAbsoluteEqual(value, 2., .05)) + self.assertTrue(bottom_bound(value, 2., .05)) class QUICKLimiter(unittest.TestCase): def test(self): @@ -142,7 +144,7 @@ def test(self): slope_precision=1) fig.save('quick-limiter.png') for label,value in fig.label_to_slope.items(): - self.assertTrue(fuzzyAbsoluteEqual(value, 2., .05)) + self.assertTrue(bottom_bound(value, 2., .05)) class KTLimitedCD(unittest.TestCase): def test(self): @@ -156,7 +158,7 @@ def test(self): slope_precision=1) fig.save('kt-cd-limiter.png') for key,value in fig.label_to_slope.items(): - self.assertTrue(fuzzyAbsoluteEqual(value, 2., .05)) + self.assertTrue(bottom_bound(value, 2., .05)) print("%s slope, %f" % (key, value)) class KTLimitedUpwind(unittest.TestCase): @@ -171,7 +173,7 @@ def test(self): slope_precision=1) fig.save('kt-upwind-limiter.png') for key,value in fig.label_to_slope.items(): - self.assertTrue(fuzzyAbsoluteEqual(value, 1., .05)) + self.assertTrue(bottom_bound(value, 1., .05)) print("%s slope, %f" % (key, value)) class KTLimitedVanLeer(unittest.TestCase): @@ -186,5 +188,5 @@ def test(self): slope_precision=1) fig.save('kt-van-leer-limiter.png') for key,value in fig.label_to_slope.items(): - self.assertTrue(fuzzyAbsoluteEqual(value, 2.5, .05)) + self.assertTrue(bottom_bound(value, 2.5, .05)) print("%s slope, %f" % (key, value)) diff --git a/test/tests/fvkernels/vector-interpolation/gold/vanLeer.e b/test/tests/fvkernels/vector-interpolation/gold/vanLeer.e index ec501e9ade36..adf6790caef8 100644 Binary files a/test/tests/fvkernels/vector-interpolation/gold/vanLeer.e and b/test/tests/fvkernels/vector-interpolation/gold/vanLeer.e differ diff --git a/unit/src/ContainerFunctors.C b/unit/src/ContainerFunctors.C index 743cdbcd524d..a24ce11b647e 100644 --- a/unit/src/ContainerFunctors.C +++ b/unit/src/ContainerFunctors.C @@ -92,7 +92,7 @@ TEST(ContainerFunctors, Test) for (const auto * const face : faces) { const auto elem = vector_functor.isInternalFace(*face) ? nullptr : face->elemPtr(); - const auto face_arg = Moose::FaceArg({face, limiter_type, true, false, elem}); + const auto face_arg = Moose::FaceArg({face, limiter_type, true, false, elem, nullptr}); const auto current_time = Moose::currentState(); EXPECT_TRUE(vector_functor(face_arg, current_time)[0] == 1.); EXPECT_TRUE(array_functor(face_arg, current_time)[0] == 1.); diff --git a/unit/src/MooseFunctorTest.C b/unit/src/MooseFunctorTest.C index 3d74948d36e7..4e6e6fdff87a 100644 --- a/unit/src/MooseFunctorTest.C +++ b/unit/src/MooseFunctorTest.C @@ -132,7 +132,7 @@ TEST(MooseFunctorTest, testArgs) QGauss qrule(1, CONSTANT); auto elem_arg = ElemArg{elem.get(), false}; - auto face = FaceArg({&fi, LimiterType::CentralDifference, true, false, nullptr}); + auto face = FaceArg({&fi, LimiterType::CentralDifference, true, false, nullptr, nullptr}); auto elem_qp = ElemQpArg({elem.get(), 0, &qrule, Point(0)}); auto elem_side_qp = ElemSideQpArg({elem.get(), 0, 0, &qrule, Point(0)}); auto elem_point = ElemPointArg({elem.get(), Point(0), false}); @@ -355,7 +355,8 @@ TEST(MooseFunctorTest, testArgs) if (!mesh_fi.neighborPtr()) continue; - auto vec_face_arg = FaceArg({&mesh_fi, LimiterType::CentralDifference, true, false, nullptr}); + auto vec_face_arg = + FaceArg({&mesh_fi, LimiterType::CentralDifference, true, false, nullptr, nullptr}); const auto vec_elem_arg = vec_face_arg.makeElem(); const auto vec_neighbor_arg = vec_face_arg.makeNeighbor(); zero_gradient_and_grad_dot_test(vec_comp, vec_elem_arg); diff --git a/unit/src/ParsedFunctionTest.C b/unit/src/ParsedFunctionTest.C index 9fd1f5e31e68..928152efe075 100644 --- a/unit/src/ParsedFunctionTest.C +++ b/unit/src/ParsedFunctionTest.C @@ -100,7 +100,8 @@ TEST_F(ParsedFunctionTest, basicConstructor) // Test face overloads _mesh->buildFiniteVolumeInfo(); const FaceInfo * const fi = _mesh->faceInfo(elem, side); - auto face = Moose::FaceArg({fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr}); + auto face = Moose::FaceArg( + {fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr}); f_traditional = f.value(0, fi->faceCentroid()); f_functor = f_wrapped(face, Moose::currentState()); gradient_traditional = f.gradient(0, fi->faceCentroid());