diff --git a/ChemistryLib/ChemicalSolverInterface.h b/ChemistryLib/ChemicalSolverInterface.h
index 17b1464d2484130833d0fb16a41d4a365e676208..21d53e9fd140591667e8aaec3cb7517ccecc76a0 100644
--- a/ChemistryLib/ChemicalSolverInterface.h
+++ b/ChemistryLib/ChemicalSolverInterface.h
@@ -36,6 +36,6 @@ public:
virtual ~ChemicalSolverInterface() = default;
public:
- std::vector<std::vector<GlobalIndexType>> chemical_system_index_map;
+ std::vector<GlobalIndexType> chemical_system_index_map;
};
} // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp
index 4f27d58b915e7507b10cfa277061f3091f26bf1d..45af3e26944053b7e5f13436ec4f3d2e0dfd90c8 100644
--- a/ChemistryLib/PhreeqcIO.cpp
+++ b/ChemistryLib/PhreeqcIO.cpp
@@ -104,11 +104,7 @@ PhreeqcIO::PhreeqcIO(std::string const project_file_name,
void PhreeqcIO::initialize()
{
- _num_chemical_systems = std::accumulate(
- begin(chemical_system_index_map), end(chemical_system_index_map), 0,
- [](auto result, auto const& chemical_system_index) {
- return result + chemical_system_index.size();
- });
+ _num_chemical_systems = chemical_system_index_map.size();
_chemical_system->initialize(_num_chemical_systems);
@@ -248,83 +244,73 @@ std::ostream& operator<<(std::ostream& os, PhreeqcIO const& phreeqc_io)
os << reaction_rates << "\n";
}
- for (std::size_t mesh_item_id = 0;
- mesh_item_id < phreeqc_io.chemical_system_index_map.size();
- ++mesh_item_id)
+ for (std::size_t chemical_system_id = 0;
+ chemical_system_id < phreeqc_io._num_chemical_systems;
+ ++chemical_system_id)
{
- auto const& local_chemical_system =
- phreeqc_io.chemical_system_index_map[mesh_item_id];
- for (std::size_t ip = 0; ip < local_chemical_system.size(); ++ip)
- {
- auto const& chemical_system_id = local_chemical_system[ip];
-
- os << "SOLUTION " << chemical_system_id + 1 << "\n";
- phreeqc_io._chemical_system->aqueous_solution->print(
- os, 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& dump = phreeqc_io._dump;
+ if (dump)
+ {
+ auto const& aqueous_solutions_prev = dump->aqueous_solutions_prev;
+ if (!aqueous_solutions_prev.empty())
{
- auto const& aqueous_solutions_prev =
- dump->aqueous_solutions_prev;
- if (!aqueous_solutions_prev.empty())
- {
- os << aqueous_solutions_prev[chemical_system_id] << "\n\n";
- }
+ 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";
+ os << "USE solution none"
+ << "\n";
+ os << "END"
+ << "\n\n";
- auto const& equilibrium_reactants =
- phreeqc_io._chemical_system->equilibrium_reactants;
- if (!equilibrium_reactants.empty())
- {
- os << "EQUILIBRIUM_PHASES " << chemical_system_id + 1 << "\n";
- for (auto const& equilibrium_reactant : equilibrium_reactants)
- {
- equilibrium_reactant.print(os, chemical_system_id);
- }
- os << "\n";
- }
+ os << "USE solution " << chemical_system_id + 1 << "\n\n";
- auto const& kinetic_reactants =
- phreeqc_io._chemical_system->kinetic_reactants;
- if (!kinetic_reactants.empty())
+ auto const& equilibrium_reactants =
+ phreeqc_io._chemical_system->equilibrium_reactants;
+ if (!equilibrium_reactants.empty())
+ {
+ os << "EQUILIBRIUM_PHASES " << chemical_system_id + 1 << "\n";
+ for (auto const& equilibrium_reactant : equilibrium_reactants)
{
- 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";
+ equilibrium_reactant.print(os, chemical_system_id);
}
+ os << "\n";
+ }
- auto const& surface = phreeqc_io._surface;
- if (!surface.empty())
+ 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)
{
- 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";
- os << "-sites_units DENSITY"
- << "\n";
- os << surface << "\n";
- os << "SAVE solution " << chemical_system_id + 1 << "\n";
+ kinetic_reactant.print(os, chemical_system_id);
}
+ os << "-steps " << phreeqc_io._dt << "\n"
+ << "\n";
+ }
- os << "END"
- << "\n\n";
+ auto const& surface = phreeqc_io._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";
+ os << "-sites_units DENSITY"
+ << "\n";
+ os << surface << "\n";
+ os << "SAVE solution " << chemical_system_id + 1 << "\n";
}
+
+ os << "END"
+ << "\n\n";
}
auto const& dump = phreeqc_io._dump;
@@ -387,186 +373,145 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
auto& equilibrium_reactants =
phreeqc_io._chemical_system->equilibrium_reactants;
- for (auto& equilibrium_reactant : equilibrium_reactants)
- {
- std::fill(equilibrium_reactant.amount_avg->begin(),
- equilibrium_reactant.amount_avg->end(), 0.);
- }
-
auto& kinetic_reactants = phreeqc_io._chemical_system->kinetic_reactants;
- for (auto& kinetic_reactant : kinetic_reactants)
- {
- std::fill(kinetic_reactant.amount_avg->begin(),
- kinetic_reactant.amount_avg->end(), 0.);
- }
- for (std::size_t mesh_item_id = 0;
- mesh_item_id < phreeqc_io.chemical_system_index_map.size();
- ++mesh_item_id)
+ for (std::size_t chemical_system_id = 0;
+ chemical_system_id < phreeqc_io._num_chemical_systems;
+ ++chemical_system_id)
{
- auto const& local_chemical_system =
- phreeqc_io.chemical_system_index_map[mesh_item_id];
- auto const num_local_chemical_system = local_chemical_system.size();
- for (std::size_t ip = 0; ip < num_local_chemical_system; ++ip)
+ // Skip equilibrium calculation result of initial solution
+ for (int i = 0; i < num_skipped_lines; ++i)
{
- auto const& chemical_system_id = local_chemical_system[ip];
-
- // 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');
- }
+ 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);
- }
+ // 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)
+ 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())
{
- if (std::find(dropped_item_ids.begin(), dropped_item_ids.end(),
- item_id) == dropped_item_ids.end())
+ double value;
+ try
{
- 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);
+ value = std::stod(items[item_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;
- 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)
+ catch (const std::invalid_argument& e)
{
- case ItemType::pH:
- {
- // Update pH value
- aqueous_solution->pH->set(
- chemical_system_id,
- 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& component = BaseLib::findElementOrError(
- components.begin(), components.end(),
- compare_by_name,
- "Could not find component '" + item_name + "'.");
- component.amount->set(chemical_system_id,
- accepted_items[item_id]);
- break;
- }
- case ItemType::EquilibriumReactant:
- {
- // Update amounts of equilibrium reactant
- auto& equilibrium_reactant =
- BaseLib::findElementOrError(
- equilibrium_reactants.begin(),
- equilibrium_reactants.end(), compare_by_name,
- "Could not find equilibrium reactant '" +
- item_name + "'.");
- (*equilibrium_reactant.amount)[chemical_system_id] =
- accepted_items[item_id];
- (*equilibrium_reactant.amount_avg)[mesh_item_id] +=
- accepted_items[item_id];
- break;
- }
- case ItemType::KineticReactant:
- {
- // Update amounts of kinetic reactants
- auto& kinetic_reactant = BaseLib::findElementOrError(
- kinetic_reactants.begin(), kinetic_reactants.end(),
- compare_by_name,
- "Could not find kinetic reactant '" + item_name +
- "'.");
- (*kinetic_reactant.amount)[chemical_system_id] =
- accepted_items[item_id];
- (*kinetic_reactant.amount_avg)[mesh_item_id] +=
- accepted_items[item_id];
- break;
- }
- case ItemType::SecondaryVariable:
- {
- assert(user_punch);
- auto& secondary_variables =
- user_punch->secondary_variables;
- // Update values of secondary variables
- auto& secondary_variable = BaseLib::findElementOrError(
- secondary_variables.begin(),
- secondary_variables.end(), compare_by_name,
- "Could not find secondary variable '" + item_name +
- "'.");
- (*secondary_variable.value)[chemical_system_id] =
- accepted_items[item_id];
- break;
- }
+ 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());
- for (auto& equilibrium_reactant : equilibrium_reactants)
+ auto& aqueous_solution = phreeqc_io._chemical_system->aqueous_solution;
+ auto& components = aqueous_solution->components;
+ auto& user_punch = phreeqc_io._user_punch;
+ for (int item_id = 0; item_id < static_cast<int>(accepted_items.size());
+ ++item_id)
{
- (*equilibrium_reactant.amount_avg)[mesh_item_id] /=
- num_local_chemical_system;
- }
+ auto const& accepted_item = output.accepted_items[item_id];
+ auto const& item_name = accepted_item.name;
- for (auto& kinetic_reactant : kinetic_reactants)
- {
- (*kinetic_reactant.amount_avg)[mesh_item_id] /=
- num_local_chemical_system;
+ 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
+ aqueous_solution->pH->set(
+ chemical_system_id,
+ 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& component = BaseLib::findElementOrError(
+ components.begin(), components.end(), compare_by_name,
+ "Could not find component '" + item_name + "'.");
+ component.amount->set(chemical_system_id,
+ accepted_items[item_id]);
+ break;
+ }
+ case ItemType::EquilibriumReactant:
+ {
+ // Update amounts of equilibrium reactant
+ auto& equilibrium_reactant = BaseLib::findElementOrError(
+ equilibrium_reactants.begin(),
+ equilibrium_reactants.end(), compare_by_name,
+ "Could not find equilibrium reactant '" + item_name +
+ "'.");
+ (*equilibrium_reactant.amount)[chemical_system_id] =
+ accepted_items[item_id];
+ break;
+ }
+ case ItemType::KineticReactant:
+ {
+ // Update amounts of kinetic reactants
+ auto& kinetic_reactant = BaseLib::findElementOrError(
+ kinetic_reactants.begin(), kinetic_reactants.end(),
+ compare_by_name,
+ "Could not find kinetic reactant '" + item_name + "'.");
+ (*kinetic_reactant.amount)[chemical_system_id] =
+ accepted_items[item_id];
+ break;
+ }
+ case ItemType::SecondaryVariable:
+ {
+ assert(user_punch);
+ auto& secondary_variables = user_punch->secondary_variables;
+ // Update values of secondary variables
+ auto& secondary_variable = BaseLib::findElementOrError(
+ secondary_variables.begin(), secondary_variables.end(),
+ compare_by_name,
+ "Could not find secondary variable '" + item_name +
+ "'.");
+ (*secondary_variable.value)[chemical_system_id] =
+ accepted_items[item_id];
+ break;
+ }
+ }
}
}
diff --git a/ChemistryLib/PhreeqcIO.h b/ChemistryLib/PhreeqcIO.h
index ab15c9d061074e7534d34e75f9d8fd9729abcdae..6df10d8efbe8ba66468bc6fca197b0af6629fd12 100644
--- a/ChemistryLib/PhreeqcIO.h
+++ b/ChemistryLib/PhreeqcIO.h
@@ -92,7 +92,7 @@ private:
Knobs const _knobs;
double _dt = std::numeric_limits<double>::quiet_NaN();
const int phreeqc_instance_id = 0;
- int _num_chemical_systems = -1;
+ std::size_t _num_chemical_systems = -1;
};
} // namespace PhreeqcIOData
} // namespace ChemistryLib