OGS
ChemistryLib Namespace Reference

Namespaces

namespace  PhreeqcIOData
namespace  PhreeqcKernelData
namespace  SelfContainedSolverData

Classes

class  ChemicalSolverInterface

Enumerations

enum class  ChemicalSolver { Phreeqc , PhreeqcKernel , SelfContained }
enum class  ChargeBalance { pH , pe , Unspecified }

Functions

ChargeBalance createChargeBalance (BaseLib::ConfigTree const &config)
template<typename ReactionRate>
std::vector< ReactionRate > createReactionRates (std::optional< BaseLib::ConfigTree > const &config)
template std::vector< PhreeqcIOData::ReactionRatecreateReactionRates< PhreeqcIOData::ReactionRate > (std::optional< BaseLib::ConfigTree > const &config)
template std::vector< PhreeqcKernelData::ReactionRatecreateReactionRates< PhreeqcKernelData::ReactionRate > (std::optional< BaseLib::ConfigTree > const &config)
template<>
std::unique_ptr< ChemicalSolverInterfacecreateChemicalSolverInterface< ChemicalSolver::Phreeqc > (std::vector< std::unique_ptr< MeshLib::Mesh > > const &meshes, std::map< std::string, std::unique_ptr< GlobalLinearSolver > > const &linear_solvers, BaseLib::ConfigTree const &config, std::string const &output_directory)
template<>
std::unique_ptr< ChemicalSolverInterfacecreateChemicalSolverInterface< ChemicalSolver::PhreeqcKernel > (std::vector< std::unique_ptr< MeshLib::Mesh > > const &meshes, std::map< std::string, std::unique_ptr< GlobalLinearSolver > > const &linear_solvers, BaseLib::ConfigTree const &config, std::string const &)
template<>
std::unique_ptr< ChemicalSolverInterfacecreateChemicalSolverInterface< ChemicalSolver::SelfContained > (std::vector< std::unique_ptr< MeshLib::Mesh > > const &meshes, std::map< std::string, std::unique_ptr< GlobalLinearSolver > > const &linear_solvers, BaseLib::ConfigTree const &config, std::string const &)
template<ChemicalSolver chemical_solver>
std::unique_ptr< ChemicalSolverInterfacecreateChemicalSolverInterface (std::vector< std::unique_ptr< MeshLib::Mesh > > const &meshes, std::map< std::string, std::unique_ptr< GlobalLinearSolver > > const &linear_solvers, BaseLib::ConfigTree const &config, std::string const &output_directory)

Enumeration Type Documentation

◆ ChargeBalance

enum class ChemistryLib::ChargeBalance
strong
Enumerator
pH 
pe 
Unspecified 

Definition at line 8 of file ChargeBalance.h.

◆ ChemicalSolver

enum class ChemistryLib::ChemicalSolver
strong
Enumerator
Phreeqc 
PhreeqcKernel 
SelfContained 

Definition at line 8 of file ChemicalSolverType.h.

Function Documentation

◆ createChargeBalance()

ChargeBalance ChemistryLib::createChargeBalance ( BaseLib::ConfigTree const & config)
Input File Parameter
prj__chemical_system__solution__charge_balance

Definition at line 12 of file CreateChargeBalance.cpp.

13{
14 auto const charge_balance =
16 config.getConfigParameterOptional<std::string>("charge_balance");
17
18 if (!charge_balance)
19 {
21 }
22 if (*charge_balance == "pH")
23 {
24 return ChargeBalance::pH;
25 }
26 if (*charge_balance == "pe")
27 {
28 return ChargeBalance::pe;
29 }
30
32 "Error in specifying the way in which charge balance is achieved. "
33 "Adjusting pH value or pe value are supported at this moment.");
34}
#define OGS_FATAL(...)
Definition Error.h:19

References BaseLib::ConfigTree::getConfigParameterOptional(), OGS_FATAL, pe, pH, and Unspecified.

Referenced by ChemistryLib::PhreeqcIOData::createAqueousSolution(), and ChemistryLib::PhreeqcKernelData::createInitialAqueousSolution().

◆ createChemicalSolverInterface()

