Refinement

Refinement covers subdomains with grids of higher resolution. Refinement is probably the most algorithmically complex operation in Berger's AMR strategy. Like integration, it is implemented recursively: a refinement at level l first refines level l+1, and so on recursively to the finest level. The reason is that for some problems, there are features which appear only in the finer level and not in the coarser level. If we start refining the grid from the coarser level as in the local uniform grid refinement (LUGR) method (Verwer and Blom [17]), those features cannot be captured. A pseudo-code description of the refinement algorithm is presented in Figure  3.2 (see also [12]).

  

Figure: Pseudo-code description of the basic refinement algorithm

We will explain the algorithm in more detail.

Except for the initial mesh refinement, the refining should start from the finest level grid. For the initial mesh refinement, since we have only the base grid, the refining starts from the base grid; the solutions on the finer grid are obtained from the initial conditions instead of from the coarser grid.

The selection algorithm flags the bad points and produces the appropriate regions of interest to be refined. A variety of selection criteria can be applied according to different monitor functions. We can construct the monitor functions from three sources: operator related monitor m_o, solution related monitor m_s and regularization related monitor m_r. These three kinds of monitor functions have been discussed in detail by Hyman and Naughton[6]. In the AMR method, because we use the monitor only to flag the bad points, not to move or distribute the grid as in the moving mesh method and conventional static rezoning, we do not need to worry about mesh crossing and mesh tangling. The monitor function can be very relaxed. Here we adopt a very simple monitor proposed by Verwer et al. [17] for local uniform grid refinement (LUGR). For each grid, we compute the spatial monitor

$$\mbox{SPCMON}(i) := \displaystyle\max_{ic}\mbox{SPCTOL}(ic)(\vert\Delta x^2u_{xx}^{ic}(i)\vert)$$

where

$$\mbox{SPCTOL}(ic) := \d{\mbox{SPCWGT(ic)}}{\mbox{UMAX}(ic)\cdot \mbox{TOLS}}.$$

The parameters in the above equations are user-specified quantities. $0< \mbox{SPCWGT} \leq 1$ is a weighting factor for the relative importance of a PDE component in the monitor. $\mbox{UMAX}$ denotes the approximate maximum absolute value for each component, and $\mbox{TOLS}$ is the spatial error tolerance. The second order derivatives of the solution are approximated by second order finite differences. A next level of refinement is required if some $\mbox{SPCMON}(i) > 1.0$, and then all grid points with $\mbox{SPCMON}(i) > 0.5$ are flagged. In order to ensure proper nesting, if the current level grid has a grandparent, those points should be flagged. Berger et al.([1], [2], etc.) adopted a seemingly more rigorous, and ultimately more expensive approach: Richardson truncation error estimation. Difficulties with this approach have been discussed by Neeman [12]. One difficulty is that this approach requires access to the time integration routine, which may be supplied by the user and may need recoding for this purpose. Expansion serves to ensure coverage of the grids on the current level for the next k time steps. To reduce the overhead of the refinement, we do not perform the refinement every time step. Instead we perform the refinement every k time steps. Therefore we must add enough buffer zones to the flagged points so that a discontinuity will not escape onto the coarse grid. The choice of k is dependent on the number of buffer zones and on the current CFL number. There is a trade-off between the number of buffer zones and the size of k. In our 1D AMR system, we tried these two approaches. In the first approach, we performed refinement only after integration on the whole grid was completed, i.e., perform Refine(1). The number of buffer zones is different for different levels in this approach, because finer levels have more integration steps so that they need more buffer zones. Usually, we assume the discontinuity moves no more than one buffer zone within one time integration. So we choose the number of buffer zones to be equivalent to one base grid. (Actually, in our tests the number of buffer zones can be much less than that). In the second approach, we perform refinement for the current level after it has advanced a fixed number of time steps (usually this is equal to factor). The number of buffer zones can be fixed and the same for all the grid levels. The finer levels have more refining. One of the advantages of this approach is that the number of buffer zones can be much smaller than in the first approach, so that the patches for finer levels can be much smaller. This saves computation time when integrating the finer level. However, the overhead due to the refinement is much larger than for the first approach. Another reason why we consider the first approach is that if the refinement is input by the user and not adapted using our monitor function, it is much easier to implement via the first approach, i.e., the refinement will not change during a number of time steps of the base grid. Following expansion, the flagged points at level l are clustered. The clustering algorithm is the most sensitive part of Berger's 2D AMR system. However, for 1D problems, it is very simple. We can just search from left to right for each coarse patch and find the adjacent flagged points and group them together. If two patches are too close (within two buffer zones), we join them as one.

Next the l+1 level is regridded. The regridding includes computing the physical locations for each fine grid and copying or injecting the solution from the old grid to the new grid. The physical positions are easy to compute by linear interpolation. The solution needs more attention. Although the solution can be obtained from the coarse grid by injection or interpolation, a more accurate solution is obtained from the old grid at the same level which partially overlays with the new grid. One of the secrets to the success of the AMR algorithm is that flow discontinuities always fall within these regions of overlay between the new grid and old grid. Thus the adaptation process cannot introduce further errors in these problem regions. Consequently, near flow discontinuities the numerical solution largely exhibits the properties of the numerical scheme used for the single mesh integration process. The method used to interpolate the field solution from the coarse grid needs to be chosen with care. In regions near a discontinuity, linear interpolation does not produce a good result. Conservative interpolation should be used. In our implementation, we use the Van Leer limiter to compute the interpolated values. The coarse grid solution is assumed to be piecewise linear. The slopes for each grid are found by applying a MinMod limiter function (see V. Leer[7]) to the forward and backward slopes between cell centers. So, for a coarse cell (i),

$$\phi_{i+\frac12} - \phi_{i-\frac12} = MinMod(\phi_{i+1} - \phi_i, \phi_i-\phi_{i-1}),$$

where

$$MinMod(a,b) = \cases{0, \quad \quad \quad \mbox{if}\ ab<0 \cr... ...\mbox{min}(\vert a\vert,\vert b\vert), \ \mbox{otherwise.}\cr}$$