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 15 of file ChargeBalance.h.

◆ ChemicalSolver

enum class ChemistryLib::ChemicalSolver
strong
Enumerator
Phreeqc 
PhreeqcKernel 
SelfContained 

Definition at line 15 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 19 of file CreateChargeBalance.cpp.

20{
21 auto const charge_balance =
23 config.getConfigParameterOptional<std::string>("charge_balance");
24
25 if (!charge_balance)
26 {
27 return ChargeBalance::Unspecified;
28 }
29 if (*charge_balance == "pH")
30 {
31 return ChargeBalance::pH;
32 }
33 if (*charge_balance == "pe")
34 {
35 return ChargeBalance::pe;
36 }
37
39 "Error in specifying the way in which charge balance is achieved. "
40 "Adjusting pH value or pe value are supported at this moment.");
41}
#define OGS_FATAL(...)
Definition Error.h:26

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__knobs
Input File Parameter
prj__chemical_system__user_punch
Input File Parameter
prj__chemical_system__use_high_precision

Definition at line 44 of file CreateChemicalSolverInterface.cpp.

73{
74 auto mesh_name =
76 config.getConfigParameter<std::string>("mesh");
77
78 // Find and extract mesh from the list of meshes.
79 auto const& mesh = MeshLib::findMeshByName(meshes, mesh_name);
80
81 assert(mesh.getID() != 0);
82 DBUG("Found mesh '{:s}' with id {:d}.", mesh.getName(), mesh.getID());
83
84 auto const ls_name =
86 config.getConfigParameter<std::string>("linear_solver");
87 auto const& linear_solver = BaseLib::getOrError(
88 linear_solvers, ls_name,
89 "A linear solver with the given name does not exist.");
90
91 auto path_to_database = parseDatabasePath(config);
92
93 // chemical system
94 auto chemical_system =
95 PhreeqcIOData::createChemicalSystem(config, *meshes[0]);
96
97 // rates
98 auto reaction_rates = createReactionRates<PhreeqcIOData::ReactionRate>(
100 config.getConfigSubtreeOptional("rates"));
101
102 // surface
103 auto const& surface = chemical_system->surface;
104
105 // exchange
106 auto const& exchangers = chemical_system->exchangers;
107
108 // dump
109 auto const project_file_name = BaseLib::joinPaths(
110 output_directory,
111 BaseLib::extractBaseNameWithoutExtension(config.getProjectFileName()));
112
113 auto dump = surface.empty() && exchangers.empty()
114 ? nullptr
115 : std::make_unique<PhreeqcIOData::Dump>(project_file_name);
116
117 // knobs
118 auto knobs = PhreeqcIOData::createKnobs(
120 config.getConfigSubtree("knobs"));
121
122 // user punch
123 auto user_punch = PhreeqcIOData::createUserPunch(
125 config.getConfigSubtreeOptional("user_punch"), *meshes[0]);
126
127 // output
128 auto const use_high_precision =
130 config.getConfigParameter<bool>("use_high_precision", true);
131 auto output = PhreeqcIOData::createOutput(
132 *chemical_system, user_punch, use_high_precision, project_file_name);
133
134 return std::make_unique<PhreeqcIOData::PhreeqcIO>(
135 mesh, *linear_solver, std::move(project_file_name),
136 std::move(path_to_database), std::move(chemical_system),
137 std::move(reaction_rates), std::move(user_punch), std::move(output),
138 std::move(dump), std::move(knobs));
139}
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
Map::mapped_type & getOrError(Map &map, Key const &key, std::string const &error_message)
Definition Algorithm.h:113
std::string extractBaseNameWithoutExtension(std::string const &pathname)
std::string joinPaths(std::string const &pathA, std::string const &pathB)
Mesh & findMeshByName(std::vector< std::unique_ptr< Mesh > > const &meshes, std::string_view const name)
Definition Mesh.cpp:363
std::string parseDatabasePath(BaseLib::ConfigTree const &config)

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

◆ 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 44 of file CreateChemicalSolverInterface.cpp.

