Demo - Mixed bases for the Helmholtz problem#
Mikael Mortensen (email: mikaem@math.uio.no), Department of Mathematics, University of Oslo.
Date: March 22, 2021
Summary. This demo shows how to solve the Helmholtz equation using different bases for test and trial spaces. The use of different bases leads for some optimal combinations to highly sparse and well-conditioned coefficient matrices.
The Helmholtz problem#
We will consider Helmholtz equation with homogeneous Dirichlet boundary conditions
where \(\alpha \in \mathbb{R^+}\). The relevant function space for the Dirichlet problem is
and the Chebyshev-Galerkin (CG) method is to find \(u_N \in V_N\) such that
where \((u,v)_{\omega^{\sigma}}=\int_{{I}}uv\omega^{\sigma} dx\) is the scalar product in the weighted space \(L^2_{\omega^{\sigma}}({I})\).
Shenfun has implemented three different Chebyshev Dirichlet basis functions
These three bases are all linearly dependent and they are all bases for \(V_N\).
Implementation#
We can get all three function spaces as
from shenfun import *
N = 40
V0 = FunctionSpace(N, 'C', basis='ShenDirichlet')
V1 = FunctionSpace(N, 'C', basis='Heinrichs')
V2 = FunctionSpace(N, 'U', basis='CompactDirichlet')
where \(V0 = \text{span}\{\phi_k\}_{k=0}^{N-3}\),
\(V1 = \text{span}\{\varphi_k\}_{k=0}^{N-3}\) and
\(V2 = \text{span}\{\psi_k\}_{k=0}^{N-3}\). Now, to solve the Helmholtz problem we simply need to choose
test and trial bases. Shen’s original method is using
V0 for both. To assemble the stiffness and mass matrices
for this choice do
u = TrialFunction(V0)
v = TestFunction(V0)
A = inner(v, div(grad(u)))
B = inner(v, u)
A manufactured solution can be chosen using Sympy We choose
implemented as
import sympy as sp
x = sp.Symbol('x', real=True)
ue = sp.sin(2*sp.pi*sp.cos(2*sp.pi*x))
The right hand side \(f\) of Helmholtz equation is
alpha = 1
f = sp.simplify(alpha*ue-ue.diff(x, 2))
To solve the problem we can do
fj = Array(V0, buffer=f) # Get f on quadrature mesh
f_hat = inner(v, fj) # Compute right hand side
M = alpha*B - A # Get coefficient matrix
u_hat = Function(V0) # Container for the solution
sol = la.Solver(M) # Solver
u_hat = sol(f_hat, u_hat) # Solve
Compare with exact solution.
uj = Array(V0, buffer=ue)
error = inner(1, (u_hat.backward()-uj)**2)
print('Error =', error)
Now that was the solution for test and trial bases from the same basis (4). Let us create a function that takes any test and any trial basis, any manufactured solution and any \(\alpha\) in the Helmholtz equation. We let the function return either the L2-error norm, the condition number of the Helmholtz coefficient matrix, or the matrix itself.
def main(N, test, trial, alpha=1, method=0, ue=sp.sin(2*sp.pi*sp.cos(2*sp.pi*x))):
"""Solve Helmholtz problem and return L2-error, condition number or matrix
Parameters
----------
N : int
Number of quadrature points
test, trial : int
Test and trial functions.
0 = :math:`T_k-T_{k+2}`
1 = :math:`(1-x^2)T_k`
2 = :math:`U_k-\frac{k+1}{k+3}U_{k+2}`
alpha : Helmholtz parameter
method : int
0 = Return L2-error norm
1 = Return condition number of matrix
2 = Return Helmholtz matrix
ue : Sympy function, optional
The manufactured solution with homogeneous boundary conditions.
Note
----
Inhomogeneous boundary conditions require a small rewrite, but is
not difficult.
"""
bases = {0: ('C', 'ShenDirichlet'), 1: ('C', 'Heinrichs'), 2: ('U', 'CompactDirichlet')}
test = FunctionSpace(N, bases[test][0], basis=bases[test][1])
trial= FunctionSpace(N, bases[trial][0], basis=bases[trial][1])
# Check that boundary conditions are homogeneous
assert abs(ue.subs(x, -1)) < 1e-8 and abs(ue.subs(x, 1)) < 1e-8
u = TrialFunction(trial)
v = TestFunction(test)
f = sp.simplify(alpha*ue-ue.diff(x, 2))
fj = Array(test, buffer=f) # Get f on quadrature mesh
f_hat = inner(v, fj) # Compute right hand side
B = inner(v, u)
A = inner(v, div(grad(u)))
M = alpha*B-A
if method == 1:
return np.linalg.cond(M.diags('csr').toarray())
if method == 2:
return M
u_hat = Function(trial)
sol = la.Solver(M)
u_hat = sol(f_hat, u_hat)
uj = Array(trial, buffer=ue)
error = np.sqrt(inner(1, (u_hat.backward()-uj)**2))
return error
Let us first try basis (4) as test function and (5) as trial function. Use otherwise default parameters.
error = main(100, 0, 1)
print(error)
So the error is small in deed. Perhaps more interesting, let’s look at the sparsity pattern of the coefficient matrix
M = main(100, 0, 1, method=2)
import plotly.express as px
z = np.where(abs(M.diags().toarray()) > 1e-6, 0, 1).astype(bool)
fig = px.imshow(z, binary_string=True)
fig.show()
#plt.spy(M.diags(), markersize=0.2) # or use matplotlib
The coefficient matrix has 4 non-zero diagonals. You can now experiment with different test and trial functions, but you will not get a better result than that. Try basis (5) for both test and trial function, and you’ll get 5 nonzero diagonals.
To see the convergence rate call main for a range of
different mesh sizes
error = []
N = (2**4, 2**6, 2**8, 2**10)
for n in N:
error.append(main(n, 0, 1))
fig = px.line(x=N, y=error, log_y=True)
fig.update_layout(yaxis=dict(showexponent='all', exponentformat='e'))
fig.show()