prep advection

bench/dolo/dolo_adv.R

@@ -0,0 +1,191 @@
+## Time-stamp: "Last modified 2023-08-17 16:28:57 mluebke"
+
+database <- normalizePath("../share/poet/bench/dolo/phreeqc_kin.dat")
+input_script <- normalizePath("../share/poet/bench/dolo/dolo_inner.pqi")
+
+#################################################################
+##                          Section 1                          ##
+##                     Grid initialization                     ##
+#################################################################
+
+n <- 50
+m <- 50
+
+types <- c("scratch", "phreeqc", "rds")
+
+init_cell <- list(
+  "H" = 110.683,
+  "O" = 55.3413,
+  "Charge" = -5.0822e-19,
+  "C(4)" = 1.2279E-4,
+  "Ca" = 1.2279E-4,
+  "Cl" = 0,
+  "Mg" = 0,
+  "O2g" = 0.499957,
+  "Calcite" = 2.07e-4,
+  "Dolomite" = 0
+)
+
+grid <- list(
+  n_cells = c(n, m),
+  s_cells = c(1, 1),
+  type = types[1]
+)
+
+##################################################################
+##                          Section 2                           ##
+##                     Advection parameters                     ##
+##################################################################
+
+## initial conditions
+init_adv <- c(
+  "H" = 110.683,
+  "O" = 55.3413,
+  "Charge" = -5.0822e-19,
+  "C(4)" = 1.2279E-4,
+  "Ca" = 1.2279E-4,
+  "Cl" = 0,
+  "Mg" = 0
+)
+
+## list of boundary conditions/inner nodes
+vecinj_adv <- list(
+  list(
+    "H" = 110.683,
+    "O" = 55.3413,
+    "Charge" = 1.90431e-16,
+    "C(4)" = 0,
+    "Ca" = 0,
+    "Cl" = 0.002,
+    "Mg" = 0.001
+  )
+)
+
+vecinj_inner <- list(
+  l1 = c(1, 1, 1)
+)
+
+# Create a list to store grid cell information
+grid_list <- vector("list", n * m)
+
+# Function to get the index of a grid cell given its row and column
+get_index <- function(row, col) {
+  return((row - 1) * m + col)
+}
+
+# Loop through each row and column to populate the grid_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()
+
+    # 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))
+    }
+
+    # Store the connections in the grid_list
+    grid_list[[index]] <- connections
+  }
+}
+
+vecinj <- do.call(rbind.data.frame, vecinj_adv)
+names(vecinj) <- names(init_adv)
+
+advection <- list(
+  init = init_adv,
+  vecinj = vecinj,
+  vecinj_inner = vecinj_inner,
+  grid = grid_list
+)
+
+#################################################################
+##                          Section 3                          ##
+##                  Chemistry module (Phreeqc)                 ##
+#################################################################
+
+
+## # Needed when using DHT
+dht_species <- c(
+  "H" = 10,
+  "O" = 10,
+  "Charge" = 3,
+  "C(4)" = 5,
+  "Ca" = 5,
+  "Cl" = 5,
+  "Mg" = 5,
+  "Calcite" = 5,
+  "Dolomite" = 5
+)
+
+check_sign_cal_dol_dht <- function(old, new) {
+  if ((old["Calcite"] == 0) != (new["Calcite"] == 0)) {
+    return(TRUE)
+  }
+  if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) {
+    return(TRUE)
+  }
+  return(FALSE)
+}
+
+fuzz_input_dht_keys <- function(input) {
+  return(input[names(dht_species)])
+}
+
+check_sign_cal_dol_interp <- function(to_interp, data_set) {
+  data_set <- as.data.frame(do.call(rbind, data_set), check.names = FALSE, optional = TRUE)
+  names(data_set) <- names(dht_species)
+  cal <- (data_set$Calcite == 0) == (to_interp["Calcite"] == 0)
+  dol <- (data_set$Dolomite == 0) == (to_interp["Dolomite"] == 0)
+
+  cal_dol_same_sig <- cal == dol
+  return(rev(which(!cal_dol_same_sig)))
+}
+
+check_neg_cal_dol <- function(result) {
+  neg_sign <- (result["Calcite"] <- 0) || (result["Dolomite"] < 0)
+  return(any(neg_sign))
+}
+
+hooks <- list(
+  dht_fill = check_sign_cal_dol_dht,
+  dht_fuzz = fuzz_input_dht_keys,
+  interp_pre_func = check_sign_cal_dol_interp,
+  interp_post_func = check_neg_cal_dol
+)
+
+chemistry <- list(
+  database = database,
+  input_script = input_script,
+  dht_species = dht_species,
+  hooks = hooks
+)
+
+#################################################################
+##                          Section 4                          ##
+##              Putting all those things together              ##
+#################################################################
+
+
+iterations <- 10
+dt <- 200
+
+setup <- list(
+  grid = grid,
+  advection = advection,
+  chemistry = chemistry,
+  iterations = iterations,
+  timesteps = rep(dt, iterations),
+  store_result = TRUE
+)

src/CMakeLists.txt

