# Non-hydrostatic mesoscale atmospheric modeling by the anisotropic mesh
adaptive discontinuous Galerkin method111This work was supported by grant No.
20-01074S of the Czech Science Foundation. The author acknowledges also the
membership in the Nečas Center for Mathematical Modeling
ncmm.karlin.mff.cuni.cz.

Vít Dolejší<EMAIL_ADDRESS>Charles University, Faculty of
Mathematics and Physics, Sokolovská 83, 186 75 Praha, Czech Republic

###### Abstract

We deal with non-hydrostatic mesoscale atmospheric modeling using the fully
implicit space-time discontinuous Galerkin method in combination with the
anisotropic $hp$-mesh adaptation technique. The time discontinuous
approximation allows the treatment of different meshes at different time
levels in a natural way which can significantly reduce the number of degrees
of freedom. The presented approach generates a sequence of triangular meshes
consisting of possible anisotropic elements and varying polynomial
approximation degrees such that the interpolation error is below the given
tolerance and the number of degrees of freedom at each time step is minimal.
We describe the discretization of the problem together with several
implementation issues related to the treatment of boundary conditions,
algebraic solver and adaptive choice of the size of the time steps. The
computational performance of the proposed method is demonstrated on several
benchmark problems.

###### keywords:

non-hydrostatic mesoscale atmospheric problems , space-time discontinuous
Galerkin method , anisotropic $hp$-mesh adaptation , numerical study

###### MSC:

65M60 , 65M50 , 76M10

##  1 Introduction

The numerical solution of the compressible Euler and/or Navier-Stokes
equations has become the standard tool for the simulation of the atmosphere in
the last decades, see e.g. the reviews in [30] or [14]. Due to the large scale
simulation, it is advantageous to use higher order discretization schemes
since they significantly reduce the number of degrees of freedom
($\mathrm{DoF}$) needed to achieve the required accuracy. It is possible to
use the spectral element (SE) method [16, 18] or the finite volume (FV)
approach with higher-order reconstruction schemes as ENO, WENO or TENO, cf.
[23] and the references cited therein. In the recent years, large attention
has been devoted to the discontinuous Galerkin (DG) methods which employ the
advantageous properties of SE and FV methods, cf. [26, 33, 3, 32, 21] for
example. DG methods employ piecewise polynomial but discontinuous
approximation, they provide high accuracy and flexibility, local conservation,
stability, and parallel efficiency. However, the disadvantage is a large
number of $\mathrm{DoF}$ in comparison to FV and SE methods.

Due to the high-order spatial approximation, the use of an explicit time
discretization is very inefficient, since the size of the time steps is
strongly restricted by the stability condition, cf. [17]. Therefore, several
variants of semi-implicit time discretization techniques have been developed,
e.g., the implicit-explicit (IMEX) formulation in [26, 15], its recent
improvement in [25], the IMEX evolution Galerkin method in [22, 33], the
positivity-preserving well-balanced central discontinuous Galerkin method in
[21] or the techniques based on the Eulerian-Lagrangian approach in [3].

In the last two decades, we have developed, analyzed and tested the
discontinuous Galerkin method for solving the compressible Navier-Stokes
equations with applications in aerodynamics, including supersonic, transonic,
subsonic, and low-Mach number flows, cf. [8]. The prominent scheme seems to us
to be the fully implicit space-time discontinuous Galerkin (STDG) method,
which uses a discontinuous piecewise polynomial approximation with respect to
both spatial and temporal coordinates. In our experience, the time
discontinuous Galerkin method is more stable than IMEX techniques and it can
naturally handle different meshes at different time levels. However, these
advantages are compensated by a higher number of $\mathrm{DoF}$ compared to
multi-step or diagonally implicit Runge-Kutta methods.

A fully implicit time discretization leads to the need to solve a nonlinear
algebraic system at each time step. We employ an iterative Newton-like method
whose first step corresponds to our earlier semi-implicit scheme from [6].
Therefore, the additional iterative steps improve the accuracy (and stability)
of the computation. In practice, only a few iterative steps are necessary to
reduce the algebraic errors.

In the last decade, we have dealt with the anisotropic $hp$-mesh adaptation
($hp$-AMA) method which constructs meshes consisting of triangles of different
size, shape and orientation together with varying polynomial approximation
degrees for each mesh element. The idea is to minimize the number of
$\mathrm{DoF}$ such that an estimate of the interpolation error is below the
given tolerance. This anisotropic mesh adaptation technique is more general
than the approach from [32] which allows the movement of grid nodes but not
local refinement/coarsening. On the other hand, the adaptive technique from
[24] offers the local mesh refinement/de-refinement but the initial mesh
structure is fixed and then the capturing of phenomena that are not parallel
to the boundary domain is not optimal. The theoretical as well as practical
issues of the $hp$-AMA method are summarized in [9].

The STDG discretization can naturally handle grids generated by the $hp$-AMA
method since it is a one-step method and the approximate solutions at
different time levels are weakly connected together by a penalty term. In
[11], we developed the $hp$-AMA approach for solving time-dependent problems
and demonstrated its performance by solving the Burgers equation, the
isentropic vortex propagation, the Kelvin-Helmholtz instability and single-
ring infiltration in porous media flow. The computational benefit of
anisotropic meshes for the numerical solution of unsteady flow problems has
been demonstrated in many works, e.g., [2, 1, 13, 19] and the references cited
therein.

The goal of this paper is to apply the STDG and $hp$-AMA methods to the
numerical simulation of non-hydrostatic mesoscale atmospheric problems. The
presented algorithm controls the errors arising from an inexact solution of
nonlinear algebraic systems and chooses adaptively the size of the time steps
whereas the polynomial approximation degree with respect to time is kept
fixed. Hence, we can speak of $hp\tau$-adaptation. Moreover, at each time
step, the algorithm evaluates an estimate of the interpolation error. If the
error estimate exceeds the prescribed tolerance, a completely new (possibly
anisotropic) mesh with varying polynomial approximation degrees with respect
to the space is generated. Attention is paid to the realization of the
boundary conditions and the choice of the physical quantity to be controlled.
The particular aspects of the presented numerical scheme are justified by
numerical experiments.

The outline of the rest of the paper is as follows. In Section 2, we present
the governing equations with the constitutive relations and the
initial/boundary conditions. In Section 3, we introduce the discretization of
the problem including the realization of the boundary conditions. Furthermore,
Section 4 presents the solution strategy of the arising algebraic system with
an adaptive stopping criterion and the adaptive choice of the time step.
Section 5 contains a description of the $hp$-AMA technique and its application
to atmospheric modeling. The performance of the proposed method is
demonstrated in Section 6 by numerical experiments for several benchmark
examples from [14]. Finally, we give several concluding remarks in Section 7.

##  2 Problem formulation

Let ${\Omega}\subset\mathbb{R}^{2}$ be the computational domain, $T>0$ the
physical time to be reached, and $Q_{T}:={\Omega}\times(0,T)$. We use the
following notation: $\rho$ is the density, $\mathrm{p}$ is the pressure,
$\bm{v}=(v_{1},v_{2})^{\mathrm{T}}$ is the velocity vector, $\theta$ is the
absolute temperature and $e=\rho c_{\mathrm{v}}\theta+\rho|\bm{v}|^{2}/2$ is
the energy per unit volume including the interior and kinetic energies (and
excluding the gravitational energy). Moreover, $c_{\mathrm{v}}>0$ and
$c_{\mathrm{p}}>0$ are the specific heat capacities at constant volume and
pressure, respectively, $\kappa=c_{\mathrm{p}}/c_{\mathrm{v}}>1$ is the
Poisson adiabatic constant, $R=c_{\mathrm{p}}-c_{\mathrm{v}}$ is the gas
constant, $\mu$ is the dynamic viscosity,
$g=9.81\,\mathrm{m}\cdot\mathrm{s}^{-2}$ is the gravity constant,
${\bm{k}}=(k_{1},k_{2})^{\mathrm{T}}$ is the upward pointing unit vector,
$\Pr$ is the Prandtl number, and $\tau_{ij},\ i,j=1,2$ denote the components
of the viscous part of the stress tensor. The symbols $\partial_{t}{}$ and
$\partial_{x_{i}}{}$ denote the partial derivatives with respect to $t$ and
$x_{i},\ i=1,2$, respectively.

We consider the system of the compressible Navier-Stokes equations with the
gravity forces written in the conservative form [14, 24]

$\displaystyle\partial_{t}{{\bm{w}}}+\nabla\cdot\vec{{\bm{f}}}({\bm{w}})-\nabla\cdot\vec{\bm{R}}({\bm{w}},\nabla{\bm{w}})={\bm{g}}\quad\mbox{in}\
Q_{T},$ (1)

where $\nabla$ and $\nabla\cdot$ denote the gradient and divergence operators,
respectively. The unknown state vector ${\bm{w}}$, the convective flux
$\vec{{\bm{f}}}=(\bm{f}_{1},\bm{f}_{2})$, the diffusive one
$\vec{\bm{R}}=(\bm{R}_{1},\bm{R}_{2})$ and the exterior forces ${\bm{g}}$ are
given by

$\displaystyle{\bm{w}}$ $\displaystyle=(\rho,\,\rho v_{1},\rho
v_{2},\,e)^{\mathrm{T}},\qquad\bm{f}_{i}=(\rho v_{i},\,\rho
v_{i}v_{1}+\delta_{i1}\mathrm{p},\,\rho
v_{i}v_{2}+\delta_{i2}\mathrm{p},\,(e+\mathrm{p})\,v_{i})^{\mathrm{T}},\qquad
i=1,2,$ (2) $\displaystyle\bm{R}_{i}$
$\displaystyle=\left(0,\,\tau_{i1},\,\tau_{i2},\sum\nolimits_{k=1}^{2}\tau_{ik}v_{i}+\frac{\mu\,c_{\mathrm{p}}}{\Pr}\partial_{x_{i}}{\theta}\right)^{\mathrm{T}},\
i=1,2,\qquad{\bm{g}}=\left(0,\,-\rho gk_{1},\,-\rho gk_{2},\,-\rho
g{\bm{k}}\cdot\bm{v}\right)^{\mathrm{T}},$

where $\delta_{ij}$ is the Kronecker delta and the symbol $\cdot$ denotes the
scalar product. The relations (1) – (2) are accompanied by the constitutive
relations of a perfect gas and the viscous part of the stress tensor as

$\displaystyle\mathrm{p}=R\rho\theta\qquad\mbox{and}\qquad\tau_{ij}=\mu\left(\partial_{x_{j}}{v_{i}}+\partial_{x_{i}}{v_{j}}\right)-\tfrac{2}{3}\delta_{ij}\nabla\cdot\bm{v},\quad
i,j=1,2,$ (3)

respectively. In atmospheric applications, the potential temperature $\Theta$
is the preferred physical quantity given as

$\displaystyle\Theta={\theta}/{\mathrm{P}},\qquad\mbox{where}\quad\mathrm{P}=\left({\mathrm{p}}/{\mathrm{p}_{0}}\right)^{(\kappa-1)/\kappa}$
(4)

is the Exner pressure and $\mathrm{p}_{0}=10^{5}\,\,\mathrm{Pa}$ is the
reference pressure.

The viscous fluxes can be written in the form

$\displaystyle\vec{\bm{R}}({\bm{w}},\nabla{\bm{w}})=\widehat{\bm{K}}({\bm{w}})\nabla{\bm{w}}\qquad\Longleftrightarrow\qquad\bm{R}_{i}({\bm{w}},\nabla{\bm{w}})=\sum\nolimits_{j=1}^{2}\bm{K}_{i,j}({\bm{w}})\partial_{x_{j}}{{\bm{w}}},\quad
i=1,2$ (5)

