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 "GeochemistryIonicStrength.h"
      11             : #include "libmesh/utility.h"
      12             : 
      13        1276 : GeochemistryIonicStrength::GeochemistryIonicStrength(Real max_ionic_strength,
      14             :                                                      Real max_stoichiometric_ionic_strength,
      15             :                                                      bool use_only_basis_molality,
      16        1276 :                                                      bool use_only_Cl_molality)
      17        1276 :   : _max_ionic_strength(max_ionic_strength),
      18        1276 :     _max_stoichiometric_ionic_strength(max_stoichiometric_ionic_strength),
      19        1276 :     _use_only_basis_molality(use_only_basis_molality),
      20        1276 :     _use_only_Cl_molality(use_only_Cl_molality)
      21             : {
      22        1276 : }
      23             : 
      24             : Real
      25       63341 : GeochemistryIonicStrength::ionicStrength(const ModelGeochemicalDatabase & mgd,
      26             :                                          const std::vector<Real> & basis_species_molality,
      27             :                                          const std::vector<Real> & eqm_species_molality,
      28             :                                          const std::vector<Real> & kin_species_molality) const
      29             : {
      30       63341 :   const unsigned num_basis = mgd.basis_species_charge.size();
      31       63341 :   const unsigned num_eqm = mgd.eqm_species_charge.size();
      32       63341 :   const unsigned num_kin = mgd.kin_species_charge.size();
      33       63341 :   if (num_basis != basis_species_molality.size())
      34           1 :     mooseError("Ionic strength calculation: Number of basis species in mgd not equal to the size "
      35             :                "of basis_species_molality");
      36       63340 :   if (num_eqm != eqm_species_molality.size())
      37           1 :     mooseError("Ionic strength calculation: Number of equilibrium species in mgd not equal to the "
      38             :                "size of eqm_species_molality");
      39       63339 :   if (num_kin != kin_species_molality.size())
      40           1 :     mooseError("Ionic strength calculation: Number of kinetic species in mgd not equal to the size "
      41             :                "of kin_species_molality");
      42             : 
      43             :   Real ionic_strength = 0.0;
      44      488518 :   for (unsigned i = 0; i < num_basis; ++i)
      45      425180 :     ionic_strength += Utility::pow<2>(mgd.basis_species_charge[i]) * basis_species_molality[i];
      46       63338 :   if (!_use_only_basis_molality)
      47             :   {
      48     1744240 :     for (unsigned i = 0; i < num_eqm; ++i)
      49     1680906 :       if (!mgd.surface_sorption_related[i])
      50     1619394 :         ionic_strength += Utility::pow<2>(mgd.eqm_species_charge[i]) * eqm_species_molality[i];
      51       71967 :     for (unsigned i = 0; i < num_kin; ++i)
      52        8633 :       ionic_strength += Utility::pow<2>(mgd.kin_species_charge[i]) * kin_species_molality[i] /
      53             :                         basis_species_molality[0]; // kin_species_molality is actually the number of
      54             :                                                    // moles of the kinetic species
      55             :   }
      56             : 
      57      119328 :   return std::max(0.0, std::min(_max_ionic_strength, 0.5 * ionic_strength));
      58             : }
      59             : 
      60             : Real
      61       55320 : GeochemistryIonicStrength::stoichiometricIonicStrength(
      62             :     const ModelGeochemicalDatabase & mgd,
      63             :     const std::vector<Real> & basis_species_molality,
      64             :     const std::vector<Real> & eqm_species_molality,
      65             :     const std::vector<Real> & kin_species_molality) const
      66             : {
      67       55320 :   const unsigned num_basis = mgd.basis_species_charge.size();
      68       55320 :   const unsigned num_eqm = mgd.eqm_species_charge.size();
      69       55320 :   const unsigned num_kin = mgd.kin_species_charge.size();
      70       55320 :   if (num_basis != basis_species_molality.size())
      71           1 :     mooseError("Stoichiometric ionic strength calculation: Number of basis species in mgd not "
      72             :                "equal to the size of basis_species_molality");
      73       55319 :   if (num_eqm != eqm_species_molality.size())
      74           1 :     mooseError("Stoichiometric ionic strength calculation: Number of equilibrium species in mgd "
      75             :                "not equal to the size of eqm_species_molality");
      76       55318 :   if (num_kin != kin_species_molality.size())
      77           1 :     mooseError("Stoichiometric ionic strength calculation: Number of kinetic species in mgd not "
      78             :                "equal to the size of kin_species_molality");
      79             : 
      80       55317 :   Real ionic_strength = 0.0;
      81       55317 :   if (_use_only_Cl_molality)
      82             :   {
      83       55446 :     if (mgd.basis_species_index.count("Cl-"))
      84       55440 :       ionic_strength = basis_species_molality[mgd.basis_species_index.at("Cl-")];
      85           6 :     else if (mgd.eqm_species_index.count("Cl-"))
      86           2 :       ionic_strength = eqm_species_molality[mgd.eqm_species_index.at("Cl-")];
      87           4 :     else if (mgd.kin_species_index.count("Cl-"))
      88           2 :       ionic_strength = kin_species_molality[mgd.kin_species_index.at("Cl-")];
      89             :     else
      90           1 :       mooseError("GeochemistryIonicStrength: attempting to compute stoichiometric ionic strength "
      91             :                  "using only the Cl- molality, but Cl- does not appear in the geochemical system");
      92       81212 :     return std::max(0.0, std::min(_max_stoichiometric_ionic_strength, ionic_strength));
      93             :   }
      94             : 
      95      122349 :   for (unsigned i = 0; i < num_basis; ++i)
      96       94755 :     ionic_strength += Utility::pow<2>(mgd.basis_species_charge[i]) * basis_species_molality[i];
      97       27594 :   if (!_use_only_basis_molality)
      98             :   {
      99      132516 :     for (unsigned i = 0; i < num_eqm; ++i)
     100      104926 :       if (!mgd.surface_sorption_related[i])
     101             :       {
     102             : 
     103       98286 :         if (mgd.eqm_species_charge[i] != 0.0)
     104       53298 :           ionic_strength += Utility::pow<2>(mgd.eqm_species_charge[i]) * eqm_species_molality[i];
     105             :         else
     106             :         {
     107      349749 :           for (unsigned j = 0; j < num_basis; ++j)
     108      304761 :             ionic_strength += Utility::pow<2>(mgd.basis_species_charge[j]) *
     109      304761 :                               eqm_species_molality[i] * mgd.eqm_stoichiometry(i, j);
     110             :         }
     111             :       }
     112       29671 :     for (unsigned i = 0; i < num_kin;
     113             :          ++i) // kin_species_molality is actually the number of moles, not molality
     114        2081 :       if (mgd.kin_species_charge[i] != 0.0)
     115           7 :         ionic_strength += Utility::pow<2>(mgd.kin_species_charge[i]) * kin_species_molality[i] /
     116             :                           basis_species_molality[0];
     117             :   }
     118             : 
     119       80867 :   return std::max(0.0, std::min(_max_stoichiometric_ionic_strength, 0.5 * ionic_strength));
     120             : }
     121             : 
     122             : void
     123       52211 : GeochemistryIonicStrength::setMaxIonicStrength(Real max_ionic_strength)
     124             : {
     125       52211 :   _max_ionic_strength = max_ionic_strength;
     126       52211 : }
     127             : 
     128             : Real
     129           2 : GeochemistryIonicStrength::getMaxIonicStrength() const
     130             : {
     131           2 :   return _max_ionic_strength;
     132             : }
     133             : 
     134             : void
     135       52212 : GeochemistryIonicStrength::setMaxStoichiometricIonicStrength(Real max_stoichiometric_ionic_strength)
     136             : {
     137       52212 :   _max_stoichiometric_ionic_strength = max_stoichiometric_ionic_strength;
     138       52212 : }
     139             : 
     140             : Real
     141           2 : GeochemistryIonicStrength::getMaxStoichiometricIonicStrength() const
     142             : {
     143           2 :   return _max_stoichiometric_ionic_strength;
     144             : }
     145             : 
     146             : void
     147           3 : GeochemistryIonicStrength::setUseOnlyBasisMolality(bool use_only_basis_molality)
     148             : {
     149           3 :   _use_only_basis_molality = use_only_basis_molality;
     150           3 : }
     151             : 
     152             : Real
     153           7 : GeochemistryIonicStrength::getUseOnlyBasisMolality() const
     154             : {
     155           7 :   return _use_only_basis_molality;
     156             : }
     157             : 
     158             : void
     159           1 : GeochemistryIonicStrength::setUseOnlyClMolality(bool use_only_Cl_molality)
     160             : {
     161           1 :   _use_only_Cl_molality = use_only_Cl_molality;
     162           1 : }
     163             : 
     164             : Real
     165           7 : GeochemistryIonicStrength::getUseOnlyClMolality() const
     166             : {
     167           7 :   return _use_only_Cl_molality;
     168             : }