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 "INSMass.h"
      11             : #include "Function.h"
      12             : 
      13             : registerMooseObject("NavierStokesApp", INSMass);
      14             : 
      15             : InputParameters
      16        1078 : INSMass::validParams()
      17             : {
      18        1078 :   InputParameters params = INSBase::validParams();
      19             : 
      20        1078 :   params.addClassDescription("This class computes the mass equation residual and Jacobian "
      21             :                              "contributions for the incompressible Navier-Stokes momentum "
      22             :                              "equation.");
      23        2156 :   params.addParam<bool>(
      24        2156 :       "pspg", false, "Whether to perform PSPG stabilization of the mass equation");
      25        2156 :   params.addParam<FunctionName>("x_vel_forcing_func", 0, "The x-velocity mms forcing function.");
      26        2156 :   params.addParam<FunctionName>("y_vel_forcing_func", 0, "The y-velocity mms forcing function.");
      27        2156 :   params.addParam<FunctionName>("z_vel_forcing_func", 0, "The z-velocity mms forcing function.");
      28        1078 :   return params;
      29           0 : }
      30             : 
      31         569 : INSMass::INSMass(const InputParameters & parameters)
      32             :   : INSBase(parameters),
      33         569 :     _pspg(getParam<bool>("pspg")),
      34         569 :     _x_ffn(getFunction("x_vel_forcing_func")),
      35         569 :     _y_ffn(getFunction("y_vel_forcing_func")),
      36        1138 :     _z_ffn(getFunction("z_vel_forcing_func"))
      37             : 
      38             : {
      39         569 : }
      40             : 
      41             : Real
      42    30654791 : INSMass::computeQpResidual()
      43             : {
      44             :   // (div u) * q
      45             :   // Note: we (arbitrarily) multiply this term by -1 so that it matches the -p(div v)
      46             :   // term in the momentum equation.  Not sure if that is really important?
      47    30654791 :   Real r = -(_grad_u_vel[_qp](0) + _grad_v_vel[_qp](1) + _grad_w_vel[_qp](2)) * _test[_i][_qp];
      48             : 
      49    30654791 :   if (_pspg)
      50     6597656 :     r += computeQpPGResidual();
      51             : 
      52    30654791 :   return r;
      53             : }
      54             : 
      55             : Real
      56     6597656 : INSMass::computeQpPGResidual()
      57             : {
      58             :   RealVectorValue viscous_term =
      59     6597656 :       _laplace ? strongViscousTermLaplace() : strongViscousTermTraction();
      60             :   RealVectorValue transient_term =
      61     6597656 :       _transient_term ? timeDerivativeTerm() : RealVectorValue(0, 0, 0);
      62     6597656 :   RealVectorValue convective_term = _convective_term ? convectiveTerm() : RealVectorValue(0, 0, 0);
      63     6597656 :   Real r = -1. / _rho[_qp] * tau() * _grad_test[_i][_qp] *
      64     6597656 :            (strongPressureTerm() + gravityTerm() + viscous_term + convective_term + transient_term -
      65     6597656 :             RealVectorValue(_x_ffn.value(_t, _q_point[_qp]),
      66     6597656 :                             _y_ffn.value(_t, _q_point[_qp]),
      67     6597656 :                             _z_ffn.value(_t, _q_point[_qp])));
      68             : 
      69     6597656 :   return r;
      70             : }
      71             : 
      72             : Real
      73    59843349 : INSMass::computeQpJacobian()
      74             : {
      75             :   // Derivative wrt to p is zero
      76             :   Real r = 0;
      77             : 
      78             :   // Unless we are doing PSPG stabilization
      79    59843349 :   if (_pspg)
      80    16222230 :     r += computeQpPGJacobian();
      81             : 
      82    59843349 :   return r;
      83             : }
      84             : 
      85             : Real
      86    16222230 : INSMass::computeQpPGJacobian()
      87             : {
      88    16222230 :   return -1. / _rho[_qp] * tau() * _grad_test[_i][_qp] * dStrongPressureDPressure();
      89             : }
      90             : 
      91             : Real
      92   240931242 : INSMass::computeQpOffDiagJacobian(const unsigned int jvar)
      93             : {
      94   240931242 :   if (jvar == _u_vel_var_number)
      95             :   {
      96   114125904 :     Real jac = -_grad_phi[_j][_qp](0) * _test[_i][_qp];
      97   114125904 :     if (_pspg)
      98    16719030 :       jac += computeQpPGOffDiagJacobian(0);
      99   114125904 :     return jac;
     100             :   }
     101             : 
     102   126805338 :   else if (jvar == _v_vel_var_number)
     103             :   {
     104   114125904 :     Real jac = -_grad_phi[_j][_qp](1) * _test[_i][_qp];
     105   114125904 :     if (_pspg)
     106    16719030 :       jac += computeQpPGOffDiagJacobian(1);
     107   114125904 :     return jac;
     108             :   }
     109             : 
     110    12679434 :   else if (jvar == _w_vel_var_number)
     111             :   {
     112       39366 :     Real jac = -_grad_phi[_j][_qp](2) * _test[_i][_qp];
     113       39366 :     if (_pspg)
     114       39366 :       jac += computeQpPGOffDiagJacobian(2);
     115       39366 :     return jac;
     116             :   }
     117             : 
     118             :   else
     119             :     return 0.0;
     120             : }
     121             : 
     122             : Real
     123    33477426 : INSMass::computeQpPGOffDiagJacobian(const unsigned int comp)
     124             : {
     125    33477426 :   RealVectorValue convective_term = _convective_term ? convectiveTerm() : RealVectorValue(0, 0, 0);
     126             :   RealVectorValue d_convective_term_d_u_comp =
     127    33477426 :       _convective_term ? dConvecDUComp(comp) : RealVectorValue(0, 0, 0);
     128             :   RealVectorValue viscous_term =
     129    33477426 :       _laplace ? strongViscousTermLaplace() : strongViscousTermTraction();
     130             :   RealVectorValue d_viscous_term_d_u_comp =
     131    33477426 :       _laplace ? dStrongViscDUCompLaplace(comp) : dStrongViscDUCompTraction(comp);
     132             :   RealVectorValue transient_term =
     133    33477426 :       _transient_term ? timeDerivativeTerm() : RealVectorValue(0, 0, 0);
     134             :   RealVectorValue d_transient_term_d_u_comp =
     135    33477426 :       _transient_term ? dTimeDerivativeDUComp(comp) : RealVectorValue(0, 0, 0);
     136             : 
     137    33477426 :   return -1. / _rho[_qp] * tau() * _grad_test[_i][_qp] *
     138             :              (d_convective_term_d_u_comp + d_viscous_term_d_u_comp + d_transient_term_d_u_comp) -
     139    33477426 :          1. / _rho[_qp] * dTauDUComp(comp) * _grad_test[_i][_qp] *
     140    33477426 :              (convective_term + viscous_term + transient_term + strongPressureTerm() +
     141    33477426 :               gravityTerm() -
     142    33477426 :               RealVectorValue(_x_ffn.value(_t, _q_point[_qp]),
     143    33477426 :                               _y_ffn.value(_t, _q_point[_qp]),
     144    33477426 :                               _z_ffn.value(_t, _q_point[_qp])));
     145             : }