where $\widehat{\bm{K}}$ is the “matrix of matrices” (fourth-order tensor)
whose entries are $4\times 4$ matrices $\bm{K}_{i,j}(\cdot)$ dependent on
${\bm{w}}$ (and independent of $\nabla{\bm{w}}$). Their explicit form of
$\bm{K}_{i,j}(\cdot)$, $i,j=1,2$ can be found e.g., in [6] or [8, Section
9.1].

The term ${\bm{g}}$ representing the gravity forces can be expressed as

$\displaystyle{\bm{g}}=\bm{B}{\bm{w}},\qquad\bm{B}=\begin{pmatrix}0&0&0&0\\\
0&0&0&0\\\ -g&0&0&0\\\ 0&0&-g&0\\\ \end{pmatrix},$ (6)

where $g$ is the gravity constant (note that
${\bm{k}}=(k_{1},k_{2})^{\mathrm{T}}=(0,1)^{\mathrm{T}}$).

The above equations are accompanied by the initial and boundary conditions.
Particularly, we set ${\bm{w}}(x,0)={\bm{w}}_{0}(x)$ in ${\Omega}$, where
${\bm{w}}_{0}$ is the given function, typically as a perturbation in terms of
the potential temperature of a steady-state flow where the mean values
$\bar{\Theta}$ and $\bar{\mathrm{P}}$ of the Exner pressure and the potential
temperature, respectively, are in a hydrostatic balance, cf. [14], i.e.,

$\displaystyle
c_{\mathrm{p}}\bar{\Theta}\frac{\mathrm{d}\bar{\mathrm{P}}}{\mathrm{d}x_{2}}=-g.$
(7)

Moreover, we consider several types of boundary conditions on the boundary
${\Gamma}:=\partial{\Omega}$:

  * •

the no-flux or reflecting boundary condition $\bm{v}\cdot\vec{n}=0$, where
$\vec{n}$ is the outer normal to the boundary,

  * •

the nonreﬂecting boundary conditions representing an artificial boundary of an
infinite domain,

  * •

the periodic boundary conditions,

  * •

the symmetric boundary conditions.

Their practical realization in the context of the discontinuous Galerkin
method is described in Section 3.2.

##  3 Discontinuous Galerkin discretization

We solve numerically the system of equations (1) – (3) by the space-time
discontinuous Galerkin method. We refer to [8, Chapters 8–9] for the detailed
discretization of the Navier-Stokes equations without the external gravity
forces. Here, only the final necessary formulas are presented.

###  3.1 Functional spaces

Let $0=t_{0}<t_{1}<\ldots<t_{r}=T$ be a partition of $(0,T)$. We set
$I_{m}=(t_{m-1},t_{m})$ and $\tau_{m}=t_{m}-t_{m-1}$ for $m=1,\ldots,r$. For
each $t_{m},\ m=0,\ldots,r$, we consider a space partition
${{\mathscr{T}}_{h,m}}$ consisting of a finite number of closed triangles $K$
with mutually disjoint interiors and covering $\overline{{\Omega}}$, the
meshes ${{\mathscr{T}}_{h,m}}$ can vary for $m=0,\dots,r$, in general. By
${\vec{n}_{\\!\\!\;\scriptscriptstyle K}}$ we denote the unit outer normal to
the element boundary $\partial K$, $K\in{{\mathscr{T}}_{h,m}}$.

For each $K\in{{\mathscr{T}}_{h,m}}$, $m=0,\dots,r$, we assign a positive
integer number $p_{K}$ denoting the polynomial approximation degree on $K$.
Then we define the spaces of discontinuous piecewise-polynomial functions by

$\displaystyle S_{m,h,{\bm{p}}}=\\{\varphi:\Omega\to\mathbb{R};\
\varphi(x)|_{K}\in P_{p_{K}}(K)\
\forall\,K\in{{\mathscr{T}}_{h,m}}\\},\quad\bm{S}_{m,h,{\bm{p}}}=[S_{m,h,{\bm{p}}}]^{4},\qquad
m=0,\dots,r,$ (8)

where $P_{p_{K}}(K)$ denotes the space of all polynomials on $K$ of degree
$\leq p_{K}$.

The set of all edges of elements of mesh ${{\mathscr{T}}_{h,m}}$ is denoted by
${{\mathcal{F}}_{h,m}}$, the symbols ${{\mathcal{F}}_{h,m}^{I}}$ and
${{\mathcal{F}}_{h,m}^{B}}$ denote its subsets consisting of all interior and
boundary edges of ${{\mathscr{T}}_{h,m}}$, respectively. Moreover,
$\boldsymbol{\psi}|_{\partial K}^{\scriptscriptstyle(+)}$ and
$\boldsymbol{\psi}|_{\partial K}^{\scriptscriptstyle(-)}$ are the interior and
exterior traces of $\boldsymbol{\psi}\in\bm{S}_{m,h,{\bm{p}}}$ on $\partial
K,\ K\in{{\mathscr{T}}_{h,m}}$, respectively. The symbols
$\left\langle{\boldsymbol{\psi}}\right\rangle$ and
$[\\![{\boldsymbol{\psi}}]\\!]$ denote the mean value and the jump of
functions $\boldsymbol{\psi}\in\bm{S}_{m,h,{\bm{p}}}$ on edges
$\gamma\subset\partial K,\ K\in{{\mathscr{T}}_{h,m}}$ given by

$\displaystyle\begin{array}[]{llll}&\left\langle{\boldsymbol{\psi}}\right\rangle_{\gamma}=\tfrac{1}{2}\left(\boldsymbol{\psi}|_{\partial
K}^{\scriptscriptstyle(+)}+\boldsymbol{\psi}|_{\partial
K}^{\scriptscriptstyle(-)}\right),&[\\![{\boldsymbol{\psi}}]\\!]_{\gamma}=\left(\boldsymbol{\psi}|_{\partial
K}^{\scriptscriptstyle(+)}-\boldsymbol{\psi}|_{\partial
K}^{\scriptscriptstyle(-)}\right){\vec{n}_{\\!\\!\;\scriptscriptstyle
K}}&\mbox{for }\gamma\in{{\mathcal{F}}_{h,m}^{I}},\\\\[5.69054pt]
&\left\langle{\boldsymbol{\psi}}\right\rangle_{\gamma}=\boldsymbol{\psi}|_{\partial
K}^{\scriptscriptstyle(+)},&[\\![{\boldsymbol{\psi}}]\\!]_{\gamma}=\boldsymbol{\psi}|_{\partial
K}^{\scriptscriptstyle(+)}{\vec{n}_{\\!\\!\;\scriptscriptstyle K}}&\mbox{for
}\gamma\in{{\mathcal{F}}_{h,m}^{B}}.\\\ \end{array}$ (11)

Furthermore, let $q\geq 0$ be the polynomial approximation degree with respect
to time. We define the space of discontinuous piecewise polynomial functions
on the space-time domain $Q_{T}$ by

$\displaystyle\bm{S}_{h,{\bm{p}}}^{\tau,q}=$
$\displaystyle\left\\{\boldsymbol{\psi}:{\Omega}\times(0,T)\to\mathbb{R}^{4};\
\boldsymbol{\psi}(x,t)|_{{\Omega}\times
I_{m}}=\sum\nolimits_{l=0}^{q}t^{l}\,\boldsymbol{\varphi}_{k}(x),\ \
\boldsymbol{\varphi}_{k}\in\bm{S}_{m,h,{\bm{p}}},\ l=0,\dots,q,\
m=0,\dots,r\right\\}.$ (12)

Since $\boldsymbol{\psi}\in\bm{S}_{h,{\bm{p}}}^{\tau,q}$ is discontinuous with
respect to the time variable at $t_{m},\ m=1,\dots,r-1$, we define the jump
with respect to time on the time level $t_{m},\ m=0,\dots,r-1$ by

$\displaystyle\left\\{{\boldsymbol{\psi}}\right\\}_{m}:=\boldsymbol{\psi}|_{m}^{+}-\boldsymbol{\psi}|_{m}^{-},\qquad\boldsymbol{\psi}|_{m}^{\pm}:=\lim_{\delta\rightarrow
0\pm}\boldsymbol{\psi}(t_{m}+\delta),$ (13)

where $\boldsymbol{\psi}|_{0}^{-}$ is typically taken from the initial
condition.

###  3.2 Space semi-discretization

The space semi-discretization of the system (1) employs the following
properties and relations:

  * •

the convective fluxes fulfill
$\bm{f}_{i}({\bm{w}})=\bm{A}_{i}({\bm{w}}){\bm{w}}$, where $\bm{A}_{i}$ are
the Jacobian matrices of $\bm{f}_{i}({\bm{w}})$, $i=1,2$,

  * •

the physical inviscid flux through element boundary $\partial K$,
$K\in{{\mathscr{T}}_{h,m}}$ having the normal vector
${\vec{n}_{\\!\\!\;\scriptscriptstyle K}}=(n_{1},n_{2})$ is approximated by
the numerical flux

$\displaystyle\sum\nolimits_{i=1}^{2}\bm{f}_{i}({\bm{w}})n_{i}|_{\partial
K}\approx\bm{P}^{\scriptscriptstyle+}({\bm{w}},{\vec{n}_{\\!\\!\;\scriptscriptstyle
K}}){\bm{w}}|_{\gamma}^{\scriptscriptstyle(+)}+\bm{P}^{\scriptscriptstyle-}({\bm{w}},{\vec{n}_{\\!\\!\;\scriptscriptstyle
K}}){\bm{w}}|_{\gamma}^{\scriptscriptstyle(-)},$ (14)

where $\bm{P}^{\scriptscriptstyle\pm}$ are matrices from a suitable
decomposition of the matrix $\sum_{i=1}^{2}\bm{A}_{i}({\bm{w}})n_{i}$, e.g.,
Lax-Friedrichs or Vijayasundaram numerical fluxes can be written in this form
[12],

  * •

the diffusive terms in (1) are discretized by the symmetric interior penalty
variant of DG method.

By symbol $\left({\bm{u}},{\boldsymbol{\psi}}\right)_{M}$ we denote the
$L^{2}(M)$-scalar product of vector-valued functions
$\bm{u},\boldsymbol{\psi}:Q_{T}\to\mathbb{R}^{4}$ over one- or two-dimensional
set $M$. In order to simplify the notation, we introduce the symbol
$\vec{\bm{A}}=(\bm{A}_{1},\bm{A}_{2})$ denoting the “vector of matrices”
(third-order tensor). Then

$\displaystyle\left({\vec{\bm{A}}({\bm{w}}){\bm{w}}},{\nabla\boldsymbol{\psi}}\right)_{M}=\sum\nolimits_{i=1}^{2}\left({\bm{A}_{i}({\bm{w}}){\bm{w}}},{\partial_{x_{i}}{\boldsymbol{\psi}}}\right)_{M}.$
(15)

Similarly, in virtue of (5) and (11), we use the notation

$\displaystyle\left({\widehat{\bm{K}}({\bm{w}})\nabla{\bm{w}}},{\nabla\boldsymbol{\psi}}\right)_{M}=\sum\nolimits_{i,j=1}^{2}\left({\bm{K}_{ij}({\bm{w}})\partial_{x_{j}}{{\bm{w}}}},{\partial_{x_{i}}{\boldsymbol{\psi}}}\right)_{M}\
\ \mbox{and}\ \
\left({\left\langle{\widehat{\bm{K}}(\bar{{{\bm{w}}}})\nabla{\bm{w}}}\right\rangle},{[\\![{\boldsymbol{\psi}}]\\!]}\right)_{\gamma}=\sum\nolimits_{i,j=1}^{2}\left({\bm{K}_{ij}(\bar{{{\bm{w}}}})\partial_{x_{j}}{{\bm{w}}}},{(\boldsymbol{\psi}|_{\gamma}^{\scriptscriptstyle(+)}-\boldsymbol{\psi}|_{\gamma}^{\scriptscriptstyle(-)})}\right)_{\gamma}n_{i},$

