PorousFlowPorosity.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 "PorousFlowPorosity.h"
 
registerMooseObject("PorousFlowApp", PorousFlowPorosity);
registerMooseObject("PorousFlowApp", ADPorousFlowPorosity);
 
template <bool is_ad>
InputParameters
PorousFlowPorosityTempl<is_ad>::validParams()
{
  InputParameters params = PorousFlowPorosityExponentialBaseTempl<is_ad>::validParams();
  params.addParam<bool>(
      "mechanical", false, "If true, porosity will be a function of total volumetric strain");
  params.addParam<bool>(
      "fluid", false, "If true, porosity will be a function of effective porepressure");
  params.addParam<bool>("thermal", false, "If true, porosity will be a function of temperature");
  params.addParam<bool>("chemical", false, "If true, porosity will be a function of precipitate");
  params.addRequiredCoupledVar("porosity_zero",
                               "The porosity at zero volumetric strain and "
                               "reference temperature and reference effective "
                               "porepressure and reference chemistry.  This must be a real number "
                               "or a constant monomial variable (not a linear lagrange or other "
                               "type of variable)");
  params.addParam<Real>("thermal_expansion_coeff",
                        "Volumetric thermal expansion coefficient of the drained porous solid "
                        "skeleton (only used if thermal=true)");
  params.addRangeCheckedParam<Real>(
      "biot_coefficient", 1, "biot_coefficient>=0 & biot_coefficient<=1", "Biot coefficient");
  params.addParam<Real>("biot_coefficient_prime",
                        "Biot coefficient that appears in the term (biot_coefficient_prime - 1) * "
                        "(P - reference_porepressure) / solid_bulk.  If not provided, this "
                        "defaults to the standard biot_coefficient");
  params.addParam<MooseFunctorName>(
      "solid_bulk", "Bulk modulus of the drained porous solid skeleton (only used if fluid=true)");
  params.addCoupledVar(
      "reference_temperature", 0.0, "Reference temperature (only used if thermal=true)");
  params.addCoupledVar(
      "reference_porepressure", 0.0, "Reference porepressure (only used if fluid=true)");
  params.addCoupledVar("reference_chemistry",
                       "Reference values of the solid mineral concentrations "
                       "(m^3(precipitate)/m^3(porous material)), entered as "
                       "a vector (one value per mineral).  (Only used if chemical=true)");
  params.addCoupledVar(
      "initial_mineral_concentrations",
      "Initial mineral concentrations (m^3(precipitate)/m^3(porous material)), entered as "
      "a vector (one value per mineral).  (Only used if chemical=true)");
  params.addParam<std::vector<Real>>("chemical_weights",
                                     "When chemical=true, porosity is a linear combination of the "
                                     "solid mineral concentrations multiplied by these weights.  "
                                     "Default=1 for all minerals.");
  params.addClassDescription("This Material calculates the porosity PorousFlow simulations");
  return params;
}
 
