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             : #include "INSADMomentumViscous.h"
      11             : #include "Assembly.h"
      12             : #include "SystemBase.h"
      13             : #include "INSADObjectTracker.h"
      14             : 
      15             : #include "metaphysicl/raw_type.h"
      16             : 
      17             : using MetaPhysicL::raw_value;
      18             : 
      19             : registerMooseObject("NavierStokesApp", INSADMomentumViscous);
      20             : 
      21             : InputParameters
      22        2896 : INSADMomentumViscous::validParams()
      23             : {
      24        2896 :   InputParameters params = ADVectorKernel::validParams();
      25        2896 :   params.addClassDescription("Adds the viscous term to the INS momentum equation");
      26        5792 :   params.addParam<MaterialPropertyName>(
      27             :       "mu_name", "mu", "The name of the viscosity material property");
      28        5792 :   MooseEnum viscous_form("traction laplace", "laplace");
      29        5792 :   params.addParam<MooseEnum>("viscous_form",
      30             :                              viscous_form,
      31             :                              "The form of the viscous term. Options are 'traction' or 'laplace'");
      32        2896 :   return params;
      33        2896 : }
      34             : 
      35        1577 : INSADMomentumViscous::INSADMomentumViscous(const InputParameters & parameters)
      36             :   : ADVectorKernel(parameters),
      37        1577 :     _mu(getADMaterialProperty<Real>("mu_name")),
      38        1577 :     _coord_sys(_assembly.coordSystem()),
      39        3154 :     _form(getParam<MooseEnum>("viscous_form")),
      40        3154 :     _rz_radial_coord(_mesh.getAxisymmetricRadialCoord())
      41             : {
      42             :   auto & obj_tracker = const_cast<INSADObjectTracker &>(
      43        1577 :       _fe_problem.getUserObject<INSADObjectTracker>("ins_ad_object_tracker"));
      44        3176 :   for (const auto block_id : blockIDs())
      45        3198 :     obj_tracker.set("viscous_form", _form, block_id);
      46        1577 : }
      47             : 
      48             : ADRealTensorValue
      49    31269077 : INSADMomentumViscous::qpViscousTerm()
      50             : {
      51    31269077 :   if (_form == "laplace")
      52    52348414 :     return _mu[_qp] * _grad_u[_qp];
      53             :   else
      54    10189740 :     return _mu[_qp] * (_grad_u[_qp] + _grad_u[_qp].transpose());
      55             : }
      56             : 
      57             : ADRealVectorValue
      58     6328092 : INSADMomentumViscous::qpAdditionalRZTerm()
      59             : {
      60             :   ADRealVectorValue ret;
      61     6328092 :   auto & extra_term = ret(_rz_radial_coord);
      62             : 
      63             :   // Add the u_r / r^2 term. There is an extra factor of 2 for the traction form
      64     6328092 :   const auto & r = _ad_q_point[_qp](_rz_radial_coord);
      65    12656184 :   extra_term = _mu[_qp] * _u[_qp](_rz_radial_coord) / (r * r);
      66     6328092 :   if (_form == "traction")
      67     4924770 :     extra_term *= 2.;
      68             : 
      69     6328092 :   return ret;
      70             : }
      71             : 
      72             : void
      73     5516117 : INSADMomentumViscous::computeResidual()
      74             : {
      75             :   // When computing the residual we use the regular version of all quantities, e.g. JxW, coord,
      76             :   // test, grad_test, because we do not need to track the derivatives
      77             : 
      78     5516117 :   prepareVectorTag(_assembly, _var.number());
      79             : 
      80     5516117 :   const unsigned int n_test = _grad_test.size();
      81             : 
      82    31039668 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
      83             :   {
      84    25523551 :     const RealTensorValue value = raw_value(qpViscousTerm()) * _JxW[_qp] * _coord[_qp];
      85   281187725 :     for (_i = 0; _i < n_test; _i++) // target for auto vectorization
      86   511328348 :       _local_re(_i) += MathUtils::dotProduct(value, _regular_grad_test[_i][_qp]);
      87             : 
      88             :     // If we're in RZ, then we need to add an additional term
      89    25523551 :     if (_coord_sys == Moose::COORD_RZ)
      90             :     {
      91     5384034 :       const RealVectorValue rz_value = raw_value(qpAdditionalRZTerm()) * _JxW[_qp] * _coord[_qp];
      92             : 
      93    44933886 :       for (_i = 0; _i < n_test; _i++) // target for auto vectorization
      94    39549852 :         _local_re(_i) += rz_value * _test[_i][_qp];
      95             :     }
      96             :   }
      97             : 
      98     5516117 :   accumulateTaggedLocalResidual();
      99             : 
     100     5516117 :   if (_has_save_in)
     101             :   {
     102             :     Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
     103           0 :     for (unsigned int i = 0; i < _save_in.size(); i++)
     104           0 :       _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
     105             :   }
     106     5516117 : }
     107             : 
     108             : void
     109     1038649 : INSADMomentumViscous::computeResidualsForJacobian()
     110             : {
     111             :   mooseAssert(_test.size() == _grad_test.size(),
     112             :               "How are the there are different number of test and grad_test objects?");
     113     1038649 :   const auto n_test = _test.size();
     114             : 
     115     1038649 :   if (_residuals.size() != n_test)
     116        1302 :     _residuals.resize(n_test, 0);
     117    12358161 :   for (auto & r : _residuals)
     118    11319512 :     r = 0;
     119             : 
     120     1038649 :   if (_use_displaced_mesh)
     121      248560 :     for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     122             :     {
     123             :       // This will also compute the derivative with respect to all dofs
     124      204992 :       const ADRealTensorValue value = qpViscousTerm() * _ad_JxW[_qp] * _ad_coord[_qp];
     125     2732928 :       for (_i = 0; _i < _grad_test.size(); _i++)
     126     2527936 :         _residuals[_i] += MathUtils::dotProduct(value, _grad_test[_i][_qp]);
     127             : 
     128             :       // If we're in RZ, then we need to add an additional term
     129      204992 :       if (_coord_sys == Moose::COORD_RZ)
     130             :       {
     131       47520 :         const ADRealVectorValue rz_value = qpAdditionalRZTerm() * _ad_JxW[_qp] * _ad_coord[_qp];
     132             : 
     133      332640 :         for (_i = 0; _i < n_test; _i++)
     134      285120 :           _residuals[_i] += rz_value * _test[_i][_qp];
     135             :       }
     136             :     }
     137             :   else
     138     6535615 :     for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     139             :     {
     140             :       // Compute scalar quanitities up front to reduce DualNumber operations
     141    11081068 :       const ADRealTensorValue value = _JxW[_qp] * _coord[_qp] * qpViscousTerm();
     142    75608886 :       for (_i = 0; _i < _grad_test.size(); _i++)
     143    70068352 :         _residuals[_i] += MathUtils::dotProduct(value, _regular_grad_test[_i][_qp]);
     144             : 
     145             :       // If we're in RZ, then we need to add an additional term
     146     5540534 :       if (_coord_sys == Moose::COORD_RZ)
     147             :       {
     148             :         // Compute scalar quanitities up front to reduce DualNumber operations
     149     1793076 :         const ADRealVectorValue rz_value = _JxW[_qp] * _coord[_qp] * qpAdditionalRZTerm();
     150    11289882 :         for (_i = 0; _i < n_test; _i++)
     151    10393344 :           _residuals[_i] += rz_value * _test[_i][_qp];
     152             :       }
     153             :     }
     154     1038649 : }
     155             : 
     156             : ADReal
     157           0 : INSADMomentumViscous::computeQpResidual()
     158             : {
     159           0 :   mooseError("computeQpResidual is not used in the INSADMomentumViscous class");
     160             : }