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             : 
      20             : InputParameters
      21       13736 : PorousFlowMassTimeDerivative::validParams()
      22             : {
      23       13736 :   InputParameters params = TimeKernel::validParams();
      24       27472 :   params.addParam<bool>("strain_at_nearest_qp",
      25       27472 :                         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       27472 :   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       27472 :   params.addParam<bool>(
      40             :       "multiply_by_density",
      41       27472 :       true,
      42             :       "If true, then this Kernel represents the time derivative of the fluid mass.  If false, then "
      43             :       "this Kernel represents the time derivative of the fluid volume (care must then be taken "
      44             :       "when using other PorousFlow objects, such as the PorousFlowFluidMass postprocessor).");
      45       27472 :   params.addParam<unsigned int>(
      46       27472 :       "fluid_component", 0, "The index corresponding to the component for this kernel");
      47       27472 :   params.addRequiredParam<UserObjectName>(
      48             :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names.");
      49       13736 :   params.set<bool>("use_displaced_mesh") = false;
      50       13736 :   params.suppressParameter<bool>("use_displaced_mesh");
      51       13736 :   params.addClassDescription("Derivative of fluid-component mass with respect to time.  Mass "
      52             :                              "lumping to the nodes is used.");
      53       13736 :   return params;
      54           0 : }
      55             : 
      56        7586 : PorousFlowMassTimeDerivative::PorousFlowMassTimeDerivative(const InputParameters & parameters)
      57             :   : TimeKernel(parameters),
      58        7586 :     _fluid_component(getParam<unsigned int>("fluid_component")),
      59        7586 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      60       15172 :     _var_is_porflow_var(_dictator.isPorousFlowVariable(_var.number())),
      61        7586 :     _num_phases(_dictator.numPhases()),
      62       15172 :     _strain_at_nearest_qp(getParam<bool>("strain_at_nearest_qp")),
      63       15172 :     _multiply_by_density(getParam<bool>("multiply_by_density")),
      64       15172 :     _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
      65        7586 :     _has_total_strain(hasMaterialProperty<RankTwoTensor>(_base_name + "total_strain")),
      66       15172 :     _total_strain_old(_has_total_strain
      67        8259 :                           ? &getMaterialPropertyOld<RankTwoTensor>(_base_name + "total_strain")
      68             :                           : nullptr),
      69       15172 :     _porosity(getMaterialProperty<Real>("PorousFlow_porosity_nodal")),
      70       15172 :     _porosity_old(getMaterialPropertyOld<Real>("PorousFlow_porosity_nodal")),
      71       15172 :     _dporosity_dvar(getMaterialProperty<std::vector<Real>>("dPorousFlow_porosity_nodal_dvar")),
      72        7586 :     _dporosity_dgradvar(
      73        7586 :         getMaterialProperty<std::vector<RealGradient>>("dPorousFlow_porosity_nodal_dgradvar")),
      74       15172 :     _nearest_qp(_strain_at_nearest_qp
      75        7641 :                     ? &getMaterialProperty<unsigned int>("PorousFlow_nearestqp_nodal")
      76             :                     : nullptr),
      77       14985 :     _fluid_density(_multiply_by_density ? &getMaterialProperty<std::vector<Real>>(
      78             :                                               "PorousFlow_fluid_phase_density_nodal")
      79             :                                         : nullptr),
      80       14985 :     _fluid_density_old(_multiply_by_density ? &getMaterialPropertyOld<std::vector<Real>>(
      81             :                                                   "PorousFlow_fluid_phase_density_nodal")
      82             :                                             : nullptr),
      83       15172 :     _dfluid_density_dvar(_multiply_by_density
      84       14985 :                              ? &getMaterialProperty<std::vector<std::vector<Real>>>(
      85             :                                    "dPorousFlow_fluid_phase_density_nodal_dvar")
      86             :                              : nullptr),
      87       15172 :     _fluid_saturation_nodal(getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_nodal")),
      88        7586 :     _fluid_saturation_nodal_old(
      89        7586 :         getMaterialPropertyOld<std::vector<Real>>("PorousFlow_saturation_nodal")),
      90        7586 :     _dfluid_saturation_nodal_dvar(
      91        7586 :         getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_saturation_nodal_dvar")),
      92       15172 :     _mass_frac(getMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal")),
      93        7586 :     _mass_frac_old(
      94        7586 :         getMaterialPropertyOld<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal")),
      95       15172 :     _dmass_frac_dvar(getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
      96        7586 :         "dPorousFlow_mass_frac_nodal_dvar"))
      97             : {
      98        7586 :   if (_fluid_component >= _dictator.numComponents())
      99           0 :     paramError(
     100             :         "fluid_component",
     101             :         "The Dictator proclaims that the maximum fluid component index in this simulation is ",
     102           0 :         _dictator.numComponents() - 1,
     103             :         " whereas you have used ",
     104           0 :         _fluid_component,
     105             :         ". Remember that indexing starts at 0. The Dictator does not take such mistakes lightly.");
     106        7586 : }
     107             : 
     108             : Real
     109   133015892 : PorousFlowMassTimeDerivative::computeQpResidual()
     110             : {
     111             :   Real mass = 0.0;
     112             :   Real mass_old = 0.0;
     113   274723224 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
     114             :   {
     115   141707332 :     const Real dens = (_multiply_by_density ? (*_fluid_density)[_i][ph] : 1.0);
     116   141707332 :     mass += dens * _fluid_saturation_nodal[_i][ph] * _mass_frac[_i][ph][_fluid_component];
     117   141707332 :     const Real dens_old = (_multiply_by_density ? (*_fluid_density_old)[_i][ph] : 1.0);
     118   141707332 :     mass_old +=
     119   141707332 :         dens_old * _fluid_saturation_nodal_old[_i][ph] * _mass_frac_old[_i][ph][_fluid_component];
     120             :   }
     121   133015892 :   const Real strain = (_has_total_strain ? (*_total_strain_old)[_qp].trace() : 0.0);
     122             : 
     123   133015892 :   return _test[_i][_qp] * (1.0 + strain) * (_porosity[_i] * mass - _porosity_old[_i] * mass_old) /
     124   133015892 :          _dt;
     125             : }
     126             : 
     127             : Real
     128   600316880 : PorousFlowMassTimeDerivative::computeQpJacobian()
     129             : {
     130             :   // If the variable is not a PorousFlow variable (very unusual), the diag Jacobian terms are 0
     131   600316880 :   if (!_var_is_porflow_var)
     132             :     return 0.0;
     133   600316272 :   return computeQpJac(_dictator.porousFlowVariableNum(_var.number()));
     134             : }
     135             : 
     136             : Real
     137   417909560 : PorousFlowMassTimeDerivative::computeQpOffDiagJacobian(unsigned int jvar)
     138             : {
     139             :   // If the variable is not a PorousFlow variable, the OffDiag Jacobian terms are 0
     140   417909560 :   if (_dictator.notPorousFlowVariable(jvar))
     141             :     return 0.0;
     142   417908344 :   return computeQpJac(_dictator.porousFlowVariableNum(jvar));
     143             : }
     144             : 
     145             : Real
     146  1018224616 : PorousFlowMassTimeDerivative::computeQpJac(unsigned int pvar)
     147             : {
     148  1018224616 :   const unsigned nearest_qp = (_strain_at_nearest_qp ? (*_nearest_qp)[_i] : _i);
     149             : 
     150  1018224616 :   const Real strain = (_has_total_strain ? (*_total_strain_old)[_qp].trace() : 0.0);
     151             : 
     152             :   // porosity is dependent on variables that are lumped to the nodes,
     153             :   // but it can depend on the gradient
     154             :   // of variables, which are NOT lumped to the nodes, hence:
     155             :   Real dmass = 0.0;
     156  2080545720 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
     157             :   {
     158  1062321104 :     const Real dens = (_multiply_by_density ? (*_fluid_density)[_i][ph] : 1.0);
     159  1062321104 :     dmass += dens * _fluid_saturation_nodal[_i][ph] * _mass_frac[_i][ph][_fluid_component] *
     160  1062321104 :              _dporosity_dgradvar[_i][pvar] * _grad_phi[_j][nearest_qp];
     161             :   }
     162             : 
     163  1018224616 :   if (_i != _j)
     164   876945332 :     return _test[_i][_qp] * (1.0 + strain) * dmass / _dt;
     165             : 
     166             :   // As the fluid mass is lumped to the nodes, only non-zero terms are for _i==_j
     167   295187436 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
     168             :   {
     169   153908152 :     if (_multiply_by_density)
     170   153859112 :       dmass += (*_dfluid_density_dvar)[_i][ph][pvar] * _fluid_saturation_nodal[_i][ph] *
     171   153859112 :                _mass_frac[_i][ph][_fluid_component] * _porosity[_i];
     172   153908152 :     const Real dens = (_multiply_by_density ? (*_fluid_density)[_i][ph] : 1.0);
     173   153908152 :     dmass += dens * _dfluid_saturation_nodal_dvar[_i][ph][pvar] *
     174   153908152 :              _mass_frac[_i][ph][_fluid_component] * _porosity[_i];
     175   153908152 :     dmass += dens * _fluid_saturation_nodal[_i][ph] *
     176   153908152 :              _dmass_frac_dvar[_i][ph][_fluid_component][pvar] * _porosity[_i];
     177   153908152 :     dmass += dens * _fluid_saturation_nodal[_i][ph] * _mass_frac[_i][ph][_fluid_component] *
     178   153908152 :              _dporosity_dvar[_i][pvar];
     179             :   }
     180   141279284 :   return _test[_i][_qp] * (1.0 + strain) * dmass / _dt;
     181             : }