refactor: cleanup of code, renaming of chemsitry module function and some output added to the diffusion module

app/poet.cpp

@@ -158,8 +158,6 @@ int main(int argc, char *argv[]) {
       cout << "CPP: R's $iter: " << ((uint32_t)(R.parseEval("mysetup$iter")))
            << ". Iteration" << endl;
 
-      cout << "CPP: Calling Advection" << endl;
-
       /* run transport */
       // TODO: transport to diffusion
       diffusion.simulate(dt);
@@ -169,11 +167,11 @@ int main(int argc, char *argv[]) {
       /* Fallback for sequential execution */
       // TODO: use new grid
       if (world_size == 1) {
-        master.ChemSim::run(dt);
+        master.ChemSim::simulate(dt);
       }
       /* otherwise run parallel */
       else {
-        master.run(dt);
+        master.simulate(dt);
       }
 
       // MDL master_iteration_end just writes on disk state_T and

include/poet/ChemSim.hpp

@@ -80,7 +80,7 @@ public:
    * @todo change function name. Maybe 'slave' to 'seq'.
    *
    */
-  virtual void run(double dt);
+  virtual void simulate(double dt);
 
   /**
    * @brief End simulation
@@ -250,7 +250,7 @@ public:
    * The main tasks are instrumented with time measurements.
    *
    */
-  void run(double dt) override;
+  void simulate(double dt) override;
 
   /**
    * @brief End chemistry simulation.

src/ChemMaster.cpp

@@ -50,16 +50,6 @@ ChemMaster::ChemMaster(SimParams &params, RRuntime &R_, Grid &grid_)
       (wp_size * this->prop_names.size()) + BUFFER_OFFSET, sizeof(double));
   mpi_buffer = (double *)calloc(grid_size, sizeof(double));
 
-  // R.parseEvalQ("wp_ids <- distribute_work_packages(len=length(mysetup$prop),
-  // "
-  //              "package_size=work_package_size)");
-
-  // // we only sort once the vector
-  // R.parseEvalQ("ordered_ids <- order(wp_ids)");
-  // R.parseEvalQ("wp_sizes_vector <- compute_wp_sizes(wp_ids)");
-  // R.parseEval("stat_wp_sizes(wp_sizes_vector)");
-  // wp_sizes_vector = as<std::vector<int>>(R["wp_sizes_vector"]);
-
   /* calculate distribution of work packages */
   uint32_t mod_pkgs = grid_size % this->wp_size;
   uint32_t n_packages = (uint32_t)(grid.getTotalCellCount() / this->wp_size) +
@@ -77,7 +67,7 @@ ChemMaster::~ChemMaster() {
   free(workerlist);
 }
 
-void ChemMaster::run(double dt) {
+void ChemMaster::simulate(double dt) {
   /* declare most of the needed variables here */
   double chem_a, chem_b;
   double seq_a, seq_b, seq_c, seq_d;

src/ChemSim.cpp

@@ -46,7 +46,7 @@ ChemSim::ChemSim(SimParams &params, RRuntime &R_, Grid &grid_)
   this->n_cells_per_prop = this->grid.getTotalCellCount();
   this->state =
       this->grid.registerState(poet::CHEMISTRY_MODULE_NAME, this->prop_names);
-  auto &field = this->state->mem;
+  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++) {
@@ -58,7 +58,7 @@ ChemSim::ChemSim(SimParams &params, RRuntime &R_, Grid &grid_)
   }
 }
 
-void ChemSim::run(double dt) {
+void ChemSim::simulate(double dt) {
   double chem_a, chem_b;
 
   /* start time measuring */

src/ChemWorker.cpp

@@ -190,7 +190,6 @@ void ChemWorker::doWork(MPI_Status &probe_status) {
   }
 
   /* Convert grid to R runtime */
-  // grid.importWP(mpi_buffer, wp_size);
   size_t rowCount = local_work_package_size;
   size_t colCount = this->prop_names.size();
 

src/DiffusionModule.cpp

@@ -25,6 +25,7 @@
 #include <Rcpp.h>
 #include <algorithm>
 #include <cstdint>
+#include <ostream>
 #include <poet/DiffusionModule.hpp>
 #include <poet/Grid.hpp>
 
@@ -78,19 +79,18 @@ void DiffusionModule::initialize(poet::DiffusionParams args) {
 
   this->state =
       this->grid.registerState(DIFFUSION_MODULE_NAME, this->prop_names);
-  auto &field = this->state->mem;
 
-  // get alphas - we cannot assume alphas are provided in same order as initial
-  // input
-  std::vector<double> local_alpha(this->prop_count);
-  for (uint32_t i = 0; i < this->prop_count; i++) {
-    local_alpha[i] = args.alpha[this->prop_names[i]];
-  }
-  this->alpha = local_alpha;
+  std::vector<double> &field = this->state->mem;
 
-  // initialize field
+  // initialize field and alphas
   field.resize(this->n_cells_per_prop * this->prop_count);
+  this->alpha.reserve(this->prop_count);
+
   for (uint32_t i = 0; i < this->prop_count; i++) {
+    // get alphas - we cannot assume alphas are provided in same order as
+    // initial input
+    this->alpha.push_back(args.alpha[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));
@@ -104,6 +104,30 @@ void DiffusionModule::initialize(poet::DiffusionParams args) {
     }
   }
 
+  // apply boundary conditions to each ghost node
+  uint8_t bc_size = (this->dim == this->DIM_1D ? 2 : 4);
+  for (uint8_t i = 0; i < bc_size; i++) {
+    const auto &side = borders[i];
+    std::vector<uint32_t> vecinj_i = Rcpp::as<std::vector<uint32_t>>(
+        args.vecinj_index[convert_bc_to_config_file(side)]);
+    for (int i = 0; i < this->prop_count; i++) {
+      std::vector<double> bc_vec = args.vecinj[this->prop_names[i]];
+      tug::bc::BoundaryCondition &curr_bc = *(this->bc_vec.begin() + i);
+      for (int j = 0; j < vecinj_i.size(); j++) {
+        if (vecinj_i[j] == 0) {
+          continue;
+        }
+        if (this->dim == this->DIM_2D) {
+          curr_bc(side, j) = {tug::bc::BC_TYPE_CONSTANT,
+                              bc_vec[vecinj_i[j] - 1]};
+        }
+        if (this->dim == this->DIM_1D) {
+          curr_bc(side) = {tug::bc::BC_TYPE_CONSTANT, bc_vec[vecinj_i[j] - 1]};
+        }
+      }
+    }
+  }
+
   // apply inner grid constant cells
   // NOTE: opening a scope here for distinguish variable names
   if (args.vecinj_inner.rows() != 0) {
@@ -134,30 +158,6 @@ void DiffusionModule::initialize(poet::DiffusionParams args) {
       }
     }
   }
-
-  // apply boundary conditions to each ghost node
-  uint8_t bc_size = (this->dim == this->DIM_1D ? 2 : 4);
-  for (uint8_t i = 0; i < bc_size; i++) {
-    const auto &side = borders[i];
-    std::vector<uint32_t> vecinj_i = Rcpp::as<std::vector<uint32_t>>(
-        args.vecinj_index[convert_bc_to_config_file(side)]);
-    for (int i = 0; i < this->prop_count; i++) {
-      std::vector<double> bc_vec = args.vecinj[this->prop_names[i]];
-      tug::bc::BoundaryCondition &curr_bc = *(this->bc_vec.begin() + i);
-      for (int j = 0; j < vecinj_i.size(); j++) {
-        if (vecinj_i[j] == 0) {
-          continue;
-        }
-        if (this->dim == this->DIM_2D) {
-          curr_bc(side, j) = {tug::bc::BC_TYPE_CONSTANT,
-                              bc_vec[vecinj_i[j] - 1]};
-        }
-        if (this->dim == this->DIM_1D) {
-          curr_bc(side) = {tug::bc::BC_TYPE_CONSTANT, bc_vec[vecinj_i[j] - 1]};
-        }
-      }
-    }
-  }
 }
 
 void DiffusionModule::simulate(double dt) {
@@ -165,8 +165,13 @@ void DiffusionModule::simulate(double dt) {
 
   sim_b_transport = MPI_Wtime();
 
+  std::cout << "DiffusionModule::simulate(): Starting diffusion ..."
+            << std::flush;
+
   std::vector<double> &curr_field = this->state->mem;
 
+  // copy output of another module (in this case there is only the chemistry
+  // module) as input for diffusion
   for (uint32_t i = 0; i < this->prop_names.size(); i++) {
     try {
       std::vector<double> t_prop_vec = this->grid.getSpeciesByName(
@@ -175,16 +180,16 @@ void DiffusionModule::simulate(double dt) {
       std::copy(t_prop_vec.begin(), t_prop_vec.end(),
                 curr_field.begin() + (i * this->n_cells_per_prop));
     } catch (...) {
+      // TODO: there might be something wrong ...
       continue;
     }
   }
 
-  auto *field = curr_field.data();
-
   this->diff_input.setTimestep(dt);
 
+  double *field = curr_field.data();
   for (int i = 0; i < this->prop_count; i++) {
-    auto *in_field = &field[i * this->n_cells_per_prop];
+    double *in_field = &field[i * this->n_cells_per_prop];
     std::vector<double> in_alpha(this->n_cells_per_prop, this->alpha[i]);
     this->diff_input.setBoundaryCondition(this->bc_vec[i]);
 
@@ -195,6 +200,8 @@ void DiffusionModule::simulate(double dt) {
     }
   }
 
+  std::cout << " done!\n";
+
   sim_a_transport = MPI_Wtime();
 
   transport_t += sim_a_transport - sim_b_transport;