diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index de6d45adaa8c46c7ad637ac92f0e7d2bf29ae164..eda8957dedad4d50459ba217404e26cbd5b9e82b 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -38,6 +38,10 @@
#include "ProcessLib/CreateJacobianAssembler.h"
#include "ProcessLib/DeactivatedSubdomain.h"
+// PhreeqcIO
+#include "ChemistryLib/CreatePhreeqcIO.h"
+#include "ProcessLib/ComponentTransport/ComponentTransportProcess.h"
+
// FileIO
#include "GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.h"
#include "MeshLib/IO/readMeshFromFile.h"
@@ -338,6 +342,11 @@ ProjectData::ProjectData(BaseLib::ConfigTree const& project_config,
//! \ogs_file_param{prj__nonlinear_solvers}
parseNonlinearSolvers(project_config.getConfigSubtree("nonlinear_solvers"));
+ parseChemicalSystem(
+ //! \ogs_file_param{prj__chemical_system}
+ project_config.getConfigSubtreeOptional("chemical_system"),
+ output_directory);
+
//! \ogs_file_param{prj__time_loop}
parseTimeLoop(project_config.getConfigSubtree("time_loop"),
output_directory);
@@ -923,8 +932,9 @@ void ProjectData::parseTimeLoop(BaseLib::ConfigTree const& config,
{
DBUG("Reading time loop configuration.");
- _time_loop = ProcessLib::createTimeLoop(
- config, output_directory, _processes, _nonlinear_solvers, _mesh_vec);
+ _time_loop = ProcessLib::createTimeLoop(config, output_directory,
+ _processes, _nonlinear_solvers,
+ _mesh_vec, _chemical_system);
if (!_time_loop)
{
@@ -996,3 +1006,52 @@ void ProjectData::parseCurves(
"The curve name is not unique.");
}
}
+
+void ProjectData::parseChemicalSystem(
+ boost::optional<BaseLib::ConfigTree> const& config,
+ std::string const& output_directory)
+{
+ if (!config)
+ return;
+
+ INFO(
+ "Ready for initializing interface to a chemical solver for water "
+ "chemistry calculation.");
+
+ auto const& process = _processes.begin()->second;
+ if (auto const* component_transport_process = dynamic_cast<
+ ProcessLib::ComponentTransport::ComponentTransportProcess const*>(
+ process.get()))
+ {
+ auto const chemical_solver =
+ //! \ogs_file_attr{prj__chemical_system__chemical_solver}
+ config->getConfigAttribute<std::string>("chemical_solver");
+
+ if (boost::iequals(chemical_solver, "Phreeqc"))
+ {
+ INFO(
+ "Configuring phreeqc interface for water chemistry "
+ "calculation.");
+
+ auto const& process_id_to_component_name_map =
+ component_transport_process->getProcessIDToComponentNameMap();
+
+ _chemical_system = ChemistryLib::createPhreeqcIO(
+ _mesh_vec[0]->getNumberOfBaseNodes(),
+ process_id_to_component_name_map, config, output_directory);
+ }
+ else
+ {
+ OGS_FATAL(
+ "Unknown chemical solver. Please specify Phreeqc as the solver "
+ "for water chemistry calculation instead.");
+ }
+ }
+ else
+ {
+ OGS_FATAL(
+ "The specified type of the process to be solved is not component "
+ "transport process so that water chemistry calculation could not "
+ "be activated.");
+ }
+}
diff --git a/Applications/ApplicationsLib/ProjectData.h b/Applications/ApplicationsLib/ProjectData.h
index b319c3341ff2c285030afd1bccaccea57e29c29e..a488bc280e9965d5d7710aebe2d28c037f2a872a 100644
--- a/Applications/ApplicationsLib/ProjectData.h
+++ b/Applications/ApplicationsLib/ProjectData.h
@@ -21,6 +21,7 @@
#include "MaterialLib/MPL/Medium.h"
#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
+#include "ChemistryLib/PhreeqcIO.h"
#include "ParameterLib/CoordinateSystem.h"
#include "ParameterLib/Parameter.h"
#include "ProcessLib/Process.h"
@@ -112,6 +113,9 @@ private:
void parseCurves(boost::optional<BaseLib::ConfigTree> const& config);
+ void parseChemicalSystem(boost::optional<BaseLib::ConfigTree> const& config,
+ const std::string& output_directory);
+
std::unique_ptr<MaterialPropertyLib::Medium> _medium;
std::vector<std::unique_ptr<MeshLib::Mesh>> _mesh_vec;
std::map<std::string, std::unique_ptr<ProcessLib::Process>> _processes;
@@ -134,4 +138,6 @@ private:
std::map<std::string,
std::unique_ptr<MathLib::PiecewiseLinearInterpolation>>
_curves;
+
+ std::unique_ptr<ChemistryLib::PhreeqcIO> _chemical_system;
};
diff --git a/ProcessLib/CreateTimeLoop.cpp b/ProcessLib/CreateTimeLoop.cpp
index 9a0de21c25335c649a53a6a9f3bd8f242d6edb09..71bbe024db3f7efd8e962efeec52cf03a5821b12 100644
--- a/ProcessLib/CreateTimeLoop.cpp
+++ b/ProcessLib/CreateTimeLoop.cpp
@@ -25,7 +25,8 @@ std::unique_ptr<TimeLoop> createTimeLoop(
const std::map<std::string, std::unique_ptr<Process>>& processes,
const std::map<std::string, std::unique_ptr<NumLib::NonlinearSolverBase>>&
nonlinear_solvers,
- std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes)
+ std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes,
+ std::unique_ptr<ChemistryLib::PhreeqcIO>& phreeqc_io)
{
auto const& coupling_config
//! \ogs_file_param{prj__time_loop__global_process_coupling}
@@ -93,6 +94,7 @@ std::unique_ptr<TimeLoop> createTimeLoop(
return std::make_unique<TimeLoop>(
std::move(output), std::move(per_process_data), max_coupling_iterations,
- std::move(global_coupling_conv_criteria), start_time, end_time);
+ std::move(global_coupling_conv_criteria), std::move(phreeqc_io),
+ start_time, end_time);
}
} // namespace ProcessLib
diff --git a/ProcessLib/CreateTimeLoop.h b/ProcessLib/CreateTimeLoop.h
index 86f6d9d89316e045d4b6210da5114a0f5620c647..167e0e68ccbdff2b6f2759890989b7c46b5233cb 100644
--- a/ProcessLib/CreateTimeLoop.h
+++ b/ProcessLib/CreateTimeLoop.h
@@ -31,6 +31,11 @@ namespace NumLib
class NonlinearSolverBase;
}
+namespace ChemistryLib
+{
+class PhreeqcIO;
+}
+
namespace ProcessLib
{
class TimeLoop;
@@ -45,6 +50,7 @@ std::unique_ptr<TimeLoop> createTimeLoop(
std::map<std::string, std::unique_ptr<Process>> const& processes,
std::map<std::string, std::unique_ptr<NumLib::NonlinearSolverBase>> const&
nonlinear_solvers,
- std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes);
+ std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes,
+ std::unique_ptr<ChemistryLib::PhreeqcIO>& phreeqc_io);
} // namespace ProcessLib
diff --git a/ProcessLib/TimeLoop.cpp b/ProcessLib/TimeLoop.cpp
index 116280db9a04926fb299ebefe9a9c07a7de39e48..efd338064c41049e84609ea894d55f7a48ca85da 100644
--- a/ProcessLib/TimeLoop.cpp
+++ b/ProcessLib/TimeLoop.cpp
@@ -13,9 +13,12 @@
#include "BaseLib/Error.h"
#include "BaseLib/RunTime.h"
+#include "ChemistryLib/CreatePhreeqcIO.h"
#include "MathLib/LinAlg/LinAlg.h"
#include "NumLib/ODESolver/ConvergenceCriterionPerComponent.h"
#include "NumLib/ODESolver/TimeDiscretizedODESystem.h"
+#include "ProcessLib/CreateProcessData.h"
+#include "ProcessLib/Output/CreateOutput.h"
#include "CoupledSolutionsForStaggeredScheme.h"
#include "ProcessData.h"
@@ -207,13 +210,15 @@ TimeLoop::TimeLoop(std::unique_ptr<Output>&& output,
const int global_coupling_max_iterations,
std::vector<std::unique_ptr<NumLib::ConvergenceCriterion>>&&
global_coupling_conv_crit,
+ std::unique_ptr<ChemistryLib::PhreeqcIO>&& chemical_system,
const double start_time, const double end_time)
: _output(std::move(output)),
_per_process_data(std::move(per_process_data)),
_start_time(start_time),
_end_time(end_time),
_global_coupling_max_iterations(global_coupling_max_iterations),
- _global_coupling_conv_crit(std::move(global_coupling_conv_crit))
+ _global_coupling_conv_crit(std::move(global_coupling_conv_crit)),
+ _chemical_system(std::move(chemical_system))
{
}
diff --git a/ProcessLib/TimeLoop.h b/ProcessLib/TimeLoop.h
index 1aca083d9d4d028e538e078e0faad20aea92b6d3..3597dad777aa5fce738434303a4df8432b5a8d51 100644
--- a/ProcessLib/TimeLoop.h
+++ b/ProcessLib/TimeLoop.h
@@ -27,6 +27,11 @@ namespace NumLib
class ConvergenceCriterion;
}
+namespace ChemistryLib
+{
+class PhreeqcIO;
+}
+
namespace ProcessLib
{
struct ProcessData;
@@ -40,6 +45,7 @@ public:
const int global_coupling_max_iterations,
std::vector<std::unique_ptr<NumLib::ConvergenceCriterion>>&&
global_coupling_conv_crit,
+ std::unique_ptr<ChemistryLib::PhreeqcIO>&& chemical_system,
const double start_time, const double end_time);
bool loop();
@@ -124,6 +130,7 @@ private:
std::vector<std::unique_ptr<NumLib::ConvergenceCriterion>>
_global_coupling_conv_crit;
+ std::unique_ptr<ChemistryLib::PhreeqcIO> _chemical_system;
/**
* Vector of solutions of the coupled processes.
* Each vector element stores the references of the solution vectors