Skip to content
Snippets Groups Projects
Commit 198b41dd authored by renchao.lu's avatar renchao.lu
Browse files

[CL] Replace with chemical system index map.

parent f95c1cfb
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
ChemistryLib/PhreeqcIO.cpp
+ 197
− 188
View file @ 198b41dd
......@@ -242,81 +242,83 @@ std::ostream& operator<<(std::ostream& os, PhreeqcIO const& phreeqc_io)
os << reaction_rates << "\n";
}
std::size_t const num_chemical_systems =
phreeqc_io._mesh.getNumberOfBaseNodes();
auto const chemical_system_map =
*phreeqc_io._mesh.getProperties()
.template getPropertyVector<std::size_t>(
"bulk_node_ids", MeshLib::MeshItemType::Node, 1);
for (std::size_t local_id = 0; local_id < num_chemical_systems; ++local_id)
for (std::size_t mesh_item_id = 0;
mesh_item_id < phreeqc_io.chemical_system_index_map.size();
++mesh_item_id)
{
auto const global_id = chemical_system_map[local_id];
os << "SOLUTION " << global_id + 1 << "\n";
phreeqc_io._chemical_system->aqueous_solution->print(os, local_id);
auto const& dump = phreeqc_io._dump;
if (dump)
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& aqueous_solutions_prev = dump->aqueous_solutions_prev;
if (!aqueous_solutions_prev.empty())
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);
auto const& dump = phreeqc_io._dump;
if (dump)
{
os << aqueous_solutions_prev[local_id] << "\n\n";
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 none" << "\n";
os << "END" << "\n\n";
os << "USE solution " << global_id + 1 << "\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())
{
os << "EQUILIBRIUM_PHASES " << global_id + 1 << "\n";
for (auto const& equilibrium_reactant : equilibrium_reactants)
auto const& equilibrium_reactants =
phreeqc_io._chemical_system->equilibrium_reactants;
if (!equilibrium_reactants.empty())
{
equilibrium_reactant.print(os, global_id);
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 << "\n";
}
auto const& kinetic_reactants =
phreeqc_io._chemical_system->kinetic_reactants;
if (!kinetic_reactants.empty())
{
os << "KINETICS " << global_id + 1 << "\n";
for (auto const& kinetic_reactant : kinetic_reactants)
auto const& kinetic_reactants =
phreeqc_io._chemical_system->kinetic_reactants;
if (!kinetic_reactants.empty())
{
kinetic_reactant.print(os, global_id);
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";
}
os << "-steps " << phreeqc_io._dt << "\n" << "\n";
}
auto const& surface = phreeqc_io._surface;
if (!surface.empty())
{
os << "SURFACE " << global_id + 1 << "\n";
std::size_t aqueous_solution_id =
dump->aqueous_solutions_prev.empty()
? global_id + 1
: num_chemical_systems + global_id + 1;
os << "-equilibrate with solution " << aqueous_solution_id << "\n";
os << "-sites_units DENSITY"
<< "\n";
os << surface << "\n";
os << "SAVE solution " << global_id + 1 << "\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";
os << "END" << "\n\n";
}
}
auto const& dump = phreeqc_io._dump;
if (dump)
{
dump->print(os, num_chemical_systems);
dump->print(os, phreeqc_io._num_chemical_systems);
}
return os;
......@@ -371,150 +373,157 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
auto const& surface = phreeqc_io._surface;
int const num_skipped_lines = surface.empty() ? 1 : 2;
std::size_t const num_chemical_systems =
phreeqc_io._mesh.getNumberOfBaseNodes();
auto const chemical_system_map =
*phreeqc_io._mesh.getProperties()
.template getPropertyVector<std::size_t>(
"bulk_node_ids", MeshLib::MeshItemType::Node, 1);
for (std::size_t local_id = 0; local_id < num_chemical_systems; ++local_id)
for (std::size_t mesh_item_id = 0;
mesh_item_id < phreeqc_io.chemical_system_index_map.size();
++mesh_item_id)
{
auto const global_id = chemical_system_map[local_id];
// Skip equilibrium calculation result of initial solution
for (int i = 0; i < num_skipped_lines; ++i)
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)
{
in.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
}
auto const& chemical_system_id = local_chemical_system[ip];
// 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.",
global_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');
}
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())
// Get calculation result of the solution after the reaction
if (!std::getline(in, line))
{
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], global_id, item_id, e.what());
}
catch (const std::out_of_range& e)
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())
{
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], global_id, item_id, e.what());
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);
}
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& equilibrium_reactants =
phreeqc_io._chemical_system->equilibrium_reactants;
auto& kinetic_reactants =
phreeqc_io._chemical_system->kinetic_reactants;
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;
assert(accepted_items.size() == output.accepted_items.size());
auto& aqueous_solution =
phreeqc_io._chemical_system->aqueous_solution;
auto& components = aqueous_solution->components;
auto& equilibrium_reactants =
phreeqc_io._chemical_system->equilibrium_reactants;
auto& kinetic_reactants =
phreeqc_io._chemical_system->kinetic_reactants;
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;
};
auto compare_by_name = [&item_name](auto const& item) {
return item.name == item_name;
};
switch (accepted_item.item_type)
{
case ItemType::pH:
switch (accepted_item.item_type)
{
// Update pH value
aqueous_solution->pH->set(local_id,
case ItemType::pH:
{
// Update pH value
aqueous_solution->pH->set(chemical_system_id,
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::pe:
{
// Update pe value
(*aqueous_solution->pe)[global_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(local_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)[global_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)[global_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)[global_id] =
accepted_items[item_id];
break;
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;
}
}
}
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment