SMR and AMR

Static Mesh Refinement

Static Mesh Refinement means that some regions of the computational domain are covered by finer grids but they are fixed and not adaptive. This is not as flexible as AMR, but it is faster (and easier to analyze) than AMR and is often useful enough in problems in astrophysics. We recommend using SMR instead of AMR if the regions of interest are known in advance. For the details of the grid structure, see Static Mesh Refinement.

Let us go back to the blast wave test in Cartesian coordinates again, and say we want to resolve the region near the center of the explosion. (Note: This is just an example and not physically motivated - you probably want to resolve the shock front. See below for such applications with AMR.)

First of all, because SMR involves a complicated grid structure, HDF5 output is strongly recommended. To configure the code with HDF5,

> python configure.py --prob blast --flux hlld -b -hdf5

This configuration includes magnetic fields. To enable MPI parallelization, simply add the -mpi option. Note that no additional option is required for SMR and AMR.

Let us place grids near the center where the explosion occurs. First, you have to set the <mesh> refinement flag to static in the input file. Then the computational domain must be decomposed into smaller MeshBlocks, because a MeshBlock is the unit of refinement. As smaller MeshBlocks are more flexible but computationally inefficient, the size of Meshblocks must be set carefully. The following example produces a root grid of 64³ cells and MeshBlocks with 16³ cells:

...
<mesh>
nx1        = 64         # Number of zones in X1-direction
x1min      = -1.0       # minimum value of X1
x1max      =  1.0       # maximum value of X1
ix1_bc     = periodic   # inner-X1 boundary flag
ox1_bc     = periodic   # outer-X1 boundary flag

nx2        = 64         # Number of zones in X2-direction
x2min      = -1.0       # minimum value of X2
x2max      =  1.0       # maximum value of X2
ix2_bc     = periodic   # inner-X2 boundary flag
ox2_bc     = periodic   # outer-X2 boundary flag

nx3        = 64         # Number of zones in X3-direction
x3min      = -1.0       # minimum value of X3
x3max      =  1.0       # maximum value of X3
ix3_bc     = periodic   # inner-X3 boundary flag
ox3_bc     = periodic   # outer-X3 boundary flag

num_threads = 1         # Number of OpenMP threads per process
refinement = static

<meshblock>
nx1    =  16
nx2    =  16
nx3    =  16
...

Then we have to specify regions to be refined. In order to get double resolution near the center, add the following block:

<refinement1>
x1min=-0.5
x1max= 0.5
x2min=-0.5
x2max= 0.5
x3min=-0.5
x3max= 0.5
level=1

You can add more blocks like refinement2, refinement3, ... It is OK for the refinement regions to overlap each other - the code generates the minimally refined grid satisfying all the conditions. If you want to get even higher resolution, simply set a higher level.

Before actually running the simulation, it is highly recommended to check the grid structure by running the code with the -m option.

> ~/athena/bin/athena -i athinput.example -m 1
Root grid = 4 x 4 x 4 MeshBlocks
Total number of MeshBlocks = 120
Number of physical refinement levels = 1
Number of logical  refinement levels = 3
  Physical level = 0 (logical level = 2): 56 MeshBlocks, cost = 56
  Physical level = 1 (logical level = 3): 64 MeshBlocks, cost = 64
Number of parallel ranks = 1
  Rank = 0: 120 MeshBlocks, cost = 120
Load Balancing:
  Minimum cost = 1, Maximum cost = 1, Average cost = 1

See the 'mesh_structure.dat' file for a complete list of MeshBlocks.
Use 'python ../vis/python/plot_mesh.py' or gnuplot to visualize the mesh structure.

This outputs statistics about MeshBlocks created, as well as load balance if you specify the number of processes you want to use (replace the 1 in -m 1). This option also generates "mesh_structure.dat", which can be read using gnuplot.

> gnuplot
> splot "mesh_structure.dat" w l

Please note that even though there are 4³ MeshBlocks in the root level and 4³ MeshBlocks in the refined level, the total number of MeshBlocks is actually 120. This is because Athena++ does not solve overlapping MeshBlocks, and replaces one coarse MeshBlock with 8 (in 3D, 4 in 2D) fine MeshBlocks. This often results in a strange number (not power of 2) of MeshBlocks but you do not have to worry about it - it is not wasting computational resources, but actually is saving. You just need to assign a reasonable number of processes for parallel calculations.

If the grid structure looks as expected, you can run it and visualize the data as usual.

You can visualize the mesh structure using the "Mesh" plot on VisIt. However this plot can be too crowded because it draws all the cells. The "IndexSelect" operator is useful for extracting the borders between MeshBlocks.

You can extract borders (surfaces) of MeshBlocks by setting Incr I = nx1, J = nx2, and K = nx3.

Adaptive Mesh Refinement

The last topic covered in this Tutorial is Adaptive Mesh Refinement. With this feature, the code can create and destroy finer MeshBlocks adaptively based on the user-defined refinement condition function. For details, see Adaptive Mesh Refinement. Here we demonstrate how to implement the refinement condition and run the AMR calculation.

The refinement condition can be anything. For example, the density or pressure gradient (or curvature) is often used to detect discontinuities like shocks. For star formation calculations, the Jeans condition is a standard technique in which finer grids are generated so that the local Jeans length is resolved with a certain number of cells. In order to resolve the shock front in the blast wave test, let us add the following code in src/pgen/blast.cpp:

int RefinementCondition(MeshBlock *pmb);

void Mesh::InitUserMeshData(ParameterInput *pin)
{
  if(adaptive==true)
    EnrollUserRefinementCondition(RefinementCondition);

  return;
}

// refinement condition: density and pressure curvature
int RefinementCondition(MeshBlock *pmb)
{
  AthenaArray<Real> &w = pmb->phydro->w;
  Real maxeps=0.0;
  for(int k=pmb->ks; k<=pmb->ke; k++) {
    for(int j=pmb->js; j<=pmb->je; j++) {
      for(int i=pmb->is; i<=pmb->ie; i++) {
        Real epsr= (std::abs(w(IDN,k,j,i+1)-2.0*w(IDN,k,j,i)+w(IDN,k,j,i-1))
                   +std::abs(w(IDN,k,j+1,i)-2.0*w(IDN,k,j,i)+w(IDN,k,j-1,i))
                   +std::abs(w(IDN,k+1,j,i)-2.0*w(IDN,k,j,i)+w(IDN,k-1,j,i)))/w(IDN,k,j,i);
        Real epsp= (std::abs(w(IPR,k,j,i+1)-2.0*w(IPR,k,j,i)+w(IPR,k,j,i-1))
                   +std::abs(w(IPR,k,j+1,i)-2.0*w(IPR,k,j,i)+w(IPR,k,j-1,i))
                   +std::abs(w(IPR,k+1,j,i)-2.0*w(IPR,k,j,i)+w(IPR,k-1,j,i)))/w(IPR,k,j,i);
        Real eps = std::max(epsr, epsp);
        maxeps = std::max(maxeps, eps);
      }
    }
  }
  // refine : curvature > 0.01
  if(maxeps > 0.01) return 1;
  // derefinement: curvature < 0.005
  if(maxeps < 0.005) return -1;
  // otherwise, stay
  return 0;
}

This is a three-dimensional version of the criterion used in the double Mach reflection test (src/pgen/dmr.cpp). When a MeshBlock needs to be refined, this function returns 1. If this MeshBlock can be derefined, it returns -1. Otherwise it returns 0. This RefinementCondition is set using the EnrollUserRefinementCondition function. Then rebuild the code, and the input file should be modified to activate AMR as follows:

...
<mesh>
nx1        = 64         # Number of zones in X1-direction
x1min      = -1.0       # minimum value of X1
x1max      =  1.0       # maximum value of X1
ix1_bc     = periodic   # inner-X1 boundary flag
ox1_bc     = periodic   # outer-X1 boundary flag

nx2        = 64         # Number of zones in X2-direction
x2min      = -1.0       # minimum value of X2
x2max      =  1.0       # maximum value of X2
ix2_bc     = periodic   # inner-X2 boundary flag
ox2_bc     = periodic   # outer-X2 boundary flag

nx3        = 64         # Number of zones in X3-direction
x3min      = -1.0       # minimum value of X3
x3max      =  1.0       # maximum value of X3
ix3_bc     = periodic   # inner-X3 boundary flag
ox3_bc     = periodic   # outer-X3 boundary flag

num_threads = 1         # Number of OpenMP threads per process
refinement  = adaptive  # Adaptive Mesh Refinement
numlevel    = 3         # Number of AMR levels
derefine_count = 5      # Number of Steps

<meshblock>
nx1 = 8
nx2 = 8
nx3 = 8

# remove the <refinement1> block as it is not needed any more
...

The refinement flag is set to adaptive, and numlevel is 3, which means that the finer levels can be generated up to three levels including the root level. In other words, the maximum resolution is 2² = 4 times finer than the root level. The derefine_count parameter means that a newly created MeshBlock is not derefined for at least derefine_count time steps (5 steps in this case). This suppresses unnecessary creation and destruction of MeshBlocks.

In order to increase the flexibility, here we use MeshBlocks with 8³ cells. Again, the size of MeshBlocks must be set considering balance between flexibility and performance, and generally speaking it is recommended to use 16³ cells or larger for performance. Note that this is a pretty time consuming calculation, so it is better to use MPI parallelization. Alternatively, you can do this in 2D, by setting nx3 to 1 in <mesh> and <meshblock>.

At the beginning of the simulation, the code checks the initial condition and generates finer levels iteratively to make sure the initial grid structure satisfies the refinement condition. Often the code will show the following warning if the initialization creates a lot of MeshBlocks. It is OK to ignore this message if it is expected.

### Warning in Mesh::Initialize
The number of MeshBlocks increased more than twice during initialization.
More computing power than you expected may be required.

While the simulation is running, AMR will create and destroy finer MeshBlocks according to the refinement condition provided.

The visualization can be done in the same way as the SMR case using VisIt. However, perhaps because of the XDMF reader of VisIt, there will be some trouble when you seek with the time bar. When the number of MeshBlocks increases between the time steps, VisIt simply visualizes as many MeshBlocks as in the previous time step. In this case, simply click the "ReOpen" button to reload the XDMF file. When the number of MeshBlocks decreases, VisIt will crash but do not panic. Again, click the "ReOpen" button, and then the "Draw" button to recover the failed plots. We are looking for better ways of visualization.

---

Continue tutorial with Special and General Relativity