@@ -11,6 +11,7 @@ add_library(poetlib
   Chemistry/SurrogateModels/ProximityHashTable.cpp
   DataStructures/Field.cpp
   Transport/DiffusionModule.cpp
+  Transport/AdvectionModule.cpp
 )
 
 target_link_libraries(poetlib PUBLIC

src/Transport/AdvectionModule.cpp

@@ -0,0 +1,50 @@
+#include "AdvectionModule.hpp"
+
+#include <cstddef>
+#include <cstdint>
+#include <string>
+#include <vector>
+
+#include <Rcpp.h>
+
+namespace poet {
+
+void simulate(double dt);
+
+void AdvectionModule::initializeParams(RInsidePOET &R) {
+  const std::uint32_t n = this->n_cells[0] =
+      R.parseEval("mysetup$grid$n_cells[1]");
+  const std::uint32_t m = this->n_cells[1] =
+      R.parseEval("mysetup$grid$n_cells[2]");
+
+  const std::uint32_t field_size{n * m};
+
+  const auto init_vec =
+      Rcpp::as<Rcpp::NumericVector>(R.parseEval("mysetup$advection$init"));
+
+  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());
+
+  for (const auto &val : init_vec) {
+    init_field.push_back(std::vector<double>(field_size, val));
+  }
+
+  const Rcpp::DataFrame vecinj = R.parseEval("mysetup$advection$vecinj");
+  const Rcpp::DataFrame inner_vecinj =
+      R.parseEval("mysetup$advection$vecinj_inner");
+
+  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);
+    for (std::size_t prop_i = 0; prop_i < prop_names.size(); prop_i++) {
+      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->t_field = Field(field_size, init_field, prop_names);
+}
+
+} // namespace poet

src/Transport/AdvectionModule.hpp

@@ -0,0 +1,90 @@
+/*
+** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
+** Potsdam)
+**
+** Copyright (C) 2018-2022 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
+** Foundation; either version 2 of the License, or (at your option) any later
+** version.
+**
+** POET is distributed in the hope that it will be useful, but WITHOUT ANY
+** WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+** A PARTICULAR PURPOSE. See the GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License along with
+** this program; if not, write to the Free Software Foundation, Inc., 51
+** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+*/
+
+#ifndef ADVECTION_MODULE_H
+#define ADVECTION_MODULE_H
+
+#include "../Base/RInsidePOET.hpp"
+#include "../DataStructures/DataStructures.hpp"
+
+#include <Rcpp.h>
+#include <array>
+#include <cstdint>
+#include <vector>
+
+namespace poet {
+/**
+ * @brief Class describing transport simulation using advection
+ *
+ * Quick and dirty implementation, assuming homogenous and constant porisity,
+ * pressure and temperature. Open boundary conditions are assumed.
+ *
+ */
+
+class AdvectionModule {
+public:
+  /**
+   * @brief Construct a new TransportSim object
+   *
+   * The instance will only be initialized with given R object.
+   *
+   * @param R RRuntime object
+   */
+  AdvectionModule(RInsidePOET &R) { initializeParams(R); }
+
+  /**
+   * @brief Run simulation for one iteration
+   *
+   * This will simply call the R function 'master_advection'
+   *
+   */
+  void simulate(double dt);
+
+  /**
+   * @brief Get the transport time
+   *
+   * @return double time spent in transport
+   */
+  double getTransportTime() const { return this->transport_t; }
+
+  /**
+   * \brief Returns the current diffusion field.
+   *
+   * \return Reference to the diffusion field.
+   */
+  const Field &getField() const { return this->t_field; }
+
+private:
+  void initializeParams(RInsidePOET &R);
+
+  std::array<std::uint32_t, 2> n_cells;
+  std::set<std::uint32_t> cells_const;
+
+  Field t_field;
+
+  /**
+   * @brief time spent for transport
+   *
+   */
+  double transport_t = 0.f;
+};
+} // namespace poet
+
+#endif // ADVECTION_MODULE_H

test/CMakeLists.txt

@@ -7,6 +7,7 @@ target_link_libraries(testPOET doctest poetlib)
 target_include_directories(testPOET PRIVATE "${PROJECT_SOURCE_DIR}/src")
 
 get_filename_component(TEST_RInsideSourceFile "RInsidePOET_funcs.R" REALPATH)
+get_filename_component(TEST_SampleInputScript "${PROJECT_SOURCE_DIR}/bench/dolo/dolo_adv.R" REALPATH)
 configure_file(testDataStructures.hpp.in testDataStructures.hpp)
 target_include_directories(testPOET PRIVATE "${CMAKE_CURRENT_BINARY_DIR}")
 

test/testAdvection.cpp

@@ -0,0 +1,15 @@
+#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);
+}

test/testDataStructures.hpp.in

@@ -2,5 +2,6 @@
 #define TESTRINSIDEPOET_H_
 
 inline const char *RInside_source_file = "@TEST_RInsideSourceFile@";
+inline const char *SampleInputScript = "@TEST_SampleInputScript@";
 
 #endif // TESTRINSIDEPOET_H_