Update dependencies and refactor code

This commit is contained in:
Max Lübke 2024-03-28 14:23:20 +00:00
parent 0bfc95bb52
commit 2c51cd90ef
21 changed files with 309 additions and 502 deletions

@ -1 +1 @@
Subproject commit d6e7130746f9823e0bc5f711179e676e27ba3737
Subproject commit f695619875fcf43e411af67d592c05b199f41a15

@ -1 +1 @@
Subproject commit 8c0687a6cd4a10a79c7a554083a35eda11cc55f0
Subproject commit 3ffa0ef6242d467a269be83751c857cb24ae532e

View File

@ -29,7 +29,6 @@
#include <cstdint>
#include <string>
#include <vector>
#include "Base/RInsidePOET.hpp"
#include "argh.hpp" // Argument handler https://github.com/adishavit/argh

View File

@ -47,36 +47,37 @@ enum { PARSER_OK, PARSER_ERROR, PARSER_HELP };
* @brief Defining all simulation parameters
*
*/
struct RuntimeParameters {
// struct RuntimeParameters {
/** Count of processes in MPI_COMM_WORLD */
int world_size;
/** rank of proces in MPI_COMM_WORLD */
int world_rank;
/** indicates if DHT should be used */
bool dht_enabled;
/** apply logarithm to key before rounding */
bool dht_log;
/** indicates if timestep dt differs between iterations */
bool dt_differ;
/** Indicates, when a DHT snapshot should be written */
int dht_snaps;
/** <b>not implemented</b>: How a DHT is distributed over processes */
int dht_strategy;
/** Size of DHt per process in byter */
unsigned int dht_size_per_process;
/** Default significant digit for rounding */
int dht_significant_digits;
/** Default work package size */
unsigned int wp_size;
/** indicates if resulting grid should be stored after every iteration */
bool store_result;
/** indicating whether the progress bar during chemistry simulation should be
* printed or not */
bool print_progressbar;
// /** Count of processes in MPI_COMM_WORLD */
// int world_size;
// /** rank of proces in MPI_COMM_WORLD */
// int world_rank;
// /** indicates if DHT should be used */
// bool dht_enabled;
// /** apply logarithm to key before rounding */
// bool dht_log;
// /** indicates if timestep dt differs between iterations */
// bool dt_differ;
// /** Indicates, when a DHT snapshot should be written */
// int dht_snaps;
// /** <b>not implemented</b>: How a DHT is distributed over processes */
// int dht_strategy;
// /** Size of DHt per process in byter */
// unsigned int dht_size_per_process;
// /** Default significant digit for rounding */
// int dht_significant_digits;
// /** Default work package size */
// unsigned int wp_size;
// /** indicates if resulting grid should be stored after every iteration */
// bool store_result;
// /** indicating whether the progress bar during chemistry simulation should
// be
// * printed or not */
// bool print_progressbar;
bool interp_enabled;
};
// bool interp_enabled;
// };
// using GridParams = struct s_GridParams {
// std::array<uint32_t, 2> n_cells;

View File

@ -8,6 +8,14 @@ target_sources(POETLib
Init/ChemistryInit.cpp
DataStructures/Field.cpp
Transport/DiffusionModule.cpp
Chemistry/ChemistryModule.cpp
Chemistry/MasterFunctions.cpp
Chemistry/WorkerFunctions.cpp
Chemistry/SurrogateModels/DHT_Wrapper.cpp
Chemistry/SurrogateModels/DHT.c
Chemistry/SurrogateModels/HashFunctions.cpp
Chemistry/SurrogateModels/InterpolationModule.cpp
Chemistry/SurrogateModels/ProximityHashTable.cpp
)
target_include_directories(POETLib PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}")
@ -16,6 +24,7 @@ target_link_libraries(
PUBLIC RRuntime
PUBLIC IPhreeqcPOET
PUBLIC tug
PUBLIC MPI::MPI_C
)
# add_library(poetlib

View File

@ -0,0 +1,23 @@
#pragma once
#include <cstdint>
#include <vector>
namespace poet {
enum DHT_PROP_TYPES { DHT_TYPE_DEFAULT, DHT_TYPE_CHARGE, DHT_TYPE_TOTAL };
enum CHEMISTRY_OUT_SOURCE { CHEM_PQC, CHEM_DHT, CHEM_INTERP };
struct WorkPackage {
std::size_t size;
std::vector<std::vector<double>> input;
std::vector<std::vector<double>> output;
std::vector<std::uint8_t> mapping;
WorkPackage(std::size_t _size) : size(_size) {
input.resize(size);
output.resize(size);
mapping.resize(size, CHEM_PQC);
}
};
} // namespace poet

View File

