ChemistryLib::PhreeqcIOData Namespace Reference

Namespaces
namespace	anonymous_namespace{PhreeqcIO.cpp}

Classes
struct	AqueousSolution
class	BasicOutputSetups
	Controls which built-in PHREEQC columns appear in the output file. More...
struct	ChemicalSystem
	Complete description of one local reactive system passed to PHREEQC. More...
struct	ClampingStats
struct	Component
struct	DensityBasedSurfaceSite
struct	Dump
struct	EquilibriumReactant
struct	ExchangeSite
struct	KineticReactant
struct	Knobs
struct	MoleBasedSurfaceSite
struct	Output
	Specification of which PHREEQC output columns are imported into OpenGeoSys. More...
struct	OutputItem
	One PHREEQC output column that OpenGeoSys will keep. More...
class	PhreeqcInstancePool
class	PhreeqcIO
	Drives the chemistry step in the operator-split reactive transport loop. More...
struct	ReactionRate
struct	SecondaryVariable
struct	UserPunch

Enumerations
enum class	ItemType {   pH , pe , Component , EquilibriumReactant ,   KineticReactant , SecondaryVariable }
	Category used to interpret a PHREEQC output column. More...

Functions
std::string	specifyFileName (std::string const &project_file_name, std::string const &file_extension)
std::ostream &	operator<< (std::ostream &os, PhreeqcIO const &phreeqc_io)
std::istream &	operator>> (std::istream &in, PhreeqcIO &phreeqc_io)
ClampingStats	setAqueousSolution (std::vector< double > const &concentrations, std::size_t const chemical_system_id, AqueousSolution &aqueous_solution, double const warning_threshold)
std::unique_ptr< AqueousSolution >	createAqueousSolution (BaseLib::ConfigTree const &config, MeshLib::Mesh &mesh)
std::unique_ptr< ChemicalSystem >	createChemicalSystem (BaseLib::ConfigTree const &config, MeshLib::Mesh &mesh)
std::vector< EquilibriumReactant >	createEquilibriumReactants (std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
std::vector< ExchangeSite >	createExchange (std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
std::vector< KineticReactant >	createKineticReactants (std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
Knobs	createKnobs (BaseLib::ConfigTree const &config)
std::unique_ptr< Output >	createOutput (ChemicalSystem const &chemical_system, std::unique_ptr< UserPunch > const &user_punch, bool const use_high_precision, std::string const &project_file_name)
std::vector< Component >	createSolutionComponents (BaseLib::ConfigTree const &config)
std::vector< std::variant< DensityBasedSurfaceSite, MoleBasedSurfaceSite > >	createSurface (std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh &mesh)
std::unique_ptr< UserPunch >	createUserPunch (std::optional< BaseLib::ConfigTree > const &config, MeshLib::Mesh const &mesh)
std::ostream &	operator<< (std::ostream &os, Knobs const &knobs)
std::ostream &	operator<< (std::ostream &os, BasicOutputSetups const &basic_output_setups)
std::ostream &	operator<< (std::ostream &os, Output const &output)
std::ostream &	operator<< (std::ostream &os, ReactionRate const &reaction_rate)
std::ostream &	operator<< (std::ostream &os, UserPunch const &user_punch)

Enumeration Type Documentation

ItemType

enum class ChemistryLib::PhreeqcIOData::ItemType

strong

Category used to interpret a PHREEQC output column.

Distinguishes pH / pe, total component amounts, equilibrium-phase amounts, kinetic-phase amounts, and user-defined secondary variables.

Enumerator
pH
pe
Component
EquilibriumReactant
KineticReactant
SecondaryVariable

Definition at line 71 of file ChemistryLib/PhreeqcIOData/Output.h.

{
    pH,
    pe,
    Component,
    EquilibriumReactant,
    KineticReactant,
    SecondaryVariable
};

Function Documentation

createAqueousSolution()

std::unique_ptr< AqueousSolution > ChemistryLib::PhreeqcIOData::createAqueousSolution	(	BaseLib::ConfigTree const &	config,
		MeshLib::Mesh &	mesh )

Input File Parameter

prj__chemical_system__solution__fixing_pe

Input File Parameter

prj__chemical_system__solution__temperature

Input File Parameter

prj__chemical_system__solution__pressure

Input File Parameter

prj__chemical_system__solution__pe

Definition at line 17 of file PhreeqcIOData/CreateAqueousSolution.cpp.

{
    auto const fixing_pe = config.getConfigAttribute<bool>("fixing_pe", false);

    auto const temperature = config.getConfigParameter<double>("temperature");

    auto const pressure = config.getConfigParameter<double>("pressure");

    auto const pe0 = config.getConfigParameter<double>("pe");
 
    auto pe = MeshLib::getOrCreateMeshProperty<double>(
        mesh, "pe", MeshLib::MeshItemType::IntegrationPoint, 1);
 
    auto components = createSolutionComponents(config);
 
    auto charge_balance = createChargeBalance(config);
 
    return std::make_unique<AqueousSolution>(fixing_pe, temperature, pressure,
                                             pe, pe0, std::move(components),
                                             charge_balance);
}

References ChemistryLib::createChargeBalance(), createSolutionComponents(), BaseLib::ConfigTree::getConfigAttribute(), BaseLib::ConfigTree::getConfigParameter(), MeshLib::getOrCreateMeshProperty(), MeshLib::IntegrationPoint, and pe.

Referenced by createChemicalSystem().

createChemicalSystem()

std::unique_ptr< ChemicalSystem > ChemistryLib::PhreeqcIOData::createChemicalSystem	(	BaseLib::ConfigTree const &	config,
		MeshLib::Mesh &	mesh )

Input File Parameter

prj__chemical_system__solution

Input File Parameter

prj__chemical_system__kinetic_reactants

Input File Parameter

prj__chemical_system__equilibrium_reactants

Input File Parameter

prj__chemical_system__exchangers

Input File Parameter

prj__chemical_system__surface

Definition at line 22 of file CreateChemicalSystem.cpp.

{
    // solution
    auto aqueous_solution = createAqueousSolution(
        config.getConfigSubtree("solution"), mesh);
 
    // kinetic reactants
    auto kinetic_reactants = createKineticReactants(
        config.getConfigSubtreeOptional("kinetic_reactants"), mesh);
 
    // equilibrium reactants
    auto equilibrium_reactants = createEquilibriumReactants(
        config.getConfigSubtreeOptional("equilibrium_reactants"), mesh);
 
    // exchangers
    auto exchangers = createExchange(
        config.getConfigSubtreeOptional("exchangers"), mesh);
 
    // surface
    auto surface = createSurface(
        config.getConfigSubtreeOptional("surface"), mesh);
 
    return std::make_unique<ChemicalSystem>(std::move(aqueous_solution),
                                            std::move(kinetic_reactants),
                                            std::move(equilibrium_reactants),
                                            std::move(exchangers),
                                            std::move(surface));
}

References createAqueousSolution(), createEquilibriumReactants(), createExchange(), createKineticReactants(), createSurface(), BaseLib::ConfigTree::getConfigSubtree(), and BaseLib::ConfigTree::getConfigSubtreeOptional().

createEquilibriumReactants()

std::vector< EquilibriumReactant > ChemistryLib::PhreeqcIOData::createEquilibriumReactants	(	std::optional< BaseLib::ConfigTree > const &	config,
		MeshLib::Mesh &	mesh )

Input File Parameter

prj__chemical_system__equilibrium_reactants__phase_component

Input File Parameter

prj__chemical_system__equilibrium_reactants__phase_component__name

Input File Parameter

prj__chemical_system__equilibrium_reactants__phase_component__saturation_index

Input File Parameter

prj__chemical_system__equilibrium_reactants__phase_component__reaction_irreversibility

Definition at line 17 of file PhreeqcIOData/CreateEquilibriumReactants.cpp.

{
    if (!config)
    {
        return {};
    }
 
    std::vector<EquilibriumReactant> equilibrium_reactants;
    for (
        auto const& equilibrium_reactant_config :
        config->getConfigSubtreeList("phase_component"))
    {
        auto name =
            equilibrium_reactant_config.getConfigParameter<std::string>("name");
 
        double const saturation_index =
            equilibrium_reactant_config.getConfigParameter<double>(
                "saturation_index");
 
        auto reaction_irreversibility =
            equilibrium_reactant_config.getConfigParameter<std::string>(
                "reaction_irreversibility", "");
 
        if (!reaction_irreversibility.empty() &&
            (reaction_irreversibility != "dissolve_only" &&
             reaction_irreversibility != "precipitate_only"))
        {
            OGS_FATAL(
                "{:s}: reaction direction only allows `dissolve_only` or "
                "`precipitate_only`",
                name);
        }
 
        auto molality = MeshLib::getOrCreateMeshProperty<double>(
            mesh, name, MeshLib::MeshItemType::IntegrationPoint, 1);
 
        auto molality_prev = MeshLib::getOrCreateMeshProperty<double>(
            mesh, name + "_prev", MeshLib::MeshItemType::IntegrationPoint, 1);
 
        auto volume_fraction = MeshLib::getOrCreateMeshProperty<double>(
            mesh, "phi_" + name, MeshLib::MeshItemType::IntegrationPoint, 1);
 
        auto volume_fraction_prev = MeshLib::getOrCreateMeshProperty<double>(
            mesh,
            "phi_" + name + "_prev",
            MeshLib::MeshItemType::IntegrationPoint,
            1);
 
        auto mesh_prop_molality = MeshLib::getOrCreateMeshProperty<double>(
            mesh, name + "_avg", MeshLib::MeshItemType::Cell, 1);
        mesh_prop_molality->resize(mesh.getNumberOfElements());
 
        equilibrium_reactants.emplace_back(std::move(name),
                                           molality,
                                           molality_prev,
                                           volume_fraction,
                                           volume_fraction_prev,
                                           mesh_prop_molality,
                                           saturation_index,
                                           std::move(reaction_irreversibility));
    }
 
    return equilibrium_reactants;
}

References MeshLib::Cell, MeshLib::Mesh::getNumberOfElements(), MeshLib::getOrCreateMeshProperty(), MeshLib::IntegrationPoint, and OGS_FATAL.

Referenced by createChemicalSystem().

createExchange()

std::vector< ExchangeSite > ChemistryLib::PhreeqcIOData::createExchange	(	std::optional< BaseLib::ConfigTree > const &	config,
		MeshLib::Mesh &	mesh )

Input File Parameter

prj__chemical_system__exchangers__exchange_site

Input File Parameter

prj__chemical_system__exchangers__exchange_site__ion_exchanging_species

Definition at line 15 of file CreateExchange.cpp.

{
    if (!config)
    {
        return {};
    }
 
    std::vector<ExchangeSite> exchangers;
    for (auto const& site_config :
         config->getConfigSubtreeList("exchange_site"))
    {
        auto name = site_config.getConfigParameter<std::string>(
            "ion_exchanging_species");
 
        auto const molality = MeshLib::getOrCreateMeshProperty<double>(
            mesh, name, MeshLib::MeshItemType::IntegrationPoint, 1);
 
        exchangers.emplace_back(std::move(name), molality);
    }
 
    return exchangers;
}

References MeshLib::getOrCreateMeshProperty(), and MeshLib::IntegrationPoint.

Referenced by createChemicalSystem().

createKineticReactants()

std::vector< KineticReactant > ChemistryLib::PhreeqcIOData::createKineticReactants	(	std::optional< BaseLib::ConfigTree > const &	config,
		MeshLib::Mesh &	mesh )

Input File Parameter

prj__chemical_system__kinetic_reactants__kinetic_reactant

Input File Parameter

prj__chemical_system__kinetic_reactants__kinetic_reactant__name

Input File Parameter

prj__chemical_system__kinetic_reactants__kinetic_reactant__chemical_formula

Input File Parameter

prj__chemical_system__kinetic_reactants__kinetic_reactant__parameters

Input File Parameter

prj__chemical_system__kinetic_reactants__kinetic_reactant__fix_amount

Definition at line 17 of file PhreeqcIOData/CreateKineticReactant.cpp.

{
    if (!config)
    {
        return {};
    }
 
    std::vector<KineticReactant> kinetic_reactants;
    for (
        auto const& reactant_config :
        config->getConfigSubtreeList("kinetic_reactant"))
    {
        auto name = reactant_config.getConfigParameter<std::string>("name");
 
        auto chemical_formula =
            reactant_config.getConfigParameter<std::string>("chemical_formula",
                                                            "");
 
        auto parameters =
            reactant_config.getConfigParameter<std::vector<double>>(
                "parameters", {});
 
        bool const fix_amount =
            reactant_config.getConfigParameter<bool>("fix_amount", false);
 
        auto molality = MeshLib::getOrCreateMeshProperty<double>(
            mesh, name, MeshLib::MeshItemType::IntegrationPoint, 1);
 
        auto molality_prev = MeshLib::getOrCreateMeshProperty<double>(
            mesh, name + "_prev", MeshLib::MeshItemType::IntegrationPoint, 1);
 
        auto volume_fraction = MeshLib::getOrCreateMeshProperty<double>(
            mesh, "phi_" + name, MeshLib::MeshItemType::IntegrationPoint, 1);
 
        auto volume_fraction_prev = MeshLib::getOrCreateMeshProperty<double>(
            mesh,
            "phi_" + name + "_prev",
            MeshLib::MeshItemType::IntegrationPoint,
            1);
 
        auto mesh_prop_molality = MeshLib::getOrCreateMeshProperty<double>(
            mesh, name + "_avg", MeshLib::MeshItemType::Cell, 1);
        assert(mesh_prop_molality->size() == mesh.getNumberOfElements());
 
        kinetic_reactants.emplace_back(std::move(name),
                                       std::move(chemical_formula),
                                       molality,
                                       molality_prev,
                                       volume_fraction,
                                       volume_fraction_prev,
                                       mesh_prop_molality,
                                       std::move(parameters),
                                       fix_amount);
    }
 
    return kinetic_reactants;
}

References MeshLib::Cell, MeshLib::Mesh::getNumberOfElements(), MeshLib::getOrCreateMeshProperty(), and MeshLib::IntegrationPoint.

Referenced by createChemicalSystem().

createKnobs()

Knobs ChemistryLib::PhreeqcIOData::createKnobs

(

BaseLib::ConfigTree const &

config

)

Input File Parameter

prj__chemical_system__knobs__max_iter

Input File Parameter

prj__chemical_system__knobs__relative_convergence_tolerance

Input File Parameter

prj__chemical_system__knobs__tolerance

Input File Parameter

prj__chemical_system__knobs__step_size

Input File Parameter

prj__chemical_system__knobs__scaling

Definition at line 13 of file CreateKnobs.cpp.

{
    auto const max_iterations =
        config.getConfigParameter<int>("max_iter");
 
    auto const relative_convergence_tolerance =
        config.getConfigParameter<double>("relative_convergence_tolerance");
 
    auto const tolerance =
        config.getConfigParameter<double>("tolerance");
 
    auto const step_size =
        config.getConfigParameter<int>("step_size");
 
    auto const scaling =
        config.getConfigParameter<bool>("scaling");
 
    return {max_iterations, relative_convergence_tolerance, tolerance,
            step_size, scaling};
}

References BaseLib::ConfigTree::getConfigParameter().

createOutput()

std::unique_ptr< Output > ChemistryLib::PhreeqcIOData::createOutput	(	ChemicalSystem const &	chemical_system,
		std::unique_ptr< UserPunch > const &	user_punch,
		bool const	use_high_precision,
		std::string const &	project_file_name )

Definition at line 18 of file ChemistryLib/PhreeqcIOData/CreateOutput.cpp.

{
    auto const& components = chemical_system.aqueous_solution->components;
    auto const& equilibrium_reactants = chemical_system.equilibrium_reactants;
    auto const& kinetic_reactants = chemical_system.kinetic_reactants;
    // Mark which phreeqc output items will be held.
    std::vector<OutputItem> accepted_items{{"pH", ItemType::pH},
                                           {"pe", ItemType::pe}};
 
    auto accepted_item = [](auto const& item)
    { return OutputItem(item.name, item.item_type); };
    std::transform(components.begin(), components.end(),
                   std::back_inserter(accepted_items), accepted_item);
    std::transform(equilibrium_reactants.begin(), equilibrium_reactants.end(),
                   std::back_inserter(accepted_items), accepted_item);
    for (auto const& kinetic_reactant : kinetic_reactants)
    {
        if (kinetic_reactant.fix_amount)
        {
            continue;
        }
        accepted_items.emplace_back(kinetic_reactant.name,
                                    kinetic_reactant.item_type);
    }
 
    if (user_punch)
    {
        auto const& secondary_variables = user_punch->secondary_variables;
        accepted_items.reserve(accepted_items.size() +
                               secondary_variables.size());
        std::transform(secondary_variables.begin(), secondary_variables.end(),
                       std::back_inserter(accepted_items), accepted_item);
    }
 
    // Record ids of which phreeqc output items will be dropped.
    BasicOutputSetups basic_output_setups(project_file_name,
                                          use_high_precision);
    auto const num_dropped_basic_items =
        basic_output_setups.getNumberOfDroppedItems();
    std::vector<int> dropped_item_ids(num_dropped_basic_items);
    std::iota(dropped_item_ids.begin(), dropped_item_ids.end(), 0);
 
    auto const num_dvalue_equilibrium_reactants = equilibrium_reactants.size();
    auto const num_dvalue_kinetic_reactants = kinetic_reactants.size();
    int const num_dvalue_items =
        num_dvalue_equilibrium_reactants + num_dvalue_kinetic_reactants;
 
    auto const num_basic_items =
        basic_output_setups.getNumberOfItemsInDisplay();
    auto const num_components = components.size();
    auto dvalue_item_id = num_basic_items + num_components + 1;
    for (int i = 0; i < num_dvalue_items; ++i)
    {
        dropped_item_ids.push_back(dvalue_item_id);
        dvalue_item_id += 2;
    }
 
    return std::make_unique<Output>(std::move(basic_output_setups),
                                    std::move(accepted_items),
                                    std::move(dropped_item_ids));
}

References ChemistryLib::PhreeqcIOData::ChemicalSystem::aqueous_solution, ChemistryLib::PhreeqcIOData::ChemicalSystem::equilibrium_reactants, ChemistryLib::PhreeqcIOData::BasicOutputSetups::getNumberOfDroppedItems(), ChemistryLib::PhreeqcIOData::BasicOutputSetups::getNumberOfItemsInDisplay(), ChemistryLib::PhreeqcIOData::ChemicalSystem::kinetic_reactants, pe, and pH.

createSolutionComponents()

std::vector< Component > ChemistryLib::PhreeqcIOData::createSolutionComponents

(

BaseLib::ConfigTree const &

config

)

Input File Parameter

prj__chemical_system__solution__components

Input File Parameter

prj__chemical_system__solution__components__component

Input File Parameter

prj__chemical_system__solution__components__component__chemical_formula

Definition at line 13 of file CreateSolutionComponent.cpp.

{
    std::vector<Component> components;
    auto components_config = config.getConfigSubtree("components");
 
    for (
        auto const& component_config :
        components_config.getConfigSubtreeList("component"))
    {
        auto const component_name = component_config.getValue<std::string>();
        auto const chemical_formula =
            component_config.getConfigAttribute<std::string>("chemical_formula",
                                                             "");
        components.emplace_back(component_name, chemical_formula);
    }
 
    return components;
}

References BaseLib::ConfigTree::getConfigSubtree().

Referenced by createAqueousSolution().

createSurface()

std::vector< std::variant< DensityBasedSurfaceSite, MoleBasedSurfaceSite > > ChemistryLib::PhreeqcIOData::createSurface	(	std::optional< BaseLib::ConfigTree > const &	config,
		MeshLib::Mesh &	mesh )

Input File Parameter

prj__chemical_system__surface__site_unit

Input File Parameter

prj__chemical_system__surface__site

Input File Parameter

prj__chemical_system__surface__site__name

Input File Parameter

prj__chemical_system__surface__site__site_density

Input File Parameter

prj__chemical_system__surface__site__specific_surface_area

Input File Parameter

prj__chemical_system__surface__site__mass

Input File Parameter

prj__chemical_system__surface__site

Input File Parameter

prj__chemical_system__surface__site__name

Definition at line 16 of file CreateSurface.cpp.

{
    if (!config)
    {
        return {};
    }
 
    std::vector<std::variant<DensityBasedSurfaceSite, MoleBasedSurfaceSite>>
        surface;
 
    auto const surface_site_unit =
        config->getConfigAttribute<std::string>("site_unit", "mole");
 
    if (surface_site_unit == "density")
    {
        for (auto const& site_config :
             config->getConfigSubtreeList("site"))
        {
            auto name = site_config.getConfigParameter<std::string>("name");
 
            auto const site_density =
                site_config.getConfigParameter<double>("site_density");
 
            auto const specific_surface_area =
                site_config.getConfigParameter<double>("specific_surface_area");
 
            auto const mass =
                site_config.getConfigParameter<double>("mass");
 
            surface.emplace_back(DensityBasedSurfaceSite(
                std::move(name), site_density, specific_surface_area, mass));
        }
 
        return surface;
    }
 
    if (surface_site_unit == "mole")
    {
        for (auto const& site_config :
             config->getConfigSubtreeList("site"))
        {
            auto name = site_config.getConfigParameter<std::string>("name");
 
            auto const molality = MeshLib::getOrCreateMeshProperty<double>(
                mesh, name, MeshLib::MeshItemType::IntegrationPoint, 1);
 
            surface.emplace_back(
                MoleBasedSurfaceSite(std::move(name), molality));
        }
 
        return surface;
    }
 
    OGS_FATAL("Surface site unit should be either of 'density' or 'mole'.");
}

References MeshLib::getOrCreateMeshProperty(), MeshLib::IntegrationPoint, and OGS_FATAL.

Referenced by createChemicalSystem().

createUserPunch()

std::unique_ptr< UserPunch > ChemistryLib::PhreeqcIOData::createUserPunch	(	std::optional< BaseLib::ConfigTree > const &	config,
		MeshLib::Mesh const &	mesh )

Input File Parameter

prj__chemical_system__user_punch__headline

Input File Parameter

prj__chemical_system__user_punch__statement

Definition at line 15 of file CreateUserPunch.cpp.

{
    if (!config)
    {
        return nullptr;
    }
 
    std::vector<SecondaryVariable> secondary_variables;
    for (auto const& variable_name :
         config->getConfigParameter<std::vector<std::string>>("headline"))
    {
        auto value = MeshLib::getOrCreateMeshProperty<double>(
            const_cast<MeshLib::Mesh&>(mesh),
            variable_name,
            MeshLib::MeshItemType::IntegrationPoint,
            1);
        std::fill(std::begin(*value),
                  std::end(*value),
                  std::numeric_limits<double>::quiet_NaN());
 
        secondary_variables.emplace_back(variable_name, value);
    }
 
    std::vector<std::string> statements;
    for (auto const& statement :
         config->getConfigParameterList<std::string>("statement"))
    {
        statements.emplace_back(statement);
    }
 
    return std::make_unique<UserPunch>(std::move(secondary_variables),
                                       std::move(statements));
}

References MeshLib::getOrCreateMeshProperty(), and MeshLib::IntegrationPoint.

operator<<() [1/6]

std::ostream & ChemistryLib::PhreeqcIOData::operator<<	(	std::ostream &	os,
		BasicOutputSetups const &	basic_output_setups )

Definition at line 39 of file ChemistryLib/PhreeqcIOData/Output.cpp.

{
    os << "SELECTED_OUTPUT"
       << "\n";
    os << "-file " << basic_output_setups.output_file << "\n";
    os << "-high_precision " << std::boolalpha
       << basic_output_setups.use_high_precision << "\n";
    os << "-simulation " << std::boolalpha
       << BasicOutputSetups::display_simulation_id << "\n";
    os << "-state " << std::boolalpha << BasicOutputSetups::display_state
       << "\n";
    os << "-distance " << std::boolalpha << BasicOutputSetups::display_distance
       << "\n";
    os << "-time " << std::boolalpha << BasicOutputSetups::display_current_time
       << "\n";
    os << "-step " << std::boolalpha << BasicOutputSetups::display_time_step
       << "\n";
 
    return os;
}

operator<<() [2/6]

std::ostream & ChemistryLib::PhreeqcIOData::operator<<	(	std::ostream &	os,
		Knobs const &	knobs )

Definition at line 12 of file Knobs.cpp.

{
    os << "KNOBS"
       << "\n";
    os << "-iterations " << knobs.max_iterations << "\n";
    os << "-convergence_tolerance " << knobs.relative_convergence_tolerance
       << "\n";
    os << "-tolerance " << knobs.tolerance << "\n";
    os << "-step_size " << knobs.step_size << "\n";
    os << "-diagonal_scale " << knobs.scaling << "\n";
 
    return os;
}

operator<<() [3/6]

std::ostream & ChemistryLib::PhreeqcIOData::operator<<	(	std::ostream &	os,
		Output const &	output )

Definition at line 61 of file ChemistryLib/PhreeqcIOData/Output.cpp.

{
    os << output.basic_output_setups;
 
    auto const component_items =
        output.getOutputItemsByItemType(ItemType::Component);
    os << "-totals";
    for (auto const& component_item : component_items)
    {
        os << " " << component_item.name;
    }
    os << "\n";
 
    auto const equilibrium_phase_items =
        output.getOutputItemsByItemType(ItemType::EquilibriumReactant);
    if (!equilibrium_phase_items.empty())
    {
        os << "-equilibrium_phases";
        for (auto const& equilibrium_phase_item : equilibrium_phase_items)
        {
            os << " " << equilibrium_phase_item.name;
        }
        os << "\n";
    }
 
    auto const kinetic_reactant_items =
        output.getOutputItemsByItemType(ItemType::KineticReactant);
    if (!kinetic_reactant_items.empty())
    {
        os << "-kinetic_reactants";
        for (auto const& kinetic_reactant_item : kinetic_reactant_items)
        {
            os << " " << kinetic_reactant_item.name;
        }
        os << "\n";
    }
 
    return os;
}

operator<<() [4/6]

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

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

operator<<() [5/6]

std::ostream & ChemistryLib::PhreeqcIOData::operator<<	(	std::ostream &	os,
		ReactionRate const &	reaction_rate )

Definition at line 12 of file PhreeqcIOData/ReactionRate.cpp.

{
    os << reaction_rate.kinetic_reactant << "\n";
    os << "-start"
       << "\n";
    int line_number = 1;
    for (auto const& expression_statement : reaction_rate.expression_statements)
    {
        os << line_number << " " << expression_statement << "\n";
        ++line_number;
    }
    os << "-end"
       << "\n";
 
    return os;
}

operator<<() [6/6]

std::ostream & ChemistryLib::PhreeqcIOData::operator<<	(	std::ostream &	os,
		UserPunch const &	user_punch )

Definition at line 20 of file UserPunch.cpp.

{
    os << "USER_PUNCH"
       << "\n";
    os << "-headings ";
    auto const& secondary_variables = user_punch.secondary_variables;
    for (auto const& secondary_variable : secondary_variables)
    {
        os << secondary_variable.name << " ";
    }
    os << "\n";
 
    os << "-start"
       << "\n";
    int line_number = 1;
    for (auto const& statement : user_punch.statements)
    {
        os << line_number << " " << statement << "\n";
        ++line_number;
    }
    os << "-end"
       << "\n";
 
    return os;
}

operator>>()

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

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

setAqueousSolution()

ClampingStats ChemistryLib::PhreeqcIOData::setAqueousSolution	(	std::vector< double > const &	concentrations,
		std::size_t	chemical_system_id,
		AqueousSolution &	aqueous_solution,
		double	warning_threshold )

Writes concentrations into aqueous_solution for the given chemical_system_id and returns the clamping statistics.

Negative component concentrations are clamped to zero: values in \([\text{warning\_threshold}, 0)\) are treated as floating-point noise (e.g. from MPI partitioning or non-monotone advection/diffusion schemes) and clamped silently, whereas values below warning_threshold (which is negative) are clamped and counted as "severe" in the returned ClampingStats. The caller is responsible for reporting these (e.g. as an aggregated warning across all chemical systems). The last entry of concentrations is the H+ activity \(10^{-\text{pH}}\) for the pH "component".

Definition at line 55 of file PhreeqcIOData/AqueousSolution.cpp.

{
    ClampingStats stats;
    auto& components = aqueous_solution.components;
    for (unsigned component_id = 0; component_id < components.size();
         ++component_id)
    {
        components[component_id].amount[chemical_system_id] =
            stats.clamp(concentrations[component_id], warning_threshold,
                        components[component_id].name);
    }
    // clamp() collects value-level counters; the per-cell counters are this
    // one cell's contribution to the run-level accumulator.
    stats.n_cells = stats.n_values > 0 ? 1 : 0;
    stats.n_severe_cells = stats.n_severe_values > 0 ? 1 : 0;
 
    // The transport process carries the H+ activity 10^-pH as the last entry
    // of the concentrations vector for the pH "component". Unlike component
    // concentrations it is not clamped: PHREEQC assumes 1 kg of water, and
    // water autoprotolysis always yields some H+, so a non-positive H+ activity
    // is physically impossible. Clamping it to zero would make
    // pH = -log10(activity) NaN/inf, so a non-positive value here signals a
    // transport-side bug rather than floating-point noise.
    double const h_plus_activity = concentrations.back();
    if (h_plus_activity <= 0.0)
    {
        OGS_FATAL(
            "H+ activity (10^-pH) at chemical system {:d} is {:g} <= 0, which "
            "is physically impossible (PHREEQC's 1 kg-water assumption "
            "guarantees a positive H+ activity).",
            chemical_system_id, h_plus_activity);
    }
    aqueous_solution.H_plus_activity[chemical_system_id] = h_plus_activity;
    return stats;
}

References ChemistryLib::PhreeqcIOData::ClampingStats::clamp(), ChemistryLib::PhreeqcIOData::AqueousSolution::components, ChemistryLib::PhreeqcIOData::AqueousSolution::H_plus_activity, ChemistryLib::PhreeqcIOData::ClampingStats::n_cells, ChemistryLib::PhreeqcIOData::ClampingStats::n_severe_cells, ChemistryLib::PhreeqcIOData::ClampingStats::n_severe_values, ChemistryLib::PhreeqcIOData::ClampingStats::n_values, and OGS_FATAL.

Referenced by ChemistryLib::PhreeqcIOData::PhreeqcIO::initializeChemicalSystemConcrete(), and ChemistryLib::PhreeqcIOData::PhreeqcIO::setChemicalSystemConcrete().

specifyFileName()

std::string ChemistryLib::PhreeqcIOData::specifyFileName ( std::string const & project_file_name, std::string const & file_extension )

extern

Definition at line 18 of file ChemistryLib/PhreeqcIOData/Output.cpp.

{
#ifdef USE_PETSC
    int mpi_rank;
 
    MPI_Comm_rank(BaseLib::MPI::OGS_COMM_WORLD, &mpi_rank);
    return project_file_name + "_phreeqc_pid_" + std::to_string(mpi_rank) +
           file_extension;
#endif
 
    return project_file_name + "_phreeqc" + file_extension;
}

References BaseLib::MPI::OGS_COMM_WORLD.

Referenced by ChemistryLib::PhreeqcIOData::BasicOutputSetups::BasicOutputSetups(), ChemistryLib::PhreeqcIOData::Dump::Dump(), and ChemistryLib::PhreeqcIOData::PhreeqcIO::PhreeqcIO().