mirror of
https://git.gfz-potsdam.de/naaice/poet.git
synced 2025-12-15 12:28:22 +01:00
Added "bench" dir with dolo_diffu_inner
This commit is contained in:
Marco De Lucia
Max Luebke
parent
5ef228c27f
commit
57af37483f
51
bench/dolo_diffu_inner/Eval.R
Normal file
51
bench/dolo_diffu_inner/Eval.R
Normal file
@ -0,0 +1,51 @@
|
||||
## Time-stamp: "Last modified 2022-12-16 20:26:03 delucia"
|
||||
|
||||
source("../../../util/data_evaluation/RFun_Eval.R")
|
||||
|
||||
sd <- ReadRTSims("naaice_2d")
|
||||
|
||||
sd <- ReadRTSims("Sim2D")
|
||||
|
||||
|
||||
sd <- ReadRTSims("inner")
|
||||
|
||||
tim <- readRDS("inner/timings.rds")
|
||||
|
||||
|
||||
simtimes <- sapply(sd, "[","simtime")
|
||||
|
||||
## workhorse function to be used with package "animation"
|
||||
PlotAn <- function(tot, prop, grid, breaks) {
|
||||
for (step in seq(1, length(tot))) {
|
||||
snap <- tot[[step]]$C
|
||||
time <- tot[[step]]$simtime/3600/24
|
||||
ind <- match(prop, colnames(snap))
|
||||
Plot2DCellData(snap[,ind], grid=grid, contour=FALSE, breaks=breaks, nlevels=length(breaks), scale=TRUE, main=paste0(prop," after ", time, "days"))
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
options(width=110)
|
||||
library(viridis)
|
||||
|
||||
Plot2DCellData(sd$iter_050$C$Cl, nx=1/100, ny=1/100, contour = TRUE,
|
||||
nlevels = 12, palette = "heat.colors",
|
||||
rev.palette = TRUE, scale = TRUE, main="Cl")
|
||||
|
||||
Plot2DCellData(sd$iter_050$C$Dolomite, nx=100, ny=100, contour = FALSE,
|
||||
nlevels = 12, palette = "heat.colors",
|
||||
rev.palette = TRUE, scale = TRUE, )
|
||||
|
||||
cairo_pdf("naaice_inner_Dolo.pdf", width=8, height = 6, family="serif")
|
||||
Plot2DCellData(sd$iter_100$C$Dolomite, nx=100, ny=100, contour = FALSE,
|
||||
nlevels = 12, palette = "viridis",
|
||||
rev.palette = TRUE, scale = TRUE, plot.axes = FALSE,
|
||||
main="2D Diffusion - Dolomite after 2E+4 s (100 iterations)")
|
||||
dev.off()
|
||||
|
||||
cairo_pdf("naaice_inner_Mg.pdf", width=8, height = 6, family="serif")
|
||||
Plot2DCellData(sd$iter_100$C$Mg, nx=100, ny=100, contour = FALSE,
|
||||
nlevels = 12, palette = "terrain.colors",
|
||||
rev.palette = TRUE, scale = TRUE, plot.axes=FALSE,
|
||||
main="2D Diffusion - Mg after 2E+4 s (100 iterations)")
|
||||
dev.off()
|
||||
144
bench/dolo_diffu_inner/dolo_diffu_inner.R
Normal file
144
bench/dolo_diffu_inner/dolo_diffu_inner.R
Normal file
@ -0,0 +1,144 @@
|
||||
## Time-stamp: "Last modified 2023-01-10 13:51:40 delucia"
|
||||
|
||||
database <- normalizePath("../data/phreeqc_kin.dat")
|
||||
input_script <- normalizePath("./dol_inner.pqi")
|
||||
|
||||
#################################################################
|
||||
## Section 1 ##
|
||||
## Grid initialization ##
|
||||
#################################################################
|
||||
|
||||
n <- 100
|
||||
m <- 100
|
||||
|
||||
types <- c("scratch", "phreeqc", "rds")
|
||||
|
||||
init_cell <- list(
|
||||
"H" = 110.683,
|
||||
"O" = 55.3413,
|
||||
"Charge" = -5.0822e-19,
|
||||
"C" = 1.2279E-4,
|
||||
"Ca" = 1.2279E-4,
|
||||
"Cl" = 0,
|
||||
"Mg" = 0,
|
||||
"O2g" = 0.499957,
|
||||
"Calcite" = 2.07e-4,
|
||||
"Dolomite" = 0
|
||||
)
|
||||
|
||||
grid <- list(
|
||||
n_cells = c(n, m),
|
||||
s_cells = c(1, 1),
|
||||
type = types[1],
|
||||
init_cell = as.data.frame(init_cell),
|
||||
props = names(init_cell),
|
||||
database = database,
|
||||
input_script = input_script
|
||||
)
|
||||
|
||||
|
||||
##################################################################
|
||||
## Section 2 ##
|
||||
## Diffusion parameters and boundary conditions ##
|
||||
##################################################################
|
||||
|
||||
## initial conditions
|
||||
init_diffu <- c(
|
||||
"H" = 110.683,
|
||||
"O" = 55.3413,
|
||||
"Charge" = -5.0822e-19,
|
||||
"C" = 1.2279E-4,
|
||||
"Ca" = 1.2279E-4,
|
||||
"Cl" = 0,
|
||||
"Mg" = 0
|
||||
)
|
||||
|
||||
## diffusion coefficients
|
||||
alpha_diffu <- c(
|
||||
"H" = 1E-6,
|
||||
"O" = 1E-6,
|
||||
"Charge" = 1E-6,
|
||||
"C" = 1E-6,
|
||||
"Ca" = 1E-6,
|
||||
"Cl" = 1E-6,
|
||||
"Mg" = 1E-6
|
||||
)
|
||||
|
||||
## list of boundary conditions/inner nodes
|
||||
vecinj_diffu <- list(
|
||||
list(
|
||||
"H" = 110.683,
|
||||
"O" = 55.3413,
|
||||
"Charge" = 1.90431e-16,
|
||||
"C" = 0,
|
||||
"Ca" = 0,
|
||||
"Cl" = 0.002,
|
||||
"Mg" = 0.001
|
||||
),
|
||||
list(
|
||||
"H" = 110.683,
|
||||
"O" = 55.3413,
|
||||
"Charge" = 1.90431e-16,
|
||||
"C" = 0,
|
||||
"Ca" = 0.0,
|
||||
"Cl" = 0.004,
|
||||
"Mg" = 0.002
|
||||
)
|
||||
)
|
||||
|
||||
vecinj_inner <- list(
|
||||
l1 = c(1,20,20),
|
||||
l2 = c(2,80,80),
|
||||
l3 = c(2,60,80)
|
||||
)
|
||||
|
||||
boundary <- list(
|
||||
# "N" = c(1, rep(0, n-1)),
|
||||
"N" = rep(0, n),
|
||||
"E" = rep(0, n),
|
||||
"S" = rep(0, n),
|
||||
"W" = rep(0, n)
|
||||
)
|
||||
|
||||
diffu_list <- names(alpha_diffu)
|
||||
|
||||
diffusion <- list(
|
||||
init = init_diffu,
|
||||
vecinj = do.call(rbind.data.frame, vecinj_diffu),
|
||||
vecinj_inner = vecinj_inner,
|
||||
vecinj_index = boundary,
|
||||
alpha = alpha_diffu
|
||||
)
|
||||
|
||||
#################################################################
|
||||
## Section 3 ##
|
||||
## Chemistry module (Phreeqc) ##
|
||||
#################################################################
|
||||
|
||||
|
||||
## # 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)
|
||||
|
||||
chemistry <- list(
|
||||
database = database,
|
||||
input_script = input_script
|
||||
)
|
||||
|
||||
#################################################################
|
||||
## Section 4 ##
|
||||
## Putting all those things together ##
|
||||
#################################################################
|
||||
|
||||
|
||||
iterations <- 100
|
||||
dt <- 200
|
||||
|
||||
setup <- list(
|
||||
grid = grid,
|
||||
diffusion = diffusion,
|
||||
chemistry = chemistry,
|
||||
iterations = iterations,
|
||||
timesteps = rep(dt, iterations)
|
||||
)
|
||||
35
bench/dolo_diffu_inner/dolo_inner.R
Normal file
35
bench/dolo_diffu_inner/dolo_inner.R
Normal file
@ -0,0 +1,35 @@
|
||||
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
|
||||
|
||||
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.1675 10
|
||||
KINETICS 1
|
||||
Calcite
|
||||
-m 0.00020
|
||||
-parms 0.05
|
||||
-tol 1e-10
|
||||
Dolomite
|
||||
-m 0.0
|
||||
-parms 0.005
|
||||
-tol 1e-10
|
||||
END
|
||||
@ -1,6 +1,6 @@
|
||||
## Simple library of functions to assess and visualize the results of the coupled simulations
|
||||
|
||||
## Time-stamp: "Last modified 2020-02-04 23:21:37 delucia"
|
||||
## Time-stamp: "Last modified 2022-12-16 17:38:12 delucia"
|
||||
|
||||
require(RedModRphree)
|
||||
require(Rmufits) ## essentially for PlotCartCellData
|
||||
@ -178,7 +178,7 @@ ExportToParaview <- function(vtu, nameout, results) {
|
||||
## "breaks" for color coding of 2D simulations
|
||||
Plot2DCellData <- function (data, grid, nx, ny, contour = TRUE,
|
||||
nlevels = 12, breaks, palette = "heat.colors",
|
||||
rev.palette = TRUE, scale = TRUE, ...) {
|
||||
rev.palette = TRUE, scale = TRUE, plot.axes=TRUE, ...) {
|
||||
if (!missing(grid)) {
|
||||
xc <- unique(sort(grid$cell$XCOORD))
|
||||
yc <- unique(sort(grid$cell$YCOORD))
|
||||
@ -207,11 +207,14 @@ Plot2DCellData <- function (data, grid, nx, ny, contour = TRUE,
|
||||
1))
|
||||
}
|
||||
par(las = 1, mar = c(5, 5, 3, 1))
|
||||
image(xc, yc, pp, xlab = "X [m]", ylab = "Y[m]", las = 1,
|
||||
asp = 1, breaks = breaks, col = colors, axes = FALSE,
|
||||
...)
|
||||
axis(1)
|
||||
axis(2)
|
||||
image(xc, yc, pp, xlab = "X [m]", ylab = "Y[m]", las = 1, asp = 1,
|
||||
breaks = breaks, col = colors, axes = FALSE, ann=plot.axes,
|
||||
...)
|
||||
|
||||
if (plot.axes) {
|
||||
axis(1)
|
||||
axis(2)
|
||||
}
|
||||
if (contour)
|
||||
contour(unique(sort(xc)), unique(sort(yc)), pp, breaks = breaks,
|
||||
add = TRUE)
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user