Compare commits
1 Commits
| Author | SHA1 | Date | |
|---|---|---|---|
|
fcfa45169d |
@ -1,8 +1,8 @@
|
||||
## Time-stamp: "Last modified 2024-04-10 22:38:14 delucia"
|
||||
## Time-stamp: "Last modified 2024-04-11 15:28:06 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")
|
||||
@ -10,11 +10,60 @@ dyn.load("SSlib.so")
|
||||
## Load the functions
|
||||
source("SS_Fun_sim.R")
|
||||
source("SS_Fun_vario.R")
|
||||
setwd(here)
|
||||
setwd("/home/work/simR/Rphree/NAAICE_tug_input/build")
|
||||
|
||||
library(lattice)
|
||||
library(PoetUtils)
|
||||
|
||||
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)
|
||||
}
|
||||
|
||||
|
||||
## 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, ...
|
||||
@ -72,54 +121,6 @@ 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)
|
||||
|
||||
@ -322,3 +323,125 @@ PlotField2((out$Cl), nx=nx, ny=ny, contour = FALSE, nlevels=10,
|
||||
|
||||
dev.off()
|
||||
|
||||
##### debug, grid 200x100
|
||||
|
||||
nx <- 200
|
||||
ny <- 100
|
||||
|
||||
alphamat <- matrix(1.1e-11, nrow = ny, ncol = nx)
|
||||
data.table::fwrite(alphamat, "debug_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, "debug_init.csv", col.names = TRUE)
|
||||
|
||||
|
||||
dput(vals[2,])
|
||||
|
||||
out <- fread("debug_output.csv")
|
||||
|
||||
x11()
|
||||
|
||||
## cairo_pdf("debug_field_U6.pdf", width = 10, height = 6, family="serif")
|
||||
# out <- fread("./debug_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("debug_field_Na.pdf", width = 10, height = 6, family="serif")
|
||||
out2 <- fread("./debug_output.csv")
|
||||
PlotField2(out2$Na, nx=nx, ny=ny, contour = FALSE, nlevels=10,
|
||||
palette = "topo.colors", plot.axes = FALSE, rev.palette = TRUE)
|
||||
## dev.off()
|
||||
out2$Na[seq(1, 200)]
|
||||
|
||||
range(out2$Na)
|
||||
|
||||
|
||||
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(out2$Fe3, nx=nx, ny=ny, contour = FALSE, nlevels=10,
|
||||
palette = "terrain.colors", plot.axes = FALSE)
|
||||
|
||||
out2$Na[seq(20*200, 21*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()
|
||||
|
||||
library(Rcpp)
|
||||
|
||||
cppFunction("static const std::vector<int> gen_seq(int from, int to) {
|
||||
int vsize = to - from + 1;
|
||||
std::vector<int> vec(vsize);
|
||||
for (int i = 0; i < vsize; i++) {
|
||||
vec[i] = i+from;
|
||||
}
|
||||
return vec;
|
||||
}")
|
||||
|
||||
gen_seq(24, 49)
|
||||
|
||||
cppFunction("static const std::vector<int> gen_vec(int elems) {
|
||||
std::vector<int> vec(elems);
|
||||
for (int i = 0; i < elems; i++) {
|
||||
vec[i] = i;
|
||||
}
|
||||
return vec;
|
||||
}")
|
||||
|
||||
|
||||
gen_vec(20)
|
||||
|
||||
@ -41,14 +41,12 @@ static const std::vector<int> gen_vec(int elems) {
|
||||
return vec;
|
||||
}
|
||||
|
||||
static const std::vector<int> gen_vec_rev(int from, int to) {
|
||||
int size = to - from+1;
|
||||
std::vector<int> vec(size);
|
||||
for (int i = from; i < to; i++) {
|
||||
vec[i] = i;
|
||||
std::cout << "i: " << i << "\n";
|
||||
static const std::vector<int> gen_seq(int from, int to) {
|
||||
int vsize = to - from + 1;
|
||||
std::vector<int> vec(vsize);
|
||||
for (int i = 0; i < vsize; i++) {
|
||||
vec[i] = i+from;
|
||||
}
|
||||
|
||||
return vec;
|
||||
}
|
||||
|
||||
@ -87,6 +85,7 @@ static bench_input barite_large_input = {
|
||||
.iterations = 5
|
||||
};
|
||||
|
||||
|
||||
static bench_input surfex_input = {
|
||||
.csv_file_init = "@SURFEX_BENCH@/surfex_init.csv",
|
||||
.csv_alpha_x = "@SURFEX_BENCH@/surfex_alpha.csv",
|
||||
@ -124,30 +123,41 @@ static bench_input surfex_input = {
|
||||
.iterations = 20
|
||||
};
|
||||
|
||||
static bench_input debug_input = {
|
||||
.csv_file_init = "@SURFEX_BENCH@/surfex_init.csv",
|
||||
.csv_alpha_x = "@SURFEX_BENCH@/surfex_alpha.csv",
|
||||
.csv_alpha_y = "@SURFEX_BENCH@/surfex_alpha.csv",
|
||||
.ncols = 200,
|
||||
.nrows = 100,
|
||||
.s_x = 0.02,
|
||||
.s_y = 0.01,
|
||||
.boundary = {.north_const = gen_vec(100),
|
||||
.south_const = gen_seq(1, 19),
|
||||
.east_const = {},
|
||||
.west_const = gen_seq(50, 99),
|
||||
.values = {120.0, // H
|
||||
55.1, // O
|
||||
8.0e-17, // Charge
|
||||
2.0e-15, // CH4
|
||||
0.2, // C
|
||||
0.03, // Ca
|
||||
0.5, // Cl
|
||||
0.0002, // Fe2
|
||||
2.0e-08, // Fe3
|
||||
2.0e-11, // H0
|
||||
1.0e-05, // K
|
||||
0.2, // Mg
|
||||
0.3, // Na
|
||||
0, // HS2
|
||||
8.3e-12, // S2
|
||||
5.1e-14, // S4
|
||||
0.026, // S6
|
||||
0.045, // Sr
|
||||
2.5e-08, // U4
|
||||
1.6e-10, // U5
|
||||
1.0e-05}}, // U6
|
||||
.timestep = 86400,
|
||||
.iterations = 20
|
||||
};
|
||||
|
||||
#endif // _BENCH_DEFS_HPP
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
43
src/main.cpp
43
src/main.cpp
@ -23,29 +23,40 @@ int main(int argc, char **argv) {
|
||||
}
|
||||
|
||||
if (with_output) {
|
||||
std::cout << "BARITE_200 started with output... ";
|
||||
timings[0] = run_bench(barite_200_input, "barite_200_output.csv");
|
||||
std::cout << "Runtime: " << timings[0] << " [ms]\n";
|
||||
// std::cout << "BARITE_200 started with output... ";
|
||||
// timings[0] = run_bench(barite_200_input, "barite_200_output.csv");
|
||||
// std::cout << "Runtime: " << timings[0] << " [ms]\n";
|
||||
|
||||
std::cout << "BARITE_LARGE started with output... ";
|
||||
timings[1] = run_bench(barite_large_input, "barite_large_output.csv");
|
||||
std::cout << "Runtime: " << timings[1] << " [ms]\n";
|
||||
// std::cout << "BARITE_LARGE started with output... ";
|
||||
// timings[1] = run_bench(barite_large_input, "barite_large_output.csv");
|
||||
// std::cout << "Runtime: " << timings[1] << " [ms]\n";
|
||||
|
||||
std::cout << "SURFEX started with output... ";
|
||||
timings[2] = run_bench(surfex_input, "surfex_output.csv");
|
||||
// std::cout << "SURFEX started with output... ";
|
||||
// timings[2] = run_bench(surfex_input, "surfex_output.csv");
|
||||
// std::cout << "Runtime: " << timings[2] << " [ms]\n";
|
||||
|
||||
std::cout << "DEBUG started with output... ";
|
||||
timings[2] = run_bench(debug_input, "debug_output.csv");
|
||||
std::cout << "Runtime: " << timings[2] << " [ms]\n";
|
||||
|
||||
} else {
|
||||
std::cout << "BARITE_200 started with no output... ";
|
||||
timings[0] = run_bench(barite_200_input);
|
||||
std::cout << " Runtime: " << timings[0] << " [ms]\n";
|
||||
|
||||
std::cout << "BARITE_LARGE started with no output... ";
|
||||
timings[1] = run_bench(barite_large_input);
|
||||
std::cout << "Runtime: " << timings[1] << " [ms]\n";
|
||||
// std::cout << "BARITE_200 started with no output... ";
|
||||
// timings[0] = run_bench(barite_200_input);
|
||||
// std::cout << " Runtime: " << timings[0] << " [ms]\n";
|
||||
|
||||
std::cout << "SURFEX started with no output... ";
|
||||
timings[2] = run_bench(surfex_input);
|
||||
// std::cout << "BARITE_LARGE started with no output... ";
|
||||
// timings[1] = run_bench(barite_large_input);
|
||||
// std::cout << "Runtime: " << timings[1] << " [ms]\n";
|
||||
|
||||
// std::cout << "SURFEX started with no output... ";
|
||||
// timings[2] = run_bench(surfex_input);
|
||||
// std::cout << "Runtime: " << timings[2] << " [ms]\n";
|
||||
|
||||
std::cout << "DEBUG started with no output... ";
|
||||
timings[2] = run_bench(debug_input);
|
||||
std::cout << "Runtime: " << timings[2] << " [ms]\n";
|
||||
|
||||
}
|
||||
|
||||
return EXIT_SUCCESS;
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user