max/poet
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/poet.git synced 2025-12-15 12:28:22 +01:00
Code Issues Packages Projects Releases Wiki Activity

feat: Add PinnableMemory class for allocating pinned memory

Browse Source 
refactor: Field class now relies on Pinnable memory
This commit is contained in:
Max Lübke 2024-05-07 08:06:17 +00:00
parent 3da5b8eb20
commit d5744191a1
17 changed files with 445 additions and 633 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/CMakeLists.txt
View File

@ -6,7 +6,7 @@ target_sources(POETLib
Init/GridInit.cpp
Init/DiffusionInit.cpp
Init/ChemistryInit.cpp
DataStructures/Field.cpp
# DataStructures/Field.cpp
Transport/DiffusionModule.cpp
Chemistry/ChemistryModule.cpp
Chemistry/MasterFunctions.cpp

38
src/Chemistry/ChemistryModule.cpp
View File

@ -334,40 +334,4 @@ void poet::ChemistryModule::initializeInterp(
interp->setInterpolationFunction(inverseDistanceWeighting);
}
}
std::vector<double>
poet::ChemistryModule::shuffleField(const std::vector<double> &in_field,
uint32_t size_per_prop, uint32_t prop_count,
uint32_t wp_count) {
std::vector<double> out_buffer(in_field.size());
uint32_t write_i = 0;
for (uint32_t i = 0; i < wp_count; i++) {
for (uint32_t j = i; j < size_per_prop; j += wp_count) {
for (uint32_t k = 0; k < prop_count; k++) {
out_buffer[(write_i * prop_count) + k] =
in_field[(k * size_per_prop) + j];
}
write_i++;
}
}
return out_buffer;
}
void poet::ChemistryModule::unshuffleField(const std::vector<double> &in_buffer,
uint32_t size_per_prop,
uint32_t prop_count,
uint32_t wp_count,
std::vector<double> &out_field) {
uint32_t read_i = 0;
for (uint32_t i = 0; i < wp_count; i++) {
for (uint32_t j = i; j < size_per_prop; j += wp_count) {
for (uint32_t k = 0; k < prop_count; k++) {
out_field[(k * size_per_prop) + j] =
in_buffer[(read_i * prop_count) + k];
}
read_i++;
}
}
}
}

12
src/Chemistry/ChemistryModule.hpp
View File

