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 "PorousFlowPreDis.h"
      11             : 
      12             : registerMooseObject("PorousFlowApp", PorousFlowPreDis);
      13             : registerMooseObject("PorousFlowApp", ADPorousFlowPreDis);
      14             : 
      15             : template <bool is_ad>
      16             : InputParameters
      17         670 : PorousFlowPreDisTempl<is_ad>::validParams()
      18             : {
      19             :   InputParameters params = GenericKernel<is_ad>::validParams();
      20        1340 :   params.set<MultiMooseEnum>("vector_tags") = "time";
      21        1340 :   params.set<MultiMooseEnum>("matrix_tags") = "system time";
      22        1340 :   params.addRequiredParam<std::vector<Real>>(
      23             :       "mineral_density",
      24             :       "Density (kg(precipitate)/m^3(precipitate)) of each secondary species in the "
      25             :       "aqueous precipitation-dissolution reaction system");
      26        1340 :   params.addRequiredParam<UserObjectName>(
      27             :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names.");
      28        1340 :   params.addRequiredParam<std::vector<Real>>("stoichiometry",
      29             :                                              "A vector of stoichiometric coefficients for the "
      30             :                                              "primary species that is the Variable of this Kernel: "
      31             :                                              "one for each precipitation-dissolution reaction "
      32             :                                              "(these are one columns of the 'reactions' matrix)");
      33         670 :   params.addClassDescription("Precipitation-dissolution of chemical species");
      34         670 :   return params;
      35           0 : }
      36             : 
      37             : template <bool is_ad>
      38         357 : PorousFlowPreDisTempl<is_ad>::PorousFlowPreDisTempl(const InputParameters & parameters)
      39             :   : PorousFlowLumpedKernelBaseTempl<is_ad>(parameters),
      40         357 :     _mineral_density(this->template getParam<std::vector<Real>>("mineral_density")),
      41         357 :     _dictator(this->template getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      42         357 :     _aq_ph(_dictator.aqueousPhaseNumber()),
      43         714 :     _porosity_old(this->template getMaterialPropertyOld<Real>("PorousFlow_porosity_nodal")),
      44         714 :     _saturation(this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
      45             :         "PorousFlow_saturation_nodal")),
      46         357 :     _dsaturation_dvar(is_ad ? nullptr
      47         357 :                             : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
      48             :                                   "dPorousFlow_saturation_nodal_dvar")),
      49         714 :     _reaction_rate(this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
      50             :         "PorousFlow_mineral_reaction_rate_nodal")),
      51         357 :     _dreaction_rate_dvar(is_ad
      52             :                              ? nullptr
      53         357 :                              : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
      54             :                                    "dPorousFlow_mineral_reaction_rate_nodal_dvar")),
      55        1428 :     _stoichiometry(this->template getParam<std::vector<Real>>("stoichiometry"))
      56             : {
      57             :   /* Not needed due to PorousFlow_mineral_reaction_rate already checking this condition
      58             :   if (_dictator.numPhases() < 1)
      59             :     mooseError("PorousFlowPreDis: The number of fluid phases must not be zero");
      60             :   */
      61             : 
      62         357 :   if (_mineral_density.size() != _dictator.numAqueousKinetic())
      63           2 :     this->paramError(
      64             :         "mineral_density",
      65             :         "The Dictator proclaims that the number of precipitation-dissolution secondary "
      66             :         "species in this simulation is ",
      67           2 :         _dictator.numAqueousKinetic(),
      68             :         " whereas you have provided ",
      69             :         _mineral_density.size(),
      70             :         ". The Dictator does not take such mistakes lightly");
      71             : 
      72         355 :   if (_stoichiometry.size() != _dictator.numAqueousKinetic())
      73           2 :     this->paramError(
      74             :         "stoichiometry",
      75             :         "The Dictator proclaims that the number of precipitation-dissolution secondary "
      76             :         "species in this simulation is ",
      77           2 :         _dictator.numAqueousKinetic(),
      78             :         " whereas you have provided ",
      79             :         _stoichiometry.size(),
      80             :         ". The Dictator does not take such mistakes lightly");
      81         353 : }
      82             : 
      83             : template <bool is_ad>
      84             : GenericReal<is_ad>
      85      930776 : PorousFlowPreDisTempl<is_ad>::computeQpResidual()
      86             : {
      87             :   /*
      88             :    *
      89             :    * Note the use of the OLD value of porosity here.
      90             :    * This strategy, which breaks the cyclic dependency between porosity
      91             :    * and mineral concentration, is used in
      92             :    * Kernel: PorousFlowPreDis
      93             :    * Material: PorousFlowPorosity
      94             :    * Material: PorousFlowAqueousPreDisChemistry
      95             :    * Material: PorousFlowAqueousPreDisMineral
      96             :    *
      97             :    */
      98           0 :   GenericReal<is_ad> res = 0.0;
      99     1866736 :   for (unsigned r = 0; r < _dictator.numAqueousKinetic(); ++r)
     100      935960 :     res += _stoichiometry[r] * _mineral_density[r] * _reaction_rate[_i][r];
     101      930776 :   return _test[_i][_qp] * res * _porosity_old[_i] * _saturation[_i][_aq_ph];
     102             : }
     103             : 
     104             : template <bool is_ad>
     105             : Real
     106     1476576 : PorousFlowPreDisTempl<is_ad>::computeQpJacobian()
     107             : {
     108             :   if constexpr (!is_ad)
     109             :   {
     110             :     /// If the variable is not a PorousFlow variable (very unusual), the diag Jacobian terms are 0
     111     1476576 :     if (_dictator.notPorousFlowVariable(_var.number()))
     112             :       return 0.0;
     113     1476576 :     return computeQpJac(_dictator.porousFlowVariableNum(_var.number()));
     114             :   }
     115           0 :   return 0.0;
     116             : }
     117             : 
     118             : template <bool is_ad>
     119             : Real
     120     3386784 : PorousFlowPreDisTempl<is_ad>::computeQpOffDiagJacobian(unsigned int jvar)
     121             : {
     122             :   if constexpr (!is_ad)
     123             :   {
     124             :     /// If the variable is not a PorousFlow variable, the OffDiag Jacobian terms are 0
     125     3386784 :     if (_dictator.notPorousFlowVariable(jvar))
     126             :       return 0.0;
     127     3386784 :     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     4863360 : PorousFlowPreDisTempl<is_ad>::computeQpJac(unsigned int pvar)
     137             : {
     138             :   if constexpr (!is_ad)
     139             :   {
     140     4863360 :     if (_i != _j)
     141             :       return 0.0;
     142             : 
     143             :     Real res = 0.0;
     144             :     Real dres = 0.0;
     145     4056128 :     for (unsigned r = 0; r < _dictator.numAqueousKinetic(); ++r)
     146             :     {
     147     2028928 :       dres += _stoichiometry[r] * _mineral_density[r] * (*_dreaction_rate_dvar)[_i][r][pvar];
     148     2028928 :       res += _stoichiometry[r] * _mineral_density[r] * _reaction_rate[_i][r];
     149             :     }
     150             : 
     151     2027200 :     return _test[_i][_qp] *
     152     2027200 :            (dres * _saturation[_i][_aq_ph] + res * (*_dsaturation_dvar)[_i][_aq_ph][pvar]) *
     153     2027200 :            _porosity_old[_i];
     154             :   }
     155             :   else
     156             :     libmesh_ignore(pvar);
     157           0 :   return 0.0;
     158             : }
     159             : 
     160             : template class PorousFlowPreDisTempl<false>;
     161             : template class PorousFlowPreDisTempl<true>;