School of Computing: assessment brief
   Module title
 Scientific Computation
  Module code
 COMP5930M
  Assignment title
 Coursework 2
  Assignment type and description
 Coursework assignment
  Rationale
 TBA
  Weighting
 20% of total mark
  Submission dead- line
 December 14th 2023 at 10:00
  Submission method
 Turnitin submission through Minerva
  Feedback provision
 Feedback provided on Minerva
  Learning outcomes assessed
 (i) Formulate and solve systems of nonlinear equations to solve challenging real-world problems arising from en- gineering and computational science; (ii) Implement al- gorithmic solutions to solve computational differential equation problems based on mathematical theory; (iii) Analyse computational linear algebra problems to iden- tify and implement the most efficient and scalable solu- tion algorithm to apply for large problems.
  Module lead
 Dr Toni Lassila
           1

1. Assignment guidance
Provide answers to the two exercises below. Answer both exercises.
2. Assessment tasks
Exercise 1: The Burgers’ equation models the propagation of a pres- sure wave in shock tube. It is a nonlinear partial-differential equation in one spatial dimension to find u(x, t) s.t.
∂u + u∂u = ν ∂2u, (1) ∂t ∂x ∂x2
where the boundary conditions u(a, t) = ua and u(b, t) = ub for all t, and the initial condition u(x, 0) = u0(x) need to be prescribed in order to obtain a well-posed problem. Here ν is the kinematic viscosity of the fluid. For ν = 0 we have the inviscid Burgers’ equation, and for ν > 0 we have the viscous Burgers’ equation.
(a) Applying the central difference formula to the second order deriva- tive in space, the upwind difference formula
􏰀Uk−Uk 􏰁
i−1
using implicit Euler’s method leads to the discrete formulation: Uk −Uk−1 􏰀Uk −Uk 􏰁 􏰀Uk −2Uk +Uk 􏰁
Fi(U)= i i +Uik i i−1 −ν i+1 i i−1 =0 ∆t h h2
(2) for i = 2,3,...,m−1 where the interval has been discretised with
m uniformly distributed nodes and a spatial grid size h. Implement the function F as a python subroutine fun burgers.py
        def fun_burgers( uk, ukp, dt, h, nu, ua, ub )
where uk is the vector Uk of size m, ukp is the previous time-step solution vector Uk−1, dt is the time-step ∆t, h is the spatial grid size parameter h, and nu is the kinematic viscosity ν. Include the boundary conditions ua and ub in the implementation. [6 marks]
2
Uik i
to the first order derivative in space, and discretising (1) in time
 h
   
(b) Derive the analytical formulas for the nonzero elements on row i of the Jacobian matrix for (2): [4 marks]
∂Fi , ∂Fi, ∂Fi . ∂Ui−1 ∂Ui ∂Ui+1
(c) Solve problem (2) numerically using your fun burgers.py and the PDE solver template solver burgers.py provided in the course- work folder. Use the viscosity value ν = 0.01, the time-step ∆t=0.01,thegridsizeh=0.01,andafinaltimeofT =1. The initial solution u(x, 0) should be taken as a unit step located at x = 0.1 (see below) and the boundary conditions as: u(0, t) = 1 and u(1, t) = 0.
   Figure 1: Initial condition u0(x) for the Burgers’ equation (1)
Plot the solution u(x, T ) at the final time step T = 1 and include it in your report. Also report the total number of Newton iterations required for the numerical solution (sum of Newton iterations over all time steps). [2 marks]
3

(d) The solution of Burgers’ equation (1) can be shown to be a (decay- ing) wavefront that travels from left to right at a constant velocity v. What is the approximate value of the numerical wavefront ve- locity vnum for ν = 0.01, ∆t = 0.01, and h = 0.01? Measure the approximate location of the wavefront using the point where the solution u(xmid) ≈ 0.5. [1 mark]
(e) Replace the discretisation of the nonlinear convection term with the downwind difference formula
􏰀Uk − Uk 􏰁
i (3)
and solve the problem with same parameters as in (c). Plot the solution u(x,T) at the final time step T = 1 and include it in your report. Also report the total number of Newton iterations required for the numerical solution (sum of Newton iterations over all time steps). What is the numerical wavefront velocity vnum in this case?
Now set ν = 0.001 and solve the problem again using the down- wind difference formula. What do you observe? Now solve the problem with ν = 0.001 using the original upwind difference for- mula and compare the results. What is the numerical wavefront velocity vnum in this case? [7 marks]
Uik i+1
h
 4

Exercise 2: Consider the anisotropic diffusion equation to find u(x, y) s.t.
􏰀 ∂2u ∂2u􏰁
− μx∂x2 +μy∂y2 =f(x,y), (x,y)∈(0,1)×(0,1), (4)
and the boundary condition u = 0 on Γ (the boundary of the unit square), where u is a scalar function that models the temperature of a heat-conducting object modelled here as a unit square and f(x,y) is a function modelling a heat source. The heat conductivity coefficients, μx > 0 and μy > 0, can have different magnitudes (anisotropy).
(a) Discretising the problem (4) using the second-order finite differ- ence formulas
∂2u ≈ ui,j−1 − 2ui,j + ui,j+1 .
Write the second-order finite difference stencil (similarly as in Tu- torial 7)
∂2u ≈ ui−1,j − 2ui,j + ui+1,j , ∂x2 h2
  ∂y2
