Compare commits

...

2 Commits

Author SHA1 Message Date
Hannes Signer
b9fe09d937 remove non relevant benchmarks 2025-12-05 00:24:10 +01:00
Hannes Martin Signer
698621097f add conda environment 2025-12-05 00:01:20 +01:00
35 changed files with 149 additions and 4951 deletions

View File

@ -40,4 +40,3 @@ add_custom_target(${BENCHTARGET} ALL)
add_subdirectory(barite)
add_subdirectory(dolo)
add_subdirectory(surfex)

View File

@ -1,12 +1,10 @@
# Create a list of files
set(bench_files
barite_200.R
barite_het.R
)
set(runtime_files
barite_200_rt.R
barite_het_rt.R
)
# add_custom_target(barite_bench DEPENDS ${bench_files} ${runtime_files})

View File

@ -47,8 +47,7 @@ dht_species <- c(
)
chemistry_setup <- list(
dht_species = dht_species,
ai_surrogate_input_script = "./barite_200ai_surrogate_input_script.R"
dht_species = dht_species
)
# Define a setup list for simulation configuration

View File

@ -1,48 +0,0 @@
## 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() {
init_model <- normalizePath(paste0(ai_surrogate_base_path,
"model_min_max_float64.keras"))
return(load_model_tf(init_model))
}
scale_min_max <- function(x, min, max, backtransform) {
if (backtransform) {
return((x * (max - min)) + min)
} else {
return((x - min) / (max - min))
}
}
preprocess <- function(df, backtransform = FALSE, outputs = FALSE) {
minmax_file <- normalizePath(paste0(ai_surrogate_base_path,
"min_max_bounds.rds"))
global_minmax <- readRDS(minmax_file)
for (column in colnames(df)) {
df[column] <- lapply(df[column],
scale_min_max,
global_minmax$min[column],
global_minmax$max[column],
backtransform)
}
return(df)
}
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 <- 3e-5
mb <- mass_balance(predictors, prediction)
msgm("Mass balance mean:", mean(mb))
msgm("Mass balance variance:", var(mb))
msgm("Rows where mass balance meets threshold", epsilon, ":",
sum(mb < epsilon))
return(mb < epsilon)
}

View File

@ -1,60 +0,0 @@
## Time-stamp: "Last modified 2024-05-30 13:34:14 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" = 3,
"O" = 3,
"Charge" = 3,
"Ba" = 6,
"Cl" = 6,
"S" = 6,
"Sr" = 6,
"Barite" = 5,
"Celestite" = 5
)
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
)

Binary file not shown.

View File

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

View File

@ -1,90 +0,0 @@
## Time-stamp: "Last modified 2024-05-30 13:27:06 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("Starting incremental training:")
## x <- as.matrix(predictor)
## y <- as.matrix(target[colnames(x)])
history <- model %>% keras3::fit(x = data.matrix(predictor),
y = data.matrix(target),
epochs = 10, verbose=1)
keras3::save_model(model,
filepath = paste0(out_dir, "/current_model.keras"),
overwrite=TRUE)
return(model)
}

View File

@ -1,32 +0,0 @@
grid_def <- matrix(c(2, 3), nrow = 2, ncol = 5)
# Define grid configuration for POET model
grid_setup <- list(
pqc_in_file = "./barite_het.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(ncol(grid_def), nrow(grid_def)), # Size of the grid in meters
constant_cells = c() # IDs of cells with constant concentration
)
diffusion_setup <- list(
boundaries = list(
"W" = list(
"type" = rep("constant", nrow(grid_def)),
"sol_id" = rep(4, nrow(grid_def)),
"cell" = seq_len(nrow(grid_def))
)
),
alpha_x = 1e-6,
alpha_y = matrix(runif(10, 1e-8, 1e-7),
nrow = nrow(grid_def),
ncol = ncol(grid_def)
)
)
# 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 = list()
)

View File

@ -1,80 +0,0 @@
## Initial: everywhere equilibrium with Celestite NB: The aqueous
## solution *resulting* from this calculation is to be used as initial
## state everywhere in the domain
SOLUTION 1
units mol/kgw
water 1
temperature 25
pH 7
pe 4
S(6) 1e-12
Sr 1e-12
Ba 1e-12
Cl 1e-12
PURE 1
Celestite 0.0 1
SAVE SOLUTION 2 # <- phreeqc keyword to store and later reuse these results
END
RUN_CELLS
-cells 1
COPY solution 1 2-3
## Here a 5x2 domain:
|---+---+---+---+---|
-> | 2 | 2 | 2 | 2 | 2 |
4 |---+---+---+---+---|
-> | 3 | 3 | 3 | 3 | 3 |
|---+---+---+---+---|
## East boundary: "injection" of solution 4. North, W, S boundaries: closed
## Here the two distinct zones: nr 2 with kinetics Celestite (initial
## amount is 0, is then allowed to precipitate) and nr 3 with kinetic
## Celestite and Barite (both initially > 0) where the actual
## replacement takes place
#USE SOLUTION 2 <- PHREEQC keyword to reuse the results from previous calculation
KINETICS 2
Celestite
-m 0 # Allowed to precipitate
-parms 10.0
-tol 1e-9
END
#USE SOLUTION 2 <- PHREEQC keyword to reuse the results from previous calculation
KINETICS 3
Barite
-m 0.001
-parms 50.
-tol 1e-9
Celestite
-m 1
-parms 10.0
-tol 1e-9
END
## A BaCl2 solution (nr 4) is "injected" from the left boundary:
SOLUTION 4
units mol/kgw
pH 7
water 1
temp 25
Ba 0.1
Cl 0.2
END
## NB: again, the *result* of the SOLUTION 4 script defines the
## concentration of all elements (+charge, tot H, tot O)
## Ideally, in the initial state SOLUTION 1 we should not have to
## define the 4 elemental concentrations (S(6), Sr, Ba and Cl) but
## obtain them having run once the scripts with the aqueous solution
## resulting from SOLUTION 1 once with KINETICS 2 and once with
## KINETICS 3.
RUN_CELLS
-cells 2-4

