Assignment 3, PHYD57

Due date and time: 8 Nov (Fri) at 11am

Problem 3.1 [25 pts] Research the history of GPGPU (i.e., GPU computing)

Problem 3.2 [25 pts] Stencils, stability and CFL numbers

Introduction and the definition of the CFL number

The CFL (Courant-Fredrichs-Levy) condition connects $\mathrm{\Delta }x$, $\mathrm{\Delta }y$, ... and $\mathrm{\Delta }t$ in wave propagation and fluid flow problems, governed by the so-called hyperbolic equations like Navier-Stokes equation, ideal gasdynamics equations called Euler equations, various water wave equations (deep water, shallow water equations). Essentially, the meaning of such conditions is to specify when a given computational stencil, i.e. pattern of combining the data from previous time step to produce the updated value at a given location in the next time step, works ok, and when it produces instabilities. Instabilities are not just truncation errors, meaning - inaccuracies caused by numerical approximations to math oprations like the derivatives and integrals. Instabilities introduce expopnentially growing errors quickly making the solution non-physical.
There are so-called implicit schemes, which use stencils transferring information about the future (next timestep) state of the spatial neighbors of a given cell, in order to compute the future state of that cell, and many of them don't have a CFL-like condition, because they always are stable, no matter how large $dt=\mathrm{\Delta }t$ one uses. Such schemes are said to be unconditionally stable. The problem is that they require information (about neighbors) which is not yet available. In practice, one can device iterative scheme, which will derive better and better approximations to all the future values of variables on the grid, but there is a price to pay for this: implicit schemes are computationally expensive. This forces the user to adopt larger timesteps than in explicit schemes, with adverse effect on accuracy. Some of the best available CFD algorthms (Computational Fluid Dynamics) are therefore explicit.

You have learned much about the approximators of the 1st and 2nd derivatives on a grid. Now is the time to see whether they work in the context of some popular PDEs. One example that we consdered in Lectures is the von Neumann analysis of error propagation in the heat (or diffusion) equation ${\partial }_{t}T=D{\partial }_{xx}T$, when the first order, forward time differencing is used for time, and 2nd order symmetric formula for the r.h.s. 2nd derivative is used. The stencil looks like this

 
    *
    |
o___o___o

The physical meaning of $(\mathrm{\Delta }x{)}^2/D$ is the diffusion time after which a local peak or trough of $T$ of width $\mathrm{\Delta }x$ gets significantly leveled by diffusion. Timestep is a couple of times smaller than diffusion timescale across one cell, for the algorithm to be able to adequately describe the diffusion or temperature changes.

Stability of schemes for advection and Schrödinger equation

The time-dependent Schrödinger equation describes the time evolution of $\psi (\overrightarrow{r},t)$, the quantum-mechanical wave function: $$i\hslash {\partial }_t\psi (\overrightarrow{r},t)=[-\frac{{\hslash }^2}{2m}{\mathrm{\nabla }}^2+V(\overrightarrow{r})]\psi (\overrightarrow{r},t).$$ where $m$ is the mass of the particle represented by the wave function, $\hslash =h/(2\pi )$ is reduced Planck's constant, and ${\mathrm{\nabla }}^2$ is Laplace operator. In 1-D, ${\mathrm{\nabla }}^2={\partial }_{xx}$, in 2-D it is ${\mathrm{\nabla }}^2={\partial }_{xx}+{\partial }_{yy}$, and so on. We are interested in the simple case of free-space propagation of the wave function, $V(\overrightarrow{r})=0$.

Show that discretization of either of the above equations using the same numerical approximations to ${\partial }_t$ and ${\mathrm{\nabla }}^2$ as we discussed above in case of the diffusion equation, is unconditionally unstable. That is, there is no CFL relation guaranteeng stability, because the numerical solution always blows up.

So what do we do in such a situation? We have to employ more sophisticated explicit methods, or implicit methods.

The Klein-Gordon and the wave equation. CFL limit

One attempt was to formulate a wave equation patterned after the relativistic energy triangle equation $$E^2=(m{c}^2{)}^2+(pc{)}^2$$ where $m$ is the rest mass of a particle, $p$ its relativistic momentum (in general not equal to $mv$ since that product would be limited in value because of $v\le c$ requirement), and $E$ total energy. This gave the Klein-Gordon equation which, as we now know, describes the spin-zero particles such as the Higgs boson.

Let's consider the application to massles particles first. If $m=0$, the Klein-Gordon equation looks exactly like the familiar wave equation of clasical physics, although the function $\psi$ is the same quantum-mechanical wave function in Schrödinger equation, and is in general complex-valued. $${\partial }_{tt}\psi =c^2{\mathrm{\nabla }}^2\psi$$

Use von Neuman's analysis, where spatial and time dependence of error is assumed to be split into Fourier components that depend on position $\overrightarrow{r}$ and time like $e^{i(\overrightarrow{k}\cdot \overrightarrow{r}-\omega t)}$, the $\overrightarrow{k}$ being the wave vector and $\omega$ the wave's frequency. Show that any wave equation (including the homogeneous, mass-less Klein-Gordon equation) can be successfully solved numerically with the 2nd derivative approximator (2nd order accuracy) used in heat equation. Derive the critical CFL number. The explicit numerical stencil has the form

 
    *
    |
o___o___o
    |
    o

HINT: Follow the philosophy that was used for diffusion equation in L5, use the trick $e^{iz}+e^{-iz}-2=-4{\sin }^{2}z/2$. Judge at which values of $q$ there are no real-valued $\omega$'s that satisfy the equation, assuming a real-valued $k$ (or $\overrightarrow{k}$ if dimension $n>1$).

