Dedalus: A flexible framework for numerical simulations with spectral methods

Numerical solutions of partial differential equations enable a broad range of scientific research. The Dedalus project is a flexible, open-source, parallelized computational framework for solving general partial differential equations using spectral methods. Dedalus translates plain-text strings describing partial differential equations into efficient solvers. This paper details the numerical method that enables this translation, describes the design and implementation of the codebase, and illustrates its capabilities with a variety of example problems. The numerical method is a first-order generalized tau formulation that discretizes equations into banded matrices. This method is implemented with an object-oriented design. Classes for spectral bases and domains manage the discretization and automatic parallel distribution of variables. Discretized fields and mathematical operators are symbolically manipulated with a basic computer algebra system. Initial value, boundary value, and eigenvalue problems are efficiently solved using high-performance linear algebra, transform, and parallel communication libraries. Custom analysis outputs can also be specified in plain text and stored in self-describing portable formats. The performance of the code is evaluated with a parallel scaling benchmark and a comparison to a finite-volume code. The features and flexibility of the codebase are illustrated by solving several examples: the nonlinear Schrödinger equation on a graph, a supersonic magnetohydrodynamic vortex, quasigeostrophic flow, Stokes flow in a cylindrical annulus, normal modes of a radiative atmosphere, and diamagnetic levitation.

Figure 1

Chebyshev polynomials (a) can be viewed as cosine functions (b) drawn on a cylinder and projected onto the bisecting plane (c).

Figure 2

Derivative matrices using (a) Chebyshev-$T$ and (b) Chebyshev-$U$ polynomials as test functions. Using different families of test and trial functions allows general differential operators to be represented with sparse matrices.

Figure 3

Conceptual stages in matrix construction for the Chebyshev discretization of Poisson's equation with Dirichlet and Neumann boundary conditions (matrix entries colored by sign). (1) Original (T-to-T) matrices with block columns corresponding to $u$ and $u_x$, and block rows corresponding to the LHS of ${\partial }_x{u}_x=0$ and ${\partial }_{x}u-u_x=f\left(x\right)$. The T-to-T differentiation matrices are dense upper triangular. (2) After T-to-U conversion and the addition of boundary conditions via the tau method. The derivatives are sparse and the identity-block bandwidths increase slightly. The boundary conditions involve all coefficients of $u$ and $u_x$. (3) After Dirichlet recombination. Boundary rows are sparse and the equation-block bandwidths increase slightly. (4) After reversing the Kronecker product to group by modes rather than variables. This final matrix is highly sparse and completely banded.

Figure 4

Example collocation grids for the implemented bases. The $\mathtt{Fourier}$ basis grid is evenly spaced on a periodic interval. The $\mathtt{Chebyshev}$ and $\mathtt{Legendre}$ basis grids cluster near the ends of a finite interval. The $\mathtt{Hermite}$ basis grid spreads out over the real line. The $\mathtt{Laguerre}$ basis grid spreads one way over the half line. The $\mathtt{Compound}$ basis concatenates other bases on adjascent segments; e.g., three Chebyshev segments.

Figure 5

The parallel data distribution for 3D data over a process mesh of shape $(4,2)$. The global data is depicted as being split into the portions that are local to each process. The dimensions are labeled e.g., $k_x$ when the corresponding dimension is in coefficient space, and, e.g., $x$ in grid space. The transforms (TF) and transpositions (TP) stepping between layouts are indicated.

Figure 6

The effective data redistribution that occurs during the distributed transposition between layouts 1 and 2 of the example shown in Fig. 5. This transposition is switching the second mesh axis with $M_2=2$ from distributing $k_y$ to distributing $z$.

Figure 7

An operator tree representing the expression $\mathtt{2}*(\mathtt{10}+\mathtt{f})+\mathtt{f}*\mathtt{dx}\left(\mathtt{f}\right)$.

Figure 8

