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        4370 : INSADTauMaterialTempl<T>::validParams()
     143             : {
     144        4370 :   InputParameters params = T::validParams();
     145        4370 :   params.addClassDescription(
     146             :       "This is the material class used to compute the stabilization parameter tau.");
     147        8740 :   params.addParam<Real>("alpha", 1., "Multiplicative factor on the stabilization parameter tau.");
     148        4370 :   return params;
     149           0 : }
     150             : 
     151             : template <typename T>
     152        3417 : INSADTauMaterialTempl<T>::INSADTauMaterialTempl(const InputParameters & parameters)
     153             :   : T(parameters),
     154        3417 :     _alpha(this->template getParam<Real>("alpha")),
     155        3417 :     _tau(this->template declareADProperty<Real>("tau")),
     156        3417 :     _viscous_strong_residual(
     157        3417 :         this->template declareADProperty<RealVectorValue>("viscous_strong_residual")),
     158        3417 :     _momentum_strong_residual(
     159        6834 :         this->template declareADProperty<RealVectorValue>("momentum_strong_residual")),
     160        3417 :     _velocity_var(getVectorVar("velocity", 0)),
     161        3417 :     _scalar_lagrange_fe(
     162        3417 :         _assembly.getFE(FEType(_velocity_var->feType().order, LAGRANGE), _mesh.dimension())),
     163        3417 :     _vel_number(_velocity_var->number()),
     164        6834 :     _vel_sys_number(_velocity_var->sys().number())
     165             : {
     166        3417 :   _scalar_lagrange_fe->get_d2phi();
     167        3417 : }
     168             : 
     169             : template <typename T>
     170             : bool
     171    44754628 : INSADTauMaterialTempl<T>::doVelocityDerivatives() const
     172             : {
     173    51387120 :   return ADReal::do_derivatives &&
     174     6632492 :          (_vel_sys_number == _fe_problem.currentNonlinearSystem().number());
     175             : }
     176             : 
     177             : template <typename T>
     178             : void
     179     5682135 : INSADTauMaterialTempl<T>::computeHMax()
     180             : {
     181     5682135 :   if (_disp_x_num == libMesh::invalid_uint || !ADReal::do_derivatives)
     182             :   {
     183     5623776 :     _hmax = _current_elem->hmax();
     184     5623776 :     return;
     185             :   }
     186             : 
     187       58359 :   _hmax = 0;
     188       58359 :   std::array<unsigned int, 3> disps = {_disp_x_num, _disp_y_num, _disp_z_num};
     189       58359 :   std::array<unsigned int, 3> disp_sys_nums = {_disp_x_sys_num, _disp_y_sys_num, _disp_z_sys_num};
     190             : 
     191      351091 :   for (unsigned int n_outer = 0; n_outer < _current_elem->n_vertices(); n_outer++)
     192      969014 :     for (unsigned int n_inner = n_outer + 1; n_inner < _current_elem->n_vertices(); n_inner++)
     193             :     {
     194      676282 :       VectorValue<ADReal> diff = (_current_elem->point(n_outer) - _current_elem->point(n_inner));
     195     2705128 :       for (const auto i : index_range(disps))
     196             :       {
     197     2028846 :         const auto disp_num = disps[i];
     198     2028846 :         if (disp_num == libMesh::invalid_uint)
     199      261210 :           continue;
     200     1767636 :         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     1767636 :         diff(i).derivatives().insert(
     206     1767636 :             _current_elem->node_ref(n_outer).dof_number(sys_num, disp_num, 0)) = 1.;
     207     1767636 :         diff(i).derivatives().insert(
     208     1767636 :             _current_elem->node_ref(n_inner).dof_number(sys_num, disp_num, 0)) = -1.;
     209             :       }
     210             : 
     211      676282 :       _hmax = std::max(_hmax, diff.norm_sq());
     212             :     }
     213             : 
     214       58359 :   _hmax = std::sqrt(_hmax);
     215             : }
     216             : 
     217             : template <typename T>
     218             : void
     219     5682135 : INSADTauMaterialTempl<T>::computeProperties()
     220             : {
     221     5682135 :   computeHMax();
     222     5682135 :   computeViscousStrongResidual();
     223             : 
     224     5682135 :   T::computeProperties();
     225     5682135 : }
     226             : 
     227             : template <typename T>
     228             : void
     229     5682135 : INSADTauMaterialTempl<T>::computeViscousStrongResidual()
     230             : {
     231    39774945 :   auto resize_and_zero = [this](auto & d2vel)
     232             :   {
     233    17046405 :     d2vel.resize(_qrule->n_points());
     234    86552499 :     for (auto & d2qp : d2vel)
     235             :       d2qp = 0;
     236             :   };
     237     5682135 :   resize_and_zero(_d2u);
     238     5682135 :   resize_and_zero(_d2v);
     239     5682135 :   resize_and_zero(_d2w);
     240             : 
     241    50436763 :   auto get_d2 = [this](const auto i) -> std::vector<ADRealTensorValue> &
     242             :   {
     243    44754628 :     switch (i)
     244             :     {
     245    22199390 :       case 0:
     246    22199390 :         return _d2u;
     247    22199390 :       case 1:
     248    22199390 :         return _d2v;
     249      355848 :       case 2:
     250      355848 :         return _d2w;
     251           0 :       default:
     252           0 :         mooseError("invalid value of 'i'");
     253             :     }
     254             :   };
     255             : 
     256             :   // libMesh does not yet have the capability for computing second order spatial derivatives of
     257             :   // vector bases. Lacking that capability, we can compute the second order spatial derivatives "by
     258             :   // hand" using the scalar field version of the vector basis, e.g. LAGRANGE instead of
     259             :   // LAGRANGE_VEC. Adding this implementation allows us to be fully consistent with results from a
     260             :   // scalar velocity field component implementation of Navier-Stokes
     261     5682135 :   const auto & vel_dof_indices = _velocity_var->dofIndices();
     262    50436763 :   for (const auto i : index_range(vel_dof_indices))
     263             :   {
     264             :     // This may not work if the element and spatial dimensions are different
     265             :     mooseAssert(_current_elem->dim() == _mesh.dimension(),
     266             :                 "Below logic only applicable if element and mesh dimension are the same");
     267    44754628 :     const auto dimensional_component = i % _mesh.dimension();
     268    44754628 :     auto & d2vel = get_d2(dimensional_component);
     269    44754628 :     const auto dof_index = vel_dof_indices[i];
     270    44754628 :     ADReal dof_value = (*_velocity_var->sys().currentSolution())(dof_index);
     271    44754628 :     if (doVelocityDerivatives())
     272     6632492 :       dof_value.derivatives().insert(dof_index) = 1;
     273    44754628 :     const auto scalar_i_component = i / _mesh.dimension();
     274   232524896 :     for (const auto qp : make_range(_qrule->n_points()))
     275   375540536 :       d2vel[qp] += dof_value * _scalar_lagrange_fe->get_d2phi()[scalar_i_component][qp];
     276             :   }
     277             : 
     278             :   // Now that we have the second order spatial derivatives of velocity, we can compute the strong
     279             :   // form of the viscous residual for use in our stabilization calculations
     280    28850833 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
     281             :   {
     282    46337396 :     _viscous_strong_residual[_qp](0) = -_mu[_qp] * _d2u[_qp].tr();
     283    23168698 :     _viscous_strong_residual[_qp](1) =
     284    69506094 :         _mesh.dimension() >= 2 ? -_mu[_qp] * _d2v[_qp].tr() : ADReal(0);
     285    23168698 :     _viscous_strong_residual[_qp](2) =
     286    46619524 :         _mesh.dimension() == 3 ? -_mu[_qp] * _d2w[_qp].tr() : ADReal(0);
     287    23168698 :     if (_viscous_form == NS::ViscousForm::Traction)
     288             :     {
     289     4717626 :       _viscous_strong_residual[_qp] -= _mu[_qp] * _d2u[_qp].row(0);
     290     4717626 :       if (_mesh.dimension() >= 2)
     291    14152878 :         _viscous_strong_residual[_qp] -= _mu[_qp] * _d2v[_qp].row(1);
     292     4717626 :       if (_mesh.dimension() == 3)
     293           0 :         _viscous_strong_residual[_qp] -= _mu[_qp] * _d2w[_qp].row(2);
     294             :     }
     295    23168698 :     if (_coord_sys == Moose::COORD_RZ)
     296     5250618 :       viscousTermRZ();
     297             :   }
     298     5682135 : }
     299             : 
     300             : template <typename T>
     301             : void
     302     5250618 : INSADTauMaterialTempl<T>::viscousTermRZ()
     303             : {
     304             :   // To understand the code immediately below, visit
     305             :   // https://en.wikipedia.org/wiki/Del_in_cylindrical_and_spherical_coordinates.
     306             :   // The u_r / r^2 term comes from the vector Laplacian. The -du_i/dr * 1/r term comes from
     307             :   // the scalar Laplacian. The scalar Laplacian in axisymmetric cylindrical coordinates is
     308             :   // equivalent to the Cartesian Laplacian plus a 1/r * du_i/dr term. And of course we are
     309             :   // applying a minus sign here because the strong form is -\nabala^2 * \vec{u}
     310             :   //
     311             :   // Another note: libMesh implements grad(v) as dvi/dxj or more obviously:
     312             :   //
     313             :   // grad(v) = (vx_x vx_y)
     314             :   //           (vy_x vy_y)
     315             :   //
     316             :   // so, the gradient of the velocity with respect to the radial coordinate will correspond to a
     317             :   // *column* slice
     318             : 
     319     5250618 :   const auto & r = _ad_q_point[_qp](_rz_radial_coord);
     320             : 
     321             :   {
     322             :     // Do the "Laplace" form. This will be present in *both* Laplace and Traction forms
     323             :     ADRealVectorValue rz_term;
     324    15751854 :     for (const auto i : make_range((unsigned int)2))
     325    21002472 :       rz_term(i) = -_mu[_qp] * _grad_velocity[_qp](i, _rz_radial_coord) / r;
     326    10501236 :     rz_term(_rz_radial_coord) += _mu[_qp] * _velocity[_qp](_rz_radial_coord) / (r * r);
     327     5250618 :     _viscous_strong_residual[_qp] += rz_term;
     328             :   }
     329     5250618 :   if (_viscous_form == NS::ViscousForm::Traction)
     330             :   {
     331             :     ADRealVectorValue rz_term;
     332    13642578 :     for (const auto i : make_range((unsigned int)2))
     333             :       // This is the transpose of the above
     334    18190104 :       rz_term(i) = -_mu[_qp] * _grad_velocity[_qp](_rz_radial_coord, i) / r;
     335             :     // This is the same as above (since the transpose of the diagonal is the diagonal)
     336     9095052 :     rz_term(_rz_radial_coord) += _mu[_qp] * _velocity[_qp](_rz_radial_coord) / (r * r);
     337     4547526 :     _viscous_strong_residual[_qp] += rz_term;
     338             :   }
     339     5250618 : }
     340             : 
     341             : template <typename T>
     342             : void
     343    23168698 : INSADTauMaterialTempl<T>::computeQpProperties()
     344             : {
     345    23168698 :   T::computeQpProperties();
     346             : 
     347    23168698 :   const auto nu = _mu[_qp] / _rho[_qp];
     348    23168698 :   const auto transient_part = _has_transient ? 4. / (_dt * _dt) : 0.;
     349    23168698 :   _speed = NS::computeSpeed<ADReal>(_relative_velocity[_qp]);
     350    92674792 :   _tau[_qp] = _alpha / std::sqrt(transient_part + (2. * _speed / _hmax) * (2. * _speed / _hmax) +
     351    92674792 :                                  9. * (4. * nu / (_hmax * _hmax)) * (4. * nu / (_hmax * _hmax)));
     352             : 
     353    46337396 :   _momentum_strong_residual[_qp] =
     354    23168698 :       _advective_strong_residual[_qp] + _viscous_strong_residual[_qp] + _grad_p[_qp];
     355             : 
     356    23168698 :   if (_has_transient)
     357     2190754 :     _momentum_strong_residual[_qp] += _td_strong_residual[_qp];
     358             : 
     359    23168698 :   if (_has_gravity)
     360    11759808 :     _momentum_strong_residual[_qp] += _gravity_strong_residual[_qp];
     361             : 
     362    23168698 :   if (_has_boussinesq)
     363    11695296 :     _momentum_strong_residual[_qp] += _boussinesq_strong_residual[_qp];
     364             : 
     365    23168698 :   if (_has_advected_mesh)
     366      690184 :     _momentum_strong_residual[_qp] += _advected_mesh_strong_residual[_qp];
     367             : 
     368    23168698 :   if (_has_coupled_force)
     369     1668432 :     _momentum_strong_residual[_qp] += _coupled_force_strong_residual[_qp];
     370             : 
     371             :   // // Future addition
     372             :   // if (_object_tracker->hasMMS())
     373             :   //   _momentum_strong_residual[_qp] += _mms_function_strong_residual[_qp];
     374    23168698 : }