Skip to content

Commit

Permalink
Merge pull request #2851 from renchao-lu/DistributeIndependentMeshFor…
Browse files Browse the repository at this point in the history 
…ChemicalSystem

[CL] Distribute an independent mesh for chemical system
  • Loading branch information
@endJunction
endJunction authored Mar 13, 2020
2 parents 39ab6cd + 12f5052 commit 487943a
Show file tree
Hide file tree
Showing 40 changed files with 279 additions and 183 deletions.
4 changes: 2 additions & 2 deletions Applications/ApplicationsLib/ProjectData.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1082,7 +1082,7 @@ void ProjectData::parseChemicalSystem(

_chemical_system = ChemistryLib::createChemicalSolverInterface<
ChemistryLib::ChemicalSolver::Phreeqc>(
*_mesh_vec[0], process_id_to_component_name_map, *config,
_mesh_vec, process_id_to_component_name_map, *config,
output_directory);
}
else if (boost::iequals(chemical_solver, "PhreeqcKernel"))
Expand All @@ -1093,7 +1093,7 @@ void ProjectData::parseChemicalSystem(

_chemical_system = ChemistryLib::createChemicalSolverInterface<
ChemistryLib::ChemicalSolver::PhreeqcKernel>(
*_mesh_vec[0], process_id_to_component_name_map, *config,
_mesh_vec, process_id_to_component_name_map, *config,
output_directory);
}
else
Expand Down
44 changes: 34 additions & 10 deletions ChemistryLib/CreateChemicalSolverInterface.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,11 +63,27 @@ namespace ChemistryLib
template <>
std::unique_ptr<ChemicalSolverInterface>
createChemicalSolverInterface<ChemicalSolver::Phreeqc>(
MeshLib::Mesh const& mesh,
std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes,
std::vector<std::pair<int, std::string>> const&
process_id_to_component_name_map,
BaseLib::ConfigTree const& config, std::string const& output_directory)
{
auto mesh_name =
//! \ogs_file_param{prj__chemical_system__mesh}
config.getConfigParameter<std::string>("mesh");

// Find and extract mesh from the list of meshes.
auto const& mesh = *BaseLib::findElementOrError(
std::begin(meshes), std::end(meshes),
[&mesh_name](auto const& mesh) {
assert(mesh != nullptr);
return mesh->getName() == mesh_name;
},
"Required mesh with name '" + mesh_name + "' not found.");

assert(mesh.getID() != 0);
DBUG("Found mesh '%s' with id %d.", mesh.getName().c_str(), mesh.getID());

auto path_to_database = parseDatabasePath(config);

// solution
Expand All @@ -77,9 +93,14 @@ createChemicalSolverInterface<ChemicalSolver::Phreeqc>(
process_id_to_component_name_map);

// kinetic reactants
auto chemical_system_map =
*mesh.getProperties().template getPropertyVector<std::size_t>(
"bulk_node_ids", MeshLib::MeshItemType::Node, 1);

auto kinetic_reactants = PhreeqcIOData::createKineticReactants(
//! \ogs_file_param{prj__chemical_system__kinetic_reactants}
config.getConfigSubtreeOptional("kinetic_reactants"), mesh);
config.getConfigSubtreeOptional("kinetic_reactants"), *meshes[0],
chemical_system_map);

// rates
auto reaction_rates = createReactionRates<PhreeqcIOData::ReactionRate>(
Expand All @@ -89,7 +110,8 @@ createChemicalSolverInterface<ChemicalSolver::Phreeqc>(
// equilibrium phases
auto equilibrium_phases = PhreeqcIOData::createEquilibriumPhases(
//! \ogs_file_param{prj__chemical_system__equilibrium_phases}
config.getConfigSubtreeOptional("equilibrium_phases"), mesh);
config.getConfigSubtreeOptional("equilibrium_phases"), *meshes[0],
chemical_system_map);

// surface
auto surface = PhreeqcIOData::createSurface(
Expand All @@ -113,7 +135,7 @@ createChemicalSolverInterface<ChemicalSolver::Phreeqc>(
// user punch
auto user_punch = PhreeqcIOData::createUserPunch(
//! \ogs_file_param{prj__chemical_system__user_punch}
config.getConfigSubtreeOptional("user_punch"), mesh);
config.getConfigSubtreeOptional("user_punch"), *meshes[0]);

// output
auto const& components = aqueous_solution.components;
Expand All @@ -129,21 +151,23 @@ createChemicalSolverInterface<ChemicalSolver::Phreeqc>(
num_chemical_systems, aqueous_solution);

return std::make_unique<PhreeqcIOData::PhreeqcIO>(
std::move(project_file_name), std::move(path_to_database),
std::move(aqueous_solutions), std::move(equilibrium_phases),
std::move(kinetic_reactants), std::move(reaction_rates),
std::move(surface), std::move(user_punch), std::move(output),
std::move(dump), std::move(knobs), process_id_to_component_name_map);
std::move(project_file_name), *meshes[mesh.getID()],
std::move(path_to_database), std::move(aqueous_solutions),
std::move(equilibrium_phases), std::move(kinetic_reactants),
std::move(reaction_rates), std::move(surface), std::move(user_punch),
std::move(output), std::move(dump), std::move(knobs),
process_id_to_component_name_map);
}