where $(n_{1},n_{2})$ is the outer normal to an edge $\gamma\subset\partial
K$.

We define the forms representing the space semi-discretization of the
convective, diffusive and gravity terms in (1) by discontinuous Galerkin
method as

$\displaystyle{a_{h,m}^{\scriptscriptstyle\mathrm{L}}}(\bar{{{\bm{w}}}},{\bm{w}},\boldsymbol{\psi}):=-{\sum_{K\in{{\mathscr{T}}_{h,m}}}\left({\vec{\bm{A}}(\bar{{{\bm{w}}}}){\bm{w}}},{\nabla\boldsymbol{\psi}}\right)_{K}}+{\sum_{K\in{{\mathscr{T}}_{h,m}}}}\left(\left({\bm{P}^{\scriptscriptstyle+}(\bar{{{\bm{w}}}},{\vec{n}_{\\!\\!\;\scriptscriptstyle
K}}){\bm{w}}|_{\partial
K}^{\scriptscriptstyle(+)}},{\boldsymbol{\psi}}\right)_{\partial
K}+\left({\bm{P}^{\scriptscriptstyle-}(\bar{{{\bm{w}}}},{\vec{n}_{\\!\\!\;\scriptscriptstyle
K}}){\bm{w}}|_{\partial
K}^{\scriptscriptstyle(-)}},{\boldsymbol{\psi}}\right)_{\partial
K\setminus{\Gamma}}\right)-{\sum_{K\in{{\mathscr{T}}_{h,m}}}\left({\bm{B}{\bm{w}}},{\boldsymbol{\psi}}\right)_{K}}$
$\displaystyle\qquad+{\sum_{K\in{{\mathscr{T}}_{h,m}}}\left({\widehat{\bm{K}}(\bar{{{\bm{w}}}})\nabla{\bm{w}}},{\nabla\boldsymbol{\psi}}\right)_{K}}-{\sum_{\gamma\in{{\mathcal{F}}_{h,m}}}\left({\left\langle{\widehat{\bm{K}}(\bar{{{\bm{w}}}})\nabla{\bm{w}}}\right\rangle},{[\\![{\boldsymbol{\psi}}]\\!]}\right)_{\gamma}}-{\sum_{\gamma\in{{\mathcal{F}}_{h,m}}}\left({\left\langle{\widehat{\bm{K}}(\bar{{{\bm{w}}}})\nabla\boldsymbol{\psi}}\right\rangle},{[\\![{{\bm{w}}}]\\!]}\right)_{\gamma}}+{\sum_{\gamma\in{{\mathcal{F}}_{h,m}}}\left({\sigma[\\![{{\bm{w}}}]\\!]},{[\\![{\boldsymbol{\psi}}]\\!]}\right)_{\gamma}},$
(16)
$\displaystyle{\tilde{a}_{h,m}}(\bar{{{\bm{w}}}},\boldsymbol{\psi}):=-{\sum_{K\in{\mathscr{T}}_{h}}\left({\bm{P}^{\scriptscriptstyle-}(\bar{{{\bm{w}}}},{\vec{n}_{\\!\\!\;\scriptscriptstyle
K}}){\bm{w}}_{D}},{\boldsymbol{\psi}}\right)_{\partial
K\cap{\Gamma}}}-{\sum_{\gamma\in{{\mathcal{F}}_{h,m}^{B}}}\left({{\widehat{\bm{K}}(\bar{{{\bm{w}}}})\nabla\boldsymbol{\psi}}},{{{\bm{w}}_{D}{\vec{n}_{\\!\\!\;\scriptscriptstyle
K}}}}\right)_{\gamma}}+{\sum_{\gamma\in{{\mathcal{F}}_{h,m}^{B}}}\left({\sigma{\bm{w}}_{D}},{\boldsymbol{\psi}}\right)_{\gamma}},$
$\displaystyle{a_{h,m}}({\bm{w}},\boldsymbol{\psi}):={a_{h,m}^{\scriptscriptstyle\mathrm{L}}}({\bm{w}},{\bm{w}},\boldsymbol{\psi})-{\tilde{a}_{h,m}}({\bm{w}},\boldsymbol{\psi}),\qquad\bar{{{\bm{w}}}},{\bm{w}},\boldsymbol{\psi}\in\bm{S}_{m,h,{\bm{p}}},\quad
m=1,\dots,r,$

where $\sigma\sim p_{K}^{2}/\mbox{diam}(K)$ is the penalty parameter and
${\bm{w}}_{D}$ is the state vector following from the boundary conditions. We
note that the form ${a_{h,m}^{\scriptscriptstyle\mathrm{L}}}$ from (3.2) is
linear with respect to its second and third arguments, which is employed in
the iterative solution of the arising nonlinear algebraic systems in Section
4.

The boundary state vector ${\bm{w}}_{D}$ is set in the following way.

  * •

No-flux boundary conditions: the vector ${\bm{w}}_{D}$ has the same density,
energy and the tangential component of the velocity as ${\bm{w}}_{\partial K}$
(= the trace of ${\bm{w}}$ on $\partial K\subset{\Gamma}$) but the opposite
normal component of the velocity, i.e., we set

$\displaystyle{\bm{w}}_{D}|_{\partial K}={\mathscr{M}}({\bm{w}}|_{\partial
K},{\vec{n}_{\\!\\!\;\scriptscriptstyle
K}}),\quad\mbox{where}\quad{\mathscr{M}}({\bm{w}},\vec{n})=\left(\rho,\rho(\bm{v}-2(\bm{v}\cdot\vec{n})\vec{n}),e\right)^{\mathrm{T}}\quad\mbox{for}\quad{\bm{w}}=(\rho,\rho\bm{v},e)^{\mathrm{T}}.$
(17)

  * •

Periodic boundary conditions: we set
${\bm{w}}_{D}|_{\gamma^{+}}={\bm{w}}|_{\gamma^{-}}$ and
${\bm{w}}_{D}|_{\gamma^{-}}={\bm{w}}|_{\gamma^{+}}$, where $\gamma^{+}$ and
$\gamma^{-}$ denote the pair of the corresponding periodic edges.

  * •

Symmetric boundary conditions: ${\bm{w}}_{D}$ is set similarly as for the
periodic boundary conditions, depending on the line of symmetry. E.g., if the
symmetry is prescribed along the axes $x_{1}=0$, the vector ${\bm{w}}_{D}$ has
the same density, energy and the vertical component of the velocity $v_{2}$
but the opposite horizontal component of the velocity $v_{1}$.

  * •

Non-reflecting boundary conditions: we employ the approach from [8, Section
8.3.2], where ${\bm{w}}_{D}$ is the solution of the nonlinear Riemann problem
${\mathscr{B}}^{\rm RP}({\bm{w}}_{\Gamma},{\bm{w}}_{\mathrm{init}})$ defined
for the Euler equations ($\vec{\bm{R}}=0$ in (1)) without the source (gravity)
terms. Here, ${\bm{w}}_{\partial K}$ is the trace of ${\bm{w}}$ on $\partial
K\subset{\Gamma}$ and ${\bm{w}}_{\mathrm{init}}$ is typically given by the
initial (unperturbed) flow field.

We are aware that the aforementioned realization of the non-reflecting
boundary conditions is rather heuristic since the treated Riemann problem does
not involve the gravity term. However, this approach works very well provided
that high polynomial approximation degrees are employed, which is documented
by the numerical experiments in Section 3.3.

###  3.3 Performance of the non-reflecting boundary conditions

We consider the steady-state solution of hydrostatic flow in the computational
domain ${\Omega}=(-25\,600,25\,600)\times(0,6\,400)\,\mathrm{m}$ with the
constant potential temperature $\bar{\Theta}=300\,\mathrm{K}$ and the velocity
$\bm{v}=0$. The non-reflecting boundary conditions are prescribed on the
vertical parts of the boundary domain and the no-flux boundary conditions are
treated on the horizontal ones. We carried out computation till
$T=10\,\mathrm{s}$ and investigated the ability of the proposed numerical
scheme to keep the initial steady-state solution. We employed $P_{1}$, $P_{3}$
and $P_{5}$ polynomial approximation with respect to space using grids having
$\\#{\mathscr{T}}_{h}=117$, 406, 2 728 and 10 916 triangles. In Table 1, we
show the quantities

$\displaystyle\min{\delta\Theta}=\min_{x\in{\Omega}}\left(\Theta(x)-\bar{\Theta}\right)/\bar{\Theta},\qquad\max{\delta\Theta}=\max_{x\in{\Omega}}\left(\Theta(x)-\bar{\Theta}\right)/\bar{\Theta},\qquad\Delta{\delta\Theta}=\max{\delta\Theta}-\min{\delta\Theta}.$
(18)

We observe the large violation of the constant flow field for $P_{1}$
approximation even for relatively very fine grids. The decrease of
$\Delta{\delta\Theta}$ with respect to the mesh size step $h$ is slow. On the
other hand, $P_{5}$ approximation gives almost negligible values of
$\Delta{\delta\Theta}$ even for the coarsest grid. This effect underlines the
importance of the use of higher order approximation and it nicely correlates
with the used mesh adaptive algorithm, described in Section 5 which generates
typically meshes having large elements with high polynomial approximation
degrees for unperturbed (smooth) flows.

degree | $\\#{\mathscr{T}}_{h}$ | $\mathrm{DoF}$ | $\min{\delta\Theta}$ | $\max{\delta\Theta}$ | $\Delta({\delta\Theta})$
---|---|---|---|---|---
$P_{1}$ | 117 |  351 | -1.0672E-02 | 3.9217E-05 | 1.0711E-02
$P_{1}$ | 703 |  2 109 | -2.2365E-03 | 3.3900E-05 | 2.2704E-03
$P_{1}$ |  2 728 |  8 184 | -4.6223E-04 | 1.3397E-06 | 4.6357E-04
$P_{1}$ |  10 916 |  32 748 | -1.2871E-04 | 1.1181E-06 | 1.2983E-04
$P_{3}$ | 117 |  1 170 | -3.3737E-06 | 4.8027E-07 | 3.8539E-06
$P_{3}$ | 703 |  7 030 | -1.6309E-07 | 2.6487E-08 | 1.8957E-07
$P_{3}$ |  2 728 |  27 280 | -8.3347E-09 | 1.2776E-09 | 9.6122E-09
$P_{5}$ | 117 |  2 457 | -2.1223E-10 | 1.7873E-11 | 2.3010E-10
$P_{5}$ | 703 |  14 763 | -2.6707E-12 | 6.0177E-13 | 3.2725E-12
Table 1: Performance of the non-reflecting boundary conditions, the violence
of the constant potential temperature field, cf. (18). The value
$\mathrm{DoF}$ means the number of degrees of freedom of the space semi-
discretization per one equation.

###  3.4 Full space-time discretization

To define the full space-time DG discretization of (1), we introduce the forms

