Implement advection into POET

bench/dolo/dolo_adv.R

@@ -1,4 +1,4 @@
-## Time-stamp: "Last modified 2023-08-18 16:13:13 mluebke"
+## Time-stamp: "Last modified 2023-08-21 12:08:54 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 <- 5
-m <- 5
+n <- 100
+m <- 100
 
 types <- c("scratch", "phreeqc", "rds")
 
@@ -28,7 +28,7 @@ init_cell <- list(
 
 grid <- list(
   n_cells = c(n, m),
-  s_cells = c(n, m),
+  s_cells = c(n*10, m*10),
   type = types[1]
 )
 
@@ -38,12 +38,15 @@ grid <- list(
 ##################################################################
 
 ## initial conditions
+
+## HACK: We need the chemical initialization here, as chem module initialization
+## depends on transport until now. This will change in the future.
 init_adv <- c(
   "H" = 110.124,
   "O" = 55.0622,
-  "Charge" = -1.217e-09,
-  "C(4)" = 0,
-  "Ca" = 0,
+  "Charge" = -1.216307659761E-09,
+  "C(4)" = 1.230067028174E-04,
+  "Ca" = 1.230067028174E-04,
   "Cl" = 0,
   "Mg" = 0
 )
@@ -88,6 +91,8 @@ get_index <- function(row, col) {
   return(index)
 }
 
+flux_val <- 0.1
+
 # Loop through each row and column to populate the flux_list
 for (row in 1:n) {
   for (col in 1:m) {
@@ -98,16 +103,16 @@ for (row in 1:n) {
 
     # Check and add connections to the east, south, west, and north cells
     # east
-    flux <- c(flux, -1)
+    flux <- c(flux, -flux_val)
 
     # south
-    flux <- c(flux, -1)
+    flux <- c(flux, -flux_val)
 
     # west
-    flux <- c(flux, 1)
+    flux <- c(flux, flux_val)
 
     # north
-    flux <- c(flux, 1)
+    flux <- c(flux, flux_val)
 
     # Store the connections in the flux_list
     flux_list[[index]] <- flux
@@ -193,7 +198,7 @@ chemistry <- list(
 
 
 iterations <- 10
-dt <- 200
+dt <- 500
 
 setup <- list(
   grid = grid,

src/Transport/AdvectionModule.cpp

@@ -4,6 +4,7 @@
 #include <cstdint>
 #include <cstdlib>
 #include <limits>
+#include <mpi.h>
 #include <string>
 #include <vector>
 
@@ -90,6 +91,8 @@ AdvectionModule::CFLTimeVec(double req_dt,
 }
 
 void AdvectionModule::simulate(double dt) {
+  double start_t = MPI_Wtime();
+
   const auto &n = this->n_cells[0];
   const auto &m = this->n_cells[1];
 
@@ -133,7 +136,7 @@ void AdvectionModule::simulate(double dt) {
   for (std::size_t species_i = 0; species_i < field_vec.size(); species_i++) {
     auto &species = field_vec[species_i];
     std::vector<double> spec_copy = species;
-    for (const double &dt : time_vec) {
+    for (const double &curr_dt : time_vec) {
       for (std::size_t cell_i = 0; cell_i < n * m; cell_i++) {
 
         // if inactive cell -> skip
@@ -146,13 +149,15 @@ void AdvectionModule::simulate(double dt) {
             calcDeltaConc(cell_i, this->boundary_condition[species_i], species,
                           flux[species_i]);
         spec_copy[cell_i] +=
-            (dt * delta_conc) / (porosity[cell_i] * cell_volume);
+            (curr_dt * delta_conc) / (porosity[cell_i] * cell_volume);
       }
       species = spec_copy;
     }
   }
 
   t_field = field_vec;
+
+  this->transport_t += MPI_Wtime() - start_t;
 }
 
 void AdvectionModule::initializeParams(RInsidePOET &R) {

src/Transport/AdvectionModule.hpp

@@ -71,7 +71,7 @@ public:
    *
    * \return Reference to the diffusion field.
    */
-  const Field &getField() const { return this->t_field; }
+  Field &getField() { return this->t_field; }
 
 private:
   void initializeParams(RInsidePOET &R);

src/poet.cpp

@@ -114,11 +114,14 @@ void set_chem_parameters(poet::ChemistryModule &chem, uint32_t wp_size,
   chem.LoadDatabase(database_path);
 }
 
-inline double RunMasterLoop(SimParams &params, RInside &R,
+inline double RunMasterLoop(SimParams &params, RInsidePOET &R,
                             const GridParams &g_params, uint32_t nxyz_master) {
 
-  DiffusionParams d_params{R};
-  DiffusionModule diffusion(d_params, g_params);
+  // DiffusionParams d_params{R};
+  // DiffusionModule diffusion(d_params, g_params);
+
+  AdvectionModule advection(R);
+
   /* Iteration Count is dynamic, retrieving value from R (is only needed by
    * master for the following loop) */
   uint32_t maxiter = R.parseEval("mysetup$iterations");
@@ -131,8 +134,7 @@ inline double RunMasterLoop(SimParams &params, RInside &R,
   set_chem_parameters(chem, nxyz_master, params.getChemParams().database_path);
   chem.RunInitFile(params.getChemParams().input_script);
 
-  poet::ChemistryModule::SingleCMap init_df = DFToHashMap(d_params.initial_t);
-  chem.initializeField(diffusion.getField());
+  chem.initializeField(advection.getField());
 
   if (params.getNumParams().print_progressbar) {
     chem.setProgressBarPrintout(true);
@@ -161,17 +163,17 @@ inline double RunMasterLoop(SimParams &params, RInside &R,
 
     /* run transport */
     // TODO: transport to diffusion
-    diffusion.simulate(dt);
+    advection.simulate(dt);
 
-    chem.getField().update(diffusion.getField());
+    chem.getField().update(advection.getField());
 
     MSG("Chemistry step");
 
     chem.SetTimeStep(dt);
     chem.RunCells();
 
-    writeFieldsToR(R, diffusion.getField(), chem.GetField());
-    diffusion.getField().update(chem.GetField());
+    writeFieldsToR(R, advection.getField(), chem.GetField());
+    advection.getField().update(chem.GetField());
 
     R["req_dt"] = dt;
     R["simtime"] = (sim_time += dt);
@@ -214,7 +216,7 @@ inline double RunMasterLoop(SimParams &params, RInside &R,
   R["phreeqc_time"] = Rcpp::wrap(chem.GetWorkerPhreeqcTimings());
   R.parseEvalQ("profiling$phreeqc <- phreeqc_time");
 
-  R["simtime_transport"] = diffusion.getTransportTime();
+  R["simtime_transport"] = advection.getTransportTime();
   R.parseEvalQ("profiling$simtime_transport <- simtime_transport");
 
   // R["phreeqc_count"] = phreeqc_counts;
@@ -247,7 +249,6 @@ inline double RunMasterLoop(SimParams &params, RInside &R,
   }
 
   chem.MasterLoopBreak();
-  diffusion.end();
 
   return dSimTime;
 }

test/advection/testAdvection.cpp

@@ -6,6 +6,7 @@
 using namespace poet;
 
 constexpr std::size_t MAX_ITER = 10;
+constexpr double DT = 200;
 
 int main(int argc, char **argv) {
   auto &R = RInsidePOET::getInstance();
@@ -19,7 +20,7 @@ int main(int argc, char **argv) {
   R.parseEval("saveRDS(field, 'adv_0.rds')");
 
   for (std::size_t i = 0; i < MAX_ITER; i++) {
-    adv.simulate(1);
+    adv.simulate(DT);
     const std::string save_q =
         "saveRDS(field, 'adv_" + std::to_string(i + 1) + ".rds')";
     R["field"] = adv.getField().asSEXP();