[Paper Review] Shenfun -- automating the spectral Galerkin method

Mikael Mortensen|arXiv (Cornell University)|Aug 10, 2017

Scientific Research and Discoveries5 references3 citations

TL;DR

Shenfun is a Python-based, high-level software framework that automates the spectral Galerkin method for solving partial differential equations (PDEs) on simple tensor product domains with one non-periodic and any number of periodic directions. It uses a FEniCS-like variational form interface, enables efficient parallel computation via MPI and FFTW, and achieves high accuracy and performance—demonstrated by solving linear and nonlinear PDEs such as the Ginzburg-Landau equation with de-aliasing and symmetry preservation.

ABSTRACT

With the shenfun Python module (github.com/spectralDNS/shenfun) an effort is made towards automating the implementation of the spectral Galerkin method for simple tensor product domains, consisting of (currently) one non-periodic and any number of periodic directions. The user interface to shenfun is intentionally made very similar to FEniCS (fenicsproject.org). Partial Differential Equations are represented through weak variational forms and solved using efficient direct solvers where available. MPI decomposition is achieved through the {mpi4py-fft} module (bitbucket.org/mpi4py/mpi4py-fft), and all developed solver may, with no additional effort, be run on supercomputers using thousands of processors. Complete solvers are shown for the linear Poisson and biharmonic problems, as well as the nonlinear and time-dependent Ginzburg-Landau equation.

Motivation & Objective

• To simplify the implementation of the spectral Galerkin method for researchers working on PDEs in simple geometries.
• To provide a high-level, user-friendly interface similar to FEniCS but tailored for spectral methods.
• To enable efficient, scalable parallel computation of PDEs using MPI and optimized FFTs on supercomputers.
• To support accurate solution of complex PDEs such as the nonlinear time-dependent Ginzburg-Landau equation with de-aliasing.
• To bridge the gap between high-level Python scripting and high-performance computing in spectral methods, avoiding low-level language constraints.

Proposed method

• The spectral Galerkin method is implemented using weak variational forms, with PDEs expressed in a FEniCS-like syntax.
• Basis functions (e.g., Chebyshev, Legendre, Fourier) are used with global support and weighted inner products for weak forms.
• The package uses Numpy for vectorized operations and mpi4py-fft for MPI decomposition across thousands of processors.
• Padded spectral transforms (3/2 rule) are applied to handle nonlinear terms and prevent aliasing in spectral space.
• A fourth-order Runge-Kutta time integrator is used for time-dependent problems, with de-aliased nonlinear terms computed via padded transforms.
• The framework leverages pyFFTW and FFTW for high-performance FFTs, matching low-level C/Fortran performance in Python.

Experimental results

Research questions

• RQ1 Can a high-level Python interface be designed to automate the spectral Galerkin method for tensor product domains?
• RQ2 Can such a framework achieve performance comparable to hand-optimized low-level codes on supercomputers?
• RQ3 How effectively can de-aliasing be integrated into time-dependent nonlinear PDE solvers using spectral methods?
• RQ4 Can symmetry-preserving solutions be reliably computed in spectral methods when using different FFT algorithms?
• RQ5 To what extent can the FEniCS-like interface be adapted to spectral methods without sacrificing computational efficiency?

Key findings

• Shenfun enables the implementation of spectral Galerkin solvers for PDEs with minimal code, closely mirroring mathematical formulations.
• The framework successfully solves the 2D time-dependent Ginzburg-Landau equation with high accuracy and symmetry preservation.
• De-aliasing via 3/2-padded transforms is essential to maintain solution symmetry; without it, symmetry is lost due to aliasing.
• The code achieves performance comparable to low-level C++ or Fortran implementations when optimized with Cython and FFTW.
• The framework has been used to implement Navier-Stokes solvers in the spectralDNS project, running on up to 2048³ meshes on the Shaheen II supercomputer.
• Different FFTW algorithms can affect numerical symmetry, highlighting the sensitivity of spectral solvers to implementation details.