Assignment 6: Using BLAS for the 1d wave

In this assignment, you will be taking the one-dimensional wave equation code and re-implementing the one_time_step function using linear algebra calls to BLAS (we could take the two-dimensional case too, but it would be a bit more involved).

Background

First, let us review the derivation of the current numerical procedure to solve the equation

$$\frac{\partial^2 \rho}{\partial t^2}= c^2\frac{\partial^2 \rho}{\partial x^2}-\frac{1}{\tau}\frac{\partial \rho}{\partial t}.$$

We discretize space as \(x_j=j\Delta x\) and time as \(t_s=s\Delta t\), denote \(\rho_{s,j}=\rho(s\Delta t, j\Delta x)\), and replace the derivatives with discretized forms:

$$\frac{\partial^2 \rho}{\partial t^2}\to\frac{\rho_{s+1,j}+\rho_{s-1,j}-2\rho_{sj}}{\Delta t^2}$$
$$\frac{\partial^2 \rho}{\partial x^2} \to\frac{\rho_{s,j+1}+\rho_{s,j-1}-2\rho_{sj}}{\Delta x^2}$$
$$\frac{\partial \rho}{\partial t} \to\frac{\rho_{s,j}-\rho_{s-1,j}}{\Delta t}$$

so we get:

$$\frac{\rho_{s+1,j}+\rho_{s-1,j}-2\rho_{sj}}{\Delta t^2}= c^2\frac{\rho_{s,j+1}+\rho_{s,j-1}-2\rho_{sj}}{\Delta x^2}- \frac{\rho_{s,j}-\rho_{s-1,j}}{\tau\Delta t}.$$

Since we are solving an initial value problem, \(\rho_s\) and \(\rho_{s-1}\) will be known, and the above equation needs to be solved for \(\rho_{s+1}\), giving

$$\rho_{s+1,j} = 2\rho_{sj} - \rho_{s-1,j} +  \Delta t^2\underbrace{\frac{c^2}{\Delta x^2} (\rho_{s,j+1}+\rho_{s,j-1}-2\rho_{sj})}_{\mbox{laplacian}}- \Delta t\underbrace{\frac{\rho_{s,j}-\rho_{s-1,j}}{\tau}}_{\mbox{friction}}.$$

In this one-timestep rule, we have indicated the variable names used in the code for parts of the right hand side. Note that in the code, only 3 time steps are kept at any moment, such that \(\rho_s\) is called rho, \(\rho_{s-1}\) is called rho_prev and \(\rho_{s+1}\) is called rho_next.

Linear Algebra Formulation

To write this as a linear algebra equation, we take together the terms involving \(\rho_{s,j}\) and \(\rho_{s-1,j}\):

$$\rho_{s+1,j} = \left(2 -2\frac{\Delta t^2c^2}{\Delta x^2}- \frac{\Delta t}{\tau}\right)\rho_{sj}  +\frac{\Delta t^2c^2}{\Delta x^2} (\rho_{s,j+1}+\rho_{s,j-1}) + \left(\frac{\Delta t}{\tau}-1\right) \rho_{s-1,j}$$

which is of the form

$$\rho_{s+1}=\mathbf{A}\rho_s+\beta \rho_{s-1}.$$

where \(\rho_{s-1},\rho_{s}\) and \(\rho_{s+1}\) are seen as vectors, and \(\mathbf{A}\) and \(\beta\) are a tri-diagonal matrix and a scalar, respectively. A tri-diagonal matrix is a special case of a banded matrix that has non-zero elements only on the diagonal and one adjacent off-diagonal above and below the diagonal.

Your assignment

• Install OpenBLAS on your account on the Teach cluster (see slides of lecture 14).
• Get the code, which is the initial wave1d code from Assignment 4, except it uses rarrays and explicitly enforces boundary conditions:

    git clone /scinet/course/phy1610/a6_wave1d_blas

• Add an integrated test to the Makefile.
• Rewrite the one_time_step function so it calls a cblas matrix-vector multiplication routine from the openblas library.
  For this, you will need to select the appropriate BLAS routine to call, look up how the matrix \(\mathbf{A}\) should be stored for that routine, and determine the precise values of \(\beta\) and the matrix elements of \(\mathbf{A}\).
• Explain your blas routine selection in comments in the code.
• Modify the Makefile to link to openblas.
• Run the integrated test to make sure it still gives the right answer.

For documentation on (C)BLAS, check the links at the bottom of the course site.

As in previous assignments, we expect a number of meaningful commits and comments in the code.  Submit the git repo as a zip file (using git2zip) by March 15, 2023, 11:59 PM.