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

\[ \begin{equation} \alpha u - u^{''} = f \quad \text{in} \, {I}=(-1, 1), \quad u(\pm 1) = 0, \label{_auto1} \tag{1} \end{equation} \]

where \(\alpha \in \mathbb{R^+}\). The relevant function space for the Dirichlet problem is

\[ \begin{equation} S_N=\text{span}\{T_k\}_{k=0}^{N-1}, \quad V_{N} = \{v \in {S}_N\,|\, v(\pm 1) = 0\}, \label{_auto2} \tag{2} \end{equation} \]

and the Chebyshev-Galerkin (CG) method is to find \(u_N \in V_N\) such that

\[ \begin{equation} \label{eq:dirGalerkin} \tag{3} \alpha (u_N, v)_{\omega^{\sigma}} -(u^{''}_N, v)_{\omega^{\sigma}} = (f, v)_{\omega^{\sigma}}, \forall \, v \in V_N, \end{equation} \]

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

\[ \label{eq:shen} \tag{4} \phi_k = T_k-T_{k+2}, \quad k=0,1, \ldots, N-3, \]

\[ \label{eq:heinrichs} \tag{5} \varphi_k = (1-x^2)T_k, \quad k=0,1, \ldots, N-3, \]

\[ \label{eq:dirichletU} \tag{6} \psi_k = U_k-\frac{k+1}{k+3}U_{k+2}, \quad k=0,1, \ldots, N-3. \]

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

\[ \begin{equation} u(x) = \sin \left( 2 \pi \cos \left( 2 \pi x \right) \right) \label{_auto3} \tag{7} \end{equation} \]

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()