CreateChemicalSolverInterface.cpp

Go to the documentation of this file.

 
#include "CreateChemicalSolverInterface.h"
 
#include <phreeqcpp/cxxKinetics.h>
 
#include "BaseLib/ConfigTree.h"
#include "BaseLib/FileTools.h"
#include "Common/CreateReactionRate.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"
 
namespace
{
std::string parseDatabasePath(BaseLib::ConfigTree const& config)
{
    // database
    auto const database = config.getConfigParameter<std::string>("database");
    auto path_to_database =
        BaseLib::joinPaths(BaseLib::getProjectDirectory(), 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 = *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 =
        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 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 surface = PhreeqcIOData::createSurface(
        config.getConfigSubtreeOptional("surface"));
 
    // exchange
    auto const& exchangers = chemical_system->exchangers;
 
    // dump
    auto const project_file_name = BaseLib::joinPaths(
        output_directory,
        BaseLib::extractBaseNameWithoutExtension(config.getProjectFileName()));
 
    if (!surface.empty() && !exchangers.empty())
    {
        OGS_FATAL(
            "Using surface and exchange reactions simultaneously is not "
            "supported at the moment");
    }
 
    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);
 
    return std::make_unique<PhreeqcIOData::PhreeqcIO>(
        *linear_solver, std::move(project_file_name),
        std::move(path_to_database), std::move(chemical_system),
        std::move(reaction_rates), std::move(surface), std::move(user_punch),
        std::move(output), std::move(dump), std::move(knobs));
}
 
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& 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>(
        *linear_solver, mesh.getNumberOfBaseNodes(),
        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));
}
}  // namespace ChemistryLib