A High-Order Discontinuous Galerkin Solver with Dynamic Adaptive Mesh Reﬁnement to Simulate Cloud Formation Processes

We present a high-order discontinuous Galerkin ( dg ) solver of the compressible Navier-Stokes equations for cloud formation processes. The scheme exploits an underlying parallelized implementation of the ader-dg method with dynamic adaptive mesh reﬁnement. We improve our method by a pde -independent general reﬁnement criterion, based on the local total variation of the numerical solution. Our generic scheme shows competitive results for both classical cfd and stratiﬁed scenarios. While established methods use numerics tailored towards the speciﬁc simulation, our scheme works scenario independent. All together, our method can be seen as a perfect candidate for large scale cloud simulation runs on current and future super computers.


Introduction
In this paper we address the resolution of basic cloud formation processes on modern super computer systems.The simulation of cloud formations, as part of convective processes, is expected to play an important role in future numerical weather prediction [1].This requires both suitable physical models and effective computational realizations.Here we focus on the simulation of simple benchmark scenarios [10].They contain relatively small scale effects which are well approximated with the compressible Navier-Stokes equations.We use the aderdg method of [5], which allows us to simulate the Navier-Stokes equations with a space-time-discretization of arbitrary high order.In contrast to Runge-Kutta time integrators or semi-implicit methods, an increase of the order of ader-dg only results in larger computational kernels and does not affect the complexity of the scheme.Additionally, ader-dg is a communication avoiding scheme and reduces the overhead on larger scale.We see our scheme in the regime of already established methods for cloud simulations, as seen for example in [10,12,13].
Due to the viscous components of the Navier-Stokes equations, it is not straightforward to apply the ader-dg formalism of [5], which addresses hyperbolic systems of partial differentials equations (pdes) in first-order formulation.
To include viscosity, we use the numerical flux for the compressible Navier-Stokes equations of Gassner, Lörcher, and Munz [8].This flux has already been applied to the ader-dg method in [4].In contrast to this paper, we focus on the simulation of complex flows with a gravitational source term and a realistic background atmosphere.Additionally, we use adaptive mesh refinement (amr) to increase the spatial resolution in areas of interest.This has been shown to work well for the simulation of cloud dynamics [12].Regarding the issue of limiting in highorder dg methods, we note that viscosity not only models the correct physics of the problem but also smooths oscillations and discontinuities, thus stabilizing the simulation.
We base our work on the ExaHyPE-Engine (www.exahype.eu),which is a framework that can solve arbitrary hyperbolic pde systems.A user of the engine is provided with a simple code interface which mirrors the parts required to formulate a well-posed Cauchy problem for a system of hyperbolic pdes of first order.The underlying ader-dg method, parallelization techniques and dynamic adaptive mesh refinement are available for simulations while the implementations are left as a black box to the user.An introduction to the communicationavoiding implementation of the whole numerical scheme can be found in [3].
To summarize, we make the following contributions in this paper: -We extend the ExaHyPE-Engine to allow viscous terms following the aforementioned numerical scheme.-We thus provide an implementation of the compressible Navier-Stokes equations.In addition, we tailor the equation set to stratified flows with gravitational source term.We emphasize that we use a standard formulation of the Navier-Stokes equations as seen in the field of computational fluid mechanics and only use small modifications of the governing equations, in contrast to a equation set that is tailored exactly to the application area.-We present a general amr-criterion that is based on the detection of outlier cells w.r.t.their total variation.Furthermore, we show how to utilize this criterion for stratified flows.-We evaluate our implementation with standard cfd scenarios and atmospheric flows and inspect the effectiveness of our proposed amr-criterion.
We thus inspect, whether our proposed general implementation can achieve results that are competitive with the state-of-the-art models that rely on heavily specified equations and numerics.

