commit 4cb75df8cdcd34f47801c908914a4f226f04d506
parent 5d5b96bf3b9c7754242aeb1f15be1a8077317822
Author: Ivan Gankevich <igankevich@ya.ru>
Date: Mon, 29 May 2017 11:14:09 +0300
Insert verification section.
Diffstat:
| R/common.R | | | 168 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| R/transform.R | | | 23 | +++++++++++++++++++++++ |
| R/velocity-potentials.R | | | 105 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| R/waves.R | | | 69 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| arma.org | | | 114 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
5 files changed, 479 insertions(+), 0 deletions(-)
diff --git a/R/common.R b/R/common.R
@@ -0,0 +1,168 @@
+source(file.path("R", "waves.R"))
+source(file.path("R", "transform.R"))
+
+arma.qqplot_grid <- function (dir, params, titles, ...) {
+ wave_params <- arma.load_wave_parameters(dir, params)
+ i <- 1
+ for (name in names(wave_params)) {
+ arma.qqplot(wave_params[[name]], 100, titles[[i]], ...)
+ i <- i + 1
+ }
+}
+
+arma.wavy_plot <- function (data, t, ...) {
+ slice <- data[data$t == t,]
+ x <- unique(slice$x)
+ y <- unique(slice$y)
+ z <- with(slice, {
+ n <- sqrt(length(z))
+ out <- matrix(nrow=n, ncol=n)
+ out[cbind(x, y)] <- z
+ out
+ })
+ nrz <- nrow(z)
+ ncz <- ncol(z)
+ # Create a function interpolating colors in the range of specified colors
+ jet.colors <- colorRampPalette( c("blue", "green") )
+ # Generate the desired number of colors from this palette
+ nbcol <- 100
+ color <- jet.colors(nbcol)
+ # Compute the z-value at the facet centres
+ zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
+ # Recode facet z-values into color indices
+ facetcol <- cut(zfacet, nbcol)
+ persp(x, y, z, phi=30, theta=30, col=color[facetcol], ...)
+}
+
+arma.skew_normal_1_plot <- function(x, params) {
+ data <- mapply(
+ function (s, k) arma.skew_normal_1(x, s, k),
+ params$skewness,
+ params$kurtosis
+ )
+ plot.new()
+ plot.window(xlim=range(x),ylim=range(data))
+ axis(1); axis(2); box()
+ for (i in seq_len(ncol(data))) {
+ d <- data[,i]
+ lines(x, d, lty=paste(params$linetypes[[i]]))
+ }
+ title(xlab="x", ylab="y")
+}
+
+
+arma.skew_normal_2_plot <- function(x, params) {
+ data <- mapply(
+ function (a) arma.skew_normal_2(x, a),
+ params$alpha
+ )
+ plot.new()
+ plot.window(xlim=range(x),ylim=range(data))
+ axis(1); axis(2); box()
+ for (i in seq_len(ncol(data))) {
+ d <- data[,i]
+ lines(x, d, lty=paste(params$linetypes[[i]]))
+ }
+ title(xlab="x", ylab="y")
+}
+
+arma.fmt <- function(x, ndigits) {
+ format(round(x, ndigits), nsmall=ndigits)
+}
+
+arma.plot_partitions <- function() {
+ library("rgl")
+ part_alpha <- 0.2
+ part_col <- "grey"
+ part_size <- 2
+ sc <- 0.8
+ off <- part_size * (1-sc)/2
+ bcol <- "grey"
+ balpha <- 1.0
+ sc2 <- 1.0 - sc
+ off2 <- part_size * sc/2
+ ccol <- "red"
+ calpha <- 1.0
+ # whole parts
+ a1 <- cube3d(color=part_col, alpha=part_alpha)
+ a2 <- cube3d(color=part_col, alpha=part_alpha)
+ a3 <- cube3d(color=part_col, alpha=part_alpha)
+ shade3d(translate3d(a1, 0*part_size, 0, 0))
+ shade3d(translate3d(a2, 1*part_size, 0, 0))
+ shade3d(translate3d(a3, 2*part_size, 0, 0))
+ # stripped parts
+ b1 <- scale3d(cube3d(color=bcol, alpha=balpha), sc, sc, sc)
+ b2 <- scale3d(cube3d(color=bcol, alpha=balpha), sc, sc, sc)
+ b3 <- scale3d(cube3d(color=bcol, alpha=balpha), sc, sc, sc)
+ shade3d(translate3d(b1, 0 + off, off, off))
+ shade3d(translate3d(b2, 2 + off, off, off))
+ shade3d(translate3d(b3, 4 + off, off, off))
+ # overlap intervals
+ c1 <- scale3d(cube3d(color=ccol, alpha=calpha), sc2, 1, 1)
+ c2 <- scale3d(cube3d(color=ccol, alpha=calpha), sc2, 1, 1)
+ c3 <- scale3d(cube3d(color=ccol, alpha=calpha), sc2, 1, 1)
+ shade3d(translate3d(c1, 0 + off2, 0, 0))
+ shade3d(translate3d(c2, 2 + off2, 0, 0))
+ shade3d(translate3d(c3, 4 + off2, 0, 0))
+}
+
+arma.plot_ramp_up_interval <- function(label="Ramp-up interval") {
+ zeta <- read.csv(file.path("build", "standing_wave", "zeta.csv"))
+ t <- round(mean(zeta$t))
+ res <- arma.wavy_plot(zeta, t, scale=FALSE)
+ library("grDevices")
+ ax <- 7
+ ay <- 7
+ my <- max(zeta$y)
+ lines(trans3d(
+ c(0, 0, ax, ax, 0),
+ c(0, my, my, 0, 0),
+ c(0, 0, 0, 0, 0),
+ pmat=res
+ ), col="red", lwd=3)
+ text(trans3d(0, my/2, max(zeta$z)*1.5, pmat=res), label, col="red", font=2)
+ from <- trans3d(0, my/2, max(zeta$z)*1.4, pmat=res)
+ to <- trans3d(0, my/2, max(zeta$z)*0.05, pmat=res)
+ arrows(from$x, from$y, to$x, to$y, lwd=2, angle=10, length=0.1, col="red")
+}
+
+arma.plot_factory_vs_openmp <- function(...) {
+ args <- list(...)
+ perf <- read.csv(file.path("data", "performance", "factory-vs-openmp.csv"))
+ scale <- 10 ** args$power
+ x <- perf$nt * perf$nx * perf$ny / scale
+ plot.new()
+ plot.window(xlim=range(x),ylim=range(perf[c("openmp", "factory")]))
+ pts <- pretty(x)
+ axis(1, at=pts, labels=sapply(pts, function(x) {as.expression(bquote(.(x) %.% 10 ** .(args$power)))}))
+ axis(2)
+ box()
+ lines(x, perf$openmp, lty="solid")
+ lines(x, perf$factory, lty="dashed")
+ title(xlab=args$xlab, ylab=args$ylab)
+}
+
+arma.plot_factory_vs_openmp_overlap <- function(...) {
+ args <- list(...)
+ openmp <- read.csv(file.path("data", "performance", "overlap-openmp.csv"), na.strings="")
+ factory <- read.csv(file.path("data", "performance", "overlap-factory.csv"), na.strings="")
+ openmp$t <- (openmp$t - min(openmp$t)) / args$scale
+ factory$t <- (factory$t - min(factory$t)) / args$scale
+ plot.new()
+ plot.window(xlim=range(c(factory$t, openmp$t)),ylim=range(0, 5))
+ axis(1)
+ axis(2, at=c(1, 3), labels=args$labels, las=1, hadj=1)
+ # OpenMP
+ lines(openmp$t, rep.int(3, length(openmp$t)))
+ openmp_pts <- openmp[!is.na(openmp$mark),]
+ openmp_y <- rep.int(3, length(openmp_pts$t))
+ points(openmp_pts$t, openmp_y)
+ text(openmp_pts$t, openmp_y, labels=openmp_pts$mark, pos=c(3, 3, 1, 1))
+ # Factory
+ lines(factory$t, rep.int(1, length(factory$t)))
+ factory_pts <- factory[!is.na(factory$mark),]
+ factory_y <- rep.int(1, length(factory_pts$t))
+ points(factory_pts$t, factory_y)
+ text(factory_pts$t, factory_y, labels=factory_pts$mark, pos=c(3, 1, 3, 1))
+ title(xlab=args$xlab)
+}
diff --git a/R/transform.R b/R/transform.R
@@ -0,0 +1,23 @@
+arma.skew_normal_1 <- function(z, skewness, kurtosis) {
+ (24 + 4*z*(-3 + z**2)* skewness + (3 - 6*z**2 + z**4)* kurtosis)/(24*(exp((z**2)/2)*sqrt(2*pi)))
+}
+
+arma.bits.erf <- function(x) {
+ 2 * pnorm(x * sqrt(2)) - 1
+}
+
+arma.bits.erfc <- function(x) {
+ 2 * pnorm(x * sqrt(2), lower = FALSE)
+}
+
+arma.bits.skewness_2 <- function(alpha) {
+ (sqrt(2) * (4 - pi) * (alpha**3)) / (sqrt(pi + (pi-2) * (alpha**2))**3)
+}
+
+arma.bits.kurtosis_2 <- function(alpha) {
+ (8*(-3 + pi)*(alpha**4))/ ((pi + (-2 + pi)*(alpha**2))**2)
+}
+
+arma.skew_normal_2 <- function(z, alpha) {
+ arma.bits.erfc( -((z*alpha)/sqrt(2)) ) / (exp((z**2)/2)*sqrt(2*pi))
+}
diff --git a/R/velocity-potentials.R b/R/velocity-potentials.R
@@ -0,0 +1,105 @@
+arma.middle_element <- function (x) {
+ x[round(length(x) + 0.5)/2]
+}
+
+arma.plot_velocity_potential_field <- function (dir, ...) {
+ args <- list(...)
+ phi <- read.csv(file.path(dir, 'phi.csv'))
+ left_top_x <- 0.1
+ right_top_x <- max(phi$x)
+ # slice time and Y ranges through the center
+ slice_t <- arma.middle_element(unique(phi$t))
+ slice_y <- arma.middle_element(unique(phi$y))
+ print(paste('Middle elements (TY) = ', slice_t, slice_y))
+ phi_slice <- phi[phi$t == slice_t & phi$y == slice_y & phi$x >= left_top_x & phi$z >= -8,]
+ x <- unique(phi_slice$x)
+ z <- unique(phi_slice$z)
+ left_top_z <- max(phi_slice$z)
+ right_top_z <- left_top_z
+ print(paste('Velocity field size (XZ) = ', length(x), length(z)))
+
+ # convert data frame to matrix
+ seq_x <- seq_along(x)
+ seq_z <- seq_along(z)
+ indices <- data.matrix(expand.grid(seq_x, seq_z))
+ u <- with(phi_slice, {
+ out <- matrix(nrow=length(seq_x), ncol=length(seq_z))
+ out[indices] <- phi
+ out
+ })
+
+ # get wave profile
+ zeta <- read.csv(file.path(dir, 'zeta.csv'))
+ zeta_slice <- zeta[zeta$t == slice_t & zeta$y == slice_y & zeta$x >= left_top_x,]
+
+ plot.new()
+ plot.window(xlim=range(x),ylim=range(z),asp=1)
+ axis(1); axis(2); box()
+
+ .filled.contour(
+ x, z, u,
+ levels=args$levels,
+ col=args$col
+ )
+
+
+ contour(
+ x, z, u,
+ levels=args$levels,
+ asp=1,
+ drawlabels=TRUE,
+ add=TRUE
+ )
+
+ top_area_x <- c(left_top_x*0.99, zeta_slice$x, right_top_x*1.01)
+ top_area_z <- c(left_top_z*1.10, zeta_slice$z, right_top_z*1.10)
+ polygon(top_area_x, top_area_z, lwd=4, border=NA, col='white')
+ lines(zeta_slice$x, zeta_slice$z, lwd=4)
+ box()
+ title(xlab="x", ylab="z")
+}
+
+arma.plot_velocity_potential_field_legend <- function (...) {
+ args <- list(...)
+ levels <- args$levels
+ col <- args$col
+ plot.new()
+ plot.window(
+ xlim = c(0, 1),
+ ylim = range(levels),
+ xaxs = "i",
+ yaxs = "i"
+ )
+ rect(0, levels[-length(levels)], 1, levels[-1L], col = col)
+ axis(4)
+ box()
+}
+
+arma.plot_velocity <- function(file1, file2, ...) {
+ args <- list(...)
+ data1 <- read.table(file1)
+ data2 <- read.table(file2)
+ ylim <- range(c(data1, data2))
+ if ("ylim" %in% names(args)) {
+ ylim <- args$ylim
+ }
+ plot.new()
+ plot.window(
+ xlim = c(0, nrow(data1)-1),
+ ylim = ylim
+ )
+ lines(data1, lty=args$linetypes[[1]])
+ lines(data2, lty=args$linetypes[[2]])
+ axis(1)
+ axis(2)
+ box()
+ title(xlab="x",ylab="u(x)")
+ legend(
+ "topleft",
+ c(
+ as.expression(bquote(u[1](x))),
+ as.expression(bquote(u[2](x)))
+ ),
+ lty = paste(args$linetypes)
+ )
+}
diff --git a/R/waves.R b/R/waves.R
@@ -0,0 +1,69 @@
+#!/usr/bin/Rscript
+
+arma.bits.qweibull <- function (x, param, sh) {
+ sh <- 2.0
+ sc = param$mean / gamma(1.0 + 1.0 / sh)
+ qweibull(x, shape=sh, scale=sc)
+}
+arma.bits.qnorm <- function (x, param) {
+ qnorm(x, mean=param$mean, sd=param$sd)
+}
+
+arma.QFUNCTIONS <- list(
+ elevation = arma.bits.qnorm,
+ heights_x = function (x, param) { arma.bits.qweibull(x, param, 2.0) },
+ heights_y = function (x, param) { arma.bits.qweibull(x, param, 2.0) },
+ lengths_x = function (x, param) { arma.bits.qweibull(x, param, 2.3) },
+ lengths_y = function (x, param) { arma.bits.qweibull(x, param, 2.3) },
+ periods = function (x, param) { arma.bits.qweibull(x, param, 3.0) }
+)
+
+arma.ALL_WAVE_CHARACTERISTICS <- list(
+ "elevation",
+ "heights_x",
+ "heights_y",
+ "lengths_x",
+ "lengths_y",
+ "periods"
+)
+
+arma.load_wave_parameters <- function (
+ dir = "output",
+ wave_params = arma.ALL_WAVE_CHARACTERISTICS,
+ qfuncs = arma.QFUNCTIONS
+) {
+ sapply(wave_params, function (filename) {
+ param <- list()
+ param$filename = filename
+ param$data <- read.table(file.path(dir, filename))[[1]]
+ param$mean <- mean(param$data)
+ param$sd <- sd(param$data)
+ param$qfunc <- function (x) { qfuncs[[filename]](x, param) }
+ param$min <- param$qfunc(.01)
+ param$max <- param$qfunc(.99)
+ elem <- list()
+ elem[[filename]] <- param
+ elem
+ })
+}
+
+arma.qqplot <- function (param, nsamples=100, ttl, ...) {
+ qdata <- param$qfunc(ppoints(nsamples))
+ qqplot(
+ qdata,
+ param$data,
+ asp=1,
+ xlim=c(param$min, param$max),
+ ylim=c(param$min, param$max),
+ ...
+ )
+ title(ttl, line=-2)
+ qqline(param$data, distribution=param$qfunc)
+}
+
+#wave_params <- arma.load_wave_parameters()
+#par(pty="s", mfrow=c(3, 3))
+#for (name in names(wave_params)) {
+# print(name)
+# arma.qqplot(wave_params[[name]])
+#}
diff --git a/arma.org b/arma.org
@@ -5,6 +5,7 @@
#+LATEX_CLASS_OPTIONS:
#+LATEX_HEADER_EXTRA: \input{preamble}
#+OPTIONS: H:5 num:0 todo:nil toc:nil
+#+PROPERTY: header-args:R :results graphics :exports results
#+begin_abstract
Simulation of sea waves is a problem appearing in the framework of developing
@@ -542,6 +543,119 @@ steps.
trigonometric identities to shift the phase.
** Evaluation
+In\nbsp{}cite:degtyarev2011modelling,degtyarev2013synoptic,boukhanovsky1997thesis AR
+model the following items are verified experimentally:
+- probability distributions of different wave characteristics (wave heights,
+ lengths, crests, periods, slopes, three-dimensionality),
+- dispersion relation,
+- retention of integral characteristics for mixed wave sea state.
+In this work both AR and MA model are verified by comparing probability
+distributions of different wave characteristics.
+
+*** Verification of wavy surface integral characteristics
+In\nbsp{}cite:рожков1990вероятностные the authors show that several ocean wave
+characteristics (listed in table\nbsp{}[[tab-weibull-shape]]) have Weibull distribution,
+and wavy surface elevation has Gaussian distribution. In order to verify that
+distributions corresponding to generated realisation are correct,
+quantile-quantile plots are used (plots where analytic quantile values are used
+for \(OX\) axis and estimated quantile values for \(OY\) axis). If the estimated
+distribution matches analytic then the graph has the form of the straight line.
+Tails of the graph may diverge from the straight line, because they can not be
+reliably estimated from the finite-size realisation. Different methods of
+extracting waves from realisation produce variations in quantile function tails,
+it is probably impractical to extract every possible wave from realisation since
+they may (and often) overlap.
+
+#+name: tab-weibull-shape
+#+caption: Values of Weibull shape parameter for different wave characteristics.
+#+attr_latex: :booktabs t
+| Characteristic | Weibull shape (\(k\)) |
+|----------------------+---------------------|
+| Wave height | 2 |
+| Wave length | 2.3 |
+| Crest length | 2.3 |
+| Wave period | 3 |
+| Wave slope | 2.5 |
+| Three-dimensionality | 2.5 |
+
+Verification was performed for standing and propagating waves. The corresponding
+ACFs and quantile-quantile plots of wave characteristics distributions are shown
+in
+fig.\nbsp{}[[propagating-wave-distributions]],\nbsp{}[[standing-wave-distributions]],\nbsp{}[[acf-slices]].
+
+#+name: propagating-wave-distributions
+#+begin_src R :file build/propagating-wave-qqplots.pdf
+source(file.path("R", "common.R"))
+par(pty="s", mfrow=c(2, 2))
+arma.qqplot_grid(
+ file.path("build", "propagating_wave"),
+ c("elevation", "heights_y", "lengths_y", "periods"),
+ c("elevation", "height Y", "length Y", "period"),
+ xlab="x",
+ ylab="y"
+)
+#+end_src
+
+#+caption: Quantile-quantile plots for propagating waves.
+#+label: propagating-wave-distributions
+#+RESULTS: propagating-wave-distributions
+[[file:build/propagating-wave-qqplots.pdf]]
+
+#+name: standing-wave-distributions
+#+begin_src R :file build/standing-wave-qqplots.pdf
+source(file.path("R", "common.R"))
+par(pty="s", mfrow=c(2, 2))
+arma.qqplot_grid(
+ file.path("build", "standing_wave"),
+ c("elevation", "heights_y", "lengths_y", "periods"),
+ c("elevation", "height Y", "length Y", "period"),
+ xlab="x",
+ ylab="y"
+)
+#+end_src
+
+#+caption: Quantile-quantile plots for standing waves.
+#+label: standing-wave-distributions
+#+RESULTS: standing-wave-distributions
+[[file:build/standing-wave-qqplots.pdf]]
+
+#+name: acf-slices
+#+header: :width 6 :height 9
+#+begin_src R :file build/acf-slices.pdf
+source(file.path("R", "common.R"))
+propagating_acf <- read.csv(file.path("build", "propagating_wave", "acf.csv"))
+standing_acf <- read.csv(file.path("build", "standing_wave", "acf.csv"))
+par(mfrow=c(5, 2), mar=c(0,0,0,0))
+for (i in seq(0, 4)) {
+ arma.wavy_plot(standing_acf, i, zlim=c(-5,5))
+ arma.wavy_plot(propagating_acf, i, zlim=c(-5,5))
+}
+#+end_src
+
+#+caption: Time slices of ACF for standing (left column) and propagating waves (right column).
+#+label: acf-slices
+#+RESULTS: acf-slices
+[[file:build/acf-slices.pdf]]
+
+Graph tails in fig.\nbsp{}[[propagating-wave-distributions]] deviate from original
+distribution for individual wave characteristics, because every wave have to be
+extracted from the resulting wavy surface to measure its length, period and
+height. There is no algorithm that guarantees correct extraction of all waves,
+because they may and often overlap each other. Weibull distribution right tail
+represents infrequently occurring waves, so it deviates more than left tail.
+
+Correspondence rate for standing waves (fig.\nbsp{}[[standing-wave-distributions]])
+is lower for height and length, roughly the same for surface elevation and
+higher for wave period distribution tails. Lower correspondence degree for
+length and height may be attributed to the fact that Weibull distributions were
+obtained empirically for ocean waves which are typically propagating, and
+distributions may be different for standings waves. Higher correspondence degree
+for wave periods is attributed to the fact that wave periods of standing waves
+are extracted more precisely as the waves do not move outside simulated wavy
+surface region. The same correspondence degree for wave elevation is obtained,
+because this is the characteristic of the wavy surface (and corresponding AR or
+MA process) and is not affected by the type of waves.
+
** Discussion
* Determining wave pressures for discretely given wavy surface
** Two-dimensional case
