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	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	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)
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.push_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.push_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 711 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;
}

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 1113 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);
        }
 
        auto accepted_items = parseAndFilterChemicalData(line, dropped_item_ids,
                                                         chemical_system_id);
        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;
}

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().