ChemistryLib Namespace Reference

Namespaces
namespace	PhreeqcIOData
namespace	PhreeqcKernelData
namespace	SelfContainedSolverData

Classes
class	ChemicalSolverInterface

Enumerations
enum class	ChemicalSolver { Phreeqc , PhreeqcKernel , SelfContained }
enum class	ChargeBalance { pH , pe , Unspecified }

Functions
ChargeBalance	createChargeBalance (BaseLib::ConfigTree const &config)
template<typename ReactionRate>
std::vector< ReactionRate >	createReactionRates (std::optional< BaseLib::ConfigTree > const &config)
template std::vector< PhreeqcIOData::ReactionRate >	createReactionRates< PhreeqcIOData::ReactionRate > (std::optional< BaseLib::ConfigTree > const &config)
template std::vector< PhreeqcKernelData::ReactionRate >	createReactionRates< PhreeqcKernelData::ReactionRate > (std::optional< BaseLib::ConfigTree > const &config)
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)
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 &)
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 &)
template<ChemicalSolver chemical_solver>
std::unique_ptr< ChemicalSolverInterface >	createChemicalSolverInterface (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)

Enumeration Type Documentation

ChargeBalance

enum class ChemistryLib::ChargeBalance

strong

Enumerator
pH
pe
Unspecified

Definition at line 8 of file ChargeBalance.h.

{
    pH,
    pe,
    Unspecified
};

ChemicalSolver

enum class ChemistryLib::ChemicalSolver

strong

Enumerator
Phreeqc
PhreeqcKernel
SelfContained

Definition at line 8 of file ChemicalSolverType.h.

{
    Phreeqc,
    PhreeqcKernel,
    SelfContained
};

Function Documentation

createChargeBalance()

ChargeBalance ChemistryLib::createChargeBalance

(

BaseLib::ConfigTree const &

config

)

Input File Parameter

prj__chemical_system__solution__charge_balance

Definition at line 12 of file CreateChargeBalance.cpp.

{
    auto const charge_balance =
        config.getConfigParameterOptional<std::string>("charge_balance");
 
    if (!charge_balance)
    {
        return ChargeBalance::Unspecified;
    }
    if (*charge_balance == "pH")
    {
        return ChargeBalance::pH;
    }
    if (*charge_balance == "pe")
    {
        return ChargeBalance::pe;
    }
 
    OGS_FATAL(
        "Error in specifying the way in which charge balance is achieved.  "
        "Adjusting pH value or pe value are supported at this moment.");
}

References BaseLib::ConfigTree::getConfigParameterOptional(), OGS_FATAL, pe, pH, and Unspecified.

Referenced by ChemistryLib::PhreeqcIOData::createAqueousSolution(), and ChemistryLib::PhreeqcKernelData::createInitialAqueousSolution().

createChemicalSolverInterface()

template<ChemicalSolver chemical_solver>

std::unique_ptr< ChemicalSolverInterface > ChemistryLib::createChemicalSolverInterface	(	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 )

Referenced by ProjectData::parseChemicalSolverInterface().

createChemicalSolverInterface< ChemicalSolver::Phreeqc >()

template<>

std::unique_ptr< ChemicalSolverInterface > ChemistryLib::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 )

Input File Parameter

prj__chemical_system__mesh

Input File Parameter

prj__chemical_system__linear_solver

Input File Parameter

prj__chemical_system__rates

Input File Parameter

prj__chemical_system__use_stream_for_data_exchange

Input File Parameter

prj__chemical_system__chemistry_threads

Input File Parameter

prj__chemical_system__concentration_warning_threshold

Input File Parameter

prj__chemical_system__knobs

Input File Parameter

prj__chemical_system__user_punch

Input File Parameter

prj__chemical_system__use_high_precision

Definition at line 39 of file CreateChemicalSolverInterface.cpp.

{
    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);
}

References BaseLib::ConfigTree::getConfigParameter(), INFO(), BaseLib::IsFileExisting(), BaseLib::joinPaths(), OGS_FATAL, and BaseLib::ConfigTree::projectDirectory().

createChemicalSolverInterface< ChemicalSolver::PhreeqcKernel >()

template<>

std::unique_ptr< ChemicalSolverInterface > ChemistryLib::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 &	)

Input File Parameter

prj__chemical_system__linear_solver

Input File Parameter

prj__chemical_system__solution

Input File Parameter

prj__chemical_system__kinetic_reactants

Input File Parameter

prj__chemical_system__rates

Input File Parameter

prj__chemical_system__equilibrium_reactants

Definition at line 39 of file CreateChemicalSolverInterface.cpp.

{
    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));
}

createChemicalSolverInterface< ChemicalSolver::SelfContained >()

template<>

std::unique_ptr< ChemicalSolverInterface > ChemistryLib::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 &	)

Input File Parameter

prj__chemical_system__mesh

Input File Parameter

prj__chemical_system__linear_solver

Input File Parameter

prj__chemical_system__chemical_reactions

Input File Parameter

prj__chemical_system__number_of_components

Definition at line 39 of file CreateChemicalSolverInterface.cpp.

{
    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));
}

createReactionRates()

template<typename ReactionRate>

std::vector< ReactionRate > ChemistryLib::createReactionRates

(

std::optional< BaseLib::ConfigTree > const &

config

)

Input File Parameter

prj__chemical_system__rates__rate

Input File Parameter

prj__chemical_system__rates__rate__kinetic_reactant

Input File Parameter

prj__chemical_system__rates__rate__expression

Input File Parameter

prj__chemical_system__rates__rate__expression__statement

Definition at line 13 of file CreateReactionRate.cpp.

{
    if (!config)
    {
        return {};
    }
 
    std::vector<ReactionRate> reaction_rates;
    for (auto const& rate_config :
         config->getConfigSubtreeList("rate"))
    {
        auto kinetic_reactant =
            rate_config.getConfigParameter<std::string>("kinetic_reactant");
 
        auto const expression_config =
            rate_config.getConfigSubtree("expression");
        auto const& statements =
            expression_config.getConfigParameterList<std::string>("statement");
 
        std::vector<std::string> expression_statements;
        expression_statements.reserve(statements.size());
        std::copy(begin(statements),
                  end(statements),
                  back_inserter(expression_statements));
 
        reaction_rates.emplace_back(std::move(kinetic_reactant),
                                    std::move(expression_statements));
    }
 
    return reaction_rates;
}

createReactionRates< PhreeqcIOData::ReactionRate >()

template std::vector< PhreeqcIOData::ReactionRate > ChemistryLib::createReactionRates< PhreeqcIOData::ReactionRate >

(

std::optional< BaseLib::ConfigTree > const &

config

)

createReactionRates< PhreeqcKernelData::ReactionRate >()

template std::vector< PhreeqcKernelData::ReactionRate > ChemistryLib::createReactionRates< PhreeqcKernelData::ReactionRate >

(

std::optional< BaseLib::ConfigTree > const &

config

)