template<ChemicalSolver chemical_solver>
std::unique_ptr< ChemicalSolverInterface > ChemistryLib::createChemicalSolverInterface ( std::vector< std::unique_ptr< MeshLib::Mesh > > const & meshes,
std::map< std::string, std::unique_ptr< GlobalLinearSolver > > const & linear_solvers,
BaseLib::ConfigTree const & config,
std::string const & output_directory )

◆ createChemicalSolverInterface< ChemicalSolver::Phreeqc >()

template<>
std::unique_ptr< ChemicalSolverInterface > ChemistryLib::createChemicalSolverInterface< ChemicalSolver::Phreeqc > ( std::vector< std::unique_ptr< MeshLib::Mesh > > const & meshes,
std::map< std::string, std::unique_ptr< GlobalLinearSolver > > const & linear_solvers,
BaseLib::ConfigTree const & config,
std::string const & output_directory )
Input File Parameter
prj__chemical_system__mesh
Input File Parameter
prj__chemical_system__linear_solver
Input File Parameter
prj__chemical_system__rates
Input File Parameter
prj__chemical_system__use_stream_for_data_exchange
Input File Parameter
prj__chemical_system__knobs
Input File Parameter
prj__chemical_system__user_punch
Input File Parameter
prj__chemical_system__use_high_precision

Definition at line 37 of file CreateChemicalSolverInterface.cpp.

66{
67 auto mesh_name =
69 config.getConfigParameter<std::string>("mesh");
70
71 // Find and extract mesh from the list of meshes.
72 auto const& mesh = MeshLib::findMeshByName(meshes, mesh_name);
73
74 assert(mesh.getID() != 0);
75 DBUG("Found mesh '{:s}' with id {:d}.", mesh.getName(), mesh.getID());
76
77 auto const ls_name =
79 config.getConfigParameter<std::string>("linear_solver");
80 auto const& linear_solver = BaseLib::getOrError(
81 linear_solvers, ls_name,
82 "A linear solver with the given name does not exist.");
83
84 auto path_to_database = parseDatabasePath(config);
85
86 // chemical system
87 auto chemical_system =
88 PhreeqcIOData::createChemicalSystem(config, *meshes[0]);
89
90 // rates
93 config.getConfigSubtreeOptional("rates"));
94
95 // surface
96 auto const& surface = chemical_system->surface;
97
98 // exchange
99 auto const& exchangers = chemical_system->exchangers;
100
101 // whether to use stream for data exchange
102 auto use_stream_for_data_exchange =
104 config.getConfigParameter<bool>("use_stream_for_data_exchange", true);
105
106 // dump
107 auto const project_file_name =
108 BaseLib::joinPaths(output_directory,
110 config.projectFilePath().string()));
111 // old dump file is deleted if it exists
112 auto dump = surface.empty() && exchangers.empty()
113 ? nullptr
114 : std::make_unique<PhreeqcIOData::Dump>(project_file_name);
115
116 // knobs
117 auto knobs = PhreeqcIOData::createKnobs(
119 config.getConfigSubtree("knobs"));
120
121 // user punch
122 auto user_punch = PhreeqcIOData::createUserPunch(
124 config.getConfigSubtreeOptional("user_punch"), *meshes[0]);
125
126 // output
127 auto const use_high_precision =
129 config.getConfigParameter<bool>("use_high_precision", true);
130 auto output = PhreeqcIOData::createOutput(
131 *chemical_system, user_punch, use_high_precision, project_file_name);
132
133 if (use_stream_for_data_exchange)
134 {
135 INFO("PhreeqcIO will use stringstream for data exchange.");
136 }
137 else
138 {
139 INFO("PhreeqcIO will use file for data exchange.");
140 }
141
142 return std::make_unique<PhreeqcIOData::PhreeqcIO>(
143 mesh, *linear_solver, std::move(project_file_name),
144 std::move(path_to_database), std::move(chemical_system),
145 std::move(reaction_rates), std::move(user_punch), std::move(output),
146 std::move(dump), std::move(knobs), use_stream_for_data_exchange);
147}
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:28
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:22
std::string extractBaseNameWithoutExtension(std::string const &pathname)
std::string joinPaths(std::string const &pathA, std::string const &pathB)
OGS_NO_DANGLING Map::mapped_type & getOrError(Map &map, Key const &key, std::string const &error_message)
Definition Algorithm.h:111
std::unique_ptr< UserPunch > createUserPunch(std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh const &mesh)
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::unique_ptr< ChemicalSystem > createChemicalSystem(BaseLib::ConfigTree const &config, MeshLib::Mesh &mesh)
Knobs createKnobs(BaseLib::ConfigTree const &config)
template std::vector< PhreeqcIOData::ReactionRate > createReactionRates< PhreeqcIOData::ReactionRate >(std::optional< BaseLib::ConfigTree > const &config)
Mesh & findMeshByName(std::vector< std::unique_ptr< Mesh > > const &meshes, std::string_view const name)
Definition Mesh.cpp:354

