diff --git a/data/SimDol1D_diffu.R b/data/SimDol1D_diffu.R index c30f5adb8..6217aec3d 100644 --- a/data/SimDol1D_diffu.R +++ b/data/SimDol1D_diffu.R @@ -1,28 +1,117 @@ -## chemical database -# NOTE: not needed now +################################################################# +## 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) -## only the directory -# demodir <- system.file("extdata", "demo_rtwithmufits", package="Rmufits") - -n <- 10 -m <- 10 - -prop <- c("C", "Ca", "Cl", "Mg", "pH", "pe", "O2g", "Calcite", "Dolomite") - -signif_vector <- c(7, 7, 7, 7, 7, 7, 7, 5, 5) -prop_type <- c("act", "act", "act", "act", "logact", "logact", "ignore", "act", "act") - -grid <- list( - n_cells = c(n), - s_cells = c(n) -) - - # NOTE: This won't be needed in the future either. Could also be done in a. pqi # file base <- c( @@ -56,116 +145,11 @@ selout <- c( "-pH", "-pe", "-totals C Ca Cl Mg", "-kinetic_reactants Calcite Dolomite", "-equilibrium O2g" ) -# -# 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") - - - -init_diffu <- c( - "C" = 1.2279E-4, - "Ca" = 1.2279E-4, - "Cl" = 0, - "Mg" = 0, - "pe" = 4, - "pH" = 7 -) - -alpha_diffu <- c( - "C" = 1E-6, - "Ca" = 1E-6, - "Cl" = 1E-6, - "Mg" = 1E-6, - "pe" = 1E-8, - "pH" = 1E-8 -) - -vecinj_diffu <- list( - list( - "C" = 0, - "Ca" = 0, - "Cl" = 0.002, - "Mg" = 0.001, - "pe" = 4, - "pH" = 7 - ) -) - -inner_index <- c(5) -inner_vecinj_index <- rep(1, 3) - -vecinj_inner <- cbind(inner_index, inner_vecinj_index) - -# constant_boundary <- list( -# n_cell = seq(1, n), -# i_vecinj = rep(1, n) -# ) -# -# closed_boundary <- list( -# n_cell = seq(1, n), -# i_vecinj = rep(0, n) -# ) - -boundary <- list( - "N" = rep(1, n), - "E" = rep(0, n), - "S" = rep(0, n), - "W" = rep(0, n) -) - -diffusion <- list( - init = init_diffu, - vecinj = do.call(rbind.data.frame, vecinj_diffu), - vecinj_inner = vecinj_inner, - boundary = boundary, - alpha = alpha_diffu -) - -# cbound <- seq(from = 0, to = n - 1) - -## setup boundary conditions for transport - we have already read the -## GRID with the following code: -## grid <- Rmufits::ReadGrid(paste0(demodir,"/d2ascii.run.GRID.SUM")) -## cbound <- which(grid$cell$ACTNUM == 2) -## dput(cbound) -# cbound <- c( -# 1L, 50L, 100L, 150L, 200L, 250L, 300L, 350L, 400L, 450L, 500L, -# 550L, 600L, 650L, 700L, 750L, 800L, 850L, 900L, 950L, 1000L, -# 1050L, 1100L, 1150L, 1200L, 1250L, 1300L, 1350L, 1400L, 1450L, -# 1500L, 1550L, 1600L, 1650L, 1700L, 1750L, 1800L, 1850L, 1900L, -# 1950L, 2000L, 2050L, 2100L, 2150L, 2200L, 2250L, 2300L, 2350L, -# 2400L, 2450L, 2451L, 2452L, 2453L, 2454L, 2455L, 2456L, 2457L, -# 2458L, 2459L, 2460L, 2461L, 2462L, 2463L, 2464L, 2465L, 2466L, -# 2467L, 2468L, 2469L, 2470L, 2471L, 2472L, 2473L, 2474L, 2475L, -# 2476L, 2477L, 2478L, 2479L, 2480L, 2481L, 2482L, 2483L, 2484L, -# 2485L, 2486L, 2487L, 2488L, 2489L, 2490L, 2491L, 2492L, 2493L, -# 2494L, 2495L, 2496L, 2497L, 2498L, 2499L, 2500L -# ) - -# boundinit <- matrix(rep(init[-c(7, 8)], length(cbound)), byrow = TRUE, nrow = length(cbound)) -# myboundmat <- cbind(cbound, boundinit) -# myboundmat[cbound == 1, c(2:7)] <- vecinj -# colnames(myboundmat) <- c("cbound", names(vecinj)) # TODO: dt and iterations +iterations <- 500 + setup <- list( # bound = myboundmat, base = base, @@ -174,7 +158,7 @@ setup <- list( # Cf = 1, grid = grid, diffusion = diffusion, - prop = prop, + prop = names(init_cell), immobile = c(7, 8, 9), kin = c(8, 9), ann = list(O2g = -0.1675), @@ -188,9 +172,10 @@ setup <- list( # cinj = c(0,1), # boundary = boundary, # injections = FALSE, - s_grid = c(n), - n_grid = c(n), - dt = 1, - iterations = 1, - timesteps = rep(1, 10) + iterations = iterations, + timesteps = rep(10, iterations) ) + +# not needed yet +# signif_vector <- c(7, 7, 7, 7, 7, 7, 7, 5, 5) +# prop_type <- c("act", "act", "act", "act", "logact", "logact", "ignore", "act", "act") diff --git a/data/SimDol2D_diffu.R b/data/SimDol2D_diffu.R index aff12a472..8e12c234d 100644 --- a/data/SimDol2D_diffu.R +++ b/data/SimDol2D_diffu.R @@ -156,7 +156,7 @@ selout <- c( # TODO: dt and iterations -iterations <- 10 +iterations <- 500 setup <- list( # bound = myboundmat, @@ -180,11 +180,8 @@ setup <- list( # cinj = c(0,1), # boundary = boundary, # injections = FALSE, - s_grid = c(n, m), - n_grid = c(n, m), - dt = 1, iterations = iterations, - timesteps = rep(1, iterations) + timesteps = rep(10, iterations) ) # not needed yet