MDL: added barite_50ai bench

bench/barite/barite_50ai.R

@@ -0,0 +1,60 @@
+## Time-stamp: "Last modified 2024-05-29 10:51:35 delucia"
+cols <- 50
+rows <- 50
+
+s_cols <- 0.25
+s_rows <- 0.25
+
+grid_def <- matrix(2, nrow = rows, ncol = cols)
+
+# Define grid configuration for POET model
+grid_setup <- list(
+  pqc_in_file = "./barite.pqi",
+  pqc_db_file = "./db_barite.dat", # Path to the database file for Phreeqc
+  grid_def = grid_def, # Definition of the grid, containing IDs according to the Phreeqc input script
+  grid_size = c(s_rows, s_cols), # Size of the grid in meters
+  constant_cells = c() # IDs of cells with constant concentration
+)
+
+bound_length <- 2
+
+bound_def <- list(
+  "type" = rep("constant", bound_length),
+  "sol_id" = rep(3, bound_length),
+  "cell" = seq(1, bound_length)
+)
+
+homogenous_alpha <- 1e-8
+
+diffusion_setup <- list(
+  boundaries = list(
+    "W" = bound_def,
+    "N" = bound_def
+  ),
+  alpha_x = homogenous_alpha,
+  alpha_y = homogenous_alpha
+)
+
+dht_species <- c(
+  "H"         = 4,
+  "O"         = 9,
+  "Charge"    = 4,
+  "Ba"        = 4,
+  "Cl"        = 4,
+  "S(6)"      = 7,
+  "Sr"        = 4,
+  "Barite"    = 7,
+  "Celestite" = 7
+)
+
+chemistry_setup <- list(
+  dht_species = dht_species,
+  ai_surrogate_input_script = "./barite_50ai_surr_mdl.R"
+)
+
+# Define a setup list for simulation configuration
+setup <- list(
+  Grid = grid_setup, # Parameters related to the grid structure
+  Diffusion = diffusion_setup, # Parameters related to the diffusion process
+  Chemistry = chemistry_setup
+)

bench/barite/barite_50ai_rt.R

@@ -0,0 +1,9 @@
+iterations <- 1000
+
+dt <- 200
+
+list(
+    timesteps = rep(dt, iterations),
+    store_result = TRUE,
+    out_save = c(1, 5, seq(20, iterations, by=20))
+)

bench/barite/barite_50ai_surr_mdl.R

@@ -0,0 +1,85 @@
+## Time-stamp: "Last modified 2024-05-30 11:16:57 delucia"
+
+## load a pretrained model from tensorflow file
+## Use the global variable "ai_surrogate_base_path" when using file paths
+## relative to the input script
+initiate_model <- function() {
+    require(keras3)
+    require(tensorflow)
+    init_model <- normalizePath(paste0(ai_surrogate_base_path,
+                                       "barite_50ai_all.keras"))
+    Model <- keras3::load_model(init_model)
+    msgm("Loaded model:")
+    print(str(Model))
+    return(Model)
+}
+
+scale_min_max <- function(x, min, max, backtransform) {
+    if (backtransform) {
+        return((x * (max - min)) + min)
+    } else {
+        return((x - min) / (max - min))
+    }
+}
+
+minmax <- list(min = c(H = 111.012433592824, O = 55.5062185549492, Charge = -3.1028354471876e-08, 
+                       Ba = 1.87312878574393e-141, Cl = 0, `S(6)` = 4.24227510643685e-07, 
+                       Sr = 0.00049382996130541, Barite = 0.000999542409828586, Celestite = 0.244801877115968),
+               max = c(H = 111.012433679682, O = 55.5087003521685, Charge = 5.27666636082035e-07, 
+                       Ba = 0.0908849779513762, Cl = 0.195697626449355, `S(6)` = 0.000620774752665846, 
+                       Sr = 0.0558680070692722, Barite = 0.756779139057097, Celestite = 1.00075422160624
+                       ))
+
+preprocess <- function(df) {
+    if (!is.data.frame(df))
+        df <- as.data.frame(df, check.names = FALSE)
+    
+    as.data.frame(lapply(colnames(df),
+                         function(x) scale_min_max(x=df[x],
+                                                   min=minmax$min[x],
+                                                   max=minmax$max[x],
+                                                   backtransform=FALSE)),
+                  check.names = FALSE)
+}
+
+postprocess <- function(df) {
+    if (!is.data.frame(df))
+        df <- as.data.frame(df, check.names = FALSE)
+    
+    as.data.frame(lapply(colnames(df),
+                         function(x) scale_min_max(x=df[x],
+                                                   min=minmax$min[x],
+                                                   max=minmax$max[x],
+                                                   backtransform=TRUE)),
+                  check.names = FALSE)
+}
+
+mass_balance <- function(predictors, prediction) {
+    dBa <- abs(prediction$Ba + prediction$Barite -
+               predictors$Ba - predictors$Barite)
+    dSr <- abs(prediction$Sr + prediction$Celestite -
+               predictors$Sr - predictors$Celestite)
+    return(dBa + dSr)
+}
+
+validate_predictions <- function(predictors, prediction) {
+    epsilon <- 1E-7
+    mb <- mass_balance(predictors, prediction)
+    msgm("Mass balance mean:", mean(mb))
+    msgm("Mass balance variance:", var(mb))
+    ret <- mb < epsilon
+    msgm("Rows where mass balance meets threshold", epsilon, ":",
+         sum(ret))
+    return(ret)
+}
+
+training_step <- function(model, predictor, target, validity) {
+  msgm("Training:")
+
+  x <- as.matrix(predictor)
+  y <- as.matrix(target[colnames(x)])
+
+  model %>% keras3::fit(x, y)
+
+  model %>% keras3::save_model(paste0(out_dir, "/current_model.keras"), overwrite=TRUE)
+}