HPC Project¶
General Remarks: This project consists of two stages. The deadline for the first stage is Friday 8th March. The deadline for the second stage is Friday 29 March. There will be no extension apart from approved extenuating circumstances. Each part of the project counts as 30% of the total module mark. The module can only be passed if at least 50% of the possible project marks have been achieved.
Stage 1¶
We are considering the poisson equation $$-\nabla \cdot \left(\sigma(x, y)\nabla\right)u(x, y) = f(x, y)$$ for $(x, y)\in [0, 1]\times [0, 1]$ with boundary conditions $u(x, y) = 0$.
The goal of Stage 1 is to gain experience with iterative solvers. The broad strategy is as follows:
1. Implement a finite difference discretisation of the above equation in OpenCL
2. Demonstrate that your implementation works by comparing with FEniCS.
3. Test the performance of iterative solvers as you refine the discretisation.
For the first point you will use a standard centered finite difference scheme. Let $u_{i, j} := u(x_i, y_j)$ be the solution on given grid point. We can approximate the application of the left-hand side operator $\nabla \cdot \left(\sigma(x, y)\nabla\right)$ through $$ \nabla \cdot \left(\sigma(x, y)\nabla\right)u \approx \frac{\left(\sigma_{i+1/2, j}\frac{\left(u_{i+1,j} – u_{i,j}\right)}{h}\right) – \left(\sigma_{i-1/2, j}\frac{\left(u_{i,j} – u_{i-1,j}\right)}{h}\right)}{h} + \frac{\left(\sigma_{i, j+1/2}\frac{\left(u_{i,j+1} – u_{i,j}\right)}{h}\right) – \left(\sigma_{i, j-1/2}\frac{\left(u_{i,j} – u_{i,j-1}\right)}{h}\right)}{h} $$ (be careful with the additional $-$ sign in the PDE).
Here, we approximate $\sigma_{i+1/2,j}\approx \frac{1}{2}\left(\sigma_{i+1,j} + \sigma_{i, j}\right)$ and $\sigma_{i,j+1/2}\approx \frac{1}{2}\left(\sigma_{i,j+1} + \sigma_{i, j}\right)$ and similarly with the other values of $\sigma$.
You should implement this iterative scheme in a matrix free way, meaning that your OpenCL implementation should involve a function that takes your vector of $u$ values and returns the discrete application of the above operator onto this vector of values. Turn this function into a linear operator class using the LinearOperator object from scipy. This allows you to apply the various iterative solvers from scipy.
To validate your implementation you should compare your own solution for a given set of right-hand side data with a FEniCS implementation of this problem. Validation is important and it is good practice in Scientific Computing to always validate model problems against existing codes.
Once you know that your implementation works experiment with different iterative solvers. In particular try gmres and bicgstab. You can use the Scipy implementations. Demonstrate the convergence by plotting the residual against the iteration count (using a log10 scale on the y axis (see semilogy plots in Matplotlib). How does the convergence change if you increase the number of discretisation points in each space dimension? Try to experimentally obtain a relationship between number of discretization points and the rate of convergence.
To make matters interesting, for your convergence plots consider a random distribution for $\sigma$, that is use $\sigma(x, y) = e^{-S(x, y)}$, where $S(x, y)$ is field of normally distributed random numbers. Your right-hand side does not matter much for convergence plots and you can use for example the function $f(x) = 1$ here.
Stage 2¶
For Stage 2 we now consider the parabolic equation $$ u_t = \nabla \cdot \left(\sigma(x, y)\nabla\right)u $$
with the same space domain as in Stage 1 and time $t$ starting at $0$. Initial conditions are given as $u(x, y, 0) = g(x, y)$ at time $t=0$ and $u(x, y, t) = 0$ on the boundary for all $t\geq 0$. For stage 2 you are asked to develop an OpenCL implementation of this parabolic PDE. For the time discretisation you can use a simple forward difference scheme given as $$ u_t(x, y)\approx \frac{u(x, y, t+\Delta t) – u(x, y, t)}{\Delta t}. $$
Demonstrate the time-domain solver for some interesting choices of initial conditions and a given random field for $\sigma$ as in Stage 1. This is a realistic simulation of a time-dependent diffusion process in two dimensions. a good initial condition is to put some weight into the middle of the domain and observe how it slowly diffuses in time.
For Stage 2 performance is relevant. Your solution must contain a directly runnable command that measures the time of performing 10 forward time steps and returns the results in seconds. For this write a Python function that takes the current iterate and performs one forward time step with your OpenCL implementation. The function should then be called 10 times in your performance measurements. I need to be able to directly run your benchmarking code.
Submission
For each Stage your submission must consist of two parts, a Jupyter Notebook that demonstrates all the codes and a pdf document that explains the implementation, demonstrates experiments, gives performance results, etc. For each Stage your pdf document must be no more than 5 pages in length. The pdf for Stage 2 can be a continuation of the Stage 1 document. You need not repeat the introduction to your code.
The Jupyter notebook for each Stage must be executable by pressing “Restart and Run All” from a fresh Notebook. Syntax errors, missing variables or other things that make the Notebook not execute are automatically marked with 0 marks for the code.
The codes will be tested in a Python 3.6 environment. Each Notebook should run in a reasonably short amount of time (around 2 minutes max). If it runs too long I will interrupt it and not mark it. Hence, not every configuration from your report needs to be in the demonstration code. But it should demonstrate all experiments that you perform, e.g. when running convergence tests only a smaller faster running configuration needs to be in the Jupyter notebook while your report also shows larger sizes.
You have almost three weeks for Stage 1 and another 2 weeks for Stage 2. Use the time and do not start too late!!
Marking: For Stage 1 marks will be as follows: Quality of the report: 50%, Code quality: 50%. For Stage 2 it will be: Report: 40%, Code quality: 50%, Performance: 10%.
In [ ]: