INSFVMixingLengthReynoldsStress kernel modification

[aiskhak]

I am modifying the INSFVMixingLengthReynoldsStress kernel in a way that I read Reynolds stresses from a .csv file (supplied as auxiliary variables). Below is the code. While the main part is sort of clear, RC interpolation is vague to me. In the original kernel the coefficients were computing using autograd feature. What should I do to make it correctly in my case?

#include "INSFV_rs_csv_ReynoldsStress.h"
#include "INSFVVelocityVariable.h"
#include "NS.h"
#include "SystemBase.h"

registerMooseObject("NavierStokesApp", INSFV_rs_csv_ReynoldsStress);

InputParameters
INSFV_rs_csv_ReynoldsStress::validParams()
{
  InputParameters params = INSFVFluxKernel::validParams();
  params.addClassDescription(
      "Computes the force due to the Reynolds stress term in the incompressible"
      " Reynolds-averaged Navier-Stokes equations.");
  params.addRequiredCoupledVar("u", "The velocity in the x direction.");
  params.addCoupledVar("v", "The velocity in the y direction.");
  params.addCoupledVar("w", "The velocity in the z direction.");
  params.addRequiredParam<MooseFunctorName>(NS::density, "fluid density");
  params.addRequiredParam<MooseFunctorName>("rs_xx", "Reynolds stress xx from CSV file.");
  params.addParam<MooseFunctorName>("rs_xy", "Reynolds stress xy from CSV file.");
  params.addParam<MooseFunctorName>("rs_yy", "Reynolds stress yy from CSV file.");
  params.addParam<MooseFunctorName>("rs_xz", "Reynolds stress xz from CSV file.");
  params.addParam<MooseFunctorName>("rs_yz", "Reynolds stress yz from CSV file.");
  params.addParam<MooseFunctorName>("rs_zz", "Reynolds stress zz from CSV file.");
  params.addParam<Real>("rs_m_xx", 1.0, "RS multiplier");
  params.addParam<Real>("rs_m_xy", 1.0, "RS multiplier");
  params.addParam<Real>("rs_m_yy", 1.0, "RS multiplier");
  params.addParam<Real>("rs_m_xz", 1.0, "RS multiplier");
  params.addParam<Real>("rs_m_yz", 1.0, "RS multiplier");
  params.addParam<Real>("rs_m_zz", 1.0, "RS multiplier");
  MooseEnum momentum_component("x=0 y=1 z=2");
  params.addRequiredParam<MooseEnum>("momentum_component",
                                      momentum_component, "The component of the momentum equation that this kernel applies to.");
  // We assume the worst, e.g. we are doing Rhie-Chow. In that case we need three layers. An 'a'
  // coefficient evaluation at a face will necessitate evaluation of *this* object at every face of
  // the adjoining element, necessitating a face gradient evaluation at those faces, necessitating a
  // cell gradient evaluation in neighboring elements, necessitating cell value evaluations in
  // neighbors of those neighbor elements
  params.set<unsigned short>("ghost_layers") = 3;
  return params;
}

