OGS
ChemistryLib::PhreeqcIOData Namespace Reference

Namespaces

namespace  anonymous_namespace{PhreeqcIO.cpp}

Classes

struct  AqueousSolution
class  BasicOutputSetups
 Controls which built-in PHREEQC columns appear in the output file. More...
struct  ChemicalSystem
 Complete description of one local reactive system passed to PHREEQC. More...
struct  Component
struct  DensityBasedSurfaceSite
struct  Dump
struct  EquilibriumReactant
struct  ExchangeSite
struct  KineticReactant
struct  Knobs
struct  MoleBasedSurfaceSite
struct  Output
 Specification of which PHREEQC output columns are imported into OpenGeoSys. More...
struct  OutputItem
 One PHREEQC output column that OpenGeoSys will keep. More...
class  PhreeqcIO
 Drives the chemistry step in the operator-split reactive transport loop. More...
struct  ReactionRate
struct  SecondaryVariable
struct  UserPunch

Enumerations

enum class  ItemType {
  pH , pe , Component , EquilibriumReactant ,
  KineticReactant , SecondaryVariable
}
 Category used to interpret a PHREEQC output column. More...

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

Category used to interpret a PHREEQC output column.

Distinguishes pH / pe, total component amounts, equilibrium-phase amounts, kinetic-phase amounts, and user-defined secondary variables.

Enumerator
pH 
pe 
Component 
EquilibriumReactant 
KineticReactant 
SecondaryVariable 

Definition at line 71 of file ChemistryLib/PhreeqcIOData/Output.h.

72{
73 pH,
74 pe,
75 Component,
76 EquilibriumReactant,
77 KineticReactant,
78 SecondaryVariable
79};

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 17 of file PhreeqcIOData/CreateAqueousSolution.cpp.

19{
21 auto const fixing_pe = config.getConfigAttribute<bool>("fixing_pe", false);
22
24 auto const temperature = config.getConfigParameter<double>("temperature");
25
27 auto const pressure = config.getConfigParameter<double>("pressure");
28
30 auto const pe0 = config.getConfigParameter<double>("pe");
31
32 auto pe = MeshLib::getOrCreateMeshProperty<double>(
33 mesh, "pe", MeshLib::MeshItemType::IntegrationPoint, 1);
34
35 auto components = createSolutionComponents(config);
36
37 auto charge_balance = createChargeBalance(config);
38
39 return std::make_unique<AqueousSolution>(fixing_pe, temperature, pressure,
40 pe, pe0, std::move(components),
41 charge_balance);
42}
ChemistryLib::PhreeqcIOData::createSolutionComponents
std::vector< Component > createSolutionComponents(BaseLib::ConfigTree const &config)
Definition CreateSolutionComponent.cpp:13
ChemistryLib::createChargeBalance
ChargeBalance createChargeBalance(BaseLib::ConfigTree const &config)
Definition CreateChargeBalance.cpp:12
MeshLib::getOrCreateMeshProperty
PropertyVector< T > * getOrCreateMeshProperty(Mesh &mesh, std::string const &property_name, MeshItemType const item_type, int const number_of_components)
Definition getOrCreateMeshProperty.h:17

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

