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 "PorousFlowFullySaturatedDarcyBase.h"
      11             : 
      12             : #include "MooseVariable.h"
      13             : 
      14             : registerMooseObject("PorousFlowApp", PorousFlowFullySaturatedDarcyBase);
      15             : registerMooseObject("PorousFlowApp", ADPorousFlowFullySaturatedDarcyBase);
      16             : 
      17             : template <bool is_ad>
      18             : InputParameters
      19        2687 : PorousFlowFullySaturatedDarcyBaseTempl<is_ad>::validParams()
      20             : {
      21             :   InputParameters params = GenericKernel<is_ad>::validParams();
      22        2758 :   params.addRequiredParam<RealVectorValue>("gravity",
      23             :                                            "Gravitational acceleration vector downwards (m/s^2)");
      24        2758 :   params.addParam<bool>("multiply_by_density",
      25        2758 :                         true,
      26             :                         "If true, then this Kernel is the fluid mass flux.  If false, then this "
      27             :                         "Kernel is the fluid volume flux (which is common in poro-mechanics)");
      28        2758 :   params.addRequiredParam<UserObjectName>(
      29             :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
      30        1379 :   params.addClassDescription(
      31             :       "Darcy flux suitable for models involving a fully-saturated, single phase, single component "
      32             :       "fluid.  No upwinding is used.  Templated on is_ad: the AD version requires no hand-coded "
      33             :       "Jacobian.");
      34        1379 :   return params;
      35           0 : }
      36             : 
      37             : template <bool is_ad>
      38         746 : PorousFlowFullySaturatedDarcyBaseTempl<is_ad>::PorousFlowFullySaturatedDarcyBaseTempl(
      39             :     const InputParameters & parameters)
      40             :   : GenericKernel<is_ad>(parameters),
      41         746 :     _multiply_by_density(this->template getParam<bool>("multiply_by_density")),
      42        1492 :     _permeability(this->template getGenericMaterialProperty<RealTensorValue, is_ad>(
      43             :         "PorousFlow_permeability_qp")),
      44         746 :     _dpermeability_dvar(is_ad ? nullptr
      45         708 :                               : &this->template getMaterialProperty<std::vector<RealTensorValue>>(
      46             :                                     "dPorousFlow_permeability_qp_dvar")),
      47         746 :     _dpermeability_dgradvar(
      48             :         is_ad ? nullptr
      49         708 :               : &this->template getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
      50             :                     "dPorousFlow_permeability_qp_dgradvar")),
      51        1492 :     _density(this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
      52             :         "PorousFlow_fluid_phase_density_qp")),
      53         746 :     _ddensity_dvar(is_ad ? nullptr
      54         708 :                          : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
      55             :                                "dPorousFlow_fluid_phase_density_qp_dvar")),
      56        1492 :     _viscosity(this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
      57             :         "PorousFlow_viscosity_qp")),
      58         746 :     _dviscosity_dvar(is_ad ? nullptr
      59         708 :                            : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
      60             :                                  "dPorousFlow_viscosity_qp_dvar")),
      61        1492 :     _pp(this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
      62             :         "PorousFlow_porepressure_qp")),
      63        1492 :     _grad_p(this->template getGenericMaterialProperty<std::vector<RealGradient>, is_ad>(
      64             :         "PorousFlow_grad_porepressure_qp")),
      65         746 :     _dgrad_p_dgrad_var(is_ad ? nullptr
      66         708 :                              : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
      67             :                                    "dPorousFlow_grad_porepressure_qp_dgradvar")),
      68         746 :     _dgrad_p_dvar(is_ad
      69             :                       ? nullptr
      70         708 :                       : &this->template getMaterialProperty<std::vector<std::vector<RealGradient>>>(
      71             :                             "dPorousFlow_grad_porepressure_qp_dvar")),
      72         746 :     _dictator(this->template getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      73        1492 :     _gravity(this->template getParam<RealVectorValue>("gravity")),
      74        1492 :     _perm_derivs(_dictator.usePermDerivs())
      75             : {
      76         746 :   if (_dictator.numPhases() != 1)
      77           0 :     mooseError(
      78             :         "PorousFlowFullySaturatedDarcyBase should not be used for multi-phase scenarios as it "
      79             :         "does no upwinding and does not include relative-permeability effects");
      80         746 : }
      81             : 
      82             : template <bool is_ad>
      83             : GenericReal<is_ad>
      84    20506308 : PorousFlowFullySaturatedDarcyBaseTempl<is_ad>::computeQpResidual()
      85             : {
      86             :   const unsigned ph = 0;
      87    20506308 :   const auto mob = mobility();
      88    20506308 :   const auto flow = _permeability[_qp] * (_grad_p[_qp][ph] - _density[_qp][ph] * _gravity);
      89    20506308 :   return _grad_test[_i][_qp] * mob * flow;
      90             : }
      91             : 
      92             : template <bool is_ad>
      93             : Real
      94    93649504 : PorousFlowFullySaturatedDarcyBaseTempl<is_ad>::computeQpJacobian()
      95             : {
      96             :   if constexpr (!is_ad)
      97    93649504 :     return computeQpOffDiagJacobian(_var.number());
      98           0 :   return 0.0;
      99             : }
     100             : 
     101             : template <bool is_ad>
     102             : Real
     103   265900544 : PorousFlowFullySaturatedDarcyBaseTempl<is_ad>::computeQpOffDiagJacobian(unsigned int jvar)
     104             : {
     105             :   if constexpr (!is_ad)
     106             :   {
     107   265900544 :     if (_dictator.notPorousFlowVariable(jvar))
     108             :       return 0.0;
     109             : 
     110             :     const unsigned ph = 0;
     111   265900544 :     const unsigned pvar = _dictator.porousFlowVariableNum(jvar);
     112             : 
     113   265900544 :     const Real mob = mobility();
     114   265900544 :     const Real dmob = dmobility(pvar) * _phi[_j][_qp];
     115             : 
     116             :     const RealVectorValue flow =
     117   265900544 :         _permeability[_qp] * (_grad_p[_qp][ph] - _density[_qp][ph] * _gravity);
     118             : 
     119             :     RealVectorValue dflow =
     120   265900544 :         _permeability[_qp] * (_grad_phi[_j][_qp] * (*_dgrad_p_dgrad_var)[_qp][ph][pvar] -
     121   265900544 :                               _phi[_j][_qp] * (*_ddensity_dvar)[_qp][ph][pvar] * _gravity);
     122   265900544 :     dflow += _permeability[_qp] * ((*_dgrad_p_dvar)[_qp][ph][pvar] * _phi[_j][_qp]);
     123             : 
     124   265900544 :     if (_perm_derivs)
     125             :     {
     126           0 :       dflow += (*_dpermeability_dvar)[_qp][pvar] * _phi[_j][_qp] *
     127           0 :                (_grad_p[_qp][ph] - _density[_qp][ph] * _gravity);
     128           0 :       for (const auto i : make_range(Moose::dim))
     129           0 :         dflow += (*_dpermeability_dgradvar)[_qp][i][pvar] * _grad_phi[_j][_qp](i) *
     130           0 :                  (_grad_p[_qp][ph] - _density[_qp][ph] * _gravity);
     131             :     }
     132             : 
     133   265900544 :     return _grad_test[_i][_qp] * (dmob * flow + mob * dflow);
     134             :   }
     135             :   else
     136             :     libmesh_ignore(jvar);
     137           0 :   return 0.0;
     138             : }
     139             : 
     140             : template <bool is_ad>
     141             : GenericReal<is_ad>
     142   381961252 : PorousFlowFullySaturatedDarcyBaseTempl<is_ad>::mobility() const
     143             : {
     144             :   const unsigned ph = 0;
     145   381961252 :   GenericReal<is_ad> mob = 1.0 / _viscosity[_qp][ph];
     146   381961252 :   if (_multiply_by_density)
     147   232363332 :     mob *= _density[_qp][ph];
     148   381961252 :   return mob;
     149             : }
     150             : 
     151             : template <bool is_ad>
     152             : Real
     153   265900544 : PorousFlowFullySaturatedDarcyBaseTempl<is_ad>::dmobility(unsigned int pvar) const
     154             : {
     155             :   if constexpr (!is_ad)
     156             :   {
     157             :     const unsigned ph = 0;
     158   265900544 :     Real dmob = -(*_dviscosity_dvar)[_qp][ph][pvar] / std::pow(_viscosity[_qp][ph], 2);
     159   265900544 :     if (_multiply_by_density)
     160   126516064 :       dmob = _density[_qp][ph] * dmob + (*_ddensity_dvar)[_qp][ph][pvar] / _viscosity[_qp][ph];
     161             :     return dmob;
     162             :   }
     163             :   else
     164             :     libmesh_ignore(pvar);
     165           0 :   return 0.0;
     166             : }
     167             : 
     168             : template class PorousFlowFullySaturatedDarcyBaseTempl<false>;
     169             : template class PorousFlowFullySaturatedDarcyBaseTempl<true>;