INSFV_rs_csv_ReynoldsStress::INSFV_rs_csv_ReynoldsStress(const InputParameters & params)
  : INSFVFluxKernel(params),
    _dim(_subproblem.mesh().dimension()),
    _axis_index(getParam<MooseEnum>("momentum_component")),
    _u_var(dynamic_cast<const INSFVVelocityVariable *>(getFieldVar("u", 0))),
    _v_var(params.isParamValid("v") ? dynamic_cast<const INSFVVelocityVariable *>(getFieldVar("v", 0))
                                    : nullptr),
    _w_var(params.isParamValid("w") ? dynamic_cast<const INSFVVelocityVariable *>(getFieldVar("w", 0))
                                    : nullptr),
    _rho(getFunctor<ADReal>(NS::density)),
    _rs_xx(getFunctor<ADReal>("rs_xx")),
    _rs_xy(params.isParamValid("v") ? getFunctor<ADReal>("rs_xy") 
                                    : getFunctor<ADReal>("rs_xx")), // assign rs_xx if 1D (not needed)
    _rs_yy(params.isParamValid("v") ? getFunctor<ADReal>("rs_yy") 
                                    : getFunctor<ADReal>("rs_xx")), // assign rs_xx if 1D (not needed)
    _rs_xz(params.isParamValid("w") ? getFunctor<ADReal>("rs_xz") 
                                    : getFunctor<ADReal>("rs_xx")), // assign rs_xx if 2D (not needed)
    _rs_yz(params.isParamValid("w") ? getFunctor<ADReal>("rs_yz") 
                                    : getFunctor<ADReal>("rs_xx")), // assign rs_xx if 2D (not needed)
    _rs_zz(params.isParamValid("w") ? getFunctor<ADReal>("rs_zz") 
                                    : getFunctor<ADReal>("rs_xx")), // assign rs_xx if 2D (not needed)
    _rs_m_xx(getParam<Real>("rs_m_xx")), // can nulify/scale the RS from input file
	_rs_m_xy(getParam<Real>("rs_m_xy")),
	_rs_m_yy(getParam<Real>("rs_m_yy")),
	_rs_m_xz(getParam<Real>("rs_m_xz")),
	_rs_m_yz(getParam<Real>("rs_m_yz")),
	_rs_m_zz(getParam<Real>("rs_m_zz"))
{
#ifndef MOOSE_GLOBAL_AD_INDEXING
  mooseError("INSFV is not supported by local AD indexing. In order to use INSFV, please run the "
             "configure script in the root MOOSE directory with the configure option "
             "'--with-ad-indexing-type=global'");
#endif

  if (!_u_var)
    paramError("u", "the u velocity must be an INSFVVelocityVariable.");

  if (_dim >= 2 && !_v_var)
    paramError("v",
               "In two or more dimensions, the v velocity must be supplied and it must be an "
               "INSFVVelocityVariable.");

  if (_dim >= 3 && !_w_var)
    paramError("w",
               "In three-dimensions, the w velocity must be supplied and it must be an "
               "INSFVVelocityVariable.");
}

