First advection with fixed timestep and volume

bench/dolo/dolo_adv.R

@@ -1,4 +1,4 @@
-## Time-stamp: "Last modified 2023-08-17 16:28:57 mluebke"
+## Time-stamp: "Last modified 2023-08-18 13:13:46 mluebke"
 
 database <- normalizePath("../share/poet/bench/dolo/phreeqc_kin.dat")
 input_script <- normalizePath("../share/poet/bench/dolo/dolo_inner.pqi")
@@ -8,8 +8,8 @@ input_script <- normalizePath("../share/poet/bench/dolo/dolo_inner.pqi")
 ##                     Grid initialization                     ##
 #################################################################
 
-n <- 50
-m <- 50
+n <- 5
+m <- 5
 
 types <- c("scratch", "phreeqc", "rds")
 
@@ -66,37 +66,42 @@ vecinj_inner <- list(
 )
 
 # Create a list to store grid cell information
-grid_list <- vector("list", n * m)
+flux_list <- list()
 
 # Function to get the index of a grid cell given its row and column
 get_index <- function(row, col) {
-  return((row - 1) * m + col)
+  index <- (row - 1) * m + col
+  if (index < 1) {
+    index <- -1
+  } else if (index > n * m) {
+    index <- -1
+  }
+  return(index)
 }
 
-# Loop through each row and column to populate the grid_list
+# Loop through each row and column to populate the flux_list
 for (row in 1:n) {
   for (col in 1:m) {
     index <- get_index(row, col)
 
     # Initialize the connections for the current cell
-    connections <- c()
+    flux <- c()
 
     # Check and add connections to the east, south, west, and north cells
-    if (col < m) {
-      connections <- c(connections, get_index(row, col + 1))
-    }
-    if (row < n) {
-      connections <- c(connections, get_index(row + 1, col))
-    }
-    if (col > index) {
-      connections <- c(connections, get_index(row, col - 1))
-    }
-    if (row > index) {
-      connections <- c(connections, get_index(row - 1, col))
-    }
+    # east
+    flux <- c(flux, -0.1)
 
-    # Store the connections in the grid_list
-    grid_list[[index]] <- connections
+    # south
+    flux <- c(flux, -0.1)
+
+    # west
+    flux <- c(flux, 0.1)
+
+    # north
+    flux <- c(flux, 0.1)
+
+    # Store the connections in the flux_list
+    flux_list[[index]] <- flux
   }
 }
 
@@ -107,7 +112,7 @@ advection <- list(
   init = init_adv,
   vecinj = vecinj,
   vecinj_inner = vecinj_inner,
-  grid = grid_list
+  const_flux = flux_list
 )
 
 #################################################################

src/Transport/AdvectionModule.cpp

