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 "PorousFlowEnergyTimeDerivative.h"
      11             : 
      12             : #include "MooseVariable.h"
      13             : 
      14             : registerMooseObject("PorousFlowApp", PorousFlowEnergyTimeDerivative);
      15             : registerMooseObject("PorousFlowApp", ADPorousFlowEnergyTimeDerivative);
      16             : 
      17             : template <bool is_ad>
      18             : InputParameters
      19        2076 : PorousFlowEnergyTimeDerivativeTempl<is_ad>::validParams()
      20             : {
      21             :   InputParameters params = GenericKernel<is_ad>::validParams();
      22        2102 :   params.set<MultiMooseEnum>("vector_tags") = "time";
      23        2102 :   params.set<MultiMooseEnum>("matrix_tags") = "system time";
      24        2102 :   params.addParam<bool>("strain_at_nearest_qp",
      25        2102 :                         false,
      26             :                         "When calculating nodal porosity that depends on strain, use the strain at "
      27             :                         "the nearest quadpoint.  This adds a small extra computational burden, and "
      28             :                         "is not necessary for simulations involving only linear lagrange elements. "
      29             :                         " If you set this to true, you will also want to set the same parameter to "
      30             :                         "true for related Kernels and Materials");
      31        2102 :   params.addParam<std::string>(
      32             :       "base_name",
      33             :       "For mechanically-coupled systems, this Kernel will depend on the volumetric strain.  "
      34             :       "base_name should almost always be the same base_name as given to the TensorMechanics object "
      35             :       "that computes strain.  Supplying a base_name to this Kernel but not defining an associated "
      36             :       "TensorMechanics strain calculator means that this Kernel will not depend on volumetric "
      37             :       "strain.  That could be useful when models contain solid mechanics that is not coupled to "
      38             :       "porous flow, for example");
      39        2102 :   params.addRequiredParam<UserObjectName>(
      40             :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names.");
      41        1051 :   params.set<bool>("use_displaced_mesh") = false;
      42        1051 :   params.suppressParameter<bool>("use_displaced_mesh");
      43        1051 :   params.addClassDescription("Derivative of heat-energy-density wrt time");
      44        1051 :   return params;
      45           0 : }
      46             : 
      47             : template <bool is_ad>
      48         564 : PorousFlowEnergyTimeDerivativeTempl<is_ad>::PorousFlowEnergyTimeDerivativeTempl(
      49             :     const InputParameters & parameters)
      50             :   : PorousFlowLumpedKernelBaseTempl<is_ad>(parameters),
      51         564 :     _dictator(this->template getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      52         564 :     _var_is_porflow_var(_dictator.isPorousFlowVariable(_var.number())),
      53         564 :     _num_phases(_dictator.numPhases()),
      54         564 :     _fluid_present(_num_phases > 0),
      55        1128 :     _strain_at_nearest_qp(this->template getParam<bool>("strain_at_nearest_qp")),
      56        1188 :     _base_name(this->isParamValid("base_name")
      57         564 :                    ? this->template getParam<std::string>("base_name") + "_"
      58             :                    : ""),
      59         564 :     _has_total_strain(
      60         564 :         this->template hasMaterialProperty<RankTwoTensor>(_base_name + "total_strain")),
      61         697 :     _total_strain_old(_has_total_strain ? &this->template getMaterialPropertyOld<RankTwoTensor>(
      62             :                                               _base_name + "total_strain")
      63             :                                         : nullptr),
      64        1128 :     _porosity(this->template getGenericMaterialProperty<Real, is_ad>("PorousFlow_porosity_nodal")),
      65        1128 :     _porosity_old(this->template getMaterialPropertyOld<Real>("PorousFlow_porosity_nodal")),
      66         564 :     _dporosity_dvar(is_ad ? nullptr
      67         550 :                           : &this->template getMaterialProperty<std::vector<Real>>(
      68             :                                 "dPorousFlow_porosity_nodal_dvar")),
      69         564 :     _dporosity_dgradvar(is_ad ? nullptr
      70         550 :                               : &this->template getMaterialProperty<std::vector<RealGradient>>(
      71             :                                     "dPorousFlow_porosity_nodal_dgradvar")),
      72         571 :     _nearest_qp(_strain_at_nearest_qp ? &this->template getMaterialProperty<unsigned int>(
      73             :                                             "PorousFlow_nearestqp_nodal")
      74             :                                       : nullptr),
      75        1128 :     _rock_energy_nodal(this->template getGenericMaterialProperty<Real, is_ad>(
      76             :         "PorousFlow_matrix_internal_energy_nodal")),
      77         564 :     _rock_energy_nodal_old(
      78         564 :         this->template getMaterialPropertyOld<Real>("PorousFlow_matrix_internal_energy_nodal")),
      79         564 :     _drock_energy_nodal_dvar(is_ad ? nullptr
      80         550 :                                    : &this->template getMaterialProperty<std::vector<Real>>(
      81             :                                          "dPorousFlow_matrix_internal_energy_nodal_dvar")),
      82        1128 :     _fluid_density(_fluid_present
      83        1074 :                        ? &this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
      84             :                              "PorousFlow_fluid_phase_density_nodal")
      85             :                        : nullptr),
      86        1074 :     _fluid_density_old(_fluid_present ? &this->template getMaterialPropertyOld<std::vector<Real>>(
      87             :                                             "PorousFlow_fluid_phase_density_nodal")
      88             :                                       : nullptr),
      89        1114 :     _dfluid_density_dvar(_fluid_present && !is_ad
      90        1053 :                              ? &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
      91             :                                    "dPorousFlow_fluid_phase_density_nodal_dvar")
      92             :                              : nullptr),
      93         564 :     _fluid_saturation_nodal(
      94        1074 :         _fluid_present ? &this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
      95             :                              "PorousFlow_saturation_nodal")
      96             :                        : nullptr),
      97        1128 :     _fluid_saturation_nodal_old(_fluid_present
      98        1074 :                                     ? &this->template getMaterialPropertyOld<std::vector<Real>>(
      99             :                                           "PorousFlow_saturation_nodal")
     100             :                                     : nullptr),
     101         564 :     _dfluid_saturation_nodal_dvar(
     102         550 :         _fluid_present && !is_ad
     103        1053 :             ? &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
     104             :                   "dPorousFlow_saturation_nodal_dvar")
     105             :             : nullptr),
     106        1128 :     _energy_nodal(_fluid_present
     107        1074 :                       ? &this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
     108             :                             "PorousFlow_fluid_phase_internal_energy_nodal")
     109             :                       : nullptr),
     110        1074 :     _energy_nodal_old(_fluid_present ? &this->template getMaterialPropertyOld<std::vector<Real>>(
     111             :                                            "PorousFlow_fluid_phase_internal_energy_nodal")
     112             :                                      : nullptr),
     113        1114 :     _denergy_nodal_dvar(_fluid_present && !is_ad
     114        1053 :                             ? &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
     115             :                                   "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")
     116         564 :                             : nullptr)
     117             : {
     118         564 : }
     119             : 
     120             : template <bool is_ad>
     121             : GenericReal<is_ad>
     122     6022744 : PorousFlowEnergyTimeDerivativeTempl<is_ad>::computeQpResidual()
     123             : {
     124             :   /// Porous matrix heat energy
     125     6022744 :   GenericReal<is_ad> energy = (1.0 - _porosity[_i]) * _rock_energy_nodal[_i];
     126     6022744 :   Real energy_old = (1.0 - _porosity_old[_i]) * _rock_energy_nodal_old[_i];
     127             : 
     128             :   /// Add the fluid heat energy
     129     6022744 :   if (_fluid_present)
     130    12375576 :     for (unsigned ph = 0; ph < _num_phases; ++ph)
     131             :     {
     132     6383696 :       energy += (*_fluid_density)[_i][ph] * (*_fluid_saturation_nodal)[_i][ph] *
     133     6381520 :                 (*_energy_nodal)[_i][ph] * _porosity[_i];
     134     6381520 :       energy_old += (*_fluid_density_old)[_i][ph] * (*_fluid_saturation_nodal_old)[_i][ph] *
     135     6381520 :                     (*_energy_nodal_old)[_i][ph] * _porosity_old[_i];
     136             :     }
     137     6022744 :   const Real strain = (_has_total_strain ? (*_total_strain_old)[_qp].trace() : 0.0);
     138             : 
     139     6026328 :   return _test[_i][_qp] * (1.0 + strain) * (energy - energy_old) / _dt;
     140             : }
     141             : 
     142             : template <bool is_ad>
     143             : Real
     144    20227116 : PorousFlowEnergyTimeDerivativeTempl<is_ad>::computeQpJacobian()
     145             : {
     146             :   if constexpr (!is_ad)
     147             :   {
     148             :     // If the variable is not a PorousFlow variable (very unusual), the diag Jacobian terms are 0
     149    20227116 :     if (!_var_is_porflow_var)
     150             :       return 0.0;
     151    20227116 :     return computeQpJac(_dictator.porousFlowVariableNum(_var.number()));
     152             :   }
     153           0 :   return 0.0;
     154             : }
     155             : 
     156             : template <bool is_ad>
     157             : Real
     158    22121196 : PorousFlowEnergyTimeDerivativeTempl<is_ad>::computeQpOffDiagJacobian(unsigned int jvar)
     159             : {
     160             :   if constexpr (!is_ad)
     161             :   {
     162             :     // If the variable is not a PorousFlow variable, the OffDiag Jacobian terms are 0
     163    22121196 :     if (_dictator.notPorousFlowVariable(jvar))
     164             :       return 0.0;
     165    22121196 :     return computeQpJac(_dictator.porousFlowVariableNum(jvar));
     166             :   }
     167             :   else
     168             :     libmesh_ignore(jvar);
     169           0 :   return 0.0;
     170             : }
     171             : 
     172             : template <bool is_ad>
     173             : Real
     174    42348312 : PorousFlowEnergyTimeDerivativeTempl<is_ad>::computeQpJac(unsigned int pvar) const
     175             : {
     176             :   if constexpr (!is_ad)
     177             :   {
     178    42348312 :     const unsigned nearest_qp = (_strain_at_nearest_qp ? (*_nearest_qp)[_i] : _i);
     179             : 
     180    42348312 :     const Real strain = (_has_total_strain ? (*_total_strain_old)[_qp].trace() : 0.0);
     181             : 
     182             :     // porosity is dependent on variables that are lumped to the nodes,
     183             :     // but it can depend on the gradient
     184             :     // of variables, which are NOT lumped to the nodes, hence:
     185    42348312 :     Real denergy = -(*_dporosity_dgradvar)[_i][pvar] * _grad_phi[_j][_i] * _rock_energy_nodal[_i];
     186    85882656 :     for (unsigned ph = 0; ph < _num_phases; ++ph)
     187    43534344 :       denergy += (*_fluid_density)[_i][ph] * (*_fluid_saturation_nodal)[_i][ph] *
     188    43534344 :                  (*_energy_nodal)[_i][ph] * (*_dporosity_dgradvar)[_i][pvar] *
     189    43534344 :                  _grad_phi[_j][nearest_qp];
     190             : 
     191    42348312 :     if (_i != _j)
     192    33964504 :       return _test[_i][_qp] * (1.0 + strain) * denergy / _dt;
     193             : 
     194             :     // As the fluid energy is lumped to the nodes, only non-zero terms are for _i==_j
     195     8383808 :     denergy += -(*_dporosity_dvar)[_i][pvar] * _rock_energy_nodal[_i];
     196     8383808 :     denergy += (1.0 - _porosity[_i]) * (*_drock_energy_nodal_dvar)[_i][pvar];
     197    17262392 :     for (unsigned ph = 0; ph < _num_phases; ++ph)
     198             :     {
     199     8878584 :       denergy += (*_dfluid_density_dvar)[_i][ph][pvar] * (*_fluid_saturation_nodal)[_i][ph] *
     200     8878584 :                  (*_energy_nodal)[_i][ph] * _porosity[_i];
     201     8878584 :       denergy += (*_fluid_density)[_i][ph] * (*_dfluid_saturation_nodal_dvar)[_i][ph][pvar] *
     202     8878584 :                  (*_energy_nodal)[_i][ph] * _porosity[_i];
     203     8878584 :       denergy += (*_fluid_density)[_i][ph] * (*_fluid_saturation_nodal)[_i][ph] *
     204     8878584 :                  (*_denergy_nodal_dvar)[_i][ph][pvar] * _porosity[_i];
     205     8878584 :       denergy += (*_fluid_density)[_i][ph] * (*_fluid_saturation_nodal)[_i][ph] *
     206     8878584 :                  (*_energy_nodal)[_i][ph] * (*_dporosity_dvar)[_i][pvar];
     207             :     }
     208     8383808 :     return _test[_i][_qp] * (1.0 + strain) * denergy / _dt;
     209             :   }
     210             :   else
     211             :     libmesh_ignore(pvar);
     212           0 :   return 0.0;
     213             : }
     214             : 
     215             : template class PorousFlowEnergyTimeDerivativeTempl<false>;
     216             : template class PorousFlowEnergyTimeDerivativeTempl<true>;