ChemistryLib::PhreeqcIOData::PhreeqcIO Class Reference

Detailed Description

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

Per timestep:

1. 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.
2. 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.
3. 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.
4. 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).
• In stream mode, executeSpeciationCalculationParallel() handles both serial (pool size 1) and parallel cases via the PhreeqcInstancePool.
• In file mode, writeInputsToFile(), callPhreeqc(), readOutputsFromFile() implement PHREEQC I/O.

Definition at line 114 of file PhreeqcIO.h.

#include <PhreeqcIO.h>

Inheritance diagram for ChemistryLib::PhreeqcIOData::PhreeqcIO:

Collaboration diagram for ChemistryLib::PhreeqcIOData::PhreeqcIO:

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, bool use_stream_mode, int num_chemistry_threads, double concentration_warning_threshold)
	~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	writeInputHeader (std::ostream &os) const
void	writeSystemBlock (std::ostream &os, std::size_t chemical_system_id, std::size_t solution_id, std::size_t prev_solution_id, double dt) const
void	updateSystemFromOutputLine (std::string_view line, std::size_t chemical_system_id)
void	setAqueousSolutionsPrevFromDumpString (std::string_view dump_content)
void	callPhreeqc () const
void	readOutputsFromFile ()
void	executeSpeciationCalculationParallel (double const dt)
std::string	generateInputForSystem (std::size_t chemical_system_id, double const dt) const
void	parseOutputForSystem (std::string_view output_content, std::size_t chemical_system_id)
void	updateChemicalSystemFromOutput (std::vector< double > const &accepted_items, std::size_t chemical_system_id)
PhreeqcIO &	operator<< (double const dt)

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

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,
		bool	use_stream_mode,
		int	num_chemistry_threads,
		double	concentration_warning_threshold )

Definition at line 347 of file PhreeqcIO.cpp.

    : ChemicalSolverInterface(mesh, linear_solver),
      _phreeqc_input_file(specifyFileName(project_file_name, ".inp")),
      _database(std::move(database)),
      _knobs(std::move(knobs)),
      _reaction_rates(std::move(reaction_rates)),
      _chemical_system(std::move(chemical_system)),
      _user_punch(std::move(user_punch)),
      _output(std::move(output)),
      _dump(std::move(dump)),
      _concentration_warning_threshold(concentration_warning_threshold),
      num_chemistry_threads_(num_chemistry_threads),
      _use_stream_mode(use_stream_mode)
{
    INFO("Chemistry threads per MPI rank: {}.", num_chemistry_threads_);
 
    if (_use_stream_mode)
    {
        // Stream mode runs exclusively on the instance pool (even for a single
        // thread) so there is one unified execution path. The standalone
        // phreeqc_instance_id is not used here and stays -1.
        if (num_chemistry_threads_ > 1)
        {
            INFO(
                "Parallel chemistry enabled: {} threads will be used for "
                "PHREEQC calculations.",
                num_chemistry_threads_);
        }
        INFO(
            "PhreeqcIO is configured for stream-based data exchange: input "
            "and output will be exchanged via in-memory strings.");
        instance_pool_ = std::make_unique<PhreeqcInstancePool>(
            _database, std::max(1, num_chemistry_threads_));
    }
    else
    {
        // File mode: create and load the standalone PHREEQC instance and
        // enable file-based selected output. PhreeqcInstancePool::
        // createInstance() OGS_FATALs on failure, so any returned id is valid.
        phreeqc_instance_id =
            PhreeqcIOData::PhreeqcInstancePool::createInstance(_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)
        {
            SetDumpFileOn(phreeqc_instance_id, 1);
        }
    }
}

References ChemistryLib::ChemicalSolverInterface::ChemicalSolverInterface(), _chemical_system, _concentration_warning_threshold, _database, _dump, _knobs, _output, _phreeqc_input_file, _reaction_rates, _use_stream_mode, _user_punch, ChemistryLib::PhreeqcIOData::PhreeqcInstancePool::createInstance(), INFO(), instance_pool_, ChemistryLib::ChemicalSolverInterface::linear_solver, num_chemistry_threads_, OGS_FATAL, phreeqc_instance_id, and ChemistryLib::PhreeqcIOData::specifyFileName().

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

