CreateChemicalSolverInterface.cpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#include "CreateChemicalSolverInterface.h"
 
#include <phreeqcpp/cxxKinetics.h>
 
#include "BaseLib/ConfigTree.h"
#include "BaseLib/Error.h"
#include "BaseLib/FileTools.h"
#include "BaseLib/OgsChemThreads.h"
#include "Common/CreateReactionRate.h"
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
#include "MeshLib/Mesh.h"
#include "PhreeqcIO.h"
#include "PhreeqcIOData/ChemicalSystem.h"
#include "PhreeqcIOData/CreateChemicalSystem.h"
#include "PhreeqcIOData/CreateExchange.h"
#include "PhreeqcIOData/CreateKnobs.h"
#include "PhreeqcIOData/CreateOutput.h"
#include "PhreeqcIOData/CreateSurface.h"
#include "PhreeqcIOData/CreateUserPunch.h"
#include "PhreeqcIOData/Dump.h"
#include "PhreeqcIOData/Knobs.h"
#include "PhreeqcIOData/ReactionRate.h"
#include "PhreeqcIOData/Surface.h"
#include "PhreeqcIOData/UserPunch.h"
#include "PhreeqcKernel.h"
#include "PhreeqcKernelData/AqueousSolution.h"
#include "PhreeqcKernelData/CreateAqueousSolution.h"
#include "PhreeqcKernelData/CreateEquilibriumReactants.h"
#include "PhreeqcKernelData/CreateKineticReactant.h"
#include "PhreeqcKernelData/ReactionRate.h"
#include "SelfContainedSolverData/CreateChemicalReactionData.h"
#include "SelfContainedSolverData/SelfContainedSolver.h"
 
namespace
{
std::string parseDatabasePath(BaseLib::ConfigTree const& config)
{
    // database
    auto const database = config.getConfigParameter<std::string>("database");
    auto path_to_database =
        BaseLib::joinPaths(config.projectDirectory().string(), database);
 
    if (!BaseLib::IsFileExisting(path_to_database))
    {
        OGS_FATAL("Not found the specified thermodynamicdatabase: {:s}",
                  path_to_database);
    }
 
    INFO("Found the specified thermodynamic database: {:s}", path_to_database);
 
    return path_to_database;
}
}  // namespace
 