template <bool is_ad>
PorousFlowPorosityTempl<is_ad>::PorousFlowPorosityTempl(const InputParameters & parameters)
  : PorousFlowPorosityExponentialBaseTempl<is_ad>(parameters),
 
    _mechanical(this->template getParam<bool>("mechanical")),
    _fluid(this->template getParam<bool>("fluid")),
    _thermal(this->template getParam<bool>("thermal")),
    _chemical(this->template getParam<bool>("chemical")),
    _phi0(coupledValue("porosity_zero")),
    _biot(this->template getParam<Real>("biot_coefficient")),
    _exp_coeff(this->isParamValid("thermal_expansion_coeff")
                   ? this->template getParam<Real>("thermal_expansion_coeff")
                   : 0.0),
    _solid_bulk(this->isParamValid("solid_bulk") ? &(this->template getFunctor<Real>("solid_bulk"))
                                                 : nullptr),
    _coeff(this->isParamValid("biot_coefficient_prime")
               ? (this->template getParam<Real>("biot_coefficient_prime") - 1.0)
               : (_biot - 1.0)),
 
    _t_reference(nodalOrQpValue("reference_temperature")),
    _p_reference(nodalOrQpValue("reference_porepressure")),
    _num_c_ref(coupledComponents("reference_chemistry")),
    _c_reference(_num_c_ref),
    _num_initial_c(coupledComponents("initial_mineral_concentrations")),
    _initial_c(_num_initial_c),
    _c_weights(this->isParamValid("chemical_weights")
                   ? this->template getParam<std::vector<Real>>("chemical_weights")
                   : std::vector<Real>(_num_c_ref, 1.0)),
 
    _porosity_old(
        _chemical ? (_nodal_material
                         ? &this->template getMaterialPropertyOld<Real>("PorousFlow_porosity_nodal")
                         : &this->template getMaterialPropertyOld<Real>("PorousFlow_porosity_qp"))
                  : nullptr),
    _vol_strain_qp(_mechanical ? &this->template getGenericMaterialProperty<Real, is_ad>(
                                     "PorousFlow_total_volumetric_strain_qp")
                               : nullptr),
    _dvol_strain_qp_dvar((_mechanical && !is_ad)
                             ? &this->template getMaterialProperty<std::vector<RealGradient>>(
                                   "dPorousFlow_total_volumetric_strain_qp_dvar")
                             : nullptr),
 
    _pf(_fluid ? (_nodal_material ? &this->template getGenericMaterialProperty<Real, is_ad>(
                                        "PorousFlow_effective_fluid_pressure_nodal")
                                  : &this->template getGenericMaterialProperty<Real, is_ad>(
                                        "PorousFlow_effective_fluid_pressure_qp"))
               : nullptr),
    _dpf_dvar((_fluid && !is_ad)
                  ? (_nodal_material ? &this->template getMaterialProperty<std::vector<Real>>(
                                           "dPorousFlow_effective_fluid_pressure_nodal_dvar")
                                     : &this->template getMaterialProperty<std::vector<Real>>(
                                           "dPorousFlow_effective_fluid_pressure_qp_dvar"))
                  : nullptr),
 
    _temperature(_thermal
                     ? (_nodal_material ? &this->template getGenericMaterialProperty<Real, is_ad>(
                                              "PorousFlow_temperature_nodal")
                                        : &this->template getGenericMaterialProperty<Real, is_ad>(
                                              "PorousFlow_temperature_qp"))
                     : nullptr),
    _dtemperature_dvar((_thermal && !is_ad)
                           ? (_nodal_material
                                  ? &this->template getMaterialProperty<std::vector<Real>>(
                                        "dPorousFlow_temperature_nodal_dvar")
                                  : &this->template getMaterialProperty<std::vector<Real>>(
                                        "dPorousFlow_temperature_qp_dvar"))
                           : nullptr),
 
    _mineral_conc_old(_chemical ? (_nodal_material
                                       ? &this->template getMaterialPropertyOld<std::vector<Real>>(
                                             "PorousFlow_mineral_concentration_nodal")
                                       : &this->template getMaterialPropertyOld<std::vector<Real>>(
                                             "PorousFlow_mineral_concentration_qp"))
                                : nullptr),
    _reaction_rate(_chemical
                       ? (_nodal_material ? &this->template getMaterialProperty<std::vector<Real>>(
                                                "PorousFlow_mineral_reaction_rate_nodal")
                                          : &this->template getMaterialProperty<std::vector<Real>>(
                                                "PorousFlow_mineral_reaction_rate_qp"))
                       : nullptr),
    _dreaction_rate_dvar(
        _chemical ? (_nodal_material
                         ? &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
                               "dPorousFlow_mineral_reaction_rate_nodal_dvar")
                         : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
                               "dPorousFlow_mineral_reaction_rate_qp_dvar"))
                  : nullptr),
    _aq_ph(_dictator.aqueousPhaseNumber()),
    _saturation(_chemical
                    ? (_nodal_material ? &this->template getMaterialProperty<std::vector<Real>>(
                                             "PorousFlow_saturation_nodal")
                                       : &this->template getMaterialProperty<std::vector<Real>>(
                                             "PorousFlow_saturation_qp"))
                    : nullptr),
    _dsaturation_dvar(
        _chemical ? (_nodal_material
                         ? &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
                               "dPorousFlow_saturation_nodal_dvar")
                         : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
                               "dPorousFlow_saturation_qp_dvar"))
                  : nullptr)
{
  // PorousFlowAqueousPreDisChemistry and PorousFlowAqueousPreDisMineral are not templated on is_ad,
  // so the mineral concentration and reaction rate consumed below exist only as non-AD properties
  if constexpr (is_ad)
    if (_chemical)
      this->mooseError(
          "ADPorousFlowPorosity does not support chemical=true: the aqueous "
          "precipitation/dissolution materials are not AD.  Use the non-AD PorousFlowPorosity for "
          "chemistry-coupled porosity");
 
  if (_thermal && !this->isParamValid("thermal_expansion_coeff"))
    this->mooseError(
        "PorousFlowPorosity: When thermal=true you must provide a thermal_expansion_coeff");
  if (_fluid && !_solid_bulk)
    this->mooseError("PorousFlowPorosity: When fluid=true you must provide a solid_bulk");
  if (_chemical && _num_c_ref != _dictator.numAqueousKinetic())
    this->mooseError(
        "PorousFlowPorosity: When chemical=true you must provide the reference_chemistry "
        "values.  The Dictator proclaims there should be ",
        _dictator.numAqueousKinetic(),
        " of these");
  if (_chemical && _num_initial_c != _dictator.numAqueousKinetic())
    this->mooseError("PorousFlowPorosity: When chemical=true you must provide the "
                     "initial_mineral_concentrations.  "
                     "The Dictator proclaims there should be ",
                     _dictator.numAqueousKinetic(),
                     " of these");
  if (_chemical && _c_weights.size() != _dictator.numAqueousKinetic())
    this->mooseError(
        "PorousFlowPorosity: When chemical=true you must provde the correct number of "
        "chemical_weights (which the Dictator knows is ",
        _dictator.numAqueousKinetic(),
        ") or do not provide any chemical_weights and use the default value of 1 for each mineral");
 
  for (unsigned i = 0; i < _num_c_ref; ++i)
  {
    _c_reference[i] = &nodalOrQpValue("reference_chemistry", i);
    _initial_c[i] = &nodalOrQpValue("initial_mineral_concentrations", i);
  }
}
 
