Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #pragma once
      11             : 
      12             : #include "InputParameters.h"
      13             : #include "NonlinearSystemBase.h"
      14             : #include "FEProblemBase.h"
      15             : #include "MaterialProperty.h"
      16             : #include "MooseArray.h"
      17             : #include "INSADMaterial.h"
      18             : #include "NavierStokesMethods.h"
      19             : #include "Assembly.h"
      20             : #include "MooseVariableFE.h"
      21             : #include "MooseMesh.h"
      22             : 
      23             : #include "libmesh/elem.h"
      24             : #include "libmesh/node.h"
      25             : #include "libmesh/fe_type.h"
      26             : 
      27             : #include <vector>
      28             : 
      29             : class INSADMaterial;
      30             : 
      31             : template <typename T>
      32             : class INSADTauMaterialTempl : public T
      33             : {
      34             : public:
      35             :   static InputParameters validParams();
      36             : 
      37             :   INSADTauMaterialTempl(const InputParameters & parameters);
      38             : 
      39             : protected:
      40             :   virtual void computeProperties() override;
      41             :   virtual void computeQpProperties() override;
      42             : 
      43             :   /**
      44             :    * Compute the maximum dimension of the current element
      45             :    */
      46             :   void computeHMax();
      47             : 
      48             :   /**
      49             :    * Compute the viscous strong residual
      50             :    */
      51             :   void computeViscousStrongResidual();
      52             : 
      53             :   /**
      54             :    * Compute the strong form corresponding to RZ pieces of the viscous term
      55             :    */
      56             :   void viscousTermRZ();
      57             : 
      58             :   /**
      59             :    * Whether to seed with respect to velocity derivatives
      60             :    */
      61             :   bool doVelocityDerivatives() const;
      62             : 
      63             :   const Real _alpha;
      64             :   ADMaterialProperty<Real> & _tau;
      65             : 
      66             :   /// Strong residual corresponding to the momentum viscous term. This is only used by stabilization
      67             :   /// kernels
      68             :   ADMaterialProperty<RealVectorValue> & _viscous_strong_residual;
      69             : 
      70             :   /// The strong residual of the momentum equation
      71             :   ADMaterialProperty<RealVectorValue> & _momentum_strong_residual;
      72             : 
      73             :   ADReal _hmax;
      74             : 
      75             :   /// The velocity variable
      76             :   const VectorMooseVariable * const _velocity_var;
      77             : 
      78             :   /// A scalar Lagrange FE data member to compute the velocity second derivatives since
      79             :   /// they're currently not supported for vector FE types
      80             :   const FEBase * const & _scalar_lagrange_fe;
      81             : 
      82             :   /// Containers to hold the matrix of second spatial derivatives of velocity
      83             :   std::vector<ADRealTensorValue> _d2u;
      84             :   std::vector<ADRealTensorValue> _d2v;
      85             :   std::vector<ADRealTensorValue> _d2w;
      86             : 
      87             :   /// The velocity variable number
      88             :   const unsigned int _vel_number;
      89             : 
      90             :   /// The velocity system number
      91             :   const unsigned int _vel_sys_number;
      92             : 
      93             :   /// The speed of the medium. This is the norm of the relative velocity, e.g. the velocity minus
      94             :   /// the mesh velocity, at the current _qp
      95             :   ADReal _speed;
      96             : 
      97             :   using T::_ad_q_point;
      98             :   using T::_advected_mesh_strong_residual;
      99             :   using T::_advective_strong_residual;
     100             :   using T::_assembly;
     101             :   using T::_boussinesq_strong_residual;
     102             :   using T::_coord_sys;
     103             :   using T::_coupled_force_strong_residual;
     104             :   using T::_current_elem;
     105             :   using T::_disp_x_num;
     106             :   using T::_disp_x_sys_num;
     107             :   using T::_disp_y_num;
     108             :   using T::_disp_y_sys_num;
     109             :   using T::_disp_z_num;
     110             :   using T::_disp_z_sys_num;
     111             :   using T::_dt;
     112             :   using T::_fe_problem;
     113             :   using T::_grad_p;
     114             :   using T::_grad_velocity;
     115             :   using T::_gravity_strong_residual;
     116             :   using T::_has_advected_mesh;
     117             :   using T::_has_boussinesq;
     118             :   using T::_has_coupled_force;
     119             :   using T::_has_gravity;
     120             :   using T::_has_transient;
     121             :   using T::_mesh;
     122             :   using T::_mu;
     123             :   using T::_object_tracker;
     124             :   using T::_q_point;
     125             :   using T::_qp;
     126             :   using T::_qrule;
     127             :   using T::_relative_velocity;
     128             :   using T::_rho;
     129             :   using T::_rz_axial_coord;
     130             :   using T::_rz_radial_coord;
     131             :   using T::_td_strong_residual;
     132             :   using T::_use_displaced_mesh;
     133             :   using T::_velocity;
     134             :   using T::_viscous_form;
     135             :   using T::getVectorVar;
     136             : };
     137             : 
     138             : typedef INSADTauMaterialTempl<INSADMaterial> INSADTauMaterial;
     139             : 
     140             : template <typename T>
     141             : InputParameters
     142        2150 : INSADTauMaterialTempl<T>::validParams()
     143             : {
     144        2150 :   InputParameters params = T::validParams();
     145        2150 :   params.addClassDescription(
     146             :       "This is the material class used to compute the stabilization parameter tau.");
     147        4300 :   params.addParam<Real>("alpha", 1., "Multiplicative factor on the stabilization parameter tau.");
     148        2150 :   return params;
     149           0 : }
     150             : 
     151             : template <typename T>
     152        1665 : INSADTauMaterialTempl<T>::INSADTauMaterialTempl(const InputParameters & parameters)
     153             :   : T(parameters),
     154        1665 :     _alpha(this->template getParam<Real>("alpha")),
     155        1665 :     _tau(this->template declareADProperty<Real>("tau")),
     156        1665 :     _viscous_strong_residual(
     157        1665 :         this->template declareADProperty<RealVectorValue>("viscous_strong_residual")),
     158        1665 :     _momentum_strong_residual(
     159        3330 :         this->template declareADProperty<RealVectorValue>("momentum_strong_residual")),
     160        1665 :     _velocity_var(getVectorVar("velocity", 0)),
     161        1665 :     _scalar_lagrange_fe(
     162        1665 :         _assembly.getFE(FEType(_velocity_var->feType().order, LAGRANGE), _mesh.dimension())),
     163        1665 :     _vel_number(_velocity_var->number()),
     164        3330 :     _vel_sys_number(_velocity_var->sys().number())
     165             : {
     166        1665 :   _scalar_lagrange_fe->get_d2phi();
     167        1665 : }
     168             : 
     169             : template <typename T>
     170             : bool
     171    34594080 : INSADTauMaterialTempl<T>::doVelocityDerivatives() const
     172             : {
     173    39357268 :   return ADReal::do_derivatives &&
     174     4763188 :          (_vel_sys_number == _fe_problem.currentNonlinearSystem().number());
     175             : }
     176             : 
     177             : template <typename T>
     178             : void
     179     4414333 : INSADTauMaterialTempl<T>::computeHMax()
     180             : {
     181     4414333 :   if (_disp_x_num == libMesh::invalid_uint || !ADReal::do_derivatives)
     182             :   {
     183     4375078 :     _hmax = _current_elem->hmax();
     184     4375078 :     return;
     185             :   }
     186             : 
     187       39255 :   _hmax = 0;
     188       39255 :   std::array<unsigned int, 3> disps = {_disp_x_num, _disp_y_num, _disp_z_num};
     189       39255 :   std::array<unsigned int, 3> disp_sys_nums = {_disp_x_sys_num, _disp_y_sys_num, _disp_z_sys_num};
     190             : 
     191      235963 :   for (unsigned int n_outer = 0; n_outer < _current_elem->n_vertices(); n_outer++)
     192      650522 :     for (unsigned int n_inner = n_outer + 1; n_inner < _current_elem->n_vertices(); n_inner++)
     193             :     {
     194      453814 :       VectorValue<ADReal> diff = (_current_elem->point(n_outer) - _current_elem->point(n_inner));
     195     1815256 :       for (const auto i : index_range(disps))
     196             :       {
     197     1361442 :         const auto disp_num = disps[i];
     198     1361442 :         if (disp_num == libMesh::invalid_uint)
     199      175998 :           continue;
     200     1185444 :         const auto sys_num = disp_sys_nums[i];
     201             : 
     202             :         // Here we insert derivatives of the difference in nodal positions with respect to the
     203             :         // displacement degrees of freedom. From above, diff = outer_node_position -
     204             :         // inner_node_position
     205     1185444 :         diff(i).derivatives().insert(
     206     1185444 :             _current_elem->node_ref(n_outer).dof_number(sys_num, disp_num, 0)) = 1.;
     207     1185444 :         diff(i).derivatives().insert(
     208     1185444 :             _current_elem->node_ref(n_inner).dof_number(sys_num, disp_num, 0)) = -1.;
     209             :       }
     210             : 
     211      907628 :       _hmax = std::max(_hmax, diff.norm_sq());
     212             :     }
     213             : 
     214             :   using std::sqrt;
     215       39255 :   _hmax = sqrt(_hmax);
     216             : }
     217             : 
     218             : template <typename T>
     219             : void
     220     4414333 : INSADTauMaterialTempl<T>::computeProperties()
     221             : {
     222     4414333 :   computeHMax();
     223     4414333 :   computeViscousStrongResidual();
     224             : 
     225     4414333 :   T::computeProperties();
     226     4414333 : }
     227             : 
     228             : template <typename T>
     229             : void
     230     4414333 : INSADTauMaterialTempl<T>::computeViscousStrongResidual()
     231             : {
     232    30900331 :   auto resize_and_zero = [this](auto & d2vel)
     233             :   {
     234    13242999 :     d2vel.resize(_qrule->n_points());
     235    67120497 :     for (auto & d2qp : d2vel)
     236             :       d2qp = 0;
     237             :   };
     238     4414333 :   resize_and_zero(_d2u);
     239     4414333 :   resize_and_zero(_d2v);
     240     4414333 :   resize_and_zero(_d2w);
     241             : 
     242    39008413 :   auto get_d2 = [this](const auto i) -> std::vector<ADRealTensorValue> &
     243             :   {
     244    34594080 :     switch (i)
     245             :     {
     246    17178096 :       case 0:
     247    17178096 :         return _d2u;
     248    17178096 :       case 1:
     249    17178096 :         return _d2v;
     250      237888 :       case 2:
     251      237888 :         return _d2w;
     252           0 :       default:
     253           0 :         mooseError("invalid value of 'i'");
     254             :     }
     255             :   };
     256             : 
     257             :   // libMesh does not yet have the capability for computing second order spatial derivatives of
     258             :   // vector bases. Lacking that capability, we can compute the second order spatial derivatives "by
     259             :   // hand" using the scalar field version of the vector basis, e.g. LAGRANGE instead of
     260             :   // LAGRANGE_VEC. Adding this implementation allows us to be fully consistent with results from a
     261             :   // scalar velocity field component implementation of Navier-Stokes
     262     4414333 :   const auto & vel_dof_indices = _velocity_var->dofIndices();
     263    39008413 :   for (const auto i : index_range(vel_dof_indices))
     264             :   {
     265             :     // This may not work if the element and spatial dimensions are different
     266             :     mooseAssert(_current_elem->dim() == _mesh.dimension(),
     267             :                 "Below logic only applicable if element and mesh dimension are the same");
     268    34594080 :     const auto dimensional_component = i % _mesh.dimension();
     269    34594080 :     auto & d2vel = get_d2(dimensional_component);
     270    34594080 :     const auto dof_index = vel_dof_indices[i];
     271    34594080 :     ADReal dof_value = (*_velocity_var->sys().currentSolution())(dof_index);
     272    34594080 :     if (doVelocityDerivatives())
     273     4763188 :       dof_value.derivatives().insert(dof_index) = 1;
     274    34594080 :     const auto scalar_i_component = i / _mesh.dimension();
     275   178959384 :     for (const auto qp : make_range(_qrule->n_points()))
     276   288730608 :       d2vel[qp] += dof_value * _scalar_lagrange_fe->get_d2phi()[scalar_i_component][qp];
     277             :   }
     278             : 
     279             :   // Now that we have the second order spatial derivatives of velocity, we can compute the strong
     280             :   // form of the viscous residual for use in our stabilization calculations
     281    22373499 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
     282             :   {
     283    35918332 :     _viscous_strong_residual[_qp](0) = -_mu[_qp] * _d2u[_qp].tr();
     284    17959166 :     _viscous_strong_residual[_qp](1) =
     285    53877498 :         _mesh.dimension() >= 2 ? -_mu[_qp] * _d2v[_qp].tr() : ADReal(0);
     286    17959166 :     _viscous_strong_residual[_qp](2) =
     287    36106980 :         _mesh.dimension() == 3 ? -_mu[_qp] * _d2w[_qp].tr() : ADReal(0);
     288    17959166 :     if (_viscous_form == NS::ViscousForm::Traction)
     289             :     {
     290     3698352 :       _viscous_strong_residual[_qp] -= _mu[_qp] * _d2u[_qp].row(0);
     291     3698352 :       if (_mesh.dimension() >= 2)
     292    11095056 :         _viscous_strong_residual[_qp] -= _mu[_qp] * _d2v[_qp].row(1);
     293     3698352 :       if (_mesh.dimension() == 3)
     294           0 :         _viscous_strong_residual[_qp] -= _mu[_qp] * _d2w[_qp].row(2);
     295             :     }
     296    17959166 :     if (_coord_sys == Moose::COORD_RZ)
     297     4069452 :       viscousTermRZ();
     298             :   }
     299     4414333 : }
     300             : 
     301             : template <typename T>
     302             : void
     303     4069452 : INSADTauMaterialTempl<T>::viscousTermRZ()
     304             : {
     305             :   // To understand the code immediately below, visit
     306             :   // https://en.wikipedia.org/wiki/Del_in_cylindrical_and_spherical_coordinates.
     307             :   // The u_r / r^2 term comes from the vector Laplacian. The -du_i/dr * 1/r term comes from
     308             :   // the scalar Laplacian. The scalar Laplacian in axisymmetric cylindrical coordinates is
     309             :   // equivalent to the Cartesian Laplacian plus a 1/r * du_i/dr term. And of course we are
     310             :   // applying a minus sign here because the strong form is -\nabala^2 * \vec{u}
     311             :   //
     312             :   // Another note: libMesh implements grad(v) as dvi/dxj or more obviously:
     313             :   //
     314             :   // grad(v) = (vx_x vx_y)
     315             :   //           (vy_x vy_y)
     316             :   //
     317             :   // so, the gradient of the velocity with respect to the radial coordinate will correspond to a
     318             :   // *column* slice
     319             : 
     320     4069452 :   const auto & r = _ad_q_point[_qp](_rz_radial_coord);
     321             : 
     322             :   {
     323             :     // Do the "Laplace" form. This will be present in *both* Laplace and Traction forms
     324             :     ADRealVectorValue rz_term;
     325    12208356 :     for (const auto i : make_range((unsigned int)2))
     326    16277808 :       rz_term(i) = -_mu[_qp] * _grad_velocity[_qp](i, _rz_radial_coord) / r;
     327     8138904 :     rz_term(_rz_radial_coord) += _mu[_qp] * _velocity[_qp](_rz_radial_coord) / (r * r);
     328     4069452 :     _viscous_strong_residual[_qp] += rz_term;
     329             :   }
     330     4069452 :   if (_viscous_form == NS::ViscousForm::Traction)
     331             :   {
     332             :     ADRealVectorValue rz_term;
     333    10754856 :     for (const auto i : make_range((unsigned int)2))
     334             :       // This is the transpose of the above
     335    14339808 :       rz_term(i) = -_mu[_qp] * _grad_velocity[_qp](_rz_radial_coord, i) / r;
     336             :     // This is the same as above (since the transpose of the diagonal is the diagonal)
     337     7169904 :     rz_term(_rz_radial_coord) += _mu[_qp] * _velocity[_qp](_rz_radial_coord) / (r * r);
     338     3584952 :     _viscous_strong_residual[_qp] += rz_term;
     339             :   }
     340     4069452 : }
     341             : 
     342             : template <typename T>
     343             : void
     344    17959166 : INSADTauMaterialTempl<T>::computeQpProperties()
     345             : {
     346    17959166 :   T::computeQpProperties();
     347             : 
     348             :   using std::sqrt;
     349             : 
     350    17959166 :   const auto nu = _mu[_qp] / _rho[_qp];
     351    17959166 :   const auto transient_part = _has_transient ? 4. / (_dt * _dt) : 0.;
     352    17959166 :   _speed = NS::computeSpeed<ADReal>(_relative_velocity[_qp]);
     353    71836664 :   _tau[_qp] = _alpha / sqrt(transient_part + (2. * _speed / _hmax) * (2. * _speed / _hmax) +
     354    71836664 :                             9. * (4. * nu / (_hmax * _hmax)) * (4. * nu / (_hmax * _hmax)));
     355             : 
     356    35918332 :   _momentum_strong_residual[_qp] =
     357    17959166 :       _advective_strong_residual[_qp] + _viscous_strong_residual[_qp] + _grad_p[_qp];
     358             : 
     359    17959166 :   if (_has_transient)
     360     1483634 :     _momentum_strong_residual[_qp] += _td_strong_residual[_qp];
     361             : 
     362    17959166 :   if (_has_gravity)
     363     9439296 :     _momentum_strong_residual[_qp] += _gravity_strong_residual[_qp];
     364             : 
     365    17959166 :   if (_has_boussinesq)
     366     9396288 :     _momentum_strong_residual[_qp] += _boussinesq_strong_residual[_qp];
     367             : 
     368    17959166 :   if (_has_advected_mesh)
     369      463018 :     _momentum_strong_residual[_qp] += _advected_mesh_strong_residual[_qp];
     370             : 
     371    17959166 :   if (_has_coupled_force)
     372     1210688 :     _momentum_strong_residual[_qp] += _coupled_force_strong_residual[_qp];
     373             : 
     374             :   // // Future addition
     375             :   // if (_object_tracker->hasMMS())
     376             :   //   _momentum_strong_residual[_qp] += _mms_function_strong_residual[_qp];
     377    17959166 : }