The SpectralDNS project:
Shenfun  automating the spectral Galerkin method
Mikael Mortensen, Miroslav Kuchta
Department of Mathematics, University of Oslo
and
Lisandro Dalcin
King Abdullah University of Science and Technology
International Conference on Computational Science and Engineering, Simula 23/10  2017
Background to
the spectralDNS/shenfun project
Started out as a toy project for
myself and Hans Petter
Using Fourier transformed Navier Stokes equations
We knew Direct Numerical Simulations (DNS) of turbulence were performed with lowlevel Fortran or C
Would it be possible to write a competitive, state of the art turbulence solver in Python?
From scratch?
(Not just wrapping a Fortran solver)
Using only standard Python modules numpy/scipy/mpi4py?
Can we make this solver run on thousands of cores ?
With billions of unknowns?
As fast as C/fortran?
Long story short  Yes
 Vectorize everything (Numpy)
 Using FFTW we implement parallel (MPI) Fast Fourier Transforms (FFT)
 100 lines of pure Python / Numpy ~30% slower than optimized C++ code
 Using Cython, we were par with C++
65,536 cores!
Gist
Concept proven
and Hans Petter is off to write 4 new books,
but with new collaborators, and fast transforms in place we can do so much more
mpi4pyfft
https://bitbucket.org/mpi4py/mpi4pyfft
With Lisandro Dalcin at KAUST
~ 500 lines of Python code
With the newly released
we can do any transform from FFTW in parallel
Highly configurable and fast
We can now implement
the spectral Galerkin method
on nonperiodic domains, with different types
of boundary conditions (Dirichlet/Neumann)
Not just periodic Fourier
The spectral Galerkin method employs global basis functions with the Galerkin approximation

Choose global basis (test) functions satisfying the correct boundary conditions. Usually modified Chebyshev or Legendre polynomials, or Fourier exponentials

Transform PDEs to weak variational forms using the method of weighted residuals

Assemble and solve system of equations
Solution
Expansion coefficients
(unknowns)
Global basis functions
and it works in three steps:
Spectral Galerkin
Poisson example
1) Create a basis that satisfies boundary conditions, e.g.,
(Legendre polynomial order k+2)
2) Transform PDE to variational form using test function
3) Assemble and solve linear system of equations
Trial function
Test function
With quadrature notation
3 steps  Sounds generic and easy?
Let's automate!
from fenics import *
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'CG', 1)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(6.0)
a = inner(grad(u), grad(v))*dx
L = f*v*dx
 Automates the finite element method (Galerkin)
 Basis functions with local support

Finite accuracy (not spectral)
How?
Shenfun
https://github.com/spectralDNS/shenfun
Python package (~3k loc) for automating the spectral Galerkin
method for high performance computing on tensor product grids
 A FEniCS like interface with a simple form language
 MPI  invisible to the user: mpi4pyfft (with Lisandro Dalcin)
 Fourier  Chebyshev  Legendre modified basis functions
 Extremely fast (!), order optimal, direct solvers for Poisson, Helmholtz and Biharmonic problems (with Miroslav Kuchta)