Equation Set
The compressible Navier-Stokes equations in the conservative form are given as with the vector of conserved quantities Q, flux F (Q, ∇Q) and source S(Q).Note that the flux can be split into a hyperbolic part F h (Q), which is identical to the flux of the Euler equations, and a viscous part F v (Q, ∇Q).The conserved quantities Q are the density ρ, the two or three-dimensional momentum ρv and the energy density ρE.The rows of Eq. ( 1) are the conservation of mass, the conservation of momentum and the conservation of energy.
The pressure p is given by the equation of state of an ideal gas The term gz is the geopotential height with the gravity of Earth g [10].The temperature T relates to the pressure by the thermal equation of state where R is the specific gas constant of a fluid.We model the diffusivity by the stress tensor with constant viscosity µ.The heat diffusion is governed by the coefficient where the ratio of specific heats γ, the heat capacity at constant pressure c p and the Prandtl number Pr depend on the fluid.Many realistic atmospheric flows can be described by a perturbation over a background state that is in hydrostatic equilibrium i.e. a state, where the pressure gradient is exactly in balance with the gravitational source term S ρv = −kρg.The vector k is the unit vector pointing in z-direction.The momentum equation is dominated by the background flow in this case.Because this can lead to numerical instabilities, problems of this kind are challenging and require some care.To lessen the impact of this, we split the pressure p = p + p into a sum of the background pressure p(z) and perturbation p (x, t).We split the density ρ = ρ + ρ in the same manner and arrive at Note that a similar and more complex splitting is performed in [12,10].In contrast to this, we use the true compressible Navier-Stokes equations with minimal modifications.

Numerics
The ExaHyPE-Engine implements an ader-dg-scheme and a muscl-Hancock finite volume method.Both can be considered as instances of the more general PnPm schemes of [5].We use a Rusanov-style flux that is adapted to pdes with viscous terms [8,7].The finite volume scheme is stabilized with the van Albada limiter [15].The user can state dynamic amr rules by supplying custom criteria that are evaluated point-wise.Our criterion uses an element-local error estimate based on the total variation of the numerical solution.We exploit the fact that the total variation of a numerical solution is a perfect indicator for edges of a wavefront.Let f (x) : R Nvars → R be a sufficiently smooth function that maps the discrete solution at a point x to an arbitrary indicator variable.The total variation (tv) of this function is defined by for each cell.The operator • 1 denotes the discrete L 1 norm in this equation.
We compute the integral efficiently with Gaussian quadrature over the collocated quadrature points.How can we decide whether a cell is important or not?To resolve this conundrum, we compute the mean and the population standard deviation of the total variation of all cells.It is important that we use the method of [2] to compute the modes in a parallel and numerical stable manner.A cell is then considered to contain significant information if its deviates from the mean more than a given threshold.This criterion can be described formally by The parameters T refine > T coarsen can be chosen freely.Chebyshev's inequality with probability P guarantees that we neither mark all cells for refinement nor for coarsening.This inequality holds for arbitrary distributions under the weak assumption that they have a finite mean µ and a finite standard deviation σ [16].Note that subcells are coarsened only if all subcells belonging to the coarse cell are marked for coarsening.In contrast to already published criteria which are either designed solely for the simulation of clouds [12] or computationally expensive [7], our criterion works for arbitrary pdes and yet, is easy to compute and intuitive.

