OGS
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
MeshLib::getOrCreateMeshProperty
PropertyVector< T > * getOrCreateMeshProperty(Mesh &mesh, std::string const &property_name, MeshItemType const item_type, int const number_of_components)
Definition getOrCreateMeshProperty.h:23

References ChemistryLib::createChargeBalance(), createSolutionComponents(), BaseLib::ConfigTree::getConfigAttribute(), BaseLib::ConfigTree::getConfigParameter(), MeshLib::getOrCreateMeshProperty(), 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::createSurface
std::vector< std::variant< DensityBasedSurfaceSite, MoleBasedSurfaceSite > > createSurface(std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
Definition CreateSurface.cpp:23
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:99

References MeshLib::Cell, MeshLib::Mesh::getNumberOfElements(), MeshLib::getOrCreateMeshProperty(), 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::getOrCreateMeshProperty(), and 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 assert(mesh_prop_molality->size() == 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(), MeshLib::getOrCreateMeshProperty(), 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().

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::getOrCreateMeshProperty(), 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::getOrCreateMeshProperty(), and MeshLib::IntegrationPoint.

◆ operator<<() [1/6]

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

Definition at line 46 of file Output.cpp.

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

◆ 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 68 of file Output.cpp.

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

◆ 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 // Print previous aqueous solution if it exists.
525 if (!aqueous_solutions_prev.empty())
526 {
527 os << aqueous_solutions_prev[chemical_system_id] << "\n\n";
528 }
529 }
530
531 os << "USE solution none\n";
532 os << "END\n\n";
533
534 os << "USE solution " << chemical_system_id + 1 << "\n\n";
535
536 auto const& equilibrium_reactants =
537 phreeqc_io._chemical_system->equilibrium_reactants;
538 if (!equilibrium_reactants.empty() || fixing_pe)
539 {
540 os << "EQUILIBRIUM_PHASES " << chemical_system_id + 1 << "\n";
541 for (auto const& equilibrium_reactant : equilibrium_reactants)
542 {
543 equilibrium_reactant.print(os, chemical_system_id);
544 }
545 fixing_pe
546 ? os << "Fix_pe "
547 << -phreeqc_io._chemical_system->aqueous_solution->pe0
548 << " O2(g)\n\n"
549 : os << "\n";
550 }
551
552 auto const& kinetic_reactants =
553 phreeqc_io._chemical_system->kinetic_reactants;
554 if (!kinetic_reactants.empty())
555 {
556 os << "KINETICS " << chemical_system_id + 1 << "\n";
557 for (auto const& kinetic_reactant : kinetic_reactants)
558 {
559 kinetic_reactant.print(os, chemical_system_id);
560 }
561 os << "-steps " << phreeqc_io._dt << "\n\n";
562 }
563
564 auto const& surface = phreeqc_io._chemical_system->surface;
565 if (!surface.empty())
566 {
567 // To get the amount of surface species from the previous time step,
568 // an equilibration calculation with the previous aqueous solution
569 // is needed. The previous aqueous solution is saved using the
570 // PHREEQC keyword "DUMP" in the dump file and then stored as
571 // SOLUTION_RAW within aqueous_solutions_prev. Along with the
572 // PHREEQC keyword 'SURFACE', distinguish between the current and
573 // previous solutions, the previous solution number is added to the
574 // total number of chemical systems (_num_chemical_systems).
575 os << "SURFACE " << chemical_system_id + 1 << "\n";
576 std::size_t aqueous_solution_id =
577 dump->aqueous_solutions_prev.empty()
578 ? chemical_system_id + 1
579 : phreeqc_io._num_chemical_systems + chemical_system_id + 1;
580 os << "-equilibrate with solution " << aqueous_solution_id << "\n";
581
582 // print unit
583 if (std::holds_alternative<DensityBasedSurfaceSite>(
584 surface.front()))
585 {
586 os << "-sites_units density\n";
587 }
588 else
589 {
590 os << "-sites_units absolute\n";
591 }
592
593 for (auto const& surface_site : surface)
594 {
595 std::visit(
596 overloaded{[&os](DensityBasedSurfaceSite const& s)
597 {
598 os << s.name << " " << s.site_density << " "
599 << s.specific_surface_area << " "
600 << s.mass << "\n";
601 },
602 [&os, chemical_system_id](
603 MoleBasedSurfaceSite const& s) {
604 os << s.name << " "
605 << (*s.molality)[chemical_system_id]
606 << "\n";
607 }},
608 surface_site);
609 }
610
611 // overlook the effect of the buildup of charges onto the surface
612 if (std::holds_alternative<MoleBasedSurfaceSite>(surface.front()))
613 {
614 os << "-no_edl\n";
615 }
616 // Here the PHREEQC "save" command is printed to save the solution
617 // internally in PHREEQC. The solution composition is saved after
618 // equilibration with the surface has been completed. The solution
619 // will not be saved for use in the next PHREEQC run.
620 os << "SAVE solution " << chemical_system_id + 1 << "\n";
621 }
622
623 auto const& exchangers = phreeqc_io._chemical_system->exchangers;
624 if (!exchangers.empty())
625 {
626 os << "EXCHANGE " << chemical_system_id + 1 << "\n";
627 std::size_t const aqueous_solution_id =
628 dump->aqueous_solutions_prev.empty()
629 ? chemical_system_id + 1
630 : phreeqc_io._num_chemical_systems + chemical_system_id + 1;
631 os << "-equilibrate with solution " << aqueous_solution_id << "\n";
632 for (auto const& exchanger : exchangers)
633 {
634 exchanger.print(os, chemical_system_id);
635 }
636 os << "SAVE solution " << chemical_system_id + 1 << "\n";
637 }
638
639 os << "END\n\n";
640 }
641
642 auto const& dump = phreeqc_io._dump;
643 if (dump)
644 {
645 dump->print(os, phreeqc_io._num_chemical_systems);
646 }
647
648 return os;
649}

◆ 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 688 of file PhreeqcIO.cpp.

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

27{
28#ifdef USE_PETSC
29 int mpi_rank;
30
31 MPI_Comm_rank(BaseLib::MPI::OGS_COMM_WORLD, &mpi_rank);
32 return project_file_name + "_phreeqc_pid_" + std::to_string(mpi_rank) +
33 file_extension;
34#endif
35
36 return project_file_name + "_phreeqc" + file_extension;
37}
BaseLib::MPI::OGS_COMM_WORLD
MPI_Comm OGS_COMM_WORLD
Definition MPI.cpp:15

References BaseLib::MPI::OGS_COMM_WORLD.

Referenced by ChemistryLib::PhreeqcIOData::BasicOutputSetups::BasicOutputSetups(), ChemistryLib::PhreeqcIOData::Dump::Dump(), and ChemistryLib::PhreeqcIOData::PhreeqcIO::PhreeqcIO().