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 "ConvectedMesh.h"
      11             : 
      12             : registerMooseObject("FsiApp", ConvectedMesh);
      13             : 
      14             : InputParameters
      15          38 : ConvectedMesh::validParams()
      16             : {
      17          38 :   InputParameters params = INSBase::validParams();
      18          38 :   params.addClassDescription(
      19             :       "Corrects the convective derivative for situations in which the fluid mesh is dynamic.");
      20          76 :   params.addRequiredCoupledVar("disp_x", "The x displacement");
      21          76 :   params.addCoupledVar("disp_y", "The y displacement");
      22          76 :   params.addCoupledVar("disp_z", "The z displacement");
      23          76 :   params.addParam<MaterialPropertyName>("rho_name", "rho", "The name of the density");
      24          76 :   params.addParam<bool>(
      25          76 :       "supg", false, "Whether to perform SUPG stabilization of the momentum residuals");
      26          38 :   return params;
      27           0 : }
      28             : 
      29          20 : ConvectedMesh::ConvectedMesh(const InputParameters & parameters)
      30             :   : INSBase(parameters),
      31          20 :     _disp_x_dot(coupledDot("disp_x")),
      32          20 :     _d_disp_x_dot(coupledDotDu("disp_x")),
      33          20 :     _disp_x_id(coupled("disp_x")),
      34          40 :     _disp_y_dot(isCoupled("disp_y") ? coupledDot("disp_y") : _zero),
      35          40 :     _d_disp_y_dot(isCoupled("disp_y") ? coupledDotDu("disp_y") : _zero),
      36          20 :     _disp_y_id(coupled("disp_y")),
      37          20 :     _disp_z_dot(isCoupled("disp_z") ? coupledDot("disp_z") : _zero),
      38          20 :     _d_disp_z_dot(isCoupled("disp_z") ? coupledDotDu("disp_z") : _zero),
      39          20 :     _disp_z_id(coupled("disp_z")),
      40          40 :     _rho(getMaterialProperty<Real>("rho_name")),
      41          60 :     _supg(getParam<bool>("supg"))
      42             : {
      43          20 :   if (_var.number() == _u_vel_var_number)
      44          10 :     _component = 0;
      45          10 :   else if (_var.number() == _v_vel_var_number)
      46          10 :     _component = 1;
      47           0 :   else if (_var.number() == _w_vel_var_number)
      48           0 :     _component = 2;
      49             :   else
      50           0 :     paramError("variable", "The variable must match one of the velocity variables.");
      51          20 : }
      52             : 
      53             : Real
      54    22494080 : ConvectedMesh::strongResidual()
      55             : {
      56    22494080 :   return -_rho[_qp] * RealVectorValue(_disp_x_dot[_qp], _disp_y_dot[_qp], _disp_z_dot[_qp]) *
      57    22494080 :          _grad_u[_qp];
      58             : }
      59             : 
      60             : Real
      61    21938560 : ConvectedMesh::computeQpResidual()
      62             : {
      63    21938560 :   auto test = _test[_i][_qp];
      64    21938560 :   const auto U = relativeVelocity();
      65    21938560 :   if (_supg)
      66    21938560 :     test += tau() * _grad_test[_i][_qp] * U;
      67    21938560 :   return test * strongResidual();
      68             : }
      69             : 
      70             : Real
      71      555520 : ConvectedMesh::computePGVelocityJacobian(const unsigned short component)
      72             : {
      73      555520 :   const auto U = relativeVelocity();
      74      555520 :   return strongResidual() * ((dTauDUComp(component) * _grad_test[_i][_qp] * U) +
      75      555520 :                              (tau() * _grad_test[_i][_qp](component) * _phi[_j][_qp]));
      76             : }
      77             : 
      78             : Real
      79      277760 : ConvectedMesh::computeQpJacobian()
      80             : {
      81      277760 :   auto test = _test[_i][_qp];
      82      277760 :   const auto U = relativeVelocity();
      83      277760 :   if (_supg)
      84      277760 :     test += tau() * _grad_test[_i][_qp] * U;
      85      277760 :   auto jac = test * -_rho[_qp] *
      86      277760 :              RealVectorValue(_disp_x_dot[_qp], _disp_y_dot[_qp], _disp_z_dot[_qp]) *
      87      277760 :              _grad_phi[_j][_qp];
      88      277760 :   if (_supg)
      89      277760 :     jac += computePGVelocityJacobian(_component);
      90             : 
      91      277760 :   return jac;
      92             : }
      93             : 
      94             : Real
      95     1111040 : ConvectedMesh::computeQpOffDiagJacobian(unsigned int jvar)
      96             : {
      97             :   mooseAssert(jvar != _var.number(), "Making sure I understand how old hand-coded Jacobians work.");
      98             : 
      99     1111040 :   auto test = _test[_i][_qp];
     100     1111040 :   const auto U = relativeVelocity();
     101     1111040 :   if (_supg)
     102     1111040 :     test += tau() * _grad_test[_i][_qp] * U;
     103             : 
     104     1111040 :   if (jvar == _disp_x_id)
     105      277760 :     return test * -_rho[_qp] * _phi[_j][_qp] * _d_disp_x_dot[_qp] * _grad_u[_qp](0);
     106      833280 :   else if (jvar == _disp_y_id)
     107      277760 :     return test * -_rho[_qp] * _phi[_j][_qp] * _d_disp_y_dot[_qp] * _grad_u[_qp](1);
     108      555520 :   else if (jvar == _disp_z_id)
     109           0 :     return test * -_rho[_qp] * _phi[_j][_qp] * _d_disp_z_dot[_qp] * _grad_u[_qp](2);
     110      555520 :   else if (jvar == _u_vel_var_number)
     111             :   {
     112      138880 :     if (_supg)
     113      138880 :       return computePGVelocityJacobian(0);
     114             :     else
     115             :       return 0;
     116             :   }
     117      416640 :   else if (jvar == _v_vel_var_number)
     118             :   {
     119      138880 :     if (_supg)
     120      138880 :       return computePGVelocityJacobian(1);
     121             :     else
     122             :       return 0;
     123             :   }
     124      277760 :   else if (jvar == _w_vel_var_number)
     125             :   {
     126           0 :     if (_supg)
     127           0 :       return computePGVelocityJacobian(2);
     128             :     else
     129             :       return 0;
     130             :   }
     131             :   else
     132             :     return 0.0;
     133             : }