Dedalus performance (mode-iterations per core-second) vs parallel scaling (pencils per core) from 32 to 2048 cores for 3D Boussinesq (a) hydrodynamics and (b) magnetohydrodynamics. Colors correspond to different 3D (Fourier-Fourier-Chebyshev) resolutions. For each resolution, the open circle indicates the run on 512 cores, with the core count increasing by factors of two to the right. Efficient strong scaling is seen down to roughly 8 pencils per core. Weak scaling shows the expected $\propto 1/{log}_{2}N$ efficiency for FFT-based computations. Execution timescales linearly with the number of problem variables.

Figure 9

Snapshots of a moderate Mach-number Kelvin-Helmholtz instability test problem simulated at various resolutions with a finite-volume code (Athena) and Dedalus. The finite-volume method introduces small errors which trigger unphysical secondary instabilities in the vortex rolls. These spurious instabilities disappear as the simulation resolution is increased. Quantitative comparisons show comparable accuracy between the finite-volume method with ${16384}^2$ degrees of freedom and Dedalus with ${2048}^2$ degrees of freedom. Figure adapted from Lecoanet et al. [24].

Figure 10

Evolution of the nonlinear Schrödinger equation on a network, simulated by coupling the boundaries of different fields on a 1D Chebyshev segment. A soliton initially isolated to one segment scatters at the vertices and fills the network over time.

Figure 11

Density at $t=0.5$ in the Orszag-Tang vortex test. The sharp density jumps visualize the MHD shocks. The structures in the upper left and lower right corners are due to shock-shock interactions, which are correctly treated in Dedalus.

Figure 12

A profile of the density at $t=0.5$ and $y=0.3125$ in the Orszag-Tang vortex test. The inset shows the density profile around $x=0.7$, with crosses denoting grid points (using the nondealiased grid). The shock is resolved by several grid points.

Figure 13

Surface and lateral slices of the PV (top) and buoyancy perturbation $\theta$ (bottom) in a 3D quasigeostrophic flow. Both images are at $t=200$ in the statistically saturated state.

Figure 14

Evolution of the tracer fields in a reversible Taylor-Couette Stokes flow. The three separate tracer fields are visualized as red, green, and blue in a single RGB image. As time progresses, we rotate the inner cylinder first counter clockwise and then clockwise. By the end of the simulation, the flow has returned to its initial conditions, modified only by mass diffusion of the dye.

Figure 15

(a) Structure of a balanced atmosphere with an Eddington closure for radiation transport and a Kramer-like opacity law with $a=1,b=0$. At depth, the atmosphere is nearly an adiabatic polytrope with constant $s/C_P$ and $N^2\approx 0$. The upper atmosphere is nearly isothermal with exponentially decaying $\rho$ and $P$. (b) The upper atmosphere is a resonant cavity for internal gravity waves with $N^2>0$. Frequencies are normalized by the Brunt-Väisälä frequency $N$ in the isothermal layer.

Figure 16

Eigenmode (a) frequencies and (b) periods for a radiative atmosphere. The mode frequencies $\omega$ are normalized by the Brunt-Väisälä frequency $N$ in the isothermal layer. Acoustic modes (“ac”, blue), gravity modes (“gw”, red), and f-modes (“f”, black) are identified by their frequencies. Acoustic waves are expected to have nearly equal frequency spacing, while gravity waves have nearly equal period spacing, as is observed. Colored circles highlight modes shown in Fig. 17.

Figure 17

Vertical velocity eigenfunctions for a radiative atmosphere, scaled by $\sqrt{\rho }$, including (a) acoustic modes and (b) gravity modes. Gravity modes are almost entirely confined to the isothermal portion of the atmosphere and are evanescent in the adiabatically stratified deep interior. Acoustic modes propagate in both regions, with high frequency waves propagating deeper.

Figure 18

Evolution of a diamagnetic object levitating in a magnetic field. The purple ellipse shows the solid-body region. The black lines show magnetic field lines the same way iron filings trace the magnetic field from a strong bar magnet. Green dots show the time-dependent trajectory of the ellipse's center of mass at equally spaced time intervals.