Adaptive Mesh Refinement

Introduction

This document briefly describes the AMR-related features in AGILE. The different options can be set in the meshlist part (see Setting parameters) of the agile.par file. For a more extensive description, you can read the article 'Parallel, grid-adaptive approaches for relativistic hydro and magnetohydrodynamics', R. Keppens, Z. Meliani, A.J. van Marle, P. Delmont, A. Vlasis, & B. van der Holst, 2011, JCP.* doi:10.1016/j.jcp.2011.01.020.

AGILE uses a standard block-based, octree AMR scheme, where we have blocks of user-controlled dimension (set by the block_nx1, block_nx2, block_nx3) in a hierarchically nested manner. To simplify the parallelization, we gave up flexibility to allow different sized refinement ratios between grid levels, fixing it to 2. Also, we now use the same time step for all levels. A generic skeleton code, generic enough to hold for any AMR code having similar restrictions, is shown below.

timeloop : do
   exit time loop by user-defined criteria
   compute timestep constraint dt over all grids
   conditionally save data
   advance all grids for one time step dt
      distinguish source- and flux-addition strategies
         select numerical scheme, geometric sources, ...
         collect fluxes at fine/coarse boundaries
      after every parallel update, on all processors
         fill all ghost cells for all grid blocks
      at the final temporal update, fix for conservation
   regrid
      quantify the error estimator on all grids
      build the new grid structure, ensure proper nesting
      load-balance the new grid structure
   update time counters
end do timeloop

Some more info follows on the different aspects involved.

Important (global) parameters

Some important global parameters are in the module src/mod_global_parameters.fpp. In particular, note that the maximum number of levels is the parameter nlevelshi. The latter is default set to 20. If your want to run with more levels, you need to change their value and recompile. The number of levels set in the par file as refine_max_level must always be smaller or equal to nlevelshi. This nlevelshi also returns in those parameters that are defined per level, such as limiter which needs to be set for all (default) 20 levels.

The parameter max_blocks controls how many blocks can be maximally in use per MPI task. It is set at runtime and can be steered by setting it in the meshlist of the parameter file.
Since AGILE pre-allocates a lot of its GPU datastructures, max_blocks is the knob to control AGILE's memory consumption. If you run out of memory, check how many blocks you expect and set max_blocks accordingly.

AMR criteria

This in essence describes the module src/amr/mod_errest.fpp, or at least its most essential aspects.

The block-tree nature implies that a decision for refining/coarsening is to be made on a block-by-block basis. This automated block-based regridding procedure involves 3 steps:

    (1) consider all blocks at level 1 < l < refine_max_level , with refine_max_level the maximal grid level selected;
    (2) quantify the local error $E_x$ at each gridpoint in a certain grid block;
    (3) if any point has this error exceeding a user-set threshold refine_threshold(l), refine this block (and ensure proper nesting);
    (4) if all points have their error below a user-set fraction of the threshold derefine_ratio(l) used in the previous step, coarsen the block (for l>1).

The local error estimator can be one of three options, selected by refine_criterion, each possibly augmented with user-defined criteria. For refine_criterion=1, only refinement based on usr_refine_grid is active. Any of the other 3 estimators use a user-selected subset of the conserved or auxiliary variables (or even variables that are computed dynamically at the time of regridding), through the formula $E_{\mathbf{x}} = \sum_{iw}\sigma_{iw}\,E^{\mathrm{Rel}}_{\mathbf{x},iw}$.

The indices included are user-identified with the w_refine_weight array, where the sum of weights of all variables should equal to 1. The estimated error is a weighted sum of contributions from all variables with non-zero weight.

For refine_criterion=3, we select a Lohner type [R. Lohner, An adaptive finite element scheme for transient problems in CFD, Comp. Meth. App. Mech. Eng. 61, 323 (1987)] prescription as also used in the PARAMESH library or the FLASH3 code. In our experience, it is computationally efficient as it employs only instantaneous values from t^n. It in essence discretizes a weighted second derivative in each grid point.

The Lohner prescription writes in formulae as

$$\begin{eqnarray} E^{\mathrm{Rel}}_{\mathbf{x},iw} &=& \sqrt{\frac{N_{iw}}{\max\!\left(D_{iw},\,\epsilon\right)}} \,,\\[2pt] N_{iw} &=& \sum_{i_1}\sum_{i_2}\left[\,\Delta_{i_1}\!\left(\Delta_{i_2} w_{iw}\right)\right]^2 \,,\\[2pt] D_{iw} &=& \sum_{i_1}\sum_{i_2}\left[\,S_{i_1}\!\left(\left|\Delta_{i_1} w_{iw}\right|\right) + f_{\mathrm{wave},1}\,S_{i_2}\!\left(S_{i_1}\!\left(\left|w_{iw}\right|\right)\right)\right]^2 \,, \end{eqnarray}$$

