Assignment

This is simultaneously a mid-term assignment and the final take home exam. The assignment should exercise you in putting together some of the ideas related to pushing particles that we have talked about and give you an opportunity to try your hands on MPI too.

The assignment is due by Friday 14th of May 1999 or earlier if you're so inclined. The assignment results should be placed in your AFS directory, preferably in a separate subdirectory. There should be an accompanying README file, which explains the content of other files in that directory. You should notify me by e-mail, when you have submitted your assignment.

This assignment is for your benefit only, i.e., it is not related to your final grade, that being based on your class attendance, participation in laboratory exercises, questions asked during the course, and so on$\ldots$

Remember: as this is a PhD level course, you're now colleagues, and not pupils. And one doesn't grade one's colleagues. But, of course, one is perfectly at ease rejecting their papers and grant applications.

So, here's the assignment itself:

Exercise 1

Write an MPI (C or F90) or a High Performance Fortran program which traces trajectories of a swarm of dust particles in a central gravitational field. The force that acts on every dust particle is

$$\boldsymbol{F} = - \frac{M m}{r^2} \frac{\boldsymbol{r}}{r}$$ (8.1)

where r is the guiding vector of the dust particle, and the source of the force is located at r = 0. The mass associated with the source is M and the mass of the dust particle is m. We have also assumed that the gravitational constant G = 1.

Because

m a = F (8.2)

the mass of the dust particle cancels out and we get:

$$\frac{\textrm{d}^2 \boldsymbol{r}} {\textrm{d} t^2} = - \frac{M}{r^2} \frac{\boldsymbol{r}}{r}$$ (8.3)

For simplicity assume that M = 1, so that

$$\frac{\textrm{d}^2 \boldsymbol{r}} {\textrm{d} t^2} = - \frac{1}{r^2} \frac{\boldsymbol{r}}{r}$$ (8.4)

initial condition

The potential energy of a dust particle in the gravitational field is given by

$$U(\boldsymbol{r}) = - \frac{Mm}{r}$$ (8.5)

Observe that U(r) becomes $-\infty$ for r=0 and converges to 0 as $r \rightarrow \infty$. Gravitational energy is always negative.

The kinetic energy of the particle is

$$K(\boldsymbol{v}) = \frac{mv^2}{2}$$ (8.6)

The kinetic energy of dust particles is always positive.

The total energy of a particle is

$$E(\boldsymbol{r}, \boldsymbol{v}) = \frac{mv^2}{2} - \frac{Mm}{r}$$ (8.7)

If the total energy is negative the particle is trapped in the gravitational well.

Generate 10,000 particle positions and velocities at random ensuring that

1.

each particle is confined within a sphere of radius 1, where the centre of the sphere coincides with the source of the central force;

2.

the energy of each particle is negative.

evolution

Evolve the system using either Runge-Kutta method or Bulirsch Stoer method for at least 3 revolutions of the slowest particle around the centre of the system.

measurements

For each particle perform the following measurements:

1.

measure the total energy of the particle: does the total energy remain constant? It should.

2.

measure the angular momentum of the particle:

$$\boldsymbol{L} = m \boldsymbol{r}\times\boldsymbol{v}$$ (8.8)

where $\times$ is the vector product (also called the cross product) of two vectors. L should stay constant both with regards to its length and with regards to its direction. Does it?

3.

as the particle moves around the central point of the system it assumes positions $\boldsymbol{r}(t_1), \boldsymbol{r}(t_2), \ldots$. If t₂ - t₁ is very small, then the area of the sector whose apex is located at the centre of the system and whose two sides are given by r(t₁) and r(t₂) is

$${\cal A} = \vert\boldsymbol{r}(t_1) \times \boldsymbol{r}(t_2)\vert$$ (8.9)

The time $\Delta t = t_2 - t_1$ taken to traverse from r(t₁) to r(t₂) should be proportional to ${\cal A}$ along the whole trajectory of an individual dust particle. Is it?

4.

measure the time it takes for the particle to travel the full circle around the centre of the system. You may have to resort to interpolation in order to measure that time exactly. Let us call that time T.

5.

measure the smallest distance between the particle and the centre of the system attained throughout a single revolution around the centre of the system. You may have to resort to interpolation in order to measure that distance exactly. Let us call that distance $r_{\textrm{{\scriptsize min}}}$

6.

evaluate $T^2/r_{\textrm{{\scriptsize min}}}^3$ for all particles. The number should be the same for all of them. Is it?

plots

Use GNUplot to plot the following quantities for selected particles:

• the total energy, L, and ${\cal A} / \Delta t$ in function of time.
• the trajectory, i.e., r in function of time
• $T^2/r_{\textrm{{\scriptsize min}}}^3$ in function of $r_{\textrm{{\scriptsize min}}}$

Exercise 2

This time assume that all particles of the system interact with each other as follows:

$$m_i \frac{\textrm{d}^2{\boldsymbol{r}}_i} {\textrm{d} t^2} ... ...symbol{r}_j} {\vert\boldsymbol{r}_i - \boldsymbol{r}_j\vert}$$ (8.10)

Observe that if r_j = 0 then the equation of motion is the same as for the central force. The direction of the vector r_i - r_j points from particle j to particle i, so - (r_i - r_j) points from particle i to particle j, therefore particle i is attracted to particle j.

initial condition

The total energy of the system is the sum of all kinetic and potential energies of individual particles. The kinetic energy of an individual particle is, as before,

$$K_i = \frac{m_i v_i^2}{2}$$ (8.11)

The potential energy of an individual particle is

$$U_i = \sum_{j \neq i} - \frac{m_i m_j} {\vert\boldsymbol{r}_i - \boldsymbol{r}_j\vert}$$ (8.12)

And the total energy is

$$E = \sum_i K_i + U_i$$ (8.13)

• generate 10,000 random positions for the particles within the sphere of radius 1.
• Assume all masses to be 1.
• For each particle evaluate its potential energy and generate random velocity so that the abolute value of its kinetic energy is less than the abolute value of its potential energy. This guarantees that the total energy of the ensemble is negative.

evolution

evolve the system for a sufficient time to observe significant changes, e.g., some of the particles should traverse across the volume of the system a couple of times

measurements

the total energy of the system E and the total angular momentum:

$$\sum_i m_i \boldsymbol{r}_i \times\boldsymbol{v}_i$$ (8.14)

should stay constant. Do they?

observations

you may notice some particles acquiring a positive total energy at a cost of some other particles cooling down. That is how the ensemble evaporates. At every step look for a particle with the lowest total energy and for a particle with the highest total energy.

plots

• Use GNUplot to plot a number of selected trajectories.
• Use GNUplot to plot the lowest and the highest observed total energies in function of time.
• Use GNUplot to plot the total angular momentum and the total energy of the system in function of time.

Exercise 3

Repeat Exercise 2, but this time assume that half of the particles have m_i = -1, and the other half m_i = +1. Distribute negative and positive values of m at random over the ensemble.