Saturday, July 31, 2010

Week July 26 - 30

This week I have been wrapping up my work at LLNL and trying to generate some comparisons between different approaches to adapt to multiple eigenvectors at once.

It turns out that most of the issues are in the way we create the new mesh for the next adaptivity iteration, in particular, how do we choose which candidates to refine. I have tried to converge to the lowest eigenvector (as well as to any other eigenvector too), to the sum of the eigenvectors, to each eigenvector individually and taking the union of the meshes and so on. I have also implemented the uniform p-FEM as well as tried p-FEM and hp-FEM using the approaches I mentioned above.

It seems to be crucial to have a good initial mesh, at least for p-FEM. If I use a good mesh and p-adaptivity, I am able to get the best results so far. In principle hp adaptivity should be at least as good, but our current approach doesn't show it yet. Hopefully we'll manage to make it work.

Besides that I have also implemented H1 norms for the Function class (based on Fekete points) by calculating coefficients with regards to a FE basis and some other little things.

Friday, July 23, 2010

Week July 19 - 23

This week I have learned how projections work in detail and wrote it up here:

including all proofs (that orthogonal projection finds the closest vector and so on). At the end of the notes I have calculated some examples in 1D, so that one can see that indeed it doesn't depend on the basis and that the basis doesn't even have to be orthogonal. Then one has to use this approach to calculate the coefficients:

I have then implemented it in h1d. Due to the lack of time, I am now developing everything in my private branch. I'll obtain the permission in about 2 or 3 weeks, so then I'll push it into the master.

Besides that I have implemented Chebyshev points for orders greater than 48, for which I don't have exact Fekete points anymore (it'd be just a matter of running my sympy script longer, but I was hitting some accuracy issues when solving those large polynomials numerically -- one needs to obtain all roots, so Chebyshev points are ok for now). So I can now represent arbitrary polynomials in 1D.

I have implemented powering of the discrete function, it automatically determines which polynomial order it has to use and creates a new discrete function (the power) on the new mesh. I wrote lots of tests for that, and I hit an interesting bug, that my naive comparison code:

assert abs(x-y) < eps

was not good enough anymore for larger numbers and I had to read some documentation, and implement the following function:

@vectorize(0, 1)
def feq(a, b, max_relative_error=1e-12, max_absolute_error=1e-12):
a = float(a)
b = float(b)
# if the numbers are close enough (absolutely), then they are equal
if abs(a-b) < max_absolute_error:
return True
# if not, they can still be equal if their relative error is small
if abs(b) > abs(a):
relative_error = abs((a-b)/b)
relative_error = abs((a-b)/a)
return relative_error <= max_relative_error

Then I implemented the global H1 and L2 projections, so far the projected function is hardwired, I still need to allow the user to specify any discrete function to be projected. I need to precalculate it and so on.

I wrote bunch of tests for the projections and powers and I always discovered some bugs by writing more tests, so the progress is slow, but at least I can trust the code that is tested.

I also helped Pavel to fix couple segfaults, as well as some other things.

Monday, July 19, 2010

Theoretical Physics Reference Book

Today I fulfilled my old dream --- I just created my first book! Here is how it looks like:

More images

Here is the source code of the book:, the repository contains a branch 'master' with the code and 'gh-pages' with the generated html pages, that are hosted at github, at the url

Then I published the book at Lulu:, I wanted a hardcover book, so I setup a project at Lulu, used some Lulu templates for the cover and that was it. Lulu's price for the book is $19.20 (166 black & white pages, hardcover), then I can set my own price and the rest of the money probably goes to me. I set the price to $20, because Lulu has free shipping for orders $20 or more. You can also download the pdf (for free) at the above link (or just use my git repository). So far this didn't cost me anything.

I have then ordered the book myself (just like anybody else would, at the above address) and it arrived today. It's a regular hardcover book. Beautiful, you can browse the pictures above. It smells deliciously (that you have to believe me). And all that it cost me was $19.20.

As for the contents itself, you can browse it online at, essentially it's most of my physics notes, that I collected over the years. I'd like to treat books like software --- release early release often. This is my first release and I would call it beta. The main purpose of it was to see if everything goes through, how long it takes (the date inside the book is July 4, 2010, I created and ordered it on July 5, got the physical book on July 19) and what the quality is (excellent). I also wanted to see how the pages correspond to the pdf (you can see for yourself on the photos, click on the picasa link above).

