Skip to content

Commit

Permalink
Added strong viscous term support to INSAD
Browse files Browse the repository at this point in the history 
Now comparing equal with INS with all stabilization terms for
PSPG
  • Loading branch information
@lindsayad
lindsayad committed Aug 23, 2023
1 parent f18bcb4 commit 33a13dc
Show file tree
Hide file tree
Showing 3 changed files with 204 additions and 91 deletions.
9 changes: 0 additions & 9 deletions modules/navier_stokes/include/materials/INSADMaterial.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,11 +25,6 @@ class INSADMaterial : public Material
protected:
virtual void computeQpProperties() override;

/**
* compute the strong form corresponding to RZ pieces of the viscous term
*/
void viscousTermRZ();

/// velocity
const ADVectorVariableValue & _velocity;

Expand All @@ -54,10 +49,6 @@ class INSADMaterial : public Material
/// Strong residual corresponding to the momentum advective term
ADMaterialProperty<RealVectorValue> & _advective_strong_residual;

/// Strong residual corresponding to the momentum viscous term. This is only used by stabilization
/// kernels
ADMaterialProperty<RealVectorValue> & _viscous_strong_residual;

/// Strong residual corresponding to the momentum transient term
ADMaterialProperty<RealVectorValue> & _td_strong_residual;

Expand Down
212 changes: 204 additions & 8 deletions modules/navier_stokes/include/materials/INSADTauMaterial.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,9 +16,13 @@
#include "MooseArray.h"
#include "INSADMaterial.h"
#include "NavierStokesMethods.h"
#include "Assembly.h"
#include "MooseVariableFE.h"
#include "MooseMesh.h"

#include "libmesh/elem.h"
#include "libmesh/node.h"
#include "libmesh/fe_type.h"

#include <vector>

Expand All @@ -35,17 +39,60 @@ class INSADTauMaterialTempl : public T
protected:
virtual void computeProperties() override;
virtual void computeQpProperties() override;

/**
* Compute the maximum dimension of the current element
*/
void computeHMax();

/**
* Compute the viscuous strong residual
*/
void computeViscousStrongResidual();

/**
* Compute the strong form corresponding to RZ pieces of the viscous term
*/
void viscousTermRZ();

/**
* Whether to seed with respect to velocity derivatives
*/
bool doVelocityDerivatives() const;

const Real _alpha;
ADMaterialProperty<Real> & _tau;

/// Strong residual corresponding to the momentum viscous term. This is only used by stabilization
/// kernels
ADMaterialProperty<RealVectorValue> & _viscous_strong_residual;

/// The strong residual of the momentum equation
ADMaterialProperty<RealVectorValue> & _momentum_strong_residual;

ADReal _hmax;

/// The velocity variable
const VectorMooseVariable * const _velocity_var;

/// A scalar Lagrange FE data member for using to compute the velocity second derivatives since
/// they're currently not supported for vector FE types
const FEBase * const & _scalar_lagrange_fe;

/// Containers to hold second derivatives with respect to velocities
std::vector<ADRealTensorValue> _d2u;
std::vector<ADRealTensorValue> _d2v;
std::vector<ADRealTensorValue> _d2w;

/// The velocity variable number
const unsigned int _vel_number;

/// The velocity system number
const unsigned int _vel_sys_number;

using T::_ad_q_point;
using T::_advective_strong_residual;
using T::_assembly;
using T::_boussinesq_strong_residual;
using T::_convected_mesh_strong_residual;
using T::_coord_sys;
Expand All @@ -60,19 +107,27 @@ class INSADTauMaterialTempl : public T
using T::_dt;
using T::_fe_problem;
using T::_grad_p;
using T::_grad_velocity;
using T::_gravity_strong_residual;
using T::_has_boussinesq;
using T::_has_convected_mesh;
using T::_has_coupled_force;
using T::_has_gravity;
using T::_has_transient;
using T::_mesh;
using T::_mu;
using T::_object_tracker;
using T::_q_point;
using T::_qp;
using T::_qrule;
using T::_rho;
using T::_rz_axial_coord;
using T::_rz_radial_coord;
using T::_td_strong_residual;
using T::_use_displaced_mesh;
using T::_velocity;
using T::_viscous_strong_residual;
using T::_viscous_form;
using T::getVectorVar;
};

