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         926 : PorousFlowFullySaturatedMassTimeDerivative::validParams()
      18             : {
      19         926 :   InputParameters params = TimeKernel::validParams();
      20        1852 :   MooseEnum coupling_type("Hydro ThermoHydro HydroMechanical ThermoHydroMechanical", "Hydro");
      21        1852 :   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        2778 :   params.addRangeCheckedParam<Real>(
      28        1852 :       "biot_coefficient", 1.0, "biot_coefficient>=0 & biot_coefficient<=1", "Biot coefficient");
      29        1852 :   params.addParam<bool>("multiply_by_density",
      30        1852 :                         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        1852 :   params.addRequiredParam<UserObjectName>(
      35             :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names.");
      36         926 :   params.addClassDescription("Fully-saturated version of the single-component, single-phase fluid "
      37             :                              "mass derivative wrt time");
      38         926 :   return params;
      39         926 : }
      40             : 
      41         497 : PorousFlowFullySaturatedMassTimeDerivative::PorousFlowFullySaturatedMassTimeDerivative(
      42         497 :     const InputParameters & parameters)
      43             :   : TimeKernel(parameters),
      44         497 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      45         497 :     _var_is_porflow_var(_dictator.isPorousFlowVariable(_var.number())),
      46         994 :     _multiply_by_density(getParam<bool>("multiply_by_density")),
      47         994 :     _coupling_type(getParam<MooseEnum>("coupling_type").getEnum<CouplingTypeEnum>()),
      48         497 :     _includes_thermal(_coupling_type == CouplingTypeEnum::ThermoHydro ||
      49             :                       _coupling_type == CouplingTypeEnum::ThermoHydroMechanical),
      50         497 :     _includes_mechanical(_coupling_type == CouplingTypeEnum::HydroMechanical ||
      51             :                          _coupling_type == CouplingTypeEnum::ThermoHydroMechanical),
      52         994 :     _biot_coefficient(getParam<Real>("biot_coefficient")),
      53         994 :     _biot_modulus(getMaterialProperty<Real>("PorousFlow_constant_biot_modulus_qp")),
      54         564 :     _thermal_coeff(_includes_thermal ? &getMaterialProperty<Real>(
      55             :                                            "PorousFlow_constant_thermal_expansion_coefficient_qp")
      56             :                                      : nullptr),
      57         596 :     _fluid_density(_multiply_by_density ? &getMaterialProperty<std::vector<Real>>(
      58             :                                               "PorousFlow_fluid_phase_density_qp")
      59             :                                         : nullptr),
      60         994 :     _dfluid_density_dvar(_multiply_by_density
      61         596 :                              ? &getMaterialProperty<std::vector<std::vector<Real>>>(
      62             :                                    "dPorousFlow_fluid_phase_density_qp_dvar")
      63             :                              : nullptr),
      64         994 :     _pp(getMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp")),
      65         994 :     _pp_old(getMaterialPropertyOld<std::vector<Real>>("PorousFlow_porepressure_qp")),
      66         497 :     _dpp_dvar(
      67         497 :         getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_porepressure_qp_dvar")),
      68         564 :     _temperature(_includes_thermal ? &getMaterialProperty<Real>("PorousFlow_temperature_qp")
      69             :                                    : nullptr),
      70         564 :     _temperature_old(_includes_thermal ? &getMaterialPropertyOld<Real>("PorousFlow_temperature_qp")
      71             :                                        : nullptr),
      72         564 :     _dtemperature_dvar(_includes_thermal ? &getMaterialProperty<std::vector<Real>>(
      73             :                                                "dPorousFlow_temperature_qp_dvar")
      74             :                                          : nullptr),
      75         994 :     _strain_rate(_includes_mechanical
      76         816 :                      ? &getMaterialProperty<Real>("PorousFlow_volumetric_strain_rate_qp")
      77             :                      : nullptr),
      78         816 :     _dstrain_rate_dvar(_includes_mechanical ? &getMaterialProperty<std::vector<RealGradient>>(
      79             :                                                   "dPorousFlow_volumetric_strain_rate_qp_dvar")
      80         497 :                                             : nullptr)
      81             : {
      82         497 :   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         497 : }
      88             : 
      89             : Real
      90    14254112 : PorousFlowFullySaturatedMassTimeDerivative::computeQpResidual()
      91             : {
      92             :   const unsigned phase = 0;
      93    14254112 :   Real volume = (_pp[_qp][phase] - _pp_old[_qp][phase]) / _dt / _biot_modulus[_qp];
      94    14254112 :   if (_includes_thermal)
      95     1844800 :     volume -= (*_thermal_coeff)[_qp] * ((*_temperature)[_qp] - (*_temperature_old)[_qp]) / _dt;
      96    14254112 :   if (_includes_mechanical)
      97     8218560 :     volume += _biot_coefficient * (*_strain_rate)[_qp];
      98    14254112 :   if (_multiply_by_density)
      99     2058496 :     return _test[_i][_qp] * (*_fluid_density)[_qp][phase] * volume;
     100    12195616 :   return _test[_i][_qp] * volume;
     101             : }
     102             : 
     103             : Real
     104    67788768 : PorousFlowFullySaturatedMassTimeDerivative::computeQpJacobian()
     105             : {
     106             :   // If the variable is not a PorousFlow variable (very unusual), the diag Jacobian terms are 0
     107    67788768 :   if (!_var_is_porflow_var)
     108             :     return 0.0;
     109    67788768 :   return computeQpJac(_dictator.porousFlowVariableNum(_var.number()));
     110             : }
     111             : 
     112             : Real
     113   121073152 : PorousFlowFullySaturatedMassTimeDerivative::computeQpOffDiagJacobian(unsigned int jvar)
     114             : {
     115             :   // If the variable is not a PorousFlow variable, the OffDiag Jacobian terms are 0
     116   121073152 :   if (_dictator.notPorousFlowVariable(jvar))
     117             :     return 0.0;
     118   121073152 :   return computeQpJac(_dictator.porousFlowVariableNum(jvar));
     119             : }
     120             : 
     121             : Real
     122   188861920 : PorousFlowFullySaturatedMassTimeDerivative::computeQpJac(unsigned int pvar)
     123             : {
     124             :   const unsigned phase = 0;
     125   188861920 :   Real volume = (_pp[_qp][phase] - _pp_old[_qp][phase]) / _dt / _biot_modulus[_qp];
     126   188861920 :   Real dvolume = _dpp_dvar[_qp][phase][pvar] / _dt / _biot_modulus[_qp] * _phi[_j][_qp];
     127   188861920 :   if (_includes_thermal)
     128             :   {
     129    19343360 :     volume -= (*_thermal_coeff)[_qp] * ((*_temperature)[_qp] - (*_temperature_old)[_qp]) / _dt;
     130    19343360 :     dvolume -= (*_thermal_coeff)[_qp] * (*_dtemperature_dvar)[_qp][pvar] / _dt * _phi[_j][_qp];
     131             :   }
     132   188861920 :   if (_includes_mechanical)
     133             :   {
     134   149857536 :     volume += _biot_coefficient * (*_strain_rate)[_qp];
     135   149857536 :     dvolume += _biot_coefficient * (*_dstrain_rate_dvar)[_qp][pvar] * _grad_phi[_j][_qp];
     136             :   }
     137   188861920 :   if (_multiply_by_density)
     138    46775584 :     return _test[_i][_qp] * ((*_fluid_density)[_qp][phase] * dvolume +
     139    46775584 :                              (*_dfluid_density_dvar)[_qp][phase][pvar] * _phi[_j][_qp] * volume);
     140   142086336 :   return _test[_i][_qp] * dvolume;
     141             : }