diff --git a/bin/dolo_fgcs_3.R b/bin/dolo_fgcs_3.R new file mode 100644 index 000000000..eb57e65bb --- /dev/null +++ b/bin/dolo_fgcs_3.R @@ -0,0 +1,135 @@ +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 <- 15000 +dt <- 200 +checkpoint_interval <- 100 +control_interval <- 100 +mape_threshold <- rep(0.1, 13) +mape_threshold[5] <- 1 #Charge +out_save <- seq(1000, iterations, by = 1000) +#out_save = c(seq(1, 10), seq(10, 100, by= 10), seq(200, iterations, by=100)) + + +list( + timesteps = rep(dt, iterations), + store_result = TRUE, + out_save = out_save, + checkpoint_interval = checkpoint_interval, + control_interval = control_interval, + mape_threshold = mape_threshold +) \ No newline at end of file diff --git a/bin/dolo_fgcs_3.qs2 b/bin/dolo_fgcs_3.qs2 new file mode 100644 index 000000000..23ed75af5 Binary files /dev/null and b/bin/dolo_fgcs_3.qs2 differ diff --git a/bin/phreeqc_kin.dat b/bin/phreeqc_kin.dat new file mode 100644 index 000000000..01d5bf516 --- /dev/null +++ b/bin/phreeqc_kin.dat @@ -0,0 +1,1307 @@ +### This is the standard "phreeqc.dat" stripped of EXCHANGE and +### SURFACE and with the RATES for Calcite and Dolomite to use with +### RedModRphree + +### Time-stamp: "Last modified 2023-05-23 10:35:56 mluebke" + +# PHREEQC.DAT for calculating pressure dependence of reactions, with +# molal volumina of aqueous species and of minerals, and +# critical temperatures and pressures of gases used in Peng-Robinson's EOS. +# Details are given at the end of this file. + +SOLUTION_MASTER_SPECIES +# +#element species alk gfw_formula element_gfw +# +H H+ -1.0 H 1.008 +H(0) H2 0 H +H(1) H+ -1.0 0 +E e- 0 0.0 0 +O H2O 0 O 16.0 +O(0) O2 0 O +O(-2) H2O 0 0 +Ca Ca+2 0 Ca 40.08 +Mg Mg+2 0 Mg 24.312 +Na Na+ 0 Na 22.9898 +K K+ 0 K 39.102 +Fe Fe+2 0 Fe 55.847 +Fe(+2) Fe+2 0 Fe +Fe(+3) Fe+3 -2.0 Fe +Mn Mn+2 0 Mn 54.938 +Mn(+2) Mn+2 0 Mn +Mn(+3) Mn+3 0 Mn +Al Al+3 0 Al 26.9815 +Ba Ba+2 0 Ba 137.34 +Sr Sr+2 0 Sr 87.62 +Si H4SiO4 0 SiO2 28.0843 +Cl Cl- 0 Cl 35.453 +C CO3-2 2.0 HCO3 12.0111 +C(+4) CO3-2 2.0 HCO3 +C(-4) CH4 0 CH4 +Alkalinity CO3-2 1.0 Ca0.5(CO3)0.5 50.05 +S SO4-2 0 SO4 32.064 +S(6) SO4-2 0 SO4 +S(-2) HS- 1.0 S +N NO3- 0 N 14.0067 +N(+5) NO3- 0 N +N(+3) NO2- 0 N +N(0) N2 0 N +N(-3) NH4+ 0 N 14.0067 +#Amm AmmH+ 0 AmmH 17.0 +B H3BO3 0 B 10.81 +P PO4-3 2.0 P 30.9738 +F F- 0 F 18.9984 +Li Li+ 0 Li 6.939 +Br Br- 0 Br 79.904 +Zn Zn+2 0 Zn 65.37 +Cd Cd+2 0 Cd 112.4 +Pb Pb+2 0 Pb 207.19 +Cu Cu+2 0 Cu 63.546 +Cu(+2) Cu+2 0 Cu +Cu(+1) Cu+1 0 Cu +# redox-uncoupled gases +Hdg Hdg 0 Hdg 2.016 # H2 gas +Oxg Oxg 0 Oxg 32 # O2 gas +Mtg Mtg 0 Mtg 16.032 # CH4 gas +Sg H2Sg 1.0 H2Sg 34.08 +Ntg Ntg 0 Ntg 28.0134 # N2 gas + +SOLUTION_SPECIES +H+ = H+ + -gamma 9.0 0 + -dw 9.31e-9 +e- = e- +H2O = H2O +Ca+2 = Ca+2 + -gamma 5.0 0.1650 + -dw 0.793e-9 + -Vm -0.3456 -7.252 6.149 -2.479 1.239 5 1.60 -57.1 -6.12e-3 1 # ref. 1 +Mg+2 = Mg+2 + -gamma 5.5 0.20 + -dw 0.705e-9 + -Vm -1.410 -8.6 11.13 -2.39 1.332 5.5 1.29 -32.9 -5.86e-3 1 # ref. 1 +Na+ = Na+ + -gamma 4.0 0.075 + -gamma 4.08 0.082 # halite solubility + -dw 1.33e-9 + -Vm 2.28 -4.38 -4.1 -0.586 0.09 4 0.3 52 -3.33e-3 0.566 # ref. 1 +# for calculating densities (rho) when I > 3... + # -Vm 2.28 -4.38 -4.1 -0.586 0.09 4 0.3 52 -3.33e-3 0.45 +K+ = K+ + -gamma 3.5 0.015 + -dw 1.96e-9 + -Vm 3.322 -1.473 6.534 -2.712 9.06e-2 3.5 0 29.7 0 1 # ref. 1 +Fe+2 = Fe+2 + -gamma 6.0 0 + -dw 0.719e-9 + -Vm -0.3255 -9.687 1.536 -2.379 0.3033 6 -4.21e-2 39.7 0 1 # ref. 1 +Mn+2 = Mn+2 + -gamma 6.0 0 + -dw 0.688e-9 + -Vm -1.10 -8.03 4.08 -2.45 1.4 6 8.07 0 -1.51e-2 0.118 # ref. 2 +Al+3 = Al+3 + -gamma 9.0 0 + -dw 0.559e-9 + -Vm -2.28 -17.1 10.9 -2.07 2.87 9 0 0 5.5e-3 1 # ref. 2 and Barta and Hepler, 1986, Can. J.C. 64, 353. +Ba+2 = Ba+2 + -gamma 5.0 0 + -gamma 4.0 0.153 # Barite solubility + -dw 0.848e-9 + -Vm 2.063 -10.06 1.9534 -2.36 0.4218 5 1.58 -12.03 -8.35e-3 1 # ref. 1 +Sr+2 = Sr+2 + -gamma 5.260 0.121 + -dw 0.794e-9 + -Vm -1.57e-2 -10.15 10.18 -2.36 0.860 5.26 0.859 -27.0 -4.1e-3 1.97 # ref. 1 +H4SiO4 = H4SiO4 + -dw 1.10e-9 + -Vm 10.5 1.7 20 -2.7 0.1291 # supcrt + 2*H2O in a1 +Cl- = Cl- + -gamma 3.5 0.015 + -gamma 3.63 0.017 # cf. pitzer.dat + -dw 2.03e-9 + -Vm 4.465 4.801 4.325 -2.847 1.748 0 -0.331 20.16 0 1 # ref. 1 +CO3-2 = CO3-2 + -gamma 5.4 0 + -dw 0.955e-9 + -Vm 5.95 0 0 -5.67 6.85 0 1.37 106 -0.0343 1 # ref. 1 +SO4-2 = SO4-2 + -gamma 5.0 -0.04 + -dw 1.07e-9 + -Vm 8.0 2.3 -46.04 6.245 3.82 0 0 0 0 1 # ref. 1 +NO3- = NO3- + -gamma 3.0 0 + -dw 1.9e-9 + -Vm 6.32 6.78 0 -3.06 0.346 0 0.93 0 -0.012 1 # ref. 1 +#AmmH+ = AmmH+ +# -gamma 2.5 0 +# -dw 1.98e-9 +# -Vm 4.837 2.345 5.522 -2.88 1.096 3 -1.456 75.0 7.17e-3 1 # ref. 1 +H3BO3 = H3BO3 + -dw 1.1e-9 + -Vm 7.0643 8.8547 3.5844 -3.1451 -.2000 # supcrt +PO4-3 = PO4-3 + -gamma 4.0 0 + -dw 0.612e-9 + -Vm 1.24 -9.07 9.31 -2.4 5.61 0 0 0 -1.41e-2 1 # ref. 2 +F- = F- + -gamma 3.5 0 + -dw 1.46e-9 + -Vm 0.928 1.36 6.27 -2.84 1.84 0 0 -0.318 0 1 # ref. 2 +Li+ = Li+ + -gamma 6.0 0 + -dw 1.03e-9 + -Vm -0.419 -0.069 13.16 -2.78 0.416 0 0.296 -12.4 -2.74e-3 1.26 # ref. 2 and Ellis, 1968, J. Chem. Soc. A, 1138 +Br- = Br- + -gamma 3.0 0 + -dw 2.01e-9 + -Vm 6.72 2.85 4.21 -3.14 1.38 0 -9.56e-2 7.08 -1.56e-3 1 # ref. 2 +Zn+2 = Zn+2 + -gamma 5.0 0 + -dw 0.715e-9 + -Vm -1.96 -10.4 14.3 -2.35 1.46 5 -1.43 24 1.67e-2 1.11 # ref. 2 +Cd+2 = Cd+2 + -dw 0.717e-9 + -Vm 1.63 -10.7 1.01 -2.34 1.47 5 0 0 0 1 # ref. 2 +Pb+2 = Pb+2 + -dw 0.945e-9 + -Vm -.0051 -7.7939 8.8134 -2.4568 1.0788 4.5 # supcrt +Cu+2 = Cu+2 + -gamma 6.0 0 + -dw 0.733e-9 + -Vm -1.13 -10.5 7.29 -2.35 1.61 6 9.78e-2 0 3.42e-3 1 # ref. 2 +# redox-uncoupled gases +Hdg = Hdg # H2 + -dw 5.13e-9 + -Vm 6.52 0.78 0.12 # supcrt +Oxg = Oxg # O2 + -dw 2.35e-9 + -Vm 5.7889 6.3536 3.2528 -3.0417 -0.3943 # supcrt +Mtg = Mtg # CH4 + -dw 1.85e-9 + -Vm 7.7 # CH4 solubility, 25-100C, 1-700atm +Ntg = Ntg # N2 + -dw 1.96e-9 + -Vm 7 # Pray et al., 1952, IEC 44. 1146 +H2Sg = H2Sg # H2S + -dw 2.1e-9 + -Vm 7.81 2.96 -0.46 # supcrt +# aqueous species +H2O = OH- + H+ + -analytic 293.29227 0.1360833 -10576.913 -123.73158 0 -6.996455e-5 + -gamma 3.5 0 + -dw 5.27e-9 + -Vm -9.66 28.5 80.0 -22.9 1.89 0 1.09 0 0 1 # ref. 1 +2 H2O = O2 + 4 H+ + 4 e- + -log_k -86.08 + -delta_h 134.79 kcal + -dw 2.35e-9 + -Vm 5.7889 6.3536 3.2528 -3.0417 -0.3943 # supcrt +2 H+ + 2 e- = H2 + -log_k -3.15 + -delta_h -1.759 kcal + -dw 5.13e-9 + -Vm 6.52 0.78 0.12 # supcrt +CO3-2 + H+ = HCO3- + -log_k 10.329 + -delta_h -3.561 kcal + -analytic 107.8871 0.03252849 -5151.79 -38.92561 563713.9 + -gamma 5.4 0 + -dw 1.18e-9 + -Vm 8.472 0 -11.5 0 1.56 0 0 146 3.16e-3 1 # ref. 1 +CO3-2 + 2 H+ = CO2 + H2O + -log_k 16.681 + -delta_h -5.738 kcal + -analytic 464.1965 0.09344813 -26986.16 -165.75951 2248628.9 + -dw 1.92e-9 + -Vm 20.85 -46.93 -79.0 27.9 -0.193 # ref. 1 +CO3-2 + 10 H+ + 8 e- = CH4 + 3 H2O + -log_k 41.071 + -delta_h -61.039 kcal + -dw 1.85e-9 + -Vm 7.7 +SO4-2 + H+ = HSO4- + -log_k 1.988 + -delta_h 3.85 kcal + -analytic -56.889 0.006473 2307.9 19.8858 + -dw 1.33e-9 + -Vm 8.2 9.2590 2.1108 -3.1618 1.1748 0 -0.3 15 0 1 # ref. 1 +HS- = S-2 + H+ + -log_k -12.918 + -delta_h 12.1 kcal + -gamma 5.0 0 + -dw 0.731e-9 +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-9 + -Vm 5.0119 4.9799 3.4765 -2.9849 1.4410 # supcrt +HS- + H+ = H2S + -log_k 6.994 + -delta_h -5.30 kcal + -analytical -11.17 0.02386 3279.0 + -dw 2.1e-9 + -Vm 7.81 2.96 -0.46 # supcrt +H2Sg = HSg- + H+ + -log_k -6.994 + -delta_h 5.30 kcal + -analytical 11.17 -0.02386 -3279.0 + -dw 2.1e-9 + -Vm 5.0119 4.9799 3.4765 -2.9849 1.4410 # supcrt +NO3- + 2 H+ + 2 e- = NO2- + H2O + -log_k 28.570 + -delta_h -43.760 kcal + -gamma 3.0 0 + -dw 1.91e-9 + -Vm 5.5864 5.8590 3.4472 -3.0212 1.1847 # supcrt +2 NO3- + 12 H+ + 10 e- = N2 + 6 H2O + -log_k 207.08 + -delta_h -312.130 kcal + -dw 1.96e-9 + -Vm 7 # Pray et al., 1952, IEC 44. 1146 +NO3- + 10 H+ + 8 e- = NH4+ + 3 H2O + -log_k 119.077 + -delta_h -187.055 kcal + -gamma 2.5 0 + -dw 1.98e-9 + -Vm 4.837 2.345 5.522 -2.88 1.096 3 -1.456 75.0 7.17e-3 1 # ref. 1 +#AmmH+ = AmmH+ +# -gamma 2.5 0.0 +# -dw 1.98e-9 +# -Vm 4.837 2.345 5.522 -2.88 1.096 3 -1.456 75.0 7.17e-3 1 # supcrt modified +NH4+ = NH3 + H+ + -log_k -9.252 + -delta_h 12.48 kcal + -analytic 0.6322 -0.001225 -2835.76 + -dw 2.28e-9 + -Vm 6.69 2.8 3.58 -2.88 1.43 # ref. 2 +#AmmH+ = Amm + H+ +# -log_k -9.252 +# -delta_h 12.48 kcal +# -analytic 0.6322 -0.001225 -2835.76 +# -dw 2.28e-9 +# -Vm 6.69 2.8 3.58 -2.88 1.43 # ref. 2 +NH4+ + SO4-2 = NH4SO4- + log_k 1.11 + -Vm 14.0 0 -35.2 0 0 0 12.3 0 -0.141 1 # ref. 2 +#AmmH+ + SO4-2 = AmmHSO4- +# -log_k 1.11 +# -Vm 14.0 0 -35.2 0 0 0 12.3 0 -0.141 1 # ref. 2 +H3BO3 = H2BO3- + H+ + -log_k -9.24 + -delta_h 3.224 kcal +H3BO3 + F- = BF(OH)3- + -log_k -0.4 + -delta_h 1.850 kcal +H3BO3 + 2 F- + H+ = BF2(OH)2- + H2O + -log_k 7.63 + -delta_h 1.618 kcal +H3BO3 + 2 H+ + 3 F- = BF3OH- + 2 H2O + -log_k 13.67 + -delta_h -1.614 kcal +H3BO3 + 3 H+ + 4 F- = BF4- + 3 H2O + -log_k 20.28 + -delta_h -1.846 kcal +PO4-3 + H+ = HPO4-2 + -log_k 12.346 + -delta_h -3.530 kcal + -gamma 5.0 0 + -dw 0.69e-9 + -Vm 3.52 1.09 8.39 -2.82 3.34 0 0 0 0 1 # ref. 2 +PO4-3 + 2 H+ = H2PO4- + -log_k 19.553 + -delta_h -4.520 kcal + -gamma 5.4 0 + -dw 0.846e-9 + -Vm 5.58 8.06 12.2 -3.11 1.3 0 0 0 1.62e-2 1 # ref. 2 +PO4-3 + 3H+ = H3PO4 + log_k 21.721 # log_k and delta_h from minteq.v4.dat, NIST46.3 + delta_h -10.1 kJ + -Vm 7.47 12.4 6.29 -3.29 0 # ref. 2 +H+ + F- = HF + -log_k 3.18 + -delta_h 3.18 kcal + -analytic -2.033 0.012645 429.01 + -Vm 3.4753 .7042 5.4732 -2.8081 -.0007 # supcrt +H+ + 2 F- = HF2- + -log_k 3.76 + -delta_h 4.550 kcal + -Vm 5.2263 4.9797 3.7928 -2.9849 1.2934 # supcrt +Ca+2 + H2O = CaOH+ + H+ + -log_k -12.78 +Ca+2 + CO3-2 = CaCO3 + -log_k 3.224 + -delta_h 3.545 kcal + -analytic -1228.732 -0.299440 35512.75 485.818 + -dw 4.46e-10 # complexes: calc'd with the Pikal formula + -Vm -.2430 -8.3748 9.0417 -2.4328 -.0300 # supcrt +Ca+2 + CO3-2 + H+ = CaHCO3+ + -log_k 11.435 + -delta_h -0.871 kcal + -analytic 1317.0071 0.34546894 -39916.84 -517.70761 563713.9 + -gamma 6.0 0 + -dw 5.06e-10 + -Vm 3.1911 .0104 5.7459 -2.7794 .3084 5.4 # supcrt +Ca+2 + SO4-2 = CaSO4 + -log_k 2.25 + -delta_h 1.325 kcal + -dw 4.71e-10 + -Vm 2.7910 -.9666 6.1300 -2.7390 -.0010 # supcrt +Ca+2 + HSO4- = CaHSO4+ + -log_k 1.08 +Ca+2 + PO4-3 = CaPO4- + -log_k 6.459 + -delta_h 3.10 kcal + -gamma 5.4 0.0 +Ca+2 + HPO4-2 = CaHPO4 + -log_k 2.739 + -delta_h 3.3 kcal +Ca+2 + H2PO4- = CaH2PO4+ + -log_k 1.408 + -delta_h 3.4 kcal + -gamma 5.4 0.0 +# Ca+2 + F- = CaF+ + # -log_k 0.94 + # -delta_h 4.120 kcal + # -gamma 5.5 0.0 + # -Vm .9846 -5.3773 7.8635 -2.5567 .6911 5.5 # supcrt +Mg+2 + H2O = MgOH+ + H+ + -log_k -11.44 + -delta_h 15.952 kcal + -gamma 6.5 0 +Mg+2 + CO3-2 = MgCO3 + -log_k 2.98 + -delta_h 2.713 kcal + -analytic 0.9910 0.00667 + -dw 4.21e-10 + -Vm -.5837 -9.2067 9.3687 -2.3984 -.0300 # supcrt +Mg+2 + H+ + CO3-2 = MgHCO3+ + -log_k 11.399 + -delta_h -2.771 kcal + -analytic 48.6721 0.03252849 -2614.335 -18.00263 563713.9 + -gamma 4.0 0 + -dw 4.78e-10 + -Vm 2.7171 -1.1469 6.2008 -2.7316 .5985 4 # supcrt +Mg+2 + SO4-2 = MgSO4 + -log_k 2.37 + -delta_h 4.550 kcal + -dw 4.45e-10 + -Vm 2.4 -0.97 6.1 -2.74 # est'd +Mg+2 + PO4-3 = MgPO4- + -log_k 6.589 + -delta_h 3.10 kcal + -gamma 5.4 0 +Mg+2 + HPO4-2 = MgHPO4 + -log_k 2.87 + -delta_h 3.3 kcal +Mg+2 + H2PO4- = MgH2PO4+ + -log_k 1.513 + -delta_h 3.4 kcal + -gamma 5.4 0 +Mg+2 + F- = MgF+ + -log_k 1.82 + -delta_h 3.20 kcal + -gamma 4.5 0 + -Vm .6494 -6.1958 8.1852 -2.5229 .9706 4.5 # supcrt +Na+ + OH- = NaOH + -log_k -10 # remove this complex +Na+ + CO3-2 = NaCO3- + -log_k 1.27 + -delta_h 8.91 kcal + -dw 5.85e-10 + -Vm 3.89 -8.23e-4 20 -9.44 3.02 9.05e-3 3.07 0 0.0233 1 # ref. 1 +Na+ + HCO3- = NaHCO3 + -log_k -0.25 + -delta_h -1 kcal + -dw 6.73e-10 + -Vm 0.431 # ref. 1 +Na+ + SO4-2 = NaSO4- + -log_k 0.7 + -delta_h 1.120 kcal + -gamma 5.4 0 + -dw 6.18e-10 + -Vm 1e-5 16.4 -0.0678 -1.05 4.14 0 6.86 0 0.0242 0.53 # ref. 1 +Na+ + HPO4-2 = NaHPO4- + -log_k 0.29 + -gamma 5.4 0 + -Vm 5.2 8.1 13 -3 0.9 0 0 1.62e-2 1 # ref. 2 +Na+ + F- = NaF + -log_k -0.24 + -Vm 2.7483 -1.0708 6.1709 -2.7347 -.030 # supcrt +K+ + SO4-2 = KSO4- + -log_k 0.85 + -delta_h 2.250 kcal + -analytical 3.106 0.0 -673.6 + -gamma 5.4 0 + -dw 7.46e-10 + -Vm 6.8 7.06 3.0 -2.07 1.1 0 0 0 0 1 # ref. 1 +K+ + HPO4-2 = KHPO4- + -log_k 0.29 + -gamma 5.4 0 + -Vm 5.4 8.1 19 -3.1 0.7 0 0 0 1.62e-2 1 # ref. 2 +Fe+2 + H2O = FeOH+ + H+ + -log_k -9.5 + -delta_h 13.20 kcal + -gamma 5.0 0 +Fe+2 + 3H2O = Fe(OH)3- + 3H+ + -log_k -31.0 + -delta_h 30.3 kcal + -gamma 5.0 0 +Fe+2 + Cl- = FeCl+ + -log_k 0.14 +Fe+2 + CO3-2 = FeCO3 + -log_k 4.38 +Fe+2 + HCO3- = FeHCO3+ + -log_k 2.0 +Fe+2 + SO4-2 = FeSO4 + -log_k 2.25 + -delta_h 3.230 kcal + -Vm -13 0 123 # ref. 2 +Fe+2 + HSO4- = FeHSO4+ + -log_k 1.08 +Fe+2 + 2HS- = Fe(HS)2 + -log_k 8.95 +Fe+2 + 3HS- = Fe(HS)3- + -log_k 10.987 +Fe+2 + HPO4-2 = FeHPO4 + -log_k 3.6 +Fe+2 + H2PO4- = FeH2PO4+ + -log_k 2.7 + -gamma 5.4 0 +Fe+2 + F- = FeF+ + -log_k 1.0 +Fe+2 = Fe+3 + e- + -log_k -13.02 + -delta_h 9.680 kcal + -gamma 9.0 0 +Fe+3 + H2O = FeOH+2 + H+ + -log_k -2.19 + -delta_h 10.4 kcal + -gamma 5.0 0 +Fe+3 + 2 H2O = Fe(OH)2+ + 2 H+ + -log_k -5.67 + -delta_h 17.1 kcal + -gamma 5.4 0 +Fe+3 + 3 H2O = Fe(OH)3 + 3 H+ + -log_k -12.56 + -delta_h 24.8 kcal +Fe+3 + 4 H2O = Fe(OH)4- + 4 H+ + -log_k -21.6 + -delta_h 31.9 kcal + -gamma 5.4 0 +Fe+2 + 2H2O = Fe(OH)2 + 2H+ + -log_k -20.57 + -delta_h 28.565 kcal +2 Fe+3 + 2 H2O = Fe2(OH)2+4 + 2 H+ + -log_k -2.95 + -delta_h 13.5 kcal +3 Fe+3 + 4 H2O = Fe3(OH)4+5 + 4 H+ + -log_k -6.3 + -delta_h 14.3 kcal +Fe+3 + Cl- = FeCl+2 + -log_k 1.48 + -delta_h 5.6 kcal + -gamma 5.0 0 +Fe+3 + 2 Cl- = FeCl2+ + -log_k 2.13 + -gamma 5.0 0 +Fe+3 + 3 Cl- = FeCl3 + -log_k 1.13 +Fe+3 + SO4-2 = FeSO4+ + -log_k 4.04 + -delta_h 3.91 kcal + -gamma 5.0 0 +Fe+3 + HSO4- = FeHSO4+2 + -log_k 2.48 +Fe+3 + 2 SO4-2 = Fe(SO4)2- + -log_k 5.38 + -delta_h 4.60 kcal +Fe+3 + HPO4-2 = FeHPO4+ + -log_k 5.43 + -delta_h 5.76 kcal + -gamma 5.0 0 +Fe+3 + H2PO4- = FeH2PO4+2 + -log_k 5.43 + -gamma 5.4 0 +Fe+3 + F- = FeF+2 + -log_k 6.2 + -delta_h 2.7 kcal + -gamma 5.0 0 +Fe+3 + 2 F- = FeF2+ + -log_k 10.8 + -delta_h 4.8 kcal + -gamma 5.0 0 +Fe+3 + 3 F- = FeF3 + -log_k 14.0 + -delta_h 5.4 kcal +Mn+2 + H2O = MnOH+ + H+ + -log_k -10.59 + -delta_h 14.40 kcal + -gamma 5.0 0 +Mn+2 + 3H2O = Mn(OH)3- + 3H+ + -log_k -34.8 + -gamma 5.0 0 +Mn+2 + Cl- = MnCl+ + -log_k 0.61 + -gamma 5.0 0 + -Vm 7.25 -1.08 -25.8 -2.73 3.99 5 0 0 0 1 # ref. 2 +Mn+2 + 2 Cl- = MnCl2 + -log_k 0.25 + -Vm 1e-5 0 144 # ref. 2 +Mn+2 + 3 Cl- = MnCl3- + -log_k -0.31 + -gamma 5.0 0 + -Vm 11.8 0 0 0 2.4 0 0 0 3.6e-2 1 # ref. 2 +Mn+2 + CO3-2 = MnCO3 + -log_k 4.9 +Mn+2 + HCO3- = MnHCO3+ + -log_k 1.95 + -gamma 5.0 0 +Mn+2 + SO4-2 = MnSO4 + -log_k 2.25 + -delta_h 3.370 kcal + -Vm -1.31 -1.83 62.3 -2.7 # ref. 2 +Mn+2 + 2 NO3- = Mn(NO3)2 + -log_k 0.6 + -delta_h -0.396 kcal + -Vm 6.16 0 29.4 0 0.9 # ref. 2 +Mn+2 + F- = MnF+ + -log_k 0.84 + -gamma 5.0 0 +Mn+2 = Mn+3 + e- + -log_k -25.51 + -delta_h 25.80 kcal + -gamma 9.0 0 +Al+3 + H2O = AlOH+2 + H+ + -log_k -5.0 + -delta_h 11.49 kcal + -analytic -38.253 0.0 -656.27 14.327 + -gamma 5.4 0 + -Vm -1.46 -11.4 10.2 -2.31 1.67 5.4 0 0 0 1 # ref. 2 and Barta and Hepler, 1986, Can. J. Chem. 64, 353. +Al+3 + 2 H2O = Al(OH)2+ + 2 H+ + -log_k -10.1 + -delta_h 26.90 kcal + -gamma 5.4 0 + -analytic 88.50 0.0 -9391.6 -27.121 +Al+3 + 3 H2O = Al(OH)3 + 3 H+ + -log_k -16.9 + -delta_h 39.89 kcal + -analytic 226.374 0.0 -18247.8 -73.597 +Al+3 + 4 H2O = Al(OH)4- + 4 H+ + -log_k -22.7 + -delta_h 42.30 kcal + -analytic 51.578 0.0 -11168.9 -14.865 + -gamma 4.5 0 +Al+3 + SO4-2 = AlSO4+ + -log_k 3.5 + -delta_h 2.29 kcal + -gamma 4.5 0 +Al+3 + 2SO4-2 = Al(SO4)2- + -log_k 5.0 + -delta_h 3.11 kcal + -gamma 4.5 0 +Al+3 + HSO4- = AlHSO4+2 + -log_k 0.46 +Al+3 + F- = AlF+2 + -log_k 7.0 + -delta_h 1.060 kcal + -gamma 5.4 0 +Al+3 + 2 F- = AlF2+ + -log_k 12.7 + -delta_h 1.980 kcal + -gamma 5.4 0 +Al+3 + 3 F- = AlF3 + -log_k 16.8 + -delta_h 2.160 kcal +Al+3 + 4 F- = AlF4- + -log_k 19.4 + -delta_h 2.20 kcal + -gamma 4.5 0 +# Al+3 + 5 F- = AlF5-2 + # log_k 20.6 + # delta_h 1.840 kcal +# Al+3 + 6 F- = AlF6-3 + # log_k 20.6 + # delta_h -1.670 kcal +H4SiO4 = H3SiO4- + H+ + -log_k -9.83 + -delta_h 6.12 kcal + -analytic -302.3724 -0.050698 15669.69 108.18466 -1119669.0 + -gamma 4 0 + -Vm 7.94 1.0881 5.3224 -2.8240 1.4767 # supcrt + H2O in a1 +H4SiO4 = H2SiO4-2 + 2 H+ + -log_k -23.0 + -delta_h 17.6 kcal + -analytic -294.0184 -0.072650 11204.49 108.18466 -1119669.0 + -gamma 5.4 0 +H4SiO4 + 4 H+ + 6 F- = SiF6-2 + 4 H2O + -log_k 30.18 + -delta_h -16.260 kcal + -gamma 5.0 0 + -Vm 8.5311 13.0492 .6211 -3.3185 2.7716 # supcrt +Ba+2 + H2O = BaOH+ + H+ + -log_k -13.47 + -gamma 5.0 0 +Ba+2 + CO3-2 = BaCO3 + -log_k 2.71 + -delta_h 3.55 kcal + -analytic 0.113 0.008721 + -Vm .2907 -7.0717 8.5295 -2.4867 -.0300 # supcrt +Ba+2 + HCO3- = BaHCO3+ + -log_k 0.982 + -delta_h 5.56 kcal + -analytic -3.0938 0.013669 +Ba+2 + SO4-2 = BaSO4 + -log_k 2.7 +Sr+2 + H2O = SrOH+ + H+ + -log_k -13.29 + -gamma 5.0 0 +Sr+2 + CO3-2 + H+ = SrHCO3+ + -log_k 11.509 + -delta_h 2.489 kcal + -analytic 104.6391 0.04739549 -5151.79 -38.92561 563713.9 + -gamma 5.4 0 +Sr+2 + CO3-2 = SrCO3 + -log_k 2.81 + -delta_h 5.22 kcal + -analytic -1.019 0.012826 + -Vm -.1787 -8.2177 8.9799 -2.4393 -.0300 # supcrt +Sr+2 + SO4-2 = SrSO4 + -log_k 2.29 + -delta_h 2.08 kcal + -Vm 6.7910 -.9666 6.1300 -2.7390 -.0010 # celestite solubility +Li+ + SO4-2 = LiSO4- + -log_k 0.64 + -gamma 5.0 0 +Cu+2 + e- = Cu+ + -log_k 2.72 + -delta_h 1.65 kcal + -gamma 2.5 0 +Cu+ + 2Cl- = CuCl2- + -log_k 5.50 + -delta_h -0.42 kcal + -gamma 4.0 0 +Cu+ + 3Cl- = CuCl3-2 + -log_k 5.70 + -delta_h 0.26 kcal + -gamma 5.0 0.0 +Cu+2 + CO3-2 = CuCO3 + -log_k 6.73 +Cu+2 + 2CO3-2 = Cu(CO3)2-2 + -log_k 9.83 +Cu+2 + HCO3- = CuHCO3+ + -log_k 2.7 +Cu+2 + Cl- = CuCl+ + -log_k 0.43 + -delta_h 8.65 kcal + -gamma 4.0 0 + -Vm -4.19 0 30.4 0 0 4 0 0 1.94e-2 1 # ref. 2 +Cu+2 + 2Cl- = CuCl2 + -log_k 0.16 + -delta_h 10.56 kcal + -Vm 26.8 0 -136 # ref. 2 +Cu+2 + 3Cl- = CuCl3- + -log_k -2.29 + -delta_h 13.69 kcal + -gamma 4.0 0 +Cu+2 + 4Cl- = CuCl4-2 + -log_k -4.59 + -delta_h 17.78 kcal + -gamma 5.0 0 +Cu+2 + F- = CuF+ + -log_k 1.26 + -delta_h 1.62 kcal +Cu+2 + H2O = CuOH+ + H+ + -log_k -8.0 + -gamma 4.0 0 +Cu+2 + 2 H2O = Cu(OH)2 + 2 H+ + -log_k -13.68 +Cu+2 + 3 H2O = Cu(OH)3- + 3 H+ + -log_k -26.9 +Cu+2 + 4 H2O = Cu(OH)4-2 + 4 H+ + -log_k -39.6 +2Cu+2 + 2H2O = Cu2(OH)2+2 + 2H+ + -log_k -10.359 + -delta_h 17.539 kcal + -analytical 2.497 0.0 -3833.0 +Cu+2 + SO4-2 = CuSO4 + -log_k 2.31 + -delta_h 1.220 kcal + -Vm 5.21 0 -14.6 # ref. 2 +Cu+2 + 3HS- = Cu(HS)3- + -log_k 25.9 +Zn+2 + H2O = ZnOH+ + H+ + -log_k -8.96 + -delta_h 13.4 kcal +Zn+2 + 2 H2O = Zn(OH)2 + 2 H+ + -log_k -16.9 +Zn+2 + 3 H2O = Zn(OH)3- + 3 H+ + -log_k -28.4 +Zn+2 + 4 H2O = Zn(OH)4-2 + 4 H+ + -log_k -41.2 +Zn+2 + Cl- = ZnCl+ + -log_k 0.43 + -delta_h 7.79 kcal + -gamma 4.0 0 + -Vm 14.8 -3.91 -105.7 -2.62 0.203 4 0 0 -5.05e-2 1 # ref. 2 +Zn+2 + 2 Cl- = ZnCl2 + -log_k 0.45 + -delta_h 8.5 kcal + -Vm -10.1 4.57 241 -2.97 -1e-3 # ref. 2 +Zn+2 + 3Cl- = ZnCl3- + -log_k 0.5 + -delta_h 9.56 kcal + -gamma 4.0 0 + -Vm 0.772 15.5 -0.349 -3.42 1.25 0 -7.77 0 0 1 # ref. 2 +Zn+2 + 4Cl- = ZnCl4-2 + -log_k 0.2 + -delta_h 10.96 kcal + -gamma 5.0 0 + -Vm 28.42 28 -5.26 -3.94 2.67 0 0 0 4.62e-2 1 # ref. 2 +Zn+2 + H2O + Cl- = ZnOHCl + H+ + -log_k -7.48 +Zn+2 + 2HS- = Zn(HS)2 + -log_k 14.94 +Zn+2 + 3HS- = Zn(HS)3- + -log_k 16.1 +Zn+2 + CO3-2 = ZnCO3 + -log_k 5.3 +Zn+2 + 2CO3-2 = Zn(CO3)2-2 + -log_k 9.63 +Zn+2 + HCO3- = ZnHCO3+ + -log_k 2.1 +Zn+2 + SO4-2 = ZnSO4 + -log_k 2.37 + -delta_h 1.36 kcal + -Vm 2.51 0 18.8 # ref. 2 +Zn+2 + 2SO4-2 = Zn(SO4)2-2 + -log_k 3.28 + -Vm 10.9 0 -98.7 0 0 0 24 0 -0.236 1 # ref. 2 +Zn+2 + Br- = ZnBr+ + -log_k -0.58 +Zn+2 + 2Br- = ZnBr2 + -log_k -0.98 +Zn+2 + F- = ZnF+ + -log_k 1.15 + -delta_h 2.22 kcal +Cd+2 + H2O = CdOH+ + H+ + -log_k -10.08 + -delta_h 13.1 kcal +Cd+2 + 2 H2O = Cd(OH)2 + 2 H+ + -log_k -20.35 +Cd+2 + 3 H2O = Cd(OH)3- + 3 H+ + -log_k -33.3 +Cd+2 + 4 H2O = Cd(OH)4-2 + 4 H+ + -log_k -47.35 +2Cd+2 + H2O = Cd2OH+3 + H+ + -log_k -9.39 + -delta_h 10.9 kcal +Cd+2 + H2O + Cl- = CdOHCl + H+ + -log_k -7.404 + -delta_h 4.355 kcal +Cd+2 + NO3- = CdNO3+ + -log_k 0.4 + -delta_h -5.2 kcal + -Vm 5.95 0 -1.11 0 2.67 7 0 0 1.53e-2 1 # ref. 2 +Cd+2 + Cl- = CdCl+ + -log_k 1.98 + -delta_h 0.59 kcal + -Vm 5.69 0 -30.2 0 0 6 0 0 0.112 1 # ref. 2 +Cd+2 + 2 Cl- = CdCl2 + -log_k 2.6 + -delta_h 1.24 kcal + -Vm 5.53 # ref. 2 +Cd+2 + 3 Cl- = CdCl3- + -log_k 2.4 + -delta_h 3.9 kcal + -Vm 4.6 0 83.9 0 0 0 0 0 0 1 # ref. 2 +Cd+2 + CO3-2 = CdCO3 + -log_k 2.9 +Cd+2 + 2CO3-2 = Cd(CO3)2-2 + -log_k 6.4 +Cd+2 + HCO3- = CdHCO3+ + -log_k 1.5 +Cd+2 + SO4-2 = CdSO4 + -log_k 2.46 + -delta_h 1.08 kcal + -Vm 10.4 0 57.9 # ref. 2 +Cd+2 + 2SO4-2 = Cd(SO4)2-2 + -log_k 3.5 + -Vm -6.29 0 -93 0 9.5 7 0 0 0 1 # ref. 2 +Cd+2 + Br- = CdBr+ + -log_k 2.17 + -delta_h -0.81 kcal +Cd+2 + 2Br- = CdBr2 + -log_k 2.9 +Cd+2 + F- = CdF+ + -log_k 1.1 +Cd+2 + 2F- = CdF2 + -log_k 1.5 +Cd+2 + HS- = CdHS+ + -log_k 10.17 +Cd+2 + 2HS- = Cd(HS)2 + -log_k 16.53 +Cd+2 + 3HS- = Cd(HS)3- + -log_k 18.71 +Cd+2 + 4HS- = Cd(HS)4-2 + -log_k 20.9 +Pb+2 + H2O = PbOH+ + H+ + -log_k -7.71 +Pb+2 + 2 H2O = Pb(OH)2 + 2 H+ + -log_k -17.12 +Pb+2 + 3 H2O = Pb(OH)3- + 3 H+ + -log_k -28.06 +Pb+2 + 4 H2O = Pb(OH)4-2 + 4 H+ + -log_k -39.7 +2 Pb+2 + H2O = Pb2OH+3 + H+ + -log_k -6.36 +Pb+2 + Cl- = PbCl+ + -log_k 1.6 + -delta_h 4.38 kcal + -Vm 2.8934 -.7165 6.0316 -2.7494 .1281 6 # supcrt +Pb+2 + 2 Cl- = PbCl2 + -log_k 1.8 + -delta_h 1.08 kcal + -Vm 6.5402 8.1879 2.5318 -3.1175 -.0300 # supcrt +Pb+2 + 3 Cl- = PbCl3- + -log_k 1.7 + -delta_h 2.17 kcal + -Vm 11.0396 19.1743 -1.7863 -3.5717 .7356 # supcrt +Pb+2 + 4 Cl- = PbCl4-2 + -log_k 1.38 + -delta_h 3.53 kcal + -Vm 16.4150 32.2997 -6.9452 -4.1143 2.3118 # supcrt +Pb+2 + CO3-2 = PbCO3 + -log_k 7.24 +Pb+2 + 2 CO3-2 = Pb(CO3)2-2 + -log_k 10.64 +Pb+2 + HCO3- = PbHCO3+ + -log_k 2.9 +Pb+2 + SO4-2 = PbSO4 + -log_k 2.75 +Pb+2 + 2 SO4-2 = Pb(SO4)2-2 + -log_k 3.47 +Pb+2 + 2HS- = Pb(HS)2 + -log_k 15.27 +Pb+2 + 3HS- = Pb(HS)3- + -log_k 16.57 +3Pb+2 + 4H2O = Pb3(OH)4+2 + 4H+ + -log_k -23.88 + -delta_h 26.5 kcal +Pb+2 + NO3- = PbNO3+ + -log_k 1.17 +Pb+2 + Br- = PbBr+ + -log_k 1.77 + -delta_h 2.88 kcal +Pb+2 + 2Br- = PbBr2 + -log_k 1.44 +Pb+2 + F- = PbF+ + -log_k 1.25 +Pb+2 + 2F- = PbF2 + -log_k 2.56 +Pb+2 + 3F- = PbF3- + -log_k 3.42 +Pb+2 + 4F- = PbF4-2 + -log_k 3.1 + +PHASES +Calcite + CaCO3 = CO3-2 + Ca+2 + -log_k -8.48 + -delta_h -2.297 kcal + -analytic -171.9065 -0.077993 2839.319 71.595 + -Vm 36.9 cm3/mol # MW (100.09 g/mol) / rho (2.71 g/cm3) +Aragonite + CaCO3 = CO3-2 + Ca+2 + -log_k -8.336 + -delta_h -2.589 kcal + -analytic -171.9773 -0.077993 2903.293 71.595 + -Vm 34.04 +Dolomite + CaMg(CO3)2 = Ca+2 + Mg+2 + 2 CO3-2 + -log_k -17.09 + -delta_h -9.436 kcal + -Vm 64.5 +Siderite + FeCO3 = Fe+2 + CO3-2 + -log_k -10.89 + -delta_h -2.480 kcal + -Vm 29.2 +Rhodochrosite + MnCO3 = Mn+2 + CO3-2 + -log_k -11.13 + -delta_h -1.430 kcal + -Vm 31.1 +Strontianite + SrCO3 = Sr+2 + CO3-2 + -log_k -9.271 + -delta_h -0.400 kcal + -analytic 155.0305 0.0 -7239.594 -56.58638 + -Vm 39.69 +Witherite + BaCO3 = Ba+2 + CO3-2 + -log_k -8.562 + -delta_h 0.703 kcal + -analytic 607.642 0.121098 -20011.25 -236.4948 + -Vm 46 +Gypsum + CaSO4:2H2O = Ca+2 + SO4-2 + 2 H2O + -log_k -4.58 + -delta_h -0.109 kcal + -analytic 68.2401 0.0 -3221.51 -25.0627 + -Vm 73.9 # 172.18 / 2.33 (Vm H2O = 13.9 cm3/mol) +Anhydrite + CaSO4 = Ca+2 + SO4-2 + -log_k -4.36 + -delta_h -1.710 kcal + -analytic 84.90 0 -3135.12 -31.79 # 50 - 160oC, 1 - 1e3 atm, anhydrite dissolution, Blount and Dickson, 1973, Am. Mineral. 58, 323. + -Vm 46.1 # 136.14 / 2.95 +Celestite + SrSO4 = Sr+2 + SO4-2 + -log_k -6.63 + -delta_h -4.037 kcal +# -analytic -14805.9622 -2.4660924 756968.533 5436.3588 -40553604.0 + -analytic -7.14 6.11e-3 75 0 0 -1.79e-5 # Howell et al., 1992, JCED 37, 464. + -Vm 46.4 +Barite + BaSO4 = Ba+2 + SO4-2 + -log_k -9.97 + -delta_h 6.35 kcal + -analytic 136.035 0.0 -7680.41 -48.595 + -Vm 51.9 +Hydroxyapatite + Ca5(PO4)3OH + 4 H+ = H2O + 3 HPO4-2 + 5 Ca+2 + -log_k -3.421 + -delta_h -36.155 kcal + -Vm 128.9 +Fluorite + CaF2 = Ca+2 + 2 F- + -log_k -10.6 + -delta_h 4.69 kcal + -analytic 66.348 0.0 -4298.2 -25.271 + -Vm 15.7 +SiO2(a) + SiO2 + 2 H2O = H4SiO4 + -log_k -2.71 + -delta_h 3.340 kcal + -analytic -0.26 0.0 -731.0 +Chalcedony + SiO2 + 2 H2O = H4SiO4 + -log_k -3.55 + -delta_h 4.720 kcal + -analytic -0.09 0.0 -1032.0 + -Vm 23.1 +Quartz + SiO2 + 2 H2O = H4SiO4 + -log_k -3.98 + -delta_h 5.990 kcal + -analytic 0.41 0.0 -1309.0 + -Vm 22.67 +Gibbsite + Al(OH)3 + 3 H+ = Al+3 + 3 H2O + -log_k 8.11 + -delta_h -22.800 kcal +Al(OH)3(a) + Al(OH)3 + 3 H+ = Al+3 + 3 H2O + -log_k 10.8 + -delta_h -26.500 kcal +Kaolinite + Al2Si2O5(OH)4 + 6 H+ = H2O + 2 H4SiO4 + 2 Al+3 + -log_k 7.435 + -delta_h -35.300 kcal +Albite + NaAlSi3O8 + 8 H2O = Na+ + Al(OH)4- + 3 H4SiO4 + -log_k -18.002 + -delta_h 25.896 kcal +Anorthite + CaAl2Si2O8 + 8 H2O = Ca+2 + 2 Al(OH)4- + 2 H4SiO4 + -log_k -19.714 + -delta_h 11.580 kcal +K-feldspar + KAlSi3O8 + 8 H2O = K+ + Al(OH)4- + 3 H4SiO4 + -log_k -20.573 + -delta_h 30.820 kcal +K-mica + KAl3Si3O10(OH)2 + 10 H+ = K+ + 3 Al+3 + 3 H4SiO4 + -log_k 12.703 + -delta_h -59.376 kcal +Chlorite(14A) + Mg5Al2Si3O10(OH)8 + 16H+ = 5Mg+2 + 2Al+3 + 3H4SiO4 + 6H2O + -log_k 68.38 + -delta_h -151.494 kcal +Ca-Montmorillonite + Ca0.165Al2.33Si3.67O10(OH)2 + 12 H2O = 0.165Ca+2 + 2.33 Al(OH)4- + 3.67 H4SiO4 + 2 H+ + -log_k -45.027 + -delta_h 58.373 kcal +Talc + Mg3Si4O10(OH)2 + 4 H2O + 6 H+ = 3 Mg+2 + 4 H4SiO4 + -log_k 21.399 + -delta_h -46.352 kcal +Illite + K0.6Mg0.25Al2.3Si3.5O10(OH)2 + 11.2H2O = 0.6K+ + 0.25Mg+2 + 2.3Al(OH)4- + 3.5H4SiO4 + 1.2H+ + -log_k -40.267 + -delta_h 54.684 kcal +Chrysotile + Mg3Si2O5(OH)4 + 6 H+ = H2O + 2 H4SiO4 + 3 Mg+2 + -log_k 32.2 + -delta_h -46.800 kcal + -analytic 13.248 0.0 10217.1 -6.1894 +Sepiolite + Mg2Si3O7.5OH:3H2O + 4 H+ + 0.5H2O = 2 Mg+2 + 3 H4SiO4 + -log_k 15.760 + -delta_h -10.700 kcal +Sepiolite(d) + Mg2Si3O7.5OH:3H2O + 4 H+ + 0.5H2O = 2 Mg+2 + 3 H4SiO4 + -log_k 18.66 +Hematite + Fe2O3 + 6 H+ = 2 Fe+3 + 3 H2O + -log_k -4.008 + -delta_h -30.845 kcal +Goethite + FeOOH + 3 H+ = Fe+3 + 2 H2O + -log_k -1.0 + -delta_h -14.48 kcal +Fe(OH)3(a) + Fe(OH)3 + 3 H+ = Fe+3 + 3 H2O + -log_k 4.891 +Pyrite + FeS2 + 2 H+ + 2 e- = Fe+2 + 2 HS- + -log_k -18.479 + -delta_h 11.300 kcal +FeS(ppt) + FeS + H+ = Fe+2 + HS- + -log_k -3.915 +Mackinawite + FeS + H+ = Fe+2 + HS- + -log_k -4.648 +Sulfur + S + 2H+ + 2e- = H2S + -log_k 4.882 + -delta_h -9.5 kcal +Vivianite + Fe3(PO4)2:8H2O = 3 Fe+2 + 2 PO4-3 + 8 H2O + -log_k -36.0 +Pyrolusite # H2O added for surface calc's + MnO2:H2O + 4 H+ + 2 e- = Mn+2 + 3 H2O + -log_k 41.38 + -delta_h -65.110 kcal +Hausmannite + Mn3O4 + 8 H+ + 2 e- = 3 Mn+2 + 4 H2O + -log_k 61.03 + -delta_h -100.640 kcal +Manganite + MnOOH + 3 H+ + e- = Mn+2 + 2 H2O + -log_k 25.34 +Pyrochroite + Mn(OH)2 + 2 H+ = Mn+2 + 2 H2O + -log_k 15.2 +Halite + NaCl = Cl- + Na+ + log_k 1.570 + -delta_h 1.37 + #-analytic -713.4616 -.1201241 37302.21 262.4583 -2106915. + -Vm 27.1 +Sylvite + KCl = K+ + Cl- + log_k 0.900 + -delta_h 8.5 + # -analytic 3.984 0.0 -919.55 + Vm 37.5 +CO2(g) + CO2 = CO2 + -log_k -1.468 + -delta_h -4.776 kcal + -analytic 109.534 1.9913e-2 -6986.04 -40.83 669370 + -T_c 304.2 # critical T, K + -P_c 72.86 # critical P, atm + -Omega 0.225 # acentric factor +H2O(g) + H2O = H2O + -log_k 1.506; delta_h -44.03 kJ + -T_c 647.3 + -P_c 217.60 + -Omega 0.344 + -analytic -16.5066 -2.0013E-3 2710.7 3.7646 0 2.24E-6 + +# Gases from LLNL... +O2(g) + O2 = O2 + -log_k -2.8983 + -analytic -7.5001 7.8981e-3 0.0 0.0 2.0027e5 + -T_c 154.6 + -P_c 49.80 + -Omega 0.021 +### MDL species added just for syntax - without parenthesis +O2g + O2 = O2 + log_k -2.8983 + -analytic -7.5001 7.8981e-3 0.0 0.0 2.0027e+5 + -T_c 154.6 + -P_c 49.80 + -Omega 0.021 +H2(g) + H2 = H2 + -log_k -3.1050 + -delta_h -4.184 kJ + -analytic -9.3114 4.6473e-3 -49.335 1.4341 1.2815e5 + -T_c 33.2 + -P_c 12.80 + -Omega -0.225 +N2(g) + N2 = N2 + -log_k -3.1864 + -analytic -58.453 1.818e-3 3199 17.909 -27460 + -T_c 126.2 + -P_c 33.50 + -Omega 0.039 +H2S(g) + H2S = H+ + HS- + -log_k -7.9759 + -analytic -97.354 -3.1576e-2 1.8285e3 37.44 28.56 + -T_c 373.2 + -P_c 88.20 + -Omega 0.1 +CH4(g) + CH4 = CH4 + -log_k -2.8502 + -analytic -24.027 4.7146e-3 372.27 6.4264 2.3362e5 + -T_c 190.6 + -P_c 45.40 + -Omega 0.008 +NH3(g) + NH3 = NH3 + -log_k 1.7966 + -analytic -18.758 3.3670e-4 2.5113e3 4.8619 39.192 + -T_c 405.6 + -P_c 111.3 + -Omega 0.25 +#Amm(g) +# Amm = Amm +# -log_k 1.7966 +# -analytic -18.758 3.3670e-4 2.5113e3 4.8619 39.192 +# -T_c 405.6 +# -P_c 111.3 +# -Omega 0.25 +# redox-uncoupled gases +Oxg(g) + Oxg = Oxg + -analytic -7.5001 7.8981e-3 0.0 0.0 2.0027e5 + -T_c 154.6 ; -P_c 49.80 ; -Omega 0.021 +Hdg(g) + Hdg = Hdg + -analytic -9.3114 4.6473e-3 -49.335 1.4341 1.2815e5 + -T_c 33.2 ; -P_c 12.80 ; -Omega -0.225 +Ntg(g) + Ntg = Ntg + -analytic -58.453 1.81800e-3 3199 17.909 -27460 + T_c 126.2 ; -P_c 33.50 ; -Omega 0.039 +Mtg(g) + Mtg = Mtg + -analytic -24.027 4.7146e-3 3.7227e2 6.4264 2.3362e5 + -T_c 190.6 ; -P_c 45.40 ; -Omega 0.008 +H2Sg(g) + H2Sg = H+ + HSg- + -analytic -97.354 -3.1576e-2 1.8285e3 37.44 28.56 + -T_c 373.2 ; -P_c 88.20 ; -Omega 0.1 +Melanterite + FeSO4:7H2O = 7 H2O + Fe+2 + SO4-2 + -log_k -2.209 + -delta_h 4.910 kcal + -analytic 1.447 -0.004153 0.0 0.0 -214949.0 +Alunite + KAl3(SO4)2(OH)6 + 6 H+ = K+ + 3 Al+3 + 2 SO4-2 + 6H2O + -log_k -1.4 + -delta_h -50.250 kcal +Jarosite-K + KFe3(SO4)2(OH)6 + 6 H+ = 3 Fe+3 + 6 H2O + K+ + 2 SO4-2 + -log_k -9.21 + -delta_h -31.280 kcal +Zn(OH)2(e) + Zn(OH)2 + 2 H+ = Zn+2 + 2 H2O + -log_k 11.5 +Smithsonite + ZnCO3 = Zn+2 + CO3-2 + -log_k -10.0 + -delta_h -4.36 kcal +Sphalerite + ZnS + H+ = Zn+2 + HS- + -log_k -11.618 + -delta_h 8.250 kcal +Willemite 289 + Zn2SiO4 + 4H+ = 2Zn+2 + H4SiO4 + -log_k 15.33 + -delta_h -33.37 kcal +Cd(OH)2 + Cd(OH)2 + 2 H+ = Cd+2 + 2 H2O + -log_k 13.65 +Otavite 315 + CdCO3 = Cd+2 + CO3-2 + -log_k -12.1 + -delta_h -0.019 kcal +CdSiO3 328 + CdSiO3 + H2O + 2H+ = Cd+2 + H4SiO4 + -log_k 9.06 + -delta_h -16.63 kcal +CdSO4 329 + CdSO4 = Cd+2 + SO4-2 + -log_k -0.1 + -delta_h -14.74 kcal +Cerrusite 365 + PbCO3 = Pb+2 + CO3-2 + -log_k -13.13 + -delta_h 4.86 kcal +Anglesite 384 + PbSO4 = Pb+2 + SO4-2 + -log_k -7.79 + -delta_h 2.15 kcal +Pb(OH)2 389 + Pb(OH)2 + 2H+ = Pb+2 + 2H2O + -log_k 8.15 + -delta_h -13.99 kcal + +RATES +Calcite +-start + 10 moles=0 + 20 IF ((M<=0) and (SI("Calcite")<0)) then goto 200 + 30 R=8.314462 # in J*K-1*mol-1 + 40 deltaT=1/TK-1/298.15 # wird für 40°C berechnet; TK is temp in Kelvin + 50 e=2.718282 # Eulersche Zahl + ## mechanism 1 (acid) + 60 Ea=14400 # Aktivierungsenergie in J/mol => 65.0 in KJ/mol + 70 logK25=-0.3 # Reaktionskonstante 25C mol/m2/s + 90 mech_a=(10^logK25)*(e^(-Ea/R*deltaT))*ACT("H+") ## removed exponent + ## base term (neutral mechanism) + 100 Ea=23500 + 110 logK25=-5.81 + 120 mech_c=(10^logK25)*(e^(-Ea/R*deltaT)) + 130 rate=mech_a+mech_c + 140 IF (SI("Calcite")<0 AND M>0) then moles=parm(1)*rate*(1-SR("Calcite")) # dissolution + ## 145 IF SI("Calcite")>0 then moles=parm(1)*M*rate*(-1+SR("Calcite")) # precipitation + ## 150 moles=parm(1)*rate*(1-SR("Calcite")) # precipitation + 200 save moles*time +-end + +Dolomite +-start + 10 moles=0 + 20 IF ((M<=0) and (SI("Dolomite")<0)) then goto 200 + 30 R=8.314462 # in J*K-1*mol-1 + 40 deltaT=1/TK-1/298.15 # wird für 40°C berechnet; TK is temp in Kelvin + 50 e=2.718282 # Eulersche Zahl + ## mechanism 1 (acid) + 60 Ea=36100 # Aktivierungsenergie in J/mol => 65.0 in KJ/mol + 70 logK25=-3.19 # Reaktionskonstante 25C mol/m2/s + 90 mech_a=(10^logK25)*(e^(-Ea/R*deltaT))*ACT("H+")^0.5 ## removed exponent + ## base term (neutral mechanism) + 100 Ea=52200 + 110 logK25=-7.53 + 120 mech_c=(10^logK25)*(e^(-Ea/R*deltaT)) + 130 rate=mech_a+mech_c + ## 140 IF SI("Dolomite")<0 then moles=parm(1)*rate*(1-SR("Dolomite")) # dissolution + ## 140 IF SI("Dolomite")<0 then moles=parm(1)*rate*(1-SR("Dolomite")) # dissolution + 150 moles=parm(1)*rate*(1-SR("Dolomite")) # precipitation + 200 save moles*time +-end + +END diff --git a/bin/run_poet.sh b/bin/run_poet.sh new file mode 100644 index 000000000..27e590a5c --- /dev/null +++ b/bin/run_poet.sh @@ -0,0 +1,19 @@ +#!/bin/bash +#SBATCH --job-name=dolo_proto1_eps01_no_zeroabs +#SBATCH --output=dolo_proto1_eps01_no_zeroabs_%j.out +#SBATCH --error=dolo_proto1_eps01_no_zeroabs_%j.err +#SBATCH --partition=long +#SBATCH --nodes=6 +#SBATCH --ntasks-per-node=24 +#SBATCH --ntasks=144 +#SBATCH --exclusive +#SBATCH --time=3-00:00:00 + + +source /etc/profile.d/modules.sh +module purge +module load cmake gcc openmpi + +#mpirun -n 144 ./poet dolo_fgcs_3.R dolo_fgcs_3.qs2 dolo_only_pqc +mpirun -n 144 ./poet --interp dolo_fgcs_3_rt.R dolo_fgcs_3.qs2 dolo_proto1_eps01_no_zeroabs +#mpirun -n 144 ./poet --interp barite_fgcs_4_new/barite_fgcs_4_new_rt.R barite_fgcs_4_new/barite_fgcs_4_new.qs2 barite \ No newline at end of file diff --git a/share/poet/barite/barite_het.qs2 b/share/poet/barite/barite_het.qs2 index 9c0e0c843..72bb5d3e6 100644 Binary files a/share/poet/barite/barite_het.qs2 and b/share/poet/barite/barite_het.qs2 differ diff --git a/share/poet/surfex/PoetEGU_surfex_500.qs2 b/share/poet/surfex/PoetEGU_surfex_500.qs2 index 135c62b06..d27b5a257 100644 Binary files a/share/poet/surfex/PoetEGU_surfex_500.qs2 and b/share/poet/surfex/PoetEGU_surfex_500.qs2 differ diff --git a/src/Chemistry/ChemistryModule.hpp b/src/Chemistry/ChemistryModule.hpp index 746e8e06c..2f1cae78e 100644 --- a/src/Chemistry/ChemistryModule.hpp +++ b/src/Chemistry/ChemistryModule.hpp @@ -102,6 +102,7 @@ public: this->ai_surrogate_enabled = setup.ai_surrogate_enabled; this->base_totals = setup.base_totals; + this->ctrl_file_out_dir = setup.dht_out_dir; if (this->dht_enabled || this->interp_enabled) { this->initializeDHT(setup.dht_size_mb, this->params.dht_species, @@ -258,7 +259,7 @@ public: std::vector ai_surrogate_validity_vector; - void SetControlModule(poet::ControlModule *ctrl) { control_module = ctrl; } + void SetControlModule(poet::ControlModule *ctrl) { control = ctrl; } void SetDhtEnabled(bool enabled) { this->dht_enabled = enabled; } bool GetDhtEnabled() const { return this->dht_enabled; } @@ -266,7 +267,8 @@ public: void SetInterpEnabled(bool enabled) { this->interp_enabled = enabled; } bool GetInterpEnabled() const { return interp_enabled; } - void SetWarmupEnabled(bool enabled) { this->warmup_enabled = enabled; } + void SetStabEnabled(bool enabled) { this->stab_enabled = enabled; } + bool GetStabEnabled() const { return stab_enabled; } void SetControlCellIds(const std::vector &ids) { this->ctrl_cell_ids = std::unordered_set(ids.begin(), ids.end()); @@ -285,13 +287,13 @@ protected: enum { CHEM_FIELD_INIT, - // CHEM_DHT_ENABLE, + CHEM_DHT_ENABLE, + CHEM_IP_ENABLE, + CHEM_CTRL_ENABLE, + CHEM_CTRL_FLAGS, 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, @@ -300,6 +302,9 @@ protected: CHEM_AI_BCAST_VALIDITY }; + /* broadcasted only every control iteration */ + enum { DHT_ENABLE = 1u << 0, IP_ENABLE = 1u << 1, STAB_ENABLE = 1u << 2 }; + enum { LOOP_WORK, LOOP_END, LOOP_CTRL }; enum { @@ -358,7 +363,7 @@ protected: void WorkerDoWork(MPI_Status &probe_status, int double_count, struct worker_s &timings); - void WorkerPostIter(MPI_Status &prope_status, uint32_t iteration); + void WorkerPostIter(MPI_Status &probe_status, uint32_t iteration); void WorkerPostSim(uint32_t iteration); void WorkerWriteDHTDump(uint32_t iteration); @@ -370,10 +375,6 @@ protected: void WorkerRunWorkPackage(WorkPackage &work_package, double dSimTime, double dTimestep); - void ProcessControlWorkPackage(std::vector> &input, - double current_sim_time, double dt, - struct worker_s &timings); - std::vector CalculateWPSizesVector(uint32_t n_cells, uint32_t wp_size) const; std::vector shuffleField(const std::vector &in_field, @@ -388,9 +389,25 @@ protected: void BCastStringVec(std::vector &io); - int packResultsIntoBuffer(std::vector &mpi_buffer, int base_count, - const WorkPackage &wp, - const WorkPackage &wp_control); + void processCtrlPkgs(std::vector> &input, + double current_sim_time, double dt, + struct worker_s &timings); + + void copyPkgs(const WorkPackage &wp, std::vector &mpi_buffer); + + inline int buildCtrlFlags(bool dht, bool interp, bool stab) { + int flags = 0; + + if (dht) + flags |= DHT_ENABLE; + if (interp) + flags |= IP_ENABLE; + if (stab) + flags |= STAB_ENABLE; + return flags; + } + + inline bool hasFlag(int flags, int type) { return (flags & type) != 0; } int comm_size, comm_rank; MPI_Comm group_comm; @@ -410,7 +427,7 @@ protected: bool ai_surrogate_enabled{false}; - static constexpr uint32_t BUFFER_OFFSET = 6; + static constexpr uint32_t BUFFER_OFFSET = 5; inline void ChemBCast(void *buf, int count, MPI_Datatype datatype) const { MPI_Bcast(buf, count, datatype, 0, this->group_comm); @@ -447,13 +464,14 @@ protected: std::unique_ptr pqc_runner; - poet::ControlModule *control_module = nullptr; + std::string ctrl_file_out_dir; + + poet::ControlModule *control = nullptr; std::vector mpi_surr_buffer; bool control_enabled{false}; - bool warmup_enabled{false}; - + bool stab_enabled{false}; std::unordered_set ctrl_cell_ids; std::vector> control_batch; }; diff --git a/src/Chemistry/MasterFunctions.cpp b/src/Chemistry/MasterFunctions.cpp index 3afa092d9..0b4bef0e3 100644 --- a/src/Chemistry/MasterFunctions.cpp +++ b/src/Chemistry/MasterFunctions.cpp @@ -280,11 +280,13 @@ inline void poet::ChemistryModule::MasterSendPkgs( // current work package start location in field send_buffer[end_of_wp + 4] = wp_start_index; // 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(flags); - + */ /* ATTENTION Worker p has rank p+1 */ // MPI_Send(send_buffer, end_of_wp + BUFFER_OFFSET, MPI_DOUBLE, p + 1, // LOOP_WORK, this->group_comm); @@ -441,6 +443,14 @@ void poet::ChemistryModule::MasterRunParallel(double dt) { MPI_INT); } + if (control->shouldBcastFlags()) { + int ftype = CHEM_CTRL_FLAGS; + PropagateFunctionType(ftype); + uint32_t ctrl_flags = buildCtrlFlags( + this->dht_enabled, this->interp_enabled, this->stab_enabled); + ChemBCast(&ctrl_flags, 1, MPI_INT); + } + ftype = CHEM_WORK_LOOP; PropagateFunctionType(ftype); @@ -512,6 +522,9 @@ void poet::ChemistryModule::MasterRunParallel(double dt) { chem_field = out_vec; /* do master stuff */ + std::cout << "[DEBUG] control_batch.size() = " + << this->control_batch.size() << std::endl; + if (!this->control_batch.empty()) { std::cout << "[Master] Processing " << this->control_batch.size() << " control cells for comparison." << std::endl; @@ -536,8 +549,11 @@ void poet::ChemistryModule::MasterRunParallel(double dt) { } metrics_a = MPI_Wtime(); - control_module->computeSpeciesErrorMetrics(this->control_batch, - surrogate_batch); + control->computeErrorMetrics(this->control_batch, surrogate_batch, + prop_names); + + control->writeErrorMetrics(ctrl_file_out_dir, prop_names); + metrics_b = MPI_Wtime(); this->metrics_t += metrics_b - metrics_a; @@ -553,12 +569,20 @@ void poet::ChemistryModule::MasterRunParallel(double dt) { this->seq_t += seq_d - seq_c; /* end time measurement of whole chemistry simulation */ + std::optional target = control->getRollbackTarget(prop_names); + int flush = target.has_value() ? 1 : 0; /* advise workers to end chemistry iteration */ for (int i = 1; i < this->comm_size; i++) { - MPI_Send(NULL, 0, MPI_DOUBLE, i, LOOP_END, this->group_comm); + MPI_Send(&flush, 1, MPI_INT, i, LOOP_END, this->group_comm); } + /* + if (flush) { + control->clearFlushRequest(); + } + */ + this->simtime += dt; iteration++; } diff --git a/src/Chemistry/WorkerFunctions.cpp b/src/Chemistry/WorkerFunctions.cpp index 9bff153ff..54fbcb69f 100644 --- a/src/Chemistry/WorkerFunctions.cpp +++ b/src/Chemistry/WorkerFunctions.cpp @@ -59,6 +59,15 @@ void poet::ChemistryModule::WorkerLoop() { MPI_INT, 0, this->group_comm); break; } + case CHEM_CTRL_FLAGS: { + int flags = 0; + ChemBCast(&flags, 1, MPI_INT); + this->dht_enabled = hasFlag(flags, DHT_ENABLE); + this->interp_enabled = hasFlag(flags, IP_ENABLE); + this->stab_enabled = hasFlag(flags, STAB_ENABLE); + break; + } + case CHEM_WORK_LOOP: { WorkerProcessPkgs(timings, iteration); break; @@ -116,7 +125,16 @@ void poet::ChemistryModule::WorkerProcessPkgs(struct worker_s &timings, } } -void poet::ChemistryModule::ProcessControlWorkPackage( +void poet::ChemistryModule::copyPkgs(const WorkPackage &wp, + std::vector &mpi_buffer) { + + for (std::size_t wp_i = 0; wp_i < wp.size; wp_i++) { + std::copy(wp.output[wp_i].begin(), wp.output[wp_i].end(), + mpi_buffer.begin() + this->prop_count * wp_i); + } +} + +void poet::ChemistryModule::processCtrlPkgs( std::vector> &input, double current_sim_time, double dt, struct worker_s &timings) { @@ -137,10 +155,7 @@ void poet::ChemistryModule::ProcessControlWorkPackage( timings.ctrl_phreeqc_t += phreeqc_end - phreeqc_start; - for (std::size_t wp_i = 0; wp_i < control_wp.size; wp_i++) { - std::copy(control_wp.output[wp_i].begin(), control_wp.output[wp_i].end(), - mpi_buffer.begin() + this->prop_count * wp_i); - } + copyPkgs(control_wp, mpi_buffer); MPI_Request send_req; MPI_Isend(mpi_buffer.data(), mpi_buffer.size(), MPI_DOUBLE, 0, LOOP_CTRL, @@ -194,13 +209,16 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status, wp_start_index = mpi_buffer[count + 4]; // read packed control flags + /* flags = static_cast(mpi_buffer[count + 5]); this->interp_enabled = (flags & 1) != 0; this->dht_enabled = (flags & 2) != 0; this->warmup_enabled = (flags & 4) != 0; this->control_enabled = (flags & 8) != 0; + */ - /*std::cout << "warmup_enabled is " << warmup_enabled << ", control_enabled is " + /*std::cout << "warmup_enabled is " << warmup_enabled << ", control_enabled is + " << control_enabled << ", dht_enabled is " << dht_enabled << ", interp_enabled is " << interp_enabled << std::endl;*/ @@ -212,7 +230,7 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status, } // std::cout << this->comm_rank << ":" << counter++ << std::endl; - if (dht_enabled || interp_enabled || warmup_enabled) { + if (dht_enabled || interp_enabled || stab_enabled) { dht->prepareKeys(s_curr_wp.input, dt); } @@ -241,18 +259,19 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status, for (std::size_t wp_i = 0; wp_i < s_curr_wp.size; wp_i++) { uint32_t cell_id = s_curr_wp.input[wp_i][0]; - - bool is_control_cell = this->ctrl_cell_ids.find(cell_id) != this->ctrl_cell_ids.end(); - bool used_surrogate = s_curr_wp.mapping[wp_i] != CHEM_PQC; - - if (is_control_cell && used_surrogate) { + + bool is_ctrl_cell = + this->ctrl_cell_ids.find(cell_id) != this->ctrl_cell_ids.end(); + bool used_sur = s_curr_wp.mapping[wp_i] != CHEM_PQC; + + if (is_ctrl_cell && used_sur) { control_batch.push_back(s_curr_wp.input[wp_i]); control_cells_processed++; if (control_batch.size() == s_curr_wp.size || control_cells_processed == this->ctrl_cell_ids.size()) { - ProcessControlWorkPackage(control_batch, current_sim_time, dt, timings); + processCtrlPkgs(control_batch, current_sim_time, dt, timings); control_batch.clear(); control_cells_processed = 0; } @@ -265,23 +284,20 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status, phreeqc_time_end = MPI_Wtime(); - for (std::size_t wp_i = 0; wp_i < s_curr_wp.size; wp_i++) { - std::copy(s_curr_wp.output[wp_i].begin(), s_curr_wp.output[wp_i].end(), - mpi_buffer.begin() + this->prop_count * wp_i); - } + copyPkgs(s_curr_wp, mpi_buffer); /* send results to master */ MPI_Request send_req; MPI_Isend(mpi_buffer.data(), count, MPI_DOUBLE, 0, LOOP_WORK, MPI_COMM_WORLD, &send_req); - if (dht_enabled || interp_enabled || warmup_enabled) { + if (dht_enabled || interp_enabled || stab_enabled) { /* write results to DHT */ dht_fill_start = MPI_Wtime(); dht->fillDHT(s_curr_wp); dht_fill_end = MPI_Wtime(); - if (interp_enabled || warmup_enabled) { + if (interp_enabled || stab_enabled) { interp->writePairs(); } timings.dht_fill += dht_fill_end - dht_fill_start; @@ -291,10 +307,18 @@ void poet::ChemistryModule::WorkerDoWork(MPI_Status &probe_status, MPI_Wait(&send_req, MPI_STATUS_IGNORE); } -void poet::ChemistryModule::WorkerPostIter(MPI_Status &prope_status, +void poet::ChemistryModule::WorkerPostIter(MPI_Status &probe_status, uint32_t iteration) { - MPI_Recv(NULL, 0, MPI_DOUBLE, 0, LOOP_END, this->group_comm, - MPI_STATUS_IGNORE); + int size, flush_request = 0; + MPI_Get_count(&probe_status, MPI_INT, &size); + + if (size == 1) { + MPI_Recv(&flush_request, size, MPI_INT, probe_status.MPI_SOURCE, LOOP_END, + this->group_comm, MPI_STATUS_IGNORE); + } else { + MPI_Recv(NULL, 0, MPI_INT, probe_status.MPI_SOURCE, LOOP_END, + this->group_comm, MPI_STATUS_IGNORE); + } if (this->dht_enabled) { dht_hits.push_back(dht->getHits()); @@ -320,7 +344,7 @@ void poet::ChemistryModule::WorkerPostIter(MPI_Status &prope_status, const auto max_mean_idx = DHT_get_used_idx_factor(this->interp->getDHTObject(), 1); - if (max_mean_idx >= 2) { + if (max_mean_idx >= 2 || flush_request) { DHT_flush(this->interp->getDHTObject()); DHT_flush(this->dht->getDHT()); if (this->comm_rank == 2) { diff --git a/src/Control/ControlModule.cpp b/src/Control/ControlModule.cpp index efef103b9..4753f4289 100644 --- a/src/Control/ControlModule.cpp +++ b/src/Control/ControlModule.cpp @@ -4,9 +4,13 @@ #include "IO/StatsIO.hpp" #include -void poet::ControlModule::updateControlIteration(const uint32_t &iter, - const bool &dht_enabled, - const bool &interp_enabled) { +poet::ControlModule::ControlModule(const ControlConfig &config_) + : config(config_) {} + +void poet::ControlModule::beginIteration(ChemistryModule &chem, + 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*/ @@ -14,81 +18,49 @@ void poet::ControlModule::updateControlIteration(const uint32_t &iter, 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)); - } - */ + updateStabilizationPhase(chem, dht_enabled, interp_enabled); 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)); +void poet::ControlModule::updateStabilizationPhase(ChemistryModule &chem, + bool dht_enabled, + bool interp_enabled) { + if (rollback_enabled) { + if (disable_surr_counter > 0) { + --disable_surr_counter; + flush_request = false; + MSG("Rollback counter: " + std::to_string(disable_surr_counter)); } else { rollback_enabled = false; } - + } + + bool prev_stab_state = chem.GetStabEnabled(); + + // user requested DHT/INTEP? keep them disabled but enable warmup-phase so + if (global_iteration <= config.stab_interval || rollback_enabled) { + chem.SetStabEnabled(true); + chem.SetDhtEnabled(false); + chem.SetInterpEnabled(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) + "."); + + chem.SetStabEnabled(false); + chem.SetDhtEnabled(dht_enabled); + chem.SetInterpEnabled(interp_enabled); + + // Mark that we need to broadcast flags if stab phase just ended + if (prev_stab_state && !chem.GetStabEnabled()) { + stab_phase_ended = true; } } -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)); +void poet::ControlModule::writeCheckpoint(DiffusionModule &diffusion, + uint32_t &iter, + const std::string &out_dir) { + double w_check_a, w_check_b; w_check_a = MPI_Wtime(); write_checkpoint(out_dir, "checkpoint" + std::to_string(iter) + ".hdf5", @@ -96,88 +68,107 @@ void poet::ControlModule::writeCheckpointAndMetrics(DiffusionModule &diffusion, w_check_b = MPI_Wtime(); this->w_check_t += w_check_b - w_check_a; + last_checkpoint_written = iter; +} + +void poet::ControlModule::readCheckpoint(DiffusionModule &diffusion, + uint32_t ¤t_iter, + uint32_t rollback_iter, + const std::string &out_dir) { + double r_check_a, r_check_b; + + 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); + current_iter = checkpoint_read.iteration; + r_check_b = MPI_Wtime(); + r_check_t += r_check_b - r_check_a; +} + +void poet::ControlModule::writeErrorMetrics( + const std::string &out_dir, const std::vector &species) { + double stats_a, stats_b; + stats_a = MPI_Wtime(); - writeStatsToCSV(metricsHistory, species_names, out_dir, "stats_overview"); + writeStatsToCSV(metrics_history, species, out_dir, "metrics_overview"); stats_b = MPI_Wtime(); this->stats_t += stats_b - stats_a; } -bool poet::ControlModule::checkAndRollback(DiffusionModule &diffusion, - uint32_t &iter) { +uint32_t poet::ControlModule::getRollbackIter() { + + uint32_t last_iter = ((global_iteration - 1) / config.checkpoint_interval) * + config.checkpoint_interval; + + uint32_t rollback_iter = (last_iter <= last_checkpoint_written) + ? last_iter + : last_checkpoint_written; + return rollback_iter; +} + +std::optional poet::ControlModule::getRollbackTarget( + const std::vector &species) { double r_check_a, r_check_b; - if (global_iteration < stabilization_interval) { - return false; + if (metrics_history.empty()) { + MSG("No error history yet, skipping rollback check."); + rollback_enabled = false; + // flush_request = false; + return std::nullopt; } - if (metricsHistory.empty()) { - MSG("No error history yet; skipping rollback check."); - return false; - } - - const auto &mape = metricsHistory.back().mape; - + const auto &mape = metrics_history.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++) { + for (size_t col = 0; col < species.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; + if (mape[row][col] > config.mape_threshold[col]) { - 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)); + if (last_checkpoint_written == 0) { + MSG(" Threshold exceeded but no checkpoint exists yet."); + return std::nullopt; + } + rollback_enabled = true; + flush_request = true; - 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("Threshold exceeded " + species[col] + " has MAPE = " + + std::to_string(mape[row][col]) + " exceeding threshold = " + + std::to_string(config.mape_threshold[col])); + + return getRollbackIter(); } } } - MSG("All species are within their MAPE thresholds."); - - return false; + rollback_enabled = false; + // flush_request = false; + return std::nullopt; } -void poet::ControlModule::computeSpeciesErrorMetrics( +void poet::ControlModule::computeErrorMetrics( std::vector> &reference_values, - std::vector> &surrogate_values) { + std::vector> &surrogate_values, + const std::vector &species) { - const uint32_t num_cells = reference_values.size(); - const uint32_t species_count = this->species_names.size(); + const uint32_t n_cells = reference_values.size(); - std::cout << "[DEBUG] computeSpeciesErrorMetrics: num_cells=" << num_cells - << ", species_count=" << species_count << std::endl; - - SpeciesErrorMetrics metrics(num_cells, species_count, global_iteration, + SpeciesErrorMetrics metrics(n_cells, species.size(), 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++) { + for (size_t cell_i = 0; cell_i < n_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++) { + for (size_t sp_i = 0; sp_i < species.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; + const double ZERO_ABS = config.zero_abs; if (std::isnan(ref_value) || std::isnan(sur_value)) { metrics.mape[cell_i][sp_i] = 0.0; @@ -200,8 +191,45 @@ void poet::ControlModule::computeSpeciesErrorMetrics( } } } - + std::cout << "[DEBUG] metrics.id.size()=" << metrics.id.size() << std::endl; - metricsHistory.push_back(metrics); - std::cout << "[DEBUG] metricsHistory.size()=" << metricsHistory.size() << std::endl; + metrics_history.push_back(metrics); + std::cout << "[DEBUG] metricsHistory.size()=" << metrics_history.size() + << std::endl; } + +void poet::ControlModule::processCheckpoint( + DiffusionModule &diffusion, uint32_t ¤t_iter, + const std::string &out_dir, const std::vector &species) { + + if (flush_request) { + uint32_t target = getRollbackIter(); + readCheckpoint(diffusion, current_iter, target, out_dir); + + rollback_enabled = true; + rollback_count++; + disable_surr_counter = config.stab_interval; + + MSG("Restored checkpoint " + std::to_string(target) + + ", surrogates disabled for " + std::to_string(config.stab_interval)); + } else { + writeCheckpoint(diffusion, global_iteration, out_dir); + } +} + +bool poet::ControlModule::shouldBcastFlags() { + if (global_iteration == 1) { + return true; + } + + if (stab_phase_ended) { + stab_phase_ended = false; + return true; + } + + if (flush_request) { + return true; + } + + return false; +} \ No newline at end of file diff --git a/src/Control/ControlModule.hpp b/src/Control/ControlModule.hpp index 91be4f984..008919fba 100644 --- a/src/Control/ControlModule.hpp +++ b/src/Control/ControlModule.hpp @@ -7,6 +7,7 @@ #include "poet.hpp" #include +#include #include #include @@ -15,107 +16,105 @@ namespace poet { class ChemistryModule; class DiffusionModule; +struct ControlConfig { + uint32_t stab_interval = 0; + uint32_t checkpoint_interval = 0; + double zero_abs = 0.0; + std::vector mape_threshold; +}; + +struct SpeciesErrorMetrics { + std::vector id; + std::vector> mape; + std::vector> rrmse; + uint32_t iteration = 0; + uint32_t rollback_count = 0; + + SpeciesErrorMetrics(uint32_t n_cells, uint32_t n_species, uint32_t iter, + uint32_t rb_count) + : mape(n_cells, std::vector(n_species, 0.0)), + rrmse(n_cells, std::vector(n_species, 0.0)), iteration(iter), + rollback_count(rb_count) {} +}; + class ControlModule { public: - /* Control configuration*/ + explicit ControlModule(const ControlConfig &config); - // std::uint32_t global_iter = 0; - // std::uint32_t sur_disabled_counter = 0; - // std::uint32_t rollback_counter = 0; + void beginIteration(ChemistryModule &chem, const uint32_t &iter, + const bool &dht_enabled, const bool &interp_enaled); - void updateControlIteration(const uint32_t &iter, const bool &dht_enabled, - const bool &interp_enaled); + void writeErrorMetrics(const std::string &out_dir, + const std::vector &species); - void initiateWarmupPhase(bool dht_enabled, bool interp_enabled); + std::optional getRollbackTarget(); - bool checkAndRollback(DiffusionModule &diffusion, uint32_t &iter); + void computeErrorMetrics(std::vector> &reference_values, + std::vector> &surrogate_values, + const std::vector &species); - struct SpeciesErrorMetrics { - std::vector id; - std::vector> mape; - std::vector> rrmse; - uint32_t iteration; // iterations in simulation after rollbacks - uint32_t rollback_count; + void processCheckpoint(DiffusionModule &diffusion, uint32_t ¤t_iter, + const std::string &out_dir, + const std::vector &species); - SpeciesErrorMetrics(uint32_t num_cells, uint32_t species_count, - uint32_t iter, uint32_t counter) - : mape(num_cells, std::vector(species_count, 0.0)), - rrmse(num_cells, std::vector(species_count, 0.0)), - iteration(iter), rollback_count(counter) {} - }; + std::optional + getRollbackTarget(const std::vector &species); - void computeSpeciesErrorMetrics( - std::vector> &reference_values, - std::vector> &surrogate_values); + + bool shouldBcastFlags(); - std::vector metricsHistory; - - struct ControlSetup { - std::string out_dir; - std::uint32_t checkpoint_interval; - std::uint32_t penalty_interval; - std::uint32_t stabilization_interval; - std::vector species_names; - std::vector mape_threshold; - std::vector 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); + bool getFlushRequest() const { return flush_request; } + void clearFlushRequest() { flush_request = false; } auto getGlobalIteration() const noexcept { return global_iteration; } - void setChemistryModule(poet::ChemistryModule *c) { chem = c; } + // void setChemistryModule(poet::ChemistryModule *c) { chem = c; } - std::vector getMapeThreshold() const { return this->mape_threshold; } + std::vector getMapeThreshold() const { + return this->config.mape_threshold; + } std::vector 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; } + auto getUpdateCtrlLogicTime() const { return prep_t; } + auto getWriteCheckpointTime() const { return w_check_t; } + auto getReadCheckpointTime() const { return r_check_t; } + auto getWriteMetricsTime() const { return stats_t; } private: - bool rollback_enabled = false; + void updateStabilizationPhase(ChemistryModule &chem, bool dht_enabled, + bool interp_enabled); - poet::ChemistryModule *chem = nullptr; + void readCheckpoint(DiffusionModule &diffusion, uint32_t ¤t_iter, + uint32_t rollback_iter, const std::string &out_dir); + void writeCheckpoint(DiffusionModule &diffusion, uint32_t &iter, + const std::string &out_dir); + + uint32_t getRollbackIter(); + + ControlConfig config; - 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 mape_threshold; + std::uint32_t disable_surr_counter = 0; std::vector ctrl_cell_ids; + std::uint32_t last_checkpoint_written = 0; + std::uint32_t penalty_interval = 0; - std::vector species_names; - std::string out_dir; + bool rollback_enabled = false; + bool flush_request = false; + bool stab_phase_ended = false; + + bool bcast_flags = false; + + std::vector metrics_history; 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 sur_shuffled; }; } // namespace poet diff --git a/src/IO/StatsIO.cpp b/src/IO/StatsIO.cpp index 3eeddaa50..f335855d6 100644 --- a/src/IO/StatsIO.cpp +++ b/src/IO/StatsIO.cpp @@ -7,7 +7,7 @@ namespace poet { - void writeStatsToCSV(const std::vector &all_stats, + void writeStatsToCSV(const std::vector &all_stats, const std::vector &species_names, const std::string &out_dir, const std::string &filename) diff --git a/src/IO/StatsIO.hpp b/src/IO/StatsIO.hpp index 5333c4fd8..512ab0d4c 100644 --- a/src/IO/StatsIO.hpp +++ b/src/IO/StatsIO.hpp @@ -3,7 +3,7 @@ namespace poet { - void writeStatsToCSV(const std::vector &all_stats, + void writeStatsToCSV(const std::vector &all_stats, const std::vector &species_names, const std::string &out_dir, const std::string &filename); diff --git a/src/poet.cpp b/src/poet.cpp index 5fbbbf6a1..5740c0f0b 100644 --- a/src/poet.cpp +++ b/src/poet.cpp @@ -252,10 +252,10 @@ int parseInitValues(int argc, char **argv, RuntimeParameters ¶ms) { Rcpp::as>(global_rt_setup->operator[]("timesteps")); params.checkpoint_interval = Rcpp::as(global_rt_setup->operator[]("checkpoint_interval")); - params.stabilization_interval = - Rcpp::as(global_rt_setup->operator[]("stabilization_interval")); - params.penalty_interval = - Rcpp::as(global_rt_setup->operator[]("penalty_interval")); + params.stab_interval = + Rcpp::as(global_rt_setup->operator[]("stab_interval")); + params.zero_abs = + Rcpp::as(global_rt_setup->operator[]("zero_abs")); params.mape_threshold = Rcpp::as>( global_rt_setup->operator[]("mape_threshold")); params.ctrl_cell_ids = Rcpp::as>( @@ -305,7 +305,7 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, RuntimeParameters ¶ms, double dSimTime{0}; for (uint32_t iter = 1; iter < maxiter + 1; iter++) { - control.updateControlIteration(iter, params.use_dht, params.use_interp); + control.beginIteration(chem, iter, params.use_dht, params.use_interp); double start_t = MPI_Wtime(); @@ -415,7 +415,9 @@ static Rcpp::List RunMasterLoop(RInsidePOET &R, RuntimeParameters ¶ms, MSG("End of *coupling* iteration " + std::to_string(iter) + "/" + std::to_string(maxiter)); - control.applyControlLogic(diffusion, iter); + control.processCheckpoint(diffusion, iter, params.out_dir, + chem.getField().GetProps()); + // MSG(); } // END SIMULATION LOOP @@ -622,9 +624,6 @@ 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), @@ -646,16 +645,6 @@ int main(int argc, char *argv[]) { 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 { @@ -697,7 +686,14 @@ int main(int argc, char *argv[]) { DiffusionModule diffusion(init_list.getDiffusionInit(), init_list.getInitialGrid()); + ControlConfig config(run_params.stab_interval, + run_params.checkpoint_interval, run_params.zero_abs, + run_params.mape_threshold); + + ControlModule control(config); + chemistry.masterSetField(init_list.getInitialGrid()); + chemistry.SetControlModule(&control); Rcpp::List profiling = RunMasterLoop(R, run_params, diffusion, chemistry, control); diff --git a/src/poet.hpp.in b/src/poet.hpp.in index ff06eeaf0..fe7ccf79a 100644 --- a/src/poet.hpp.in +++ b/src/poet.hpp.in @@ -50,12 +50,11 @@ struct RuntimeParameters { std::string out_ext; bool print_progress = false; - std::uint32_t penalty_interval = 0; - std::uint32_t stabilization_interval = 0; + std::uint32_t stab_interval = 0; std::uint32_t checkpoint_interval = 0; - std::uint32_t control_interval = 0; + double zero_abs = 0.0; std::vector mape_threshold; - std::vector ctrl_cell_ids; + std::vector ctrl_cell_ids; static constexpr std::uint32_t WORK_PACKAGE_SIZE_DEFAULT = 32; std::uint32_t work_package_size = WORK_PACKAGE_SIZE_DEFAULT;