~PhreeqcIO()

ChemistryLib::PhreeqcIOData::PhreeqcIO::~PhreeqcIO

(

)

Definition at line 414 of file PhreeqcIO.cpp.

{
    // The standalone instance is only created in file mode; in stream mode
    // phreeqc_instance_id stays -1 and there is nothing to destroy.
    if (phreeqc_instance_id >= 0)
    {
        DestroyIPhreeqc(phreeqc_instance_id);
    }
}

References phreeqc_instance_id.

Member Function Documentation

callPhreeqc()

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

private

Definition at line 803 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 1016 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 512 of file PhreeqcIO.cpp.

{
    // Report the negative-concentration clamping accumulated over all chemical
    // systems in this step as a single aggregated message, then reset the
    // accumulator for the next step.
    _clamping_totals.report(_concentration_warning_threshold);
    _clamping_totals = {};
 
    if (_use_stream_mode)
    {
        // Stream mode always uses the pool (serial or parallel), giving one
        // unified execution path. The pool was created with max(1, threads).
        executeSpeciationCalculationParallel(dt);
        return;
    }
 
    // File-based data exchange.
    DBUG("Executing speciation with file-based data exchange.");
    writeInputsToFile(dt);
    callPhreeqc();
    readOutputsFromFile();
}

References _clamping_totals, _concentration_warning_threshold, _use_stream_mode, callPhreeqc(), DBUG(), executeSpeciationCalculationParallel(), readOutputsFromFile(), and writeInputsToFile().

executeSpeciationCalculationParallel()

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

private

Definition at line 1092 of file PhreeqcIO.cpp.

{
    INFO("Phreeqc: Executing parallel chemical calculation with {} threads.",
         num_chemistry_threads_);
 
    // 1. Pre-generate inputs for all chemical systems (sequential)
    // This must be sequential because it reads from shared data structures
    std::vector<std::string> inputs(_num_chemical_systems);
    for (std::size_t i = 0; i < _num_chemical_systems; ++i)
    {
        inputs[i] = generateInputForSystem(i, dt);
    }
 
    // 2. Storage for outputs
    std::vector<std::string> outputs(_num_chemical_systems);
    std::vector<std::string> dumps(_num_chemical_systems);
 
    // 3. Parallel execution
    std::vector<std::size_t> failed_systems;
    std::mutex failed_mutex;
 
#pragma omp parallel num_threads(num_chemistry_threads_)
    {
#ifdef _OPENMP
        int const thread_id = omp_get_thread_num();
#else
        int const thread_id = 0;
#endif
        int const phreeqc_id = instance_pool_->getInstanceForThread(thread_id);
 
#pragma omp for schedule(dynamic)
        for (std::ptrdiff_t i = 0;
             i < static_cast<std::ptrdiff_t>(_num_chemical_systems);
             ++i)
        {
            if (RunString(phreeqc_id, inputs[i].c_str()) != IPQ_OK)
            {
                OutputErrorString(phreeqc_id);
                std::lock_guard<std::mutex> guard(failed_mutex);
                failed_systems.push_back(i);
                continue;
            }
 
            // Retrieve output string
            const char* output_ptr = GetSelectedOutputString(phreeqc_id);
            if (output_ptr)
            {
                outputs[i] = output_ptr;
            }
 
            // Retrieve dump string if needed
            if (_dump)
            {
                const char* dump_ptr = GetDumpString(phreeqc_id);
                if (dump_ptr)
                {
                    dumps[i] = dump_ptr;
                }
            }
        }
    }
 
    std::exception_ptr local_error;
    if (!failed_systems.empty())
    {
        std::sort(failed_systems.begin(), failed_systems.end());
        std::string ids;
        for (auto const id : failed_systems)
        {
            if (!ids.empty())
            {
                ids += ", ";
            }
            ids += std::to_string(id);
        }
        local_error = std::make_exception_ptr(std::runtime_error(
            "Failed in performing speciation calculation for "
            "chemical system(s) " +
            ids + "."));
    }
    BaseLib::MPI::allRanksThrowOrNone<std::runtime_error>(local_error);
 
    // 4. Parse results (sequential, updates shared state)
    for (std::size_t i = 0; i < _num_chemical_systems; ++i)
    {
        parseOutputForSystem(outputs[i], i);
    }
 
    // 5. Handle dump data if surfaces/exchangers exist.
    // Each pool instance dumped solution 1; remap to the correct system ID.
    if (_dump && !dumps.empty())
    {
        _dump->aqueous_solutions_prev.resize(_num_chemical_systems);
        for (std::size_t i = 0; i < _num_chemical_systems; ++i)
        {
            if (!dumps[i].empty())
            {
                _dump->readDumpFromStringForSystem(dumps[i], i,
                                                   _num_chemical_systems);
            }
        }
    }
}