24{
25 // solution
26 auto aqueous_solution = createAqueousSolution(
28 config.getConfigSubtree("solution"), mesh);
29
30 // kinetic reactants
31 auto kinetic_reactants = createKineticReactants(
33 config.getConfigSubtreeOptional("kinetic_reactants"), mesh);
34
35 // equilibrium reactants
36 auto equilibrium_reactants = createEquilibriumReactants(
38 config.getConfigSubtreeOptional("equilibrium_reactants"), mesh);
39
40 // exchangers
41 auto exchangers = createExchange(
43 config.getConfigSubtreeOptional("exchangers"), mesh);
44
45 // surface
46 auto surface = createSurface(
48 config.getConfigSubtreeOptional("surface"), 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::move(surface));
55}
ChemistryLib::PhreeqcIOData::createSurface
std::vector< std::variant< DensityBasedSurfaceSite, MoleBasedSurfaceSite > > createSurface(std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
Definition CreateSurface.cpp:16
ChemistryLib::PhreeqcIOData::createAqueousSolution
std::unique_ptr< AqueousSolution > createAqueousSolution(BaseLib::ConfigTree const &config, MeshLib::Mesh &mesh)
Definition PhreeqcIOData/CreateAqueousSolution.cpp:17
ChemistryLib::PhreeqcIOData::createExchange
std::vector< ExchangeSite > createExchange(std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
Definition CreateExchange.cpp:15
ChemistryLib::PhreeqcIOData::createEquilibriumReactants
std::vector< EquilibriumReactant > createEquilibriumReactants(std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
Definition PhreeqcIOData/CreateEquilibriumReactants.cpp:17
ChemistryLib::PhreeqcIOData::createKineticReactants
std::vector< KineticReactant > createKineticReactants(std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
Definition PhreeqcIOData/CreateKineticReactant.cpp:17

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 17 of file PhreeqcIOData/CreateEquilibriumReactants.cpp.

19{
20 if (!config)
21 {
22 return {};
23 }
24
25 std::vector<EquilibriumReactant> equilibrium_reactants;
26 for (
27 auto const& equilibrium_reactant_config :
29 config->getConfigSubtreeList("phase_component"))
30 {
31 auto name =
33 equilibrium_reactant_config.getConfigParameter<std::string>("name");
34
35 double const saturation_index =
37 equilibrium_reactant_config.getConfigParameter<double>(
38 "saturation_index");
39
40 auto reaction_irreversibility =
42 equilibrium_reactant_config.getConfigParameter<std::string>(
43 "reaction_irreversibility", "");
44
45 if (!reaction_irreversibility.empty() &&
46 (reaction_irreversibility != "dissolve_only" &&
47 reaction_irreversibility != "precipitate_only"))
48 {
49 OGS_FATAL(
50 "{:s}: reaction direction only allows `dissolve_only` or "
51 "`precipitate_only`",
52 name);
53 }
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 equilibrium_reactants.emplace_back(std::move(name),
75 molality,
76 molality_prev,
77 volume_fraction,
78 volume_fraction_prev,
79 mesh_prop_molality,
80 saturation_index,
81 std::move(reaction_irreversibility));
82 }
83
84 return equilibrium_reactants;
85}
OGS_FATAL
#define OGS_FATAL(...)
Definition Error.h:19
MeshLib::Mesh::getNumberOfElements
std::size_t getNumberOfElements() const
Get the number of elements.
Definition Mesh.h:88

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

17{
18 if (!config)
19 {
20 return {};
21 }
22
23 std::vector<ExchangeSite> exchangers;
24 for (auto const& site_config :
26 config->getConfigSubtreeList("exchange_site"))
27 {
29 auto name = site_config.getConfigParameter<std::string>(
30 "ion_exchanging_species");
31
32 auto const molality = MeshLib::getOrCreateMeshProperty<double>(
33 mesh, name, MeshLib::MeshItemType::IntegrationPoint, 1);
34
35 exchangers.emplace_back(std::move(name), molality);
36 }
37
38 return exchangers;
39}

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 17 of file PhreeqcIOData/CreateKineticReactant.cpp.

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

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 13 of file CreateKnobs.cpp.

14{
15 auto const max_iterations =
17 config.getConfigParameter<int>("max_iter");
18
19 auto const relative_convergence_tolerance =
21 config.getConfigParameter<double>("relative_convergence_tolerance");
22
23 auto const tolerance =
25 config.getConfigParameter<double>("tolerance");
26
27 auto const step_size =
29 config.getConfigParameter<int>("step_size");
30
31 auto const scaling =
33 config.getConfigParameter<bool>("scaling");
34
35 return {max_iterations, relative_convergence_tolerance, tolerance,
36 step_size, scaling};
37}

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 18 of file ChemistryLib/PhreeqcIOData/CreateOutput.cpp.

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

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 13 of file CreateSolutionComponent.cpp.

15{
16 std::vector<Component> components;
18 auto components_config = config.getConfigSubtree("components");
19
20 for (
21 auto const& component_config :
23 components_config.getConfigSubtreeList("component"))
24 {
25 auto const component_name = component_config.getValue<std::string>();
26 auto const chemical_formula =
28 component_config.getConfigAttribute<std::string>("chemical_formula",
29 "");
30 components.emplace_back(component_name, chemical_formula);
31 }
32
33 return components;
34}

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 16 of file CreateSurface.cpp.

18{
19 if (!config)
20 {
21 return {};
22 }
23
24 std::vector<std::variant<DensityBasedSurfaceSite, MoleBasedSurfaceSite>>
25 surface;
26
27 auto const surface_site_unit =
29 config->getConfigAttribute<std::string>("site_unit", "mole");
30
31 if (surface_site_unit == "density")
32 {
33 for (auto const& site_config :
35 config->getConfigSubtreeList("site"))
36 {
38 auto name = site_config.getConfigParameter<std::string>("name");
39
40 auto const site_density =
42 site_config.getConfigParameter<double>("site_density");
43
44 auto const specific_surface_area =
46 site_config.getConfigParameter<double>("specific_surface_area");
47
48 auto const mass =
50 site_config.getConfigParameter<double>("mass");
51
52 surface.push_back(DensityBasedSurfaceSite(
53 std::move(name), site_density, specific_surface_area, mass));
54 }
55
56 return surface;
57 }
58
59 if (surface_site_unit == "mole")
60 {
61 for (auto const& site_config :
63 config->getConfigSubtreeList("site"))
64 {
66 auto name = site_config.getConfigParameter<std::string>("name");
67
68 auto const molality = MeshLib::getOrCreateMeshProperty<double>(
69 mesh, name, MeshLib::MeshItemType::IntegrationPoint, 1);
70
71 surface.push_back(MoleBasedSurfaceSite(std::move(name), molality));
72 }
73
74 return surface;
75 }
76
77 OGS_FATAL("Surface site unit should be either of 'density' or 'mole'.");
78}

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

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

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 39 of file ChemistryLib/PhreeqcIOData/Output.cpp.

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

◆ operator<<() [2/6]

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

Definition at line 12 of file Knobs.cpp.

13{
14 os << "KNOBS"
15 << "\n";
16 os << "-iterations " << knobs.max_iterations << "\n";
17 os << "-convergence_tolerance " << knobs.relative_convergence_tolerance
18 << "\n";
19 os << "-tolerance " << knobs.tolerance << "\n";
20 os << "-step_size " << knobs.step_size << "\n";
21 os << "-diagonal_scale " << knobs.scaling << "\n";
22
23 return os;
24}

◆ operator<<() [3/6]

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

Definition at line 61 of file ChemistryLib/PhreeqcIOData/Output.cpp.

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

◆ operator<<() [4/6]

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

Definition at line 476 of file PhreeqcIO.cpp.

477{
478 bool const fixing_pe =
479 phreeqc_io._chemical_system->aqueous_solution->fixing_pe;
480 if (fixing_pe)
481 {
482 os << "PHASES\n"
483 << "Fix_pe\n"
484 << "e- = e-\n"
485 << "log_k 0.0\n\n";
486 }
487
488 os << phreeqc_io._knobs << "\n";
489
490 os << *phreeqc_io._output << "\n";
491
492 auto const& user_punch = phreeqc_io._user_punch;
493 if (user_punch)
494 {
495 os << *user_punch << "\n";
496 }
497
498 auto const& reaction_rates = phreeqc_io._reaction_rates;
499 if (!reaction_rates.empty())
500 {
501 os << "RATES\n";
502 os << reaction_rates << "\n";
503 }
504
505 for (std::size_t chemical_system_id = 0;
506 chemical_system_id < phreeqc_io._num_chemical_systems;
507 ++chemical_system_id)
508 {
509 os << "SOLUTION " << chemical_system_id + 1 << "\n";
510 phreeqc_io._chemical_system->aqueous_solution->print(
511 os, chemical_system_id);
512
513 auto const& dump = phreeqc_io._dump;
514 if (dump)
515 {
516 auto const& aqueous_solutions_prev = dump->aqueous_solutions_prev;
517 // Print previous aqueous solution if it exists.
518 if (!aqueous_solutions_prev.empty())
519 {
520 os << aqueous_solutions_prev[chemical_system_id] << "\n\n";
521 }
522 }
523
524 os << "USE solution none\n";
525 os << "END\n\n";
526
527 os << "USE solution " << chemical_system_id + 1 << "\n\n";
528
529 auto const& equilibrium_reactants =
530 phreeqc_io._chemical_system->equilibrium_reactants;
531 if (!equilibrium_reactants.empty() || fixing_pe)
532 {
533 os << "EQUILIBRIUM_PHASES " << chemical_system_id + 1 << "\n";
534 for (auto const& equilibrium_reactant : equilibrium_reactants)
535 {
536 equilibrium_reactant.print(os, chemical_system_id);
537 }
538 fixing_pe
539 ? os << "Fix_pe "
540 << -phreeqc_io._chemical_system->aqueous_solution->pe0
541 << " O2(g)\n\n"
542 : os << "\n";
543 }
544
545 auto const& kinetic_reactants =
546 phreeqc_io._chemical_system->kinetic_reactants;
547 if (!kinetic_reactants.empty())
548 {
549 os << "KINETICS " << chemical_system_id + 1 << "\n";
550 for (auto const& kinetic_reactant : kinetic_reactants)
551 {
552 kinetic_reactant.print(os, chemical_system_id);
553 }
554 os << "-steps " << phreeqc_io._dt << "\n\n";
555 }
556
557 auto const& surface = phreeqc_io._chemical_system->surface;
558 if (!surface.empty())
559 {
560 // To get the amount of surface species from the previous time step,
561 // an equilibration calculation with the previous aqueous solution
562 // is needed. The previous aqueous solution is saved using the
563 // PHREEQC keyword "DUMP" in the dump file and then stored as
564 // SOLUTION_RAW within aqueous_solutions_prev. Along with the
565 // PHREEQC keyword 'SURFACE', distinguish between the current and
566 // previous solutions, the previous solution number is added to the
567 // total number of chemical systems (_num_chemical_systems).
568 os << "SURFACE " << chemical_system_id + 1 << "\n";
569 std::size_t aqueous_solution_id =
570 dump->aqueous_solutions_prev.empty()
571 ? chemical_system_id + 1
572 : phreeqc_io._num_chemical_systems + chemical_system_id + 1;
573 os << "-equilibrate with solution " << aqueous_solution_id << "\n";
574
575 // print unit
576 if (std::holds_alternative<DensityBasedSurfaceSite>(
577 surface.front()))
578 {
579 os << "-sites_units density\n";
580 }
581 else
582 {
583 os << "-sites_units absolute\n";
584 }
585
586 for (auto const& surface_site : surface)
587 {
588 std::visit(
589 overloaded{[&os](DensityBasedSurfaceSite const& s)
590 {
591 os << s.name << " " << s.site_density << " "
592 << s.specific_surface_area << " "
593 << s.mass << "\n";
594 },
595 [&os, chemical_system_id](
596 MoleBasedSurfaceSite const& s) {
597 os << s.name << " "
598 << (*s.molality)[chemical_system_id]
599 << "\n";
600 }},
601 surface_site);
602 }
603
604 // overlook the effect of the buildup of charges onto the surface
605 if (std::holds_alternative<MoleBasedSurfaceSite>(surface.front()))
606 {
607 os << "-no_edl\n";
608 }
609 // Here the PHREEQC "save" command is printed to save the solution
610 // internally in PHREEQC. The solution composition is saved after
611 // equilibration with the surface has been completed. The solution
612 // will not be saved for use in the next PHREEQC run.
613 os << "SAVE solution " << chemical_system_id + 1 << "\n";
614 }
615
616 auto const& exchangers = phreeqc_io._chemical_system->exchangers;
617 if (!exchangers.empty())
618 {
619 os << "EXCHANGE " << chemical_system_id + 1 << "\n";
620 std::size_t const aqueous_solution_id =
621 dump->aqueous_solutions_prev.empty()
622 ? chemical_system_id + 1
623 : phreeqc_io._num_chemical_systems + chemical_system_id + 1;
624 os << "-equilibrate with solution " << aqueous_solution_id << "\n";
625 for (auto const& exchanger : exchangers)
626 {
627 exchanger.print(os, chemical_system_id);
628 }
629 os << "SAVE solution " << chemical_system_id + 1 << "\n";
630 }
631
632 os << "END\n\n";
633 }
634
635 auto const& dump = phreeqc_io._dump;
636 if (dump)
637 {
638 dump->print(os, phreeqc_io._num_chemical_systems);
639 }
640
641 return os;
642}

◆ operator<<() [5/6]

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

Definition at line 12 of file PhreeqcIOData/ReactionRate.cpp.

13{
14 os << reaction_rate.kinetic_reactant << "\n";
15 os << "-start"
16 << "\n";
17 int line_number = 1;
18 for (auto const& expression_statement : reaction_rate.expression_statements)
19 {
20 os << line_number << " " << expression_statement << "\n";
21 ++line_number;
22 }
23 os << "-end"
24 << "\n";
25
26 return os;
27}

◆ operator<<() [6/6]

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

Definition at line 20 of file UserPunch.cpp.

21{
22 os << "USER_PUNCH"
23 << "\n";
24 os << "-headings ";
25 auto const& secondary_variables = user_punch.secondary_variables;
26 for (auto const& secondary_variable : secondary_variables)
27 {
28 os << secondary_variable.name << " ";
29 }
30 os << "\n";
31
32 os << "-start"
33 << "\n";
34 int line_number = 1;
35 for (auto const& statement : user_punch.statements)
36 {
37 os << line_number << " " << statement << "\n";
38 ++line_number;
39 }
40 os << "-end"
41 << "\n";
42
43 return os;
44}

◆ operator>>()

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

Definition at line 681 of file PhreeqcIO.cpp.

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

◆ specifyFileName()

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

Definition at line 18 of file ChemistryLib/PhreeqcIOData/Output.cpp.

20{
21#ifdef USE_PETSC
22 int mpi_rank;
23
24 MPI_Comm_rank(BaseLib::MPI::OGS_COMM_WORLD, &mpi_rank);
25 return project_file_name + "_phreeqc_pid_" + std::to_string(mpi_rank) +
26 file_extension;
27#endif
28
29 return project_file_name + "_phreeqc" + file_extension;
30}
BaseLib::MPI::OGS_COMM_WORLD
MPI_Comm OGS_COMM_WORLD
Definition MPI.cpp:9

References BaseLib::MPI::OGS_COMM_WORLD.

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