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__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 // dump
102 auto const project_file_name =
103 BaseLib::joinPaths(output_directory,
105 config.projectFilePath().string()));
106 // old dump file is deleted if it exists
107 auto dump = surface.empty() && exchangers.empty()
108 ? nullptr
109 : std::make_unique<PhreeqcIOData::Dump>(project_file_name);
110
111 // knobs
112 auto knobs = PhreeqcIOData::createKnobs(
114 config.getConfigSubtree("knobs"));
115
116 // user punch
117 auto user_punch = PhreeqcIOData::createUserPunch(
119 config.getConfigSubtreeOptional("user_punch"), *meshes[0]);
120
121 // output
122 auto const use_high_precision =
124 config.getConfigParameter<bool>("use_high_precision", true);
125 auto output = PhreeqcIOData::createOutput(
126 *chemical_system, user_punch, use_high_precision, project_file_name);
127
128 return std::make_unique<PhreeqcIOData::PhreeqcIO>(
129 mesh, *linear_solver, std::move(project_file_name),
130 std::move(path_to_database), std::move(chemical_system),
131 std::move(reaction_rates), std::move(user_punch), std::move(output),
132 std::move(dump), std::move(knobs));
133}
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.

142{
143 auto mesh = *meshes[0];
144
145 auto const ls_name =
147 config.getConfigParameter<std::string>("linear_solver");
148 auto const& linear_solver = BaseLib::getOrError(
149 linear_solvers, ls_name,
150 "A linear solver with the given name does not exist.");
151
152 auto path_to_database = parseDatabasePath(config);
153
154 // TODO (renchao): remove mapping process id to component name.
155 std::vector<std::pair<int, std::string>> process_id_to_component_name_map;
156 // solution
157 auto aqueous_solution = PhreeqcKernelData::createAqueousSolution(
159 config.getConfigSubtree("solution"),
160 process_id_to_component_name_map);
161
162 // kinetic reactants
163 auto kinetic_reactants = PhreeqcKernelData::createKineticReactants(
165 config.getConfigSubtreeOptional("kinetic_reactants"), mesh);
166
167 // rates
170 config.getConfigSubtreeOptional("rates"));
171
172 // equilibrium reactants
173 auto equilibrium_reactants = PhreeqcKernelData::createEquilibriumReactants(
175 config.getConfigSubtreeOptional("equilibrium_reactants"), mesh);
176
177 return std::make_unique<PhreeqcKernelData::PhreeqcKernel>(
178 mesh, *linear_solver, mesh.computeNumberOfBaseNodes(),
179 process_id_to_component_name_map, std::move(path_to_database),
180 std::move(aqueous_solution), std::move(equilibrium_reactants),
181 std::move(kinetic_reactants), std::move(reaction_rates));
182}
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.

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