Results
We use a mix of various benchmarking scenarios that consist of standard cfd scenarios and atmospheric flows.For the standard cfd testing scenarios, we use -0.5 -0.25 0 0.25 0.5 x, z 0.0 0.5 (a) Our approximation (solid lines) of the lid-driven cavity flow vs. reference solution (crosses) of [9].The respective other coordinate is held constant at a value of 0.
(b) Our result (markers) of the Taylor-Green vortex vs. the analytical solution (lines) Eq. ( 12).The plot shows two velocity slices, the respective other coordinate is held constant at a value of π.
The constant C = 100 /γ governs the speed of sound and thus the Mach number Ma = 0.1 [6].The viscosity is set to µ = 0.1.We simulate on the domain [0, 2π] 2 and impose the analytical solution at the boundary.A comparison at time t = 10.0 of the analytical solution for the pressure with our approximation (fig.1b) shows excellent agreement.We used an ader-dg-scheme of order 5 with a grid of 25 2 cells.
The Arnold-Beltrami-Childress (abc) flow is similar to the Taylor-Green vortex but possesses an analytical solution for three-dimensions as well [14].It is defined in the domain [−π, π]  13)).All other axes are held constant at a value of 0. We show every 6th value.again under the assumption of incompressibility.The constant C = 100 /γ is chosen as before.We use a viscosity of µ = 0.01 and analytical boundary conditions.Our results (fig.2) show a good agreement between the analytical solution and our approximation with an ader-dg-scheme of order 3 with a mesh consisting of 27 3 cells at time t = 0.1 s.As a final example of standard flow scenarios, we consider the lid-driven cavity flow where the fluid is initially at rest, with ρ = 1 and p(x) = 100 /γ.We consider a domain of size 1 m × 1 m which is surrounded by no-slip walls.The flow is driven entirely by the upper wall which has a velocity of v x = 1 m/s.The simulation runs for 10 s.Again, our results (fig.1a) have an excellent agreement with the reference solution of [9].We used an ader-dg-method of order 3 with a mesh of size 27 2 .
With the constants γ = 1.4,Pr = 0.71, R = 287.058,p 0 = 10 × 10 5 Pa, g = 9.8 m/s 2 , ( all following stratified flow scenarios are described in terms of the potential temperature with reference pressure p 0 [12,10].We compute the initial background density and pressure by inserting the assumption of a constant background energy in Eq. ( 6).The background atmosphere is then perturbed.We set the density and energy at the boundary such that it corresponds to the background atmosphere.Furthermore, to ensure that the atmosphere stays in hydrostatic balance, we need to impose the viscous heat flux at the boundary [10].In this equation, T (z) is the background temperature at position z, which can be computed from Eqs. ( 2) and (6).Our first scenario is the colliding bubbles scenario [12].We use perturbations of the form where s is the decay rate and r is the radius to the center i.e., r denotes the Euclidean distance between the spatial positions x = (x, z) and the center of a bubble x c = (x c , z c ) -for three-dimensional scenarios x and x c also contain a y coordinate.
We have two bubbles, with constants warm: (19) Similar to [12], we use a constant viscosity of µ = 0.01 to regularize the solution.Note that we use a different implementation of viscosity than [12].Hence, it is difficult to compare the parametrization directly.We ran this scenario twice: once without amr and a mesh of size 1000/81 m = 12.35 m and once with amr with two adaptive refinement levels and parameters T refine = 2.5 and T coarsen = −0.5.We specialize the amr-criterion (Eq.9) to our stratified flows by using the potential temperature.This resulted in a mesh with cell-size lengths of approx.111.1 m, 37.04 m, and 12.34 m.The resulting mesh can be seen in fig. 4. We observe that the L 2 difference between the potential temperature of the amr run, which uses 1953 cells, and the one of the fully refined run with 6561 cells, is only 1.87.The relative error is 6.69 × 10 −6 .We further emphasize that our amr-criterion accurately tracks the position of the edges of the cloud instead of only its position.This is the main advantage of our gradient-based method in contrast to methods working directly with the value of the solution, as for example [12].Overall, our result for this benchmark shows an excellent agreement to the previous solutions of [12].
In addition, we simulated the same scenario with our muscl-Hancock method, using 7 2 patches with 90 2 finite volume cells each.As we use limiting, we do not need any viscosity.The results of this method (fig.3) also agree with the reference but contain fewer details.Note that the numerical dissipativity of the finite volume scheme has a smoothing effect that is similar to the smoothing caused by viscosity.
For our second scenario, the cosine bubble, we use a perturbation of the form For the three-dimensional bubble, we set y c = x c = 500 m.This corresponds to the parameters used in [11] 1 .For the 2D case, we use a constant viscosity of µ = 0.01 and an ader-dg-method of order 5 with two levels of dynamic amr, resulting again in cell sizes of roughly 111.1 m, 37.04 m, 12.34 m.We use slightly different amr parameters of T refine = 1.5 and T coarsen = −0.5 and let the simulation run for 600 s.Note that, as seen in fig.5a, our amr-criterion tracks the wavefront of the cloud accurately.This result shows an excellent agreement to the ones achieved in [10,12].For the 3D case, we use an ader-dg-scheme of order 3 with a static mesh with cell sizes of 40 m and a shorter simulation duration of 400 s.Due to the relatively coarse resolution and the hence increased aliasing errors, we need to increase the viscosity to µ = 0.05.This corresponds to a larger amount of smoothing.Our results (fig.5b) capture the dynamics of the scenario well and agree with the reference solution of [11].

Conclusion
We presented an implementation of an ader-dg-method and a muscl-Hancockscheme with amr for the Navier-Stokes equations, based on the ExaHyPE-Engine.Our implementation is capable of simulating different scenarios: We evaluated our method for standard cfd scenarios and achieved competitive results for all used scenarios.This is true for both two-dimensional scenarios (Taylor-Green vortex and lid-driven cavity) and for the three-dimensional abc-flow.Furthermore, our method allows us to simulate flows in hydrostatic equilibrium correctly, as our results for the cosine and colliding bubble scenarios showed.We showed that our amr-criterion is able to vastly reduce the number of grid cells while preserving the quality of the results.

Figure 2 :
Figure 2: Our approximation (markers) of the abc-flow vs. analytical solution (lines, Eq. (13)).All other axes are held constant at a value of 0. We show every 6th value.