typedef INSADTauMaterialTempl<INSADMaterial> INSADTauMaterial;
Expand All @@ -93,9 +148,25 @@ INSADTauMaterialTempl<T>::INSADTauMaterialTempl(const InputParameters & paramete
: T(parameters),
_alpha(this->template getParam<Real>("alpha")),
_tau(this->template declareADProperty<Real>("tau")),
_viscous_strong_residual(
this->template declareADProperty<RealVectorValue>("viscous_strong_residual")),
_momentum_strong_residual(
this->template declareADProperty<RealVectorValue>("momentum_strong_residual"))
this->template declareADProperty<RealVectorValue>("momentum_strong_residual")),
_velocity_var(getVectorVar("velocity", 0)),
_scalar_lagrange_fe(
_assembly.getFE(FEType(_velocity_var->feType().order, LAGRANGE), _mesh.dimension())),
_vel_number(_velocity_var->number()),
_vel_sys_number(_velocity_var->sys().number())
{
_scalar_lagrange_fe->get_d2phi();
}

template <typename T>
bool
INSADTauMaterialTempl<T>::doVelocityDerivatives() const
{
return ADReal::do_derivatives &&
(_vel_sys_number == _fe_problem.currentNonlinearSystem().number());
}

template <typename T>
Expand Down Expand Up @@ -140,10 +211,139 @@ void
INSADTauMaterialTempl<T>::computeProperties()
{
computeHMax();
computeViscousStrongResidual();

T::computeProperties();
}

template <typename T>
void
INSADTauMaterialTempl<T>::computeViscousStrongResidual()
{
auto resize_and_zero = [this](auto & d2vel)
{
d2vel.resize(_qrule->n_points());
for (auto & d2qp : d2vel)
d2qp = 0;
};
resize_and_zero(_d2u);
resize_and_zero(_d2v);
resize_and_zero(_d2w);

auto get_d2 = [this](const auto i) -> std::vector<ADRealTensorValue> &
{
switch (i)
{
case 0:
return _d2u;
case 1:
return _d2v;
case 2:
return _d2w;
default:
mooseError("invalid value of 'i'");
}
};

const auto & vel_dof_indices = _velocity_var->dofIndices();
for (const auto i : index_range(vel_dof_indices))
{
// This may not work if the element and spatial dimensions are different
const auto dimensional_component = i % _mesh.dimension();
auto & d2vel = get_d2(dimensional_component);
const auto dof_index = vel_dof_indices[i];
ADReal dof_value = (*_velocity_var->sys().currentSolution())(dof_index);
if (doVelocityDerivatives())
dof_value.derivatives().insert(dof_index) = 1;
const auto scalar_i_component = i / _mesh.dimension();
for (const auto qp : make_range(_qrule->n_points()))
d2vel[qp] += dof_value * _scalar_lagrange_fe->get_d2phi()[scalar_i_component][qp];
}

for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
{
_viscous_strong_residual[_qp](0) = -_mu[_qp] * _d2u[_qp].tr();
_viscous_strong_residual[_qp](1) =
_mesh.dimension() >= 2 ? -_mu[_qp] * _d2v[_qp].tr() : ADReal(0);
_viscous_strong_residual[_qp](2) =
_mesh.dimension() == 3 ? -_mu[_qp] * _d2w[_qp].tr() : ADReal(0);
if (_viscous_form == "traction")
{
_viscous_strong_residual[_qp] -= _mu[_qp] * _d2u[_qp].row(0);
if (_mesh.dimension() >= 2)
_viscous_strong_residual[_qp] -= _mu[_qp] * _d2v[_qp].row(1);
if (_mesh.dimension() == 3)
_viscous_strong_residual[_qp] -= _mu[_qp] * _d2w[_qp].row(2);
}
if (_coord_sys == Moose::COORD_RZ)
viscousTermRZ();
}
}

