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 "PorousFlowMassTimeDerivative.h"
      11             : 
      12             : #include "MooseVariable.h"
      13             : 
      14             : #include "libmesh/quadrature.h"
      15             : 
      16             : #include <limits>
      17             : 
      18             : registerMooseObject("PorousFlowApp", PorousFlowMassTimeDerivative);
      19             : registerMooseObject("PorousFlowApp", ADPorousFlowMassTimeDerivative);
      20             : 
      21             : template <bool is_ad>
      22             : InputParameters
      23        6819 : PorousFlowMassTimeDerivativeTempl<is_ad>::validParams()
      24             : {
      25             :   InputParameters params = GenericKernel<is_ad>::validParams();
      26       13638 :   params.set<MultiMooseEnum>("vector_tags") = "time";
      27       13638 :   params.set<MultiMooseEnum>("matrix_tags") = "system time";
      28       13638 :   params.addParam<bool>("strain_at_nearest_qp",
      29       13638 :                         false,
      30             :                         "When calculating nodal porosity that depends on strain, use the strain at "
      31             :                         "the nearest quadpoint.  This adds a small extra computational burden, and "
      32             :                         "is not necessary for simulations involving only linear lagrange elements. "
      33             :                         " If you set this to true, you will also want to set the same parameter to "
      34             :                         "true for related Kernels and Materials");
      35       13638 :   params.addParam<std::string>(
      36             :       "base_name",
      37             :       "For mechanically-coupled systems, this Kernel will depend on the volumetric strain.  "
      38             :       "base_name should almost always be the same base_name as given to the TensorMechanics object "
      39             :       "that computes strain.  Supplying a base_name to this Kernel but not defining an associated "
      40             :       "TensorMechanics strain calculator means that this Kernel will not depend on volumetric "
      41             :       "strain.  That could be useful when models contain solid mechanics that is not coupled to "
      42             :       "porous flow, for example");
      43       13638 :   params.addParam<bool>(
      44             :       "multiply_by_density",
      45       13638 :       true,
      46             :       "If true, then this Kernel represents the time derivative of the fluid mass.  If false, then "
      47             :       "this Kernel represents the time derivative of the fluid volume (care must then be taken "
      48             :       "when using other PorousFlow objects, such as the PorousFlowFluidMass postprocessor).");
      49       13638 :   params.addParam<unsigned int>(
      50       13638 :       "fluid_component", 0, "The index corresponding to the component for this kernel");
      51       13638 :   params.addRequiredParam<UserObjectName>(
      52             :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names.");
      53        6819 :   params.set<bool>("use_displaced_mesh") = false;
      54        6819 :   params.suppressParameter<bool>("use_displaced_mesh");
      55        6819 :   params.addClassDescription("Derivative of fluid-component mass with respect to time.  Mass "
      56             :                              "lumping to the nodes is used.");
      57        6819 :   return params;
      58           0 : }
      59             : 
      60             : template <bool is_ad>
      61        3684 : PorousFlowMassTimeDerivativeTempl<is_ad>::PorousFlowMassTimeDerivativeTempl(
      62             :     const InputParameters & parameters)
      63             :   : PorousFlowLumpedKernelBaseTempl<is_ad>(parameters),
      64        3682 :     _fluid_component(this->template getParam<unsigned int>("fluid_component")),
      65        3682 :     _dictator(this->template getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      66        3682 :     _var_is_porflow_var(_dictator.isPorousFlowVariable(_var.number())),
      67        3682 :     _num_phases(_dictator.numPhases()),
      68        7364 :     _strain_at_nearest_qp(this->template getParam<bool>("strain_at_nearest_qp")),
      69        7364 :     _multiply_by_density(this->template getParam<bool>("multiply_by_density")),
      70        7364 :     _base_name(this->isParamValid("base_name")
      71        3682 :                    ? this->template getParam<std::string>("base_name") + "_"
      72             :                    : ""),
      73        3682 :     _has_total_strain(
      74        3682 :         this->template hasMaterialProperty<RankTwoTensor>(_base_name + "total_strain")),
      75        4003 :     _total_strain_old(_has_total_strain ? &this->template getMaterialPropertyOld<RankTwoTensor>(
      76             :                                               _base_name + "total_strain")
      77             :                                         : nullptr),
      78        7364 :     _porosity(this->template getGenericMaterialProperty<Real, is_ad>("PorousFlow_porosity_nodal")),
      79        7364 :     _porosity_old(this->template getMaterialPropertyOld<Real>("PorousFlow_porosity_nodal")),
      80        3682 :     _dporosity_dvar(is_ad ? nullptr
      81        3661 :                           : &this->template getMaterialProperty<std::vector<Real>>(
      82             :                                 "dPorousFlow_porosity_nodal_dvar")),
      83        3682 :     _dporosity_dgradvar(is_ad ? nullptr
      84        3661 :                               : &this->template getMaterialProperty<std::vector<RealGradient>>(
      85             :                                     "dPorousFlow_porosity_nodal_dgradvar")),
      86        3713 :     _nearest_qp(_strain_at_nearest_qp ? &this->template getMaterialProperty<unsigned int>(
      87             :                                             "PorousFlow_nearestqp_nodal")
      88             :                                       : nullptr),
      89        7364 :     _fluid_density(_multiply_by_density
      90        7269 :                        ? &this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
      91             :                              "PorousFlow_fluid_phase_density_nodal")
      92             :                        : nullptr),
      93        7364 :     _fluid_density_old(_multiply_by_density
      94        7269 :                            ? &this->template getMaterialPropertyOld<std::vector<Real>>(
      95             :                                  "PorousFlow_fluid_phase_density_nodal")
      96             :                            : nullptr),
      97        7343 :     _dfluid_density_dvar(_multiply_by_density && !is_ad
      98        7227 :                              ? &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
      99             :                                    "dPorousFlow_fluid_phase_density_nodal_dvar")
     100             :                              : nullptr),
     101        7364 :     _fluid_saturation_nodal(this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
     102             :         "PorousFlow_saturation_nodal")),
     103        3682 :     _fluid_saturation_nodal_old(
     104        3682 :         this->template getMaterialPropertyOld<std::vector<Real>>("PorousFlow_saturation_nodal")),
     105        3682 :     _dfluid_saturation_nodal_dvar(
     106             :         is_ad ? nullptr
     107        3661 :               : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
     108             :                     "dPorousFlow_saturation_nodal_dvar")),
     109        7364 :     _mass_frac(this->template getGenericMaterialProperty<std::vector<std::vector<Real>>, is_ad>(
     110             :         "PorousFlow_mass_frac_nodal")),
     111        7364 :     _mass_frac_old(this->template getMaterialPropertyOld<std::vector<std::vector<Real>>>(
     112             :         "PorousFlow_mass_frac_nodal")),
     113        3682 :     _dmass_frac_dvar(
     114             :         is_ad ? nullptr
     115        3661 :               : &this->template getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
     116        3684 :                     "dPorousFlow_mass_frac_nodal_dvar"))
     117             : {
     118        3682 :   if (_fluid_component >= _dictator.numComponents())
     119           0 :     this->paramError(
     120             :         "fluid_component",
     121             :         "The Dictator proclaims that the maximum fluid component index in this simulation is ",
     122           0 :         _dictator.numComponents() - 1,
     123             :         " whereas you have used ",
     124           0 :         _fluid_component,
     125             :         ". Remember that indexing starts at 0. The Dictator does not take such mistakes lightly.");
     126        3682 : }
     127             : 
     128             : template <bool is_ad>
     129             : GenericReal<is_ad>
     130    91678288 : PorousFlowMassTimeDerivativeTempl<is_ad>::computeQpResidual()
     131             : {
     132        3264 :   GenericReal<is_ad> mass = 0.0;
     133             :   Real mass_old = 0.0;
     134   189398108 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
     135             :   {
     136        3264 :     const GenericReal<is_ad> dens =
     137    97719820 :         (_multiply_by_density ? (*_fluid_density)[_i][ph] : GenericReal<is_ad>(1.0));
     138    97723084 :     mass += dens * _fluid_saturation_nodal[_i][ph] * _mass_frac[_i][ph][_fluid_component];
     139    97719820 :     const Real dens_old = (_multiply_by_density ? (*_fluid_density_old)[_i][ph] : 1.0);
     140    97719820 :     mass_old +=
     141    97719820 :         dens_old * _fluid_saturation_nodal_old[_i][ph] * _mass_frac_old[_i][ph][_fluid_component];
     142             :   }
     143    91678288 :   const Real strain = (_has_total_strain ? (*_total_strain_old)[_qp].trace() : 0.0);
     144             : 
     145    91678288 :   return _test[_i][_qp] * (1.0 + strain) * (_porosity[_i] * mass - _porosity_old[_i] * mass_old) /
     146    91681552 :          _dt;
     147             : }
     148             : 
     149             : template <bool is_ad>
     150             : Real
     151   415803164 : PorousFlowMassTimeDerivativeTempl<is_ad>::computeQpJacobian()
     152             : {
     153             :   if constexpr (!is_ad)
     154             :   {
     155   415803164 :     if (!_var_is_porflow_var)
     156             :       return 0.0;
     157   415802556 :     return computeQpJac(_dictator.porousFlowVariableNum(_var.number()));
     158             :   }
     159           0 :   return 0.0;
     160             : }
     161             : 
     162             : template <bool is_ad>
     163             : Real
     164   283435148 : PorousFlowMassTimeDerivativeTempl<is_ad>::computeQpOffDiagJacobian(unsigned int jvar)
     165             : {
     166             :   if constexpr (!is_ad)
     167             :   {
     168   283435148 :     if (_dictator.notPorousFlowVariable(jvar))
     169             :       return 0.0;
     170   283433932 :     return computeQpJac(_dictator.porousFlowVariableNum(jvar));
     171             :   }
     172             :   else
     173             :     libmesh_ignore(jvar);
     174           0 :   return 0.0;
     175             : }
     176             : 
     177             : template <bool is_ad>
     178             : Real
     179   699236488 : PorousFlowMassTimeDerivativeTempl<is_ad>::computeQpJac(unsigned int pvar)
     180             : {
     181             :   if constexpr (!is_ad)
     182             :   {
     183   699236488 :     const unsigned nearest_qp = (_strain_at_nearest_qp ? (*_nearest_qp)[_i] : _i);
     184             : 
     185   699236488 :     const Real strain = (_has_total_strain ? (*_total_strain_old)[_qp].trace() : 0.0);
     186             : 
     187             :     Real dmass = 0.0;
     188  1428177272 :     for (unsigned ph = 0; ph < _num_phases; ++ph)
     189             :     {
     190   728940784 :       const Real dens = (_multiply_by_density ? (*_fluid_density)[_i][ph] : 1.0);
     191   728940784 :       dmass += dens * _fluid_saturation_nodal[_i][ph] * _mass_frac[_i][ph][_fluid_component] *
     192   728940784 :                (*_dporosity_dgradvar)[_i][pvar] * _grad_phi[_j][nearest_qp];
     193             :     }
     194             : 
     195   699236488 :     if (_i != _j)
     196   602770160 :       return _test[_i][_qp] * (1.0 + strain) * dmass / _dt;
     197             : 
     198   201430260 :     for (unsigned ph = 0; ph < _num_phases; ++ph)
     199             :     {
     200   104963932 :       if (_multiply_by_density)
     201   104931452 :         dmass += (*_dfluid_density_dvar)[_i][ph][pvar] * _fluid_saturation_nodal[_i][ph] *
     202   104931452 :                  _mass_frac[_i][ph][_fluid_component] * _porosity[_i];
     203   104963932 :       const Real dens = (_multiply_by_density ? (*_fluid_density)[_i][ph] : 1.0);
     204   104963932 :       dmass += dens * (*_dfluid_saturation_nodal_dvar)[_i][ph][pvar] *
     205   104963932 :                _mass_frac[_i][ph][_fluid_component] * _porosity[_i];
     206   104963932 :       dmass += dens * _fluid_saturation_nodal[_i][ph] *
     207   104963932 :                (*_dmass_frac_dvar)[_i][ph][_fluid_component][pvar] * _porosity[_i];
     208   104963932 :       dmass += dens * _fluid_saturation_nodal[_i][ph] * _mass_frac[_i][ph][_fluid_component] *
     209   104963932 :                (*_dporosity_dvar)[_i][pvar];
     210             :     }
     211    96466328 :     return _test[_i][_qp] * (1.0 + strain) * dmass / _dt;
     212             :   }
     213             :   else
     214             :     libmesh_ignore(pvar);
     215           0 :   return 0.0;
     216             : }
     217             : 
     218             : template class PorousFlowMassTimeDerivativeTempl<false>;
     219             : template class PorousFlowMassTimeDerivativeTempl<true>;