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             : 
      16             : InputParameters
      17        1308 : PorousFlowFullySaturatedDarcyBase::validParams()
      18             : {
      19        1308 :   InputParameters params = Kernel::validParams();
      20        2616 :   params.addRequiredParam<RealVectorValue>("gravity",
      21             :                                            "Gravitational acceleration vector downwards (m/s^2)");
      22        2616 :   params.addParam<bool>("multiply_by_density",
      23        2616 :                         true,
      24             :                         "If true, then this Kernel is the fluid mass "
      25             :                         "flux.  If false, then this Kernel is the "
      26             :                         "fluid volume flux (which is common in "
      27             :                         "poro-mechanics)");
      28        2616 :   params.addRequiredParam<UserObjectName>(
      29             :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
      30        1308 :   params.addClassDescription("Darcy flux suitable for models involving a fully-saturated, single "
      31             :                              "phase, single component fluid.  No upwinding is used");
      32        1308 :   return params;
      33           0 : }
      34             : 
      35         708 : PorousFlowFullySaturatedDarcyBase::PorousFlowFullySaturatedDarcyBase(
      36         708 :     const InputParameters & parameters)
      37             :   : Kernel(parameters),
      38         708 :     _multiply_by_density(getParam<bool>("multiply_by_density")),
      39        1416 :     _permeability(getMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
      40         708 :     _dpermeability_dvar(
      41         708 :         getMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar")),
      42        1416 :     _dpermeability_dgradvar(getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
      43             :         "dPorousFlow_permeability_qp_dgradvar")),
      44        1416 :     _density(getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
      45        1416 :     _ddensity_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
      46             :         "dPorousFlow_fluid_phase_density_qp_dvar")),
      47        1416 :     _viscosity(getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_qp")),
      48         708 :     _dviscosity_dvar(
      49         708 :         getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_qp_dvar")),
      50        1416 :     _pp(getMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp")),
      51        1416 :     _grad_p(getMaterialProperty<std::vector<RealGradient>>("PorousFlow_grad_porepressure_qp")),
      52        1416 :     _dgrad_p_dgrad_var(getMaterialProperty<std::vector<std::vector<Real>>>(
      53             :         "dPorousFlow_grad_porepressure_qp_dgradvar")),
      54        1416 :     _dgrad_p_dvar(getMaterialProperty<std::vector<std::vector<RealGradient>>>(
      55             :         "dPorousFlow_grad_porepressure_qp_dvar")),
      56         708 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      57        1416 :     _gravity(getParam<RealVectorValue>("gravity")),
      58        1416 :     _perm_derivs(_dictator.usePermDerivs())
      59             : {
      60         708 :   if (_dictator.numPhases() != 1)
      61           0 :     mooseError("PorousFlowFullySaturatedDarcyBase should not be used for multi-phase scenarios as "
      62             :                "it does no upwinding and does not include relative-permeability effects");
      63         708 : }
      64             : 
      65             : Real
      66    20469940 : PorousFlowFullySaturatedDarcyBase::computeQpResidual()
      67             : {
      68             :   const unsigned ph = 0;
      69    20469940 :   const Real mob = mobility();
      70             :   const RealVectorValue flow =
      71    20469940 :       _permeability[_qp] * (_grad_p[_qp][ph] - _density[_qp][ph] * _gravity);
      72    20469940 :   return _grad_test[_i][_qp] * mob * flow;
      73             : }
      74             : 
      75             : Real
      76    93649504 : PorousFlowFullySaturatedDarcyBase::computeQpJacobian()
      77             : {
      78    93649504 :   return computeQpOffDiagJacobian(_var.number());
      79             : }
      80             : 
      81             : Real
      82   265900544 : PorousFlowFullySaturatedDarcyBase::computeQpOffDiagJacobian(unsigned int jvar)
      83             : {
      84   265900544 :   if (_dictator.notPorousFlowVariable(jvar))
      85             :     return 0.0;
      86             : 
      87             :   const unsigned ph = 0;
      88   265900544 :   const unsigned pvar = _dictator.porousFlowVariableNum(jvar);
      89             : 
      90   265900544 :   const Real mob = mobility();
      91   265900544 :   const Real dmob = dmobility(pvar) * _phi[_j][_qp];
      92             : 
      93             :   const RealVectorValue flow =
      94   265900544 :       _permeability[_qp] * (_grad_p[_qp][ph] - _density[_qp][ph] * _gravity);
      95             : 
      96             :   RealVectorValue dflow =
      97   265900544 :       _permeability[_qp] * (_grad_phi[_j][_qp] * _dgrad_p_dgrad_var[_qp][ph][pvar] -
      98   265900544 :                             _phi[_j][_qp] * _ddensity_dvar[_qp][ph][pvar] * _gravity);
      99   265900544 :   dflow += _permeability[_qp] * (_dgrad_p_dvar[_qp][ph][pvar] * _phi[_j][_qp]);
     100             : 
     101   265900544 :   if (_perm_derivs)
     102             :   {
     103           0 :     dflow += _dpermeability_dvar[_qp][pvar] * _phi[_j][_qp] *
     104           0 :              (_grad_p[_qp][ph] - _density[_qp][ph] * _gravity);
     105           0 :     for (const auto i : make_range(Moose::dim))
     106           0 :       dflow += _dpermeability_dgradvar[_qp][i][pvar] * _grad_phi[_j][_qp](i) *
     107           0 :                (_grad_p[_qp][ph] - _density[_qp][ph] * _gravity);
     108             :   }
     109             : 
     110   265900544 :   return _grad_test[_i][_qp] * (dmob * flow + mob * dflow);
     111             : }
     112             : 
     113             : Real
     114   381924884 : PorousFlowFullySaturatedDarcyBase::mobility() const
     115             : {
     116             :   const unsigned ph = 0;
     117   381924884 :   Real mob = 1.0 / _viscosity[_qp][ph];
     118   381924884 :   if (_multiply_by_density)
     119   232326964 :     mob *= _density[_qp][ph];
     120   381924884 :   return mob;
     121             : }
     122             : 
     123             : Real
     124   265900544 : PorousFlowFullySaturatedDarcyBase::dmobility(unsigned pvar) const
     125             : {
     126             :   const unsigned ph = 0;
     127   265900544 :   Real dmob = -_dviscosity_dvar[_qp][ph][pvar] / std::pow(_viscosity[_qp][ph], 2);
     128   265900544 :   if (_multiply_by_density)
     129   126516064 :     dmob = _density[_qp][ph] * dmob + _ddensity_dvar[_qp][ph][pvar] / _viscosity[_qp][ph];
     130   265900544 :   return dmob;
     131             : }