View File

@ -1,4 +0,0 @@
list(
timesteps = rep(50, 100),
store_result = TRUE
)

View File

@ -1,10 +1,8 @@
set(bench_files
dolo_inner_large.R
dolo_interp.R
)
set(runtime_files
dolo_inner_large_rt.R
dolo_interp_rt.R
)

Binary file not shown.

View File

@ -1,115 +0,0 @@
rows <- 2000
cols <- 1000
grid_def <- matrix(2, nrow = rows, ncol = cols)
# Define grid configuration for POET model
grid_setup <- list(
pqc_in_file = "./dol.pqi",
pqc_db_file = "./phreeqc_kin.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(cols, rows) / 100, # Size of the grid in meters
constant_cells = c() # IDs of cells with constant concentration
)
bound_size <- 2
diffusion_setup <- list(
inner_boundaries = list(
"row" = c(400, 1400, 1600),
"col" = c(200, 800, 800),
"sol_id" = c(3, 4, 4)
),
alpha_x = 1e-6,
alpha_y = 1e-6
)
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) {
dht_species <- c(
"H" = 3,
"O" = 3,
"Charge" = 3,
"C(4)" = 6,
"Ca" = 6,
"Cl" = 3,
"Mg" = 5,
"Calcite" = 4,
"Dolomite" = 4
)
return(input[names(dht_species)])
}
check_sign_cal_dol_interp <- function(to_interp, data_set) {
dht_species <- c(
"H" = 3,
"O" = 3,
"Charge" = 3,
"C(4)" = 6,
"Ca" = 6,
"Cl" = 3,
"Mg" = 5,
"Calcite" = 4,
"Dolomite" = 4
)
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(neg_sign)
}
# Optional when using Interpolation (example with less key species and custom
# significant digits)
pht_species <- c(
"C(4)" = 3,
"Ca" = 3,
"Mg" = 2,
"Calcite" = 2,
"Dolomite" = 2
)
chemistry_setup <- list(
dht_species = c(
"H" = 3,
"O" = 3,
"Charge" = 3,
"C(4)" = 6,
"Ca" = 6,
"Cl" = 3,
"Mg" = 5,
"Calcite" = 4,
"Dolomite" = 4
),
pht_species = pht_species,
hooks = list(
dht_fill = check_sign_cal_dol_dht,
dht_fuzz = fuzz_input_dht_keys,
interp_pre = check_sign_cal_dol_interp,
interp_post = check_neg_cal_dol
)
)
# 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 # Parameters related to the chemistry process
)

View File

@ -1,10 +0,0 @@
iterations <- 500
dt <- 50
out_save <- seq(5, iterations, by = 5)
list(
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = out_save
)

View File

@ -1,12 +0,0 @@
iterations <- 2000
dt <- 200
out_save <- c(1, 10, 20, seq(40, iterations, by = 40))
list(
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = out_save
)

View File

