Stop tracking bin/dolo/ and update .gitignore

updated .gitignore
Store per-cell MAPE/RRMSE and add stabilization interval
2025-12-15 12:28:22 +01:00 · 2025-11-13 16:00:47 +01:00 · 2025-11-12 10:25:58 +01:00 · 2025-11-12 10:20:28 +01:00 · 2025-11-07 15:55:26 +01:00 · 2025-11-07 15:52:39 +01:00
47 changed files with 1482 additions and 3478 deletions
--- a/.gitignore
+++ b/.gitignore
@ -143,4 +143,23 @@ build/
 /.cache/

 .vscode
-.codechecker
+.codechecker
+
+# Prevent upload of local installations
+bin/*
+share/
+lib/
+include/
+
+# But keep these specific files
+!bin/plot/
+!bin/barite_fgcs_3.pqi
+!bin/barite_fgcs_4_rt.R
+!bin/barite_fgcs_4.R
+!bin/barite_fgcs_4.qs2
+!bin/db_barite.dat
+!bin/dolo_fgcs_3.qs2
+!bin/dolo_fgcs_3.R
+!bin/dolo_fgcs_3.pqi
+!bin/phreeqc_kin.dat
+!bin/run_poet.shbin/dolo/
--- a/bin/barite_fgcs_2.R
+++ b/bin/barite_fgcs_2.R
@ -1,121 +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
-)
-
-
-iterations <- 200
-dt <- 1000
-control_iteration <- 50
-species_epsilon <- rep(0.001, length(dht_species))
-out_save <- seq(50, iterations, by = 50)
-
-
-list(
-    timesteps = rep(dt, iterations),
-    store_result = TRUE,
-    out_save = out_save,
-    control_iteration = control_iteration,
-    species_epsilon = species_epsilon
-)
--- a/bin/barite_fgcs_2.pqi
+++ b/bin/barite_fgcs_2.pqi
@ -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
--- a/bin/barite_het.qs2
+++ b/bin/barite_het.qs2
--- a/bin/barite_het_rt.R
+++ b/bin/barite_het_rt.R
@ -1,6 +0,0 @@
-list(
-    timesteps = rep(50, 100),
-    store_result = TRUE,
-    control_iteration <- 20,
-    species_epsilon <- c(0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001)
-)
--- a/bin/db_barite.dat
+++ b/bin/db_barite.dat
@ -1,195 +0,0 @@
-DATABASE
-	###########################
-SOLUTION_MASTER_SPECIES
-	H      H+     -1            H         1.008  # phreeqc/
-	H(0)   H2      0            H                # phreeqc/
-	H(1)   H+     -1     0.0                     # phreeqc/
-	E      e-      0     0.0                0.0  # phreeqc/
-	O      H2O     0            O          16.0  # phreeqc/
-	O(0)   O2      0            O                # phreeqc/
-	O(-2)  H2O     0     0.0                     # phreeqc/
-	Na     Na+     0            Na      22.9898  # phreeqc/
-	Ba     Ba+2    0            Ba       137.34  # phreeqc/
-	Sr     Sr+2    0            Sr        87.62  # phreeqc/
-	Cl     Cl-     0            Cl       35.453  # phreeqc/
-	S      SO4-2   0            SO4      32.064  # phreeqc/
-	S(6)   SO4-2   0            SO4              # phreeqc/
-	S(-2)  HS-     1            S                # phreeqc/
-SOLUTION_SPECIES
-	H+ = H+
-		-gamma 9 0
-		-dw 9.31e-09
-		# source: phreeqc
-	e- = e-
-		# source: phreeqc
-	H2O = H2O
-		# source: phreeqc
-	Na+ = Na+
-		-gamma 4.08 0.082
-		-dw 1.33e-09
-		-Vm 2.28 -4.38 -4.1 -0.586 0.09 4 0.3 52 -0.00333 0.566
-		# source: phreeqc
-	Ba+2 = Ba+2
-		-gamma 4 0.153
-		-dw 8.48e-10
-		-Vm 2.063 -10.06 1.9534 -2.36 0.4218 5 1.58 -12.03 -0.00835 1
-		# source: phreeqc
-	Sr+2 = Sr+2
-		-gamma 5.26 0.121
-		-dw 7.94e-10
-		-Vm -0.0157 -10.15 10.18 -2.36 0.86 5.26 0.859 -27 -0.0041 1.97
-		# source: phreeqc
-	Cl- = Cl-
-		-gamma 3.63 0.017
-		-dw 2.03e-09
-		-Vm 4.465 4.801 4.325 -2.847 1.748 0 -0.331 20.16 0 1
-		# source: phreeqc
-	SO4-2 = SO4-2
-		-gamma 5 -0.04
-		-dw 1.07e-09
-		-Vm 8 2.3 -46.04 6.245 3.82 0 0 0 0 1
-		# source: phreeqc
-	H2O = OH- + H+
-		-analytical_expression 293.29227 0.1360833 -10576.913 -123.73158 0 -6.996455e-05
-		-gamma 3.5 0
-		-dw 5.27e-09
-		-Vm -9.66 28.5 80 -22.9 1.89 0 1.09 0 0 1
-		# source: phreeqc
-	2 H2O = O2 + 4 H+ + 4 e-
-		-log_k -86.08
-		-delta_h 134.79 kcal
-		-dw 2.35e-09
-		-Vm 5.7889 6.3536 3.2528 -3.0417 -0.3943     
-		# source: phreeqc
-	2 H+ + 2 e- = H2
-		-log_k -3.15
-		-delta_h -1.759 kcal
-		-dw 5.13e-09
-		-Vm 6.52 0.78 0.12       
-		# source: phreeqc
-	SO4-2 + H+ = HSO4-
-		-log_k 1.988
-		-delta_h 3.85 kcal
-		-analytical_expression -56.889 0.006473 2307.9 19.8858  
-		-dw 1.33e-09
-		-Vm 8.2 9.259 2.1108 -3.1618 1.1748 0 -0.3 15 0 1
-		# source: phreeqc
-	HS- = S-2 + H+
-		-log_k -12.918
-		-delta_h 12.1 kcal
-		-gamma 5 0
-		-dw 7.31e-10
-		# source: phreeqc
-	SO4-2 + 9 H+ + 8 e- = HS- + 4 H2O
-		-log_k 33.65
-		-delta_h -60.140 kcal
-		-gamma 3.5 0
-		-dw 1.73e-09
-		-Vm 5.0119 4.9799 3.4765 -2.9849 1.441     
-		# source: phreeqc
-	HS- + H+ = H2S
-		-log_k 6.994
-		-delta_h -5.30 kcal
-		-analytical_expression -11.17 0.02386 3279   
-		-dw 2.1e-09
-		-Vm 7.81 2.96 -0.46       
-		# source: phreeqc
-	Na+ + OH- = NaOH
-		-log_k -10
-		# source: phreeqc
-	Na+ + SO4-2 = NaSO4-
-		-log_k 0.7
-		-delta_h 1.120 kcal
-		-gamma 5.4 0
-		-dw 1.33e-09
-		-Vm 1e-05 16.4 -0.0678 -1.05 4.14 0 6.86 0 0.0242 0.53
-		# source: phreeqc
-	Ba+2 + H2O = BaOH+ + H+
-		-log_k -13.47
-		-gamma 5 0
-		# source: phreeqc
-	Ba+2 + SO4-2 = BaSO4
-		-log_k 2.7
-		# source: phreeqc
-	Sr+2 + H2O = SrOH+ + H+
-		-log_k -13.29
-		-gamma 5 0
-		# source: phreeqc
-	Sr+2 + SO4-2 = SrSO4
-		-log_k 2.29
-		-delta_h 2.08 kcal
-		-Vm 6.791 -0.9666 6.13 -2.739 -0.001     
-		# source: phreeqc
-PHASES
-	Barite
-		BaSO4 = Ba+2 + SO4-2
-		-log_k -9.97
-		-delta_h 6.35 kcal
-		-analytical_expression -282.43 -0.08972 5822 113.08  
-		-Vm 52.9
-		# source: phreeqc
-		# comment: 
-	Celestite
-		SrSO4 = Sr+2 + SO4-2
-		-log_k -6.63
-		-delta_h -4.037 kcal
-		-analytical_expression -7.14 0.00611 75 0 0 -1.79e-05
-		-Vm 46.4
-		# source: phreeqc
-		# comment: 
-RATES
-	Celestite  # Palandri & Kharaka 2004<--------------------------------change me
-# PARM(1): reactive surface area
-# am: acid mechanism, nm: neutral mechanism, bm: base mechanism
-    -start
-        10 sr_i = SR("Celestite")  # saturation ratio, (-)<----------change me
-        20 moles = 0  # init target variable, (mol)
-        30 IF ((M <= 0) AND (sr_i < 1)) OR (sr_i = 1.0) THEN GOTO 310
-        40 sa = PARM(1)  # reactive surface area, (m2)
-
-       100 r = 8.314462  # gas constant, (J K-1 mol-1)
-       110 dTi = (1 / TK) - (1 / 298.15)  # (K-1)
-       120 ea_am = 23800  # activation energy am, (J mol-1)<-----------change me
-       130 ea_nm = 0  # activation energy nm, (J mol-1)<-----------change me
-       140 ea_bm = 0  # activation energy bm, (J mol-1)<-----------change me
-       150 log_k_am = -5.66  # reaction constant am<-------------------change me
-       rem log_k_nm = -99  # reaction constant nm<-------------------change me
-       rem log_k_bm = -99  # reaction constant bm<-------------------change me
-       180 n_am = 0.109  # H+ reaction order am<-----------------------change me
-       rem n_bm = 0  # H+ reaction order bm<-----------------------change me
-       200 am = (10 ^ log_k_am) * EXP(-ea_am * dTi / r) * ACT("H+") ^ n_am
-       rem nm = (10 ^ log_k_nm) * EXP(-ea_nm * dTi / r)
-       rem bm = (10 ^ log_k_bm) * EXP(-ea_bm * dTi / r) * ACT("H+") ^ n_bm
-
-       300 moles = sa * (am) * (1 - sr_i)
-       310 save moles * time
-    -end
-
-	Barite  # Palandri & Kharaka 2004<-----------------------------------change me
-# PARM(1): reactive surface area
-# am: acid mechanism, nm: neutral mechanism, bm: base mechanism
-    -start
-        10 sr_i = SR("Barite")  # saturation ratio, (-)<----------change me
-        20 moles = 0  # init target variable, (mol)
-        30 IF ((M <= 0) AND (sr_i < 1)) OR (sr_i = 1.0) THEN GOTO 310
-        40 sa = PARM(1)  # reactive surface area, (m2)
-
-       100 r = 8.314462  # gas constant, (J K-1 mol-1)
-       110 dTi = (1 / TK) - (1 / 298.15)  # (K-1)
-       120 ea_am = 30800  # activation energy am, (J mol-1)<---------change me
-       130 ea_nm = 30800  # activation energy nm, (J mol-1)<---------change me
-       rem ea_bm = 0  # activation energy bm, (J mol-1)<---------change me
-       150 log_k_am = -6.90  # reaction constant am<-----------------change me
-       160 log_k_nm = -7.90  # reaction constant nm<-----------------change me
-       rem log_k_bm = -99  # reaction constant bm<-------------------change me
-       180 n_am = 0.22  # H+ reaction order am<----------------------change me
-       rem n_bm = 0  # H+ reaction order bm<-----------------------change me
-       200 am = (10 ^ log_k_am) * EXP(-ea_am * dTi / r) * ACT("H+") ^ n_am
-       210 nm = (10 ^ log_k_nm) * EXP(-ea_nm * dTi / r)
-       rem bm = (10 ^ log_k_bm) * EXP(-ea_bm * dTi / r) * ACT("H+") ^ n_bm
-
-       300 moles = sa * (am + nm) * (1 - sr_i)
-       310 save moles * time
-    -end
-
-END
--- a/bin/dolo/iter_00.qs2
+++ b/bin/dolo/iter_00.qs2
--- a/bin/dolo/iter_000.qs2
+++ b/bin/dolo/iter_000.qs2
--- a/bin/dolo/iter_050.qs2
+++ b/bin/dolo/iter_050.qs2
--- a/bin/dolo/iter_100.qs2
+++ b/bin/dolo/iter_100.qs2
--- a/bin/dolo/iter_150.qs2
+++ b/bin/dolo/iter_150.qs2
--- a/bin/dolo/iter_200.qs2
+++ b/bin/dolo/iter_200.qs2
--- a/bin/dolo/iter_250.qs2
+++ b/bin/dolo/iter_250.qs2
--- a/bin/dolo/iter_300.qs2
+++ b/bin/dolo/iter_300.qs2
--- a/bin/dolo/iter_350.qs2
+++ b/bin/dolo/iter_350.qs2
--- a/bin/dolo/iter_400.qs2
+++ b/bin/dolo/iter_400.qs2
--- a/bin/dolo/iter_450.qs2
+++ b/bin/dolo/iter_450.qs2
--- a/bin/dolo/iter_500.qs2
+++ b/bin/dolo/iter_500.qs2
--- a/bin/dolo/timings.qs2
+++ b/bin/dolo/timings.qs2
--- a/bin/dolo_fgcs.pqi
+++ b/bin/dolo_fgcs.pqi
@ -1,48 +0,0 @@
-SOLUTION 1
-    units	mol/kgw
-    water	1
-    temperature	25
-    pH	7
-    pe	4
-PURE 1
-    Calcite 0.0 1
-END
-
-RUN_CELLS
-	-cells 1
-END
-
-COPY solution 1 2
-
-#PURE 2
-#    O2g -0.1675 10
-KINETICS 2
-    Calcite
-    -m      0.00207
-    -parms  0.05
-    -tol    1e-10
-    Dolomite
-    -m      0.0
-    -parms  0.01
-    -tol    1e-10
-END
-
-SOLUTION 3
-    units	mol/kgw
-    water 1
-    temp 25 
-    Mg 0.001
-    Cl 0.002
-END
-
-SOLUTION 4
-    units mol/kgw 
-    water 1
-    temp 25
-    Mg 0.002
-    Cl 0.004
-END
-
-RUN_CELLS 
-    -cells 2-4
-END
--- a/bin/dolo_fgcs_3.R
+++ b/bin/dolo_fgcs_3.R
@ -1,131 +0,0 @@
-rows <- 400
-cols <- 400
-
-grid_def <- matrix(2, nrow = rows, ncol = cols)
-
-# Define grid configuration for POET model
-grid_setup <- list(
-    pqc_in_file = "./dolo_fgcs.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(5, 5), # Size of the grid in meters
-    constant_cells = c() # IDs of cells with constant concentration
-)
-
-bound_def_we <- list(
-    "type" = rep("constant", rows),
-    "sol_id" = rep(1, rows),
-    "cell" = seq(1, rows)
-)
-
-bound_def_ns <- list(
-    "type" = rep("constant", cols),
-    "sol_id" = rep(1, cols),
-    "cell" = seq(1, cols)
-)
-
-diffusion_setup <- list(
-    boundaries = list(
-        "W" = bound_def_we,
-        "E" = bound_def_we,
-        "N" = bound_def_ns,
-        "S" = bound_def_ns
-    ),
-    inner_boundaries = list(
-        "row" = floor(rows / 2),
-        "col" = floor(cols / 2),
-        "sol_id" = c(3)
-    ),
-    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)
-}
-
-check_sign_cal_dol_interp <- function(to_interp, data_set) {
-    dht_species <- c(
-        "H" = 3,
-        "O" = 3,
-        "C" = 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"  = 3,
-    "Ca" = 3,
-    "Mg" = 3,
-    "Cl" = 3,
-    "Calcite" = 3,
-    "Dolomite" = 3
-)
-
-
-dht_species <- c(
-        "H" = 3,
-        "O" = 3,
-        "C" = 6,
-        "Ca" = 6,
-        "Cl" = 3,
-        "Mg" = 5,
-        "Calcite" = 4,
-        "Dolomite" = 4)
-
-chemistry_setup <- list(
-    dht_species = dht_species,
-    pht_species = pht_species,
-    hooks = list(
-        dht_fill = check_sign_cal_dol_dht,
-        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
-)
-
-iterations <- 20
-dt <- 100
-control_iteration <- 10
-species_epsilon <- rep(0.001, length(dht_species))
-#out_save <- seq(50, iterations, by = 50)
-
-
-list(
-    timesteps = rep(dt, iterations),
-    store_result = TRUE,
-    #out_save = out_save,
-    control_iteration = control_iteration,
-    species_epsilon = species_epsilon
-)
--- a/bin/dolo_fgcs_3.qs2
+++ b/bin/dolo_fgcs_3.qs2
--- a/bin/dolo_interp.qs2
+++ b/bin/dolo_interp.qs2
--- a/bin/dolo_interp_rt.R
+++ b/bin/dolo_interp_rt.R
@ -1,17 +0,0 @@
-iterations <- 500
-dt <- 100
-control_iteration <- 20
-species_epsilon <- c(.0, 1e-12, 1e-1, 1e-1, 1e-2, 1e-2, 1e+1, 1e+1, 1e-3, 0, 1e-1, .0, 0.12, .0)
-out_save <- seq(50, iterations, by = 50)
-penalty_iteration <- 40
-max_penalty_iteration <- 80
-
-list(
-    timesteps = rep(dt, iterations),
-    store_result = TRUE,
-    out_save = out_save,
-    control_iteration = control_iteration,
-    species_epsilon = species_epsilon,
-    penalty_iteration = penalty_iteration,
-    max_penalty_iteration = max_penalty_iteration
-)
--- a/bin/phreeqc_kin.dat
+++ b/bin/phreeqc_kin.dat
--- a/bin/poet
+++ b/bin/poet
--- a/bin/poet_init
+++ b/bin/poet_init
--- a/bin/stats_overview
+++ b/bin/stats_overview
@ -1,376 +0,0 @@
-Iteration, Species, MAPE, RRSME 
-20,	ID,	0,	0
-20,	H,	1.41445e-14,	2.35911e-16
-20,	O,	0.0553419,	0.000567653
-20,	Charge,	0.0182496,	0.000187253
-20,	C,	9.41889e-05,	1.47295e-05
-20,	Ca,	9.42715e-05,	1.47294e-05
-20,	Cl,	6.0345,	0.0869993
-20,	Mg,	6.09838,	0.0882693
-20,	Calcite,	5.63369e-05,	8.82505e-06
-20,	Calcite_p1,	0,	0
-20,	Dolomite,	4.13256e-09,	1.16886e-08
-20,	Dolomite_p1,	0,	0
-20,	O2g_eq,	0.145118,	0.0014895
-20,	O2g_si,	0,	0
-
-40,	ID,	0,	0
-40,	H,	3.88971e-14,	4.87626e-16
-40,	O,	0.056556,	0.000586363
-40,	Charge,	0.0175789,	0.000183408
-40,	C,	0.000325979,	3.81983e-05
-40,	Ca,	0.000326053,	3.8197e-05
-40,	Cl,	3.04641,	0.058268
-40,	Mg,	3.09965,	0.0603568
-40,	Calcite,	0.000196233,	2.30896e-05
-40,	Calcite_p1,	0,	0
-40,	Dolomite,	8.48003e-06,	8.69316e-06
-40,	Dolomite_p1,	0,	0
-40,	O2g_eq,	0.148463,	0.00154008
-40,	O2g_si,	0,	0
-
-60,	ID,	0,	0
-60,	H,	0,	0
-60,	O,	0,	0
-60,	Charge,	0,	0
-60,	C,	0,	0
-60,	Ca,	0,	0
-60,	Cl,	0,	0
-60,	Mg,	0,	0
-60,	Calcite,	0,	0
-60,	Calcite_p1,	0,	0
-60,	Dolomite,	0,	0
-60,	Dolomite_p1,	0,	0
-60,	O2g_eq,	0,	0
-60,	O2g_si,	0,	0
-
-80,	ID,	0,	0
-80,	H,	0,	0
-80,	O,	0,	0
-80,	Charge,	0,	0
-80,	C,	0,	0
-80,	Ca,	0,	0
-80,	Cl,	0,	0
-80,	Mg,	0,	0
-80,	Calcite,	0,	0
-80,	Calcite_p1,	0,	0
-80,	Dolomite,	0,	0
-80,	Dolomite_p1,	0,	0
-80,	O2g_eq,	0,	0
-80,	O2g_si,	0,	0
-
-100,	ID,	0,	0
-100,	H,	2.52054e-14,	4.54572e-16
-100,	O,	0.0499718,	0.000533381
-100,	Charge,	0.0157197,	0.000172845
-100,	C,	0.000521522,	3.64097e-05
-100,	Ca,	0.000521537,	3.64034e-05
-100,	Cl,	1.79663,	0.0479452
-100,	Mg,	1.89434,	0.050751
-100,	Calcite,	0.000314689,	2.20679e-05
-100,	Calcite_p1,	0,	0
-100,	Dolomite,	0.000611171,	0.000558382
-100,	Dolomite_p1,	0,	0
-100,	O2g_eq,	0.130299,	0.0013932
-100,	O2g_si,	0,	0
-
-120,	ID,	0,	0
-120,	H,	2.45742e-14,	4.82097e-16
-120,	O,	0.0511434,	0.000552565
-120,	Charge,	0.0151224,	0.000169403
-120,	C,	0.00055729,	3.68181e-05
-120,	Ca,	0.000557564,	3.68497e-05
-120,	Cl,	1.26629,	0.0384465
-120,	Mg,	1.33603,	0.0414565
-120,	Calcite,	0.000344085,	2.33848e-05
-120,	Calcite_p1,	0,	0
-120,	Dolomite,	0.00565502,	0.0070734
-120,	Dolomite_p1,	0,	0
-120,	O2g_eq,	0.133569,	0.00144506
-120,	O2g_si,	0,	0
-
-140,	ID,	0,	0
-140,	H,	0,	0
-140,	O,	0,	0
-140,	Charge,	0,	0
-140,	C,	0,	0
-140,	Ca,	0,	0
-140,	Cl,	0,	0
-140,	Mg,	0,	0
-140,	Calcite,	0,	0
-140,	Calcite_p1,	0,	0
-140,	Dolomite,	0,	0
-140,	Dolomite_p1,	0,	0
-140,	O2g_eq,	0,	0
-140,	O2g_si,	0,	0
-
-160,	ID,	0,	0
-160,	H,	0,	0
-160,	O,	0,	0
-160,	Charge,	0,	0
-160,	C,	0,	0
-160,	Ca,	0,	0
-160,	Cl,	0,	0
-160,	Mg,	0,	0
-160,	Calcite,	0,	0
-160,	Calcite_p1,	0,	0
-160,	Dolomite,	0,	0
-160,	Dolomite_p1,	0,	0
-160,	O2g_eq,	0,	0
-160,	O2g_si,	0,	0
-
-180,	ID,	0,	0
-180,	H,	4.21264e-14,	6.43691e-16
-180,	O,	0.0451131,	0.000507232
-180,	Charge,	0.0133282,	0.000158663
-180,	C,	0.000689696,	3.58042e-05
-180,	Ca,	0.000689595,	3.57904e-05
-180,	Cl,	0.920225,	0.0359217
-180,	Mg,	0.986528,	0.0403564
-180,	Calcite,	0.000420255,	2.2125e-05
-180,	Calcite_p1,	0,	0
-180,	Dolomite,	0.0148754,	0.00869066
-180,	Dolomite_p1,	0,	0
-180,	O2g_eq,	0.117169,	0.00131922
-180,	O2g_si,	0,	0
-
-200,	ID,	0,	0
-200,	H,	3.09518e-14,	5.96967e-16
-200,	O,	0.0460042,	0.000524667
-200,	Charge,	0.0127803,	0.000155455
-200,	C,	0.000624372,	3.16039e-05
-200,	Ca,	0.000624798,	3.16382e-05
-200,	Cl,	0.711112,	0.0301469
-200,	Mg,	0.766742,	0.0335527
-200,	Calcite,	0.000388114,	2.04472e-05
-200,	Calcite_p1,	0,	0
-200,	Dolomite,	0.00418652,	0.00305248
-200,	Dolomite_p1,	0,	0
-200,	O2g_eq,	0.119619,	0.00136598
-200,	O2g_si,	0,	0
-
-220,	ID,	0,	0
-220,	H,	3.38653e-14,	6.03322e-16
-220,	O,	0.0448395,	0.000519484
-220,	Charge,	0.0122,	0.000151847
-220,	C,	0.000602414,	3.03438e-05
-220,	Ca,	0.000602795,	3.0383e-05
-220,	Cl,	0.57974,	0.026681
-220,	Mg,	0.630163,	0.0304884
-220,	Calcite,	0.000371929,	1.92166e-05
-220,	Calcite_p1,	0,	0
-220,	Dolomite,	0.027979,	0.0156879
-220,	Dolomite_p1,	0,	0
-220,	O2g_eq,	0.116433,	0.00135073
-220,	O2g_si,	0,	0
-
-240,	ID,	0,	0
-240,	H,	4.01928e-14,	6.63038e-16
-240,	O,	0.0437732,	0.00051452
-240,	Charge,	0.0116593,	0.000148471
-240,	C,	0.000622541,	3.15577e-05
-240,	Ca,	0.000623301,	3.1637e-05
-240,	Cl,	0.486397,	0.0242533
-240,	Mg,	0.532535,	0.0278207
-240,	Calcite,	0.000388216,	2.05683e-05
-240,	Calcite_p1,	0,	0
-240,	Dolomite,	0.0318458,	0.0173306
-240,	Dolomite_p1,	0,	0
-240,	O2g_eq,	0.113552,	0.00133602
-240,	O2g_si,	0,	0
-
-260,	ID,	0,	0
-260,	H,	3.49348e-14,	6.34319e-16
-260,	O,	0.0426325,	0.000508416
-260,	Charge,	0.011194,	0.000145704
-260,	C,	0.000560489,	2.68571e-05
-260,	Ca,	0.000561128,	2.69044e-05
-260,	Cl,	0.416622,	0.023058
-260,	Mg,	0.459164,	0.026071
-260,	Calcite,	0.000347105,	1.71572e-05
-260,	Calcite_p1,	0,	0
-260,	Dolomite,	0.0209413,	0.012389
-260,	Dolomite_p1,	0,	0
-260,	O2g_eq,	0.110485,	0.00131874
-260,	O2g_si,	0,	0
-
-280,	ID,	0,	0
-280,	H,	3.60996e-14,	6.6192e-16
-280,	O,	0.0413622,	0.000500831
-280,	Charge,	0.0107693,	0.000142997
-280,	C,	0.000577356,	2.6604e-05
-280,	Ca,	0.000578462,	2.6685e-05
-280,	Cl,	0.360292,	0.0214539
-280,	Mg,	0.401591,	0.0253684
-280,	Calcite,	0.000364216,	1.81639e-05
-280,	Calcite_p1,	0,	0
-280,	Dolomite,	0.0221824,	0.0141538
-280,	Dolomite_p1,	0,	0
-280,	O2g_eq,	0.107019,	0.00129705
-280,	O2g_si,	0,	0
-
-300,	ID,	0,	0
-300,	H,	3.56994e-14,	6.52017e-16
-300,	O,	0.0399861,	0.000491851
-300,	Charge,	0.0103422,	0.000140284
-300,	C,	0.000604904,	2.83669e-05
-300,	Ca,	0.000605783,	2.84371e-05
-300,	Cl,	0.313853,	0.0201635
-300,	Mg,	0.353232,	0.0235415
-300,	Calcite,	0.000386924,	2.08332e-05
-300,	Calcite_p1,	0,	0
-300,	Dolomite,	0.0226333,	0.0141679
-300,	Dolomite_p1,	0,	0
-300,	O2g_eq,	0.103311,	0.00127183
-300,	O2g_si,	0,	0
-
-320,	ID,	0,	0
-320,	H,	3.40886e-14,	6.62584e-16
-320,	O,	0.0387726,	0.000484965
-320,	Charge,	0.00993204,	0.000137438
-320,	C,	0.000604705,	2.84329e-05
-320,	Ca,	0.000605988,	2.85747e-05
-320,	Cl,	0.274623,	0.0191566
-320,	Mg,	0.313395,	0.0223996
-320,	Calcite,	0.000376625,	1.86093e-05
-320,	Calcite_p1,	0,	0
-320,	Dolomite,	0.0134875,	0.0101849
-320,	Dolomite_p1,	0,	0
-320,	O2g_eq,	0.100005,	0.00125156
-320,	O2g_si,	0,	0
-
-340,	ID,	0,	0
-340,	H,	3.70111e-14,	6.80177e-16
-340,	O,	0.0376343,	0.00047866
-340,	Charge,	0.00956443,	0.000134982
-340,	C,	0.000601699,	2.71597e-05
-340,	Ca,	0.000603066,	2.72861e-05
-340,	Cl,	0.238462,	0.018479
-340,	Mg,	0.280258,	0.0219955
-340,	Calcite,	0.000385246,	1.94428e-05
-340,	Calcite_p1,	0,	0
-340,	Dolomite,	0.0318,	0.0161574
-340,	Dolomite_p1,	0,	0
-340,	O2g_eq,	0.0969628,	0.00123385
-340,	O2g_si,	0,	0
-
-360,	ID,	0,	0
-360,	H,	3.60799e-14,	6.63514e-16
-360,	O,	0.0365169,	0.000472697
-360,	Charge,	0.00920089,	0.000132444
-360,	C,	0.000589346,	2.56004e-05
-360,	Ca,	0.000590568,	2.56932e-05
-360,	Cl,	0.193642,	0.0171236
-360,	Mg,	0.248585,	0.0212179
-360,	Calcite,	0.000377203,	1.84335e-05
-360,	Calcite_p1,	0,	0
-360,	Dolomite,	0.0218774,	0.0141493
-360,	Dolomite_p1,	0,	0
-360,	O2g_eq,	0.0939613,	0.00121694
-360,	O2g_si,	0,	0
-
-380,	ID,	0,	0
-380,	H,	3.55796e-14,	6.8648e-16
-380,	O,	0.0354729,	0.000466356
-380,	Charge,	0.0088669,	0.000130032
-380,	C,	0.00063697,	2.96211e-05
-380,	Ca,	0.000638083,	2.97091e-05
-380,	Cl,	0.140909,	0.0150374
-380,	Mg,	0.216507,	0.0225187
-380,	Calcite,	0.000406558,	2.07877e-05
-380,	Calcite_p1,	0,	0
-380,	Dolomite,	0.047141,	0.021219
-380,	Dolomite_p1,	0,	0
-380,	O2g_eq,	0.0911081,	0.00119849
-380,	O2g_si,	0,	0
-
-400,	ID,	0,	0
-400,	H,	3.6232e-14,	7.02432e-16
-400,	O,	0.0345795,	0.0004623
-400,	Charge,	0.00852493,	0.000127455
-400,	C,	0.000603678,	2.58639e-05
-400,	Ca,	0.000605169,	2.59922e-05
-400,	Cl,	0.10123,	0.014209
-400,	Mg,	0.175309,	0.020241
-400,	Calcite,	0.000381044,	1.74899e-05
-400,	Calcite_p1,	0,	0
-400,	Dolomite,	0.0368562,	0.018712
-400,	Dolomite_p1,	0,	0
-400,	O2g_eq,	0.0887188,	0.0011865
-400,	O2g_si,	0,	0
-
-420,	ID,	0,	0
-420,	H,	3.83382e-14,	6.68783e-16
-420,	O,	0.0335484,	0.000455385
-420,	Charge,	0.00817171,	0.000124832
-420,	C,	0.000603529,	2.54134e-05
-420,	Ca,	0.000605869,	2.56068e-05
-420,	Cl,	0.067397,	0.0126332
-420,	Mg,	0.132868,	0.0190658
-420,	Calcite,	0.000382948,	1.74596e-05
-420,	Calcite_p1,	0,	0
-420,	Dolomite,	0.0226142,	0.0141506
-420,	Dolomite_p1,	0,	0
-420,	O2g_eq,	0.0859497,	0.00116688
-420,	O2g_si,	0,	0
-
-440,	ID,	0,	0
-440,	H,	3.67274e-14,	6.59047e-16
-440,	O,	0.0326937,	0.000451041
-440,	Charge,	0.0078609,	0.000122411
-440,	C,	0.0006084,	2.44131e-05
-440,	Ca,	0.000610448,	2.45656e-05
-440,	Cl,	0.03675,	0.0100792
-440,	Mg,	0.0979748,	0.0177069
-440,	Calcite,	0.000390181,	1.72703e-05
-440,	Calcite_p1,	0,	0
-440,	Dolomite,	0.0246411,	0.0142021
-440,	Dolomite_p1,	0,	0
-440,	O2g_eq,	0.0836881,	0.0011549
-440,	O2g_si,	0,	0
-
-460,	ID,	0,	0
-460,	H,	3.51522e-14,	6.28657e-16
-460,	O,	0.0317917,	0.000445065
-460,	Charge,	0.00753923,	0.000119952
-460,	C,	0.000598644,	2.42974e-05
-460,	Ca,	0.000600299,	2.44311e-05
-460,	Cl,	0.0143973,	0.00637336
-460,	Mg,	0.0677921,	0.0158669
-460,	Calcite,	0.0003902,	1.9257e-05
-460,	Calcite_p1,	0,	0
-460,	Dolomite,	0.0128875,	0.0100163
-460,	Dolomite_p1,	0,	0
-460,	O2g_eq,	0.0812584,	0.00113768
-460,	O2g_si,	0,	0
-
-480,	ID,	0,	0
-480,	H,	3.58576e-14,	6.04745e-16
-480,	O,	0.0309867,	0.000441054
-480,	Charge,	0.00722698,	0.000117467
-480,	C,	0.000607134,	2.35285e-05
-480,	Ca,	0.00060929,	2.37256e-05
-480,	Cl,	0.00374346,	0.00196574
-480,	Mg,	0.04042,	0.0134234
-480,	Calcite,	0.000388391,	1.64682e-05
-480,	Calcite_p1,	0,	0
-480,	Dolomite,	0.00312571,	0.000640164
-480,	Dolomite_p1,	0,	0
-480,	O2g_eq,	0.0791289,	0.00112619
-480,	O2g_si,	0,	0
-
-500,	ID,	0,	0
-500,	H,	3.64637e-14,	6.36339e-16
-500,	O,	0.0302192,	0.000436687
-500,	Charge,	0.00693139,	0.000114918
-500,	C,	0.000609576,	2.32842e-05
-500,	Ca,	0.000612185,	2.34914e-05
-500,	Cl,	0.000865377,	0.000464931
-500,	Mg,	0.017595,	0.00765859
-500,	Calcite,	0.000387095,	1.60962e-05
-500,	Calcite_p1,	0,	0
-500,	Dolomite,	0.00378159,	0.000835746
-500,	Dolomite_p1,	0,	0
-500,	O2g_eq,	0.0770378,	0.00111302
-500,	O2g_si,	0,	0
-
--- a/share/poet/barite/barite_200.qs2
+++ b/share/poet/barite/barite_200.qs2
--- a/share/poet/barite/barite_het.qs2
+++ b/share/poet/barite/barite_het.qs2
--- a/share/poet/dolo/dolo_inner_large.qs2
+++ b/share/poet/dolo/dolo_inner_large.qs2
--- a/share/poet/dolo/dolo_interp.qs2
+++ b/share/poet/dolo/dolo_interp.qs2
--- a/share/poet/surfex/PoetEGU_surfex_500.qs2
+++ b/share/poet/surfex/PoetEGU_surfex_500.qs2
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@ -33,6 +33,7 @@ add_library(POETLib
  Chemistry/SurrogateModels/HashFunctions.cpp
  Chemistry/SurrogateModels/InterpolationModule.cpp
  Chemistry/SurrogateModels/ProximityHashTable.cpp
+  Control/ControlModule.cpp
 )

 set(POET_TUG_APPROACH "Implicit" CACHE STRING "tug numerical approach to use")
--- a/src/Chemistry/ChemistryModule.hpp
+++ b/src/Chemistry/ChemistryModule.hpp
@ -2,455 +2,461 @@
 #ifndef CHEMISTRYMODULE_H_
 #define CHEMISTRYMODULE_H_

+#include "ChemistryDefs.hpp"
+#include "Control/ControlModule.hpp"
 #include "DataStructures/Field.hpp"
 #include "DataStructures/NamedVector.hpp"
-
-#include "ChemistryDefs.hpp"
-
 #include "Init/InitialList.hpp"
 #include "NameDouble.h"
+#include "PhreeqcRunner.hpp"
 #include "SurrogateModels/DHT_Wrapper.hpp"
 #include "SurrogateModels/Interpolation.hpp"

-#include "poet.hpp"
-
-#include "PhreeqcRunner.hpp"
 #include <array>
 #include <cstdint>
 #include <map>
 #include <memory>
 #include <mpi.h>
 #include <string>
+#include <unordered_set>
 #include <vector>

-namespace poet
-{
+namespace poet {
+class ControlModule;
+/**
+ * \brief Wrapper around PhreeqcRM to provide POET specific parallelization with
+ * easy access.
+ */
+class ChemistryModule {
+public:
  /**
-   * \brief Wrapper around PhreeqcRM to provide POET specific parallelization with
-   * easy access.
+   * Creates a new instance of Chemistry module with given grid cell count, work
+   * package size and communicator.
+   *
+   * This constructor shall only be called by the master. To create workers, see
+   * ChemistryModule::createWorker .
+   *
+   * When the use of parallelization is intended, the nxyz value shall be set to
+   * 1 to save memory and only one node is needed for initialization.
+   *
+   * \param nxyz Count of grid cells to allocate and initialize for each
+   * process. For parellel use set to 1 at the master.
+   * \param wp_size Count of grid cells to fill each work package at maximum.
+   * \param communicator MPI communicator to distribute work in.
   */
-  class ChemistryModule
-  {
-  public:
-    /**
-     * Creates a new instance of Chemistry module with given grid cell count, work
-     * package size and communicator.
-     *
-     * This constructor shall only be called by the master. To create workers, see
-     * ChemistryModule::createWorker .
-     *
-     * When the use of parallelization is intended, the nxyz value shall be set to
-     * 1 to save memory and only one node is needed for initialization.
-     *
-     * \param nxyz Count of grid cells to allocate and initialize for each
-     * process. For parellel use set to 1 at the master.
-     * \param wp_size Count of grid cells to fill each work package at maximum.
-     * \param communicator MPI communicator to distribute work in.
-     */
-    ChemistryModule(uint32_t wp_size,
-                    const InitialList::ChemistryInit chem_params,
-                    MPI_Comm communicator);
-
-    /**
-     * Deconstructor, which frees DHT data structure if used.
-     */
-    ~ChemistryModule();
-
-    void masterSetField(Field field);
-    /**
-     * Run the chemical simulation with parameters set.
-     */
-    void simulate(double dt);
-
-    /**
-     * Returns all known species names, including not only aqueous species, but
-     * also equilibrium, exchange, surface and kinetic reactants.
-     */
-    // auto GetPropNames() const { return this->prop_names; }
-
-    /**
-     * Return the accumulated runtime in seconds for chemical simulation.
-     */
-    auto GetChemistryTime() const { return this->chem_t; }
-
-    void setFilePadding(std::uint32_t maxiter)
-    {
-      this->file_pad =
-          static_cast<std::uint8_t>(std::ceil(std::log10(maxiter + 1)));
-    }
-
-    struct SurrogateSetup
-    {
-      std::vector<std::string> prop_names;
-      std::array<double, 2> base_totals;
-      bool has_het_ids;
-
-      bool dht_enabled;
-      std::uint32_t dht_size_mb;
-      int dht_snaps;
-      std::string dht_out_dir;
-
-      bool interp_enabled;
-      std::uint32_t interp_bucket_size;
-      std::uint32_t interp_size_mb;
-      std::uint32_t interp_min_entries;
-      bool ai_surrogate_enabled;
-    };
-
-    void masterEnableSurrogates(const SurrogateSetup &setup)
-    {
-      // FIXME: This is a hack to get the prop_names and prop_count from the setup
-      this->prop_names = setup.prop_names;
-      this->prop_count = setup.prop_names.size();
-
-      this->dht_enabled = setup.dht_enabled;
-      this->interp_enabled = setup.interp_enabled;
-      this->ai_surrogate_enabled = setup.ai_surrogate_enabled;
-
-      this->base_totals = setup.base_totals;
-
-      if (this->dht_enabled || this->interp_enabled)
-      {
-        this->initializeDHT(setup.dht_size_mb, this->params.dht_species,
-                            setup.has_het_ids);
-
-        if (setup.dht_snaps != DHT_SNAPS_DISABLED)
-        {
-          this->setDHTSnapshots(setup.dht_snaps, setup.dht_out_dir);
-        }
-      }
-
-      if (this->interp_enabled)
-      {
-        this->initializeInterp(setup.interp_bucket_size, setup.interp_size_mb,
-                               setup.interp_min_entries,
-                               this->params.interp_species);
-      }
-    }
-
-    /**
-     * Intended to alias input parameters for grid initialization with a single
-     * value per species.
-     */
-    using SingleCMap = std::unordered_map<std::string, double>;
-
-    /**
-     * Intended to alias input parameters for grid initialization with mutlitple
-     * values per species.
-     */
-    using VectorCMap = std::unordered_map<std::string, std::vector<double>>;
-
-    /**
-     * Enumerating DHT file options
-     */
-    enum
-    {
-      DHT_SNAPS_DISABLED = 0, //!< disabled file output
-      DHT_SNAPS_SIMEND,       //!< only output of snapshot after simulation
-      DHT_SNAPS_ITEREND       //!< output snapshots after each iteration
-    };
-
-    /**
-     * **Only called by workers!** Start the worker listening loop.
-     */
-    void WorkerLoop();
-
-    /**
-     * **Called by master** Advise the workers to break the loop.
-     */
-    void MasterLoopBreak();
-
-    /**
-     * **Master only** Return count of grid cells.
-     */
-    auto GetNCells() const { return this->n_cells; }
-    /**
-     * **Master only** Return work package size.
-     */
-    auto GetWPSize() const { return this->wp_size; }
-    /**
-     * **Master only** Return the time in seconds the master spent waiting for any
-     * free worker.
-     */
-    auto GetMasterIdleTime() const { return this->idle_t; }
-    /**
-     * **Master only** Return the time in seconds the master spent in sequential
-     * part of the simulation, including times for shuffling/unshuffling field
-     * etc.
-     */
-    auto GetMasterSequentialTime() const { return this->seq_t; }
-    /**
-     * **Master only** Return the time in seconds the master spent in the
-     * send/receive loop.
-     */
-    auto GetMasterLoopTime() const { return this->send_recv_t; }
-
-    /**
-     * **Master only** Collect and return all accumulated timings recorded by
-     * workers to run Phreeqc simulation.
-     *
-     * \return Vector of all accumulated Phreeqc timings.
-     */
-    std::vector<double> GetWorkerPhreeqcTimings() const;
-    /**
-     * **Master only** Collect and return all accumulated timings recorded by
-     * workers to get values from the DHT.
-     *
-     * \return Vector of all accumulated DHT get times.
-     */
-    std::vector<double> GetWorkerDHTGetTimings() const;
-    /**
-     * **Master only** Collect and return all accumulated timings recorded by
-     * workers to write values to the DHT.
-     *
-     * \return Vector of all accumulated DHT fill times.
-     */
-    std::vector<double> GetWorkerDHTFillTimings() const;
-    /**
-     * **Master only** Collect and return all accumulated timings recorded by
-     * workers waiting for work packages from the master.
-     *
-     * \return Vector of all accumulated waiting times.
-     */
-    std::vector<double> GetWorkerIdleTimings() const;
-
-    /**
-     * **Master only** Collect and return DHT hits of all workers.
-     *
-     * \return Vector of all count of DHT hits.
-     */
-    std::vector<uint32_t> GetWorkerDHTHits() const;
-
-    /**
-     * **Master only** Collect and return DHT evictions of all workers.
-     *
-     * \return Vector of all count of DHT evictions.
-     */
-    std::vector<uint32_t> GetWorkerDHTEvictions() const;
-
-    /**
-     * **Master only** Returns the current state of the chemical field.
-     *
-     * \return Reference to the chemical field.
-     */
-    Field &getField() { return this->chem_field; }
-
-    /**
-     * **Master only** Enable/disable progress bar.
-     *
-     * \param enabled True if print progressbar, false if not.
-     */
-    void setProgressBarPrintout(bool enabled)
-    {
-      this->print_progessbar = enabled;
-    };
-
-    /**
-     *  **Master only** Set the ai surrogate validity vector from R
-     */
-    void set_ai_surrogate_validity_vector(std::vector<int> r_vector);
-
-    std::vector<uint32_t> GetWorkerInterpolationCalls() const;
-
-    std::vector<double> GetWorkerInterpolationWriteTimings() const;
-    std::vector<double> GetWorkerInterpolationReadTimings() const;
-    std::vector<double> GetWorkerInterpolationGatherTimings() const;
-    std::vector<double> GetWorkerInterpolationFunctionCallTimings() const;
-
-    std::vector<uint32_t> GetWorkerPHTCacheHits() const;
-
-    std::vector<int> ai_surrogate_validity_vector;
-
-    RuntimeParameters *runtime_params = nullptr;
-    uint32_t control_iteration_counter = 0;
-
-    struct error_stats
-    {
-      std::vector<double> mape;
-      std::vector<double> rrsme;
-      uint32_t iteration;
-
-      error_stats(size_t species_count, size_t iter)
-      : mape(species_count, 0.0), rrsme(species_count, 0.0), iteration(iter) {}
-    };
-
-    std::vector<error_stats> error_stats_history;
-
-    static void computeStats(const std::vector<double> &pqc_vector,
-                      const std::vector<double> &sur_vector,
-                      uint32_t size_per_prop, uint32_t species_count,
-                      error_stats &stats);
-                      
-  protected:
-    void initializeDHT(uint32_t size_mb,
-                       const NamedVector<std::uint32_t> &key_species,
-                       bool has_het_ids);
-    void setDHTSnapshots(int type, const std::string &out_dir);
-    void setDHTReadFile(const std::string &input_file);
-
-    void initializeInterp(std::uint32_t bucket_size, std::uint32_t size_mb,
-                          std::uint32_t min_entries,
-                          const NamedVector<std::uint32_t> &key_species);
-
-    enum
-    {
-      CHEM_FIELD_INIT,
-      CHEM_DHT_ENABLE,
-      CHEM_DHT_SIGNIF_VEC,
-      CHEM_DHT_SNAPS,
-      CHEM_DHT_READ_FILE,
-      CHEM_INTERP,
-      CHEM_IP_ENABLE,
-      CHEM_IP_MIN_ENTRIES,
-      CHEM_IP_SIGNIF_VEC,
-      CHEM_WORK_LOOP,
-      CHEM_PERF,
-      CHEM_BREAK_MAIN_LOOP,
-      CHEM_AI_BCAST_VALIDITY
-    };
-
-    enum
-    {
-      LOOP_WORK,
-      LOOP_END,
-      LOOP_CTRL
-    };
-
-    enum
-    {
-      WORKER_PHREEQC,
-      WORKER_DHT_GET,
-      WORKER_DHT_FILL,
-      WORKER_IDLE,
-      WORKER_IP_WRITE,
-      WORKER_IP_READ,
-      WORKER_IP_GATHER,
-      WORKER_IP_FC,
-      WORKER_DHT_HITS,
-      WORKER_DHT_EVICTIONS,
-      WORKER_PHT_CACHE_HITS,
-      WORKER_IP_CALLS
-    };
-
-    std::vector<uint32_t> interp_calls;
-    std::vector<uint32_t> dht_hits;
-    std::vector<uint32_t> dht_evictions;
-
-    struct worker_s
-    {
-      double phreeqc_t = 0.;
-      double dht_get = 0.;
-      double dht_fill = 0.;
-      double idle_t = 0.;
-    };
-
-    struct worker_info_s
-    {
-      char has_work = 0;
-      double *send_addr;
-      double *surrogate_addr;
-    };
-
-    using worker_list_t = std::vector<struct worker_info_s>;
-    using workpointer_t = std::vector<double>::iterator;
-
-    void MasterRunParallel(double dt);
-    void MasterRunSequential();
-
-    void MasterSendPkgs(worker_list_t &w_list, workpointer_t &work_pointer, workpointer_t &sur_pointer,
-                        int &pkg_to_send, int &count_pkgs, int &free_workers,
-                        double dt, uint32_t iteration, uint32_t control_iteration,
-                        const std::vector<uint32_t> &wp_sizes_vector);
-    void MasterRecvPkgs(worker_list_t &w_list, int &pkg_to_recv, bool to_send,
-                        int &free_workers);
-
-    std::vector<double> MasterGatherWorkerTimings(int type) const;
-    std::vector<uint32_t> MasterGatherWorkerMetrics(int type) const;
-
-    void WorkerProcessPkgs(struct worker_s &timings, uint32_t &iteration);
-
-    void WorkerDoWork(MPI_Status &probe_status, int double_count,
-                      struct worker_s &timings);
-    void WorkerPostIter(MPI_Status &prope_status, uint32_t iteration);
-    void WorkerPostSim(uint32_t iteration);
-
-    void WorkerWriteDHTDump(uint32_t iteration);
-    void WorkerReadDHTDump(const std::string &dht_input_file);
-
-    void WorkerPerfToMaster(int type, const struct worker_s &timings);
-    void WorkerMetricsToMaster(int type);
-
-    void WorkerRunWorkPackage(WorkPackage &work_package, double dSimTime,
-                              double dTimestep);
-
-    std::vector<uint32_t> CalculateWPSizesVector(uint32_t n_cells,
-                                                 uint32_t wp_size) const;
-    std::vector<double> shuffleField(const std::vector<double> &in_field,
-                                     uint32_t size_per_prop, uint32_t prop_count,
-                                     uint32_t wp_count);
-    void unshuffleField(const std::vector<double> &in_buffer,
-                        uint32_t size_per_prop, uint32_t prop_count,
-                        uint32_t wp_count, std::vector<double> &out_field);
-    std::vector<std::int32_t>
-    parseDHTSpeciesVec(const NamedVector<std::uint32_t> &key_species,
-                       const std::vector<std::string> &to_compare) const;
-
-    void BCastStringVec(std::vector<std::string> &io);
-
-    int comm_size, comm_rank;
-    MPI_Comm group_comm;
-
-    bool is_sequential;
-    bool is_master;
-
-    uint32_t wp_size;
-    bool dht_enabled{false};
-    int dht_snaps_type{DHT_SNAPS_DISABLED};
-    std::string dht_file_out_dir;
-
-    poet::DHT_Wrapper *dht = nullptr;
-
-    bool interp_enabled{false};
-    std::unique_ptr<poet::InterpolationModule> interp;
-
-    bool ai_surrogate_enabled{false};
-
-    static constexpr uint32_t BUFFER_OFFSET = 6;
-
-    inline void ChemBCast(void *buf, int count, MPI_Datatype datatype) const
-    {
-      MPI_Bcast(buf, count, datatype, 0, this->group_comm);
-    }
-
-    inline void PropagateFunctionType(int &type) const
-    {
-      ChemBCast(&type, 1, MPI_INT);
-    }
-    double simtime = 0.;
-    double idle_t = 0.;
-    double seq_t = 0.;
-    double send_recv_t = 0.;
-
-    std::array<double, 2> base_totals{0};
-
-    bool print_progessbar{false};
-
-    std::uint8_t file_pad{1};
-
-    double chem_t{0.};
-
-    uint32_t n_cells = 0;
-    uint32_t prop_count = 0;
+  ChemistryModule(uint32_t wp_size,
+                  const InitialList::ChemistryInit chem_params,
+                  MPI_Comm communicator);
+
+  /**
+   * Deconstructor, which frees DHT data structure if used.
+   */
+  ~ChemistryModule();
+
+  void masterSetField(Field field);
+  /**
+   * Run the chemical simulation with parameters set.
+   */
+  void simulate(double dt);
+
+  /**
+   * Returns all known species names, including not only aqueous species, but
+   * also equilibrium, exchange, surface and kinetic reactants.
+   */
+  // auto GetPropNames() const { return this->prop_names; }
+
+  /**
+   * Return the accumulated runtime in seconds for chemical simulation.
+   */
+  auto GetChemistryTime() const { return this->chem_t; }
+
+  void setFilePadding(std::uint32_t maxiter) {
+    this->file_pad =
+        static_cast<std::uint8_t>(std::ceil(std::log10(maxiter + 1)));
+  }
+
+  struct SurrogateSetup {
    std::vector<std::string> prop_names;
+    std::array<double, 2> base_totals;
+    bool has_het_ids;

-    Field chem_field;
+    bool dht_enabled;
+    std::uint32_t dht_size_mb;
+    int dht_snaps;
+    std::string dht_out_dir;

-    const InitialList::ChemistryInit params;
-
-    std::unique_ptr<PhreeqcRunner> pqc_runner;
-
-    std::vector<double> sur_shuffled;
+    bool interp_enabled;
+    std::uint32_t interp_bucket_size;
+    std::uint32_t interp_size_mb;
+    std::uint32_t interp_min_entries;
+    bool ai_surrogate_enabled;
  };
+
+  void masterEnableSurrogates(const SurrogateSetup &setup) {
+    // FIXME: This is a hack to get the prop_names and prop_count from the setup
+    this->prop_names = setup.prop_names;
+    this->prop_count = setup.prop_names.size();
+
+    this->dht_enabled = setup.dht_enabled;
+    this->interp_enabled = setup.interp_enabled;
+    this->ai_surrogate_enabled = setup.ai_surrogate_enabled;
+
+    this->base_totals = setup.base_totals;
+
+    if (this->dht_enabled || this->interp_enabled) {
+      this->initializeDHT(setup.dht_size_mb, this->params.dht_species,
+                          setup.has_het_ids);
+
+      if (setup.dht_snaps != DHT_SNAPS_DISABLED) {
+        this->setDHTSnapshots(setup.dht_snaps, setup.dht_out_dir);
+      }
+    }
+
+    if (this->interp_enabled) {
+      this->initializeInterp(setup.interp_bucket_size, setup.interp_size_mb,
+                             setup.interp_min_entries,
+                             this->params.interp_species);
+    }
+  }
+
+  /**
+   * Intended to alias input parameters for grid initialization with a single
+   * value per species.
+   */
+  using SingleCMap = std::unordered_map<std::string, double>;
+
+  /**
+   * Intended to alias input parameters for grid initialization with mutlitple
+   * values per species.
+   */
+  using VectorCMap = std::unordered_map<std::string, std::vector<double>>;
+
+  /**
+   * Enumerating DHT file options
+   */
+  enum {
+    DHT_SNAPS_DISABLED = 0, //!< disabled file output
+    DHT_SNAPS_SIMEND,       //!< only output of snapshot after simulation
+    DHT_SNAPS_ITEREND       //!< output snapshots after each iteration
+  };
+
+  /**
+   * **Only called by workers!** Start the worker listening loop.
+   */
+  void WorkerLoop();
+
+  /**
+   * **Called by master** Advise the workers to break the loop.
+   */
+  void MasterLoopBreak();
+
+  /**
+   * **Master only** Return count of grid cells.
+   */
+  auto GetNCells() const { return this->n_cells; }
+  /**
+   * **Master only** Return work package size.
+   */
+  auto GetWPSize() const { return this->wp_size; }
+  /**
+   * **Master only** Return the time in seconds the master spent waiting for any
+   * free worker.
+   */
+  auto GetMasterIdleTime() const { return this->idle_t; }
+  /**
+   * **Master only** Return the time in seconds the master spent in sequential
+   * part of the simulation, including times for shuffling/unshuffling field
+   * etc.
+   */
+  auto GetMasterSequentialTime() const { return this->seq_t; }
+  /**
+   * **Master only** Return the time in seconds the master spent in the
+   * send/receive loop.
+   */
+  auto GetMasterLoopTime() const { return this->send_recv_t; }
+
+  auto GetMasterRecvCtrlDataTime() const { return this->recv_ctrl_t; }
+
+  auto GetMasterUnshuffleTime() const { return this->shuf_t; }
+
+  auto GetMasterCtrlMetricsTime() const { return this->metrics_t; }
+
+  /**
+   * **Master only** Collect and return all accumulated timings recorded by
+   * workers to run Phreeqc simulation.
+   *
+   * \return Vector of all accumulated Phreeqc timings.
+   */
+  std::vector<double> GetWorkerPhreeqcTimings() const;
+  /**
+   * **Master only** Collect and return all accumulated timings recorded by
+   * workers to get values from the DHT.
+   *
+   * \return Vector of all accumulated DHT get times.
+   */
+  std::vector<double> GetWorkerDHTGetTimings() const;
+  /**
+   * **Master only** Collect and return all accumulated timings recorded by
+   * workers to write values to the DHT.
+   *
+   * \return Vector of all accumulated DHT fill times.
+   */
+  std::vector<double> GetWorkerDHTFillTimings() const;
+  /**
+   * **Master only** Collect and return all accumulated timings recorded by
+   * workers waiting for work packages from the master.
+   *
+   * \return Vector of all accumulated waiting times.
+   */
+  std::vector<double> GetWorkerIdleTimings() const;
+
+  std::vector<double> GetWorkerControlTimings() const;
+
+  /**
+   * **Master only** Collect and return DHT hits of all workers.
+   *
+   * \return Vector of all count of DHT hits.
+   */
+  std::vector<uint32_t> GetWorkerDHTHits() const;
+
+  /**
+   * **Master only** Collect and return DHT evictions of all workers.
+   *
+   * \return Vector of all count of DHT evictions.
+   */
+  std::vector<uint32_t> GetWorkerDHTEvictions() const;
+
+  /**
+   * **Master only** Returns the current state of the chemical field.
+   *
+   * \return Reference to the chemical field.
+   */
+  Field &getField() { return this->chem_field; }
+
+  /**
+   * **Master only** Enable/disable progress bar.
+   *
+   * \param enabled True if print progressbar, false if not.
+   */
+  void setProgressBarPrintout(bool enabled) {
+    this->print_progessbar = enabled;
+  };
+
+  /**
+   *  **Master only** Set the ai surrogate validity vector from R
+   */
+  void set_ai_surrogate_validity_vector(std::vector<int> r_vector);
+
+  std::vector<uint32_t> GetWorkerInterpolationCalls() const;
+
+  std::vector<double> GetWorkerInterpolationWriteTimings() const;
+  std::vector<double> GetWorkerInterpolationReadTimings() const;
+  std::vector<double> GetWorkerInterpolationGatherTimings() const;
+  std::vector<double> GetWorkerInterpolationFunctionCallTimings() const;
+
+  std::vector<uint32_t> GetWorkerPHTCacheHits() const;
+
+  std::vector<int> ai_surrogate_validity_vector;
+
+  void SetControlModule(poet::ControlModule *ctrl) { control_module = ctrl; }
+
+  void SetDhtEnabled(bool enabled) { this->dht_enabled = enabled; }
+  bool GetDhtEnabled() const { return this->dht_enabled; }
+
+  void SetInterpEnabled(bool enabled) { this->interp_enabled = enabled; }
+  bool GetInterpEnabled() const { return interp_enabled; }
+
+  void SetWarmupEnabled(bool enabled) { this->warmup_enabled = enabled; }
+
+  void SetControlCellIds(const std::vector<uint32_t> &ids) {
+    this->ctrl_cell_ids = std::unordered_set<uint32_t>(ids.begin(), ids.end());
+  }
+
+protected:
+  void initializeDHT(uint32_t size_mb,
+                     const NamedVector<std::uint32_t> &key_species,
+                     bool has_het_ids);
+  void setDHTSnapshots(int type, const std::string &out_dir);
+  void setDHTReadFile(const std::string &input_file);
+
+  void initializeInterp(std::uint32_t bucket_size, std::uint32_t size_mb,
+                        std::uint32_t min_entries,
+                        const NamedVector<std::uint32_t> &key_species);
+
+  enum {
+    CHEM_FIELD_INIT,
+    // CHEM_DHT_ENABLE,
+    CHEM_DHT_SIGNIF_VEC,
+    CHEM_DHT_SNAPS,
+    CHEM_DHT_READ_FILE,
+    // CHEM_WARMUP_PHASE,  // Control flag
+    // CHEM_CTRL_ENABLE, // Control flag
+    // CHEM_IP_ENABLE,
+    CHEM_IP_MIN_ENTRIES,
+    CHEM_IP_SIGNIF_VEC,
+    CHEM_WORK_LOOP,
+    CHEM_PERF,
+    CHEM_BREAK_MAIN_LOOP,
+    CHEM_AI_BCAST_VALIDITY
+  };
+
+  enum { LOOP_WORK, LOOP_END, LOOP_CTRL };
+
+  enum {
+    WORKER_PHREEQC,
+    WORKER_CTRL_ITER,
+    WORKER_DHT_GET,
+    WORKER_DHT_FILL,
+    WORKER_IDLE,
+    WORKER_IP_WRITE,
+    WORKER_IP_READ,
+    WORKER_IP_GATHER,
+    WORKER_IP_FC,
+    WORKER_DHT_HITS,
+    WORKER_DHT_EVICTIONS,
+    WORKER_PHT_CACHE_HITS,
+    WORKER_IP_CALLS
+  };
+
+  std::vector<uint32_t> interp_calls;
+  std::vector<uint32_t> dht_hits;
+  std::vector<uint32_t> dht_evictions;
+
+  struct worker_s {
+    double phreeqc_t = 0.;
+    double dht_get = 0.;
+    double dht_fill = 0.;
+    double idle_t = 0.;
+    double ctrl_t = 0.;
+    double ctrl_phreeqc_t = 0.;
+  };
+
+  struct worker_info_s {
+    char has_work = 0;
+    double *send_addr;
+    double *surrogate_addr;
+  };
+
+  using worker_list_t = std::vector<struct worker_info_s>;
+  using workpointer_t = std::vector<double>::iterator;
+
+  void MasterRunParallel(double dt);
+  void MasterRunSequential();
+
+  void MasterSendPkgs(worker_list_t &w_list, workpointer_t &work_pointer,
+                      workpointer_t &sur_pointer, int &pkg_to_send,
+                      int &count_pkgs, int &free_workers, double dt,
+                      uint32_t iteration,
+                      const std::vector<uint32_t> &wp_sizes_vector);
+  void MasterRecvPkgs(worker_list_t &w_list, int &pkg_to_recv, bool to_send,
+                      int &free_workers);
+
+  std::vector<double> MasterGatherWorkerTimings(int type) const;
+  std::vector<uint32_t> MasterGatherWorkerMetrics(int type) const;
+
+  void WorkerProcessPkgs(struct worker_s &timings, uint32_t &iteration);
+
+  void WorkerDoWork(MPI_Status &probe_status, int double_count,
+                    struct worker_s &timings);
+  void WorkerPostIter(MPI_Status &prope_status, uint32_t iteration);
+  void WorkerPostSim(uint32_t iteration);
+
+  void WorkerWriteDHTDump(uint32_t iteration);
+  void WorkerReadDHTDump(const std::string &dht_input_file);
+
+  void WorkerPerfToMaster(int type, const struct worker_s &timings);
+  void WorkerMetricsToMaster(int type);
+
+  void WorkerRunWorkPackage(WorkPackage &work_package, double dSimTime,
+                            double dTimestep);
+
+  void ProcessControlWorkPackage(std::vector<std::vector<double>> &input,
+                                 double current_sim_time, double dt,
+                                 struct worker_s &timings);
+
+  std::vector<uint32_t> CalculateWPSizesVector(uint32_t n_cells,
+                                               uint32_t wp_size) const;
+  std::vector<double> shuffleField(const std::vector<double> &in_field,
+                                   uint32_t size_per_prop, uint32_t prop_count,
+                                   uint32_t wp_count);
+  void unshuffleField(const std::vector<double> &in_buffer,
+                      uint32_t size_per_prop, uint32_t prop_count,
+                      uint32_t wp_count, std::vector<double> &out_field);
+  std::vector<std::int32_t>
+  parseDHTSpeciesVec(const NamedVector<std::uint32_t> &key_species,
+                     const std::vector<std::string> &to_compare) const;
+
+  void BCastStringVec(std::vector<std::string> &io);
+
+  int packResultsIntoBuffer(std::vector<double> &mpi_buffer, int base_count,
+                            const WorkPackage &wp,
+                            const WorkPackage &wp_control);
+
+  int comm_size, comm_rank;
+  MPI_Comm group_comm;
+
+  bool is_sequential;
+  bool is_master;
+
+  uint32_t wp_size;
+  bool dht_enabled{false};
+  int dht_snaps_type{DHT_SNAPS_DISABLED};
+  std::string dht_file_out_dir;
+
+  poet::DHT_Wrapper *dht = nullptr;
+
+  bool interp_enabled{false};
+  std::unique_ptr<poet::InterpolationModule> interp;
+
+  bool ai_surrogate_enabled{false};
+
+  static constexpr uint32_t BUFFER_OFFSET = 6;
+
+  inline void ChemBCast(void *buf, int count, MPI_Datatype datatype) const {
+    MPI_Bcast(buf, count, datatype, 0, this->group_comm);
+  }
+
+  inline void PropagateFunctionType(int &type) const {
+    ChemBCast(&type, 1, MPI_INT);
+  }
+
+  double simtime = 0.;
+  double idle_t = 0.;
+  double seq_t = 0.;
+  double send_recv_t = 0.;
+
+  double recv_ctrl_t = 0.;
+  double shuf_t = 0.;
+  double metrics_t = 0.;
+
+  std::array<double, 2> base_totals{0};
+
+  bool print_progessbar{false};
+
+  std::uint8_t file_pad{1};
+
+  double chem_t{0.};
+
+  uint32_t n_cells = 0;
+  uint32_t prop_count = 0;
+  std::vector<std::string> prop_names;
+
+  Field chem_field;
+
+  const InitialList::ChemistryInit params;
+
+  std::unique_ptr<PhreeqcRunner> pqc_runner;
+
+  poet::ControlModule *control_module = nullptr;
+
+  std::vector<double> mpi_surr_buffer;
+
+  bool control_enabled{false};
+  bool warmup_enabled{false};
+
+  std::unordered_set<uint32_t> ctrl_cell_ids;
+  std::vector<std::vector<double>> control_batch;
+};
 } // namespace poet

 #endif // CHEMISTRYMODULE_H_