@@ -9,7 +9,82 @@
 
 namespace poet {
 
-void simulate(double dt);
+// inline std::vector<std::vector<double>> getFlux(const Rcpp::List &flux_list)
+// {
+//   std::vector<std::vector<double>> fluxes(flux_lisVt.size());
+//   for (std::size_t i = 0; i < flux_list.size(); i++) {
+//     const auto flux_2d = Rcpp::as<std::vector<double>>(flux_list[i + 1]);
+//     fluxes[i] = flux_2d;
+//   }
+
+//   return fluxes;
+// }
+
+inline std::array<std::int32_t, 4> getIndices(std::int32_t index,
+                                              std::int32_t n, std::int32_t m) {
+  std::array<std::int32_t, 4> indices;
+
+  // east index
+  indices[0] = (index % n == n - 1 ? -1 : index + 1);
+  // south index
+  indices[1] = (index + n >= m * n ? -1 : index + n);
+  // west index
+  indices[2] = (index % n == 0 ? -1 : index - 1);
+  // north index
+  indices[3] = (index - n < 0 ? -1 : index - n);
+
+  return indices;
+}
+
+inline double getFluxApplyConc(bool inbound, std::uint32_t curr_index,
+                               std::int32_t neighbor_index,
+                               const std::vector<double> &conc) {
+  if (inbound && neighbor_index >= 0) {
+    return conc[neighbor_index];
+  }
+
+  // outbound flow only affects current cell with given index
+  return conc[curr_index];
+}
+
+void AdvectionModule::simulate(double dt) {
+  const auto &n = this->n_cells[0];
+  const auto &m = this->n_cells[1];
+
+  // HACK: constant flux for this moment imported from R runtime
+
+  RInsidePOET &R = RInsidePOET::getInstance();
+
+  const auto flux_list =
+      Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$advection$const_flux"));
+  std::vector<std::vector<double>> flux(flux_list.size());
+  for (std::size_t i = 0; i < flux_list.size(); i++) {
+    const auto flux_2d = Rcpp::as<std::vector<double>>(flux_list[i]);
+    flux[i] = flux_2d;
+  }
+
+  auto field_vec = t_field.As2DVector();
+  for (auto &species : field_vec) {
+    std::vector<double> spec_copy = species;
+    for (std::size_t cell_i = 0; cell_i < n * m; cell_i++) {
+      if (this->inactive_cells.find(cell_i) != this->inactive_cells.end()) {
+        continue;
+      }
+      const auto indices = getIndices(cell_i, n, m);
+      double ds{.0};
+      for (std::size_t neighbor_i = 0; neighbor_i < indices.size();
+           neighbor_i++) {
+        const double &curr_flux = flux[cell_i][neighbor_i];
+        const double flux_apply_val = getFluxApplyConc(
+            curr_flux > 0, cell_i, indices[neighbor_i], species);
+        ds += curr_flux * flux_apply_val;
+      }
+      spec_copy[cell_i] += ds;
+    }
+    species = spec_copy;
+  }
+  t_field = field_vec;
+}
 
 void AdvectionModule::initializeParams(RInsidePOET &R) {
   const std::uint32_t n = this->n_cells[0] =
@@ -22,6 +97,7 @@ void AdvectionModule::initializeParams(RInsidePOET &R) {
   const auto init_vec =
       Rcpp::as<Rcpp::NumericVector>(R.parseEval("mysetup$advection$init"));
 
+  // Initialize prop names and according values
   const auto prop_names = Rcpp::as<std::vector<std::string>>(init_vec.names());
   std::vector<std::vector<double>> init_field;
   init_field.reserve(prop_names.size());
@@ -30,10 +106,12 @@ void AdvectionModule::initializeParams(RInsidePOET &R) {
     init_field.push_back(std::vector<double>(field_size, val));
   }
 
+  // Set inner constant cells. Is it needed?
   const Rcpp::DataFrame vecinj = R.parseEval("mysetup$advection$vecinj");
   const Rcpp::DataFrame inner_vecinj =
       R.parseEval("mysetup$advection$vecinj_inner");
 
+  // same logic as for diffusion module applied
   for (std::size_t i = 0; i < inner_vecinj.size(); i++) {
     const Rcpp::NumericVector tuple = inner_vecinj[i];
     const std::uint32_t cell_index = (tuple[2] - 1) + ((tuple[1] - 1) * n);
@@ -41,7 +119,7 @@ void AdvectionModule::initializeParams(RInsidePOET &R) {
       const Rcpp::NumericVector curr_prop_vec_inj = vecinj[prop_names[prop_i]];
       init_field[prop_i][cell_index] = curr_prop_vec_inj[tuple[0] - 1];
     }
-    this->cells_const.insert(cell_index);
+    this->inactive_cells.insert(cell_index);
   }
 
   this->t_field = Field(field_size, init_field, prop_names);

src/Transport/AdvectionModule.hpp

@@ -27,6 +27,8 @@
 #include <Rcpp.h>
 #include <array>
 #include <cstdint>
+#include <map>
+#include <utility>
 #include <vector>
 
 namespace poet {
@@ -75,7 +77,7 @@ private:
   void initializeParams(RInsidePOET &R);
 
   std::array<std::uint32_t, 2> n_cells;
-  std::set<std::uint32_t> cells_const;
+  std::set<std::uint32_t> inactive_cells;
 
   Field t_field;
 

test/CMakeLists.txt

@@ -15,3 +15,7 @@ add_custom_target(check
     COMMAND $<TARGET_FILE:testPOET>
     DEPENDS testPOET
 )
+
+add_executable(advection advection/testAdvection.cpp)
+target_include_directories(advection PRIVATE "${CMAKE_CURRENT_BINARY_DIR}")
+target_link_libraries(advection poet_lib)

test/advection/testAdvection.cpp

@@ -0,0 +1,26 @@
+#include "testDataStructures.hpp"
+#include <cstddef>
+#include <poet/AdvectionModule.hpp>
+
+using namespace poet;
+
+constexpr std::size_t MAX_ITER = 100;
+
+int main(int argc, char **argv) {
+  auto &R = RInsidePOET::getInstance();
+
+  R["input_script"] = SampleInputScript;
+  R.parseEvalQ("source(input_script)");
+  R.parseEval("mysetup <- setup");
+
+  AdvectionModule adv(R);
+  R["field"] = adv.getField().asSEXP();
+  R.parseEval("saveRDS(field, 'adv_0.rds')");
+
+  for (std::size_t i = 0; i < MAX_ITER; i++) {
+    adv.simulate(1);
+  }
+
+  R["field"] = adv.getField().asSEXP();
+  R.parseEval("saveRDS(field, 'adv_10.rds')");
+}

test/testAdvection.cpp

@@ -1,15 +0,0 @@
-#include "testDataStructures.hpp"
-#include <poet/AdvectionModule.hpp>
-#include <doctest/doctest.h>
-
-using namespace poet;
-
-TEST_CASE("Advection Module") {
-    auto &R = RInsidePOET::getInstance();
-
-    R["input_script"] = SampleInputScript;
-    R.parseEvalQ("source(input_script)");
-    R.parseEval("mysetup <- setup");
-
-    AdvectionModule adv(R);
-}