References BaseLib::ConfigTree::getConfigParameter(), INFO(), BaseLib::IsFileExisting(), BaseLib::joinPaths(), OGS_FATAL, and BaseLib::ConfigTree::projectDirectory().

◆ createChemicalSolverInterface< ChemicalSolver::PhreeqcKernel >()

template<>
std::unique_ptr< ChemicalSolverInterface > ChemistryLib::createChemicalSolverInterface< ChemicalSolver::PhreeqcKernel > ( std::vector< std::unique_ptr< MeshLib::Mesh > > const & meshes,
std::map< std::string, std::unique_ptr< GlobalLinearSolver > > const & linear_solvers,
BaseLib::ConfigTree const & config,
std::string const &  )
Input File Parameter
prj__chemical_system__linear_solver
Input File Parameter
prj__chemical_system__solution
Input File Parameter
prj__chemical_system__kinetic_reactants
Input File Parameter
prj__chemical_system__rates
Input File Parameter
prj__chemical_system__equilibrium_reactants

Definition at line 37 of file CreateChemicalSolverInterface.cpp.

156{
157 auto mesh = *meshes[0];
158
159 auto const ls_name =
161 config.getConfigParameter<std::string>("linear_solver");
162 auto const& linear_solver = BaseLib::getOrError(
163 linear_solvers, ls_name,
164 "A linear solver with the given name does not exist.");
165
166 auto path_to_database = parseDatabasePath(config);
167
168 // TODO (renchao): remove mapping process id to component name.
169 std::vector<std::pair<int, std::string>> process_id_to_component_name_map;
170 // solution
171 auto aqueous_solution = PhreeqcKernelData::createAqueousSolution(
173 config.getConfigSubtree("solution"),
174 process_id_to_component_name_map);
175
176 // kinetic reactants
177 auto kinetic_reactants = PhreeqcKernelData::createKineticReactants(
179 config.getConfigSubtreeOptional("kinetic_reactants"), mesh);
180
181 // rates
184 config.getConfigSubtreeOptional("rates"));
185
186 // equilibrium reactants
187 auto equilibrium_reactants = PhreeqcKernelData::createEquilibriumReactants(
189 config.getConfigSubtreeOptional("equilibrium_reactants"), mesh);
190
191 return std::make_unique<PhreeqcKernelData::PhreeqcKernel>(
192 mesh, *linear_solver, mesh.computeNumberOfBaseNodes(),
193 process_id_to_component_name_map, std::move(path_to_database),
194 std::move(aqueous_solution), std::move(equilibrium_reactants),
195 std::move(kinetic_reactants), std::move(reaction_rates));
196}
std::unique_ptr< EquilibriumReactants > createEquilibriumReactants(std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh const &mesh)
std::unique_ptr< Kinetics > createKineticReactants(std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh const &mesh)
AqueousSolution createAqueousSolution(BaseLib::ConfigTree const &config, std::vector< std::pair< int, std::string > > const &process_id_to_component_name_map)
template std::vector< PhreeqcKernelData::ReactionRate > createReactionRates< PhreeqcKernelData::ReactionRate >(std::optional< BaseLib::ConfigTree > const &config)

◆ createChemicalSolverInterface< ChemicalSolver::SelfContained >()

