Mesh Mess

FEM, CFD, CAE, and all that.

Month: March, 2014

Moving to ggplot2 from gnuplot

Finally decided to move to ggplot2, after years’ use of gnuplot. I love gnuplot’s simple and straightforward, though its syntax could be a pain sometimes. All the plots in my previous latex docs are generated from gnuplot, using its nice latex terminal for publication (this is still unmatched by ggplot2, about this in a minute).

One thing that motivates me to use ggplot is Hadley’s book. The idea of grammar of graphics looks very sweet. Another reason is that by moving in R ecosystem I can use sweave to produce dynamic docs for some previosu docs (new docs will be written in org-mode so that’s not necessary). As to this,  another option is to use org-mode (babel) + matplotlib + latex. I love the elegant idea of writing whole thing in org,  and there are plenty examples to check out, see here and here. On the other hand, sweave is just that simple but has been having some long-time problems (some are addressed by knitr). As to matplotlib, it does the work alright, but in my opinion the plots are not as nice looking as those from ggplot.

However, gnuplot and matplotlib have one great advantage against ggplot, which is latex text rendering in plots. Latex rendering is supported by activating the latex option in matplotlib, similar to using latex terminal in gnuplot. While in ggplot, there is no such thing out of box. Workaround includes using tikzdevice (seems the project is in stall) or extrafont, but neither is satisfactory to me. That’s why org-mode gives a straightforward (though not satisfactory) solution by allowing me use gnuplot when needed, at least for a while, because better solution appears.

POD snippet in R

I’ve been working on proper orthogonal decomposition (POD) for a while. Though it’s nice to have packages such as modred out there, I cannot help myself writing a snippet using R, since lots of my previosu data processing was in that language. Without all the theory chitchat  like K-L transformation, from the the implementation-on-discrete-data perspective the (snapshot) POD method boils down to just singular value decomposition (SVD), and it’s out of box in R, which makes the task mostly straightforward. That’s why a POD function in R essentially calls svd in R (and it further calls LAPACK routines DGESDD and ZGESVD.):

pod <- function(arr, nmod=dim(arr)[1], nv=nmod){
  pod <- svd(arr, nmod, nv)
  return(pod)
}

Here “nmod” is the number of POD modes to be retrieved. Though it’s default to be the length of each snapshot, it’s rarely necessary this large. Mostly of the time a reasonable correlated set of snapshots requires less than 10 modes to obtain most of the system energy.

With pod of data “arr” obtained, we can do several thing immediately. We  can obtain the modes by retrieving the right singular vector:

pod_u<- function(pod, nnod){
  u <- matrix(unlist(pod["u"]), nrow = nnod )
  return(u)
}

or get the singular values:

pod_sv<- function(pod){
  sv <- as.vector(unlist(pod["d"]))
  return(sv)
}

or the coefficients of the pod modes for the approximation:

pod_coef<- function(pod, nshot, nmod){
  v <- matrix(unlist(pod["v"]), nrow = nshot )
  sv <- pod_sv(pod)[1:nmod]
  coef <- diag(sv) %*% t(v)
  return(coef)
}

and eventually the approximation itself

pod_approx<- function(pod, nshot, nmod, nnod){
  u <- pod_u(pod, nnod)
  coef <- pod_coef(pod, nshot, nmod)
  approx <- u %*% coef
  return(approx)
}