diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp index b0bf49b8158c5a15f6090f53e701672aa8679fb7..5446b0fa99ba2c0bd8c8627d0ae79eaf2cb4059e 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 1c0f7a880c48a8fed2e35316859fb095d2e12637..7ee84dc15ce19e6c509795d82f80f49da7f78571 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 5448c2c2c42a162e9b6c3551e43da10be3927e21..4a34925d754cab809035d5196287f526ebce02f5 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 0b72c904eb00c9efe40e6557a452c1cab328234e..87058f69db360e52bd5b843481ef9de60fc87240 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 0000000000000000000000000000000000000000..51876fe33cd8b70dd85e5841efe90d58df62f293 --- /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