Shenfun implementation of 1D Poisson problem
import numpy as np
from shenfun import TestFunction, TrialFunction, inner, div, grad, chebyshev
N = 20
V = chebyshev.ShenDirichletBasis(N, plan=True)
v = TestFunction(V)
u = TrialFunction(V)
# Assemble stiffness matrix
A = inner(v, div(grad(u)))
# Some consistent right hand side
f = Function(V, False)
f[:] = 2
# Assemble right hand side
f_hat = inner(v, f)
# Solve equations and transform spectral solution to real space
f_hat = A.solve(f_hat)
fj = V.backward(f_hat)
x = V.mesh(N)
assert np.allclose(fj, 1x**2)
Very similar to FEniCS!
just with spectral accuracy
(1) Create basis
(2) Variational forms
(3) Solve system
2D/3D/4D.. Tensor product spaces
Cartesian products of Fourier and Chebyshev/Legendre bases
Any Cartesian combination of spaces (ultimately)
With solutions of the form:
Unknown expansion coefficients
Tensor product spaces
Implementation
from shenfun import *
from mpi4py import MPI
comm = MPI.COMM_WORLD
N = (20, 30)
# Create two bases
V0 = chebyshev.ShenDirichletBasis(N[0])
V1 = fourier.R2CBasis(N[1])
# Create tensor product space
T = TensorProductSpace(comm, (V0, V1))
v = TestFunction(T)
u = TrialFunction(T)
A = inner(v, div(grad(u)))
# Or 3D
V2 = fourier.C2CBasis(40)
T3 = TensorProductSpace(comm, (V0, V2, V1))
v = TestFunction(T3)
u = TrialFunction(T3)
A = inner(v, div(grad(u)))
MPI under the hood
using mpi4pyfft
on tensor product grids like:
Slab
Pencil
A simple form language
influenced heavily by UFL/FEniCS
Mathematical operators act on test, trial, or functions with data
V = TensorProductSpace(comm, (V0, V1, V2))
# Form arguments:
v = TestFunction(V)
u = TrialFunction(V)
u_ = Function(V, False)
# Forms
du = grad(u_)
d2u = div(du)
cu = curl(du)
dudx = Dx(u_, 0)
Inner products and projections
A = inner(v, u) # Bilinear form, assembles to matrices
f = inner(v, u_) # Linear form, assembles to vector
project(div(grad(u_)), V) # projection to V
Fast transforms
u_hat = V.forward(u, u_hat) u = V.backward(u_hat, u)
Move back and forth between and using fast transforms
Periodic basis using fast Fourier transforms
Nonperiodic basis using fast Chebyshev transforms (not for Legendre)
Nonlinearities
all treated with pseudospectral techniques
(1) Compute inverse
u = T.backward(u_hat)
with dealiasing using
3/2rule (padding) or 2/3rule
(3) Transform back to spectral space
w = u*u
(2) Compute nonlinearity
w_hat = T.forward(w)
T = TensorProductSpace(comm, (V0, V1, V2))
Fast direct solvers
Fourier: A is diagonal
Modified Legendre: A is diagonal
Modified Chebyshev: A sparse  solved as tridiagonal matrix
Fast direct solvers
2 nonperiodic directions
Matrix form:
Sparse matrices A (stiffness) and B (mass) are 1D matrices of basis
Very fast direct solvers implemented in parallel using matrix decomposition method
Linear algebra:
Complex GinzburgLandau
example
2D periodic (x, y) \( \in [50, 50]^2 \)
Use Fourier in both directions
Weak form, only spatial discretization:
Complex GinzburgLandau
example
2D periodic (x, y) \( \in [50, 50]^2 \)
import sympy
from shenfun.fourier.bases import R2CBasis, C2CBasis
from shenfun import inner, grad, TestFunction, TrialFunction, \
Function, TensorProductSpace
# Use sympy to set up initial condition
x, y = sympy.symbols("x,y")
ue = (x + y)*sympy.exp(0.03*(x**2+y**2))
N = (200, 200)
K0 = C2CBasis(N[0], domain=(50., 50.))
K1 = C2CBasis(N[1], domain=(50., 50.))
T = TensorProductSpace(comm, (K0, K1))
u = TrialFunction(T)
v = TestFunction(T)
U = Function(T, False)
def compute_rhs(rhs, u_hat, U, T):
rhs = inner(grad(v), grad(U), output_array=rhs, uh_hat=u_hat)
rhs /= (2*np.pi)**2
rhs += u_hat
rhs = T.forward((1+1.5j)*U*abs(U)**2, w0)
return rhs
# And, e.g., RungeKutta integrator
#
#
Key code:
Navier Stokes
Turbulent channel flow
3D tensor product space
2 Fourier bases +
Modified Chebyshev bases for wall normal direction
Kim and Moin modified version of NavierStokes:
is convection
Wall normal velocity
Wall normal vorticity
Navier Stokes
Turbulent channel flow
First fully spectral DNS with
3D Tensor product mesh:
Shaheen II at KAUST
Some more demos prepared with Hans Petter's doconce:
Thank you for your time!
http://github.com/spectralDNS/shenfun
Questions?
Presented at the International Conference on Computational Science and Engineering 23/102017
By Mikael Mortensen
Presented at the International Conference on Computational Science and Engineering 23/102017
 334