{
"cells": [
{
"cell_type": "markdown",
"id": "5451ef11-f683-4995-bda8-c9d87abaec49",
"metadata": {},
"source": [
"# Particle Dynamics with Python\n",
"<div class=\"dateauthor\">\n",
"07 June 2021 | Jan H. Meinke\n",
"</div>"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5822f3b3-bc03-4e2f-85f1-57cb246e3a05",
"metadata": {},
"outputs": [],
"source": [
"import math\n",
"import random\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f7d1939b-7d73-4c0c-9d8a-d6ea39d48b49",
"metadata": {},
"outputs": [],
"source": [
"%matplotlib widget"
]
},
{
"cell_type": "markdown",
"id": "19819e70-b42c-405a-958f-70c05a972ee6",
"metadata": {},
"source": [
"Particle dynamics simulations are common in various scientific fields. They are used to simulate \n",
"the formation of galaxies and the movements of molecules in a cell. Particles can have different\n",
"properties such as mass and charge and interact in different ways.\n",
"\n",
"A classical particle dynamics code solves Newton's equation of motion:\n",
"\n",
"$$\\mathbf F = m \\mathbf a,$$\n",
"\n",
"where $\\mathbf F$ is the force, $m$ the mass, and $\\mathbf a$ the acceleration. $\\mathbf F$ and \n",
"$\\mathbf a$ are vectors.\n",
"\n",
"In general, this problem is only solvable analytically for two particles. If there are more \n",
"particles, we have to look for a numerical solution.\n",
"\n",
"You may remember that you can calculate the velocity $\\mathbf v$ of a particle as\n",
"\n",
"$$\\mathbf v(t + dt) = \\mathbf v(t) + \\mathbf a(t) dt$$\n",
"\n",
"and the position $\\mathbf r$ as\n",
"\n",
"$$\\mathbf r(t + dt) = \\mathbf r(t) + \\mathbf v(t)dt + \\frac 1 2 \\mathbf a(t) dt^2.$$\n",
"\n",
"If we know all the positions, velocities and masses at time $t$ and can calculate the forces, we \n",
"can follow the motion of the particles over time."
]
},
{
"cell_type": "markdown",
"id": "50ad1731-c5b0-4922-adc8-14e507a7b6b8",
"metadata": {},
"source": [
"## Gravitational force\n",
"Let's assume our particles only interact via gravity. Then the force between two particles is given \n",
"by\n",
"\n",
"$$\\mathbf F_{ij}(t) = G\\frac{m_i m_j}{r_{ij}^2(t)} \\mathbf {\\hat r}_{ij}(t),$$\n",
"\n",
"where $\\mathbf F_{ij}(t)$ is the force on particle $i$ due to particle $j$. $r_{ij}(t)$ is the \n",
"distance between particles $i$ and $j$, and $\\mathbf {\\hat r}_{ij}(t)$ is the unit vector pointing\n",
"from $i$ to $j$.\n",
"\n",
"To get the total force on particle $i$, we need to sum over all $j \\neq i$:\n",
"\n",
"$$\\mathbf F_{i}(t) = \\sum_{j\\neq i} \\mathbf F_{ij}(t).$$"
]
},
{
"cell_type": "markdown",
"id": "32f7c975-ed21-4c70-9168-5b7bfa5ca276",
"metadata": {},
"source": [
"## The algorithm"
]
},
{
"cell_type": "markdown",
"id": "4288de12-8bf3-41b2-96ca-5c3c47fc0d84",
"metadata": {},
"source": [
"1. Calculate the force on each particle by summing up all the forces acting on it.\n",
"2. Integrate the equation of motion\n",
"\n",
" a) Calculate the position of each particle after a time step *dt*\n",
" \n",
" b) Calculate the velocity of each particle after a time step *dt*"
]
},
{
"cell_type": "markdown",
"id": "539c2d60-df7b-471b-a438-d9b4efb51781",
"metadata": {},
"source": [
"## (Parallel) Patterns\n",
"In Think Vector, we got to know some patterns. Let's see how we can apply them here:\n",
"\n",
"(i, j) -> $\\mathbf F_{ij}$:\n",
" This is a map\n",
" \n",
"$\\mathbf F_{ij}$ -> $\\mathbf F_{i}$:\n",
" This is a reduction\n",
" \n",
"Calulate the new velocity:\n",
" This is map\n",
" \n",
"Calculate the new position:\n",
" This is a map, too.\n",
" \n",
"Now, let's try to express this in code.\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b4525c8a-378a-45b7-b1e2-b67f5f07d397",
"metadata": {},
"outputs": [],
"source": [
"# Initialize positions and velocities\n",
"N = 50\n",
"L2 = 5 # half the length of the box size\n",
"epsilon = 1e-9 # softening factor\n",
"dt = 0.1 # time step\n",
"G = 1 # For simplicity we set the universal graviational constant to 1\n",
"m = 1 # This corresponds to 150 x 10^9 kg\n",
"x = [random.uniform(-L2, L2) for i in range(N)]\n",
"y = [random.uniform(-L2, L2) for i in range(N)]\n",
"z = [random.uniform(-L2, L2) for i in range(N)]\n",
"vx = [0 for i in range(N)]\n",
"vy = [0 for i in range(N)]\n",
"vz = [0 for i in range(N)]"
]
},