Detailed Description

Definition at line 25 of file PhreeqcKernel.h.

#include <PhreeqcKernel.h>

Inheritance diagram for ChemistryLib::PhreeqcKernelData::PhreeqcKernel:

Collaboration diagram for ChemistryLib::PhreeqcKernelData::PhreeqcKernel:

Public Member Functions
	PhreeqcKernel (MeshLib::Mesh const &mesh, GlobalLinearSolver &linear_solver, std::size_t const num_chemical_systems, std::vector< std::pair< int, std::string > > const &process_id_to_component_name_map, std::string const &database, AqueousSolution aqueous_solution, std::unique_ptr< EquilibriumReactants > &&equilibrium_reactants, std::unique_ptr< Kinetics > &&kinetic_reactants, std::vector< ReactionRate > &&reaction_rates)
void	executeSpeciationCalculation (double const dt) override
void	setAqueousSolutions (std::vector< GlobalVector * > const &process_solutions)
void	callPhreeqc (std::vector< GlobalVector * > &process_solutions)
void	updateNodalProcessSolutions (std::vector< GlobalVector * > const &process_solutions, std::size_t const node_id)
Public Member Functions inherited from ChemistryLib::ChemicalSolverInterface
	ChemicalSolverInterface (MeshLib::Mesh const &mesh, GlobalLinearSolver &linear_solver_)
std::vector< std::size_t > const &	activeElementIDs () const
virtual void	initialize ()
virtual void	initializeChemicalSystemConcrete (std::vector< double > const &, GlobalIndexType const &, MaterialPropertyLib::Medium const &, ParameterLib::SpatialPosition const &, double const)
virtual void	setChemicalSystemConcrete (std::vector< double > const &, GlobalIndexType const &, MaterialPropertyLib::Medium const *, MaterialPropertyLib::VariableArray const &, ParameterLib::SpatialPosition const &, double const, double const)
virtual void	setAqueousSolutionsPrevFromDumpFile ()
virtual double	getConcentration (int const, GlobalIndexType const) const
virtual std::vector< std::string > const	getComponentList () const
virtual Eigen::SparseMatrix< double > const *	getStoichiometricMatrix () const
virtual double	getKineticPrefactor (std::size_t reaction_id) const
virtual void	updateVolumeFractionPostReaction (GlobalIndexType const &, MaterialPropertyLib::Medium const &, ParameterLib::SpatialPosition const &, double const, double const, double const)
virtual void	updatePorosityPostReaction (GlobalIndexType const &, MaterialPropertyLib::Medium const &, double &)
virtual void	computeSecondaryVariable (std::size_t const, std::vector< GlobalIndexType > const &)
virtual	~ChemicalSolverInterface ()=default

Private Member Functions
void	initializePhreeqcGeneralSettings ()
void	tidyEquilibriumReactants (EquilibriumReactants const &equilibrium_reactants)
void	loadDatabase (std::string const &database)
void	reinitializeRates ()
void	setConvergenceTolerance ()
void	configureOutputSettings ()
cxxISolution *	getOrCreateInitialAqueousSolution (cxxSolution &aqueous_solution)
void	setTimeStepSize (double const dt)
void	reset (std::size_t const chemical_system_id)

Private Attributes
std::map< int, struct master * >	_process_id_to_master_map
std::unique_ptr< cxxISolution const >	_initial_aqueous_solution
std::unique_ptr< cxxSolution const >	_aqueous_solution
std::vector< ReactionRate > const	_reaction_rates

Additional Inherited Members
Public Attributes inherited from ChemistryLib::ChemicalSolverInterface
std::vector< GlobalIndexType >	chemical_system_index_map
MeshLib::Mesh const &	_mesh
GlobalLinearSolver &	linear_solver

Constructor & Destructor Documentation

◆ PhreeqcKernel()

ChemistryLib::PhreeqcKernelData::PhreeqcKernel::PhreeqcKernel	(	MeshLib::Mesh const &	mesh,
		GlobalLinearSolver &	linear_solver,
		std::size_t const	num_chemical_systems,
		std::vector< std::pair< int, std::string > > const &	process_id_to_component_name_map,
		std::string const &	database,
		AqueousSolution	aqueous_solution,
		std::unique_ptr< EquilibriumReactants > &&	equilibrium_reactants,
		std::unique_ptr< Kinetics > &&	kinetic_reactants,
		std::vector< ReactionRate > &&	reaction_rates )

