POD snippet in R

by Yi Zhang

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)
}

 

Advertisements