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 "PorousFlowMassRadioactiveDecay.h"
      11             : 
      12             : #include "MooseVariable.h"
      13             : 
      14             : #include "libmesh/quadrature.h"
      15             : 
      16             : #include <limits>
      17             : 
      18             : registerMooseObject("PorousFlowApp", PorousFlowMassRadioactiveDecay);
      19             : 
      20             : InputParameters
      21          41 : PorousFlowMassRadioactiveDecay::validParams()
      22             : {
      23          41 :   InputParameters params = TimeKernel::validParams();
      24          82 :   params.addParam<bool>("strain_at_nearest_qp",
      25          82 :                         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          82 :   params.addParam<unsigned int>(
      32          82 :       "fluid_component", 0, "The index corresponding to the fluid component for this kernel");
      33          82 :   params.addRequiredParam<Real>("decay_rate",
      34             :                                 "The decay rate (units 1/time) for the fluid component");
      35          82 :   params.addRequiredParam<UserObjectName>(
      36             :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names.");
      37          41 :   params.addClassDescription("Radioactive decay of a fluid component");
      38          41 :   return params;
      39           0 : }
      40             : 
      41          22 : PorousFlowMassRadioactiveDecay::PorousFlowMassRadioactiveDecay(const InputParameters & parameters)
      42             :   : TimeKernel(parameters),
      43          22 :     _decay_rate(getParam<Real>("decay_rate")),
      44          44 :     _fluid_component(getParam<unsigned int>("fluid_component")),
      45          22 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      46          22 :     _var_is_porflow_var(_dictator.isPorousFlowVariable(_var.number())),
      47          22 :     _num_phases(_dictator.numPhases()),
      48          44 :     _strain_at_nearest_qp(getParam<bool>("strain_at_nearest_qp")),
      49          44 :     _porosity(getMaterialProperty<Real>("PorousFlow_porosity_nodal")),
      50          44 :     _dporosity_dvar(getMaterialProperty<std::vector<Real>>("dPorousFlow_porosity_nodal_dvar")),
      51          22 :     _dporosity_dgradvar(
      52          22 :         getMaterialProperty<std::vector<RealGradient>>("dPorousFlow_porosity_nodal_dgradvar")),
      53          44 :     _nearest_qp(_strain_at_nearest_qp
      54          22 :                     ? &getMaterialProperty<unsigned int>("PorousFlow_nearestqp_nodal")
      55             :                     : nullptr),
      56          44 :     _fluid_density(getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal")),
      57          44 :     _dfluid_density_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
      58             :         "dPorousFlow_fluid_phase_density_nodal_dvar")),
      59          44 :     _fluid_saturation_nodal(getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_nodal")),
      60          22 :     _dfluid_saturation_nodal_dvar(
      61          22 :         getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_saturation_nodal_dvar")),
      62          44 :     _mass_frac(getMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal")),
      63          44 :     _dmass_frac_dvar(getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
      64          22 :         "dPorousFlow_mass_frac_nodal_dvar"))
      65             : {
      66          22 :   if (_fluid_component >= _dictator.numComponents())
      67           0 :     paramError(
      68             :         "fluid_component",
      69             :         "The Dictator proclaims that the maximum fluid component index in this simulation is ",
      70           0 :         _dictator.numComponents() - 1,
      71             :         " whereas you have used ",
      72           0 :         _fluid_component,
      73             :         ". Remember that indexing starts at 0. The Dictator does not take such mistakes lightly.");
      74          22 : }
      75             : 
      76             : Real
      77         720 : PorousFlowMassRadioactiveDecay::computeQpResidual()
      78             : {
      79             :   Real mass = 0.0;
      80        1440 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
      81         720 :     mass += _fluid_density[_i][ph] * _fluid_saturation_nodal[_i][ph] *
      82         720 :             _mass_frac[_i][ph][_fluid_component];
      83             : 
      84         720 :   return _test[_i][_qp] * _decay_rate * _porosity[_i] * mass;
      85             : }
      86             : 
      87             : Real
      88        1200 : PorousFlowMassRadioactiveDecay::computeQpJacobian()
      89             : {
      90             :   // If the variable is not a PorousFlow variable (very unusual), the diag Jacobian terms are 0
      91        1200 :   if (!_var_is_porflow_var)
      92             :     return 0.0;
      93        1200 :   return computeQpJac(_dictator.porousFlowVariableNum(_var.number()));
      94             : }
      95             : 
      96             : Real
      97           0 : PorousFlowMassRadioactiveDecay::computeQpOffDiagJacobian(unsigned int jvar)
      98             : {
      99             :   // If the variable is not a PorousFlow variable, the OffDiag Jacobian terms are 0
     100           0 :   if (_dictator.notPorousFlowVariable(jvar))
     101             :     return 0.0;
     102           0 :   return computeQpJac(_dictator.porousFlowVariableNum(jvar));
     103             : }
     104             : 
     105             : Real
     106        1200 : PorousFlowMassRadioactiveDecay::computeQpJac(unsigned int pvar)
     107             : {
     108        1200 :   const unsigned nearest_qp = (_strain_at_nearest_qp ? (*_nearest_qp)[_i] : _i);
     109             : 
     110             :   // porosity is dependent on variables that are lumped to the nodes,
     111             :   // but it can depend on the gradient
     112             :   // of variables, which are NOT lumped to the nodes, hence:
     113             :   Real dmass = 0.0;
     114        2400 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
     115        1200 :     dmass += _fluid_density[_i][ph] * _fluid_saturation_nodal[_i][ph] *
     116        1200 :              _mass_frac[_i][ph][_fluid_component] * _dporosity_dgradvar[_i][pvar] *
     117        1200 :              _grad_phi[_j][nearest_qp];
     118             : 
     119        1200 :   if (_i != _j)
     120         600 :     return _test[_i][_qp] * _decay_rate * dmass;
     121             : 
     122             :   // As the fluid mass is lumped to the nodes, only non-zero terms are for _i==_j
     123        1200 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
     124             :   {
     125         600 :     dmass += _dfluid_density_dvar[_i][ph][pvar] * _fluid_saturation_nodal[_i][ph] *
     126         600 :              _mass_frac[_i][ph][_fluid_component] * _porosity[_i];
     127         600 :     dmass += _fluid_density[_i][ph] * _dfluid_saturation_nodal_dvar[_i][ph][pvar] *
     128         600 :              _mass_frac[_i][ph][_fluid_component] * _porosity[_i];
     129         600 :     dmass += _fluid_density[_i][ph] * _fluid_saturation_nodal[_i][ph] *
     130         600 :              _dmass_frac_dvar[_i][ph][_fluid_component][pvar] * _porosity[_i];
     131         600 :     dmass += _fluid_density[_i][ph] * _fluid_saturation_nodal[_i][ph] *
     132         600 :              _mass_frac[_i][ph][_fluid_component] * _dporosity_dvar[_i][pvar];
     133             :   }
     134         600 :   return _test[_i][_qp] * _decay_rate * dmass;
     135             : }