References _dump, _num_chemical_systems, BaseLib::MPI::allRanksThrowOrNone(), generateInputForSystem(), INFO(), instance_pool_, num_chemistry_threads_, and parseOutputForSystem().

Referenced by executeSpeciationCalculation().

generateInputForSystem()

std::string ChemistryLib::PhreeqcIOData::PhreeqcIO::generateInputForSystem ( std::size_t chemical_system_id, double const dt ) const

private

Definition at line 1035 of file PhreeqcIO.cpp.

{
    std::ostringstream os;
    os << std::scientific
       << std::setprecision(std::numeric_limits<double>::max_digits10);
 
    writeInputHeader(os);
 
    // Each pool instance processes exactly one system with solution ID 1.
    // The prev_solution_id follows the same remapping convention as in file
    // mode: num_chemical_systems + chemical_system_id + 1.
    std::size_t const prev_solution_id =
        _num_chemical_systems + chemical_system_id + 1;
    writeSystemBlock(os, chemical_system_id, 1, prev_solution_id, dt);
 
    // Each pool instance dumps its single solution (ID 1) for the next
    // timestep's dump-restore cycle.
    if (_dump)
    {
        os << "DUMP\n";
        os << "-solution 1\n";
        os << "END\n";
    }
 
    return os.str();
}

References _dump, _num_chemical_systems, writeInputHeader(), and writeSystemBlock().

Referenced by executeSpeciationCalculationParallel().

getComponentList()

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

overridevirtual

Reimplemented from ChemistryLib::ChemicalSolverInterface.

Definition at line 966 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.emplace_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 535 of file PhreeqcIO.cpp.

{
    auto const& aqueous_solution = *_chemical_system->aqueous_solution;
    auto const& components = aqueous_solution.components;
    auto const& H_plus_activity = aqueous_solution.H_plus_activity;
 
    if (component_id < static_cast<int>(components.size()))
    {
        return components[component_id].amount[chemical_system_id];
    }
 
    if (component_id != static_cast<int>(components.size()))
    {
        OGS_FATAL(
            "Invalid component_id {:d}: must be in [0, {:d}] "
            "(the last index represents H+ activity).",
            component_id, components.size());
    }
 
    // H+ activity 10^-pH (the transport-side state vector stores this
    // activity for the pH "component" — see setAqueousSolution).
    return H_plus_activity[chemical_system_id];
}

References _chemical_system, and OGS_FATAL.

initialize()

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

overridevirtual

Reimplemented from ChemistryLib::ChemicalSolverInterface.

Definition at line 424 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 436 of file PhreeqcIO.cpp.

