Detailed Description

Drives the chemistry step in the operator-split reactive transport loop.

Per timestep:

initializeChemicalSystemConcrete(...) / setChemicalSystemConcrete(...): For each chemical_system_id, store the transported totals \(c_{T\alpha}\), the current porous-medium state (porosity, mineral fractions, reactive surface area), and the process variables (e.g. temperature \(T\), pressure \(p\), time \(t\), and timestep \(\Delta t\)) in the corresponding ChemicalSystem. Each chemical_system_id is treated as one local batch reactor.
executeSpeciationCalculation(dt): Build PHREEQC input for all chemical_system_id, run PHREEQC via IPhreeqc, and obtain the reacted state after advancing chemistry over \(\Delta t\):
- updated aqueous composition,
- reaction/source terms \(R_{\alpha}\) for each transported component,
- updated mineral amounts.
updateVolumeFractionPostReaction(...), updatePorosityPostReaction(...): Apply precipitation / dissolution back into OpenGeoSys by updating mineral / solid volume fractions and porosity for each local system. These updates feed into flow / transport / mechanics in the next global step.
computeSecondaryVariable(...): Optionally expose derived quantities (e.g. pH fields, mineral fractions, surface loading) for output.

Accessors:

getConcentration(component_id, chemical_system_id): Reacted total \(c_{T\alpha}\) of a transported component for that local system after the last chemistry step. Used as the starting composition for the next transport solve.
getComponentList(): List of transported components in the order expected by the transport process.

Internal data:

_chemical_system holds the per-chemical_system_id definition (aqueous solution, kinetic and equilibrium reactants, exchangers, surface sites).
writeInputsToFile(), callPhreeqc(), readOutputsFromFile() implement PHREEQC I/O.

Definition at line 110 of file PhreeqcIO.h.

#include <PhreeqcIO.h>

Inheritance diagram for ChemistryLib::PhreeqcIOData::PhreeqcIO:

[legend]

Collaboration diagram for ChemistryLib::PhreeqcIOData::PhreeqcIO:

[legend]

Public Member Functions
	PhreeqcIO (MeshLib::Mesh const &mesh, GlobalLinearSolver &linear_solver, std::string const &project_file_name, std::string &&database, std::unique_ptr< ChemicalSystem > &&chemical_system, std::vector< ReactionRate > &&reaction_rates, std::unique_ptr< UserPunch > &&user_punch, std::unique_ptr< Output > &&output, std::unique_ptr< Dump > &&dump, Knobs &&knobs)
	~PhreeqcIO ()
void	initialize () override
void	initializeChemicalSystemConcrete (std::vector< double > const &concentrations, GlobalIndexType const &chemical_system_id, MaterialPropertyLib::Medium const &medium, ParameterLib::SpatialPosition const &pos, double const t) override
void	setChemicalSystemConcrete (std::vector< double > const &concentrations, GlobalIndexType const &chemical_system_id, MaterialPropertyLib::Medium const *medium, MaterialPropertyLib::VariableArray const &vars, ParameterLib::SpatialPosition const &pos, double const t, double const dt) override
void	setAqueousSolutionsPrevFromDumpFile () override
void	executeSpeciationCalculation (double const dt) override
double	getConcentration (int const component_id, GlobalIndexType const chemical_system_id) const override
void	updateVolumeFractionPostReaction (GlobalIndexType const &chemical_system_id, MaterialPropertyLib::Medium const &medium, ParameterLib::SpatialPosition const &pos, double const porosity, double const t, double const dt) override
void	updatePorosityPostReaction (GlobalIndexType const &chemical_system_id, MaterialPropertyLib::Medium const &medium, double &porosity) override
void	computeSecondaryVariable (std::size_t const ele_id, std::vector< GlobalIndexType > const &chemical_system_indices) override
std::vector< std::string > const	getComponentList () const override
Public Member Functions inherited from ChemistryLib::ChemicalSolverInterface
	ChemicalSolverInterface (MeshLib::Mesh const &mesh, GlobalLinearSolver &linear_solver_)
std::vector< std::size_t > const &	activeElementIDs () const
virtual Eigen::SparseMatrix< double > const *	getStoichiometricMatrix () const
virtual double	getKineticPrefactor (std::size_t reaction_id) const
virtual	~ChemicalSolverInterface ()=default

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

Private Member Functions
void	writeInputsToFile (double const dt)
void	callPhreeqc () const
void	readOutputsFromFile ()
PhreeqcIO &	operator<< (double const dt)