Definition at line 19 of file PhreeqcKernel.cpp.

    : ChemicalSolverInterface(mesh, linear_solver),
      _initial_aqueous_solution(aqueous_solution.getInitialAqueousSolution()),
      _aqueous_solution(aqueous_solution.castToBaseClassNoninitialized()),
      _reaction_rates(std::move(reaction_rates))
{
    initializePhreeqcGeneralSettings();
 
    loadDatabase(database);
 
    // solution
    for (std::size_t chemical_system_id = 0;
         chemical_system_id < num_chemical_systems;
         ++chemical_system_id)
    {
        aqueous_solution.setChemicalSystemID(chemical_system_id);
        Rxn_solution_map.emplace(chemical_system_id,
                                 *aqueous_solution.castToBaseClass());
    }
    use.Set_solution_in(true);
 
    // equilibrium reactants
    if (equilibrium_reactants)
    {
        tidyEquilibriumReactants(*equilibrium_reactants);
 
        for (std::size_t chemical_system_id = 0;
             chemical_system_id < num_chemical_systems;
             ++chemical_system_id)
        {
            equilibrium_reactants->setChemicalSystemID(chemical_system_id);
            Rxn_pp_assemblage_map.emplace(
                chemical_system_id, *equilibrium_reactants->castToBaseClass());
        }
        // explicitly release and delete equilibrium_reactants
        equilibrium_reactants.reset(nullptr);
 
        use.Set_pp_assemblage_in(true);
    }
 
    // kinetic reactants
    if (kinetic_reactants)
    {
        for (std::size_t chemical_system_id = 0;
             chemical_system_id < num_chemical_systems;
             ++chemical_system_id)
        {
            kinetic_reactants->setChemicalSystemID(chemical_system_id);
            Rxn_kinetics_map.emplace(chemical_system_id,
                                     *kinetic_reactants->castToBaseClass());
        }
        // explicitly release and delete kinetic_reactants
        kinetic_reactants.reset(nullptr);
 
        use.Set_kinetics_in(true);
    }
 
    // rates
    reinitializeRates();
 
    setConvergenceTolerance();
 
    configureOutputSettings();
 
    for (auto const& map_pair : process_id_to_component_name_map)
    {
        auto const transport_process_id = map_pair.first;
        auto const& transport_process_variable = map_pair.second;
 
        auto master_species =
            master_bsearch(transport_process_variable.c_str());
 
        _process_id_to_master_map[transport_process_id] = master_species;
    }
}

References ChemistryLib::ChemicalSolverInterface::ChemicalSolverInterface(), _aqueous_solution, _initial_aqueous_solution, _process_id_to_master_map, _reaction_rates, ChemistryLib::PhreeqcKernelData::AqueousSolution::castToBaseClass(), configureOutputSettings(), initializePhreeqcGeneralSettings(), ChemistryLib::ChemicalSolverInterface::linear_solver, loadDatabase(), reinitializeRates(), ChemistryLib::PhreeqcKernelData::AqueousSolution::setChemicalSystemID(), setConvergenceTolerance(), and tidyEquilibriumReactants().

Member Function Documentation

◆ callPhreeqc()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::callPhreeqc ( std::vector< GlobalVector * > & process_solutions )

Definition at line 257 of file PhreeqcKernel.cpp.

{
    if (process_solutions.empty())
    {
        OGS_FATAL("About to access an empty process solutions vector.");
    }
    std::size_t const num_chemical_systems = process_solutions[0]->size();
    for (std::size_t chemical_system_id = 0;
         chemical_system_id < num_chemical_systems;
         ++chemical_system_id)
    {
        Rxn_new_solution.insert(chemical_system_id);
        new_solution = 1;
        use.Set_n_solution_user(chemical_system_id);
 
        if (!Rxn_pp_assemblage_map.empty())
        {
            Rxn_new_pp_assemblage.insert(chemical_system_id);
            use.Set_pp_assemblage_in(true);
            use.Set_n_pp_assemblage_user(chemical_system_id);
        }
 
        if (!Rxn_kinetics_map.empty())
        {
            use.Set_kinetics_in(true);
            use.Set_n_kinetics_user(chemical_system_id);
        }
 
        initial_solutions(false);
 
        reactions();
 
        updateNodalProcessSolutions(process_solutions, chemical_system_id);
 
        reset(chemical_system_id);
    }
}