OPTIONAL PART: Finally, if you dare, derive the stability criterion for the Klein-Gordon equation with $m>0$, discretized as above. This is not necessary to get full points for this problem! But if you do it, some mistakes in the earlier parts of this problem will be forgiven.

Problem 3.3 [40 pts] N-body problem: Gravitational collapse

Consider an initially motionless, uniform-density, spherical cloud of randomly placed $N=1000$ points, extending to radius $R=1$. The total mass is equal $M=1$, and each body's mass is $m=M/N$. The gravitational constant is $G=1$. [You are free to use the actual physical constants and dimensional values. I'll not omit $M$ ($m$) or $G$ in the formulae below. The task you're facing will not depend on what values of $R$ and $M$ you use. Only the $dt$ and the simulated times will have different values, but the number of steps taken or wall clock time should not be much different.]

The acceleration acting on particle $i$ is equal $${\overrightarrow{f}}_i=Gm\sum _{j\ne i}\frac{{\overrightarrow{r}}_{ij}}{(r_{ij}^2+{ϵ}^2{)}^{3/2}}$$ where ${\overrightarrow{r}}_{ij}={\overrightarrow{r}}_j-{\overrightarrow{r}}_i$, $r_{ij}=|{\overrightarrow{r}}_{ij}|$ , and the sum is done over all particles other than $i$. It's not required, but you may use non-zero softening parameter $ϵ\approx {10}^{-5}$, effectively limiting the strength of gravity force at distances $r<ϵ$. It's rare that particles are that close to each other, but then their trajectories will be computed more reasonably (as if gravity was actually non-Newtonian at small distances, or if gravity was Newtonian but particles were really fuzzy balls of radius $\sim ϵ$ of much smaller components. In astronomy, such objects with potential $-G{m}^2/\sqrt{r^2+{ϵ}^2}$ are known as Plummer spheres.)

Monitor the relative change of the total mechanical energy $(E-E_0)/E_0$ and the instantaneous total momentum of the system divided by $N$, where $E_0=E(t=0)$ and $$E=m\sum _i[\frac{v_i^2}{2}-Gm\sum _{j<i}\frac{1}{\sqrt{r_{ij}^2+{ϵ}^2}}],\overrightarrow{p}=-\frac{m}{N}\sum _i{\overrightarrow{v}}_i.$$

Choose a constant timestep $dt$ that assures good conservations of energy and nearly zero total momentum during the simulation. State the absolute value of momentum divided by $M$. Try to keep the relative energy error below ${10}^{-5}$. Consider any results where the relative error is larger than 0.01 garbage, affected by a very close encounter of particles. If the running time is too long, kill the program by CTRL-C after 10 or 15 minutes and run it again with much smaller $N$, it should then finish quickly. Establish the biggest N that gives reasonably short runtime by trial and error guided by the method of bisection. Hopefully you won't have to decrease N, but who knows..

You will finish the computation once the following happens. Write a subroutine or procedure, which you will call once in a while, not after every timestep, since that would slow down the working of the code unnecessarily, but maybe every 10 or more steps. It will return some information about the system, such as $E$ and $\overrightarrow{p}$, as well as the mean radius of the system $R_m=(1/N)\sum _i|\overrightarrow{r_i}|$ and the mean velocity, $V_m=(1/N)\sum _i|\overrightarrow{v_i}|$. The ratio $t_c=R_m/V_m$ tells you the average crossing time, which is similar to the typical infall timescale (initially) and the orbital period of a particle at the end. The code should finish when the time exceeds $20t_c$, by which time we will consider that the system settled into equilibrium.

The main goal of your code is to find by what factor the mean radius of the system in final equilibrium is smaller than the initial mean radius. (The virial theorem can be used to make a prediction for that ratio, and you can try to make the prediction, but that is not required.)

Some particles may escape. Don't include in the average radius those, which escaped beyond 10 times the initial ball's radius, and don't plot them (you can choose to plot only the particles within 3 initial radii). But study wherether they escaped because of code inaccuracies- by varying dt. If there is a noticeble change in results when making energy more exactly conserved (smaller dt) then it might mean that even dE/E ~ 1e-5 is not sufficient to converge to a reliable answer. Consider adding a small number called softening parameter ε² ~ 10^-10 to the square of every particle-particle distance to limit acceleration and decrease errors due to close encounters. Cf. the ${\overrightarrow{f}}_i$ for a proper way to modify the equation.
Cluster evaporation is a physical phenomenon, but the theoretical timescale of it (to lose significant part of the cluster) much exceeds 20 crossing times.

Code details need to be described in comments, and in a separate text that will also include the results and their discussion.

Problem 3.4 [10 pts] A variant of code in 3.3

What about the following idea for another optimization. The method, as described above, computes a given pair $ij$ twice, once while summing the acceleration on particle $i$, and another time doing summation for particle $j$. Meanwhile, Newton's 3rd law states that the second calculation of force is redundant, because the force is equal and opposite to the one already computed (both acceleration can be updated simultaneously). One can only consider pairs that satisfy $i<j$ or $j<i$, to reduce the number of 3-D distances calculated. But will that work well with OMP? You can try and tell us what happened. (We are looking forward to your bright ideas on code optimization. E.g., you might check whether 9 by N or N by 9 arrays work faster, whether putting them on stack or heap, by palcing the declaration in or outside a procedure or fortran module, works better. Just remember to increase OS stack size if you get error messages having to do with memory allocations, or segmentation faults, which is the typical symptom of those troubles. This is just a suggestion. Describe your efforts even if they didn't bring much benefit. We ask you to describe your solution to 3.3 and 3.4 in a text file, points will be cut for submitting only the codes.)