
{
"cell_type": "markdown",
"id": "8fd053d2-8c88-4666-82ed-0316fe21ac34",
"metadata": {},
"source": [
"### Calculating forces"
]
},
{
"cell_type": "markdown",
"id": "41861767-e08d-45b8-802a-28b269e3f7ee",
"metadata": {},
"source": [
"To calculate the force, we need the distance vector first. These are actually 3 maps (one for each component). The result is a distance matrix. As mentioned before maps are expressed as list generators:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "338142b6-f973-4f7a-b5a4-77e76f3b758f",
"metadata": {},
"outputs": [],
"source": [
"Dxx = [(i - j) for j in x for i in x]\n",
"Dyy = [(i - j) for j in y for i in y]\n",
"Dzz = [(i - j) for j in z for i in z]\n",
"D = [math.sqrt(i * i + j * j + k * k) for i, j, k in zip(Dxx, Dyy, Dzz)]"
]
},
{
"cell_type": "markdown",
"id": "d0156a2d-13ae-46dd-b3a8-cb7eb1aca0bf",
"metadata": {},
"source": [
"Now that we have the vector components and the magnitude of the vector, we can calculate the forces."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e841a076-504d-445b-b006-b931e3cb0bc2",
"metadata": {},
"outputs": [],
"source": [
"Fxx = [G * m * m * i / (d * d * d + epsilon) for i, d in zip(Dxx, D)] # epsilon prevents a zero in the dominator.\n",
"Fyy = [G * m * m * i / (d * d * d + epsilon) for i, d in zip(Dyy, D)]\n",
"Fzz = [G * m * m * i / (d * d * d + epsilon) for i, d in zip(Dzz, D)]\n",
"Fx = [sum(Fxx[i * N: (i + 1) * N]) for i in range(N)]\n",
"Fy = [sum(Fyy[i * N: (i + 1) * N]) for i in range(N)]\n",
"Fz = [sum(Fzz[i * N: (i + 1) * N]) for i in range(N)]"
]
},
{
"cell_type": "markdown",
"id": "3de052ac-7591-4477-8285-cc15c0019a7a",
"metadata": {},
"source": [
"Let's visualize the forces on the particles:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1133b4bb-111b-4aca-9326-22a7c29c8522",
"metadata": {},
"outputs": [],
"source": [
"ax = plt.figure(figsize=(6, 6)).add_subplot(projection='3d')\n",
"ax.scatter3D(x, y, z)\n",
"ax.quiver(x, y, z, Fx, Fy, Fz)"
]
},
{
"cell_type": "markdown",
"id": "ccea23e5-4f4b-4ff6-b379-8d45e3fe15f4",
"metadata": {},
"source": [
"### Integrating the equation of motion"
]
},
{
"cell_type": "markdown",
"id": "dba27f9b-350e-4e65-9f42-e3615ee30a84",
"metadata": {},
"source": [
"We are ready to update the positions and velocities of our particles:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5ddc24f9-eaf3-491c-bf81-232efa584c1c",
"metadata": {},
"outputs": [],
"source": [
"x = [i + v * dt + 0.5 * f / m * dt * dt for i, v, f in zip(x, vx, Fx)]\n",
"y = [i + v * dt + 0.5 * f / m * dt * dt for i, v, f in zip(y, vy, Fy)]\n",
"z = [i + v * dt + 0.5 * f / m * dt * dt for i, v, f in zip(z, vz, Fz)]"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2266d4e8-8f67-4979-ae47-abf8508673a4",
"metadata": {},
"outputs": [],
"source": [
"vx = [v + f / m * dt for v, f in zip(vx, Fx)]\n",
"vy = [v + f / m * dt for v, f in zip(vy, Fy)]\n",
"vz = [v + f / m * dt for v, f in zip(vz, Fz)]"
]
},
{
"cell_type": "markdown",
"id": "e4cff076-759c-477c-9758-41bb730cd606",
"metadata": {},
"source": [
"Let's take a look at the particle positions and velocities:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bf48e0d0-34f6-47ba-8a30-ba0c1e19489d",
"metadata": {},
"outputs": [],
"source": [
"ax = plt.figure(figsize=(6, 6)).add_subplot(projection='3d')\n",
"ax.scatter3D(x, y, z)\n",
"ax.quiver(x, y, z, vx, vy, vz)"
]
},
{
"cell_type": "markdown",
"id": "65984f53-4b54-4f6d-aaa1-6de391150539",
"metadata": {},
"source": [
"That's it. By going back to the [calculation of the forces](#Calculating-forces), we can follow the motion of the particles over time."
]
},
{
"cell_type": "markdown",
"id": "f1f30004-a9c3-4499-84e0-976937b9f8a8",
"metadata": {},
"source": [
"## Exercise\n",
"Rewrite the program in a vectorized manner using `ndarray`s."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "039819a6-698f-43a6-a4f0-4f7b8852fbb1",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"id": "8cb45f43-29e2-49df-a976-bf7790fe5a44",
"metadata": {},
"source": [
"### Solution:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1f236119-af8c-499d-86cf-1d6b98f9e5fd",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "HPC Python 2021",
"language": "python",
"name": "hpcpy21"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.5"
}
},
"nbformat": 4,
"nbformat_minor": 5
}