In [1]:
from shenfun import *
print('hello world')
approximates solutions $u(x)$ using global trial functions $\phi_k(x)$ and unknown expansion coefficients $\hat{u}_k$
$$ u(x) = \sum_{k=0}^{N-1}\hat{u}_k \phi_k(x) $$Multidimensional solutions are formed from outer products of 1D bases
$$ u(x, y) = \sum_{k=0}^{N_0-1}\sum_{l=0}^{N_1-1}\hat{u}_{kl} \phi_{kl}(x, y)\quad \text{ or }\quad u(x, y, z) = \sum_{k=0}^{N_0-1}\sum_{l=0}^{N_1-1} \sum_{m=0}^{N_2-1}\hat{u}_{klm} \phi_{klm}(x, y, z) $$where, for example
$$ \begin{align} \phi_{kl}(x, y) &= T_k(x) L_l(y)\\ \phi_{klm}(x, y, z) &= T_k(x) L_l(y) \exp(\text{i}mz) \end{align} $$$T_k$ and $L_k$ are Chebyshev and Legendre polynomials.
solves PDEs, like Poisson's equation
\begin{align} \nabla^2 u(x) &= f(x), \quad x \in [-1, 1] \\ u(\pm 1) &= 0 \end{align}using variational forms by the method of weighted residuals. I.e., multiply PDE by a test function $v$ and integrate over the domain. For Poisson this leads to the problem:
Find $u \in H^1_0$ such that
$$(\nabla u, \nabla v)_w^N = -(f, v)_w^N \quad \forall v \in H^1_0$$Here $(u, v)_w^{N}$ is a weighted inner product and $v(=\phi_j)$ is a test function. Note that test and trial functions are the same for the Galerkin method.
The weighted inner product is defined as
$$ (u, v)_w = \int_{\Omega} u \overline{v} w \, d\Omega, $$where $w(\mathbf{x})$ is a weight associated with the chosen basis (different bases have different weights). The overline represents a complex conjugate (for Fourier).
$\Omega$ is a tensor product domain spanned by the chosen 1D bases.
1D with Chebyshev basis:
$$ (u, v)_w ^N = \sum_{i=0}^{N-1} u(x_i) v(x_i) \omega_i \approx \int_{-1}^1 \frac{u v}{\sqrt{1-x^2}} \, {dx}, $$where $\{\omega_i\}_{i=0}^{N-1}$ are the quadrature weights associated with the chosen basis and quadrature rule. The associated quadrature points are denoted as $\{x_i\}_{i=0}^{N-1}$.
2D with mixed Chebyshev-Fourier:
$$ (u, v)_w^N = \int_{-1}^1\int_{0}^{2\pi} \frac{u \overline{v}}{2\pi\sqrt{1-x^2}} \, {dxdy} \approx \sum_{i=0}^{N_0-1}\sum_{j=0}^{N_1-1} u(x_i, y_j) \overline{v}(x_i, y_j) \omega^{(x)}_i \omega_j^{(y)} , $$
In [2]:
from shenfun import *
N = 8
C = FunctionSpace(N, 'Chebyshev', quad='GC', domain=[-2, 2])
L = FunctionSpace(N, 'Legendre')
x, w = C.points_and_weights()
print(L.points_and_weights())
| family | Basis | Boundary condition |
|---|---|---|
| Chebyshev | $$\{T_k-T_{k+2}\}_{k=0}^{N-3}$$ | $$u(\pm 1) = 0$$ |
| Legendre | $$\{L_k-L_{k+2}\}_{k=0}^{N-3}$$ | $$u(\pm 1) = 0$$ |
| Hermite | $$\exp(-x^2)\{H_k\}_{k=0}^{N-1}$$ | $$u(\pm \infty) = 0$$ |
| Laguerre | $$\exp(-x/2)\{La_k-La_{k+1}\}_{k=0}^{N-2}$$ | $$u(0) = u(\infty) = 0$$ |
In [3]:
C0 = FunctionSpace(N, 'Chebyshev', bc=(0, 0))
L0 = FunctionSpace(N, 'Legendre', bc=(0, 0))
H0 = FunctionSpace(N, 'Hermite')
La = FunctionSpace(N, 'Laguerre', bc=(0, 0))
In [4]:
CN = FunctionSpace(N, 'Chebyshev', bc='Neumann')
LN = FunctionSpace(N, 'Legendre', bc='Neumann')
In [5]:
CB = FunctionSpace(N, 'Chebyshev', bc='Biharmonic')
LB = FunctionSpace(N, 'Legendre', bc='Biharmonic')
In [6]:
L0 = FunctionSpace(N, 'Legendre', bc=(0, 0))
C0 = FunctionSpace(N, 'Chebyshev', bc=(0, 0))
L1 = FunctionSpace(N, 'Legendre')
LL = TensorProductSpace(comm, (L0, L1)) # comm is MPI.COMM_WORLD
CL = TensorProductSpace(comm, (C0, L1))
V = VectorTensorProductSpace(LL) # For vector valued functions
f = Array(LL)
print(f.shape)
In [7]:
L0 = FunctionSpace(N, 'Legendre', bc=(0, 0))
L1 = FunctionSpace(N, 'Chebyshev', bc='Neumann')
F2 = FunctionSpace(N, 'Fourier', dtype='d')
T = TensorProductSpace(comm, (L0, L1, F2))
f = Array(T)
print(f.dtype)
help(project)
In [8]:
L0 = FunctionSpace(N, 'Legendre', bc=(0, 0))
L1 = FunctionSpace(N, 'Legendre')
u = TrialFunction(L0)
v = TestFunction(L0)
uh = Function(L0)
g = Array(L0)
du = grad(u) # vector valued expression
h = div(du) # scalar valued expression
A = inner(Dx(u, 0, 3), v)
print(A.diags().todense())
Function represents the solutionuh = Function(L0)
The function evaluated for all quadrature points, $\{x_j\}_{j=0}^{N-1}$, is an Array
uj = Array(L0)
There is a (fast) backward transform for moving from Function to Array, and a forward transform to go the other way.
In [9]:
uj = Array(L0)
uj = uh.backward(uj)
uh = uj.forward(uh)
Project $g(\mathbf{x})$ to $V$:
Find $u$ in $V$ such that:
$$(u, v)_w = (Ig, v)_w \quad \text{for} \, v \in V $$where $Ig$ is $\{g(x_j)\}_{j=0}^{N-1}$, i.e., $g(x)$ evaluated on the quadrature mesh.
Works if $g(x)$ is
Array, which is exactly a Function evaluated on the meshFunction, like div(grad(uh))sympy expression, like sin(x)
In [10]:
dudx = project(Dx(uh, 0, 1), L1) # Compute du/dx
wh = project(uj, L1)
import sympy as sp
x, y = sp.symbols('x,y')
ws = project(sp.sin(4*x), L1)
In [11]:
F0 = FunctionSpace(N, 'F', dtype='D')
F1 = FunctionSpace(N, 'F', dtype='D')
F2 = FunctionSpace(N, 'F', dtype='d')
FF = TensorProductSpace(comm, (F0, F1, F2))
uh = Function(FF)
ua = Array(FF)
ua[:] = np.random.random(ua.shape)
uh = ua.forward(uh)
du = div(grad(uh))
V = VectorTensorProductSpace(FF)
vv = Function(V)
c = curl(vv)
ch = project(c, V)
In [12]:
A = inner(grad(u), grad(v))
In [13]:
print(A)
In [14]:
print(A.diags().todense())
A diagonal stiffness matrix!
In [15]:
# Solve Poisson's equation
import matplotlib.pyplot as plt
from sympy import symbols, sin, cos, lambdify
from shenfun import *
# Use sympy to compute manufactured solution
x, y = symbols("x,y")
ue = sin(6*np.pi*x)*(1-x**2) # `ue` is the manufactured solution
fe = ue.diff(x, 2) # `fe` is Poisson's right hand side for `ue`
SD = FunctionSpace(20, 'L', bc=(0, 0))
u = TrialFunction(SD)
v = TestFunction(SD)
b = inner(v, Array(SD, buffer=fe)) # Array is initialized with `fe`
A = inner(v, div(grad(u)))
uh = Function(SD)
uh = A.solve(b, uh) # Very fast solver due to Jie Shen
print(uh.backward()-Array(SD, buffer=ue))
In [16]:
plt.plot(SD.mesh(), uh.backward())
Out[16]:
In [17]:
L0 = FunctionSpace(N, 'Legendre', bc=(0, 0))
F1 = FunctionSpace(N, 'Fourier', dtype='d')
TP = TensorProductSpace(comm, (L0, F1))
u = TrialFunction(TP)
v = TestFunction(TP)
A = inner(grad(u), grad(v))
In [18]:
print(A)
TPMatrix is a Tensor Product matrixA TPMatrix is the outer product of smaller matrices (2 in 2D, 3 in 3D etc).
Consider the inner product:
$$ \begin{align} (\nabla u, \nabla v)_w &= \frac{1}{2\pi}\int_{-1}^{1}\int_{0}^{2\pi} \left(\frac{\partial u}{\partial x}, \frac{\partial u}{\partial y}\right) \cdot \left(\frac{\partial \overline{v}}{\partial x}, \frac{\partial \overline{v}}{\partial y}\right) {dxdy} \\ (\nabla u, \nabla v)_w &= \frac{1}{2\pi}\int_{-1}^1 \int_{0}^{2\pi} \frac{\partial u}{\partial x}\frac{\partial \overline{v}}{\partial x} {dxdy} + \int_{-1}^1 \int_{0}^{2\pi} \frac{\partial u}{\partial y}\frac{\partial \overline{v}}{\partial y} {dxdy} \end{align} $$which, like A, is a sum of two terms. These two terms are the two TPMatrixes returned by inner above.
Now each one of these two terms can be written as the outer product of two smaller matrices.
Consider the first, inserting for test and trial functions
$$ \begin{align} v &= \phi_{kl} = (L_k(x)-L_{k+2}(x))\exp(\text{i}ly) \\ u &= \phi_{mn} \end{align} $$The first term becomes
$$ \small \begin{align} \int_{-1}^1 \int_{0}^{2\pi} \frac{\partial u}{\partial x}\frac{\partial \overline{v}}{\partial x} \frac{dxdy}{2\pi} &= \underbrace{\int_{-1}^1 \frac{\partial (L_m-L_{m+2})}{\partial x}\frac{\partial (L_k-L_{k+2})}{\partial x} {dx}}_{a_{km}} \underbrace{\int_{0}^{2\pi} \exp(iny) \exp(-ily) \frac{dy}{2\pi}}_{\delta_{ln}} \\ &= a_{km} \delta_{ln} \end{align} $$and the second
$$ \small \begin{align} \int_{-1}^1 \int_{0}^{2\pi} \frac{\partial u}{\partial y}\frac{\partial \overline{v}}{\partial y} \frac{dxdy}{2\pi} &= \underbrace{\int_{-1}^1 (L_m-L_{m+2})(L_k-L_{k+2}) {dx}}_{b_{km}} \underbrace{\int_{0}^{2\pi} ln \exp(iny) \exp(-ily)\frac{dy}{2\pi}}_{l^2\delta_{ln}} \\ &= l^2 b_{km} \delta_{ln} \end{align} $$All in all:
$$ (\nabla u, \nabla v)_w = \left(a_{km} \delta_{ln} + l^2 b_{km} \delta_{ln}\right) $$
In [19]:
A = inner(grad(u), grad(v)) # <- list of two TPMatrices
print(A[0].mats)
print('Or as dense matrices:')
for mat in A[0].mats:
print(mat.diags().todense())
In [20]:
print(A[1].mats[0].diags().todense())
#print(A[1].scale) # l^2
In [21]:
from sympy import symbols, sin, cos, lambdify
from shenfun import *
# Use sympy to compute manufactured solution
x, y, z = symbols("x,y,z")
ue = (cos(4*x) + sin(2*y) + sin(4*z))*(1-x**2)
fe = ue.diff(x, 2) + ue.diff(y, 2) + ue.diff(z, 2)
C0 = FunctionSpace(32, 'Chebyshev', bc=(0, 0))
F1 = FunctionSpace(32, 'Fourier', dtype='D')
F2 = FunctionSpace(32, 'Fourier', dtype='d')
T = TensorProductSpace(comm, (C0, F1, F2))
u = TrialFunction(T)
v = TestFunction(T)
# Assemble left and right hand
f_hat = inner(v, Array(T, buffer=fe))
A = inner(v, div(grad(u)))
# Solve
solver = chebyshev.la.Helmholtz(*A) # Very fast solver due to Jie Shen
u_hat = Function(T)
u_hat = solver(u_hat, f_hat)
assert np.linalg.norm(u_hat.backward()-Array(T, buffer=ue)) < 1e-12
print(u_hat.shape)
In [22]:
X = T.local_mesh()
ua = u_hat.backward()
plt.contourf(X[2][0, 0, :], X[0][:, 0, 0], ua[:, 2], 100)
plt.colorbar()
Out[22]:
In [23]:
import subprocess
subprocess.check_output('mpirun -np 4 python poisson3D.py', shell=True)
Out[23]:
Note that Fourier bases are especially attractive because of features easily handled with MPI:
- diagonal matrices
- fast transformsAll treated with pseudo-spectral techniques
$$ \begin{align} \hat{w}_k &= \widehat{u^2}_k \end{align} $$That is, transform Functions to real space Arrays, perform the nonlinear operation there and transform the nonlinear product back to spectral space (to a Function).
3/2-rule or 2/3-rule is possible for dealiasing with Fourier. Not for the remaining bases.
In [24]:
uj = Array(SD)
#uj[:] = np.random.random(uj.shape)
uj = uh.backward(uj)
wh = Function(SD)
wh = SD.forward(uj*uj, wh)
$$ \begin{align} \nabla^2 u(x) + \alpha |\nabla u(x)|^2 &= -1, \quad x \in [-1, 1]\\ u(\pm 1) = 0 &\text{ and } \alpha \in \mathbb{R^+} \end{align} $$
Note that if $\alpha=0$ then obviously $u(x) = a_0+a_1x-0.5x^2$. This solution can be realized with basis function $T_0(x)-T_2(x)$, which equals $2(1-x^2)$, and which is the first basis function in the space
$$ V = \text{span}\{T_k-T_{k+2}\}_{k=0}^{N-3} $$In other words, for $u$ in $V$ the solution is $\hat{u}_0 = 0.25$ and $\hat{u}_k = 0$ for all $k=1, 2, \ldots$.
$$ \begin{align} u(x) &= \sum_{k=0}^{N-3}\hat{u}_k (T_k-T_{k+2}) = 0.25(1-(2x^2-1)) \\ &= 0.5(1-x^2) \end{align} $$This is a good initial guess to find the nonlinear solution when $\alpha>0$! What is a good initial guess for Legendre? You need to iterate to find the solution when $\alpha > 0$!
$$ \begin{align*} \nabla^2 \mathbf{u} - \nabla p &= \mathbf{f} \quad \text{in } \Omega, \quad \quad \Omega = [-1, 1]\times[-1, 1]\\ \nabla \cdot \mathbf{u} &= h \quad \text{in } \Omega \\ \int_{\Omega} p dx &= 0 \\ \mathbf{u}(\pm 1, y) = \mathbf{u}(x, -1) = (0, 0) &\text{ and }\mathbf{u}(x, 1) = (1, 0) \text{ or } ((1-x^2)(1+x^2), 0) \end{align*} $$
Given appropriate spaces $V$ and $Q$ a variational form reads: find $(\mathbf{u}, p) \in V \times Q$ such that
$$ \begin{equation} a((\mathbf{u}, p), (\mathbf{v}, q)) = L((\mathbf{v}, q)) \quad \forall (\mathbf{v}, q) \in V \times Q \end{equation} $$where bilinear and linear forms are, respectively
$$ \begin{equation} a((\mathbf{u}, p), (\mathbf{v}, q)) = \int_{\Omega} (\nabla^2 \mathbf{u} - \nabla p) \cdot {\mathbf{v}} \, dx_w + \int_{\Omega} \nabla \cdot \mathbf{u} \, {q} \, dx_w, \end{equation} $$$$ \begin{equation} L((\mathbf{v}, q)) = \int_{\Omega} \mathbf{f} \cdot {\mathbf{v}}\, dx_w + \int_{\Omega} h {q} \, dx_w \end{equation} $$Using integration by parts for Legendre
$$ \begin{equation} a((\mathbf{u}, p), (\mathbf{v}, q)) = -\int_{\Omega} \nabla \mathbf{u} \cdot \nabla{\mathbf{v}} \, dx_w + \int_{\Omega} \nabla \cdot \mathbf{v} \, {p} \, dx_w + \int_{\Omega} \nabla \cdot \mathbf{u} \, {q} \, dx_w, \end{equation} $$
In [25]:
N = (40, 40)
family = 'Legendre'
D0X = FunctionSpace(N[0], 'Legendre', bc=(0, 0))
#D1Y = FunctionSpace(N[1], 'Legendre', bc=(1, 0)) # Regular lid
D1Y = FunctionSpace(N[1], 'Legendre', bc=(0, (1-x)**2*(1+x)**2)) # Regularized lid
D0Y = FunctionSpace(N[1], 'Legendre', bc=(0, 0))
PX = FunctionSpace(N[0], 'Legendre')
PY = FunctionSpace(N[1], 'Legendre')
PX.slice = lambda: slice(0, PX.N-2)
PY.slice = lambda: slice(0, PY.N-2)
# All required spaces
V1 = TensorProductSpace(comm, (D0X, D1Y))
V0 = TensorProductSpace(comm, (D0X, D0Y))
Q = TensorProductSpace(comm, (PX, PY))
V = MixedTensorProductSpace([V1, V0])
W = MixedTensorProductSpace([V0, V0])
VQ = MixedTensorProductSpace([V, Q])
# All required test and trial functions
up = TrialFunction(VQ)
vq = TestFunction(VQ)
u, p = up
v, q = vq
In [26]:
# Assemble matrices
A = inner(grad(v), -grad(u))
G = inner(div(v), p)
D = inner(q, div(u))
# Extract the boundary matrices
bc_mats = extract_bc_matrices([A, G, D])
# Create Block matrix
M = BlockMatrix(A+G+D)
BM = BlockMatrix(bc_mats)
# Add contribution to rhs from inhomogeneous basis
up_hat = Function(VQ).set_boundary_dofs()
fh_hat = Function(VQ)
fh_hat = BM.matvec(-up_hat, fh_hat)
# Solve Stokes problem. Note constraint for pressure
up_hat = M.solve(fh_hat, u=up_hat, constraints=((2, 0, 0),))
# Move solution to Array in real space
up = up_hat.backward()
u_, p_ = up
In [27]:
X = Q.local_mesh(True)
plt.quiver(X[0], X[1], u_[0], u_[1])
Out[27]:
In [28]:
%matplotlib notebook
plt.spy(M.diags(), markersize=0.5)
Out[28]:
where $D = G^T$ for the Legendre basis, making $M$ symmetric. For Chebyshev $M$ will not be symmetric.
Solver through scipy.sparse.linalg
For Navier-Stokes of the lid-driven cavity, see https://github.com/spectralDNS/shenfun/blob/master/demo/NavierStokesDrivenCavity.py