OGS
Namespaces | Classes | Enumerations | Functions
ChemistryLib::PhreeqcIOData Namespace Reference

Namespaces

namespace  anonymous_namespace{PhreeqcIO.cpp}
 

Classes

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

Enumerations

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

Functions

std::string specifyFileName (std::string const &project_file_name, std::string const &file_extension)
 
std::ostream & operator<< (std::ostream &os, PhreeqcIO const &phreeqc_io)
 
std::istream & operator>> (std::istream &in, PhreeqcIO &phreeqc_io)
 
std::unique_ptr< AqueousSolution > createAqueousSolution (BaseLib::ConfigTree const &config, MeshLib::Mesh &mesh)
 
std::unique_ptr< ChemicalSystem > createChemicalSystem (BaseLib::ConfigTree const &config, MeshLib::Mesh &mesh)
 
std::vector< EquilibriumReactant > createEquilibriumReactants (std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
 
std::vector< ExchangeSite > createExchange (std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
 
std::vector< KineticReactant > createKineticReactants (std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
 
Knobs createKnobs (BaseLib::ConfigTree const &config)
 
std::unique_ptr< Output > createOutput (ChemicalSystem const &chemical_system, std::unique_ptr< UserPunch > const &user_punch, bool const use_high_precision, std::string const &project_file_name)
 
std::vector< Component > createSolutionComponents (BaseLib::ConfigTree const &config)
 
std::vector< std::variant< DensityBasedSurfaceSite, MoleBasedSurfaceSite > > createSurface (std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
 
std::unique_ptr< UserPunch > createUserPunch (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, UserPunch const &user_punch)
 

Enumeration Type Documentation

◆ ItemType

enum class ChemistryLib::PhreeqcIOData::ItemType
strong
Enumerator
pH 
pe 
Component 
EquilibriumReactant 
KineticReactant 
SecondaryVariable 

Definition at line 59 of file Output.h.

60{
61 pH,
62 pe,
63 Component,
64 EquilibriumReactant,
65 KineticReactant,
66 SecondaryVariable
67};

Function Documentation

◆ createAqueousSolution()

std::unique_ptr< AqueousSolution > ChemistryLib::PhreeqcIOData::createAqueousSolution ( BaseLib::ConfigTree const & config,
MeshLib::Mesh & mesh )
Input File Parameter
prj__chemical_system__solution__fixing_pe
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 24 of file CreateAqueousSolution.cpp.

26{
28 auto const fixing_pe = config.getConfigAttribute<bool>("fixing_pe", false);
29
31 auto const temperature = config.getConfigParameter<double>("temperature");
32
34 auto const pressure = config.getConfigParameter<double>("pressure");
35
37 auto const pe0 = config.getConfigParameter<double>("pe");
38
39 auto pe = MeshLib::getOrCreateMeshProperty<double>(
40 mesh, "pe", MeshLib::MeshItemType::IntegrationPoint, 1);
41
42 auto components = createSolutionComponents(config);
43
44 auto charge_balance = createChargeBalance(config);
45
46 return std::make_unique<AqueousSolution>(fixing_pe, temperature, pressure,
47 pe, pe0, std::move(components),
48 charge_balance);
49}
ChemistryLib::PhreeqcIOData::createSolutionComponents
std::vector< Component > createSolutionComponents(BaseLib::ConfigTree const &config)
Definition CreateSolutionComponent.cpp:20
ChemistryLib::createChargeBalance
ChargeBalance createChargeBalance(BaseLib::ConfigTree const &config)
Definition CreateChargeBalance.cpp:19

References ChemistryLib::createChargeBalance(), createSolutionComponents(), BaseLib::ConfigTree::getConfigAttribute(), BaseLib::ConfigTree::getConfigParameter(), MeshLib::IntegrationPoint, and 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
Input File Parameter
prj__chemical_system__surface

Definition at line 29 of file CreateChemicalSystem.cpp.

31{
32 // solution
33 auto aqueous_solution = createAqueousSolution(
35 config.getConfigSubtree("solution"), mesh);
36
37 // kinetic reactants
38 auto kinetic_reactants = createKineticReactants(
40 config.getConfigSubtreeOptional("kinetic_reactants"), mesh);
41
42 // equilibrium reactants
43 auto equilibrium_reactants = createEquilibriumReactants(
45 config.getConfigSubtreeOptional("equilibrium_reactants"), mesh);
46
47 // exchangers
48 auto exchangers = createExchange(
50 config.getConfigSubtreeOptional("exchangers"), mesh);
51
52 // surface
53 auto surface = createSurface(
55 config.getConfigSubtreeOptional("surface"), mesh);
56
57 return std::make_unique<ChemicalSystem>(std::move(aqueous_solution),
58 std::move(kinetic_reactants),
59 std::move(equilibrium_reactants),
60 std::move(exchangers),
61 std::move(surface));
62}
ChemistryLib::PhreeqcIOData::createAqueousSolution
std::unique_ptr< AqueousSolution > createAqueousSolution(BaseLib::ConfigTree const &config, MeshLib::Mesh &mesh)
Definition CreateAqueousSolution.cpp:24
ChemistryLib::PhreeqcIOData::createExchange
std::vector< ExchangeSite > createExchange(std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
Definition CreateExchange.cpp:22
ChemistryLib::PhreeqcIOData::createEquilibriumReactants
std::vector< EquilibriumReactant > createEquilibriumReactants(std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
Definition CreateEquilibriumReactants.cpp:24
ChemistryLib::PhreeqcIOData::createKineticReactants
std::vector< KineticReactant > createKineticReactants(std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
Definition CreateKineticReactant.cpp:24

References createAqueousSolution(), createEquilibriumReactants(), createExchange(), createKineticReactants(), createSurface(), 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 24 of file CreateEquilibriumReactants.cpp.

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

References MeshLib::Cell, MeshLib::Mesh::getNumberOfElements(), MeshLib::IntegrationPoint, and OGS_FATAL.

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 22 of file CreateExchange.cpp.

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

References MeshLib::IntegrationPoint.

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 24 of file CreateKineticReactant.cpp.

26{
27 if (!config)
28 {
29 return {};
30 }
31
32 std::vector<KineticReactant> kinetic_reactants;
33 for (
34 auto const& reactant_config :
36 config->getConfigSubtreeList("kinetic_reactant"))
37 {
39 auto name = reactant_config.getConfigParameter<std::string>("name");
40
41 auto chemical_formula =
43 reactant_config.getConfigParameter<std::string>("chemical_formula",
44 "");
45
46 auto parameters =
48 reactant_config.getConfigParameter<std::vector<double>>(
49 "parameters", {});
50
51 bool const fix_amount =
53 reactant_config.getConfigParameter<bool>("fix_amount", false);
54
55 auto molality = MeshLib::getOrCreateMeshProperty<double>(
56 mesh, name, MeshLib::MeshItemType::IntegrationPoint, 1);
57
58 auto molality_prev = MeshLib::getOrCreateMeshProperty<double>(
59 mesh, name + "_prev", MeshLib::MeshItemType::IntegrationPoint, 1);
60
61 auto volume_fraction = MeshLib::getOrCreateMeshProperty<double>(
62 mesh, "phi_" + name, MeshLib::MeshItemType::IntegrationPoint, 1);
63
64 auto volume_fraction_prev = MeshLib::getOrCreateMeshProperty<double>(
65 mesh,
66 "phi_" + name + "_prev",
67 MeshLib::MeshItemType::IntegrationPoint,
68 1);
69
70 auto mesh_prop_molality = MeshLib::getOrCreateMeshProperty<double>(
71 mesh, name + "_avg", MeshLib::MeshItemType::Cell, 1);
72 mesh_prop_molality->resize(mesh.getNumberOfElements());
73
74 kinetic_reactants.emplace_back(std::move(name),
75 std::move(chemical_formula),
76 molality,
77 molality_prev,
78 volume_fraction,
79 volume_fraction_prev,
80 mesh_prop_molality,
81 std::move(parameters),
82 fix_amount);
83 }
84
85 return kinetic_reactants;
86}

References MeshLib::Cell, MeshLib::Mesh::getNumberOfElements(), and MeshLib::IntegrationPoint.

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}

References BaseLib::ConfigTree::getConfigParameter().

◆ 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< std::variant< DensityBasedSurfaceSite, MoleBasedSurfaceSite > > ChemistryLib::PhreeqcIOData::createSurface ( std::optional< BaseLib::ConfigTree > const & config,
MeshLib::Mesh & mesh )
Input File Parameter
prj__chemical_system__surface__site_unit
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
Input File Parameter
prj__chemical_system__surface__site
Input File Parameter
prj__chemical_system__surface__site__name

Definition at line 23 of file CreateSurface.cpp.

25{
26 if (!config)
27 {
28 return {};
29 }
30
31 std::vector<std::variant<DensityBasedSurfaceSite, MoleBasedSurfaceSite>>
32 surface;
33
34 auto const surface_site_unit =
36 config->getConfigAttribute<std::string>("site_unit", "mole");
37
38 if (surface_site_unit == "density")
39 {
40 for (auto const& site_config :
42 config->getConfigSubtreeList("site"))
43 {
45 auto name = site_config.getConfigParameter<std::string>("name");
46
47 auto const site_density =
49 site_config.getConfigParameter<double>("site_density");
50
51 auto const specific_surface_area =
53 site_config.getConfigParameter<double>("specific_surface_area");
54
55 auto const mass =
57 site_config.getConfigParameter<double>("mass");
58
59 surface.push_back(DensityBasedSurfaceSite(
60 std::move(name), site_density, specific_surface_area, mass));
61 }
62
63 return surface;
64 }
65
66 if (surface_site_unit == "mole")
67 {
68 for (auto const& site_config :
70 config->getConfigSubtreeList("site"))
71 {
73 auto name = site_config.getConfigParameter<std::string>("name");
74
75 auto const molality = MeshLib::getOrCreateMeshProperty<double>(
76 mesh, name, MeshLib::MeshItemType::IntegrationPoint, 1);
77
78 surface.push_back(MoleBasedSurfaceSite(std::move(name), molality));
79 }
80
81 return surface;
82 }
83
84 OGS_FATAL("Surface site unit should be either of 'density' or 'mole'.");
85}

References MeshLib::IntegrationPoint, and OGS_FATAL.

Referenced by createChemicalSystem().

◆ 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 22 of file CreateUserPunch.cpp.

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

References MeshLib::IntegrationPoint.

◆ operator<<() [1/6]

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

Definition at line 44 of file Output.cpp.

46{
47 os << "SELECTED_OUTPUT"
48 << "\n";
49 os << "-file " << basic_output_setups.output_file << "\n";
50 os << "-high_precision " << std::boolalpha
51 << basic_output_setups.use_high_precision << "\n";
52 os << "-simulation " << std::boolalpha
53 << BasicOutputSetups::display_simulation_id << "\n";
54 os << "-state " << std::boolalpha << BasicOutputSetups::display_state
55 << "\n";
56 os << "-distance " << std::boolalpha << BasicOutputSetups::display_distance
57 << "\n";
58 os << "-time " << std::boolalpha << BasicOutputSetups::display_current_time
59 << "\n";
60 os << "-step " << std::boolalpha << BasicOutputSetups::display_time_step
61 << "\n";
62
63 return os;
64}

◆ operator<<() [2/6]

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/6]

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

Definition at line 66 of file Output.cpp.

67{
68 os << output.basic_output_setups;
69
70 auto const component_items =
71 output.getOutputItemsByItemType(ItemType::Component);
72 os << "-totals";
73 for (auto const& component_item : component_items)
74 {
75 os << " " << component_item.name;
76 }
77 os << "\n";
78
79 auto const equilibrium_phase_items =
80 output.getOutputItemsByItemType(ItemType::EquilibriumReactant);
81 if (!equilibrium_phase_items.empty())
82 {
83 os << "-equilibrium_phases";
84 for (auto const& equilibrium_phase_item : equilibrium_phase_items)
85 {
86 os << " " << equilibrium_phase_item.name;
87 }
88 os << "\n";
89 }
90
91 auto const kinetic_reactant_items =
92 output.getOutputItemsByItemType(ItemType::KineticReactant);
93 if (!kinetic_reactant_items.empty())
94 {
95 os << "-kinetic_reactants";
96 for (auto const& kinetic_reactant_item : kinetic_reactant_items)
97 {
98 os << " " << kinetic_reactant_item.name;
99 }
100 os << "\n";
101 }
102
103 return os;
104}

◆ operator<<() [4/6]

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

Definition at line 483 of file PhreeqcIO.cpp.

484{
485 bool const fixing_pe =
486 phreeqc_io._chemical_system->aqueous_solution->fixing_pe;
487 if (fixing_pe)
488 {
489 os << "PHASES\n"
490 << "Fix_pe\n"
491 << "e- = e-\n"
492 << "log_k 0.0\n\n";
493 }
494
495 os << phreeqc_io._knobs << "\n";
496
497 os << *phreeqc_io._output << "\n";
498
499 auto const& user_punch = phreeqc_io._user_punch;
500 if (user_punch)
501 {
502 os << *user_punch << "\n";
503 }
504
505 auto const& reaction_rates = phreeqc_io._reaction_rates;
506 if (!reaction_rates.empty())
507 {
508 os << "RATES\n";
509 os << reaction_rates << "\n";
510 }
511
512 for (std::size_t chemical_system_id = 0;
513 chemical_system_id < phreeqc_io._num_chemical_systems;
514 ++chemical_system_id)
515 {
516 os << "SOLUTION " << chemical_system_id + 1 << "\n";
517 phreeqc_io._chemical_system->aqueous_solution->print(
518 os, chemical_system_id);
519
520 auto const& dump = phreeqc_io._dump;
521 if (dump)
522 {
523 auto const& aqueous_solutions_prev = dump->aqueous_solutions_prev;
524 if (!aqueous_solutions_prev.empty())
525 {
526 os << aqueous_solutions_prev[chemical_system_id] << "\n\n";
527 }
528 }
529
530 os << "USE solution none\n";
531 os << "END\n\n";
532
533 os << "USE solution " << chemical_system_id + 1 << "\n\n";
534
535 auto const& equilibrium_reactants =
536 phreeqc_io._chemical_system->equilibrium_reactants;
537 if (!equilibrium_reactants.empty() || fixing_pe)
538 {
539 os << "EQUILIBRIUM_PHASES " << chemical_system_id + 1 << "\n";
540 for (auto const& equilibrium_reactant : equilibrium_reactants)
541 {
542 equilibrium_reactant.print(os, chemical_system_id);
543 }
544 fixing_pe
545 ? os << "Fix_pe "
546 << -phreeqc_io._chemical_system->aqueous_solution->pe0
547 << " O2(g)\n\n"
548 : os << "\n";
549 }
550
551 auto const& kinetic_reactants =
552 phreeqc_io._chemical_system->kinetic_reactants;
553 if (!kinetic_reactants.empty())
554 {
555 os << "KINETICS " << chemical_system_id + 1 << "\n";
556 for (auto const& kinetic_reactant : kinetic_reactants)
557 {
558 kinetic_reactant.print(os, chemical_system_id);
559 }
560 os << "-steps " << phreeqc_io._dt << "\n\n";
561 }
562
563 auto const& surface = phreeqc_io._chemical_system->surface;
564 if (!surface.empty())
565 {
566 os << "SURFACE " << chemical_system_id + 1 << "\n";
567 std::size_t aqueous_solution_id =
568 dump->aqueous_solutions_prev.empty()
569 ? chemical_system_id + 1
570 : phreeqc_io._num_chemical_systems + chemical_system_id + 1;
571 os << "-equilibrate with solution " << aqueous_solution_id << "\n";
572
573 // print unit
574 if (std::holds_alternative<DensityBasedSurfaceSite>(
575 surface.front()))
576 {
577 os << "-sites_units density\n";
578 }
579 else
580 {
581 os << "-sites_units absolute\n";
582 }
583
584 for (auto const& surface_site : surface)
585 {
586 std::visit(
587 overloaded{[&os](DensityBasedSurfaceSite const& s)
588 {
589 os << s.name << " " << s.site_density << " "
590 << s.specific_surface_area << " "
591 << s.mass << "\n";
592 },
593 [&os, chemical_system_id](
594 MoleBasedSurfaceSite const& s) {
595 os << s.name << " "
596 << (*s.molality)[chemical_system_id]
597 << "\n";
598 }},
599 surface_site);
600 }
601
602 // overlook the effect of the buildup of charges onto the surface
603 if (std::holds_alternative<MoleBasedSurfaceSite>(surface.front()))
604 {
605 os << "-no_edl\n";
606 }
607
608 os << "SAVE solution " << chemical_system_id + 1 << "\n";
609 }
610
611 auto const& exchangers = phreeqc_io._chemical_system->exchangers;
612 if (!exchangers.empty())
613 {
614 os << "EXCHANGE " << chemical_system_id + 1 << "\n";
615 std::size_t const aqueous_solution_id =
616 dump->aqueous_solutions_prev.empty()
617 ? chemical_system_id + 1
618 : phreeqc_io._num_chemical_systems + chemical_system_id + 1;
619 os << "-equilibrate with solution " << aqueous_solution_id << "\n";
620 for (auto const& exchanger : exchangers)
621 {
622 exchanger.print(os, chemical_system_id);
623 }
624 os << "SAVE solution " << chemical_system_id + 1 << "\n";
625 }
626
627 os << "END\n\n";
628 }
629
630 auto const& dump = phreeqc_io._dump;
631 if (dump)
632 {
633 dump->print(os, phreeqc_io._num_chemical_systems);
634 }
635
636 return os;
637}

◆ operator<<() [5/6]

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/6]

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 const& 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 676 of file PhreeqcIO.cpp.

677{
678 // Skip the headline
679 in.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
680
681 std::string line;
682 auto const& output = *phreeqc_io._output;
683 auto const& dropped_item_ids = output.dropped_item_ids;
684
685 auto const& surface = phreeqc_io._chemical_system->surface;
686 auto const& exchangers = phreeqc_io._chemical_system->exchangers;
687
688 int const num_skipped_lines =
689 1 + (!surface.empty() ? 1 : 0) + (!exchangers.empty() ? 1 : 0);
690
691 auto& equilibrium_reactants =
692 phreeqc_io._chemical_system->equilibrium_reactants;
693 auto& kinetic_reactants = phreeqc_io._chemical_system->kinetic_reactants;
694
695 for (std::size_t chemical_system_id = 0;
696 chemical_system_id < phreeqc_io._num_chemical_systems;
697 ++chemical_system_id)
698 {
699 // Skip equilibrium calculation result of initial solution
700 for (int i = 0; i < num_skipped_lines; ++i)
701 {
702 in.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
703 }
704
705 // Get calculation result of the solution after the reaction
706 if (!std::getline(in, line))
707 {
708 OGS_FATAL(
709 "Error when reading calculation result of Solution {:d} "
710 "after the reaction.",
711 chemical_system_id);
712 }
713
714 std::vector<double> accepted_items;
715 std::vector<std::string> items;
716 boost::trim_if(line, boost::is_any_of("\t "));
717 boost::algorithm::split(items, line, boost::is_any_of("\t "),
718 boost::token_compress_on);
719 for (int item_id = 0; item_id < static_cast<int>(items.size());
720 ++item_id)
721 {
722 if (std::find(dropped_item_ids.begin(), dropped_item_ids.end(),
723 item_id) == dropped_item_ids.end())
724 {
725 double value;
726 try
727 {
728 value = std::stod(items[item_id]);
729 }
730 catch (const std::invalid_argument& e)
731 {
732 OGS_FATAL(
733 "Invalid argument. Could not convert string '{:s}' to "
734 "double for chemical system {:d}, column {:d}. "
735 "Exception '{:s}' was thrown.",
736 items[item_id], chemical_system_id + 1, item_id,
737 e.what());
738 }
739 catch (const std::out_of_range& e)
740 {
741 OGS_FATAL(
742 "Out of range error. Could not convert string "
743 "'{:s}' to double for chemical system {:d}, column "
744 "{:d}. Exception '{:s}' was thrown.",
745 items[item_id], chemical_system_id + 1, item_id,
746 e.what());
747 }
748 accepted_items.push_back(value);
749 }
750 }
751 assert(accepted_items.size() == output.accepted_items.size());
752
753 auto& aqueous_solution = phreeqc_io._chemical_system->aqueous_solution;
754 auto& components = aqueous_solution->components;
755 auto& user_punch = phreeqc_io._user_punch;
756
757 GlobalIndexType const offset = aqueous_solution->pH->getRangeBegin();
758 GlobalIndexType const global_index = offset + chemical_system_id;
759
760 for (int item_id = 0; item_id < static_cast<int>(accepted_items.size());
761 ++item_id)
762 {
763 auto const& accepted_item = output.accepted_items[item_id];
764 auto const& item_name = accepted_item.name;
765
766 auto compare_by_name = [&item_name](auto const& item)
767 { return item.name == item_name; };
768
769 switch (accepted_item.item_type)
770 {
771 case ItemType::pH:
772 {
773 // Update pH value
774 MathLib::LinAlg::setLocalAccessibleVector(
775 *aqueous_solution->pH);
776 aqueous_solution->pH->set(
777 global_index, std::pow(10, -accepted_items[item_id]));
778 break;
779 }
780 case ItemType::pe:
781 {
782 // Update pe value
783 (*aqueous_solution->pe)[chemical_system_id] =
784 accepted_items[item_id];
785 break;
786 }
787 case ItemType::Component:
788 {
789 // Update component concentrations
790 auto const& component = BaseLib::findElementOrError(
791 components, compare_by_name,
792 [&]() {
793 OGS_FATAL("Could not find component '{:s}'.",
794 item_name);
795 });
796 MathLib::LinAlg::setLocalAccessibleVector(
797 *component.amount);
798 component.amount->set(global_index,
799 accepted_items[item_id]);
800 break;
801 }
802 case ItemType::EquilibriumReactant:
803 {
804 // Update amounts of equilibrium reactant
805 auto const& equilibrium_reactant =
806 BaseLib::findElementOrError(
807 equilibrium_reactants, compare_by_name,
808 [&]()
809 {
810 OGS_FATAL(
811 "Could not find equilibrium reactant "
812 "'{:s}'",
813 item_name);
814 });
815 (*equilibrium_reactant.molality)[chemical_system_id] =
816 accepted_items[item_id];
817 break;
818 }
819 case ItemType::KineticReactant:
820 {
821 // Update amounts of kinetic reactants
822 auto const& kinetic_reactant = BaseLib::findElementOrError(
823 kinetic_reactants, compare_by_name,
824 [&]() {
825 OGS_FATAL("Could not find kinetic reactant '{:s}'.",
826 item_name);
827 });
828 (*kinetic_reactant.molality)[chemical_system_id] =
829 accepted_items[item_id];
830 break;
831 }
832 case ItemType::SecondaryVariable:
833 {
834 assert(user_punch);
835 auto const& secondary_variables =
836 user_punch->secondary_variables;
837 // Update values of secondary variables
838 auto const& secondary_variable =
839 BaseLib::findElementOrError(
840 secondary_variables, compare_by_name,
841 [&]() {
842 OGS_FATAL(
843 "Could not find secondary variable '{:s}'.",
844 item_name);
845 });
846 (*secondary_variable.value)[chemical_system_id] =
847 accepted_items[item_id];
848 break;
849 }
850 }
851 }
852 }
853
854 return in;
855}
GlobalIndexType
GlobalMatrix::IndexType GlobalIndexType
Definition GlobalMatrixVectorTypes.h:51
ChemistryLib::PhreeqcIOData::PhreeqcIO::_chemical_system
std::unique_ptr< ChemicalSystem > _chemical_system
Definition PhreeqcIO.h:111
ChemistryLib::PhreeqcIOData::PhreeqcIO::_user_punch
std::unique_ptr< UserPunch > _user_punch
Definition PhreeqcIO.h:113
ChemistryLib::PhreeqcIOData::PhreeqcIO::_output
std::unique_ptr< Output > const _output
Definition PhreeqcIO.h:114
BaseLib::findElementOrError
ranges::range_reference_t< Range > findElementOrError(Range &range, std::predicate< ranges::range_reference_t< Range > > auto &&predicate, std::invocable auto error_callback)
Definition Algorithm.h:76
MathLib::LinAlg::setLocalAccessibleVector
void setLocalAccessibleVector(PETScVector const &x)
Definition LinAlg.cpp:27

◆ specifyFileName()

std::string ChemistryLib::PhreeqcIOData::specifyFileName ( std::string const & project_file_name,
std::string const & file_extension )
extern

Definition at line 23 of file Output.cpp.

25{
26#ifdef USE_PETSC
27 int mpi_rank;
28
29 MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
30 return project_file_name + "_phreeqc_pid_" + std::to_string(mpi_rank) +
31 file_extension;
32#endif
33
34 return project_file_name + "_phreeqc" + file_extension;
35}