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 "ConvectedMeshPSPG.h"
      11             : 
      12             : registerMooseObject("FsiApp", ConvectedMeshPSPG);
      13             : 
      14             : InputParameters
      15          41 : ConvectedMeshPSPG::validParams()
      16             : {
      17          41 :   InputParameters params = INSBase::validParams();
      18          41 :   params.addClassDescription(
      19             :       "Corrects the convective derivative for situations in which the fluid mesh is dynamic.");
      20          82 :   params.addRequiredCoupledVar("disp_x", "The x displacement");
      21          82 :   params.addCoupledVar("disp_y", "The y displacement");
      22          82 :   params.addCoupledVar("disp_z", "The z displacement");
      23          82 :   params.addParam<MaterialPropertyName>("rho_name", "rho", "The name of the density");
      24          41 :   return params;
      25           0 : }
      26             : 
      27          22 : ConvectedMeshPSPG::ConvectedMeshPSPG(const InputParameters & parameters)
      28             :   : INSBase(parameters),
      29          22 :     _disp_x_dot(coupledDot("disp_x")),
      30          22 :     _d_disp_x_dot(coupledDotDu("disp_x")),
      31          22 :     _disp_x_id(coupled("disp_x")),
      32          44 :     _disp_y_dot(isCoupled("disp_y") ? coupledDot("disp_y") : _zero),
      33          44 :     _d_disp_y_dot(isCoupled("disp_y") ? coupledDotDu("disp_y") : _zero),
      34          22 :     _disp_y_id(coupled("disp_y")),
      35          22 :     _disp_z_dot(isCoupled("disp_z") ? coupledDot("disp_z") : _zero),
      36          22 :     _d_disp_z_dot(isCoupled("disp_z") ? coupledDotDu("disp_z") : _zero),
      37          22 :     _disp_z_id(coupled("disp_z")),
      38          66 :     _rho(getMaterialProperty<Real>("rho_name"))
      39             : {
      40          22 : }
      41             : 
      42             : RealVectorValue
      43    16854880 : ConvectedMeshPSPG::strongResidual()
      44             : {
      45             :   const auto minus_rho_ddisp_dt =
      46    16854880 :       -_rho[_qp] * RealVectorValue(_disp_x_dot[_qp], _disp_y_dot[_qp], _disp_z_dot[_qp]);
      47    16854880 :   return RealVectorValue(minus_rho_ddisp_dt * _grad_u_vel[_qp],
      48    16854880 :                          minus_rho_ddisp_dt * _grad_v_vel[_qp],
      49    16854880 :                          minus_rho_ddisp_dt * _grad_w_vel[_qp]);
      50             : }
      51             : 
      52             : RealVectorValue
      53      412160 : ConvectedMeshPSPG::dStrongResidualDDisp(const unsigned short component)
      54             : {
      55     1236480 :   const auto & ddisp_dot = [&]() -> const VariableValue &
      56             :   {
      57      412160 :     switch (component)
      58             :     {
      59      206080 :       case 0:
      60      206080 :         return _d_disp_x_dot;
      61      206080 :       case 1:
      62      206080 :         return _d_disp_y_dot;
      63           0 :       case 2:
      64           0 :         return _d_disp_z_dot;
      65           0 :       default:
      66           0 :         mooseError("Invalid component");
      67             :     }
      68      412160 :   }();
      69             : 
      70             :   // Only non-zero component will be from 'component'
      71             :   RealVectorValue ddisp_dt;
      72      412160 :   ddisp_dt(component) = _phi[_j][_qp] * ddisp_dot[_qp];
      73             : 
      74      412160 :   const auto minus_rho_ddisp_dt = -_rho[_qp] * ddisp_dt;
      75      412160 :   return RealVectorValue(minus_rho_ddisp_dt * _grad_u_vel[_qp],
      76      412160 :                          minus_rho_ddisp_dt * _grad_v_vel[_qp],
      77      412160 :                          minus_rho_ddisp_dt * _grad_w_vel[_qp]);
      78             : }
      79             : 
      80             : RealVectorValue
      81      412160 : ConvectedMeshPSPG::dStrongResidualDVel(const unsigned short component)
      82             : {
      83             :   const auto minus_rho_ddisp_dt =
      84      412160 :       -_rho[_qp] * RealVectorValue(_disp_x_dot[_qp], _disp_y_dot[_qp], _disp_z_dot[_qp]);
      85             : 
      86             :   // Only non-zero component will be from 'component'
      87             :   RealVectorValue ret;
      88      412160 :   ret(component) = minus_rho_ddisp_dt * _grad_phi[_j][_qp];
      89      412160 :   return ret;
      90             : }
      91             : 
      92             : Real
      93    16442720 : ConvectedMeshPSPG::computeQpResidual()
      94             : {
      95    16442720 :   return -tau() / _rho[_qp] * _grad_test[_i][_qp] * strongResidual();
      96             : }
      97             : 
      98             : Real
      99      206080 : ConvectedMeshPSPG::computeQpJacobian()
     100             : {
     101             :   // No derivative with respect to pressure
     102      206080 :   return 0;
     103             : }
     104             : 
     105             : Real
     106      824320 : ConvectedMeshPSPG::computeQpOffDiagJacobian(unsigned int jvar)
     107             : {
     108             :   mooseAssert(jvar != _var.number(), "Making sure I understand how old hand-coded Jacobians work.");
     109             : 
     110      824320 :   if (jvar == _disp_x_id)
     111      206080 :     return -tau() / _rho[_qp] * _grad_test[_i][_qp] * dStrongResidualDDisp(0);
     112      618240 :   else if (jvar == _disp_y_id)
     113      206080 :     return -tau() / _rho[_qp] * _grad_test[_i][_qp] * dStrongResidualDDisp(1);
     114      412160 :   else if (jvar == _disp_z_id)
     115           0 :     return -tau() / _rho[_qp] * _grad_test[_i][_qp] * dStrongResidualDDisp(2);
     116      412160 :   else if (jvar == _u_vel_var_number)
     117      412160 :     return -dTauDUComp(0) / _rho[_qp] * _grad_test[_i][_qp] * strongResidual() -
     118      206080 :            tau() / _rho[_qp] * _grad_test[_i][_qp] * dStrongResidualDVel(0);
     119      206080 :   else if (jvar == _v_vel_var_number)
     120      412160 :     return -dTauDUComp(1) / _rho[_qp] * _grad_test[_i][_qp] * strongResidual() -
     121      206080 :            tau() / _rho[_qp] * _grad_test[_i][_qp] * dStrongResidualDVel(1);
     122           0 :   else if (jvar == _w_vel_var_number)
     123           0 :     return -dTauDUComp(2) / _rho[_qp] * _grad_test[_i][_qp] * strongResidual() -
     124           0 :            tau() / _rho[_qp] * _grad_test[_i][_qp] * dStrongResidualDVel(2);
     125             :   else
     126             :     return 0.0;
     127             : }