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             : 
      16             : InputParameters
      17        2046 : PorousFlowEnergyTimeDerivative::validParams()
      18             : {
      19        2046 :   InputParameters params = TimeKernel::validParams();
      20        4092 :   params.addParam<bool>("strain_at_nearest_qp",
      21        4092 :                         false,
      22             :                         "When calculating nodal porosity that depends on strain, use the strain at "
      23             :                         "the nearest quadpoint.  This adds a small extra computational burden, and "
      24             :                         "is not necessary for simulations involving only linear lagrange elements. "
      25             :                         " If you set this to true, you will also want to set the same parameter to "
      26             :                         "true for related Kernels and Materials");
      27        4092 :   params.addParam<std::string>(
      28             :       "base_name",
      29             :       "For mechanically-coupled systems, this Kernel will depend on the volumetric strain.  "
      30             :       "base_name should almost always be the same base_name as given to the TensorMechanics object "
      31             :       "that computes strain.  Supplying a base_name to this Kernel but not defining an associated "
      32             :       "TensorMechanics strain calculator means that this Kernel will not depend on volumetric "
      33             :       "strain.  That could be useful when models contain solid mechanics that is not coupled to "
      34             :       "porous flow, for example");
      35        4092 :   params.addRequiredParam<UserObjectName>(
      36             :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names.");
      37        2046 :   params.set<bool>("use_displaced_mesh") = false;
      38        2046 :   params.suppressParameter<bool>("use_displaced_mesh");
      39        2046 :   params.addClassDescription("Derivative of heat-energy-density wrt time");
      40        2046 :   return params;
      41           0 : }
      42             : 
      43        1120 : PorousFlowEnergyTimeDerivative::PorousFlowEnergyTimeDerivative(const InputParameters & parameters)
      44             :   : TimeKernel(parameters),
      45        1120 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      46        1120 :     _var_is_porflow_var(_dictator.isPorousFlowVariable(_var.number())),
      47        1120 :     _num_phases(_dictator.numPhases()),
      48        1120 :     _fluid_present(_num_phases > 0),
      49        2240 :     _strain_at_nearest_qp(getParam<bool>("strain_at_nearest_qp")),
      50        2372 :     _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
      51        1120 :     _has_total_strain(hasMaterialProperty<RankTwoTensor>(_base_name + "total_strain")),
      52        2240 :     _total_strain_old(_has_total_strain
      53        1385 :                           ? &getMaterialPropertyOld<RankTwoTensor>(_base_name + "total_strain")
      54             :                           : nullptr),
      55        2240 :     _porosity(getMaterialProperty<Real>("PorousFlow_porosity_nodal")),
      56        2240 :     _porosity_old(getMaterialPropertyOld<Real>("PorousFlow_porosity_nodal")),
      57        2240 :     _dporosity_dvar(getMaterialProperty<std::vector<Real>>("dPorousFlow_porosity_nodal_dvar")),
      58        1120 :     _dporosity_dgradvar(
      59        1120 :         getMaterialProperty<std::vector<RealGradient>>("dPorousFlow_porosity_nodal_dgradvar")),
      60        2240 :     _nearest_qp(_strain_at_nearest_qp
      61        1131 :                     ? &getMaterialProperty<unsigned int>("PorousFlow_nearestqp_nodal")
      62             :                     : nullptr),
      63        2240 :     _rock_energy_nodal(getMaterialProperty<Real>("PorousFlow_matrix_internal_energy_nodal")),
      64        2240 :     _rock_energy_nodal_old(getMaterialPropertyOld<Real>("PorousFlow_matrix_internal_energy_nodal")),
      65        1120 :     _drock_energy_nodal_dvar(
      66        1120 :         getMaterialProperty<std::vector<Real>>("dPorousFlow_matrix_internal_energy_nodal_dvar")),
      67        2141 :     _fluid_density(_fluid_present ? &getMaterialProperty<std::vector<Real>>(
      68             :                                         "PorousFlow_fluid_phase_density_nodal")
      69             :                                   : nullptr),
      70        2141 :     _fluid_density_old(_fluid_present ? &getMaterialPropertyOld<std::vector<Real>>(
      71             :                                             "PorousFlow_fluid_phase_density_nodal")
      72             :                                       : nullptr),
      73        2141 :     _dfluid_density_dvar(_fluid_present ? &getMaterialProperty<std::vector<std::vector<Real>>>(
      74             :                                               "dPorousFlow_fluid_phase_density_nodal_dvar")
      75             :                                         : nullptr),
      76        1120 :     _fluid_saturation_nodal(
      77        2141 :         _fluid_present ? &getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_nodal")
      78             :                        : nullptr),
      79        1120 :     _fluid_saturation_nodal_old(
      80        2141 :         _fluid_present ? &getMaterialPropertyOld<std::vector<Real>>("PorousFlow_saturation_nodal")
      81             :                        : nullptr),
      82        2240 :     _dfluid_saturation_nodal_dvar(_fluid_present
      83        2141 :                                       ? &getMaterialProperty<std::vector<std::vector<Real>>>(
      84             :                                             "dPorousFlow_saturation_nodal_dvar")
      85             :                                       : nullptr),
      86        2141 :     _energy_nodal(_fluid_present ? &getMaterialProperty<std::vector<Real>>(
      87             :                                        "PorousFlow_fluid_phase_internal_energy_nodal")
      88             :                                  : nullptr),
      89        2141 :     _energy_nodal_old(_fluid_present ? &getMaterialPropertyOld<std::vector<Real>>(
      90             :                                            "PorousFlow_fluid_phase_internal_energy_nodal")
      91             :                                      : nullptr),
      92        2141 :     _denergy_nodal_dvar(_fluid_present ? &getMaterialProperty<std::vector<std::vector<Real>>>(
      93             :                                              "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")
      94        1120 :                                        : nullptr)
      95             : {
      96        1120 : }
      97             : 
      98             : Real
      99     8380204 : PorousFlowEnergyTimeDerivative::computeQpResidual()
     100             : {
     101             :   /// Porous matrix heat energy
     102     8380204 :   Real energy = (1.0 - _porosity[_i]) * _rock_energy_nodal[_i];
     103     8380204 :   Real energy_old = (1.0 - _porosity_old[_i]) * _rock_energy_nodal_old[_i];
     104             : 
     105             :   /// Add the fluid heat energy
     106     8380204 :   if (_fluid_present)
     107    17212004 :     for (unsigned ph = 0; ph < _num_phases; ++ph)
     108             :     {
     109     8872168 :       energy += (*_fluid_density)[_i][ph] * (*_fluid_saturation_nodal)[_i][ph] *
     110     8872168 :                 (*_energy_nodal)[_i][ph] * _porosity[_i];
     111     8872168 :       energy_old += (*_fluid_density_old)[_i][ph] * (*_fluid_saturation_nodal_old)[_i][ph] *
     112     8872168 :                     (*_energy_nodal_old)[_i][ph] * _porosity_old[_i];
     113             :     }
     114     8380204 :   const Real strain = (_has_total_strain ? (*_total_strain_old)[_qp].trace() : 0.0);
     115             : 
     116     8380204 :   return _test[_i][_qp] * (1.0 + strain) * (energy - energy_old) / _dt;
     117             : }
     118             : 
     119             : Real
     120    27083880 : PorousFlowEnergyTimeDerivative::computeQpJacobian()
     121             : {
     122             :   // If the variable is not a PorousFlow variable (very unusual), the diag Jacobian terms are 0
     123    27083880 :   if (!_var_is_porflow_var)
     124             :     return 0.0;
     125    27083880 :   return computeQpJac(_dictator.porousFlowVariableNum(_var.number()));
     126             : }
     127             : 
     128             : Real
     129    29953960 : PorousFlowEnergyTimeDerivative::computeQpOffDiagJacobian(unsigned int jvar)
     130             : {
     131             :   // If the variable is not a PorousFlow variable, the OffDiag Jacobian terms are 0
     132    29953960 :   if (_dictator.notPorousFlowVariable(jvar))
     133             :     return 0.0;
     134    29953960 :   return computeQpJac(_dictator.porousFlowVariableNum(jvar));
     135             : }
     136             : 
     137             : Real
     138    57037840 : PorousFlowEnergyTimeDerivative::computeQpJac(unsigned int pvar) const
     139             : {
     140    57037840 :   const unsigned nearest_qp = (_strain_at_nearest_qp ? (*_nearest_qp)[_i] : _i);
     141             : 
     142    57037840 :   const Real strain = (_has_total_strain ? (*_total_strain_old)[_qp].trace() : 0.0);
     143             : 
     144             :   // porosity is dependent on variables that are lumped to the nodes,
     145             :   // but it can depend on the gradient
     146             :   // of variables, which are NOT lumped to the nodes, hence:
     147    57037840 :   Real denergy = -_dporosity_dgradvar[_i][pvar] * _grad_phi[_j][_i] * _rock_energy_nodal[_i];
     148   115767392 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
     149    58729552 :     denergy += (*_fluid_density)[_i][ph] * (*_fluid_saturation_nodal)[_i][ph] *
     150    58729552 :                (*_energy_nodal)[_i][ph] * _dporosity_dgradvar[_i][pvar] * _grad_phi[_j][nearest_qp];
     151             : 
     152    57037840 :   if (_i != _j)
     153    45189840 :     return _test[_i][_qp] * (1.0 + strain) * denergy / _dt;
     154             : 
     155             :   // As the fluid energy is lumped to the nodes, only non-zero terms are for _i==_j
     156    11848000 :   denergy += -_dporosity_dvar[_i][pvar] * _rock_energy_nodal[_i];
     157    11848000 :   denergy += (1.0 - _porosity[_i]) * _drock_energy_nodal_dvar[_i][pvar];
     158    24418624 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
     159             :   {
     160    12570624 :     denergy += (*_dfluid_density_dvar)[_i][ph][pvar] * (*_fluid_saturation_nodal)[_i][ph] *
     161    12570624 :                (*_energy_nodal)[_i][ph] * _porosity[_i];
     162    12570624 :     denergy += (*_fluid_density)[_i][ph] * (*_dfluid_saturation_nodal_dvar)[_i][ph][pvar] *
     163    12570624 :                (*_energy_nodal)[_i][ph] * _porosity[_i];
     164    12570624 :     denergy += (*_fluid_density)[_i][ph] * (*_fluid_saturation_nodal)[_i][ph] *
     165    12570624 :                (*_denergy_nodal_dvar)[_i][ph][pvar] * _porosity[_i];
     166    12570624 :     denergy += (*_fluid_density)[_i][ph] * (*_fluid_saturation_nodal)[_i][ph] *
     167    12570624 :                (*_energy_nodal)[_i][ph] * _dporosity_dvar[_i][pvar];
     168             :   }
     169    11848000 :   return _test[_i][_qp] * (1.0 + strain) * denergy / _dt;
     170             : }