poet/data/SimDol1D_diffu.R
2023-07-27 15:23:05 +02:00

181 lines
4.2 KiB
R

#################################################################
## Section 1 ##
## Grid initialization ##
#################################################################
n <- 5
m <- 5
types <- c("scratch", "phreeqc", "rds")
# initsim <- c("SELECTED_OUTPUT", "-high_precision true",
# "-reset false",
# "USER_PUNCH",
# "-head C Ca Cl Mg pH pe O2g Calcite Dolomite",
# "10 PUNCH TOT(\"C\"), TOT(\"Ca\"), TOT(\"Cl\"), TOT(\"Mg\"), -LA(\"H+\"), -LA(\"e-\"), EQUI(\"O2g\"), EQUI(\"Calcite\"), EQUI(\"Dolomite\")",
# "SOLUTION 1",
# "units mol/kgw",
# "temp 25.0", "water 1",
# "pH 9.91 charge",
# "pe 4.0",
# "C 1.2279E-04",
# "Ca 1.2279E-04",
# "Cl 1E-12",
# "Mg 1E-12",
# "PURE 1",
# "O2g -0.6788 10.0",
# "Calcite 0.0 2.07E-3",
# "Dolomite 0.0 0.0",
# "END")
# needed if init type is set to "scratch"
# prop <- c("C", "Ca", "Cl", "Mg", "pH", "pe", "O2g", "Calcite", "Dolomite")
init_cell <- list(
"C" = 1.2279E-4,
"Ca" = 1.2279E-4,
"Cl" = 0,
"Mg" = 0,
"pH" = 9.91,
"pe" = 4,
"O2g" = 10,
"Calcite" = 2.07e-4,
"Dolomite" = 0
)
grid <- list(
n_cells = n,
s_cells = n,
type = types[1],
init_cell = as.data.frame(init_cell),
props = names(init_cell),
init_script = NULL
)
##################################################################
## Section 2 ##
## Diffusion parameters and boundary conditions ##
##################################################################
init_diffu <- c(
"C" = 1.2279E-4,
"Ca" = 1.2279E-4,
"Cl" = 0,
"Mg" = 0,
"pe" = 4,
"pH" = 7
)
alpha_diffu <- c(
"C" = 1E-4,
"Ca" = 1E-4,
"Cl" = 1E-4,
"Mg" = 1E-4,
"pe" = 1E-4,
"pH" = 1E-4
)
vecinj_diffu <- list(
list(
"C" = 0,
"Ca" = 0,
"Cl" = 0.002,
"Mg" = 0.001,
"pe" = 4,
"pH" = 9.91
)
)
boundary <- list(
"E" = 0,
"W" = 1
)
diffu_list <- names(alpha_diffu)
diffusion <- list(
init = init_diffu,
vecinj = do.call(rbind.data.frame, vecinj_diffu),
vecinj_index = boundary,
alpha = alpha_diffu
)
##################################################################
## Section 3 ##
## Phreeqc simulation ##
##################################################################
db <- RPhreeFile(system.file("extdata", "phreeqc_kin.dat",
package = "RedModRphree"
), is.db = TRUE)
phreeqc::phrLoadDatabaseString(db)
# NOTE: This won't be needed in the future either. Could also be done in a. pqi
# file
base <- c(
"SOLUTION 1",
"units mol/kgw",
"temp 25.0",
"water 1",
"pH 9.91 charge",
"pe 4.0",
"C 1.2279E-04",
"Ca 1.2279E-04",
"Mg 0.001",
"Cl 0.002",
"PURE 1",
"O2g -0.1675 10",
"KINETICS 1",
"-steps 100",
"-step_divide 100",
"-bad_step_max 2000",
"Calcite", "-m 0.000207",
"-parms 0.0032",
"Dolomite",
"-m 0.0",
"-parms 0.00032",
"END"
)
selout <- c(
"SELECTED_OUTPUT", "-high_precision true", "-reset false",
"-time", "-soln", "-temperature true", "-water true",
"-pH", "-pe", "-totals C Ca Cl Mg",
"-kinetic_reactants Calcite Dolomite", "-equilibrium O2g"
)
# Needed when using DHT
signif_vector <- c(7, 7, 7, 7, 7, 7, 7, 5, 5)
prop_type <- c("act", "act", "act", "act", "logact", "logact", "ignore", "act", "act")
prop <- names(init_cell)
iterations <- 500
setup <- list(
# bound = myboundmat,
base = base,
first = selout,
# initsim = initsim,
# Cf = 1,
grid = grid,
diffusion = diffusion,
prop = prop,
immobile = c(7, 8, 9),
kin = c(8, 9),
ann = list(O2g = -0.1675),
# phase = "FLUX1",
# density = "DEN1",
reduce = FALSE,
# snapshots = demodir, ## directory where we will read MUFITS SUM files
# gridfile = paste0(demodir, "/d2ascii.run.GRID.SUM")
# init = init,
# vecinj = vecinj,
# cinj = c(0,1),
# boundary = boundary,
# injections = FALSE,
iterations = iterations,
timesteps = rep(10, iterations)
)