template<>
std::unique_ptr< ChemicalSolverInterface > ChemistryLib::createChemicalSolverInterface< ChemicalSolver::SelfContained > ( std::vector< std::unique_ptr< MeshLib::Mesh > > const & meshes,
std::map< std::string, std::unique_ptr< GlobalLinearSolver > > const & linear_solvers,
BaseLib::ConfigTree const & config,
std::string const &  )
Input File Parameter
prj__chemical_system__mesh
Input File Parameter
prj__chemical_system__linear_solver
Input File Parameter
prj__chemical_system__chemical_reactions
Input File Parameter
prj__chemical_system__number_of_components

Definition at line 37 of file CreateChemicalSolverInterface.cpp.

205{
206 auto mesh_name =
208 config.getConfigParameter<std::string>("mesh");
209
210 // Find and extract mesh from the list of meshes.
211 auto const& mesh = MeshLib::findMeshByName(meshes, mesh_name);
212
213 assert(mesh.getID() != 0);
214 DBUG("Found mesh '{:s}' with id {:d}.", mesh.getName(), mesh.getID());
215
216 auto const ls_name =
218 config.getConfigParameter<std::string>("linear_solver");
219 auto const& linear_solver = BaseLib::getOrError(
220 linear_solvers, ls_name,
221 "A linear solver with the given name does not exist.");
222
223 auto chemical_reaction_data =
226 config.getConfigSubtree("chemical_reactions"));
227
228 // create sparse stoichiometric matrix
229 std::vector<double> stoichiometric_matrix_vec;
230 for (auto const& per_chemical_reaction_data : chemical_reaction_data)
231 {
232 stoichiometric_matrix_vec.insert(
233 stoichiometric_matrix_vec.end(),
234 per_chemical_reaction_data->stoichiometric_vector.begin(),
235 per_chemical_reaction_data->stoichiometric_vector.end());
236 }
237
238 auto const num_components =
240 config.getConfigParameter<int>("number_of_components");
241
242 Eigen::MatrixXd const stoichiometric_matrix = MathLib::toMatrix(
243 stoichiometric_matrix_vec, num_components, num_components);
244
245 Eigen::SparseMatrix<double> sparse_stoichiometric_matrix;
246 sparse_stoichiometric_matrix = stoichiometric_matrix.sparseView();
247
248 return std::make_unique<SelfContainedSolverData::SelfContainedSolver>(
249 mesh, *linear_solver, sparse_stoichiometric_matrix,
250 std::move(chemical_reaction_data));
251}
std::vector< std::unique_ptr< ChemicalReaction > > createChemicalReactionData(BaseLib::ConfigTree const &config)
Eigen::Map< const Matrix > toMatrix(std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)

◆ createReactionRates()

template<typename ReactionRate>
std::vector< ReactionRate > ChemistryLib::createReactionRates ( std::optional< BaseLib::ConfigTree > const & config)
Input File Parameter
prj__chemical_system__rates__rate
Input File Parameter
prj__chemical_system__rates__rate__kinetic_reactant
Input File Parameter
prj__chemical_system__rates__rate__expression
Input File Parameter
prj__chemical_system__rates__rate__expression__statement

Definition at line 13 of file CreateReactionRate.cpp.

15{
16 if (!config)
17 {
18 return {};
19 }
20
21 std::vector<ReactionRate> reaction_rates;
22 for (auto const& rate_config :
24 config->getConfigSubtreeList("rate"))
25 {
26 auto kinetic_reactant =
28 rate_config.getConfigParameter<std::string>("kinetic_reactant");
29
30 auto const expression_config =
32 rate_config.getConfigSubtree("expression");
33 auto const& statements =
35 expression_config.getConfigParameterList<std::string>("statement");
36
37 std::vector<std::string> expression_statements;
38 expression_statements.reserve(statements.size());
39 std::copy(begin(statements),
40 end(statements),
41 back_inserter(expression_statements));
42
43 reaction_rates.emplace_back(std::move(kinetic_reactant),
44 std::move(expression_statements));
45 }
46
47 return reaction_rates;
48}

◆ createReactionRates< PhreeqcIOData::ReactionRate >()

◆ createReactionRates< PhreeqcKernelData::ReactionRate >()