ADReal
INSFV_rs_csv_ReynoldsStress::computeStrongResidual()
{
#ifdef MOOSE_GLOBAL_AD_INDEXING
  
  // interpolate density to the face
  const auto face = Moose::FV::makeCDFace(*_face_info, faceArgSubdomains());
  const ADReal rho = _rho(face);
  
  // compute Reynolds stresses contribution fo each direction
  ADReal rs = 0.0;	// initialize
  
  // 1D case
  if (_dim == 1)
  {
    if (_index == 0)
    {
      ADReal rs_xx = _rs_xx(face);
      rs += rs_xx*_normal(0)*_rs_m_xx; // x direction
    }
      else
    {
      mooseError("if case is 1D use x direction");
    }
  }
  
  // 2D case
  else if (_dim == 2) 
  {
    if (_index == 0) 
    {
      ADReal rs_xx = _rs_xx(face)*_rs_m_xx;
      ADReal rs_xy = _rs_xy(face)*_rs_m_xy;
      rs += rs_xx*_normal(0) + rs_xy*_normal(1); // x direction
    }
    else if (_index == 1)
    {
      ADReal rs_xy = _rs_xy(face)*_rs_m_xy;
      ADReal rs_yy = _rs_yy(face)*_rs_m_yy;
      rs += rs_xy*_normal(0) + rs_yy*_normal(1); // y direction
    }
    else 
    {
      mooseError("if case is 2D use x and y directions"); 
    }
  }
 
  // 3D case
  else if (_dim == 3) 
  {
    if (_index == 0) 
    {
      ADReal rs_xx = _rs_xx(face)*_rs_m_xx;
      ADReal rs_xy = _rs_xy(face)*_rs_m_xy;
      ADReal rs_xz = _rs_xz(face)*_rs_m_xz;
      rs += rs_xx*_normal(0) + rs_xy*_normal(1) + rs_xz*_normal(2); // x direction
    }
    else if (_index == 1)
    {
      ADReal rs_xy = _rs_xy(face)*_rs_m_xy;
      ADReal rs_yy = _rs_yy(face)*_rs_m_yy;
      ADReal rs_yz = _rs_yz(face)*_rs_m_yz;
      rs += rs_xy*_normal(0) + rs_yy*_normal(1) + rs_yz*_normal(2); // y direction
    }
    else if (_index == 2)
    {
      ADReal rs_xz = _rs_xz(face)*_rs_m_xz;
      ADReal rs_yz = _rs_yz(face)*_rs_m_yz;
      ADReal rs_zz = _rs_zz(face)*_rs_m_zz;
      rs += rs_xz*_normal(0) + rs_yz*_normal(1) + rs_zz*_normal(2); // z direction
    }
    else
    {
      mooseError("wrong index"); 
    }
  }
  
  // tally contribution of face flux to Rhie-Chow coefficient 
  // of element (ae) and neighbor (an) of the face
  if (_face_type == FaceInfo::VarFaceNeighbors::ELEM || _face_type == FaceInfo::VarFaceNeighbors::BOTH)
  {
    const auto dof_number = _face_info->elem().dof_number(_sys.number(), _var.number(), 0);
    // norm_strain_rate is a linear combination of degrees of freedom so it's safe to straight-up
    // index into the derivatives vector at the dof we care about
    //_ae = norm_strain_rate.derivatives()[dof_number];
    //_ae *= -rho * eddy_diff;
    _ae = 0.0; //rho*rs.derivatives()[dof_number];
  }
  if (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR || _face_type == FaceInfo::VarFaceNeighbors::BOTH)
  {
    const auto dof_number = _face_info->neighbor().dof_number(_sys.number(), _var.number(), 0);
    //_an = norm_strain_rate.derivatives()[dof_number];
    //_an *= rho * eddy_diff;
    _an = 0.0; //rho*rs.derivatives()[dof_number];
  }

  // Return the turbulent stress contribution to the momentum equation
  //return -1 * rho * eddy_diff * norm_strain_rate;
  return rho * rs;

#else
  return 0;
#endif
}

// Move tallied contribution of Rhie-Chow coefficient into the actual coefficient, 
// held by the RC Interpolator user object. 
// The contribution is weighted by the area of the face (coord() is the RZ radius)
void
INSFV_rs_csv_ReynoldsStress::gatherRCData(const FaceInfo & fi)
{
  if (skipForBoundary(fi))
    return;

  _face_info = &fi;
  _normal = fi.normal();
  _face_type = fi.faceType(_var.name());

  processResidualAndJacobian(computeStrongResidual() * (fi.faceArea() * fi.faceCoord()));

  if (_face_type == FaceInfo::VarFaceNeighbors::ELEM || _face_type == FaceInfo::VarFaceNeighbors::BOTH)
    _rc_uo.addToA(&fi.elem(), _index, _ae * (fi.faceArea() * fi.faceCoord()));
  if (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR || _face_type == FaceInfo::VarFaceNeighbors::BOTH)
    _rc_uo.addToA(fi.neighborPtr(), _index, _an * (fi.faceArea() * fi.faceCoord()));
}

[GiudGiud]

@lindsayad I m curious about your thoughts on a flux kernel with the flux that does not depend on velocity (reynolds stresses preset from a CSV file), whether they should contribute to the RC coefficients.

[lindsayad]

If the stresses do not depend on the velocity, then they should not contribute to the RC 'a' coefficients

[GiudGiud]

I think the real Reynolds stresses should, but here they dont because they are plugged in functions from CSV

[lindsayad]

Well this is the user's choice as to how they want to model it. If in their model the Reynolds stress does not change as a function of velocity, than in the linearized form of the discretized momentum equation, the Reynolds stresses would not be in the 'a' coefficients as they multiply velocity.