diff --git a/framework/include/utils/MathFVUtils.h b/framework/include/utils/MathFVUtils.h index 8bcdee962ea1..bf657855f675 100644 --- a/framework/include/utils/MathFVUtils.h +++ b/framework/include/utils/MathFVUtils.h @@ -90,6 +90,23 @@ faceArgSubdomains(const SubdomainRestrictable & obj, const FaceInfo & fi) makeSidedFace(obj, fi.neighborPtr(), fi).sub_id); } +inline FaceArg +makeFace(const FaceInfo & fi, + const LimiterType limiter_type, + const bool elem_is_upwind, + const std::pair & subs, + const bool skewness_correction = false, + const bool apply_gradient_correction = false) +{ + return {&fi, + limiter_type, + elem_is_upwind, + skewness_correction, + apply_gradient_correction, + subs.first, + subs.second}; +} + /** * Make a functor face argument with a central differencing limiter, e.g. compose a face argument * that will tell functors to perform (possibly skew-corrected) linear interpolations from cell @@ -111,13 +128,12 @@ makeCDFace(const FaceInfo & fi, const bool skewness_correction = false, const bool apply_gradient_correction = false) { - return {&fi, - Moose::FV::LimiterType::CentralDifference, - true, - skewness_correction, - apply_gradient_correction, - subs.first, - subs.second}; + return makeFace(fi, + LimiterType::CentralDifference, + true, + subs, + skewness_correction, + apply_gradient_correction); } /** @@ -482,11 +498,6 @@ interpCoeffs(const Limiter & limiter, return std::make_pair(1. - g, g); } -template -struct GradientType -{ -}; - /** * Interpolates with a limiter */ @@ -598,6 +609,105 @@ interpolate(const Limiter & limiter, return ret; } +/** + * 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 + */ +template ::value>::type> +std::pair, std::pair> +interpCoeffsAndAdvected(const FunctorBase & functor, const FaceArg & face) +{ + typedef typename FunctorBase::GradientType GradientType; + static const GradientType zero(0); + + mooseAssert(face.fi, "this must be non-null"); + const auto limiter = Limiter::value_type>::build(face.limiter_type); + + const auto upwind_arg = face.elem_is_upwind ? face.elemFromFace() : face.neighborFromFace(); + const auto downwind_arg = face.elem_is_upwind ? face.neighborFromFace() : face.elemFromFace(); + auto phi_upwind = functor(upwind_arg); + auto phi_downwind = functor(downwind_arg); + + std::pair interp_coeffs; + 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); + else + { + const auto grad_phi_upwind = functor.gradient(upwind_arg); + interp_coeffs = interpCoeffs( + *limiter, phi_upwind, phi_downwind, &grad_phi_upwind, *face.fi, face.elem_is_upwind); + } + + if (face.elem_is_upwind) + return std::make_pair(std::move(interp_coeffs), + std::make_pair(std::move(phi_upwind), std::move(phi_downwind))); + else + return std::make_pair( + std::make_pair(std::move(interp_coeffs.second), std::move(interp_coeffs.first)), + std::make_pair(std::move(phi_downwind), std::move(phi_upwind))); +} + +template ::value>::type> +T +interpolate(const FunctorBase & functor, const FaceArg & face) +{ + const auto [interp_coeffs, advected] = interpCoeffsAndAdvected(functor, face); + return interp_coeffs.first * advected.first + interp_coeffs.second * advected.second; +} + +template +VectorValue +interpolate(const FunctorBase> & functor, const FaceArg & face) +{ + static const VectorValue grad_zero(0); + + mooseAssert(face.fi, "this must be non-null"); + const auto limiter = Limiter::value_type>::build(face.limiter_type); + + const auto upwind_arg = face.elem_is_upwind ? face.elemFromFace() : face.neighborFromFace(); + const auto downwind_arg = face.elem_is_upwind ? face.neighborFromFace() : face.elemFromFace(); + auto phi_upwind = functor(upwind_arg); + auto phi_downwind = functor(downwind_arg); + + VectorValue ret; + T coeff_upwind, coeff_downwind; + + if (face.limiter_type == LimiterType::Upwind || + face.limiter_type == LimiterType::CentralDifference) + { + for (unsigned int i = 0; i < LIBMESH_DIM; ++i) + { + const auto &component_upwind = phi_upwind(i), component_downwind = phi_downwind(i); + std::tie(coeff_upwind, coeff_downwind) = interpCoeffs(*limiter, + component_upwind, + component_downwind, + &grad_zero, + *face.fi, + face.elem_is_upwind); + ret(i) = coeff_upwind * component_upwind + coeff_downwind * component_downwind; + } + } + else + { + const auto grad_phi_upwind = functor.gradient(upwind_arg); + 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; + } + } + + return ret; +} + /** * Return whether the supplied face is on a boundary of the \p object's execution */ diff --git a/framework/include/utils/PiecewiseByBlockLambdaFunctor.h b/framework/include/utils/PiecewiseByBlockLambdaFunctor.h index cd8b918e382f..aa62935ea5bc 100644 --- a/framework/include/utils/PiecewiseByBlockLambdaFunctor.h +++ b/framework/include/utils/PiecewiseByBlockLambdaFunctor.h @@ -78,6 +78,8 @@ class PiecewiseByBlockLambdaFunctor : public Moose::FunctorBase using Moose::FunctorBase::evaluateGradient; GradientType evaluateGradient(const Moose::ElemArg & elem_arg, unsigned int) const override; GradientType evaluateGradient(const Moose::FaceArg & face_arg, unsigned int) const override; + GradientType evaluateGradient(const Moose::ElemFromFaceArg & elem_from_face_arg, + unsigned int) const override; private: /** @@ -228,36 +230,12 @@ PiecewiseByBlockLambdaFunctor::evaluate(const Moose::SingleSidedFaceArg & fac template typename PiecewiseByBlockLambdaFunctor::ValueType -PiecewiseByBlockLambdaFunctor::evaluate(const Moose::FaceArg & face, unsigned int state) const +PiecewiseByBlockLambdaFunctor::evaluate(const Moose::FaceArg & face, + unsigned int libmesh_dbg_var(state)) const { using namespace Moose::FV; - - mooseAssert(face.fi, - "We must have a non-null face_info in order to prepare our ElemFromFace tuples"); - static const GradientType example_gradient; - - switch (face.limiter_type) - { - case LimiterType::CentralDifference: - case LimiterType::Upwind: - { - const auto elem_from_face = face.elemFromFace(); - const auto neighbor_from_face = face.neighborFromFace(); - const auto & upwind_elem = face.elem_is_upwind ? elem_from_face : neighbor_from_face; - const auto & downwind_elem = face.elem_is_upwind ? neighbor_from_face : elem_from_face; - auto limiter_obj = - Limiter::value_type>::build(face.limiter_type); - return interpolate(*limiter_obj, - evaluate(upwind_elem, state), - evaluate(downwind_elem, state), - &example_gradient, - *face.fi, - face.elem_is_upwind); - } - - default: - mooseError("Unsupported limiter type in PiecewiseByBlockLambdaFunctor::evaluate(FaceArg)"); - } + mooseAssert(state == 0, "Only current time state supported."); + return interpolate(*this, face); } template @@ -301,3 +279,15 @@ PiecewiseByBlockLambdaFunctor::evaluateGradient(const Moose::FaceArg & face_a { return Moose::FV::greenGaussGradient(face_arg, *this, true, _mesh); } + +template +typename PiecewiseByBlockLambdaFunctor::GradientType +PiecewiseByBlockLambdaFunctor::evaluateGradient( + const Moose::ElemFromFaceArg & elem_from_face_arg, unsigned int) const +{ + const auto elem_arg = elem_from_face_arg.makeElem(); + if (elem_arg.elem) + return Moose::FV::greenGaussGradient(elem_arg, *this, true, _mesh); + else + return GradientType(0); +} diff --git a/framework/src/limiters/Limiter.C b/framework/src/limiters/Limiter.C index d8660c3bafa5..da7661ad8d94 100644 --- a/framework/src/limiters/Limiter.C +++ b/framework/src/limiters/Limiter.C @@ -66,6 +66,18 @@ limiterType(const InterpMethod interp_method) case InterpMethod::Upwind: return LimiterType::Upwind; + case InterpMethod::VanLeer: + return LimiterType::VanLeer; + + case InterpMethod::MinMod: + return LimiterType::MinMod; + + case InterpMethod::SOU: + return LimiterType::SOU; + + case InterpMethod::QUICK: + return LimiterType::QUICK; + default: mooseError("Unrecognized interpolation method type."); } diff --git a/modules/navier_stokes/include/fvkernels/INSFVMomentumAdvection.h b/modules/navier_stokes/include/fvkernels/INSFVMomentumAdvection.h index 1064b0538aee..628d2821785d 100644 --- a/modules/navier_stokes/include/fvkernels/INSFVMomentumAdvection.h +++ b/modules/navier_stokes/include/fvkernels/INSFVMomentumAdvection.h @@ -11,6 +11,7 @@ #include "INSFVAdvectionKernel.h" #include "INSFVMomentumResidualObject.h" +#include "PiecewiseByBlockLambdaFunctor.h" /** * An advection kernel that implements interpolation schemes specific to Navier-Stokes flow @@ -30,6 +31,9 @@ class INSFVMomentumAdvection : public INSFVAdvectionKernel, public INSFVMomentum /// Density const Moose::Functor & _rho; + /// Our local momentum functor + const PiecewiseByBlockLambdaFunctor _rho_u; + /// The a coefficient for the element ADReal _ae = 0; diff --git a/modules/navier_stokes/src/fvkernels/INSFVAdvectionKernel.C b/modules/navier_stokes/src/fvkernels/INSFVAdvectionKernel.C index e8a9a1422956..7f49547664f4 100644 --- a/modules/navier_stokes/src/fvkernels/INSFVAdvectionKernel.C +++ b/modules/navier_stokes/src/fvkernels/INSFVAdvectionKernel.C @@ -10,12 +10,14 @@ #include "INSFVAdvectionKernel.h" #include "NS.h" #include "MooseVariableFV.h" +#include "RelationshipManager.h" InputParameters INSFVAdvectionKernel::validParams() { InputParameters params = FVFluxKernel::validParams(); - MooseEnum advected_interp_method("average upwind skewness-corrected", "upwind"); + MooseEnum advected_interp_method("average upwind sou min_mod vanLeer quick skewness-corrected", + "upwind"); params.addParam( "advected_interp_method", advected_interp_method, @@ -53,8 +55,46 @@ INSFVAdvectionKernel::INSFVAdvectionKernel(const InputParameters & params) else if (advected_interp_method == "upwind") _advected_interp_method = InterpMethod::Upwind; else - mooseError("Unrecognized interpolation type ", - static_cast(advected_interp_method)); + { + if (advected_interp_method == "sou") + _advected_interp_method = InterpMethod::SOU; + else if (advected_interp_method == "min_mod") + _advected_interp_method = InterpMethod::MinMod; + else if (advected_interp_method == "vanLeer") + _advected_interp_method = InterpMethod::VanLeer; + else if (advected_interp_method == "quick") + _advected_interp_method = InterpMethod::QUICK; + else + mooseError("Unrecognized interpolation type ", + static_cast(advected_interp_method)); + + if (_tid == 0) + { + auto & factory = _app.getFactory(); + + auto rm_params = factory.getValidParams("ElementSideNeighborLayers"); + + rm_params.set("for_whom") = name(); + rm_params.set("mesh") = &const_cast(_mesh); + rm_params.set("rm_type") = + Moose::RelationshipManagerType::GEOMETRIC | Moose::RelationshipManagerType::ALGEBRAIC | + Moose::RelationshipManagerType::COUPLING; + FVKernel::setRMParams( + _pars, + rm_params, + std::max((unsigned short)(3), _pars.get("ghost_layers"))); + mooseAssert(rm_params.areAllRequiredParamsValid(), + "All relationship manager parameters should be valid."); + + auto rm_obj = factory.create( + "ElementSideNeighborLayers", name() + "_skew_correction", rm_params); + + // Delete the resources created on behalf of the RM if it ends up not being added to the + // App. + if (!_app.addRelationshipManager(rm_obj)) + factory.releaseSharedObjects(*rm_obj); + } + } const auto & velocity_interp_method = getParam("velocity_interp_method"); if (velocity_interp_method == "average") diff --git a/modules/navier_stokes/src/fvkernels/INSFVEnergyAdvection.C b/modules/navier_stokes/src/fvkernels/INSFVEnergyAdvection.C index 1bbac4425cc1..b09e051b7442 100644 --- a/modules/navier_stokes/src/fvkernels/INSFVEnergyAdvection.C +++ b/modules/navier_stokes/src/fvkernels/INSFVEnergyAdvection.C @@ -35,18 +35,12 @@ INSFVEnergyAdvection::INSFVEnergyAdvection(const InputParameters & params) ADReal INSFVEnergyAdvection::computeQpResidual() { - ADReal adv_quant_interface; - - const auto elem_face = elemFromFace(); - const auto neighbor_face = neighborFromFace(); - const auto v = _rc_vel_provider.getVelocity(_velocity_interp_method, *_face_info, _tid); - Moose::FV::interpolate(_advected_interp_method, - adv_quant_interface, - _adv_quant(elem_face), - _adv_quant(neighbor_face), - v, - *_face_info, - true); + const auto adv_quant_interface = + Moose::FV::interpolate(_adv_quant, + Moose::FV::makeFace(*_face_info, + limiterType(_advected_interp_method), + MetaPhysicL::raw_value(v) * _normal > 0, + faceArgSubdomains())); return _normal * v * adv_quant_interface; } diff --git a/modules/navier_stokes/src/fvkernels/INSFVMassAdvection.C b/modules/navier_stokes/src/fvkernels/INSFVMassAdvection.C index 8b96ff546ee9..1af153d6a0c4 100644 --- a/modules/navier_stokes/src/fvkernels/INSFVMassAdvection.C +++ b/modules/navier_stokes/src/fvkernels/INSFVMassAdvection.C @@ -34,18 +34,12 @@ INSFVMassAdvection::INSFVMassAdvection(const InputParameters & params) ADReal INSFVMassAdvection::computeQpResidual() { - ADReal rho_interface; - - const auto elem_face = elemFromFace(); - const auto neighbor_face = neighborFromFace(); - const auto v = _rc_vel_provider.getVelocity(_velocity_interp_method, *_face_info, _tid); - Moose::FV::interpolate(_advected_interp_method, - rho_interface, - _rho(elem_face), - _rho(neighbor_face), - v, - *_face_info, - true); + const auto rho_interface = + Moose::FV::interpolate(_rho, + Moose::FV::makeFace(*_face_info, + limiterType(_advected_interp_method), + MetaPhysicL::raw_value(v) * _normal > 0, + faceArgSubdomains())); return _normal * v * rho_interface; } diff --git a/modules/navier_stokes/src/fvkernels/INSFVMomentumAdvection.C b/modules/navier_stokes/src/fvkernels/INSFVMomentumAdvection.C index 2596e3267a48..416c90711161 100644 --- a/modules/navier_stokes/src/fvkernels/INSFVMomentumAdvection.C +++ b/modules/navier_stokes/src/fvkernels/INSFVMomentumAdvection.C @@ -27,7 +27,13 @@ INSFVMomentumAdvection::validParams() INSFVMomentumAdvection::INSFVMomentumAdvection(const InputParameters & params) : INSFVAdvectionKernel(params), INSFVMomentumResidualObject(*this), - _rho(getFunctor(NS::density)) + _rho(getFunctor(NS::density)), + _rho_u( + "rho_u", + [this](const auto & r, const auto & t) -> ADReal { return _rho(r, t) * _var(r, t); }, + std::set({EXEC_ALWAYS}), + _mesh, + this->blockIDs()) { #ifndef MOOSE_GLOBAL_AD_INDEXING mooseError("INSFV is not supported by local AD indexing. In order to use INSFV, please run the " @@ -39,16 +45,27 @@ INSFVMomentumAdvection::INSFVMomentumAdvection(const InputParameters & params) ADReal INSFVMomentumAdvection::computeQpResidual() { + using namespace Moose::FV; + + const auto v = _rc_vel_provider.getVelocity(_velocity_interp_method, *_face_info, _tid); + const auto [interp_coeffs, advected] = + interpCoeffsAndAdvected(_rho_u, + makeFace(*_face_info, + limiterType(_advected_interp_method), + MetaPhysicL::raw_value(v) * _normal > 0, + faceArgSubdomains())); + const auto elem_face = elemFromFace(); const auto neighbor_face = neighborFromFace(); - const auto v = _rc_vel_provider.getVelocity(_velocity_interp_method, *_face_info, _tid); - const auto interp_coeffs = Moose::FV::interpCoeffs(_advected_interp_method, *_face_info, true, v); - _ae = _normal * v * _rho(elem_face) * interp_coeffs.first; + const auto rho_elem = _rho(elem_face), rho_neighbor = _rho(neighbor_face); + const auto var_elem = advected.first / rho_elem, var_neighbor = advected.second / rho_neighbor; + + _ae = _normal * v * rho_elem * interp_coeffs.first; // Minus sign because we apply a minus sign to the residual in computeResidual - _an = -_normal * v * _rho(neighbor_face) * interp_coeffs.second; + _an = -_normal * v * rho_neighbor * interp_coeffs.second; - return _ae * _var(elem_face) - _an * _var(neighbor_face); + return _ae * var_elem - _an * var_neighbor; } void