@ -51,7 +51,7 @@ public:
*/
~ChemistryModule();
void masterSetField(Field field);
void masterSetField(Field<TugType> &&field);
/**
* Run the chemical simulation with parameters set.
*/
@ -208,7 +208,7 @@ public:
*
* \return Reference to the chemical field.
*/
Field &getField() { return this->chem_field; }
Field<TugType> &getField() { return *this->chem_field; }
/**
* **Master only** Enable/disable progress bar.
@ -320,12 +320,6 @@ protected:
std::vector<uint32_t> CalculateWPSizesVector(uint32_t n_cells,
uint32_t wp_size) const;
std::vector<double> shuffleField(const std::vector<double> &in_field,
uint32_t size_per_prop, uint32_t prop_count,
uint32_t wp_count);
void unshuffleField(const std::vector<double> &in_buffer,
uint32_t size_per_prop, uint32_t prop_count,
uint32_t wp_count, std::vector<double> &out_field);
std::vector<std::int32_t>
parseDHTSpeciesVec(const NamedVector<std::uint32_t> &key_species,
const std::vector<std::string> &to_compare) const;
@ -374,7 +368,7 @@ protected:
uint32_t prop_count = 0;
std::vector<std::string> prop_names;
Field chem_field;
std::unique_ptr<Field<TugType>> chem_field;
const InitialList::ChemistryInit params;

33
src/Chemistry/MasterFunctions.cpp
View File

@ -1,9 +1,13 @@
#include "ChemistryModule.hpp"
#include "DataStructures/Field.hpp"
#include "Init/InitialList.hpp"
#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <memory>
#include <mpi.h>
#include <span>
#include <vector>
std::vector<uint32_t>
@ -159,10 +163,9 @@ std::vector<uint32_t> poet::ChemistryModule::GetWorkerPHTCacheHits() const {
return ret;
}
inline std::vector<double> shuffleField(const std::vector<double> &in_field,
uint32_t size_per_prop,
uint32_t prop_count,
uint32_t wp_count) {
static inline std::vector<double>
shuffleField(const std::span<poet::TugType> &in_field, uint32_t size_per_prop,
uint32_t prop_count, uint32_t wp_count) {
std::vector<double> out_buffer(in_field.size());
uint32_t write_i = 0;
for (uint32_t i = 0; i < wp_count; i++) {
@ -177,9 +180,10 @@ inline std::vector<double> shuffleField(const std::vector<double> &in_field,
return out_buffer;
}
inline void unshuffleField(const std::vector<double> &in_buffer,
uint32_t size_per_prop, uint32_t prop_count,
uint32_t wp_count, std::vector<double> &out_field) {
static inline void unshuffleField(const std::vector<double> &in_buffer,
uint32_t size_per_prop, uint32_t prop_count,
uint32_t wp_count,
std::span<poet::TugType> &&out_field) {
uint32_t read_i = 0;
for (uint32_t i = 0; i < wp_count; i++) {
@ -368,9 +372,8 @@ void poet::ChemistryModule::MasterRunParallel(double dt) {
/* shuffle grid */
// grid.shuffleAndExport(mpi_buffer);
std::vector<double> mpi_buffer =
shuffleField(chem_field.AsVector(), this->n_cells, this->prop_count,
wp_sizes_vector.size());
std::vector<double> mpi_buffer = shuffleField(
*chem_field, this->n_cells, this->prop_count, wp_sizes_vector.size());
/* setup local variables */
pkg_to_send = wp_sizes_vector.size();
@ -418,10 +421,8 @@ void poet::ChemistryModule::MasterRunParallel(double dt) {
/* unshuffle grid */
// grid.importAndUnshuffle(mpi_buffer);
std::vector<double> out_vec{mpi_buffer};
unshuffleField(mpi_buffer, this->n_cells, this->prop_count,
wp_sizes_vector.size(), out_vec);
chem_field = out_vec;
wp_sizes_vector.size(), *chem_field);
/* do master stuff */
@ -464,9 +465,9 @@ 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();
void poet::ChemistryModule::masterSetField(Field<TugType> &&field) {
this->chem_field = std::make_unique<Field<TugType>>(field);
this->prop_count = chem_field->size();
int ftype = CHEM_FIELD_INIT;
PropagateFunctionType(ftype);

3
src/DataStructures/CMakeLists.txt
View File

@ -1,3 +0,0 @@
target_sources(POETLib PRIVATE
Field.cpp
)

201
src/DataStructures/Field.cpp
View File

@ -1,201 +0,0 @@
#include "Field.hpp"
#include <Rcpp.h>
#include <Rcpp/DataFrame.h>
#include <Rcpp/exceptions.h>
#include <Rcpp/vector/instantiation.h>
#include <Rinternals.h>
#include <cstddef>
#include <cstdint>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
poet::Field::Field(std::uint32_t vec_s,
const std::vector<std::vector<double>> &input,
const std::vector<std::string> &prop_vec, bool row_major)
: req_vec_size(vec_s), props(prop_vec) {
if (row_major) {
if (this->props.size() != input[0].size()) {
throw std::runtime_error("Prop vector shall name all elements.");
}
const std::size_t n_input = input.size();
for (std::size_t i = 0; i < this->props.size(); i++) {
std::vector<double> curr_col(n_input);
for (std::size_t j = 0; j < n_input; j++) {
curr_col[j] = input[j][i];
}
this->insert({this->props[i], std::move(curr_col)});
}
} else {
if (this->props.size() != input.size()) {
throw std::runtime_error("Prop vector shall name all elements.");
}
auto name_it = this->props.begin();
for (const auto &in_vec : input) {
if (in_vec.size() != req_vec_size) {
throw std::runtime_error(
"Input vector doesn't match expected vector size.");
}
this->insert({*(name_it++), in_vec});
}
}
}
poet::Field::Field(std::uint32_t vec_s, const std::vector<double> &input,
const std::vector<std::string> &prop_vec)
: Field(vec_s) {
const uint32_t expected_size = prop_vec.size() * req_vec_size;
if (expected_size != input.size()) {
throw std::runtime_error(
"Input vector have more (or less) elements than expected.");
}
auto name_it = prop_vec.begin();
for (uint32_t i = 0; i < expected_size; i += req_vec_size) {
auto input_pair = std::make_pair(
*(name_it++), std::vector<double>(&input[i], &input[i + req_vec_size]));
this->insert(input_pair);
}
this->props = prop_vec;
}
auto poet::Field::AsVector() const -> poet::FieldColumn {
const uint32_t min_size = req_vec_size * this->size();
poet::FieldColumn output;
output.reserve(min_size);
for (const auto &elem : props) {
const auto map_it = this->find(elem);
const auto start = map_it->second.begin();
const auto end = map_it->second.end();
output.insert(output.end(), start, end);
}
return output;
}
void poet::Field::update(const poet::Field &input) {
for (const auto &input_elem : input) {
auto it_self = this->find(input_elem.first);
if (it_self == this->end()) {
continue;
}
it_self->second = input_elem.second;
}
}
auto poet::Field::As2DVector() const -> std::vector<poet::FieldColumn> {
std::vector<poet::FieldColumn> output;
output.reserve(this->size());
for (const auto &elem : props) {
const auto map_it = this->find(elem);
output.push_back(map_it->second);
}
return output;
}
poet::FieldColumn &poet::Field::operator[](const std::string &key) {
if (this->find(key) == this->end()) {
props.push_back(key);
}
return std::unordered_map<std::string, FieldColumn>::operator[](key);
}
SEXP poet::Field::asSEXP() const {
Rcpp::List output;
for (const auto &elem : this->props) {
const auto map_it = this->find(elem);
output[elem] = Rcpp::wrap(map_it->second);
}
return output;
}
poet::Field &poet::Field::operator=(const FieldColumn &cont_field) {
if (cont_field.size() != this->size() * req_vec_size) {
throw std::runtime_error(
"Field::SetFromVector: vector does not match expected size");
}
uint32_t vec_p = 0;
for (const auto &elem : props) {
const auto start = cont_field.begin() + vec_p;
const auto end = start + req_vec_size;
const auto map_it = this->find(elem);
map_it->second = FieldColumn(start, end);
vec_p += req_vec_size;
}
return *this;
}
poet::Field &
poet::Field::operator=(const std::vector<FieldColumn> &cont_field) {
if (cont_field.size() != this->size()) {
throw std::runtime_error(
"Input field contains more or less elements than this container.");
}
auto in_field_it = cont_field.begin();
for (const auto &elem : props) {
if (in_field_it->size() != req_vec_size) {
throw std::runtime_error(
"One vector contains more or less elements than expected.");
}
const auto map_it = this->find(elem);
map_it->second = *(in_field_it++);
}
return *this;
}
void poet::Field::fromSEXP(const SEXP &s_rhs) {
this->clear();
Rcpp::List in_list;
try {
in_list = Rcpp::List(s_rhs);
} catch (const Rcpp::exception &e) {
throw std::runtime_error("Input cannot be casted as list.");
}
if (in_list.size() == 0) {
return;
}
this->props = Rcpp::as<std::vector<std::string>>(in_list.names());
const Rcpp::NumericVector &in_vec = in_list[0];
this->req_vec_size = static_cast<std::uint32_t>(in_vec.size());
for (const auto &elem : this->props) {
const auto values = Rcpp::as<std::vector<double>>(in_list[elem]);
this->insert({elem, values});
}
}

288
src/DataStructures/Field.hpp
View File

@ -1,182 +1,154 @@
#ifndef DATASTRUCTURES_H_
#define DATASTRUCTURES_H_
#pragma once
#include <Rcpp.h>
#include <Rcpp/vector/Matrix.h>
#include <Rcpp/vector/instantiation.h>
#include <Rinternals.h>
#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <string>
#include <cstring>
#include <iterator>
#include <optional>
#include <span>
#include <utility>
#include <vector>
#include "DataStructures/PinnableMemory.hpp"
namespace poet {
using FieldColumn = std::vector<double>;
/**
* Stores data for input/output of a module. The class keeps track of all
* defined properties with name and import/export to 1D and 2D STL vectors.
* Also, it eases the update process of a field with an input from another
* field.
*
* It can be seen as an R-like data frame, but with less access to the members.
* Species values are stored as "columns", where column is a STL vector.
*/
class Field : private std::unordered_map<std::string, FieldColumn> {
template <class T> class Field {
public:
Field(){};
Field() = delete;
Field(Field &&) = delete;
Field &operator=(Field &&) = delete;
/**
* Creates a new instance of a field with fixed expected vector size.
*
* \param vec_s expected vector size of each component/column.
*/
Field(std::uint32_t vec_s) : req_vec_size(vec_s){};
/**
* Initializes instance with a 2D vector and according names for each columnn.
* There is no check if names were given multiple times. The order of name
* vector also defines the ordering of the output.
*
* \param vec_s expected vector size of each component/column.
* \param input 2D vector using STL semantic describing the current state of
* the field.
* \param prop_vec Name of each vector in the input. Shall match
* the count of vectors.
*
* \exception std::runtime_error If prop_vec size doesn't match input vector
* size or column vectors size doesn't match expected vector size.
*/
Field(std::uint32_t vec_s, const std::vector<std::vector<double>> &input,
const std::vector<std::string> &prop_vec, bool row_major = false);
/**
* Initializes instance with a 1D continious memory vector and according names
* for each columnn. There is no check if names were given multiple times. The
* order of name vector also defines the ordering of the output.
*
* \param vec_s expected vector size of each component/column.
* \param input 1D vector using STL semantic describing the current state of
* the field storing each column starting at index *i times requested vector
* size*.
* \param prop_vec Name of each vector in the input. Shall match the
* count of vectors.
*
* \exception std::runtime_error If prop_vec size doesn't match input vector
* size or column vectors size doesn't match expected vector size.
*/
Field(std::uint32_t vec_s, const std::vector<double> &input,
const std::vector<std::string> &prop_vec);
Field(const SEXP &s_init) { fromSEXP(s_init); }
Field &operator=(const Field &rhs) {
this->req_vec_size = rhs.req_vec_size;
this->props = rhs.props;
this->clear();
for (const auto &pair : rhs) {
this->insert(pair);
Field(const SEXP &s_rhs) {
Rcpp::List in_list;
try {
in_list = Rcpp::List(s_rhs);
} catch (const Rcpp::not_compatible &) {
throw std::runtime_error("Input is not a list.");
}
if (in_list.size() == 0) {
throw std::runtime_error("Input list is empty.");
}
Rcpp::NumericVector first_vec = Rcpp::NumericVector(in_list[0]);
_m_vec_size = first_vec.size();
// check if all elements are vectors of the same size
for (const Rcpp::NumericVector &vec : in_list) {
if (vec.size() != _m_vec_size) {
throw std::runtime_error("Input vectors have different sizes.");
}
}
_m_data.emplace(in_list.size() * _m_vec_size);
// copy data
for (std::size_t i = 0; i < in_list.size(); i++) {
Rcpp::NumericVector curr_vec = Rcpp::NumericVector(in_list[i]);
for (std::size_t j = 0; j < _m_vec_size; j++) {
_m_data.value()[i * _m_vec_size + j] = curr_vec[j];
}
}
// get the column names
Rcpp::CharacterVector colnames = in_list.names();
if (colnames.size() != in_list.size()) {
throw std::runtime_error(
"Column names do not match the number of columns.");
}
_m_colnames.resize(colnames.size());
for (std::size_t i = 0; i < colnames.size(); i++) {
_m_colnames[i] = (Rcpp::as<std::string>(colnames[i]));
}
}
Field &operator=(Field &rhs) {
if (this == &rhs) {
return *this;
}
if (_m_vec_size != rhs._m_vec_size) {
throw std::runtime_error("Vector sizes do not match.");
}
for (const auto &colname : _m_colnames) {
auto it =
std::find(rhs._m_colnames.begin(), rhs._m_colnames.end(), colname);
if (it == rhs._m_colnames.end()) {
continue;
}
std::span<T> rhs_col = rhs[colname];
std::span<T> this_col = this->operator[](colname);
std::memcpy(this_col.data(), rhs_col.data(), rhs_col.size_bytes());
}
return *this;
}
/**
* Read in a (previously exported) 1D vector. It has to have the same
* dimensions as the current column count times the requested vector size of
* this instance. Import order is given by the species name vector.
*
* \param cont_field 1D field as vector.
*
* \exception std::runtime_error Input vector does not match the expected
* size.
*/
Field &operator=(const FieldColumn &cont_field);
/**
* Read in a (previously exported) 2D vector. It has to have the same
* dimensions as the current column count of this instance and each vector
* must have the size of the requested vector size. Import order is given by
* the species name vector.
*
* \param cont_field 2D field as vector.
*
* \exception std::runtime_error Input vector has more or less elements than
* the instance or a column vector does not match expected vector size.
*/
Field &operator=(const std::vector<FieldColumn> &cont_field);
Field &operator=(const SEXP &s_rhs) {
fromSEXP(s_rhs);
return *this;
Field(Field &rhs)
: _m_vec_size(rhs._m_vec_size), _m_colnames(rhs._m_colnames) {
_m_data.emplace(rhs._m_data.value().size());
std::memcpy(_m_data.value().data(), rhs._m_data.value().data(),
rhs._m_data.value().size_bytes());
}
/**
* Returns a reference to the column vector with given name. Creates a new
* vector if prop was not found. The prop name will be added to the end of the
* list.
*
* \param key Name of the prop.
*
* \return Reference to the column vector.
*/
FieldColumn &operator[](const std::string &key);
std::span<T> operator[](const std::string &colname) {
auto it = std::find(_m_colnames.begin(), _m_colnames.end(), colname);
if (it == _m_colnames.end()) {
throw std::runtime_error("Column name not found.");
}
/**
* Returns the names of all species defined in the instance.
*
* \return Vector containing all species names in output order.
*/
auto GetProps() const -> std::vector<std::string> { return this->props; }
const std::size_t col_idx = std::distance(_m_colnames.begin(), it);
/**
* Return the requested vector size.
*
* \return Requested vector size set in the instanciation of the object.
*/
auto GetRequestedVecSize() const -> std::uint32_t {
return this->req_vec_size;
};
return this->operator[](col_idx);
}
/**
* Updates all species with values from another field. If one element of the
* input field doesn't match the names of the calling instance it will get
* skipped.
*
* \param input Field to update the current instance's columns.
*/
void update(const Field &input);
std::span<T> operator[](std::size_t col_idx) {
if (col_idx >= _m_colnames.size()) {
throw std::runtime_error("Column index out of bounds.");
}
/**
* Builds a new 1D vector with the current state of the instance. The output
* order is given by the given species name vector set earlier and/or added
* values using the [] operator.
*
* \return 1D STL vector stored each column one after another.
*/
auto AsVector() const -> FieldColumn;
return std::span<T>(&_m_data.value()[col_idx * _m_vec_size], _m_vec_size);
}
/**
* Builds a new 2D vector with the current state of the instance. The output
* order is given by the given species name vector set earlier and/or added
* values using the [] operator.
*
* \return 2D STL vector stored each column one after anothe in a new vector
* element.
*/
auto As2DVector() const -> std::vector<FieldColumn>;
operator std::span<T>() { return std::span<T>(_m_data.value()); }
SEXP asSEXP() const;
operator SEXP() {
Rcpp::List out_list(_m_colnames.size());
std::size_t cols() const { return this->props.size(); }
for (std::size_t i = 0; i < _m_colnames.size(); i++) {
Rcpp::NumericVector vec(_m_vec_size);
std::memcpy(vec.begin(), _m_data.value().data() + i * _m_vec_size,
_m_vec_size * sizeof(T));
out_list[i] = vec;
}
out_list.names() =
Rcpp::CharacterVector(_m_colnames.begin(), _m_colnames.end());
return std::move(out_list);
}
std::size_t size() const { return _m_colnames.size(); }
std::size_t rows() const { return _m_vec_size; }
T *data() { return _m_data.value().data(); }
const std::vector<std::string> &colnames() const { return _m_colnames; }
private:
std::uint32_t req_vec_size{0};
std::vector<std::string> props{{}};
void fromSEXP(const SEXP &s_rhs);
std::optional<PinnableMemory<T>> _m_data;
std::uint32_t _m_vec_size;
std::vector<std::string> _m_colnames;
};
} // namespace poet
#endif // DATASTRUCTURES_H_
} // namespace poet

14
src/DataStructures/PinnableMemory.hpp
View File

@ -33,6 +33,20 @@ public:
*/
T *data() { return _m_data.get(); }
/**
* @brief Return the count of elements of type T allocated
*
* @return std::size_t Count of elements allocated
*/
std::size_t size() const { return _m_count; }
/**
* @brief Return the size of the allocated memory in bytes
*
* @return std::size_t Size of the allocated memory in bytes
*/
std::size_t size_bytes() const { return _m_count * sizeof(T); }
/**
* @brief Returns the element at the given index
*

7
src/Init/DiffusionInit.cpp
View File

@ -379,7 +379,8 @@ RcppListToInnerBoundaryMap(const std::vector<std::string> &trans_names,
}
InitialList::DiffusionInit InitialList::getDiffusionInit() const {
DiffusionInit diff_init;
DiffusionInit diff_init = {.alpha_x = Field<TugType>(this->alpha_x),
.alpha_y = Field<TugType>(this->alpha_y)};
diff_init.dim = this->dim;
@ -397,8 +398,8 @@ InitialList::DiffusionInit InitialList::getDiffusionInit() const {
diff_init.inner_boundaries =
RcppListToInnerBoundaryMap(this->transport_names, this->inner_boundaries,
this->n_cols, this->n_rows);
diff_init.alpha_x = Field(this->alpha_x);
diff_init.alpha_y = Field(this->alpha_y);
// diff_init.alpha_x = Field(this->alpha_x);
// diff_init.alpha_y = Field(this->alpha_y);
return diff_init;
}

8
src/Init/InitialList.cpp
View File

@ -12,6 +12,10 @@ namespace poet {
void InitialList::initializeFromList(const Rcpp::List &setup) {
prepareGrid(setup[grid_key]);
initDiffusion(setup[diffusion_key]);
this->grid_colnames =
Rcpp::as<std::vector<std::string>>(this->initial_grid.names());
initChemistry(setup[chemistry_key]);
}
@ -39,6 +43,9 @@ void InitialList::importList(const Rcpp::List &setup, bool minimal) {
Rcpp::List(setup[static_cast<int>(ExportList::GRID_INITIAL)]);
}
this->grid_colnames = Rcpp::as<std::vector<std::string>>(
setup[static_cast<int>(ExportList::GRID_COLNAMES)]);
this->transport_names = Rcpp::as<std::vector<std::string>>(
setup[static_cast<int>(ExportList::DIFFU_TRANSPORT)]);
if (!minimal) {
@ -96,6 +103,7 @@ Rcpp::List InitialList::exportList() {
Rcpp::wrap(this->constant_cells);
out[static_cast<int>(ExportList::GRID_POROSITY)] = Rcpp::wrap(this->porosity);
out[static_cast<int>(ExportList::GRID_INITIAL)] = this->initial_grid;
out[static_cast<int>(ExportList::GRID_COLNAMES)] = this->grid_colnames;
out[static_cast<int>(ExportList::DIFFU_TRANSPORT)] =
Rcpp::wrap(this->transport_names);

15
src/Init/InitialList.hpp
View File

@ -35,7 +35,13 @@ public:
void importList(const Rcpp::List &setup, bool minimal = false);
Rcpp::List exportList();
Field getInitialGrid() const { return Field(this->initial_grid); }
Field<TugType> getInitialGrid() const {
return Field<TugType>(this->initial_grid);
}
std::vector<std::string> getGridColnames() const {
return this->grid_colnames;
}
private:
RInside &R;
@ -47,6 +53,7 @@ private:
GRID_CONSTANT,
GRID_POROSITY,
GRID_INITIAL,
GRID_COLNAMES,
DIFFU_TRANSPORT,
DIFFU_BOUNDARIES,
DIFFU_INNER_BOUNDARIES,
@ -113,6 +120,8 @@ private:
Rcpp::List initial_grid;
std::vector<std::string> grid_colnames;
// No export
Rcpp::NumericMatrix phreeqc_mat;
@ -139,8 +148,8 @@ public:
BoundaryMap boundaries;
InnerBoundaryMap inner_boundaries;
Field alpha_x;
Field alpha_y;
Field<TugType> alpha_x;
Field<TugType> alpha_y;
};
DiffusionInit getDiffusionInit() const;

20
src/Transport/DiffusionModule.cpp
View File

@ -40,22 +40,18 @@
using namespace poet;
static inline std::vector<TugType>
MatrixToVec(const Eigen::MatrixX<TugType> &mat) {
std::vector<TugType> vec(mat.rows() * mat.cols());
static inline void MatrixToVec(const Eigen::MatrixX<TugType> &mat,
std::span<TugType> &out_vec) {
for (std::uint32_t i = 0; i < mat.cols(); i++) {
for (std::uint32_t j = 0; j < mat.rows(); j++) {
vec[j * mat.cols() + i] = mat(j, i);
out_vec[j * mat.cols() + i] = mat(j, i);
}
}
return vec;
}
static inline Eigen::MatrixX<TugType>
VecToMatrix(const std::vector<TugType> &vec, std::uint32_t n_rows,
std::uint32_t n_cols) {
static inline Eigen::MatrixX<TugType> VecToMatrix(const std::span<TugType> &vec,
std::uint32_t n_rows,
std::uint32_t n_cols) {
Eigen::MatrixX<TugType> mat(n_rows, n_cols);
for (std::uint32_t i = 0; i < n_cols; i++) {
@ -85,7 +81,7 @@ void DiffusionModule::simulate(double requested_dt) {
sim.setIterations(1);
for (const auto &sol_name : this->param_list.transport_names) {
auto &species_conc = this->transport_field[sol_name];
auto species_conc = this->transport_field[sol_name];
Eigen::MatrixX<TugType> conc = VecToMatrix(species_conc, n_rows, n_cols);
Eigen::MatrixX<TugType> alpha_x =
@ -106,7 +102,7 @@ void DiffusionModule::simulate(double requested_dt) {
sim.run();
species_conc = MatrixToVec(grid.getConcentrations());
MatrixToVec(grid.getConcentrations(), species_conc);
}
const auto end_diffusion_t = std::chrono::high_resolution_clock::now();

7
src/Transport/DiffusionModule.hpp
View File

@ -46,7 +46,8 @@ public:
*
* @param R RRuntime object
*/
DiffusionModule(const InitialList::DiffusionInit &init_list, Field field)
DiffusionModule(InitialList::DiffusionInit &&init_list,
Field<TugType> &&field)
: param_list(init_list), transport_field(field){};
/**
@ -69,7 +70,7 @@ public:
*
* \return Reference to the diffusion field.
*/
Field &getField() { return this->transport_field; }
Field<TugType> &getField() { return this->transport_field; }
private:
/**
@ -79,7 +80,7 @@ private:
InitialList::DiffusionInit param_list;
Field transport_field;
Field<TugType> transport_field;
/**
* @brief time spent for transport

31
src/poet.cpp
View File

@ -36,6 +36,7 @@
#include <memory>
#include <mpi.h>
#include <optional>
#include <span>
#include <string>
#include "Base/argh.hpp"
@ -225,18 +226,17 @@ ParseRet parseInitValues(char **argv, RuntimeParameters &params) {
// 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 call_master_iter_end(RInside &R, const Field &trans, const Field &chem) {
R["TMP"] = Rcpp::wrap(trans.AsVector());
R["TMP_PROPS"] = Rcpp::wrap(trans.GetProps());
void call_master_iter_end(RInside &R, Field<TugType> &trans,
Field<TugType> &chem) {
R["TMP"] = Rcpp::wrap(static_cast<std::span<double>>(trans));
R["TMP_PROPS"] = Rcpp::wrap(trans.colnames());
R.parseEval(std::string("state_T <- setNames(data.frame(matrix(TMP, nrow=" +
std::to_string(trans.GetRequestedVecSize()) +
")), TMP_PROPS)"));
std::to_string(trans.rows()) + ")), TMP_PROPS)"));
R["TMP"] = Rcpp::wrap(chem.AsVector());
R["TMP_PROPS"] = Rcpp::wrap(chem.GetProps());
R["TMP"] = Rcpp::wrap(static_cast<std::span<double>>(chem));
R["TMP_PROPS"] = Rcpp::wrap(chem.colnames());
R.parseEval(std::string("state_C <- setNames(data.frame(matrix(TMP, nrow=" +
std::to_string(chem.GetRequestedVecSize()) +
")), TMP_PROPS)"));
std::to_string(chem.rows()) + ")), TMP_PROPS)"));
R["setup"] = *global_rt_setup;
R.parseEval("setup <- master_iteration_end(setup, state_T, state_C)");
*global_rt_setup = R["setup"];
@ -271,7 +271,8 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters &params,
/* run transport */
diffusion.simulate(dt);
chem.getField().update(diffusion.getField());
chem.getField() = diffusion.getField();
;
// chem.getfield().update(diffusion.getfield());
@ -285,7 +286,7 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, const RuntimeParameters &params,
// store_result is TRUE)
call_master_iter_end(R, diffusion.getField(), chem.getField());
diffusion.getField().update(chem.getField());
diffusion.getField() = chem.getField();
MSG("End of *coupling* iteration " + std::to_string(iter) + "/" +
std::to_string(maxiter));
@ -378,7 +379,7 @@ int main(int argc, char *argv[]) {
init_list.getChemistryInit(), MPI_COMM_WORLD);
const ChemistryModule::SurrogateSetup surr_setup = {
init_list.getInitialGrid().GetProps(),
init_list.getGridColnames(),
run_params.use_dht,
run_params.dht_size,
run_params.use_interp,
@ -397,9 +398,9 @@ int main(int argc, char *argv[]) {
// R.parseEvalQ("mysetup <- setup");
// // if (MY_RANK == 0) { // get timestep vector from
// // grid_init function ... //
*global_rt_setup =
master_init_R.value()(*global_rt_setup, run_params.out_dir,
init_list.getInitialGrid().asSEXP());
*global_rt_setup = master_init_R.value()(
*global_rt_setup, run_params.out_dir,
static_cast<std::span<double>>(init_list.getInitialGrid()));
// MDL: store all parameters
// MSG("Calling R Function to store calling parameters");

2
test/CMakeLists.txt
View File

@ -15,4 +15,4 @@ target_include_directories(testPOET PRIVATE "${CMAKE_CURRENT_BINARY_DIR}")
add_custom_target(check
COMMAND $<TARGET_FILE:testPOET>
DEPENDS testPOET
)
)

2
test/RInsidePOET_funcs.R
View File

@ -14,7 +14,7 @@ bool_named_vec <- function(input) {
simple_field <- function(field) {
field <- as.data.frame(field, check.names = FALSE)
field$Na <- 0
field$C <- 0
return(field)
}

395
test/testField.cpp
View File

@ -1,4 +1,5 @@
#include <Rcpp.h>
#include <Rcpp/vector/instantiation.h>
#include <cstddef>
#include <cstdint>
#include <doctest/doctest.h>
@ -18,194 +19,248 @@ using namespace poet;
}
TEST_CASE("Field") {
constexpr uint32_t VEC_SIZE = 5;
constexpr uint32_t VEC_COUNT = 3;
constexpr double INIT_VAL = 1;
std::vector<std::string> names = {"C", "Ca", "Na"};
std::vector<FieldColumn> init_values(names.size(),
FieldColumn(VEC_SIZE, INIT_VAL));
RInsidePOET &R = RInsidePOET::getInstance();
R["sourcefile"] = RInside_source_file;
R.parseEval("source(sourcefile)");
SUBCASE("Initialize field with 2D vector") {
Field dut(VEC_SIZE, init_values, names);
R.parseEval("mylist <- list(C = c(1, 2, 3), Ca = c(4, 5, 6))");
auto props = dut.GetProps();
Rcpp::List in_list = R.parseEval("mylist");
CHECK_EQ(names.size(), props.size());
Field<double> dut(in_list);
const auto res = dut["C"];
SUBCASE("operator[](string)") {
auto res = dut["C"];
CHECK_EQ(res.size(), VEC_SIZE);
CHECK_EQ(res.size(), 3);
for (const auto &elem : res) {
CHECK_AND_FAIL_LOOP(elem, INIT_VAL);
CHECK_EQ(res[0], 1);
CHECK_EQ(res[1], 2);
CHECK_EQ(res[2], 3);
CHECK_THROWS(dut["D"]);
}
SUBCASE("operator std::span<T>") {
std::span<double> res = dut;
CHECK_EQ(res.size(), 6);
for (std::size_t i = 0; i < res.size(); i++) {
CHECK_EQ(res[i], i + 1);
}
}
SUBCASE("Initialize field with 2D vector") {
std::vector<double> init_values(VEC_SIZE * VEC_COUNT, 1);
Field dut(VEC_SIZE, init_values, names);
SUBCASE("operator =(Field)") {
R.parseEval(
"mylist2 <- list(C = c(1, 1, 1), Ca = c(2, 2, 2), H = c(7, 8, 9))");
Rcpp::List new_in_list = R.parseEval("mylist2");
auto props = dut.GetProps();
Field<double> other_field(new_in_list);
CHECK_EQ(names.size(), props.size());
CHECK_NOTHROW(dut = other_field);
const auto res = dut["C"];
auto res = dut["Ca"];
CHECK_EQ(res.size(), VEC_SIZE);
CHECK_EQ(res.size(), 3);
for (const auto &elem : res) {
CHECK_AND_FAIL_LOOP(elem, INIT_VAL);
}
}
Field dut(VEC_SIZE, init_values, names);
SUBCASE("Return vector") {
std::vector<double> output = dut.AsVector();
CHECK(output.size() == VEC_SIZE * VEC_COUNT);
}
SUBCASE("SEXP return") {
R["test"] = dut.asSEXP();
Field field_constructed = R.parseEval("test");
Rcpp::DataFrame R_df = R.parseEval("test");
Rcpp::CharacterVector R_names = R_df.names();
for (std::size_t i = 0; i < VEC_COUNT; i++) {
const auto r_values = Rcpp::as<std::vector<double>>(R_df[i]);
CHECK_EQ(r_values, dut[names[i]]);
CHECK_EQ(R_names[i], names[i]);
CHECK_EQ(field_constructed[names[i]], dut[names[i]]);
}
}
SUBCASE("Apply R function (set Na to zero)") {
RHookFunction<Field> to_call(R, "simple_field");
Field field_proc = to_call(dut.asSEXP());
CHECK_EQ(field_proc["Na"], FieldColumn(dut.GetRequestedVecSize(), 0));
}
SUBCASE("Apply R function (add two fields)") {
RHookFunction<Field> to_call(R, "extended_field");
Field field_proc = to_call(dut.asSEXP(), dut.asSEXP());
CHECK_EQ(field_proc["Na"],
FieldColumn(dut.GetRequestedVecSize(), INIT_VAL + INIT_VAL));
}
constexpr double NEW_VAL = 2.;
std::vector<double> new_val_vec(VEC_SIZE, NEW_VAL);
dut["C"] = new_val_vec;
SUBCASE("Check manipulation of column") {
auto test_it = new_val_vec.begin();
for (const auto &val : dut["C"]) {
CHECK_EQ(val, *test_it);
test_it++;
}
}
SUBCASE("Check correctly manipulated values of 1D vector") {
auto out_res = dut.AsVector();
auto out_it = out_res.begin();
for (uint32_t i = 0; i < VEC_SIZE; i++, out_it++) {
CHECK_AND_FAIL_LOOP(NEW_VAL, *out_it);
}
for (; out_it != out_res.end(); out_it++) {
CHECK_AND_FAIL_LOOP(INIT_VAL, *out_it);
}
}
std::vector<double> new_field(VEC_SIZE * VEC_COUNT);
for (uint32_t i = 0; i < VEC_COUNT; i++) {
for (uint32_t j = 0; j < VEC_SIZE; j++) {
new_field[j + (i * VEC_SIZE)] = (double)(i + 1) / (double)(j + 1);
}
}
dut = new_field;
SUBCASE("SetFromVector return correct field vector") {
auto ret_vec = dut.AsVector();
auto ret_it = ret_vec.begin();
auto new_it = new_field.begin();
for (; ret_it != ret_vec.end(); ret_it++, new_it++) {
CHECK_AND_FAIL_LOOP(*ret_it, *new_it);
}
}
SUBCASE("Get single column with new values") {
auto new_it = new_field.begin() + 2 * VEC_SIZE;
auto ret_vec = dut["Na"];
auto ret_it = ret_vec.begin();
CHECK_EQ(ret_vec.size(), VEC_SIZE);
for (; ret_it != ret_vec.end(); ret_it++, new_it++) {
CHECK_AND_FAIL_LOOP(*ret_it, *new_it);
}
}
std::vector<double> mg_vec = {1, 2, 3, 4, 5};
dut["Mg"] = mg_vec;
SUBCASE("Operator creates new element and places it at the end") {
auto new_props = dut.GetProps();
CHECK_EQ(new_props.size(), 4);
CHECK_EQ(new_props[3], "Mg");
auto ret_vec = dut.As2DVector();
auto mg_vec_it = mg_vec.begin();
for (const auto &val : ret_vec[3]) {
CHECK_AND_FAIL_LOOP(val, *(mg_vec_it++));
}
}
// reset field
names = dut.GetProps();
dut = std::vector<FieldColumn>(names.size(), FieldColumn(VEC_SIZE, INIT_VAL));
constexpr double SOME_OTHER_VAL = -0.5;
std::vector<std::string> some_other_props = {"Na", "Cl"};
std::vector<double> some_other_values(VEC_SIZE * some_other_props.size(),
SOME_OTHER_VAL);
Field some_other_field(VEC_SIZE, some_other_values, some_other_props);
SUBCASE("Update existing field from another field") {
dut.update(some_other_field);
auto ret_vec = dut.As2DVector();
auto ret_it = ret_vec.begin();
for (const auto &prop : names) {
const auto &curr_vec = *(ret_it++);
for (const auto &val : curr_vec) {
CHECK_AND_FAIL_LOOP((prop == "Na" ? SOME_OTHER_VAL : INIT_VAL), val);
}
for (std::size_t i = 0; i < res.size(); i++) {
CHECK_EQ(res[i], 2);
}
}
}
// TEST_CASE("Field") {
// constexpr uint32_t VEC_SIZE = 5;
// constexpr uint32_t VEC_COUNT = 3;
// constexpr double INIT_VAL = 1;
// std::vector<std::string> names = {"C", "Ca", "Na"};
// std::vector<FieldColumn> init_values(names.size(),
// FieldColumn(VEC_SIZE, INIT_VAL));
// RInsidePOET &R = RInsidePOET::getInstance();
// R["sourcefile"] = RInside_source_file;
// R.parseEval("source(sourcefile)");
// SUBCASE("Initialize field with 2D vector") {
// Field dut(VEC_SIZE, init_values, names);
// auto props = dut.GetProps();
// CHECK_EQ(names.size(), props.size());
// const auto res = dut["C"];
// CHECK_EQ(res.size(), VEC_SIZE);
// for (const auto &elem : res) {
// CHECK_AND_FAIL_LOOP(elem, INIT_VAL);
// }
// }
// SUBCASE("Initialize field with 2D vector") {
// std::vector<double> init_values(VEC_SIZE * VEC_COUNT, 1);
// Field dut(VEC_SIZE, init_values, names);
// auto props = dut.GetProps();
// CHECK_EQ(names.size(), props.size());
// const auto res = dut["C"];
// CHECK_EQ(res.size(), VEC_SIZE);
// for (const auto &elem : res) {
// CHECK_AND_FAIL_LOOP(elem, INIT_VAL);
// }
// }
// Field dut(VEC_SIZE, init_values, names);
// SUBCASE("Return vector") {
// std::vector<double> output = dut.AsVector();
// CHECK(output.size() == VEC_SIZE * VEC_COUNT);
// }
// SUBCASE("SEXP return") {
// R["test"] = dut.asSEXP();
// Field field_constructed = R.parseEval("test");
// Rcpp::DataFrame R_df = R.parseEval("test");
// Rcpp::CharacterVector R_names = R_df.names();
// for (std::size_t i = 0; i < VEC_COUNT; i++) {
// const auto r_values = Rcpp::as<std::vector<double>>(R_df[i]);
// CHECK_EQ(r_values, dut[names[i]]);
// CHECK_EQ(R_names[i], names[i]);
// CHECK_EQ(field_constructed[names[i]], dut[names[i]]);
// }
// }
// SUBCASE("Apply R function (set Na to zero)") {
// RHookFunction<Field> to_call(R, "simple_field");
// Field field_proc = to_call(dut.asSEXP());
// CHECK_EQ(field_proc["Na"], FieldColumn(dut.GetRequestedVecSize(), 0));
// }
// SUBCASE("Apply R function (add two fields)") {
// RHookFunction<Field> to_call(R, "extended_field");
// Field field_proc = to_call(dut.asSEXP(), dut.asSEXP());
// CHECK_EQ(field_proc["Na"],
// FieldColumn(dut.GetRequestedVecSize(), INIT_VAL + INIT_VAL));
// }
// constexpr double NEW_VAL = 2.;
// std::vector<double> new_val_vec(VEC_SIZE, NEW_VAL);
// dut["C"] = new_val_vec;
// SUBCASE("Check manipulation of column") {
// auto test_it = new_val_vec.begin();
// for (const auto &val : dut["C"]) {
// CHECK_EQ(val, *test_it);
// test_it++;
// }
// }
// SUBCASE("Check correctly manipulated values of 1D vector") {
// auto out_res = dut.AsVector();
// auto out_it = out_res.begin();
// for (uint32_t i = 0; i < VEC_SIZE; i++, out_it++) {
// CHECK_AND_FAIL_LOOP(NEW_VAL, *out_it);
// }
// for (; out_it != out_res.end(); out_it++) {
// CHECK_AND_FAIL_LOOP(INIT_VAL, *out_it);
// }
// }
// std::vector<double> new_field(VEC_SIZE * VEC_COUNT);
// for (uint32_t i = 0; i < VEC_COUNT; i++) {
// for (uint32_t j = 0; j < VEC_SIZE; j++) {
// new_field[j + (i * VEC_SIZE)] = (double)(i + 1) / (double)(j + 1);
// }
// }
// dut = new_field;
// SUBCASE("SetFromVector return correct field vector") {
// auto ret_vec = dut.AsVector();
// auto ret_it = ret_vec.begin();
// auto new_it = new_field.begin();
// for (; ret_it != ret_vec.end(); ret_it++, new_it++) {
// CHECK_AND_FAIL_LOOP(*ret_it, *new_it);
// }
// }
// SUBCASE("Get single column with new values") {
// auto new_it = new_field.begin() + 2 * VEC_SIZE;
// auto ret_vec = dut["Na"];
// auto ret_it = ret_vec.begin();
// CHECK_EQ(ret_vec.size(), VEC_SIZE);
// for (; ret_it != ret_vec.end(); ret_it++, new_it++) {
// CHECK_AND_FAIL_LOOP(*ret_it, *new_it);
// }
// }
// std::vector<double> mg_vec = {1, 2, 3, 4, 5};
// dut["Mg"] = mg_vec;
// SUBCASE("Operator creates new element and places it at the end") {
// auto new_props = dut.GetProps();
// CHECK_EQ(new_props.size(), 4);
// CHECK_EQ(new_props[3], "Mg");
// auto ret_vec = dut.As2DVector();
// auto mg_vec_it = mg_vec.begin();
// for (const auto &val : ret_vec[3]) {
// CHECK_AND_FAIL_LOOP(val, *(mg_vec_it++));
// }
// }
// // reset field
// names = dut.GetProps();
// dut = std::vector<FieldColumn>(names.size(), FieldColumn(VEC_SIZE,
// INIT_VAL));
// constexpr double SOME_OTHER_VAL = -0.5;
// std::vector<std::string> some_other_props = {"Na", "Cl"};
// std::vector<double> some_other_values(VEC_SIZE * some_other_props.size(),
// SOME_OTHER_VAL);
// Field some_other_field(VEC_SIZE, some_other_values, some_other_props);
// SUBCASE("Update existing field from another field") {
// dut.update(some_other_field);
// auto ret_vec = dut.As2DVector();
// auto ret_it = ret_vec.begin();
// for (const auto &prop : names) {
// const auto &curr_vec = *(ret_it++);
// for (const auto &val : curr_vec) {
// CHECK_AND_FAIL_LOOP((prop == "Na" ? SOME_OTHER_VAL : INIT_VAL), val);
// }
// }
// }
// }