-
Notifications
You must be signed in to change notification settings - Fork 42
/
biharmonic2D.py
80 lines (61 loc) · 2.41 KB
/
biharmonic2D.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
r"""
Solve Biharmonic equation in 2D with periodic bcs in one direction
and homogeneous Dirichlet and Neumann in the other
\nabla^4 u = f,
Use Fourier basis for the periodic direction and Shen's Biharmonic
basis for the non-periodic direction.
"""
import os
import sympy as sp
import numpy as np
from shenfun import inner, div, grad, TestFunction, TrialFunction, Array, \
Function, TensorProductSpace, FunctionSpace, comm, la, chebyshev
# Use sympy to compute a rhs, given an analytical solution
x, y = sp.symbols("x,y", real=True)
ue = x**4*sp.sin(2*y)
fe = ue.diff(x, 4) + ue.diff(y, 4) + 2*ue.diff(x, 2, y, 2)
def main(N, family):
bcs = {'left': {'D': ue.subs(x, -1), 'N': ue.diff(x, 1).subs(x, -1)},
'right': {'D': ue.subs(x, 1), 'N': ue.diff(x, 1).subs(x, 1)}}
SD = FunctionSpace(N, family=family, bc=bcs)
K1 = FunctionSpace(N, family='F')
T = TensorProductSpace(comm, (SD, K1), axes=(0, 1))
u = TrialFunction(T)
v = TestFunction(T)
# Get f on quad points
fj = Array(T, buffer=fe)
# Compute right hand side of biharmonic equation
f_hat = inner(v, fj)
# Get left hand side of biharmonic equation
matrices = inner(v, div(grad(div(grad(u)))))
u_hat = Function(T) # Solution spectral space
# Create linear algebra solver
BiharmonicSolver = chebyshev.la.Biharmonic if SD.family() == 'chebyshev' and N % 2 == 0 else la.SolverGeneric1ND
H = BiharmonicSolver(matrices)
# Solve and transform to real space
u_hat = H(f_hat, u_hat)
uq = u_hat.backward()
# Compare with analytical solution
uj = Array(T, buffer=ue)
error = np.sqrt(inner(1, (uj-uq)**2))
if comm.Get_rank() == 0:
print(f"biharmonic2D {SD.family():14s} L2 error = {error:2.6e}")
if 'pytest' not in os.environ:
import matplotlib.pyplot as plt
plt.figure()
X = T.local_mesh(True) # With broadcasting=True the shape of X is local_shape, even though the number of datapoints are still the same as in 1D
plt.contourf(X[0], X[1], uq)
plt.colorbar()
plt.figure()
plt.contourf(X[0], X[1], uj)
plt.colorbar()
plt.figure()
plt.contourf(X[0], X[1], uq-uj)
plt.colorbar()
plt.title('Error')
plt.show()
else:
assert error < 1e-6
if __name__ == '__main__':
for family in ('legendre', 'chebyshev', 'chebyshevu', 'ultraspherical', 'jacobi'):
main(24, family)