🔢 mathematics

An IMEX-DG solver with non-conforming mesh refinement for atmospheric dynamics with rotation

This paper presents a high-order implicit-explicit discontinuous Galerkin solver on adaptive non-conforming meshes for simulating rotating atmospheric flows, featuring a novel discrete treatment of rotational source terms and demonstrating accuracy and stability through convergence studies and 3D stratified flow simulations.

Original authors: Letizia Bottani, Tommaso Benacchio, Giuseppe Orlando, Luca Bonaventura, Allan Peter Engsig-Karup

Published 2026-05-29

📖 4 min read🧠 Deep dive

CC BY 4.0

Original authors: Letizia Bottani, Tommaso Benacchio, Giuseppe Orlando, Luca Bonaventura, Allan Peter Engsig-Karup

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine the Earth's atmosphere as a giant, invisible ocean of air. Scientists try to predict how this ocean moves using super-computers. But there's a catch: this "ocean" has two very different types of waves happening at the same time.

First, there are fast, jittery sound waves (like the pop of a balloon) that move incredibly quickly. Second, there are slow, rolling weather waves (like the gentle rise and fall of a tide) that carry the actual weather patterns we care about.

The Problem: The "Speed Trap"

In the past, computer models had to take tiny, tiny steps to catch those fast sound waves. If they took a step too big, the math would explode. This meant the computer spent 99% of its time calculating sound waves that don't really affect the weather, and only 1% of its time figuring out the actual wind and rain. It was like trying to watch a slow-motion movie of a glacier, but your camera was forced to take 1,000 pictures every second just to keep the lens from shaking.

The Solution: The "IMEX" Magic Trick

The authors of this paper built a new type of solver (a math engine) called IMEX-DG. Think of it as a smart camera that knows how to handle the two speeds differently:

The "IM" (Implicit) part: It handles the fast, jittery sound waves with a special "time-travel" math trick that lets it take big steps without the picture blurring.
The "EX" (Explicit) part: It handles the slow, important weather waves normally, step-by-step.

This allows the computer to zoom through time quickly, focusing on the weather that matters, while still keeping the math stable.

The New Challenge: The Spinning Earth

The Earth isn't just sitting still; it's spinning. This spin creates a force called the Coriolis effect (which makes hurricanes swirl and winds curve). Adding this spin to the math is tricky. It's like trying to calculate the path of a ball thrown on a merry-go-round while the ball is also bouncing off walls.

The authors had to figure out how to add this "spin" into their fast-and-slow math engine without breaking it. They tested two different ways to do this:

Method R1: Treat the spin as a "side note" that gets calculated after the main math is done.
Method R2: Bake the spin directly into the main math engine so it's solved at the same time as the pressure and wind.

The Test: The "Mountain" and the "Wave"

To see if their new engine worked, they ran two types of tests:

The Ripple Test: They created a perfect, theoretical wave in a computer box and watched it travel. They checked if the wave kept its shape and speed over time. Both methods (R1 and R2) worked perfectly, proving the math was accurate.
The Mountain Test: They simulated wind blowing over a 3D mountain.
- Without Spin: The wind waves went straight over the mountain, looking like a mirror image on both sides.
- With Spin: The wind waves got twisted! Because the Earth is spinning, the waves didn't look the same on the left and right sides. They curved and tilted.

The authors' new solver successfully captured this "twist." It showed that the wind behaves differently when the Earth spins, just like real life, without creating any weird, fake ripples or crashing the computer.

The "Lego" Advantage

One of the coolest features of their tool is that it works on non-conforming meshes. Imagine building a model out of Lego bricks. Usually, you need bricks of the exact same size to make a smooth wall. But this new solver is like having a Lego set where you can use giant bricks for the open sky and tiny, detailed bricks right around the mountain peak, all without the wall falling apart. This saves a massive amount of computer power because it only does the hard math where it's actually needed.

The Bottom Line

The authors have successfully upgraded their weather simulation engine to include the Earth's spin. They proved that:

