refactor: make Grid class a bit more generic

This commit is contained in:
Max Luebke 2023-01-23 14:43:25 +01:00
parent fc20d7a7f0
commit 689c55bf32
9 changed files with 220 additions and 187 deletions

View File

@ -179,11 +179,11 @@ int main(int argc, char *argv[]) {
// TODO: Grid anpassen
GridParams grid_params = GridParams(R);
Grid grid(grid_params);
// grid.init_from_R();
Grid grid;
params.initVectorParams(R, grid.getSpeciesCount());
grid.InitModuleFromParams(GridParams(R));
params.initVectorParams(R, grid.GetSpeciesCount());
// MDL: store all parameters
if (world_rank == 0) {

@ -1 +1 @@
Subproject commit a9e3d9fa355160ef7efdb5e421076e72c0e58893
Subproject commit e7769e1a718c9993c1500c5df13dc64b39362b9a

@ -1 +1 @@
Subproject commit 0abd4e04ae40724f9ddb62da9f7c6ebd234330bb
Subproject commit d4e3ab8544628c72921bf60d4462d6d131b0ffb5

View File

@ -2,7 +2,7 @@
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam)
** Copyright (C) 2018-2023 Marco De Lucia, Max Luebke (GFZ Potsdam)
**
** POET is free software; you can redistribute it and/or modify it under the
** terms of the GNU General Public License as published by the Free Software
@ -21,10 +21,10 @@
#ifndef GRID_H
#define GRID_H
#include "PhreeqcRM.h"
#include "poet/SimParams.hpp"
#include <Rcpp.h>
#include <array>
#include <cstddef>
#include <cstdint>
#include <map>
#include <string>
@ -53,65 +53,54 @@ constexpr const char *GRID_MODULE_NAME = "grid_init";
*
*/
class Grid {
private:
using prop_vec = std::vector<std::string>;
public:
/**
* @brief Construct a new Grid object
*
* This will call the default constructor and initializes private RRuntime
* with given R runtime.
*
* @param R
*/
Grid(poet::GridParams grid_args);
Grid();
~Grid();
/**
* @brief Init the grid
*
* At this moment init will only declare and define a variable inside the R
* runtime called grid_tmp since the whole Grid initialization and management
* is done by the R runtime. This may change in the future.
*
*/
void init_from_R();
auto getGridDimension() -> uint8_t;
auto getTotalCellCount() -> uint32_t;
auto getGridCellsCount(uint8_t direction) -> uint32_t;
auto getDomainSize(uint8_t dircetion) -> uint32_t;
void InitModuleFromParams(const poet::GridParams &grid_args);
StateMemory *registerState(std::string module_name,
std::vector<std::string> props);
StateMemory *getStatePointer(std::string module_name);
void deregisterState(std::string module_name);
void SetGridDimension(uint8_t dim);
void SetGridCellCount(uint32_t n_x, uint32_t n_y = 0, uint32_t n_z = 0);
void SetGridSize(double s_x, double s_y = 0., double s_z = 0.);
void SetPropNames(const std::vector<std::string> &prop_names);
auto getSpeciesCount() -> uint32_t;
auto getPropNames() -> prop_vec;
void InitGridFromScratch(const Rcpp::DataFrame &init_cell);
void InitGridFromRDS();
void InitGridFromPhreeqc();
auto getSpeciesByName(std::string name,
std::string module_name = poet::GRID_MODULE_NAME)
auto GetGridDimension() const -> uint8_t;
auto GetTotalCellCount() const -> uint32_t;
auto GetGridCellsCount(uint8_t direction) const -> uint32_t;
auto GetGridSize(uint8_t direction) const -> uint32_t;
auto RegisterState(std::string module_name, std::vector<std::string> props)
-> StateMemory *;
auto GetStatePointer(std::string module_name) -> StateMemory *;
auto GetInitialGrid() const -> StateMemory *;
auto GetSpeciesCount() const -> uint32_t;
auto GetPropNames() const -> std::vector<std::string>;
auto GetSpeciesByName(std::string name,
std::string module_name = poet::GRID_MODULE_NAME) const
-> std::vector<double>;
private:
PhreeqcRM *initPhreeqc(const poet::GridParams &params);
// auto InitPhreeqc(const poet::GridParams &params) -> PhreeqcRM *;
std::uint8_t dim;
std::array<std::uint32_t, MAX_DIM> d_spatial;
std::uint8_t dim = 0;
std::array<double, MAX_DIM> grid_size;
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;
StateMemory *grid_init = std::nullptr_t();
prop_vec prop_names;
std::vector<std::string> prop_names;
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);
// PhreeqcRM *phreeqc_rm = std::nullptr_t();
};
} // namespace poet
#endif // GRID_H

