From 608926ea0297284db9a5240d3b6f9c3236cfe53b Mon Sep 17 00:00:00 2001 From: Marco De Lucia Date: Wed, 10 Apr 2024 23:15:25 +0200 Subject: [PATCH] doc/SpectralSim --- doc/SpectralSim_NAAICE_tug.R | 324 +++++++++++++++++++++++++++++++++++ 1 file changed, 324 insertions(+) create mode 100644 doc/SpectralSim_NAAICE_tug.R diff --git a/doc/SpectralSim_NAAICE_tug.R b/doc/SpectralSim_NAAICE_tug.R new file mode 100644 index 0000000..5497c91 --- /dev/null +++ b/doc/SpectralSim_NAAICE_tug.R @@ -0,0 +1,324 @@ +## Time-stamp: "Last modified 2024-04-10 22:38:14 delucia" + +## Compile the fortran code. On Linux, if a fortran compiler is +## installed, it should work out of the box +here <- getwd() +setwd("/home/work/simR/Rphree/SpecSim/") +system("R CMD SHLIB SS_Sim.f SS_Vario.f -o SSlib.so") +dyn.load("SSlib.so") + +## Load the functions +source("SS_Fun_sim.R") +source("SS_Fun_vario.R") +setwd(here) + +library(lattice) +library(PoetUtils) + + +## generate a grid. At the moment it must be a cartesian square grid +## with number of cells a power of two, e.g. 128, 256, ... +n <- 256 ## we want a 256x256 grid + +## Simple structure describing the simulation grid +gr <- SS.Grid(nodes=n, side=100) + +## fix the seed of RNG for reproducibility +set.seed(64276) + +## Extract the random numbers (uniform distribution between 0 and 1) +## once and for all +randoms <- runif(n**2) + +## Anisotropic SpectralSimulations with 2 different isotropic correlation lengths +## using the same random numbers. Here "range" is in meters! +a1 <- SS.Sim(vec.rn = randoms, nugget=0., sill=1, type="spherical", range=20, ratio=8, angle = -30, grid=gr, normalize=TRUE) + +## Compute the experimental variograms ONLY in X and Y direction. +## "lag" refers to the nodes on one side of the grid, not meters! +variog1 <- SS.Vario(a1, grid=gr, lag=2, nstep=72) + +## Visualize + +# cairo_pdf("Spherical_40_80.pdf", height=9, width=9) +par(mfrow=c(1,2)) +SS.Plot(a1, grid=gr, main="Range 20 m, anisotropy ratio 4", axes=FALSE) +axis(1, at=c(-100, 0, 100)) +axis(2, at=c(-100, 0, 100)) +SS.VarioPlot(vario=variog1, main="Experimental Variograms", lwd=2) +lines(variog1$htrue, SS.Spherical(variog1$htrue, range=20, sill=1), col="blue") +lines(variog1$htrue, SS.Spherical(variog1$htrue, range=20/4, sill=1), col="red") + + +sim1 <- matrix(a1, n, n) +levelplot(sim1) + +cut1 <- sim1[1:200, 1:200] + +levelplot(cut1) +sd(cut1) +mean(cut1) + +ax <- 1E-6*(cut1 - min(cut1))/max(cut1)+1E-7 +ay <- 1E-7*(cut1 - min(cut1))/max(cut1)+1E-7 + +sd(ax) +mean(ax) + +sd(ay) +mean(ay) + +range(ax) +range(ay) + + +PlotField2 <- function (data, grid, nx, ny, contour = TRUE, nlevels = 12, breaks, + palette = "heat.colors", rev.palette = TRUE, scale = TRUE, + plot.axes = TRUE, ...) +{ + if (!missing(grid)) { + xc <- unique(sort(grid$cell$XCOORD)) + yc <- unique(sort(grid$cell$YCOORD)) + nx <- length(xc) + ny <- length(yc) + if (!length(data) == nx * ny) + stop("Wrong nx, ny or grid") + } + else { + xc <- seq(1, nx) + yc <- seq(1, ny) + } + z <- matrix(round(data, 6), ncol = nx, nrow = ny, byrow = TRUE) + pp <- t(z[rev(seq(1, nrow(z))), ]) + if (missing(breaks)) { + breaks <- pretty(data, n = nlevels) + } + breakslen <- length(breaks) + colors <- do.call(palette, list(n = breakslen - 1)) + if (rev.palette) + colors <- rev(colors) + if (scale) { + par(mfrow = c(1, 2)) + nf <- layout(matrix(c(1, 2), 1, 2, byrow = TRUE), widths = c(4, + 1)) + } + par(las = 1, mar = c(5, 5, 3, 1)) + image(xc, yc, pp, 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) + if (scale) { + par(las = 1, mar = c(5, 1, 5, 5)) + PlotImageScale(data, breaks = breaks, add.axis = FALSE, + axis.pos = 4, col = colors) + axis(4, at = breaks) + } + invisible(pp) +} + +PlotField2(log10(ax), nx=200, ny=200, contour = FALSE, xaxs="i", yaxs="i", nlevels=6, plot.axes = FALSE) + +title(expression(log[10](alpha[x]))) + +levelplot(log10(t(ax))) + +levelplot(log10(ay)) + +data.table::fwrite(ax, "alpha_x.csv", col.names = FALSE) +data.table::fwrite(ay, "alpha_y.csv", col.names = FALSE) + +init <- c(110.01240000000783,55.50867875567298,-1.2162968960604385e-9,4.455329999159673e-7,2.0000000000005577e-12,0.0006151616266928334,0.0006147160946983908) + +names <- c("H","O","Charge","Ba","Cl","S_6_","Sr") + +inmat <- matrix(rep(init, 200*200), nrow = 200*200, byrow = TRUE) + +colnames(inmat) <- names + +data.table::fwrite(inmat, "barite_200_init.csv", col.names = TRUE) + +library(PoetUtils) + +out <- fread("../../build/barite_200_output.csv") + +cairo_pdf("barite_200_field_alphax.pdf", width = 9, height = 6, family="serif") +PlotField2(log10(ax), nx=200, ny=200, contour = FALSE, nlevels=10, + palette = "heat.colors", rev.palette = FALSE, + plot.axes = FALSE, main=expression(log[10](Sr))) +dev.off() + +cairo_pdf("barite_200_field_Ba.pdf", width = 9, height = 6, family="serif") +PlotField2(log10(out$Ba), nx=200, ny=200, contour = FALSE, nlevels=10, + palette = "terrain.colors", plot.axes = FALSE, main=expression(log[10](Sr))) +dev.off() + +PlotField2(log10(out$H), nx=200, ny=200, contour = FALSE, palette = "terrain.colors", nlevels=10, plot.axes = FALSE, main=expression(log[10](Sr))) + +PlotField2(out$S_6_, nx=200, ny=200, contour = FALSE, nlevels=10, plot.axes = FALSE, main=expression(log[10](Ba))) + +###################### barite large + +alphamat <- matrix(1e-6, nrow = 1000, ncol = 1000) +data.table::fwrite(alphamat, "barite_large_alpha.csv", col.names = FALSE) + + +init <- c(110.0124, 55.508678, -1.216296896e-9, 4.45533e-7, 0, 0.0006151616, 0.0006147161) +nams <- c("H","O","Charge","Ba","Cl","S_6_","Sr") +names(init) <- nams + +n <- 1000 +inmat <- matrix(rep(init, 1000*1000-n), nrow = 1000*1000 - n, byrow = TRUE) +colnames(inmat) <- nams + +bounds <- c(111.0124, 55.50622, -3.0E-07, 1, 2, 0.01, 0.001) +names(bounds) <- nams +bmat <- matrix(rep(bounds, n), nrow = n, byrow = TRUE) +colnames(bmat) <- nams + +set.seed(9342) +inds <- sample(seq_len(1E6), size=1E6) + +dim(inmat) +dim(bmat) + +totmat <- rbind(bmat, inmat) +dim(totmat) + +fmat <- totmat[inds, ] + + +xc <- rep(seq(1, 1000), times = 1000) +yc <- rep(seq(1, 1000), each = 1000) +iplot <- which(fmat[,"Cl"]==2) + +plot(xc[iplot], 1000-yc[iplot], pch=4, cex=0.5, xaxs="i",yaxs="i") + + +cairo_pdf("barite_large_init_locs.pdf", width = 6, height = 6, family="serif") +par(mar=c(1,1,0.5,0.5)) +plot(xc[iplot], yc[iplot], pch=4, cex=0.5, las=1, xaxs="i", yaxs="i", + xlab="Easting", ylab="Northig", asp=1, axes = FALSE) +box() +dev.off() + +data.table::fwrite(fmat, "barite_large_init.csv", col.names = TRUE) + + +out <- fread("../../build/barite_large_output.csv") + +x11() + +cairo_pdf("barite_large_field_Ba.pdf", width = 8, height = 6, family="serif") +##par(mar=c(4,4,0.5,0.5)) +PlotField2(log10(out$Ba), nx=1000, ny=1000, contour = FALSE, nlevels=10, + palette = "terrain.colors", plot.axes = FALSE) +box() +dev.off() + + +init + +round(rbind(init, bounds), 3) + + + +##### surfex, grid 500x250 + +nx <- 200 +ny <- 100 + +alphamat <- matrix(1.1e-12, nrow = ny, ncol = nx) +data.table::fwrite(alphamat, "surfex_alpha.csv", col.names = FALSE) + +a <- read.table( + textConnection( +"H 1.11e+02 120.0, +O 5.55e+01 55.1, +Charge -2.0e-13 8.0e-17, +C(-4) 2.0e-16 2.0e-15, +C(4) 2.0e-03 0.2, +Ca 2.0e-01 0.03, +Cl 3.0e-01 0.5, +Fe(2) 1.4e-04 0.0002, +Fe(3) 1.3e-09 2.0e-08, +H(0) 6.0e-12 2.0e-11, +K 2.0e-03 1.0e-05, +Mg 1.0e-02 0.2, +Na 2.0e-01 0.3, +S(-2) 5.9e-10 0, +S(2) 8.3e-15 8.3e-12, +S(4) 2.1e-14 5.1e-14, +S(6) 1.6e-02 0.026, +Sr 4.5e-04 0.045, +U(4) 2.5e-09 2.5e-08, +U(5) 1.6e-10 1.6e-10, +U(6) 2.3e-07 1.0e-05" +)) +vals <- t(a[, 2:3]) + +nams <- c("H", "O", "Charge", "CH4", "C", "Ca", "Cl", "Fe2", "Fe3", + "H0", "K", "Mg", "Na", "HS2", "S2", "S4", "S6","Sr", "U4", + "U5", "U6") +colnames(vals) <- nams +rownames(vals) <- c("IC", "BC") +vals + +length(nams) + +dim(vals) +t(vals) + +ic <- matrix(rep(vals[1,], nx*ny), nrow = nx*ny, byrow = TRUE) + +colnames(ic) <- nams + +data.table::fwrite(ic, "surfex_init.csv", col.names = TRUE) + + +dput(vals[2,]) + +out <- fread("../../build/surfex_output.csv") + +x11() + +cairo_pdf("surfex_field_U6.pdf", width = 10, height = 6, family="serif") +out <- fread("../../build/surfex_output.csv") +PlotField2(log10(out$U6), nx=nx, ny=ny, contour = FALSE, nlevels=10, + palette = "cm.colors", plot.axes = FALSE, rev.palette = FALSE) +dev.off() + +cairo_pdf("surfex_field_Na.pdf", width = 10, height = 6, family="serif") +out <- fread("../../build/surfex_output.csv") +PlotField2(out$Na, nx=nx, ny=ny, contour = FALSE, nlevels=10, + palette = "topo.colors", plot.axes = FALSE, rev.palette = TRUE) +dev.off() + + + +PlotField2(log10(out$U4), nx=nx, ny=ny, contour = FALSE, nlevels=10, + palette = "terrain.colors", plot.axes = FALSE) +out$U4[seq(50*200, 51*200)] + + + + +PlotField2(out$Na, nx=nx, ny=ny, contour = FALSE, nlevels=10, + palette = "terrain.colors", plot.axes = FALSE) +out$Na[seq(50*200, 51*200)] +out$Cl[seq(50*200, 51*200)] + + +PlotField2(log10(out$Fe2), nx=nx, ny=ny, contour = FALSE, nlevels=12, + palette = "terrain.colors", plot.axes = FALSE) + +PlotField2((out$Cl), nx=nx, ny=ny, contour = FALSE, nlevels=10, + palette = "terrain.colors", plot.axes = FALSE) + + +dev.off() +