148{
149 auto mesh = *meshes[0];
150
151 auto const ls_name =
153 config.getConfigParameter<std::string>("linear_solver");
154 auto const& linear_solver = BaseLib::getOrError(
155 linear_solvers, ls_name,
156 "A linear solver with the given name does not exist.");
157
158 auto path_to_database = parseDatabasePath(config);
159
160 // TODO (renchao): remove mapping process id to component name.
161 std::vector<std::pair<int, std::string>> process_id_to_component_name_map;
162 // solution
163 auto aqueous_solution = PhreeqcKernelData::createAqueousSolution(
165 config.getConfigSubtree("solution"),
166 process_id_to_component_name_map);
167
168 // kinetic reactants
169 auto kinetic_reactants = PhreeqcKernelData::createKineticReactants(
171 config.getConfigSubtreeOptional("kinetic_reactants"), mesh);
172
173 // rates
174 auto reaction_rates = createReactionRates<PhreeqcKernelData::ReactionRate>(
176 config.getConfigSubtreeOptional("rates"));
177
178 // equilibrium reactants
179 auto equilibrium_reactants = PhreeqcKernelData::createEquilibriumReactants(
181 config.getConfigSubtreeOptional("equilibrium_reactants"), mesh);
182
183 return std::make_unique<PhreeqcKernelData::PhreeqcKernel>(
184 mesh, *linear_solver, mesh.computeNumberOfBaseNodes(),
185 process_id_to_component_name_map, std::move(path_to_database),
186 std::move(aqueous_solution), std::move(equilibrium_reactants),
187 std::move(kinetic_reactants), std::move(reaction_rates));
188}

◆ 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 44 of file CreateChemicalSolverInterface.cpp.

197{
198 auto mesh_name =
200 config.getConfigParameter<std::string>("mesh");
201
202 // Find and extract mesh from the list of meshes.
203 auto const& mesh = MeshLib::findMeshByName(meshes, mesh_name);
204
205 assert(mesh.getID() != 0);
206 DBUG("Found mesh '{:s}' with id {:d}.", mesh.getName(), mesh.getID());
207
208 auto const ls_name =
210 config.getConfigParameter<std::string>("linear_solver");
211 auto const& linear_solver = BaseLib::getOrError(
212 linear_solvers, ls_name,
213 "A linear solver with the given name does not exist.");
214
215 auto chemical_reaction_data =
216 SelfContainedSolverData::createChemicalReactionData(
218 config.getConfigSubtree("chemical_reactions"));
219
220 // create sparse stoichiometric matrix
221 std::vector<double> stoichiometric_matrix_vec;
222 for (auto const& per_chemical_reaction_data : chemical_reaction_data)
223 {
224 stoichiometric_matrix_vec.insert(
225 stoichiometric_matrix_vec.end(),
226 per_chemical_reaction_data->stoichiometric_vector.begin(),
227 per_chemical_reaction_data->stoichiometric_vector.end());
228 }
229
230 auto const num_components =
232 config.getConfigParameter<int>("number_of_components");
233
234 Eigen::MatrixXd const stoichiometric_matrix = MathLib::toMatrix(
235 stoichiometric_matrix_vec, num_components, num_components);
236
237 Eigen::SparseMatrix<double> sparse_stoichiometric_matrix;
238 sparse_stoichiometric_matrix = stoichiometric_matrix.sparseView();
239
240 return std::make_unique<SelfContainedSolverData::SelfContainedSolver>(
241 mesh, *linear_solver, sparse_stoichiometric_matrix,
242 std::move(chemical_reaction_data));
243}
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 20 of file CreateReactionRate.cpp.

22{
23 if (!config)
24 {
25 return {};
26 }
27
28 std::vector<ReactionRate> reaction_rates;
29 for (auto const& rate_config :
31 config->getConfigSubtreeList("rate"))
32 {
33 auto kinetic_reactant =
35 rate_config.getConfigParameter<std::string>("kinetic_reactant");
36
37 auto const expression_config =
39 rate_config.getConfigSubtree("expression");
40 auto const& statements =
42 expression_config.getConfigParameterList<std::string>("statement");
43
44 std::vector<std::string> expression_statements;
45 expression_statements.reserve(statements.size());
46 std::copy(begin(statements),
47 end(statements),
48 back_inserter(expression_statements));
49
50 reaction_rates.emplace_back(std::move(kinetic_reactant),
51 std::move(expression_statements));
52 }
53
54 return reaction_rates;
55}

◆ createReactionRates< PhreeqcIOData::ReactionRate >()

◆ createReactionRates< PhreeqcKernelData::ReactionRate >()