Adaptive mesh refinement (AMR) is a technique that has been featured only sporadically in atmospheric science literature. This paper aims to demonstrate the utility of AMR for simulating atmospheric flows. Several test cases are implemented in a 2D shallow-water model on the sphere using the Chombo-AMR dynamical core. This high-order finite-volume model implements adaptive refinement in both space and time on a cubed-sphere grid using a mapped-multiblock mesh technique. The tests consist of the passive advection of a tracer around moving vortices, a steady-state geostrophic flow, an unsteady solid-body rotation, a gravity wave impinging on a mountain, and the interaction of binary vortices. Both static and dynamic refinements are analyzed to determine the strengths and weaknesses of AMR in both complex flows with small-scale features and large-scale smooth flows. The different test cases required different AMR criteria, such as vorticity or height-gradient based thresholds, in order to achieve the best accuracy for cost. The simulations show that the model can accurately resolve key local features without requiring global high-resolution grids. The adaptive grids are able to track features of interest reliably without inducing noise or visible distortions at the coarse–fine interfaces. Furthermore, the AMR grids keep any degradations of the large-scale smooth flows to a minimum.
Global climate models have become vital tools for simulating the present and future climate and for predicting important climate trends. However, current global models are limited in their ability to represent many multiscale aspects of atmospheric flows. Their resolutions, limited by computational costs, are too coarse to accurately represent key processes that span a wide range of temporal and spatial scales. High-resolution simulations are essential for capturing these scale interactions and for accurately representing local and regional phenomena such as convection, orographically induced precipitation, mesoscale storm systems, and tropical cyclones. Such events have large regional impacts as well as broader feedbacks onto the large-scale climate system. Today, high-resolution general circulation models (GCMs) used for global climate simulations can utilize uniform grid spacings down to 10 km as documented by Manganello et al. (2012) or Kinter et al. (2013). However, they are computationally very expensive and still unable to represent key processes such as clouds explicitly. Exceptions are the cloud-permitting, partly cloud-resolving, global simulations by Miura et al. (2007), Putman and Suarez (2011) or Miyamoto et al. (2013) that were run for short, multiday time periods with grid spacings in the 0.9–3.5-km range. To employ such high resolutions over longer climate time periods, modelers typically use limited-area models (LAMs) that focus the computational resources in areas of interest. A major drawback is that LAMs require their lateral boundaries to be externally forced. These boundary conditions are typically derived from much coarser GCMs, which use different numerical schemes and physical parameterizations, thereby introducing possible biases or numerical discrepancies. In addition, it is an open question how well LAMs can capture teleconnections of global large-scale dynamics and localized features particularly for tropical cyclones and other phenomena that have feedbacks onto the larger climate system.
Variable-resolution GCMs can utilize static or dynamic grid refinements, which are promising options to bridge the gap between global and regional climate modeling. Application examples for static (nonmoving) mesh adaptations are provided in Zarzycki et al. (2014), Rauscher and Ringler (2014), Zarzycki and Jablonowski (2015), and Huang et al. (2016) (see also further references therein). Our paper focuses on dynamically adaptive grids, which track features of interest during the model simulation by locally adding or removing grid points as needed. Adaptation criteria based on error estimates (e.g., Skamarock et al. 1989; Behrens 1998; Blaise and St-Cyr 2012) or flow characteristics (e.g., Hubbard and Nikiforakis 2003; Jablonowski et al. 2006, 2009; St-Cyr et al. 2008) can be used to determine where the high-resolution mesh should be placed. By increasing resolution only locally, dynamic refinement significantly decreases the total number of degrees of freedom for the simulation. However, since dynamic refinement is used within a global model, it also eliminates the need for forced boundary conditions and solves the local high-resolution area and global flow using the same dynamical core and physics package. Global models allow for a better representation of global and synoptic-scale phenomena and permit them to interact better with meso- and small-scale features that can be resolved in the model. Additionally, a key advantage of dynamic refinement compared to a static refinement setup is that the location of the refined area does not have to be determined a priori.
Adaptive refinement methods are well established in many areas of computational fluid dynamics like aerospace engineering or space weather modeling (e.g., Tóth et al. 2012). In atmospheric science though, the costs and benefits of AMR methods have mostly been evaluated in idealized simulations and simplified models so far. Some of the first adaptive atmospheric models were developed by Skamarock et al. (1989), Skamarock and Klemp (1993), and Dietachmayer and Droegemeier (1992). In general, the grid refinement strategies can be categorized into three overarching types: r refinement (or r adaptivity), h refinement, and p refinement.
In r refinement, the number of grid cells remains unchanged; instead, the cells are dynamically redistributed to increase resolution in parts of the grid while decreasing it everywhere else. This dynamic grid adaptation creates smoother transition regions between resolutions but requires a complex global remapping of the mesh to move the location of the high resolution. Dietachmayer and Droegemeier (1992) use this global grid redistribution technique to increase resolution in areas where the estimated solution error is high. Giraldo (2000) and Iselin et al. (2002) have also applied this type of dynamic adaptation for the 2D shallow-water equations and advection problems, respectively. More recently, Kühnlein et al. (2012) implemented r adaptivity within a 3D Cartesian framework, Bauer et al. (2014) used r-refinement grids guided by error estimates in a shallow-water model, and Weller et al. (2016) demonstrated r-refinement use on the sphere.
Adaptive mesh refinement (AMR), another term for h refinement, increases resolution locally either by adding cells within the grid structure or by overlaying additional cells of finer resolution on top of the grid without changing the base grid structure. Skamarock et al. (1989) and Skamarock and Klemp (1993) implemented AMR by placing finer-resolution meshes over the coarse grid in areas which had large truncation error estimates. The gridcell solutions and boundary conditions between the higher-resolution meshes and the base grid are continually updated. In a more recent example, Chen et al. (2011a) use AMR that overlays high-resolution meshes in areas of interest in a shallow-water model on a cubed-sphere grid. Examples of AMR techniques that locally add and remove cells to the base grid for the shallow-water equations on the sphere have been presented by Behrens et al. (2005), Läuter et al. (2007), St-Cyr et al. (2008), and Marras et al. (2015). Both Behrens et al. (2005) and Läuter et al. (2007) use conformal unstructured finite-element meshes. In conformal grids, each cell shares an edge with exactly one other cell, while on a nonconforming grid, cells can share an edge with more than one neighboring element. The two AMR models described in St-Cyr et al. (2008)—a block-structured finite-volume method on a latitude–longitude grid and a spectral-element method on a cubed-sphere grid—use nonconforming meshes and a quad-tree based AMR method with gradient- or vorticity-based refinement criteria. Marras et al. (2015) compared the use of an AMR approach on several structured and unstructured nonconformal grids. Use of AMR in a regional model paired with a physical parameterization package was presented in Bacon et al. (2000). Furthermore, AMR methods have also been investigated for 2D flow fields in Cartesian geometry. Recent examples include Müller et al. (2013) and Kopera and Giraldo (2014) who analyzed a tree-structured AMR algorithm for nonhydrostatic dynamical cores in the x–z plane, and Hendricks et al. (2016) who explored static and dynamic AMR for tropical cyclone–like vortices in a shallow-water model on an f plane with a constant Coriolis parameter f.
The third type of dynamic refinement, p refinement, holds the grid spacing fixed but changes the order of the polynomial approximation within each grid element to increase local resolution. Use of p refinement with discontinuous Galerkin (DG) methods for the shallow-water equations were described by Kubatko et al. (2009) and Tumolo and Bonaventura (2015). A hybrid refinement method that combines the h- and p-refinements methods was demonstrated by Eskilsson (2011), and Blaise and St-Cyr (2012) used an hp-adaptive DG method to model the shallow-water equations on the sphere for global tsunami simulations. Recently, Aechtner et al. (2015) implemented a new adaptive wavelet approach for local dynamic refinement with the 2D shallow-water equations on the sphere.
The purpose of this study is to demonstrate the pros and cons of using AMR, h refinement, for simulating atmospheric flows. It assesses the effectiveness of AMR, employing a nonconformal grid, in achieving similar results as uniform-grid simulations while reducing computational cost. Furthermore, it is revealed that AMR can be implemented without harming or degrading the large-scale smooth flows or inducing numerical noise and wavelike reflections at AMR boundaries. The 2D shallow-water equations serve as a useful test bed for an AMR model as they exhibit many of the complexities present in a full 3D model. We utilize the cubed-sphere fourth-order finite-volume AMR model presented in McCorquodale et al. (2015) for the 2D shallow-water equations on the sphere. The model implements a mapped-multiblock AMR technique that overlays refined patches on the coarser grid. Our work tests the model’s ability to track and refine over dynamic small-scale features of interest and to evaluate refinement criteria. We investigate various refinement criteria, such as thresholds for the height gradient or relative vorticity, which guide the locations of refinement patches. In addition, we shed light on factors that may limit AMR applications including the size of the refinement ratios between grid levels and the total number of levels. Last, we examine the effect of AMR on large relatively smooth flows where extra refinement is unnecessary. Specifically we focus on how the coarse–fine interfaces of the AMR patches influence the overall flow and error.
The paper is organized into three main sections. A brief description of the model and a discussion of the multiblock AMR techniques are provided in section 2. Section 3 discusses the results of the numerical tests. One advection test and four shallow-water tests are presented. The tests are the moving vortices advection test by Nair and Jablonowski (2008), a steady-state geostrophic flow [test case 2 from Williamson et al. (1992)], the unsteady solid-body rotation test of Läuter et al. (2005), a shallow-water test consisting of a gravity wave impinging on an isolated mountain, and last a test that assesses the interaction of idealized binary vortices. Conclusions are provided in section 4.
2. Model description
The Chombo-AMR dynamical core (dycore) is a new model that is built upon an unstaggered high-order finite-volume (FV) multiblock approach with a classical fourth-order Runge-Kutta (RK4) time discretization scheme. A detailed description of the model setup to solve the shallow-water equations in conservative flux form can be found in McCorquodale et al. (2015). In addition, the Chombo-AMR library is described in Adams et al. (2015). The finite-volume approach is implemented on an equiangular cubed-sphere grid. This grid consists of six separate panels that are projected onto the surface of the sphere. The mesh thereby eliminates the singularities due to converging meridians at the poles found in spherical latitude–longitude grids. Additionally, the equiangular cubed sphere leads to a quasi-uniform spherical grid with grid cells of similar size across the sphere. The discrete resolution of the cubed-sphere grid is represented as , where denotes the number of grid cells in each direction on the six panels. A list of properties of the equiangular cubed-sphere grid, including the approximate grid spacings, is given in Table 1 for several resolutions. The finite-volume method for the spatial discretization uses a fourth-order accurate discretization to compute flux averages on the faces. The central difference operators used to obtain the fluxes are smoothed by an explicitly added sixth-order diffusive operator that maintains the fourth-order accuracy of the scheme [see McCorquodale et al. (2015) for details]. No additional limiters or filters are implemented. The numerical scheme is mass conserving to machine precision and energy conserving up to the temporal truncation order, when used without limiters or explicit dissipation. The total-energy conservation properties for the model with added dissipation are demonstrated in Fig. 11 of McCorquodale et al. (2015).
Since high-order FV schemes make use of neighboring elements, a mapped-multiblock approach is used to coordinate the remapping of element values that are needed for the flux calculations across panel boundaries on the cubed sphere. Though the cells at panel edges are conformal with neighboring cells across panel boundaries, the transition between the panels is not smooth due to the separate mapping on each panel. To preserve the order of accuracy of the fluxes, the domain is expanded at the panel edges with the addition of three layers of ghost cells to perform the FV calculations on each panel. As a result of different mappings, ghost cells of one panel will not have the same shape as cells on the neighboring panel. Therefore, the values in the ghost cells are set by least squares interpolation from a stencil of surrounding cells that are within the domain of the ghost cell’s panel and on neighboring panels [section 3.4 in McCorquodale et al. (2015) describes in detail the interpolation process]. Additionally, flux values for the cell faces that lie on a panel edge are calculated separately for each panel, and the mean of the two fluxes is taken as the value for that face to ensure conservation. Thus, communication between the separate domains for each panel is limited to the fluxes at the domain boundary and the neighboring cell values needed to interpolate the solution to the ghost cell regions. The block-structured AMR method allows for further subdivision of the computational domain of each panel into rectangular regions of grid cells called patches, which allow the calculation to be distributed efficiently on parallel computing platforms.
Adaptive mesh refinement
Our 2D AMR shallow-water model uses the strategies within the Chombo library (Adams et al. 2015) that has been designed for parallel computing architectures. AMR calculations are performed on a hierarchy of nested meshes, called levels, which have a defined refinement ratio between them. This refinement ratio must be a power of 2. The finer levels are overlaid on top of the coarser levels and are organized in the block structure described in the previous section. Figure 1a provides an example of the AMR grid structure with two refinement levels. Whenever new cells are created, they are initialized via interpolations from the coarser level and ghost cells are used to calculated the fluxes at patch boundaries. At these coarse–fine interfaces, as at panel boundaries, the space–time accuracy drops from fourth order to third order due to a lack of error cancellation that would normally occur with the FV method. Intermediate levels must have a sufficient number of cells separating a finer level and a coarser one. This ensures that the finer level is properly nested so that the interpolation to fill ghost cells on the finer level can be performed using cells from only one level.
Figure 1b depicts an overview of the time-stepping and subcycling process. Rather than having a uniform time step dependent on the smallest gridcell size, Chombo-AMR subcycles the refined levels in time maintaining a constant Courant number. As noted by Ullrich (2014) the single-wave-mode characteristics of a numerical method often have an unexpected and nonlinear dependence on the Courant number. Therefore, a constant Courant number helps ensure consistent dispersive properties across the refinement levels. The typical work flow for advancing an AMR grid level l in time is as follows:
Regrid levels finer than l if required:
Evaluate the refinement criterion and mark (tag) all cells that should be included in finer levels. In these regions, new blocks of cells at levels are overlaid. The new cell values are interpolated from the coarser level using a fourth-order least squares algorithm that maintains conservation as described in section 4.1 of McCorquodale et al. (2015).
Advance level l one time step using the RK4 time-stepping method.
Interpolate values to the ghost cells that surround level using the least squares algorithm also implemented for ghost cells at panel boundaries. Three layers of ghost cells are required. The interpolation does not need to be conservative as the ghost cell values are only used to reconstruct the flux on the faces of the level cells. Figure 2 depicts the location of ghost cells for two of the three layers and the stencil of coarse grid cells used for their interpolation.
Perform previous steps for level . Level is advanced using refined time steps (subcycling) as depicted in Fig. 1b. A temporal interpolation closely related to the RK4 method is used to update the values in the level ghost cells from cells on level l at the intermediate time steps (McCorquodale et al. 2015).
After the subcycling in time is completed, average the solution from level and sum up the fluxes to update the values on the coarse grid. Corrections are applied to the fluxes at coarse–fine interfaces to ensure conservation (Berger and Colella 1989).
As mentioned in the description of the first step, the additional levels are placed in locations that have been tagged by the refinement criterion of the model. The size of the refined grid that is added over tagged cells is determined by three aspects: 1) the need to ensure proper nesting of finer levels, 2) parameters from the Chombo library designed for efficient parallelization, and 3) a user-defined buffer parameter. Refinement (tagging) criteria are based on thresholds for user-selected flow properties, like relative vorticity, that indicate where refinement should be placed. The tagging strategies can be based on a variety of properties including tracer values, vorticity thresholds, gradients, or a combination of several criteria. The type of refinement criteria and the threshold values are set independently for each simulation. The threshold values can be uniform across all refinement levels or designed to scale with increasing refinement. For example, the relative vorticity threshold can be set so that it increases with increasing resolution. With the idealized test cases presented in this paper, the selection of refinement criteria and thresholds was problem dependent, with the simplicity of the tests offering only a few possible options. A range of refinement thresholds for the AMR test cases was explored. Here, we present various thresholds for the gravity wave and binary vortices tests to demonstrate how changes in refinement affect the solutions and grids. For the moving-vortices advection test we only present a single threshold value that represents a compromise between reducing the error and increasing the computational run time.
The Chombo-AMR dycore is designed to have multiple refinement levels (up to 10), and the maximum number of levels is set for each simulation. In our simulations, we explore the effects of adding up to two refinement levels (called two-level AMR). In addition, we explore the grid-resolution refinement ratios between successive levels, and present selected results for three powers of 2: ×2, ×4, and ×8. As an example, a grid with a base level of c32 (2.8°) resolution and two levels of ×4 refinement has a maximum resolution of c512 (0.18°).
3. Results of the numerical experiments
For our assessment of the Chombo-AMR shallow-water model, we select five test cases: one advection test and four shallow-water tests. The test cases are divided into two categories: large-scale smooth flows and simulations with either sharp gradients or strong, nonlinear, localized flows. The first category comprises the following two shallow-water tests:
These large-scale flows, which have no realistic use for AMR, serve as “do no harm” tests. They are used to check the model’s ability to preserve the characteristics of smooth flows as they cross the AMR patches. We measure the impact that the refinement ratios, the number of AMR levels, and the location of the refinement patches have on the solution. Convergence tests are also performed with these test cases that both have analytical solutions.
The second test category with localized flow features consists of three tests for which AMR could improve the solution, and we seek to evaluate how effectively it is able to do so. These tests are
the moving vortices advection test by Nair and Jablonowski (2008) (with analytical solution),
a gravity wave impinging on an idealized mountain shallow-water test, and
a binary vortices test case in which two vortices interact.
The bottom two tests do not have analytical solutions and the evaluations rely on high-resolution reference solutions. A variety of refinement criteria are used with these test cases to demonstrate the AMR’s ability to track, adapt to, and resolve these localized features accurately. The model results are presented using normalized and error norms as defined in Williamson et al. (1992). Additionally, the total number of grid cells quoted for AMR runs include the sum of all valid grid cells from all refinement levels (not including ghost cells) since the finer levels overlay the coarser grids beneath them. The total number of grid cells can serve as a rough benchmark of computational cost when comparing AMR runs to uniform runs.
a. Moving-vortices advection test
The moving-vortices test is a challenging deformational-flow advection test proposed by Nair and Jablonowski (2008). The test represents the roll-up of an initially smooth tracer into tight spiral bands. The roll-up creates steep gradients that provide a good test for the AMR. In this test, a pair of vortices is generated on diametrically opposite sides of the sphere. The wind field is the summation of a solid-body rotation and a deformational flow such that the two vortices move along a great circle and an exact solution is known at all times [see Nair and Jablonowski (2008) for details]. A 12-day time period is simulated that advects the spiraling vortices once around the sphere. The background flow is prescribed with a rotation angle of so that the two vortices are advected through the corners of the cubed sphere (located at ). Figures 3a–d depict the analytical solution for the roll-up of the tracer at days 0, 4, 8, and 12.
Numerical tests were carried out with uniform grids and AMR grids with two different tagging criteria: a tracer-gradient tagging and a combination tagging. In tracer-gradient tagging, the model refines in regions where , where q is the unitless value of the passive tracer. Combination tagging combines the tracer-gradient tagging criterion with a relative-vorticity tagging that refines in areas where , so that the model refines over regions in which either threshold is reached. The gradient threshold value demonstrates a balance between reducing error and limiting computational costs. A lower threshold increases run time without significantly reducing error and a higher one results in a large increase in error. The vorticity threshold is set to maximize initial refinement around the vorticity patches without being too low as to refine on the background vorticity. With both criteria, a ×4 refinement ratio between levels was used. The evolution of the grids can be seen in Fig. 3, which shows the tracer error for the two AMR runs. Figures 3e–h depict the two-level AMR run (c32/c128/c512), which has a base grid with c32 resolution, a c128 intermediate AMR level, and a finest c512 AMR level, using the tracer-gradient tagging. Figures 3i–l depict the same two-level AMR setup but with the combination tagging. The adaptive block structure is clearly successful at tracking the evolving vortices. The block structure of patches, not the individual grids cells, is outlined by the black lines in Fig. 3. Each patch contains a user-selected maximum number of grid cells. However, patches might also be split into smaller blocks by the Chombo library. This optimizes the load balancing of the model on parallel computing architectures and increases the efficiency of the code. Therefore, various patch sizes are present in Fig. 3.
For the tracer-gradient tagging scheme (Figs. 3e–h), the error is largest near the center of the spiral. The combination tagging scheme (Figs. 3i–l) significantly reduces the error near the center, and the largest error is now toward the outer edges of the spiral just beyond the coarse–fine grid boundaries. When using only the tracer-gradient tagging, refinement does not begin until after day 1, so errors accumulate on the coarse grid until that time and remain higher even after being overlaid with finer levels. With the combination tagging, the central region is already refined (see Fig. 3i), which limits the error growth. This effect can also be seen in Fig. 4a where a time series of the normalized tracer errors with respect to the analytic solution are depicted for 15 different model configurations. In particular, Fig. 4a compares the tracer errors of the five uniform-resolution simulations c32, c64, c128, c256, and c512 (see Table 1 for the associated grid spacings) to various one- and two-level AMR experiments with either the tracer-gradient or combination tagging. The time evolution of the associated total number of grid cells, including the underlying coarse cells, is shown in Fig. 4b. In general, the error in the combination tagging simulations reaches the lower error values of the uniform runs (with the same resolution as the finest AMR level) more quickly than the tracer-gradient tagging ones, despite having a significant difference in gridcell count for only the first few days (see Fig. 4b). The errors in all single-level AMR runs in Fig. 4a using both tagging criteria converge to the error of uniform runs at the highest resolution (e.g., the c32/c128 AMR run and the uniform c128 run). This result is achieved using significantly fewer grid cells than the uniform runs. AMR limits the normalized error growth. AMR runs begin with global errors that are near the level of uniform runs with a grid resolution of their coarsest grid. However, their global error increases at a much slower rate than that of the uniform runs until the AMR error approaches the error level of the matching high-resolution uniform run. In the simulations presented here the two-level AMR runs see a diminishing effectiveness at decreasing the global error as the errors from the coarse section of the grid dominate.
The error at day 12 as a function of total number of grid cells for both uniform and AMR runs is plotted in Fig. 5. AMR runs that have error values plotted to the left of the errors of the uniform runs (line of hollow black squares) achieve a lower error than a uniform run with a comparable number of grid cells. While the one-level AMR runs show a decrease in error while using fewer grid cells, the two-level AMR runs generally result in a higher error for a comparable number of grid cells than the uniform runs. Only the two-level AMR run with a c16 base using the combination tagging has a slightly improved error compared to the uniform runs (see the leftmost magenta star). These results demonstrate the reduced improvement to the global error from additional levels of AMR with very coarse base meshes. However, AMR still slows down the error growth over time in comparison to uniform runs, and even asymptotes to the finest uniform mesh errors (e.g., c128/c512). Extending the run time to 18 days, the error from the c32 base two-level AMR run with the combination tagging (orange triangle in Fig. 5) is now slightly lower than the error from the uniform runs at day 18 (red downward-pointing triangles). Figure 5 also confirms that the convergence rate for the uniform run is fourth order in the normalized error.
Nair and Jablonowski (2008) applied this test using an α = 0° setup (not α = 45° as we present here) in an FV-AMR model using a coarse 5° base grid and one to three levels of ×2 refinement, which were guided by a tracer-gradient threshold. Their AMR errors were generally the same or lower than the errors of their comparable uniform runs for coarse grids with one or two levels of refinement. However, for their three-level AMR run (finest resolution 0.625°) errors were slightly higher than their uniform 0.625° run, agreeing with what we observe for multiple levels of AMR. Additionally, our error measures for uniform-resolution runs are comparable to other higher-order models. Our uniform c64 (~1.4°) run with α = 45° has a normalized error of 8.9 × 10−3 after 12 days, which is comparable to the results from the 1.25° grid using a multimoment method in Chen et al. (2011b) and slightly higher than the c80 (1.125°) run in Lauritzen et al. (2010).
We also include the total run time versus number of grid cells in Table 2, for some of the 12-day runs in Fig. 5. In this case, the wall-clock time (as a % of the finest uniform run) is closely related to the total number of grid points, because most of the time steps and grid cells are at the finest level. For AMR runs, coarser levels must be completed before fine levels. There is, therefore, a slight performance penalty for the two-level AMR (c32/c128/c512) run, which actually increases the total number of finest-level points. This configuration only has 24 boxes (each cells) at the c32 level to be distributed across 32 processors on the Yellowstone computing platform [operated by the National Center for Atmospheric Research (NCAR)]. This means that some processors run idle when the solution on the coarsest grid is computed, which slightly lessens the parallel performance of this AMR run. Overall, a good heuristic for this test is that the run time is approximately proportional to the total number of grid cells at the finest level.
b. Global steady-state geostrophic flow
This zonal steady-state flow is the second test case from the Williamson et al. (1992) test case suite for the shallow-water equations on the sphere. The test consists of a solid-body rotation along an axis that differs from the polar axis by angle α and a corresponding balanced height field. We use the more challenging case, which means that the flow travels over the cubed-sphere corner points at a 45° angle. Since the flow is initialized in a gradient-wind balance, any changes from the initial conditions (which serve as the analytical solution) are considered errors. No topography is present. The test was designed to measure how well the model can maintain a large-scale smooth balanced flow. Thus, we expect little benefit from AMR refinement. We use the test primarily to assess the sensitivity of the flow to the grid structure and abrupt changes in the grid resolution along coarse–fine mesh interfaces.
We implement two static refinement configurations, which can be seen in the bottom two panels of Fig. 6. The first configuration (Fig. 6b) consists of a statically refined patch centered at 0° latitude and longitude (fully contained within an equatorial cubed-sphere panel) over an area with strong gradients in the height field. In the second configuration (Fig. 6c), we place static patches over the locations of high relative vorticity, refining where . This criterion results in two midlatitudinal patches that transect polar-equatorial panel boundaries, a challenging location for the cubed-sphere grid. Using these static refinement configurations, we ran simulations that have one and two levels of refinement using ×2, ×4, and ×8 refinement ratios. Increasing the refinement ratio permits us to test how abruptly resolution can increase without harming accuracy or causing spurious numerical noise at the boundary between the coarse and fine regions. We compare the height errors, characterized as the difference between the analytic initial condition and the numerical solution for the height at day 5, of uniform-resolution simulations and simulations using the equatorial and midlatitudinal patch configurations. The addition of refined patches should ideally result in no additional error if the flow is well resolved, or a small decrease in global height error if it is not.
Results in Table 3 show the normalized and height errors for the uniform c32 (2.8°) run and c32 base runs with the two refinement configurations after 5 days. The errors for the runs with the equatorial patch are essentially unchanged from that of the uniform c32 run. Even the extreme cases of a ×8 refinement ratio or multiple levels of refinement increase the height error by less than . The results for midlatitudinal patch runs show a reduction in global error with even a ×8 refinement ratio reducing the height error by more than and the error by at least compared to the c32 uniform run. The height error plots at day 5 for the uniform c32 run and the runs for the two refinement configuration using two levels of ×4 additional refinement (Figs. 6a–c) depict a similar result. The equatorial patch run (Fig. 6b) has essentially the same height error profile as the uniform run (Fig. 6a). The height error for a base c32 run with the midlatitudinal refinement patches (Fig. 6c) shows a clear improvement in error since the coarse base resolution does not fully resolve the flow and the refined patches cover areas of high error. The refined patches do not create any spurious wave reflections or lead to an increase in error along the coarse–fine boundary.
Figure 7 depicts a comparison of the day-5 normalized and height errors for runs with the midlatitudinal refinement patches and uniform runs for base resolutions of c16 (5.6°) to c256 (0.35°). At coarse resolutions, we see a slight improvement in the error for runs with refinement compared to the uniform runs. However, at higher resolutions when the flow is well resolved the change in error is indistinguishable. The figure also shows that fourth-order convergence is maintained in the runs with the static refinement patches. Results for runs with the equatorial refinement patch at other higher resolutions followed a similar pattern as the c32 base runs in Table 3 and also demonstrated fourth-order convergence (not shown).
The steady-state geostrophic flow test case has been used in other AMR and static refinement studies. Similar refined grid locations were used with the finite-volume AMR models in Chen et al. (2011a) (on the cubed sphere) and St-Cyr et al. (2008) (on a latitude–longitude grid). In both models, the introduction of refined patches led to increases in error when compared with the uniform runs, with significantly larger increases in the height error for configurations in which the refinement patch was placed over strong height gradients. The error increased by in Chen et al. (2011a) and a factor of 2.5 in St-Cyr et al. (2008). However, for the higher-order spectral-element method (SEM) in St-Cyr et al. (2008), the error was considerably reduced with the addition of a refined patch. Weller et al. (2009) tested a number of grid geometries with variable resolutions using a ×2 refinement ratio and found increases in error when a refinement patch was added. Harris and Lin (2013) used a nested-grid FV model with a ×3 refinement ratio. After 5 days their height errors roughly doubled in comparison to their uniform run, though their errors were nearly unchanged. Our results with static refinements are very competitive as they show almost no increase in error or even in some cases an improvement in the error. Thus, our model preserves the large-scale flow and limits the errors at the refinement patches very effectively.
c. Unsteady solid-body rotation
The time-dependent zonal flow test proposed in Läuter et al. (2005) (example 3) consists of an unsteady solid-body rotation (USBR), which is forced by topography. It possesses an analytical solution. The large-scale flow and topography are smooth, zonally symmetric, and somewhat artificial. In particular, the topography is zero at the equator and rises to its maximum (around 11 km) at both poles. As with the previous test case, we expect little benefit from AMR given the smooth characteristics. The benefit of the USBR test is that it has the added complication of moving features that can be tracked with AMR while still having an analytic solution to determine the error. One can observe how well the flow is maintained and whether numerical artifacts, if any, are created by the resolution change at grid boundaries and by the AMR regridding process. Using the setup described in Läuter et al. (2005), we have set the parameter to let the flow field pass over the corners of the cubed sphere at a 45° angle. To force the initial condition to repeat itself after exactly 1 day for better comparison of the results, Earth’s angular velocity is slightly modified to be based on a solar day instead of sidereal day, so that the angular velocity is s−1 s−1.
We conduct a series of simulations over a range of base resolutions with either one or two levels of refinement using three refinement configurations:
Statically refined patch used in section 3b now centered at 0°, 90°E.
Dynamic AMR refinement with height tag, the threshold for the free surface height is m.
Dynamic AMR refinement with relative vorticity tag, the threshold is .
The three grid configurations can be seen in Figs. 8b–d, which show the patches at the identical initial and final (day 5) positions. The second and third configurations (Figs. 8c,d) provide for the AMR tracking of a moving feature, so the effects of a moving mesh and regridding can be observed. The vorticity tag provides a more challenging test since it bisects a cubed-sphere panel edge. We see little deviation in error among the ×2, ×4, and ×8 refinement ratio simulations, so we only discuss runs with the ×4 refinement ratio. Whenever the dynamic AMR grids are moved, the underlying topography is reinitialized with the analytical formulation given in Läuter et al. (2005).
Simulations were run for 5 days. The normalized global and height errors after 5 days are shown for simulations with a c32 base grid in Table 4. Errors are calculated by comparing runs to the analytic solution of the test case. The uniform grid results are compared with one- and two-level refinement runs using the three grid configurations. For the static refinement simulations, the height error increases by approximately in comparison to the uniform c32 run. In the two-level height-tag AMR run, we observe that the height error increases by about , while the two-level vorticity-tag AMR run decreases the error by roughly . Figure 8 depicts the USBR height errors at day 5 for the c32 uniform run and c32 base grid runs with the three grid configurations. For the static equatorial patch run (Fig. 8b), the height errors remain nearly the same as for the uniform run (Fig. 8a), with only very slight increases in the large error areas on the polar panels. Along the coarse–fine boundary, no spurious grid-induced error is observed. In the height-tag AMR run (Fig. 8c) the errors along the polar-equatorial panel boundaries increase, while in the vorticity-tag AMR run (Fig. 8d), we observe that the large errors on the polar panels are reduced due to the addition of refinement over that area. With the height tagging, we do not see a similar improvement because the refined patches are over the low error areas on the equatorial panel.
After 5 days, the uniform c32 run (~313-km grid) had a normalized height error of 2.6 × 10−6. For comparison, the second-order icosahedral model by Düben et al. (2012) obtained a normalized height error of 1.05 × 10−4 at day 5 using a uniform grid with an average edge length of 240 km. Other investigations focused on results at 12 h, after which the flow features have progressed only halfway around the sphere. Pudykiewicz (2011) showed a normalized height error of ~6 × 10−6 after 12 h using a second-order icosahedral geodesic model with a 2° (~220 km) grid resolution, while our uniform c32 (2.8°) run produced a normalized height error of 3.5 × 10−7 at 12 h. These results are comparable with those obtained by the fourth-order multimoment model on a Yin–Yang grid with an effective resolution of 1.875° in Li et al. (2015). They reported a normalized height error of ~3 × 10−7 after 12 h. Additionally, the third-order multimoment method on a cubed-sphere grid with a N = 40 (2.25°) grid resolution in Chen et al. (2014) had a normalized height error of ~6.5 × 10−7 at 12 h. We are unaware of other results that use the USBR test with AMR applications.
We performed the USBR test in simulations with increasing base resolutions of up to c128 (~78-km spacing) with the two dynamic grid configurations. The normalized and height errors after 5 days are plotted in Figs. 9a and 9b, respectively. At higher base resolutions the slight improvements in the height error no longer occur for the vorticity-tag AMR as the large-scale flow features are well resolved (Fig. 9a). The fourth-order convergence is maintained for all grid configurations. In the height error plot (Fig. 9b), we observe the same slight decrease in error for the coarse c16 and c32 base resolutions for the vorticity-tag AMR and the slight increase in the height-tag AMR errors as observed earlier in Figs. 8c and 8d. However, at higher base resolutions we find a large increase in the error for the vorticity-tag AMR runs. While fourth-order convergence is maintained at all resolutions for the uniform, static refinement, and height-tag AMR configurations across all resolutions, the convergence rate drops to between 3 and 2.5 for the vorticity-tag AMR runs at higher base resolutions. The increased error is due to the regridding of the AMR patches. The maximum errors occur in cells bordering the coarse–fine boundary of the AMR patch and the base grid when that boundary intersects an edge of the cubed sphere. This point-source-like artifact of the AMR grid occurs in both the height-tag and vorticity-tag AMR simulations. In the height-tag runs, the artifact is triggered only when the AMR grid passes over the corners of the cubed sphere, resulting in the slight error increase observed in Fig. 8c. In the vorticity-tag AMR runs, the refined grid bisects the polar-equatorial panel edge during the entire run, thus this small error is compounded at each regridding step. This results in the sharp increase in the error seen in Fig. 9b for the vorticity-tag runs with c64 and c128 base resolutions. Overall, though, the error is small and very localized at the cells where the AMR patches intersect the cubed-sphere edge. It is, therefore, only obvious in the strict error measure and only at high horizontal base resolutions. At lower base resolutions the magnitudes of other errors are bigger which then dominate the global error measure.
d. Isolated mountain gravity wave
This shallow-water test was developed to assess AMR when topography is present. In the test a gravity wave, which is triggered by an unbalanced initial height perturbation in a quiescent background environment, passes over an idealized mountain. The change in topography deforms the structure of the gravity wave. The bottom topography consists of a cosine mountain and is defined by
where and . Outside the radius R the topography is set to zero. The peak height of the mountain is , and it is centered at in the longitudinal and latitudinal direction, respectively. The initial velocity field is set to zero and the initial free surface height has a constant background value of with a local Gaussian dip perturbation. Thus, the initial free surface height field (above the reference sphere at sea level) is given as
The maximum depth of the perturbation is set to , is a width parameter, and τ is the great-circle distance from point to the dip’s center such that
where is the average radius of Earth. The Gaussian dip is centered at in the longitudinal and latitudinal direction, respectively. Figure 10 depicts the initial height field and its distance from the mountain.
The simulation is run for a period of 12 h so that the gravity wave has moved halfway around the sphere. Figures 11 and 12 show the perturbation height (defined as deviations from , top panels) and the height difference from a uniform c1024 (~10 km) reference solution (bottom panels) at hours 6 and 12, respectively, for a uniform c128 run and a c32 base one-level AMR run (c32/c128) tagged on a height-gradient threshold of . After 6 h (Fig. 11), the gravity wave has just passed over the mountain and the distortion to the wave from the mountain is clearly visible. The presence of the mountain breaks the symmetry of the circular, outward-propagating gravity wave. This propagation is captured by the AMR refinement criterion as indicated by the overlaid block structure in Figs. 11b and 11d. The height differences for the c32/c128 AMR run (Figs. 11d and 12d) are similar in position and magnitude to the uniform c128 run errors (Fig. 11c) within the refined AMR domain. The areas of larger error at the borders of the AMR region seen in Figs. 11d and 12d are located over the leading and trailing edges of the gravity wave, which are not fully covered by the AMR refinement criterion. The location and magnitude of these larger errors correlate with the error at the leading edge of the gravity wave observed in the uniform c32 run. At hour 6, the AMR refined grids are over the mountain and by hour 12, the mountain is once again covered by only the coarse grid. The AMR reinitializes the topography [using Eq. (1)] whenever adaptations are triggered. No spurious errors appear as the AMR refines and coarsens over the topography region (Figs. 11d and 12d).
Results in Fig. 13 show the normalized height error (Fig. 13a) and the total number of grid cells (Fig. 13b) as a function of forecast hour for uniform runs from c32 to c512 and one-level AMR runs using height-gradient tagging with thresholds of , , and , labeled T1, T2, and T3, respectively. The normalized error metrics are determined with respect to the uniform c1024 simulation, which serves as the reference solution. For uniform runs, the solution error converges to fourth order in both the normalized and height error norms. The AMR runs have improved error but do not reach the error of the uniform run with the same resolution as the highest refinement level. The c32 base one-level AMR T3 run (c32/c128) has a maximum number of grid cells roughly equivalent to the uniform c64 run, but its error is nearly an order of magnitude smaller than the uniform c64 run. Reducing the gradient threshold in AMR runs so that more area around the gravity wave is covered by the AMR patches improves the solution and reduces error. A c32/c128 AMR run with a refinement threshold of , lower than the T3 criterion, results in the AMR grid covering of the globe by hour 12 and a normalized height error of 3.957 × 10−6. In comparison, the uniform c128 run has an error of 3.051 × 10−6 and the T3 run has an error of 5.806 × 10−6 with AMR blocks covering only of the area. As more of the leading edge of the gravity wave is refined with lower tagging thresholds, the error is decreased further, though at a diminishing rate.
e. Binary–vortices interaction
The binary–vortices interaction test demonstrates the AMR benefits and its effectiveness in studying an important and more realistic problem. The interaction of two neighboring tropical cyclones (TCs) often alters the structures of the two, leads to complex tracks for the storms, and in some instances results in a merger of the two cyclones. These interactions were first studied by Fujiwhara (Fujiwhara 1921) and are commonly called the Fujiwhara effect. Idealized binary–vortex interactions have been extensively investigated using 2D idealized models by Melander et al. (1988), Waugh (1992), Ritchie and Holland (1993), Prieto et al. (2003), and Shin et al. (2006). A majority of the research has been conducted on two-dimensional Cartesian systems using a constant Coriolis force. These studies have used a variety of initial vortex profiles featuring discrete (Ritchie and Holland 1993) or continuous (Bauer et al. 2014) vortices in both symmetric and asymmetric pairs (Dritschel and Waugh 1992). They demonstrate that slight changes in initial conditions will cause widely diverging results. The vortices will either merge or repel each other depending on the strength, size, and separation distances of the vortices, and the postinteraction shapes of the vortices will be vastly different. Bauer et al. (2014) used an r-adaptive shallow-water model to demonstrate that the vortices’ tracks are sensitive to initial conditions and to initial grid resolution. Given the sensitivity to resolution, binary–vortex interactions are a well-suited test of AMR. The steep gradients, localized areas of high vorticity, and complex flow fields around the vortices are transient and resolution-dependent, mimicking the multiscale nature of tropical cyclones. With this test, we can evaluate the AMR’s ability to refine and track these features of interest and measure the AMR’s accuracy in resolving the vortex interaction. It can assess how well the results and errors in AMR runs match the results of uniform high-resolution runs and can determine the sensitivity of the vortex tracks to the changing grid resolutions.
In our binary tropical cyclone–like vortices test, we use the full shallow-water equations on a spherical grid with a changing Coriolis parameter, whereas most other published studies use a nondivergent barotropic model on an f plane. We also restrict our study to only the symmetric case so that the two vortices are identical in size and strength. The two vortices are initialized near each other and are allowed to interact over a simulation period of several days. Two variations of this setup are presented. In the separation case, the two vortices orbit around each other and then slowly drift apart. In the merger case, the vortices merge. Our initializations of the continuous vortex profiles were inspired by the definition of the initial state in Holland and Dietachmayer (1993). The initial wind and height profiles are derived from the shallow-water equations in cylindrical coordinates using an f-plane approximation. The vortex structure is depicted as a radial perturbation in the geopotential field and is given by
where is the background geopotential with the background height m and Earth’s gravity m s−2, denotes the geopotential perturbation, symbolizes the maximum geopotential perturbation, is the radius of maximum wind, b stands for a scaling parameter set to 1.5, and r is the great-circle distance from point to the vortex center [see also Eq. (3)]. The values of , , , and are provided later.
A balanced tangential wind field is then found by using the steady-state shallow-water momentum equations in cylindrical coordinates. In particular, the tangential wind:
where f is the constant Coriolis parameter for an f-plane approximation at the latitude of the vortex center (specified later). The last initialization step is to project the tangential velocity onto the sphere with the zonal u and meridional υ spherical wind components given by
The threshold value prevents division by zero. The topography is set to zero.
This initialization technique represents a perfect balance for a single vortex in cylindrical coordinates and leads to a very good balance in spherical coordinates. Note that the perfect balance is broken on the sphere since f varies in the spherical domain and is held constant for the purpose of the initialization. In addition, an analytically derived balance is not fully balanced in a numerical (discrete) sense and the overlap region of two vortices is not strictly balanced either. However, the initial imbalances for our separation and merger test cases are very minor and the resulting small gravity waves do not interfere with our AMR analysis. The same initialization technique was also used for idealized tropical cyclone simulations in 3D GCMs (Reed and Jablonowski 2011).
The radial cross sections of the initial relative vorticity, height field, and tangential wind for a single vortex as a function of the great-circle distance from the center are depicted in Fig. 14a. Here, all magnitudes are normalized to 1 and are provided below for each test case. The initial relative vorticity for the merger test case with two initial vortices is depicted in Fig. 14b. It shows that the tangential wind profile from Eq. (7) results in a relative vorticity profile with a core of positive relative vorticity in the center of each storm surrounded by a broad ring of negative vorticity with relatively small magnitudes. The two vortices slightly overlap with very minor magnitudes of u, υ, and [Eq. (5)]. Here, we use the sum of the u, υ, and values of both vortices to initialize our shallow-water system.
1) Separation case
In the separation case, the two vortices are centered at N, which defines the constant Coriolis parameter with Earth’s angular velocity s−1. The maximum geopotential perturbation is set to with the maximum height perturbation m. The radius of maximum wind is set to km. This results in a maximum tangential wind of 64 m s−1 and a maximum relative vorticity of 9.4 × 10−4 s−1. The centers of the two vortices are 13.5° apart ( km) in the longitudinal direction, 6 times the radius of maximum wind, so that their negative vorticity regions still overlap. In particular, the two vortex centers are located at and with the midway point between the two cyclones at 90°W.
This scenario is sensitive to variations in initial conditions, making it desirable for testing adaptive grids. Decreasing the separation distance by 20 km results in the merger of the two vortices. During the first 3 days of the simulation, the two vortices make one complete orbit around each other as the beta-drift steers them toward the northwest, after which the two vortices then drift apart. In that time, the negative vorticity is stretched as it is advected around the pair of positive cores before being spun off behind the pair as an anticyclone. We also observed the growth of a large-scale wave train that forms in the lee of the orbiting pair. The time evolution of the flow can be seen in Fig. 15, which shows the vorticity field of the cyclone pair at days 1, 2, 3, 4, and 6 for several uniform-resolution runs. These serve as references for the AMR simulations. We ran uniform runs with resolutions from c32 to c1024, of which the c128, c256, and c1024 runs are depicted in Fig. 15. Results vary significantly with resolution, though results do converge with increasing resolution. Runs with coarser resolution than c128 (not shown) have very weak vortices that merge instead of drifting apart. In the c128 run (Figs. 15k–o), the vortices start separating earlier and at the end of the run are in markedly different positions. The c256 simulation (Figs. 15f–j) more closely resembles the highest-resolution c1024 run (Figs. 15a–e), but there are still significant differences. However, the c512 run (not shown as a time series) is nearly indistinguishable from c1024 with only slight differences in the center of the vortex cores and in the finescale vorticity filaments. A comparison of the uniform c512 run and c1024 run at day 6 can be seen in Figs. 16e and 16i.
To assess the AMR performance, we ran the model using relative vorticity refinement criteria with one and two levels of AMR and ×4 refinement on base grid resolutions from c16 to c256. Samples of the resulting relative vorticity fields at day 6 using two different relative vorticity thresholds of and are displayed in Fig. 16. They are divided into columns that share the same highest resolutions; for example, in the leftmost column, Fig. 16a is the uniform c256 run, while Figs. 16b–d are for one- or two-level AMR runs that have a finest grid resolution of c256. The AMR runs agree well with the uniform solutions having the same resolution as the finest AMR level. The AMR blocks accurately capture the positions of the two vortices and the shape of their high-vorticity cores. The AMR runs also effectively reproduce the overall shapes of the anticyclonic filaments and patches around the cores, and the wave train developing to the lee of the vortices with only minor differences in the anticyclone filament between the two vortices in a few runs. Given the nonlinear sensitivity to initial conditions and grid resolution, unrefined areas and the coarse–fine boundaries near the vortices in AMR runs may cause divergent solutions. In none of our AMR simulations did this occur. The results of the AMR runs differed only slightly from the uniform-resolution runs.
2) Merging vortices
In the merging-vortices case, the two vortices have the same maximum height perturbation with m, but the radius of maximum wind is increased to km so that the maximum tangential wind is 61 m s−1 and the maximum relative vorticity is around 5.7 × 10−4 s−1. The vortices are centered at N where the constant Coriolis parameter f is evaluated for this test case. The vortex centers are 15.65° apart (~1700 km) in the longitudinal direction. In particular, they are located at and with the same midway point as before at 90°W. This separation distance is about 4.3 times the radius of maximum wind, so that the edges of the positive vorticity cores slightly overlap.
Figure 17 depicts the evolution of these vortices as they merge over the course of four simulation days. As in the previous case, though small changes in initial condition lead to very different results, we do observe a slow convergence with increasing resolution. We ran the test case with several configurations using one and two levels of AMR with criteria based on a relative vorticity or a height gradient threshold. The AMR run shown in Fig. 17 has a c64 base resolution with two levels of ×4 refinement so that its finest level has a c1024 resolution. The AMR is triggered when the absolute value of the relative vorticity exceeds 2.3 × 10−5 s−1. With that relatively low refinement threshold, the AMR captures not only the main vortex cores, but also the small-scale anticyclonic filament that extends far south of the merged vortex, and the small-magnitude wave train that develops by day 4. Figure 18 shows a column comparison of the vorticity field at day 4 for uniform-resolution runs and AMR runs using three refinement criteria. The first column contains the uniform c1024 run and 3 two-level AMR runs with ×4 refinement on c64 base resolution that have a finest resolution of c1024. The second column has the c512 uniform run and 3 two-level AMR runs with a c32 base. The last column has the uniform c256 run and 3 two-level AMR runs with a c16 base. The three AMR refinement criteria in Fig. 18 are as follows:
Large-height-gradient-tag AMR: tag where the absolute value of the height gradient is (second row of Fig. 18);
Small-height-gradient-tag AMR: tag where (third row);
Vorticity-tag AMR: tag where the absolute value of the relative vorticity is s−1 (fourth row).
Locally refining the grid resolution with AMR effectively achieves a similar result in the refined areas as the corresponding high-resolution uniform runs. Even the large-height-gradient refinement threshold used in Figs. 18b and 18f, which results in very little refinement, is still able to produce a very similar vortex structure and position within the refined area demonstrating little to no negative effects from the coarse–fine boundaries surrounding the vortex. The lower refinement thresholds are further able to capture the anticyclonic filament wrapping around the new vortex and extending down from it as well as the development of the secondary leeside wave train.
The convergence to the error of the uniform-resolution runs can be observed in the normalized vorticity error seen in Fig. 19. The error norm is computed by comparing runs to the uniform c1024 run, which serves as a reference solution. The figure depicts the normalized relative vorticity error (Fig. 19a) and total number of grid cells (Fig. 19b) as a function of forecast days for uniform-resolution runs and one- and two-level AMR runs using the large-height-gradient tag or the vorticity tag with ×4 refinement ratios. The vorticity-tag AMR simulations (both one- and two-level AMR) have nearly the same error as the uniform runs with the highest resolution while the gradient-tag runs have slightly higher error. This agrees with the fact that the gradient-tag runs use fewer grid cells and only cover the merged vortex core (see Figs. 18b,f,j). Although the large-scale shape and locations of the two merging vortices and the postmerger vortex appear visually to converge to a solution with increasing resolution, we do not observe a large reduction in the global errors with increasing resolution. The source of this error is the small differences that occur in the core of the vortices, caused by small-scale nonlinear features in the high-vorticity filaments, as well as slight variations in the beta drift created by small changes in the vorticity magnitude for the different resolution runs. These small differences in these features lead to localized large-magnitude errors in the vorticity. As in the separation case, the AMR improves the solution using fewer grid cells. Even when the AMR patch over the vortex is small and the coarse–fine boundary is near the high vorticity cores, the solution is not negatively distorted, showing the robustness of the model given the sensitivity of the test case to grid resolution and slight changes in initial conditions.
The computing run times versus number of grid cells for a 4-day simulation with vorticity tagging is presented in Table 5. The table thereby represents some of the runs from Fig. 19b. Eight processors on one node of NCAR’s Yellowstone computing platform are used. We see the approximate 8× reduction in cost between the c512 and c256 uniform runs, as expected for a doubling of the horizontal resolution and a halving of the time step. For this test, the wall-clock run time for AMR runs is closer to for 2× resolution changes, demonstrating some of the overhead of the AMR algorithm. The total wall-clock time roughly correlates with the total number of grid cells, as in the moving-vortices advection test, even for the coarsest two-level AMR runs (c32/c128/c512 and c16/c64/c256).
In this paper, we utilized a fourth-order finite-volume model on a cubed-sphere grid, which is adaptive in both space and time, to demonstrate the effectiveness of the AMR in resolving and tracking chosen features of interest while maintaining large-scale smooth flows. Using selected shallow-water and advection test cases, we evaluated the AMR’s ability to track and resolve features of interest without creating distortions or numerical noise in the large-scale smooth flows at the interfaces between meshes. A variety of static and dynamic refinement criteria and strategies are implemented to assess the strengths and weaknesses of the AMR method. With the large-scale smooth “do no harm” tests, one and two levels of static and adaptive refinement meshes with several refinement ratios were placed at several locations on the cubed-sphere grid. The results confirmed that multiple levels of refinement and abrupt ×4 or ×8 refinement ratios between levels still allowed flows to move smoothly through the refined areas. There was little induced noise and numerical error at the refinement boundaries. For coarse resolutions, the refinement improved global errors slightly, and the errors remained nearly unchanged when refinement was added to a higher-resolution base grid for the two shallow-water tests. Only for high resolutions in the USBR tests when a moving AMR grid transected the cubed-sphere panel boundaries did we see a noticeable increase in error. This error, however, was very localized and only becomes apparent because the base global error is so low in the uniform-resolution simulations for this smooth idealized test. In the coarser runs for the USBR and in the other more complex shallow-water tests with larger expected global errors, this was not observed.
With the three AMR test cases, we demonstrated that AMR is able to track the features of interest and closely reproduces the results of uniform high-resolution runs using fewer grid cells. AMR was implemented using tracer and height-field gradients as well as relative-vorticity magnitude as tagging criteria with multiple refinement levels and a range of thresholds. The AMR grids are added and removed in time without creating significant distortions or noise at the mesh interfaces. In the tracer advection test, fourth-order convergence was maintained while using AMR, and the error of AMR runs with one level of refinement was comparable to the error of uniform runs having the same fine-level resolution. The test showed the importance of refining early, as errors developed from the coarse grid propagate through the fine grids for the rest of the simulation. The gravity wave impinging on the mountain test demonstrated the use of AMR with topography. Refinement was added and removed from the areas with topography without creating additional negative impacts. In the binary vortices test, AMR improved accuracy of the position of the vortices as they interacted and their structures when compared with uniform-resolution runs. Even stringent criteria with high threshold values, which did not create a large buffer of high resolution around the vortices, still produced accurate results and improved the solution in the refined patches. Additional refinement, though, significantly improved the representation of the vorticity filaments that extended well away from the central vortices and the developing Rossby wave train.
All three test cases demonstrated that a variety of AMR criteria and thresholds lead to improvements in the results, though to maximize that improvement, the refinement criteria needed careful tailoring. Several conditions increased the effectiveness of AMR; however, there was no clear strategy for establishing the best general refinement criteria. Having initial refinement or refining early in the run before errors developed on the coarse grid was one of the key strategies for improving accuracy. When there is no initial refinement, the benefits of AMR are limited by the coarseness of the base mesh; AMR with two levels was ineffective due to large errors introduced by the coarsest base meshes early in the calculation. This speaks to the need for a sufficiently refined base mesh to avoid contaminating finer levels. Using more than one level of refinement and effective tagging strategies resulted in better-resolved features of interest, but at a diminishing rate of return of improvement. Our conclusion is that the benefit of AMR does not come automatically from the computational savings of a very coarse base mesh. However, there may still be benefits of two or more levels over uniform-resolution calculations that otherwise would not be computationally feasible without AMR. In a realistic climate simulation, tropical cyclones could, perhaps, be effectively captured early by criteria that place high resolution over the cyclogenesis region and then refine over and track emerging storms to ensure continued accuracy. For other, more complex or moist flows, more advanced criteria than just a simple relative-vorticity threshold need to be investigated. They could be based, for example, on combinations of physics-based properties (like rainfall), thresholds of vorticity, or gradients. Future work will explore such refinement criteria in the 3D nonhydrostatic version of the Chombo-AMR model with and without a variety of physical parameterization schemes. The addition of physical parameterizations will also allow us to test the scale awareness of the physics schemes.
Support for this work has been provided by the Office of Science, U.S. Department of Energy, Award DE-SC0003990 and by the Director, Office of Science, Office of Advanced Scientific Computing Research of the U.S. Department of Energy under Contract DE-AC02-05CH11231 as part of their Mathematical, Computational, and Computer Sciences Research/Computational Partnerships Program. We would like to acknowledge high-performance computing support from Yellowstone provided by NCAR’s Computational and Information Systems Laboratory, sponsored by the National Science Foundation. We thank the reviewers for their helpful comments and suggestions.