max/poet
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/poet.git synced 2025-12-16 12:54:50 +01:00
Code Issues Packages Projects Releases Wiki Activity

Control component with minimum features

Browse Source
This commit is contained in:
rastogi 2025-10-15 10:15:21 +02:00
parent 1a1054faa4
commit 96567f9074
39 changed files with 941 additions and 895 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

43
bin/dol.pqi Normal file
View File

@ -0,0 +1,43 @@
SOLUTION 1
units mol/kgw
water 1
temperature 25
pH 7
pe 4
PURE 1
Calcite 0.0 1
END
RUN_CELLS
-cells 1
COPY solution 1 2
PURE 2
O2g -0.1675 10
KINETICS 2
Calcite
-m 0.000207
-parms 0.05
-tol 1e-10
Dolomite
-m 0.0
-parms 0.005
-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

BIN
bin/dolo/iter_00.qs2
View File

Binary file not shown.

BIN
bin/dolo/iter_000.qs2
View File

Binary file not shown.

BIN
bin/dolo/iter_050.qs2
View File

Binary file not shown.

BIN
bin/dolo/iter_100.qs2
View File

Binary file not shown.

BIN
bin/dolo/iter_150.qs2
View File

Binary file not shown.

BIN
bin/dolo/iter_200.qs2
View File

Binary file not shown.

BIN
bin/dolo/iter_250.qs2
View File

Binary file not shown.

BIN
bin/dolo/iter_300.qs2
View File

Binary file not shown.

BIN
bin/dolo/iter_350.qs2
View File

Binary file not shown.

BIN
bin/dolo/iter_400.qs2
View File

Binary file not shown.

BIN
bin/dolo/iter_450.qs2
View File

Binary file not shown.

BIN
bin/dolo/iter_500.qs2
View File

Binary file not shown.

BIN
bin/dolo/timings.qs2
View File

Binary file not shown.

BIN
bin/dolo_inner.rds Normal file
View File

Binary file not shown.

131
bin/dolo_interp.R Normal file
View File

@ -0,0 +1,131 @@
rows <- 400
cols <- 200
grid_def <- matrix(2, nrow = rows, ncol = cols)
# Define grid configuration for POET model
grid_setup <- list(
pqc_in_file = "./dol.pqi",
pqc_db_file = "./phreeqc_kin.dat", # Path to the database file for Phreeqc
grid_def = grid_def, # Definition of the grid, containing IDs according to the Phreeqc input script
grid_size = c(2.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)
}
# fuzz_input_dht_keys <- function(input) {
# dht_species <- c(
# "H" = 3,
# "O" = 3,
# "Charge" = 3,
# "C" = 6,
# "Ca" = 6,
# "Cl" = 3,
# "Mg" = 5,
# "Calcite" = 4,
# "Dolomite" = 4
# )
# return(input[names(dht_species)])
# }
check_sign_cal_dol_interp <- function(to_interp, data_set) {
dht_species <- c(
"H" = 3,
"O" = 3,
"Charge" = 3,
"C" = 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" = 2,
"Calcite" = 2,
"Dolomite" = 2
)
chemistry_setup <- list(
dht_species = c(
"H" = 3,
"O" = 3,
"Charge" = 3,
"C" = 6,
"Ca" = 6,
"Cl" = 3,
"Mg" = 5,
"Calcite" = 4,
"Dolomite" = 4
),
pht_species = pht_species,
hooks = list(
dht_fill = check_sign_cal_dol_dht,
# dht_fuzz = fuzz_input_dht_keys,
interp_pre = check_sign_cal_dol_interp,
interp_post = check_neg_cal_dol
)
)
# Define a setup list for simulation configuration
setup <- list(
Grid = grid_setup, # Parameters related to the grid structure
Diffusion = diffusion_setup, # Parameters related to the diffusion process
Chemistry = chemistry_setup # Parameters related to the chemistry process
)

BIN
bin/dolo_interp.qs2
View File

Binary file not shown.

BIN
bin/dolo_interp/iter_00.qs2 Normal file
View File

Binary file not shown.

BIN
bin/dolo_interp/iter_50.qs2 Normal file
View File

Binary file not shown.

BIN
bin/dolo_interp/timings.qs2 Normal file
View File

Binary file not shown.

16
bin/dolo_interp_rt.R
View File

@ -1,17 +1,15 @@
iterations <- 500
iterations <- 50
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)
checkpoint_interval <- 30
control_interval <- 40
mape_threshold <- 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
checkpoint_interval = checkpoint_interval,
control_interval = control_interval,
mape_threshold = mape_threshold
)

79
bin/error.inp Normal file
View File

@ -0,0 +1,79 @@
SOLUTION_RAW 1 Solution after simulation 2.
-temp 25
-pressure 1
-potential 0
-total_h 2
-total_o 111.01243359386
-cb 55.506585818252
-density 1
-viscosity 1
-viscos_0 0.89002391829221
-totals
Ca 0.00012300690436983
Cl 0.00012300690436997
-pH 7
-pe 4
-mu 1e-07
-ah2o 1
-mass_water 1.9984236470542
-soln_vol 1
-total_alkalinity 0
-activities
-gammas
KINETICS_RAW 1
# KINETICS_MODIFY candidate identifiers #
-step_divide 1
-rk 3
-bad_step_max 500
-use_cvode 0
-cvode_steps 100
-cvode_order 5
-component Calcite
# KINETICS_MODIFY candidate identifiers #
-tol 1e-10
-m 1.4255219292182e-208
-m0 0.000207
-namecoef
Calcite 1
-d_params
0.000207
# KineticsComp workspace variables #
-moles 0
-initial_moles 0
-component Dolomite
# KINETICS_MODIFY candidate identifiers #
-tol 1e-10
-m 0.05
-m0 0
-namecoef
Dolomite 1
-d_params
0
# KineticsComp workspace variables #
-moles 0
-initial_moles 0
-equal_increments 0
-count 0
-steps
1
# KINETICS workspace variables #
-totals
EQUILIBRIUM_PHASES_RAW 1
# EXCHANGE_MODIFY candidates; use new_def=true #
-new_def 0
-component O2g
# EQUILIBRIUM_PHASES_MODIFY candidate identifiers #
-si 10
-moles 0.005
-force_equality 0
-dissolve_only 0
-precipitate_only 0
# PPassemblage workspace variables #
-si_org -0.1675
-delta 0
-initial_moles 0
-totals
-eltList # List of all elements in phases and alternate reactions
O 2
# PPassemblage workspace variables #
-assemblage_totals

BIN
bin/poet
View File

Binary file not shown.

BIN
bin/poet_init
View File

Binary file not shown.

392
bin/stats_overview
View File

@ -1,376 +1,18 @@
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
Iteration Rollback Species MAPE RRSME
---------------------------------------------------------------------------
40 0 cell_ID 0 0
40 0 ID 0 0
40 0 H 1.98104e-14 3.20223e-16
40 0 O 0.108665 0.00113056
40 0 Charge 0.0118752 0.000124408
40 0 C 0.000321187 3.76695e-05
40 0 Ca 0.000321235 3.76683e-05
40 0 Cl 2.96333 0.0577691
40 0 Mg 3.01563 0.0600123
40 0 Calcite 0.000193346 2.27729e-05
40 0 Calcite_p1 0 0
40 0 Dolomite 3.52007e-06 5.97467e-06
40 0 Dolomite_p1 0 0
40 0 O2g_eq 0.292809 0.00304717
40 0 O2g_si 0 0

BIN
share/poet/barite/barite_200.qs2
View File

Binary file not shown.

BIN
share/poet/barite/barite_het.qs2
View File

Binary file not shown.

BIN
share/poet/dolo/dolo_inner_large.qs2
View File

Binary file not shown.

