Merge branch 'fix-copy-fields' into 'main'

fix: do pre-copy of fields in grid class

See merge request sec34/port!13
This commit is contained in:
Max Lübke 2023-01-23 16:41:12 +01:00
commit 6a911962db
6 changed files with 48 additions and 123 deletions

View File

@ -56,6 +56,7 @@ inline double RunMasterLoop(SimParams &params, RInside &R, Grid &grid,
double dStartTime = MPI_Wtime();
for (uint32_t iter = 1; iter < maxiter + 1; iter++) {
uint32_t tick = 0;
// cout << "CPP: Evaluating next time step" << endl;
// R.parseEvalQ("mysetup <- master_iteration_setup(mysetup)");
@ -72,9 +73,12 @@ inline double RunMasterLoop(SimParams &params, RInside &R, Grid &grid,
// TODO: transport to diffusion
diffusion.simulate(dt);
grid.PreModuleFieldCopy(tick++);
cout << "CPP: Chemistry" << endl;
C.Simulate(dt);
grid.PreModuleFieldCopy(tick++);
R["req_dt"] = dt;
R["simtime"] = (sim_time += dt);
@ -182,6 +186,10 @@ int main(int argc, char *argv[]) {
Grid grid;
grid.InitModuleFromParams(GridParams(R));
grid.PushbackModuleFlow(poet::DIFFUSION_MODULE_NAME,
poet::BaseChemModule::CHEMISTRY_MODULE_NAME);
grid.PushbackModuleFlow(poet::BaseChemModule::CHEMISTRY_MODULE_NAME,
poet::DIFFUSION_MODULE_NAME);
params.initVectorParams(R, grid.GetSpeciesCount());

View File

@ -26,6 +26,7 @@
#include <array>
#include <cstddef>
#include <cstdint>
#include <list>
#include <map>
#include <string>
#include <vector>
@ -41,6 +42,11 @@ using StateMemory = struct s_StateMemory {
std::vector<std::string> props;
};
using FlowInputOutputInfo = struct inOut_info {
std::string input_field;
std::string output_field;
};
constexpr const char *GRID_MODULE_NAME = "grid_init";
/**
@ -66,6 +72,9 @@ public:
void SetGridSize(double s_x, double s_y = 0., double s_z = 0.);
void SetPropNames(const std::vector<std::string> &prop_names);
void PushbackModuleFlow(const std::string &input, const std::string &output);
void PreModuleFieldCopy(uint32_t tick);
void InitGridFromScratch(const Rcpp::DataFrame &init_cell);
void InitGridFromRDS();
void InitGridFromPhreeqc();
@ -89,7 +98,7 @@ public:
-> std::vector<double>;
private:
// auto InitPhreeqc(const poet::GridParams &params) -> PhreeqcRM *;
std::vector<FlowInputOutputInfo> flow_vec;
std::uint8_t dim = 0;
std::array<double, MAX_DIM> grid_size;
@ -99,8 +108,6 @@ private:
StateMemory *grid_init = std::nullptr_t();
std::vector<std::string> prop_names;
// PhreeqcRM *phreeqc_rm = std::nullptr_t();
};
} // namespace poet
#endif // GRID_H

View File

@ -103,18 +103,6 @@ void ChemMaster::Simulate(double dt) {
std::vector<double> &field = this->state->mem;
for (uint32_t i = 0; i < this->prop_names.size(); i++) {
try {
std::vector<double> t_prop_vec = this->grid.GetSpeciesByName(
this->prop_names[i], poet::DIFFUSION_MODULE_NAME);
std::copy(t_prop_vec.begin(), t_prop_vec.end(),
field.begin() + (i * this->n_cells_per_prop));
} catch (...) {
continue;
}
}
// HACK: transfer the field into R data structure serving as input for phreeqc
R["TMP_T"] = field;

View File

@ -72,18 +72,6 @@ void ChemSeq::Simulate(double dt) {
std::vector<double> &field = this->state->mem;
for (uint32_t i = 0; i < this->prop_names.size(); i++) {
try {
std::vector<double> t_prop_vec = this->grid.GetSpeciesByName(
this->prop_names[i], poet::DIFFUSION_MODULE_NAME);
std::copy(t_prop_vec.begin(), t_prop_vec.end(),
field.begin() + (i * this->n_cells_per_prop));
} catch (...) {
continue;
}
}
// HACK: transfer the field into R data structure serving as input for phreeqc
R["TMP_T"] = field;

View File

@ -167,21 +167,6 @@ void DiffusionModule::simulate(double dt) {
std::vector<double> &curr_field = this->state->mem;
// copy output of another module (in this case there is only the chemistry
// module) as input for diffusion
for (uint32_t i = 0; i < this->prop_names.size(); i++) {
try {
std::vector<double> t_prop_vec = this->grid.GetSpeciesByName(
this->prop_names[i], poet::BaseChemModule::CHEMISTRY_MODULE_NAME);
std::copy(t_prop_vec.begin(), t_prop_vec.end(),
curr_field.begin() + (i * this->n_cells_per_prop));
} catch (...) {
// TODO: there might be something wrong ...
continue;
}
}
this->diff_input.setTimestep(dt);
double *field = curr_field.data();

View File

@ -101,6 +101,36 @@ void poet::Grid::SetPropNames(const std::vector<std::string> &prop_names) {
this->prop_names = prop_names;
}
void poet::Grid::PushbackModuleFlow(const std::string &input,
const std::string &output) {
FlowInputOutputInfo element = {input, output};
this->flow_vec.push_back(element);
}
void poet::Grid::PreModuleFieldCopy(uint32_t tick) {
FlowInputOutputInfo curr_element = this->flow_vec.at(tick);
const std::string input_module_name = curr_element.input_field;
StateMemory *out_state = this->GetStatePointer(curr_element.output_field);
std::vector<std::string> &mod_props = out_state->props;
std::vector<double> &mod_field = out_state->mem;
// copy output of another module as input for another module
for (uint32_t i = 0; i < mod_props.size(); i++) {
try {
std::vector<double> t_prop_vec =
this->GetSpeciesByName(mod_props[i], input_module_name);
std::copy(t_prop_vec.begin(), t_prop_vec.end(),
mod_field.begin() + (i * this->GetTotalCellCount()));
} catch (...) {
// TODO: there might be something wrong ...
continue;
}
}
}
Grid::Grid() {
this->n_cells.fill(0);
this->grid_size.fill(0);
@ -110,9 +140,6 @@ Grid::~Grid() {
for (auto &[key, val] : this->state_register) {
delete val;
}
// if (this->phreeqc_rm != nullptr) {
// delete this->phreeqc_rm;
// }
}
void poet::Grid::InitGridFromScratch(const Rcpp::DataFrame &init_cell) {
@ -168,20 +195,6 @@ poet::StateMemory *Grid::GetStatePointer(std::string module_name) {
return it->second;
}
// void Grid::DeregisterState(std::string module_name) {
// const auto it = this->state_register.find(module_name);
// if (it == this->state_register.end()) {
// throw std::range_error(
// std::string("Requested module " + module_name + " was not found"));
// }
// auto *p = it->second;
// delete p;
// this->state_register.erase(it);
// }
auto Grid::GetGridSize(uint8_t direction) const -> uint32_t {
return this->grid_size.at(direction);
}
@ -217,20 +230,6 @@ auto poet::Grid::GetSpeciesByName(std::string name,
std::string module_name) const
-> std::vector<double> {
// if module name references to initial grid
// if (module_name == poet::GRID_MODULE_NAME) {
// auto const it = this->grid_init.find(name);
// if (it == this->grid_init.end()) {
// throw std::range_error(
// std::string("Species " + name + " was not found for initial
// grid"));
// }
// return it->second;
// }
// else find module with name
auto const it = this->state_register.find(module_name);
if (it == this->state_register.end()) {
@ -254,53 +253,3 @@ auto poet::Grid::GetSpeciesByName(std::string name,
return std::vector<double>(module_memory->mem.begin() + begin_vec,
module_memory->mem.begin() + end_vec);
}
// PhreeqcRM *poet::Grid::InitPhreeqc(const poet::GridParams &params) {
// PhreeqcRM *chem_model = new PhreeqcRM(this->GetTotalCellCount(), 1);
// // FIXME: hardcoded options ...
// chem_model->SetErrorHandlerMode(1);
// chem_model->SetComponentH2O(false);
// chem_model->SetRebalanceFraction(0.5);
// chem_model->SetRebalanceByCell(true);
// chem_model->UseSolutionDensityVolume(false);
// chem_model->SetPartitionUZSolids(false);
// // Set concentration units
// // 1, mg/L; 2, mol/L; 3, kg/kgs
// chem_model->SetUnitsSolution(2);
// // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
// chem_model->SetUnitsPPassemblage(1);
// // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
// chem_model->SetUnitsExchange(1);
// // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
// chem_model->SetUnitsSurface(1);
// // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
// chem_model->SetUnitsGasPhase(1);
// // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
// chem_model->SetUnitsSSassemblage(1);
// // 0, mol/L cell; 1, mol/L water; 2 mol/L rock
// chem_model->SetUnitsKinetics(1);
// // Set representative volume
// std::vector<double> rv;
// rv.resize(this->GetTotalCellCount(), 1.0);
// chem_model->SetRepresentativeVolume(rv);
// // Set initial porosity
// std::vector<double> por;
// por.resize(this->GetTotalCellCount(), 0.05);
// chem_model->SetPorosity(por);
// // Set initial saturation
// std::vector<double> sat;
// sat.resize(this->GetTotalCellCount(), 1.0);
// chem_model->SetSaturation(sat);
// // Load database
// chem_model->LoadDatabase(params.database_path);
// chem_model->RunFile(true, true, true, params.input_script);
// chem_model->RunString(true, false, true, "DELETE; -all");
// return chem_model;
// }