MDL: suppressed warnings (hard coded for now); added 1k*2k simulation bench/dolo_diffu_inner_large

bench/dolo_diffu_inner/dolo_diffu_inner_large.R

@@ -0,0 +1,144 @@
+## Time-stamp: "Last modified 2023-01-10 13:51:40 delucia"
+
+database <- normalizePath("../share/poet/examples/phreeqc_kin.dat")
+input_script <- normalizePath("../share/poet/bench/dolo_inner.pqi")
+
+#################################################################
+##                          Section 1                          ##
+##                     Grid initialization                     ##
+#################################################################
+
+n <- 2000
+m <- 1000
+
+types <- c("scratch", "phreeqc", "rds")
+
+init_cell <- list(
+  "H" = 110.683,
+  "O" = 55.3413,
+  "Charge" = -5.0822e-19,
+  "C" = 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],
+  init_cell = as.data.frame(init_cell),
+  props = names(init_cell),
+  database = database,
+  input_script = input_script
+)
+
+
+##################################################################
+##                          Section 2                           ##
+##         Diffusion parameters and boundary conditions         ##
+##################################################################
+
+## initial conditions
+init_diffu <- c(
+  "H" = 110.683,
+  "O" = 55.3413,
+  "Charge" = -5.0822e-19,
+  "C" = 1.2279E-4,
+  "Ca" = 1.2279E-4,
+  "Cl" = 0,
+  "Mg" = 0
+)
+
+## diffusion coefficients
+alpha_diffu <- c(
+  "H" = 1E-6,
+  "O" = 1E-6,
+  "Charge" = 1E-6,
+  "C" = 1E-6,
+  "Ca" = 1E-6,
+  "Cl" = 1E-6,
+  "Mg" = 1E-6
+)
+
+## list of boundary conditions/inner nodes
+vecinj_diffu <- list(
+    list(
+        "H" = 110.683,
+        "O" = 55.3413,
+        "Charge" = 1.90431e-16,
+        "C" = 0,
+        "Ca" = 0,
+        "Cl" = 0.002,
+        "Mg" = 0.001
+    ),
+    list(
+        "H" = 110.683,
+        "O" = 55.3413,
+        "Charge" = 1.90431e-16,
+        "C" = 0,
+        "Ca" = 0.0,
+        "Cl" = 0.004,
+        "Mg" = 0.002
+    )
+)
+
+vecinj_inner <- list(
+  l1 = c(1,400,200),
+  l2 = c(2,1600,800),
+  l3 = c(2,1600,800)
+)
+
+boundary <- list(
+#  "N" = c(1, rep(0, n-1)),
+  "N" = rep(0, n),
+  "E" = rep(0, m),
+  "S" = rep(0, n),
+  "W" = rep(0, m)
+)
+
+diffu_list <- names(alpha_diffu)
+
+diffusion <- list(
+  init = init_diffu,
+  vecinj = do.call(rbind.data.frame, vecinj_diffu),
+  vecinj_inner = vecinj_inner,
+  vecinj_index = boundary,
+  alpha = alpha_diffu
+)
+
+#################################################################
+##                          Section 3                          ##
+##                  Chemistry module (Phreeqc)                 ##
+#################################################################
+
+
+## # Needed when using DHT
+signif_vector <- c(10, 10, 2, 5, 5, 5, 5, 0, 5, 5)
+prop_type <- c("", "", "", "act", "act", "act", "act", "ignore", "", "")
+prop <- names(init_cell)
+
+chemistry <- list(
+  database = database,
+  input_script = input_script
+)
+
+#################################################################
+##                          Section 4                          ##
+##              Putting all those things together              ##
+#################################################################
+
+
+iterations <- 10
+dt <- 200
+
+setup <- list(
+  grid = grid,
+  diffusion = diffusion,
+  chemistry = chemistry,
+  iterations = iterations,
+  timesteps = rep(dt, iterations)
+)

src/PhreeqcWrapper.cpp

@@ -66,7 +66,8 @@ void poet::PhreeqcWrapper::SetupAndLoadDB(
 
 void poet::PhreeqcWrapper::InitFromFile(const std::string &strInputFile) {
   this->RunFile(true, true, false, strInputFile);
-  this->RunString(true, false, true, "DELETE; -all");
+  // MDL: this is run only by the workers
+  this->RunString(true, false, true, "DELETE; -all; PRINT; -warnings 0;");
 
   this->FindComponents();
 

util/data_evaluation/RFun_Eval.R

@@ -1,6 +1,6 @@
 ## Simple library of functions to assess and visualize the results of the coupled simulations
 
-## Time-stamp: "Last modified 2022-12-16 17:38:12 delucia"
+## Time-stamp: "Last modified 2023-01-17 19:12:13 delucia"
 
 require(RedModRphree)
 require(Rmufits)  ## essentially for PlotCartCellData 
@@ -24,7 +24,7 @@ ReadRTSims <- function(dir) {
 }
 
 ## function which reads all successive DHT stored in a given directory
-ReadAllDHT <- function(dir, new_scheme = T) {
+ReadAllDHT <- function(dir, new_scheme = TRUE) {
     files_full <- list.files(dir, pattern="iter.*dht", full.names=TRUE)
     files_name <- list.files(dir, pattern="iter.*dht", full.names=FALSE)
     res <- lapply(files_full, ReadDHT, new_scheme = new_scheme)
@@ -33,7 +33,7 @@ ReadAllDHT <- function(dir, new_scheme = T) {
 }
 
 ## function which reads one .dht file and gives a matrix
-ReadDHT <- function(file, new_scheme = T) {
+ReadDHT <- function(file, new_scheme = TRUE) {
     conn <- file(file, "rb")  ## open for reading in binary mode
     if (!isSeekable(conn))
         stop("Connection not seekable")
@@ -91,27 +91,27 @@ PlotScatter <- function(sam1, sam2, which=NULL, labs=c("NO DHT", "DHT"), pch="."
 ##### Some metrics for relative comparison
 
 ## Using range as norm
-RranRMSE <- function (pred, obs)
+RranRMSE <- function(pred, obs)
     sqrt(mean((pred - obs)^2))/abs(max(pred) - min(pred))
 
 ## Using max val as norm
-RmaxRMSE <- function (pred, obs)
+RmaxRMSE <- function(pred, obs)
     sqrt(mean((pred - obs)^2)/abs(max(pred)))
 
 ## Using sd as norm
-RsdRMSE <- function (pred, obs)
+RsdRMSE <- function(pred, obs)
     sqrt(mean((pred - obs)^2))/sd(pred)
 
 ## Using mean as norm
-RmeanRMSE <- function (pred, obs)
+RmeanRMSE <- function(pred, obs)
     sqrt(mean((pred - obs)^2))/mean(pred)
 
 ## Using mean as norm
-RAEmax <- function (pred, obs)
+RAEmax <- function(pred, obs)
     mean(abs(pred - obs))/max(pred)
 
 ## Max absolute error
-MAE <- function (pred, obs)
+MAE <- function(pred, obs)
     max(abs(pred - obs))
 
 ## workhorse function for ComputeErrors and its use with mapply