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