Merge branch '17-new-barite-benchmarks-for-fgcs-journal' into 'main'

Draft: Resolve "New barite benchmarks for FGCS journal"

Closes #17

See merge request naaice/poet!44

bench/dolo/dolo_interp.R

@@ -41,36 +41,36 @@ diffusion_setup <- list(
 )
 
 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)
-    }
+    # 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)])
-}
+# fuzz_input_dht_keys <- function(input) {
+#     dht_species <- c(
+#         "H" = 3,
+#         "O" = 3,
+#         "Charge" = 3,
+#         "C" = 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,
+        "C" = 6,
         "Ca" = 6,
         "Cl" = 3,
         "Mg" = 5,
@@ -95,7 +95,7 @@ check_neg_cal_dol <- function(result) {
 # significant digits)
 
 pht_species <- c(
-    "C(4)" = 3,
+    "C" = 3,
     "Ca" = 3,
     "Mg" = 2,
     "Calcite" = 2,
@@ -107,7 +107,7 @@ chemistry_setup <- list(
         "H" = 3,
         "O" = 3,
         "Charge" = 3,
-        "C(4)" = 6,
+        "C" = 6,
         "Ca" = 6,
         "Cl" = 3,
         "Mg" = 5,
@@ -117,7 +117,7 @@ chemistry_setup <- list(
     pht_species = pht_species,
     hooks = list(
         dht_fill = check_sign_cal_dol_dht,
-        dht_fuzz = fuzz_input_dht_keys,
+        # dht_fuzz = fuzz_input_dht_keys,
         interp_pre = check_sign_cal_dol_interp,
         interp_post = check_neg_cal_dol
     )

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 <delucia@gfz.de>}
+\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:

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

bench/fgcs/README.org

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

bench/fgcs/barite_fgcs_2.R

@@ -0,0 +1,105 @@
+## 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
+)

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

bench/fgcs/barite_fgcs_2_rt.R

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

src/CMakeLists.txt

@@ -1,7 +1,4 @@
-add_library(POETLib)
-
-target_sources(POETLib 
-  PRIVATE
+add_library(POETLib
   Init/InitialList.cpp 
   Init/GridInit.cpp 
   Init/DiffusionInit.cpp

src/Chemistry/SurrogateModels/DHT.c

@@ -105,6 +105,8 @@ DHT *DHT_create(MPI_Comm comm, uint64_t size, unsigned int data_size,
   object->index_count = 9 - (index_bytes / 8);
   object->index = (uint64_t *)malloc((object->index_count) * sizeof(uint64_t));
   object->mem_alloc = mem_alloc;
+  object->sum_idx = 0;
+  object->cnt_idx = 0;
 
   // if set, initialize dht_stats
 #ifdef DHT_STATISTICS
@@ -189,6 +191,9 @@ int DHT_write_accumulate(DHT *table, const void *send_key, int data_size,
     }
   }
 
+  table->cnt_idx += 1;
+  table->sum_idx += (i + 1);
+
   if (result == DHT_WRITE_SUCCESS_WITH_EVICTION) {
     memset((char *)table->send_entry + 1 + table->key_size, '\0',
            table->data_size);
@@ -279,6 +284,9 @@ int DHT_write(DHT *table, void *send_key, void *send_data, uint32_t *proc,
     }
   }
 
+  table->cnt_idx += 1;
+  table->sum_idx += (i + 1);
+
   // put data to DHT (with last selected index by value i)
   if (MPI_Put(table->send_entry, 1 + table->data_size + table->key_size,
               MPI_BYTE, dest_rank, table->index[i],
@@ -550,6 +558,41 @@ int DHT_free(DHT *table, int *eviction_counter, int *readerror_counter) {
   return DHT_SUCCESS;
 }
 
+float DHT_get_used_idx_factor(DHT *table, int with_reset) {
+  int rank;
+  MPI_Comm_rank(table->communicator, &rank);
+
+  float my_avg_idx = (float)table->sum_idx / (float)table->cnt_idx;
+
+  float max_mean_index;
+
+  MPI_Reduce(&my_avg_idx, &max_mean_index, 1, MPI_FLOAT, MPI_MAX, 0,
+             table->communicator);
+
+  MPI_Bcast(&max_mean_index, 1, MPI_FLOAT, 0, table->communicator);
+
+  if (!!with_reset) {
+    table->sum_idx = 0;
+    table->cnt_idx = 0;
+  }
+
+  return max_mean_index;
+}
+
+int DHT_flush(DHT *table) {
+  // make sure all processes are synchronized
+  MPI_Barrier(table->communicator);
+
+  // wipe local memory with zeros
+  memset(table->mem_alloc, '\0',
+         table->table_size * (1 + table->data_size + table->key_size));
+
+  table->sum_idx = 0;
+  table->cnt_idx = 0;
+
+  return DHT_SUCCESS;
+}
+
 int DHT_print_statistics(DHT *table) {
 #ifdef DHT_STATISTICS
   int *written_buckets;

src/Chemistry/SurrogateModels/DHT.h

@@ -117,6 +117,9 @@ typedef struct {
   unsigned int index_count;
 
   int (*accumulate_callback)(int, void *, int, void *);
+
+  size_t sum_idx;
+  size_t cnt_idx;
 #ifdef DHT_STATISTICS
   /** Detailed statistics of the usage of the DHT. */
   DHT_stats *stats;
@@ -125,10 +128,11 @@ typedef struct {
 
 extern void DHT_set_accumulate_callback(DHT *table,
                                         int (*callback_func)(int, void *, int,
-                                                              void *));
+                                                             void *));
 
 extern int DHT_write_accumulate(DHT *table, const void *key, int send_size,
-                                void *data, uint32_t *proc, uint32_t *index, int *callback_ret);
+                                void *data, uint32_t *proc, uint32_t *index,
+                                int *callback_ret);
 
 /**
  * @brief Create a DHT.
@@ -284,44 +288,8 @@ extern int DHT_free(DHT *table, int *eviction_counter, int *readerror_counter);
  */
 extern int DHT_print_statistics(DHT *table);
 
-/**
- * @brief Determine destination rank and index.
- *
- * This is done by looping over all possbile indices. First of all, set a
- * temporary index to zero and copy count of bytes for each index into the
- * memory area of the temporary index. After that the current index is
- * calculated by the temporary index modulo the table size. The destination rank
- * of the process is simply determined by hash modulo the communicator size.
- *
- * @param hash Calculated 64 bit hash.
- * @param comm_size Communicator size.
- * @param table_size Count of buckets per process.
- * @param dest_rank Reference to the destination rank variable.
- * @param index Pointer to the array index.
- * @param index_count Count of possible indeces.
- */
-static void determine_dest(uint64_t hash, int comm_size,
-                           unsigned int table_size, unsigned int *dest_rank,
-                           uint64_t *index, unsigned int index_count);
+extern float DHT_get_used_idx_factor(DHT *table, int with_reset);
 
-/**
- * @brief Set the occupied flag.
- *
- * This will set the first bit of a bucket to 1.
- *
- * @param flag_byte First byte of a bucket.
- */
-static void set_flag(char *flag_byte);
-
-/**
- * @brief Get the occupied flag.
- *
- * This function determines whether the occupied flag of a bucket was set or
- * not.
- *
- * @param flag_byte First byte of a bucket.
- * @return int Returns 1 for true or 0 for false.
- */
-static int read_flag(char flag_byte);
+extern int DHT_flush(DHT *table);
 
 #endif /* DHT_H */

src/Chemistry/SurrogateModels/DHT_Wrapper.cpp

@@ -133,6 +133,10 @@ void DHT_Wrapper::fillDHT(const WorkPackage &work_package) {
       continue;
     }
 
+    if (work_package.input[i][0] != 2) {
+      continue;
+    }
+
     // check if calcite or dolomite is absent and present, resp.n and vice
     // versa in input/output. If this is the case -> Do not write to DHT!
     // HACK: hardcoded, should be fixed!
@@ -140,7 +144,7 @@ void DHT_Wrapper::fillDHT(const WorkPackage &work_package) {
       NamedVector<double> old_values(output_names, work_package.input[i]);
       NamedVector<double> new_values(output_names, work_package.output[i]);
 
-      if (hooks.dht_fill(old_values, new_values)) {
+      if (Rcpp::as<bool>(hooks.dht_fill(old_values, new_values))) {
         continue;
       }
     }

src/Chemistry/SurrogateModels/Interpolation.hpp

@@ -166,6 +166,8 @@ public:
 
   enum result_status { RES_OK, INSUFFICIENT_DATA, NOT_NEEDED };
 
+  DHT *getDHTObject() { return this->pht->getDHTObject(); }
+
   struct InterpolationResult {
     std::vector<std::vector<double>> results;
     std::vector<result_status> status;

src/Chemistry/SurrogateModels/InterpolationModule.cpp

@@ -76,6 +76,11 @@ void InterpolationModule::tryInterpolation(WorkPackage &work_package) {
   const auto dht_results = this->dht_instance.getDHTResults();
 
   for (int wp_i = 0; wp_i < work_package.size; wp_i++) {
+    if (work_package.input[wp_i][0] != 2) {
+      interp_result.status[wp_i] = INSUFFICIENT_DATA;
+      continue;
+    }
+
     if (work_package.mapping[wp_i] != CHEM_PQC) {
       interp_result.status[wp_i] = NOT_NEEDED;
       continue;
@@ -140,7 +145,7 @@ void InterpolationModule::tryInterpolation(WorkPackage &work_package) {
 
     if (hooks.interp_post.isValid()) {
       NamedVector<double> nv_result(this->out_names, work_package.output[wp_i]);
-      if (hooks.interp_post(nv_result)) {
+      if (Rcpp::as<bool>(hooks.interp_post(nv_result))) {
         interp_result.status[wp_i] = INSUFFICIENT_DATA;
         continue;
       }

src/Chemistry/WorkerFunctions.cpp

@@ -9,7 +9,6 @@
 #include <cstdint>
 #include <iomanip>
 #include <iostream>
-#include <map>
 #include <mpi.h>
 #include <stdexcept>
 #include <string>
@@ -247,6 +246,17 @@ void poet::ChemistryModule::WorkerPostIter(MPI_Status &prope_status,
           << std::setw(this->file_pad) << iteration << ".pht";
       interp->dumpPHTState(out.str());
     }
+
+    const auto max_mean_idx =
+        DHT_get_used_idx_factor(this->interp->getDHTObject(), 1);
+
+    if (max_mean_idx >= 2) {
+      DHT_flush(this->interp->getDHTObject());
+      DHT_flush(this->dht->getDHT());
+      if (this->comm_rank == 2) {
+        std::cout << "Flushed both DHT and PHT!\n\n";
+      }
+    }
   }
 
   RInsidePOET::getInstance().parseEvalQ("gc()");