Private Attributes
std::string const	_database
std::unique_ptr< ChemicalSystem >	_chemical_system
std::vector< ReactionRate > const	_reaction_rates
std::unique_ptr< UserPunch >	_user_punch
std::unique_ptr< Output > const	_output
std::unique_ptr< Dump > const	_dump
Knobs const	_knobs
double	_dt = std::numeric_limits<double>::quiet_NaN()
const int	phreeqc_instance_id = 0
std::size_t	_num_chemical_systems = -1

Friends
std::ostream &	operator<< (std::ostream &os, PhreeqcIO const &phreeqc_io)
std::istream &	operator>> (std::istream &in, PhreeqcIO &phreeqc_io)

Constructor & Destructor Documentation

◆ PhreeqcIO()

ChemistryLib::PhreeqcIOData::PhreeqcIO::PhreeqcIO	(	MeshLib::Mesh const &	mesh,
		GlobalLinearSolver &	linear_solver,
		std::string const &	project_file_name,
		std::string &&	database,
		std::unique_ptr< ChemicalSystem > &&	chemical_system,
		std::vector< ReactionRate > &&	reaction_rates,
		std::unique_ptr< UserPunch > &&	user_punch,
		std::unique_ptr< Output > &&	output,
		std::unique_ptr< Dump > &&	dump,
		Knobs &&	knobs )

Definition at line 258 of file PhreeqcIO.cpp.

    : ChemicalSolverInterface(mesh, linear_solver),
      _phreeqc_input_file(specifyFileName(project_file_name, ".inp")),
      _database(std::move(database)),
      _chemical_system(std::move(chemical_system)),
      _reaction_rates(std::move(reaction_rates)),
      _user_punch(std::move(user_punch)),
      _output(std::move(output)),
      _dump(std::move(dump)),
      _knobs(std::move(knobs))
{
    // initialize phreeqc instance
    if (CreateIPhreeqc() != phreeqc_instance_id)
    {
        OGS_FATAL(
            "Failed to initialize phreeqc instance, due to lack of memory.");
    }
 
    // load specified thermodynamic database
    if (LoadDatabase(phreeqc_instance_id, _database.c_str()) != IPQ_OK)
    {
        OGS_FATAL(
            "Failed in loading the specified thermodynamic database file: "
            "{:s}.",
            _database);
    }
 
    if (SetSelectedOutputFileOn(phreeqc_instance_id, 1) != IPQ_OK)
    {
        OGS_FATAL(
            "Failed to fly the flag for the specified file {:s} where phreeqc "
            "will write output.",
            _output->basic_output_setups.output_file);
    }
 
    if (_dump)
    {
        // Chemical composition of the aqueous solution of last time step will
        // be written into .dmp file once the second function argument is set to
        // one.
        SetDumpFileOn(phreeqc_instance_id, 1);
    }
}

References ChemistryLib::ChemicalSolverInterface::ChemicalSolverInterface(), _chemical_system, _database, _dump, _knobs, _output, _phreeqc_input_file, _reaction_rates, _user_punch, ChemistryLib::ChemicalSolverInterface::linear_solver, OGS_FATAL, phreeqc_instance_id, and ChemistryLib::PhreeqcIOData::specifyFileName().

Referenced by operator<<(), operator<<, and operator>>.

◆ ~PhreeqcIO()

ChemistryLib::PhreeqcIOData::PhreeqcIO::~PhreeqcIO ( )

Definition at line 311 of file PhreeqcIO.cpp.

{
    DestroyIPhreeqc(phreeqc_instance_id);
}

References phreeqc_instance_id.

Member Function Documentation

◆ callPhreeqc()

void ChemistryLib::PhreeqcIOData::PhreeqcIO::callPhreeqc ( ) const

private

Definition at line 644 of file PhreeqcIO.cpp.

{
    INFO("Phreeqc: Executing chemical calculation.");
    if (RunFile(phreeqc_instance_id, _phreeqc_input_file.c_str()) != IPQ_OK)
    {
        OutputErrorString(phreeqc_instance_id);
        OGS_FATAL(
            "Failed in performing speciation calculation with the generated "
            "phreeqc input file '{:s}'.",
            _phreeqc_input_file);
    }
}

References _phreeqc_input_file, INFO(), OGS_FATAL, and phreeqc_instance_id.

Referenced by executeSpeciationCalculation().

◆ computeSecondaryVariable()

void ChemistryLib::PhreeqcIOData::PhreeqcIO::computeSecondaryVariable	(	std::size_t const	ele_id,
		std::vector< GlobalIndexType > const &	chemical_system_indices )

overridevirtual

Reimplemented from ChemistryLib::ChemicalSolverInterface.

Definition at line 912 of file PhreeqcIO.cpp.