BIN
share/poet/dolo/dolo_interp.qs2
View File

Binary file not shown.

BIN
share/poet/surfex/PoetEGU_surfex_500.qs2
View File

Binary file not shown.

46
src/Chemistry/ChemistryModule.hpp
View File

@ -185,6 +185,13 @@ namespace poet
*/
auto GetMasterLoopTime() const { return this->send_recv_t; }
auto GetMasterCtrlLogicTime() const { return this->ctrl_t; }
auto GetMasterCtrlBcastTime() const { return this->bcast_ctrl_t; }
auto GetMasterRecvCtrlLogicTime() const { return this->recv_ctrl_t; }
/**
* **Master only** Collect and return all accumulated timings recorded by
* workers to run Phreeqc simulation.
@ -214,6 +221,8 @@ namespace poet
*/
std::vector<double> GetWorkerIdleTimings() const;
std::vector<double> GetWorkerControlTimings() const;
/**
* **Master only** Collect and return DHT hits of all workers.
*
@ -262,25 +271,29 @@ namespace poet
std::vector<int> ai_surrogate_validity_vector;
RuntimeParameters *runtime_params = nullptr;
uint32_t control_iteration_counter = 0;
struct error_stats
struct SimulationErrorStats
{
std::vector<double> mape;
std::vector<double> rrsme;
uint32_t iteration;
std::vector<double> rrmse;
uint32_t iteration; // iterations in simulation after rollbacks
uint32_t rollback_count;
error_stats(size_t species_count, size_t iter)
: mape(species_count, 0.0), rrsme(species_count, 0.0), iteration(iter) {}
SimulationErrorStats(size_t species_count, uint32_t iter, uint32_t counter)
: mape(species_count, 0.0),
rrmse(species_count, 0.0),
iteration(iter),
rollback_count(counter){}
};
std::vector<error_stats> error_stats_history;
std::vector<SimulationErrorStats> error_history;
static void computeSpeciesErrors(const std::vector<double> &reference_values,
const std::vector<double> &surrogate_values,
uint32_t size_per_prop,
uint32_t species_count,
SimulationErrorStats &species_error_stats);
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,
@ -319,6 +332,7 @@ namespace poet
enum
{
WORKER_PHREEQC,
WORKER_CTRL_ITER,
WORKER_DHT_GET,
WORKER_DHT_FILL,
WORKER_IDLE,
@ -342,6 +356,7 @@ namespace poet
double dht_get = 0.;
double dht_fill = 0.;
double idle_t = 0.;
double ctrl_t = 0.;
};
struct worker_info_s
@ -410,6 +425,7 @@ namespace poet
poet::DHT_Wrapper *dht = nullptr;
bool dht_fill_during_rollback{false};
bool interp_enabled{false};
std::unique_ptr<poet::InterpolationModule> interp;
@ -431,6 +447,10 @@ namespace poet
double seq_t = 0.;
double send_recv_t = 0.;
double ctrl_t = 0.;
double bcast_ctrl_t = 0.;
double recv_ctrl_t = 0.;
std::array<double, 2> base_totals{0};
bool print_progessbar{false};
@ -449,7 +469,7 @@ namespace poet
std::unique_ptr<PhreeqcRunner> pqc_runner;
std::vector<double> sur_shuffled;
std::vector<double> sur_shuffled;
};
} // namespace poet

128
src/Chemistry/MasterFunctions.cpp
View File

@ -41,6 +41,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,35 +166,35 @@ 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,
void poet::ChemistryModule::computeSpeciesErrors(const std::vector<double> &reference_values,
const std::vector<double> &surrogate_values,
uint32_t size_per_prop,
uint32_t species_count,
error_stats &stats) {
SimulationErrorStats &species_error_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];
const double ref_value = reference_values[base_idx + j];
const double sur_value = surrogate_values[base_idx + j];
if (pqc_value == 0.0) {
if (ref_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);
double alpha = 1.0 - (sur_value / ref_value);
err_sum += std::abs(alpha);
sqr_err_sum += alpha * alpha;
}
}
stats.mape[i] = 100.0 * (err_sum / size_per_prop);
stats.rrsme[i] =
species_error_stats.mape[i] = 100.0 * (err_sum / size_per_prop);
species_error_stats.rrmse[i] =
(size_per_prop > 0) ? std::sqrt(sqr_err_sum / size_per_prop) : 0.0;
}
}
@ -264,7 +270,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) {
uint32_t control_interval, const std::vector<uint32_t> &wp_sizes_vector) {
/* declare variables */
int local_work_package_size;
@ -305,7 +311,7 @@ inline void poet::ChemistryModule::MasterSendPkgs(
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;
send_buffer[end_of_wp + 5] = control_interval;
/* ATTENTION Worker p has rank p+1 */
// MPI_Send(send_buffer, end_of_wp + BUFFER_OFFSET, MPI_DOUBLE, p + 1,
@ -329,6 +335,7 @@ inline void poet::ChemistryModule::MasterRecvPkgs(worker_list_t &w_list,
int need_to_receive = 1;
double idle_a, idle_b;
int p, size;
double recv_a, recv_b;
MPI_Status probe_status;
// master_recv_a = MPI_Wtime();
@ -361,6 +368,7 @@ inline void poet::ChemistryModule::MasterRecvPkgs(worker_list_t &w_list,
free_workers++;
}
if (probe_status.MPI_TAG == LOOP_CTRL) {
recv_a = MPI_Wtime();
MPI_Get_count(&probe_status, MPI_DOUBLE, &size);
// layout of buffer is [phreeqc][surrogate]
@ -378,6 +386,8 @@ inline void poet::ChemistryModule::MasterRecvPkgs(worker_list_t &w_list,
w_list[p - 1].has_work = 0;
pkg_to_recv -= 1;
free_workers++;
recv_b = MPI_Wtime();
this->recv_ctrl_t += recv_b - recv_a;
}
}
}
@ -432,6 +442,10 @@ void poet::ChemistryModule::MasterRunParallel(double dt) {
int i_pkgs;
int ftype;
double ctrl_a, ctrl_b;
double worker_ctrl_a, worker_ctrl_b;
double ctrl_bcast_a, ctrl_bcast_b;
const std::vector<uint32_t> wp_sizes_vector =
CalculateWPSizesVector(this->n_cells, this->wp_size);
@ -445,28 +459,44 @@ void poet::ChemistryModule::MasterRunParallel(double dt) {
MPI_INT);
}
ftype = CHEM_WORK_LOOP;
PropagateFunctionType(ftype);
/* start time measurement of broadcasting interpolation status */
ctrl_bcast_a = MPI_Wtime();
ftype = CHEM_INTERP;
PropagateFunctionType(ftype);
if(this->runtime_params->rollback_simulation){
int interp_flag = 0;
int dht_fill_flag = 0;
if(this->runtime_params->rollback_enabled){
this->interp_enabled = false;
int interp_flag = 0;
ChemBCast(&interp_flag, 1, MPI_INT);
} else {
this->dht_fill_during_rollback = true;
interp_flag = 0;
dht_fill_flag = 1;
}
else {
this->interp_enabled = true;
int interp_flag = 1;
ChemBCast(&interp_flag, 1, MPI_INT);
this->dht_fill_during_rollback = false;
interp_flag = 1;
dht_fill_flag = 0;
}
ChemBCast(&interp_flag, 1, MPI_INT);
ChemBCast(&dht_fill_flag, 1, MPI_INT);
/* end time measurement of broadcasting interpolation status */
ctrl_bcast_b = MPI_Wtime();
this->bcast_ctrl_t += ctrl_bcast_b - ctrl_bcast_a;
ftype = CHEM_WORK_LOOP;
PropagateFunctionType(ftype);
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) {
uint32_t control_logic_enabled = this->runtime_params->control_interval_enabled ? 1 : 0;
if (control_logic_enabled) {
sur_shuffled.clear();
sur_shuffled.reserve(this->n_cells * this->prop_count);
}
@ -512,7 +542,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,
i_pkgs, free_workers, dt, iteration, control_logic_enabled,
wp_sizes_vector);
}
// ... and try to receive them from workers who has finished their work
@ -522,39 +552,43 @@ void poet::ChemistryModule::MasterRunParallel(double dt) {
// Just to complete the progressbar
std::cout << std::endl;
/* stop time measurement of chemistry time needed for send/recv loop */
worker_chemistry_b = MPI_Wtime();
this->send_recv_t += worker_chemistry_b - worker_chemistry_a;
/* stop time measurement of chemistry time needed for send/recv loop */
worker_chemistry_b = MPI_Wtime();
this->send_recv_t += worker_chemistry_b - worker_chemistry_a;
/* start time measurement of sequential part */
seq_c = MPI_Wtime();
/* start time measurement of sequential part */
seq_c = MPI_Wtime();
/* unshuffle grid */
// grid.importAndUnshuffle(mpi_buffer);
std::vector<double> out_vec{mpi_buffer};
unshuffleField(mpi_buffer, this->n_cells, this->prop_count,
wp_sizes_vector.size(), out_vec);
chem_field = out_vec;
/* unshuffle grid */
// grid.importAndUnshuffle(mpi_buffer);
std::vector<double> out_vec{mpi_buffer};
unshuffleField(mpi_buffer, this->n_cells, this->prop_count,
wp_sizes_vector.size(), out_vec);
chem_field = out_vec;
/* do master stuff */
/* do master stuff */
if (control_iteration) {
control_iteration_counter++;
/* start time measurement of control logic */
ctrl_a = MPI_Wtime();
std::vector<double> sur_unshuffled{sur_shuffled};
if (control_logic_enabled && !this->runtime_params->rollback_enabled) {
unshuffleField(sur_shuffled, this->n_cells, this->prop_count,
wp_sizes_vector.size(), sur_unshuffled);
std::vector<double> sur_unshuffled{sur_shuffled};;
error_stats stats(this->prop_count, control_iteration_counter *
runtime_params->control_iteration);
unshuffleField(sur_shuffled, this->n_cells, this->prop_count,
wp_sizes_vector.size(), sur_unshuffled);
computeStats(out_vec, sur_unshuffled, this->n_cells, this->prop_count,
stats);
error_stats_history.push_back(stats);
SimulationErrorStats stats(this->prop_count, this->runtime_params->global_iter, this->runtime_params->rollback_counter);
computeSpeciesErrors(out_vec, sur_unshuffled, this->n_cells, this->prop_count, stats);
error_history.push_back(stats);
}
/* end time measurement of control logic */
ctrl_b = MPI_Wtime();
this->ctrl_t += ctrl_b - ctrl_a;
// to do: control values to epsilon
}
/* start time measurement of master chemistry */
sim_e_chemistry = MPI_Wtime();

538
src/Chemistry/SurrogateModels/DHT_Wrapper.cpp
View File

@ -36,315 +36,357 @@
using namespace std;
namespace poet {
namespace poet
{
DHT_Wrapper::DHT_Wrapper(MPI_Comm dht_comm, std::uint64_t dht_size,
const NamedVector<std::uint32_t> &key_species,
const std::vector<std::int32_t> &key_indices,
const std::vector<std::string> &_output_names,
const InitialList::ChemistryHookFunctions &_hooks,
uint32_t data_count, bool _with_interp,
bool _has_het_ids)
: key_count(key_indices.size()), data_count(data_count),
input_key_elements(key_indices), communicator(dht_comm),
key_species(key_species), output_names(_output_names), hooks(_hooks),
with_interp(_with_interp), has_het_ids(_has_het_ids) {
// initialize DHT object
// key size = count of key elements + timestep
uint32_t key_size = (key_count + 1) * sizeof(Lookup_Keyelement);
uint32_t data_size =
(data_count + (with_interp ? input_key_elements.size() : 0)) *
sizeof(double);
uint32_t buckets_per_process =
static_cast<std::uint32_t>(dht_size / (data_size + key_size));
dht_object = DHT_create(dht_comm, buckets_per_process, data_size, key_size,
&poet::Murmur2_64A);
DHT_Wrapper::DHT_Wrapper(MPI_Comm dht_comm, std::uint64_t dht_size,
const NamedVector<std::uint32_t> &key_species,
const std::vector<std::int32_t> &key_indices,
const std::vector<std::string> &_output_names,
const InitialList::ChemistryHookFunctions &_hooks,
uint32_t data_count, bool _with_interp,
bool _has_het_ids)
: key_count(key_indices.size()), data_count(data_count),
input_key_elements(key_indices), communicator(dht_comm),
key_species(key_species), output_names(_output_names), hooks(_hooks),
with_interp(_with_interp), has_het_ids(_has_het_ids)
{
// initialize DHT object
// key size = count of key elements + timestep
uint32_t key_size = (key_count + 1) * sizeof(Lookup_Keyelement);
uint32_t data_size =
(data_count + (with_interp ? input_key_elements.size() : 0)) *
sizeof(double);
uint32_t buckets_per_process =
static_cast<std::uint32_t>(dht_size / (data_size + key_size));
dht_object = DHT_create(dht_comm, buckets_per_process, data_size, key_size,
&poet::Murmur2_64A);
dht_signif_vector = key_species.getValues();
dht_signif_vector = key_species.getValues();
// this->dht_signif_vector.resize(key_size, DHT_KEY_SIGNIF_DEFAULT);
// this->dht_signif_vector.resize(key_size, DHT_KEY_SIGNIF_DEFAULT);
this->dht_prop_type_vector.resize(key_count, DHT_TYPE_DEFAULT);
this->dht_prop_type_vector.resize(key_count, DHT_TYPE_DEFAULT);
const auto key_names = key_species.getNames();
const auto key_names = key_species.getNames();
auto tot_h = std::find(key_names.begin(), key_names.end(), "H");
if (tot_h != key_names.end()) {
this->dht_prop_type_vector[tot_h - key_names.begin()] = DHT_TYPE_TOTAL;
}
auto tot_o = std::find(key_names.begin(), key_names.end(), "O");
if (tot_o != key_names.end()) {
this->dht_prop_type_vector[tot_o - key_names.begin()] = DHT_TYPE_TOTAL;
}
auto charge = std::find(key_names.begin(), key_names.end(), "Charge");
if (charge != key_names.end()) {
this->dht_prop_type_vector[charge - key_names.begin()] = DHT_TYPE_CHARGE;
}
}
DHT_Wrapper::~DHT_Wrapper() {
// free DHT
DHT_free(dht_object, NULL, NULL);
}
auto DHT_Wrapper::checkDHT(WorkPackage &work_package)
-> const DHT_ResultObject & {
const auto length = work_package.size;
std::vector<double> bucket_writer(
this->data_count + (with_interp ? input_key_elements.size() : 0));
// loop over every grid cell contained in work package
for (int i = 0; i < length; i++) {
// point to current grid cell
auto &key_vector = dht_results.keys[i];
// overwrite input with data from DHT, IF value is found in DHT
int res =
DHT_read(this->dht_object, key_vector.data(), bucket_writer.data());
switch (res) {
case DHT_SUCCESS:
work_package.output[i] =
(with_interp
? inputAndRatesToOutput(bucket_writer, work_package.input[i])
: bucket_writer);
work_package.mapping[i] = CHEM_DHT;
this->dht_hits++;
break;
case DHT_READ_MISS:
break;
auto tot_h = std::find(key_names.begin(), key_names.end(), "H");
if (tot_h != key_names.end())
{
this->dht_prop_type_vector[tot_h - key_names.begin()] = DHT_TYPE_TOTAL;
}
auto tot_o = std::find(key_names.begin(), key_names.end(), "O");
if (tot_o != key_names.end())
{
this->dht_prop_type_vector[tot_o - key_names.begin()] = DHT_TYPE_TOTAL;
}
auto charge = std::find(key_names.begin(), key_names.end(), "Charge");
if (charge != key_names.end())
{
this->dht_prop_type_vector[charge - key_names.begin()] = DHT_TYPE_CHARGE;
}
}
return dht_results;
}
DHT_Wrapper::~DHT_Wrapper()
{
// free DHT
DHT_free(dht_object, NULL, NULL);
}
auto DHT_Wrapper::checkDHT(WorkPackage &work_package)
-> const DHT_ResultObject &
{
void DHT_Wrapper::fillDHT(const WorkPackage &work_package) {
const auto length = work_package.size;
const auto length = work_package.size;
std::vector<double> bucket_writer(
this->data_count + (with_interp ? input_key_elements.size() : 0));
// loop over every grid cell contained in work package
dht_results.locations.resize(length);
dht_results.filledDHT = std::vector<bool>(length, false);
for (int i = 0; i < length; i++) {
// If true grid cell was simulated, needs to be inserted into dht
if (work_package.mapping[i] != CHEM_PQC) {
continue;
}
// loop over every grid cell contained in work package
for (int i = 0; i < length; i++)
{
// point to current grid cell
auto &key_vector = dht_results.keys[i];
if (work_package.input[i][0] != 2) {
continue;
}
// overwrite input with data from DHT, IF value is found in DHT
int res =
DHT_read(this->dht_object, key_vector.data(), bucket_writer.data());
// check if calcite or dolomite is absent and present, resp.n and vice
// versa in input/output. If this is the case -> Do not write to DHT!
// HACK: hardcoded, should be fixed!
if (hooks.dht_fill.isValid()) {
NamedVector<double> old_values(output_names, work_package.input[i]);
NamedVector<double> new_values(output_names, work_package.output[i]);
if (Rcpp::as<bool>(hooks.dht_fill(old_values, new_values))) {
continue;
switch (res)
{
case DHT_SUCCESS:
work_package.output[i] =
(with_interp
? inputAndRatesToOutput(bucket_writer, work_package.input[i])
: bucket_writer);
work_package.mapping[i] = CHEM_DHT;
this->dht_hits++;
break;
case DHT_READ_MISS:
break;
}
}
uint32_t proc, index;
auto &key = dht_results.keys[i];
const auto data =
(with_interp ? outputToInputAndRates(work_package.input[i],
work_package.output[i])
: work_package.output[i]);
// void *data = (void *)&(work_package[i * this->data_count]);
// fuzz data (round, logarithm etc.)
return dht_results;
}
// insert simulated data with fuzzed key into DHT
int res = DHT_write(this->dht_object, key.data(),
const_cast<double *>(data.data()), &proc, &index);
void DHT_Wrapper::fillDHT(const WorkPackage &work_package)
{
dht_results.locations[i] = {proc, index};
const auto length = work_package.size;
// if data was successfully written ...
if ((res != DHT_SUCCESS) && (res == DHT_WRITE_SUCCESS_WITH_EVICTION)) {
dht_evictions++;
// loop over every grid cell contained in work package
dht_results.locations.resize(length);
dht_results.filledDHT = std::vector<bool>(length, false);
for (int i = 0; i < length; i++)
{
// If true grid cell was simulated, needs to be inserted into dht
if (work_package.mapping[i] != CHEM_PQC)
{
continue;
}
if (work_package.input[i][1] != 2)
{
continue;
}
// check if calcite or dolomite is absent and present, resp.n and vice
// versa in input/output. If this is the case -> Do not write to DHT!
// HACK: hardcoded, should be fixed!
if (hooks.dht_fill.isValid())
{
NamedVector<double> old_values(output_names, work_package.input[i]);
NamedVector<double> new_values(output_names, work_package.output[i]);
if (Rcpp::as<bool>(hooks.dht_fill(old_values, new_values)))
{
continue;
}
}
uint32_t proc, index;
auto &key = dht_results.keys[i];
const auto data =
(with_interp ? outputToInputAndRates(work_package.input[i],
work_package.output[i])
: work_package.output[i]);
// void *data = (void *)&(work_package[i * this->data_count]);
// fuzz data (round, logarithm etc.)
// insert simulated data with fuzzed key into DHT
int res = DHT_write(this->dht_object, key.data(),
const_cast<double *>(data.data()), &proc, &index);
dht_results.locations[i] = {proc, index};
// if data was successfully written ...
if ((res != DHT_SUCCESS) && (res == DHT_WRITE_SUCCESS_WITH_EVICTION))
{
dht_evictions++;
}
dht_results.filledDHT[i] = true;
}
}
inline std::vector<double>
DHT_Wrapper::outputToInputAndRates(const std::vector<double> &old_results,
const std::vector<double> &new_results)
{
const int prefix_size = this->input_key_elements.size();
std::vector<double> output(prefix_size + this->data_count);
std::copy(new_results.begin(), new_results.end(),
output.begin() + prefix_size);
for (int i = 0; i < prefix_size; i++)
{
const int data_elem_i = input_key_elements[i];
output[i] = old_results[data_elem_i];
output[prefix_size + data_elem_i] -= old_results[data_elem_i];
}
dht_results.filledDHT[i] = true;
}
}
inline std::vector<double>
DHT_Wrapper::outputToInputAndRates(const std::vector<double> &old_results,
const std::vector<double> &new_results) {
const int prefix_size = this->input_key_elements.size();
std::vector<double> output(prefix_size + this->data_count);
std::copy(new_results.begin(), new_results.end(),
output.begin() + prefix_size);
for (int i = 0; i < prefix_size; i++) {
const int data_elem_i = input_key_elements[i];
output[i] = old_results[data_elem_i];
output[prefix_size + data_elem_i] -= old_results[data_elem_i];
return output;
}
return output;
}
inline std::vector<double>
DHT_Wrapper::inputAndRatesToOutput(const std::vector<double> &dht_data,
const std::vector<double> &input_values)
{
const int prefix_size = this->input_key_elements.size();
inline std::vector<double>
DHT_Wrapper::inputAndRatesToOutput(const std::vector<double> &dht_data,
const std::vector<double> &input_values) {
const int prefix_size = this->input_key_elements.size();
std::vector<double> output(input_values);
std::vector<double> output(input_values);
for (int i = 0; i < prefix_size; i++)
{
const int data_elem_i = input_key_elements[i];
output[data_elem_i] += dht_data[i];
}
for (int i = 0; i < prefix_size; i++) {
const int data_elem_i = input_key_elements[i];
output[data_elem_i] += dht_data[i];
return output;
}
return output;
}
inline std::vector<double>
DHT_Wrapper::outputToRates(const std::vector<double> &old_results,
const std::vector<double> &new_results)
{
std::vector<double> output(new_results);
inline std::vector<double>
DHT_Wrapper::outputToRates(const std::vector<double> &old_results,
const std::vector<double> &new_results) {
std::vector<double> output(new_results);
for (const auto &data_elem_i : input_key_elements)
{
output[data_elem_i] -= old_results[data_elem_i];
}
for (const auto &data_elem_i : input_key_elements) {
output[data_elem_i] -= old_results[data_elem_i];
return output;
}
return output;
}
inline std::vector<double>
DHT_Wrapper::ratesToOutput(const std::vector<double> &dht_data,
const std::vector<double> &input_values)
{
std::vector<double> output(input_values);
inline std::vector<double>
DHT_Wrapper::ratesToOutput(const std::vector<double> &dht_data,
const std::vector<double> &input_values) {
std::vector<double> output(input_values);
for (const auto &data_elem_i : input_key_elements)
{
output[data_elem_i] += dht_data[data_elem_i];
}
for (const auto &data_elem_i : input_key_elements) {
output[data_elem_i] += dht_data[data_elem_i];
return output;
}
return output;
}
// void DHT_Wrapper::resultsToWP(std::vector<double> &work_package) {
// for (int i = 0; i < dht_results.length; i++) {
// if (!dht_results.needPhreeqc[i]) {
// std::copy(dht_results.results[i].begin(), dht_results.results[i].end(),
// work_package.begin() + (data_count * i));
// }
// }
// }
// void DHT_Wrapper::resultsToWP(std::vector<double> &work_package) {
// for (int i = 0; i < dht_results.length; i++) {
// if (!dht_results.needPhreeqc[i]) {
// std::copy(dht_results.results[i].begin(), dht_results.results[i].end(),
// work_package.begin() + (data_count * i));
// }
// }
// }
int DHT_Wrapper::tableToFile(const char *filename) {
int res = DHT_to_file(dht_object, filename);
return res;
}
int DHT_Wrapper::fileToTable(const char *filename) {
int res = DHT_from_file(dht_object, filename);
if (res != DHT_SUCCESS)
int DHT_Wrapper::tableToFile(const char *filename)
{
int res = DHT_to_file(dht_object, filename);
return res;
}
int DHT_Wrapper::fileToTable(const char *filename)
{
int res = DHT_from_file(dht_object, filename);
if (res != DHT_SUCCESS)
return res;
#ifdef DHT_STATISTICS
DHT_print_statistics(dht_object);
DHT_print_statistics(dht_object);
#endif
return DHT_SUCCESS;
}
void DHT_Wrapper::printStatistics() {
int res;
res = DHT_print_statistics(dht_object);
if (res != DHT_SUCCESS) {
// MPI ERROR ... WHAT TO DO NOW?
// RUNNING CIRCLES WHILE SCREAMING
return DHT_SUCCESS;
}
}
LookupKey DHT_Wrapper::fuzzForDHT_R(const std::vector<double> &cell,
double dt) {
const auto c_zero_val = std::pow(10, AQUEOUS_EXP);
void DHT_Wrapper::printStatistics()
{
int res;
NamedVector<double> input_nv(this->output_names, cell);
res = DHT_print_statistics(dht_object);
const std::vector<double> eval_vec =
Rcpp::as<std::vector<double>>(hooks.dht_fuzz(input_nv));
assert(eval_vec.size() == this->key_count);
LookupKey vecFuzz(this->key_count + 1 + has_het_ids, {.0});
if (res != DHT_SUCCESS)
{
// MPI ERROR ... WHAT TO DO NOW?
// RUNNING CIRCLES WHILE SCREAMING
}
}
DHT_Rounder rounder;
LookupKey DHT_Wrapper::fuzzForDHT_R(const std::vector<double> &cell,
double dt)
{
const auto c_zero_val = std::pow(10, AQUEOUS_EXP);
int totals_i = 0;
// introduce fuzzing to allow more hits in DHT
// loop over every variable of grid cell
for (std::uint32_t i = 0; i < eval_vec.size(); i++) {
double curr_key = eval_vec[i];
if (curr_key != 0) {
if (this->dht_prop_type_vector[i] == DHT_TYPE_TOTAL) {
curr_key -= base_totals[totals_i++];
NamedVector<double> input_nv(this->output_names, cell);
const std::vector<double> eval_vec =
Rcpp::as<std::vector<double>>(hooks.dht_fuzz(input_nv));
assert(eval_vec.size() == this->key_count);
LookupKey vecFuzz(this->key_count + 1 + has_het_ids, {.0});
DHT_Rounder rounder;
int totals_i = 0;
// introduce fuzzing to allow more hits in DHT
// loop over every variable of grid cell
for (std::uint32_t i = 0; i < eval_vec.size(); i++)
{
double curr_key = eval_vec[i];
if (curr_key != 0)
{
if (this->dht_prop_type_vector[i] == DHT_TYPE_TOTAL)
{
curr_key -= base_totals[totals_i++];
}
vecFuzz[i] =
rounder.round(curr_key, dht_signif_vector[i],
this->dht_prop_type_vector[i] == DHT_TYPE_TOTAL);
}
vecFuzz[i] =
rounder.round(curr_key, dht_signif_vector[i],
this->dht_prop_type_vector[i] == DHT_TYPE_TOTAL);
}
}
// add timestep to the end of the key as double value
vecFuzz[this->key_count].fp_element = dt;
if (has_het_ids) {
vecFuzz[this->key_count + 1].fp_element = cell[0];
// add timestep to the end of the key as double value
vecFuzz[this->key_count].fp_element = dt;
if (has_het_ids)
{
vecFuzz[this->key_count + 1].fp_element = cell[0];
}
return vecFuzz;
}
return vecFuzz;
}
LookupKey DHT_Wrapper::fuzzForDHT(const std::vector<double> &cell, double dt)
{
const auto c_zero_val = std::pow(10, AQUEOUS_EXP);
LookupKey DHT_Wrapper::fuzzForDHT(const std::vector<double> &cell, double dt) {
const auto c_zero_val = std::pow(10, AQUEOUS_EXP);
LookupKey vecFuzz(this->key_count + 1 + has_het_ids, {.0});
DHT_Rounder rounder;
LookupKey vecFuzz(this->key_count + 1 + has_het_ids, {.0});
DHT_Rounder rounder;
int totals_i = 0;
// introduce fuzzing to allow more hits in DHT
// loop over every variable of grid cell
for (std::uint32_t i = 0; i < input_key_elements.size(); i++) {
if (input_key_elements[i] == DHT_KEY_INPUT_CUSTOM) {
continue;
}
double curr_key = cell[input_key_elements[i]];
if (curr_key != 0) {
if (curr_key < c_zero_val &&
this->dht_prop_type_vector[i] == DHT_TYPE_DEFAULT) {
int totals_i = 0;
// introduce fuzzing to allow more hits in DHT
// loop over every variable of grid cell
for (std::uint32_t i = 0; i < input_key_elements.size(); i++)
{
if (input_key_elements[i] == DHT_KEY_INPUT_CUSTOM)
{
continue;
}
if (this->dht_prop_type_vector[i] == DHT_TYPE_TOTAL) {
curr_key -= base_totals[totals_i++];
double curr_key = cell[input_key_elements[i]];
if (curr_key != 0)
{
if (curr_key < c_zero_val &&
this->dht_prop_type_vector[i] == DHT_TYPE_DEFAULT)
{
continue;
}
if (this->dht_prop_type_vector[i] == DHT_TYPE_TOTAL)
{
curr_key -= base_totals[totals_i++];
}
vecFuzz[i] =
rounder.round(curr_key, dht_signif_vector[i],
this->dht_prop_type_vector[i] == DHT_TYPE_TOTAL);
}
vecFuzz[i] =
rounder.round(curr_key, dht_signif_vector[i],
this->dht_prop_type_vector[i] == DHT_TYPE_TOTAL);
}
}
// add timestep to the end of the key as double value
vecFuzz[this->key_count].fp_element = dt;
if (has_het_ids) {
vecFuzz[this->key_count + 1].fp_element = cell[0];
// add timestep to the end of the key as double value
vecFuzz[this->key_count].fp_element = dt;
if (has_het_ids)
{
vecFuzz[this->key_count + 1].fp_element = cell[0];
}
return vecFuzz;
}
return vecFuzz;
}
void poet::DHT_Wrapper::SetSignifVector(std::vector<uint32_t> signif_vec)
{
if (signif_vec.size() != this->key_count)
{
throw std::runtime_error(
"Significant vector size mismatches count of key elements.");
}
void poet::DHT_Wrapper::SetSignifVector(std::vector<uint32_t> signif_vec) {
if (signif_vec.size() != this->key_count) {
throw std::runtime_error(
"Significant vector size mismatches count of key elements.");
this->dht_signif_vector = signif_vec;
}
this->dht_signif_vector = signif_vec;
}
} // namespace poet

255
src/Chemistry/SurrogateModels/InterpolationModule.cpp
View File

@ -25,152 +25,175 @@
#include <utility>
#include <vector>
extern "C" {
extern "C"
{
#include "DHT.h"
}
namespace poet {
namespace poet
{
InterpolationModule::InterpolationModule(
std::uint32_t entries_per_bucket, std::uint64_t size_per_process,
std::uint32_t min_entries_needed, DHT_Wrapper &dht,
const NamedVector<std::uint32_t> &interp_key_signifs,
const std::vector<std::int32_t> &dht_key_indices,
const std::vector<std::string> &_out_names,
const InitialList::ChemistryHookFunctions &_hooks)
: dht_instance(dht), key_signifs(interp_key_signifs),
key_indices(dht_key_indices), min_entries_needed(min_entries_needed),
dht_names(dht.getKeySpecies().getNames()), out_names(_out_names),
hooks(_hooks) {
InterpolationModule::InterpolationModule(
std::uint32_t entries_per_bucket, std::uint64_t size_per_process,
std::uint32_t min_entries_needed, DHT_Wrapper &dht,
const NamedVector<std::uint32_t> &interp_key_signifs,
const std::vector<std::int32_t> &dht_key_indices,
const std::vector<std::string> &_out_names,
const InitialList::ChemistryHookFunctions &_hooks)
: dht_instance(dht), key_signifs(interp_key_signifs),
key_indices(dht_key_indices), min_entries_needed(min_entries_needed),
dht_names(dht.getKeySpecies().getNames()), out_names(_out_names),
hooks(_hooks)
{
initPHT(this->key_signifs.size(), entries_per_bucket, size_per_process,
dht.getCommunicator());
initPHT(this->key_signifs.size(), entries_per_bucket, size_per_process,
dht.getCommunicator());
pht->setSourceDHT(dht.getDHT());
}
pht->setSourceDHT(dht.getDHT());
}
void InterpolationModule::initPHT(std::uint32_t key_count,
std::uint32_t entries_per_bucket,
std::uint32_t size_per_process,
MPI_Comm communicator) {
uint32_t key_size = key_count * sizeof(Lookup_Keyelement) + sizeof(double);
uint32_t data_size = sizeof(DHT_Location);
void InterpolationModule::initPHT(std::uint32_t key_count,
std::uint32_t entries_per_bucket,
std::uint32_t size_per_process,
MPI_Comm communicator)
{
uint32_t key_size = key_count * sizeof(Lookup_Keyelement) + sizeof(double);
uint32_t data_size = sizeof(DHT_Location);
pht = std::make_unique<ProximityHashTable>(
key_size, data_size, entries_per_bucket, size_per_process, communicator);
}
pht = std::make_unique<ProximityHashTable>(
key_size, data_size, entries_per_bucket, size_per_process, communicator);
}
void InterpolationModule::writePairs() {
const auto in = this->dht_instance.getDHTResults();
for (int i = 0; i < in.filledDHT.size(); i++) {
if (in.filledDHT[i]) {
const auto coarse_key = roundKey(in.keys[i]);
pht->writeLocationToPHT(coarse_key, in.locations[i]);
void InterpolationModule::writePairs()
{
const auto in = this->dht_instance.getDHTResults();
for (int i = 0; i < in.filledDHT.size(); i++)
{
if (in.filledDHT[i])
{
const auto coarse_key = roundKey(in.keys[i]);
pht->writeLocationToPHT(coarse_key, in.locations[i]);
}
}
}
}
void InterpolationModule::tryInterpolation(WorkPackage &work_package) {
interp_result.status.resize(work_package.size, NOT_NEEDED);
void InterpolationModule::tryInterpolation(WorkPackage &work_package)
{
interp_result.status.resize(work_package.size, NOT_NEEDED);
const auto dht_results = this->dht_instance.getDHTResults();
const auto dht_results = this->dht_instance.getDHTResults();
for (int wp_i = 0; wp_i < work_package.size; wp_i++) {
if (work_package.input[wp_i][0] != 2) {
interp_result.status[wp_i] = INSUFFICIENT_DATA;
continue;
}
if (work_package.mapping[wp_i] != CHEM_PQC) {
interp_result.status[wp_i] = NOT_NEEDED;
continue;
}
const auto rounded_key = roundKey(dht_results.keys[wp_i]);
auto pht_result =
pht->query(rounded_key, this->min_entries_needed,
dht_instance.getInputCount(), dht_instance.getOutputCount());
if (pht_result.size < this->min_entries_needed) {
interp_result.status[wp_i] = INSUFFICIENT_DATA;
continue;
}
if (hooks.interp_pre.isValid()) {
NamedVector<double> nv_in(this->out_names, work_package.input[wp_i]);
std::vector<int> rm_indices = Rcpp::as<std::vector<int>>(
hooks.interp_pre(nv_in, pht_result.in_values));
pht_result.size -= rm_indices.size();
if (pht_result.size < this->min_entries_needed) {
for (int wp_i = 0; wp_i < work_package.size; wp_i++)
{
if (work_package.input[wp_i][1] != 2)
{
interp_result.status[wp_i] = INSUFFICIENT_DATA;
continue;
}
for (const auto &index : rm_indices) {
pht_result.in_values.erase(
std::next(pht_result.in_values.begin(), index - 1));
pht_result.out_values.erase(
std::next(pht_result.out_values.begin(), index - 1));
if (work_package.mapping[wp_i] != CHEM_PQC)
{
interp_result.status[wp_i] = NOT_NEEDED;
continue;
}
}
#ifdef POET_PHT_ADD
this->pht->incrementReadCounter(roundKey(rounded_key));
#endif
const auto rounded_key = roundKey(dht_results.keys[wp_i]);
const int cell_id = static_cast<int>(work_package.input[wp_i][0]);
auto pht_result =
pht->query(rounded_key, this->min_entries_needed,
dht_instance.getInputCount(), dht_instance.getOutputCount());
if (!to_calc_cache.contains(cell_id)) {
const std::vector<std::int32_t> &to_calc = dht_instance.getKeyElements();
std::vector<std::int32_t> keep_indices;
if (pht_result.size < this->min_entries_needed)
{
interp_result.status[wp_i] = INSUFFICIENT_DATA;
continue;
}
for (std::size_t i = 0; i < to_calc.size(); i++) {
if (!std::isnan(work_package.input[wp_i][to_calc[i]])) {
keep_indices.push_back(to_calc[i]);
if (hooks.interp_pre.isValid())
{
NamedVector<double> nv_in(this->out_names, work_package.input[wp_i]);
std::vector<int> rm_indices = Rcpp::as<std::vector<int>>(
hooks.interp_pre(nv_in, pht_result.in_values));
pht_result.size -= rm_indices.size();
if (pht_result.size < this->min_entries_needed)
{
interp_result.status[wp_i] = INSUFFICIENT_DATA;
continue;
}
for (const auto &index : rm_indices)
{
pht_result.in_values.erase(
std::next(pht_result.in_values.begin(), index - 1));
pht_result.out_values.erase(
std::next(pht_result.out_values.begin(), index - 1));
}
}
to_calc_cache[cell_id] = keep_indices;
#ifdef POET_PHT_ADD
this->pht->incrementReadCounter(roundKey(rounded_key));
#endif
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();
std::vector<std::int32_t> keep_indices;
for (std::size_t i = 0; i < to_calc.size(); i++)
{
if (!std::isnan(work_package.input[wp_i][to_calc[i]]))
{
keep_indices.push_back(to_calc[i]);
}
}
to_calc_cache[cell_id] = keep_indices;
}
double start_fc = MPI_Wtime();
work_package.output[wp_i] =
f_interpolate(to_calc_cache[cell_id], work_package.input[wp_i],
pht_result.in_values, pht_result.out_values);
if (hooks.interp_post.isValid())
{
NamedVector<double> nv_result(this->out_names, work_package.output[wp_i]);
if (Rcpp::as<bool>(hooks.interp_post(nv_result)))
{
interp_result.status[wp_i] = INSUFFICIENT_DATA;
continue;
}
}
// interp_result.results[i][0] = mean_water;
this->interp_t += MPI_Wtime() - start_fc;
this->interpolations++;
work_package.mapping[wp_i] = CHEM_INTERP;
interp_result.status[wp_i] = RES_OK;
}
}
double start_fc = MPI_Wtime();
work_package.output[wp_i] =
f_interpolate(to_calc_cache[cell_id], work_package.input[wp_i],
pht_result.in_values, pht_result.out_values);
if (hooks.interp_post.isValid()) {
NamedVector<double> nv_result(this->out_names, work_package.output[wp_i]);
if (Rcpp::as<bool>(hooks.interp_post(nv_result))) {
interp_result.status[wp_i] = INSUFFICIENT_DATA;
continue;
void InterpolationModule::resultsToWP(std::vector<double> &work_package) const
{
for (uint32_t i = 0; i < interp_result.status.size(); i++)
{
if (interp_result.status[i] == RES_OK)
{
const std::size_t length =
interp_result.results[i].end() - interp_result.results[i].begin();
std::copy(interp_result.results[i].begin(),
interp_result.results[i].end(),
work_package.begin() + (length * i));
}
}
// interp_result.results[i][0] = mean_water;
this->interp_t += MPI_Wtime() - start_fc;
this->interpolations++;
work_package.mapping[wp_i] = CHEM_INTERP;
interp_result.status[wp_i] = RES_OK;
}
}
void InterpolationModule::resultsToWP(std::vector<double> &work_package) const {
for (uint32_t i = 0; i < interp_result.status.size(); i++) {
if (interp_result.status[i] == RES_OK) {
const std::size_t length =
interp_result.results[i].end() - interp_result.results[i].begin();
std::copy(interp_result.results[i].begin(),
interp_result.results[i].end(),
work_package.begin() + (length * i));
}
}
}
} // namespace poet

40
src/Chemistry/WorkerFunctions.cpp
View File

@ -69,9 +69,14 @@ namespace poet
}
case CHEM_INTERP:
{
int interp_flag;
int interp_flag = 0;
int dht_fill_flag = 0;
ChemBCast(&interp_flag, 1, MPI_INT);
ChemBCast(&dht_fill_flag, 1, MPI_INT);
this->interp_enabled = (interp_flag == 1);
this->dht_fill_during_rollback = (dht_fill_flag == 1);
break;
}
case CHEM_WORK_LOOP:
@ -150,13 +155,14 @@ namespace poet
double dht_get_start, dht_get_end;
double phreeqc_time_start, phreeqc_time_end;
double dht_fill_start, dht_fill_end;
double ctrl_time_c, ctrl_time_d;
uint32_t iteration;
double dt;
double current_sim_time;
uint32_t wp_start_index;
int count = double_count;
bool control_iteration_active = false;
bool control_logic_enabled = false;
std::vector<double> mpi_buffer(count);
/* receive */
@ -183,7 +189,7 @@ namespace poet
// current work package start location in field
wp_start_index = mpi_buffer[count + 4];
control_iteration_active = (mpi_buffer[count + 5] == 1);
control_logic_enabled = (mpi_buffer[count + 5] == 1);
for (std::size_t wp_i = 0; wp_i < s_curr_wp.size; wp_i++)
{
@ -229,7 +235,7 @@ namespace poet
poet::WorkPackage s_curr_wp_control = s_curr_wp;
if (control_iteration_active)
if (control_logic_enabled)
{
for (std::size_t wp_i = 0; wp_i < s_curr_wp_control.size; wp_i++)
{
@ -240,12 +246,15 @@ namespace poet
phreeqc_time_start = MPI_Wtime();
WorkerRunWorkPackage(control_iteration_active ? s_curr_wp_control : s_curr_wp, current_sim_time, dt);
WorkerRunWorkPackage(control_logic_enabled ? s_curr_wp_control : s_curr_wp, current_sim_time, dt);
phreeqc_time_end = MPI_Wtime();
if (control_iteration_active)
{
if (control_logic_enabled)
{
/* start time measurement for copying control workpackage */
ctrl_time_c = MPI_Wtime();
std::size_t sur_wp_offset = s_curr_wp.size * this->prop_count;
mpi_buffer.resize(count + sur_wp_offset);
@ -275,6 +284,10 @@ namespace poet
}
count += sur_wp_offset;
/* end time measurement for copying control workpackage */
ctrl_time_d = MPI_Wtime();
timings.ctrl_t += ctrl_time_d - ctrl_time_c;
}
else
{
@ -288,14 +301,14 @@ namespace poet
/* send results to master */
MPI_Request send_req;
int mpi_tag = control_iteration_active ? LOOP_CTRL : LOOP_WORK;
int mpi_tag = control_logic_enabled ? LOOP_CTRL : LOOP_WORK;
MPI_Isend(mpi_buffer.data(), count, MPI_DOUBLE, 0, mpi_tag, MPI_COMM_WORLD, &send_req);
if (dht_enabled || interp_enabled)
if (dht_enabled || interp_enabled || dht_fill_during_rollback)
{
/* write results to DHT */
dht_fill_start = MPI_Wtime();
dht->fillDHT(control_iteration_active ? s_curr_wp_control : s_curr_wp);
dht->fillDHT(control_logic_enabled ? s_curr_wp_control : s_curr_wp);
dht_fill_end = MPI_Wtime();
if (interp_enabled)
@ -306,7 +319,6 @@ namespace poet
}
timings.phreeqc_t += phreeqc_time_end - phreeqc_time_start;
MPI_Wait(&send_req, MPI_STATUS_IGNORE);
}
@ -460,6 +472,12 @@ namespace poet
this->group_comm);
break;
}
case WORKER_CTRL_ITER:
{
MPI_Gather(&timings.ctrl_t, 1, MPI_DOUBLE, NULL, 1, MPI_DOUBLE, 0,
this->group_comm);
break;
}
case WORKER_DHT_GET:
{
MPI_Gather(&timings.dht_get, 1, MPI_DOUBLE, NULL, 1, MPI_DOUBLE, 0,

28
src/IO/StatsIO.cpp
View File

@ -2,10 +2,11 @@
#include <fstream>
#include <iostream>
#include <string>
#include <iomanip> // for std::setw and std::setprecision
namespace poet
{
void writeStatsToCSV(const std::vector<ChemistryModule::error_stats> &all_stats,
void writeStatsToCSV(const std::vector<ChemistryModule::SimulationErrorStats> &all_stats,
const std::vector<std::string> &species_names,
const std::string &filename)
{
@ -17,21 +18,32 @@ namespace poet
}
// header
out << "Iteration, Species, MAPE, RRSME \n";
out << std::left << std::setw(15) << "Iteration"
<< std::setw(15) << "Rollback"
<< std::setw(15) << "Species"
<< std::setw(15) << "MAPE"
<< std::setw(15) << "RRSME" << "\n";
out << std::string(75, '-') << "\n"; // separator line
// data rows
for (size_t i = 0; i < all_stats.size(); ++i)
{
for (size_t j = 0; j < species_names.size(); ++j)
{
out << all_stats[i].iteration << ",\t"
<< species_names[j] << ",\t"
<< all_stats[i].mape[j] << ",\t"
<< all_stats[i].rrsme[j] << "\n";
out << std::left
<< std::setw(15) << all_stats[i].iteration
<< std::setw(15) << all_stats[i].rollback_count
<< std::setw(15) << species_names[j]
<< std::setw(15) << all_stats[i].mape[j]
<< std::setw(15) << all_stats[i].rrmse[j]
<< "\n";
}
out << std::endl;
out << "\n"; // blank line between iterations
}
out.close();
std::cout << "Stats written to " << filename << "\n";
}
} // namespace poet
}
// namespace poet

2
src/IO/StatsIO.hpp
View File

@ -3,7 +3,7 @@
namespace poet
{
void writeStatsToCSV(const std::vector<ChemistryModule::error_stats> &all_stats,
void writeStatsToCSV(const std::vector<ChemistryModule::SimulationErrorStats> &all_stats,
const std::vector<std::string> &species_names,
const std::string &filename);
} // namespace poet

118
src/poet.cpp
View File

@ -253,7 +253,6 @@ int parseInitValues(int argc, char **argv, RuntimeParameters &params)
try
{
Rcpp::List init_params_(ReadRObj_R(init_file));
params.init_params = init_params_;
@ -266,14 +265,12 @@ 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"));
params.control_interval =
Rcpp::as<uint32_t>(global_rt_setup->operator[]("control_interval"));
params.checkpoint_interval =
Rcpp::as<uint32_t>(global_rt_setup->operator[]("checkpoint_interval"));
params.mape_threshold =
Rcpp::as<std::vector<double>>(global_rt_setup->operator[]("mape_threshold"));
}
catch (const std::exception &e)
{
@ -304,53 +301,38 @@ 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)
bool triggerRollbackIfExceeded(ChemistryModule &chem, RuntimeParameters &params, uint32_t &current_iteration)
{
const std::vector<double> &latest_mape = chem.error_stats_history.back().mape;
const std::vector<double> &mape_values = chem.error_history.back().mape;
for (uint32_t j = 0; j < params.species_epsilon.size(); j++)
for (uint32_t i = 0; i < params.mape_threshold.size(); i++)
{
if (params.species_epsilon[j] < latest_mape[j] && latest_mape[j] != 0)
// Skip if no meaningful MAPE value
if(mape_values[i] == 0){
continue;
}
if (mape_values[i] > params.mape_threshold[i])
{
uint32_t rollback_iter = iter - (iter % params.control_iteration);
uint32_t rollback_iteration = ((current_iteration - 1) / params.checkpoint_interval) * params.checkpoint_interval;
std::cout << chem.getField().GetProps()[j] << " with a MAPE value of " << latest_mape[j] << " exceeds epsilon of "
<< params.species_epsilon[j] << "! " << std::endl;
MSG("[THRESHOLD EXCEEDED] " + chem.getField().GetProps()[i] + " has MAPE = " +
std::to_string(mape_values[i]) + " exceeding threshold = " + std::to_string(params.mape_threshold[i]) +
" → rolling back to iteration " + std::to_string(rollback_iteration));
Checkpoint_s checkpoint_read{.field = chem.getField()};
read_checkpoint("checkpoint" + std::to_string(rollback_iter) + ".hdf5", checkpoint_read);
iter = checkpoint_read.iteration;
read_checkpoint("checkpoint" + std::to_string(rollback_iteration) + ".hdf5", checkpoint_read);
current_iteration = checkpoint_read.iteration;
// Rollback happend
return true;
}
}
}
MSG("All spezies are below their threshold values");
MSG("All species are within their error thresholds.");
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,
@ -367,21 +349,25 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, RuntimeParameters &params,
}
R["TMP_PROPS"] = Rcpp::wrap(chem.getField().GetProps());
params.next_penalty_check = params.penalty_iteration;
/* SIMULATION LOOP */
double dSimTime{0};
double chkTime = 0.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;
// Rollback countdowm
if (params.rollback_enabled) {
if (params.sur_disabled_counter > 0) {
--params.sur_disabled_counter;
MSG("Rollback counter: " + std::to_string(params.sur_disabled_counter));
} else {
params.rollback_enabled = false;
}
}
params.control_iteration_active = (iter % params.control_iteration == 0 /* && iter != 0 */);
params.global_iter = iter;
params.control_interval_enabled = (iter % params.control_interval == 0);
double start_t = MPI_Wtime();
@ -495,20 +481,27 @@ 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");
double chk_start = MPI_Wtime();
if(iter % params.checkpoint_interval == 0){
MSG("Writing checkpoint of iteration " + std::to_string(iter));
write_checkpoint("checkpoint" + std::to_string(iter) + ".hdf5",
{.field = chem.getField(), .iteration = iter});
}
if (iter == params.next_penalty_check)
if (params.control_interval_enabled && !params.rollback_enabled)
{
bool roolback_happend = checkAndRollback(chem, params, iter);
updatePenaltyLogic(params, roolback_happend);
writeStatsToCSV(chem.error_history, chem.getField().GetProps(), "stats_overview");
params.next_penalty_check = iter + params.penalty_iteration;
if(triggerRollbackIfExceeded(chem, params, iter)){
params.rollback_enabled = true;
params.rollback_counter ++;
params.sur_disabled_counter = params.control_interval;
MSG("Interpolation disabled for the next " + std::to_string(params.control_interval) + ".");
}
}
double chk_end = MPI_Wtime();
chkTime += chk_end - chk_start;
// MSG();
} // END SIMULATION LOOP
@ -526,6 +519,13 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, RuntimeParameters &params,
Rcpp::List diffusion_profiling;
diffusion_profiling["simtime"] = diffusion.getTransportTime();
Rcpp::List ctrl_profiling;
ctrl_profiling["checkpointing_time"] = chkTime;
ctrl_profiling["ctrl_logic_master"] = chem.GetMasterCtrlLogicTime();
ctrl_profiling["bcast_ctrl_logic_master"] = chem.GetMasterCtrlBcastTime();
ctrl_profiling["recv_ctrl_logic_maser"] = chem.GetMasterRecvCtrlLogicTime();
ctrl_profiling["ctrl_logic_worker"] = Rcpp::wrap(chem.GetWorkerControlTimings());
if (params.use_dht)
{
chem_profiling["dht_hits"] = Rcpp::wrap(chem.GetWorkerDHTHits());
@ -554,6 +554,8 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, RuntimeParameters &params,
profiling["simtime"] = dSimTime;
profiling["chemistry"] = chem_profiling;
profiling["diffusion"] = diffusion_profiling;
profiling["ctrl_logic"] = ctrl_profiling;
chem.MasterLoopBreak();

20
src/poet.hpp.in
View File

@ -38,10 +38,10 @@ static const inline std::string ai_surrogate_r_library =
R"(@R_AI_SURROGATE_LIB@)";
static const inline std::string r_runtime_parameters = "mysetup";
struct RuntimeParameters {
struct RuntimeParameters
{
std::string out_dir;
std::vector<double> timesteps;
std::vector<double> species_epsilon;
Rcpp::List init_params;
@ -52,13 +52,15 @@ struct RuntimeParameters {
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;
bool rollback_enabled = false;
bool control_interval_enabled = false;
std::uint32_t global_iter = 0;
std::uint32_t sur_disabled_counter = 0;
std::uint32_t rollback_counter = 0;
std::uint32_t checkpoint_interval = 0;
std::uint32_t control_interval = 0;
std::vector<double> mape_threshold;
std::vector<double> rrmse_threshold;
static constexpr std::uint32_t WORK_PACKAGE_SIZE_DEFAULT = 32;
std::uint32_t work_package_size = WORK_PACKAGE_SIZE_DEFAULT;