Methods in Numerical Relativity

Question

I have been reading the book of Masaru Shibata Numerical Relativity to grasp some ideas on the methods used. I see that at the heart of the method the system of differential equations is "converted" in a set of algebraic equations, by the means of finite difference method and this works fairly well.

On the other hand, comparing with other areas of physics and engineering, I know from my own experience that the finite element method has many advantages over the finite difference method, and, in fact, that is one of the reasons why in many applications it is used much more.

I wonder why in numerical relativity this the case is not also. I suspect that it is the following difficulty, in two and three dimensions there are very good algorithms to divide the domain where the system of differential equations is solved in a partition or set of "finite elements", but that in 4D this collides with certain difficulties. I understand that if these practical issues were solved, simulations in numerical relativity could be done using the finite element method and take advantage of all its benefits.

I would like to know what people who know the subject think and see if they share my suspicion that the difficulty of constructing a proper mesh in 4D is part of the problem.

Answer 1

I see that at the heart of the method the system of differential equations is "converted" in a set of algebraic equations, by the means of finite difference method and this works fairly well.

I've not read the book, but most of the "more advanced" codes out there actually use finite volume instead of finite difference method. While they are similar in structure, based on discretizing space on a domain and using differences between those discrete cells, the finite volume method approximates an integral whereas the finite difference approximates the PDE.

I know from my own experience that the finite element method has many advantages over the finite difference method, and, in fact, that is one of the reasons why in many applications it is used much more.

Were any of this statement true, then no CFD code written since 1950 would have used FVM, it all would have been FEM. Since that is decidedly not true, then the whole statement is wrong.

There are some advantages of FEM, e.g. when studying stresses and strains of materials, but those disappear when you are considering a compressible, turbulent fluid where shocks may arise, which is fairly common in astrophysical systems. FVM is the better method for handling such scenarios.
Of course, this is not to say that FEM cannot handle such scenarios. Of course it can, since the two methods are fully equivalent. It is just that the mathematical work needed to make FEM handle such cases vs. what is needed for FVM is not generally worth the additional efforts from a software development perspective.

As to FEM being used "much more," this is almost certainly your own bias. My experience is quite the opposite, but I spent my (astro)physics career developing/maintaining FVM codes and never looked into FEM. I would naively assume roughly equal proportions of FEM vs FVM codes out there (though I could see a bias towards FVM due to many people choosing to use existing software (COMSOL, OpenFOAM, etc) rather than developing a new code).

I understand that if these practical issues were solved, simulations in numerical relativity could be done using the finite element method and take advantage of all its benefits.

This assumes there has been no research or developments into FEM GR, which is not true. See for instance arXiv:0509099, Meier (1999), Radice & Rezzolla 2011 (NB: PDF) or this dissertation (NB: PDF), etc.

I would like to know what people who know the subject think and see if they share my suspicion that the difficulty of constructing a proper mesh in 4D is part of the problem.

According to a former contributor to this site, the most challenging aspect of numerical relativity is converting the primitive variables (density, velocity, pressure) from the conserved variables (density, momentum, energy);¹ the difficulty being that the conserved variables are nonlinearly related to the primitive variables with no simple inverses. In addition to that, you have physical constraints on primitive variables (e.g., strictly positive densities & pressures, avoiding superluminal velocities) that may not cleanly be matched to the conserved variables (or cause problems with the Riemann solver)--I believe numerical root finding is required for each cell to find the correct conversion.

^{1. The conserved density differs from the primitive density by a factor of $\gamma$.}

Answer 2

Finite elements or actually, spectral methods (if you agree for this exchange), are indeed used in modern numerical relativity codes - see SPEC or SPECtre codes. I think what a lot of these things boil down to is the challenge of covering so many magnitudes of scale - and using appropriate mesh refinement techiques to handle it not only efficiently, but scalably!

The groups that develop codes that use spectral methods at some point encounter problems that arise in the mathematical formulation - in simulations of binary black holes, the singularities have to be excised from the computational domain, and there has to be some clever tracking of all the domains, the element cells, etc. If I'm not mistaken, successful dealing with this exact problem was one of the missing pieces in the code developed by the group of Frans Pretorius (BBH simulations).

On the other hand, in codes that use FDM, the handling of the singularities is done via the so-called 'moving puncture' method/gauge (Campanelli, Baker, around 2006), which derives from the interpretation and methods of the Schwarzschild solution in isotropic coordinates. The $r=0$ singularity is considered to be second asymptotic infinity. The conformal factor $\psi$ is then evolved as $W=\psi^{-2}$, $\chi=\psi^{-4}$ or $\phi=\ln\psi$ to avoid finite differencing of very large numbers around $r=0$. In this method of handling of the singularity, suitable gauge conditions (most commonly, the 1+log evolution equation for the lapse and Gamma driver evolution equation for the shift) are employed and no excision is necessary.

I think one struggle of the finite element method is its inability to capture hydrodynamical shocks accurately (or, perhaps, the difficulties associated with it) and preserve fluxes.

Finite volume methods are used for high-resolution shock capturing methods in general relativistic hydrodynamical codes, such as WhiskyTHC, or IllinoisGRMHD, which are both distributed as part of the open-source Einstein Toolkit.

For the question about the cell dimensionality - as a matter of fact, all the above codes will employ some version of the 3+1 decomposition, where you evolve fields on a 3-dimensional hypersurface. Therefore what you really require in FEM in these scenarios is a distribution of 3-dimensional cells.

Let me add that nowadays there are codes developed that aim to use a hybrid approach - for example spectral methods far-away, where there is either no matter and/or small curvature - and where the FD method is computationally less efficient - and FD in large curvature or highly dynamical regions, where you need greater finesse.

Answer 3

One reason is inertia. FDM has worked (and continues to work) for over 50 years, and a lot of development has been done in that time to optimize and build on the early codes. Is there any indication FEM will be better than what we already have? Not that I'm aware, so it would be a huge task to get back to the cutting edge starting from scratch with FEM, and there is no guarantee it will be any better. You risk reinventing the wheel and making it square

You could argue that spectral methods and FEM are the same thing, and spectral methods are quite popular. But any "new" method has to fight against the proven success of finite differences