{
    for (auto const& kinetic_reactant : _chemical_system->kinetic_reactants)
    {
        (*kinetic_reactant.mesh_prop_molality)[ele_id] =
            averageReactantMolality(kinetic_reactant, chemical_system_indices);
    }
 
    for (auto const& equilibrium_reactant :
         _chemical_system->equilibrium_reactants)
    {
        (*equilibrium_reactant.mesh_prop_molality)[ele_id] =
            averageReactantMolality(equilibrium_reactant,
                                    chemical_system_indices);
    }
}

References _chemical_system.

◆ executeSpeciationCalculation()

void ChemistryLib::PhreeqcIOData::PhreeqcIO::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 396 of file PhreeqcIO.cpp.

{
    writeInputsToFile(dt);
 
    callPhreeqc();
 
    readOutputsFromFile();
}

References callPhreeqc(), readOutputsFromFile(), and writeInputsToFile().

◆ getComponentList()

std::vector< std::string > const ChemistryLib::PhreeqcIOData::PhreeqcIO::getComponentList ( ) const

overridevirtual

Reimplemented from ChemistryLib::ChemicalSolverInterface.

Definition at line 862 of file PhreeqcIO.cpp.

{
    std::vector<std::string> component_names;
    auto const& components = _chemical_system->aqueous_solution->components;
    std::transform(components.begin(), components.end(),
                   std::back_inserter(component_names),
                   [](auto const& c) { return c.name; });
 
    component_names.push_back("H");
 
    return component_names;
}

References _chemical_system.

◆ getConcentration()

double ChemistryLib::PhreeqcIOData::PhreeqcIO::getConcentration	(	int const	component_id,
		GlobalIndexType const	chemical_system_id ) const

overridevirtual

Reimplemented from ChemistryLib::ChemicalSolverInterface.

Definition at line 405 of file PhreeqcIO.cpp.

{
    auto const& aqueous_solution = *_chemical_system->aqueous_solution;
    auto& components = aqueous_solution.components;
    auto const& pH = *aqueous_solution.pH;
 
    if (component_id < static_cast<int>(components.size()))
    {
        MathLib::LinAlg::setLocalAccessibleVector(
            *components[component_id].amount);
 
        return components[component_id].amount->get(chemical_system_id);
    }
 
    // pH
    MathLib::LinAlg::setLocalAccessibleVector(*aqueous_solution.pH);
 
    return pH.get(chemical_system_id);
}

References _chemical_system, ChemistryLib::PhreeqcIOData::pH, and MathLib::LinAlg::setLocalAccessibleVector().

◆ initialize()

void ChemistryLib::PhreeqcIOData::PhreeqcIO::initialize ( )

overridevirtual

Reimplemented from ChemistryLib::ChemicalSolverInterface.

Definition at line 316 of file PhreeqcIO.cpp.

{
    _num_chemical_systems = chemical_system_index_map.size();
 
    _chemical_system->initialize(_num_chemical_systems);
 
    if (_user_punch)
    {
        _user_punch->initialize(_num_chemical_systems);
    }
}

References _chemical_system, _num_chemical_systems, _user_punch, and ChemistryLib::ChemicalSolverInterface::chemical_system_index_map.

◆ initializeChemicalSystemConcrete()

void ChemistryLib::PhreeqcIOData::PhreeqcIO::initializeChemicalSystemConcrete	(	std::vector< double > const &	concentrations,
		GlobalIndexType const &	chemical_system_id,
		MaterialPropertyLib::Medium const &	medium,
		ParameterLib::SpatialPosition const &	pos,
		double const	t )

overridevirtual

Reimplemented from ChemistryLib::ChemicalSolverInterface.

Definition at line 328 of file PhreeqcIO.cpp.

{
    setAqueousSolution(concentrations, chemical_system_id,
                       *_chemical_system->aqueous_solution);
 
    auto const& solid_phase = medium.phase("Solid");
    auto const& liquid_phase = medium.phase("AqueousLiquid");
 
    for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
    {
        initializeReactantMolality(kinetic_reactant, chemical_system_id,
                                   solid_phase, liquid_phase, medium, pos, t);
    }
 
    for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants)
    {
        initializeReactantMolality(equilibrium_reactant, chemical_system_id,
                                   solid_phase, liquid_phase, medium, pos, t);
    }
 
    for (auto& exchanger : _chemical_system->exchangers)
    {
        initializeSiteMolality(exchanger, chemical_system_id, solid_phase, pos,
                               t);
    }
 
    for (auto& surface_site : _chemical_system->surface)
    {
        if (auto const surface_site_ptr =
                std::get_if<MoleBasedSurfaceSite>(&surface_site))
        {
            initializeSiteMolality(*surface_site_ptr, chemical_system_id,
                                   solid_phase, pos, t);
        }
    }
}