{
    setAqueousSolution(concentrations, chemical_system_id,
                       *_chemical_system->aqueous_solution,
                       _concentration_warning_threshold);
 
    auto const& solid_phase =
        medium.phase(MaterialPropertyLib::PhaseName::Solid);
    auto const& liquid_phase =
        medium.phase(MaterialPropertyLib::PhaseName::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, _concentration_warning_threshold, MaterialPropertyLib::AqueousLiquid, MaterialPropertyLib::Medium::phase(), ChemistryLib::PhreeqcIOData::setAqueousSolution(), and MaterialPropertyLib::Solid.

operator<<()

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

inlineprivate

Definition at line 230 of file PhreeqcIO.h.

    {
        _dt = dt;
        return *this;
    }

References PhreeqcIO(), and _dt.

parseOutputForSystem()

void ChemistryLib::PhreeqcIOData::PhreeqcIO::parseOutputForSystem ( std::string_view output_content, std::size_t chemical_system_id )

private

Definition at line 1063 of file PhreeqcIO.cpp.

{
    if (output_content.empty())
    {
        OGS_FATAL("Empty output for chemical system {}.", chemical_system_id);
    }
 
    StringViewLineIterator line_iter(output_content);
    std::string_view line;
 
    line_iter.getline(line);  // skip headline
 
    int const num_skipped_lines =
        1 + (!_chemical_system->surface.empty() ? 1 : 0) +
        (!_chemical_system->exchangers.empty() ? 1 : 0);
    line_iter.skip(num_skipped_lines);
 
    if (!line_iter.getline(line))
    {
        OGS_FATAL(
            "Error when reading calculation result of Solution {} after the "
            "reaction.",
            chemical_system_id);
    }
 
    updateSystemFromOutputLine(line, chemical_system_id);
}

References _chemical_system, OGS_FATAL, and updateSystemFromOutputLine().

Referenced by executeSpeciationCalculationParallel().

readOutputsFromFile()

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

private

Definition at line 816 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 560 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.

setAqueousSolutionsPrevFromDumpString()

void ChemistryLib::PhreeqcIOData::PhreeqcIO::setAqueousSolutionsPrevFromDumpString ( std::string_view dump_content )

private

Definition at line 586 of file PhreeqcIO.cpp.

{
    if (!_dump)
    {
        return;
    }
 
    if (dump_content.empty())
    {
        // return if dump content is empty. This happens in
        // the first time step when no dump data is available.
        DBUG(
            "Dump content is empty, skipping aqueous solutions initialization "
            "from dump.");
        return;
    }
 
    _dump->readDumpFromString(dump_content, _num_chemical_systems);
}

References _dump, _num_chemical_systems, and DBUG().

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 481 of file PhreeqcIO.cpp.

{
    // Accumulate clamping over all chemical systems; reported in a single
    // aggregated message before speciation (see executeSpeciationCalculation).
    _clamping_totals += setAqueousSolution(concentrations, chemical_system_id,
                                           *_chemical_system->aqueous_solution,
                                           _concentration_warning_threshold);
 
    auto const& solid_phase =
        medium->phase(MaterialPropertyLib::PhaseName::Solid);
    auto const& liquid_phase =
        medium->phase(MaterialPropertyLib::PhaseName::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, _clamping_totals, _concentration_warning_threshold, MaterialPropertyLib::AqueousLiquid, MaterialPropertyLib::Medium::phase(), ChemistryLib::PhreeqcIOData::setAqueousSolution(), and MaterialPropertyLib::Solid.

updateChemicalSystemFromOutput()

void ChemistryLib::PhreeqcIOData::PhreeqcIO::updateChemicalSystemFromOutput ( std::vector< double > const & accepted_items, std::size_t chemical_system_id )

private

Definition at line 840 of file PhreeqcIO.cpp.

{
    auto const& output = *_output;
    auto& aqueous_solution = _chemical_system->aqueous_solution;
    auto& components = aqueous_solution->components;
    auto& equilibrium_reactants = _chemical_system->equilibrium_reactants;
    auto& kinetic_reactants = _chemical_system->kinetic_reactants;
 
    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:
            {
                aqueous_solution->H_plus_activity[chemical_system_id] =
                    std::pow(10, -accepted_items[item_id]);
                break;
            }
            case ItemType::pe:
            {
                (*aqueous_solution->pe)[chemical_system_id] =
                    accepted_items[item_id];
                break;
            }
            case ItemType::Component:
            {
                auto& component = BaseLib::findElementOrError(
                    components, compare_by_name,
                    [&]()
                    {
                        OGS_FATAL("Could not find component '{:s}'.",
                                  item_name);
                    });
                component.amount[chemical_system_id] = accepted_items[item_id];
                break;
            }
            case ItemType::EquilibriumReactant:
            {
                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:
            {
                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;
                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;
            }
        }
    }
}

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

Referenced by updateSystemFromOutputLine().

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 998 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.

updateSystemFromOutputLine()

void ChemistryLib::PhreeqcIOData::PhreeqcIO::updateSystemFromOutputLine ( std::string_view line, std::size_t chemical_system_id )

private

Parse one line of PHREEQC selected output and update the chemical system state. Used by both output-reading code paths (file-mode operator>> and stream-mode parseOutputForSystem()).

Definition at line 770 of file PhreeqcIO.cpp.

{
    auto const& output = *_output;
    auto accepted_items = parseAndFilterChemicalData(
        line, output.dropped_item_ids, chemical_system_id);
    assert(accepted_items.size() == output.accepted_items.size());
    updateChemicalSystemFromOutput(accepted_items, chemical_system_id);
}

References _output, and updateChemicalSystemFromOutput().

Referenced by operator>>, and parseOutputForSystem().

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 979 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.

writeInputHeader()

void ChemistryLib::PhreeqcIOData::PhreeqcIO::writeInputHeader ( std::ostream & os ) const

private

Write the global input header (PHASES, KNOBS, SELECTED_OUTPUT, USER_PUNCH, RATES) shared by both file and stream mode.

Definition at line 632 of file PhreeqcIO.cpp.

{
    bool const fixing_pe = _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 << _knobs << "\n";
    os << *_output << "\n";
 
    if (_user_punch)
    {
        os << *_user_punch << "\n";
    }
 
    if (!_reaction_rates.empty())
    {
        os << "RATES\n";
        os << _reaction_rates << "\n";
    }
}

References _chemical_system, _knobs, _output, _reaction_rates, and _user_punch.

Referenced by generateInputForSystem(), and operator<<.

writeInputsToFile()

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

private

Definition at line 607 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);
    *this << dt;
    out << *this;
 
    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().

writeSystemBlock()

void ChemistryLib::PhreeqcIOData::PhreeqcIO::writeSystemBlock ( std::ostream & os, std::size_t chemical_system_id, std::size_t solution_id, std::size_t prev_solution_id, double dt ) const

private

Write the per-system input block (SOLUTION, EQUILIBRIUM_PHASES, KINETICS, SURFACE, EXCHANGE).

Parameters

os	Output stream to write the PHREEQC input to.
chemical_system_id	Zero-based index of the chemical system (mesh node/cell) being written.
solution_id	PHREEQC solution number for this system (1 in stream mode, chemical_system_id+1 in file mode).
prev_solution_id	Solution number of the previous-timestep aqueous solution used for SURFACE/EXCHANGE equilibration.
dt	Time step size in seconds.

Definition at line 658 of file PhreeqcIO.cpp.

{
    bool const fixing_pe = _chemical_system->aqueous_solution->fixing_pe;
 
    os << "SOLUTION " << solution_id << "\n";
    _chemical_system->aqueous_solution->print(os, chemical_system_id);
 
    if (_dump && !_dump->aqueous_solutions_prev.empty())
    {
        os << _dump->aqueous_solutions_prev[chemical_system_id] << "\n\n";
    }
 
    os << "USE solution none\n";
    os << "END\n\n";
 
    os << "USE solution " << solution_id << "\n\n";
 
    auto const& equilibrium_reactants = _chemical_system->equilibrium_reactants;
    if (!equilibrium_reactants.empty() || fixing_pe)
    {
        os << "EQUILIBRIUM_PHASES " << solution_id << "\n";
        for (auto const& r : equilibrium_reactants)
        {
            r.print(os, chemical_system_id);
        }
        fixing_pe ? os << "Fix_pe " << -_chemical_system->aqueous_solution->pe0
                       << " O2(g)\n\n"
                  : os << "\n";
    }
 
    auto const& kinetic_reactants = _chemical_system->kinetic_reactants;
    if (!kinetic_reactants.empty())
    {
        os << "KINETICS " << solution_id << "\n";
        for (auto const& k : kinetic_reactants)
        {
            k.print(os, chemical_system_id);
        }
        os << "-steps " << dt << "\n\n";
    }
 
    auto const& surface = _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" and stored as SOLUTION_RAW within
        // aqueous_solutions_prev. Along with the PHREEQC keyword 'SURFACE',
        // distinguish between the current and previous solutions via
        // prev_solution_id.
        os << "SURFACE " << solution_id << "\n";
        std::size_t const aq_id =
            (_dump && !_dump->aqueous_solutions_prev.empty()) ? prev_solution_id
                                                              : solution_id;
        os << "-equilibrate with solution " << aq_id << "\n";
 
        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);
        }
 
        if (std::holds_alternative<MoleBasedSurfaceSite>(surface.front()))
        {
            os << "-no_edl\n";
        }
        os << "SAVE solution " << solution_id << "\n";
    }
 
    auto const& exchangers = _chemical_system->exchangers;
    if (!exchangers.empty())
    {
        os << "EXCHANGE " << solution_id << "\n";
        std::size_t const aq_id =
            (_dump && !_dump->aqueous_solutions_prev.empty()) ? prev_solution_id
                                                              : solution_id;
        os << "-equilibrate with solution " << aq_id << "\n";
        for (auto const& exchanger : exchangers)
        {
            exchanger.print(os, chemical_system_id);
        }
        os << "SAVE solution " << solution_id << "\n";
    }
 
    os << "END\n\n";
}