References OGS_FATAL, reset(), and updateNodalProcessSolutions().

Referenced by executeSpeciationCalculation().

◆ configureOutputSettings()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::configureOutputSettings ( )

inlineprivate

Definition at line 60 of file PhreeqcKernel.h.

60{ pr.all = false; }

Referenced by PhreeqcKernel().

◆ executeSpeciationCalculation()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::executeSpeciationCalculation ( double const dt )

overridevirtual

Run PHREEQC for all local chemical systems for the current timestep.

Uses the state previously provided by initializeChemicalSystemConcrete() / setChemicalSystemConcrete() for each chemical_system_id. Each local system is advanced as an isolated batch reactor over \(\Delta t\): no mass is exchanged between different systems inside PHREEQC.

PHREEQC returns, for each system:

updated aqueous composition,
reaction / source terms \(R_{\alpha}\) for each transported component,
updated mineral amounts and saturation state.

After this call completes, these reacted values are available for porosity update and for assembling the reaction/source term in the transport equation.

Parameters

dt	Timestep size \(\Delta t\) used for kinetic reactions.

Reimplemented from ChemistryLib::ChemicalSolverInterface.

Definition at line 167 of file PhreeqcKernel.cpp.

{
    std::vector<GlobalVector*> process_solutions;
 
    setAqueousSolutions(process_solutions);
 
    setTimeStepSize(dt);
 
    callPhreeqc(process_solutions);
}

References callPhreeqc(), setAqueousSolutions(), and setTimeStepSize().

◆ getOrCreateInitialAqueousSolution()

cxxISolution * ChemistryLib::PhreeqcKernelData::PhreeqcKernel::getOrCreateInitialAqueousSolution ( cxxSolution & aqueous_solution )

private

Definition at line 235 of file PhreeqcKernel.cpp.

{
    if (!aqueous_solution.Get_initial_data())
    {
        aqueous_solution.Set_initial_data(_initial_aqueous_solution.get());
        aqueous_solution.Set_new_def(true);
    }
 
    return aqueous_solution.Get_initial_data();
}

References _initial_aqueous_solution.

Referenced by setAqueousSolutions().

◆ initializePhreeqcGeneralSettings()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::initializePhreeqcGeneralSettings ( )

inlineprivate

Definition at line 49 of file PhreeqcKernel.h.

49{ do_initialize(); }

Referenced by PhreeqcKernel().

◆ loadDatabase()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::loadDatabase ( std::string const & database )

private

Definition at line 126 of file PhreeqcKernel.cpp.

{
    std::ifstream in(database);
    if (!in)
    {
        OGS_FATAL("Unable to open database file '{:s}'.", database);
    }
    assert(phrq_io->get_istream() == nullptr);
    phrq_io->push_istream(&in, false);
    read_database();
}

References OGS_FATAL.

Referenced by PhreeqcKernel().

◆ reinitializeRates()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::reinitializeRates ( )

private

Definition at line 138 of file PhreeqcKernel.cpp.

{
    std::vector<struct rate> rates(_reaction_rates.size());
    int rate_id = 0;
    for (auto const& reaction_rate : _reaction_rates)
    {
        rates[rate_id].name = reaction_rate.kinetic_reactant.data();
 
        // Commands' strings are freed by Phreeqc dtor.
        rates[rate_id].commands = static_cast<char*>(
            malloc(sizeof(char) * reaction_rate.commands().size() + 1));
        if (rates[rate_id].commands == nullptr)
        {
            OGS_FATAL("Could not allocate memory for rate[{:d}] commands.",
                      rate_id);
        }
        reaction_rate.commands().copy(rates[rate_id].commands,
                                      std::string::npos);
        rates[rate_id].commands[reaction_rate.commands().size()] = '\0';
 
        rates[rate_id].new_def = 1;
        rates[rate_id].linebase = nullptr;
        rates[rate_id].varbase = nullptr;
        rates[rate_id].loopbase = nullptr;
 
        ++rate_id;
    };
}