References _chemical_system, and MaterialPropertyLib::Medium::phase().

◆ operator<<()

PhreeqcIO & ChemistryLib::PhreeqcIOData::PhreeqcIO::operator<< ( double const dt )

inlineprivate

Definition at line 181 of file PhreeqcIO.h.

    {
        _dt = dt;
        return *this;
    }

References PhreeqcIO(), and _dt.

◆ readOutputsFromFile()

void ChemistryLib::PhreeqcIOData::PhreeqcIO::readOutputsFromFile ( )

private

Definition at line 657 of file PhreeqcIO.cpp.

{
    auto const& basic_output_setups = _output->basic_output_setups;
    auto const& phreeqc_result_file = basic_output_setups.output_file;
    DBUG("Reading phreeqc results from file '{:s}'.", phreeqc_result_file);
    std::ifstream in(phreeqc_result_file);
 
    if (!in)
    {
        OGS_FATAL("Could not open phreeqc result file '{:s}'.",
                  phreeqc_result_file);
    }
 
    in >> *this;
 
    if (!in)
    {
        OGS_FATAL("Error when reading phreeqc result file '{:s}'",
                  phreeqc_result_file);
    }
 
    in.close();
}

References _output, DBUG(), and OGS_FATAL.

Referenced by executeSpeciationCalculation().

◆ setAqueousSolutionsPrevFromDumpFile()

void ChemistryLib::PhreeqcIOData::PhreeqcIO::setAqueousSolutionsPrevFromDumpFile ( )

overridevirtual

Reimplemented from ChemistryLib::ChemicalSolverInterface.

Definition at line 426 of file PhreeqcIO.cpp.

{
    if (!_dump)
    {
        return;
    }
 
    auto const& dump_file = _dump->dump_file;
    std::ifstream in(dump_file);
    if (!in)
    {
        // return if phreeqc dump file doesn't exist. This happens in
        // the first time step when no dump file is provided by the user.
        return;
    }
 
    _dump->readDumpFile(in, _num_chemical_systems);
 
    if (!in)
    {
        OGS_FATAL("Error when reading phreeqc dump file '{:s}'", dump_file);
    }
 
    in.close();
}

References _dump, _num_chemical_systems, and OGS_FATAL.

◆ setChemicalSystemConcrete()

void ChemistryLib::PhreeqcIOData::PhreeqcIO::setChemicalSystemConcrete	(	std::vector< double > const &	concentrations,
		GlobalIndexType const &	chemical_system_id,
		MaterialPropertyLib::Medium const *	medium,
		MaterialPropertyLib::VariableArray const &	vars,
		ParameterLib::SpatialPosition const &	pos,
		double const	t,
		double const	dt )

overridevirtual

Reimplemented from ChemistryLib::ChemicalSolverInterface.

Definition at line 370 of file PhreeqcIO.cpp.

{
    setAqueousSolution(concentrations, chemical_system_id,
                       *_chemical_system->aqueous_solution);
 
    auto const& solid_phase = medium->phase("Solid");
    auto const& liquid_phase = medium->phase("AqueousLiquid");
 
    for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
    {
        setReactantMolality(kinetic_reactant, chemical_system_id, solid_phase,
                            liquid_phase, vars, pos, t, dt);
    }
 
    for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants)
    {
        setReactantMolality(equilibrium_reactant, chemical_system_id,
                            solid_phase, liquid_phase, vars, pos, t, dt);
    }
}

References _chemical_system, and MaterialPropertyLib::Medium::phase().

◆ updatePorosityPostReaction()

void ChemistryLib::PhreeqcIOData::PhreeqcIO::updatePorosityPostReaction	(	GlobalIndexType const &	,
		MaterialPropertyLib::Medium const &	,
		double &	)

overridevirtual

Update porosity after chemical reactions.

This applies porosity changes caused by precipitation / dissolution in the local chemical system. The updated porosity can then be used to update flow parameters (e.g. permeability) in the next timestep.

Reimplemented from ChemistryLib::ChemicalSolverInterface.

Definition at line 894 of file PhreeqcIO.cpp.

{
    for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
    {
        setPorosityPostReaction(kinetic_reactant, chemical_system_id, medium,
                                porosity);
    }
 
    for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants)
    {
        setPorosityPostReaction(equilibrium_reactant, chemical_system_id,
                                medium, porosity);
    }
}

References _chemical_system.

◆ updateVolumeFractionPostReaction()

void ChemistryLib::PhreeqcIOData::PhreeqcIO::updateVolumeFractionPostReaction	(	GlobalIndexType const &	,
		MaterialPropertyLib::Medium const &	,
		ParameterLib::SpatialPosition const &	,
		double const	,
		double const	,
		double const	)

