Skip to content

Commit

Permalink
Move to functor material property for mu in (P)INSFV
Browse files Browse the repository at this point in the history
  • Loading branch information
@lindsayad
lindsayad committed Jul 20, 2021
1 parent aa0a407 commit 22e5abd
Show file tree
Hide file tree
Showing 5 changed files with 81 additions and 78 deletions.
11 changes: 4 additions & 7 deletions modules/navier_stokes/include/fvkernels/INSFVMomentumAdvection.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<Real> & _mu_elem;

/// The dynamic viscosity on the FaceInfo neighbor
const ADMaterialProperty<Real> & _mu_neighbor;
/// The dynamic viscosity
const FunctorMaterialProperty<ADReal> & _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<ADReal> & rcCoeff(const Elem & elem, const ADReal & mu) const;
const VectorValue<ADReal> & 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<ADReal> coeffCalculator(const Elem & elem, const ADReal & mu) const;
virtual VectorValue<ADReal> coeffCalculator(const Elem & elem) const;

/**
* Clear the RC 'a' coefficient cache
Expand Down
2 changes: 1 addition & 1 deletion modules/navier_stokes/include/fvkernels/PINSFVMomentumAdvection.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ class PINSFVMomentumAdvection : public INSFVMomentumAdvection
const ADRealVectorValue & neighbor_v) override;

virtual ADReal computeQpResidual() override;
VectorValue<ADReal> coeffCalculator(const Elem & elem, const ADReal & mu) const override;
VectorValue<ADReal> coeffCalculator(const Elem & elem) const override;

/// porosity variable to compute gradients
const MooseVariableFV<Real> * const _eps_var;
Expand Down
71 changes: 36 additions & 35 deletions modules/navier_stokes/src/fvkernels/INSFVMomentumAdvection.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,8 +67,7 @@ INSFVMomentumAdvection::validParams()

INSFVMomentumAdvection::INSFVMomentumAdvection(const InputParameters & params)
: FVMatAdvection(params),
_mu_elem(getADMaterialProperty<Real>("mu")),
_mu_neighbor(getNeighborADMaterialProperty<Real>("mu")),
_mu(getFunctorMaterialProperty<ADReal>("mu")),
_p_var(dynamic_cast<const INSFVPressureVariable *>(getFieldVar(NS::pressure, 0))),
_u_var(dynamic_cast<const INSFVVelocityVariable *>(getFieldVar("u", 0))),
_v_var(params.isParamValid("v")
Expand Down Expand Up @@ -233,7 +232,7 @@ INSFVMomentumAdvection::skipForBoundary(const FaceInfo & fi) const
}

const VectorValue<ADReal> &
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(),
Expand All @@ -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");

Expand All @@ -260,7 +259,7 @@ INSFVMomentumAdvection::rcCoeff(const Elem & elem, const ADReal & mu) const

#ifdef MOOSE_GLOBAL_AD_INDEXING
VectorValue<ADReal>
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.
Expand Down Expand Up @@ -296,7 +295,6 @@ INSFVMomentumAdvection::coeffCalculator(const Elem & elem, const ADReal & mu) co

auto action_functor = [&coeff,
&elem_velocity,
&mu,
#ifndef NDEBUG
&elem,
#endif
Expand All @@ -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
Expand All @@ -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));

Expand Down Expand Up @@ -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))
Expand All @@ -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);

Expand All @@ -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));
Expand All @@ -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))
Expand Down Expand Up @@ -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<ADReal> & grad_p = _p_var->adGradSln(*_face_info);
Expand All @@ -482,17 +492,15 @@ INSFVMomentumAdvection::interpolate(Moose::FV::InterpMethod m,
const VectorValue<ADReal> & 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<ADReal> & elem_a = rcCoeff(*elem, elem_mu);
const VectorValue<ADReal> & 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;
Expand All @@ -509,26 +517,19 @@ INSFVMomentumAdvection::interpolate(Moose::FV::InterpMethod m,

VectorValue<ADReal> face_D;

if (neighbor && this->hasBlocks(neighbor->subdomain_id()))
{
const auto & neighbor_mu = _mu_neighbor[_qp];

const VectorValue<ADReal> & neighbor_a = rcCoeff(*neighbor, neighbor_mu);
const VectorValue<ADReal> & 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<ADReal> 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<ADReal> 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))
Expand All @@ -537,7 +538,7 @@ INSFVMomentumAdvection::interpolate(Moose::FV::InterpMethod m,
#else

VectorValue<ADReal>
INSFVMomentumAdvection::coeffCalculator(const Elem &, const ADReal &) const
INSFVMomentumAdvection::coeffCalculator(const Elem &) const
{
mooseError("INSFVMomentumAdvection only works with global AD indexing");
}
Expand Down
64 changes: 32 additions & 32 deletions modules/navier_stokes/src/fvkernels/PINSFVMomentumAdvection.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ PINSFVMomentumAdvection::PINSFVMomentumAdvection(const InputParameters & params)

#ifdef MOOSE_GLOBAL_AD_INDEXING
VectorValue<ADReal>
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.
Expand Down Expand Up @@ -76,7 +76,6 @@ PINSFVMomentumAdvection::coeffCalculator(const Elem & elem, const ADReal & mu) c

auto action_functor = [&coeff,
&elem_velocity,
&mu,
#ifndef NDEBUG
&elem,
#endif
Expand Down Expand Up @@ -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
Expand All @@ -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));

Expand Down Expand Up @@ -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))
Expand All @@ -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);

Expand All @@ -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(
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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<ADReal> & grad_p = _p_var->adGradSln(*_face_info);
Expand All @@ -283,17 +292,15 @@ PINSFVMomentumAdvection::interpolate(Moose::FV::InterpMethod m,
const VectorValue<ADReal> & 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<ADReal> & elem_a = rcCoeff(*elem, elem_mu);
const VectorValue<ADReal> & 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;
Expand All @@ -310,26 +317,19 @@ PINSFVMomentumAdvection::interpolate(Moose::FV::InterpMethod m,

VectorValue<ADReal> face_D;

if (neighbor && this->hasBlocks(neighbor->subdomain_id()))
{
const auto & neighbor_mu = _mu_neighbor[_qp];

const VectorValue<ADReal> & neighbor_a = rcCoeff(*neighbor, neighbor_mu);
const VectorValue<ADReal> & 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<ADReal> 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<ADReal> 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(
Expand All @@ -342,7 +342,7 @@ PINSFVMomentumAdvection::interpolate(Moose::FV::InterpMethod m,
#else

VectorValue<ADReal>
PINSFVMomentumAdvection::coeffCalculator(const Elem &, const ADReal &) const
PINSFVMomentumAdvection::coeffCalculator(const Elem &) const
{
mooseError("PINSFVMomentumAdvection only works with global AD indexing");
}
Expand Down
Loading

0 comments on commit 22e5abd

Please sign in to comment.