@ -1,90 +0,0 @@
## Time-stamp: "Last modified 2024-12-11 23:21:25 delucia"
library(PoetUtils)
library(viridis)
res <- ReadPOETSims("./res_fgcs2_96/")
pp <- PlotField(res$iter_200$C$Barite, rows = 200, cols = 200, contour = FALSE,
nlevels=12, palette=terrain.colors)
cairo_pdf("fgcs_Celestite_init.pdf", family="serif")
par(mar=c(0,0,0,0))
pp <- PlotField((res$iter_000$Celestite), rows = 200, cols = 200,
contour = FALSE, breaks=c(-0.5,0.5,1.5),
palette = grey.colors, plot.axes = FALSE, scale = FALSE,
main="Initial Celestite crystals")
dev.off()
cairo_pdf("fgcs_Ba_init.pdf", family="serif")
par(mar=c(0,0,0,0))
pp <- PlotField(log10(res$iter_001$C$Cl), rows = 200, cols = 200,
contour = FALSE,
palette = terrain.colors, plot.axes = FALSE, scale = FALSE,
main="log10(Ba)")
dev.off()
pp <- PlotField(log10(res$iter_002$C$Ba), rows = 200, cols = 200,
contour = FALSE, palette = viridis, rev.palette = FALSE,
main = "log10(Ba) after 5 iterations")
pp <- PlotField(log10(res$iter_200$C$`S(6)`), rows = 200, cols = 200, contour = FALSE)
str(res$iter_00)
res$iter_178$C$Barite
pp <- res$iter_043$C$Barite
breaks <- pretty(pp, n = 5)
br <- c(0, 0.0005, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1)
pp <- PlotField(res$iter_200$C$Barite, rows = 200, cols = 200, contour = FALSE,
breaks = br, palette=terrain.colors)
cairo_pdf("fgcs_Barite_200.pdf", family="serif")
pp <- PlotField(log10(res$iter_200$C$Barite), rows = 200, cols = 200,
contour = FALSE, palette = terrain.colors, plot.axes = FALSE,
rev.palette = FALSE, main = "log10(Barite) after 200 iter")
dev.off()
ref <- ReadPOETSims("./res_fgcs_2_ref")
rei <- ReadPOETSims("./res_fgcs_2_interp1/")
timref <- ReadRObj("./res_fgcs_2_ref/timings.qs")
timint <- ReadRObj("./res_fgcs_2_interp1/timings.qs")
timref
timint
wch <- c("H","O", "Ba", "Sr","Cl", "S(6)")
rf <- data.matrix(ref$iter_001$C[, wch])
r1 <- data.matrix(rei$iter_001$C[, wch])
r1[is.nan(r1)] <- NA
rf[is.nan(rf)] <- NA
cairo_pdf("fgcs_interp_1.pdf", family="serif", width = 10, height = 7)
PlotScatter(rf, r1, which = wch, labs = c("ref", "interp"), cols = 3, log="", las = 1, pch=4)
dev.off()
head(r1)
head(rf)
rf$O
r1$O

View File

@ -1,2 +0,0 @@
* Refer to the LaTeX file (and pdf) for more information

View File

@ -1,105 +0,0 @@
## Time-stamp: "Last modified 2024-12-11 16:08:11 delucia"
cols <- 1000
rows <- 1000
dim_cols <- 50
dim_rows <- 50
ncirc <- 20 ## number of crystals
rmax <- cols / 10 ## max radius (in nodes)
set.seed(22933)
centers <- cbind(sample(seq_len(cols), ncirc), sample(seq_len(rows), ncirc))
radii <- sample(seq_len(rmax), ncirc, replace = TRUE)
mi <- matrix(rep(seq_len(cols), rows), byrow = TRUE, nrow = rows)
mj <- matrix(rep(seq_len(cols), each = rows), byrow = TRUE, nrow = rows)
tmpl <- lapply(seq_len(ncirc), function(x) which((mi - centers[x, 1])^2 + (mj - centers[x, 2])^2 < radii[x]^2, arr.ind = TRUE))
inds <- do.call(rbind, tmpl)
grid <- matrix(1, nrow = rows, ncol = cols)
grid[inds] <- 2
alpha <- matrix(1e-5, ncol = cols, nrow = rows)
alpha[inds] <- 1e-7
## image(grid, asp=1)
## Define grid configuration for POET model
grid_setup <- list(
pqc_in_file = "./barite_fgcs_2.pqi",
pqc_db_file = "../barite/db_barite.dat", ## database file
grid_def = grid, ## grid definition, IDs according to the Phreeqc input
grid_size = c(dim_cols, dim_rows), ## grid size in meters
constant_cells = c() ## IDs of cells with constant concentration
)
bound_length <- cols / 10
bound_N <- list(
"type" = rep("constant", bound_length),
"sol_id" = rep(3, bound_length),
"cell" = seq(1, bound_length)
)
bound_W <- list(
"type" = rep("constant", bound_length),
"sol_id" = rep(3, bound_length),
"cell" = seq(1, bound_length)
)
bound_E <- list(
"type" = rep("constant", bound_length),
"sol_id" = rep(4, bound_length),
"cell" = seq(rows - bound_length + 1, rows)
)
bound_S <- list(
"type" = rep("constant", bound_length),
"sol_id" = rep(4, bound_length),
"cell" = seq(cols - bound_length + 1, cols)
)
diffusion_setup <- list(
boundaries = list(
"W" = bound_W,
"N" = bound_N,
"E" = bound_E,
"S" = bound_S
),
alpha_x = alpha,
alpha_y = alpha
)
dht_species <- c(
"H" = 7,
"O" = 7,
"Ba" = 7,
"Cl" = 7,
"S" = 7,
"Sr" = 7,
"Barite" = 4,
"Celestite" = 4
)
pht_species <- c(
"Ba" = 4,
"Cl" = 3,
"S" = 3,
"Sr" = 3,
"Barite" = 0,
"Celestite" = 0
)
chemistry_setup <- list(
dht_species = dht_species,
pht_species = pht_species
)
## 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
)

View File