$\displaystyle{A_{h,m}}({\bm{w}},\boldsymbol{\psi})$
$\displaystyle:=\int_{I_{m}}\Bigl{(}\left({\partial_{t}{{\bm{w}}}},{\boldsymbol{\psi}}\right)_{\Omega}+{a_{h,m}}({\bm{w}},\psi)\Bigr{)}{\,\mathrm{d}t}+\left({\left\\{{{\bm{w}}}\right\\}_{m-1}},{\boldsymbol{\psi}|_{m-1}^{+}}\right)_{{\Omega}},$
(19)
$\displaystyle{A_{h,m}^{\scriptscriptstyle\mathrm{L}}}(\bar{{{\bm{w}}}},{\bm{w}},\boldsymbol{\psi})$
$\displaystyle:=\int_{I_{m}}\left(\left({\partial_{t}{{\bm{w}}}},{\boldsymbol{\psi}}\right)_{\Omega}+{a_{h,m}^{\scriptscriptstyle\mathrm{L}}}(\bar{{{\bm{w}}}},{\bm{w}},\psi)\right){\,\mathrm{d}t}+\left({\left\\{{{\bm{w}}}\right\\}_{m-1}},{\boldsymbol{\psi}|_{m-1}^{+}}\right)_{{\Omega}},\quad\bar{{{\bm{w}}}},{\bm{w}},\boldsymbol{\psi}\in\bm{S}_{h,{\bm{p}}}^{\tau,q},\quad
m=1,\dots,r.$

We note that form ${A_{h,m}^{\scriptscriptstyle\mathrm{L}}}$ is linear with
respect to its second and third arguments.

######  Definition 3.1.

We say that ${{\bm{w}}}_{h\tau}\in\bm{S}_{h,{\bm{p}}}^{\tau,q}$ is the space-
time discontinuous Galerkin solution of (1) if

$\displaystyle{A_{h,m}}({{\bm{w}}}_{h\tau},\boldsymbol{\psi}_{h})=0\qquad\forall\boldsymbol{\psi}_{h}\in\bm{S}_{h,{\bm{p}}}^{\tau,q},\quad
m=1,\dots,r,$ (20)

where we set ${{\bm{w}}}_{h\tau}|_{0}^{-}:={\bm{w}}_{0}$ (= the initial
condition).

The term
$\left({\left\\{{{{\bm{w}}}_{h\tau}}\right\\}_{m-1}},{\boldsymbol{\psi}|_{m-1}^{+}}\right)_{{\Omega}}$
in (19) represents the time penalty which joins together the approximate
solution on two subsequent time layers $I_{m-1}$ and $I_{m}$ in a weak sense.
This term is applicable also for varying spaces $\bm{S}_{m-1,h,{\bm{p}}}$ and
$\bm{S}_{m,h,{\bm{p}}}$. Indeed, due to (13), we have

$\displaystyle\left({\left\\{{{{\bm{w}}}_{h\tau}}\right\\}_{m-1}},{\boldsymbol{\psi}|_{m-1}^{+}}\right)_{{\Omega}}=\left({{{\bm{w}}}_{h\tau}|_{m-1}^{+}},{\boldsymbol{\psi}|_{m-1}^{+}}\right)_{{\Omega}}-\left({{{\bm{w}}}_{h\tau}|_{m-1}^{-}},{\boldsymbol{\psi}|_{m-1}^{+}}\right)_{{\Omega}}.$
(21)

If $\bm{S}_{m-1,h,{\bm{p}}}\not=\bm{S}_{m,h,{\bm{p}}}$ then the last term in
(21) exhibits the integration of the product of two functions which are
piecewise polynomial on different meshes ${{\mathscr{T}}_{h,m-1}}$ and
${{\mathscr{T}}_{h,m}}$. In Section 5, we introduce the $hp$-AMA algorithm
which generates completely different meshes for two sequential time steps, see
Figure 1 for an illustration. In this case, the evaluation of the last
integral in (21) is a little delicate since we integrate
${{\bm{w}}}_{h\tau}|_{m-1}^{-}$ over $K\in{{\mathscr{T}}_{h,m}}$ but
${{\bm{w}}}_{h\tau}|_{m-1}^{-}$ is discontinuous on $K$. In [11], we presented
a strategy, which employs a non-adaptive composite quadrature rule. The
evaluation of this integral exhibits another source of computational errors.
However, the numerical experiments in Section 5.1 show that the inaccuracy of
the mass and energy conservation is on the same level as the computation on a
fixed mesh. On the other hand, if
${{\mathscr{T}}_{h,m-1}}={{\mathscr{T}}_{h,m}}$, the numerical integration of
the terms in (21) can be carried out exactly (up to rounding errors) since the
integrated functions are polynomial.

Figure 1: Illustration of meshes ${{\mathscr{T}}_{h,m-1}}$ (red) and
${{\mathscr{T}}_{h,m}}$ (blue), total view (left) and a detail (right).

The STDG discretization (20) represents the system of ${N_{h,m}}$ algebraic
equations for each $m=1,\dots,r$, where

$\displaystyle{N_{h,m}}=(q+1){M_{h,m}}\quad\mbox{with}\quad{M_{h,m}}=\sum_{K\in{{\mathscr{T}}_{h,m}}}\nolimits(p_{K}+1)(p_{K}+2)/2.$
(22)

On the contrary, the diagonally implicit Runge-Kutta methods give $q$ systems
of size ${M_{h,m}}$ and multi-step methods give only one system of size
${M_{h,m}}$, hence their solution is cheaper. This is the main disadvantage of
the time discontinuous Galerkin discretization. We note that if $q=0$
(piecewise constant in time) then the time DG scheme is identical to the
backward (implicit) Euler method.

On the other hand, the STDG method has the super-convergence property in the
time nodes $t_{m},\ m=1,\dots,r$, namely $O(\tau^{2q+1})$. Therefore, the use
of piecewise linear approximation ($q=1$) seems to be sufficient for the
majority of numerical experiments since it provides a third-order method at
time nodes. Moreover, whereas multi-step of Runge-Kutta methods give the
approximate solution only at time nodes $t_{m},\ m=1,\dots,r$, the STDG
discretization provides an approximate solution on the whole interval $I_{m}$,
$m=1,\dots,r$ directly. It is advantageous in the definition of the space and
time discretization errors, cf. Section 4.2.

##  4 Solution strategy

In this section, we describe the strategy for solving the discrete problem
(20). We omit the details and we refer to [10], [8, Chapters 8-9] for a
complete description.

###  4.1 Solution of nonlinear algebraic systems

Let $m=1,\dots,r$ be an arbitrary but fixed index of the time level. The
nonlinear algebraic system (20) is solved iteratively by a Newton-like method
where we do not differentiate the form ${A_{h,m}}$ but employ its
linearization ${A_{h,m}^{\scriptscriptstyle\mathrm{L}}}$, cf. (19).
Particularly, let ${{\bm{w}}}_{h\tau}^{m}:={{\bm{w}}}_{h\tau}|_{{\Omega}\times
I_{m}}$ denote the restriction of ${{\bm{w}}}_{h\tau}$ on the space-time layer
${\Omega}\times I_{m}$. Then we define the sequence
${{\bm{w}}}_{h\tau}^{m,k}$, $k=0,1,\dots,$ approximating
${{\bm{w}}}_{h\tau}^{m}$ by

$\displaystyle{{\bm{w}}}_{h\tau}^{m,0}$
$\displaystyle={\mathscr{E}}({{\bm{w}}}_{h\tau}^{m-1}),$ (23a)
$\displaystyle{{\bm{w}}}_{h\tau}^{m,k}$
$\displaystyle={{\bm{w}}}_{h\tau}^{m,k-1}+\lambda_{k}{\bm{u}}_{h\tau}^{k},\qquad
k=1,2,\dots$ (23b) $\displaystyle\quad\mbox{where
}{\bm{u}}_{h\tau}^{k}\in\bm{S}_{h,{\bm{p}}}^{\tau,q}:\ \
{A_{h,m}^{\scriptscriptstyle\mathrm{L}}}({{\bm{w}}}_{h\tau}^{m,k-1},{\bm{u}}_{h\tau}^{k},\boldsymbol{\psi}_{h})=-{A_{h,m}}({{\bm{w}}}_{h\tau}^{m,k-1},\boldsymbol{\psi}_{h})\quad\forall\boldsymbol{\psi}_{h}\in\bm{S}_{h,{\bm{p}}}^{\tau,q}.$
(23c)

The symbol ${\mathscr{E}}$ in formula (23a) means that the initial guess
${{\bm{w}}}_{h\tau}^{m,0}$ is extrapolated either from the previous time layer
or from the initial condition. If
$\bm{S}_{m-1,h,{\bm{p}}}\not=\bm{S}_{m,h,{\bm{p}}}$, we have to project
${{\bm{w}}}_{h\tau}^{m-1}\in\bm{S}_{m-1,h,{\bm{p}}}$ to the new space
$\bm{S}_{m,h,{\bm{p}}}$ but we do not emphasize it explicitly. Moreover, the
function ${\bm{u}}_{h\tau}^{k}$ in (23b) is the update in the $k$-th iteration
and it is obtained from the linear algebraic system (23c). The convergence is
improved by the suitable choice of the damping factor $\lambda_{k}\in(0,1]$.
Its choice is controlled by the monitoring factor [4]

$\displaystyle\zeta_{k}=\frac{\left\|\left\\{{A_{h,m}}({{\bm{w}}}_{h\tau}^{m,k},\boldsymbol{\psi}_{i})\right\\}_{i=1}^{{N_{h,m}}}\right\|_{\ell^{2}}}{\left\|\left\\{{A_{h,m}}({{\bm{w}}}_{h\tau}^{m,k-1},\boldsymbol{\psi}_{i})\right\\}_{i=1}^{{N_{h,m}}}\right\|_{\ell^{2}}},$
(24)

where $\boldsymbol{\psi}_{i}$, $i=1,\dots,{N_{h,m}}$ are the basis functions
on ${\Omega}\times I_{m}$ such that
$\left\|\boldsymbol{\psi}_{i}\right\|_{L^{2}({\Omega}\times I_{m})}=1$ and
$\left\|\cdot\right\|_{\ell^{2}}$ is the Euclidean discrete norm. Hence, the
damping factor $\lambda_{k}$ has to be chosen such that $\zeta_{k}<1$ for each
iterates $k=1,2,\dots$.

It is advantageous to solve the linear algebraic system (23c) iteratively
since ${\bm{u}}_{h\tau}^{k}=0$ is a suitable initial guess and only a few
iterative steps are necessary. We employ the GMRES method with block ILU(0)
preconditioner but other solvers are possible. The stopping criterion for the
iterative process (23) is presented in Section 4.2. We note that it is not
necessary to update the matrix on the left-hand side of (23c) for each $k$. It
means that in the first argument of
${A_{h,m}^{\scriptscriptstyle\mathrm{L}}}$, we can employ a former
approximation ${{\bm{w}}}_{h\tau}^{m,k-s}$ for some $s>1$. Particularly, if
the monitoring function $\zeta_{k}$ is sufficiently smaller than 1, we do not
update the matrix on the left-hand side of (23c).

The presented approach can be treated also in the context of the IMEX
formulation. Let ${{\bm{w}}}_{h\tau}$ be the approximate solution given by
(20) and $\lambda\in(0,1]$ then we have the identity

$\displaystyle{A_{h,m}^{\scriptscriptstyle\mathrm{L}}}({{\bm{w}}}_{h\tau},{{\bm{w}}}_{h\tau},\boldsymbol{\psi}_{h})={A_{h,m}^{\scriptscriptstyle\mathrm{L}}}({{\bm{w}}}_{h\tau},{{\bm{w}}}_{h\tau},\boldsymbol{\psi}_{h})-\lambda{A_{h,m}}({{\bm{w}}}_{h\tau},\boldsymbol{\psi}_{h})\qquad\forall\boldsymbol{\psi}_{h}\in\bm{S}_{h,{\bm{p}}}^{\tau,q},\quad
m=1,\dots,r.$ (25)

