Friday, August 13, 2010

Week Aug 9 - 13

On Monday I learned fwrap (excellent piece of software btw), there were a few minor technical issues, that I communicated with Kurt on the fwrap mailinglist and I also send him a simple patch, so that it works fine for my fortran code (functions returning the value itself, instead of a tuple of length 1).

Then I took my old fortran based shooting method solvers that I wrote couple years ago and wrapped them using fwrap and run couple simulations against my FE solver.

On Tuesday we had a lunch with all the advisors and students and the llnl director and a little presentation about what we did.

On Wednesday I run shooting method calculations for 50 states of silver, both for selfconsistent DFT potential and Z/r potential. I then also run the FE solver for the same DFT potential and compared results. There are lots of small technical issues, for example I had to use cubic splines to interpolate the potential, play with the mesh for the shooting method and so on.

However, the shooting method and FE agrees to every single printed digit, after making sure that the mesh is ok for both methods. For all potentials that I tried. That's very cool.

In the process of it, I also wrote a patch to SymPy to calculate exact energies for the Hydrogen atom, both from Schroedinger and Dirac equations. I still need to polish it a bit.

On Thursday I run couple more calculations and setup a poster and had a poster session, it was two hours, and I think around 7 people (not counting other students and people from our group) stopped by and talked with me about it, so I was very happy. Being able to solve radial Schroedinger and especially Dirac equations robustly is something that several people in the lab would really need.

Today I talked little bit (finally) about some Green functions in QM and QFT with a postdoc in the Quantum Simulations group, that I always wanted to, but didn't have time before, then packed my things and went back to Reno.

My plan for the next week(s) is to wrap up what I did and put it into articles. I already have enough material for some articles, so it has to be done. In parallel, I'd like to finish the FE Dirac solver, the coding is done, but now I need to play with adaptivity and also investigate if we are getting the spurious states, that other people are getting when using b-splines.

Friday, August 6, 2010

Week Aug 2 - 6

This week I essentially only worked on my LLNL poster, which I finally finished about two hours ago. I have created a web page for it:

where you can download pdf, sources, I also put there some relevant info and links.

It turned out to be a lot more work than I expected (well, as usual), but we were very thorough with John and in the process I discovered several bugs in my program, so I am glad we did it. I used to generate all the plots by hand, by manually adjusting all the parameters in the Python script (like atomic number, mesh parameters, element orders, adaptivity parameters, error tolerance and so on). Essentially I had to remember all these parameters for each of the plots (about 10 of them). Then I settled to have a Python dictionary, that holds all the parameters, and then I just pass them to a radial_schroedinger_equation_adapt(params, error_tol=1e-8) function.

Here are example of the parameters:

params_hydrogen_p_L = dict(l=0, Z=1, a=0, b=100, el_num=4, el_order=1,
eig_num=3, mesh_uniform=False, mesh_par1=20, adapt_type="p",
params_hydrogen_p_U = dict(l=0, Z=1, a=0, b=100, el_num=4, el_order=2,
eig_num=3, mesh_uniform=True, adapt_type="p", eqn_type="R")
params_hydrogen_hp_U = dict(l=0, Z=1, a=0, b=100, el_num=4, el_order=2,
eig_num=3, mesh_uniform=True, adapt_type="hp", eqn_type="R")
params_hydrogen_h_U = dict(l=0, Z=1, a=0, b=100, el_num=4, el_order=6,
eig_num=3, mesh_uniform=True, adapt_type="romanowski",

params_silver_p_L = dict(l=0, Z=47, a=0, b=150, el_num=4, el_order=13,
eig_num=50, mesh_uniform=False, mesh_par1=35, adapt_type="p",
params_silver_hp_L = dict(l=0, Z=47, a=0, b=150, el_num=4, el_order=13,
eig_num=50, mesh_uniform=False, mesh_par1=35, adapt_type="hp",

I mean, this is kind of obvious if you think about it, but for some reason I didn't do that at all at the beginning, because I thought --- I'll just run this once and I am done with it. But I had to run it like 20x, e.g. regenerating he plots, then creating a table about meshes, then redoing the table after changing the error tolerance, and so on.

Besides that I also got permission to release my code, so I'll go over it in the coming days and generate nice patches against Hermes1D.

Also in the process of creating the poster, I played a lot with p-FEM, uniform-p-FEM, hp-FEM and h-FEM and I will keep playing with that. It's clear to me now, that our current Hermes1D is not optimal. Especially the convergence of hp-FEM and p-FEM (as it is implemented right now) greatly depends on the initial mesh.

Nevertheless, even with the above limitations, hp-FEM seems to be really good if you don't a-priori know anything about the problem/mesh. One should not make any deep conclusions in 1D (it might be a bit different in 2D and 3D, and also I only did couple test problems), but from my experience so far, hp-FEM is a really good choice, if you just want to solve the problem and get a decent convergence (way better than h-FEM, and in general about the same as uniform-p-FEM with optimized mesh).

Another conclusion is that uniform-p-FEM (also called spectral element method), if you optimize the mesh for the problem, is very fast. All you have to do is increase the polynomial order and it goes very straight on the convergence graphs, it's very hard to beat. Also, and that I would like to write in the coming days, the algorithm for optimizing the mesh is really simple: just solve it with high "p", then play with the mesh parameters (for logarithmic mesh, there are only 2 parameters --- number of elements, and a ratio of first vs. last element), so that the eigenvalues (that one is interested in) are converged (with given accuracy) and optimize it wrt DOFs. The algorithm can also "look" at the convergence graphs and make sure it's steep enough. For atomic problems, my experience shows that the logarithmic mesh is good enough (as long as you optimize it). The advantage is that you do this once, and then (for close enough potentials in the Schroedinger equation), you just increase "p", and it's very robust and fast (no need for reference mesh, or trial refinements and so on).

When I get back to Reno, we'll do more research on hp-FEM with Pavel and I think this is not the last word to say. We need to review how we choose the candidates for eigenvectors, especially "p" vs "hp" and make it more robust. We'll see.