overridevirtual

Apply mineral precipitation / dissolution to the solid fraction in OGS.

Called after executeSpeciationCalculation(). For a given chemical_system_id, this function updates the solid / mineral inventory and related volume fractions based on the PHREEQC results for the last timestep.

Typical effect:

change in mineral volume fraction at this location,
updated solid composition available to mechanics / flow.

Reimplemented from ChemistryLib::ChemicalSolverInterface.

Definition at line 875 of file PhreeqcIO.cpp.

{
    for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
    {
        updateReactantVolumeFraction(kinetic_reactant, chemical_system_id,
                                     medium, pos, porosity, t, dt);
    }
 
    for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants)
    {
        updateReactantVolumeFraction(equilibrium_reactant, chemical_system_id,
                                     medium, pos, porosity, t, dt);
    }
}

References _chemical_system.

◆ writeInputsToFile()

void ChemistryLib::PhreeqcIOData::PhreeqcIO::writeInputsToFile ( double const dt )

private

Definition at line 452 of file PhreeqcIO.cpp.

{
    DBUG("Writing phreeqc inputs into file '{:s}'.", _phreeqc_input_file);
    std::ofstream out(_phreeqc_input_file, std::ofstream::out);
 
    if (!out)
    {
        OGS_FATAL("Could not open file '{:s}' for writing phreeqc inputs.",
                  _phreeqc_input_file);
    }
 
    out << std::scientific
        << std::setprecision(std::numeric_limits<double>::max_digits10);
    out << (*this << dt);
 
    if (!out)
    {
        OGS_FATAL("Failed in generating phreeqc input file '{:s}'.",
                  _phreeqc_input_file);
    }
 
    out.close();
}

References _phreeqc_input_file, DBUG(), and OGS_FATAL.

Referenced by executeSpeciationCalculation().

◆ operator<<

std::ostream & operator<<	(	std::ostream &	os,
		PhreeqcIO const &	phreeqc_io )

friend

Definition at line 476 of file PhreeqcIO.cpp.

