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 "GeochemicalDatabaseReader.h"
      11             : #include "GeochemicalDatabaseValidator.h"
      12             : 
      13             : #include "MooseUtils.h"
      14             : #include "Conversion.h"
      15             : #include "string"
      16             : #include <fstream>
      17             : 
      18        1002 : GeochemicalDatabaseReader::GeochemicalDatabaseReader(
      19             :     const FileName filename,
      20             :     const bool reexpress_free_electron,
      21             :     const bool use_piecewise_interpolation,
      22        1002 :     const bool remove_all_extrapolated_secondary_species)
      23        1002 :   : _filename(filename)
      24             : {
      25        2004 :   read(_filename);
      26        1002 :   validate(_filename, _root);
      27             : 
      28        1001 :   if (reexpress_free_electron)
      29         999 :     reexpressFreeElectron();
      30             : 
      31        1001 :   if (use_piecewise_interpolation && _root["Header"].contains("logk model"))
      32         933 :     _root["Header"]["logk model"] = "piecewise-linear";
      33             : 
      34        1001 :   if (remove_all_extrapolated_secondary_species)
      35          10 :     removeExtrapolatedSecondarySpecies();
      36             : 
      37        1001 :   setTemperatures();
      38        1001 :   setDebyeHuckel();
      39        1001 :   setNeutralSpeciesActivity();
      40        1005 : }
      41             : 
      42             : void
      43        1002 : GeochemicalDatabaseReader::read(const FileName filename)
      44             : {
      45        1002 :   MooseUtils::checkFileReadable(filename);
      46             : 
      47             :   // Read the JSON database
      48        1002 :   std::ifstream jsondata(filename);
      49        1002 :   jsondata >> _root;
      50        1002 : }
      51             : 
      52             : void
      53        1002 : GeochemicalDatabaseReader::validate(const FileName filename, const nlohmann::json & db)
      54             : {
      55             :   // Validate the JSON database so that we don't have to check array sizes,
      56             :   // check for conversion issues, etc when extracting data using get methods
      57        1002 :   GeochemicalDatabaseValidator dbv(filename, db);
      58        1002 :   dbv.validate();
      59        1001 : }
      60             : 
      61             : void
      62         999 : GeochemicalDatabaseReader::reexpressFreeElectron()
      63             : {
      64        1808 :   if (!_root.contains("free electron") || !_root["free electron"].contains("e-") ||
      65         809 :       !_root["free electron"]["e-"]["species"].contains("O2(g)"))
      66         202 :     return;
      67         807 :   if (!_root.contains("basis species") || !_root["basis species"].contains("O2(aq)"))
      68           9 :     return;
      69        2392 :   if (!_root.contains("gas species") || !_root["gas species"].contains("O2(g)") ||
      70        1595 :       !_root["gas species"]["O2(g)"]["species"].contains("O2(aq)") ||
      71         797 :       (_root["gas species"]["O2(g)"]["species"].size() != 1))
      72           1 :     return;
      73             : 
      74             :   // remove O2(g) in the "e-" and replace with O2(aq)
      75             :   const std::string stoi_o2g =
      76         797 :       nlohmann::to_string(_root["free electron"]["e-"]["species"]["O2(g)"]);
      77         797 :   _root["free electron"]["e-"]["species"].erase("O2(g)");
      78         797 :   _root["free electron"]["e-"]["species"]["O2(aq)"] = stoi_o2g;
      79         797 :   const Real stoi = getReal(stoi_o2g);
      80             : 
      81             :   // alter equilibrium constants
      82         797 :   if (!_root["Header"].contains("temperatures"))
      83             :     return;
      84        7173 :   for (unsigned i = 0; i < _root["Header"]["temperatures"].size(); ++i)
      85             :   {
      86        6376 :     const Real logk_e = getReal(_root["free electron"]["e-"]["logk"][i]);
      87        6376 :     const Real logk_o2 = getReal(_root["gas species"]["O2(g)"]["logk"][i]);
      88        6376 :     const Real newk = logk_e + stoi * logk_o2;
      89       12752 :     _root["free electron"]["e-"]["logk"][i] = std::to_string(newk);
      90             :   }
      91             : }
      92             : 
      93             : void
      94          10 : GeochemicalDatabaseReader::removeExtrapolatedSecondarySpecies()
      95             : {
      96          10 :   if (_root.contains("secondary species"))
      97             :   {
      98             :     std::set<std::string> remove; // items to remove
      99        4986 :     for (const auto & item : _root["secondary species"].items())
     100        4966 :       if (item.value().contains("note"))
     101        2683 :         remove.insert(item.key());
     102             : 
     103        2693 :     for (const auto & name : remove)
     104        2683 :       _root["secondary species"].erase(name);
     105             :   }
     106          10 : }
     107             : 
     108             : std::string
     109        2270 : GeochemicalDatabaseReader::getActivityModel() const
     110             : {
     111        2270 :   return _root["Header"]["activity model"];
     112             : }
     113             : 
     114             : std::string
     115           1 : GeochemicalDatabaseReader::getFugacityModel() const
     116             : {
     117           1 :   return _root["Header"]["fugacity model"];
     118             : }
     119             : 
     120             : std::string
     121       18754 : GeochemicalDatabaseReader::getLogKModel() const
     122             : {
     123       18754 :   return _root["Header"]["logk model"];
     124             : }
     125             : 
     126             : void
     127        1001 : GeochemicalDatabaseReader::setTemperatures()
     128             : {
     129        1001 :   if (_root["Header"].contains("temperatures"))
     130             :   {
     131        1001 :     auto temperatures = _root["Header"]["temperatures"];
     132        1001 :     _temperature_points.resize(temperatures.size());
     133        9009 :     for (unsigned int i = 0; i < temperatures.size(); ++i)
     134        8008 :       _temperature_points[i] = getReal(temperatures[i]);
     135             :   }
     136        1001 : }
     137             : 
     138             : const std::vector<Real> &
     139       22551 : GeochemicalDatabaseReader::getTemperatures() const
     140             : {
     141       22551 :   return _temperature_points;
     142             : }
     143             : 
     144             : std::vector<Real>
     145           1 : GeochemicalDatabaseReader::getPressures()
     146             : {
     147             :   // Read pressure points
     148           1 :   if (_root["Header"].contains("pressures"))
     149             :   {
     150           1 :     auto pressures = _root["Header"]["pressures"];
     151           1 :     _pressure_points.resize(pressures.size());
     152           9 :     for (unsigned int i = 0; i < pressures.size(); ++i)
     153           8 :       _pressure_points[i] = getReal(pressures[i]);
     154             :   }
     155             : 
     156           1 :   return _pressure_points;
     157             : }
     158             : 
     159             : void
     160        1001 : GeochemicalDatabaseReader::setDebyeHuckel()
     161             : {
     162        2002 :   if (getActivityModel() == "debye-huckel")
     163             :   {
     164        1000 :     if (_root["Header"].contains("adh"))
     165             :     {
     166        1000 :       std::vector<Real> adhvals(_root["Header"]["adh"].size());
     167        9000 :       for (unsigned int i = 0; i < _root["Header"]["adh"].size(); ++i)
     168        8000 :         adhvals[i] = getReal(_root["Header"]["adh"][i]);
     169        1000 :       _debye_huckel.adh = adhvals;
     170             :     }
     171             : 
     172        1000 :     if (_root["Header"].contains("bdh"))
     173             :     {
     174        1000 :       std::vector<Real> bdhvals(_root["Header"]["bdh"].size());
     175        9000 :       for (unsigned int i = 0; i < _root["Header"]["bdh"].size(); ++i)
     176        8000 :         bdhvals[i] = getReal(_root["Header"]["bdh"][i]);
     177        1000 :       _debye_huckel.bdh = bdhvals;
     178             :     }
     179             : 
     180        1000 :     if (_root["Header"].contains("bdot"))
     181             :     {
     182        1000 :       std::vector<Real> bdotvals(_root["Header"]["bdot"].size());
     183        9000 :       for (unsigned int i = 0; i < _root["Header"]["bdot"].size(); ++i)
     184        8000 :         bdotvals[i] = getReal(_root["Header"]["bdot"][i]);
     185        1000 :       _debye_huckel.bdot = bdotvals;
     186             :     }
     187             :   }
     188        1001 : }
     189             : 
     190             : const GeochemistryDebyeHuckel &
     191        1267 : GeochemicalDatabaseReader::getDebyeHuckel() const
     192             : {
     193        2534 :   if (getActivityModel() != "debye-huckel")
     194           1 :     mooseError("Attempted to get Debye-Huckel activity parameters but the activity model is ",
     195           1 :                getActivityModel());
     196        1266 :   return _debye_huckel;
     197             : }
     198             : 
     199             : std::map<std::string, GeochemistryElements>
     200           1 : GeochemicalDatabaseReader::getElements()
     201             : {
     202           1 :   if (_root.contains("elements"))
     203             :   {
     204           4 :     for (auto & el : _root["elements"].items())
     205             :     {
     206           6 :       _elements[el.key()].name = el.value()["name"];
     207           2 :       _elements[el.key()].molecular_weight = getReal(el.value()["molecular weight"]);
     208             :     }
     209             :   }
     210             : 
     211           1 :   return _elements;
     212             : }
     213             : 
     214             : std::map<std::string, GeochemistryBasisSpecies>
     215        6552 : GeochemicalDatabaseReader::getBasisSpecies(const std::vector<std::string> & names)
     216             : {
     217             :   // Parse the basis species specified in names
     218       13105 :   for (const auto & species : names)
     219        6554 :     if (_root["basis species"].contains(species))
     220             :     {
     221             :       GeochemistryBasisSpecies dbs;
     222             : 
     223       13106 :       auto basis_species = _root["basis species"][species];
     224             :       dbs.name = species;
     225        6553 :       dbs.radius = getReal(basis_species["radius"]);
     226        6553 :       dbs.charge = getReal(basis_species["charge"]);
     227        6553 :       dbs.molecular_weight = getReal(basis_species["molecular weight"]);
     228             : 
     229             :       std::map<std::string, Real> elements;
     230       23231 :       for (auto & el : basis_species["elements"].items())
     231       10125 :         elements[el.key()] = getReal(el.value());
     232             : 
     233             :       dbs.elements = elements;
     234             : 
     235        6553 :       _basis_species[species] = dbs;
     236        6553 :     }
     237             :     else
     238           1 :       mooseError(species + " does not exist in database " + _filename);
     239             : 
     240        6551 :   return _basis_species;
     241             : }
     242             : 
     243             : std::map<std::string, GeochemistryEquilibriumSpecies>
     244      319305 : GeochemicalDatabaseReader::getEquilibriumSpecies(const std::vector<std::string> & names)
     245             : {
     246             :   // Parse the secondary species specified in names
     247      638613 :   for (const auto & species : names)
     248      319309 :     if (_root["secondary species"].contains(species) or _root["free electron"].contains(species))
     249             :     {
     250             :       GeochemistryEquilibriumSpecies dbs;
     251             : 
     252      319308 :       auto sec_species = _root["secondary species"].contains(species)
     253      637935 :                              ? _root["secondary species"][species]
     254      638616 :                              : _root["free electron"][species];
     255             :       dbs.name = species;
     256      319308 :       dbs.radius = getReal(sec_species["radius"]);
     257      319308 :       dbs.charge = getReal(sec_species["charge"]);
     258      319308 :       dbs.molecular_weight = getReal(sec_species["molecular weight"]);
     259             : 
     260      319308 :       std::vector<Real> eq_const(sec_species["logk"].size());
     261     2873772 :       for (unsigned int i = 0; i < sec_species["logk"].size(); ++i)
     262     2554464 :         eq_const[i] = getReal(sec_species["logk"][i]);
     263             : 
     264      319308 :       dbs.equilibrium_const = eq_const;
     265             : 
     266             :       std::map<std::string, Real> basis_species;
     267     1443595 :       for (auto & bs : sec_species["species"].items())
     268      804979 :         basis_species[bs.key()] = getReal(bs.value());
     269             : 
     270             :       dbs.basis_species = basis_species;
     271             : 
     272      319308 :       _equilibrium_species[species] = dbs;
     273      319308 :     }
     274             :     else
     275           1 :       mooseError(species + " does not exist in database " + _filename);
     276             : 
     277      319304 :   return _equilibrium_species;
     278             : }
     279             : 
     280             : std::map<std::string, GeochemistryMineralSpecies>
     281        1972 : GeochemicalDatabaseReader::getMineralSpecies(const std::vector<std::string> & names)
     282             : {
     283             :   // Parse the mineral species specified in names
     284        3947 :   for (const auto & species : names)
     285        1978 :     if (_root["mineral species"].contains(species))
     286             :     {
     287             :       GeochemistryMineralSpecies dbs;
     288             : 
     289        3950 :       auto mineral_species = _root["mineral species"][species];
     290             :       dbs.name = species;
     291        1975 :       dbs.molecular_weight = getReal(mineral_species["molecular weight"]);
     292        1975 :       dbs.molecular_volume = getReal(mineral_species["molar volume"]);
     293             : 
     294        1975 :       std::vector<Real> eq_const(mineral_species["logk"].size());
     295       17775 :       for (unsigned int i = 0; i < mineral_species["logk"].size(); ++i)
     296       15800 :         eq_const[i] = getReal(mineral_species["logk"][i]);
     297             : 
     298        1975 :       dbs.equilibrium_const = eq_const;
     299             : 
     300             :       std::map<std::string, Real> basis_species;
     301       10681 :       for (auto & bs : mineral_species["species"].items())
     302        6731 :         basis_species[bs.key()] = getReal(bs.value());
     303             : 
     304             :       dbs.basis_species = basis_species;
     305             : 
     306             :       // recover sorption information, if any
     307             :       std::map<std::string, Real> species_and_sorbing_density;
     308        1975 :       dbs.surface_area = 0.0;
     309        1975 :       if (_root["sorbing minerals"].contains(species))
     310             :       {
     311         120 :         auto sorbing_mineral = _root["sorbing minerals"][species];
     312         120 :         dbs.surface_area = getReal(sorbing_mineral["surface area"]);
     313             : 
     314         480 :         for (auto & site : sorbing_mineral["sorbing sites"].items())
     315         240 :           species_and_sorbing_density[site.key()] = getReal(site.value());
     316             :       }
     317             :       dbs.sorption_sites = species_and_sorbing_density;
     318             : 
     319        1975 :       _mineral_species[species] = dbs;
     320        1975 :     }
     321             :     else
     322           3 :       mooseError(species + " does not exist in database " + _filename);
     323             : 
     324        1969 :   return _mineral_species;
     325             : }
     326             : 
     327             : std::map<std::string, GeochemistryGasSpecies>
     328         351 : GeochemicalDatabaseReader::getGasSpecies(const std::vector<std::string> & names)
     329             : {
     330             :   // Parse the gas species specified in names
     331         701 :   for (const auto & species : names)
     332         352 :     if (_root["gas species"].contains(species))
     333             :     {
     334             :       GeochemistryGasSpecies dbs;
     335             : 
     336         700 :       auto gas_species = _root["gas species"][species];
     337             :       dbs.name = species;
     338         350 :       dbs.molecular_weight = getReal(gas_species["molecular weight"]);
     339             : 
     340         350 :       std::vector<Real> eq_const(gas_species["logk"].size());
     341        3150 :       for (unsigned int i = 0; i < gas_species["logk"].size(); ++i)
     342        2800 :         eq_const[i] = getReal(gas_species["logk"][i]);
     343             : 
     344         350 :       dbs.equilibrium_const = eq_const;
     345             : 
     346             :       std::map<std::string, Real> basis_species;
     347        1420 :       for (auto & bs : gas_species["species"].items())
     348         720 :         basis_species[bs.key()] = getReal(bs.value());
     349             : 
     350             :       dbs.basis_species = basis_species;
     351             : 
     352             :       // Optional fugacity coefficients
     353         350 :       if (gas_species.contains("chi"))
     354             :       {
     355         191 :         std::vector<Real> chi(gas_species["chi"].size());
     356        1337 :         for (unsigned int i = 0; i < gas_species["chi"].size(); ++i)
     357        1146 :           chi[i] = getReal(gas_species["chi"][i]);
     358             : 
     359         191 :         dbs.chi = chi;
     360             :       }
     361             : 
     362         350 :       if (gas_species.contains("Pcrit"))
     363         350 :         dbs.Pcrit = getReal(gas_species["Pcrit"]);
     364             : 
     365         350 :       if (gas_species.contains("Tcrit"))
     366         350 :         dbs.Tcrit = getReal(gas_species["Tcrit"]);
     367             : 
     368         350 :       if (gas_species.contains("omega"))
     369         350 :         dbs.omega = getReal(gas_species["omega"]);
     370             : 
     371         350 :       _gas_species[species] = dbs;
     372         350 :     }
     373             :     else
     374           2 :       mooseError(species + " does not exist in database " + _filename);
     375             : 
     376         349 :   return _gas_species;
     377             : }
     378             : 
     379             : std::map<std::string, GeochemistryRedoxSpecies>
     380       31597 : GeochemicalDatabaseReader::getRedoxSpecies(const std::vector<std::string> & names)
     381             : {
     382             :   // Parse the redox species specified in names
     383       63192 :   for (const auto & species : names)
     384       31597 :     if (_root["redox couples"].contains(species))
     385             :     {
     386             :       GeochemistryRedoxSpecies dbs;
     387             : 
     388       63190 :       auto redox_species = _root["redox couples"][species];
     389             :       dbs.name = species;
     390       31595 :       dbs.radius = getReal(redox_species["radius"]);
     391       31595 :       dbs.charge = getReal(redox_species["charge"]);
     392       31595 :       dbs.molecular_weight = getReal(redox_species["molecular weight"]);
     393             : 
     394       31595 :       std::vector<Real> eq_const(redox_species["logk"].size());
     395      284355 :       for (unsigned int i = 0; i < redox_species["logk"].size(); ++i)
     396      252760 :         eq_const[i] = getReal(redox_species["logk"][i]);
     397             : 
     398       31595 :       dbs.equilibrium_const = eq_const;
     399             : 
     400             :       std::map<std::string, Real> basis_species;
     401      182636 :       for (auto & bs : redox_species["species"].items())
     402      119446 :         basis_species[bs.key()] = getReal(bs.value());
     403             : 
     404             :       dbs.basis_species = basis_species;
     405             : 
     406       31595 :       _redox_species[species] = dbs;
     407       31595 :     }
     408             :     else
     409           2 :       mooseError(species + " does not exist in database " + _filename);
     410             : 
     411       31595 :   return _redox_species;
     412             : }
     413             : 
     414             : std::map<std::string, GeochemistryOxideSpecies>
     415           2 : GeochemicalDatabaseReader::getOxideSpecies(const std::vector<std::string> & names)
     416             : {
     417             :   // Parse the oxide species specified in names
     418           3 :   for (auto & species : names)
     419           2 :     if (_root["oxides"].contains(species))
     420             :     {
     421             :       GeochemistryOxideSpecies dbs;
     422             : 
     423           2 :       auto oxide_species = _root["oxides"][species];
     424             :       dbs.name = species;
     425           1 :       dbs.molecular_weight = getReal(oxide_species["molecular weight"]);
     426             : 
     427             :       std::map<std::string, Real> basis_species;
     428           5 :       for (auto & bs : oxide_species["species"].items())
     429           3 :         basis_species[bs.key()] = getReal(bs.value());
     430             : 
     431             :       dbs.basis_species = basis_species;
     432             : 
     433           1 :       _oxide_species[species] = dbs;
     434           1 :     }
     435             :     else
     436           1 :       mooseError(species + " does not exist in database " + _filename);
     437             : 
     438           1 :   return _oxide_species;
     439             : }
     440             : 
     441             : std::map<std::string, GeochemistrySurfaceSpecies>
     442       11797 : GeochemicalDatabaseReader::getSurfaceSpecies(const std::vector<std::string> & names)
     443             : {
     444             :   // Parse the secondary species specified in names
     445       23592 :   for (const auto & species : names)
     446       11797 :     if (_root["surface species"].contains(species))
     447             :     {
     448             :       GeochemistrySurfaceSpecies dbs;
     449             : 
     450       23590 :       auto surface_species = _root["surface species"][species];
     451             :       dbs.name = species;
     452       11795 :       dbs.charge = getReal(surface_species["charge"]);
     453       11795 :       dbs.molecular_weight = getReal(surface_species["molecular weight"]);
     454       11795 :       dbs.log10K = getReal(surface_species["log K"]);
     455       11795 :       dbs.dlog10KdT = getReal(surface_species["dlogK/dT"]);
     456             : 
     457             :       std::map<std::string, Real> basis_species;
     458       58341 :       for (auto & bs : surface_species["species"].items())
     459       34751 :         basis_species[bs.key()] = getReal(bs.value());
     460             : 
     461             :       dbs.basis_species = basis_species;
     462             : 
     463       11795 :       _surface_species[species] = dbs;
     464       11795 :     }
     465             :     else
     466           2 :       mooseError(species + " does not exist in database " + _filename);
     467             : 
     468       11795 :   return _surface_species;
     469             : }
     470             : 
     471             : void
     472        1001 : GeochemicalDatabaseReader::setNeutralSpeciesActivity()
     473             : {
     474        1001 :   if (_root["Header"].contains("neutral species"))
     475             :   {
     476        1000 :     auto neutral_species = _root["Header"]["neutral species"];
     477        3999 :     for (auto & ns : neutral_species.items())
     478             :     {
     479             :       std::vector<std::vector<Real>> coeffs;
     480             : 
     481       12557 :       for (auto & nsac : ns.value().items())
     482             :       {
     483       17118 :         if (nsac.key() == "note")
     484         804 :           continue;
     485        7755 :         std::vector<Real> coeffvec(nsac.value().size());
     486             : 
     487       69795 :         for (unsigned int i = 0; i < coeffvec.size(); ++i)
     488       62040 :           coeffvec[i] = getReal(nsac.value()[i]);
     489             : 
     490        7755 :         coeffs.push_back(coeffvec);
     491             :       }
     492             : 
     493             :       // GeochemistryNeutralSpeciesActivity expects four vectos, so
     494             :       // add empty vectors if coeffs.size() != 4
     495        1999 :       coeffs.resize(4, {});
     496             : 
     497        1999 :       GeochemistryNeutralSpeciesActivity nsa(coeffs);
     498             : 
     499        1999 :       _neutral_species_activity[ns.key()] = nsa;
     500        1999 :     }
     501             :   }
     502        1001 : }
     503             : 
     504             : const std::map<std::string, GeochemistryNeutralSpeciesActivity> &
     505        5060 : GeochemicalDatabaseReader::getNeutralSpeciesActivity() const
     506             : {
     507        5060 :   if (!_root["Header"].contains("neutral species"))
     508           1 :     mooseError("No neutral species activity coefficients in database");
     509        5059 :   return _neutral_species_activity;
     510             : }
     511             : 
     512             : std::vector<std::string>
     513           2 : GeochemicalDatabaseReader::equilibriumReactions(const std::vector<std::string> & names) const
     514             : {
     515           2 :   std::vector<std::map<std::string, Real>> basis_species(names.size());
     516             : 
     517           7 :   for (unsigned int i = 0; i < names.size(); ++i)
     518             :   {
     519           6 :     auto species = names[i];
     520             : 
     521           6 :     if (_root["secondary species"].contains(species))
     522             :     {
     523           5 :       auto sec_species = _root["secondary species"][species];
     524             : 
     525             :       // The basis species in this reaction
     526             :       std::map<std::string, Real> this_basis_species;
     527          23 :       for (auto & bs : sec_species["species"].items())
     528          13 :         this_basis_species[bs.key()] = getReal(bs.value());
     529             : 
     530             :       basis_species[i] = this_basis_species;
     531             :     }
     532             :     else
     533           2 :       mooseError(species + " does not exist in database " + _filename);
     534             :   }
     535             : 
     536           1 :   auto reactions = printReactions(names, basis_species);
     537             : 
     538           1 :   return reactions;
     539           2 : }
     540             : 
     541             : std::vector<std::string>
     542           2 : GeochemicalDatabaseReader::mineralReactions(const std::vector<std::string> & names) const
     543             : {
     544           2 :   std::vector<std::map<std::string, Real>> basis_species(names.size());
     545             : 
     546           3 :   for (unsigned int i = 0; i < names.size(); ++i)
     547             :   {
     548           2 :     const auto species = names[i];
     549             : 
     550           2 :     if (_root["mineral species"].contains(species))
     551             :     {
     552           1 :       auto min_species = _root["mineral species"][species];
     553             : 
     554             :       // The basis species in this reaction
     555             :       std::map<std::string, Real> this_basis_species;
     556           5 :       for (auto & bs : min_species["species"].items())
     557           3 :         this_basis_species[bs.key()] = getReal(bs.value());
     558             : 
     559             :       basis_species[i] = this_basis_species;
     560             :     }
     561             :     else
     562           2 :       mooseError(species + " does not exist in database " + _filename);
     563             :   }
     564             : 
     565           1 :   auto reactions = printReactions(names, basis_species);
     566             : 
     567           1 :   return reactions;
     568           2 : }
     569             : 
     570             : std::vector<std::string>
     571           2 : GeochemicalDatabaseReader::gasReactions(const std::vector<std::string> & names) const
     572             : {
     573           2 :   std::vector<std::map<std::string, Real>> basis_species(names.size());
     574             : 
     575           4 :   for (unsigned int i = 0; i < names.size(); ++i)
     576             :   {
     577           3 :     const auto species = names[i];
     578             : 
     579           3 :     if (_root["gas species"].contains(species))
     580             :     {
     581           2 :       auto gas_species = _root["gas species"][species];
     582             : 
     583             :       // The basis species in this reaction
     584             :       std::map<std::string, Real> this_basis_species;
     585           6 :       for (auto & bs : gas_species["species"].items())
     586           2 :         this_basis_species[bs.key()] = getReal(bs.value());
     587             : 
     588             :       basis_species[i] = this_basis_species;
     589             :     }
     590             :     else
     591           2 :       mooseError(species + " does not exist in database " + _filename);
     592             :   }
     593             : 
     594           1 :   auto reactions = printReactions(names, basis_species);
     595             : 
     596           1 :   return reactions;
     597           2 : }
     598             : 
     599             : std::vector<std::string>
     600           2 : GeochemicalDatabaseReader::redoxReactions(const std::vector<std::string> & names) const
     601             : {
     602           2 :   std::vector<std::map<std::string, Real>> basis_species(names.size());
     603             : 
     604           4 :   for (unsigned int i = 0; i < names.size(); ++i)
     605             :   {
     606           3 :     const auto species = names[i];
     607             : 
     608           3 :     if (_root["redox couples"].contains(species))
     609             :     {
     610           2 :       auto redox_species = _root["redox couples"][species];
     611             : 
     612             :       // The basis species in this reaction
     613             :       std::map<std::string, Real> this_basis_species;
     614          12 :       for (auto & bs : redox_species["species"].items())
     615           8 :         this_basis_species[bs.key()] = getReal(bs.value());
     616             : 
     617             :       basis_species[i] = this_basis_species;
     618             :     }
     619             :     else
     620           2 :       mooseError(species + " does not exist in database " + _filename);
     621             :   }
     622             : 
     623           1 :   auto reactions = printReactions(names, basis_species);
     624             : 
     625           1 :   return reactions;
     626           2 : }
     627             : 
     628             : std::vector<std::string>
     629           2 : GeochemicalDatabaseReader::oxideReactions(const std::vector<std::string> & names) const
     630             : {
     631           2 :   std::vector<std::map<std::string, Real>> basis_species(names.size());
     632             : 
     633           3 :   for (unsigned int i = 0; i < names.size(); ++i)
     634             :   {
     635           2 :     const auto species = names[i];
     636             : 
     637           2 :     if (_root["oxides"].contains(species))
     638             :     {
     639           1 :       auto oxide_species = _root["oxides"][species];
     640             : 
     641             :       // The basis species in this reaction
     642             :       std::map<std::string, Real> this_basis_species;
     643           5 :       for (auto & bs : oxide_species["species"].items())
     644           3 :         this_basis_species[bs.key()] = getReal(bs.value());
     645             : 
     646             :       basis_species[i] = this_basis_species;
     647             :     }
     648             :     else
     649           2 :       mooseError(species + " does not exist in database " + _filename);
     650             :   }
     651             : 
     652           1 :   auto reactions = printReactions(names, basis_species);
     653             : 
     654           1 :   return reactions;
     655           2 : }
     656             : 
     657             : std::vector<std::string>
     658           5 : GeochemicalDatabaseReader::printReactions(
     659             :     const std::vector<std::string> & names,
     660             :     const std::vector<std::map<std::string, Real>> & basis_species) const
     661             : {
     662             :   std::vector<std::string> reactions;
     663             : 
     664          16 :   for (unsigned int i = 0; i < names.size(); ++i)
     665             :   {
     666          11 :     std::string reaction = "";
     667          40 :     for (auto & bs : basis_species[i])
     668             :     {
     669          29 :       if (bs.second < 0.0)
     670             :       {
     671          10 :         if (bs.second == -1.0)
     672          12 :           reaction += " - " + bs.first;
     673             :         else
     674           8 :           reaction += " " + Moose::stringify(bs.second) + bs.first;
     675             :       }
     676             :       else
     677             :       {
     678          19 :         if (bs.second == 1.0)
     679          30 :           reaction += " + " + bs.first;
     680             :         else
     681           8 :           reaction += " + " + Moose::stringify(bs.second) + bs.first;
     682             :       }
     683             :     }
     684             : 
     685             :     // Trim off leading +
     686          11 :     if (reaction.size() > 1 && reaction[1] == '+')
     687           9 :       reaction.erase(1, 2);
     688             : 
     689          22 :     reactions.push_back(names[i] + " =" + reaction);
     690             :   }
     691             : 
     692           5 :   return reactions;
     693           0 : }
     694             : 
     695             : std::vector<std::string>
     696           2 : GeochemicalDatabaseReader::mineralSpeciesNames() const
     697             : {
     698             :   std::vector<std::string> names;
     699           2 :   if (_root.contains("mineral species"))
     700          17 :     for (auto & item : _root["mineral species"].items())
     701          13 :       names.push_back(item.key());
     702           2 :   return names;
     703           0 : }
     704             : 
     705             : std::vector<std::string>
     706        1896 : GeochemicalDatabaseReader::secondarySpeciesNames() const
     707             : {
     708             :   std::vector<std::string> names;
     709        1896 :   if (_root.contains("secondary species"))
     710      505578 :     for (auto & item : _root["secondary species"].items())
     711      501788 :       names.push_back(item.key());
     712             : 
     713        1896 :   if (_root.contains("free electron"))
     714        5121 :     for (const auto & nm : _root["free electron"].items())
     715        1707 :       names.push_back(nm.key());
     716        1896 :   return names;
     717           0 : }
     718             : 
     719             : std::vector<std::string>
     720        1217 : GeochemicalDatabaseReader::redoxCoupleNames() const
     721             : {
     722             :   std::vector<std::string> names;
     723        1217 :   if (_root.contains("redox couples"))
     724       33504 :     for (const auto & item : _root["redox couples"].items())
     725       31070 :       names.push_back(item.key());
     726        1217 :   return names;
     727           0 : }
     728             : 
     729             : std::vector<std::string>
     730        1217 : GeochemicalDatabaseReader::surfaceSpeciesNames() const
     731             : {
     732             :   std::vector<std::string> names;
     733        1217 :   if (_root.contains("surface species"))
     734       13060 :     for (const auto & item : _root["surface species"].items())
     735       11792 :       names.push_back(item.key());
     736        1217 :   return names;
     737           0 : }
     738             : 
     739             : const FileName &
     740           2 : GeochemicalDatabaseReader::filename() const
     741             : {
     742           2 :   return _filename;
     743             : }
     744             : 
     745             : bool
     746        7274 : GeochemicalDatabaseReader::isBasisSpecies(const std::string & name) const
     747             : {
     748        7274 :   return _root["basis species"].contains(name);
     749             : }
     750             : 
     751             : bool
     752        3870 : GeochemicalDatabaseReader::isRedoxSpecies(const std::string & name) const
     753             : {
     754        3870 :   return _root["redox couples"].contains(name);
     755             : }
     756             : 
     757             : bool
     758          13 : GeochemicalDatabaseReader::isSorbingMineral(const std::string & name) const
     759             : {
     760          13 :   return _root["sorbing minerals"].contains(name);
     761             : }
     762             : 
     763             : bool
     764          28 : GeochemicalDatabaseReader::isSecondarySpecies(const std::string & name) const
     765             : {
     766          28 :   return _root["secondary species"].contains(name) || _root["free electron"].contains(name);
     767             : }
     768             : 
     769             : bool
     770          13 : GeochemicalDatabaseReader::isGasSpecies(const std::string & name) const
     771             : {
     772          13 :   return _root["gas species"].contains(name);
     773             : }
     774             : 
     775             : bool
     776          13 : GeochemicalDatabaseReader::isMineralSpecies(const std::string & name) const
     777             : {
     778          13 :   return _root["mineral species"].contains(name);
     779             : }
     780             : 
     781             : bool
     782          13 : GeochemicalDatabaseReader::isOxideSpecies(const std::string & name) const
     783             : {
     784          13 :   return _root["oxides"].contains(name);
     785             : }
     786             : 
     787             : bool
     788          13 : GeochemicalDatabaseReader::isSurfaceSpecies(const std::string & name) const
     789             : {
     790          13 :   return _root["surface species"].contains(name);
     791             : }
     792             : 
     793             : std::string
     794           2 : GeochemicalDatabaseReader::getSpeciesData(const std::string name) const
     795             : {
     796             :   std::string output;
     797          26 :   for (auto & item : _root.items())
     798          22 :     if (_root[item.key()].contains(name))
     799             :     {
     800           1 :       std::ostringstream os;
     801           2 :       os << item.value()[name].dump(4);
     802           1 :       output = os.str();
     803           1 :     }
     804             : 
     805           2 :   if (output.empty())
     806           1 :     mooseError(name + " is not a species in the database");
     807             : 
     808           3 :   return name + ":\n" + output;
     809             : }
     810             : 
     811             : Real
     812     5036620 : GeochemicalDatabaseReader::getReal(const nlohmann::json & node)
     813             : {
     814     5036620 :   if (node.is_string())
     815      285542 :     return MooseUtils::convert<Real>(node);
     816     4893849 :   return node;
     817             : }