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         949 : GeochemicalDatabaseReader::GeochemicalDatabaseReader(
      19             :     const FileName filename,
      20             :     const bool reexpress_free_electron,
      21             :     const bool use_piecewise_interpolation,
      22         949 :     const bool remove_all_extrapolated_secondary_species)
      23         949 :   : _filename(filename)
      24             : {
      25        1898 :   read(_filename);
      26         949 :   validate(_filename, _root);
      27             : 
      28         948 :   if (reexpress_free_electron)
      29         946 :     reexpressFreeElectron();
      30             : 
      31         948 :   if (use_piecewise_interpolation && _root["Header"].contains("logk model"))
      32         903 :     _root["Header"]["logk model"] = "piecewise-linear";
      33             : 
      34         948 :   if (remove_all_extrapolated_secondary_species)
      35           9 :     removeExtrapolatedSecondarySpecies();
      36             : 
      37         948 :   setTemperatures();
      38         948 :   setDebyeHuckel();
      39         948 :   setNeutralSpeciesActivity();
      40         952 : }
      41             : 
      42             : void
      43         949 : GeochemicalDatabaseReader::read(const FileName filename)
      44             : {
      45         949 :   MooseUtils::checkFileReadable(filename);
      46             : 
      47             :   // Read the JSON database
      48         949 :   std::ifstream jsondata(filename);
      49         949 :   jsondata >> _root;
      50         949 : }
      51             : 
      52             : void
      53         949 : 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         949 :   GeochemicalDatabaseValidator dbv(filename, db);
      58         949 :   dbv.validate();
      59         948 : }
      60             : 
      61             : void
      62         946 : GeochemicalDatabaseReader::reexpressFreeElectron()
      63             : {
      64        1717 :   if (!_root.contains("free electron") || !_root["free electron"].contains("e-") ||
      65         771 :       !_root["free electron"]["e-"]["species"].contains("O2(g)"))
      66         186 :     return;
      67         769 :   if (!_root.contains("basis species") || !_root["basis species"].contains("O2(aq)"))
      68           8 :     return;
      69        2281 :   if (!_root.contains("gas species") || !_root["gas species"].contains("O2(g)") ||
      70        1521 :       !_root["gas species"]["O2(g)"]["species"].contains("O2(aq)") ||
      71         760 :       (_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         760 :       nlohmann::to_string(_root["free electron"]["e-"]["species"]["O2(g)"]);
      77         760 :   _root["free electron"]["e-"]["species"].erase("O2(g)");
      78         760 :   _root["free electron"]["e-"]["species"]["O2(aq)"] = stoi_o2g;
      79         760 :   const Real stoi = getReal(stoi_o2g);
      80             : 
      81             :   // alter equilibrium constants
      82         760 :   if (!_root["Header"].contains("temperatures"))
      83             :     return;
      84        6840 :   for (unsigned i = 0; i < _root["Header"]["temperatures"].size(); ++i)
      85             :   {
      86        6080 :     const Real logk_e = getReal(_root["free electron"]["e-"]["logk"][i]);
      87        6080 :     const Real logk_o2 = getReal(_root["gas species"]["O2(g)"]["logk"][i]);
      88        6080 :     const Real newk = logk_e + stoi * logk_o2;
      89       12160 :     _root["free electron"]["e-"]["logk"][i] = std::to_string(newk);
      90             :   }
      91             : }
      92             : 
      93             : void
      94           9 : GeochemicalDatabaseReader::removeExtrapolatedSecondarySpecies()
      95             : {
      96           9 :   if (_root.contains("secondary species"))
      97             :   {
      98             :     std::set<std::string> remove; // items to remove
      99        4433 :     for (const auto & item : _root["secondary species"].items())
     100        4415 :       if (item.value().contains("note"))
     101        2385 :         remove.insert(item.key());
     102             : 
     103        2394 :     for (const auto & name : remove)
     104        2385 :       _root["secondary species"].erase(name);
     105             :   }
     106           9 : }
     107             : 
     108             : std::string
     109        2161 : GeochemicalDatabaseReader::getActivityModel() const
     110             : {
     111        2161 :   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       17834 : GeochemicalDatabaseReader::getLogKModel() const
     122             : {
     123       17834 :   return _root["Header"]["logk model"];
     124             : }
     125             : 
     126             : void
     127         948 : GeochemicalDatabaseReader::setTemperatures()
     128             : {
     129         948 :   if (_root["Header"].contains("temperatures"))
     130             :   {
     131         948 :     auto temperatures = _root["Header"]["temperatures"];
     132         948 :     _temperature_points.resize(temperatures.size());
     133        8532 :     for (unsigned int i = 0; i < temperatures.size(); ++i)
     134        7584 :       _temperature_points[i] = getReal(temperatures[i]);
     135             :   }
     136         948 : }
     137             : 
     138             : const std::vector<Real> &
     139       21461 : GeochemicalDatabaseReader::getTemperatures() const
     140             : {
     141       21461 :   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         948 : GeochemicalDatabaseReader::setDebyeHuckel()
     161             : {
     162        1896 :   if (getActivityModel() == "debye-huckel")
     163             :   {
     164         947 :     if (_root["Header"].contains("adh"))
     165             :     {
     166         947 :       std::vector<Real> adhvals(_root["Header"]["adh"].size());
     167        8523 :       for (unsigned int i = 0; i < _root["Header"]["adh"].size(); ++i)
     168        7576 :         adhvals[i] = getReal(_root["Header"]["adh"][i]);
     169         947 :       _debye_huckel.adh = adhvals;
     170             :     }
     171             : 
     172         947 :     if (_root["Header"].contains("bdh"))
     173             :     {
     174         947 :       std::vector<Real> bdhvals(_root["Header"]["bdh"].size());
     175        8523 :       for (unsigned int i = 0; i < _root["Header"]["bdh"].size(); ++i)
     176        7576 :         bdhvals[i] = getReal(_root["Header"]["bdh"][i]);
     177         947 :       _debye_huckel.bdh = bdhvals;
     178             :     }
     179             : 
     180         947 :     if (_root["Header"].contains("bdot"))
     181             :     {
     182         947 :       std::vector<Real> bdotvals(_root["Header"]["bdot"].size());
     183        8523 :       for (unsigned int i = 0; i < _root["Header"]["bdot"].size(); ++i)
     184        7576 :         bdotvals[i] = getReal(_root["Header"]["bdot"][i]);
     185         947 :       _debye_huckel.bdot = bdotvals;
     186             :     }
     187             :   }
     188         948 : }
     189             : 
     190             : const GeochemistryDebyeHuckel &
     191        1211 : GeochemicalDatabaseReader::getDebyeHuckel() const
     192             : {
     193        2422 :   if (getActivityModel() != "debye-huckel")
     194           1 :     mooseError("Attempted to get Debye-Huckel activity parameters but the activity model is ",
     195           1 :                getActivityModel());
     196        1210 :   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        6278 : GeochemicalDatabaseReader::getBasisSpecies(const std::vector<std::string> & names)
     216             : {
     217             :   // Parse the basis species specified in names
     218       12557 :   for (const auto & species : names)
     219        6280 :     if (_root["basis species"].contains(species))
     220             :     {
     221             :       GeochemistryBasisSpecies dbs;
     222             : 
     223       12558 :       auto basis_species = _root["basis species"][species];
     224             :       dbs.name = species;
     225        6279 :       dbs.radius = getReal(basis_species["radius"]);
     226        6279 :       dbs.charge = getReal(basis_species["charge"]);
     227        6279 :       dbs.molecular_weight = getReal(basis_species["molecular weight"]);
     228             : 
     229             :       std::map<std::string, Real> elements;
     230       22281 :       for (auto & el : basis_species["elements"].items())
     231        9723 :         elements[el.key()] = getReal(el.value());
     232             : 
     233             :       dbs.elements = elements;
     234             : 
     235        6279 :       _basis_species[species] = dbs;
     236        6279 :     }
     237             :     else
     238           1 :       mooseError(species + " does not exist in database " + _filename);
     239             : 
     240        6277 :   return _basis_species;
     241             : }
     242             : 
     243             : std::map<std::string, GeochemistryEquilibriumSpecies>
     244      299007 : GeochemicalDatabaseReader::getEquilibriumSpecies(const std::vector<std::string> & names)
     245             : {
     246             :   // Parse the secondary species specified in names
     247      598017 :   for (const auto & species : names)
     248      299011 :     if (_root["secondary species"].contains(species) or _root["free electron"].contains(species))
     249             :     {
     250             :       GeochemistryEquilibriumSpecies dbs;
     251             : 
     252      299010 :       auto sec_species = _root["secondary species"].contains(species)
     253      597349 :                              ? _root["secondary species"][species]
     254      598020 :                              : _root["free electron"][species];
     255             :       dbs.name = species;
     256      299010 :       dbs.radius = getReal(sec_species["radius"]);
     257      299010 :       dbs.charge = getReal(sec_species["charge"]);
     258      299010 :       dbs.molecular_weight = getReal(sec_species["molecular weight"]);
     259             : 
     260      299010 :       std::vector<Real> eq_const(sec_species["logk"].size());
     261     2691090 :       for (unsigned int i = 0; i < sec_species["logk"].size(); ++i)
     262     2392080 :         eq_const[i] = getReal(sec_species["logk"][i]);
     263             : 
     264      299010 :       dbs.equilibrium_const = eq_const;
     265             : 
     266             :       std::map<std::string, Real> basis_species;
     267     1351822 :       for (auto & bs : sec_species["species"].items())
     268      753802 :         basis_species[bs.key()] = getReal(bs.value());
     269             : 
     270             :       dbs.basis_species = basis_species;
     271             : 
     272      299010 :       _equilibrium_species[species] = dbs;
     273      299010 :     }
     274             :     else
     275           1 :       mooseError(species + " does not exist in database " + _filename);
     276             : 
     277      299006 :   return _equilibrium_species;
     278             : }
     279             : 
     280             : std::map<std::string, GeochemistryMineralSpecies>
     281        1900 : GeochemicalDatabaseReader::getMineralSpecies(const std::vector<std::string> & names)
     282             : {
     283             :   // Parse the mineral species specified in names
     284        3803 :   for (const auto & species : names)
     285        1906 :     if (_root["mineral species"].contains(species))
     286             :     {
     287             :       GeochemistryMineralSpecies dbs;
     288             : 
     289        3806 :       auto mineral_species = _root["mineral species"][species];
     290             :       dbs.name = species;
     291        1903 :       dbs.molecular_weight = getReal(mineral_species["molecular weight"]);
     292        1903 :       dbs.molecular_volume = getReal(mineral_species["molar volume"]);
     293             : 
     294        1903 :       std::vector<Real> eq_const(mineral_species["logk"].size());
     295       17127 :       for (unsigned int i = 0; i < mineral_species["logk"].size(); ++i)
     296       15224 :         eq_const[i] = getReal(mineral_species["logk"][i]);
     297             : 
     298        1903 :       dbs.equilibrium_const = eq_const;
     299             : 
     300             :       std::map<std::string, Real> basis_species;
     301       10292 :       for (auto & bs : mineral_species["species"].items())
     302        6486 :         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        1903 :       dbs.surface_area = 0.0;
     309        1903 :       if (_root["sorbing minerals"].contains(species))
     310             :       {
     311         113 :         auto sorbing_mineral = _root["sorbing minerals"][species];
     312         113 :         dbs.surface_area = getReal(sorbing_mineral["surface area"]);
     313             : 
     314         452 :         for (auto & site : sorbing_mineral["sorbing sites"].items())
     315         226 :           species_and_sorbing_density[site.key()] = getReal(site.value());
     316             :       }
     317             :       dbs.sorption_sites = species_and_sorbing_density;
     318             : 
     319        1903 :       _mineral_species[species] = dbs;
     320        1903 :     }
     321             :     else
     322           3 :       mooseError(species + " does not exist in database " + _filename);
     323             : 
     324        1897 :   return _mineral_species;
     325             : }
     326             : 
     327             : std::map<std::string, GeochemistryGasSpecies>
     328         344 : GeochemicalDatabaseReader::getGasSpecies(const std::vector<std::string> & names)
     329             : {
     330             :   // Parse the gas species specified in names
     331         687 :   for (const auto & species : names)
     332         345 :     if (_root["gas species"].contains(species))
     333             :     {
     334             :       GeochemistryGasSpecies dbs;
     335             : 
     336         686 :       auto gas_species = _root["gas species"][species];
     337             :       dbs.name = species;
     338         343 :       dbs.molecular_weight = getReal(gas_species["molecular weight"]);
     339             : 
     340         343 :       std::vector<Real> eq_const(gas_species["logk"].size());
     341        3087 :       for (unsigned int i = 0; i < gas_species["logk"].size(); ++i)
     342        2744 :         eq_const[i] = getReal(gas_species["logk"][i]);
     343             : 
     344         343 :       dbs.equilibrium_const = eq_const;
     345             : 
     346             :       std::map<std::string, Real> basis_species;
     347        1389 :       for (auto & bs : gas_species["species"].items())
     348         703 :         basis_species[bs.key()] = getReal(bs.value());
     349             : 
     350             :       dbs.basis_species = basis_species;
     351             : 
     352             :       // Optional fugacity coefficients
     353         343 :       if (gas_species.contains("chi"))
     354             :       {
     355         186 :         std::vector<Real> chi(gas_species["chi"].size());
     356        1302 :         for (unsigned int i = 0; i < gas_species["chi"].size(); ++i)
     357        1116 :           chi[i] = getReal(gas_species["chi"][i]);
     358             : 
     359         186 :         dbs.chi = chi;
     360             :       }
     361             : 
     362         343 :       if (gas_species.contains("Pcrit"))
     363         343 :         dbs.Pcrit = getReal(gas_species["Pcrit"]);
     364             : 
     365         343 :       if (gas_species.contains("Tcrit"))
     366         343 :         dbs.Tcrit = getReal(gas_species["Tcrit"]);
     367             : 
     368         343 :       if (gas_species.contains("omega"))
     369         343 :         dbs.omega = getReal(gas_species["omega"]);
     370             : 
     371         343 :       _gas_species[species] = dbs;
     372         343 :     }
     373             :     else
     374           2 :       mooseError(species + " does not exist in database " + _filename);
     375             : 
     376         342 :   return _gas_species;
     377             : }
     378             : 
     379             : std::map<std::string, GeochemistryRedoxSpecies>
     380       29730 : GeochemicalDatabaseReader::getRedoxSpecies(const std::vector<std::string> & names)
     381             : {
     382             :   // Parse the redox species specified in names
     383       59458 :   for (const auto & species : names)
     384       29730 :     if (_root["redox couples"].contains(species))
     385             :     {
     386             :       GeochemistryRedoxSpecies dbs;
     387             : 
     388       59456 :       auto redox_species = _root["redox couples"][species];
     389             :       dbs.name = species;
     390       29728 :       dbs.radius = getReal(redox_species["radius"]);
     391       29728 :       dbs.charge = getReal(redox_species["charge"]);
     392       29728 :       dbs.molecular_weight = getReal(redox_species["molecular weight"]);
     393             : 
     394       29728 :       std::vector<Real> eq_const(redox_species["logk"].size());
     395      267552 :       for (unsigned int i = 0; i < redox_species["logk"].size(); ++i)
     396      237824 :         eq_const[i] = getReal(redox_species["logk"][i]);
     397             : 
     398       29728 :       dbs.equilibrium_const = eq_const;
     399             : 
     400             :       std::map<std::string, Real> basis_species;
     401      171844 :       for (auto & bs : redox_species["species"].items())
     402      112388 :         basis_species[bs.key()] = getReal(bs.value());
     403             : 
     404             :       dbs.basis_species = basis_species;
     405             : 
     406       29728 :       _redox_species[species] = dbs;
     407       29728 :     }
     408             :     else
     409           2 :       mooseError(species + " does not exist in database " + _filename);
     410             : 
     411       29728 :   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       10927 : GeochemicalDatabaseReader::getSurfaceSpecies(const std::vector<std::string> & names)
     443             : {
     444             :   // Parse the secondary species specified in names
     445       21852 :   for (const auto & species : names)
     446       10927 :     if (_root["surface species"].contains(species))
     447             :     {
     448             :       GeochemistrySurfaceSpecies dbs;
     449             : 
     450       21850 :       auto surface_species = _root["surface species"][species];
     451             :       dbs.name = species;
     452       10925 :       dbs.charge = getReal(surface_species["charge"]);
     453       10925 :       dbs.molecular_weight = getReal(surface_species["molecular weight"]);
     454       10925 :       dbs.log10K = getReal(surface_species["log K"]);
     455       10925 :       dbs.dlog10KdT = getReal(surface_species["dlogK/dT"]);
     456             : 
     457             :       std::map<std::string, Real> basis_species;
     458       54006 :       for (auto & bs : surface_species["species"].items())
     459       32156 :         basis_species[bs.key()] = getReal(bs.value());
     460             : 
     461             :       dbs.basis_species = basis_species;
     462             : 
     463       10925 :       _surface_species[species] = dbs;
     464       10925 :     }
     465             :     else
     466           2 :       mooseError(species + " does not exist in database " + _filename);
     467             : 
     468       10925 :   return _surface_species;
     469             : }
     470             : 
     471             : void
     472         948 : GeochemicalDatabaseReader::setNeutralSpeciesActivity()
     473             : {
     474         948 :   if (_root["Header"].contains("neutral species"))
     475             :   {
     476         947 :     auto neutral_species = _root["Header"]["neutral species"];
     477        3787 :     for (auto & ns : neutral_species.items())
     478             :     {
     479             :       std::vector<std::vector<Real>> coeffs;
     480             : 
     481       11884 :       for (auto & nsac : ns.value().items())
     482             :       {
     483       16196 :         if (nsac.key() == "note")
     484         767 :           continue;
     485        7331 :         std::vector<Real> coeffvec(nsac.value().size());
     486             : 
     487       65979 :         for (unsigned int i = 0; i < coeffvec.size(); ++i)
     488       58648 :           coeffvec[i] = getReal(nsac.value()[i]);
     489             : 
     490        7331 :         coeffs.push_back(coeffvec);
     491             :       }
     492             : 
     493             :       // GeochemistryNeutralSpeciesActivity expects four vectos, so
     494             :       // add empty vectors if coeffs.size() != 4
     495        1893 :       coeffs.resize(4, {});
     496             : 
     497        1893 :       GeochemistryNeutralSpeciesActivity nsa(coeffs);
     498             : 
     499        1893 :       _neutral_species_activity[ns.key()] = nsa;
     500        1893 :     }
     501             :   }
     502         948 : }
     503             : 
     504             : const std::map<std::string, GeochemistryNeutralSpeciesActivity> &
     505        4836 : GeochemicalDatabaseReader::getNeutralSpeciesActivity() const
     506             : {
     507        4836 :   if (!_root["Header"].contains("neutral species"))
     508           1 :     mooseError("No neutral species activity coefficients in database");
     509        4835 :   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        1833 : GeochemicalDatabaseReader::secondarySpeciesNames() const
     707             : {
     708             :   std::vector<std::string> names;
     709        1833 :   if (_root.contains("secondary species"))
     710      479654 :     for (auto & item : _root["secondary species"].items())
     711      475990 :       names.push_back(item.key());
     712             : 
     713        1833 :   if (_root.contains("free electron"))
     714        4977 :     for (const auto & nm : _root["free electron"].items())
     715        1659 :       names.push_back(nm.key());
     716        1833 :   return names;
     717           0 : }
     718             : 
     719             : std::vector<std::string>
     720        1164 : GeochemicalDatabaseReader::redoxCoupleNames() const
     721             : {
     722             :   std::vector<std::string> names;
     723        1164 :   if (_root.contains("redox couples"))
     724       31532 :     for (const auto & item : _root["redox couples"].items())
     725       29204 :       names.push_back(item.key());
     726        1164 :   return names;
     727           0 : }
     728             : 
     729             : std::vector<std::string>
     730        1164 : GeochemicalDatabaseReader::surfaceSpeciesNames() const
     731             : {
     732             :   std::vector<std::string> names;
     733        1164 :   if (_root.contains("surface species"))
     734       12160 :     for (const auto & item : _root["surface species"].items())
     735       10922 :       names.push_back(item.key());
     736        1164 :   return names;
     737           0 : }
     738             : 
     739             : const FileName &
     740           2 : GeochemicalDatabaseReader::filename() const
     741             : {
     742           2 :   return _filename;
     743             : }
     744             : 
     745             : bool
     746        6984 : GeochemicalDatabaseReader::isBasisSpecies(const std::string & name) const
     747             : {
     748        6984 :   return _root["basis species"].contains(name);
     749             : }
     750             : 
     751             : bool
     752        3796 : GeochemicalDatabaseReader::isRedoxSpecies(const std::string & name) const
     753             : {
     754        3796 :   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     4720437 : GeochemicalDatabaseReader::getReal(const nlohmann::json & node)
     813             : {
     814     4720437 :   if (node.is_string())
     815      263510 :     return MooseUtils::convert<Real>(node);
     816     4588682 :   return node;
     817             : }