−μx h2 −μy h2 = fi,j.
h2 􏰀ui−1,j − 2ui,j + ui+1,j 􏰁 􏰀ui,j−1 − 2ui,j + ui,j+1 􏰁
leads to the discretised form
  ?**7;
s11 s12 s13 ?**8; ?**8;
S=s s s?**8;  21 22 23?**8;
?**8; s s s?**9;
corresponding to this finite difference scheme. [4 marks] (b) Implement a python function source function.py
    def source_function( x, y, h )
that returns the right-hand side by evaluating the function:
f(x,y) :=
⭺**;1, ifx≥0.1andx≤0.3andy≥0.1andy≤0.3 0, otherwise
.
Include the source code in your answer. [3 marks] 5
31 ** 33
(5)

 Figure 2: Computational domain for problem (4) and the sub-region where the heat source is located (in red).
(c) Modify the solver from Tutorial 7 to numerically solve the diffusion problem (4) for the right-hand side (5).
Solve the linear problem AU = F using the conjugate gradient method (without preconditioning) with the diffusion coefficients μx = 1 and μy = 1, stopping tolerance tol = 10−6, and maxi- mum of 1000 CG iterations. You can use the CG implementation in scipy.sparse.linalg.cg for this problem or code your own implementation.
Plot the solution surface and include the plot in your answer. How many iterations does it take for CG to converge in this case?
[2 marks]
(d) Consider now the use of a preconditioner to accelerate the con- vergence of CG. The incomplete-LU preconditioner approximates the system matrix A ≈ LincUinc by performing Gaussian elimi- nation but setting to zero any elements that are smaller than a dropoff tolerance ε chosen by the user. You can use the imple- mentation provided in scipy.sparse.linalg.spilu to compute
6

the incomplete factors Linc and Uinc.
Write a python implementation myPCG.py of the preconditioned
conjugate gradient from Lecture 18:
            def myPCG( A, b, L, U, tol, maxit )
that solves the preconditioning step for the residual, Mzi+1 = LU zi+1 = ri+1 , using appropriate solution algorithms. Include the source code as part of your answer. [4 marks]
(e) Solve the problem (4) again using your preconditioned CG imple- mentation from (d). Use a dropout tolerance of ε = 0.1 for the incomplete LU-factorisation.
How many nonzero elements (nnz) do the factors Linc and Uinc have in this case?
How many PCG iterations does the problem take to converge to tol = 10−6 now?
[2 marks]
(f) Repeat the experiment from (e) with different values of the dif- fusion coefficients. Solve the problem (4) with μx = 0.1 and μx = 0.01, while keeping the other value at μy = 1. Solve the problem using PCG with the same ILU-preconditioner as before with a dropout tolerance of ε = 0.1. Plot the two respective solu- tions and the respective number of CG iterations. What do you observe?
[5 marks]
3. General guidance and study support
The MS Teams group for COMP53**M Scientific Computation will be used for general support for this assignment. If your question would reveal parts of the answer to any problem, please send a private message to the module leader on MS Teams instead. You can also use the tutorial sessions to ask questions about coursework.
4. Assessment criteria and marking process
Assessment marks and feedback will be available on Minerva within
three weeks of the submission deadline. Late submissions are allowed 7

within 14 days of the original deadline providing that a request for an extension is submitted before the deadline. Standard late penalties apply for submissions without approved extensions.
5. Presentation and referencing
When writing mathematical formulas, use similar notation and sym- bols as during the lectures and tutorials. Hand-written sections for mathematical notation are acceptable but need to be clearly readable.
You may assume theorems and other results that have been presented during lectures and tutorials as known. Any other theorems need to be cited using standard citation practice.
6. Submission requirements
This is an individual piece of work. Submit your answers through Tur- nitin as one PDF document (generated either in Word or with LaTeX). You may use hand-written and scanned pages for mathematical formu- las, but these need to be clearly legible and the document must contain at least some typeset text or Turnitin will reject it. All submissions will be checked for academic integrity.
7. Academic misconduct and plagiarism
Academic integrity means engaging in good academic practice. This involves essential academic skills, such as keeping track of where you find ideas and information and referencing these accurately in your work.
By submitting this assignment you are confirming that the work is a true expression of your own work and ideas and that you have given credit to others where their work has contributed to yours.
8. Assessment/marking criteria grid
Total number of marks is 40, divided as follows:
Exercise 1 (One-dimensional Burgers equation): 20 marks
Exercise 2 (Anisotropic diffusion and conjugate gradient): 20 marks
請加QQ：99515681 或郵箱：99515681@qq.com   WX：codehelp