namespace ChemistryLib
{
template <>
std::unique_ptr<ChemicalSolverInterface>
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)
{
    auto mesh_name =
        config.getConfigParameter<std::string>("mesh");
 
    // Find and extract mesh from the list of meshes.
    auto const& mesh = MeshLib::findMeshByName(meshes, mesh_name);
 
    assert(mesh.getID() != 0);
    DBUG("Found mesh '{:s}' with id {:d}.", mesh.getName(), mesh.getID());
 
    auto const ls_name =
        config.getConfigParameter<std::string>("linear_solver");
    auto const& linear_solver = BaseLib::getOrError(
        linear_solvers, ls_name,
        "A linear solver with the given name does not exist.");
 
    auto path_to_database = parseDatabasePath(config);
 
    // chemical system
    auto chemical_system =
        PhreeqcIOData::createChemicalSystem(config, *meshes[0]);
 
    // rates
    auto reaction_rates = createReactionRates<PhreeqcIOData::ReactionRate>(
        config.getConfigSubtreeOptional("rates"));
 
    // surface
    auto const& surface = chemical_system->surface;
 
    // exchange
    auto const& exchangers = chemical_system->exchangers;
 
    // whether to use stream for data exchange
    auto use_stream_for_data_exchange =
        config.getConfigParameter<bool>("use_stream_for_data_exchange", true);
 
    // number of OpenMP threads for parallel chemistry execution
    auto const num_chemistry_threads_opt =
        config.getConfigParameterOptional<int>("chemistry_threads");
    if (num_chemistry_threads_opt && *num_chemistry_threads_opt <= 0)
    {
        OGS_FATAL(
            "<chemistry_threads> must be a positive integer, but got {:d}.",
            *num_chemistry_threads_opt);
    }
    int const num_chemistry_threads =
        num_chemistry_threads_opt ? *num_chemistry_threads_opt
                                  : BaseLib::getNumberOfChemistryThreads();
 
    // Negative (or zero) threshold below which clamped concentrations are
    // reported. All negative concentrations are clamped to zero unconditionally
    // (PHREEQC rejects negatives); this only sets the value below which an
    // aggregated warning is emitted. Default -1e-12 mol/kgw.
    auto const concentration_warning_threshold =
        config.getConfigParameter<double>("concentration_warning_threshold",
                                          -1e-12);
    if (concentration_warning_threshold > 0.0)
    {
        OGS_FATAL(
            "<concentration_warning_threshold> must not be positive (it is a "
            "threshold for negative concentrations), but got {:g}.",
            concentration_warning_threshold);
    }
 
    // dump
    auto const project_file_name =
        BaseLib::joinPaths(output_directory,
                           BaseLib::extractBaseNameWithoutExtension(
                               config.projectFilePath().string()));
    // old dump file is deleted if it exists
    auto dump = surface.empty() && exchangers.empty()
                    ? nullptr
                    : std::make_unique<PhreeqcIOData::Dump>(project_file_name);
 
    // knobs
    auto knobs = PhreeqcIOData::createKnobs(
        config.getConfigSubtree("knobs"));
 
    // user punch
    auto user_punch = PhreeqcIOData::createUserPunch(
        config.getConfigSubtreeOptional("user_punch"), *meshes[0]);
 
    // output
    auto const use_high_precision =
        config.getConfigParameter<bool>("use_high_precision", true);
    auto output = PhreeqcIOData::createOutput(
        *chemical_system, user_punch, use_high_precision, project_file_name);
 
    if (use_stream_for_data_exchange)
    {
        INFO("PhreeqcIO will use stringstream for data exchange.");
    }
    else
    {
        INFO("PhreeqcIO will use file for data exchange.");
    }
 
    return std::make_unique<PhreeqcIOData::PhreeqcIO>(
        mesh, *linear_solver, std::move(project_file_name),
        std::move(path_to_database), std::move(chemical_system),
        std::move(reaction_rates), std::move(user_punch), std::move(output),
        std::move(dump), std::move(knobs), use_stream_for_data_exchange,
        num_chemistry_threads, concentration_warning_threshold);
}
 
template <>
std::unique_ptr<ChemicalSolverInterface>
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& /*output_directory*/)
{
    auto mesh = *meshes[0];
 
    auto const ls_name =
        config.getConfigParameter<std::string>("linear_solver");
    auto const& linear_solver = BaseLib::getOrError(
        linear_solvers, ls_name,
        "A linear solver with the given name does not exist.");
 
    auto path_to_database = parseDatabasePath(config);
 
    // TODO (renchao): remove mapping process id to component name.
    std::vector<std::pair<int, std::string>> process_id_to_component_name_map;
    // solution
    auto aqueous_solution = PhreeqcKernelData::createAqueousSolution(
        config.getConfigSubtree("solution"),
        process_id_to_component_name_map);
 
    // kinetic reactants
    auto kinetic_reactants = PhreeqcKernelData::createKineticReactants(
        config.getConfigSubtreeOptional("kinetic_reactants"), mesh);
 
    // rates
    auto reaction_rates = createReactionRates<PhreeqcKernelData::ReactionRate>(
        config.getConfigSubtreeOptional("rates"));
 
    // equilibrium reactants
    auto equilibrium_reactants = PhreeqcKernelData::createEquilibriumReactants(
        config.getConfigSubtreeOptional("equilibrium_reactants"), mesh);
 
    return std::make_unique<PhreeqcKernelData::PhreeqcKernel>(
        mesh, *linear_solver, mesh.computeNumberOfBaseNodes(),
        process_id_to_component_name_map, std::move(path_to_database),
        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 =
        config.getConfigParameter<std::string>("mesh");
 
    // Find and extract mesh from the list of meshes.
    auto const& mesh = MeshLib::findMeshByName(meshes, mesh_name);
 
    assert(mesh.getID() != 0);
    DBUG("Found mesh '{:s}' with id {:d}.", mesh.getName(), mesh.getID());
 
    auto const ls_name =
        config.getConfigParameter<std::string>("linear_solver");
    auto const& linear_solver = BaseLib::getOrError(
        linear_solvers, ls_name,
        "A linear solver with the given name does not exist.");
 
    auto chemical_reaction_data =
        SelfContainedSolverData::createChemicalReactionData(
            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 =
        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