Contents
In this entry I benchmark and test a finite-element method (FEM) Poisson solver. The solver is implemented in Lucee/Gkeyll as the class FemPoissonStructUpdater and solves the equation
where \(s\) is the spatially dependent source term. This updater works in 1D, 2D and 3D and is agnostic of the underlying nodal basis functions, which need to be provided separately.
In this test the convergence of the 1D solver is tested with an exact solution. The source is chosen to be
where \(a=2\). With this the exact solution is
where \(c_0\) and \(c_1\) are constants of integration.
For the first test we pick the domain \(x\in [0,1]\) and Dirichlet boundary conditions \(\psi(0)=\psi(1)=0\). With these we get \(c_1=0\) and \(c_0=a/12-1/2\).
The following table shows the errors for the second-order Lobatto scheme with different cell sizes corresponding to 8, 16, 32, and 64 elements and with \(a=2\).
Grid size \(\Delta x\) | Average error | Order | Simulation |
---|---|---|---|
\(0.125\) | \(3.797 \times 10^{-4}\) | s77 | |
\(0.0625\) | \(1.01725 \times 10^{-4}\) | 1.90 | s78 |
\(0.03125\) | \(2.627\times 10^{-5}\) | 1.95 | s79 |
\(0.015625\) | \(6.675726\times 10^{-6}\) | 1.98 | s80 |
An example solution with 16 elements is shown below.
For the second test we use Dirichlet and Neumann boundary conditions \({\partial \psi}/{\partial x}=0\) at \(x=0\) and \(\psi(1)=0\). With these we get \(c_0=0\) and \(c_1=a/12-1/2\).
The following table shows the errors for the second-order Lobatto scheme with different cell sizes corresponding to 8, 16, 32, and 64 elements and with \(a=5\).
Grid size \(\Delta x\) | Average error | Order | Simulation |
---|---|---|---|
\(0.125\) | \(4.20464 \times 10^{-3}\) | s81 | |
\(0.0625\) | \(1.06812 \times 10^{-3}\) | 1.98 | s82 |
\(0.03125\) | \(2.69148\times 10^{-4}\) | 1.99 | s83 |
\(0.015625\) | \(6.75519\times 10^{-5}\) | 1.99 | s84 |
An example solution with 16 elements is shown below.
In this test the convergence of the 3rd and 4th order 1D solvers is tested with an exact solution. The source is chosen to be
where \(a=2\) and \(b=-12\). With this the exact solution is
where \(c_0\) and \(c_1\) are constants of integration. For the first test we pick the domain \(x\in [0,1]\) and Dirichlet boundary conditions \(\psi(0)=\psi(1)=0\). With these we get \(c_1=0\) and \(c_0=-(1/2+a/12+b/30)\).
The following table shows the errors for the third-order Lobatto scheme with different cell sizes corresponding to 4, 8, 16, and 32 elements and with \(a=2\).
Grid size \(\Delta x\) | Average error | Order | Simulation |
---|---|---|---|
\(0.25\) | \(7.742 \times 10^{-5}\) | s89 | |
\(0.125\) | \(5.397 \times 10^{-6}\) | 3.84 | s90 |
\(0.0625\) | \(3.535 \times 10^{-7}\) | 3.93 | s91 |
\(0.03125\) | \(2.249\times 10^{-8}\) | 3.97 | s92 |
Note that the 3rd order scheme is actually converging with 4th order accuracy. In fact, the 3rd order scheme gives the exact solution when the same source as was used in testing the 2nd order scheme is used.
The following table shows the errors for the fourth-order Lobatto scheme with different cell sizes corresponding to 4, 8, 16, and 32 elements and with \(a=2\).
Grid size \(\Delta x\) | Average error | Order | Simulation |
---|---|---|---|
\(0.5\) | \(5.536 \times 10^{-5}\) | s93 | |
\(0.25\) | \(1.847 \times 10^{-6}\) | 4.90 | s94 |
\(0.125\) | \(5.994 \times 10^{-8}\) | 4.94 | s95 |
\(0.0625\) | \(1.910\times 10^{-8}\) | 4.97 | s96 |
In this test the convergence of the 2D solver is tested with an exact solution. The exact solution is chosen to be
where
Here, I have chosen \(a=2\), \(b=5\), \(c_1=d_0=0\) and \(c_0=a/12-1/2\) and \(d_1=b/12-1/2\). This corresponds to Dirichlet boundary conditions on the left, right and top edge and a Neumann boundary condition on the bottom edge.
The following table shows the errors for the second-order Lobatto scheme with different cell sizes corresponding to \(8\times 8\), \(16\times 16\), \(32\times 32\), and \(64\times 64\) element grids.
Grid size \(\Delta x\) | Average error | Order | Simulation |
---|---|---|---|
\(0.125\) | \(2.90322 \times 10^{-5}\) | s85 | |
\(0.0625\) | \(7.8699 \times 10^{-6}\) | 1.88 | s86 |
\(0.03125\) | \(2.04355\times 10^{-6}\) | 1.95 | s87 |
\(0.015625\) | \(5.5524\times 10^{-7}\) | 1.89 | s88 |
In this test the convergence of the third-order 2D solver is tested with an exact solution. The problem setup is the same as used in testing the second-order solver, except using the third-order Serendipity basis functions.
The following table shows the errors for the third-order scheme with different cell sizes corresponding to \(8\times 8\), \(16\times 16\), \(32\times 32\), and \(64\times 64\) element grids.
Grid size \(\Delta x\) | Average error | Order | Simulation |
---|---|---|---|
\(0.25\) | \(1.156 \times 10^{-5}\) | s97 | |
\(0.125\) | \(8.767 \times 10^{-7}\) | 3.72 | s98 |
\(0.0625\) | \(6.043 \times 10^{-8}\) | 3.85 | s99 |
\(0.03125\) | \(5.200\times 10^{-9}\) | 3.53 | s100 |
The solution converges with greater than third order, as it did for the 1D solver.
The finite-element solver needs significant modification when periodic boundary conditions need to be applied. First, we need to ensure that the source integrated over the domain vanishes, i.e
To ensure this condition the updater first computes the integrated source and removes that when computing the global source vector.
Further, one needs to carefully build the stiffness matrix to take into account those nodes that the identified with each other due to the periodic boundary conditions. This needs careful book-keeping of local to global index mappings and also accounting for the same periodicity in computing the global source vector. Finally, the stiffness matrix will be singular as with periodic BCs the solution is only determined to an additive constant. To avoid problems with the linear solve one needs to set one location in the domain to an arbitrary value. I have chosen to set the bottom left corner as \(\phi(0,0) = 0\). All these modifications are not trivial and need 2X to 3X more code than needed to just support Dirichlet and/or Neumann boundary conditions.
In this section I test the convergence of the solver using an exact solution given by
where \(N\) is a normalizing factor and \(a_{mn}\) and \(b_{mn}\) are specified coefficients. The source corresponding to this solution is
The domain is \(x \in [0,2\pi]\) and \(y \in [0,2\pi]\). For the other quantities used see the Lua scripts linked below.
The following table shows the errors for the second-order scheme with different cell sizes corresponding to \(32\times 32\), \(64\times 64\) and \(128\times 128\) element grids.
Grid size \(\Delta x\) | Average error | Order | Simulation |
---|---|---|---|
\(2\pi/32\) | \(3.718 \times 10^{-3}\) | s101 | |
\(2\pi/64\) | \(9.595 \times 10^{-4}\) | 1.95 | s102 |
\(2\pi/128\) | \(2.903 \times 10^{-4}\) | 1.72 | s103 |
With these tests I am reasonably confident that the Poisson finite-element updater works correctly. I have shown the expected convergence behaviour of the 1D and 2D schemes. In 2D I also tested periodic boundary conditions and showed second-order convergence. I have not yet tested the periodic BCs in 1D or the third order 2D periodic BC solver. I will update this entry when I do so.