References _reaction_rates, and OGS_FATAL.

Referenced by PhreeqcKernel().

◆ reset()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::reset ( std::size_t const chemical_system_id )

private

Definition at line 295 of file PhreeqcKernel.cpp.

{
    // Clean up
    Rxn_new_solution.clear();
 
    // Solution
    {
        Rxn_solution_map[chemical_system_id] = *_aqueous_solution;
        Rxn_solution_map[chemical_system_id].Set_n_user_both(
            chemical_system_id);
        Rxn_solution_map[chemical_system_id].Set_pe(-s_eminus->la);
    }
 
    // Equilibrium reactants
    if (!Rxn_pp_assemblage_map.empty())
    {
        Rxn_new_pp_assemblage.clear();
        for (int j = 0; j < count_unknowns; j++)
        {
            if (x == nullptr || x[j]->type != PP)
                continue;
 
            auto& phase_components = Rxn_pp_assemblage_map[chemical_system_id]
                                         .Get_pp_assemblage_comps();
            phase_components[x[j]->pp_assemblage_comp_name].Set_moles(
                x[j]->moles);
        }
    }
 
    // Kinetics
    if (!Rxn_kinetics_map.empty())
    {
        Rxn_kinetics_map[chemical_system_id].Get_steps().clear();
    }
}

References _aqueous_solution.

Referenced by callPhreeqc().

◆ setAqueousSolutions()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::setAqueousSolutions ( std::vector< GlobalVector * > const & process_solutions )

Definition at line 183 of file PhreeqcKernel.cpp.

{
    assert(!process_solutions.empty());
    // Loop over chemical systems
    std::size_t const num_chemical_systems = process_solutions[0]->size();
    for (std::size_t chemical_system_id = 0;
         chemical_system_id < num_chemical_systems;
         ++chemical_system_id)
    {
        auto& aqueous_solution = Rxn_solution_map[chemical_system_id];
 
        auto initial_aqueous_solution =
            getOrCreateInitialAqueousSolution(aqueous_solution);
 
        auto& components = initial_aqueous_solution->Get_comps();
        // Loop over transport process id map to retrieve component
        // concentrations from process solutions
        for (auto const& map_pair : _process_id_to_master_map)
        {
            auto const transport_process_id = map_pair.first;
            auto const& master_species = map_pair.second;
 
            auto& transport_process_solution =
                process_solutions[transport_process_id];
 
            char const* const element_name = master_species->elt->name;
            auto const concentration =
                transport_process_solution->get(chemical_system_id);
            if (isHydrogen(element_name))
            {
                // Set pH value by hydrogen concentration.
                double const pH = -std::log10(concentration);
                aqueous_solution.Set_ph(pH);
 
                {
                    auto hydrogen = components.find("H(1)");
                    if (hydrogen != components.end())
                    {
                        hydrogen->second.Set_input_conc(pH);
                    }
                }
            }
            else
            {
                // Set component concentrations.
                components[element_name].Set_input_conc(concentration);
            }
        }
    }
}

References _process_id_to_master_map, getOrCreateInitialAqueousSolution(), ChemistryLib::PhreeqcKernelData::isHydrogen(), and ChemistryLib::pH.

Referenced by executeSpeciationCalculation().

◆ setConvergenceTolerance()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::setConvergenceTolerance ( )

inlineprivate

Definition at line 58 of file PhreeqcKernel.h.

58{ convergence_tolerance = 1e-12; }

Referenced by PhreeqcKernel().

◆ setTimeStepSize()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::setTimeStepSize ( double const dt )

private

Definition at line 247 of file PhreeqcKernel.cpp.

{
    // Loop over rxn kinetics map
    for (auto& map_pair : Rxn_kinetics_map)
    {
        auto& kinetics = map_pair.second;
        kinetics.Get_steps().push_back(dt);
    }
}

Referenced by executeSpeciationCalculation().

◆ tidyEquilibriumReactants()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::tidyEquilibriumReactants ( EquilibriumReactants const & equilibrium_reactants )

private

Definition at line 105 of file PhreeqcKernel.cpp.