{
    bool const fixing_pe =
        phreeqc_io._chemical_system->aqueous_solution->fixing_pe;
    if (fixing_pe)
    {
        os << "PHASES\n"
           << "Fix_pe\n"
           << "e- = e-\n"
           << "log_k 0.0\n\n";
    }
 
    os << phreeqc_io._knobs << "\n";
 
    os << *phreeqc_io._output << "\n";
 
    auto const& user_punch = phreeqc_io._user_punch;
    if (user_punch)
    {
        os << *user_punch << "\n";
    }
 
    auto const& reaction_rates = phreeqc_io._reaction_rates;
    if (!reaction_rates.empty())
    {
        os << "RATES\n";
        os << reaction_rates << "\n";
    }
 
    for (std::size_t chemical_system_id = 0;
         chemical_system_id < phreeqc_io._num_chemical_systems;
         ++chemical_system_id)
    {
        os << "SOLUTION " << chemical_system_id + 1 << "\n";
        phreeqc_io._chemical_system->aqueous_solution->print(
            os, chemical_system_id);
 
        auto const& dump = phreeqc_io._dump;
        if (dump)
        {
            auto const& aqueous_solutions_prev = dump->aqueous_solutions_prev;
            // Print previous aqueous solution if it exists.
            if (!aqueous_solutions_prev.empty())
            {
                os << aqueous_solutions_prev[chemical_system_id] << "\n\n";
            }
        }
 
        os << "USE solution none\n";
        os << "END\n\n";
 
        os << "USE solution " << chemical_system_id + 1 << "\n\n";
 
        auto const& equilibrium_reactants =
            phreeqc_io._chemical_system->equilibrium_reactants;
        if (!equilibrium_reactants.empty() || fixing_pe)
        {
            os << "EQUILIBRIUM_PHASES " << chemical_system_id + 1 << "\n";
            for (auto const& equilibrium_reactant : equilibrium_reactants)
            {
                equilibrium_reactant.print(os, chemical_system_id);
            }
            fixing_pe
                ? os << "Fix_pe "
                     << -phreeqc_io._chemical_system->aqueous_solution->pe0
                     << " O2(g)\n\n"
                : os << "\n";
        }
 
        auto const& kinetic_reactants =
            phreeqc_io._chemical_system->kinetic_reactants;
        if (!kinetic_reactants.empty())
        {
            os << "KINETICS " << chemical_system_id + 1 << "\n";
            for (auto const& kinetic_reactant : kinetic_reactants)
            {
                kinetic_reactant.print(os, chemical_system_id);
            }
            os << "-steps " << phreeqc_io._dt << "\n\n";
        }
 
        auto const& surface = phreeqc_io._chemical_system->surface;
        if (!surface.empty())
        {
            // To get the amount of surface species from the previous time step,
            // an equilibration calculation with the previous aqueous solution
            // is needed. The previous aqueous solution is saved using the
            // PHREEQC keyword "DUMP" in the dump file and then stored as
            // SOLUTION_RAW within aqueous_solutions_prev. Along with the
            // PHREEQC keyword 'SURFACE', distinguish between the current and
            // previous solutions, the previous solution number is added to the
            // total number of chemical systems (_num_chemical_systems).
            os << "SURFACE " << chemical_system_id + 1 << "\n";
            std::size_t aqueous_solution_id =
                dump->aqueous_solutions_prev.empty()
                    ? chemical_system_id + 1
                    : phreeqc_io._num_chemical_systems + chemical_system_id + 1;
            os << "-equilibrate with solution " << aqueous_solution_id << "\n";
 
            // print unit
            if (std::holds_alternative<DensityBasedSurfaceSite>(
                    surface.front()))
            {
                os << "-sites_units density\n";
            }
            else
            {
                os << "-sites_units absolute\n";
            }
 
            for (auto const& surface_site : surface)
            {
                std::visit(
                    overloaded{[&os](DensityBasedSurfaceSite const& s)
                               {
                                   os << s.name << " " << s.site_density << " "
                                      << s.specific_surface_area << " "
                                      << s.mass << "\n";
                               },
                               [&os, chemical_system_id](
                                   MoleBasedSurfaceSite const& s) {
                                   os << s.name << " "
                                      << (*s.molality)[chemical_system_id]
                                      << "\n";
                               }},
                    surface_site);
            }
 
            // overlook the effect of the buildup of charges onto the surface
            if (std::holds_alternative<MoleBasedSurfaceSite>(surface.front()))
            {
                os << "-no_edl\n";
            }
            // Here the PHREEQC "save" command is printed to save the solution
            // internally in PHREEQC. The solution composition is saved after
            // equilibration with the surface has been completed. The solution
            // will not be saved for use in the next PHREEQC run.
            os << "SAVE solution " << chemical_system_id + 1 << "\n";
        }
 
        auto const& exchangers = phreeqc_io._chemical_system->exchangers;
        if (!exchangers.empty())
        {
            os << "EXCHANGE " << chemical_system_id + 1 << "\n";
            std::size_t const aqueous_solution_id =
                dump->aqueous_solutions_prev.empty()
                    ? chemical_system_id + 1
                    : phreeqc_io._num_chemical_systems + chemical_system_id + 1;
            os << "-equilibrate with solution " << aqueous_solution_id << "\n";
            for (auto const& exchanger : exchangers)
            {
                exchanger.print(os, chemical_system_id);
            }
            os << "SAVE solution " << chemical_system_id + 1 << "\n";
        }
 
        os << "END\n\n";
    }
 
    auto const& dump = phreeqc_io._dump;
    if (dump)
    {
        dump->print(os, phreeqc_io._num_chemical_systems);
    }
 
    return os;
}

References PhreeqcIO(), _chemical_system, _dt, _dump, _knobs, _num_chemical_systems, _output, _reaction_rates, and _user_punch.

◆ operator>>

std::istream & operator>>	(	std::istream &	in,
		PhreeqcIO &	phreeqc_io )

friend

Definition at line 681 of file PhreeqcIO.cpp.