References _chemical_system, and _dump.

Referenced by generateInputForSystem(), and operator<<.

Friends And Related Symbol Documentation

operator<<

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

friend

Definition at line 780 of file PhreeqcIO.cpp.

{
    phreeqc_io.writeInputHeader(os);
 
    for (std::size_t chemical_system_id = 0;
         chemical_system_id < phreeqc_io._num_chemical_systems;
         ++chemical_system_id)
    {
        std::size_t const solution_id = chemical_system_id + 1;
        std::size_t const prev_solution_id =
            phreeqc_io._num_chemical_systems + chemical_system_id + 1;
        phreeqc_io.writeSystemBlock(os, chemical_system_id, solution_id,
                                    prev_solution_id, phreeqc_io._dt);
    }
 
    if (phreeqc_io._dump)
    {
        phreeqc_io._dump->print(os, phreeqc_io._num_chemical_systems);
    }
 
    return os;
}

References PhreeqcIO(), _dt, _dump, _num_chemical_systems, writeInputHeader(), and writeSystemBlock().

operator>>

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

friend

Definition at line 930 of file PhreeqcIO.cpp.

{
    // Skip the headline
    in.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
 
    std::string line;
 
    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);
 
    for (std::size_t chemical_system_id = 0;
         chemical_system_id < phreeqc_io._num_chemical_systems;
         ++chemical_system_id)
    {
        for (int i = 0; i < num_skipped_lines; ++i)
        {
            in.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
        }
 
        if (!std::getline(in, line))
        {
            OGS_FATAL(
                "Error when reading calculation result of Solution {:d} "
                "after the reaction.",
                chemical_system_id);
        }
 
        phreeqc_io.updateSystemFromOutputLine(line, chemical_system_id);
    }
 
    return in;
}