template <>
std::unique_ptr<ChemicalSolverInterface>
createChemicalSolverInterface<ChemicalSolver::PhreeqcKernel>(
MeshLib::Mesh const& mesh,
std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes,
std::vector<std::pair<int, std::string>> const&
process_id_to_component_name_map,
BaseLib::ConfigTree const& config, std::string const& /*output_directory*/)
{
auto mesh = *meshes[0];
auto path_to_database = parseDatabasePath(config);

// solution
Expand Down
2 changes: 1 addition & 1 deletion ChemistryLib/CreateChemicalSolverInterface.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ class ChemicalSolverInterface;

template <ChemicalSolver chemical_solver>
std::unique_ptr<ChemicalSolverInterface> createChemicalSolverInterface(
MeshLib::Mesh const& mesh,
std::vector<std::unique_ptr<MeshLib::Mesh>> const& meshes,
std::vector<std::pair<int, std::string>> const&
process_id_to_component_name_map,
BaseLib::ConfigTree const& config, std::string const& output_directory);
Expand Down
107 changes: 61 additions & 46 deletions ChemistryLib/PhreeqcIO.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@

#include "BaseLib/Algorithm.h"
#include "BaseLib/ConfigTreeUtil.h"
#include "MeshLib/Mesh.h"
#include "PhreeqcIO.h"
#include "PhreeqcIOData/AqueousSolution.h"
#include "PhreeqcIOData/Dump.h"
Expand Down Expand Up @@ -44,6 +45,7 @@ std::ostream& operator<<(std::ostream& os,
} // namespace

PhreeqcIO::PhreeqcIO(std::string const project_file_name,
MeshLib::Mesh const& mesh,
std::string&& database,
std::vector<AqueousSolution>&& aqueous_solutions,
std::vector<EquilibriumPhase>&& equilibrium_phases,
Expand All @@ -57,6 +59,7 @@ PhreeqcIO::PhreeqcIO(std::string const project_file_name,
std::vector<std::pair<int, std::string>> const&
process_id_to_component_name_map)
: _phreeqc_input_file(project_file_name + "_phreeqc.inp"),
_mesh(mesh),
_database(std::move(database)),
_aqueous_solutions(std::move(aqueous_solutions)),
_equilibrium_phases(std::move(equilibrium_phases)),
Expand Down Expand Up @@ -139,14 +142,18 @@ void PhreeqcIO::doWaterChemistryCalculation(
void PhreeqcIO::setAqueousSolutionsOrUpdateProcessSolutions(
std::vector<GlobalVector*> const& process_solutions, Status const status)
{
std::size_t const num_chemical_systems = _aqueous_solutions.size();
std::size_t const num_chemical_systems = _mesh.getNumberOfBaseNodes();

auto const chemical_system_map =
*_mesh.getProperties().template getPropertyVector<std::size_t>(
"bulk_node_ids", MeshLib::MeshItemType::Node, 1);

// Loop over chemical systems
for (std::size_t chemical_system_id = 0;
chemical_system_id < num_chemical_systems;
++chemical_system_id)
for (std::size_t local_id = 0; local_id < num_chemical_systems; ++local_id)
{
auto const global_id = chemical_system_map[local_id];
// Get chemical compostion of solution in a particular chemical system
auto& aqueous_solution = _aqueous_solutions[chemical_system_id];
auto& aqueous_solution = _aqueous_solutions[local_id];
auto& components = aqueous_solution.components;
// Loop over transport process id map to retrieve component
// concentrations from process solutions or to update process solutions
Expand Down Expand Up @@ -175,11 +182,11 @@ void PhreeqcIO::setAqueousSolutionsOrUpdateProcessSolutions(
case Status::SettingAqueousSolutions:
// Set component concentrations.
component->amount =
transport_process_solution->get(chemical_system_id);
transport_process_solution->get(global_id);
break;
case Status::UpdatingProcessSolutions:
// Update solutions of component transport processes.
transport_process_solution->set(chemical_system_id,
transport_process_solution->set(global_id,
component->amount);
break;
}
Expand All @@ -192,17 +199,16 @@ void PhreeqcIO::setAqueousSolutionsOrUpdateProcessSolutions(
case Status::SettingAqueousSolutions:
{
// Set pH value by hydrogen concentration.
aqueous_solution.pH =
-std::log10(transport_process_solution->get(
chemical_system_id));
aqueous_solution.pH = -std::log10(
transport_process_solution->get(global_id));
break;
}
case Status::UpdatingProcessSolutions:
{
// Update hydrogen concentration by pH value.
auto hydrogen_concentration =
std::pow(10, -aqueous_solution.pH);
transport_process_solution->set(chemical_system_id,
transport_process_solution->set(global_id,
hydrogen_concentration);
break;
}
Expand All @@ -226,7 +232,7 @@ void PhreeqcIO::setAqueousSolutionsPrevFromDumpFile()
OGS_FATAL("Could not open phreeqc dump file '%s'.", dump_file.c_str());
}

std::size_t const num_chemical_systems = _aqueous_solutions.size();
std::size_t const num_chemical_systems = _mesh.getNumberOfBaseNodes();
_dump->readDumpFile(in, num_chemical_systems);

if (!in)
Expand Down Expand Up @@ -280,14 +286,19 @@ std::ostream& operator<<(std::ostream& os, PhreeqcIO const& phreeqc_io)
}

std::size_t const num_chemical_systems =
phreeqc_io._aqueous_solutions.size();
for (std::size_t chemical_system_id = 0;
chemical_system_id < num_chemical_systems;
++chemical_system_id)
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)
{
auto const global_id = chemical_system_map[local_id];
auto const& aqueous_solution =
phreeqc_io._aqueous_solutions[chemical_system_id];
os << "SOLUTION " << chemical_system_id + 1 << "\n";
phreeqc_io._aqueous_solutions[local_id];
os << "SOLUTION " << global_id + 1 << "\n";
os << aqueous_solution << "\n";

auto const& dump = phreeqc_io._dump;
Expand All @@ -296,49 +307,50 @@ std::ostream& operator<<(std::ostream& os, PhreeqcIO const& phreeqc_io)
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[local_id] << "\n\n";
}
}

os << "USE solution none" << "\n";
os << "END" << "\n\n";

os << "USE solution " << chemical_system_id + 1 << "\n\n";
os << "USE solution " << global_id + 1 << "\n\n";

auto const& equilibrium_phases = phreeqc_io._equilibrium_phases;
if (!equilibrium_phases.empty())
{
os << "EQUILIBRIUM_PHASES " << chemical_system_id + 1 << "\n";
os << "EQUILIBRIUM_PHASES " << global_id + 1 << "\n";
for (auto const& equilibrium_phase : equilibrium_phases)
{
equilibrium_phase.print(os, chemical_system_id);
equilibrium_phase.print(os, global_id);
}
os << "\n";
}

auto const& kinetic_reactants = phreeqc_io._kinetic_reactants;
if (!kinetic_reactants.empty())
{
os << "KINETICS " << chemical_system_id + 1 << "\n";
os << "KINETICS " << global_id + 1 << "\n";
for (auto const& kinetic_reactant : kinetic_reactants)
{
kinetic_reactant.print(os, chemical_system_id);
kinetic_reactant.print(os, global_id);
}
os << "-steps " << phreeqc_io._dt << "\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
: 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 << "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";
}

os << "END" << "\n\n";
Expand Down Expand Up @@ -404,11 +416,16 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
int const num_skipped_lines = surface.empty() ? 1 : 2;

std::size_t const num_chemical_systems =
phreeqc_io._aqueous_solutions.size();
for (std::size_t chemical_system_id = 0;
chemical_system_id < num_chemical_systems;
++chemical_system_id)
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)
{
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)
{
Expand All @@ -421,7 +438,7 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
OGS_FATAL(
"Error when reading calculation result of Solution %u after "
"the reaction.",
chemical_system_id);
global_id);
}

std::vector<double> accepted_items;
Expand All @@ -446,25 +463,23 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
"Invalid argument. Could not convert string '%s' to "
"double for chemical system %d, column %d. Exception "
"'%s' was thrown.",
items[item_id].c_str(), chemical_system_id, item_id,
e.what());
items[item_id].c_str(), global_id, 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].c_str(), chemical_system_id, item_id,
e.what());
items[item_id].c_str(), global_id, item_id, e.what());
}
accepted_items.push_back(value);
}
}
assert(accepted_items.size() == output.accepted_items.size());

