diff --git a/modules/navier_stokes/include/fvkernels/INSFVMomentumAdvection.h b/modules/navier_stokes/include/fvkernels/INSFVMomentumAdvection.h index 6f76af8e44a0..ba13d61daa0d 100644 --- a/modules/navier_stokes/include/fvkernels/INSFVMomentumAdvection.h +++ b/modules/navier_stokes/include/fvkernels/INSFVMomentumAdvection.h @@ -46,25 +46,22 @@ class INSFVMomentumAdvection : public FVMatAdvection void residualSetup() override final { clearRCCoeffs(); } void jacobianSetup() override final { clearRCCoeffs(); } - /// The dynamic viscosity on the FaceInfo elem - const ADMaterialProperty & _mu_elem; - - /// The dynamic viscosity on the FaceInfo neighbor - const ADMaterialProperty & _mu_neighbor; + /// The dynamic viscosity + const FunctorMaterialProperty & _mu; /** * Returns the Rhie-Chow 'a' coefficient for the requested element \p elem * @param elem The elem to get the Rhie-Chow coefficient for * @param mu The dynamic viscosity */ - const VectorValue & rcCoeff(const Elem & elem, const ADReal & mu) const; + const VectorValue & rcCoeff(const Elem & elem) const; /** * method for computing the Rhie-Chow 'a' coefficients for the given elem \p elem * @param elem The elem to compute the Rhie-Chow coefficient for * @param mu The dynamic viscosity */ - virtual VectorValue coeffCalculator(const Elem & elem, const ADReal & mu) const; + virtual VectorValue coeffCalculator(const Elem & elem) const; /** * Clear the RC 'a' coefficient cache diff --git a/modules/navier_stokes/include/fvkernels/PINSFVMomentumAdvection.h b/modules/navier_stokes/include/fvkernels/PINSFVMomentumAdvection.h index 4243af3f6145..8b8d89fa5e17 100644 --- a/modules/navier_stokes/include/fvkernels/PINSFVMomentumAdvection.h +++ b/modules/navier_stokes/include/fvkernels/PINSFVMomentumAdvection.h @@ -31,7 +31,7 @@ class PINSFVMomentumAdvection : public INSFVMomentumAdvection const ADRealVectorValue & neighbor_v) override; virtual ADReal computeQpResidual() override; - VectorValue coeffCalculator(const Elem & elem, const ADReal & mu) const override; + VectorValue coeffCalculator(const Elem & elem) const override; /// porosity variable to compute gradients const MooseVariableFV * const _eps_var; diff --git a/modules/navier_stokes/src/fvkernels/INSFVMomentumAdvection.C b/modules/navier_stokes/src/fvkernels/INSFVMomentumAdvection.C index 223148fbfe19..fc74744f279e 100644 --- a/modules/navier_stokes/src/fvkernels/INSFVMomentumAdvection.C +++ b/modules/navier_stokes/src/fvkernels/INSFVMomentumAdvection.C @@ -67,8 +67,7 @@ INSFVMomentumAdvection::validParams() INSFVMomentumAdvection::INSFVMomentumAdvection(const InputParameters & params) : FVMatAdvection(params), - _mu_elem(getADMaterialProperty("mu")), - _mu_neighbor(getNeighborADMaterialProperty("mu")), + _mu(getFunctorMaterialProperty("mu")), _p_var(dynamic_cast(getFieldVar(NS::pressure, 0))), _u_var(dynamic_cast(getFieldVar("u", 0))), _v_var(params.isParamValid("v") @@ -233,7 +232,7 @@ INSFVMomentumAdvection::skipForBoundary(const FaceInfo & fi) const } const VectorValue & -INSFVMomentumAdvection::rcCoeff(const Elem & elem, const ADReal & mu) const +INSFVMomentumAdvection::rcCoeff(const Elem & elem) const { auto it = _rc_a_coeffs.find(&_app); mooseAssert(it != _rc_a_coeffs.end(), @@ -251,7 +250,7 @@ INSFVMomentumAdvection::rcCoeff(const Elem & elem, const ADReal & mu) const // Returns a pair with the first being an iterator pointing to the key-value pair and the second a // boolean denoting whether a new insertion took place - auto emplace_ret = my_map.emplace(&elem, coeffCalculator(elem, mu)); + auto emplace_ret = my_map.emplace(&elem, coeffCalculator(elem)); mooseAssert(emplace_ret.second, "We should have inserted a new key-value pair"); @@ -260,7 +259,7 @@ INSFVMomentumAdvection::rcCoeff(const Elem & elem, const ADReal & mu) const #ifdef MOOSE_GLOBAL_AD_INDEXING VectorValue -INSFVMomentumAdvection::coeffCalculator(const Elem & elem, const ADReal & mu) const +INSFVMomentumAdvection::coeffCalculator(const Elem & elem) const { // these coefficients arise from simple control volume balances of advection and diffusion. These // coefficients are the linear coefficients associated with the centroid of the control volume. @@ -296,7 +295,6 @@ INSFVMomentumAdvection::coeffCalculator(const Elem & elem, const ADReal & mu) co auto action_functor = [&coeff, &elem_velocity, - &mu, #ifndef NDEBUG &elem, #endif @@ -321,6 +319,8 @@ INSFVMomentumAdvection::coeffCalculator(const Elem & elem, const ADReal & mu) co "Let's make sure our normal is what we think it is"); #endif + const auto elem_mu = _mu(&elem); + // Unless specified otherwise, "elem" here refers to the element we're computing the // Rhie-Chow coefficient for. "neighbor" is the element across the current FaceInfo (fi) // face from the Rhie-Chow element @@ -335,7 +335,7 @@ INSFVMomentumAdvection::coeffCalculator(const Elem & elem, const ADReal & mu) co { // Need to account for viscous shear stress from wall for (const auto i : make_range(_dim)) - coeff(i) += mu * surface_vector.norm() / + coeff(i) += elem_mu * surface_vector.norm() / std::abs((fi->faceCentroid() - rc_centroid) * normal) * (1 - normal(i) * normal(i)); @@ -369,7 +369,8 @@ INSFVMomentumAdvection::coeffCalculator(const Elem & elem, const ADReal & mu) co // Moukalled 8.80, 8.82, and the orthogonal correction approach equation for E_f, // equation 8.89. So relative to the internal face viscous term, we have substituted // eqn. 8.82 for 8.78 - temp_coeff += mu * surface_vector.norm() / (fi->faceCentroid() - rc_centroid).norm(); + temp_coeff += + elem_mu * surface_vector.norm() / (fi->faceCentroid() - rc_centroid).norm(); // For flow boundaries, the coefficient addition is the same for every velocity component for (const auto i : make_range(_dim)) @@ -382,7 +383,7 @@ INSFVMomentumAdvection::coeffCalculator(const Elem & elem, const ADReal & mu) co { // Moukalled eqns. 15.154 - 15.156 for (const auto i : make_range(_dim)) - coeff(i) += 2. * mu * surface_vector.norm() / + coeff(i) += 2. * elem_mu * surface_vector.norm() / std::abs((fi->faceCentroid() - rc_centroid) * normal) * normal(i) * normal(i); @@ -399,6 +400,11 @@ INSFVMomentumAdvection::coeffCalculator(const Elem & elem, const ADReal & mu) co // Else we are on an internal face + const auto neighbor_mu = _mu(neighbor); + ADReal face_mu; + Moose::FV::interpolate( + Moose::FV::InterpMethod::Average, face_mu, elem_mu, neighbor_mu, *fi, elem_has_info); + ADRealVectorValue neighbor_velocity(_u_var->getNeighborValue(neighbor, *fi, elem_velocity(0))); if (_v_var) neighbor_velocity(1) = _v_var->getNeighborValue(neighbor, *fi, elem_velocity(1)); @@ -422,7 +428,8 @@ INSFVMomentumAdvection::coeffCalculator(const Elem & elem, const ADReal & mu) co // Now add the viscous flux. Note that this includes only the orthogonal component! See // Moukalled equations 8.80, 8.78, and the orthogonal correction approach equation for // E_f, equation 8.69 - temp_coeff += mu * surface_vector.norm() / (fi->neighborCentroid() - fi->elemCentroid()).norm(); + temp_coeff += + face_mu * surface_vector.norm() / (fi->neighborCentroid() - fi->elemCentroid()).norm(); // For internal faces the coefficient is the same for every velocity component. for (const auto i : make_range(_dim)) @@ -473,6 +480,9 @@ INSFVMomentumAdvection::interpolate(Moose::FV::InterpMethod m, if (m == Moose::FV::InterpMethod::Average) return; + mooseAssert(neighbor && this->hasBlocks(neighbor->subdomain_id()), + "We should be on an internal face..."); + // Get pressure gradient. This is the uncorrected gradient plus a correction from cell centroid // values on either side of the face const VectorValue & grad_p = _p_var->adGradSln(*_face_info); @@ -482,17 +492,15 @@ INSFVMomentumAdvection::interpolate(Moose::FV::InterpMethod m, const VectorValue & unc_grad_p = _p_var->uncorrectedAdGradSln(*_face_info); const Point & elem_centroid = _face_info->elemCentroid(); - const Point * const neighbor_centroid = neighbor ? &_face_info->neighborCentroid() : nullptr; + const Point & neighbor_centroid = _face_info->neighborCentroid(); Real elem_volume = _face_info->elemVolume(); - Real neighbor_volume = neighbor ? _face_info->neighborVolume() : 0; - const auto & elem_mu = _mu_elem[_qp]; + Real neighbor_volume = _face_info->neighborVolume(); // Now we need to perform the computations of D - const VectorValue & elem_a = rcCoeff(*elem, elem_mu); + const VectorValue & elem_a = rcCoeff(*elem); - mooseAssert(neighbor ? _subproblem.getCoordSystem(elem->subdomain_id()) == - _subproblem.getCoordSystem(neighbor->subdomain_id()) - : true, + mooseAssert(_subproblem.getCoordSystem(elem->subdomain_id()) == + _subproblem.getCoordSystem(neighbor->subdomain_id()), "Coordinate systems must be the same between the two elements"); Real coord; @@ -509,26 +517,19 @@ INSFVMomentumAdvection::interpolate(Moose::FV::InterpMethod m, VectorValue face_D; - if (neighbor && this->hasBlocks(neighbor->subdomain_id())) - { - const auto & neighbor_mu = _mu_neighbor[_qp]; - - const VectorValue & neighbor_a = rcCoeff(*neighbor, neighbor_mu); + const VectorValue & neighbor_a = rcCoeff(*neighbor); - coordTransformFactor(_subproblem, neighbor->subdomain_id(), *neighbor_centroid, coord); - neighbor_volume *= coord; + coordTransformFactor(_subproblem, neighbor->subdomain_id(), neighbor_centroid, coord); + neighbor_volume *= coord; - VectorValue neighbor_D = 0; - for (const auto i : make_range(_dim)) - { - mooseAssert(neighbor_a(i).value() != 0, "We should not be dividing by zero"); - neighbor_D(i) = neighbor_volume / neighbor_a(i); - } - Moose::FV::interpolate( - Moose::FV::InterpMethod::Average, face_D, elem_D, neighbor_D, *_face_info, true); + VectorValue neighbor_D = 0; + for (const auto i : make_range(_dim)) + { + mooseAssert(neighbor_a(i).value() != 0, "We should not be dividing by zero"); + neighbor_D(i) = neighbor_volume / neighbor_a(i); } - else - face_D = elem_D; + Moose::FV::interpolate( + Moose::FV::InterpMethod::Average, face_D, elem_D, neighbor_D, *_face_info, true); // perform the pressure correction for (const auto i : make_range(_dim)) @@ -537,7 +538,7 @@ INSFVMomentumAdvection::interpolate(Moose::FV::InterpMethod m, #else VectorValue -INSFVMomentumAdvection::coeffCalculator(const Elem &, const ADReal &) const +INSFVMomentumAdvection::coeffCalculator(const Elem &) const { mooseError("INSFVMomentumAdvection only works with global AD indexing"); } diff --git a/modules/navier_stokes/src/fvkernels/PINSFVMomentumAdvection.C b/modules/navier_stokes/src/fvkernels/PINSFVMomentumAdvection.C index 96c82d388326..7b1aafe1a1ad 100644 --- a/modules/navier_stokes/src/fvkernels/PINSFVMomentumAdvection.C +++ b/modules/navier_stokes/src/fvkernels/PINSFVMomentumAdvection.C @@ -40,7 +40,7 @@ PINSFVMomentumAdvection::PINSFVMomentumAdvection(const InputParameters & params) #ifdef MOOSE_GLOBAL_AD_INDEXING VectorValue -PINSFVMomentumAdvection::coeffCalculator(const Elem & elem, const ADReal & mu) const +PINSFVMomentumAdvection::coeffCalculator(const Elem & elem) const { // these coefficients arise from simple control volume balances of advection and diffusion. These // coefficients are the linear coefficients associated with the centroid of the control volume. @@ -76,7 +76,6 @@ PINSFVMomentumAdvection::coeffCalculator(const Elem & elem, const ADReal & mu) c auto action_functor = [&coeff, &elem_velocity, - &mu, #ifndef NDEBUG &elem, #endif @@ -105,6 +104,8 @@ PINSFVMomentumAdvection::coeffCalculator(const Elem & elem, const ADReal & mu) c // Rhie-Chow coefficient for. "neighbor" is the element across the current FaceInfo (fi) // face from the Rhie-Chow element + const auto elem_mu = _mu(&elem); + if (onBoundary(*fi)) { // Compute the face porosity @@ -118,7 +119,7 @@ PINSFVMomentumAdvection::coeffCalculator(const Elem & elem, const ADReal & mu) c { // Need to account for viscous shear stress from wall for (const auto i : make_range(_dim)) - coeff(i) += mu / eps_face * surface_vector.norm() / + coeff(i) += elem_mu / eps_face * surface_vector.norm() / std::abs((fi->faceCentroid() - rc_centroid) * normal) * (1 - normal(i) * normal(i)); @@ -154,8 +155,8 @@ PINSFVMomentumAdvection::coeffCalculator(const Elem & elem, const ADReal & mu) c // equation 8.89. So relative to the internal face viscous term, we have substituted // eqn. 8.82 for 8.78 // Note: If mu is an effective diffusivity, this should not be divided by eps_face - temp_coeff += - mu / eps_face * surface_vector.norm() / (fi->faceCentroid() - rc_centroid).norm(); + temp_coeff += elem_mu / eps_face * surface_vector.norm() / + (fi->faceCentroid() - rc_centroid).norm(); // For flow boundaries, the coefficient addition is the same for every velocity component for (const auto i : make_range(_dim)) @@ -168,7 +169,7 @@ PINSFVMomentumAdvection::coeffCalculator(const Elem & elem, const ADReal & mu) c { // Moukalled eqns. 15.154 - 15.156, adapted for porosity for (const auto i : make_range(_dim)) - coeff(i) += 2. * mu / eps_face * surface_vector.norm() / + coeff(i) += 2. * elem_mu / eps_face * surface_vector.norm() / std::abs((fi->faceCentroid() - rc_centroid) * normal) * normal(i) * normal(i); @@ -186,6 +187,11 @@ PINSFVMomentumAdvection::coeffCalculator(const Elem & elem, const ADReal & mu) c // Else we are on an internal face + const auto neighbor_mu = _mu(neighbor); + ADReal face_mu; + Moose::FV::interpolate( + Moose::FV::InterpMethod::Average, face_mu, elem_mu, neighbor_mu, *fi, elem_has_info); + // Compute the face porosity // Note: Try to be consistent with how the superficial velocity is computed in computeQpResidual const Real eps_face = MetaPhysicL::raw_value( @@ -216,7 +222,7 @@ PINSFVMomentumAdvection::coeffCalculator(const Elem & elem, const ADReal & mu) c // Now add the viscous flux. Note that this includes only the orthogonal component! See // Moukalled equations 8.80, 8.78, and the orthogonal correction approach equation for // E_f, equation 8.69 - temp_coeff += mu / eps_face * surface_vector.norm() / + temp_coeff += face_mu / eps_face * surface_vector.norm() / (fi->neighborCentroid() - fi->elemCentroid()).norm(); // For internal faces the coefficient is the same for every velocity component. @@ -274,6 +280,9 @@ PINSFVMomentumAdvection::interpolate(Moose::FV::InterpMethod m, MetaPhysicL::raw_value(_eps_var->adGradSln(neighbor)).norm() > 1e-12) return; + mooseAssert((neighbor && this->hasBlocks(neighbor->subdomain_id())), + "We should be on an internal face..."); + // Get pressure gradient. This is the uncorrected gradient plus a correction from cell centroid // values on either side of the face const VectorValue & grad_p = _p_var->adGradSln(*_face_info); @@ -283,17 +292,15 @@ PINSFVMomentumAdvection::interpolate(Moose::FV::InterpMethod m, const VectorValue & unc_grad_p = _p_var->uncorrectedAdGradSln(*_face_info); const Point & elem_centroid = _face_info->elemCentroid(); - const Point * const neighbor_centroid = neighbor ? &_face_info->neighborCentroid() : nullptr; + const Point & neighbor_centroid = _face_info->neighborCentroid(); Real elem_volume = _face_info->elemVolume(); - Real neighbor_volume = neighbor ? _face_info->neighborVolume() : 0; - const auto & elem_mu = _mu_elem[_qp]; + Real neighbor_volume = _face_info->neighborVolume(); // Now we need to perform the computations of D - const VectorValue & elem_a = rcCoeff(*elem, elem_mu); + const VectorValue & elem_a = rcCoeff(*elem); - mooseAssert(neighbor ? _subproblem.getCoordSystem(elem->subdomain_id()) == - _subproblem.getCoordSystem(neighbor->subdomain_id()) - : true, + mooseAssert(_subproblem.getCoordSystem(elem->subdomain_id()) == + _subproblem.getCoordSystem(neighbor->subdomain_id()), "Coordinate systems must be the same between the two elements"); Real coord; @@ -310,26 +317,19 @@ PINSFVMomentumAdvection::interpolate(Moose::FV::InterpMethod m, VectorValue face_D; - if (neighbor && this->hasBlocks(neighbor->subdomain_id())) - { - const auto & neighbor_mu = _mu_neighbor[_qp]; - - const VectorValue & neighbor_a = rcCoeff(*neighbor, neighbor_mu); + const VectorValue & neighbor_a = rcCoeff(*neighbor); - coordTransformFactor(_subproblem, neighbor->subdomain_id(), *neighbor_centroid, coord); - neighbor_volume *= coord; + coordTransformFactor(_subproblem, neighbor->subdomain_id(), neighbor_centroid, coord); + neighbor_volume *= coord; - VectorValue neighbor_D = 0; - for (const auto i : make_range(_dim)) - { - mooseAssert(neighbor_a(i).value() != 0, "We should not be dividing by zero"); - neighbor_D(i) = neighbor_volume / neighbor_a(i); - } - Moose::FV::interpolate( - Moose::FV::InterpMethod::Average, face_D, elem_D, neighbor_D, *_face_info, true); + VectorValue neighbor_D = 0; + for (const auto i : make_range(_dim)) + { + mooseAssert(neighbor_a(i).value() != 0, "We should not be dividing by zero"); + neighbor_D(i) = neighbor_volume / neighbor_a(i); } - else - face_D = elem_D; + Moose::FV::interpolate( + Moose::FV::InterpMethod::Average, face_D, elem_D, neighbor_D, *_face_info, true); // evaluate face porosity, see (18) in Hanimann 2021 or (11) in Nordlund 2016 const auto eps_face = MetaPhysicL::raw_value(_eps_var->getInternalFaceValue( @@ -342,7 +342,7 @@ PINSFVMomentumAdvection::interpolate(Moose::FV::InterpMethod m, #else VectorValue -PINSFVMomentumAdvection::coeffCalculator(const Elem &, const ADReal &) const +PINSFVMomentumAdvection::coeffCalculator(const Elem &) const { mooseError("PINSFVMomentumAdvection only works with global AD indexing"); } diff --git a/modules/navier_stokes/test/tests/finite_volume/ins/lid-driven/lid-driven.i b/modules/navier_stokes/test/tests/finite_volume/ins/lid-driven/lid-driven.i index 3912edc4009f..28f017b65722 100644 --- a/modules/navier_stokes/test/tests/finite_volume/ins/lid-driven/lid-driven.i +++ b/modules/navier_stokes/test/tests/finite_volume/ins/lid-driven/lid-driven.i @@ -61,7 +61,7 @@ rho=1 pressure = pressure u = u v = v - mu = ${mu} + mu = 'mu' rho = ${rho} [] [mean_zero_pressure] @@ -77,7 +77,7 @@ rho=1 pressure = pressure u = u v = v - mu = ${mu} + mu = 'mu' rho = ${rho} [] @@ -101,7 +101,7 @@ rho=1 pressure = pressure u = u v = v - mu = ${mu} + mu = 'mu' rho = ${rho} [] @@ -150,6 +150,11 @@ rho=1 pressure = 'pressure' rho = ${rho} [] + [mu] + type = ADGenericConstantFunctorMaterial + prop_names = 'mu' + prop_values = '${mu}' + [] [] [Preconditioning]