References PhreeqcIO(), _chemical_system, _num_chemical_systems, OGS_FATAL, and updateSystemFromOutputLine().

Member Data Documentation

_chemical_system

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

private

Definition at line 246 of file PhreeqcIO.h.

Referenced by PhreeqcIO(), computeSecondaryVariable(), getComponentList(), getConcentration(), initialize(), initializeChemicalSystemConcrete(), operator>>, parseOutputForSystem(), setChemicalSystemConcrete(), updateChemicalSystemFromOutput(), updatePorosityPostReaction(), updateVolumeFractionPostReaction(), writeInputHeader(), and writeSystemBlock().

_clamping_totals

ClampingStats ChemistryLib::PhreeqcIOData::PhreeqcIO::_clamping_totals

private

Negative-concentration clamping accumulated over all chemical systems within one speciation step. Folded in setChemicalSystemConcrete() and flushed (reported + reset) in executeSpeciationCalculation() so that clamping is reported once per step rather than once per cell.

Definition at line 242 of file PhreeqcIO.h.

Referenced by executeSpeciationCalculation(), and setChemicalSystemConcrete().

_concentration_warning_threshold

double ChemistryLib::PhreeqcIOData::PhreeqcIO::_concentration_warning_threshold

private

Negative (or zero) warning threshold for clamped concentrations. All negative concentrations are clamped to zero unconditionally (PHREEQC rejects negatives); this only sets the value below which a clamp is reported via an aggregated WARN. Concentrations in [threshold, 0) are treated as floating-point noise and clamped silently.

