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