Treating ${{\bm{w}}}_{h\tau}$ in the second argument of
${A_{h,m}^{\scriptscriptstyle\mathrm{L}}}$ on the left-hand side of (25)
implicitly and in the other arguments explicitly by extrapolating from the
previous time layer (cf. (23a)), we obtain the linear systems

$\displaystyle{A_{h,m}^{\scriptscriptstyle\mathrm{L}}}({\mathscr{E}}({{\bm{w}}}_{h\tau}^{m-1}),{{\bm{w}}}_{h\tau}^{m},\boldsymbol{\psi}_{h})={A_{h,m}^{\scriptscriptstyle\mathrm{L}}}({\mathscr{E}}({{\bm{w}}}_{h\tau}^{m-1}),{\mathscr{E}}({{\bm{w}}}_{h\tau}^{m-1}),\boldsymbol{\psi}_{h})-\lambda{A_{h,m}}({\mathscr{E}}({{\bm{w}}}_{h\tau}^{m-1}),\boldsymbol{\psi}_{h})\qquad\forall\boldsymbol{\psi}_{h}\in\bm{S}_{h,{\bm{p}}}^{\tau,q},\quad
m=1,\dots,r.$ (26)

The formulation (26) is equivalent to the first iterative step of method (23)
for $k=1$. So the additional iterates in (23) for $k\geq 2$ exhibits an
improvement of the accuracy which can be controlled by the stopping criterion
presented in Section 4.2. In our experience, the fully-implicit technique (23)
is more stable than the semi-implicit one (26), particularly it is possible to
use larger time steps. We demonstrate this effect by a numerical example in
Section 4.3.

###  4.2 Choice of the time step and iterative stopping criterion

The accuracy and efficiency of the whole computational process strongly depend
on the choice of the time step and the stopping criterion for the iterative
solver (23). A too large time step and/or a too weak stopping criterion cause
a loss of accuracy. On the other hand, a too small time step and/or a too
strong stopping criterion lead to an increase in the computational time
without a significant impact on the accuracy. In order to balance these two
aspects, we employ residual based estimators indicating the errors arising
from the spatial and temporal discretizations and the inaccuracy of iterative
algebraic solvers.

For each time interval $I_{m},\ m=1,\dots,r$, we define the estimators

$\displaystyle\eta^{\scriptscriptstyle\mathrm{A}}_{m}({{\bm{w}}}_{h\tau})=\sup_{0\not=\boldsymbol{\psi}_{h}\in\bm{S}_{h,{\bm{p}}}^{\tau,q}}\frac{{A_{h,m}}({{\bm{w}}}_{h\tau},\boldsymbol{\psi}_{h})}{\left\|\boldsymbol{\psi}_{h}\right\|_{X}},\
\
\eta^{\scriptscriptstyle\mathrm{S}}_{m}({{\bm{w}}}_{h\tau})=\sup_{0\not=\boldsymbol{\psi}_{h}\in\bm{S}_{h,{\bm{p}}+1}^{\tau,{q}}}\frac{{A_{h,m}}({{\bm{w}}}_{h\tau},\boldsymbol{\psi}_{h})}{\left\|\boldsymbol{\psi}_{h}\right\|_{X}},\
\
\eta^{\scriptscriptstyle\mathrm{T}}_{m}({{\bm{w}}}_{h\tau})=\sup_{0\not=\boldsymbol{\psi}_{h}\in\bm{S}_{h,{\bm{p}}}^{\tau,{q}+1}}\frac{{A_{h,m}}({{\bm{w}}}_{h\tau},\boldsymbol{\psi}_{h})}{\left\|\boldsymbol{\psi}_{h}\right\|_{X}},$
(27)

where the space of functions $\bm{S}_{h,{\bm{p}}}^{\tau,q}$ is given by (12)
and the spaces $\bm{S}_{h,{\bm{p}}+1}^{\tau,{q}}$ and
$\bm{S}_{h,{\bm{p}}}^{\tau,{q}+1}$ are given analogously by an increase of the
polynomial approximation degree by one with respect to space and time,
respectively. Particularly, we have

$\displaystyle\bm{S}_{h,{\bm{p}}+1}^{\tau,{q}}=$
$\displaystyle\left\\{\boldsymbol{\psi}:{\Omega}\times(0,T)\to\mathbb{R}^{4};\
\boldsymbol{\psi}(x,t)|_{{\Omega}\times
I_{m}}=\sum\nolimits_{l=0}^{q}t^{l}\,\boldsymbol{\varphi}_{k}(x),\ \
\boldsymbol{\varphi}_{k}\in\bm{S}_{m,h,{\bm{p}}+1},\ l=0,\dots,q,\
m=0,\dots,r\right\\},$ (28) $\displaystyle\bm{S}_{h,{\bm{p}}}^{\tau,{q}+1}=$
$\displaystyle\left\\{\boldsymbol{\psi}:{\Omega}\times(0,T)\to\mathbb{R}^{4};\
\boldsymbol{\psi}(x,t)|_{{\Omega}\times
I_{m}}=\sum\nolimits_{l=0}^{q+1}t^{l}\,\boldsymbol{\varphi}_{k}(x),\ \
\boldsymbol{\varphi}_{k}\in\bm{S}_{m,h,{\bm{p}}},\ l=0,\dots,q+1,\
m=0,\dots,r\right\\},$

where (cf. (8))

$\displaystyle\bm{S}_{m,h,{\bm{p}}+1}=[S_{m,h,{\bm{p}}+1}]^{4}\quad\mbox{with}\quad
S_{m,h,{\bm{p}}+1}=\\{\varphi:\Omega\to\mathbb{R};\ \varphi(x)|_{K}\in
P_{p_{K}+1}(K)\ \forall\,K\in{{\mathscr{T}}_{h,m}}\\},\qquad m=0,\dots,r.$
(29)

Moreover, the norm $\left\|\cdot\right\|_{X}$ is defined by

$\displaystyle\left\|\boldsymbol{\psi}_{h}\right\|_{X}=\left(\left\|\boldsymbol{\psi}_{h}\right\|_{L^{2}({\Omega}\times
I_{m})}^{2}+\left\|\nabla\boldsymbol{\psi}_{h}\right\|_{L^{2}({\Omega}\times
I_{m})}^{2}+\left\|\partial_{t}{\boldsymbol{\psi}}_{h}\right\|_{L^{2}({\Omega}\times
I_{m})}^{2}\right)^{1/2}.$ (30)

The quantity $\eta^{\scriptscriptstyle\mathrm{A}}_{m}({{\bm{w}}}_{h\tau})$
represents the inaccuracy of the solution of the algebraic system (20) since
it vanishes if (20) is solved exactly. In the same spirit,
$\eta^{\scriptscriptstyle\mathrm{S}}_{m}({{\bm{w}}}_{h\tau})$ and
$\eta^{\scriptscriptstyle\mathrm{T}}_{m}({{\bm{w}}}_{h\tau})$ represent the
inaccuracy of the spatial and temporal discretization, respectively.
Particularly, they vanish if the approximate solution is equal to the exact
one. The numerical study of these error estimators was presented in [10]. We
note that due to the choice of the norm (30), the estimators in (27) can be
evaluated as the sum of local estimators computed for each space-time element
$K\times I_{m}$, $K\in{{\mathscr{T}}_{h,m}}$ separately. These local
estimators are computable by standard tools of the constrained minimization
since the corresponding systems of equations are small and thus the evaluation
of (27) is fast.

Employing estimators (27), we define the stopping criterion of the iterative
method (23) by

$\displaystyle\eta^{\scriptscriptstyle\mathrm{A}}_{m}({{\bm{w}}}_{h\tau})\leq
c_{A}\min\left(\eta^{\scriptscriptstyle\mathrm{S}}_{m}({{\bm{w}}}_{h\tau}),\eta^{\scriptscriptstyle\mathrm{T}}_{m}({{\bm{w}}}_{h\tau})\right),\quad
m=1,\dots,r,$ (31)

where $c_{A}\in(0,1)$. It means that the algebraic error estimator has to be
bounded with respect to the spatial and temporal ones. Similarly, we choose
adaptively the time step $\tau_{m}$, $m=1,\dots,r$ such that spatial and
temporal error estimators are balanced, particularly, we set the time step
such that

$\displaystyle\eta^{\scriptscriptstyle\mathrm{T}}_{m}({{\bm{w}}}_{h\tau})\approx
c_{T}\eta^{\scriptscriptstyle\mathrm{S}}_{m}({{\bm{w}}}_{h\tau}),\quad
m=1,\dots,r,$ (32)

for some $c_{T}\in(0,1)$. In practice, we employ the values $c_{A}=0.01$ and
$c_{T}=0.2$.

###  4.3 Numerical study of the stability and adaptivity of the STDG scheme

To demonstrate the stability and the adaptive performance of the fully
implicit scheme, we consider the rising thermal bubble from Section 6.2 with
the smooth initial condition (41). We consider a fixed quasi-uniform grid
having 556 elements, $P_{2}$ approximation with respect to space and $P_{1}$
approximation with respect to time variables. Therefore, the size of the
algebraic system is ${N_{h,m}}=\,\numprint{26688}$. We compare the
computations obtained by

  * •

fully implicit scheme (23) using the adaptive choice of the time step
according to (32) with $c_{T}=0.2$ and adaptive stopping criterion (31) with
$c_{A}=0.01$,

  * •

semi-implicit scheme (26) using the adaptive choice of the time step according
to (32) with $c_{T}=0.1$ which guarantees the stability and a comparable
quality of the results as the fully implicit scheme – the larger value of
$c_{T}$ gives an unstable solution.

Figure 2 shows the size of the time step $\tau_{m}$ and the accumulated number
of GMRES iterations for each time step $m=1,2,\dots$. We observe that the
fully implicit scheme admits the use of a larger size of the time steps and
the corresponding number of GMRES iterations is smaller. We note that the size
of $\tau_{m}$ in the last time steps corresponds to $\mbox{CFL}=400$ in the
“common” stability condition [8, Section 8.4]

$\displaystyle\tau_{m}\leq\min_{\gamma\subset\partial
K}\frac{\mbox{CFL}\,|K|}{|\gamma|\,\Lambda({{\bm{w}}}_{h\tau},{\vec{n}_{\\!\\!\;\scriptscriptstyle
K}})}\qquad\forall\gamma\subset\partial K\quad\forall
K\in{{\mathscr{T}}_{h,m}},$ (33)

where $|K|$ denotes the area of $K\in{{\mathscr{T}}_{h,m}}$, $|\gamma|$ is the
length of an edge $\gamma\subset\partial K$ and
$\Lambda({{\bm{w}}}_{h\tau},{\vec{n}_{\\!\\!\;\scriptscriptstyle K}})$ is the
maximal eigenvalue of matrix
$\bm{P}=\bm{P}^{\scriptscriptstyle(+)}+\bm{P}^{\scriptscriptstyle(-)}$ arising
in the numerical flux (14). The computational time for the semi-implicit
method was $\,\numprint{2437}\,\mathrm{s}$ whereas for the implicit one
$\,\numprint{1363}\,\mathrm{s}$.

Figure 2: Comparison of the semi-implicit (blue) and fully implicit (red)
schemes: the dependence of the size of the time step $\tau_{m}$ (left) and the
accumulated number of GMRES iterations (right) with respect to the physical
time $t\in(0,700)\,\mathrm{s}$. Each dot corresponds to one time step
$m=1,\dots,r$.

Moreover, Figure 3 shows the performance of the adaptive algebraic stopping
criterion (31) and the choice of the time step (32). We plot the size of
estimators $\eta^{\scriptscriptstyle\mathrm{A}}_{m}$,
$\eta^{\scriptscriptstyle\mathrm{S}}_{m}$,
$\eta^{\scriptscriptstyle\mathrm{S}}_{m}$ (cf. (27)) and the size of the time
step $\tau_{m}$ for several starting and final time levels $m$. We observe
that the size of the time steps and algebraic errors are controlled by the
aforementioned conditions. Moreover, only 2-3 nonlinear iterations are
performed for smaller $m$ when the time steps are shorter and about $3-5$
steps at the end of the computation where the time steps are larger.

