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 "PorousFlowMassVolumetricExpansion.h"
      11             : 
      12             : #include "MooseVariable.h"
      13             : 
      14             : registerMooseObject("PorousFlowApp", PorousFlowMassVolumetricExpansion);
      15             : 
      16             : InputParameters
      17        1248 : PorousFlowMassVolumetricExpansion::validParams()
      18             : {
      19        1248 :   InputParameters params = TimeKernel::validParams();
      20        2496 :   params.addParam<bool>("strain_at_nearest_qp",
      21        2496 :                         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        2496 :   params.addParam<bool>(
      28             :       "multiply_by_density",
      29        2496 :       true,
      30             :       "If true, then this Kernel represents component_mass*rate_of_solid_volumetric_expansion.  If "
      31             :       "flase, then this Kernel represents component_volume*rate_of_solid_volumetric_expansion "
      32             :       "(care must then be taken when using other PorousFlow objects, such as the "
      33             :       "PorousFlowFluidMass postprocessor).");
      34        2496 :   params.addParam<unsigned int>(
      35        2496 :       "fluid_component", 0, "The index corresponding to the component for this kernel");
      36        2496 :   params.addRequiredParam<UserObjectName>(
      37             :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names.");
      38        1248 :   params.set<bool>("use_displaced_mesh") = false;
      39        1248 :   params.suppressParameter<bool>("use_displaced_mesh");
      40        1248 :   params.addClassDescription("Component_mass*rate_of_solid_volumetric_expansion.  This Kernel "
      41             :                              "lumps the component mass to the nodes.");
      42        1248 :   return params;
      43           0 : }
      44             : 
      45         673 : PorousFlowMassVolumetricExpansion::PorousFlowMassVolumetricExpansion(
      46         673 :     const InputParameters & parameters)
      47             :   : TimeKernel(parameters),
      48         673 :     _fluid_component(getParam<unsigned int>("fluid_component")),
      49         673 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      50         673 :     _var_is_porflow_var(!_dictator.notPorousFlowVariable(_var.number())),
      51         673 :     _num_phases(_dictator.numPhases()),
      52        1346 :     _strain_at_nearest_qp(getParam<bool>("strain_at_nearest_qp")),
      53        1346 :     _multiply_by_density(getParam<bool>("multiply_by_density")),
      54        1346 :     _porosity(getMaterialProperty<Real>("PorousFlow_porosity_nodal")),
      55        1346 :     _dporosity_dvar(getMaterialProperty<std::vector<Real>>("dPorousFlow_porosity_nodal_dvar")),
      56         673 :     _dporosity_dgradvar(
      57         673 :         getMaterialProperty<std::vector<RealGradient>>("dPorousFlow_porosity_nodal_dgradvar")),
      58        1346 :     _nearest_qp(_strain_at_nearest_qp
      59         706 :                     ? &getMaterialProperty<unsigned int>("PorousFlow_nearestqp_nodal")
      60             :                     : nullptr),
      61        1335 :     _fluid_density(_multiply_by_density ? &getMaterialProperty<std::vector<Real>>(
      62             :                                               "PorousFlow_fluid_phase_density_nodal")
      63             :                                         : nullptr),
      64        1346 :     _dfluid_density_dvar(_multiply_by_density
      65        1335 :                              ? &getMaterialProperty<std::vector<std::vector<Real>>>(
      66             :                                    "dPorousFlow_fluid_phase_density_nodal_dvar")
      67             :                              : nullptr),
      68        1346 :     _fluid_saturation(getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_nodal")),
      69         673 :     _dfluid_saturation_dvar(
      70         673 :         getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_saturation_nodal_dvar")),
      71        1346 :     _mass_frac(getMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal")),
      72        1346 :     _dmass_frac_dvar(getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
      73             :         "dPorousFlow_mass_frac_nodal_dvar")),
      74        1346 :     _strain_rate_qp(getMaterialProperty<Real>("PorousFlow_volumetric_strain_rate_qp")),
      75        1346 :     _dstrain_rate_qp_dvar(getMaterialProperty<std::vector<RealGradient>>(
      76         673 :         "dPorousFlow_volumetric_strain_rate_qp_dvar"))
      77             : {
      78         673 :   if (_fluid_component >= _dictator.numComponents())
      79           0 :     mooseError("The Dictator proclaims that the number of components in this simulation is ",
      80           0 :                _dictator.numComponents(),
      81             :                " whereas you have used the Kernel PorousFlowComponetMassVolumetricExpansion with "
      82             :                "component = ",
      83           0 :                _fluid_component,
      84             :                ".  The Dictator is watching you");
      85         673 : }
      86             : 
      87             : Real
      88    68788640 : PorousFlowMassVolumetricExpansion::computeQpResidual()
      89             : {
      90             :   Real mass = 0.0;
      91   137612800 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
      92             :   {
      93    68824160 :     const Real dens = (_multiply_by_density ? (*_fluid_density)[_i][ph] : 1.0);
      94    68824160 :     mass += dens * _fluid_saturation[_i][ph] * _mass_frac[_i][ph][_fluid_component];
      95             :   }
      96             : 
      97    68788640 :   return _test[_i][_qp] * mass * _porosity[_i] * _strain_rate_qp[_qp];
      98             : }
      99             : 
     100             : Real
     101   362304064 : PorousFlowMassVolumetricExpansion::computeQpJacobian()
     102             : {
     103   362304064 :   return computedMassQpJac(_var.number()) + computedVolQpJac(_var.number());
     104             : }
     105             : 
     106             : Real
     107   332065536 : PorousFlowMassVolumetricExpansion::computeQpOffDiagJacobian(unsigned int jvar)
     108             : {
     109   332065536 :   return computedMassQpJac(jvar) + computedVolQpJac(jvar);
     110             : }
     111             : 
     112             : Real
     113   694369600 : PorousFlowMassVolumetricExpansion::computedVolQpJac(unsigned int jvar) const
     114             : {
     115   694369600 :   if (_dictator.notPorousFlowVariable(jvar))
     116             :     return 0.0;
     117             : 
     118   694369600 :   const unsigned int pvar = _dictator.porousFlowVariableNum(jvar);
     119             : 
     120             :   Real mass = 0.0;
     121  1389200000 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
     122             :   {
     123   694830400 :     const Real dens = (_multiply_by_density ? (*_fluid_density)[_i][ph] : 1.0);
     124   694830400 :     mass += dens * _fluid_saturation[_i][ph] * _mass_frac[_i][ph][_fluid_component];
     125             :   }
     126             : 
     127   694369600 :   Real dvol = _dstrain_rate_qp_dvar[_qp][pvar] * _grad_phi[_j][_qp];
     128             : 
     129   694369600 :   return _test[_i][_qp] * mass * _porosity[_i] * dvol;
     130             : }
     131             : Real
     132   694369600 : PorousFlowMassVolumetricExpansion::computedMassQpJac(unsigned int jvar) const
     133             : {
     134   694369600 :   if (_dictator.notPorousFlowVariable(jvar))
     135             :     return 0.0;
     136             : 
     137   694369600 :   const unsigned int pvar = _dictator.porousFlowVariableNum(jvar);
     138   694369600 :   const unsigned nearest_qp = (_strain_at_nearest_qp ? (*_nearest_qp)[_i] : _i);
     139             : 
     140             :   Real dmass = 0.0;
     141  1389200000 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
     142             :   {
     143   694830400 :     const Real dens = (_multiply_by_density ? (*_fluid_density)[_i][ph] : 1.0);
     144   694830400 :     dmass += dens * _fluid_saturation[_i][ph] * _mass_frac[_i][ph][_fluid_component] *
     145   694830400 :              _dporosity_dgradvar[_i][pvar] * _grad_phi[_j][nearest_qp];
     146             :   }
     147             : 
     148   694369600 :   if (_i != _j)
     149   619931352 :     return _test[_i][_qp] * dmass * _strain_rate_qp[_qp];
     150             : 
     151   148991696 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
     152             :   {
     153    74553448 :     if (_multiply_by_density)
     154    74552168 :       dmass += (*_dfluid_density_dvar)[_i][ph][pvar] * _fluid_saturation[_i][ph] *
     155    74552168 :                _mass_frac[_i][ph][_fluid_component] * _porosity[_i];
     156    74553448 :     const Real dens = (_multiply_by_density ? (*_fluid_density)[_i][ph] : 1.0);
     157    74553448 :     dmass += dens * _dfluid_saturation_dvar[_i][ph][pvar] * _mass_frac[_i][ph][_fluid_component] *
     158    74553448 :              _porosity[_i];
     159    74553448 :     dmass += dens * _fluid_saturation[_i][ph] * _dmass_frac_dvar[_i][ph][_fluid_component][pvar] *
     160             :              _porosity[_i];
     161    74553448 :     dmass += dens * _fluid_saturation[_i][ph] * _mass_frac[_i][ph][_fluid_component] *
     162    74553448 :              _dporosity_dvar[_i][pvar];
     163             :   }
     164             : 
     165    74438248 :   return _test[_i][_qp] * dmass * _strain_rate_qp[_qp];
     166             : }