Definition at line 258 of file PhreeqcIO.h.

Referenced by PhreeqcIO(), executeSpeciationCalculation(), initializeChemicalSystemConcrete(), and setChemicalSystemConcrete().

_database

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

private

Definition at line 243 of file PhreeqcIO.h.

Referenced by PhreeqcIO().

_dt

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

private

Definition at line 251 of file PhreeqcIO.h.

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

_dump

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

private

Definition at line 249 of file PhreeqcIO.h.

Referenced by PhreeqcIO(), executeSpeciationCalculationParallel(), generateInputForSystem(), operator<<, setAqueousSolutionsPrevFromDumpFile(), setAqueousSolutionsPrevFromDumpString(), and writeSystemBlock().

_knobs

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

private

Definition at line 244 of file PhreeqcIO.h.

Referenced by PhreeqcIO(), and writeInputHeader().

_num_chemical_systems

std::size_t ChemistryLib::PhreeqcIOData::PhreeqcIO::_num_chemical_systems = 0

private

Definition at line 252 of file PhreeqcIO.h.

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

_output

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

private

Definition at line 248 of file PhreeqcIO.h.

Referenced by PhreeqcIO(), readOutputsFromFile(), updateChemicalSystemFromOutput(), updateSystemFromOutputLine(), and writeInputHeader().

_phreeqc_input_file

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

Definition at line 179 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 245 of file PhreeqcIO.h.

Referenced by PhreeqcIO(), and writeInputHeader().

_use_stream_mode

bool ChemistryLib::PhreeqcIOData::PhreeqcIO::_use_stream_mode

private

Definition at line 268 of file PhreeqcIO.h.

Referenced by PhreeqcIO(), and executeSpeciationCalculation().

_user_punch

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

private

Definition at line 247 of file PhreeqcIO.h.

Referenced by PhreeqcIO(), initialize(), updateChemicalSystemFromOutput(), and writeInputHeader().

instance_pool_

std::unique_ptr<PhreeqcInstancePool> ChemistryLib::PhreeqcIOData::PhreeqcIO::instance_pool_

private

Definition at line 250 of file PhreeqcIO.h.

Referenced by PhreeqcIO(), and executeSpeciationCalculationParallel().

num_chemistry_threads_

int ChemistryLib::PhreeqcIOData::PhreeqcIO::num_chemistry_threads_

private

Definition at line 267 of file PhreeqcIO.h.

Referenced by PhreeqcIO(), and executeSpeciationCalculationParallel().

phreeqc_instance_id

int ChemistryLib::PhreeqcIOData::PhreeqcIO::phreeqc_instance_id = -1

private

ID of the standalone file-mode PHREEQC instance (independent of the parallel pool). Created by PhreeqcInstancePool::createInstance() in the ctor only in file mode; in stream mode the pool is used instead and this stays -1. It must not be hard-coded to 0 because the IPhreeqc library hands out monotonically increasing IDs and another consumer (e.g. a unit test that constructed an IPhreeqc instance earlier) may have used 0 first.

Definition at line 266 of file PhreeqcIO.h.

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

---

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

• PhreeqcIO.h
• PhreeqcIO.cpp