Update CMakeLists.txt and add new files for barite_het simulation

CMakeLists.txt

@@ -22,8 +22,8 @@ find_package(MPI REQUIRED)
 
 find_package(RRuntime REQUIRED)
 
-add_compile_options(-fsanitize=address -fno-omit-frame-pointer)
-add_link_options(-fsanitize=address)
+# add_compile_options(-fsanitize=address -fno-omit-frame-pointer)
+# add_link_options(-fsanitize=address)
 
 add_subdirectory(src)
 add_subdirectory(bench)

bench/barite/het/barite_het.pqi

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

bench/barite/het/test.R

@@ -0,0 +1,15 @@
+grid_def <- matrix(c(2, 3), nrow = 2, ncol = 5)
+
+# Define grid configuration for POET model
+grid_setup <- list(
+    pqc_in = "./barite_het.pqi",
+    pqc_db = "./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
+)
+
+# Define a setup list for simulation configuration
+setup <- list(
+    grid = grid_setup # Parameters related to the grid structure
+)

src/Init/GridInit.cpp

@@ -99,10 +99,10 @@ void InitialList::initGrid(const Rcpp::List &grid_input) {
       pqcScriptToGrid(R, script_rp, database_rp, this->pqc_raw_dumps);
   this->initial_grid = matToGrid(R, this->phreeqc_mat, grid_def);
 
-  // R["pqc_mat"] = this->phreeqc_mat;
-  // R["grid_def"] = initial_grid;
+  R["pqc_mat"] = this->phreeqc_mat;
+  R["grid_def"] = initial_grid;
 
-  // R.parseEval("print(pqc_mat)");
-  // R.parseEval("print(grid_def)");
+  R.parseEval("print(pqc_mat)");
+  R.parseEval("print(grid_def)");
 }
 } // namespace poet