1 MATH4063
The University of Nottingham
SCHOOL OF MATHEMATICAL SCIENCES
AUTUMN SEMESTER 2022-2023
MATH4063 - SCIENTIFIC COMPUTING AND C++
Coursework 1 - Released 30th October 2023, 4pm
Your work should be submitted electronically via the MATH4063 Moodle page by 12noon on Monday 20th
November (unless you have arranged an extension). Since this work is assessed, your submission must be
entirely your own work (see the University’s policy on Academic Misconduct). Submissions up to five working
days late will be marked, but subject to a penalty of 5% of the maximum mark per working day.
The marks for each question are given by means of a figure enclosed by square brackets, eg [20]. There are
a total of 100 marks available for the coursework and it contributes 45% to the module. The marking rubric
available on Moodle will be applied to each full question to further break down this mark.
You are free to name the functions you write as you wish, but bear in mind these names should be meaningful.
Functions should be grouped together in .cpp files and accessed in other files using correspondingly named
.hpp files.
All calculations should be done in double precision.
A single zip file containing your full solution should be submitted on Moodle. This zip file should contain three
folders called main, source and include, with the following files in them:
main:
• q1d.cpp
• q2c.cpp
• q3c.cpp
• q4b.cpp
source:
• vector.cpp
• dense_matrix.cpp
• csr_matrix.cpp
• linear_algebra.cpp
• finite_volume.cpp
include:
• vector.hpp
• dense_matrix.hpp
• csr_matrix.hpp
• linear_algebra.hpp
• finite_volume.hpp
Prior to starting the coursework, please download the CW1_code.zip from Moodle and extract the files. More
information about the contents of the files included in this zip file is given in the questions below.
Hint: When using a C++ struct with header files, the whole struct needs to be defined fully in the header file,
and the header file included in the corresponding .cpp file. Include guards should also be used.
MATH4063 Turn Over
2 MATH4063
In this coursework you will build a 2D finite volume solver for the following PDE boundary value problem
−𝛥𝑢 + ∇ ⋅ (b𝑢) = 𝑓 (𝑥, 𝑦) ∈ 𝛺, (1)
𝑢 = 𝑔, (𝑥, 𝑦) ∈ 𝜕𝛺, (2)
where 𝑓 ∶ 𝛺 → ℝ, 𝑔 ∶ 𝜕𝛺 → ℝ and b ∶ 𝛺 → ℝ2
.
In order to solve this problem, you will first define a sparse matrix structure, then write functions to apply
the GMRES linear algebra solver and finally build and solve the linear system arising from the finite volume
approximation of (1)-(2).
1. Matrices arising from the discretisation of partial differential equations using, for example, finite volume
methods, are generally sparse in the sense that they have many more zero entries than nonzero ones.
We would like to avoid storing the zero entries and only store the nonzero ones.
A commonly employed sparse matrix storage format is the Compressed Sparse Row (CSR) format. Here,
the nonzero entries of an 𝑛 × 𝑛 matrix are stored in a vector matrix_entries, the vector column_no gives
the column position of the corresponding entries in matrix_entries, while the vector row_start of length
𝑛+1 is the list of indices which indicates where each row starts in matrix_entries. For example, consider
the following:
𝐴 =
⎛
⎜
⎜
⎝
8 0 0 2
0 3 1 0
0 0 4 0
6 0 0 7
⎞
⎟
⎟
⎠
⟶
matrix_entries = (8 2 3 1 4 6 7)
column_no = (0 3 1 2 2 0 3)
row_start = (0 2 4 5 7)
Note, in the above, C++ indexing has been assumed, i.e, indices begin at 0.
(a) In csr_matrix.hpp, define a C++ struct called csr_matrix to store a matrix in CSR format. In
addition to matrix_entries, column_no and row_start, you should store the number of rows of the
matrix explicitly.
(b) In csr_matrix.cpp, write a C++ function that will set up the matrix 𝐴 from above in CSR format.
Remember, if you are using dynamically allocated memory, then you should also have corresponding
functions that will deallocate the memory you have set up.
(c) In csr_matrix.cpp, write a C++ function that takes as input a matrix 𝐴 stored in CSR format and a
vector x and computes the product 𝐴x. The prototype for your function should be:
void MultiplyMatrixVector ( csr_matrix & matrix ,double* vector ,
double* productVector )
Hence, the input vector and the output productVector should be pointers to dynamically allocated
arrays. In particular, it should be assumed that productVector has been preallocated to the correct
size already.
(d) By setting a vector x = (4, −1, 3, 6)⊤, write a test program in q1d.cpp to compute and print to the
screen the product 𝐴x, where 𝐴 is the matrix given above.
[20 marks]
MATH4063
3 MATH4063
2. Suppose we wish to find x ∈ ℝ𝑛
such that
𝐴x = b, (3)
where 𝐴 is an 𝑛 × 𝑛 matrix and b ∈ ℝ𝑛
.
One algorithm for solving this problem is the (restarted) Generalised Minimal RESidual (GMRES) algorithm.
The method is too complicated to explain here, but works to quickly find approximations x𝑘 = x0 + y𝑘
where y𝑘 ∈ 𝒦𝑘 ∶= Span{𝐴q0
, 𝐴2q0 … 𝐴𝑘q0
} for 𝑘 = 1, 2, …. y𝑘 is chosen to minimise the residual
‖b − 𝐴x𝑘‖2
.
Here x0
is some initial guess vector and q0
is the normed initial residual
q0 =
b − 𝐴x0
‖b − 𝐴x0‖2
.
𝒦𝑘 is called a Krylov subspace of 𝐴.
The algorithm stops when ‖b − 𝐴x𝑘‖2 < tol for some termination tolerance tol. As the method becomes
very memory inefficient when 𝑘 is large, the method is restarted every so often and x𝑘 reset to be x0
.
An incomplete GMRES algorithm function PerformGMRESRestarted() has been written in
linear_algebra.cpp.
A key component of the GMRES algorithm is the Arnoldi iteration that seeks to find an orthonormal basis
of 𝒦𝑘. At the 𝑘th step of the iteration, the Arnoldi method constructs the following matrix decomposition
of 𝐴:
𝐴𝑄𝑘 = 𝑄𝑘+1𝐻̃
𝑘,
where the columns of 𝑄𝑘 (𝑄𝑘+1) contain the orthonormal basis of 𝒦𝑘 (𝒦𝑘+1, resp.) and 𝐻̃
𝑘 is a (𝑘+1)× 𝑘
upper Hessenberg matrix. That is, a matrix that is nearly upper triangular but has non-zero components
on the first subdiagonal.
The 𝑘th step of the Arnoldi algorithm is:
Algorithm 1 One step of the Arnoldi Iteration.
Require: 𝑘 > 0, 𝐴, 𝑄𝑘:
1: Let q𝑖 be the 𝑖th column of 𝑄𝑘.
2: Let h = {ℎ𝑖
}
𝑘+1
𝑖=1 be a vector of length 𝑘 + 1.
3: Compute q𝑘+1 = 𝐴q𝑘
4: for 𝑖 = 1, … , 𝑘 do
5: ℎ𝑖 = q𝑘+1 ⋅ q𝑖
.
6: q𝑘+1 = q𝑘+1 − ℎ𝑖q𝑖
.
7: end for
8: ℎ𝑘+1 = ‖q𝑘+1‖2
.
9: q𝑘+1 = q𝑘+1/ℎ𝑘.
10: 𝑄𝑘+1 = [𝑄𝑘, q𝑘+1].
11: return 𝑄𝑘+1 and h.
(a) In linear_algebra.cpp, write a C++ function which implements one step of the Arnoldi iteration
method defined above.
The function should have the following prototype
void PerformArnoldiIteration ( csr_matrix & matrix ,
dense_matrix & krylov_matrix , int k, double* hessenberg )
MATH4063 Turn Over
4 MATH4063
Here, matrix is 𝐴, k is the step of the iteration to perform, krylov_matrix is the matrix containing
the orthonormal basis, where each row is a basis vector. Upon entry, krylov_matrix should have 𝑘
rows and upon exit it should contain 𝑘 + 1 rows, with the new basis vector in the last row.
Finally, upon exit, hessenberg should contain h, which is the final column of 𝐻̃
𝑘. You may assume that
hessenberg has been preallocated to be of length 𝑘+1 before the call to PerformArnoldiIteration.
Your function should make use, where possible, of prewritten functions defined in dense_matrix.cpp
and vector.cpp. Your code should also make use of the matrix multiplication function from Q1.
Once you have written PerformArnoldiIteration() the GMRES function should function as intended.
Note: Storage of the basis functions in the rows of krylov_matrix, rather than in the columns,
improves efficiency of the code.
(b) In csr_matrix.cpp, write a C++ function that will read from a file a matrix already stored in CSR
format and a vector. You may assume the file structures are as in matrix1.dat and vector1.dat on
Moodle and you may use these data files to test your function.
(c) Write a test program in file q2c.cpp that will read in the matrix 𝐴 from matrix2.dat and the vector
x from vector2.dat, compute b = 𝐴x, then use PerformGMRESRestarted() with the default input
arguments to find an approximation x̂to x. At the end of the calculation, print to the screen the error
‖x − ̂ x‖2
.
[30 marks]
3. The file mesh.hpp contains a struct that defines a mesh data structure mesh for a general mesh comprising
axis-aligned rectangular cells. In particular, each cell in the mesh has an additional struct called
cell_information that contains, among other things, information about the cell neighbours. Familiarise
yourself with these data structures by looking in mesh.hpp.
mesh.cpp contains two functions that will generate meshes, they are:
• ConstructRectangularMesh() - this constructs a mesh on the rectangular domain 𝛺𝑅 = [𝑎, 𝑏] ×
[𝑐, 𝑑].
• ConstructLShapedMesh() - this constructs a mesh on the L-shaped domain 𝛺𝐿 = 𝛺𝑅\𝛺𝐶, where
𝛺𝐶 = [(𝑎 + 𝑏)/2, 𝑏] × [(𝑐 + 𝑑)/2, 𝑑].
(a) In finite_volume.cpp, write a C++ function that will create the storage for a matrix 𝐴 in CSR format
and a RHS vector F required for a cell-centred finite volume method for solving (1)-(2). You should
follow the procedure outlined in the Unit 6 lecture notes. As one of the inputs, your function should
take in a variable of type mesh.
(b) In csr_matrix.cpp, write a C++ function that will output to the screen a matrix stored in CSR format
in the same style as in matrix1.dat.
(c) In Q3c.cpp, write a program that will ask the user to supply the number of cells in each coordinate
direction of a rectangular mesh, sets up the mesh using ConstructRectangularMesh() then calls the
function from part (a) to set up the corresponding matrix and finally prints it to the screen using the
function from part (b).
[30 marks]
MATH4063
5 MATH4063
4. (a) In finite_volume.cpp, write a function that takes in a mesh, uses the function from Q3(a) to construct
𝐴 and F, then populates it with the correct entries to solve problem (1)-(2) using the cell-centred finite
volume method, as outlined in the Unit 6 notes. The function should also take as input the functions
𝑓(𝑥, 𝑦), b(𝑥, 𝑦) and the Dirichlet boundary function 𝑔(𝑥, 𝑦).
(b) In Q4b.cpp, write a main program to ask the user to select from the following problems and supply
the number of cells in each coordinate direction.
1. • Rectangular Mesh - 𝑎 = 0, 𝑏 = 1, 𝑐 = 0 and 𝑑 = 1;
• 𝑓(𝑥, 𝑦) = 1;
• 𝑔(𝑥, 𝑦) = 0;
• b = 0.
2. • L-shaped Mesh - 𝑎 = 0, 𝑏 = 1, 𝑐 = 0 and 𝑑 = 1;
• 𝑓(𝑥, 𝑦) = 8𝜋2
cos(2𝜋𝑥) cos(2𝜋𝑦);
• 𝑔(𝑥, 𝑦) = cos(2𝜋𝑥) cos(2𝜋𝑦);
• b = 0.
3. • Rectangular Mesh - 𝑎 = −1, 𝑏 = 1, 𝑐 = −1 and 𝑑 = 1;
• 𝑓(𝑥, 𝑦) = 1;
• 𝑔(𝑥, 𝑦) = 0;
• b = (10, 10)⊤.
4. • L-Shaped Mesh - 𝑎 = 0, 𝑏 = 1, 𝑐 = 0 and 𝑑 = 1;
• 𝑓(𝑥, 𝑦) = 0;
•
𝑔(𝑥, 𝑦) = {
1, 𝑥 = 0, 0.25 < 𝑦 < 0.75,
0, otherwise;
• b = (
50𝑦
√𝑥2+𝑦2
,
−50𝑥
√𝑥2+𝑦2
)
⊤
.
The code should then set up the linear system arising from the finite volume discretisation and solve
the system
𝐴uℎ = F
using PerformGMRESRestarted().
Finally, print to the screen the maximum value of uℎ.
Hint: Once you have computed uℎ you can output it to together with the mesh to a file using
OutputSolution() in mesh.cpp. plot_solution.py can then be used to plot the solution in Python.
Note, if you are unable to get the iterative solver from Q2 working, then you may create the finite volume
matrix 𝐴 as if it were a dense matrix (i.e store all the zero entries) and use the function
PerformGaussianElimination() from dense_matrix.cpp to solve the system of equations. This will incur
a small penalty. Note, an illustration of the use of PerformGaussianElimination() can be found in the
main program inside gaussian_elimination_test.cpp.
[20 marks]
MATH4063 End

请加QQ：99515681 或邮箱：99515681@qq.com   WX：codehelp