View File

@ -6,7 +6,7 @@ poet::BaseChemModule::BaseChemModule(SimParams &params, RInside &R_,
t_simparams tmp = params.getNumParams();
this->world_rank = tmp.world_rank;
this->world_size = tmp.world_size;
this->prop_names = this->grid.getPropNames();
this->prop_names = this->grid.GetPropNames();
this->n_cells_per_prop = this->grid.getTotalCellCount();
this->n_cells_per_prop = this->grid.GetTotalCellCount();
}

View File

@ -42,7 +42,7 @@ ChemMaster::ChemMaster(SimParams &params, RInside &R_, Grid &grid_)
this->wp_size = tmp.wp_size;
this->dht_enabled = tmp.dht_enabled;
uint32_t grid_size = grid.getTotalCellCount() * this->prop_names.size();
uint32_t grid_size = grid.GetTotalCellCount() * this->prop_names.size();
/* allocate memory */
this->workerlist = new worker_struct[this->world_size - 1];
@ -52,8 +52,8 @@ ChemMaster::ChemMaster(SimParams &params, RInside &R_, Grid &grid_)
this->mpi_buffer = new double[grid_size];
/* calculate distribution of work packages */
uint32_t mod_pkgs = grid.getTotalCellCount() % this->wp_size;
uint32_t n_packages = (uint32_t)(grid.getTotalCellCount() / 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, this->wp_size);
@ -64,14 +64,14 @@ ChemMaster::ChemMaster(SimParams &params, RInside &R_, Grid &grid_)
}
}
this->state = this->grid.registerState(
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]);
this->grid.GetSpeciesByName(this->prop_names[i]);
std::copy(prop_vec.begin(), prop_vec.end(),
field.begin() + (i * this->n_cells_per_prop));
@ -105,7 +105,7 @@ void ChemMaster::Simulate(double dt) {
for (uint32_t i = 0; i < this->prop_names.size(); i++) {
try {
std::vector<double> t_prop_vec = this->grid.getSpeciesByName(
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(),
@ -233,7 +233,7 @@ inline void ChemMaster::sendPkgs(int &pkg_to_send, int &count_pkgs,
workerlist[p].send_addr = work_pointer;
/* push work pointer to next work package */
end_of_wp = local_work_package_size * grid.getSpeciesCount();
end_of_wp = local_work_package_size * grid.GetSpeciesCount();
work_pointer = &(work_pointer[end_of_wp]);
// fill send buffer starting with work_package ...

View File

@ -37,14 +37,14 @@ using namespace poet;
ChemSeq::ChemSeq(SimParams &params, RInside &R_, Grid &grid_)
: BaseChemModule(params, R_, grid_) {
this->state = this->grid.registerState(
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]);
this->grid.GetSpeciesByName(this->prop_names[i]);
std::copy(prop_vec.begin(), prop_vec.end(),
field.begin() + (i * this->n_cells_per_prop));
@ -52,7 +52,6 @@ ChemSeq::ChemSeq(SimParams &params, RInside &R_, Grid &grid_)
}
poet::ChemSeq::~ChemSeq() {
this->grid.deregisterState(poet::BaseChemModule::CHEMISTRY_MODULE_NAME);
if (this->phreeqc_rm) {
delete this->phreeqc_rm;
}
@ -75,7 +74,7 @@ void ChemSeq::Simulate(double dt) {
for (uint32_t i = 0; i < this->prop_names.size(); i++) {
try {
std::vector<double> t_prop_vec = this->grid.getSpeciesByName(
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(),

View File

@ -24,13 +24,11 @@
#include <Rcpp.h>
#include <algorithm>
#include <cstdint>
#include <ostream>
#include <poet/ChemSimSeq.hpp>
#include <poet/DiffusionModule.hpp>
#include <poet/Grid.hpp>
#include <array>
#include <bits/stdint-uintn.h>
#include <cassert>
#include <mpi.h>
#include <string>
@ -58,14 +56,14 @@ inline const char *convert_bc_to_config_file(uint8_t in) {
DiffusionModule::DiffusionModule(poet::DiffusionParams diffu_args, Grid &grid_)
: grid(grid_) {
this->diff_input.setGridCellN(grid_.getGridCellsCount(GRID_X_DIR),
grid_.getGridCellsCount(GRID_Y_DIR));
this->diff_input.setDomainSize(grid_.getDomainSize(GRID_X_DIR),
grid_.getDomainSize(GRID_Y_DIR));
this->diff_input.setGridCellN(grid_.GetGridCellsCount(GRID_X_DIR),
grid_.GetGridCellsCount(GRID_Y_DIR));
this->diff_input.setDomainSize(grid_.GetGridSize(GRID_X_DIR),
grid_.GetGridSize(GRID_Y_DIR));
this->dim = grid_.getGridDimension();
this->dim = grid_.GetGridDimension();
this->n_cells_per_prop = grid_.getTotalCellCount();
this->n_cells_per_prop = grid_.GetTotalCellCount();
this->initialize(diffu_args);
}
@ -78,7 +76,7 @@ void DiffusionModule::initialize(poet::DiffusionParams args) {
this->prop_count = this->prop_names.size();
this->state =
this->grid.registerState(DIFFUSION_MODULE_NAME, this->prop_names);
this->grid.RegisterState(DIFFUSION_MODULE_NAME, this->prop_names);
std::vector<double> &field = this->state->mem;
@ -91,15 +89,15 @@ void DiffusionModule::initialize(poet::DiffusionParams args) {
// initial input
this->alpha.push_back(args.alpha[this->prop_names[i]]);
std::vector<double> prop_vec = grid.getSpeciesByName(this->prop_names[i]);
std::vector<double> prop_vec = grid.GetSpeciesByName(this->prop_names[i]);
std::copy(prop_vec.begin(), prop_vec.end(),
field.begin() + (i * this->n_cells_per_prop));
if (this->dim == this->DIM_2D) {
tug::bc::BoundaryCondition bc(this->grid.getGridCellsCount(GRID_X_DIR),
this->grid.getGridCellsCount(GRID_Y_DIR));
tug::bc::BoundaryCondition bc(this->grid.GetGridCellsCount(GRID_X_DIR),
this->grid.GetGridCellsCount(GRID_Y_DIR));
this->bc_vec.push_back(bc);
} else {
tug::bc::BoundaryCondition bc(this->grid.getGridCellsCount(GRID_X_DIR));
tug::bc::BoundaryCondition bc(this->grid.GetGridCellsCount(GRID_X_DIR));
this->bc_vec.push_back(bc);
}
}
@ -173,7 +171,7 @@ void DiffusionModule::simulate(double dt) {
// 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(
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(),

View File

@ -2,7 +2,7 @@
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam)
** Copyright (C) 2018-2023 Marco De Lucia, Max Luebke (GFZ Potsdam)
**
** POET is free software; you can redistribute it and/or modify it under the
** terms of the GNU General Public License as published by the Free Software
@ -18,7 +18,6 @@
** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include "PhreeqcRM.h"
#include "poet/SimParams.hpp"
#include <Rcpp.h>
#include <algorithm>
@ -49,62 +48,100 @@ static inline int8_t get_type_id(std::string type) {
return -1;
}
Grid::Grid(poet::GridParams grid_args) {
// get dimension of grid
this->dim = grid_args.n_cells.size();
void poet::Grid::InitModuleFromParams(const poet::GridParams &grid_args) {
assert(grid_args.n_cells.size() == grid_args.s_cells.size());
// assert that dim is 1 or 2
assert(this->dim <= MAX_DIM && this->dim > 0);
// also assert that vetor spatial discretization is equal in size of
// dimensions
assert(this->dim == grid_args.s_cells.size());
this->SetGridDimension(grid_args.n_cells.size());
this->n_cells.fill(0);
this->d_spatial.fill(0);
// store grid dimensions
std::copy(grid_args.n_cells.begin(), grid_args.n_cells.end(),
this->n_cells.begin());
std::copy(grid_args.s_cells.begin(), grid_args.s_cells.end(),
this->d_spatial.begin());
this->grid_size.begin());
auto init_type = get_type_id(grid_args.type);
assert(init_type >= 0);
uint32_t vec_size = this->getTotalCellCount();
if (init_type != INIT_PHREEQC) {
this->prop_names = grid_args.props;
switch (init_type) {
case INIT_SCRATCH:
this->SetPropNames(grid_args.props);
this->InitGridFromScratch(grid_args.init_df);
break;
case INIT_RDS:
this->SetPropNames(grid_args.props);
this->InitGridFromRDS();
break;
case INIT_PHREEQC:
this->InitGridFromPhreeqc();
}
}
for (auto const &prop : this->prop_names) {
std::vector<double> prop_vec(vec_size);
// size of the vector shall be 1
if (init_type == INIT_SCRATCH) {
double prop_val =
((Rcpp::NumericVector)grid_args.init_df[prop.c_str()])[0];
std::fill(prop_vec.begin(), prop_vec.end(), prop_val);
} else if (init_type == INIT_RDS) {
prop_vec = Rcpp::as<std::vector<double>>(grid_args.init_df[prop.c_str()]);
}
// TODO: define behaviour for phreeqc script input
void poet::Grid::SetGridDimension(uint8_t dim) {
assert(dim < MAX_DIM + 1);
// NOTE: if there are props defined more than once the first written value
// will be used
this->grid_init.insert(std::pair(prop, prop_vec));
}
this->dim = dim;
}
void poet::Grid::SetGridCellCount(uint32_t n_x, uint32_t n_y, uint32_t n_z) {
assert(this->dim < MAX_DIM + 1);
this->n_cells[0] = n_x;
this->n_cells[1] = n_y;
// this->n_cells[2] = n_z;
}
void poet::Grid::SetGridSize(double s_x, double s_y, double s_z) {
assert(this->dim < MAX_DIM + 1);
this->grid_size[0] = s_x;
this->grid_size[1] = s_y;
// this->grid_size[2] = s_z;
}
void poet::Grid::SetPropNames(const std::vector<std::string> &prop_names) {
this->prop_names = prop_names;
}
Grid::Grid() {
this->n_cells.fill(0);
this->grid_size.fill(0);
};
Grid::~Grid() {
for (auto &[key, val] : this->state_register) {
delete val;
}
if (this->phreeqc_rm != nullptr) {
delete this->phreeqc_rm;
}
// if (this->phreeqc_rm != nullptr) {
// delete this->phreeqc_rm;
// }
}
poet::StateMemory *Grid::registerState(std::string module_name,
void poet::Grid::InitGridFromScratch(const Rcpp::DataFrame &init_cell) {
const uint32_t vec_size = this->GetTotalCellCount();
StateMemory *curr_state =
this->RegisterState(poet::GRID_MODULE_NAME, this->prop_names);
curr_state->props = this->prop_names;
std::vector<double> &curr_field = curr_state->mem;
for (auto const &prop : this->prop_names) {
// size of the vector shall be 1
double prop_val = ((Rcpp::NumericVector)init_cell[prop.c_str()])[0];
std::vector<double> prop_vec(vec_size, prop_val);
curr_field.insert(curr_field.end(), prop_vec.begin(), prop_vec.end());
}
this->grid_init = curr_state;
}
void poet::Grid::InitGridFromRDS() {
// TODO
}
void poet::Grid::InitGridFromPhreeqc() {
// TODO
}
poet::StateMemory *Grid::RegisterState(std::string module_name,
std::vector<std::string> props) {
poet::StateMemory *mem = new poet::StateMemory;
@ -120,7 +157,7 @@ poet::StateMemory *Grid::registerState(std::string module_name,
return mem;
}
poet::StateMemory *Grid::getStatePointer(std::string module_name) {
poet::StateMemory *Grid::GetStatePointer(std::string module_name) {
const auto it = this->state_register.find(module_name);
if (it == this->state_register.end()) {
@ -131,27 +168,27 @@ 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);
// 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"));
}
// 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;
// auto *p = it->second;
// delete p;
this->state_register.erase(it);
// this->state_register.erase(it);
// }
auto Grid::GetGridSize(uint8_t direction) const -> uint32_t {
return this->grid_size.at(direction);
}
auto Grid::getDomainSize(uint8_t dircetion) -> uint32_t {
return this->d_spatial.at(dircetion);
}
auto Grid::GetGridDimension() const -> uint8_t { return this->dim; }
auto Grid::getGridDimension() -> uint8_t { return this->dim; }
auto Grid::getTotalCellCount() -> uint32_t {
auto Grid::GetTotalCellCount() const -> uint32_t {
uint32_t sum = 1;
for (auto const &dim : this->n_cells) {
sum *= (dim != 0 ? dim : 1);
@ -160,28 +197,38 @@ auto Grid::getTotalCellCount() -> uint32_t {
return sum;
}
auto Grid::getGridCellsCount(uint8_t direction) -> uint32_t {
auto Grid::GetGridCellsCount(uint8_t direction) const -> uint32_t {
return this->n_cells.at(direction);
}
auto Grid::getSpeciesCount() -> uint32_t { return this->prop_names.size(); }
auto poet::Grid::GetInitialGrid() const -> StateMemory * {
return this->grid_init;
}
auto Grid::getPropNames() -> prop_vec { return this->prop_names; }
auto Grid::GetSpeciesCount() const -> uint32_t {
return this->prop_names.size();
}
auto poet::Grid::getSpeciesByName(std::string name, std::string module_name)
auto Grid::GetPropNames() const -> std::vector<std::string> {
return this->prop_names;
}
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 (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"));
}
// if (it == this->grid_init.end()) {
// throw std::range_error(
// std::string("Species " + name + " was not found for initial
// grid"));
// }
return it->second;
}
// return it->second;
// }
// else find module with name
auto const it = this->state_register.find(module_name);
@ -201,59 +248,59 @@ auto poet::Grid::getSpeciesByName(std::string name, std::string module_name)
}
uint32_t prop_index = find_res - prop_vector.begin();
uint32_t begin_vec = prop_index * this->getTotalCellCount(),
end_vec = ((prop_index + 1) * this->getTotalCellCount());
uint32_t begin_vec = prop_index * this->GetTotalCellCount(),
end_vec = ((prop_index + 1) * this->GetTotalCellCount());
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);
// 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);
// // 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 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 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 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);
// // 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");
// // 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;
}
// return chem_model;
// }