It is accurate (it catches the waves correctly).
It is stable (it doesn't crash over long periods).
It is efficient (it runs fast on supercomputers, even with the extra "spin" math).
It handles both smooth skies and detailed mountain terrain equally well.

Essentially, they built a better, faster, and more realistic way to simulate how our spinning planet's atmosphere moves.

Technical Summary: An IMEX-DG Solver with Non-Conforming Mesh Refinement for Atmospheric Dynamics with Rotation

Problem Statement
Atmospheric flows are characterized by the interaction of multiple spatial and temporal scales, ranging from large-scale balanced circulations to fast acoustic and gravity waves. In fully compressible atmospheric models, the coexistence of slow geostrophic modes and fast acoustic modes introduces severe stiffness, which strongly constrains the time step size in fully explicit integration methods. While Implicit-Explicit (IMEX) time integration and high-order Discontinuous Galerkin (DG) spatial discretization have been successfully combined to address acoustic stiffness in non-rotating Cartesian configurations, extending these frameworks to include planetary rotation is essential for realistic numerical weather prediction (NWP) and standard benchmark configurations. The inclusion of the Coriolis force introduces additional complexity, particularly regarding the preservation of geostrophic balance and the efficient implementation of rotational terms within matrix-free, high-order solvers.

Methodology
The authors present an extension of an existing high-order IMEX-DG solver, implemented within the deal.II finite element library, to handle the compressible Euler equations with rotational effects.

Governing Equations: The solver addresses the inviscid, compressible Euler equations on a bounded domain, incorporating the Coriolis acceleration ( $-f \mathbf{k} \times \mathbf{u}$ ) and gravitational acceleration. The system is closed using a calorically perfect ideal gas equation of state.
Temporal Discretization: A second-order IMEX Additive Runge–Kutta (ARK) scheme (specifically the TR-BDF2 method) is employed. Stiff contributions (pressure gradients, gravity, and rotation) are treated implicitly, while non-stiff advective terms are treated explicitly. This allows the time step to be determined by meteorologically relevant flow scales rather than the speed of sound.
Spatial Discretization: A high-order DG method is used on non-conforming adaptive meshes (quadrilateral in 2D, hexahedral in 3D). The formulation utilizes a nodal basis on tensor-product elements to enable efficient matrix-free operator evaluation.
Rotational Term Treatment: Two alternative strategies for handling the rotational term within the nonlinear fixed-point (Picard) iteration are introduced and compared:
1. Strategy R1: Rotational and gravitational terms are lagged to the right-hand side, evaluated using the previous Picard iterate.
2. Strategy R2: Rotational and gravitational terms are treated implicitly on the left-hand side of the coupled momentum-energy system.
Matrix-Free Implementation: A key technical contribution is the derivation of a formulation for the inverse of the operator $(A + R)$ , where $A$ is the mass matrix and $R$ represents the Coriolis term. By exploiting the algebraic structure of the Coriolis operator (specifically that $R = \beta A J$ , where $J$ is a planar rotation), the authors derive a factorization $(A+R)^{-1} = (I + \beta J)^{-1}A^{-1}$ . This allows the application of the inverse to be computed via standard matrix-free operations followed by a local algebraic transformation, avoiding the need to solve linear systems with the full coupled operator.

Key Contributions

Extension to Rotating Flows: The successful integration of the Coriolis force into a fully compressible, high-order IMEX-DG framework within the deal.II library.
Efficient Matrix-Free Rotational Solver: The development of a specific algebraic factorization for the Coriolis operator that preserves the efficiency of matrix-free implementations, avoiding the computational cost of solving large coupled linear systems for the rotational term.
Comparison of Iterative Strategies: A systematic comparison of two nonlinear fixed-point treatments (R1 and R2) for rotational terms, evaluating them in terms of accuracy, robustness, and computational efficiency.
Validation on Adaptive Meshes: The first application of such a rotating IMEX-DG framework on non-conforming adaptive meshes, demonstrating the capability to perform local refinement near orography while maintaining stability and accuracy.

Numerical Results
The proposed formulation was validated through a series of two- and three-dimensional benchmarks:

2D Inertia-Gravity Waves: Convergence studies on a standard inertia-gravity wave benchmark showed that both R1 and R2 strategies achieve near-optimal convergence rates consistent with the theoretical order of the scheme. The results confirmed that the solver accurately captures wave propagation and the preservation of balanced dynamics.
Planetary-Scale Waves: A large-domain test (60,000 km horizontal extent) demonstrated the solver's robustness over long integration times (96 hours). The discrete geostrophic equilibrium error remained bounded and below $10^{-7}$ , indicating that the lack of an exactly well-balanced discretization did not lead to spurious oscillations or instability.
3D Flow Over Orography:
- Medium-Steep Hill: Simulations on a 60x40x16 km domain showed that the inclusion of rotation correctly breaks the axial symmetry of the wake and wave patterns, consistent with physical expectations. The results were consistent between uniform and non-conforming meshes, with the latter providing better resolution of gradients near the mountain at reduced computational cost.
- Large-Scale Hill: A new, larger-scale test (240x240x30 km) highlighted the pronounced asymmetry induced by rotation in the velocity fields, which was clearly captured without spurious oscillations.
Parallel Scalability: Strong-scaling tests on up to 2048 cores (using the MeluXina supercomputer) demonstrated that the inclusion of rotational terms does not significantly degrade the parallel efficiency of the original non-rotating solver. The solver maintained close-to-ideal scaling behavior, and the non-conforming mesh configuration showed comparable scalability to the uniform mesh case.

Significance and Claims
The paper claims that the proposed rotating IMEX-DG framework successfully extends the capabilities of existing compressible atmospheric solvers to include planetary rotation without compromising the accuracy, robustness, or scalability of the original method. The authors emphasize that the matrix-free treatment of the rotational operator is a critical enabler for efficiency, allowing the solver to handle the additional complexity of the Coriolis force with minimal overhead.

The results indicate that the framework is suitable for simulating large-scale atmospheric flows where rotation is a leading-order physical effect. The comparison between the two rotational treatments suggests they are numerically equivalent for the tested cases, providing flexibility in implementation. The authors conclude that the method correctly captures the asymmetry and wave structures induced by rotation, validating its use for both standard benchmarks and more realistic, large-scale atmospheric configurations. Future work is identified as focusing on planetary-scale 3D configurations, detailed comparisons with exactly well-balanced discretizations, and extensions to spherical geometries.