@ -1,49 +0,0 @@
SOLUTION 1
units mol/kgw
water 1
temperature 25
pH 7.008
pe 10.798
S 6.205e-04
Sr 6.205e-04
END
SOLUTION 2
units mol/kgw
water 1
temperature 25
pH 7.008
pe 10.798
S 6.205e-04
Sr 6.205e-04
KINETICS 2
Barite
-m 0.00
-parms 50. # reactive surface area
-tol 1e-9
Celestite
-m 1
-parms 10.0 # reactive surface area
-tol 1e-9
END
SOLUTION 3
units mol/kgw
water 1
temperature 25
Ba 0.1
Cl 0.2
END
SOLUTION 4
units mol/kgw
water 1
temperature 25
Ba 0.2
Cl 0.4
END
RUN_CELLS
-cells 1 2 3 4
END

View File

@ -1,7 +0,0 @@
iterations <- 200
dt <- 1000
list(
timesteps = rep(dt, iterations),
store_result = TRUE
)

View File

@ -1,20 +0,0 @@
set(bench_files
# surfex.R
# ex.R
PoetEGU_surfex_500.R
)
set(runtime_files
# surfex_rt.R
# ex_rt.R
PoetEGU_surfex_500_rt.R
)
ADD_BENCH_TARGET(
surfex_bench
bench_files
runtime_files
"surfex"
)
add_dependencies(${BENCHTARGET} surfex_bench)

View File

@ -1,63 +0,0 @@
## Time-stamp: "Last modified 2023-03-21 11:49:43 mluebke"
KNOBS
-logfile false
-iterations 10000
-convergence_tolerance 1E-12
-step_size 2
-pe_step_size 2
SELECTED_OUTPUT
-reset false
-high_precision true
-solution true
-state true
-step true
-pH true
-pe true
-ionic_strength true
-water true
SOLUTION 1
temp 13
units mol/kgw
pH 7.06355
pe -2.626517
C(4) 0.001990694
Ca 0.02172649
Cl 0.3227673 charge
Fe 0.0001434717
K 0.001902357
Mg 0.01739704
Na 0.2762882
S(6) 0.01652701
Sr 0.0004520361
U(4) 8.147792e-12
U(6) 2.237946e-09
-water 0.00133
EXCHANGE 1
-equil 1
Z 0.0012585
Y 0.0009418
END
SOLUTION 2
temp 13
units mol/kgw
C(-4) 2.92438561098248e-21
C(4) 2.65160558871092e-06
Ca 2.89001071336443e-05
Cl 0.000429291158114428
Fe(2) 1.90823391198114e-07
Fe(3) 3.10832423034763e-12
H(0) 2.7888235127385e-15
K 2.5301787e-06
Mg 2.31391999937907e-05
Na 0.00036746969
S(-2) 1.01376078438546e-14
S(2) 1.42247026981542e-19
S(4) 9.49422092568557e-18
S(6) 2.19812504654191e-05
Sr 6.01218519999999e-07
U(4) 4.82255946569383e-12
U(5) 5.49050615347901e-13
U(6) 1.32462838991902e-09
END

View File

@ -1,40 +0,0 @@
rows <- 500
cols <- 200
grid_left <- matrix(1, nrow = rows, ncol = cols/2)
grid_rght <- matrix(2, nrow = rows, ncol = cols/2)
grid_def <- cbind(grid_left, grid_rght)
# Define grid configuration for POET model
grid_setup <- list(
pqc_in_file = "./SurfexEGU.pqi",
pqc_db_file = "./SMILE_2021_11_01_TH.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(10, 4), # Size of the grid in meters
constant_cells = c() # IDs of cells with constant concentration
)
bound_def <- list(
"type" = rep("constant", cols),
"sol_id" = rep(3, cols),
"cell" = seq(1, cols)
)
diffusion_setup <- list(
boundaries = list(
"N" = bound_def
),
alpha_x = matrix(runif(rows*cols))*1e-8,
alpha_y = matrix(runif(rows*cols))*1e-9## ,1e-10
)
chemistry_setup <- list()
# 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 # Parameters related to the chemistry process
)

View File

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

View File

