This exercise aims to provide a flavor of the expertise required for the SPIN project as well as giving you an idea of the work expected of you. If you choose to apply to work with me on this project then I will ask you to work through the exercise and send in your results before the interview with me. I provide a worked out solution so you are free to look at it to verify your own solution or to copy and past from it. However please keep in mind that the exercise mostly serves to give you an idea of the skills required to avoid signing up for more than you can handle.
Background
Numerical relativity solves Einstein's equations of General Relativity using supercomputers to predict for example how two neutron stars collide and how gravitational waves such as the ones detected by LIGO are generated in such events. Computer simulations require initial data that they then evolve forwarded in time using the computer. In General Relativity such initial data has to satisfy what is called initial data constraints that are typically of the form
Δ φ =−4 πρ
in some domain Ω with a boundary condition on the edge of the domain ∂ Ω. Typically the coordinates used to discretize the problem are chosen to match with the domain, for example using spherical coordinates if Ω is a sphere. For very irregular Ω however it is convenient to stick to Cartesian coordinates x,y,z and instead embed the boundary ∂ Ω into the solution domain by turning the condition in ∂ Ω into a condition in Ω. The project aims to explore different methods of doing so. A starting point will be this paper http://www.cs.ucr.edu/~shinar/papers/2007_boundary_condition.pdf that however handles a slightly different situation but can be used as a source of references to dig down.
In this exercise we will use a very simple, square domain Ω = [0,1] × [0,1] and employ a very simple numerical method called the Jacobi method (Wikipedia) to solve the elliptic equation. Basically one starts with a guess for φ and then at each point (x,y) one replace φ(x,y) by the average of the φ values around it (which on a Cartesian grid are its 4 neighbours), then repeats until φ no longer changes. This (eventually, it is very slow) converges to the solution.
The file poisson.py solves the equation above for ρ = 0 and boundary conditions φ(x=0,y) = φ(x,y=0) = 0, φ(x=1,y) = y, φ(x,y=1) = x which has the solution φ(x,y) = x y . The Jacobi iteration itself is in these lines of code:
# update phi according to Jacobi method formula
# https://en.wikipedia.org/wiki/Jacobi_method
del_phi = phi[2:,1:-1] + phi[0:-2,1:-1] + \
phi[1:-1,2:] + phi[1:-1,0:-2]
phi_new[1:-1,1:-1] = -h**2/4. * (source[1:-1,1:-1] - del_phi/h**2)
while boundary conditions are enforced here
# enforce boundary condition sol = solution(xmesh, ymesh) phi[ 0, :] = sol[ 0, :] # x = xmin phi[-1, :] = sol[-1, :] # x = xmax phi[ :, 0] = sol[ :, 0] # x = ymin phi[ :,-1] = sol[ :,-1] # x = ymax
After 100 iterations the difference between the correct solution and the current approximation looks like this
As an exercise change the code so that instead it solves
Δ φ = -2 π² sin(π x) sin(π y)
using the same boundary conditions φ(x=0,y) = φ(x,y=0) = 0, φ(x=1,y) = y, φ(x,y=1) = x . This has a solution φ(x,y) = x y + sin(π x) sin(π y) .
After 10000 iterations the difference between the correct solution and the current approximation looks like this
You can find my solution in poisson_solved.py .