@ -1,20 +1,16 @@
#include "ChemistryModule.hpp"
#include "IPhreeqcPOET.hpp"
#include "SurrogateModels/DHT_Wrapper.hpp"
#include "SurrogateModels/Interpolation.hpp"
#include <PhreeqcRM.h>
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <map>
#include <memory>
#include <mpi.h>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
constexpr uint32_t MB_FACTOR = 1E6;
@ -158,25 +154,26 @@ inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
return results;
}
poet::ChemistryModule::ChemistryModule(uint32_t nxyz, uint32_t wp_size,
std::uint32_t maxiter,
const ChemistryParams &chem_param,
MPI_Comm communicator)
: PhreeqcRM(nxyz, 1), group_comm(communicator), wp_size(wp_size),
params(chem_param) {
poet::ChemistryModule::ChemistryModule(
uint32_t wp_size_, const InitialList::ChemistryInit &chem_params,
MPI_Comm communicator)
: params(chem_params), wp_size(wp_size_), group_comm(communicator) {
MPI_Comm_rank(communicator, &comm_rank);
MPI_Comm_size(communicator, &comm_size);
MPI_Comm_size(communicator, &this->comm_size);
MPI_Comm_rank(communicator, &this->comm_rank);
this->is_sequential = comm_size == 1;
this->is_master = comm_rank == 0;
this->is_sequential = (this->comm_size == 1);
this->is_master = (this->comm_rank == 0);
this->n_cells = chem_params.total_grid_cells;
if (!is_sequential && is_master) {
MPI_Bcast(&wp_size, 1, MPI_UINT32_T, 0, this->group_comm);
MPI_Bcast(&maxiter, 1, MPI_UINT32_T, 0, this->group_comm);
if (!is_master) {
for (std::size_t i = 0; i < chem_params.pqc_ids.size(); i++) {
this->phreeqc_instances[chem_params.pqc_ids[i]] =
std::make_unique<IPhreeqcPOET>(chem_params.database,
chem_params.pqc_scripts[i],
chem_params.pqc_sol_order, wp_size_);
}
}
this->file_pad = std::ceil(std::log10(maxiter + 1));
}
poet::ChemistryModule::~ChemistryModule() {
@ -185,206 +182,45 @@ poet::ChemistryModule::~ChemistryModule() {
}
}
poet::ChemistryModule
poet::ChemistryModule::createWorker(MPI_Comm communicator,
const ChemistryParams &chem_param) {
uint32_t wp_size;
MPI_Bcast(&wp_size, 1, MPI_UINT32_T, 0, communicator);
std::uint32_t maxiter;
MPI_Bcast(&maxiter, 1, MPI_UINT32_T, 0, communicator);
return ChemistryModule(wp_size, wp_size, maxiter, chem_param, communicator);
}
void poet::ChemistryModule::RunInitFile(const std::string &input_script_path) {
if (this->is_master) {
int f_type = CHEM_INIT;
PropagateFunctionType(f_type);
int count = input_script_path.size();
ChemBCast(&count, 1, MPI_INT);
ChemBCast(const_cast<char *>(input_script_path.data()), count, MPI_CHAR);
}
this->RunFile(true, true, false, input_script_path);
this->RunString(true, false, false, "DELETE; -all; PRINT; -warnings 0;");
this->FindComponents();
PhreeqcRM::initializePOET(this->speciesPerModule, this->prop_names);
this->prop_count = prop_names.size();
char exchange = (speciesPerModule[1] == 0 ? -1 : 1);
char kinetics = (speciesPerModule[2] == 0 ? -1 : 1);
char equilibrium = (speciesPerModule[3] == 0 ? -1 : 1);
char surface = (speciesPerModule[4] == 0 ? -1 : 1);
std::vector<int> ic1;
ic1.resize(this->nxyz * 7, -1);
// TODO: hardcoded reaction modules
for (int i = 0; i < nxyz; i++) {
ic1[i] = 1; // Solution 1
ic1[nxyz + i] = equilibrium; // Equilibrium 1
ic1[2 * nxyz + i] = exchange; // Exchange none
ic1[3 * nxyz + i] = surface; // Surface none
ic1[4 * nxyz + i] = -1; // Gas phase none
ic1[5 * nxyz + i] = -1; // Solid solutions none
ic1[6 * nxyz + i] = kinetics; // Kinetics 1
}
this->InitialPhreeqc2Module(ic1);
}
void poet::ChemistryModule::initializeField(const Field &trans_field) {
if (is_master) {
int f_type = CHEM_INIT_SPECIES;
PropagateFunctionType(f_type);
}
std::vector<std::string> essentials_backup{
prop_names.begin() + speciesPerModule[0], prop_names.end()};
std::vector<std::string> new_solution_names{
this->prop_names.begin(), this->prop_names.begin() + speciesPerModule[0]};
if (is_master) {
for (auto &prop : trans_field.GetProps()) {
if (std::find(new_solution_names.begin(), new_solution_names.end(),
prop) == new_solution_names.end()) {
int size = prop.size();
ChemBCast(&size, 1, MPI_INT);
ChemBCast(prop.data(), prop.size(), MPI_CHAR);
new_solution_names.push_back(prop);
}
}
int end = 0;
ChemBCast(&end, 1, MPI_INT);
} else {
constexpr int MAXSIZE = 128;
MPI_Status status;
int recv_size;
char recv_buffer[MAXSIZE];
while (1) {
ChemBCast(&recv_size, 1, MPI_INT);
if (recv_size == 0) {
break;
}
ChemBCast(recv_buffer, recv_size, MPI_CHAR);
recv_buffer[recv_size] = '\0';
new_solution_names.push_back(std::string(recv_buffer));
}
}
// now sort the new values
std::sort(new_solution_names.begin() + 3, new_solution_names.end());
this->SetPOETSolutionNames(new_solution_names);
this->speciesPerModule[0] = new_solution_names.size();
// and append other processes than solutions
std::vector<std::string> new_prop_names = new_solution_names;
new_prop_names.insert(new_prop_names.end(), essentials_backup.begin(),
essentials_backup.end());
std::vector<std::string> old_prop_names{this->prop_names};
this->prop_names = std::move(new_prop_names);
this->prop_count = prop_names.size();
if (is_master) {
this->n_cells = trans_field.GetRequestedVecSize();
std::vector<std::vector<double>> phreeqc_dump(this->nxyz);
this->getDumpedField(phreeqc_dump);
if (is_sequential) {
std::vector<double> init_vec;
for (std::size_t i = 0; i < n_cells; i++) {
init_vec.insert(init_vec.end(), phreeqc_dump[i].begin(),
phreeqc_dump[i].end());
}
const auto tmp_buffer{init_vec};
this->unshuffleField(tmp_buffer, n_cells, prop_count, 1, init_vec);
this->chem_field = Field(n_cells, init_vec, prop_names);
return;
}
std::vector<double> &phreeqc_init = phreeqc_dump[0];
std::vector<std::vector<double>> initial_values;
for (const auto &vec : trans_field.As2DVector()) {
initial_values.push_back(vec);
}
this->base_totals = {initial_values.at(0).at(0),
initial_values.at(1).at(0)};
ChemBCast(base_totals.data(), 2, MPI_DOUBLE);
for (int i = speciesPerModule[0]; i < phreeqc_init.size(); i++) {
std::vector<double> init(n_cells, phreeqc_init[i]);
initial_values.push_back(std::move(init));
}
this->chem_field = Field(n_cells, initial_values, prop_names);
} else {
ChemBCast(base_totals.data(), 2, MPI_DOUBLE);
}
if (this->params.use_dht || this->params.use_interp) {
initializeDHT(this->params.dht_size, this->params.dht_signifs);
setDHTSnapshots(this->params.dht_snaps, this->params.dht_outdir);
setDHTReadFile(this->params.dht_file);
this->dht_enabled = this->params.use_dht;
if (this->params.use_interp) {
initializeInterp(this->params.pht_max_entries, this->params.pht_size,
this->params.interp_min_entries,
this->params.pht_signifs);
this->interp_enabled = this->params.use_interp;
}
}
}
void poet::ChemistryModule::initializeDHT(
uint32_t size_mb, const NamedVector<std::uint32_t> &key_species) {
constexpr uint32_t MB_FACTOR = 1E6;
// constexpr uint32_t MB_FACTOR = 1E6;
this->dht_enabled = true;
// this->dht_enabled = true;
MPI_Comm dht_comm;
// MPI_Comm dht_comm;
if (this->is_master) {
MPI_Comm_split(this->group_comm, MPI_UNDEFINED, this->comm_rank, &dht_comm);
return;
}
// if (this->is_master) {
// MPI_Comm_split(this->group_comm, MPI_UNDEFINED, this->comm_rank,
// &dht_comm); return;
// }
if (!this->is_master) {
// if (!this->is_master) {
MPI_Comm_split(this->group_comm, 1, this->comm_rank, &dht_comm);
// MPI_Comm_split(this->group_comm, 1, this->comm_rank, &dht_comm);
auto map_copy = key_species;
// auto map_copy = key_species;
if (key_species.empty()) {
std::vector<std::uint32_t> default_signif(
this->prop_names.size(), DHT_Wrapper::DHT_KEY_SIGNIF_DEFAULT);
map_copy = NamedVector<std::uint32_t>(this->prop_names, default_signif);
}
// if (key_species.empty()) {
// std::vector<std::uint32_t> default_signif(
// this->prop_names.size(), DHT_Wrapper::DHT_KEY_SIGNIF_DEFAULT);
// map_copy = NamedVector<std::uint32_t>(this->prop_names,
// default_signif);
// }
auto key_indices = parseDHTSpeciesVec(key_species, this->prop_names);
// auto key_indices = parseDHTSpeciesVec(key_species, this->prop_names);
if (this->dht) {
delete this->dht;
}
// if (this->dht) {
// delete this->dht;
// }
const std::uint64_t dht_size = size_mb * MB_FACTOR;
// const std::uint64_t dht_size = size_mb * MB_FACTOR;
this->dht = new DHT_Wrapper(dht_comm, dht_size, map_copy, key_indices,
this->prop_names, params.hooks,
this->prop_count, params.use_interp);
this->dht->setBaseTotals(base_totals.at(0), base_totals.at(1));
}
// this->dht = new DHT_Wrapper(dht_comm, dht_size, map_copy, key_indices,
// this->prop_names, params.hooks,
// this->prop_count, params.use_interp);
// this->dht->setBaseTotals(base_totals.at(0), base_totals.at(1));
// }
}
inline std::vector<std::int32_t> poet::ChemistryModule::parseDHTSpeciesVec(
@ -457,42 +293,42 @@ void poet::ChemistryModule::initializeInterp(
std::uint32_t bucket_size, std::uint32_t size_mb, std::uint32_t min_entries,
const NamedVector<std::uint32_t> &key_species) {
if (!this->is_master) {
// if (!this->is_master) {
constexpr uint32_t MB_FACTOR = 1E6;
// constexpr uint32_t MB_FACTOR = 1E6;
assert(this->dht);
// assert(this->dht);
this->interp_enabled = true;
// this->interp_enabled = true;
auto map_copy = key_species;
// auto map_copy = key_species;
if (key_species.empty()) {
map_copy = this->dht->getKeySpecies();
for (std::size_t i = 0; i < map_copy.size(); i++) {
const std::uint32_t signif =
map_copy[i] - (map_copy[i] > InterpolationModule::COARSE_DIFF
? InterpolationModule::COARSE_DIFF
: 0);
map_copy[i] = signif;
}
}
// if (key_species.empty()) {
// map_copy = this->dht->getKeySpecies();
// for (std::size_t i = 0; i < map_copy.size(); i++) {
// const std::uint32_t signif =
// map_copy[i] - (map_copy[i] > InterpolationModule::COARSE_DIFF
// ? InterpolationModule::COARSE_DIFF
// : 0);
// map_copy[i] = signif;
// }
// }
auto key_indices =
parseDHTSpeciesVec(map_copy, dht->getKeySpecies().getNames());
// auto key_indices =
// parseDHTSpeciesVec(map_copy, dht->getKeySpecies().getNames());
if (this->interp) {
this->interp.reset();
}
// if (this->interp) {
// this->interp.reset();
// }
const uint64_t pht_size = size_mb * MB_FACTOR;
// const uint64_t pht_size = size_mb * MB_FACTOR;
interp = std::make_unique<poet::InterpolationModule>(
bucket_size, pht_size, min_entries, *(this->dht), map_copy, key_indices,
this->prop_names, this->params.hooks);
// interp = std::make_unique<poet::InterpolationModule>(
// bucket_size, pht_size, min_entries, *(this->dht), map_copy,
// key_indices, this->prop_names, this->params.hooks);
interp->setInterpolationFunction(inverseDistanceWeighting);
}
// interp->setInterpolationFunction(inverseDistanceWeighting);
// }
}
std::vector<double>

View File

@ -2,16 +2,17 @@
#ifndef CHEMISTRYMODULE_H_
#define CHEMISTRYMODULE_H_
#include "../Base/SimParams.hpp"
#include "../DataStructures/DataStructures.hpp"
#include "DataStructures/Field.hpp"
#include "DataStructures/NamedVector.hpp"
#include "ChemistryDefs.hpp"
#include "Init/InitialList.hpp"
#include "SurrogateModels/DHT_Wrapper.hpp"
#include "SurrogateModels/Interpolation.hpp"
#include <IrmResult.h>
#include <PhreeqcRM.h>
#include "IPhreeqcPOET.hpp"
#include <array>
#include <cstddef>
#include <cstdint>
#include <map>
#include <memory>
@ -24,7 +25,7 @@ namespace poet {
* \brief Wrapper around PhreeqcRM to provide POET specific parallelization with
* easy access.
*/
class ChemistryModule : public PhreeqcRM {
class ChemistryModule {
public:
/**
* Creates a new instance of Chemistry module with given grid cell count, work
@ -41,45 +42,40 @@ public:
* \param wp_size Count of grid cells to fill each work package at maximum.
* \param communicator MPI communicator to distribute work in.
*/
ChemistryModule(uint32_t nxyz, uint32_t wp_size, std::uint32_t maxiter,
const ChemistryParams &chem_param, MPI_Comm communicator);
ChemistryModule(uint32_t wp_size,
const InitialList::ChemistryInit &chem_params,
MPI_Comm communicator);
/**
* Deconstructor, which frees DHT data structure if used.
*/
~ChemistryModule();
/**
* Parses the input script and extract information needed during runtime.
*
* **Only run by master**.
*
* Database must be loaded beforehand.
*
* \param input_script_path Path to input script to parse.
*/
void RunInitFile(const std::string &input_script_path);
void masterSetField(Field field);
/**
* Run the chemical simulation with parameters set.
*/
void RunCells();
void simulate(double dt);
/**
* Returns the chemical field.
*/
auto GetField() const { return this->chem_field; }
auto &GetField() { return this->chem_field; }
/**
* Returns all known species names, including not only aqueous species, but
* also equilibrium, exchange, surface and kinetic reactants.
*/
auto GetPropNames() const { return this->prop_names; }
// auto GetPropNames() const { return this->prop_names; }
/**
* Return the accumulated runtime in seconds for chemical simulation.
*/
auto GetChemistryTime() const { return this->chem_t; }
void setFilePadding(std::uint32_t maxiter) {
this->file_pad = std::ceil(std::log10(maxiter + 1));
}
/**
* Create a new worker instance inside given MPI communicator.
*
@ -118,17 +114,6 @@ public:
DHT_SNAPS_ITEREND //!< output snapshots after each iteration
};
/**
* **This function has to be run!**
*
* Merge initial values from existing module with the chemistry module and set
* according internal variables.
*
* \param other Field to merge chemistry with. Most likely it is something
* like the diffusion field.
*/
void initializeField(const Field &other);
/**
* **Only called by workers!** Start the worker listening loop.
*/
@ -243,9 +228,7 @@ protected:
const NamedVector<std::uint32_t> &key_species);
enum {
CHEM_INIT,
CHEM_WP_SIZE,
CHEM_INIT_SPECIES,
CHEM_FIELD_INIT,
CHEM_DHT_ENABLE,
CHEM_DHT_SIGNIF_VEC,
CHEM_DHT_SNAPS,
@ -294,7 +277,7 @@ protected:
using worker_list_t = std::vector<struct worker_info_s>;
using workpointer_t = std::vector<double>::iterator;
void MasterRunParallel();
void MasterRunParallel(double dt);
void MasterRunSequential();
void MasterSendPkgs(worker_list_t &w_list, workpointer_t &work_pointer,
@ -320,8 +303,8 @@ protected:
void WorkerPerfToMaster(int type, const struct worker_s &timings);
void WorkerMetricsToMaster(int type);
IRM_RESULT WorkerRunWorkPackage(WorkPackage &work_package, double dSimTime,
double dTimestep);
void WorkerRunWorkPackage(WorkPackage &work_package, double dSimTime,
double dTimestep);
std::vector<uint32_t> CalculateWPSizesVector(uint32_t n_cells,
uint32_t wp_size) const;
@ -372,21 +355,18 @@ protected:
bool print_progessbar{false};
std::uint32_t file_pad;
std::uint32_t file_pad{1};
double chem_t{0.};
uint32_t n_cells = 0;
uint32_t prop_count = 0;
std::vector<std::string> prop_names;
Field chem_field;
static constexpr int MODULE_COUNT = 5;
const InitialList::ChemistryInit &params;
const ChemistryParams &params;
std::array<std::uint32_t, MODULE_COUNT> speciesPerModule{};
std::map<int, std::unique_ptr<IPhreeqcPOET>> phreeqc_instances;
};
} // namespace poet

View File

@ -1,12 +1,9 @@
#include "ChemistryModule.hpp"
#include <PhreeqcRM.h>
#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <mpi.h>
#include <stdexcept>
#include <vector>
std::vector<uint32_t>
@ -308,45 +305,46 @@ inline void poet::ChemistryModule::MasterRecvPkgs(worker_list_t &w_list,
}
}
void poet::ChemistryModule::RunCells() {
void poet::ChemistryModule::simulate(double dt) {
double start_t{MPI_Wtime()};
if (this->is_sequential) {
MasterRunSequential();
return;
}
MasterRunParallel();
MasterRunParallel(dt);
double end_t{MPI_Wtime()};
this->chem_t += end_t - start_t;
}
void poet::ChemistryModule::MasterRunSequential() {
std::vector<double> shuffled_field =
shuffleField(chem_field.AsVector(), n_cells, prop_count, 1);
// std::vector<double> shuffled_field =
// shuffleField(chem_field.AsVector(), n_cells, prop_count, 1);
std::vector<std::vector<double>> input;
for (std::size_t i = 0; i < n_cells; i++) {
input.push_back(
std::vector<double>(shuffled_field.begin() + (i * prop_count),
shuffled_field.begin() + ((i + 1) * prop_count)));
}
// std::vector<std::vector<double>> input;
// for (std::size_t i = 0; i < n_cells; i++) {
// input.push_back(
// std::vector<double>(shuffled_field.begin() + (i * prop_count),
// shuffled_field.begin() + ((i + 1) *
// prop_count)));
// }
this->setDumpedField(input);
PhreeqcRM::RunCells();
this->getDumpedField(input);
// this->setDumpedField(input);
// PhreeqcRM::RunCells();
// this->getDumpedField(input);
shuffled_field.clear();
for (std::size_t i = 0; i < n_cells; i++) {
shuffled_field.insert(shuffled_field.end(), input[i].begin(),
input[i].end());
}
// shuffled_field.clear();
// for (std::size_t i = 0; i < n_cells; i++) {
// shuffled_field.insert(shuffled_field.end(), input[i].begin(),
// input[i].end());
// }
std::vector<double> out_vec{shuffled_field};
unshuffleField(shuffled_field, n_cells, prop_count, 1, out_vec);
chem_field = out_vec;
// std::vector<double> out_vec{shuffled_field};
// unshuffleField(shuffled_field, n_cells, prop_count, 1, out_vec);
// chem_field = out_vec;
}
void poet::ChemistryModule::MasterRunParallel() {
void poet::ChemistryModule::MasterRunParallel(double dt) {
/* declare most of the needed variables here */
double seq_a, seq_b, seq_c, seq_d;
double worker_chemistry_a, worker_chemistry_b;
@ -360,7 +358,6 @@ void poet::ChemistryModule::MasterRunParallel() {
MPI_Barrier(this->group_comm);
double dt = this->PhreeqcRM::GetTimeStep();
static uint32_t iteration = 0;
/* start time measurement of sequential part */
@ -466,3 +463,13 @@ poet::ChemistryModule::CalculateWPSizesVector(uint32_t n_cells,
return wp_sizes_vector;
}
void poet::ChemistryModule::masterSetField(Field field) {
this->chem_field = field;
this->prop_count = field.GetProps().size();
int ftype = CHEM_FIELD_INIT;
PropagateFunctionType(ftype);
ChemBCast(&this->prop_count, 1, MPI_UINT32_T);
}

View File

@ -22,13 +22,14 @@
#include "DHT_Wrapper.hpp"
#include "Rounding.hpp"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <cstring>
#include <iostream>
#include <stdexcept>
#include <vector>

View File

@ -23,19 +23,18 @@
#ifndef DHT_WRAPPER_H
#define DHT_WRAPPER_H
#include "../../Base/RInsidePOET.hpp"
#include "Base/RInsidePOET.hpp"
#include "DataStructures/NamedVector.hpp"
#include "../../Base/SimParams.hpp"
#include "../../DataStructures/DataStructures.hpp"
#include "../enums.hpp"
#include "HashFunctions.hpp"
#include "Chemistry/ChemistryDefs.hpp"
#include "LookupKey.hpp"
#include "Rounding.hpp"
#include <array>
#include <cstdint>
#include <limits>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>

View File

@ -3,14 +3,12 @@
#ifndef INTERPOLATION_H_
#define INTERPOLATION_H_
#include "../../Base/SimParams.hpp"
#include "../../DataStructures/DataStructures.hpp"
#include "DataStructures/NamedVector.hpp"
#include "DHT_Wrapper.hpp"
#include "LookupKey.hpp"
#include "Rounding.hpp"
#include <cassert>
#include <iostream>
#include <list>
#include <memory>
#include <mpi.h>
@ -22,7 +20,6 @@ extern "C" {
}
#include <cstdint>
#include <functional>
#include <unordered_map>
#include <vector>

View File

@ -1,8 +1,8 @@
// Time-stamp: "Last modified 2023-08-16 17:02:31 mluebke"
#include "Interpolation.hpp"
#include "../../DataStructures/DataStructures.hpp"
#include "DHT_Wrapper.hpp"
#include "DataStructures/NamedVector.hpp"
#include "HashFunctions.hpp"
#include "LookupKey.hpp"
#include "Rounding.hpp"

View File

@ -2,10 +2,8 @@
#include "SurrogateModels/DHT_Wrapper.hpp"
#include "SurrogateModels/Interpolation.hpp"
#include <IrmResult.h>
#include <algorithm>
#include <cassert>
#include <cmath>
#include "Chemistry/ChemistryDefs.hpp"
#include <cstddef>
#include <cstdint>
#include <iomanip>
@ -18,19 +16,6 @@
namespace poet {
struct WorkPackage {
std::size_t size;
std::vector<std::vector<double>> input;
std::vector<std::vector<double>> output;
std::vector<std::uint8_t> mapping;
WorkPackage(size_t _size) : size(_size) {
input.resize(size);
output.resize(size);
mapping.resize(size, CHEM_PQC);
}
};
inline std::string get_string(int root, MPI_Comm communicator) {
int count;
MPI_Bcast(&count, 1, MPI_INT, root, communicator);
@ -59,13 +44,8 @@ void poet::ChemistryModule::WorkerLoop() {
PropagateFunctionType(func_type);
switch (func_type) {
case CHEM_INIT: {
RunInitFile(get_string(0, this->group_comm));
break;
}
case CHEM_INIT_SPECIES: {
Field dummy;
initializeField(dummy);
case CHEM_FIELD_INIT: {
ChemBCast(&this->prop_count, 1, MPI_UINT32_T);
break;
}
case CHEM_WORK_LOOP: {
@ -191,9 +171,7 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status,
phreeqc_time_start = MPI_Wtime();
if (WorkerRunWorkPackage(s_curr_wp, current_sim_time, dt) != IRM_OK) {
std::cerr << "Phreeqc error" << std::endl;
};
WorkerRunWorkPackage(s_curr_wp, current_sim_time, dt);
phreeqc_time_end = MPI_Wtime();
@ -299,42 +277,65 @@ void poet::ChemistryModule::WorkerReadDHTDump(
}
}
IRM_RESULT
poet::ChemistryModule::WorkerRunWorkPackage(WorkPackage &work_package,
double dSimTime, double dTimestep) {
void poet::ChemistryModule::WorkerRunWorkPackage(WorkPackage &work_package,
double dSimTime,
double dTimestep) {
// check if we actually need to start phreeqc
std::vector<std::uint32_t> pqc_mapping;
bool queued_cell = false;
for (std::size_t i = 0; i < work_package.size; i++) {
if (work_package.mapping[i] == CHEM_PQC) {
pqc_mapping.push_back(i);
for (std::size_t wp_id = 0; wp_id < work_package.size; wp_id++) {
if (work_package.mapping[wp_id] != CHEM_PQC) {
continue;
}
const auto &input = work_package.input[wp_id];
const auto &sol_id = input[0];
auto &phreeqc_instance = this->phreeqc_instances[sol_id];
work_package.output[wp_id] = work_package.input[wp_id];
work_package.input[wp_id].erase(work_package.input[wp_id].begin());
// remove NaNs from the input
work_package.input[wp_id].erase(
std::remove_if(work_package.input[wp_id].begin(),
work_package.input[wp_id].end(),
[](double d) { return std::isnan(d); }),
work_package.input[wp_id].end());
phreeqc_instance->queueCell(work_package.input[wp_id]);
queued_cell = true;
}
if (!queued_cell) {
return;
}
std::map<int, std::vector<std::size_t>> zone_mapping;
// run the phreeqc instances
for (const auto &[pqc_id, phreeqc_instance] : this->phreeqc_instances) {
if (zone_mapping.find(pqc_id) == zone_mapping.end()) {
continue;
}
phreeqc_instance->runQueuedCells(dTimestep);
// remap the output to the work_package
std::vector<std::vector<double>> pqc_out;
phreeqc_instance->dequeueCells(pqc_out);
std::size_t output_id = 0;
for (const auto &wp_id : zone_mapping[pqc_id]) {
std::size_t output_index = 0;
for (std::size_t i = 1; i < work_package.output[wp_id].size(); i++) {
if (!(std::isnan(work_package.output[wp_id][i]))) {
work_package.output[wp_id][i] = pqc_out[output_id][output_index++];
}
}
output_id++;
}
}
if (pqc_mapping.empty()) {
return IRM_OK;
}
IRM_RESULT result;
this->PhreeqcRM::setPOETMapping(pqc_mapping);
this->setDumpedField(work_package.input);
this->PhreeqcRM::SetTime(dSimTime);
this->PhreeqcRM::SetTimeStep(dTimestep);
result = this->PhreeqcRM::RunCells();
std::vector<std::vector<double>> output_tmp(work_package.size);
this->getDumpedField(output_tmp);
for (std::size_t i = 0; i < work_package.size; i++) {
if (work_package.mapping[i] == CHEM_PQC) {
work_package.output[i] = output_tmp[i];
}
}
return result;
}
void poet::ChemistryModule::WorkerPerfToMaster(int type,

View File

@ -1,10 +0,0 @@
#ifndef ENUMS_H_
#define ENUMS_H_
namespace poet {
enum DHT_PROP_TYPES { DHT_TYPE_DEFAULT, DHT_TYPE_CHARGE, DHT_TYPE_TOTAL };
enum CHEMISTRY_OUT_SOURCE { CHEM_PQC, CHEM_DHT, CHEM_INTERP };
} // namespace poet
#endif // ENUMS_H_

View File

@ -4,6 +4,7 @@
#include <Rcpp.h>
#include <cstddef>
#include <cstdint>
#include <string>
#include <vector>

View File

@ -4,8 +4,6 @@ namespace poet {
InitialList::ChemistryInit InitialList::getChemistryInit() const {
ChemistryInit chem_init;
chem_init.initial_grid = Field(initial_grid);
chem_init.total_grid_cells = this->n_cols * this->n_rows;
chem_init.database = database;

View File

@ -39,9 +39,9 @@ static Rcpp::NumericMatrix pqcScriptToGrid(IPhreeqcPOET &phreeqc, RInside &R) {
static inline Rcpp::List matToGrid(RInside &R, const Rcpp::NumericMatrix &mat,
const Rcpp::NumericMatrix &grid) {
Rcpp::Function pqc_to_grid("pqc_to_grid");
Rcpp::Function pqc_to_grid_R("pqc_to_grid");
return pqc_to_grid(mat, grid);
return pqc_to_grid_R(mat, grid);
}
static inline std::map<int, std::string>

View File

@ -160,8 +160,6 @@ private:
public:
struct ChemistryInit {
Field initial_grid;
uint32_t total_grid_cells;
std::string database;

View File

@ -18,7 +18,9 @@ int main(int argc, char **argv) {
R.parseEvalQ(init_r_library);
const std::string script = argv[1];
R.parseEvalQ("source('" + script + "')");
Rcpp::Function source_R("source");
source_R(script);
Rcpp::List setup = R["setup"];

View File

@ -23,6 +23,7 @@
#include "Base/Macros.hpp"
#include "Base/RInsidePOET.hpp"
// #include "Chemistry/ChemistryModule.hpp"
#include "Chemistry/ChemistryModule.hpp"
#include "DataStructures/Field.hpp"
#include "Init/InitialList.hpp"
#include "Transport/DiffusionModule.hpp"
@ -48,41 +49,16 @@ using namespace Rcpp;
static int MY_RANK = 0;
// poet::ChemistryModule::SingleCMap DFToHashMap(const Rcpp::DataFrame &df) {
// std::unordered_map<std::string, double> out_map;
// vector<string> col_names = Rcpp::as<vector<string>>(df.names());
// for (const auto &name : col_names) {
// double val = df[name.c_str()];
// out_map.insert({name, val});
// }
// return out_map;
// }
// HACK: this is a step back as the order and also the count of fields is
// predefined, but it will change in the future
void writeFieldsToR(RInside &R, const Field &trans) {
// , const Field &chem) {
void writeFieldsToR(RInside &R, const Field &trans, const Field &chem) {
Rcpp::DataFrame t_field(trans.asSEXP());
R["TMP"] = t_field;
R.parseEval("mysetup$state_T <- TMP");
R["TMP"] = chem.asSEXP();
R.parseEval("mysetup$state_C <- TMP");
// R["TMP"] = Rcpp::wrap(trans.AsVector());
// R["TMP_PROPS"] = Rcpp::wrap(trans.GetProps());
// R.parseEval(std::string(
// "mysetup$state_T <- setNames(data.frame(matrix(TMP, nrow=" +
// std::to_string(trans.GetRequestedVecSize()) + ")), TMP_PROPS)"));
// R["TMP"] = Rcpp::wrap(chem.AsVector());
// R["TMP_PROPS"] = Rcpp::wrap(chem.GetProps());
// R.parseEval(std::string(
// "mysetup$state_C <- setNames(data.frame(matrix(TMP, nrow=" +
// std::to_string(chem.GetRequestedVecSize()) + ")), TMP_PROPS)"));
}
enum ParseRet { PARSER_OK, PARSER_ERROR, PARSER_HELP };
@ -102,11 +78,12 @@ ParseRet parseInitValues(char **argv, RInsidePOET &R,
return ParseRet::PARSER_HELP;
}
// if positional arguments are missing
else if (!cmdl(2)) {
if (!cmdl(3)) {
if (MY_RANK == 0) {
ERRMSG("POET needs 2 positional arguments: ");
ERRMSG("1) the R script defining your simulation and");
ERRMSG("2) the directory prefix where to save results and profiling");
ERRMSG("POET needs 3 positional arguments: ");
ERRMSG("1) the R script defining your simulation runtime.");
ERRMSG("2) the initial .rds file generated by poet_init.");
ERRMSG("3) the directory prefix where to save results and profiling");
}
return ParseRet::PARSER_ERROR;
}
@ -203,8 +180,8 @@ ParseRet parseInitValues(char **argv, RInsidePOET &R,
std::string runtime_file;
std::string out_dir;
cmdl(1) >> init_file;
cmdl(2) >> runtime_file;
cmdl(1) >> runtime_file;
cmdl(2) >> init_file;
cmdl(3) >> out_dir;
// chem_params.dht_outdir = out_dir;
@ -300,7 +277,7 @@ ParseRet parseInitValues(char **argv, RInsidePOET &R,
// }
static double RunMasterLoop(RInside &R, const RuntimeParameters &params,
DiffusionModule &diffusion) {
DiffusionModule &diffusion, ChemistryModule &chem) {
/* Iteration Count is dynamic, retrieving value from R (is only needed by
* master for the following loop) */
@ -320,7 +297,7 @@ static double RunMasterLoop(RInside &R, const RuntimeParameters &params,
// chem.initializeField(diffusion.getField());
// if (params.getNumParams().print_progressbar) {
// chem.setProgressBarPrintout(true);
chem.setProgressBarPrintout(true);
// }
/* SIMULATION LOOP */
@ -332,7 +309,7 @@ static double RunMasterLoop(RInside &R, const RuntimeParameters &params,
// cout << "CPP: Evaluating next time step" << endl;
// R.parseEvalQ("mysetup <- master_iteration_setup(mysetup)");
double dt = params.timesteps[iter - 1];
const double &dt = params.timesteps[iter - 1];
// cout << "CPP: Next time step is " << dt << "[s]" << endl;
MSG("Next time step is " + std::to_string(dt) + " [s]");
@ -347,14 +324,15 @@ static double RunMasterLoop(RInside &R, const RuntimeParameters &params,
// TODO: transport to diffusion
diffusion.simulate(dt);
// chem.getField().update(diffusion.getField());
chem.GetField().update(diffusion.getField());
// chem.getfield().update(diffusion.getfield());
MSG("Chemistry step");
// chem.SetTimeStep(dt);
// chem.RunCells();
chem.simulate(dt);
writeFieldsToR(R, diffusion.getField());
writeFieldsToR(R, diffusion.getField(), chem.GetField());
// diffusion.getField().update(chem.GetField());
R["req_dt"] = dt;
@ -431,7 +409,7 @@ static double RunMasterLoop(RInside &R, const RuntimeParameters &params,
// R.parseEvalQ("profiling$interp_cached <- interp_cached");
// }
// chem.MasterLoopBreak();
chem.MasterLoopBreak();
return dSimTime;
}
@ -453,37 +431,6 @@ int main(int argc, char *argv[]) {
MSG("Running POET version " + std::string(poet_version));
}
if (world_rank > 0) {
{
// SimParams params(world_rank, world_size);
// int pret = params.parseFromCmdl(argv, R);
// if (pret == poet::PARSER_ERROR) {
// MPI_Finalize();
// return EXIT_FAILURE;
// } else if (pret == poet::PARSER_HELP) {
// MPI_Finalize();
// return EXIT_SUCCESS;
// }
// // ChemistryModule worker(nxyz, nxyz, MPI_COMM_WORLD);
// ChemistryModule worker = poet::ChemistryModule::createWorker(
// MPI_COMM_WORLD, params.getChemParams());
// set_chem_parameters(worker, worker.GetWPSize(),
// params.getChemParams().database_path);
// worker.WorkerLoop();
}
MPI_Barrier(MPI_COMM_WORLD);
MSG("finished, cleanup of process " + std::to_string(world_rank));
MPI_Finalize();
return EXIT_SUCCESS;
}
/*Loading Dependencies*/
// TODO: kann raus
R.parseEvalQ(kin_r_library);
@ -499,8 +446,25 @@ int main(int argc, char *argv[]) {
break;
}
InitialList init_list(R);
init_list.importList(run_params.init_params);
MSG("RInside initialized on process " + std::to_string(world_rank));
if (world_rank > 0) {
ChemistryModule worker(1, init_list.getChemistryInit(), MPI_COMM_WORLD);
worker.WorkerLoop();
MPI_Barrier(MPI_COMM_WORLD);
MSG("finished, cleanup of process " + std::to_string(world_rank));
MPI_Finalize();
return EXIT_SUCCESS;
}
// R.parseEvalQ("mysetup <- setup");
// // if (world_rank == 0) { // get timestep vector from
// // grid_init function ... //
@ -521,13 +485,14 @@ int main(int argc, char *argv[]) {
// MPI_Barrier(MPI_COMM_WORLD);
InitialList init_list(R);
init_list.importList(run_params.init_params);
DiffusionModule diffusion(init_list.getDiffusionInit(),
init_list.getInitialGrid());
dSimTime = RunMasterLoop(R, run_params, diffusion);
ChemistryModule chemistry(1, init_list.getChemistryInit(), MPI_COMM_WORLD);
chemistry.masterSetField(init_list.getInitialGrid());
dSimTime = RunMasterLoop(R, run_params, diffusion, chemistry);
MSG("finished simulation loop");