@ -1,100 +0,0 @@
#+TITLE: Description of =surfex= benchmark
#+AUTHOR: MDL <delucia@gfz-potsdam.de>
#+DATE: 2023-08-26
#+STARTUP: inlineimages
#+LATEX_CLASS_OPTIONS: [a4paper,9pt]
#+LATEX_HEADER: \usepackage{fullpage}
#+LATEX_HEADER: \usepackage{amsmath, systeme}
#+LATEX_HEADER: \usepackage{graphicx}
#+LATEX_HEADER: \usepackage{charter}
#+OPTIONS: toc:nil
* Quick start
#+begin_src sh :language sh :frame single
mpirun -np 4 ./poet ex.R ex_res
mpirun -np 4 ./poet surfex.R surfex_res
#+end_src
* List of Files
- =ex.R=: POET input script for a 100x100 simulation grid, only
exchange
- =ExBase.pqi=: PHREEQC input script for the =ex.R= model
- =surfex.R=: POET input script for a 1000x1000 simulation grid
considering both cation exchange and surface complexation
- =SurfExBase.pqi=: PHREEQC input script for the =surfex.R= model
- =SMILE_2021_11_01_TH.dat=: PHREEQC database containing the
parametrized data for Surface and Exchange, based on the SMILE
Thermodynamic Database (Version 01-November-2021)
* Chemical system
This model describes migration of Uranium radionuclide in Opalinus
clay subject to surface complexation and cation exchange on the
surface of clay minerals. These two processes account for the binding
of aqueous complexes to the surfaces of minerals, which may have a
significant impact on safety of underground nuclear waste repository.
Namely, they can act as retardation buffer for uranium complexes
entering into a natural system. The system is kindly provided by Dr.
T. Hennig and is inspired to the sandy facies BWS-A3 sample from the
Mont Terri underground lab (Hennig and Kühn, 2021).
This chemical system is highly redox-sensitive, and several elements
are defined in significant amounts in different valence states. In
total, 20 elemental concentrations and valences are transported:
C(-4), C(4), Ca, Cl, Fe(2), Fe(3), K, Mg, Na, S(-2), S(2), S(4), S(6),
Sr , U(4), U(5), U(6); plus the total H, total O and Charge implicitly
required by PHREEQC_RM.
** Exchange
The SMILE database defines thermodynamical data for exchange of all
major cations and uranyl-ions on Illite and Montmorillonite. In
PHREEQC terms:
- *Y* for Montmorillonite, with a total amount of 1.2585
milliequivalents and
- *Z* for Illite, with a total amount of 0.9418 meq
** Surface
Here we consider a Donnan diffuse double layer of 0.49 nm. Six
distinct sorption sites are defined:
- Kln_aOH (aluminol site) and Kln_siOH (silanol) for Kaolinite
- For Illite, strong and weak sites Ill_sOH and Ill_wOH respectively
- For Montmorillonite, strong and weak sites Mll_sOH and Mll_wOH
respectively
Refer to the =SurfExBase.pqi= script for the actual numerical values
of the parameters.
* POET simulations
** =ex.R=
This benchmark only considers EXCHANGE, no mineral or SURFACE
complexation is involved.
- Grid discretization: square domain of 1 \cdot 1 m^{2} discretized in
100x100 cells
- Boundary conditions: E, S and W sides of the domain are closed.
*Fixed concentrations* are fixed at the N boundary.
- Diffusion coefficients: isotropic homogeneous \alpha = 1E-06
- Time steps & iterations: 10 iterations with \Delta t of 200 s
- *DHT* is not implemented as of yet for models including SURFACE and
EXCHANGE geochemical processes *TODO*
- Hooks: no hooks defined *TODO*
** =surfex.R=
- Grid discretization: rectangular domain of 1 \cdot 1 m^{2}
discretized in 10 \times 10 cells
- Boundary conditions: E, S and W sides of the domain are closed.
*Fixed concentrations* are fixed at the N boundary.
- Diffusion coefficients: isotropic homogeneous \alpha = 1E-06
- Time steps & iterations: 10 iterations with \Delta t of 200 s
* References
- Hennig, T.; Kühn, M.Surrogate Model for Multi-Component Diffusion of
Uranium through Opalinus Clay on the Host Rock Scale. Appl. Sci.
2021, 11, 786. https://doi.org/10.3390/app11020786

File diff suppressed because it is too large Load Diff

View File

@ -1,80 +0,0 @@
## Time-stamp: "Last modified 2023-02-27 14:31:11 delucia"
KNOBS
-logfile false
-iterations 10000
-convergence_tolerance 1E-12
-step_size 2
-pe_step_size 2
SELECTED_OUTPUT
-reset false
-high_precision true
-solution true
-state true
-step true
-pH true
-pe true
-ionic_strength true
-water true
USER_PUNCH
-head total_o total_h cb C(-4) C(4) Ca Cl Fe(2) Fe(3) H(0) K Mg Na S(-2) S(2) S(4) S(6) Sr U(4) U(5) U(6) UO2(am,hyd) KdU
-start
5 w=TOT("water")
10 PUNCH TOTMOLE("O"), TOTMOLE("H"), CHARGE_BALANCE, w*TOT("C(-4)"), w*TOT("C(4)"), w*TOT("Ca"), w*TOT("Cl"), w*TOT("Fe(2)"), w*TOT("Fe(3)"), w*TOT("H(0)"), w*TOT("K"), w*TOT("Mg"), w*TOT("Na"), w*TOT("S(-2)"), w*TOT("S(2)"), w*TOT("S(4)"), w*TOT("S(6)"), w*TOT("Sr"), w*TOT("U(4)"), w*TOT("U(5)"), w*TOT("U(6)"), EQUI("UO2(am,hyd)")
20 PUNCH ((SURF("U, Ill")+SURF("U, Mll")+SURF("U, Kln")+EDL("U, Ill")+EDL("U, Mll")+EDL("U, Kln"))/((TOT("U")*1.01583)))/(0.002251896406*1000)
-end
SOLUTION 1
temp 13
units mol/kgw
pH 7.06355
pe -2.626517
C(4) 0.001990694
Ca 0.02172649
Cl 0.3227673 charge
Fe 0.0001434717
K 0.001902357
Mg 0.01739704
Na 0.2762882
S(6) 0.01652701
Sr 0.0004520361
U(4) 8.147792e-12
U(6) 2.237946e-09
-water 0.00133
SURFACE 1
-equil 1
-sites_units density
-donnan 4.9e-10
Kln_aOH 1.155 11. 5.0518
Kln_siOH 1.155
Ill_sOH 0.05 100. 5.5931
Ill_wOH 2.26
Mll_sOH 0.05 100. 1.0825
Mll_wOH 2.26
EXCHANGE 1
-equil 1
Z 0.0012585
Y 0.0009418
END
SOLUTION 2
temp 13
units mol/kgw
C(-4) 2.92438561098248e-21
C(4) 2.65160558871092e-06
Ca 2.89001071336443e-05
Cl 0.000429291158114428
Fe(2) 1.90823391198114e-07
Fe(3) 3.10832423034763e-12
H(0) 2.7888235127385e-15
K 2.5301787e-06
Mg 2.31391999937907e-05
Na 0.00036746969
S(-2) 1.01376078438546e-14
S(2) 1.42247026981542e-19
S(4) 9.49422092568557e-18
S(6) 2.19812504654191e-05
Sr 6.01218519999999e-07
U(4) 4.82255946569383e-12
U(5) 5.49050615347901e-13
U(6) 1.32462838991902e-09
END