Figure 3: Performance of the adaptive algebraic stopping criteria and the
adaptive choice of the time step: estimators
$\eta^{\scriptscriptstyle\mathrm{A}}_{m}$,
$\eta^{\scriptscriptstyle\mathrm{S}}_{m}$,
$\eta^{\scriptscriptstyle\mathrm{T}}_{m}$ (cf. (27)) and $\tau_{m}$ for
$m=1,\dots,15$ (left) and $m=141,\dots,156$ (right), each dot corresponds to
one iterative step $k$ in (23).

##  5 Mesh adaptation

Mesh adaptive methods try to reduce the number of degrees of freedom
($\mathrm{DoF}$) while the accuracy is preserved. This is beneficial namely
for the simulation of evolution problems with traveling phenomena having sharp
interfaces. The idea of the presented adaptation procedure is to generate the
sequence of the spaces $S_{m,h,{\bm{p}}}$, $m=0,\dots,r$ (i.e., the meshes
${{\mathscr{T}}_{h,m}}$ and the corresponding polynomial degrees $p_{K}$,
$K\in{{\mathscr{T}}_{h,m}}$) such that

  1. (i)

the interpolation error of a selected scalar quantity is under the given
tolerance $\omega>0$,

  2. (ii)

the number of degrees of freedom ($={N_{h,m}}$, cf. (22)) is minimal for each
$m=0,\dots,r$.

More precisely, let
${q}_{h\tau}={q}_{h\tau}({{\bm{w}}}_{h\tau}):{\Omega}\times(0,T)\to\mathbb{R}$
denote a scalar solution-dependent quantity, e.g., the density and the Mach
number are commonly used for aerodynamic problems, cf. [2, 1, 13, 19]. Then
condition (i) means

$\displaystyle\left\|\widehat{q}_{h\tau}(t)-\Pi\widehat{q}_{h\tau}(t)\right\|_{L^{2}({\Omega})}\leq\omega\quad\forall
t\in I_{m},\ m=1,\dots,r,$ (34)

where $\widehat{q}_{h\tau}(t)\in S_{m,h,{\bm{p}}+1}$, cf. (29) is the higher-
order reconstruction of the scalar piecewise polynomial function
${q}_{h\tau}(t)={q}_{h\tau}({{\bm{w}}}_{h\tau}(\cdot,t))$ for $t\in I_{m}$ and
$\Pi$ is the $L^{2}({\Omega})$-projection to $S_{m,h,{\bm{p}}}$. In practice,
we relax condition (34) and require

$\displaystyle\left\|\widehat{q}_{h\tau}(t)-\Pi\widehat{q}_{h\tau}(t)\right\|_{L^{2}({\Omega})}\leq\omega\quad\forall
t\in\mathcal{J}_{m},\ m=1,\dots,r,$ (35)

where $\mathcal{J}_{m}$ is the finite set of selected $t\in I_{m}$. Typically
$\mathcal{J}_{m}$ contains endpoints of $I_{m}$ and the time integration
nodes. The reconstruction $\widehat{q}_{h\tau}$ is carried out locally by a
least-square technique, cf. [9, Section 9.1.1].

The whole adaptive procedure is described by Algorithm 1. We start with an
initial mesh, perform one time step by solving problem (20), evaluate
${q}_{h\tau}={q}_{h\tau}({{\bm{w}}}_{h\tau})$ and reconstruct
$\widehat{q}_{h\tau}(t)$, $t\in\mathcal{J}_{m}$. If condition (35) is valid we
proceed to the next time step using the same mesh and polynomial approximation
degrees. Otherwise, we carry out the re-meshing using the anisotropic
$hp$-mesh adaptation ($hp$-AMA) technique (step 18 of Algorithm 1), whose
theoretical as well as practical issues are described in detail in the recent
monographs [9], see also [11].

1: inputs: tolerance $\omega>0$, initial time step $\tau_{1}$, initial
$hp$-mesh ${{\mathscr{T}}_{h,0}}$ and $p_{K}\geq 1$ $\forall
K\in{{\mathscr{T}}_{h,0}}$

2: set $m=1$, $t=0$

3: while $t<T$ do

4: repeat

5: repeat

6: solve ${A_{h,m}}({{\bm{w}}}_{h\tau},\boldsymbol{\psi}_{h})=0$
$\forall\boldsymbol{\psi}_{h}\in\bm{S}_{h,{\bm{p}}}^{\tau,q}$, cf. (20) until
(31) is fulfilled

7: evaluate $\eta^{\scriptscriptstyle\mathrm{S}}_{m}({{\bm{w}}}_{h\tau})$ and
$\eta^{\scriptscriptstyle\mathrm{T}}_{m}({{\bm{w}}}_{h\tau})$ by (27)

8: if
$\eta^{\scriptscriptstyle\mathrm{T}}_{m}>c_{T}\eta^{\scriptscriptstyle\mathrm{S}}_{m}$
then

9: reduce the size of the time step $\tau_{m}$

10: end if

11: until $\eta^{\scriptscriptstyle\mathrm{T}}_{m}\leq
c_{T}\eta^{\scriptscriptstyle\mathrm{S}}_{m}$

12: evaluate $\widehat{q}_{h\tau}$ by a higher-order reconstruction of the
quantity ${q}_{h\tau}({{\bm{w}}}_{h\tau})$

13: set
$\eta_{m}^{\mathrm{int}}:=\max_{t\in\mathcal{J}_{m}}\left\|\widehat{q}_{h\tau}(t)-\Pi\widehat{q}_{h\tau}(t)\right\|_{L^{2}({\Omega})}$

14: if $\eta_{m}^{\mathrm{int}}\leq\omega$ then

15: set $t:=t+\tau_{m}$, propose the new size of the time step $\tau_{m+1}$

16: set ${{\mathscr{T}}_{h,m+1}}:={{\mathscr{T}}_{h,m}}$ and the same $p_{K},\
K\in{{\mathscr{T}}_{h,m}}$, set $m:=m+1$

17: else

18: generate a new mesh ${{\mathscr{T}}_{h,m}}$ & new polynomial degrees
$p_{K},\ K\in{{\mathscr{T}}_{h,m}}$ by $hp$-AMA technique

19: reduce the size of the time step $\tau_{m}$

20: end if

21: until $\eta_{m}^{\mathrm{int}}\leq\omega$

22: end while

Algorithm 1 Space-time anisotropic $hp$-mesh adaptive algorithm.

The $hp$-AMA method sets an optimal size of mesh elements by the
equidistribution of an interpolation error estimate, i.e., the aim is to
fulfill the condition

