Assignment set #3  PHYD57  Adv. Comp. Methods for Physics
===========================================================
Due 20 Mar 2020 (23:59)
You can arbitrarily mix handwritten answers with printouts, in any proportion
you like, but please staple them or else number pages consecutively and mark
each with your name & student id (last 3 digits).
Problem 3.1  Research GPGPU (i.e., GPU computing)
___________________________________________________
We are going to talk about GPU computing during Lectures, but will not
cover in detail the history of GPU computing. Study it using numerous
resources on
(i) course web page, Links > Parallel computing  CUDA C & CUDA Fortran,
(ii) auxiliary page with nonpublic URL.
(iii) wiki
In your 2page writeup (max. 3 pages please) mention not just the history
in the sense of who did what & when, but also how GPU computing differs
from CPU computing. You may include some diagrams or short code snippets
(in Word or Latex or jupyter notebook format).
Try to capture all the essential things; it's easy to get carried away
citing various interesting but not essential information and suddenly
run out space (no more than 3 pages please). Cite more than two sources,
mention the types of applications GPU has found, and personal experience
with graphics computing, if any.
Problem 3.2  The wave equation. CFL limit
__________________________________________
Develop, using von Neuman's analysis [where spatial form of
error is taken as u = exp(i k x)], an analytical criterion of
stability for the method of solution of the homogenous,
linear wave equation
d^2 W/ dt^2 = c^2 (d^2 W/dx^2 + d^2 W/dy^2 +...)
c = wave speed = const. ~ 300 m/s.
Assume a numerical stencil of the type
X

OOO

O
(where time is running upwards and space is extending horizontally; here we
illustrate only one spatial dimention, but the equation can be either 1D, 2D
or 3D).
The criterion for the minimum time step is called CFL or CourantFridrichs
Levi criterion. Coefficients in it will depend on k=1,2,3, the number of
dimensions.
Interprete the result by comparing the numerical speed of propagation of
information on the uniform spacetime grid, and the physical speed equal c.
HINT: We derived a criterion for diffusion equation in a similar way.
Check the recordings of lectures/tutorials.
Problem 3.3  Simulate homogeneous wave equation in 1D
_________________________________________________________
This is one of the most celebrated linear equations of mathematical physics,
please read this article: https://en.wikipedia.org/wiki/Wave_equation.
Equation for deflection W(x,y,z) reads
W'' = c^2 (d^2 W/dx^2) (W'' means d^2 W/dt^2.)
We could write the parenthesis on the r.h.s. differently, using nablasquared
(Laplace operator). Nabla is a collection (a vector) of 3 partial drivatives
along 3 directions x,y,z: nabla = (dW/dx,dW/dy,dW/dz), and nabla^2 =
nabla * nabla, where * stands for dot product of vectors. In our case we'll
solve a 1D equation, so nabla*nabla operator is a second derivative
d^2 W / dx^2.
In our problem c = 1 m/s, is the speed of the wavefront. Spatial coordinate is
between 0 and 1 values of the x variable. Space is resolved into
N = 600 cells. When declaring array(s) remember to leave a little extra
space for two boundary cells, sometimes called ghost cells.
Boundary condition is closed or fixed, that is the value of
solution in the ghost cells at the beginning & end of the solution
vector W[601], will be kept constant in time and equal W[0] = W[601] = 0.
The concrete index values may differ depending on language and coding.
What is the best time step dt? Can you show what numerical trouble accompanies
a time step larger than the one given by CLF condition? All details of
how to implement the method are up to you, but it has to be done in C(++)
or Fortran.
Wave equation "only" differs from the diffusion equation in having second
instead of first time derivative on the left side.
You already know how to numerically approximate the second derivative
on a 3point stencil. Use it for time axis as well.
Set timestep to slightly less than the CLF minimum value.
As initial shape of W(x;t) choose a flat W=0 everywhere
except in the 21 central cells, where W = 1. Wave eq. is a second order partial
diff. eq, thus it requires specification not only of W but also W' in
the initial state. Assume W'(x;0) = 0 in the whole computational domain.
Make a set of plots of W(x) in the computational domain from
index 1 to 600, or overlay many W(x;t) curves on one plot,
at 4 times t in the interval [1,2]. All other details are up to you.
Problem 3.4  Double pendulum in C or Fortran (double prec.)
_____________________________________________________________
Find out, e.g. from Classical Mechanics textbook by Taylor (2004) or
wikipedia, what equations govern a frictionless double pendulum made of
two massless rigid struts with two equal masses
(e.g., m=1kg or some other value; it should not matter: do you know why?)
attached like so (x is showing a fixed but a freely rotating attachment point).
x
\
\
mm
The two arms have the same length of 0.5 m. Two degrees of freedom are
p1 and p2, two angles of the arms counted from the vertical line.
Problem: Starting from p1 = 0 (upward position of the first rod), and p2=1e4
(nearly upward second rod, a bit to the left of the first mass).
Integrate the Newtonian equations of motion from t=0s to time t_max = 8s.
Use a 4th order symplectic integrator described by Yoshida
https://planets.utsc.utoronto.ca/~pawel/PHYD57/yoshida.pdf
earlier discovered by Ruth and Forest in 1980s.
The double pendulum does not have any friction, and so its kinetic plus
gravit. potential energy should be conserved. Track the deviation
of the total energy E from the initial value, and print it at the end of
each integration, together with final angles p1 and p2. No graphic output
is required.
Run the program 3 times with these timesteps:
dt = 0.01s, 0.001s, and 0.0001s. Compare the magnitude of energy error and
position (assume that the 3rd position is precise and compare the
error of p1 and p2 w.r.t. this assumed exact result in the 0.01 and 0.001s
timestep cases. Do the errors drop by 10**(4) for a 10fold decrease
of timestep? Why (or why not)?
Finally, repeat the study integrating in the time interval t = 0...80s.
Did you get the same answers to the last questions above?