View File

@ -1,108 +0,0 @@
## Time-stamp: "Last modified 2024-04-12 10:59:59 delucia"
## KNOBS
## -logfile false
## -iterations 10000
## -convergence_tolerance 1E-12
## -step_size 2
## -pe_step_size 2
SOLUTION 1 ## Porewater composition Opalinus Clay, WITHOUT radionuclides, AFTER EQUI_PHASES
pe -2.627 ## Eh = -227 mV, Value from Bossart & Thury (2008)-> PC borehole measurement 2003, Eh still decreasing
density 1.01583 ## kg/dm³ = g/cm³
temp 13 ## mean temperature Mont Terri, Bossart & Thury (2008), calculations performed for 25°C
units mol/kgw
## Mean composition
pH 7.064
Na 2.763e-01
Cl 3.228e-01 charge
S(6) 1.653e-02 as SO4
Ca 2.173e-02
Mg 1.740e-02
K 1.902e-03
Sr 4.520e-04
Fe 1.435e-04
U 2.247e-09
SURFACE 1 Opalinus Clay, clay minerals
## calculated with rho_b=2.2903 kg/dm³, poro=0.1662
## 1 dm³ = 13.565641 kg_sed/kg_pw
-equil 1 ## equilibrate with solution 1
-sites_units density ## set unit for binding site density to sites/nm2
-donnan 4.9e-10 ## calculated after Wigger & Van Loon (2018) for ionic strength after equilibration with minerales for pCO2=2.2 log10 bar
# surface density SSA (m2/g) mass (g/kgw)
Kln_aOH 1.155 11. 3798.4 ## Kaolinite 28 wt% (aluminol and silanol sites)
Kln_siOH 1.155
Ill_sOH 0.05 100. 4205.35 ## Illite 31 wt% (weak und strong binding sites)
Ill_wOH 2.26 ## 2 % strong binding sites
Mll_sOH 0.05 100. 813.94 ## Montmorillonite = smektite = 6 wt% (weak und strong binding sites)
Mll_wOH 2.26 ## 2 % strong binding sites
EXCHANGE 1 Exchanger, only illite+montmorillonite
## Illite = 0.225 eq/kg_rock, Montmorillonit = 0.87 eq/kg_rock
-equil 1 ## equilibrate with solution 1
Z 0.9462 ## = Illite
Y 0.70813 ## = Montmorillonite
END
SOLUTION 2 ## Porewater composition Opalinus Clay, WITHOUT radionuclides, AFTER EQUI_PHASES
pe -2.627 ## Eh = -227 mV, Value from Bossart & Thury (2008)-> PC borehole measurement 2003, Eh still decreasing
density 1.01583 ## kg/dm³ = g/cm³
temp 13 ## mean temperature Mont Terri, Bossart & Thury (2008), calculations performed for 25°C
units mol/kgw
## Mean composition
pH 7.064
Na 2.763e-01
Cl 3.228e-01 charge
S(6) 1.653e-02 as SO4
Ca 2.173e-02
Mg 1.740e-02
K 1.902e-03
Sr 4.520e-04
Fe 1.435e-04
U 2.247e-09
SURFACE 2 Opalinus Clay, clay minerals
-equil 2 ## equilibrate with solution 2
-sites_units density ## set unit for binding site density to
## sites/nm2
-donnan 4.9e-10 ## calculated after Wigger & Van Loon (2018)
## for ionic strength after equilibration
## with minerales for pCO2=2.2 log10 bar
## surface density SSA (m2/g) mass (g/kgw)
Kln_aOH 1.155 11. 2798.4 ## Kaolinite 28 wt% (aluminol and silanol sites)
Kln_siOH 1.155
Ill_sOH 0.05 100. 1205.35 ## Illite 31 wt% (weak und strong binding sites)
Ill_wOH 2.26 ## 2 % strong binding sites
Mll_sOH 0.05 100. 113.94 ## Montmorillonite = smektite = 6 wt% (weak und strong binding sites)
Mll_wOH 2.26 ## 2 % strong binding sites
EXCHANGE 2 Exchanger, only illite+montmorillonite
## Illite = 0.225 eq/kg_rock, Montmorillonit = 0.87 eq/kg_rock
-equil 2 ## equilibrate with solution 1
Z 0.5 ## = Illite
Y 0.2 ## = Montmorillonite
END
SOLUTION 3
pe -2.627 ## Eh = -227 mV, Value from Bossart & Thury (2008)-> PC borehole measurement 2003, Eh still decreasing
density 1.01583 ## kg/dm³ = g/cm³
temp 13 ## mean temperature Mont Terri, Bossart & Thury (2008), calculations performed for 25°C
units mol/kgw
## Mean composition
pH 7.064
Na 3.763e-01
Cl 4.228e-01 charge
S(6) 1.653e-02 as SO4
Ca 2.173e-02
Mg 1.740e-02
K 1.902e-03
Sr 4.520e-04
Fe 1.435e-04
U 1e-6
C 1.991e-03
END
RUN_CELLS
END

