diff --git a/modules/navier_stokes/include/materials/INSADMaterial.h b/modules/navier_stokes/include/materials/INSADMaterial.h index db2e69094256..8bf6ce46fcc2 100644 --- a/modules/navier_stokes/include/materials/INSADMaterial.h +++ b/modules/navier_stokes/include/materials/INSADMaterial.h @@ -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; @@ -54,10 +49,6 @@ class INSADMaterial : public Material /// Strong residual corresponding to the momentum advective term ADMaterialProperty & _advective_strong_residual; - /// Strong residual corresponding to the momentum viscous term. This is only used by stabilization - /// kernels - ADMaterialProperty & _viscous_strong_residual; - /// Strong residual corresponding to the momentum transient term ADMaterialProperty & _td_strong_residual; diff --git a/modules/navier_stokes/include/materials/INSADTauMaterial.h b/modules/navier_stokes/include/materials/INSADTauMaterial.h index 8c2629e1c5fb..1b8cca108a4b 100644 --- a/modules/navier_stokes/include/materials/INSADTauMaterial.h +++ b/modules/navier_stokes/include/materials/INSADTauMaterial.h @@ -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 @@ -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 & _tau; + /// Strong residual corresponding to the momentum viscous term. This is only used by stabilization + /// kernels + ADMaterialProperty & _viscous_strong_residual; + /// The strong residual of the momentum equation ADMaterialProperty & _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 _d2u; + std::vector _d2v; + std::vector _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; @@ -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 INSADTauMaterial; @@ -93,9 +148,25 @@ INSADTauMaterialTempl::INSADTauMaterialTempl(const InputParameters & paramete : T(parameters), _alpha(this->template getParam("alpha")), _tau(this->template declareADProperty("tau")), + _viscous_strong_residual( + this->template declareADProperty("viscous_strong_residual")), _momentum_strong_residual( - this->template declareADProperty("momentum_strong_residual")) + this->template declareADProperty("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 +bool +INSADTauMaterialTempl::doVelocityDerivatives() const +{ + return ADReal::do_derivatives && + (_vel_sys_number == _fe_problem.currentNonlinearSystem().number()); } template @@ -140,10 +211,139 @@ void INSADTauMaterialTempl::computeProperties() { computeHMax(); + computeViscousStrongResidual(); T::computeProperties(); } +template +void +INSADTauMaterialTempl::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 & + { + 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 +void +INSADTauMaterialTempl::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 void INSADTauMaterialTempl::computeQpProperties() @@ -156,12 +356,8 @@ INSADTauMaterialTempl::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]; diff --git a/modules/navier_stokes/src/materials/INSADMaterial.C b/modules/navier_stokes/src/materials/INSADMaterial.C index 688b1b2ade2a..26c395d0085e 100644 --- a/modules/navier_stokes/src/materials/INSADMaterial.C +++ b/modules/navier_stokes/src/materials/INSADMaterial.C @@ -40,7 +40,6 @@ INSADMaterial::INSADMaterial(const InputParameters & parameters) _velocity_dot(nullptr), _mass_strong_residual(declareADProperty("mass_strong_residual")), _advective_strong_residual(declareADProperty("advective_strong_residual")), - _viscous_strong_residual(declareADProperty("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 @@ -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); - } }