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        2234 : INSMass::validParams()
      17             : {
      18        2234 :   InputParameters params = INSBase::validParams();
      19             : 
      20        2234 :   params.addClassDescription("This class computes the mass equation residual and Jacobian "
      21             :                              "contributions for the incompressible Navier-Stokes momentum "
      22             :                              "equation.");
      23        4468 :   params.addParam<bool>(
      24        4468 :       "pspg", false, "Whether to perform PSPG stabilization of the mass equation");
      25        4468 :   params.addParam<FunctionName>("x_vel_forcing_func", 0, "The x-velocity mms forcing function.");
      26        4468 :   params.addParam<FunctionName>("y_vel_forcing_func", 0, "The y-velocity mms forcing function.");
      27        4468 :   params.addParam<FunctionName>("z_vel_forcing_func", 0, "The z-velocity mms forcing function.");
      28        2234 :   return params;
      29           0 : }
      30             : 
      31        1205 : INSMass::INSMass(const InputParameters & parameters)
      32             :   : INSBase(parameters),
      33        1205 :     _pspg(getParam<bool>("pspg")),
      34        1205 :     _x_ffn(getFunction("x_vel_forcing_func")),
      35        1205 :     _y_ffn(getFunction("y_vel_forcing_func")),
      36        2410 :     _z_ffn(getFunction("z_vel_forcing_func"))
      37             : 
      38             : {
      39        1205 : }
      40             : 
      41             : Real
      42    48509042 : 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    48509042 :   Real r = -(_grad_u_vel[_qp](0) + _grad_v_vel[_qp](1) + _grad_w_vel[_qp](2)) * _test[_i][_qp];
      48             : 
      49    48509042 :   if (_pspg)
      50     9807602 :     r += computeQpPGResidual();
      51             : 
      52    48509042 :   return r;
      53             : }
      54             : 
      55             : Real
      56     9807602 : INSMass::computeQpPGResidual()
      57             : {
      58             :   RealVectorValue viscous_term =
      59     9807602 :       _laplace ? strongViscousTermLaplace() : strongViscousTermTraction();
      60             :   RealVectorValue transient_term =
      61     9807602 :       _transient_term ? timeDerivativeTerm() : RealVectorValue(0, 0, 0);
      62     9807602 :   RealVectorValue convective_term = _convective_term ? convectiveTerm() : RealVectorValue(0, 0, 0);
      63     9807602 :   Real r = -1. / _rho[_qp] * tau() * _grad_test[_i][_qp] *
      64     9807602 :            (strongPressureTerm() + gravityTerm() + viscous_term + convective_term + transient_term -
      65     9807602 :             RealVectorValue(_x_ffn.value(_t, _q_point[_qp]),
      66     9807602 :                             _y_ffn.value(_t, _q_point[_qp]),
      67     9807602 :                             _z_ffn.value(_t, _q_point[_qp])));
      68             : 
      69     9807602 :   return r;
      70             : }
      71             : 
      72             : Real
      73    85434894 : INSMass::computeQpJacobian()
      74             : {
      75             :   // Derivative wrt to p is zero
      76             :   Real r = 0;
      77             : 
      78             :   // Unless we are doing PSPG stabilization
      79    85434894 :   if (_pspg)
      80    22912362 :     r += computeQpPGJacobian();
      81             : 
      82    85434894 :   return r;
      83             : }
      84             : 
      85             : Real
      86    22912362 : INSMass::computeQpPGJacobian()
      87             : {
      88    22912362 :   return -1. / _rho[_qp] * tau() * _grad_test[_i][_qp] * dStrongPressureDPressure();
      89             : }
      90             : 
      91             : Real
      92   343822842 : INSMass::computeQpOffDiagJacobian(const unsigned int jvar)
      93             : {
      94   343822842 :   if (jvar == _u_vel_var_number)
      95             :   {
      96   163222254 :     Real jac = -_grad_phi[_j][_qp](0) * _test[_i][_qp];
      97   163222254 :     if (_pspg)
      98    23658822 :       jac += computeQpPGOffDiagJacobian(0);
      99   163222254 :     return jac;
     100             :   }
     101             : 
     102   180600588 :   else if (jvar == _v_vel_var_number)
     103             :   {
     104   163222254 :     Real jac = -_grad_phi[_j][_qp](1) * _test[_i][_qp];
     105   163222254 :     if (_pspg)
     106    23658822 :       jac += computeQpPGOffDiagJacobian(1);
     107   163222254 :     return jac;
     108             :   }
     109             : 
     110    17378334 :   else if (jvar == _w_vel_var_number)
     111             :   {
     112      118098 :     Real jac = -_grad_phi[_j][_qp](2) * _test[_i][_qp];
     113      118098 :     if (_pspg)
     114      118098 :       jac += computeQpPGOffDiagJacobian(2);
     115      118098 :     return jac;
     116             :   }
     117             : 
     118             :   else
     119             :     return 0.0;
     120             : }
     121             : 
     122             : Real
     123    47435742 : INSMass::computeQpPGOffDiagJacobian(const unsigned int comp)
     124             : {
     125    47435742 :   RealVectorValue convective_term = _convective_term ? convectiveTerm() : RealVectorValue(0, 0, 0);
     126             :   RealVectorValue d_convective_term_d_u_comp =
     127    47435742 :       _convective_term ? dConvecDUComp(comp) : RealVectorValue(0, 0, 0);
     128             :   RealVectorValue viscous_term =
     129    47435742 :       _laplace ? strongViscousTermLaplace() : strongViscousTermTraction();
     130             :   RealVectorValue d_viscous_term_d_u_comp =
     131    47435742 :       _laplace ? dStrongViscDUCompLaplace(comp) : dStrongViscDUCompTraction(comp);
     132             :   RealVectorValue transient_term =
     133    47435742 :       _transient_term ? timeDerivativeTerm() : RealVectorValue(0, 0, 0);
     134             :   RealVectorValue d_transient_term_d_u_comp =
     135    47435742 :       _transient_term ? dTimeDerivativeDUComp(comp) : RealVectorValue(0, 0, 0);
     136             : 
     137    47435742 :   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    47435742 :          1. / _rho[_qp] * dTauDUComp(comp) * _grad_test[_i][_qp] *
     140    47435742 :              (convective_term + viscous_term + transient_term + strongPressureTerm() +
     141    47435742 :               gravityTerm() -
     142    47435742 :               RealVectorValue(_x_ffn.value(_t, _q_point[_qp]),
     143    47435742 :                               _y_ffn.value(_t, _q_point[_qp]),
     144    47435742 :                               _z_ffn.value(_t, _q_point[_qp])));
     145             : }