Now I need to improve the first pages a bit, as well as the last pages, improve the index, write some foreword and so on. I also need to think how to better organize the contents itself and generally improve it. I also need to figure out some versioning scheme, so far this is version 0.1. I think I'll do edition 1, edition 2, edition 3, and so on. And whenever I feel that I have added enough new content, I'll just publish it as a new edition. So if you want to buy it, I suggest to wait for my 1.0 version, that will have the mentioned improvements.

It'd be also cool to have all the editions online somehow and create nice webpages for it (currently points directly to the book html itself).

So far the book is just text. I still need to figure out how to handle pictures and also whether or not to use program examples (in Python, using sympy, scipy, etc.). So far I am inclining not to put there any program codes, as then I don't need to maintain them.

Overall I am very pleased with the quality, up to some minor issues that I mentioned above, everything else end up just fine. I think we have come a long way from the discovery of the printing press. Anybody can now create a book for free, and if you want to hold the hardcopy in your hands, it costs around $20. You don't need to order certain amounts of books, nor partner with some publisher etc. I think that's just awesome.

Friday, July 16, 2010

Week June 12 -- July 16

This week I was implementing hp-adaptivity based on H1 projections first for the ground state and later for other states.

It took me the whole week to debug things and I am still not done. There were issues with normalizing vectors and other problems. I am now able to converge any individual vector and that seems to work fine, but I am still not able to converge multiple vectors somehow.

I spent also some time with speeding up my implementation of fekete points, and I have implemented the fast evaluation based on Lagrange interpolation polynomials. In general, I made for example the l2_norm() method about 10x faster. I think it's about as fast now as if it was written in C++ directly. The main loops are now optimized C without any Python C/API calls.

My plan is now to implement taking squares of the solutions (exactly), and so on, and projecting it back to the original mesh. I'll try converging to that. At the same time, I'll try to converge to multiple eigenvalues. Hermes1d needs some improvements, that take me quite some time to implement.

Saturday, July 10, 2010

Week June 5 -- July 9

This week I have learned how to use the new C++ support in Cython, and posted a demo project at github (, then I have refactored the Python wrappers in hermes1d and I think things are now way cleaner.

Then I have finished my reimplementation of Romanowski's algorithm for adapting eigenvalues for the radial Schroedinger equation, I've extracted his values from the graphs, calculated the same thing myself and it agrees perfectly.

I am now polishing the h1 adaptivity in h1d, I essentially reimplemented it myself, so that I can consider more candidates (all possible combinations of "p" and "h", and also soon I'll consider bisecting, trisecting and so on of the interval).

I've implemented hydrogen wavefunctions in sympy and use it to project onto a very fine mesh (12 order, lots of elements), convert to Fekete points and use it for the adaptivity. It's in pure Python, as I need to develop very fast to get some results and see what approach is the best. Everything is fast, except the selection of candidates, which I'll now rewrite into Cython and try to use hermes1d whenever possible.

Next week I'll try to merge my adaptivity with hermes1d adaptivity and make it fast. And see how it converges.

I've also spent about 10 hours with improving Pavel's hermes2d branch, as well as implementing the Vector class in hermes_common.

Saturday, July 3, 2010

Week June 28 -- July 2

On Saturday I flew to Prague, then to Pilsen to the ESCO2010 conference. I had to presentations there and I met lots of awesome people.

I discussed with Dmitri how to solve the Euler equations using FEM and I think I know how to implement it now in hermes2d. I will try to give it a shot at the end of the summer, as this is something that was bothering me a lot.

I have made a progress in my fekete points code, I also wrote analytic solutions for the hydrogen atom in sympy. I'll wrap it up over the weekend and push to h1d.

I have investigated how to use github for collaboration using git, here is how it looks like for sympy:

I have figured out how to host any domain at github, in particular this:

now runs at github. I think this is really awesome.

Besides that, I spent some time discussing how to do hp-adaptivity in h1d with Pavel and other people, as well as I helped couple people with hermes2d and so on.

My plan for the next week is to implement hp-adaptivity in h1d for eigenproblems, as well as finish my fekete machinery and make it fast using cython and use it in conjunction with the h1d adaptivity.