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 "INSMomentumBase.h"
      11             : #include "Function.h"
      12             : 
      13             : InputParameters
      14        2094 : INSMomentumBase::validParams()
      15             : {
      16        2094 :   InputParameters params = INSBase::validParams();
      17             : 
      18        4188 :   params.addRequiredParam<unsigned>("component", "The velocity component that this is applied to.");
      19        4188 :   params.addParam<bool>(
      20        4188 :       "integrate_p_by_parts", true, "Whether to integrate the pressure term by parts.");
      21        4188 :   params.addParam<bool>(
      22        4188 :       "supg", false, "Whether to perform SUPG stabilization of the momentum residuals");
      23        4188 :   params.addParam<FunctionName>("forcing_func", 0, "The mms forcing function.");
      24        2094 :   return params;
      25           0 : }
      26             : 
      27        1142 : INSMomentumBase::INSMomentumBase(const InputParameters & parameters)
      28             :   : INSBase(parameters),
      29        1142 :     _component(getParam<unsigned>("component")),
      30        2284 :     _integrate_p_by_parts(getParam<bool>("integrate_p_by_parts")),
      31        2284 :     _supg(getParam<bool>("supg")),
      32        2284 :     _ffn(getFunction("forcing_func"))
      33             : {
      34        1142 :   if (_supg && !_convective_term)
      35           0 :     mooseError("It doesn't make sense to conduct SUPG stabilization without a convective term.");
      36        1142 : }
      37             : 
      38             : Real
      39   121398118 : INSMomentumBase::computeQpResidual()
      40             : {
      41             :   Real r = 0;
      42             : 
      43             :   // viscous term
      44   121398118 :   r += computeQpResidualViscousPart();
      45             : 
      46             :   // pressure term
      47   121398118 :   if (_integrate_p_by_parts)
      48    59530728 :     r += _grad_test[_i][_qp](_component) * weakPressureTerm();
      49             :   else
      50    61867390 :     r += _test[_i][_qp] * strongPressureTerm()(_component);
      51             : 
      52             :   // body force term
      53   121398118 :   r += _test[_i][_qp] * (gravityTerm()(_component) - _ffn.value(_t, _q_point[_qp]));
      54             : 
      55             :   // convective term
      56   121398118 :   if (_convective_term)
      57   121388398 :     r += _test[_i][_qp] * convectiveTerm()(_component);
      58             : 
      59   121398118 :   if (_supg)
      60    19909368 :     r += computeQpPGResidual();
      61             : 
      62   121398118 :   return r;
      63             : }
      64             : 
      65             : Real
      66    19909368 : INSMomentumBase::computeQpPGResidual()
      67             : {
      68    19909368 :   const auto U = relativeVelocity();
      69             : 
      70    19909368 :   const auto convective_term = _convective_term ? convectiveTerm() : RealVectorValue(0, 0, 0);
      71    19909368 :   const auto viscous_term = _laplace ? strongViscousTermLaplace() : strongViscousTermTraction();
      72    19909368 :   const auto transient_term = _transient_term ? timeDerivativeTerm() : RealVectorValue(0, 0, 0);
      73             : 
      74    19909368 :   return tau() * U * _grad_test[_i][_qp] *
      75    19909368 :          ((convective_term + viscous_term + transient_term + strongPressureTerm() +
      76    39818736 :            gravityTerm())(_component)-_ffn.value(_t, _q_point[_qp]));
      77             : 
      78             :   // For GLS as opposed to SUPG stabilization, one would need to modify the test function functional
      79             :   // space to include second derivatives of the Galerkin test functions corresponding to the viscous
      80             :   // term. This would look like:
      81             :   // Real lap_test =
      82             :   //     _second_test[_i][_qp](0, 0) + _second_test[_i][_qp](1, 1) + _second_test[_i][_qp](2, 2);
      83             : 
      84             :   // Real pg_viscous_r = -_mu[_qp] * lap_test * tau() *
      85             :   //                     (convective_term + viscous_term + strongPressureTerm()(_component) +
      86             :   //                      gravityTerm())(_component);
      87             : }
      88             : 
      89             : Real
      90   470850624 : INSMomentumBase::computeQpJacobian()
      91             : {
      92             :   Real jac = 0;
      93             : 
      94             :   // viscous term
      95   470850624 :   jac += computeQpJacobianViscousPart();
      96             : 
      97             :   // convective term
      98   470850624 :   if (_convective_term)
      99   470806884 :     jac += _test[_i][_qp] * dConvecDUComp(_component)(_component);
     100             : 
     101   470850624 :   if (_supg)
     102    82420938 :     jac += computeQpPGJacobian(_component);
     103             : 
     104   470850624 :   return jac;
     105             : }
     106             : 
     107             : Real
     108   164959974 : INSMomentumBase::computeQpPGJacobian(unsigned comp)
     109             : {
     110   164959974 :   const auto U = relativeVelocity();
     111             :   RealVectorValue d_U_d_U_comp(0, 0, 0);
     112   164959974 :   d_U_d_U_comp(comp) = _phi[_j][_qp];
     113             : 
     114   164959974 :   const Real convective_term = _convective_term ? convectiveTerm()(_component) : 0;
     115   164959974 :   const Real d_convective_term_d_u_comp = _convective_term ? dConvecDUComp(comp)(_component) : 0;
     116             :   const Real viscous_term =
     117   164959974 :       _laplace ? strongViscousTermLaplace()(_component) : strongViscousTermTraction()(_component);
     118   164959974 :   const Real d_viscous_term_d_u_comp = _laplace ? dStrongViscDUCompLaplace(comp)(_component)
     119    18702747 :                                                 : dStrongViscDUCompTraction(comp)(_component);
     120   164959974 :   const Real transient_term = _transient_term ? timeDerivativeTerm()(_component) : 0;
     121             :   const Real d_transient_term_d_u_comp =
     122   164959974 :       _transient_term ? dTimeDerivativeDUComp(comp)(_component) : 0;
     123             : 
     124   164959974 :   return dTauDUComp(comp) * U * _grad_test[_i][_qp] *
     125   164959974 :              (convective_term + viscous_term + strongPressureTerm()(_component) +
     126   164959974 :               gravityTerm()(_component) + transient_term - _ffn.value(_t, _q_point[_qp])) +
     127   164959974 :          tau() * d_U_d_U_comp * _grad_test[_i][_qp] *
     128   164959974 :              (convective_term + viscous_term + strongPressureTerm()(_component) +
     129   164959974 :               gravityTerm()(_component) + transient_term - _ffn.value(_t, _q_point[_qp])) +
     130   164959974 :          tau() * U * _grad_test[_i][_qp] *
     131   164959974 :              (d_convective_term_d_u_comp + d_viscous_term_d_u_comp + d_transient_term_d_u_comp);
     132             : }
     133             : 
     134             : Real
     135   755846442 : INSMomentumBase::computeQpOffDiagJacobian(unsigned jvar)
     136             : {
     137             :   Real jac = 0;
     138   755846442 :   if (jvar == _u_vel_var_number)
     139             :   {
     140   235442115 :     Real convective_term = _convective_term ? _test[_i][_qp] * dConvecDUComp(0)(_component) : 0.;
     141   235442115 :     Real viscous_term = computeQpOffDiagJacobianViscousPart(jvar);
     142             : 
     143   235442115 :     jac += convective_term + viscous_term;
     144             : 
     145   235442115 :     if (_supg)
     146    41230152 :       jac += computeQpPGJacobian(0);
     147             : 
     148   235442115 :     return jac;
     149             :   }
     150   520404327 :   else if (jvar == _v_vel_var_number)
     151             :   {
     152   235442115 :     Real convective_term = _convective_term ? _test[_i][_qp] * dConvecDUComp(1)(_component) : 0.;
     153   235442115 :     Real viscous_term = computeQpOffDiagJacobianViscousPart(jvar);
     154             : 
     155   235442115 :     jac += convective_term + viscous_term;
     156             : 
     157   235442115 :     if (_supg)
     158    41230152 :       jac += computeQpPGJacobian(1);
     159             : 
     160   235442115 :     return jac;
     161             :   }
     162   284962212 :   else if (jvar == _w_vel_var_number)
     163             :   {
     164       78732 :     Real convective_term = _convective_term ? _test[_i][_qp] * dConvecDUComp(2)(_component) : 0.;
     165       78732 :     Real viscous_term = computeQpOffDiagJacobianViscousPart(jvar);
     166             : 
     167       78732 :     jac += convective_term + viscous_term;
     168             : 
     169       78732 :     if (_supg)
     170       78732 :       jac += computeQpPGJacobian(2);
     171             : 
     172       78732 :     return jac;
     173             :   }
     174             : 
     175   284883480 :   else if (jvar == _p_var_number)
     176             :   {
     177   228291174 :     if (_integrate_p_by_parts)
     178   175929725 :       jac += _grad_test[_i][_qp](_component) * dWeakPressureDPressure();
     179             :     else
     180    52361449 :       jac += _test[_i][_qp] * dStrongPressureDPressure()(_component);
     181             : 
     182   228291174 :     if (_supg)
     183             :     {
     184    52201458 :       const auto U = relativeVelocity();
     185    52201458 :       jac += tau() * U * _grad_test[_i][_qp] * dStrongPressureDPressure()(_component);
     186             :     }
     187             : 
     188   228291174 :     return jac;
     189             :   }
     190             : 
     191             :   else
     192             :     return 0;
     193             : }