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             : registerMooseObject("PorousFlowApp", ADPorousFlowMassRadioactiveDecay);
      20             : 
      21             : template <bool is_ad>
      22             : InputParameters
      23          32 : PorousFlowMassRadioactiveDecayTempl<is_ad>::validParams()
      24             : {
      25             :   InputParameters params = GenericKernel<is_ad>::validParams();
      26          64 :   params.set<MultiMooseEnum>("vector_tags") = "time";
      27          64 :   params.set<MultiMooseEnum>("matrix_tags") = "system time";
      28          64 :   params.addParam<bool>("strain_at_nearest_qp",
      29          64 :                         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          64 :   params.addParam<unsigned int>(
      36          64 :       "fluid_component", 0, "The index corresponding to the fluid component for this kernel");
      37          64 :   params.addRequiredParam<Real>("decay_rate",
      38             :                                 "The decay rate (units 1/time) for the fluid component");
      39          64 :   params.addRequiredParam<UserObjectName>(
      40             :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names.");
      41          32 :   params.addClassDescription("Radioactive decay of a fluid component");
      42          32 :   return params;
      43           0 : }
      44             : 
      45             : template <bool is_ad>
      46          17 : PorousFlowMassRadioactiveDecayTempl<is_ad>::PorousFlowMassRadioactiveDecayTempl(
      47             :     const InputParameters & parameters)
      48             :   : PorousFlowLumpedKernelBaseTempl<is_ad>(parameters),
      49          17 :     _decay_rate(this->template getParam<Real>("decay_rate")),
      50          34 :     _fluid_component(this->template getParam<unsigned int>("fluid_component")),
      51          17 :     _dictator(this->template getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      52          17 :     _var_is_porflow_var(_dictator.isPorousFlowVariable(_var.number())),
      53          17 :     _num_phases(_dictator.numPhases()),
      54          34 :     _strain_at_nearest_qp(this->template getParam<bool>("strain_at_nearest_qp")),
      55          34 :     _porosity(this->template getGenericMaterialProperty<Real, is_ad>("PorousFlow_porosity_nodal")),
      56          17 :     _dporosity_dvar(is_ad ? nullptr
      57          10 :                           : &this->template getMaterialProperty<std::vector<Real>>(
      58             :                                 "dPorousFlow_porosity_nodal_dvar")),
      59          17 :     _dporosity_dgradvar(is_ad ? nullptr
      60          10 :                               : &this->template getMaterialProperty<std::vector<RealGradient>>(
      61             :                                     "dPorousFlow_porosity_nodal_dgradvar")),
      62          17 :     _nearest_qp(_strain_at_nearest_qp ? &this->template getMaterialProperty<unsigned int>(
      63             :                                             "PorousFlow_nearestqp_nodal")
      64             :                                       : nullptr),
      65          34 :     _fluid_density(this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
      66             :         "PorousFlow_fluid_phase_density_nodal")),
      67          17 :     _dfluid_density_dvar(is_ad
      68             :                              ? nullptr
      69          10 :                              : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
      70             :                                    "dPorousFlow_fluid_phase_density_nodal_dvar")),
      71          34 :     _fluid_saturation_nodal(this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
      72             :         "PorousFlow_saturation_nodal")),
      73          17 :     _dfluid_saturation_nodal_dvar(
      74             :         is_ad ? nullptr
      75          10 :               : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
      76             :                     "dPorousFlow_saturation_nodal_dvar")),
      77          34 :     _mass_frac(this->template getGenericMaterialProperty<std::vector<std::vector<Real>>, is_ad>(
      78             :         "PorousFlow_mass_frac_nodal")),
      79          17 :     _dmass_frac_dvar(
      80             :         is_ad ? nullptr
      81          10 :               : &this->template getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
      82          17 :                     "dPorousFlow_mass_frac_nodal_dvar"))
      83             : {
      84          17 :   if (_fluid_component >= _dictator.numComponents())
      85           0 :     this->paramError(
      86             :         "fluid_component",
      87             :         "The Dictator proclaims that the maximum fluid component index in this simulation is ",
      88           0 :         _dictator.numComponents() - 1,
      89             :         " whereas you have used ",
      90           0 :         _fluid_component,
      91             :         ". Remember that indexing starts at 0. The Dictator does not take such mistakes lightly.");
      92          17 : }
      93             : 
      94             : template <bool is_ad>
      95             : GenericReal<is_ad>
      96        1912 : PorousFlowMassRadioactiveDecayTempl<is_ad>::computeQpResidual()
      97             : {
      98        1408 :   GenericReal<is_ad> mass = 0.0;
      99        3824 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
     100        3320 :     mass += _fluid_density[_i][ph] * _fluid_saturation_nodal[_i][ph] *
     101        1912 :             _mass_frac[_i][ph][_fluid_component];
     102             : 
     103        3320 :   return _test[_i][_qp] * _decay_rate * _porosity[_i] * mass;
     104             : }
     105             : 
     106             : template <bool is_ad>
     107             : Real
     108         840 : PorousFlowMassRadioactiveDecayTempl<is_ad>::computeQpJacobian()
     109             : {
     110             :   if constexpr (!is_ad)
     111             :   {
     112         840 :     if (!_var_is_porflow_var)
     113             :       return 0.0;
     114         840 :     return computeQpJac(_dictator.porousFlowVariableNum(_var.number()));
     115             :   }
     116           0 :   return 0.0;
     117             : }
     118             : 
     119             : template <bool is_ad>
     120             : Real
     121           0 : PorousFlowMassRadioactiveDecayTempl<is_ad>::computeQpOffDiagJacobian(unsigned int jvar)
     122             : {
     123             :   if constexpr (!is_ad)
     124             :   {
     125           0 :     if (_dictator.notPorousFlowVariable(jvar))
     126             :       return 0.0;
     127           0 :     return computeQpJac(_dictator.porousFlowVariableNum(jvar));
     128             :   }
     129             :   else
     130             :     libmesh_ignore(jvar);
     131           0 :   return 0.0;
     132             : }
     133             : 
     134             : template <bool is_ad>
     135             : Real
     136         840 : PorousFlowMassRadioactiveDecayTempl<is_ad>::computeQpJac(unsigned int pvar)
     137             : {
     138             :   if constexpr (!is_ad)
     139             :   {
     140         840 :     const unsigned nearest_qp = (_strain_at_nearest_qp ? (*_nearest_qp)[_i] : _i);
     141             : 
     142             :     // porosity can depend on grad(variables) evaluated at qps, so the grad-term is nonzero even
     143             :     // for off-node DOFs
     144             :     Real dmass = 0.0;
     145        1680 :     for (unsigned ph = 0; ph < _num_phases; ++ph)
     146         840 :       dmass += _fluid_density[_i][ph] * _fluid_saturation_nodal[_i][ph] *
     147         840 :                _mass_frac[_i][ph][_fluid_component] * (*_dporosity_dgradvar)[_i][pvar] *
     148         840 :                _grad_phi[_j][nearest_qp];
     149             : 
     150         840 :     if (_i != _j)
     151         420 :       return _test[_i][_qp] * _decay_rate * dmass;
     152             : 
     153             :     // As the fluid mass is lumped to the nodes, only non-zero terms are for _i==_j
     154         840 :     for (unsigned ph = 0; ph < _num_phases; ++ph)
     155             :     {
     156         420 :       dmass += (*_dfluid_density_dvar)[_i][ph][pvar] * _fluid_saturation_nodal[_i][ph] *
     157         420 :                _mass_frac[_i][ph][_fluid_component] * _porosity[_i];
     158         420 :       dmass += _fluid_density[_i][ph] * (*_dfluid_saturation_nodal_dvar)[_i][ph][pvar] *
     159         420 :                _mass_frac[_i][ph][_fluid_component] * _porosity[_i];
     160         420 :       dmass += _fluid_density[_i][ph] * _fluid_saturation_nodal[_i][ph] *
     161         420 :                (*_dmass_frac_dvar)[_i][ph][_fluid_component][pvar] * _porosity[_i];
     162         420 :       dmass += _fluid_density[_i][ph] * _fluid_saturation_nodal[_i][ph] *
     163         420 :                _mass_frac[_i][ph][_fluid_component] * (*_dporosity_dvar)[_i][pvar];
     164             :     }
     165         420 :     return _test[_i][_qp] * _decay_rate * dmass;
     166             :   }
     167             :   else
     168             :     libmesh_ignore(pvar);
     169           0 :   return 0.0;
     170             : }
     171             : 
     172             : template class PorousFlowMassRadioactiveDecayTempl<false>;
     173             : template class PorousFlowMassRadioactiveDecayTempl<true>;