{
    // extract a part of function body from int
    // Phreeqc::tidy_pp_assemblage(void)
    count_elts = 0;
    double coef = 1.0;
    for (auto const& phase_component :
         equilibrium_reactants.getPhaseComponents())
    {
        int phase_id;
        struct phase* phase_component_ptr =
            phase_bsearch(phase_component.first.c_str(), &phase_id, FALSE);
        add_elt_list(phase_component_ptr->next_elt, coef);
    }
 
    cxxNameDouble nd = elt_list_NameDouble();
    const_cast<cxxPPassemblage*>(equilibrium_reactants.castToBaseClass())
        ->Set_eltList(nd);
}

References ChemistryLib::PhreeqcKernelData::EquilibriumReactants::castToBaseClass(), and ChemistryLib::PhreeqcKernelData::EquilibriumReactants::getPhaseComponents().

Referenced by PhreeqcKernel().

◆ updateNodalProcessSolutions()

void ChemistryLib::PhreeqcKernelData::PhreeqcKernel::updateNodalProcessSolutions	(	std::vector< GlobalVector * > const &	process_solutions,
		std::size_t const	node_id )

Definition at line 331 of file PhreeqcKernel.cpp.

{
    for (auto const& map_pair : _process_id_to_master_map)
    {
        auto const transport_process_id = map_pair.first;
        auto const& master_species = map_pair.second;
 
        auto& transport_process_solution =
            process_solutions[transport_process_id];
 
        char const* const element_name = master_species->elt->name;
        if (isHydrogen(element_name))
        {
            // Update hydrogen concentration by pH value.
            auto const hydrogen_concentration = std::pow(10, s_hplus->la);
            transport_process_solution->set(node_id, hydrogen_concentration);
        }
        else
        {
            // Update solutions of component transport processes.
            auto const concentration =
                master_species->primary
                    ? master_species->total_primary / mass_water_aq_x
                    : master_species->total / mass_water_aq_x;
            transport_process_solution->set(node_id, concentration);
        }
    }
}

References _process_id_to_master_map, and ChemistryLib::PhreeqcKernelData::isHydrogen().

Referenced by callPhreeqc().

Member Data Documentation

◆ _aqueous_solution

std::unique_ptr<cxxSolution const> ChemistryLib::PhreeqcKernelData::PhreeqcKernel::_aqueous_solution

private

Definition at line 71 of file PhreeqcKernel.h.

Referenced by PhreeqcKernel(), and reset().

◆ _initial_aqueous_solution

std::unique_ptr<cxxISolution const> ChemistryLib::PhreeqcKernelData::PhreeqcKernel::_initial_aqueous_solution

private

Definition at line 70 of file PhreeqcKernel.h.

Referenced by PhreeqcKernel(), and getOrCreateInitialAqueousSolution().

◆ _process_id_to_master_map

std::map<int, struct master*> ChemistryLib::PhreeqcKernelData::PhreeqcKernel::_process_id_to_master_map

private

Definition at line 69 of file PhreeqcKernel.h.

Referenced by PhreeqcKernel(), setAqueousSolutions(), and updateNodalProcessSolutions().

◆ _reaction_rates

std::vector<ReactionRate> const ChemistryLib::PhreeqcKernelData::PhreeqcKernel::_reaction_rates

private

Definition at line 72 of file PhreeqcKernel.h.

Referenced by PhreeqcKernel(), and reinitializeRates().

The documentation for this class was generated from the following files:

Detailed Description

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Constructor & Destructor Documentation

◆ PhreeqcKernel()

Member Function Documentation

◆ callPhreeqc()

◆ configureOutputSettings()

◆ executeSpeciationCalculation()

◆ getOrCreateInitialAqueousSolution()

◆ initializePhreeqcGeneralSettings()

◆ loadDatabase()

◆ reinitializeRates()

◆ reset()

◆ setAqueousSolutions()

◆ setConvergenceTolerance()

◆ setTimeStepSize()

◆ tidyEquilibriumReactants()

◆ updateNodalProcessSolutions()

Member Data Documentation

◆ _aqueous_solution

◆ _initial_aqueous_solution

◆ _process_id_to_master_map

◆ _reaction_rates