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             : // Navier-Stokes inclues
      11             : #include "NS.h"
      12             : #include "NSMomentumInviscidFluxWithGradP.h"
      13             : 
      14             : registerMooseObject("NavierStokesApp", NSMomentumInviscidFluxWithGradP);
      15             : 
      16             : InputParameters
      17           0 : NSMomentumInviscidFluxWithGradP::validParams()
      18             : {
      19           0 :   InputParameters params = NSKernel::validParams();
      20           0 :   params.addClassDescription(
      21             :       "This class computes the inviscid flux with pressure gradient in the momentum equation.");
      22           0 :   params.addRequiredCoupledVar(NS::pressure, "pressure");
      23           0 :   params.addRequiredParam<unsigned int>("component", "");
      24           0 :   return params;
      25           0 : }
      26             : 
      27           0 : NSMomentumInviscidFluxWithGradP::NSMomentumInviscidFluxWithGradP(const InputParameters & parameters)
      28             :   : NSKernel(parameters),
      29           0 :     _grad_p(coupledGradient(NS::pressure)),
      30           0 :     _component(getParam<unsigned int>("component")),
      31           0 :     _pressure_derivs(*this)
      32             : {
      33             :   // Store pointers to all variable gradients in a single vector.
      34             :   // This is needed for computing pressure Hessian values with a small
      35             :   // amount of code.
      36           0 :   _gradU.resize(5);
      37           0 :   _gradU[0] = &_grad_rho;
      38           0 :   _gradU[1] = &_grad_rho_u;
      39           0 :   _gradU[2] = &_grad_rho_v;
      40           0 :   _gradU[3] = &_grad_rho_w;
      41           0 :   _gradU[4] = &_grad_rho_et;
      42           0 : }
      43             : 
      44             : Real
      45           0 : NSMomentumInviscidFluxWithGradP::computeQpResidual()
      46             : {
      47             :   // For _component = k,
      48             : 
      49             :   // (rho*u) * u_k = (rho*u_k) * u <- we write it this way
      50           0 :   RealVectorValue vec(_u[_qp] * _u_vel[_qp],  // (U_k) * u_1
      51           0 :                       _u[_qp] * _v_vel[_qp],  // (U_k) * u_2
      52           0 :                       _u[_qp] * _w_vel[_qp]); // (U_k) * u_3
      53             : 
      54             :   // -((rho*u_k) * u) * grad(phi) + dp/dx_k * phi
      55           0 :   return -(vec * _grad_test[_i][_qp]) + _grad_p[_qp](_component) * _test[_i][_qp];
      56             : }
      57             : 
      58             : Real
      59           0 : NSMomentumInviscidFluxWithGradP::computeQpJacobian()
      60             : {
      61           0 :   RealVectorValue vec(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
      62             : 
      63             :   // The component'th entry of the on-diagonal Jacobian value is 2*u_i without the pressure
      64             :   // contribution.
      65           0 :   vec(_component) = 2. * vec(_component);
      66             : 
      67             :   // The Jacobian contribution due to grad(p) for the on-diagonal
      68             :   // variable, which is equal to _component+1.
      69           0 :   Real dFdp = pressureQpJacobianHelper(_component + 1);
      70             : 
      71             :   return
      72             :       // Convective terms Jacobian
      73           0 :       -(vec * _grad_test[_i][_qp]) * _phi[_j][_qp]
      74             :       // Pressure term Jacobian
      75           0 :       + dFdp * _test[_i][_qp];
      76             : }
      77             : 
      78             : Real
      79           0 : NSMomentumInviscidFluxWithGradP::computeQpOffDiagJacobian(unsigned int jvar)
      80             : {
      81           0 :   if (isNSVariable(jvar))
      82             :   {
      83             : 
      84             :     // Map jvar into the numbering expected by this->compute_pressure_jacobain_value()
      85           0 :     unsigned int var_number = mapVarNumber(jvar);
      86             : 
      87             :     // The Jacobian contribution due to differentiating the grad(p)
      88             :     // term wrt variable var_number.
      89           0 :     Real dFdp = pressureQpJacobianHelper(var_number);
      90             : 
      91           0 :     if (jvar == _rho_var_number)
      92             :     {
      93             :       // Derivative of inviscid flux convective terms wrt density:
      94             :       // x-mom: (-u_1^2   , -u_1*u_2  , -u_1*u_3 ) * grad(phi_i) * phi_j
      95             :       // y-mom: (-u_2*u_1 , -u_2^2    , -u_2*u_3 ) * grad(phi_i) * phi_j
      96             :       // z-mom: (-u_3*u_1 , -u_3*u_2  , -u_3^2   ) * grad(phi_i) * phi_j
      97             : 
      98             :       // Start with the velocity vector
      99           0 :       RealVectorValue vec(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
     100             : 
     101             :       // Scale velocity vector by -1 * vec(_component)
     102           0 :       vec *= -vec(_component);
     103             : 
     104             :       return
     105             :           // Convective terms Jacobian
     106           0 :           -(vec * _grad_test[_i][_qp]) * _phi[_j][_qp]
     107             :           // Pressure term Jacobian
     108           0 :           + dFdp * _test[_i][_qp];
     109             :     }
     110             : 
     111             :     // Handle off-diagonal derivatives wrt momentums
     112           0 :     else if (jvar == _rhou_var_number || jvar == _rhov_var_number || jvar == _rhow_var_number)
     113             :     {
     114             :       // Start with the velocity vector
     115           0 :       RealVectorValue vel(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
     116             : 
     117             :       // Map jvar into jlocal = {0,1,2}, regardless of how Moose has numbered things.
     118             :       // Can't do a case statement here since _rhou_var_number, etc. are not constants...
     119             :       unsigned int jlocal = 0;
     120             : 
     121           0 :       if (jvar == _rhov_var_number)
     122             :         jlocal = 1;
     123           0 :       else if (jvar == _rhow_var_number)
     124             :         jlocal = 2;
     125             : 
     126             :       return
     127             :           // Convective terms Jacobian
     128           0 :           -vel(_component) * _grad_test[_i][_qp](jlocal) * _phi[_j][_qp]
     129             :           // Pressure term Jacobian
     130           0 :           + dFdp * _test[_i][_qp];
     131             :     }
     132             : 
     133           0 :     else if (jvar == _rho_et_var_number)
     134             :     {
     135             :       // Pressure term Jacobian
     136           0 :       return dFdp * _test[_i][_qp];
     137             :     }
     138             :     else
     139             :       return 0.0;
     140             :   }
     141             :   else
     142             :     return 0.0;
     143             : }
     144             : 
     145             : Real
     146           0 : NSMomentumInviscidFluxWithGradP::pressureQpJacobianHelper(unsigned var_number)
     147             : {
     148             :   // Make sure our local gradient and Hessian data
     149             :   // structures are up-to-date for this quadrature point
     150             :   //  this->recalculate_gradient_and_hessian();
     151             : 
     152             :   Real hessian_sum = 0.0;
     153           0 :   for (unsigned int n = 0; n < 5; ++n)
     154           0 :     hessian_sum += _pressure_derivs.get_hess(var_number, n) * (*_gradU[n])[_qp](_component);
     155             : 
     156             :   // Hit hessian_sum with phij, then add to dp/dU_m * dphij/dx_k, finally return the result
     157           0 :   return _pressure_derivs.get_grad(var_number) * _grad_phi[_j][_qp](_component) +
     158           0 :          hessian_sum * _phi[_j][_qp];
     159             : }