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             : 
      16             : InputParameters
      17         442 : PorousFlowFullySaturatedMassTimeDerivative::validParams()
      18             : {
      19         442 :   InputParameters params = TimeKernel::validParams();
      20         884 :   MooseEnum coupling_type("Hydro ThermoHydro HydroMechanical ThermoHydroMechanical", "Hydro");
      21         884 :   params.addParam<MooseEnum>("coupling_type",
      22             :                              coupling_type,
      23             :                              "The type of simulation.  For simulations involving Mechanical "
      24             :                              "deformations, you will need to supply the correct Biot coefficient.  "
      25             :                              "For simulations involving Thermal flows, you will need an associated "
      26             :                              "ConstantThermalExpansionCoefficient Material");
      27        1326 :   params.addRangeCheckedParam<Real>(
      28         884 :       "biot_coefficient", 1.0, "biot_coefficient>=0 & biot_coefficient<=1", "Biot coefficient");
      29         884 :   params.addParam<bool>("multiply_by_density",
      30         884 :                         true,
      31             :                         "If true, then this Kernel is the time derivative of the fluid "
      32             :                         "mass.  If false, then this Kernel is the derivative of the "
      33             :                         "fluid volume (which is common in poro-mechanics)");
      34         884 :   params.addRequiredParam<UserObjectName>(
      35             :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names.");
      36         442 :   params.addClassDescription("Fully-saturated version of the single-component, single-phase fluid "
      37             :                              "mass derivative wrt time");
      38         442 :   return params;
      39         442 : }
      40             : 
      41         233 : PorousFlowFullySaturatedMassTimeDerivative::PorousFlowFullySaturatedMassTimeDerivative(
      42         233 :     const InputParameters & parameters)
      43             :   : TimeKernel(parameters),
      44         233 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      45         233 :     _var_is_porflow_var(_dictator.isPorousFlowVariable(_var.number())),
      46         466 :     _multiply_by_density(getParam<bool>("multiply_by_density")),
      47         466 :     _coupling_type(getParam<MooseEnum>("coupling_type").getEnum<CouplingTypeEnum>()),
      48         233 :     _includes_thermal(_coupling_type == CouplingTypeEnum::ThermoHydro ||
      49             :                       _coupling_type == CouplingTypeEnum::ThermoHydroMechanical),
      50         233 :     _includes_mechanical(_coupling_type == CouplingTypeEnum::HydroMechanical ||
      51             :                          _coupling_type == CouplingTypeEnum::ThermoHydroMechanical),
      52         466 :     _biot_coefficient(getParam<Real>("biot_coefficient")),
      53         466 :     _biot_modulus(getMaterialProperty<Real>("PorousFlow_constant_biot_modulus_qp")),
      54         268 :     _thermal_coeff(_includes_thermal ? &getMaterialProperty<Real>(
      55             :                                            "PorousFlow_constant_thermal_expansion_coefficient_qp")
      56             :                                      : nullptr),
      57         280 :     _fluid_density(_multiply_by_density ? &getMaterialProperty<std::vector<Real>>(
      58             :                                               "PorousFlow_fluid_phase_density_qp")
      59             :                                         : nullptr),
      60         466 :     _dfluid_density_dvar(_multiply_by_density
      61         280 :                              ? &getMaterialProperty<std::vector<std::vector<Real>>>(
      62             :                                    "dPorousFlow_fluid_phase_density_qp_dvar")
      63             :                              : nullptr),
      64         466 :     _pp(getMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp")),
      65         466 :     _pp_old(getMaterialPropertyOld<std::vector<Real>>("PorousFlow_porepressure_qp")),
      66         233 :     _dpp_dvar(
      67         233 :         getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_porepressure_qp_dvar")),
      68         268 :     _temperature(_includes_thermal ? &getMaterialProperty<Real>("PorousFlow_temperature_qp")
      69             :                                    : nullptr),
      70         268 :     _temperature_old(_includes_thermal ? &getMaterialPropertyOld<Real>("PorousFlow_temperature_qp")
      71             :                                        : nullptr),
      72         268 :     _dtemperature_dvar(_includes_thermal ? &getMaterialProperty<std::vector<Real>>(
      73             :                                                "dPorousFlow_temperature_qp_dvar")
      74             :                                          : nullptr),
      75         466 :     _strain_rate(_includes_mechanical
      76         380 :                      ? &getMaterialProperty<Real>("PorousFlow_volumetric_strain_rate_qp")
      77             :                      : nullptr),
      78         380 :     _dstrain_rate_dvar(_includes_mechanical ? &getMaterialProperty<std::vector<RealGradient>>(
      79             :                                                   "dPorousFlow_volumetric_strain_rate_qp_dvar")
      80         233 :                                             : nullptr)
      81             : {
      82         233 :   if (_dictator.numComponents() != 1 || _dictator.numPhases() != 1)
      83           0 :     mooseError("PorousFlowFullySaturatedMassTimeDerivative is only applicable to single-phase, "
      84             :                "single-component fluid-flow problems.  The Dictator proclaims that you have more "
      85             :                "than one phase or more than one fluid component.  The Dictator does not take such "
      86             :                "mistakes lightly");
      87         233 : }
      88             : 
      89             : Real
      90    10051728 : PorousFlowFullySaturatedMassTimeDerivative::computeQpResidual()
      91             : {
      92             :   const unsigned phase = 0;
      93    10051728 :   Real volume = (_pp[_qp][phase] - _pp_old[_qp][phase]) / _dt / _biot_modulus[_qp];
      94    10051728 :   if (_includes_thermal)
      95     1487744 :     volume -= (*_thermal_coeff)[_qp] * ((*_temperature)[_qp] - (*_temperature_old)[_qp]) / _dt;
      96    10051728 :   if (_includes_mechanical)
      97     5459296 :     volume += _biot_coefficient * (*_strain_rate)[_qp];
      98    10051728 :   if (_multiply_by_density)
      99     1374800 :     return _test[_i][_qp] * (*_fluid_density)[_qp][phase] * volume;
     100     8676928 :   return _test[_i][_qp] * volume;
     101             : }
     102             : 
     103             : Real
     104    48153184 : PorousFlowFullySaturatedMassTimeDerivative::computeQpJacobian()
     105             : {
     106             :   // If the variable is not a PorousFlow variable (very unusual), the diag Jacobian terms are 0
     107    48153184 :   if (!_var_is_porflow_var)
     108             :     return 0.0;
     109    48153184 :   return computeQpJac(_dictator.porousFlowVariableNum(_var.number()));
     110             : }
     111             : 
     112             : Real
     113    81463808 : PorousFlowFullySaturatedMassTimeDerivative::computeQpOffDiagJacobian(unsigned int jvar)
     114             : {
     115             :   // If the variable is not a PorousFlow variable, the OffDiag Jacobian terms are 0
     116    81463808 :   if (_dictator.notPorousFlowVariable(jvar))
     117             :     return 0.0;
     118    81463808 :   return computeQpJac(_dictator.porousFlowVariableNum(jvar));
     119             : }
     120             : 
     121             : Real
     122   129616992 : PorousFlowFullySaturatedMassTimeDerivative::computeQpJac(unsigned int pvar)
     123             : {
     124             :   const unsigned phase = 0;
     125   129616992 :   Real volume = (_pp[_qp][phase] - _pp_old[_qp][phase]) / _dt / _biot_modulus[_qp];
     126   129616992 :   Real dvolume = _dpp_dvar[_qp][phase][pvar] / _dt / _biot_modulus[_qp] * _phi[_j][_qp];
     127   129616992 :   if (_includes_thermal)
     128             :   {
     129    15719680 :     volume -= (*_thermal_coeff)[_qp] * ((*_temperature)[_qp] - (*_temperature_old)[_qp]) / _dt;
     130    15719680 :     dvolume -= (*_thermal_coeff)[_qp] * (*_dtemperature_dvar)[_qp][pvar] / _dt * _phi[_j][_qp];
     131             :   }
     132   129616992 :   if (_includes_mechanical)
     133             :   {
     134    99047168 :     volume += _biot_coefficient * (*_strain_rate)[_qp];
     135    99047168 :     dvolume += _biot_coefficient * (*_dstrain_rate_dvar)[_qp][pvar] * _grad_phi[_j][_qp];
     136             :   }
     137   129616992 :   if (_multiply_by_density)
     138    31202304 :     return _test[_i][_qp] * ((*_fluid_density)[_qp][phase] * dvolume +
     139    31202304 :                              (*_dfluid_density_dvar)[_qp][phase][pvar] * _phi[_j][_qp] * volume);
     140    98414688 :   return _test[_i][_qp] * dvolume;
     141             : }