View File

@ -1,37 +0,0 @@
rows <- 100
cols <- 100
grid_def <- matrix(1, nrow = rows, ncol = cols)
# Define grid configuration for POET model
grid_setup <- list(
pqc_in_file = "./SurfExBase.pqi",
pqc_db_file = "./SMILE_2021_11_01_TH.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(1, 1), # Size of the grid in meters
constant_cells = c() # IDs of cells with constant concentration
)
bound_def <- list(
"type" = rep("constant", cols),
"sol_id" = rep(2, cols),
"cell" = seq(1, cols)
)
diffusion_setup <- list(
boundaries = list(
"N" = bound_def
),
alpha_x = 1e-6,
alpha_y = 1e-6
)
chemistry_setup <- list()
# 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 # Parameters related to the chemistry process
)

View File

@ -1,7 +0,0 @@
iterations <- 10
dt <- 200
list(
timesteps = rep(dt, iterations),
store_result = TRUE
)

View File

@ -1,37 +0,0 @@
rows <- 1000
cols <- 1000
grid_def <- matrix(1, nrow = rows, ncol = cols)
# Define grid configuration for POET model
grid_setup <- list(
pqc_in_file = "./SurfExBase.pqi",
pqc_db_file = "./SMILE_2021_11_01_TH.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(rows, cols) / 10, # Size of the grid in meters
constant_cells = c() # IDs of cells with constant concentration
)
bound_def <- list(
"type" = rep("constant", cols),
"sol_id" = rep(2, cols),
"cell" = seq(1, cols)
)
diffusion_setup <- list(
boundaries = list(
"N" = bound_def
),
alpha_x = 1e-6,
alpha_y = 1e-6
)
chemistry_setup <- list()
# 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 # Parameters related to the chemistry process
)

View File

@ -1,10 +0,0 @@
iterations <- 100
dt <- 200
out_save <- seq(5, iterations, by = 5)
list(
timesteps = rep(dt, iterations),
store_result = TRUE,
out_save = out_save
)

148
poet.yml Normal file
View File

