PorousFlowFullySaturatedDarcyBase.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "PorousFlowFullySaturatedDarcyBase.h"

#include "MooseVariable.h"

registerMooseObject("PorousFlowApp", PorousFlowFullySaturatedDarcyBase);
registerMooseObject("PorousFlowApp", ADPorousFlowFullySaturatedDarcyBase);

template <bool is_ad>
InputParameters
PorousFlowFullySaturatedDarcyBaseTempl<is_ad>::validParams()
{
  InputParameters params = GenericKernel<is_ad>::validParams();
  params.addRequiredParam<RealVectorValue>("gravity",
                                           "Gravitational acceleration vector downwards (m/s^2)");
  params.addParam<bool>("multiply_by_density",
                        true,
                        "If true, then this Kernel is the fluid mass flux.  If false, then this "
                        "Kernel is the fluid volume flux (which is common in poro-mechanics)");
  params.addRequiredParam<UserObjectName>(
      "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
  params.addClassDescription(
      "Darcy flux suitable for models involving a fully-saturated, single phase, single component "
      "fluid.  No upwinding is used.  Templated on is_ad: the AD version requires no hand-coded "
      "Jacobian.");
  return params;
}

template <bool is_ad>
PorousFlowFullySaturatedDarcyBaseTempl<is_ad>::PorousFlowFullySaturatedDarcyBaseTempl(
    const InputParameters & parameters)
  : GenericKernel<is_ad>(parameters),
    _multiply_by_density(this->template getParam<bool>("multiply_by_density")),
    _permeability(this->template getGenericMaterialProperty<RealTensorValue, is_ad>(
        "PorousFlow_permeability_qp")),
    _dpermeability_dvar(is_ad ? nullptr
                              : &this->template getMaterialProperty<std::vector<RealTensorValue>>(
                                    "dPorousFlow_permeability_qp_dvar")),
    _dpermeability_dgradvar(
        is_ad ? nullptr
              : &this->template getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
                    "dPorousFlow_permeability_qp_dgradvar")),
    _density(this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
        "PorousFlow_fluid_phase_density_qp")),
    _ddensity_dvar(is_ad ? nullptr
                         : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
                               "dPorousFlow_fluid_phase_density_qp_dvar")),
    _viscosity(this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
        "PorousFlow_viscosity_qp")),
    _dviscosity_dvar(is_ad ? nullptr
                           : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
                                 "dPorousFlow_viscosity_qp_dvar")),
    _pp(this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
        "PorousFlow_porepressure_qp")),
    _grad_p(this->template getGenericMaterialProperty<std::vector<RealGradient>, is_ad>(
        "PorousFlow_grad_porepressure_qp")),
    _dgrad_p_dgrad_var(is_ad ? nullptr
                             : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
                                   "dPorousFlow_grad_porepressure_qp_dgradvar")),
    _dgrad_p_dvar(is_ad
                      ? nullptr
                      : &this->template getMaterialProperty<std::vector<std::vector<RealGradient>>>(
                            "dPorousFlow_grad_porepressure_qp_dvar")),
    _dictator(this->template getUserObject<PorousFlowDictator>("PorousFlowDictator")),
    _gravity(this->template getParam<RealVectorValue>("gravity")),
    _perm_derivs(_dictator.usePermDerivs())
{
  if (_dictator.numPhases() != 1)
    mooseError(
        "PorousFlowFullySaturatedDarcyBase should not be used for multi-phase scenarios as it "
        "does no upwinding and does not include relative-permeability effects");
}

template <bool is_ad>
GenericReal<is_ad>
PorousFlowFullySaturatedDarcyBaseTempl<is_ad>::computeQpResidual()
{
  const unsigned ph = 0;
  const auto mob = mobility();
  const auto flow = _permeability[_qp] * (_grad_p[_qp][ph] - _density[_qp][ph] * _gravity);
  return _grad_test[_i][_qp] * mob * flow;
}

template <bool is_ad>
Real
PorousFlowFullySaturatedDarcyBaseTempl<is_ad>::computeQpJacobian()
{
  if constexpr (!is_ad)
    return computeQpOffDiagJacobian(_var.number());
  return 0.0;
}

template <bool is_ad>
Real
PorousFlowFullySaturatedDarcyBaseTempl<is_ad>::computeQpOffDiagJacobian(unsigned int jvar)
{
  if constexpr (!is_ad)
  {
    if (_dictator.notPorousFlowVariable(jvar))
      return 0.0;

    const unsigned ph = 0;
    const unsigned pvar = _dictator.porousFlowVariableNum(jvar);

    const Real mob = mobility();
    const Real dmob = dmobility(pvar) * _phi[_j][_qp];

    const RealVectorValue flow =
        _permeability[_qp] * (_grad_p[_qp][ph] - _density[_qp][ph] * _gravity);

    RealVectorValue dflow =
        _permeability[_qp] * (_grad_phi[_j][_qp] * (*_dgrad_p_dgrad_var)[_qp][ph][pvar] -
                              _phi[_j][_qp] * (*_ddensity_dvar)[_qp][ph][pvar] * _gravity);
    dflow += _permeability[_qp] * ((*_dgrad_p_dvar)[_qp][ph][pvar] * _phi[_j][_qp]);

    if (_perm_derivs)
    {
      dflow += (*_dpermeability_dvar)[_qp][pvar] * _phi[_j][_qp] *
               (_grad_p[_qp][ph] - _density[_qp][ph] * _gravity);
      for (const auto i : make_range(Moose::dim))
        dflow += (*_dpermeability_dgradvar)[_qp][i][pvar] * _grad_phi[_j][_qp](i) *
                 (_grad_p[_qp][ph] - _density[_qp][ph] * _gravity);
    }

    return _grad_test[_i][_qp] * (dmob * flow + mob * dflow);
  }
  else
    libmesh_ignore(jvar);
  return 0.0;
}

template <bool is_ad>
GenericReal<is_ad>
PorousFlowFullySaturatedDarcyBaseTempl<is_ad>::mobility() const
{
  const unsigned ph = 0;
  GenericReal<is_ad> mob = 1.0 / _viscosity[_qp][ph];
  if (_multiply_by_density)
    mob *= _density[_qp][ph];
  return mob;
}

template <bool is_ad>
Real
PorousFlowFullySaturatedDarcyBaseTempl<is_ad>::dmobility(unsigned int pvar) const
{
  if constexpr (!is_ad)
  {
    const unsigned ph = 0;
    Real dmob = -(*_dviscosity_dvar)[_qp][ph][pvar] / std::pow(_viscosity[_qp][ph], 2);
    if (_multiply_by_density)
      dmob = _density[_qp][ph] * dmob + (*_ddensity_dvar)[_qp][ph][pvar] / _viscosity[_qp][ph];
    return dmob;
  }
  else
    libmesh_ignore(pvar);
  return 0.0;
}

template class PorousFlowFullySaturatedDarcyBaseTempl<false>;
template class PorousFlowFullySaturatedDarcyBaseTempl<true>;