Skip to main content
SciNet
  • Home
  • All Courses
  • Calendar
  • Certificates
  • SciNet
    Main Site Documentation my.SciNet
  • CCDB
  • More
Close
Toggle search input
English
English Français
You are currently using guest access
Log in
SciNet
Home All Courses Calendar Certificates SciNet Collapse Expand
Main Site Documentation my.SciNet
CCDB
Expand all Collapse all
  1. Dashboard
  2. PHY1610 - Winter 2023
  3. Assignment 6: Using BLAS for the 1d wave

Assignment 6: Using BLAS for the 1d wave

Completion requirements
Opened: Wednesday, 8 March 2023, 12:00 AM
Due: Wednesday, 15 March 2023, 11:59 PM

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.


  • a6_wave1d_blas.zip a6_wave1d_blas.zip
    8 March 2023, 11:35 AM
Contact site support
You are currently using guest access (Log in)
Data retention summary


All content on this website is made available under the Creative Commons Attribution 4.0 International licence, with the exception of all videos which are released under the Creative Commons Attribution-NoDerivatives 4.0 International licence.
Powered by Moodle