where $\Delta_i$ and $S_i$ mean a central difference and a sum along dimension $i$, $\epsilon$ a floor guarding the division, and $f_{\mathrm{wave},1}$ the wave-filter parameter. The wave filter is set per level, and defaults as amr_wavefilter(1:nlevelshi)=1.0d-2.

Data structures

The data structures are defined in src/mod_physicaldata.fpp and src/amr/mod_forest.fpp, you can inspect them for learning more details.

We provide details on useful data structures. All of these are suited for any curvilinear coordinate system, and merely reflect the tree structure of the block-AMR. We implicitly assume a fixed refinement ratio of two. Schematic figures for a 2D Cartesian case generalize straightforwardly to higher or lower dimensionality.

The overall domain is considered 'rectangular', i.e. bounded by the coordinate pairs xprobmin1,xprobmax1, xprobmin2,xprobmax2 and xprobmin3,xprobmax3. On the lowest grid level l=1, one controls the coarsest resolution as well as a suitable domain decomposition, by specifying both the total number of level 1 grid cells domain_nx1, domain_nx2, domain_nx3 along with the individual block size per dimension block_nx1, block_nx2, block_nx3, which exclude the ghost cells. The total cell number must be an integer multiple of the block size in each dimension, so e.g. domain_nx1 = 4*block_nx1.

A hypothetical 2D domain is shown below, which corresponds to a domain where 4 by 3 blocks on level 1 are exploited in this domain decomposition, and where local refinement was activated in the top-right corner: the two level-1 blocks (3,2) and (3,3) plus (4,2) and (4,3) are resolved at level l=2, and one of those level-2 blocks is refined once more to level l=3. Global, integer grid indices are introduced per dimension, in a manner where knowledge of these grid indices, combined with AMR level knowledge, instantly allows one to localize the grid when needed. Following the figure, the grid on level l=2 indicated by global grid indices (5,3) is indeed the fifth grid block horizontally, and the third vertically, when the domain would be resolved fully with level l=2 blocks. The total amount of grid blocks per dimension, per level l, is stored in ng1(l), ng2(l), ng3(l), and the actual length of a grid block on level l, per dimension, is dg1(l), dg2(l), dg3(l).

The figure below reflects the tree representation of the same hypothetical grid hierarchy, where the presence of a grid leaf at a certain grid level is identified through a boolean variable. As indicated before, the total number of active grid leafs nleafs may change from timestep to timestep. This tree info is stored in the structure tree_root(ig1,ig2,ig3), indexed by the level-l global grid indices, whose %node knows about the global grid index through %node%ig1, %node%ig2, %node%ig3, the level through %node%level, the processor on which it resides through the integer %node%ipe, and its presence or absence through the logical %node%leaf.

Various extra indices are helpful to traverse the tree structure. Local grid indices across AMR levels are schematically given below, which are used to identify the directional neighbours, as well as the children and parent blocks. These are used to realize and facilitate the possible interprocessor communication patterns.

The directional neighbours of a grid block, and the exchange that fills a block's ghost layer from them, are shown below for the 3D case.

For parallelization, we adopted a fairly straightforward Z-order or Morton- order space filling curve (SFC). The Morton curve threads the leaf blocks in a recursive Z pattern — descending into each refined block in the same Z order — which gives every leaf a single integer Morton number. The blocks are then handed to processors in contiguous arcs of that ordering. Load-balancing is done after every timestep, following the adaptive remeshing. When exploiting N_p GPUs, our strategy for load balancing merely ensures that each GPU has at least [nleafs/N_p]int (denoting integer division) grid blocks, while the remainder increase this number by 1 for the first as many GPUs. For the hypothetical hierarchy above — 27 leaf blocks on 4 GPUs — the first three GPUs hold 7 blocks each and the fourth holds 6. The grid Morton numbers of all grids residing on processor _mype lie between Morton_start(mype),Morton_stop(mype). The global grid index, once you know the grid Morton number Morton_no is found from sfc_to_igrid(Morton_no), which gives the relation between the SFC and the global grid index igrid. The data for the conservative variables for grid igrid is then actually found from pw(igrid)%w.