I have, in the past few months, been studying numerical relativity, specifically the problem of spherical collapse of a scalar field as studied by Choptuik. I have also taken a look at some books that seem to cover extensive discussion on coordinates selection and singularity treatments (as well as the ADM formalism itself), such as

1 - 3+1 Formalism and Bases of Numerical Relativity, by Éric Gourgoulhon [gr-qc/0703035];

2 - Numerical Relativity: Solving Einstein's Equations on the Computer, by T. Baumgarte & S. Shapiro;

and many chapters of proceedings that have been written by Choptuik himself.

Still, I found that my greatest difficulty is with the so called adaptive mesh refinement (AMR) algorithm--this is the process of changing grid size on the fly, i.e. you perform some error estimation during the calculation and, if needed, you increase the number of grid pointz on a specific portion of the domain. This can (and should) be done recursively on regions that require a lot of resolution to produce meaningful results.

While it is quite intuitive to understand what is being done and why, I find it extremely difficult to apply the ideas in practice (i.e. sitting down and writing a code that works properly).

Although I am currently able to reproduce most of Choptuik's results, my errors are still pretty large because of the lack of AMR (making it impossible to study, for example, critical phenomena). Are there any resources other than the "canonical reference" by Berger & Oliger that you guys could recommend?

  • $\begingroup$ +1 for a numerical question which fits the site and one I’d like an answer to! $\endgroup$ – JamalS Nov 22 '17 at 20:48
  • $\begingroup$ A lot of the issues boil down to how the information of the grid is stored in memory and can be dynamically allocated in a linked list in such a way that (internal) boundary information is preserved. $\endgroup$ – Kyle Kanos Nov 22 '17 at 23:41
  • $\begingroup$ One of the authors of AstroBEAR wrote a good paper in it in 2010 or so (albeit for a non-rel. MHD code), but can't think of his name (pretty sure it begins with a C, as unhelpful as that may be). If I think of it, I'll find it and pass it on. $\endgroup$ – Kyle Kanos Nov 22 '17 at 23:42
  • $\begingroup$ Found it: Carroll-Nellenback et al (2013) (arXiv preprint). It's highly specific to their chosen AMR implementation, but much of the discussion should be applicable in other models (well, grid-based models anyway). $\endgroup$ – Kyle Kanos Nov 22 '17 at 23:55
  • $\begingroup$ @KyleKanos That's a good point, and whenever thinking about memory one should be wary of the CPU cache, i.e. cache misses. $\endgroup$ – JamalS Nov 23 '17 at 8:18

Part of the problem with the question asked is that there isn't one type of AMR engine. The two most common I've seen are either grid-based (or patch-based) or block-based. There's also choices of time-stepping (e.g., W-cycle or all grids updated simultaneously to $t+{\rm d}t$) and refinement strategies (though any good code would allow for at least user-specified criteria). There's even static mesh refinement that could be useful in some circumstances.

In any event, I think a good starting point for learning the basics is a deep reading of the appendix in Esquivel et al (2010). The relevant section you want is A.4.1:

The adaptive mesh...consists of a series of “root” blocks with a predetermined (set at compile time) number of cells, $n_x × n_y$ . Each block can be subdivided in a binary fashion into four blocks with the same number of cells of the parent block. We will refer to these new blocks as “siblings,” which have twice the resolution of their parent block. The maximum number of levels of refinement allowed ($n_{levs}$) is set at compile time....The result is a tree structure (of blocks or patches, not individual cells): the root blocks branch into higher resolutions up to $n_{levs}$ . The equations are only solved in “leaf” blocks, that is, only on blocks that are not further refined.

The next section or two discuss more about AMR strategies and gridding techniques for MPI, but the key is the above. Basically, the whole idea in the above case is to have the whole grid already refined and only evolve the maximum level needed based on the physical scenario (e.g., strong shocks, jeans instability, etc). You'd need to track where refinement's needed in the region, but your basic algorithm would be like,

do nlev=max_level,1,-1
   if update_flag(nlev)
   end if
   checkCoarsenRefine(nlev, U(nlev), update_flag(nlev))
end do

So this is reasonably efficient in terms of coding to get the AMR resolution needed, as it really requires a new boolean flag & an extra loop over the levels.

Hopefully you can see that a more robust algorithm would involve a linked list with details of the grid (e.g., location on the main grid, parent & children grids (if any), neighbor boundary values, actual neighbors) needing to be stored in the object/structure you create. In this method, you'd create a new object/structure when refinement is needed and store it in a linked list of objects (probably a vector of linked lists so that you can have one list at level $n$ and another at level $n+1$, etc). Carroll-Nellenback et al (2013) outline details of the distributed tree algorithm needed for efficient implementation of AMR; I suspect the references therein (and probably references to it) would be useful to look over (I've not done so).

My experience has been in 2D MHD (ended up w/o AMR w/ my research), so my knowledge of GRHD/GRMHD is very limited. The best I could tell you is that it would likely be useful for you to take a look at other papers in the field about how their AMR engines are implemented. Part of it, though, is going to be reading between the lines to see what is being done vs what is written down.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.