diff --git a/doc/SpectralSim_NAAICE_tug.R b/doc/SpectralSim_NAAICE_tug.R index 5497c91..a3d7ae0 100644 --- a/doc/SpectralSim_NAAICE_tug.R +++ b/doc/SpectralSim_NAAICE_tug.R @@ -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 gen_seq(int from, int to) { + int vsize = to - from + 1; + std::vector vec(vsize); + for (int i = 0; i < vsize; i++) { + vec[i] = i+from; + } + return vec; +}") + +gen_seq(24, 49) + +cppFunction("static const std::vector gen_vec(int elems) { + std::vector vec(elems); + for (int i = 0; i < elems; i++) { + vec[i] = i; + } + return vec; +}") + + +gen_vec(20) diff --git a/eval/bench_defs.hpp.in b/eval/bench_defs.hpp.in index 175030e..3233ae1 100644 --- a/eval/bench_defs.hpp.in +++ b/eval/bench_defs.hpp.in @@ -41,14 +41,12 @@ static const std::vector gen_vec(int elems) { return vec; } -static const std::vector gen_vec_rev(int from, int to) { - int size = to - from+1; - std::vector vec(size); - for (int i = from; i < to; i++) { - vec[i] = i; - std::cout << "i: " << i << "\n"; +static const std::vector gen_seq(int from, int to) { + int vsize = to - from + 1; + std::vector 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 - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/src/main.cpp b/src/main.cpp index 35aa8ac..470573a 100644 --- a/src/main.cpp +++ b/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;