auto& aqueous_solution =
phreeqc_io._aqueous_solutions[chemical_system_id];
phreeqc_io._aqueous_solutions[local_id];
auto& components = aqueous_solution.components;
auto& equilibrium_phases = phreeqc_io._equilibrium_phases;
auto& kinetic_reactants = phreeqc_io._kinetic_reactants;
Expand Down Expand Up @@ -510,7 +525,7 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
compare_by_name,
"Could not find equilibrium phase '" + item_name +
"'.");
(*equilibrium_phase.amount)[chemical_system_id] =
(*equilibrium_phase.amount)[global_id] =
accepted_items[item_id];
break;
}
Expand All @@ -521,7 +536,7 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
kinetic_reactants.begin(), kinetic_reactants.end(),
compare_by_name,
"Could not find kinetic reactant '" + item_name + "'.");
(*kinetic_reactant.amount)[chemical_system_id] =
(*kinetic_reactant.amount)[global_id] =
accepted_items[item_id];
break;
}
Expand All @@ -535,7 +550,7 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
compare_by_name,
"Could not find secondary variable '" + item_name +
"'.");
(*secondary_variable.value)[chemical_system_id] =
(*secondary_variable.value)[global_id] =
accepted_items[item_id];
break;
}
Expand Down
Loading

0 comments on commit 487943a

Please sign in to comment.