@ -0,0 +1,148 @@
name: poet
channels:
- defaults
- conda-forge
- https://repo.anaconda.com/pkgs/main
- https://repo.anaconda.com/pkgs/r
dependencies:
- _libgcc_mutex=0.1=conda_forge
- _openmp_mutex=4.5=2_gnu
- attr=2.5.2=h39aace5_0
- binutils_impl_linux-64=2.43=h4bf12b8_5
- binutils_linux-64=2.43=h4852527_5
- bzip2=1.0.8=h5eee18b_6
- c-ares=1.34.5=hb9d3cd8_0
- ca-certificates=2025.11.12=hbd8a1cb_0
- cached-property=1.5.2=py_0
- cmake=4.1.2=hc946e07_0
- eigen=3.4.0=h171cf75_1
- expat=2.7.1=h6a678d5_0
- gcc_impl_linux-64=13.3.0=h1e990d8_2
- gcc_linux-64=13.3.0=h6f18a23_11
- gettext=0.21.0=hedfda30_2
- gxx_impl_linux-64=13.3.0=hae580e1_2
- gxx_linux-64=13.3.0=hb14504d_11
- h5py=3.14.0=nompi_py313hfaf8fd4_100
- hdf5=1.14.6=nompi_h6e4c0c1_103
- highfive=2.10.1=he6560a2_2
- icu=73.1=h6a678d5_0
- kernel-headers_linux-64=3.10.0=he073ed8_18
- krb5=1.21.3=h143b758_0
- ld_impl_linux-64=2.43=h712a8e2_5
- libaec=1.1.4=h3f801dc_0
- libblas=3.9.0=32_h59b9bed_openblas
- libcap=2.77=h3ff7636_0
- libcblas=3.9.0=32_he106b2a_openblas
- libcurl=8.16.0=heebcbe5_0
- libedit=3.1.20230828=h5eee18b_0
- libev=4.33=h7f8727e_1
- libevent=2.1.12=hf998b51_1
- libexpat=2.7.1=hecca717_0
- libfabric=2.3.1=ha770c72_1
- libfabric1=2.3.1=h6c8fc0a_1
- libffi=3.4.4=h6a678d5_1
- libgcc=15.1.0=h767d61c_2
- libgcc-devel_linux-64=13.3.0=hc03c837_102
- libgcc-ng=15.1.0=h69a702a_2
- libgfortran=15.1.0=h69a702a_2
- libgfortran-ng=15.1.0=h69a702a_2
- libgfortran5=15.1.0=hcea5267_2
- libgomp=15.1.0=h767d61c_2
- libhwloc=2.12.1=default_h3d81e11_1000
- libiconv=1.16=h5eee18b_3
- libidn2=2.3.8=hf80d704_0
- liblapack=3.9.0=32_h7ac8fdf_openblas
- liblzma=5.8.1=hb9d3cd8_2
- liblzma-devel=5.8.1=hb9d3cd8_2
- libmpdec=4.0.0=h5eee18b_0
- libnghttp2=1.64.0=h161d5f1_0
- libnl=3.11.0=hb9d3cd8_0
- libopenblas=0.3.30=pthreads_h94d23a6_0
- libpmix=5.0.8=h4bd6b51_2
- libsanitizer=13.3.0=he8ea267_2
- libsqlite=3.50.2=h6cd9bfd_0
- libssh2=1.11.1=hcf80075_0
- libstdcxx=15.1.0=h8f9b012_2
- libstdcxx-devel_linux-64=13.3.0=hc03c837_102
- libstdcxx-ng=15.1.0=h4852527_2
- libsystemd0=257.10=hd0affe5_2
- libudev1=257.10=hd0affe5_2
- libunistring=1.3=hb25bd0a_0
- libuuid=2.38.1=h0b41bf4_0
- libuv=1.48.0=h5eee18b_0
- libxcb=1.17.0=h9b100fa_0
- libxml2=2.13.9=h2c43086_0
- libzlib=1.3.1=hb9d3cd8_2
- mpi=1.0.1=openmpi
- ncurses=6.5=h2d0b736_3
- numpy=2.3.0=py313h17eae1a_0
- openmpi=5.0.8=h2fe1745_108
- openssl=3.6.0=h26f9b46_0
- pip=25.1=pyhc872135_2
- pthread-stubs=0.3=h0ce48e5_1
- pybind11=2.13.6=py313hdb19cb5_1
- pybind11-global=2.13.6=py313hdb19cb5_1
- python=3.13.2=hf636f53_101_cp313
- python_abi=3.13=0_cp313
- rdma-core=60.0=hecca717_0
- readline=8.2=h5eee18b_0
- rhash=1.4.6=ha914fed_0
- setuptools=78.1.1=py313h06a4308_0
- sqlite=3.31.1=h7b6447c_0
- sysroot_linux-64=2.17=h0157908_18
- tk=8.6.13=noxft_hd72426e_102
- ucc=1.6.0=hb729f83_0
- ucx=1.19.0=h63b5c0b_5
- wheel=0.45.1=py313h06a4308_0
- xorg-libx11=1.8.12=h9b100fa_1
- xorg-libxau=1.0.12=h9b100fa_0
- xorg-libxdmcp=1.1.5=h9b100fa_0
- xorg-xorgproto=2024.1=h5eee18b_1
- xz=5.8.1=hbcc6ac9_2
- xz-gpl-tools=5.8.1=hbcc6ac9_2
- xz-tools=5.8.1=hb9d3cd8_2
- zlib=1.3.1=hb9d3cd8_2
- zstd=1.5.7=hb8e6e7a_2
- pip:
- absl-py==2.3.1
- ai-benchmarks-utils==0.1.0
- astunparse==1.6.3
- certifi==2025.11.12
- charset-normalizer==3.4.4
- flatbuffers==25.9.23
- gast==0.6.0
- google-pasta==0.2.0
- grpcio==1.76.0
- idna==3.11
- joblib==1.5.2
- keras==3.12.0
- libclang==18.1.1
- markdown==3.10
- markdown-it-py==4.0.0
- markupsafe==3.0.3
- mdurl==0.1.2
- ml-dtypes==0.5.4
- namex==0.1.0
- opt-einsum==3.4.0
- optree==0.18.0
- packaging==25.0
- pandas==2.3.3
- pillow==12.0.0
- protobuf==6.33.1
- pygments==2.19.2
- python-dateutil==2.9.0.post0
- pytz==2025.2
- requests==2.32.5
- rich==14.2.0
- scikit-learn==1.7.2
- scipy==1.16.3
- six==1.17.0
- tensorflow==2.20.0
- termcolor==3.2.0
- threadpoolctl==3.6.0
- typing-extensions==4.15.0
- tzdata==2025.2
- urllib3==2.5.0
- werkzeug==3.1.3
- wrapt==2.0.1
prefix: /mnt/scratch/miniconda3/envs/poet-dummy