FEM, CFD, CAE, and all that.

### Offlineimap + mu4e

So I live in emacs most of the time: writing papers, coding, and composing emails. The markup from org mode and a latex front end such as auctex make it really uncessary to resort to any other writing/typesetting systems, let alone MS word, especially when I’m doing scientific writing. So after replacing a failed hard drive and upgrading to Ubuntu 14.04, one of the first things was to set up the email for emacs, so I can assure others that I am still kicking.

To do that is not difficult at all, by using offlineimap + mu4e. No I don’t use gnus and yes I’ve tried that and it didn’t work out. In this scenario, we use offlineimap as a mail manager and retriever, and use mu4e as the front-end of emacs to do the work. Setting up offlineimap is pretty straightforward, one needs only to set up a profile file “~/.offlineimaprc”. Setting up mu4e is pretty easy too, as everything should work out of box from a mu installation and a minimum emacs setup.

On top of these two, I needed to set up the “~/.authinfo” file so that I don’t have to type in password everytime checking mails. This can be done in several ways. I opt for a simple python script and emacs’s in-house EasyPG package. Basically what I did was make the offlineimap use the passphrase encrypted in the ~/.authinfo.gpg file.

And that’s it. Now I can fire up mu4e and manage emails through the message-mode. IMHO this is how emacs’ email management should be, and I’m not the only one.

### Soft shell structure under hydrodynamical load

In the middle of constructing a reduced-order model (ROM) of a highly flexible structure interacting with surround flow, I built a structure model and loaded it with oscillating pressure at the lower part. The animation from the simulation is as follows.

There are several challenges to this type of FEA problems, mostly are related to the highly flexibility nature of the material. The structure is first subject to a quasi-static hydrostatic load, which makes it buckle to a new equilibrium position. After that, ambient oscillating pressure with a much smaller magnitude than the initial buckling load is applied to make the fabric vibrate.

This is a perfect example for using explicit-implicit switching for time integration. Since it’s usually difficult for an implicit solver to find the bulking configuration, as a result of the fact that the Newton-Ralphson method does not gaarantee  convergence, one needs to issue an explicit analysis for the initial loading stage. Some solvers provide specific “bulking simulation” controls for this type of the problem. After the initial loading the structure to a new equilibrium (passing a bifurcation point, as one clearly see in the animation), we proceed to impose the oscillating load. This type of load is supposed to be of a much smaller magnitude than the initial load, so that the structural vibration would be a perturbation near the new equilibrium.

Since explicit method suffers from a general stability problem (notorious forward Euler blow-up in every toy example in every textbook), it is only applicable to short-duration problems such as impact/contact, and it is very easy for the blow-up to happen during a cyclic loading problem, due to change of sign of the derivatives. Thus we need to turn off the explicit simulation and replace it with an implicit one. The benefit? (Unconditional) stability and hence greater allowable time step. The cost, on the other hand, is much higher since now we are actually inverting the stiffness matrix and perform in-time iterations.

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

### Size of deallocated pointer in Fortran

Something I just found out. A deallocated pointer retains its size in terms of the return  of the SIZE function. Specifically, if the pointer is first allocated with size $n$ and later deallocated, SIZE function still returns $n$. I’ve checked this on GNU and PGI Fortran and got same results, but not sure whether this is compiler dependent. Check standards later.