template <bool is_ad>
GenericReal<is_ad>
PorousFlowPorosityTempl<is_ad>::atNegInfinityQp() const
{
  /*
   *
   * Note the use of the OLD value of porosity here.
   * This strategy, which breaks the cyclic dependency between porosity
   * and mineral concentration, is used in
   * Kernel: PorousFlowPreDis
   * Material: PorousFlowPorosity
   * Material: PorousFlowAqueousPreDisChemistry
   * Material: PorousFlowAqueousPreDisMineral
   *
   */
  GenericReal<is_ad> result = _biot;
  if (_chemical)
  {
    if (_t_step == 0 && !_app.isRestarting())
      for (unsigned i = 0; i < _num_c_ref; ++i)
        result -= _c_weights[i] * (*_initial_c[i])[_qp];
    else
      for (unsigned i = 0; i < _num_c_ref; ++i)
        result -= _c_weights[i] * ((*_mineral_conc_old)[_qp][i] + _dt * (*_porosity_old)[_qp] *
                                                                      (*_saturation)[_qp][_aq_ph] *
                                                                      (*_reaction_rate)[_qp][i]);
  }
  return result;
}
 
template <bool is_ad>
Real
PorousFlowPorosityTempl<is_ad>::datNegInfinityQp(unsigned pvar) const
{
  Real result = 0.0;
  if (_chemical && (_t_step >= 1 || _app.isRestarting()))
    for (unsigned i = 0; i < _num_c_ref; ++i)
      result -= _c_weights[i] * _dt * (*_porosity_old)[_qp] *
                ((*_saturation)[_qp][_aq_ph] * (*_dreaction_rate_dvar)[_qp][i][pvar] +
                 (*_dsaturation_dvar)[_qp][_aq_ph][pvar] * (*_reaction_rate)[_qp][i]);
  return result;
}
 
template <bool is_ad>
GenericReal<is_ad>
PorousFlowPorosityTempl<is_ad>::atZeroQp() const
{
  // note the [0] below: _phi0 is a constant monomial and we use [0] regardless of _nodal_material
  GenericReal<is_ad> result = _phi0[0];
  if (_chemical)
  {
    if (_t_step == 0 && !_app.isRestarting())
      for (unsigned i = 0; i < _num_c_ref; ++i)
        result -= _c_weights[i] * ((*_initial_c[i])[_qp] - (*_c_reference[i])[_qp]);
    else
      for (unsigned i = 0; i < _num_c_ref; ++i)
        result -= _c_weights[i] * ((*_mineral_conc_old)[_qp][i] +
                                   _dt * (*_porosity_old)[_qp] * (*_saturation)[_qp][_aq_ph] *
                                       (*_reaction_rate)[_qp][i] -
                                   (*_c_reference[i])[_qp]);
  }
  return result;
}
 
template <bool is_ad>
Real
PorousFlowPorosityTempl<is_ad>::datZeroQp(unsigned pvar) const
{
  Real result = 0.0;
  if (_chemical && (_t_step >= 1 || _app.isRestarting()))
    for (unsigned i = 0; i < _num_c_ref; ++i)
      result -= _c_weights[i] * _dt * (*_porosity_old)[_qp] *
                ((*_saturation)[_qp][_aq_ph] * (*_dreaction_rate_dvar)[_qp][i][pvar] +
                 (*_dsaturation_dvar)[_qp][_aq_ph][pvar] * (*_reaction_rate)[_qp][i]);
  return result;
}
 