template <typename T>
void
INSADTauMaterialTempl<T>::viscousTermRZ()
{
// To understand the code immediately below, visit
// https://en.wikipedia.org/wiki/Del_in_cylindrical_and_spherical_coordinates.
// The u_r / r^2 term comes from the vector Laplacian. The -du_i/dr * 1/r term comes from
// the scalar Laplacian. The scalar Laplacian in axisymmetric cylindrical coordinates is
// equivalent to the Cartesian Laplacian plus a 1/r * du_i/dr term. And of course we are
// applying a minus sign here because the strong form is -\nabala^2 * \vec{u}
//
// Another note: libMesh implements grad(v) as dvi/dxj

if (_use_displaced_mesh)
{
ADReal r = _ad_q_point[_qp](_rz_radial_coord);

if (_viscous_form == "laplace")
_viscous_strong_residual[_qp] += ADRealVectorValue(
// u_r
// Additional term from vector Laplacian
_mu[_qp] * (_velocity[_qp](_rz_radial_coord) / (r * r) -
// Additional term from scalar Laplacian
_grad_velocity[_qp](_rz_radial_coord, _rz_radial_coord) / r),
// u_z
// Additional term from scalar Laplacian
-_mu[_qp] * _grad_velocity[_qp](_rz_axial_coord, _rz_radial_coord) / r,
0);
else
_viscous_strong_residual[_qp] +=
ADRealVectorValue(2. * _mu[_qp] *
(_velocity[_qp](_rz_radial_coord) / (r * r) -
_grad_velocity[_qp](_rz_radial_coord, _rz_radial_coord) / r),
-_mu[_qp] / r * (_grad_velocity[_qp](1, 0) + _grad_velocity[_qp](0, 1)),
0);
}
else
{
Real r = _q_point[_qp](_rz_radial_coord);
if (_viscous_form == "laplace")
_viscous_strong_residual[_qp] +=
// u_r
// Additional term from vector Laplacian
ADRealVectorValue(
_mu[_qp] * (_velocity[_qp](_rz_radial_coord) /
(_q_point[_qp](_rz_radial_coord) * _q_point[_qp](_rz_radial_coord)) -
// Additional term from scalar Laplacian
_grad_velocity[_qp](_rz_radial_coord, _rz_radial_coord) /
_q_point[_qp](_rz_radial_coord)),
// u_z
// Additional term from scalar Laplacian
-_mu[_qp] * _grad_velocity[_qp](_rz_axial_coord, _rz_radial_coord) /
_q_point[_qp](_rz_radial_coord),
0);
else
_viscous_strong_residual[_qp] +=
ADRealVectorValue(2. * _mu[_qp] *
(_velocity[_qp](_rz_radial_coord) / (r * r) -
_grad_velocity[_qp](_rz_radial_coord, _rz_radial_coord) / r),
-_mu[_qp] / r * (_grad_velocity[_qp](1, 0) + _grad_velocity[_qp](0, 1)),
0);
}
}

template <typename T>
void
INSADTauMaterialTempl<T>::computeQpProperties()
Expand All @@ -156,12 +356,8 @@ INSADTauMaterialTempl<T>::computeQpProperties()
_tau[_qp] = _alpha / std::sqrt(transient_part + (2. * speed / _hmax) * (2. * speed / _hmax) +
9. * (4. * nu / (_hmax * _hmax)) * (4. * nu / (_hmax * _hmax)));

_momentum_strong_residual[_qp] = _advective_strong_residual[_qp] + _grad_p[_qp];

// Since we can't current compute vector Laplacians we only have strong residual contributions
// from the viscous term in the RZ coordinate system
if (_coord_sys == Moose::COORD_RZ)
_momentum_strong_residual[_qp] += _viscous_strong_residual[_qp];
_momentum_strong_residual[_qp] =
_advective_strong_residual[_qp] + _viscous_strong_residual[_qp] + _grad_p[_qp];

