Namespaces
namespace	anonymous_namespace{PhreeqcIO.cpp}

Classes
struct	AqueousSolution

class	BasicOutputSetups

struct	ChemicalSystem

struct	Component

struct	DensityBasedSurfaceSite

struct	Dump

struct	EquilibriumReactant

struct	ExchangeSite

struct	KineticReactant

struct	Knobs

struct	MoleBasedSurfaceSite

struct	Output

struct	OutputItem

class	PhreeqcIO

struct	ReactionRate

struct	SecondaryVariable

struct	UserPunch

Enumerations
enum class	ItemType { pH , pe , Component , EquilibriumReactant , KineticReactant , SecondaryVariable }

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

Enumerator
pH
pe
Component
EquilibriumReactant
KineticReactant
SecondaryVariable

Definition at line 59 of file 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 24 of file 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 29 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 24 of file 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 22 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 24 of file 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);
        mesh_prop_molality->resize(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 20 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 25 of file 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 20 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(), and BaseLib::ConfigTree::getValue().

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 23 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 22 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 44 of file 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 19 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 66 of file 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 483 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;
            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())
        {
            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";
            }
 
            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 19 of file 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 27 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 676 of file PhreeqcIO.cpp.

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

◆ specifyFileName()

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

extern

Definition at line 23 of file Output.cpp.

{
#ifdef USE_PETSC
    int mpi_rank;
 
    MPI_Comm_rank(MPI_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;
}

Namespaces

Classes

Enumerations

Functions

Enumeration Type Documentation

◆ ItemType

Function Documentation

◆ createAqueousSolution()

◆ createChemicalSystem()

◆ createEquilibriumReactants()

◆ createExchange()

◆ createKineticReactants()

◆ createKnobs()

◆ createOutput()

◆ createSolutionComponents()

◆ createSurface()

◆ createUserPunch()

◆ operator<<() [1/6]

◆ operator<<() [2/6]

◆ operator<<() [3/6]

◆ operator<<() [4/6]

◆ operator<<() [5/6]

◆ operator<<() [6/6]

◆ operator>>()

◆ specifyFileName()