mirror of
https://git.gfz-potsdam.de/naaice/poet.git
synced 2025-12-15 20:38:23 +01:00
BREAKING CHANGE: utilize PhreeqcRM and bring back old DHT bindings
This commit is contained in:
Max Luebke
Max Lübke
parent
8d13748545
commit
ded8fbd0ae
16
app/poet.cpp
16
app/poet.cpp
@ -32,8 +32,8 @@
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
#include <poet.h>
|
||||
#include <mpi.h>
|
||||
#include <poet.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace poet;
|
||||
@ -52,7 +52,7 @@ inline double RunMasterLoop(SimParams ¶ms, RInside &R, Grid &grid,
|
||||
C.InitModule(chem_params);
|
||||
/* SIMULATION LOOP */
|
||||
|
||||
double ret_time = MPI_Wtime();
|
||||
double dStartTime = MPI_Wtime();
|
||||
for (uint32_t iter = 1; iter < maxiter + 1; iter++) {
|
||||
// cout << "CPP: Evaluating next time step" << endl;
|
||||
// R.parseEvalQ("mysetup <- master_iteration_setup(mysetup)");
|
||||
@ -93,11 +93,11 @@ inline double RunMasterLoop(SimParams ¶ms, RInside &R, Grid &grid,
|
||||
C.End();
|
||||
diffusion.end();
|
||||
|
||||
return ret_time;
|
||||
return MPI_Wtime() - dStartTime;
|
||||
}
|
||||
|
||||
int main(int argc, char *argv[]) {
|
||||
double sim_start, sim_end;
|
||||
double dSimTime, sim_end;
|
||||
|
||||
int world_size, world_rank;
|
||||
|
||||
@ -195,18 +195,16 @@ int main(int argc, char *argv[]) {
|
||||
/* THIS IS EXECUTED BY THE MASTER */
|
||||
if (world_rank == 0) {
|
||||
if (world_size == 1) {
|
||||
sim_start = RunMasterLoop<ChemSeq>(params, R, grid, chem_params);
|
||||
dSimTime = RunMasterLoop<ChemSeq>(params, R, grid, chem_params);
|
||||
} else {
|
||||
sim_start = RunMasterLoop<ChemMaster>(params, R, grid, chem_params);
|
||||
dSimTime = RunMasterLoop<ChemMaster>(params, R, grid, chem_params);
|
||||
}
|
||||
|
||||
cout << "CPP: finished simulation loop" << endl;
|
||||
|
||||
sim_end = MPI_Wtime();
|
||||
|
||||
cout << "CPP: start timing profiling" << endl;
|
||||
|
||||
R["simtime"] = sim_end - sim_start;
|
||||
R["simtime"] = dSimTime;
|
||||
R.parseEvalQ("profiling$simtime <- simtime");
|
||||
|
||||
string r_vis_code;
|
||||
|
||||
@ -164,8 +164,8 @@ selout <- c(
|
||||
|
||||
|
||||
# Needed when using DHT
|
||||
signif_vector <- c(7, 7, 7, 7, 7, 7, 7, 5, 5)
|
||||
prop_type <- c("act", "act", "act", "act", "logact", "logact", "ignore", "act", "act")
|
||||
signif_vector <- c(7, 7, 7, 7, 7, 7, 7, 5, 5, 7)
|
||||
prop_type <- c("act", "act", "act", "act", "logact", "logact", "ignore", "act", "act", "act")
|
||||
prop <- names(init_cell)
|
||||
|
||||
chemistry <- list(
|
||||
|
||||
13
data/SimDol2D_diffu.yaml
Normal file
13
data/SimDol2D_diffu.yaml
Normal file
@ -0,0 +1,13 @@
|
||||
phreeqc:
|
||||
RebalanceByCell: True
|
||||
RebalanceFraction: 0.5
|
||||
SolutionDensityVolume: False
|
||||
PartitionUZSolids: False
|
||||
Units:
|
||||
Solution: "mg/L"
|
||||
PPassamblage: "mol/L"
|
||||
Exchange: "mol/L"
|
||||
Surface: "mol/L"
|
||||
GasPhase: "mol/L"
|
||||
SSassemblage: "mol/L"
|
||||
Kinetics: "mol/L"
|
||||
@ -1 +1 @@
|
||||
Subproject commit 0c160a810aeb64c359531c7f84911b9c1df42a9d
|
||||
Subproject commit a9e3d9fa355160ef7efdb5e421076e72c0e58893
|
||||
@ -27,12 +27,10 @@
|
||||
#include "RInside.h"
|
||||
#include "SimParams.hpp"
|
||||
#include "poet/PhreeqcWrapper.hpp"
|
||||
#include <PhreeqcRM.h>
|
||||
|
||||
#include <array>
|
||||
#include <bits/stdint-uintn.h>
|
||||
#include <cstdint>
|
||||
#include <mpi.h>
|
||||
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
@ -169,7 +167,7 @@ public:
|
||||
private:
|
||||
void printProgressbar(int count_pkgs, int n_wp, int barWidth = 70);
|
||||
inline void sendPkgs(int &pkg_to_send, int &count_pkgs, int &free_workers,
|
||||
double *work_pointer, const double &dt,
|
||||
double *&work_pointer, const double &dt,
|
||||
const uint32_t iteration);
|
||||
inline void recvPkgs(int &pkg_to_recv, bool to_send, int &free_workers);
|
||||
void shuffleField(const std::vector<double> &in_field, uint32_t size_per_prop,
|
||||
@ -251,6 +249,7 @@ private:
|
||||
uint32_t wp_size;
|
||||
uint32_t dht_size_per_process;
|
||||
uint32_t iteration = 0;
|
||||
|
||||
bool dht_enabled;
|
||||
bool dt_differ;
|
||||
int dht_snaps;
|
||||
|
||||
@ -23,6 +23,7 @@
|
||||
|
||||
#include "SimParams.hpp"
|
||||
|
||||
#include <cstdint>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
@ -38,8 +39,8 @@ extern "C" {
|
||||
* Macro to round a double value by cutting every digit after significant digit
|
||||
*
|
||||
*/
|
||||
#define ROUND(value, signif) \
|
||||
(((int)(pow(10.0, (double)signif) * value)) * pow(10.0, (double)-signif))
|
||||
// #define ROUND(value, signif) \
|
||||
// (((int)(pow(10.0, (double)signif) * value)) * pow(10.0, (double)-signif))
|
||||
|
||||
namespace poet {
|
||||
/**
|
||||
@ -63,7 +64,7 @@ static uint64_t get_md5(int key_size, void *key);
|
||||
*
|
||||
*/
|
||||
class DHT_Wrapper {
|
||||
public:
|
||||
public:
|
||||
/**
|
||||
* @brief Construct a new dht wrapper object
|
||||
*
|
||||
@ -107,8 +108,9 @@ class DHT_Wrapper {
|
||||
* @param[in,out] work_package Pointer to current work package
|
||||
* @param dt Current timestep of simulation
|
||||
*/
|
||||
void checkDHT(int length, std::vector<bool> &out_result_index,
|
||||
double *work_package, double dt);
|
||||
auto checkDHT(int length, std::vector<bool> &out_result_index,
|
||||
double *work_package, double dt)
|
||||
-> std::vector<std::vector<double>>;
|
||||
|
||||
/**
|
||||
* @brief Write simulated values into DHT
|
||||
@ -181,7 +183,7 @@ class DHT_Wrapper {
|
||||
*/
|
||||
uint64_t getEvictions();
|
||||
|
||||
private:
|
||||
private:
|
||||
/**
|
||||
* @brief Transform given workpackage into DHT key
|
||||
*
|
||||
@ -209,7 +211,7 @@ class DHT_Wrapper {
|
||||
* @param key Pointer to work package handled as the key
|
||||
* @param dt Current time step of the simulation
|
||||
*/
|
||||
void fuzzForDHT(int var_count, void *key, double dt);
|
||||
std::vector<double> fuzzForDHT(int var_count, void *key, double dt);
|
||||
|
||||
/**
|
||||
* @brief DHT handle
|
||||
@ -295,6 +297,6 @@ class DHT_Wrapper {
|
||||
*/
|
||||
std::vector<std::string> dht_prop_type_vector;
|
||||
};
|
||||
} // namespace poet
|
||||
} // namespace poet
|
||||
|
||||
#endif // DHT_WRAPPER_H
|
||||
#endif // DHT_WRAPPER_H
|
||||
|
||||
@ -21,11 +21,10 @@
|
||||
#ifndef GRID_H
|
||||
#define GRID_H
|
||||
|
||||
#include "PhreeqcRM.h"
|
||||
#include "poet/SimParams.hpp"
|
||||
#include <Rcpp.h>
|
||||
#include <array>
|
||||
#include <bits/stdint-intn.h>
|
||||
#include <bits/stdint-uintn.h>
|
||||
#include <cstdint>
|
||||
#include <map>
|
||||
#include <string>
|
||||
@ -97,16 +96,22 @@ public:
|
||||
-> std::vector<double>;
|
||||
|
||||
private:
|
||||
PhreeqcRM *initPhreeqc(const poet::GridParams ¶ms);
|
||||
|
||||
std::uint8_t dim;
|
||||
std::array<std::uint32_t, MAX_DIM> d_spatial;
|
||||
std::array<std::uint32_t, MAX_DIM> n_cells;
|
||||
|
||||
std::map<std::string, StateMemory *> state_register;
|
||||
std::map<std::string, std::vector<double>> grid_init;
|
||||
|
||||
prop_vec prop_names;
|
||||
|
||||
std::map<std::string, std::vector<double>> grid_init;
|
||||
PhreeqcRM *phreeqc_rm = std::nullptr_t();
|
||||
|
||||
void initFromPhreeqc(const poet::GridParams &grid_args);
|
||||
void initFromRDS(const poet::GridParams &grid_args);
|
||||
void initFromScratch(const poet::GridParams &grid_args);
|
||||
};
|
||||
} // namespace poet
|
||||
#endif // GRID_H
|
||||
|
||||
@ -4,9 +4,10 @@
|
||||
#include "IrmResult.h"
|
||||
#include "poet/SimParams.hpp"
|
||||
#include <PhreeqcRM.h>
|
||||
#include <bits/stdint-uintn.h>
|
||||
#include <cstdint>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
namespace poet {
|
||||
class PhreeqcWrapper : public PhreeqcRM {
|
||||
public:
|
||||
@ -20,14 +21,23 @@ public:
|
||||
auto GetSpeciesNamesFull() -> const std::vector<std::string> &;
|
||||
|
||||
void SetInternalsFromWP(const std::vector<double> &vecWP,
|
||||
uint32_t iCurrWPSize);
|
||||
auto GetWPFromInternals(uint32_t iCurrWPSize) -> std::vector<double>;
|
||||
uint32_t iCurrWPSize);
|
||||
void GetWPFromInternals(std::vector<double> &vecWP, uint32_t iCurrWPSize);
|
||||
|
||||
auto ReplaceTotalsByPotentials(const std::vector<double> &vecWP,
|
||||
uint32_t iCurrWPSize) -> std::vector<double>;
|
||||
|
||||
IRM_RESULT RunWorkPackage(std::vector<double> &vecWP,
|
||||
std::vector<int32_t> &vecMapping, double dSimTime,
|
||||
double dTimestep);
|
||||
|
||||
private:
|
||||
void InitInternals();
|
||||
void SetMappingForWP(uint32_t iCurrWPSize);
|
||||
inline void SetMappingForWP(uint32_t iCurrWPSize);
|
||||
|
||||
static constexpr int MODULE_COUNT = 5;
|
||||
static constexpr uint32_t DHT_SELOUT = 100;
|
||||
|
||||
uint32_t iWPSize;
|
||||
bool bInactiveCells = false;
|
||||
uint32_t iSpeciesCount;
|
||||
|
||||
@ -18,17 +18,16 @@
|
||||
** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
||||
*/
|
||||
|
||||
#include "poet/ChemSimPar.hpp"
|
||||
#include "poet/ChemSimSeq.hpp"
|
||||
#include "poet/DiffusionModule.hpp"
|
||||
#include "poet/Grid.hpp"
|
||||
#include "poet/SimParams.hpp"
|
||||
|
||||
#include <Rcpp.h>
|
||||
|
||||
#include <array>
|
||||
#include <bits/stdint-uintn.h>
|
||||
#include <cstring>
|
||||
#include <iostream>
|
||||
|
||||
#include <poet/ChemSimPar.hpp>
|
||||
#include <poet/Grid.hpp>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
@ -47,19 +46,35 @@ ChemMaster::ChemMaster(SimParams ¶ms, RInside &R_, Grid &grid_)
|
||||
|
||||
/* allocate memory */
|
||||
this->workerlist = new worker_struct[this->world_size - 1];
|
||||
std::memset(this->workerlist, '\0', sizeof(worker_struct) * (world_size - 1));
|
||||
this->send_buffer =
|
||||
new double[this->wp_size * this->prop_names.size() + BUFFER_OFFSET];
|
||||
this->mpi_buffer = new double[grid_size];
|
||||
|
||||
/* calculate distribution of work packages */
|
||||
uint32_t mod_pkgs = grid_size % this->wp_size;
|
||||
uint32_t mod_pkgs = grid.getTotalCellCount() % this->wp_size;
|
||||
uint32_t n_packages = (uint32_t)(grid.getTotalCellCount() / this->wp_size) +
|
||||
(mod_pkgs != 0 ? 1 : 0);
|
||||
|
||||
this->wp_sizes_vector =
|
||||
std::vector<uint32_t>(n_packages - mod_pkgs, this->wp_size);
|
||||
for (uint32_t i = 0; i < mod_pkgs; i++) {
|
||||
this->wp_sizes_vector.push_back(this->wp_size - 1);
|
||||
this->wp_sizes_vector = std::vector<uint32_t>(n_packages, this->wp_size);
|
||||
if (mod_pkgs) {
|
||||
auto itEndVector = this->wp_sizes_vector.end() - 1;
|
||||
for (uint32_t i = 0; i < this->wp_size - mod_pkgs; i++) {
|
||||
*(itEndVector - i) -= 1;
|
||||
}
|
||||
}
|
||||
|
||||
this->state = this->grid.registerState(
|
||||
poet::BaseChemModule::CHEMISTRY_MODULE_NAME, this->prop_names);
|
||||
std::vector<double> &field = this->state->mem;
|
||||
|
||||
field.resize(this->n_cells_per_prop * this->prop_names.size());
|
||||
for (uint32_t i = 0; i < this->prop_names.size(); i++) {
|
||||
std::vector<double> prop_vec =
|
||||
this->grid.getSpeciesByName(this->prop_names[i]);
|
||||
|
||||
std::copy(prop_vec.begin(), prop_vec.end(),
|
||||
field.begin() + (i * this->n_cells_per_prop));
|
||||
}
|
||||
}
|
||||
|
||||
@ -116,7 +131,6 @@ void ChemMaster::Simulate(double dt) {
|
||||
/* retrieve needed data from R runtime */
|
||||
iteration = (int)R.parseEval("mysetup$iter");
|
||||
// dt = (double)R.parseEval("mysetup$requested_dt");
|
||||
current_sim_time = (double)R.parseEval("mysetup$simulation_time") - dt;
|
||||
|
||||
/* setup local variables */
|
||||
pkg_to_send = wp_sizes_vector.size();
|
||||
@ -191,10 +205,12 @@ void ChemMaster::Simulate(double dt) {
|
||||
for (int i = 1; i < world_size; i++) {
|
||||
MPI_Send(NULL, 0, MPI_DOUBLE, i, TAG_DHT_ITER, MPI_COMM_WORLD);
|
||||
}
|
||||
|
||||
this->current_sim_time += dt;
|
||||
}
|
||||
|
||||
inline void ChemMaster::sendPkgs(int &pkg_to_send, int &count_pkgs,
|
||||
int &free_workers, double *work_pointer,
|
||||
int &free_workers, double *&work_pointer,
|
||||
const double &dt, const uint32_t iteration) {
|
||||
/* declare variables */
|
||||
double master_send_a, master_send_b;
|
||||
|
||||
@ -100,8 +100,8 @@ void ChemSeq::Simulate(double dt) {
|
||||
|
||||
// HACK: we will copy resulting field into global grid field. Maybe we will
|
||||
// find a more performant way here ...
|
||||
std::vector<double> vecSimResult =
|
||||
this->phreeqc_rm->GetWPFromInternals(this->n_cells_per_prop);
|
||||
std::vector<double> vecSimResult;
|
||||
this->phreeqc_rm->GetWPFromInternals(vecSimResult, this->n_cells_per_prop);
|
||||
std::copy(vecSimResult.begin(), vecSimResult.end(), field.begin());
|
||||
|
||||
R["TMP_C"] = field;
|
||||
|
||||
@ -18,15 +18,16 @@
|
||||
** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
||||
*/
|
||||
|
||||
#include "poet/ChemSimPar.hpp"
|
||||
#include "poet/SimParams.hpp"
|
||||
#include <Rcpp.h>
|
||||
|
||||
#include <algorithm>
|
||||
#include <cstdint>
|
||||
#include <iostream>
|
||||
#include <mpi.h>
|
||||
#include <ostream>
|
||||
#include <string>
|
||||
|
||||
#include <poet/ChemSimPar.hpp>
|
||||
#include <vector>
|
||||
|
||||
using namespace poet;
|
||||
@ -57,8 +58,8 @@ ChemWorker::ChemWorker(SimParams ¶ms, RInside &R_, Grid &grid_,
|
||||
|
||||
if (this->dht_enabled) {
|
||||
int data_size = this->prop_names.size() * sizeof(double);
|
||||
int key_size =
|
||||
this->prop_names.size() * sizeof(double) + (dt_differ * sizeof(double));
|
||||
int key_size = (this->prop_names.size() - 1) * sizeof(double) +
|
||||
(dt_differ * sizeof(double));
|
||||
int dht_buckets_per_process =
|
||||
dht_size_per_process / (1 + data_size + key_size);
|
||||
|
||||
@ -78,8 +79,6 @@ ChemWorker::ChemWorker(SimParams ¶ms, RInside &R_, Grid &grid_,
|
||||
readFile();
|
||||
// set size
|
||||
dht_flags.resize(wp_size, true);
|
||||
// assign all elements to true (default)
|
||||
dht_flags.assign(wp_size, true);
|
||||
}
|
||||
|
||||
this->timing.fill(0.0);
|
||||
@ -139,6 +138,8 @@ void ChemWorker::doWork(MPI_Status &probe_status) {
|
||||
|
||||
double dt;
|
||||
|
||||
bool bNoPhreeqc = false;
|
||||
|
||||
/* get number of doubles to be received */
|
||||
MPI_Get_count(&probe_status, MPI_DOUBLE, &count);
|
||||
|
||||
@ -153,125 +154,86 @@ void ChemWorker::doWork(MPI_Status &probe_status) {
|
||||
* mpi_buffer */
|
||||
|
||||
// work_package_size
|
||||
if (mpi_buffer[count] != local_work_package_size) { // work_package_size
|
||||
local_work_package_size = mpi_buffer[count];
|
||||
R["work_package_size"] = local_work_package_size;
|
||||
R.parseEvalQ("mysetup$work_package_size <- work_package_size");
|
||||
}
|
||||
local_work_package_size = mpi_buffer[count];
|
||||
|
||||
// current iteration of simulation
|
||||
this->iteration = mpi_buffer[count + 1];
|
||||
R["iter"] = this->iteration;
|
||||
R.parseEvalQ("mysetup$iter <- iter");
|
||||
|
||||
// current timestep size
|
||||
if (mpi_buffer[count + 2] != dt) {
|
||||
dt = mpi_buffer[count + 2];
|
||||
R["dt"] = dt;
|
||||
R.parseEvalQ("mysetup$dt <- dt");
|
||||
}
|
||||
dt = mpi_buffer[count + 2];
|
||||
|
||||
// current simulation time ('age' of simulation)
|
||||
if (mpi_buffer[count + 3] != current_sim_time) {
|
||||
current_sim_time = mpi_buffer[count + 3];
|
||||
R["simulation_time"] = current_sim_time;
|
||||
R.parseEvalQ("mysetup$simulation_time <- simulation_time");
|
||||
}
|
||||
current_sim_time = mpi_buffer[count + 3];
|
||||
|
||||
/* 4th double value is currently a placeholder */
|
||||
// if (mpi_buffer[count+4] != placeholder) {
|
||||
// placeholder = mpi_buffer[count+4];
|
||||
// R["mysetup$placeholder"] = placeholder;
|
||||
// }
|
||||
// placeholder = mpi_buffer[count+4];
|
||||
|
||||
std::vector<double> vecCurrWP(
|
||||
mpi_buffer,
|
||||
mpi_buffer + (local_work_package_size * this->prop_names.size()));
|
||||
|
||||
std::vector<int32_t> vecMappingWP(this->wp_size);
|
||||
std::vector<std::vector<double>> vecDHTResults;
|
||||
std::vector<double> vecDHTKeys;
|
||||
std::vector<bool> vecNeedPhreeqc(this->wp_size, true);
|
||||
|
||||
for (uint32_t i = 0; i < local_work_package_size; i++) {
|
||||
vecMappingWP[i] = i;
|
||||
}
|
||||
|
||||
if (local_work_package_size != this->wp_size) {
|
||||
std::vector<double> vecFiller(
|
||||
(this->wp_size - local_work_package_size) * this->prop_names.size(), 0);
|
||||
vecCurrWP.insert(vecCurrWP.end(), vecFiller.begin(), vecFiller.end());
|
||||
|
||||
// set all remaining cells to inactive
|
||||
for (int i = local_work_package_size; i < this->wp_size; i++) {
|
||||
vecMappingWP[i] = -1;
|
||||
}
|
||||
}
|
||||
|
||||
if (dht_enabled) {
|
||||
/* resize helper vector dht_flags of work_package_size changes */
|
||||
if ((int)dht_flags.size() != local_work_package_size) {
|
||||
dht_flags.resize(local_work_package_size, true); // set size
|
||||
dht_flags.assign(local_work_package_size,
|
||||
true); // assign all elements to true (default)
|
||||
}
|
||||
|
||||
/* check for values in DHT */
|
||||
dht_get_start = MPI_Wtime();
|
||||
dht->checkDHT(local_work_package_size, dht_flags, mpi_buffer, dt);
|
||||
vecDHTKeys = this->phreeqc_rm->ReplaceTotalsByPotentials(
|
||||
vecCurrWP, local_work_package_size);
|
||||
vecDHTResults = dht->checkDHT(local_work_package_size, vecNeedPhreeqc,
|
||||
vecDHTKeys.data(), dt);
|
||||
dht_get_end = MPI_Wtime();
|
||||
|
||||
/* distribute dht_flags to R Runtime */
|
||||
R["dht_flags"] = as<LogicalVector>(wrap(dht_flags));
|
||||
}
|
||||
|
||||
/* Convert grid to R runtime */
|
||||
size_t rowCount = local_work_package_size;
|
||||
size_t colCount = this->prop_names.size();
|
||||
|
||||
std::vector<std::vector<double>> input(colCount);
|
||||
|
||||
for (size_t i = 0; i < rowCount; i++) {
|
||||
for (size_t j = 0; j < colCount; j++) {
|
||||
input[j].push_back(mpi_buffer[i * colCount + j]);
|
||||
uint32_t iMappingIndex = 0;
|
||||
for (uint32_t i = 0; i < local_work_package_size; i++) {
|
||||
if (vecMappingWP[i] == -1) {
|
||||
continue;
|
||||
}
|
||||
vecMappingWP[i] = (vecNeedPhreeqc[i] ? iMappingIndex++ : -1);
|
||||
}
|
||||
}
|
||||
|
||||
R["work_package_full"] = Rcpp::as<Rcpp::DataFrame>(Rcpp::wrap(input));
|
||||
phreeqc_time_start = MPI_Wtime();
|
||||
this->phreeqc_rm->RunWorkPackage(vecCurrWP, vecMappingWP, current_sim_time,
|
||||
dt);
|
||||
phreeqc_time_end = MPI_Wtime();
|
||||
|
||||
if (dht_enabled) {
|
||||
R.parseEvalQ("work_package <- work_package_full[dht_flags,]");
|
||||
} else {
|
||||
R.parseEvalQ("work_package <- work_package_full");
|
||||
}
|
||||
|
||||
R.parseEvalQ("work_package <- as.matrix(work_package)");
|
||||
|
||||
unsigned int nrows = R.parseEval("nrow(work_package)");
|
||||
|
||||
if (nrows > 0) {
|
||||
/*Single Line error Workaround*/
|
||||
if (nrows <= 1) {
|
||||
// duplicate line to enable correct simmulation
|
||||
R.parseEvalQ("work_package <- work_package[rep(1:nrow(work_package), "
|
||||
"times = 2), ]");
|
||||
}
|
||||
/* Run PHREEQC */
|
||||
phreeqc_time_start = MPI_Wtime();
|
||||
R.parseEvalQ("result <- as.data.frame(slave_chemistry(setup=mysetup, "
|
||||
"data = work_package))");
|
||||
phreeqc_time_end = MPI_Wtime();
|
||||
} else {
|
||||
// undefined behaviour, isn't it?
|
||||
}
|
||||
|
||||
phreeqc_count++;
|
||||
|
||||
if (dht_enabled) {
|
||||
R.parseEvalQ("result_full <- work_package_full");
|
||||
if (nrows > 0)
|
||||
R.parseEvalQ("result_full[dht_flags,] <- result");
|
||||
} else {
|
||||
R.parseEvalQ("result_full <- result");
|
||||
}
|
||||
|
||||
/* convert grid to C domain */
|
||||
|
||||
std::vector<std::vector<double>> output =
|
||||
Rcpp::as<std::vector<std::vector<double>>>(R.parseEval("result_full"));
|
||||
|
||||
for (size_t i = 0; i < rowCount; i++) {
|
||||
for (size_t j = 0; j < colCount; j++) {
|
||||
mpi_buffer_results[i * colCount + j] = output[j][i];
|
||||
for (uint32_t i = 0; i < local_work_package_size; i++) {
|
||||
if (!vecNeedPhreeqc[i]) {
|
||||
std::copy(vecDHTResults[i].begin(), vecDHTResults[i].end(),
|
||||
vecCurrWP.begin() + (this->prop_names.size() * i));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// grid.exportWP(mpi_buffer_results);
|
||||
/* send results to master */
|
||||
MPI_Request send_req;
|
||||
MPI_Isend(mpi_buffer_results, count, MPI_DOUBLE, 0, TAG_WORK, MPI_COMM_WORLD,
|
||||
MPI_Isend(vecCurrWP.data(), count, MPI_DOUBLE, 0, TAG_WORK, MPI_COMM_WORLD,
|
||||
&send_req);
|
||||
|
||||
if (dht_enabled) {
|
||||
/* write results to DHT */
|
||||
dht_fill_start = MPI_Wtime();
|
||||
dht->fillDHT(local_work_package_size, dht_flags, mpi_buffer,
|
||||
mpi_buffer_results, dt);
|
||||
dht->fillDHT(local_work_package_size, vecNeedPhreeqc, vecDHTKeys.data(),
|
||||
vecCurrWP.data(), dt);
|
||||
dht_fill_end = MPI_Wtime();
|
||||
|
||||
timing[1] += dht_get_end - dht_get_start;
|
||||
|
||||
@ -20,6 +20,7 @@
|
||||
|
||||
#include "poet/HashFunctions.hpp"
|
||||
#include <cstddef>
|
||||
#include <cmath>
|
||||
#include <cstdint>
|
||||
#include <openssl/evp.h>
|
||||
#include <poet/DHT_Wrapper.hpp>
|
||||
@ -27,11 +28,16 @@
|
||||
#include <math.h>
|
||||
|
||||
#include <iostream>
|
||||
#include <stdexcept>
|
||||
#include <vector>
|
||||
|
||||
using namespace poet;
|
||||
using namespace std;
|
||||
|
||||
inline double round_signif(double value, int32_t signif) {
|
||||
const double multiplier = std::pow(10.0, signif);
|
||||
return .0 + std::trunc(value * multiplier) / multiplier;
|
||||
}
|
||||
|
||||
DHT_Wrapper::DHT_Wrapper(SimParams ¶ms, MPI_Comm dht_comm,
|
||||
int buckets_per_process, int data_size, int key_size) {
|
||||
poet::initHashCtx(EVP_md5());
|
||||
@ -59,40 +65,43 @@ DHT_Wrapper::~DHT_Wrapper() {
|
||||
poet::freeHashCtx();
|
||||
}
|
||||
|
||||
void DHT_Wrapper::checkDHT(int length, std::vector<bool> &out_result_index,
|
||||
double *work_package, double dt) {
|
||||
auto DHT_Wrapper::checkDHT(int length, std::vector<bool> &out_result_index,
|
||||
double *work_package, double dt)
|
||||
-> std::vector<std::vector<double>> {
|
||||
void *key;
|
||||
int res;
|
||||
// var count -> count of variables per grid cell
|
||||
int var_count = dht_prop_type_vector.size();
|
||||
int var_count = dht_prop_type_vector.size() - 1;
|
||||
std::vector<std::vector<double>> data(length);
|
||||
// loop over every grid cell contained in work package
|
||||
for (int i = 0; i < length; i++) {
|
||||
// point to current grid cell
|
||||
key = (void *)&(work_package[i * var_count]);
|
||||
data[i].resize(dht_prop_type_vector.size());
|
||||
|
||||
// fuzz data (round, logarithm etc.)
|
||||
fuzzForDHT(var_count, key, dt);
|
||||
std::vector<double> vecFuzz = fuzzForDHT(var_count, key, dt);
|
||||
|
||||
// overwrite input with data from DHT, IF value is found in DHT
|
||||
res = DHT_read(dht_object, fuzzing_buffer, key);
|
||||
res = DHT_read(dht_object, vecFuzz.data(), data[i].data());
|
||||
|
||||
// if DHT_SUCCESS value was found ...
|
||||
if (res == DHT_SUCCESS) {
|
||||
// ... and grid cell will be marked as 'not to be simulating'
|
||||
out_result_index[i] = false;
|
||||
dht_hits++;
|
||||
|
||||
}
|
||||
// ... otherwise ...
|
||||
else if (res == DHT_READ_MISS) {
|
||||
// grid cell needs to be simulated by PHREEQC
|
||||
out_result_index[i] = true;
|
||||
dht_miss++;
|
||||
} else {
|
||||
// MPI ERROR ... WHAT TO DO NOW?
|
||||
// RUNNING CIRCLES WHILE SCREAMING
|
||||
}
|
||||
}
|
||||
|
||||
return data;
|
||||
}
|
||||
|
||||
void DHT_Wrapper::fillDHT(int length, std::vector<bool> &result_index,
|
||||
@ -104,16 +113,16 @@ void DHT_Wrapper::fillDHT(int length, std::vector<bool> &result_index,
|
||||
int var_count = dht_prop_type_vector.size();
|
||||
// loop over every grid cell contained in work package
|
||||
for (int i = 0; i < length; i++) {
|
||||
key = (void *)&(work_package[i * var_count]);
|
||||
key = (void *)&(work_package[i * (var_count - 1)]);
|
||||
data = (void *)&(results[i * var_count]);
|
||||
|
||||
// If true grid cell was simulated, needs to be inserted into dht
|
||||
if (result_index[i]) {
|
||||
// fuzz data (round, logarithm etc.)
|
||||
fuzzForDHT(var_count, key, dt);
|
||||
std::vector<double> vecFuzz = fuzzForDHT(var_count - 1, key, dt);
|
||||
|
||||
// insert simulated data with fuzzed key into DHT
|
||||
res = DHT_write(dht_object, fuzzing_buffer, data);
|
||||
res = DHT_write(dht_object, vecFuzz.data(), data);
|
||||
|
||||
// if data was successfully written ...
|
||||
if (res != DHT_SUCCESS) {
|
||||
@ -164,7 +173,10 @@ uint64_t DHT_Wrapper::getMisses() { return this->dht_miss; }
|
||||
|
||||
uint64_t DHT_Wrapper::getEvictions() { return this->dht_evictions; }
|
||||
|
||||
void DHT_Wrapper::fuzzForDHT(int var_count, void *key, double dt) {
|
||||
std::vector<double> DHT_Wrapper::fuzzForDHT(int var_count, void *key,
|
||||
double dt) {
|
||||
|
||||
std::vector<double> vecFuzz(var_count, .0);
|
||||
unsigned int i = 0;
|
||||
// introduce fuzzing to allow more hits in DHT
|
||||
// loop over every variable of grid cell
|
||||
@ -179,34 +191,37 @@ void DHT_Wrapper::fuzzForDHT(int var_count, void *key, double dt) {
|
||||
cerr << "dht_wrapper.cpp::fuzz_for_dht(): Warning! Negative value in "
|
||||
"key!"
|
||||
<< endl;
|
||||
fuzzing_buffer[i] = 0;
|
||||
vecFuzz[i] = 0;
|
||||
}
|
||||
// if variable is 0 set fuzzing buffer to 0
|
||||
else if (((double *)key)[i] == 0)
|
||||
fuzzing_buffer[i] = 0;
|
||||
vecFuzz[i] = 0;
|
||||
// otherwise ...
|
||||
else
|
||||
else {
|
||||
// round current variable value by applying log with base 10, negate
|
||||
// (since the actual values will be between 0 and 1) and cut result
|
||||
// after significant digit
|
||||
fuzzing_buffer[i] =
|
||||
ROUND(-(std::log10(((double *)key)[i])), dht_signif_vector[i]);
|
||||
vecFuzz[i] = round_signif(-(std::log10(((double *)key)[i])),
|
||||
dht_signif_vector[i]);
|
||||
}
|
||||
}
|
||||
// if log is disabled
|
||||
else {
|
||||
// just round by cutting after signifanct digit
|
||||
fuzzing_buffer[i] = ROUND((((double *)key)[i]), dht_signif_vector[i]);
|
||||
vecFuzz[i] = round_signif((((double *)key)[i]), dht_signif_vector[i]);
|
||||
}
|
||||
} else if (dht_prop_type_vector[i] == "charge") {
|
||||
vecFuzz[i] = round_signif(((double *)key)[i], dht_signif_vector[i]);
|
||||
}
|
||||
// if variable is defined as 'logact' (log was already applied e.g. pH)
|
||||
else if (dht_prop_type_vector[i] == "logact") {
|
||||
// just round by cutting after signifanct digit
|
||||
fuzzing_buffer[i] = ROUND((((double *)key)[i]), dht_signif_vector[i]);
|
||||
vecFuzz[i] = round_signif((((double *)key)[i]), dht_signif_vector[i]);
|
||||
}
|
||||
// if defined ass 'ignore' ...
|
||||
else if (dht_prop_type_vector[i] == "ignore") {
|
||||
// ... just set fuzzing buffer to 0
|
||||
fuzzing_buffer[i] = 0;
|
||||
vecFuzz[i] = 0;
|
||||
}
|
||||
// and finally, if type is not defined, print error message
|
||||
else {
|
||||
@ -217,6 +232,9 @@ void DHT_Wrapper::fuzzForDHT(int var_count, void *key, double dt) {
|
||||
}
|
||||
// if timestep differs over iterations set current current time step at the
|
||||
// end of fuzzing buffer
|
||||
if (dt_differ)
|
||||
fuzzing_buffer[var_count] = dt;
|
||||
if (dt_differ) {
|
||||
vecFuzz[var_count] = dt;
|
||||
}
|
||||
|
||||
return vecFuzz;
|
||||
}
|
||||
|
||||
57
src/Grid.cpp
57
src/Grid.cpp
@ -18,12 +18,12 @@
|
||||
** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
||||
*/
|
||||
|
||||
#include "PhreeqcRM.h"
|
||||
#include "poet/SimParams.hpp"
|
||||
#include <Rcpp.h>
|
||||
#include <algorithm>
|
||||
#include <bits/stdint-intn.h>
|
||||
#include <bits/stdint-uintn.h>
|
||||
#include <cassert>
|
||||
#include <cstdint>
|
||||
#include <new>
|
||||
#include <numeric>
|
||||
#include <poet/Grid.hpp>
|
||||
@ -99,6 +99,9 @@ Grid::~Grid() {
|
||||
for (auto &[key, val] : this->state_register) {
|
||||
delete val;
|
||||
}
|
||||
if (this->phreeqc_rm != nullptr) {
|
||||
delete this->phreeqc_rm;
|
||||
}
|
||||
}
|
||||
|
||||
poet::StateMemory *Grid::registerState(std::string module_name,
|
||||
@ -204,3 +207,53 @@ auto poet::Grid::getSpeciesByName(std::string name, std::string module_name)
|
||||
return std::vector<double>(module_memory->mem.begin() + begin_vec,
|
||||
module_memory->mem.begin() + end_vec);
|
||||
}
|
||||
|
||||
PhreeqcRM *poet::Grid::initPhreeqc(const poet::GridParams ¶ms) {
|
||||
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;
|
||||
}
|
||||
|
||||
@ -1,13 +1,18 @@
|
||||
#include "poet/PhreeqcWrapper.hpp"
|
||||
#include "IPhreeqc.hpp"
|
||||
#include "IrmResult.h"
|
||||
#include "PhreeqcRM.h"
|
||||
|
||||
#include <algorithm>
|
||||
#include <bits/stdint-uintn.h>
|
||||
#include <cstdint>
|
||||
#include <iterator>
|
||||
#include <poet/PhreeqcWrapper.hpp>
|
||||
#include <stdexcept>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
poet::PhreeqcWrapper::PhreeqcWrapper(uint32_t inxyz)
|
||||
: PhreeqcRM(inxyz, 1), iWPSize(inxyz) {
|
||||
this->vecDefMapping.reserve(inxyz);
|
||||
this->vecDefMapping.resize(inxyz);
|
||||
|
||||
for (uint32_t i = 0; i < inxyz; i++) {
|
||||
this->vecDefMapping[i] = i;
|
||||
@ -60,7 +65,7 @@ void poet::PhreeqcWrapper::SetupAndLoadDB(
|
||||
}
|
||||
|
||||
void poet::PhreeqcWrapper::InitFromFile(const std::string &strInputFile) {
|
||||
this->RunFile(true, true, true, strInputFile);
|
||||
this->RunFile(true, true, false, strInputFile);
|
||||
this->RunString(true, false, true, "DELETE; -all");
|
||||
|
||||
this->FindComponents();
|
||||
@ -81,11 +86,11 @@ void poet::PhreeqcWrapper::InitFromFile(const std::string &strInputFile) {
|
||||
this->InitialPhreeqc2Module(ic1);
|
||||
|
||||
// Initial equilibration of cells
|
||||
double time = 0.0;
|
||||
double time_step = 0.0;
|
||||
this->SetTime(time);
|
||||
this->SetTimeStep(time_step);
|
||||
this->RunCells();
|
||||
// double time = 0.0;
|
||||
// double time_step = 0.0;
|
||||
// this->SetTime(time);
|
||||
// this->SetTimeStep(time_step);
|
||||
// this->RunCells();
|
||||
|
||||
this->InitInternals();
|
||||
}
|
||||
@ -142,22 +147,25 @@ void poet::PhreeqcWrapper::SetInternalsFromWP(const std::vector<double> &vecWP,
|
||||
uint32_t iCurrWPSize) {
|
||||
uint32_t iCurrElements;
|
||||
|
||||
if (iCurrWPSize != this->iWPSize) {
|
||||
this->SetMappingForWP(iCurrWPSize);
|
||||
}
|
||||
|
||||
auto itStart = vecWP.begin();
|
||||
auto itEnd = itStart;
|
||||
|
||||
// this->SetMappingForWP(iCurrWPSize);
|
||||
|
||||
int nchem = this->GetChemistryCellCount();
|
||||
|
||||
iCurrElements = this->vecSpeciesPerModule[0];
|
||||
|
||||
itEnd += iCurrElements * iCurrWPSize;
|
||||
itEnd += iCurrElements * this->iWPSize;
|
||||
std::vector<double> out;
|
||||
this->GetConcentrations(out);
|
||||
this->SetConcentrations(std::vector<double>(itStart, itEnd));
|
||||
itStart = itEnd;
|
||||
|
||||
// Equlibirum Phases
|
||||
if ((iCurrElements = this->vecSpeciesPerModule[1]) != 0) {
|
||||
itEnd += iCurrElements * iCurrWPSize;
|
||||
itEnd += iCurrElements * this->iWPSize;
|
||||
this->GetPPhaseMoles(out);
|
||||
this->SetPPhaseMoles(std::vector<double>(itStart, itEnd));
|
||||
itStart = itEnd;
|
||||
}
|
||||
@ -168,24 +176,25 @@ void poet::PhreeqcWrapper::SetInternalsFromWP(const std::vector<double> &vecWP,
|
||||
// // BLOCK_END
|
||||
|
||||
if ((iCurrElements = this->vecSpeciesPerModule[4]) != 0) {
|
||||
itEnd += iCurrElements * iCurrWPSize;
|
||||
itEnd += iCurrElements * this->iWPSize;
|
||||
this->GetKineticsMoles(out);
|
||||
this->SetKineticsMoles(std::vector<double>(itStart, itEnd));
|
||||
itStart = itEnd;
|
||||
}
|
||||
}
|
||||
auto poet::PhreeqcWrapper::GetWPFromInternals(uint32_t iCurrWPSize)
|
||||
-> std::vector<double> {
|
||||
std::vector<double> vecOutput, vecCurrOutput;
|
||||
void poet::PhreeqcWrapper::GetWPFromInternals(std::vector<double> &vecWP,
|
||||
uint32_t iCurrWPSize) {
|
||||
std::vector<double> vecCurrOutput;
|
||||
|
||||
vecOutput.reserve(this->iSpeciesCount * iCurrWPSize);
|
||||
vecWP.clear();
|
||||
vecWP.reserve(this->iSpeciesCount * this->iWPSize);
|
||||
|
||||
this->GetConcentrations(vecCurrOutput);
|
||||
vecOutput.insert(vecOutput.end(), vecCurrOutput.begin(), vecCurrOutput.end());
|
||||
vecWP.insert(vecWP.end(), vecCurrOutput.begin(), vecCurrOutput.end());
|
||||
|
||||
if (this->vecSpeciesPerModule[1] != 0) {
|
||||
this->GetPPhaseMoles(vecCurrOutput);
|
||||
vecOutput.insert(vecOutput.end(), vecCurrOutput.begin(),
|
||||
vecCurrOutput.end());
|
||||
vecWP.insert(vecWP.end(), vecCurrOutput.begin(), vecCurrOutput.end());
|
||||
}
|
||||
|
||||
// NOTE: Block for 'Surface' and 'Exchange' is missing because of missing
|
||||
@ -195,21 +204,132 @@ auto poet::PhreeqcWrapper::GetWPFromInternals(uint32_t iCurrWPSize)
|
||||
|
||||
if (this->vecSpeciesPerModule[4] != 0) {
|
||||
this->GetKineticsMoles(vecCurrOutput);
|
||||
vecOutput.insert(vecOutput.end(), vecCurrOutput.begin(),
|
||||
vecCurrOutput.end());
|
||||
vecWP.insert(vecWP.end(), vecCurrOutput.begin(), vecCurrOutput.end());
|
||||
}
|
||||
}
|
||||
|
||||
auto poet::PhreeqcWrapper::ReplaceTotalsByPotentials(
|
||||
const std::vector<double> &vecWP, uint32_t iCurrWPSize)
|
||||
-> std::vector<double> {
|
||||
|
||||
uint32_t iPropsPerCell = this->vecSpeciesNames.size();
|
||||
|
||||
int iphreeqcResult;
|
||||
std::vector<double> vecOutput;
|
||||
vecOutput.reserve((iPropsPerCell - 1) * iCurrWPSize);
|
||||
|
||||
auto itCStart = vecWP.begin();
|
||||
auto itCEnd = itCStart + (this->vecSpeciesPerModule[0]);
|
||||
uint32_t iDiff = iPropsPerCell - this->vecSpeciesPerModule[0];
|
||||
|
||||
for (uint32_t i = 0; i < iCurrWPSize; i++) {
|
||||
std::vector<double> vecCIn(itCStart, itCEnd);
|
||||
double pH, pe;
|
||||
|
||||
// FIXME: Hardcoded temperatures and pressures here!
|
||||
IPhreeqc *util_ptr = this->Concentrations2Utility(
|
||||
vecCIn, std::vector<double>(1, 25.0), std::vector<double>(1, 1.0));
|
||||
|
||||
std::string sInput = "SELECTED_OUTPUT " + std::to_string(this->DHT_SELOUT) +
|
||||
"; -pH; -pe;RUN_CELLS; -cells 1; -time_step 0";
|
||||
iphreeqcResult = util_ptr->RunString(sInput.c_str());
|
||||
if (iphreeqcResult != 0) {
|
||||
throw std::runtime_error("IPhreeqc Utility returned non-zero value: " +
|
||||
std::to_string(iphreeqcResult));
|
||||
}
|
||||
|
||||
util_ptr->SetCurrentSelectedOutputUserNumber(this->DHT_SELOUT);
|
||||
|
||||
int vtype;
|
||||
static std::string svalue(100, '\0');
|
||||
|
||||
iphreeqcResult = util_ptr->GetSelectedOutputValue2(
|
||||
1, 0, &vtype, &pH, svalue.data(), svalue.capacity());
|
||||
if (iphreeqcResult != 0) {
|
||||
throw std::runtime_error("IPhreeqc Utility returned non-zero value: " +
|
||||
std::to_string(iphreeqcResult));
|
||||
}
|
||||
|
||||
iphreeqcResult = util_ptr->GetSelectedOutputValue2(
|
||||
1, 1, &vtype, &pe, svalue.data(), svalue.capacity());
|
||||
if (iphreeqcResult != 0) {
|
||||
throw std::runtime_error("IPhreeqc Utility returned non-zero value: " +
|
||||
std::to_string(iphreeqcResult));
|
||||
}
|
||||
|
||||
vecOutput.insert(vecOutput.end(), vecCIn.begin() + 3, vecCIn.end());
|
||||
vecOutput.push_back(pH);
|
||||
vecOutput.push_back(pe);
|
||||
vecOutput.insert(vecOutput.end(), itCEnd, itCEnd + iDiff);
|
||||
|
||||
itCStart = itCEnd + iDiff;
|
||||
itCEnd = itCStart + (this->vecSpeciesPerModule[0]);
|
||||
}
|
||||
|
||||
return vecOutput;
|
||||
}
|
||||
|
||||
void poet::PhreeqcWrapper::SetMappingForWP(uint32_t iCurrWPSize) {
|
||||
IRM_RESULT
|
||||
poet::PhreeqcWrapper::RunWorkPackage(std::vector<double> &vecWP,
|
||||
std::vector<int32_t> &vecMapping,
|
||||
double dSimTime, double dTimestep) {
|
||||
if (this->iWPSize != vecMapping.size()) {
|
||||
return IRM_INVALIDARG;
|
||||
}
|
||||
|
||||
if ((this->iWPSize * this->iSpeciesCount) != vecWP.size()) {
|
||||
return IRM_INVALIDARG;
|
||||
}
|
||||
|
||||
// check if we actually need to start phreeqc
|
||||
bool bRunPhreeqc = false;
|
||||
for (const auto &aMappingNum : vecMapping) {
|
||||
if (aMappingNum != -1) {
|
||||
bRunPhreeqc = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if (!bRunPhreeqc) {
|
||||
return IRM_OK;
|
||||
}
|
||||
|
||||
std::vector<double> vecCopy;
|
||||
|
||||
vecCopy = vecWP;
|
||||
for (uint32_t i = 0; i < this->iSpeciesCount; i++) {
|
||||
for (uint32_t j = 0; j < this->iWPSize; j++) {
|
||||
vecWP[(i * this->iWPSize) + j] = vecCopy[(j * this->iSpeciesCount) + i];
|
||||
}
|
||||
}
|
||||
|
||||
IRM_RESULT result;
|
||||
this->CreateMapping(vecMapping);
|
||||
this->SetInternalsFromWP(vecWP, this->iWPSize);
|
||||
this->SetTime(dSimTime);
|
||||
this->SetTimeStep(dTimestep);
|
||||
result = this->PhreeqcRM::RunCells();
|
||||
this->GetWPFromInternals(vecWP, this->iWPSize);
|
||||
|
||||
vecCopy = vecWP;
|
||||
for (uint32_t i = 0; i < this->iSpeciesCount; i++) {
|
||||
for (uint32_t j = 0; j < this->iWPSize; j++) {
|
||||
vecWP[(j * this->iSpeciesCount) + i] = vecCopy[(i * this->iWPSize) + j];
|
||||
}
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
inline void poet::PhreeqcWrapper::SetMappingForWP(uint32_t iCurrWPSize) {
|
||||
if (iCurrWPSize == this->iWPSize) {
|
||||
this->CreateMapping(this->vecDefMapping);
|
||||
} else {
|
||||
std::vector<int> vecCurrMapping = this->vecDefMapping;
|
||||
for (uint32_t i = iCurrWPSize; i < vecCurrMapping.size(); i++) {
|
||||
vecCurrMapping[i] = -1;
|
||||
}
|
||||
this->CreateMapping(vecCurrMapping);
|
||||
return;
|
||||
}
|
||||
|
||||
std::vector<int> vecCurrMapping(this->vecDefMapping);
|
||||
for (uint32_t i = iCurrWPSize; i < vecCurrMapping.size(); i++) {
|
||||
vecCurrMapping[i] = -1;
|
||||
}
|
||||
this->CreateMapping(vecCurrMapping);
|
||||
}
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user