diff --git a/bench/fgcs/20241211_README.pdf b/bench/fgcs/20241211_README.pdf new file mode 100644 index 000000000..2cd4332f2 Binary files /dev/null and b/bench/fgcs/20241211_README.pdf differ diff --git a/bench/fgcs/20241211_README.tex b/bench/fgcs/20241211_README.tex new file mode 100644 index 000000000..579445633 --- /dev/null +++ b/bench/fgcs/20241211_README.tex @@ -0,0 +1,102 @@ +% Created 2024-12-11 mer 23:24 +% Intended LaTeX compiler: pdflatex +\documentclass[a4paper, 9pt]{article} +\usepackage[utf8]{inputenc} +\usepackage[T1]{fontenc} +\usepackage{graphicx} +\usepackage{longtable} +\usepackage{wrapfig} +\usepackage{rotating} +\usepackage[normalem]{ulem} +\usepackage{amsmath} +\usepackage{amssymb} +\usepackage{capt-of} +\usepackage{hyperref} +\usepackage{fullpage} +\usepackage{amsmath} +\usepackage{graphicx} +\usepackage{charter} +\usepackage{listings} +\lstloadlanguages{R} +\author{MDL } +\date{2024-12-11} +\title{A \texttt{barite}-based benchmark for FGCS interpolation paper} + + + +\begin{document} + +\maketitle + +\section{Description} +\label{sec:org739879a} + +\begin{itemize} +\item \texttt{barite\_fgcs\_2.R}: POET input script with circular + "crystals" on a 200x200 nodes grid + +\item \(\alpha\): isotropic 10\textsuperscript{-5} + m\textsuperscript{2}/s outside of the crystals, + 10\textsuperscript{-7} inside +\item 200 iterations, dt = 1000 +\item \texttt{barite\_fgcs\_2.pqi}: PHREEQC input, 4 SOLUTIONS + (basically the same as in \texttt{barite} benchmark): + \begin{enumerate} + \item Equilibrium with Celestite, no mineral \(Rightarrow\) + \item Equilibrium with Celestite, KINETICS Celestite (1 mol) and + Barite (0 mol) + \item Injection of 0.1 BaCl2 from NW corner + \item Injection of 0.2 BaCl2 from SE corner +\end{enumerate} + +\item \texttt{db\_barite.dat}: PHREEQC database containing the kinetic + expressions for barite and celestite, stripped down from + \texttt{phreeqc.dat} +\end{itemize} + +\begin{figure}[htbp] + \centering + \includegraphics[width=0.48\textwidth]{./fgcs_Celestite_init.pdf} + \includegraphics[width=0.48\textwidth]{./fgcs_Barite_200.pdf} + \caption{\textbf{Left:} Initial distribution of Celestite + "crystals". \textbf{Right:} precipitated Barite} +\end{figure} + +\section{Interpolation} +\label{sec:org2a09431} + +Using the following parametrization: + +\begin{lstlisting} +dht_species <- c("H" = 7, + "O" = 7, + "Ba" = 7, + "Cl" = 7, + "S(6)" = 7, + "Sr" = 7, + "Barite" = 4, + "Celestite" = 4) + +pht_species <- c("Ba" = 4, + "Cl" = 3, + "S(6)" = 3, + "Sr" = 3, + "Barite" = 2, + "Celestite" = 2 ) +\end{lstlisting} + +Runtime goes from 1800 to 600 s (21 CPUs) but there are "suspect" +errors especially in O and H, where "suspect" means some values appear +to be multiplied by 2: +\begin{figure}[htbp] + \centering + \includegraphics[width=0.9\textwidth]{./fgcs_interp_1.pdf} + \caption{Scatterplots reference vs interpolated after 1 coupling + iteration} +\end{figure} +\end{document} + +%%% Local Variables: +%%% mode: LaTeX +%%% TeX-master: t +%%% End: diff --git a/bench/fgcs/EvalFGCS.R b/bench/fgcs/EvalFGCS.R new file mode 100644 index 000000000..5176c3e99 --- /dev/null +++ b/bench/fgcs/EvalFGCS.R @@ -0,0 +1,90 @@ +## 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 diff --git a/bench/fgcs/README.org b/bench/fgcs/README.org new file mode 100644 index 000000000..0fb5acb53 --- /dev/null +++ b/bench/fgcs/README.org @@ -0,0 +1,2 @@ +* Refer to the LaTeX file (and pdf) for more information + diff --git a/bench/fgcs/barite_fgcs_2.R b/bench/fgcs/barite_fgcs_2.R new file mode 100644 index 000000000..500b372fd --- /dev/null +++ b/bench/fgcs/barite_fgcs_2.R @@ -0,0 +1,105 @@ +## Time-stamp: "Last modified 2024-12-11 16:08:11 delucia" + +cols <- 200 +rows <- 200 + +dim_cols <- 10 +dim_rows <- 10 + +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 = "./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(6)" = 7, + "Sr" = 7, + "Barite" = 4, + "Celestite" = 4 +) + +pht_species <- c( + "Ba" = 4, + "Cl" = 3, + "S(6)" = 3, + "Sr" = 3, + "Barite" = 2, + "Celestite" = 2 +) + +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 +) diff --git a/bench/fgcs/barite_fgcs_2.pqi b/bench/fgcs/barite_fgcs_2.pqi new file mode 100644 index 000000000..f806d5aac --- /dev/null +++ b/bench/fgcs/barite_fgcs_2.pqi @@ -0,0 +1,49 @@ +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 diff --git a/bench/fgcs/barite_fgcs_2_rt.R b/bench/fgcs/barite_fgcs_2_rt.R new file mode 100644 index 000000000..373ad79dd --- /dev/null +++ b/bench/fgcs/barite_fgcs_2_rt.R @@ -0,0 +1,7 @@ +iterations <- 200 +dt <- 1000 + +list( + timesteps = rep(dt, iterations), + store_result = TRUE +)