OGS
ChemistryLib::PhreeqcIOData Namespace Reference

Namespaces

 anonymous_namespace{PhreeqcIO.cpp}
 

Classes

class  PhreeqcIO
 
struct  Component
 
struct  AqueousSolution
 
struct  ChemicalSystem
 
struct  Dump
 
struct  EquilibriumReactant
 
struct  ExchangeSite
 
struct  KineticReactant
 
struct  Knobs
 
class  BasicOutputSetups
 
struct  OutputItem
 
struct  Output
 
struct  ReactionRate
 
struct  SurfaceSite
 
struct  SecondaryVariable
 
struct  UserPunch
 

Enumerations

enum class  ItemType {
  pH , pe , Component , EquilibriumReactant ,
  KineticReactant , SecondaryVariable
}
 

Functions

std::ostream & operator<< (std::ostream &os, PhreeqcIO const &phreeqc_io)
 
std::istream & operator>> (std::istream &in, PhreeqcIO &phreeqc_io)
 
std::unique_ptr< AqueousSolutioncreateAqueousSolution (BaseLib::ConfigTree const &config, MeshLib::Mesh &mesh)
 
std::unique_ptr< ChemicalSystemcreateChemicalSystem (BaseLib::ConfigTree const &config, MeshLib::Mesh &mesh)
 
