From f9a7ef1e360f14f3b58f42e21c244f5350fac474 Mon Sep 17 00:00:00 2001
From: renchao_lu <renchao.lu@gmail.com>
Date: Sun, 10 Jul 2022 15:52:55 +0200
Subject: [PATCH] [CL] add self-contained chemical solver.
---
Applications/ApplicationsLib/ProjectData.cpp | 10 +++
ChemistryLib/CMakeLists.txt | 1 +
ChemistryLib/ChemicalSolverType.h | 3 +-
.../CreateChemicalSolverInterface.cpp | 65 +++++++++++++++++++
.../SelfContainedSolver.h | 53 +++++++++++++++
5 files changed, 131 insertions(+), 1 deletion(-)
create mode 100644 ChemistryLib/SelfContainedSolverData/SelfContainedSolver.h
diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index b0bf49b8158..5446b0fa99b 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -552,6 +552,16 @@ ProjectData::parseChemicalSolverInterface(
"the time being. Please set 'Phreeqc'' as the chemical solver for "
"reactive transport modeling.");
}
+ else if (boost::iequals(chemical_solver, "SelfContained"))
+ {
+ INFO(
+ "Use self-contained chemical solver for water chemistry "
+ "calculation.");
+
+ chemical_solver_interface = ChemistryLib::createChemicalSolverInterface<
+ ChemistryLib::ChemicalSolver::SelfContained>(
+ _mesh_vec, _linear_solvers, *config, output_directory);
+ }
else
{
OGS_FATAL(
diff --git a/ChemistryLib/CMakeLists.txt b/ChemistryLib/CMakeLists.txt
index 1c0f7a880c4..7ee84dc15ce 100644
--- a/ChemistryLib/CMakeLists.txt
+++ b/ChemistryLib/CMakeLists.txt
@@ -2,6 +2,7 @@
get_source_files(SOURCES)
append_source_files(SOURCES PhreeqcIOData)
append_source_files(SOURCES PhreeqcKernelData)
+append_source_files(SOURCES SelfContainedSolverData)
append_source_files(SOURCES Common)
# Create the library
diff --git a/ChemistryLib/ChemicalSolverType.h b/ChemistryLib/ChemicalSolverType.h
index 5448c2c2c42..4a34925d754 100644
--- a/ChemistryLib/ChemicalSolverType.h
+++ b/ChemistryLib/ChemicalSolverType.h
@@ -15,6 +15,7 @@ namespace ChemistryLib
enum class ChemicalSolver
{
Phreeqc,
- PhreeqcKernel
+ PhreeqcKernel,
+ SelfContained
};
} // namespace ChemistryLib
diff --git a/ChemistryLib/CreateChemicalSolverInterface.cpp b/ChemistryLib/CreateChemicalSolverInterface.cpp
index 0b72c904eb0..87058f69db3 100644
--- a/ChemistryLib/CreateChemicalSolverInterface.cpp
+++ b/ChemistryLib/CreateChemicalSolverInterface.cpp
@@ -15,6 +15,7 @@
#include "BaseLib/ConfigTree.h"
#include "BaseLib/FileTools.h"
#include "Common/CreateReactionRate.h"
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
#include "MeshLib/Mesh.h"
#include "PhreeqcIO.h"
#include "PhreeqcIOData/ChemicalSystem.h"
@@ -35,6 +36,8 @@
#include "PhreeqcKernelData/CreateEquilibriumReactants.h"
#include "PhreeqcKernelData/CreateKineticReactant.h"
#include "PhreeqcKernelData/ReactionRate.h"
+#include "SelfContainedSolverData/CreateChemicalReactionData.h"
+#include "SelfContainedSolverData/SelfContainedSolver.h"
namespace
{
@@ -190,4 +193,66 @@ createChemicalSolverInterface<ChemicalSolver::PhreeqcKernel>(
std::move(aqueous_solution), std::move(equilibrium_reactants),
std::move(kinetic_reactants), std::move(reaction_rates));
}
+
+template <>
+std::unique_ptr<ChemicalSolverInterface>
+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& /*output_directory*/)
+{
+ auto mesh_name =
+ //! \ogs_file_param{prj__chemical_system__mesh}
+ config.getConfigParameter<std::string>("mesh");
+
+ // Find and extract mesh from the list of meshes.
+ auto const& mesh = *BaseLib::findElementOrError(
+ std::begin(meshes), std::end(meshes),
+ [&mesh_name](auto const& mesh)
+ {
+ assert(mesh != nullptr);
+ return mesh->getName() == mesh_name;
+ },
+ "Required mesh with name '" + mesh_name + "' not found.");
+
+ assert(mesh.getID() != 0);
+ DBUG("Found mesh '{:s}' with id {:d}.", mesh.getName(), mesh.getID());
+
+ auto const ls_name =
+ //! \ogs_file_param{prj__chemical_system__linear_solver}
+ config.getConfigParameter<std::string>("linear_solver");
+ auto& linear_solver = BaseLib::getOrError(
+ linear_solvers, ls_name,
+ "A linear solver with the given name does not exist.");
+
+ auto chemical_reaction_data =
+ SelfContainedSolverData::createChemicalReactionData(
+ //! \ogs_file_param{prj__chemical_system__chemical_reactions}
+ config.getConfigSubtree("chemical_reactions"));
+
+ // create sparse stoichiometric matrix
+ std::vector<double> stoichiometric_matrix_vec;
+ for (auto const& per_chemical_reaction_data : chemical_reaction_data)
+ {
+ stoichiometric_matrix_vec.insert(
+ stoichiometric_matrix_vec.end(),
+ per_chemical_reaction_data->stoichiometric_vector.begin(),
+ per_chemical_reaction_data->stoichiometric_vector.end());
+ }
+
+ auto const num_components =
+ //! \ogs_file_param{prj__chemical_system__number_of_components}
+ config.getConfigParameter<int>("number_of_components");
+
+ Eigen::MatrixXd const stoichiometric_matrix = MathLib::toMatrix(
+ stoichiometric_matrix_vec, num_components, num_components);
+
+ Eigen::SparseMatrix<double> sparse_stoichiometric_matrix;
+ sparse_stoichiometric_matrix = stoichiometric_matrix.sparseView();
+
+ return std::make_unique<SelfContainedSolverData::SelfContainedSolver>(
+ mesh, *linear_solver, sparse_stoichiometric_matrix,
+ std::move(chemical_reaction_data));
+}
} // namespace ChemistryLib
diff --git a/ChemistryLib/SelfContainedSolverData/SelfContainedSolver.h b/ChemistryLib/SelfContainedSolverData/SelfContainedSolver.h
new file mode 100644
index 00000000000..51876fe33cd
--- /dev/null
+++ b/ChemistryLib/SelfContainedSolverData/SelfContainedSolver.h
@@ -0,0 +1,53 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2022, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <Eigen/Sparse>
+#include <memory>
+
+#include "ChemicalReaction.h"
+#include "ChemistryLib/ChemicalSolverInterface.h"
+
+namespace ChemistryLib
+{
+namespace SelfContainedSolverData
+{
+class SelfContainedSolver final : public ChemicalSolverInterface
+{
+public:
+ SelfContainedSolver(MeshLib::Mesh const& mesh,
+ GlobalLinearSolver& linear_solver,
+ Eigen::SparseMatrix<double>
+ stoichiometric_matrix,
+ std::vector<std::unique_ptr<ChemicalReaction>>
+ chemical_reactions)
+ : ChemicalSolverInterface(mesh, linear_solver),
+ _stoichiometric_matrix(stoichiometric_matrix),
+ _chemical_reactions(std::move(chemical_reactions))
+ {
+ }
+
+ Eigen::SparseMatrix<double> const* getStoichiometricMatrix() const override
+ {
+ return &_stoichiometric_matrix;
+ }
+
+ double getKineticPrefactor(unsigned int reaction_id) const override
+ {
+ return _chemical_reactions[reaction_id]->getKineticPrefactor();
+ }
+
+private:
+ Eigen::SparseMatrix<double> _stoichiometric_matrix;
+ std::vector<std::unique_ptr<ChemicalReaction>> _chemical_reactions;
+};
+} // namespace SelfContainedSolverData
+} // namespace ChemistryLib
--
GitLab