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 "RichardsFlux.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "Material.h"
      14             : #include "MooseVariable.h"
      15             : 
      16             : // C++ includes
      17             : #include <iostream>
      18             : 
      19             : registerMooseObject("RichardsApp", RichardsFlux);
      20             : 
      21             : InputParameters
      22         543 : RichardsFlux::validParams()
      23             : {
      24         543 :   InputParameters params = Kernel::validParams();
      25        1086 :   params.addParam<bool>(
      26             :       "linear_shape_fcns",
      27        1086 :       true,
      28             :       "If you are using second-order Lagrange shape functions you need to set this to false.");
      29        1086 :   params.addRequiredParam<UserObjectName>(
      30             :       "richardsVarNames_UO", "The UserObject that holds the list of Richards variable names.");
      31         543 :   return params;
      32           0 : }
      33             : 
      34         259 : RichardsFlux::RichardsFlux(const InputParameters & parameters)
      35             :   : Kernel(parameters),
      36         259 :     _richards_name_UO(getUserObject<RichardsVarNames>("richardsVarNames_UO")),
      37         259 :     _pvar(_richards_name_UO.richards_var_num(_var.number())),
      38             : 
      39             :     // This kernel gets lots of things from the material
      40         518 :     _flux(getMaterialProperty<std::vector<RealVectorValue>>("flux")),
      41         518 :     _dflux_dv(getMaterialProperty<std::vector<std::vector<RealVectorValue>>>("dflux_dv")),
      42         518 :     _dflux_dgradv(getMaterialProperty<std::vector<std::vector<RealTensorValue>>>("dflux_dgradv")),
      43         259 :     _d2flux_dvdv(
      44         259 :         getMaterialProperty<std::vector<std::vector<std::vector<RealVectorValue>>>>("d2flux_dvdv")),
      45         518 :     _d2flux_dgradvdv(getMaterialProperty<std::vector<std::vector<std::vector<RealTensorValue>>>>(
      46             :         "d2flux_dgradvdv")),
      47         518 :     _d2flux_dvdgradv(getMaterialProperty<std::vector<std::vector<std::vector<RealTensorValue>>>>(
      48             :         "d2flux_dvdgradv")),
      49             : 
      50         518 :     _second_u(getParam<bool>("linear_shape_fcns")
      51         259 :                   ? _second_zero
      52           0 :                   : (_is_implicit ? _var.secondSln() : _var.secondSlnOld())),
      53         518 :     _second_phi(getParam<bool>("linear_shape_fcns") ? _second_phi_zero : secondPhi()),
      54             : 
      55         518 :     _tauvel_SUPG(getMaterialProperty<std::vector<RealVectorValue>>("tauvel_SUPG")),
      56         259 :     _dtauvel_SUPG_dgradv(
      57         259 :         getMaterialProperty<std::vector<std::vector<RealTensorValue>>>("dtauvel_SUPG_dgradv")),
      58         259 :     _dtauvel_SUPG_dv(
      59         518 :         getMaterialProperty<std::vector<std::vector<RealVectorValue>>>("dtauvel_SUPG_dv"))
      60             : {
      61         259 : }
      62             : 
      63             : Real
      64    12356644 : RichardsFlux::computeQpResidual()
      65             : {
      66    12356644 :   Real flux_part = _grad_test[_i][_qp] * _flux[_qp][_pvar];
      67             : 
      68    12356644 :   Real supg_test = _tauvel_SUPG[_qp][_pvar] * _grad_test[_i][_qp];
      69             :   Real supg_kernel = 0.0;
      70    12356644 :   if (supg_test != 0)
      71             :     // NOTE: Libmesh does not correctly calculate grad(_grad_u) correctly, so following might not be
      72             :     // correct
      73             :     // NOTE: The following is -divergence(flux)
      74             :     // NOTE: The following must be generalised if a non PPPP formalism is used
      75             :     // NOTE: The generalisation will look like
      76             :     // supg_kernel = sum_j {-(_dflux_dgradv[_qp][_pvar][j]*_second_u[_qp][j]).tr() -
      77             :     // _dflux_dv[_qp][_pvar][j]*_grad_u[_qp][j]}
      78             :     // where _grad_u[_qp][j] is the gradient of the j^th variable at the quadpoint.
      79    10815940 :     supg_kernel = -(_dflux_dgradv[_qp][_pvar][_pvar] * _second_u[_qp]).tr() -
      80    10815940 :                   _dflux_dv[_qp][_pvar][_pvar] * _grad_u[_qp];
      81             : 
      82    12356644 :   return flux_part + supg_test * supg_kernel;
      83             : }
      84             : 
      85             : Real
      86    44352280 : RichardsFlux::computeQpJac(unsigned int wrt_num)
      87             : {
      88    44352280 :   Real flux_prime = _grad_test[_i][_qp] * (_dflux_dgradv[_qp][_pvar][wrt_num] * _grad_phi[_j][_qp] +
      89    44352280 :                                            _dflux_dv[_qp][_pvar][wrt_num] * _phi[_j][_qp]);
      90             : 
      91    44352280 :   Real supg_test = _tauvel_SUPG[_qp][_pvar] * _grad_test[_i][_qp];
      92             :   Real supg_test_prime =
      93    44352280 :       _grad_phi[_j][_qp] * (_dtauvel_SUPG_dgradv[_qp][_pvar][wrt_num] * _grad_test[_i][_qp]) +
      94    44352280 :       _phi[_j][_qp] * _dtauvel_SUPG_dv[_qp][_pvar][wrt_num] * _grad_test[_i][_qp];
      95             :   Real supg_kernel = 0.0;
      96             :   Real supg_kernel_prime = 0.0;
      97             : 
      98    44352280 :   if (supg_test != 0)
      99             :   {
     100             :     // NOTE: since Libmesh does not correctly calculate grad(_grad_u) correctly, so following might
     101             :     // not be correct
     102    40497688 :     supg_kernel = -(_dflux_dgradv[_qp][_pvar][_pvar] * _second_u[_qp]).tr() -
     103    40497688 :                   _dflux_dv[_qp][_pvar][_pvar] * _grad_u[_qp];
     104             : 
     105             :     // NOTE: just like supg_kernel, this must be generalised for non-PPPP formulations
     106    40497688 :     supg_kernel_prime =
     107    40497688 :         -(_d2flux_dvdv[_qp][_pvar][_pvar][wrt_num] * _phi[_j][_qp] * _grad_u[_qp] +
     108    40497688 :           _phi[_j][_qp] * (_d2flux_dgradvdv[_qp][_pvar][_pvar][wrt_num] * _second_u[_qp]).tr() +
     109    40497688 :           (_d2flux_dvdgradv[_qp][_pvar][_pvar][wrt_num] * _grad_u[_qp]) * _grad_phi[_j][_qp]);
     110    40497688 :     if (wrt_num == _pvar)
     111    36446792 :       supg_kernel_prime -= _dflux_dv[_qp][_pvar][_pvar] * _grad_phi[_j][_qp];
     112    40497688 :     supg_kernel_prime -= (_dflux_dgradv[_qp][_pvar][_pvar] * _second_phi[_j][_qp]).tr();
     113             :   }
     114             : 
     115    44352280 :   return flux_prime + supg_test_prime * supg_kernel + supg_test * supg_kernel_prime;
     116             : }
     117             : 
     118             : Real
     119    40172792 : RichardsFlux::computeQpJacobian()
     120             : {
     121    40172792 :   return computeQpJac(_pvar);
     122             : }
     123             : 
     124             : Real
     125     4179488 : RichardsFlux::computeQpOffDiagJacobian(unsigned int jvar)
     126             : {
     127     4179488 :   if (_richards_name_UO.not_richards_var(jvar))
     128             :     return 0.0;
     129     4179488 :   unsigned int dvar = _richards_name_UO.richards_var_num(jvar);
     130     4179488 :   return computeQpJac(dvar);
     131             : }