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 "PorousFlowFullySaturatedMassTimeDerivative.h"
      11             : 
      12             : #include "MooseVariable.h"
      13             : 
      14             : registerMooseObject("PorousFlowApp", PorousFlowFullySaturatedMassTimeDerivative);
      15             : registerMooseObject("PorousFlowApp", ADPorousFlowFullySaturatedMassTimeDerivative);
      16             : 
      17             : template <bool is_ad>
      18             : InputParameters
      19         916 : PorousFlowFullySaturatedMassTimeDerivativeTempl<is_ad>::validParams()
      20             : {
      21             :   InputParameters params = GenericKernel<is_ad>::validParams();
      22             :   // Mark residual contributions as time terms (same as TimeKernel / ADTimeKernel)
      23         948 :   params.set<MultiMooseEnum>("vector_tags") = "time";
      24         948 :   params.set<MultiMooseEnum>("matrix_tags") = "system time";
      25         948 :   MooseEnum coupling_type("Hydro ThermoHydro HydroMechanical ThermoHydroMechanical", "Hydro");
      26         948 :   params.addParam<MooseEnum>(
      27             :       "coupling_type",
      28             :       coupling_type,
      29             :       "The type of simulation.  For simulations involving Mechanical deformations, you will need "
      30             :       "to supply the correct Biot coefficient.  For simulations involving Thermal flows, you will "
      31             :       "need an associated ConstantThermalExpansionCoefficient Material");
      32        1422 :   params.addRangeCheckedParam<Real>(
      33         948 :       "biot_coefficient", 1.0, "biot_coefficient>=0 & biot_coefficient<=1", "Biot coefficient");
      34         948 :   params.addParam<bool>("multiply_by_density",
      35         948 :                         true,
      36             :                         "If true, then this Kernel is the time derivative of the fluid mass.  If "
      37             :                         "false, then this Kernel is the derivative of the fluid volume (which is "
      38             :                         "common in poro-mechanics)");
      39         948 :   params.addRequiredParam<UserObjectName>(
      40             :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names.");
      41         474 :   params.addClassDescription(
      42             :       "Fully-saturated version of the single-component, single-phase fluid mass derivative wrt "
      43             :       "time.");
      44         474 :   return params;
      45         474 : }
      46             : 
      47             : template <bool is_ad>
      48         250 : PorousFlowFullySaturatedMassTimeDerivativeTempl<
      49             :     is_ad>::PorousFlowFullySaturatedMassTimeDerivativeTempl(const InputParameters & parameters)
      50             :   : GenericKernel<is_ad>(parameters),
      51         250 :     _dictator(this->template getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      52         250 :     _var_is_porflow_var(_dictator.isPorousFlowVariable(_var.number())),
      53         500 :     _multiply_by_density(this->template getParam<bool>("multiply_by_density")),
      54         250 :     _coupling_type(
      55         250 :         this->template getParam<MooseEnum>("coupling_type").template getEnum<CouplingTypeEnum>()),
      56         250 :     _includes_thermal(_coupling_type == CouplingTypeEnum::ThermoHydro ||
      57             :                       _coupling_type == CouplingTypeEnum::ThermoHydroMechanical),
      58         250 :     _includes_mechanical(_coupling_type == CouplingTypeEnum::HydroMechanical ||
      59             :                          _coupling_type == CouplingTypeEnum::ThermoHydroMechanical),
      60         500 :     _biot_coefficient(this->template getParam<Real>("biot_coefficient")),
      61         500 :     _biot_modulus(this->template getMaterialProperty<Real>("PorousFlow_constant_biot_modulus_qp")),
      62         285 :     _thermal_coeff(_includes_thermal ? &this->template getMaterialProperty<Real>(
      63             :                                            "PorousFlow_constant_thermal_expansion_coefficient_qp")
      64             :                                      : nullptr),
      65         500 :     _fluid_density(_multiply_by_density
      66         314 :                        ? &this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
      67             :                              "PorousFlow_fluid_phase_density_qp")
      68             :                        : nullptr),
      69         483 :     _dfluid_density_dvar(_multiply_by_density && !is_ad
      70         280 :                              ? &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
      71             :                                    "dPorousFlow_fluid_phase_density_qp_dvar")
      72             :                              : nullptr),
      73         500 :     _pp(this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
      74             :         "PorousFlow_porepressure_qp")),
      75         500 :     _pp_old(this->template getMaterialPropertyOld<std::vector<Real>>("PorousFlow_porepressure_qp")),
      76         250 :     _dpp_dvar(is_ad ? nullptr
      77         233 :                     : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
      78             :                           "dPorousFlow_porepressure_qp_dvar")),
      79         285 :     _temperature(_includes_thermal ? &this->template getGenericMaterialProperty<Real, is_ad>(
      80             :                                          "PorousFlow_temperature_qp")
      81             :                                    : nullptr),
      82         500 :     _temperature_old(_includes_thermal
      83         285 :                          ? &this->template getMaterialPropertyOld<Real>("PorousFlow_temperature_qp")
      84             :                          : nullptr),
      85         483 :     _dtemperature_dvar(_includes_thermal && !is_ad
      86         268 :                            ? &this->template getMaterialProperty<std::vector<Real>>(
      87             :                                  "dPorousFlow_temperature_qp_dvar")
      88             :                            : nullptr),
      89         397 :     _strain_rate(_includes_mechanical ? &this->template getMaterialProperty<Real>(
      90             :                                             "PorousFlow_volumetric_strain_rate_qp")
      91             :                                       : nullptr),
      92         483 :     _dstrain_rate_dvar(_includes_mechanical && !is_ad
      93         380 :                            ? &this->template getMaterialProperty<std::vector<RealGradient>>(
      94             :                                  "dPorousFlow_volumetric_strain_rate_qp_dvar")
      95         250 :                            : nullptr)
      96             : {
      97         250 :   if (_dictator.numComponents() != 1 || _dictator.numPhases() != 1)
      98           0 :     mooseError("PorousFlowFullySaturatedMassTimeDerivative is only applicable to single-phase, "
      99             :                "single-component fluid-flow problems.  The Dictator proclaims that you have more "
     100             :                "than one phase or more than one fluid component.  The Dictator does not take such "
     101             :                "mistakes lightly");
     102         250 : }
     103             : 
     104             : template <bool is_ad>
     105             : GenericReal<is_ad>
     106    10080416 : PorousFlowFullySaturatedMassTimeDerivativeTempl<is_ad>::computeQpResidual()
     107             : {
     108             :   const unsigned phase = 0;
     109    10080416 :   GenericReal<is_ad> volume = (_pp[_qp][phase] - _pp_old[_qp][phase]) / _dt / _biot_modulus[_qp];
     110    10080416 :   if (_includes_thermal)
     111     1487744 :     volume -= (*_thermal_coeff)[_qp] * ((*_temperature)[_qp] - (*_temperature_old)[_qp]) / _dt;
     112    10080416 :   if (_includes_mechanical)
     113     5459296 :     volume += _biot_coefficient * (*_strain_rate)[_qp];
     114    10080416 :   if (_multiply_by_density)
     115     1432176 :     return _test[_i][_qp] * (*_fluid_density)[_qp][phase] * volume;
     116     8676928 :   return _test[_i][_qp] * volume;
     117             : }
     118             : 
     119             : template <bool is_ad>
     120             : Real
     121    48153184 : PorousFlowFullySaturatedMassTimeDerivativeTempl<is_ad>::computeQpJacobian()
     122             : {
     123             :   if constexpr (!is_ad)
     124             :   {
     125    48153184 :     if (!_var_is_porflow_var)
     126             :       return 0.0;
     127    48153184 :     return computeQpJac(_dictator.porousFlowVariableNum(_var.number()));
     128             :   }
     129           0 :   return 0.0;
     130             : }
     131             : 
     132             : template <bool is_ad>
     133             : Real
     134    81463808 : PorousFlowFullySaturatedMassTimeDerivativeTempl<is_ad>::computeQpOffDiagJacobian(unsigned int jvar)
     135             : {
     136             :   if constexpr (!is_ad)
     137             :   {
     138    81463808 :     if (_dictator.notPorousFlowVariable(jvar))
     139             :       return 0.0;
     140    81463808 :     return computeQpJac(_dictator.porousFlowVariableNum(jvar));
     141             :   }
     142             :   else
     143             :     libmesh_ignore(jvar);
     144           0 :   return 0.0;
     145             : }
     146             : 
     147             : template <bool is_ad>
     148             : Real
     149   129616992 : PorousFlowFullySaturatedMassTimeDerivativeTempl<is_ad>::computeQpJac(unsigned int pvar)
     150             : {
     151             :   if constexpr (!is_ad)
     152             :   {
     153             :     const unsigned phase = 0;
     154   129616992 :     Real volume = (_pp[_qp][phase] - _pp_old[_qp][phase]) / _dt / _biot_modulus[_qp];
     155   129616992 :     Real dvolume = (*_dpp_dvar)[_qp][phase][pvar] / _dt / _biot_modulus[_qp] * _phi[_j][_qp];
     156   129616992 :     if (_includes_thermal)
     157             :     {
     158    15719680 :       volume -= (*_thermal_coeff)[_qp] * ((*_temperature)[_qp] - (*_temperature_old)[_qp]) / _dt;
     159    15719680 :       dvolume -= (*_thermal_coeff)[_qp] * (*_dtemperature_dvar)[_qp][pvar] / _dt * _phi[_j][_qp];
     160             :     }
     161   129616992 :     if (_includes_mechanical)
     162             :     {
     163    99047168 :       volume += _biot_coefficient * (*_strain_rate)[_qp];
     164    99047168 :       dvolume += _biot_coefficient * (*_dstrain_rate_dvar)[_qp][pvar] * _grad_phi[_j][_qp];
     165             :     }
     166   129616992 :     if (_multiply_by_density)
     167    31202304 :       return _test[_i][_qp] * ((*_fluid_density)[_qp][phase] * dvolume +
     168    31202304 :                                (*_dfluid_density_dvar)[_qp][phase][pvar] * _phi[_j][_qp] * volume);
     169    98414688 :     return _test[_i][_qp] * dvolume;
     170             :   }
     171             :   else
     172             :     libmesh_ignore(pvar);
     173           0 :   return 0.0;
     174             : }
     175             : 
     176             : template class PorousFlowFullySaturatedMassTimeDerivativeTempl<false>;
     177             : template class PorousFlowFullySaturatedMassTimeDerivativeTempl<true>;