rollback implemented, triggered when MAPE exceeds epsilon

This commit is contained in:
rastogi 2025-09-12 09:47:00 +02:00
parent 733a13d9d9
commit 882dff4873
14 changed files with 1590 additions and 158 deletions

BIN
bin/checkpoint10.hdf5 Normal file

Binary file not shown.

BIN
bin/checkpoint20.hdf5 Normal file

Binary file not shown.

BIN
bin/dolo/iter_00.qs2 Normal file

Binary file not shown.

BIN
bin/dolo/iter_050.qs2 Normal file

Binary file not shown.

48
bin/dolo_fgcs.pqi Normal file
View File

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

131
bin/dolo_fgcs_3.R Normal file
View File

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

BIN
bin/dolo_fgcs_3.qs2 Normal file

Binary file not shown.

View File

@ -1,7 +1,7 @@
iterations <- 100 iterations <- 100
dt <- 200 dt <- 200
control_iteration <- 10 control_iteration <- 10
species_epsilon <- seq(0.001, 14) species_epsilon <- rep(1e-6, 14)
out_save <- seq(50, iterations, by = 50) out_save <- seq(50, iterations, by = 50)
list( list(

1307
bin/phreeqc_kin.dat Normal file

File diff suppressed because it is too large Load Diff

BIN
bin/poet

Binary file not shown.

View File

@ -1,151 +1,76 @@
Iteration, Species, MAPE, RRSME Iteration, Species, MAPE, RRSME
0, ID, 6.95625, 0.567065 10, ID, 0, 0
0, H, 10.7506, 0.321891 10, H, 2.18792e-14, 3.11625e-16
0, O, 10.6833, 0.320578 10, O, 0.0261378, 0.000268639
0, Charge, 6.72019e+10, 1.98915e+10 10, Charge, 0.0242586, 0.000249429
0, C, 1.82297e+07, 5.35302e+06 10, C, 4.72509e-05, 1.01573e-05
0, Ca, 379.805, 109.274 10, Ca, 4.7361e-05, 1.01573e-05
0, Cl, 1.30731e+06, 3.12249e+06 10, Cl, 13.2935, 0.180335
0, Mg, 1.30727e+06, 3.12241e+06 10, Mg, 14.6907, 0.277312
0, Calcite, 11.8903, 0.34551 10, Calcite, 2.81987e-05, 6.07374e-06
0, Calcite_p1, 5.27625, 0.229701 10, Calcite_p1, 0, 0
0, Dolomite, 100, 1 10, Dolomite, 0, 0
0, Dolomite_p1, 5.27625, 0.229701 10, Dolomite_p1, 0, 0
0, O2g_eq, 5.38873, 0.229706 10, O2g_eq, 0.0643273, 0.000661513
0, O2g_si, 5.27625, 0.229701 10, O2g_si, 0, 0
0, ID, 8.08, 0.405956 20, ID, 0, 0
0, H, 13.5939, 0.363642 20, H, 2.19678e-14, 3.29442e-16
0, O, 13.5705, 0.363318 20, O, 0.0280859, 0.000292534
0, Charge, 4.02457e+10, 2.1846e+10 20, Charge, 0.0233668, 0.000244248
0, C, 4.97993e+06, 2.56055e+06 20, C, 0.000537105, 6.41251e-05
0, Ca, 118.297, 56.8558 20, Ca, 0.000537197, 6.41227e-05
0, Cl, 25377, 23205.5 20, Cl, 7.24733, 0.177101
0, Mg, 25377, 23205.5 20, Mg, 7.43845, 0.196059
0, Calcite, 17.3444, 0.454722 20, Calcite, 0.000323255, 3.87549e-05
0, Calcite_p1, 7.52, 0.274226 20, Calcite_p1, 0, 0
0, Dolomite, 100, 1 20, Dolomite, 4.16664e-07, 6.48128e-07
0, Dolomite_p1, 7.52, 0.274226 20, Dolomite_p1, 0, 0
0, O2g_eq, 7.72812, 0.274241 20, O2g_eq, 0.0699382, 0.000728667
0, O2g_si, 7.52, 0.274226 20, O2g_si, 0, 0
0, ID, 11.1113, 0.407722 20, ID, 0, 0
0, H, 15.9914, 0.396067 20, H, 2.17942e-14, 3.78672e-16
0, O, 15.977, 0.395888 20, O, 0.0272958, 0.00028732
0, Charge, 3.80159e+10, 2.64461e+10 20, Charge, 0.0224743, 0.000238903
0, C, 2.64497e+06, 1.6471e+06 20, C, 0.000693403, 6.89355e-05
0, Ca, 80.799, 45.1004 20, Ca, 0.000693447, 6.89274e-05
0, Cl, 54.2385, 19.2069 20, Cl, 4.91228, 0.145942
0, Mg, 54.2385, 19.2069 20, Mg, 5.07712, 0.164729
0, Calcite, 22.1764, 0.629237 20, Calcite, 0.00041918, 4.19111e-05
0, Calcite_p1, 10.7438, 0.327777 20, Calcite_p1, 0, 0
0, Dolomite, 100, 1 20, Dolomite, 0.00175364, 0.00252754
0, Dolomite_p1, 10.7438, 0.327777 20, Dolomite_p1, 0, 0
0, O2g_eq, 11.0267, 0.327812 20, O2g_eq, 0.0678277, 0.000714771
0, O2g_si, 10.7438, 0.327777 20, O2g_si, 0, 0
0, ID, 13.885, 0.412735 20, ID, 0, 0
0, H, 19.47, 0.438484 20, H, 2.56842e-14, 4.21869e-16
0, O, 19.4621, 0.438327 20, O, 0.026514, 0.000282294
0, Charge, 2.43063e+10, 1.9873e+10 20, Charge, 0.0216724, 0.000234381
0, C, 1.30867e+06, 1.06974e+06 20, C, 0.000801173, 7.00726e-05
0, Ca, 50.2849, 25.2101 20, Ca, 0.000801423, 7.00996e-05
0, Cl, 31.7636, 4.11168 20, Cl, 3.63922, 0.123201
0, Mg, 31.7636, 4.11168 20, Mg, 3.86128, 0.144882
0, Calcite, 27.974, 0.877967 20, Calcite, 0.000487895, 4.31642e-05
0, Calcite_p1, 13.675, 0.369797 20, Calcite_p1, 0, 0
0, Dolomite, 100, 1 20, Dolomite, 0.00168115, 0.00101652
0, Dolomite_p1, 13.675, 0.369797 20, Dolomite_p1, 0, 0
0, O2g_eq, 14.1733, 0.369879 20, O2g_eq, 0.0657776, 0.000701359
0, O2g_si, 13.675, 0.369797 20, O2g_si, 0, 0
0, ID, 17.12, 0.467386 20, ID, 0, 0
0, H, 21.5302, 0.461047 20, H, 2.39005e-14, 4.35959e-16
0, O, 21.5263, 0.460881 20, O, 0.0257638, 0.000277486
0, Charge, 3.86401e+10, 2.5938e+10 20, Charge, 0.0209118, 0.000230386
0, C, 1.78992e+06, 1.19609e+06 20, C, 0.000836776, 6.68445e-05
0, Ca, 60.5235, 26.1934 20, Ca, 0.000836791, 6.68233e-05
0, Cl, 25.5508, 1.64324 20, Cl, 2.7374, 0.103349
0, Mg, 25.5508, 1.64324 20, Mg, 2.94845, 0.129046
0, Calcite, 33.1836, 1.31935 20, Calcite, 0.000508786, 4.10914e-05
0, Calcite_p1, 16.805, 0.409939 20, Calcite_p1, 0, 0
0, Dolomite, 100, 1 20, Dolomite, 0.000173459, 9.94268e-05
0, Dolomite_p1, 16.805, 0.409939 20, Dolomite_p1, 0, 0
0, O2g_eq, 17.914, 0.410276 20, O2g_eq, 0.06381, 0.000688484
0, O2g_si, 16.805, 0.409939 20, O2g_si, 0, 0
0, ID, 20.2338, 0.542414
0, H, 24.6071, 0.498266
0, O, 24.5965, 0.498068
0, Charge, 7.30486e+10, 4.29408e+10
0, C, 2.63576e+06, 1.28868e+06
0, Ca, 85.8119, 36.0012
0, Cl, 28.2702, 1.3898
0, Mg, 28.2702, 1.3898
0, Calcite, 37.9446, 1.65495
0, Calcite_p1, 19.6212, 0.442959
0, Dolomite, 100, 1
0, Dolomite_p1, 19.6212, 0.442959
0, O2g_eq, 20.4653, 0.443141
0, O2g_si, 19.6212, 0.442959
0, ID, 22.64, 0.512981
0, H, 26.5835, 0.544559
0, O, 26.5937, 0.544477
0, Charge, 2.32884e+10, 2.30432e+10
0, C, 563829, 497083
0, Ca, 43.2574, 16.53
0, Cl, 30.0228, 1.05626
0, Mg, 30.0228, 1.05626
0, Calcite, 52.2208, 2.91508
0, Calcite_p1, 22.395, 0.473234
0, Dolomite, 100, 1
0, Dolomite_p1, 22.395, 0.473234
0, O2g_eq, 24.9942, 0.474601
0, O2g_si, 22.395, 0.473234
0, ID, 25.295, 0.555045
0, H, 28.454, 0.539834
0, O, 28.4584, 0.539817
0, Charge, 6.87698e+10, 6.32776e+10
0, C, 352658, 313612
0, Ca, 71.4835, 39.6263
0, Cl, 31.4124, 2.41365
0, Mg, 31.4124, 2.41365
0, Calcite, 50.1774, 2.94726
0, Calcite_p1, 24.9275, 0.499274
0, Dolomite, 100, 1
0, Dolomite_p1, 24.9275, 0.499274
0, O2g_eq, 26.0254, 0.499577
0, O2g_si, 24.9275, 0.499274
0, ID, 27.6025, 0.588982
0, H, 30.3794, 0.594796
0, O, 30.3879, 0.594742
0, Charge, 1.03453e+11, 7.39959e+10
0, C, 719488, 1.0116e+06
0, Ca, 87.8521, 41.1368
0, Cl, 32.5071, 1.04708
0, Mg, 32.5071, 1.04708
0, Calcite, 58.4251, 3.98278
0, Calcite_p1, 27.13, 0.520865
0, Dolomite, 100, 1
0, Dolomite_p1, 27.13, 0.520865
0, O2g_eq, 30.637, 0.523121
0, O2g_si, 27.13, 0.520865
0, ID, 32.27, 0.657989
0, H, 42.1, 0.66122
0, O, 42.0933, 0.661169
0, Charge, 1.78607e+11, 1.01683e+11
0, C, 1.21472e+06, 1.5525e+06
0, Ca, 131.324, 50.8261
0, Cl, 46.4231, 0.949978
0, Mg, 46.4231, 0.949977
0, Calcite, 81.7668, 5.63658
0, Calcite_p1, 31.535, 0.56156
0, Dolomite, 100, 1
0, Dolomite_p1, 31.535, 0.56156
0, O2g_eq, 33.2615, 0.562111
0, O2g_si, 31.535, 0.56156

Binary file not shown.

View File

@ -300,6 +300,35 @@ void call_master_iter_end(RInside &R, const Field &trans, const Field &chem)
*global_rt_setup = R["setup"]; *global_rt_setup = R["setup"];
} }
bool checkAndRollback(ChemistryModule &chem, RuntimeParameters &params, uint32_t &iter)
{
for (uint32_t i = 0; i < chem.error_stats_history.size(); i++)
{
if (iter == chem.error_stats_history[i].iteration)
{
for (uint32_t j = 0; j < params.species_epsilon.size(); j++)
{
if (params.species_epsilon[j] < chem.error_stats_history[i].mape[j] && chem.error_stats_history[i].mape[j] != 0 && chem.control_iteration_counter > 1)
{
uint32_t rollback_iter = iter - params.control_iteration;
std::cout << chem.getField().GetProps()[j] << " with a MAPE value of " << chem.error_stats_history[i].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;
chem.control_iteration_counter--;
return true;
}
}
}
}
return false;
}
static Rcpp::List RunMasterLoop(RInsidePOET &R, RuntimeParameters &params, static Rcpp::List RunMasterLoop(RInsidePOET &R, RuntimeParameters &params,
DiffusionModule &diffusion, DiffusionModule &diffusion,
ChemistryModule &chem) ChemistryModule &chem)
@ -441,21 +470,15 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, RuntimeParameters &params,
MSG("End of *coupling* iteration " + std::to_string(iter) + "/" + MSG("End of *coupling* iteration " + std::to_string(iter) + "/" +
std::to_string(maxiter)); std::to_string(maxiter));
/* if (iter % params.control_iteration == 0)
if (params.control_iteration_active)
{ {
std::string file_path = "checkpoint" + std::to_string(iter) + ".hdf5"; writeStatsToCSV(chem.error_stats_history, chem.getField().GetProps(), "stats_overview");
write_checkpoint(file_path,
{.field = chem.getField(), .iteration = iter});
}
if (iter % params.control_iteration == 0 && iter != 0)
{
write_checkpoint("checkpoint" + std::to_string(iter) + ".hdf5", write_checkpoint("checkpoint" + std::to_string(iter) + ".hdf5",
{.field = chem.getField(), .iteration = iter}); {.field = chem.getField(), .iteration = iter});
checkAndRollback(chem, params, iter);
} }
*/
// MSG(); // MSG();
} // END SIMULATION LOOP } // END SIMULATION LOOP
@ -503,8 +526,6 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, RuntimeParameters &params,
chem.MasterLoopBreak(); chem.MasterLoopBreak();
writeStatsToCSV(chem.error_stats_history, chem.getField().GetProps(), "stats_overview");
return profiling; return profiling;
} }