{
    // Skip the headline
    in.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
 
    std::string line;
    auto const& output = *phreeqc_io._output;
    auto const& dropped_item_ids = output.dropped_item_ids;
 
    auto const& surface = phreeqc_io._chemical_system->surface;
    auto const& exchangers = phreeqc_io._chemical_system->exchangers;
 
    int const num_skipped_lines =
        1 + (!surface.empty() ? 1 : 0) + (!exchangers.empty() ? 1 : 0);
 
    auto& equilibrium_reactants =
        phreeqc_io._chemical_system->equilibrium_reactants;
    auto& kinetic_reactants = phreeqc_io._chemical_system->kinetic_reactants;
 
    for (std::size_t chemical_system_id = 0;
         chemical_system_id < phreeqc_io._num_chemical_systems;
         ++chemical_system_id)
    {
        // Skip equilibrium calculation result of initial solution
        for (int i = 0; i < num_skipped_lines; ++i)
        {
            in.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
        }
 
        // Get calculation result of the solution after the reaction
        if (!std::getline(in, line))
        {
            OGS_FATAL(
                "Error when reading calculation result of Solution {:d} "
                "after the reaction.",
                chemical_system_id);
        }
 
        std::vector<double> accepted_items;
        std::vector<std::string> items;
        boost::trim_if(line, boost::is_any_of("\t "));
        boost::algorithm::split(items, line, boost::is_any_of("\t "),
                                boost::token_compress_on);
        for (int item_id = 0; item_id < static_cast<int>(items.size());
             ++item_id)
        {
            if (std::find(dropped_item_ids.begin(), dropped_item_ids.end(),
                          item_id) == dropped_item_ids.end())
            {
                double value;
                try
                {
                    value = std::stod(items[item_id]);
                }
                catch (const std::invalid_argument& e)
                {
                    OGS_FATAL(
                        "Invalid argument. Could not convert string '{:s}' to "
                        "double for chemical system {:d}, column {:d}. "
                        "Exception '{:s}' was thrown.",
                        items[item_id], chemical_system_id + 1, item_id,
                        e.what());
                }
                catch (const std::out_of_range& e)
                {
                    OGS_FATAL(
                        "Out of range error. Could not convert string "
                        "'{:s}' to double for chemical system {:d}, column "
                        "{:d}. Exception '{:s}' was thrown.",
                        items[item_id], chemical_system_id + 1, item_id,
                        e.what());
                }
                accepted_items.push_back(value);
            }
        }
        assert(accepted_items.size() == output.accepted_items.size());
 
        auto& aqueous_solution = phreeqc_io._chemical_system->aqueous_solution;
        auto& components = aqueous_solution->components;
        auto& user_punch = phreeqc_io._user_punch;
 
        GlobalIndexType const offset = aqueous_solution->pH->getRangeBegin();
        GlobalIndexType const global_index = offset + chemical_system_id;
 
        for (int item_id = 0; item_id < static_cast<int>(accepted_items.size());
             ++item_id)
        {
            auto const& accepted_item = output.accepted_items[item_id];
            auto const& item_name = accepted_item.name;
 
            auto compare_by_name = [&item_name](auto const& item)
            { return item.name == item_name; };
 
            switch (accepted_item.item_type)
            {
                case ItemType::pH:
                {
                    // Update pH value
                    MathLib::LinAlg::setLocalAccessibleVector(
                        *aqueous_solution->pH);
                    aqueous_solution->pH->set(
                        global_index, std::pow(10, -accepted_items[item_id]));
                    break;
                }
                case ItemType::pe:
                {
                    // Update pe value
                    (*aqueous_solution->pe)[chemical_system_id] =
                        accepted_items[item_id];
                    break;
                }
                case ItemType::Component:
                {
                    // Update component concentrations
                    auto const& component = BaseLib::findElementOrError(
                        components, compare_by_name,
                        [&]() {
                            OGS_FATAL("Could not find component '{:s}'.",
                                      item_name);
                        });
                    MathLib::LinAlg::setLocalAccessibleVector(
                        *component.amount);
                    component.amount->set(global_index,
                                          accepted_items[item_id]);
                    break;
                }
                case ItemType::EquilibriumReactant:
                {
                    // Update amounts of equilibrium reactant
                    auto const& equilibrium_reactant =
                        BaseLib::findElementOrError(
                            equilibrium_reactants, compare_by_name,
                            [&]()
                            {
                                OGS_FATAL(
                                    "Could not find equilibrium reactant "
                                    "'{:s}'",
                                    item_name);
                            });
                    (*equilibrium_reactant.molality)[chemical_system_id] =
                        accepted_items[item_id];
                    break;
                }
                case ItemType::KineticReactant:
                {
                    // Update amounts of kinetic reactants
                    auto const& kinetic_reactant = BaseLib::findElementOrError(
                        kinetic_reactants, compare_by_name,
                        [&]() {
                            OGS_FATAL("Could not find kinetic reactant '{:s}'.",
                                      item_name);
                        });
                    (*kinetic_reactant.molality)[chemical_system_id] =
                        accepted_items[item_id];
                    break;
                }
                case ItemType::SecondaryVariable:
                {
                    assert(user_punch);
                    auto const& secondary_variables =
                        user_punch->secondary_variables;
                    // Update values of secondary variables
                    auto const& secondary_variable =
                        BaseLib::findElementOrError(
                            secondary_variables, compare_by_name,
                            [&]() {
                                OGS_FATAL(
                                    "Could not find secondary variable '{:s}'.",
                                    item_name);
                            });
                    (*secondary_variable.value)[chemical_system_id] =
                        accepted_items[item_id];
                    break;
                }
            }
        }
    }
 
    return in;
}

