doc/SpectralSim
This commit is contained in:
Marco De Lucia
parent
1c3528a53d
commit
608926ea02
324
doc/SpectralSim_NAAICE_tug.R
Normal file
324
doc/SpectralSim_NAAICE_tug.R
Normal file
@ -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()
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user