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        4358 : INSMomentumBase::validParams()
      15             : {
      16        4358 :   InputParameters params = INSBase::validParams();
      17             : 
      18        8716 :   params.addRequiredParam<unsigned>("component", "The velocity component that this is applied to.");
      19        8716 :   params.addParam<bool>(
      20        8716 :       "integrate_p_by_parts", true, "Whether to integrate the pressure term by parts.");
      21        8716 :   params.addParam<bool>(
      22        8716 :       "supg", false, "Whether to perform SUPG stabilization of the momentum residuals");
      23        8716 :   params.addParam<FunctionName>("forcing_func", 0, "The mms forcing function.");
      24        4358 :   return params;
      25           0 : }
      26             : 
      27        2422 : INSMomentumBase::INSMomentumBase(const InputParameters & parameters)
      28             :   : INSBase(parameters),
      29        2422 :     _component(getParam<unsigned>("component")),
      30        4844 :     _integrate_p_by_parts(getParam<bool>("integrate_p_by_parts")),
      31        4844 :     _supg(getParam<bool>("supg")),
      32        4844 :     _ffn(getFunction("forcing_func"))
      33             : {
      34        2422 :   if (_supg && !_convective_term)
      35           0 :     mooseError("It doesn't make sense to conduct SUPG stabilization without a convective term.");
      36        2422 : }
      37             : 
      38             : Real
      39   193933702 : INSMomentumBase::computeQpResidual()
      40             : {
      41             :   Real r = 0;
      42             : 
      43             :   // viscous term
      44   193933702 :   r += computeQpResidualViscousPart();
      45             : 
      46             :   // pressure term
      47   193933702 :   if (_integrate_p_by_parts)
      48    85195500 :     r += _grad_test[_i][_qp](_component) * weakPressureTerm();
      49             :   else
      50   108738202 :     r += _test[_i][_qp] * strongPressureTerm()(_component);
      51             : 
      52             :   // body force term
      53   193933702 :   r += _test[_i][_qp] * (gravityTerm()(_component) - _ffn.value(_t, _q_point[_qp]));
      54             : 
      55             :   // convective term
      56   193933702 :   if (_convective_term)
      57   193919122 :     r += _test[_i][_qp] * convectiveTerm()(_component);
      58             : 
      59   193933702 :   if (_supg)
      60    29744364 :     r += computeQpPGResidual();
      61             : 
      62   193933702 :   return r;
      63             : }
      64             : 
      65             : Real
      66    29744364 : INSMomentumBase::computeQpPGResidual()
      67             : {
      68    29744364 :   const auto U = relativeVelocity();
      69             : 
      70    29744364 :   const auto convective_term = _convective_term ? convectiveTerm() : RealVectorValue(0, 0, 0);
      71    29744364 :   const auto viscous_term = _laplace ? strongViscousTermLaplace() : strongViscousTermTraction();
      72    29744364 :   const auto transient_term = _transient_term ? timeDerivativeTerm() : RealVectorValue(0, 0, 0);
      73             : 
      74    29744364 :   return tau() * U * _grad_test[_i][_qp] *
      75    29744364 :          ((convective_term + viscous_term + transient_term + strongPressureTerm() +
      76    59488728 :            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   674025696 : INSMomentumBase::computeQpJacobian()
      91             : {
      92             :   Real jac = 0;
      93             : 
      94             :   // viscous term
      95   674025696 :   jac += computeQpJacobianViscousPart();
      96             : 
      97             :   // convective term
      98   674025696 :   if (_convective_term)
      99   673960086 :     jac += _test[_i][_qp] * dConvecDUComp(_component)(_component);
     100             : 
     101   674025696 :   if (_supg)
     102   119543832 :     jac += computeQpPGJacobian(_component);
     103             : 
     104   674025696 :   return jac;
     105             : }
     106             : 
     107             : Real
     108   239441958 : INSMomentumBase::computeQpPGJacobian(unsigned comp)
     109             : {
     110   239441958 :   const auto U = relativeVelocity();
     111             :   RealVectorValue d_U_d_U_comp(0, 0, 0);
     112   239441958 :   d_U_d_U_comp(comp) = _phi[_j][_qp];
     113             : 
     114   239441958 :   const Real convective_term = _convective_term ? convectiveTerm()(_component) : 0;
     115   239441958 :   const Real d_convective_term_d_u_comp = _convective_term ? dConvecDUComp(comp)(_component) : 0;
     116             :   const Real viscous_term =
     117   239441958 :       _laplace ? strongViscousTermLaplace()(_component) : strongViscousTermTraction()(_component);
     118   239441958 :   const Real d_viscous_term_d_u_comp = _laplace ? dStrongViscDUCompLaplace(comp)(_component)
     119    28367361 :                                                 : dStrongViscDUCompTraction(comp)(_component);
     120   239441958 :   const Real transient_term = _transient_term ? timeDerivativeTerm()(_component) : 0;
     121             :   const Real d_transient_term_d_u_comp =
     122   239441958 :       _transient_term ? dTimeDerivativeDUComp(comp)(_component) : 0;
     123             : 
     124   239441958 :   return dTauDUComp(comp) * U * _grad_test[_i][_qp] *
     125   239441958 :              (convective_term + viscous_term + strongPressureTerm()(_component) +
     126   239441958 :               gravityTerm()(_component) + transient_term - _ffn.value(_t, _q_point[_qp])) +
     127   239441958 :          tau() * d_U_d_U_comp * _grad_test[_i][_qp] *
     128   239441958 :              (convective_term + viscous_term + strongPressureTerm()(_component) +
     129   239441958 :               gravityTerm()(_component) + transient_term - _ffn.value(_t, _q_point[_qp])) +
     130   239441958 :          tau() * U * _grad_test[_i][_qp] *
     131   239441958 :              (d_convective_term_d_u_comp + d_viscous_term_d_u_comp + d_transient_term_d_u_comp);
     132             : }
     133             : 
     134             : Real
     135  1078175898 : INSMomentumBase::computeQpOffDiagJacobian(unsigned jvar)
     136             : {
     137             :   Real jac = 0;
     138  1078175898 :   if (jvar == _u_vel_var_number)
     139             :   {
     140   337069017 :     Real convective_term = _convective_term ? _test[_i][_qp] * dConvecDUComp(0)(_component) : 0.;
     141   337069017 :     Real viscous_term = computeQpOffDiagJacobianViscousPart(jvar);
     142             : 
     143   337069017 :     jac += convective_term + viscous_term;
     144             : 
     145   337069017 :     if (_supg)
     146    59830965 :       jac += computeQpPGJacobian(0);
     147             : 
     148   337069017 :     return jac;
     149             :   }
     150   741106881 :   else if (jvar == _v_vel_var_number)
     151             :   {
     152   337069017 :     Real convective_term = _convective_term ? _test[_i][_qp] * dConvecDUComp(1)(_component) : 0.;
     153   337069017 :     Real viscous_term = computeQpOffDiagJacobianViscousPart(jvar);
     154             : 
     155   337069017 :     jac += convective_term + viscous_term;
     156             : 
     157   337069017 :     if (_supg)
     158    59830965 :       jac += computeQpPGJacobian(1);
     159             : 
     160   337069017 :     return jac;
     161             :   }
     162   404037864 :   else if (jvar == _w_vel_var_number)
     163             :   {
     164      236196 :     Real convective_term = _convective_term ? _test[_i][_qp] * dConvecDUComp(2)(_component) : 0.;
     165      236196 :     Real viscous_term = computeQpOffDiagJacobianViscousPart(jvar);
     166             : 
     167      236196 :     jac += convective_term + viscous_term;
     168             : 
     169      236196 :     if (_supg)
     170      236196 :       jac += computeQpPGJacobian(2);
     171             : 
     172      236196 :     return jac;
     173             :   }
     174             : 
     175   403801668 :   else if (jvar == _p_var_number)
     176             :   {
     177   326562606 :     if (_integrate_p_by_parts)
     178   248922631 :       jac += _grad_test[_i][_qp](_component) * dWeakPressureDPressure();
     179             :     else
     180    77639975 :       jac += _test[_i][_qp] * dStrongPressureDPressure()(_component);
     181             : 
     182   326562606 :     if (_supg)
     183             :     {
     184    74941182 :       const auto U = relativeVelocity();
     185    74941182 :       jac += tau() * U * _grad_test[_i][_qp] * dStrongPressureDPressure()(_component);
     186             :     }
     187             : 
     188   326562606 :     return jac;
     189             :   }
     190             : 
     191             :   else
     192             :     return 0;
     193             : }