References PhreeqcIO(), _chemical_system, _num_chemical_systems, _output, _user_punch, ChemistryLib::PhreeqcIOData::Component, ChemistryLib::PhreeqcIOData::EquilibriumReactant, BaseLib::findElementOrError(), ChemistryLib::PhreeqcIOData::KineticReactant, OGS_FATAL, ChemistryLib::PhreeqcIOData::pe, ChemistryLib::PhreeqcIOData::pH, ChemistryLib::PhreeqcIOData::SecondaryVariable, and MathLib::LinAlg::setLocalAccessibleVector().

Member Data Documentation

◆ _chemical_system

std::unique_ptr<ChemicalSystem> ChemistryLib::PhreeqcIOData::PhreeqcIO::_chemical_system

private

Definition at line 188 of file PhreeqcIO.h.

Referenced by PhreeqcIO(), computeSecondaryVariable(), getComponentList(), getConcentration(), initialize(), initializeChemicalSystemConcrete(), operator<<, operator>>, setChemicalSystemConcrete(), updatePorosityPostReaction(), and updateVolumeFractionPostReaction().

◆ _database

std::string const ChemistryLib::PhreeqcIOData::PhreeqcIO::_database

private

Definition at line 187 of file PhreeqcIO.h.

Referenced by PhreeqcIO().

◆ _dt

double ChemistryLib::PhreeqcIOData::PhreeqcIO::_dt = std::numeric_limits<double>::quiet_NaN()

private

Definition at line 194 of file PhreeqcIO.h.

Referenced by operator<<(), and operator<<.

◆ _dump

std::unique_ptr<Dump> const ChemistryLib::PhreeqcIOData::PhreeqcIO::_dump

private

Definition at line 192 of file PhreeqcIO.h.

Referenced by PhreeqcIO(), operator<<, and setAqueousSolutionsPrevFromDumpFile().

◆ _knobs

Knobs const ChemistryLib::PhreeqcIOData::PhreeqcIO::_knobs

private

Definition at line 193 of file PhreeqcIO.h.

Referenced by PhreeqcIO(), and operator<<.

◆ _num_chemical_systems

std::size_t ChemistryLib::PhreeqcIOData::PhreeqcIO::_num_chemical_systems = -1

private

Definition at line 196 of file PhreeqcIO.h.

Referenced by initialize(), operator<<, operator>>, and setAqueousSolutionsPrevFromDumpFile().

◆ _output

std::unique_ptr<Output> const ChemistryLib::PhreeqcIOData::PhreeqcIO::_output

private

Definition at line 191 of file PhreeqcIO.h.

Referenced by PhreeqcIO(), operator<<, operator>>, and readOutputsFromFile().

◆ _phreeqc_input_file

std::string const ChemistryLib::PhreeqcIOData::PhreeqcIO::_phreeqc_input_file

Definition at line 172 of file PhreeqcIO.h.

Referenced by PhreeqcIO(), callPhreeqc(), and writeInputsToFile().

◆ _reaction_rates

std::vector<ReactionRate> const ChemistryLib::PhreeqcIOData::PhreeqcIO::_reaction_rates

private

Definition at line 189 of file PhreeqcIO.h.

Referenced by PhreeqcIO(), and operator<<.

◆ _user_punch

std::unique_ptr<UserPunch> ChemistryLib::PhreeqcIOData::PhreeqcIO::_user_punch

private

Definition at line 190 of file PhreeqcIO.h.

Referenced by PhreeqcIO(), initialize(), operator<<, and operator>>.

◆ phreeqc_instance_id

const int ChemistryLib::PhreeqcIOData::PhreeqcIO::phreeqc_instance_id = 0

private

Definition at line 195 of file PhreeqcIO.h.

Referenced by PhreeqcIO(), ~PhreeqcIO(), and callPhreeqc().

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

Detailed Description

Public Member Functions

Public Attributes

Private Member Functions

Private Attributes

Friends

Constructor & Destructor Documentation

◆ PhreeqcIO()

◆ ~PhreeqcIO()

Member Function Documentation

◆ callPhreeqc()

◆ computeSecondaryVariable()

◆ executeSpeciationCalculation()

◆ getComponentList()

◆ getConcentration()

◆ initialize()

◆ initializeChemicalSystemConcrete()

◆ operator<<()

◆ readOutputsFromFile()

◆ setAqueousSolutionsPrevFromDumpFile()

◆ setChemicalSystemConcrete()

◆ updatePorosityPostReaction()

◆ updateVolumeFractionPostReaction()

◆ writeInputsToFile()

Friends And Related Symbol Documentation

◆ operator<<

◆ operator>>

Member Data Documentation

◆ _chemical_system

◆ _database

◆ _dt

◆ _dump

◆ _knobs

◆ _num_chemical_systems

◆ _output

◆ _phreeqc_input_file

◆ _reaction_rates

◆ _user_punch

◆ phreeqc_instance_id