if (_has_transient)
_momentum_strong_residual[_qp] += _td_strong_residual[_qp];
Expand Down
74 changes: 0 additions & 74 deletions modules/navier_stokes/src/materials/INSADMaterial.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,6 @@ INSADMaterial::INSADMaterial(const InputParameters & parameters)
_velocity_dot(nullptr),
_mass_strong_residual(declareADProperty<Real>("mass_strong_residual")),
_advective_strong_residual(declareADProperty<RealVectorValue>("advective_strong_residual")),
_viscous_strong_residual(declareADProperty<RealVectorValue>("viscous_strong_residual")),
// We have to declare the below strong residuals for integrity check purposes even though we may
// not compute them. This may incur some unnecessary cost for a non-sparse derivative container
// since when the properties are resized the entire non-sparse derivative containers will be
Expand Down Expand Up @@ -239,77 +238,4 @@ INSADMaterial::computeQpProperties()
// _mms_function_strong_residual[_qp] = -RealVectorValue(_x_vel_fn.value(_t, _q_point[_qp]),
// _y_vel_fn.value(_t, _q_point[_qp]),
// _z_vel_fn.value(_t, _q_point[_qp]));

// // The code immediately below is fictional. E.g. there is no Moose::Laplacian nor is there
// // currently a _second_velocity member because TypeNTensor (where N = 3 in this case) math is not
// // really implemented in libMesh. Hence we cannot add this strong form contribution of the viscous
// // term at this time. Note that for linear elements this introduces no error in the consistency of
// // stabilization methods, and in general for bi-linear elements, the error introduced is small
// _viscous_strong_residual[_qp] = -_mu[_qp] * Moose::Laplacian(_second_velocity[_qp]);

if (_coord_sys == Moose::COORD_RZ)
viscousTermRZ();
}

void
INSADMaterial::viscousTermRZ()
{
// To understand the code immediately below, visit
// https://en.wikipedia.org/wiki/Del_in_cylindrical_and_spherical_coordinates.
// The u_r / r^2 term comes from the vector Laplacian. The -du_i/dr * 1/r term comes from
// the scalar Laplacian. The scalar Laplacian in axisymmetric cylindrical coordinates is
// equivalent to the Cartesian Laplacian plus a 1/r * du_i/dr term. And of course we are
// applying a minus sign here because the strong form is -\nabala^2 * \vec{u}
//
// Another note: libMesh implements grad(v) as dvi/dxj

if (_use_displaced_mesh)
{
ADReal r = _ad_q_point[_qp](_rz_radial_coord);

if (_viscous_form == "laplace")
_viscous_strong_residual[_qp] = ADRealVectorValue(
// u_r
// Additional term from vector Laplacian
_mu[_qp] * (_velocity[_qp](_rz_radial_coord) / (r * r) -
// Additional term from scalar Laplacian
_grad_velocity[_qp](_rz_radial_coord, _rz_radial_coord) / r),
// u_z
// Additional term from scalar Laplacian
-_mu[_qp] * _grad_velocity[_qp](_rz_axial_coord, _rz_radial_coord) / r,
0);
else
_viscous_strong_residual[_qp] =
ADRealVectorValue(2. * _mu[_qp] *
(_velocity[_qp](_rz_radial_coord) / (r * r) -
_grad_velocity[_qp](_rz_radial_coord, _rz_radial_coord) / r),
-_mu[_qp] / r * (_grad_velocity[_qp](1, 0) + _grad_velocity[_qp](0, 1)),
0);
}
else
{
Real r = _q_point[_qp](_rz_radial_coord);
if (_viscous_form == "laplace")
_viscous_strong_residual[_qp] =
// u_r
// Additional term from vector Laplacian
ADRealVectorValue(
_mu[_qp] * (_velocity[_qp](_rz_radial_coord) /
(_q_point[_qp](_rz_radial_coord) * _q_point[_qp](_rz_radial_coord)) -
// Additional term from scalar Laplacian
_grad_velocity[_qp](_rz_radial_coord, _rz_radial_coord) /
_q_point[_qp](_rz_radial_coord)),
// u_z
// Additional term from scalar Laplacian
-_mu[_qp] * _grad_velocity[_qp](_rz_axial_coord, _rz_radial_coord) /
_q_point[_qp](_rz_radial_coord),
0);
else
_viscous_strong_residual[_qp] =
ADRealVectorValue(2. * _mu[_qp] *
(_velocity[_qp](_rz_radial_coord) / (r * r) -
_grad_velocity[_qp](_rz_radial_coord, _rz_radial_coord) / r),
-_mu[_qp] / r * (_grad_velocity[_qp](1, 0) + _grad_velocity[_qp](0, 1)),
0);
}
}

0 comments on commit 33a13dc

Please sign in to comment.