std::vector< EquilibriumReactantcreateEquilibriumReactants (std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
 
std::vector< ExchangeSitecreateExchange (std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
 
std::vector< KineticReactantcreateKineticReactants (std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
 
Knobs createKnobs (BaseLib::ConfigTree const &config)
 
std::unique_ptr< OutputcreateOutput (ChemicalSystem const &chemical_system, std::unique_ptr< UserPunch > const &user_punch, bool const use_high_precision, std::string const &project_file_name)
 
std::vector< ComponentcreateSolutionComponents (BaseLib::ConfigTree const &config)
 
std::vector< SurfaceSitecreateSurface (std::optional< BaseLib::ConfigTree > const &config)
 
std::unique_ptr< UserPunchcreateUserPunch (std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh const &mesh)
 
std::ostream & operator<< (std::ostream &os, Knobs const &knobs)
 
std::ostream & operator<< (std::ostream &os, BasicOutputSetups const &basic_output_setups)
 
std::ostream & operator<< (std::ostream &os, Output const &output)
 
std::ostream & operator<< (std::ostream &os, ReactionRate const &reaction_rate)
 
std::ostream & operator<< (std::ostream &os, SurfaceSite const &surface_site)
 
std::ostream & operator<< (std::ostream &os, UserPunch const &user_punch)
 

Enumeration Type Documentation

◆ ItemType

Enumerator
pH 
pe 
Component 
EquilibriumReactant 
KineticReactant 
SecondaryVariable 

Definition at line 62 of file Output.h.

Function Documentation

◆ createAqueousSolution()

std::unique_ptr< AqueousSolution > ChemistryLib::PhreeqcIOData::createAqueousSolution ( BaseLib::ConfigTree const &  config,
MeshLib::Mesh mesh 
)
Input File Parameter:
prj__chemical_system__solution__temperature
Input File Parameter:
prj__chemical_system__solution__pressure
Input File Parameter:
prj__chemical_system__solution__pe

Definition at line 23 of file CreateAqueousSolution.cpp.

25 {
27  auto const temperature = config.getConfigParameter<double>("temperature");
28 
30  auto const pressure = config.getConfigParameter<double>("pressure");
31 
33  auto const pe0 = config.getConfigParameter<double>("pe");
34 
35  auto pe = MeshLib::getOrCreateMeshProperty<double>(
37 
38  auto components = createSolutionComponents(config);
39 
40  auto charge_balance = createChargeBalance(config);
41 
42  return std::make_unique<AqueousSolution>(
43  temperature, pressure, pe, pe0, std::move(components), charge_balance);
44 }
std::vector< Component > createSolutionComponents(BaseLib::ConfigTree const &config)
ChargeBalance createChargeBalance(BaseLib::ConfigTree const &config)

References ChemistryLib::createChargeBalance(), createSolutionComponents(), BaseLib::ConfigTree::getConfigParameter(), MeshLib::IntegrationPoint, and ChemistryLib::pe.

Referenced by createChemicalSystem().

◆ createChemicalSystem()

std::unique_ptr< ChemicalSystem > ChemistryLib::PhreeqcIOData::createChemicalSystem ( BaseLib::ConfigTree const &  config,
MeshLib::Mesh mesh 
)
Input File Parameter:
prj__chemical_system__solution
Input File Parameter:
prj__chemical_system__kinetic_reactants
Input File Parameter:
prj__chemical_system__equilibrium_reactants
Input File Parameter:
prj__chemical_system__exchangers

Definition at line 27 of file CreateChemicalSystem.cpp.

29 {
30  // solution
31  auto aqueous_solution = createAqueousSolution(
33  config.getConfigSubtree("solution"), mesh);
34 
35  // kinetic reactants
36  auto kinetic_reactants = createKineticReactants(
38  config.getConfigSubtreeOptional("kinetic_reactants"), mesh);
39 
40  // equilibrium reactants
41  auto equilibrium_reactants = createEquilibriumReactants(
43  config.getConfigSubtreeOptional("equilibrium_reactants"), mesh);
44 
45  // exchangers
46  auto exchangers = createExchange(
48  config.getConfigSubtreeOptional("exchangers"), mesh);
49 
50  return std::make_unique<ChemicalSystem>(std::move(aqueous_solution),
51  std::move(kinetic_reactants),
52  std::move(equilibrium_reactants),
53  std::move(exchangers));
54 }
std::unique_ptr< AqueousSolution > createAqueousSolution(BaseLib::ConfigTree const &config, MeshLib::Mesh &mesh)
std::vector< ExchangeSite > createExchange(std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
std::vector< EquilibriumReactant > createEquilibriumReactants(std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
std::vector< KineticReactant > createKineticReactants(std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)

References createAqueousSolution(), createEquilibriumReactants(), createExchange(), createKineticReactants(), BaseLib::ConfigTree::getConfigSubtree(), and BaseLib::ConfigTree::getConfigSubtreeOptional().

◆ createEquilibriumReactants()

std::vector< EquilibriumReactant > ChemistryLib::PhreeqcIOData::createEquilibriumReactants ( std::optional< BaseLib::ConfigTree > const &  config,
MeshLib::Mesh mesh 
)
Input File Parameter:
prj__chemical_system__equilibrium_reactants__phase_component
Input File Parameter:
prj__chemical_system__equilibrium_reactants__phase_component__name
Input File Parameter:
prj__chemical_system__equilibrium_reactants__phase_component__saturation_index
Input File Parameter:
prj__chemical_system__equilibrium_reactants__phase_component__reaction_irreversibility

Definition at line 23 of file CreateEquilibriumReactants.cpp.

25 {
26  if (!config)
27  {
28  return {};
29  }
30 
31  std::vector<EquilibriumReactant> equilibrium_reactants;
32  for (
33  auto const& equilibrium_reactant_config :
35  config->getConfigSubtreeList("phase_component"))
36  {
37  auto name =
39  equilibrium_reactant_config.getConfigParameter<std::string>("name");
40 
41  double const saturation_index =
43  equilibrium_reactant_config.getConfigParameter<double>(
44  "saturation_index");
45 
46  auto reaction_irreversibility =
48  equilibrium_reactant_config.getConfigParameter<std::string>(
49  "reaction_irreversibility", "");
50 
51  if (!reaction_irreversibility.empty() &&
52  (reaction_irreversibility != "dissolve_only" &&
53  reaction_irreversibility != "precipitate_only"))
54  {
55  OGS_FATAL(
56  "{:s}: reaction direction only allows `dissolve_only` or "
57  "`precipitate_only`",
58  name);
59  }
60 
61  auto molality = MeshLib::getOrCreateMeshProperty<double>(
63 
64  auto molality_prev = MeshLib::getOrCreateMeshProperty<double>(
65  mesh, name + "_prev", MeshLib::MeshItemType::IntegrationPoint, 1);
66 
67  auto volume_fraction = MeshLib::getOrCreateMeshProperty<double>(
69 
70  auto volume_fraction_prev = MeshLib::getOrCreateMeshProperty<double>(
71  mesh,
72  "phi_" + name + "_prev",
74  1);
75 
76  auto mesh_prop_molality = MeshLib::getOrCreateMeshProperty<double>(
77  mesh, name + "_avg", MeshLib::MeshItemType::Cell, 1);
78  mesh_prop_molality->resize(mesh.getNumberOfElements());
79 
80  equilibrium_reactants.emplace_back(std::move(name),
81  molality,
82  molality_prev,
84  volume_fraction_prev,
85  mesh_prop_molality,
86  saturation_index,
87  std::move(reaction_irreversibility));
88  }
89 
90  return equilibrium_reactants;
91 }
#define OGS_FATAL(...)
Definition: Error.h:26
std::size_t getNumberOfElements() const
Get the number of elements.
Definition: Mesh.h:86

References MeshLib::Cell, MeshLib::Mesh::getNumberOfElements(), MeshLib::IntegrationPoint, MaterialPropertyLib::molality, MaterialPropertyLib::name, OGS_FATAL, and MaterialPropertyLib::volume_fraction.

Referenced by createChemicalSystem().

◆ createExchange()

std::vector< ExchangeSite > ChemistryLib::PhreeqcIOData::createExchange ( std::optional< BaseLib::ConfigTree > const &  config,
MeshLib::Mesh mesh 
)
Input File Parameter:
prj__chemical_system__exchangers__exchange_site
Input File Parameter:
prj__chemical_system__exchangers__exchange_site__ion_exchanging_species

Definition at line 21 of file CreateExchange.cpp.

23 {
24  if (!config)
25  {
26  return {};
27  }
28 
29  std::vector<ExchangeSite> exchangers;
30  for (auto const& site_config :
32  config->getConfigSubtreeList("exchange_site"))
33  {
35  auto name = site_config.getConfigParameter<std::string>(
36  "ion_exchanging_species");
37 
38  auto const molality = MeshLib::getOrCreateMeshProperty<double>(
40 
41  exchangers.emplace_back(std::move(name), molality);
42  }
43 
44  return exchangers;
45 }

References MeshLib::IntegrationPoint, MaterialPropertyLib::molality, and MaterialPropertyLib::name.

Referenced by createChemicalSystem().

◆ createKineticReactants()

std::vector< KineticReactant > ChemistryLib::PhreeqcIOData::createKineticReactants ( std::optional< BaseLib::ConfigTree > const &  config,
MeshLib::Mesh mesh 
)
Input File Parameter:
prj__chemical_system__kinetic_reactants__kinetic_reactant
Input File Parameter:
prj__chemical_system__kinetic_reactants__kinetic_reactant__name
Input File Parameter:
prj__chemical_system__kinetic_reactants__kinetic_reactant__chemical_formula
Input File Parameter:
prj__chemical_system__kinetic_reactants__kinetic_reactant__parameters
Input File Parameter:
prj__chemical_system__kinetic_reactants__kinetic_reactant__fix_amount

Definition at line 23 of file CreateKineticReactant.cpp.

25 {
26  if (!config)
27  {
28  return {};
29  }
30 
31  std::vector<KineticReactant> kinetic_reactants;
32  for (
33  auto const& reactant_config :
35  config->getConfigSubtreeList("kinetic_reactant"))
36  {
38  auto name = reactant_config.getConfigParameter<std::string>("name");
39 
40  auto chemical_formula =
42  reactant_config.getConfigParameter<std::string>("chemical_formula",
43  "");
44 
45  auto parameters =
47  reactant_config.getConfigParameter<std::vector<double>>(
48  "parameters", {});
49 
50  bool const fix_amount =
52  reactant_config.getConfigParameter<bool>("fix_amount", false);
53 
54  auto molality = MeshLib::getOrCreateMeshProperty<double>(
56 
57  auto molality_prev = MeshLib::getOrCreateMeshProperty<double>(
58  mesh, name + "_prev", MeshLib::MeshItemType::IntegrationPoint, 1);
59 
60  auto volume_fraction = MeshLib::getOrCreateMeshProperty<double>(
62 
63  auto volume_fraction_prev = MeshLib::getOrCreateMeshProperty<double>(
64  mesh,
65  "phi_" + name + "_prev",
67  1);
68 
69  auto mesh_prop_molality = MeshLib::getOrCreateMeshProperty<double>(
70  mesh, name + "_avg", MeshLib::MeshItemType::Cell, 1);
71  mesh_prop_molality->resize(mesh.getNumberOfElements());
72 
73  if (chemical_formula.empty() && fix_amount)
74  {
75  OGS_FATAL(
76  "fix_amount can only be used if a chemical_formula has been "
77  "defined");
78  }
79 
80  kinetic_reactants.emplace_back(std::move(name),
81  std::move(chemical_formula),
82  molality,
83  molality_prev,
85  volume_fraction_prev,
86  mesh_prop_molality,
87  std::move(parameters),
88  fix_amount);
89  }
90 
91  return kinetic_reactants;
92 }

References MeshLib::Cell, MeshLib::Mesh::getNumberOfElements(), MeshLib::IntegrationPoint, MaterialPropertyLib::molality, MaterialPropertyLib::name, OGS_FATAL, and MaterialPropertyLib::volume_fraction.

Referenced by createChemicalSystem().

◆ createKnobs()

Knobs ChemistryLib::PhreeqcIOData::createKnobs ( BaseLib::ConfigTree const &  config)
Input File Parameter:
prj__chemical_system__knobs__max_iter
Input File Parameter:
prj__chemical_system__knobs__relative_convergence_tolerance
Input File Parameter:
prj__chemical_system__knobs__tolerance
Input File Parameter:
prj__chemical_system__knobs__step_size
Input File Parameter:
prj__chemical_system__knobs__scaling

Definition at line 20 of file CreateKnobs.cpp.

21 {
22  auto const max_iterations =
24  config.getConfigParameter<int>("max_iter");
25 
26  auto const relative_convergence_tolerance =
28  config.getConfigParameter<double>("relative_convergence_tolerance");
29 
30  auto const tolerance =
32  config.getConfigParameter<double>("tolerance");
33 
34  auto const step_size =
36  config.getConfigParameter<int>("step_size");
37 
38  auto const scaling =
40  config.getConfigParameter<bool>("scaling");
41 
42  return {max_iterations, relative_convergence_tolerance, tolerance,
43  step_size, scaling};
44 }
static double const tolerance

References BaseLib::ConfigTree::getConfigParameter(), and ParameterLib::tolerance.

◆ createOutput()

std::unique_ptr< Output > ChemistryLib::PhreeqcIOData::createOutput ( ChemicalSystem const &  chemical_system,
std::unique_ptr< UserPunch > const &  user_punch,
bool const  use_high_precision,
std::string const &  project_file_name 
)

Definition at line 25 of file CreateOutput.cpp.

30 {
31  auto const& components = chemical_system.aqueous_solution->components;
32  auto const& equilibrium_reactants = chemical_system.equilibrium_reactants;
33  auto const& kinetic_reactants = chemical_system.kinetic_reactants;
34  // Mark which phreeqc output items will be held.
35  std::vector<OutputItem> accepted_items{{"pH", ItemType::pH},
36  {"pe", ItemType::pe}};
37 
38  auto accepted_item = [](auto const& item)
39  { return OutputItem(item.name, item.item_type); };
40  std::transform(components.begin(), components.end(),
41  std::back_inserter(accepted_items), accepted_item);
42  std::transform(equilibrium_reactants.begin(), equilibrium_reactants.end(),
43  std::back_inserter(accepted_items), accepted_item);
44  for (auto const& kinetic_reactant : kinetic_reactants)
45  {
46  if (kinetic_reactant.fix_amount)
47  {
48  continue;
49  }
50  accepted_items.emplace_back(kinetic_reactant.name,
51  kinetic_reactant.item_type);
52  }
53 
54  if (user_punch)
55  {
56  auto const& secondary_variables = user_punch->secondary_variables;
57  accepted_items.reserve(accepted_items.size() +
58  secondary_variables.size());
59  std::transform(secondary_variables.begin(), secondary_variables.end(),
60  std::back_inserter(accepted_items), accepted_item);
61  }
62 
63  // Record ids of which phreeqc output items will be dropped.
64  BasicOutputSetups basic_output_setups(project_file_name,
65  use_high_precision);
66  auto const num_dropped_basic_items =
67  basic_output_setups.getNumberOfDroppedItems();
68  std::vector<int> dropped_item_ids(num_dropped_basic_items);
69  std::iota(dropped_item_ids.begin(), dropped_item_ids.end(), 0);
70 
71  auto const num_dvalue_equilibrium_reactants = equilibrium_reactants.size();
72  auto const num_dvalue_kinetic_reactants = kinetic_reactants.size();
73  int const num_dvalue_items =
74  num_dvalue_equilibrium_reactants + num_dvalue_kinetic_reactants;
75 
76  auto const num_basic_items =
77  basic_output_setups.getNumberOfItemsInDisplay();
78  auto const num_components = components.size();
79  auto dvalue_item_id = num_basic_items + num_components + 1;
80  for (int i = 0; i < num_dvalue_items; ++i)
81  {
82  dropped_item_ids.push_back(dvalue_item_id);
83  dvalue_item_id += 2;
84  }
85 
86  return std::make_unique<Output>(std::move(basic_output_setups),
87  std::move(accepted_items),
88  std::move(dropped_item_ids));
89 }

References ChemistryLib::PhreeqcIOData::ChemicalSystem::aqueous_solution, ChemistryLib::PhreeqcIOData::ChemicalSystem::equilibrium_reactants, ChemistryLib::PhreeqcIOData::BasicOutputSetups::getNumberOfDroppedItems(), ChemistryLib::PhreeqcIOData::BasicOutputSetups::getNumberOfItemsInDisplay(), ChemistryLib::PhreeqcIOData::ChemicalSystem::kinetic_reactants, pe, and pH.

◆ createSolutionComponents()

std::vector< Component > ChemistryLib::PhreeqcIOData::createSolutionComponents ( BaseLib::ConfigTree const &  config)
Input File Parameter:
prj__chemical_system__solution__components
Input File Parameter:
prj__chemical_system__solution__components__component
Input File Parameter:
prj__chemical_system__solution__components__component__chemical_formula

Definition at line 20 of file CreateSolutionComponent.cpp.

22 {
23  std::vector<Component> components;
25  auto components_config = config.getConfigSubtree("components");
26 
27  for (
28  auto const& component_config :
30  components_config.getConfigSubtreeList("component"))
31  {
32  auto const component_name = component_config.getValue<std::string>();
33  auto const chemical_formula =
35  component_config.getConfigAttribute<std::string>("chemical_formula",
36  "");
37  components.emplace_back(component_name, chemical_formula);
38  }
39 
40  return components;
41 }

References BaseLib::ConfigTree::getConfigSubtree(), and BaseLib::ConfigTree::getValue().

Referenced by createAqueousSolution().

◆ createSurface()

std::vector< SurfaceSite > ChemistryLib::PhreeqcIOData::createSurface ( std::optional< BaseLib::ConfigTree > const &  config)
Input File Parameter:
prj__chemical_system__surface__site
Input File Parameter:
prj__chemical_system__surface__site__name
Input File Parameter:
prj__chemical_system__surface__site__site_density
Input File Parameter:
prj__chemical_system__surface__site__specific_surface_area
Input File Parameter:
prj__chemical_system__surface__site__mass

Definition at line 20 of file CreateSurface.cpp.

22 {
23  if (!config)
24  {
25  return {};
26  }
27 
28  std::vector<SurfaceSite> surface;
29  for (auto const& site_config :
31  config->getConfigSubtreeList("site"))
32  {
34  auto name = site_config.getConfigParameter<std::string>("name");
35 
36  auto const site_density =
38  site_config.getConfigParameter<double>("site_density");
39 
40  auto const specific_surface_area =
42  site_config.getConfigParameter<double>("specific_surface_area");
43 
44  auto const mass =
46  site_config.getConfigParameter<double>("mass");
47 
48  surface.emplace_back(
49  std::move(name), site_density, specific_surface_area, mass);
50  }
51 
52  return surface;
53 }

References MaterialPropertyLib::name.

◆ createUserPunch()

std::unique_ptr< UserPunch > ChemistryLib::PhreeqcIOData::createUserPunch ( std::optional< BaseLib::ConfigTree > const &  config,
MeshLib::Mesh const &  mesh 
)
Input File Parameter:
prj__chemical_system__user_punch__headline
Input File Parameter:
prj__chemical_system__user_punch__statement

Definition at line 21 of file CreateUserPunch.cpp.

23 {
24  if (!config)
25  {
26  return nullptr;
27  }
28 
29  std::vector<SecondaryVariable> secondary_variables;
30  for (auto const& variable_name :
32  config->getConfigParameter<std::vector<std::string>>("headline"))
33  {
34  auto value = MeshLib::getOrCreateMeshProperty<double>(
35  const_cast<MeshLib::Mesh&>(mesh),
36  variable_name,
38  1);
39  std::fill(std::begin(*value),
40  std::end(*value),
41  std::numeric_limits<double>::quiet_NaN());
42 
43  secondary_variables.emplace_back(variable_name, value);
44  }
45 
46  std::vector<std::string> statements;
47  for (auto const& statement :
49  config->getConfigParameterList<std::string>("statement"))
50  {
51  statements.emplace_back(statement);
52  }
53 
54  return std::make_unique<UserPunch>(std::move(secondary_variables),
55  std::move(statements));
56 }

References MeshLib::IntegrationPoint.

◆ operator<<() [1/7]

std::ostream& ChemistryLib::PhreeqcIOData::operator<< ( std::ostream &  os,
BasicOutputSetups const &  basic_output_setups 
)

Definition at line 19 of file Output.cpp.

21 {
22  os << "SELECTED_OUTPUT"
23  << "\n";
24  os << "-file " << basic_output_setups.output_file << "\n";
25  os << "-high_precision " << std::boolalpha
26  << basic_output_setups.use_high_precision << "\n";
27  os << "-simulation " << std::boolalpha
28  << BasicOutputSetups::display_simulation_id << "\n";
29  os << "-state " << std::boolalpha << BasicOutputSetups::display_state
30  << "\n";
31  os << "-distance " << std::boolalpha << BasicOutputSetups::display_distance
32  << "\n";
33  os << "-time " << std::boolalpha << BasicOutputSetups::display_current_time
34  << "\n";
35  os << "-step " << std::boolalpha << BasicOutputSetups::display_time_step
36  << "\n";
37 
38  return os;
39 }

◆ operator<<() [2/7]

std::ostream& ChemistryLib::PhreeqcIOData::operator<< ( std::ostream &  os,
Knobs const &  knobs 
)

Definition at line 19 of file Knobs.cpp.

20 {
21  os << "KNOBS"
22  << "\n";
23  os << "-iterations " << knobs.max_iterations << "\n";
24  os << "-convergence_tolerance " << knobs.relative_convergence_tolerance
25  << "\n";
26  os << "-tolerance " << knobs.tolerance << "\n";
27  os << "-step_size " << knobs.step_size << "\n";
28  os << "-diagonal_scale " << knobs.scaling << "\n";
29 
30  return os;
31 }

◆ operator<<() [3/7]

std::ostream& ChemistryLib::PhreeqcIOData::operator<< ( std::ostream &  os,
Output const &  output 
)

Definition at line 41 of file Output.cpp.

42 {
43  os << output.basic_output_setups;
44 
45  auto const component_items =
46  output.getOutputItemsByItemType(ItemType::Component);
47  os << "-totals";
48  for (auto const& component_item : component_items)
49  {
50  os << " " << component_item.name;
51  }
52  os << "\n";
53 
54  auto const equilibrium_phase_items =
55  output.getOutputItemsByItemType(ItemType::EquilibriumReactant);
56  if (!equilibrium_phase_items.empty())
57  {
58  os << "-equilibrium_phases";
59  for (auto const& equilibrium_phase_item : equilibrium_phase_items)
60  {
61  os << " " << equilibrium_phase_item.name;
62  }
63  os << "\n";
64  }
65 
66  auto const kinetic_reactant_items =
67  output.getOutputItemsByItemType(ItemType::KineticReactant);
68  if (!kinetic_reactant_items.empty())
69  {
70  os << "-kinetic_reactants";
71  for (auto const& kinetic_reactant_item : kinetic_reactant_items)
72  {
73  os << " " << kinetic_reactant_item.name;
74  }
75  os << "\n";
76  }
77 
78  return os;
79 }

◆ operator<<() [4/7]

std::ostream& ChemistryLib::PhreeqcIOData::operator<< ( std::ostream &  os,
PhreeqcIO const &  phreeqc_io 
)

Definition at line 460 of file PhreeqcIO.cpp.

461 {
462  os << phreeqc_io._knobs << "\n";
463 
464  os << *phreeqc_io._output << "\n";
465 
466  auto const& user_punch = phreeqc_io._user_punch;
467  if (user_punch)
468  {
469  os << *user_punch << "\n";
470  }
471 
472  auto const& reaction_rates = phreeqc_io._reaction_rates;
473  if (!reaction_rates.empty())
474  {
475  os << "RATES"
476  << "\n";
477  os << reaction_rates << "\n";
478  }
479 
480  for (std::size_t chemical_system_id = 0;
481  chemical_system_id < phreeqc_io._num_chemical_systems;
482  ++chemical_system_id)
483  {
484  os << "SOLUTION " << chemical_system_id + 1 << "\n";
485  phreeqc_io._chemical_system->aqueous_solution->print(
486  os, chemical_system_id);
487 
488  auto const& dump = phreeqc_io._dump;
489  if (dump)
490  {
491  auto const& aqueous_solutions_prev = dump->aqueous_solutions_prev;
492  if (!aqueous_solutions_prev.empty())
493  {
494  os << aqueous_solutions_prev[chemical_system_id] << "\n\n";
495  }
496  }
497 
498  os << "USE solution none"
499  << "\n";
500  os << "END"
501  << "\n\n";
502 
503  os << "USE solution " << chemical_system_id + 1 << "\n\n";
504 
505  auto const& equilibrium_reactants =
506  phreeqc_io._chemical_system->equilibrium_reactants;
507  if (!equilibrium_reactants.empty())
508  {
509  os << "EQUILIBRIUM_PHASES " << chemical_system_id + 1 << "\n";
510  for (auto const& equilibrium_reactant : equilibrium_reactants)
511  {
512  equilibrium_reactant.print(os, chemical_system_id);
513  }
514  os << "\n";
515  }
516 
517  auto const& kinetic_reactants =
518  phreeqc_io._chemical_system->kinetic_reactants;
519  if (!kinetic_reactants.empty())
520  {
521  os << "KINETICS " << chemical_system_id + 1 << "\n";
522  for (auto const& kinetic_reactant : kinetic_reactants)
523  {
524  kinetic_reactant.print(os, chemical_system_id);
525  }
526  os << "-steps " << phreeqc_io._dt << "\n"
527  << "\n";
528  }
529 
530  auto const& surface = phreeqc_io._surface;
531  if (!surface.empty())
532  {
533  os << "SURFACE " << chemical_system_id + 1 << "\n";
534  std::size_t aqueous_solution_id =
535  dump->aqueous_solutions_prev.empty()
536  ? chemical_system_id + 1
537  : phreeqc_io._num_chemical_systems + chemical_system_id + 1;
538  os << "-equilibrate with solution " << aqueous_solution_id << "\n";
539  os << "-sites_units DENSITY"
540  << "\n";
541  os << surface << "\n";
542  os << "SAVE solution " << chemical_system_id + 1 << "\n";
543  }
544 
545  auto const& exchangers = phreeqc_io._chemical_system->exchangers;
546  if (!exchangers.empty())
547  {
548  os << "EXCHANGE " << chemical_system_id + 1 << "\n";
549  std::size_t const aqueous_solution_id =
550  dump->aqueous_solutions_prev.empty()
551  ? chemical_system_id + 1
552  : phreeqc_io._num_chemical_systems + chemical_system_id + 1;
553  os << "-equilibrate with solution " << aqueous_solution_id << "\n";
554  for (auto const& exchanger : exchangers)
555  {
556  exchanger.print(os, chemical_system_id);
557  }
558  os << "SAVE solution " << chemical_system_id + 1 << "\n";
559  }
560 
561  os << "END"
562  << "\n\n";
563  }
564 
565  auto const& dump = phreeqc_io._dump;
566  if (dump)
567  {
568  dump->print(os, phreeqc_io._num_chemical_systems);
569  }
570 
571  return os;
572 }

◆ operator<<() [5/7]

std::ostream& ChemistryLib::PhreeqcIOData::operator<< ( std::ostream &  os,
ReactionRate const &  reaction_rate 
)

Definition at line 19 of file ReactionRate.cpp.

20 {
21  os << reaction_rate.kinetic_reactant << "\n";
22  os << "-start"
23  << "\n";
24  int line_number = 1;
25  for (auto const& expression_statement : reaction_rate.expression_statements)
26  {
27  os << line_number << " " << expression_statement << "\n";
28  ++line_number;
29  }
30  os << "-end"
31  << "\n";
32 
33  return os;
34 }

◆ operator<<() [6/7]

std::ostream& ChemistryLib::PhreeqcIOData::operator<< ( std::ostream &  os,
SurfaceSite const &  surface_site 
)

Definition at line 19 of file Surface.cpp.

20 {
21  os << surface_site.name << " " << surface_site.site_density << " "
22  << surface_site.specific_surface_area << " " << surface_site.mass
23  << "\n";
24 
25  return os;
26 }

◆ operator<<() [7/7]

std::ostream& ChemistryLib::PhreeqcIOData::operator<< ( std::ostream &  os,
UserPunch const &  user_punch 
)

Definition at line 27 of file UserPunch.cpp.

28 {
29  os << "USER_PUNCH"
30  << "\n";
31  os << "-headings ";
32  auto const& secondary_variables = user_punch.secondary_variables;
33  for (auto& secondary_variable : secondary_variables)
34  {
35  os << secondary_variable.name << " ";
36  }
37  os << "\n";
38 
39  os << "-start"
40  << "\n";
41  int line_number = 1;
42  for (auto const& statement : user_punch.statements)
43  {
44  os << line_number << " " << statement << "\n";
45  ++line_number;
46  }
47  os << "-end"
48  << "\n";
49 
50  return os;
51 }

◆ operator>>()

std::istream& ChemistryLib::PhreeqcIOData::operator>> ( std::istream &  in,
PhreeqcIO phreeqc_io 
)

Definition at line 611 of file PhreeqcIO.cpp.

612 {
613  // Skip the headline
614  in.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
615 
616  std::string line;
617  auto const& output = *phreeqc_io._output;
618  auto const& dropped_item_ids = output.dropped_item_ids;
619 
620  auto const& surface = phreeqc_io._surface;
621  auto const& exchangers = phreeqc_io._chemical_system->exchangers;
622 
623  if (!surface.empty() && !exchangers.empty())
624  {
625  OGS_FATAL(
626  "Using surface and exchange reactions simultaneously is not "
627  "supported at the moment");
628  }
629 
630  int const num_skipped_lines = surface.empty() && exchangers.empty() ? 1 : 2;
631 
632  auto& equilibrium_reactants =
633  phreeqc_io._chemical_system->equilibrium_reactants;
634  auto& kinetic_reactants = phreeqc_io._chemical_system->kinetic_reactants;
635 
636  for (std::size_t chemical_system_id = 0;
637  chemical_system_id < phreeqc_io._num_chemical_systems;
638  ++chemical_system_id)
639  {
640  // Skip equilibrium calculation result of initial solution
641  for (int i = 0; i < num_skipped_lines; ++i)
642  {
643  in.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
644  }
645 
646  // Get calculation result of the solution after the reaction
647  if (!std::getline(in, line))
648  {
649  OGS_FATAL(
650  "Error when reading calculation result of Solution {:d} "
651  "after the reaction.",
652  chemical_system_id);
653  }
654 
655  std::vector<double> accepted_items;
656  std::vector<std::string> items;
657  boost::trim_if(line, boost::is_any_of("\t "));
658  boost::algorithm::split(items, line, boost::is_any_of("\t "),
659  boost::token_compress_on);
660  for (int item_id = 0; item_id < static_cast<int>(items.size());
661  ++item_id)
662  {
663  if (std::find(dropped_item_ids.begin(), dropped_item_ids.end(),
664  item_id) == dropped_item_ids.end())
665  {
666  double value;
667  try
668  {
669  value = std::stod(items[item_id]);
670  }
671  catch (const std::invalid_argument& e)
672  {
673  OGS_FATAL(
674  "Invalid argument. Could not convert string '{:s}' to "
675  "double for chemical system {:d}, column {:d}. "
676  "Exception '{:s}' was thrown.",
677  items[item_id], chemical_system_id + 1, item_id,
678  e.what());
679  }
680  catch (const std::out_of_range& e)
681  {
682  OGS_FATAL(
683  "Out of range error. Could not convert string "
684  "'{:s}' to double for chemical system {:d}, column "
685  "{:d}. Exception '{:s}' was thrown.",
686  items[item_id], chemical_system_id + 1, item_id,
687  e.what());
688  }
689  accepted_items.push_back(value);
690  }
691  }
692  assert(accepted_items.size() == output.accepted_items.size());
693 
694  auto& aqueous_solution = phreeqc_io._chemical_system->aqueous_solution;
695  auto& components = aqueous_solution->components;
696  auto& user_punch = phreeqc_io._user_punch;
697  for (int item_id = 0; item_id < static_cast<int>(accepted_items.size());
698  ++item_id)
699  {
700  auto const& accepted_item = output.accepted_items[item_id];
701  auto const& item_name = accepted_item.name;
702 
703  auto compare_by_name = [&item_name](auto const& item)
704  { return item.name == item_name; };
705 
706  switch (accepted_item.item_type)
707  {
708  case ItemType::pH:
709  {
710  // Update pH value
711  aqueous_solution->pH->set(
712  chemical_system_id,
713  std::pow(10, -accepted_items[item_id]));
714  break;
715  }
716  case ItemType::pe:
717  {
718  // Update pe value
719  (*aqueous_solution->pe)[chemical_system_id] =
720  accepted_items[item_id];
721  break;
722  }
723  case ItemType::Component:
724  {
725  // Update component concentrations
726  auto& component = BaseLib::findElementOrError(
727  components.begin(), components.end(), compare_by_name,
728  "Could not find component '" + item_name + "'.");
729  component.amount->set(chemical_system_id,
730  accepted_items[item_id]);
731  break;
732  }
733  case ItemType::EquilibriumReactant:
734  {
735  // Update amounts of equilibrium reactant
736  auto& equilibrium_reactant = BaseLib::findElementOrError(
737  equilibrium_reactants.begin(),
738  equilibrium_reactants.end(), compare_by_name,
739  "Could not find equilibrium reactant '" + item_name +
740  "'.");
741  (*equilibrium_reactant.molality)[chemical_system_id] =
742  accepted_items[item_id];
743  break;
744  }
745  case ItemType::KineticReactant:
746  {
747  // Update amounts of kinetic reactants
748  auto& kinetic_reactant = BaseLib::findElementOrError(
749  kinetic_reactants.begin(), kinetic_reactants.end(),
750  compare_by_name,
751  "Could not find kinetic reactant '" + item_name + "'.");
752  (*kinetic_reactant.molality)[chemical_system_id] =
753  accepted_items[item_id];
754  break;
755  }
756  case ItemType::SecondaryVariable:
757  {
758  assert(user_punch);
759  auto& secondary_variables = user_punch->secondary_variables;
760  // Update values of secondary variables
761  auto& secondary_variable = BaseLib::findElementOrError(
762  secondary_variables.begin(), secondary_variables.end(),
763  compare_by_name,
764  "Could not find secondary variable '" + item_name +
765  "'.");
766  (*secondary_variable.value)[chemical_system_id] =
767  accepted_items[item_id];
768  break;
769  }
770  }
771  }
772  }
773 
774  return in;
775 }
std::iterator_traits< InputIt >::reference findElementOrError(InputIt begin, InputIt end, Predicate predicate, std::string const &error="")
Definition: Algorithm.h:69