$\displaystyle\max_{t\in\mathcal{J}_{m}}\left\|\widehat{q}_{h\tau}(t)-\Pi\widehat{q}_{h\tau}(t)\right\|_{L^{2}(K)}\approx\frac{\omega}{\sqrt{\\#{\mathscr{T}}_{h}}}\quad\forall
K\in{{\mathscr{T}}_{h,m}},\ m=1,\dots,r,$ (36)

where $\\#{\mathscr{T}}_{h}$ is the number of elements of
${{\mathscr{T}}_{h,m}}$. Moreover, the shape and orientation of mesh elements
are optimized locally in such a way that for each $K\in{{\mathscr{T}}_{h,m}}$,
we fix its size and seek a new shape and orientation by minimizing the
interpolation error estimate. Similarly, we can modify the polynomial
approximation degree such that we fix the density of the number of degrees of
freedom and minimize the interpolation error estimate. This estimate depends
on the ($p_{K}+1$)-directional derivatives of the reconstructed solution
$\widehat{q}_{h\tau}$.

There is a natural question of which quantity ${q}_{h\tau}$ should be used for
the higher-order reconstruction. Solving non-hydrostatic mesoscale atmospheric
problems, the very frequent quantity of interest is the perturbation of the
potential temperature $\Theta$, cf. (4) with respect to the initial state.
Hence, it seems to be favorable to employ $\Theta$ in the $hp$-AMA algorithm.
Nevertheless, numerical experiments in Section 5.1 demonstrate that the simple
choice ${q}_{h\tau}$ as the density is also possible.

######  Remark 5.1.

We are aware that the computational time necessary for the solution of
algebraic systems (the most time-consuming part of the whole computational
process) does not depend only on the size of these systems. It is influenced
also by the shape of elements, mesh grading, polynomial approximation degrees
and the size of time steps. Nevertheless, the reduction of the number of
$\mathrm{DoF}$ typically accelerates numerical simulations.

###  5.1 Numerical study of the mass and energy conservation and the choice of
the quantity of interest

The goal of the following numerical experiments is to demonstrate that (i) the
use of different meshes at different time layers does not have a significant
impact on the accuracy of the computations and (ii) the density is suitable as
the quantity of interest in the $hp$-AMA algorithm. We consider the rising
thermal bubble test case in Section 6.2 with the initial condition (41). We
compare results achieved on

  * •

quasi-uniform grid having 2 248 elements using $P_{3}$ polynomial
approximation with respect to space and $P_{1}$ polynomial approximation with
respect to time variables,

  * •

anisotropic $hp$-meshes produced by Algorithm 1 using density as the quantity
${q}_{h\tau}$ with three tolerances $\omega$.

Figure 4 shows the final mesh (without the polynomial degree distribution),
the isolines of the potential temperature perturbation and the isolines of the
density. We observe that the potential temperature perturbation isolines give
the expected mushroom cloud with increasing accuracy for the increasing mesh
refinement. On the other hand, the density distribution is constant in the
horizontal direction and linear in the vertical one with a small perturbation
around the mushroom cloud. Even these small perturbations are sufficient for
the $hp$-AMA algorithm to detect the regions of the largest interpolation
error and using the evaluated higher order directional derivatives, the
algorithm sets the required shape and orientation of the mesh elements.

mesh

potential temperature

density

uniform |  $hp$-AMA, $\omega=0.02$ |  $hp$-AMA, $\omega=0.015$ |  $hp$-AMA, $\omega=0.01$
---|---|---|---
$\Delta_{\mathrm{mass}}$ = 9.761E-06 |  $\Delta_{\mathrm{mass}}$ = 9.425E-06 |  $\Delta_{\mathrm{mass}}$ = 9.537E-06 |  $\Delta_{\mathrm{mass}}$ = 9.649E-06
$\Delta_{\mathrm{energy}}$ = 1.207E-05 |  $\Delta_{\mathrm{energy}}$ = 1.173E-05 |  $\Delta_{\mathrm{energy}}$ = 1.179E-05 |  $\Delta_{\mathrm{energy}}$ = 1.193E-05
$\\#{\mathscr{T}}_{h}=2\,248$, $\mathrm{DoF}$ = 107 904 |  $\\#{\mathscr{T}}_{h}=247$, $\mathrm{DoF}$ =11 856 |  $\\#{\mathscr{T}}_{h}=247$, $\mathrm{DoF}$ =11 856 |  $\\#{\mathscr{T}}_{h}=319$, $\mathrm{DoF}$ =15 312

Figure 4: Smooth rising thermal bubble: meshes (first line), isolines of the
potential temperature perturbation with the contour interval of $0.05$ (second
line) and the isolines of the density with the contour interval of $0.05$
(third line); the computations using the uniform grid (left column) and
anisotropic meshes generated with the tolerances $\omega=0.02$, $0.015$ and
$0.01$ (from second to fourth columns) at $T=500\,\mathrm{s}$. Note: The
numbers of $\\#{\mathscr{T}}_{h}$ and $\mathrm{DoF}$ are the same for
$\omega=0.02$ and $\omega=0.015$ which is just a coincidence.

Moreover, we evaluate the quantities

$\displaystyle\Delta_{\mathrm{mass}}=\left|\frac{M(T)-M(0)}{M(0)}\right|\ \
\mbox{with}\ \
M(t)=\int_{{\Omega}}\rho_{h\tau}(x,t){\,\mathrm{d}x}\quad\mbox{and}\quad\Delta_{\mathrm{energy}}=\left|\frac{E(T)-E(0)}{E(0)}\right|\
\ \mbox{with}\ \ E(t)=\int_{{\Omega}}e_{h\tau}(x,t){\,\mathrm{d}x},$ (37)

where $\rho_{h\tau}$ and $e_{h\tau}$ are the approximations of the density and
the energy, respectively, given by the STDG method. The corresponding values
are given at the bottom of Figure 4 together with the number of mesh elements
${{\mathscr{T}}_{h,m}}$ and $\mathrm{DoF}={N_{h,m}}=\dim\bm{S}_{m,h,{\bm{p}}}$
for the last time layer $m=r$. We observe the same level of violation of the
mass and energy conservation for the computations using fixed and varying
meshes which indicates that the resulting numerical solutions do not suffer
from the errors arising from multiple recomputations between adapted meshes.

##  6 Numerical experiments

In this section, we present the numerical solution of several benchmarks from
[14] by the combination of the STDG and $hp$-AMA methods, cf. Algorithm 1. The
computations have been carried out by the in-house code ADGFEM [7] employing
ANGENER code [5] for the mesh adaptation. For some of the presented examples,
the $hp$-AMA method significantly reduces the computational costs in
comparison to non-refined grids. On the other hand, some test problems do not
have anisotropic features so the $hp$-AMA method generates grids that are
close to the uniform ones. Similarly as in [14], all examples except the sharp
thermal bubble in Section 6.2 and the density current flow in Section 6.3 are
treated as inviscid flows, i.e., $\vec{\bm{R}}=0$ in (1). The physical
constants appearing in (2) – (4) take the values $\kappa=1.4$, $\Pr=0.72$ and
$c_{\mathrm{v}}=718\,\mathrm{J\cdot kg^{-1}\cdot K^{-1}}$. We use the
piecewise linear approximation with respect to time ($q=1$ in (12)) for all
numerical examples except the steady-state problem in Section 6.4 when we set
$q=0$.

###  6.1 Inertia gravity wave

The inertia gravity wave benchmark (cf. [14, 29]) exhibits the evolution of a
potential temperature perturbation in the rectangular domain
${\Omega}=(0,\,\numprint{300000})\times(0,\,\numprint{10000})\,\mathrm{m}$
with the periodic boundary conditions on the left and right boundaries of
${\Omega}$ and the no-flux boundary conditions on the horizontal boundary
parts. The initial state of the atmosphere has the mean horizontal velocity
$v_{1}=20\,\mathrm{m\cdot s^{-1}}$ in a uniformly stratified atmosphere with
the Brunt-Väisälä frequency $N=0.01\,\mathrm{s^{-1}}$ where the mean potential
temperature fulfills

$\displaystyle\bar{\Theta}=\Theta_{0}\exp(N^{2}x_{2}/g),\qquad\Theta_{0}=300\,\mathrm{K}.$
(38)

The Exner pressure is then given by the hydrostatic balance equation (7),
namely

$\displaystyle\bar{\mathrm{P}}=1+\frac{g^{2}}{c_{\mathrm{p}}\,\Theta_{0}\,N^{2}}\left(\exp(-N^{2}x_{2}/g)-1\right).$
(39)

Moreover, to the mean potential temperature $\bar{\Theta}$, we add the
perturbation

$\displaystyle{\Theta}^{\prime}=\Theta_{c}\sin(\pi
x_{2}/h_{c})\left(1+((x_{1}-\bar{x}_{1})/a_{c})^{2}\right)^{-1},$ (40)

where $\Theta_{c}=0.01\,\mathrm{K}$, $h_{c}=\,\numprint{10000}\,\mathrm{m}$,
$a_{c}=\,\numprint{5000}\,\mathrm{m}$ and
$\bar{x}_{1}=\,\numprint{100000}\,\mathrm{m}$. The initial perturbation
propagates to the left and right symmetrically, but it is carried away by the
mean horizontal flow till the final time $T=3\,000\,\mathrm{s}$.

The solution of this problem is smooth, so we expect that the mesh adaptive
algorithm generates a sequence of coarse meshes with a low number of elements
and with high polynomial approximation degrees. Figure 5 presents the
potential temperature perturbation at the physical times
$t=\,\numprint{0}\,\mathrm{s}$, $t=\,\numprint{2000}\,\mathrm{s}$ and
$t=\,\numprint{3000}\,\mathrm{s}$, namely the isolines and the horizontal
distribution along $x_{2}=\,\numprint{5000}\,\mathrm{m}$. The left figures
indicate also the generated triangular meshes, the polynomial approximation
degrees $p_{K}$ for (almost) all $K\in{{\mathscr{T}}_{h,m}}$ are equal to 8
which is the maximal degree supported by our code.

Figure 5: Interior gravity wave: isolines of the potential temperature
perturbation (the contour interval of $5\cdot 10^{-4}\,\mathrm{K}$) with the
corresponding mesh (left) and the potential temperature perturbation along
$x_{2}=5\,000\,\mathrm{m}$ at $t=0\,\mathrm{s}$ (top),
$t=\,\numprint{2000}\,\mathrm{s}$ (center) and
$t=\,\numprint{3000}\,\mathrm{s}$ (bottom).

We observe a smooth resolution of the isolines although relatively coarse
meshes have been generated. The profiles are perfectly symmetric along the
vertical position $x_{1}=\bar{x}_{1}+v_{1}t$, i.e.,
$x_{1}=\,\numprint{140000}\,\mathrm{m}$ for $t=\,\numprint{2000}\,\mathrm{s}$
and $x_{1}=\,\numprint{160000}\,\mathrm{m}$ for
$t=\,\numprint{3000}\,\mathrm{s}$. The results are in full agreement with the
reference data in [14, Section 6.1].

###  6.2 Rising thermal bubbles

The rising thermal bubble test case simulates the evolution of a warm bubble
in a constant potential temperature field. Since the bubble is warmer, it
rises while deforming as a consequence of the shearing motion caused by the
velocity field gradients until it forms a mushroom cloud. According to [14,
27], we consider the computational domain
${\Omega}=(0,1000)\times(0,1000)\,\mathrm{m}$, the constant mean potential
temperature $\bar{\Theta}=300\,\mathrm{K}$ and the Exner pressure follows from
the hydrostatic balance (7) as
$\mathrm{P}=1-g/(c_{\mathrm{p}}\bar{\Theta})x_{2}$. The initial velocity is
equal to zero. Moreover, to the mean flow, we add two different potential
temperature perturbations, the smooth one

$\displaystyle\Theta^{\prime}=\begin{cases}0&\mbox{for }r>r_{c},\\\
\tfrac{1}{2}\Theta_{c}\left(1+\cos(\pi\,r/r_{c})\right)&\mbox{for }r\leq
r_{c},\\\ \end{cases}$ (41)

where $\Theta_{c}=0.5\,\mathrm{K}$, $r=|x-x_{c}|$,
$x_{c}=(500,300)\,\mathrm{m}$ and $r_{c}=250\,\mathrm{m}$, and the sharp one

$\displaystyle\Theta^{\prime}=\begin{cases}\Theta_{c}&\mbox{for }r<r_{c},\\\
\Theta_{c}\exp(-(r-r_{c})^{2}/s^{2})&\mbox{for }r\geq r_{c},\\\ \end{cases}$
(42)

where $\Theta_{c}=0.5\,\mathrm{K}$, $r=|x-x_{c}|$,
$x_{c}=(500,260)\,\mathrm{m}$, $r_{c}=150\,\mathrm{m}$ and $s=10$. We note
that in the original paper [27], the parameter $s$ in (42) was chosen as
$s=100$, so we have much sharper initial perturbation in our case. For the
smooth initial bubble, the final physical time is set to $T=700\,\mathrm{s}$.
For the sharper one, we consider $T=850\,\mathrm{s}$ and a small amount of the
dynamic viscosity $\mu=3.44\cdot 10^{-3}\,\mathrm{Pa\cdot s}$. The no-flux
boundary conditions are prescribed on the whole boundary.

Figures 6 and 7 show the distribution of the potential temperature
perturbation and the corresponding $hp$-meshes at selected time steps for both
initial settings. We observe an accurate capture of the deformed bubbles
together with the expected anisotropic mesh refinement. We note that all
figures are visualized by the library Matplotlib [20] where no additional
smoothing of the results is done.

Figure 6: Smooth rising thermal bubble: the distribution of the potential
temperature perturbation (top) and the corresponding $hp$-meshes (bottom) at
$t=0\,\mathrm{s}$, $t=350\,\mathrm{s}$ and $t=700\,\mathrm{s}$ (from left to
right).

Figure 7: Sharp rising thermal bubble: the distribution of the potential
temperature perturbation (top) and the corresponding $hp$-meshes (bottom) at
$t=0\,\mathrm{s}$, $t=425\,\mathrm{s}$ and $t=850\,\mathrm{s}$ (from left to
right).

###  6.3 Density current flow problem

The density current flow problem benchmark was proposed in [31] and it
exhibits the evolution of a cold bubble dropped in a neutrally stratified
atmosphere. Since the bubble is cold, it sinks to the ground and starts to
move along the ground forming Kelvin–Helmholtz vortexes.

In the same way as [14], we consider the computational domain
${\Omega}=(0,\,\numprint{25600})\times(0,\,\numprint{6400})\,\mathrm{m}$, the
final time $T=900\,\mathrm{s}$. We prescribe the symmetric boundary condition
along $x_{1}=0$ and the no-flux boundary conditions over the rest of the
boundary. It was discussed in [31] that a viscosity is required to obtain a
grid-converged solution, hence we set $\mu=0.1\,\mathrm{Pa\cdot s}$ in (3).
The initial flow is given by the velocity $\bm{v}=0$ and the constant mean
potential temperature $\bar{\Theta}=300\,\mathrm{K}$ perturbed by

$\displaystyle\Theta^{\prime}=\frac{\Theta_{C}}{2}\left(1+\cos(\pi
r)\right),\qquad\Theta_{C}=-15\,\mathrm{K},\quad
r=\left(\left({(x_{1}-x_{1}^{c})}/{x_{1}^{r}}\right)^{2}+\left({(x_{2}-x_{2}^{c})}/{x_{2}^{r}}\right)^{2}\right)^{1/2},$
(43)

where the center of the bubble is at
$(x_{1}^{c},x_{x}^{c})=(0,\,\numprint{3000})\,\mathrm{m}$ and its size is
$(x_{1}^{r},x_{2}^{r})=\left(\,\numprint{4000},\,\numprint{2000}\right)\,\mathrm{m}$.
The Exner pressure follows from the balance equation (7).

Figure 8 shows the distribution of the perturbation of the potential
temperature and the corresponding $hp$-meshes at selected time levels. The
development of Kelvin–Helmholtz vortexes is easily observed together with the
expected mesh adaptation. Due to high polynomial approximation degrees, the
vortexes are captured sharply even for relatively coarse meshes, the maximal
resolution is approximately $250\,\mathrm{m}$. We note that in the reference
paper [31], the larger viscosity $\mu=75\,\mathrm{Pa\cdot s}$ has been used.

Figure 8: Density current flow: the distribution of the potential temperature
perturbation (left) and the corresponding $hp$-meshes (right) at
$t=300\,\mathrm{s}$, $t=500\,\mathrm{s}$ , $t=700\,\mathrm{s}$ and
$t=\,\numprint{900}\,\mathrm{s}$ (from top to bottom).

###  6.4 Schär mountain

This test case deals with the steady-state solution of hydrostatic flow over a
five-peak mountain chain, with steady inflow and outflow boundary conditions.
In the same way as in [14, 28], the initial state of the atmosphere has a
constant mean flow $v_{1}=10\,\mathrm{m\cdot s^{-1}}$ in a uniformly
stratified atmosphere with a Brunt-Väisälä frequency of
$N=0.01\,\mathrm{s^{-1}}$. Using (38)–(39), we set the reference potential
temperature and the Exner pressure with $\Theta_{0}=280\,\mathrm{K}$. The
computational domain is defined by

$\displaystyle{\Omega}$ $\displaystyle=\left\\{(x_{1},x_{2}),\
g(x_{1})<x_{2}<\,\numprint{21000},\
\,\numprint{-25000}<x_{1}<\,\numprint{25000}\right\\},\quad
g(x_{1})=h_{c}\exp(-(x_{1}/a_{c})^{2})\cos^{2}(\pi\,x_{1}/\lambda_{c}),$ (44)

where with the parameters $h_{c}=\,\numprint{250}\,\mathrm{m}$,
$\lambda_{c}=\,\numprint{4000}\,\mathrm{m}$ and
$a_{c}=\,\numprint{5000}\,\mathrm{m}$. No-flux boundary conditions are used
along the bottom surface
$\\{(x_{1},g(x_{1})),x_{1}\in(\,\numprint{-25000},\,\numprint{25000})\\}$,
while non-reflecting boundary conditions are used on the rest of boundary
${\Gamma}$.

Since the flow is steady, we adopt Algorithm 1. In step 13, we set
$\eta_{m}^{\mathrm{int}}:=\left\|\widehat{q}_{h\tau}(t_{m})-\Pi\widehat{q}_{h\tau}(t_{m})\right\|_{L^{2}({\Omega})}$
and carry out the computation formally for $T\to\infty$. Moreover, we do not
repeat the time steps if the stopping criterion (step 14) is not met and
proceed directly to the next time step. Furthermore, the constant $c_{T}$ in
(32) is set large since we do not take care of the time discretization error,
hence we replace (31) by
$\eta^{\scriptscriptstyle\mathrm{A}}_{m}({{\bm{w}}}_{h\tau})\leq
c_{A}\eta^{\scriptscriptstyle\mathrm{S}}_{m}({{\bm{w}}}_{h\tau})$. The
computation is stopped when $\eta_{m}^{\mathrm{int}}\leq\omega$.

Figure 9 shows isolines of the horizontal and vertical components of the
velocity after 20 levels of the mesh adaptation when the tolerance
$\omega=10^{-5}$ has been achieved. The corresponding $hp$-mesh with a zoom of
the bottom boundary part is presented in Figure 10. We note that the mesh
triangles adjacent to the non-polygonal part of the boundary are treated as
the curvilinear elements, namely we use a piecewise cubic approximation. We
observe small elements near the bottom part of the boundary which are
necessary to capture its geometry. On the other hand, the mesh is mostly
quasi-uniform with large elements having the highest polynomial approximation
degree in the majority of the domain since this example does not contain any
strong anisotropic feature. Again, these results are in agreement with the
reference data in [14, Section 6.4].

Figure 9: Schär mountain: zoom of isolines of the horizontal (left) and
vertical (right) components of velocity with the contour interval of $0.2$ and
$0.0005$, respectively.

Figure 10: Schär mountain: the final $hp$-mesh, total view (left) and its zoom
around the bottom part of the boundary (right).

##  7 Conclusion

We presented the $hp\tau$-adaptive space-time discontinuous Galerkin method
for the numerical solution of non-hydrostatic mesoscale atmospheric problems.
The time discontinuous discretization with a weak coupling by the penalty
admits using different meshes at different time levels in a natural way. We
demonstrated by numerical experiments that the complete re-meshing during the
computational process does not significantly affect the accuracy of the
approximation. Moreover, measuring the interpolation error in terms of the
density is sufficient for the anisotropic mesh adaptation algorithm to
generate appropriately adapted meshes, namely to capture physically relevant
quantity as the perturbation of the potential temperature. The robustness of
the adaptive algorithm is supported by numerical simulation of several
benchmarks.

## References

  * [1] Alauzet, F., Loseille, A., Olivier, G.: Time-accurate multi-scale anisotropic mesh adaptation for unsteady flows in CFD.  J. Comput. Phys. 373, 28–63 (2018)
  * [2] Belme, A., Dervieux, A., Alauzet, F.: Time accurate anisotropic goal-oriented mesh adaptation for unsteady flows.  J. Comput. Phys. 231(19), 6323–6348 (2012)
  * [3] Busto, S., Tavelli, M., Boscheri, W., Dumbser, M.: Efficient high order accurate staggered semi-implicit discontinuous Galerkin methods for natural convection problems.  Comput. Fluids 198 (2020)
  * [4] Deuflhard, P.: Newton Methods for Nonlinear Problems, Springer Series in Computational Mathematics, Vol. 35.  Springer (2004)
  * [5] Dolejší, V.: ANGENER – Anisotropic mesh generator, in-house code.  Charles University, Prague, Faculty of Mathematics and Physics (2000).  https://msekce.karlin.mff.cuni.cz/~dolejsi/angen/
  * [6] Dolejší, V.: Semi-implicit interior penalty discontinuous Galerkin methods for viscous compressible flows.  Commun. Comput. Phys. 4(2), 231–274 (2008)
  * [7] Dolejší, V.: ADGFEM – Adaptive discontinuous Galerkin finite element method, in-house code.  Charles University, Prague, Faculty of Mathematics and Physics (2020).  https://msekce.karlin.mff.cuni.cz/~dolejsi/adgfem/
  * [8] Dolejší, V., Feistauer, M.: Discontinuous Galerkin Method – Analysis and Applications to Compressible Flow.  Springer Series in Computational Mathematics 48. Springer, Cham (2015)
  * [9] Dolejší, V., May, G.: Anisotropic $hp$-Mesh Adaptation Methods.  Birkhäuser (2022)
  * [10] Dolejší, V., Roskovec, F., Vlasák, M.: Residual based error estimates for the space-time discontinuous Galerkin method applied to the compressible flows.  Comput. Fluids 117, 304–324 (2015)
  * [11] Dolejší, V., May, G.: An anisotropic $hp$-mesh adaptation method for time-dependent problems based on interpolation error control.  J. Sci. Comput. 95(2) (2023)
  * [12] Feistauer, M., Felcman, J., Straškraba, I.: Mathematical and Computational Methods for Compressible Flow.  Clarendon Press, Oxford (2003)
  * [13] Frey, P.J., Alauzet, F.: Anisotropic mesh adaptation for CFD computations.  Comput. Methods Appl. Mech. Engrg. 194, 5068–5082 (2005)
  * [14] Giraldo, F., Restelli, M.: A study of spectral element and discontinuous Galerkin methods for the Navier–Stokes equations in nonhydrostatic mesoscale atmospheric modeling: Equation sets and test cases.  J. Comput. Phys. 227(8), 3849–3877 (2008)
  * [15] Giraldo, F., Restelli, M., Läuter, M.: Semi-implicit formulations of the Navier–Stokes equations: Application to nonhydrostatic atmospheric modeling.  SIAM J. Sci. Comput. 32(6), 3394–3425 (2010)
  * [16] Giraldo, F., Warburton, T.: A nodal triangle-based spectral element method for the shallow water equations on the sphere.  J. Comput. Phys. 207(1), 129–150 (2005)
  * [17] Giraldo, F.X.: An Introduction to Element-Based Galerkin Methods on Tensor-Product Bases - Analysis, Algorithms and Applications.  Springer (2020)
  * [18] Giraldo, F.X., Rosmond, T.E.: A scalable spectral element eulerian atmospheric model (SEE-AM) for NWP: Dynamical core tests.  Mon. Weather Rev. 132(1), 133–153 (2004)
  * [19] Habashi, W.G., Dompierre, J., Bourgault, Y., Ait-Ali-Yahia, D., Fortin, M., Vallet, M.G.: Anisotropic mesh adaptation: towards user-independent, mesh-independent and solver-independent CFD. Part I: general principles.  Int. J. Numer. Methods Fluids 32(6), 725–744 (2000)
  * [20] Hunter, J.D.: Matplotlib: A 2D graphics environment.  Computing in Science & Engineering 9(3), 90–95 (2007)
  * [21] Jiang, H., Tang, H., Wu, K.: Positivity-preserving well-balanced central discontinuous Galerkin schemes for the Euler equations under gravitational fields.  J. Comput. Phys. 463 (2022)
  * [22] Lukáčová-Medvid’ová, M., Morton, K.W., Warnecke, G.: Finite volume evolution Galerkin methods for hyperbolic systems.  SIAM J. Sci. Comput. 26(1), 1–30 (2004)
  * [23] Navas-Montilla, A., Echeverribar, I.: A family of well-balanced WENO and TENO schemes for atmospheric flows.  J. Comput. Phys. 489 (2023)
  * [24] Orlando, G., Benacchio, T., Bonaventura, L.: An IMEX-DG solver for atmospheric dynamics simulations with adaptive mesh refinement.  J. Comput. Appl. Math. 427 (2023)
  * [25] Reddy, S., Maciej Waruszewski and, F.A.d.B.A., Giraldo, F.X.: Schur complement IMplicit-EXplicit formulations for discontinuous Galerkin non-hydrostatic atmospheric models.  J. Comput. Phys. 491 (2023)
  * [26] Restelli, M., Giraldo, F.X.: A conservative discontinuous Galerkin semi-implicit formulation for the Navier–Stokes equations in nonhydrostatic mesoscale modeling.  SIAM J. Sci. Comput. 31(3), 2231–2257 (2009)
  * [27] Robert, A.: Bubble convection experiments with a semi-implicit formulation of the Euler equations.  J. Atmos. Sci. 50(13), 1865 – 1873 (1993)
  * [28] Schär, C., Leuenberger, D., Fuhrer, O., Lüthic, D., Girard, C.: A new terrain-following vertical coordinate formulation for atmospheric prediction models.  Mon. Weather Rev. 130, 2459–2480 (2002)
  * [29] Skamarock, W., Klemp, J.: Efficiency and accuracy of the Klemp-Wilhelmson time-splitting technique.  Mon. Weather Rev. 122(11), 2623–2630 (1994)
  * [30] Steppeler, J., Hess, R., Doms, G., Schättler, U., Bonaventura, L.: Review of numerical methods for nonhydrostatic weather prediction models.  Meteorol. Atmos. Phys. pp. 287–301 (2003)
  * [31] Straka, J., Wilhelmson, R.B., Wicker, L.J., Anderson, J.R., Droegemeier, K.K.: Numerical solutions of a non‐linear density current: A benchmark solution and comparisons.  Int. J. Numer. Meth. Fluids 17(1), 1– 22 (1993)
  * [32] Waruszewski, M., Kozdon, J., Wilcox, L., Gibson, T., Giraldo, F.: Entropy stable discontinuous Galerkin methods for balance laws in non-conservative form: Applications to the Euler equations with gravity.  J. Comput. Phys. 468 (2022)
  * [33] Yelash, L., Müller, A., Lukáčová-Medvid’ová, M., Giraldo, F., Wirth, V.: Adaptive discontinuous evolution Galerkin method for dry atmospheric flow.  J. Comput. Phys. 268, 106 – 133 (2014)