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 non-public URL. 
(iii) wiki
In your 2-page 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  
      |
   O--O--O
      |
      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 Courant-Fridrichs-
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 space-time 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 1-D 
_________________________________________________________

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 nabla-squared 
(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 1-D 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 3-point 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
 \
  \
   m-----m

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=1e-4
(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 10-fold 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?