--- a/src/Chemistry/MasterFunctions.cpp
+++ b/src/Chemistry/MasterFunctions.cpp
@ -3,7 +3,6 @@
 #include <algorithm>
 #include <cstddef>
 #include <cstdint>
-#include <iomanip>
 #include <mpi.h>
 #include <vector>

@ -41,6 +40,12 @@ std::vector<double> poet::ChemistryModule::GetWorkerPhreeqcTimings() const {
  return MasterGatherWorkerTimings(WORKER_PHREEQC);
 }

+std::vector<double> poet::ChemistryModule::GetWorkerControlTimings() const {
+  int type = CHEM_PERF;
+  MPI_Bcast(&type, 1, MPI_INT, 0, this->group_comm);
+  return MasterGatherWorkerTimings(WORKER_CTRL_ITER);
+}
+
 std::vector<double> poet::ChemistryModule::GetWorkerDHTGetTimings() const {
  int type = CHEM_PERF;
  MPI_Bcast(&type, 1, MPI_INT, 0, this->group_comm);
@ -160,39 +165,6 @@ std::vector<uint32_t> poet::ChemistryModule::GetWorkerPHTCacheHits() const {
  return ret;
 }

-void poet::ChemistryModule::computeStats(const std::vector<double> &pqc_vector,
-                                         const std::vector<double> &sur_vector,
-                                         uint32_t size_per_prop,
-                                         uint32_t species_count,
-                                         error_stats &stats) {
-  for (uint32_t i = 0; i < species_count; ++i) {
-    double err_sum = 0.0;
-    double sqr_err_sum = 0.0;
-    uint32_t base_idx = i * size_per_prop;
-
-    for (uint32_t j = 0; j < size_per_prop; ++j) {
-      const double pqc_value = pqc_vector[base_idx + j];
-      const double sur_value = sur_vector[base_idx + j];
-
-      if (pqc_value == 0.0) {
-        if (sur_value != 0.0) {
-          err_sum += 1.0;
-          sqr_err_sum += 1.0;
-        }
-        // Both zero: skip
-      } else {
-        double alpha = 1.0 - (sur_value / pqc_value);
-        err_sum += std::abs(alpha);
-        sqr_err_sum += alpha * alpha;
-      }
-    }
-
-    stats.mape[i] = 100.0 * (err_sum / size_per_prop);
-    stats.rrsme[i] =
-        (size_per_prop > 0) ? std::sqrt(sqr_err_sum / size_per_prop) : 0.0;
-  }
-}
-
 inline std::vector<int> shuffleVector(const std::vector<int> &in_vector,
                                      uint32_t size_per_prop,
                                      uint32_t wp_count) {
@ -264,7 +236,7 @@ inline void poet::ChemistryModule::MasterSendPkgs(
    worker_list_t &w_list, workpointer_t &work_pointer,
    workpointer_t &sur_pointer, int &pkg_to_send, int &count_pkgs,
    int &free_workers, double dt, uint32_t iteration,
-    uint32_t control_iteration, const std::vector<uint32_t> &wp_sizes_vector) {
+    const std::vector<uint32_t> &wp_sizes_vector) {
  /* declare variables */
  int local_work_package_size;

@ -278,9 +250,15 @@ inline void poet::ChemistryModule::MasterSendPkgs(
      local_work_package_size = (int)wp_sizes_vector[count_pkgs];
      count_pkgs++;

+      uint32_t wp_start_index =
+          std::accumulate(wp_sizes_vector.begin(),
+                          std::next(wp_sizes_vector.begin(), count_pkgs), 0);
+
      /* note current processed work package in workerlist */
      w_list[p].send_addr = work_pointer.base();
      w_list[p].surrogate_addr = sur_pointer.base();
+      // this->control_enabled ? sur_pointer.base() : w_list[p].surrogate_addr =
+      // nullptr;

      /* push work pointer to next work package */
      const uint32_t end_of_wp = local_work_package_size * this->prop_count;
@ -300,12 +278,12 @@ inline void poet::ChemistryModule::MasterSendPkgs(
      // current time of simulation (age) in seconds
      send_buffer[end_of_wp + 3] = this->simtime;
      // current work package start location in field
-      uint32_t wp_start_index =
-          std::accumulate(wp_sizes_vector.begin(),
-                          std::next(wp_sizes_vector.begin(), count_pkgs), 0);
      send_buffer[end_of_wp + 4] = wp_start_index;
-      // whether this iteration is a control iteration
-      send_buffer[end_of_wp + 5] = control_iteration;
+      // control flags (bitmask)
+      int flags = (this->interp_enabled ? 1 : 0) | (this->dht_enabled ? 2 : 0) |
+                  (this->warmup_enabled ? 4 : 0) |
+                  (this->control_enabled ? 8 : 0);
+      send_buffer[end_of_wp + 5] = static_cast<double>(flags);

      /* ATTENTION Worker p has rank p+1 */
      // MPI_Send(send_buffer, end_of_wp + BUFFER_OFFSET, MPI_DOUBLE, p + 1,
@ -328,8 +306,11 @@ inline void poet::ChemistryModule::MasterRecvPkgs(worker_list_t &w_list,
  /* declare most of the variables here */
  int need_to_receive = 1;
  double idle_a, idle_b;
+  double recv_ctrl_a, recv_ctrl_b;
  int p, size;
+  std::vector<double> recv_buffer;

+  recv_buffer.reserve(wp_size * prop_count * 2);
  MPI_Status probe_status;
  // master_recv_a = MPI_Wtime();
  /* start to loop as long there are packages to recv and the need to receive
@ -347,38 +328,52 @@ inline void poet::ChemistryModule::MasterRecvPkgs(worker_list_t &w_list,
      idle_b = MPI_Wtime();
      this->idle_t += idle_b - idle_a;
    }
-
+    if (!need_to_receive) {
+      continue;
+    }
    /* if need_to_receive was set to true above, so there is a message to
     * receive */
-    if (need_to_receive) {
-      p = probe_status.MPI_SOURCE;
-      if (probe_status.MPI_TAG == LOOP_WORK) {
-        MPI_Get_count(&probe_status, MPI_DOUBLE, &size);
-        MPI_Recv(w_list[p - 1].send_addr, size, MPI_DOUBLE, p, LOOP_WORK,
-                 this->group_comm, MPI_STATUS_IGNORE);
-        w_list[p - 1].has_work = 0;
-        pkg_to_recv -= 1;
-        free_workers++;
-      }
-      if (probe_status.MPI_TAG == LOOP_CTRL) {
-        MPI_Get_count(&probe_status, MPI_DOUBLE, &size);
-
-        // layout of buffer is [phreeqc][surrogate]
-        std::vector<double> recv_buffer(size);
-
-        MPI_Recv(recv_buffer.data(), size, MPI_DOUBLE, p, LOOP_CTRL,
-                 this->group_comm, MPI_STATUS_IGNORE);
-
-        std::copy(recv_buffer.begin(), recv_buffer.begin() + (size / 2),
-                  w_list[p - 1].send_addr);
-
-        std::copy(recv_buffer.begin() + (size / 2), recv_buffer.begin() + size,
-                  w_list[p - 1].surrogate_addr);
-
-        w_list[p - 1].has_work = 0;
-        pkg_to_recv -= 1;
-        free_workers++;
+    p = probe_status.MPI_SOURCE;
+    bool handled = false;
+
+    switch (probe_status.MPI_TAG) {
+    case LOOP_WORK: {
+      MPI_Get_count(&probe_status, MPI_DOUBLE, &size);
+      MPI_Recv(w_list[p - 1].send_addr, size, MPI_DOUBLE, p, LOOP_WORK,
+               this->group_comm, MPI_STATUS_IGNORE);
+      // Only LOOP_WORK completes a work package
+      w_list[p - 1].has_work = 0;
+      pkg_to_recv -= 1;
+      free_workers++;
+      break;
+    }
+    case LOOP_CTRL: {
+      recv_ctrl_a = MPI_Wtime();
+
+      MPI_Get_count(&probe_status, MPI_DOUBLE, &size);
+      recv_buffer.resize(size);
+      MPI_Recv(recv_buffer.data(), size, MPI_DOUBLE, p, LOOP_CTRL,
+               this->group_comm, MPI_STATUS_IGNORE);
+
+      recv_ctrl_b = MPI_Wtime();
+      recv_ctrl_t += recv_ctrl_b - recv_ctrl_a;
+
+      // Collect PHREEQC rows for control cells
+      const std::size_t cells_per_batch =
+          static_cast<std::size_t>(size) /
+          static_cast<std::size_t>(this->prop_count);
+      for (std::size_t i = 0; i < cells_per_batch; i++) {
+        std::vector<double> cell_output(
+            recv_buffer.begin() + this->prop_count * i,
+            recv_buffer.begin() + this->prop_count * (i + 1));
+        this->control_batch.push_back(std::move(cell_output));
      }
+      break;
+    }
+    default: {
+      throw std::runtime_error("Master received unknown MPI tag: " +
+                               std::to_string(probe_status.MPI_TAG));
+    }
    }
  }
 }
@ -431,6 +426,7 @@ void poet::ChemistryModule::MasterRunParallel(double dt) {
  int free_workers;
  int i_pkgs;
  int ftype;
+  double shuf_a, shuf_b, metrics_a, metrics_b;

  const std::vector<uint32_t> wp_sizes_vector =
      CalculateWPSizesVector(this->n_cells, this->wp_size);
@ -448,47 +444,27 @@ void poet::ChemistryModule::MasterRunParallel(double dt) {
  ftype = CHEM_WORK_LOOP;
  PropagateFunctionType(ftype);

-  ftype = CHEM_INTERP;
-  PropagateFunctionType(ftype);
-
-  if(this->runtime_params->rollback_simulation){
-    this->interp_enabled = false;
-    int interp_flag = 0;
-    ChemBCast(&interp_flag, 1, MPI_INT);
-  } else {
-    this->interp_enabled = true;
-    int interp_flag = 1;
-    ChemBCast(&interp_flag, 1, MPI_INT);
-  }
-
  MPI_Barrier(this->group_comm);

  static uint32_t iteration = 0;
-  uint32_t control_iteration = static_cast<uint32_t>(
-      this->runtime_params->control_iteration_active ? 1 : 0);
-  if (control_iteration) {
-    sur_shuffled.clear();
-    sur_shuffled.reserve(this->n_cells * this->prop_count);
-  }

  /* start time measurement of sequential part */
  seq_a = MPI_Wtime();

  /* shuffle grid */
  // grid.shuffleAndExport(mpi_buffer);
-
  std::vector<double> mpi_buffer =
      shuffleField(chem_field.AsVector(), this->n_cells, this->prop_count,
                   wp_sizes_vector.size());

-  this->sur_shuffled.resize(mpi_buffer.size());
+  // this->mpi_surr_buffer.resize(mpi_buffer.size());

  /* setup local variables */
  pkg_to_send = wp_sizes_vector.size();
  pkg_to_recv = wp_sizes_vector.size();

  workpointer_t work_pointer = mpi_buffer.begin();
-  workpointer_t sur_pointer = sur_shuffled.begin();
+  workpointer_t sur_pointer = mpi_buffer.begin();
  worker_list_t worker_list(this->comm_size - 1);

  free_workers = this->comm_size - 1;
@ -512,8 +488,7 @@ void poet::ChemistryModule::MasterRunParallel(double dt) {
    if (pkg_to_send > 0) {
      // send packages to all free workers ...
      MasterSendPkgs(worker_list, work_pointer, sur_pointer, pkg_to_send,
-                     i_pkgs, free_workers, dt, iteration, control_iteration,
-                     wp_sizes_vector);
+                     i_pkgs, free_workers, dt, iteration, wp_sizes_vector);
    }
    // ... and try to receive them from workers who has finished their work
    MasterRecvPkgs(worker_list, pkg_to_recv, pkg_to_send > 0, free_workers);
@ -537,23 +512,37 @@ void poet::ChemistryModule::MasterRunParallel(double dt) {
  chem_field = out_vec;

  /* do master stuff */
+  if (!this->control_batch.empty()) {
+    std::cout << "[Master] Processing " << this->control_batch.size()
+              << " control cells for comparison." << std::endl;

-  if (control_iteration) {
-    control_iteration_counter++;
+    /* using mpi-buffer because we need cell-major layout*/
+    std::vector<std::vector<double>> surrogate_batch;
+    surrogate_batch.reserve(this->control_batch.size());

-    std::vector<double> sur_unshuffled{sur_shuffled};
+    for (const auto &element : this->control_batch) {

-    unshuffleField(sur_shuffled, this->n_cells, this->prop_count,
-                   wp_sizes_vector.size(), sur_unshuffled);
+      for (size_t i = 0; i < this->n_cells; i++) {
+        uint32_t curr_cell_id = mpi_buffer[this->prop_count * i];

-    error_stats stats(this->prop_count, control_iteration_counter *
-                                            runtime_params->control_iteration);
+        if (curr_cell_id == element[0]) {
+          std::vector<double> surrogate_output(
+              mpi_buffer.begin() + this->prop_count * i,
+              mpi_buffer.begin() + this->prop_count * (i + 1));
+          surrogate_batch.push_back(surrogate_output);
+          break;
+        }
+      }
+    }

-    computeStats(out_vec, sur_unshuffled, this->n_cells, this->prop_count,
-                 stats);
-    error_stats_history.push_back(stats);
+    metrics_a = MPI_Wtime();
+    control_module->computeSpeciesErrorMetrics(this->control_batch,
+                                               surrogate_batch);
+    metrics_b = MPI_Wtime();
+    this->metrics_t += metrics_b - metrics_a;

-    // to do: control values to epsilon
+    // Clear for next control iteration
+    this->control_batch.clear();
  }

  /* start time measurement of master chemistry */
--- a/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp
+++ b/src/Chemistry/SurrogateModels/DHT_Wrapper.cpp
@ -133,7 +133,7 @@ void DHT_Wrapper::fillDHT(const WorkPackage &work_package) {
      continue;
    }

-    if (work_package.input[i][0] != 2) {
+    if (work_package.input[i][1] != 2) {
      continue;
    }

--- a/src/Chemistry/SurrogateModels/InterpolationModule.cpp
+++ b/src/Chemistry/SurrogateModels/InterpolationModule.cpp
@ -76,7 +76,7 @@ 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) {
+    if (work_package.input[wp_i][1] != 2) {
      interp_result.status[wp_i] = INSUFFICIENT_DATA;
      continue;
    }
@ -122,7 +122,7 @@ void InterpolationModule::tryInterpolation(WorkPackage &work_package) {
    this->pht->incrementReadCounter(roundKey(rounded_key));
 #endif

-    const int cell_id = static_cast<int>(work_package.input[wp_i][0]);
+    const int cell_id = static_cast<int>(work_package.input[wp_i][1]);

    if (!to_calc_cache.contains(cell_id)) {
      const std::vector<std::int32_t> &to_calc = dht_instance.getKeyElements();
--- a/src/Chemistry/WorkerFunctions.cpp
+++ b/src/Chemistry/WorkerFunctions.cpp
--- a/src/Control/ControlModule.cpp
+++ b/src/Control/ControlModule.cpp
@ -0,0 +1,207 @@
+#include "ControlModule.hpp"
+#include "IO/Datatypes.hpp"
+#include "IO/HDF5Functions.hpp"
+#include "IO/StatsIO.hpp"
+#include <cmath>
+
+void poet::ControlModule::updateControlIteration(const uint32_t &iter,
+                                                 const bool &dht_enabled,
+                                                 const bool &interp_enabled) {
+
+  /* dht_enabled and inter_enabled are user settings set before startig the
+   * simulation*/
+  double prep_a, prep_b;
+
+  prep_a = MPI_Wtime();
+
+  /*
+  if (control_interval == 0) {
+    control_interval_enabled = false;
+    return;
+  }
+    */
+  global_iteration = iter;
+  initiateWarmupPhase(dht_enabled, interp_enabled);
+
+  /*
+  control_interval_enabled =
+      (control_interval > 0 && iter % control_interval == 0);
+
+
+  if (control_interval_enabled) {
+    MSG("[Control] Control interval enabled at iteration " +
+        std::to_string(iter));
+  }
+  */
+  prep_b = MPI_Wtime();
+  this->prep_t += prep_b - prep_a;
+}
+
+void poet::ControlModule::initiateWarmupPhase(bool dht_enabled,
+                                              bool interp_enabled) {
+
+  // user requested DHT/INTEP? keep them disabled but enable warmup-phase so
+  if (global_iteration < stabilization_interval || rollback_enabled) {
+    chem->SetWarmupEnabled(true);
+    chem->SetDhtEnabled(false);
+    chem->SetInterpEnabled(false);
+
+    MSG("Stabilization enabled until next control interval at iteration " +
+        std::to_string(stabilization_interval) + ".");
+
+    if (sur_disabled_counter > 0) {
+      --sur_disabled_counter;
+      MSG("Rollback counter: " + std::to_string(sur_disabled_counter));
+    } else {
+      rollback_enabled = false;
+    }
+
+    return;
+  }
+
+  chem->SetWarmupEnabled(false);
+  chem->SetDhtEnabled(dht_enabled);
+  chem->SetInterpEnabled(interp_enabled);
+}
+
+void poet::ControlModule::applyControlLogic(DiffusionModule &diffusion,
+                                            uint32_t &iter) {
+
+  /*
+  if (!control_interval_enabled) {
+    return;
+  }
+  */
+  writeCheckpointAndMetrics(diffusion, iter);
+
+  if (checkAndRollback(diffusion, iter)) {
+    rollback_enabled = true;
+    rollback_count++;
+    sur_disabled_counter = stabilization_interval;
+
+    MSG("Interpolation disabled for the next " +
+        std::to_string(stabilization_interval) + ".");
+  }
+}
+
+void poet::ControlModule::writeCheckpointAndMetrics(DiffusionModule &diffusion,
+                                                    uint32_t iter) {
+
+  double w_check_a, w_check_b, stats_a, stats_b;
+  MSG("Writing checkpoint of iteration " + std::to_string(iter));
+
+  w_check_a = MPI_Wtime();
+  write_checkpoint(out_dir, "checkpoint" + std::to_string(iter) + ".hdf5",
+                   {.field = diffusion.getField(), .iteration = iter});
+  w_check_b = MPI_Wtime();
+  this->w_check_t += w_check_b - w_check_a;
+
+  stats_a = MPI_Wtime();
+  writeStatsToCSV(metricsHistory, species_names, out_dir, "stats_overview");
+  stats_b = MPI_Wtime();
+
+  this->stats_t += stats_b - stats_a;
+}
+
+bool poet::ControlModule::checkAndRollback(DiffusionModule &diffusion,
+                                           uint32_t &iter) {
+  double r_check_a, r_check_b;
+
+  if (global_iteration < stabilization_interval) {
+    return false;
+  }
+
+  if (metricsHistory.empty()) {
+    MSG("No error history yet; skipping rollback check.");
+    return false;
+  }
+
+  const auto &mape = metricsHistory.back().mape;
+
+  for (size_t row = 0; row < mape.size(); row++) {
+    for (size_t col = 0; col < species_names.size() && col < mape[row].size(); col++) {
+      if (mape[row][col] == 0) {
+        continue;
+      }
+
+      if (mape[row][col] > mape_threshold[col]) {
+        uint32_t rollback_iter =
+            ((iter - 1) / checkpoint_interval) * checkpoint_interval;
+
+        MSG("[THRESHOLD EXCEEDED] " + species_names[col] +
+            " has MAPE = " + std::to_string(mape[row][col]) +
+            " exceeding threshold = " + std::to_string(mape_threshold[col]) +
+            ", rolling back to iteration " + std::to_string(rollback_iter));
+
+        r_check_a = MPI_Wtime();
+        Checkpoint_s checkpoint_read{.field = diffusion.getField()};
+        read_checkpoint(out_dir,
+                        "checkpoint" + std::to_string(rollback_iter) + ".hdf5",
+                        checkpoint_read);
+        iter = checkpoint_read.iteration;
+        r_check_b = MPI_Wtime();
+        r_check_t += r_check_b - r_check_a;
+        return true;
+      }
+    }
+  }
+  MSG("All species are within their MAPE thresholds.");
+
+  return false;
+}
+
+void poet::ControlModule::computeSpeciesErrorMetrics(
+    std::vector<std::vector<double>> &reference_values,
+    std::vector<std::vector<double>> &surrogate_values) {
+
+  const uint32_t num_cells = reference_values.size();
+  const uint32_t species_count = this->species_names.size();
+
+  std::cout << "[DEBUG] computeSpeciesErrorMetrics: num_cells=" << num_cells 
+            << ", species_count=" << species_count << std::endl;
+
+  SpeciesErrorMetrics metrics(num_cells, species_count, global_iteration,
+                              rollback_count);
+
+  if (reference_values.size() != surrogate_values.size()) {
+    MSG(" Reference and surrogate vectors differ in size: " +
+        std::to_string(reference_values.size()) + " vs " +
+        std::to_string(surrogate_values.size()));
+    return;
+  }
+
+  for (size_t cell_i = 0; cell_i < num_cells; cell_i++) {
+
+    metrics.id.push_back(reference_values[cell_i][0]);
+
+    for (size_t sp_i = 0; sp_i < reference_values[cell_i].size(); sp_i++) {
+      const double ref_value = reference_values[cell_i][sp_i];
+      const double sur_value = surrogate_values[cell_i][sp_i];
+      const double ZERO_ABS = 1e-13;
+
+      if (std::isnan(ref_value) || std::isnan(sur_value)) {
+        metrics.mape[cell_i][sp_i] = 0.0;
+        metrics.rrmse[cell_i][sp_i] = 0.0;
+        continue;
+      }
+
+      if (std::abs(ref_value) < ZERO_ABS) {
+        if (std::abs(sur_value) >= ZERO_ABS) {
+          metrics.mape[cell_i][sp_i] = 1.0;
+          metrics.rrmse[cell_i][sp_i] = 1.0;
+        } else {
+          metrics.mape[cell_i][sp_i] = 0.0;
+          metrics.rrmse[cell_i][sp_i] = 0.0;
+        }
+      } else {
+        double alpha = 1.0 - (sur_value / ref_value);
+        metrics.mape[cell_i][sp_i] = 100.0 * std::abs(alpha);
+        metrics.rrmse[cell_i][sp_i] = alpha * alpha;
+      }
+    }
+  }
+  
+  std::cout << "[DEBUG] metrics.id.size()=" << metrics.id.size() << std::endl;
+  metricsHistory.push_back(metrics);
+  std::cout << "[DEBUG] metricsHistory.size()=" << metricsHistory.size() << std::endl;
+}
--- a/src/Control/ControlModule.hpp
+++ b/src/Control/ControlModule.hpp
@ -0,0 +1,123 @@
+#ifndef CONTROLMODULE_H_
+#define CONTROLMODULE_H_
+
+#include "Base/Macros.hpp"
+#include "Chemistry/ChemistryModule.hpp"
+#include "Transport/DiffusionModule.hpp"
+#include "poet.hpp"
+
+#include <cstdint>
+#include <string>
+#include <vector>
+
+namespace poet {
+
+class ChemistryModule;
+class DiffusionModule;
+
+class ControlModule {
+
+public:
+  /* Control configuration*/
+
+  // std::uint32_t global_iter = 0;
+  // std::uint32_t sur_disabled_counter = 0;
+  // std::uint32_t rollback_counter = 0;
+
+  void updateControlIteration(const uint32_t &iter, const bool &dht_enabled,
+                              const bool &interp_enaled);
+
+  void initiateWarmupPhase(bool dht_enabled, bool interp_enabled);
+
+  bool checkAndRollback(DiffusionModule &diffusion, uint32_t &iter);
+
+  struct SpeciesErrorMetrics {
+    std::vector<std::uint32_t> id;
+    std::vector<std::vector<double>> mape;
+    std::vector<std::vector<double>> rrmse;
+    uint32_t iteration; // iterations in simulation after rollbacks
+    uint32_t rollback_count;
+
+    SpeciesErrorMetrics(uint32_t num_cells, uint32_t species_count,
+                        uint32_t iter, uint32_t counter)
+        : mape(num_cells, std::vector<double>(species_count, 0.0)),
+          rrmse(num_cells, std::vector<double>(species_count, 0.0)),
+          iteration(iter), rollback_count(counter) {}
+  };
+
+  void computeSpeciesErrorMetrics(
+      std::vector<std::vector<double>> &reference_values,
+      std::vector<std::vector<double>> &surrogate_values);
+
+  std::vector<SpeciesErrorMetrics> metricsHistory;
+
+  struct ControlSetup {
+    std::string out_dir;
+    std::uint32_t checkpoint_interval;
+    std::uint32_t penalty_interval;
+    std::uint32_t stabilization_interval;
+    std::vector<std::string> species_names;
+    std::vector<double> mape_threshold;
+    std::vector<uint32_t> ctrl_cell_ids;
+  };
+
+  void enableControlLogic(const ControlSetup &setup) {
+    this->out_dir = setup.out_dir;
+    this->checkpoint_interval = setup.checkpoint_interval;
+    this->stabilization_interval = setup.stabilization_interval;
+    this->species_names = setup.species_names;
+    this->mape_threshold = setup.mape_threshold;
+    this->ctrl_cell_ids = setup.ctrl_cell_ids;
+  }
+
+  void applyControlLogic(DiffusionModule &diffusion, uint32_t &iter);
+
+  void writeCheckpointAndMetrics(DiffusionModule &diffusion, uint32_t iter);
+
+  auto getGlobalIteration() const noexcept { return global_iteration; }
+
+  void setChemistryModule(poet::ChemistryModule *c) { chem = c; }
+
+  std::vector<double> getMapeThreshold() const { return this->mape_threshold; }
+
+  std::vector<uint32_t> getCtrlCellIds() const { return this->ctrl_cell_ids; }
+
+  /* Profiling getters */
+
+  auto getUpdateCtrlLogicTime() const { return this->prep_t; }
+
+  auto getWriteCheckpointTime() const { return this->w_check_t; }
+
+  auto getReadCheckpointTime() const { return this->r_check_t; }
+
+  auto getWriteMetricsTime() const { return this->stats_t; }
+
+private:
+  bool rollback_enabled = false;
+
+  poet::ChemistryModule *chem = nullptr;
+
+  std::uint32_t stabilization_interval = 0;
+  std::uint32_t penalty_interval = 0;
+  std::uint32_t checkpoint_interval = 0;
+  std::uint32_t global_iteration = 0;
+  std::uint32_t rollback_count = 0;
+  std::uint32_t sur_disabled_counter = 0;
+  std::vector<double> mape_threshold;
+  std::vector<uint32_t> ctrl_cell_ids;
+
+  std::vector<std::string> species_names;
+  std::string out_dir;
+
+  double prep_t = 0.;
+  double r_check_t = 0.;
+  double w_check_t = 0.;
+  double stats_t = 0.;
+
+  /* Buffer for shuffled surrogate data */
+  std::vector<double> sur_shuffled;
+};
+
+} // namespace poet
+
+#endif // CONTROLMODULE_H_
--- a/src/IO/HDF5Functions.hpp
+++ b/src/IO/HDF5Functions.hpp
@ -3,7 +3,7 @@
 #include <string>
 #include "Datatypes.hpp"

-int write_checkpoint(const std::string &file_path, struct Checkpoint_s &&checkpoint);

-int read_checkpoint(const std::string &file_path, struct Checkpoint_s &checkpoint);
+int write_checkpoint(const std::string &dir_path, const std::string &file_name, struct Checkpoint_s &&checkpoint);

+int read_checkpoint(const std::string &dir_path, const std::string &file_name, struct Checkpoint_s &checkpoint);
--- a/src/IO/StatsIO.cpp
+++ b/src/IO/StatsIO.cpp
@ -2,36 +2,61 @@
 #include <fstream>
 #include <iostream>
 #include <string>
+#include <iomanip> 
+#include <filesystem>

 namespace poet
 {
-    void writeStatsToCSV(const std::vector<ChemistryModule::error_stats> &all_stats,
+    void writeStatsToCSV(const std::vector<ControlModule::SpeciesErrorMetrics> &all_stats,
                         const std::vector<std::string> &species_names,
+                         const std::string &out_dir,
                         const std::string &filename)
    {
-        std::ofstream out(filename);
+        std::filesystem::path full_path = std::filesystem::path(out_dir) / filename;
+
+        std::ofstream out(full_path);
        if (!out.is_open())
        {
            std::cerr << "Could not open " << filename << " !" << std::endl;
            return;
        }

-        // header
-        out << "Iteration, Species, MAPE, RRSME \n";
+        // header: CellID, Iteration, Rollback, Species, MAPE, RRMSE
+        out << std::left << std::setw(15) << "CellID"
+            << std::setw(15) << "Iteration"
+            << std::setw(15) << "Rollback"
+            << std::setw(15) << "Species"
+            << std::setw(15) << "MAPE"
+            << std::setw(15) << "RRMSE" << "\n";

-        for (size_t i = 0; i < all_stats.size(); ++i)
+        out << std::string(90, '-') << "\n"; 
+
+        // data rows: iterate over iterations
+        for (size_t iter_idx = 0; iter_idx < all_stats.size(); ++iter_idx)
        {
-            for (size_t j = 0; j < species_names.size(); ++j)
+            const auto &metrics = all_stats[iter_idx];
+            
+            // Iterate over cells
+            for (size_t cell_idx = 0; cell_idx < metrics.id.size(); ++cell_idx)
            {
-                out << all_stats[i].iteration << ",\t"
-                    << species_names[j] << ",\t"
-                    << all_stats[i].mape[j] << ",\t"
-                    << all_stats[i].rrsme[j] << "\n";
+                // Iterate over species for this cell
+                for (size_t species_idx = 0; species_idx < species_names.size(); ++species_idx)
+                {
+                    out << std::left
+                        << std::setw(15) << metrics.id[cell_idx]
+                        << std::setw(15) << metrics.iteration
+                        << std::setw(15) << metrics.rollback_count
+                        << std::setw(15) << species_names[species_idx]
+                        << std::setw(15) << metrics.mape[cell_idx][species_idx]
+                        << std::setw(15) << metrics.rrmse[cell_idx][species_idx]
+                        << "\n";
+                }
            }
-            out << std::endl;
+            out << "\n"; 
        }

        out.close();
-        std::cout << "Stats written to " << filename << "\n";
+        std::cout << "Error metrics written to " << out_dir << "/" << filename << "\n";
    }
-} // namespace poet
+}
+ // namespace poet
--- a/src/IO/StatsIO.hpp
+++ b/src/IO/StatsIO.hpp
@ -1,9 +1,10 @@
 #include <string>
-#include "Chemistry/ChemistryModule.hpp"
+#include "Control/ControlModule.hpp"

 namespace poet
 {
-    void writeStatsToCSV(const std::vector<ChemistryModule::error_stats> &all_stats,
+    void writeStatsToCSV(const std::vector<ControlModule::SpeciesErrorMetrics> &all_stats,
                         const std::vector<std::string> &species_names,
+                         const std::string &out_dir,
                         const std::string &filename);
 } // namespace poet
--- a/src/IO/checkpoint.cpp
+++ b/src/IO/checkpoint.cpp
@ -1,13 +1,20 @@
 #include "IO/Datatypes.hpp"
 #include <cstdint>
 #include <highfive/H5Easy.hpp>
+#include <filesystem>

-int write_checkpoint(const std::string &file_path, struct Checkpoint_s &&checkpoint){
+namespace fs = std::filesystem;

+int write_checkpoint(const std::string &dir_path, const std::string &file_name, struct Checkpoint_s &&checkpoint){
+    
+    if (!fs::exists(dir_path)) {
+        std::cerr << "Directory does not exist: " << dir_path << std::endl;
+        return -1;
+    }
+    fs::path file_path = fs::path(dir_path) / file_name;
    // TODO: errorhandling
    H5Easy::File file(file_path, H5Easy::File::Overwrite);

-
    H5Easy::dump(file, "/MetaParam/Iterations", checkpoint.iteration);
    H5Easy::dump(file, "/Grid/Names", checkpoint.field.GetProps());
    H5Easy::dump(file, "/Grid/Chemistry", checkpoint.field.As2DVector());
@ -15,9 +22,16 @@ int write_checkpoint(const std::string &file_path, struct Checkpoint_s &&checkpo
    return 0;
 }

-int read_checkpoint(const std::string &file_path, struct Checkpoint_s &checkpoint){
+int read_checkpoint(const std::string &dir_path, const std::string &file_name, struct Checkpoint_s &checkpoint){
+    
+    fs::path file_path = fs::path(dir_path) / file_name;

-     H5Easy::File file(file_path, H5Easy::File::ReadOnly);
+    if (!fs::exists(file_path)) {
+        std::cerr << "File does not exist: " << file_path << std::endl;
+        return -1;
+    }
+     
+    H5Easy::File file(file_path, H5Easy::File::ReadOnly);

     checkpoint.iteration = H5Easy::load<uint32_t>(file, "/MetaParam/Iterations");

--- a/src/poet.cpp
+++ b/src/poet.cpp
@ -25,10 +25,8 @@
 #include "Base/RInsidePOET.hpp"
 #include "CLI/CLI.hpp"
 #include "Chemistry/ChemistryModule.hpp"
+#include "Control/ControlModule.hpp"
 #include "DataStructures/Field.hpp"
-#include "IO/Datatypes.hpp"
-#include "IO/HDF5Functions.hpp"
-#include "IO/StatsIO.hpp"
 #include "Init/InitialList.hpp"
 #include "Transport/DiffusionModule.hpp"

@ -68,8 +66,7 @@ static poet::DEFunc ReadRObj_R;
 static poet::DEFunc SaveRObj_R;
 static poet::DEFunc source_R;

-static void init_global_functions(RInside &R)
-{
+static void init_global_functions(RInside &R) {
  R.parseEval(kin_r_library);
  master_init_R = DEFunc("master_init");
  master_iteration_end_R = DEFunc("master_iteration_end");
@ -92,15 +89,9 @@ static void init_global_functions(RInside &R)
 //   R.parseEval("mysetup$state_C <- TMP");
 // }

-enum ParseRet
-{
-  PARSER_OK,
-  PARSER_ERROR,
-  PARSER_HELP
-};
+enum ParseRet { PARSER_OK, PARSER_ERROR, PARSER_HELP };

-int parseInitValues(int argc, char **argv, RuntimeParameters &params)
-{
+int parseInitValues(int argc, char **argv, RuntimeParameters &params) {

  CLI::App app{"POET - Potsdam rEactive Transport simulator"};

@ -182,12 +173,9 @@ int parseInitValues(int argc, char **argv, RuntimeParameters &params)
                 "Output directory of the simulation")
      ->required();

-  try
-  {
+  try {
    app.parse(argc, argv);
-  }
-  catch (const CLI::ParseError &e)
-  {
+  } catch (const CLI::ParseError &e) {
    app.exit(e);
    return -1;
  }
@ -199,16 +187,14 @@ int parseInitValues(int argc, char **argv, RuntimeParameters &params)
  if (params.as_qs)
    params.out_ext = "qs";

-  if (MY_RANK == 0)
-  {
+  if (MY_RANK == 0) {
    // MSG("Complete results storage is " + BOOL_PRINT(simparams.store_result));
    MSG("Output format/extension is " + params.out_ext);
    MSG("Work Package Size: " + std::to_string(params.work_package_size));
    MSG("DHT is " + BOOL_PRINT(params.use_dht));
    MSG("AI Surrogate is " + BOOL_PRINT(params.use_ai_surrogate));

-    if (params.use_dht)
-    {
+    if (params.use_dht) {
      // MSG("DHT strategy is " + std::to_string(simparams.dht_strategy));
      // MDL: these should be outdated (?)
      // MSG("DHT key default digits (ignored if 'signif_vector' is "
@ -222,8 +208,7 @@ int parseInitValues(int argc, char **argv, RuntimeParameters &params)
      // MSG("DHT load file is " + chem_params.dht_file);
    }

-    if (params.use_interp)
-    {
+    if (params.use_interp) {
      MSG("PHT interpolation enabled: " + BOOL_PRINT(params.use_interp));
      MSG("PHT interp-size = " + std::to_string(params.interp_size));
      MSG("PHT interp-min  = " + std::to_string(params.interp_min_entries));
@ -251,8 +236,7 @@ int parseInitValues(int argc, char **argv, RuntimeParameters &params)
  // // log before rounding?
  // R["dht_log"] = simparams.dht_log;

-  try
-  {
+  try {

    Rcpp::List init_params_(ReadRObj_R(init_file));
    params.init_params = init_params_;
@ -266,17 +250,18 @@ int parseInitValues(int argc, char **argv, RuntimeParameters &params)

    params.timesteps =
        Rcpp::as<std::vector<double>>(global_rt_setup->operator[]("timesteps"));
-    params.control_iteration =
-        Rcpp::as<uint32_t>(global_rt_setup->operator[]("control_iteration"));
-    params.species_epsilon =
-        Rcpp::as<std::vector<double>>(global_rt_setup->operator[]("species_epsilon"));
-    params.penalty_iteration =
-        Rcpp::as<uint32_t>(global_rt_setup->operator[]("penalty_iteration"));
-    params.max_penalty_iteration = 
-        Rcpp::as<uint32_t>(global_rt_setup->operator[]("max_penalty_iteration"));
-  }
-  catch (const std::exception &e)
-  {
+    params.checkpoint_interval =
+        Rcpp::as<uint32_t>(global_rt_setup->operator[]("checkpoint_interval"));
+    params.stabilization_interval = 
+        Rcpp::as<uint32_t>(global_rt_setup->operator[]("stabilization_interval"));
+    params.penalty_interval =
+        Rcpp::as<uint32_t>(global_rt_setup->operator[]("penalty_interval"));
+    params.mape_threshold = Rcpp::as<std::vector<double>>(
+        global_rt_setup->operator[]("mape_threshold"));
+    params.ctrl_cell_ids = Rcpp::as<std::vector<uint32_t>>(
+        global_rt_setup->operator[]("ctrl_cell_ids"));
+
+  } catch (const std::exception &e) {
    ERRMSG("Error while parsing R scripts: " + std::string(e.what()));
    return ParseRet::PARSER_ERROR;
  }
@ -286,8 +271,7 @@ int parseInitValues(int argc, char **argv, RuntimeParameters &params)

 // HACK: this is a step back as the order and also the count of fields is
 // predefined, but it will change in the future
-void call_master_iter_end(RInside &R, const Field &trans, const Field &chem)
-{
+void call_master_iter_end(RInside &R, const Field &trans, const Field &chem) {
  R["TMP"] = Rcpp::wrap(trans.AsVector());
  R["TMP_PROPS"] = Rcpp::wrap(trans.GetProps());
  R.parseEval(std::string("state_T <- setNames(data.frame(matrix(TMP, nrow=" +
@ -304,84 +288,24 @@ void call_master_iter_end(RInside &R, const Field &trans, const Field &chem)
  *global_rt_setup = R["setup"];
 }

-bool checkAndRollback(ChemistryModule &chem, RuntimeParameters &params, uint32_t &iter)
-{
-  const std::vector<double> &latest_mape = chem.error_stats_history.back().mape;
-
-  for (uint32_t j = 0; j < params.species_epsilon.size(); j++)
-  {
-    if (params.species_epsilon[j] < latest_mape[j] && latest_mape[j] != 0)
-    {
-      uint32_t rollback_iter = iter - (iter % params.control_iteration);
-
-      std::cout << chem.getField().GetProps()[j] << " with a MAPE value of " << latest_mape[j] << " exceeds epsilon of "
-                << params.species_epsilon[j] << "! " << std::endl;
-
-      Checkpoint_s checkpoint_read{.field = chem.getField()};
-      read_checkpoint("checkpoint" + std::to_string(rollback_iter) + ".hdf5", checkpoint_read);
-      iter = checkpoint_read.iteration;
-
-      return true;
-    } 
-  }
-  MSG("All spezies are below their threshold values");
-  return false;
-}
-
-void updatePenaltyLogic(RuntimeParameters &params, bool roolback_happend)
-{
-  if (roolback_happend)
-  {
-    params.rollback_simulation = true;
-    params.penalty_counter = params.penalty_iteration;
-    std::cout << "Penalty counter reset to: " << params.penalty_counter << std::endl;
-    MSG("Rollback! Penalty phase started for " + std::to_string(params.penalty_iteration) + " iterations.");
-  }
-  else
-  {
-    if (params.rollback_simulation && params.penalty_counter == 0)
-    {
-      params.rollback_simulation = false;
-      MSG("Penalty phase ended. Interpolation re-enabled.");
-    }
-    else if (!params.rollback_simulation)
-    {
-      params.penalty_iteration = std::min(params.penalty_iteration *= 2, params.max_penalty_iteration);
-      MSG("Stable surrogate phase detected. Penalty iteration doubled to " + std::to_string(params.penalty_iteration) + " iterations.");
-    }
-  }
-}
-
 static Rcpp::List RunMasterLoop(RInsidePOET &R, RuntimeParameters &params,
                                DiffusionModule &diffusion,
-                                ChemistryModule &chem)
-{
+                                ChemistryModule &chem, ControlModule &control) {

  /* Iteration Count is dynamic, retrieving value from R (is only needed by
   * master for the following loop) */
  uint32_t maxiter = params.timesteps.size();

-  if (params.print_progress)
-  {
+  if (params.print_progress) {
    chem.setProgressBarPrintout(true);
  }
  R["TMP_PROPS"] = Rcpp::wrap(chem.getField().GetProps());

-  params.next_penalty_check = params.penalty_iteration;
-
  /* SIMULATION LOOP */

  double dSimTime{0};
-  for (uint32_t iter = 1; iter < maxiter + 1; iter++)
-  {
-    // Penalty countdown
-    if (params.rollback_simulation && params.penalty_counter > 0)
-    {
-      params.penalty_counter--;
-      std::cout << "Penalty counter: " << params.penalty_counter << std::endl;
-    }
-
-    params.control_iteration_active = (iter % params.control_iteration == 0 /* && iter != 0 */);
+  for (uint32_t iter = 1; iter < maxiter + 1; iter++) {
+    control.updateControlIteration(iter, params.use_dht, params.use_interp);

    double start_t = MPI_Wtime();

@ -398,13 +322,10 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, RuntimeParameters &params,
    /* run transport */
    diffusion.simulate(dt);

-    chem.runtime_params = &params;
-
    chem.getField().update(diffusion.getField());

    // MSG("Chemistry start");
-    if (params.use_ai_surrogate)
-    {
+    if (params.use_ai_surrogate) {
      double ai_start_t = MPI_Wtime();
      // Save current values from the tug field as predictor for the ai step
      R["TMP"] = Rcpp::wrap(chem.getField().AsVector());
@ -455,8 +376,7 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, RuntimeParameters &params,
    chem.simulate(dt);

    /* AI surrogate iterative training*/
-    if (params.use_ai_surrogate)
-    {
+    if (params.use_ai_surrogate) {
      double ai_start_t = MPI_Wtime();

      R["TMP"] = Rcpp::wrap(chem.getField().AsVector());
@ -495,21 +415,7 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, RuntimeParameters &params,
    MSG("End of *coupling* iteration " + std::to_string(iter) + "/" +
        std::to_string(maxiter));

-    if (iter % params.control_iteration == 0)
-    {
-      writeStatsToCSV(chem.error_stats_history, chem.getField().GetProps(), "stats_overview");
-      write_checkpoint("checkpoint" + std::to_string(iter) + ".hdf5",
-                       {.field = chem.getField(), .iteration = iter});
-    }
-
-    if (iter == params.next_penalty_check)
-    {
-      bool roolback_happend = checkAndRollback(chem, params, iter);
-      updatePenaltyLogic(params, roolback_happend);
-
-      params.next_penalty_check = iter + params.penalty_iteration;
-    }
-
+    control.applyControlLogic(diffusion, iter);
    // MSG();
  } // END SIMULATION LOOP

@ -526,8 +432,7 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, RuntimeParameters &params,
  Rcpp::List diffusion_profiling;
  diffusion_profiling["simtime"] = diffusion.getTransportTime();

-  if (params.use_dht)
-  {
+  if (params.use_dht) {
    chem_profiling["dht_hits"] = Rcpp::wrap(chem.GetWorkerDHTHits());
    chem_profiling["dht_evictions"] = Rcpp::wrap(chem.GetWorkerDHTEvictions());
    chem_profiling["dht_get_time"] = Rcpp::wrap(chem.GetWorkerDHTGetTimings());
@ -535,8 +440,7 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, RuntimeParameters &params,
        Rcpp::wrap(chem.GetWorkerDHTFillTimings());
  }

-  if (params.use_interp)
-  {
+  if (params.use_interp) {
    chem_profiling["interp_w"] =
        Rcpp::wrap(chem.GetWorkerInterpolationWriteTimings());
    chem_profiling["interp_r"] =
@ -560,9 +464,31 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, RuntimeParameters &params,
  return profiling;
 }

+static void getControlCellIds(std::vector<std::uint32_t> &ids, int root,
+                              MPI_Comm comm) {
+  std::uint32_t n_ids = 0;
+  int rank;
+  MPI_Comm_rank(comm, &rank);
+  bool is_master = root == rank;
+
+  if (is_master) {
+    n_ids = ids.size();
+  }
+  // broadcast size of id vector
+  MPI_Bcast(&n_ids, 1, MPI_UINT32_T, root, comm);
+
+  // worker
+  if (!is_master) {
+    ids.resize(n_ids);
+  }
+  // broadcast control cell ids
+  if (n_ids > 0) {
+    MPI_Bcast(ids.data(), n_ids, MPI_UINT32_T, root, comm);
+  }
+}
+
 std::vector<std::string> getSpeciesNames(const Field &&field, int root,
-                                         MPI_Comm comm)
-{
+                                         MPI_Comm comm) {
  std::uint32_t n_elements;
  std::uint32_t n_string_size;

@ -572,13 +498,11 @@ std::vector<std::string> getSpeciesNames(const Field &&field, int root,
  const bool is_master = root == rank;

  // first, the master sends all the species names iterative
-  if (is_master)
-  {
+  if (is_master) {
    n_elements = field.GetProps().size();
    MPI_Bcast(&n_elements, 1, MPI_UINT32_T, root, MPI_COMM_WORLD);

-    for (std::uint32_t i = 0; i < n_elements; i++)
-    {
+    for (std::uint32_t i = 0; i < n_elements; i++) {
      n_string_size = field.GetProps()[i].size();
      MPI_Bcast(&n_string_size, 1, MPI_UINT32_T, root, MPI_COMM_WORLD);
      MPI_Bcast(const_cast<char *>(field.GetProps()[i].c_str()), n_string_size,
@ -593,8 +517,7 @@ std::vector<std::string> getSpeciesNames(const Field &&field, int root,

  std::vector<std::string> species_names_out(n_elements);

-  for (std::uint32_t i = 0; i < n_elements; i++)
-  {
+  for (std::uint32_t i = 0; i < n_elements; i++) {
    MPI_Bcast(&n_string_size, 1, MPI_UINT32_T, root, MPI_COMM_WORLD);

    char recv_buf[n_string_size];
@ -607,8 +530,7 @@ std::vector<std::string> getSpeciesNames(const Field &&field, int root,
  return species_names_out;
 }

-std::array<double, 2> getBaseTotals(Field &&field, int root, MPI_Comm comm)
-{
+std::array<double, 2> getBaseTotals(Field &&field, int root, MPI_Comm comm) {
  std::array<double, 2> base_totals;

  int rank;
@ -616,8 +538,7 @@ std::array<double, 2> getBaseTotals(Field &&field, int root, MPI_Comm comm)

  const bool is_master = root == rank;

-  if (is_master)
-  {
+  if (is_master) {
    const auto h_col = field["H"];
    const auto o_col = field["O"];

@ -632,8 +553,7 @@ std::array<double, 2> getBaseTotals(Field &&field, int root, MPI_Comm comm)
  return base_totals;
 }

-bool getHasID(Field &&field, int root, MPI_Comm comm)
-{
+bool getHasID(Field &&field, int root, MPI_Comm comm) {
  bool has_id;

  int rank;
@ -641,8 +561,7 @@ bool getHasID(Field &&field, int root, MPI_Comm comm)

  const bool is_master = root == rank;

-  if (is_master)
-  {
+  if (is_master) {
    const auto ID_field = field["ID"];

    std::set<double> unique_IDs(ID_field.begin(), ID_field.end());
@ -659,8 +578,7 @@ bool getHasID(Field &&field, int root, MPI_Comm comm)
  return has_id;
 }

-int main(int argc, char *argv[])
-{
+int main(int argc, char *argv[]) {
  int world_size;

  MPI_Init(&argc, &argv);
@ -671,8 +589,7 @@ int main(int argc, char *argv[])

    RInsidePOET &R = RInsidePOET::getInstance();

-    if (MY_RANK == 0)
-    {
+    if (MY_RANK == 0) {
      MSG("Running POET version " + std::string(poet_version));
    }

@ -680,8 +597,7 @@ int main(int argc, char *argv[])

    RuntimeParameters run_params;

-    if (parseInitValues(argc, argv, run_params) != 0)
-    {
+    if (parseInitValues(argc, argv, run_params) != 0) {
      MPI_Finalize();
      return 0;
    }
@ -706,6 +622,9 @@ int main(int argc, char *argv[])

    ChemistryModule chemistry(run_params.work_package_size,
                              init_list.getChemistryInit(), MPI_COMM_WORLD);
+    ControlModule control;
+    chemistry.SetControlModule(&control);
+    control.setChemistryModule(&chemistry);

    const ChemistryModule::SurrogateSetup surr_setup = {
        getSpeciesNames(init_list.getInitialGrid(), 0, MPI_COMM_WORLD),
@ -723,12 +642,23 @@ int main(int argc, char *argv[])

    chemistry.masterEnableSurrogates(surr_setup);

-    if (MY_RANK > 0)
-    {
+    /* broadcast control cell ids before simulation starts */
+    getControlCellIds(run_params.ctrl_cell_ids, 0, MPI_COMM_WORLD);
+    chemistry.SetControlCellIds(run_params.ctrl_cell_ids);
+
+    const ControlModule::ControlSetup ctrl_setup = {
+        run_params.out_dir, // added
+        run_params.checkpoint_interval,
+        run_params.penalty_interval,
+        run_params.stabilization_interval,
+        getSpeciesNames(init_list.getInitialGrid(), 0, MPI_COMM_WORLD),
+        run_params.mape_threshold};
+
+    control.enableControlLogic(ctrl_setup);
+
+    if (MY_RANK > 0) {
      chemistry.WorkerLoop();
-    }
-    else
-    {
+    } else {
      // R.parseEvalQ("mysetup <- setup");
      // // if (MY_RANK == 0) { // get timestep vector from
      // // grid_init function ... //
@ -742,8 +672,7 @@ int main(int argc, char *argv[])
      R["out_ext"] = run_params.out_ext;
      R["out_dir"] = run_params.out_dir;

-      if (run_params.use_ai_surrogate)
-      {
+      if (run_params.use_ai_surrogate) {
        /* Incorporate ai surrogate from R */
        R.parseEvalQ(ai_surrogate_r_library);
        /* Use dht species for model input and output */
@ -770,7 +699,8 @@ int main(int argc, char *argv[])

      chemistry.masterSetField(init_list.getInitialGrid());

-      Rcpp::List profiling = RunMasterLoop(R, run_params, diffusion, chemistry);
+      Rcpp::List profiling =
+          RunMasterLoop(R, run_params, diffusion, chemistry, control);

      MSG("finished simulation loop");

@ -792,8 +722,7 @@ int main(int argc, char *argv[])

  MPI_Finalize();

-  if (MY_RANK == 0)
-  {
+  if (MY_RANK == 0) {
    MSG("done, bye!");
  }

--- a/src/poet.hpp.in
+++ b/src/poet.hpp.in
@ -41,7 +41,6 @@ static const inline std::string r_runtime_parameters = "mysetup";
 struct RuntimeParameters {
  std::string out_dir;
  std::vector<double> timesteps;
-  std::vector<double> species_epsilon;

  Rcpp::List init_params;

@ -51,14 +50,12 @@ struct RuntimeParameters {
  std::string out_ext;

  bool print_progress = false;
-
-  std::uint32_t penalty_iteration = 0;
-  std::uint32_t max_penalty_iteration = 0;
-  std::uint32_t penalty_counter = 0;
-  std::uint32_t next_penalty_check = 0;
-  bool rollback_simulation = false;
-  bool control_iteration_active = false;
-  std::uint32_t control_iteration = 1;
+  std::uint32_t penalty_interval = 0;
+  std::uint32_t stabilization_interval = 0;
+  std::uint32_t checkpoint_interval = 0;
+  std::uint32_t control_interval = 0;
+  std::vector<double> mape_threshold;
+  std::vector<uint32t_t> ctrl_cell_ids;

  static constexpr std::uint32_t WORK_PACKAGE_SIZE_DEFAULT = 32;
  std::uint32_t work_package_size = WORK_PACKAGE_SIZE_DEFAULT;
Author	SHA1	Message	Date
rastogi	f05e65779e	Stop tracking bin/dolo/ and update .gitignore	2025-11-13 16:00:47 +01:00
rastogi	93c4463784	updated .gitignore	2025-11-12 10:25:58 +01:00
rastogi	5a3fe3a888	Store per-cell MAPE/RRMSE and add stabilization interval	2025-11-12 10:20:28 +01:00
rastogi	f15f9049b8	Fix cell_ID access in DHT/interpolation and restructure control cell metric computation to use row-major layout	2025-11-07 15:55:26 +01:00
rastogi	23d0cc2dd8	Fix cell_ID access in DHT/interpolation and restructure control cell metric computation to use row-major layout	2025-11-07 15:52:39 +01:00
rastogi	1f70dc4070	feat: Add control logic for surrogate validation with selective PHREEQC re-execution, error metrics (MAPE/RRMSE), and automatic checkpoint/rollback	2025-11-06 20:55:02 +01:00
rastogi	d0faa8e6c3	feat: bcast control_cell_ids to workers	2025-11-02 22:26:14 +01:00
rastogi	9810221182	Added ctrl_cell_ids	2025-10-30 21:36:03 +01:00
rastogi	ed4fb3b1b9	Restructuring prototyp 3	2025-10-30 16:15:27 +01:00
rastogi	aaf7dc7692	Add .gitignore	2025-10-30 16:14:22 +01:00