template <bool is_ad>
GenericReal<is_ad>
PorousFlowPorosityTempl<is_ad>::decayQp() const
{
  GenericReal<is_ad> result = 0.0;
 
  if (_thermal)
    result += _exp_coeff * ((*_temperature)[_qp] - _t_reference[_qp]);
 
  if (_fluid)
  {
    Real solid_bulk;
    // Using Qp 0 can leverage the functor caching
    // TODO: Find a way to effectively use subdomain-constant-ness
    unsigned int qp_used = (_constant_option == Material::ConstantTypeEnum::NONE) ? _qp : 0;
    if (_nodal_material)
    {
      const std::set<SubdomainID> subdomain_set = {_current_elem->subdomain_id()};
      const Moose::NodeArg space_arg = {_current_elem->node_ptr(qp_used), &subdomain_set};
      solid_bulk = (*_solid_bulk)(space_arg, Moose::currentState());
    }
    else if (_bnd)
    {
      const Moose::ElemSideQpArg space_arg = {
          _current_elem, _current_side, qp_used, _qrule, _q_point[qp_used]};
      solid_bulk = (*_solid_bulk)(space_arg, Moose::currentState());
    }
    else
    {
      const Moose::ElemQpArg space_arg = {_current_elem, qp_used, _qrule, _q_point[qp_used]};
      solid_bulk = (*_solid_bulk)(space_arg, Moose::currentState());
    }
    if (solid_bulk <= 0)
      this->mooseError("PorousFlowPorosity: solid_bulk must be larger than Zero");
    result += _coeff / solid_bulk * ((*_pf)[_qp] - _p_reference[_qp]);
  }
 
  if (_mechanical)
  {
    // Note that in the following _strain[_qp] is evaluated at q quadpoint
    // So _porosity_nodal[_qp], which should be the nodal value of porosity
    // actually uses the strain at a quadpoint.  This
    // is OK for LINEAR elements, as strain is constant over the element anyway.
    const unsigned qp_to_use =
        (_nodal_material && (_bnd || _strain_at_nearest_qp) ? nearestQP(_qp) : _qp);
    result += -(*_vol_strain_qp)[qp_to_use];
  }
 
  return result;
}
 
template <bool is_ad>
Real
PorousFlowPorosityTempl<is_ad>::ddecayQp_dvar(unsigned pvar) const
{
  Real result = 0.0;
 
  if (_thermal)
    result += _exp_coeff * (*_dtemperature_dvar)[_qp][pvar];
 
  if (_fluid)
  {
    Real solid_bulk;
    // Using Qp 0 can leverage the functor caching
    // TODO: Find a way to effectively use subdomain-constant-ness
    unsigned int qp_used = (_constant_option == Material::ConstantTypeEnum::NONE) ? _qp : 0;
    if (_nodal_material)
    {
      const std::set<SubdomainID> subdomain_set = {_current_elem->subdomain_id()};
      const Moose::NodeArg space_arg = {_current_elem->node_ptr(qp_used), &subdomain_set};
      solid_bulk = (*_solid_bulk)(space_arg, Moose::currentState());
    }
    else if (_bnd)
    {
      const Moose::ElemSideQpArg space_arg = {
          _current_elem, _current_side, qp_used, _qrule, _q_point[qp_used]};
      solid_bulk = (*_solid_bulk)(space_arg, Moose::currentState());
    }
    else
    {
      const Moose::ElemQpArg space_arg = {_current_elem, qp_used, _qrule, _q_point[qp_used]};
      solid_bulk = (*_solid_bulk)(space_arg, Moose::currentState());
    }
    if (solid_bulk <= 0)
      this->mooseError("PorousFlowPorosity: solid_bulk must be larger than Zero.");
    result += _coeff / solid_bulk * (*_dpf_dvar)[_qp][pvar];
  }
 
  return result;
}
 
template <bool is_ad>
RealGradient
PorousFlowPorosityTempl<is_ad>::ddecayQp_dgradvar(unsigned pvar) const
{
  RealGradient result(0.0, 0.0, 0.0);
  if (_mechanical)
  {
    const unsigned qp_to_use =
        (_nodal_material && (_bnd || _strain_at_nearest_qp) ? nearestQP(_qp) : _qp);
    result += -(*_dvol_strain_qp_dvar)[qp_to_use][pvar];
  }
  return result;
}
 
template class PorousFlowPorosityTempl<false>;
template class PorousFlowPorosityTempl<true>;