probably faulty but running at least

nbody/checklist.md

@@ -1,18 +1,18 @@
 # N-Body project - Checklist
 
 ### Task 1
-- [ ] Compute characteristic quantities/scales
+- [x] Compute characteristic quantities/scales
 - [x] Compare analytical model and particle density distribution
-- [ ] Compute forces through nbody simulation
+- [x] Compute forces through nbody simulation
 - [x] vary softening length and compare results
 - [x] compare with the analytical expectation from Newtons 2nd law
-- [ ] compute the relaxation time
+- [x] compute the relaxation time
 
 ### Task 2 (particle mesh)
 - [ ] Choose reasonable units
-- [ ] Implement force computation on mesh
+- [~x] Implement force computation on mesh
-- [ ] Find optimal mesh size
+- [x] Find optimal mesh size
-- [ ] Compare with direct nbody simulation
+- [x] Compare with direct nbody simulation
 - [ ] Time integration for direct method AND mesh method
 
 

nbody/task1.ipynb

@@ -42,7 +42,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "16:03:58 - utils.load - Loaded 50010 rows and 10 columns from data/data.txt\n"
+      "08:59:02 - utils.load - Loaded 50010 rows and 10 columns from data/data.txt\n"
      ]
     }
    ],
@@ -104,8 +104,8 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "16:04:00 - task1 - Considering a globular cluster - total mass of particles: 4622219.258999999, maximum radius of particles: 724.689657812915\n",
+      "08:59:04 - task1 - Considering a globular cluster - total mass of particles: 4622219.258999999, maximum radius of particles: 724.689657812915\n",
-      "16:04:00 - utils.units - Set scales: M_SCALE = 0.022 solMass, R_SCALE = 0.028 pc\n"
+      "08:59:04 - utils.units - Set scales: M_SCALE = 0.022 solMass, R_SCALE = 0.028 pc\n"
      ]
     }
    ],
@@ -206,7 +206,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "16:04:00 - utils.particles - Found mean interparticle distance: 0.010402746349924056\n"
+      "08:59:05 - utils.particles - Found mean interparticle distance: 0.010402746349924056\n"
      ]
     }
    ],
@@ -235,56 +235,6 @@
    "cell_type": "code",
    "execution_count": 8,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "\"## compare the two force calculations\\n# since the forces were computed for each particle, rather than comparing them directly we compare the relative error in the magnitude and direction of the forces\\n\\n\\n# f_diff = f_nsquare_1e - f_analytical\\nf_diff = f_nsquare_2e - f_analytical\\ndiff_mag = np.linalg.norm(f_diff, axis=1)\\n\\n\\n# plot the distribution of the error\\n# create 4 stacked histograms, sharing the same x axis\\nfig, ax = plt.subplots(4, sharex=True)\\nax[0].hist(diff_mag, bins=NBINS)\\nax[0].set_title('Magnitude of the force difference')\\nax[0].set_yscale('log')\\n\\nax[1].hist(f_diff[:,0], bins=NBINS)\\nax[1].set_title('X component of the force difference')\\nax[1].set_yscale('log')\\n\\nax[2].hist(f_diff[:,1], bins=NBINS)\\nax[2].set_title('Y component of the force difference')\\nax[2].set_yscale('log')\\n\\nax[3].hist(f_diff[:,2], bins=NBINS)\\nax[3].set_title('Z component of the force difference')\\nax[3].set_yscale('log')\\n\\nplt.title('Error in forces')\\nplt.show()\\n\""
-      ]
-     },
-     "execution_count": 8,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "\"\"\"## compare the two force calculations\n",
-    "# since the forces were computed for each particle, rather than comparing them directly we compare the relative error in the magnitude and direction of the forces\n",
-    "\n",
-    "\n",
-    "# f_diff = f_nsquare_1e - f_analytical\n",
-    "f_diff = f_nsquare_2e - f_analytical\n",
-    "diff_mag = np.linalg.norm(f_diff, axis=1)\n",
-    "\n",
-    "\n",
-    "# plot the distribution of the error\n",
-    "# create 4 stacked histograms, sharing the same x axis\n",
-    "fig, ax = plt.subplots(4, sharex=True)\n",
-    "ax[0].hist(diff_mag, bins=NBINS)\n",
-    "ax[0].set_title('Magnitude of the force difference')\n",
-    "ax[0].set_yscale('log')\n",
-    "\n",
-    "ax[1].hist(f_diff[:,0], bins=NBINS)\n",
-    "ax[1].set_title('X component of the force difference')\n",
-    "ax[1].set_yscale('log')\n",
-    "\n",
-    "ax[2].hist(f_diff[:,1], bins=NBINS)\n",
-    "ax[2].set_title('Y component of the force difference')\n",
-    "ax[2].set_yscale('log')\n",
-    "\n",
-    "ax[3].hist(f_diff[:,2], bins=NBINS)\n",
-    "ax[3].set_title('Z component of the force difference')\n",
-    "ax[3].set_yscale('log')\n",
-    "\n",
-    "plt.title('Error in forces')\n",
-    "plt.show()\n",
-    "\"\"\""
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
    "outputs": [
     {
      "data": {
@@ -347,15 +297,15 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": 9,
    "metadata": {},
    "outputs": [
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "16:04:17 - task1 - Crossing time for half mass system: 1.7e-06 pc(3/2) / solMass(1/2)\n",
+      "08:59:24 - task1 - Crossing time for half mass system: 1.7e-06 pc(3/2) / solMass(1/2)\n",
-      "16:04:17 - task1 - Direct estimate of the relaxation timescale: 0.00078 pc(3/2) / solMass(1/2)\n"
+      "08:59:24 - task1 - Direct estimate of the relaxation timescale: 0.00078 pc(3/2) / solMass(1/2)\n"
      ]
     }
    ],

nbody/utils/forces_basic.py

@@ -33,13 +33,13 @@ def n_body_forces(particles: np.ndarray, G: float, softening: float = 0):
         # m is a list of scalars and displacements is a list of vectors (2D array)
         # we only want row_wise multiplication
         num = G * (masses * displacements.T).T
-        
+
         # a zero value is expected where we have the same particle
         r_adjusted[i] = 1
         num[i] = 0
-        
+
         f = np.sum((num.T / r_adjusted**1.5).T, axis=0) * m_current
-        forces[i] = -f
+        forces[i] = - f
 
         if i % 5000 == 0:
             logger.debug(f"Particle {i} done")
@@ -47,6 +47,24 @@ def n_body_forces(particles: np.ndarray, G: float, softening: float = 0):
     return forces
 
 
+def n_body_forces_basic(particles: np.ndarray, G: float, softening: float = 0):
+    if particles.shape[1] != 4:
+        raise ValueError("Particles array must have 4 columns: x, y, z, m")
+
+    x_vec = particles[:, 0:3]
+    masses = particles[:, 3]
+    n = particles.shape[0]
+    forces = np.zeros((n, 3))
+    for i in range(n):
+        for j in range(n):
+            if i == j:
+                continue # keep the value at zero
+            r_vec = x_vec[j] - x_vec[i]
+            r = np.linalg.norm(r_vec)
+            f = - G * masses[i] * masses[j] * r_vec / (r**3 + softening**3)
+            forces[i] += f
+
+    return forces
 
 def analytical_forces(particles: np.ndarray):
     """

nbody/utils/forces_mesh.py

@@ -75,7 +75,7 @@ def mesh_forces_v2(particles: np.ndarray, G: float, n_grid: int, mapping: callab
     if logger.level >= logging.DEBUG:
         show_mesh_information(phi, "Potential mesh")
         show_mesh_information(phi_grad[0], "Potential gradient")
-    logger.debug(f"Got phi_grad with: {phi_grad.shape}, {np.max(phi_grad)}")
+        logger.debug(f"Got phi_grad with: {phi_grad.shape}, {np.max(phi_grad)}")
     
     # compute the particle forces from the mesh potential
     forces = np.zeros_like(particles[:, :3])

nbody/utils/integrate.py

@@ -16,10 +16,9 @@ def ode_setup(particles: np.ndarray, force_function: callable) -> tuple[np.ndarr
     if particles.shape[1] != 7:
         raise ValueError("Particles array must have 7 columns: x, y, z, vx, vy, vz, m")
     
-    # for scipy integrators we need to flatten array which contains 7 columns for now
+    # for the integrators we need to flatten array which contains 7 columns for now
-    # we don't really care how we reshape as long as we unflatten consistently afterwards
+    # we don't really care how we reshape as long as we unflatten consistently
     particles = particles.reshape(-1, copy=False, order='A')
-    # this is consistent with the unflattening in to_particles()!
     logger.debug(f"Reshaped 7 columns into {particles.shape=}")
 
     def f(y, t):
@@ -29,33 +28,35 @@ def ode_setup(particles: np.ndarray, force_function: callable) -> tuple[np.ndarr
         """
         y = to_particles(y)
         # now y has shape (n, 7), with columns x, y, z, vx, vy, vz, m
-        
 
         forces = force_function(y[:, [0, 1, 2, -1]])
-        
+
         # compute the accelerations
         masses = y[:, -1]
         a = forces / masses[:, None]
         # the [:, None] is to force broadcasting in order to divide each row of forces by the corresponding mass
-        # a.flatten()
+
-        
+        dydt = np.zeros_like(y)
-        # replace some values in y:
         # the position columns become the velocities
         # the velocity columns become the accelerations
-        y[:, 0:3] = y[:, 3:6]
+        dydt[:, 0:3] = y[:, 3:6]
-        y[:, 3:6] = a
+        dydt[:, 3:6] = a
         # the masses remain unchanged
+        dydt[:, -1] = masses
 
         # flatten the array again
-        y = y.reshape(-1, copy=False, order='A')
+        # logger.debug(f"As particles: {y}")
-        return y
+        dydt = dydt.reshape(-1, copy=False, order='A')
+        
+        # logger.debug(f"As column: {y}")
+        return dydt
 
     return particles, f
 
 
 def to_particles(y: np.ndarray) -> np.ndarray:
     """
-    Converts the 1D array y into a 2D array IN PLACE
+    Converts the 1D array y into a 2D array
     The new shape is (n, 7) where n is the number of particles.
     The columns are x, y, z, vx, vy, vz, m
     """
@@ -64,10 +65,21 @@ def to_particles(y: np.ndarray) -> np.ndarray:
     
     n = y.size // 7
     y = y.reshape((n, 7), copy=False, order='F')
-    logger.debug(f"Unflattened array into {y.shape=}")
+    # logger.debug(f"Unflattened array into {y.shape=}")
     return y
 
 
+def to_particles_3d(y: np.ndarray) -> np.ndarray:
+    """
+    Converts the 2D sol array with one vector per timestep into a 3D array:
+    2d particles (nx7) x nsteps
+    """
+    n_steps = y.shape[0]
+    n_particles = y.shape[1] // 7
+    y = y.reshape((n_steps, n_particles, 7), copy=False, order='F')
+    # logger.debug(f"Unflattened array into {y.shape=}")
+    return y
+
 
 def runge_kutta_4(y0 : np.ndarray, t : float, f, dt : float):
     k1 = f(y0, t)

nbody/utils/logging_config.py

@@ -2,10 +2,10 @@ import logging
 
 logging.basicConfig(
     ## set logging level
-    # level=logging.INFO,
+    level = logging.DEBUG,
-    level=logging.INFO,
+    # level = logging.INFO,
-    format='%(asctime)s - %(name)s - %(message)s',
+    format = '%(asctime)s - %(name)s - %(message)s',
-    datefmt='%H:%M:%S'
+    datefmt = '%H:%M:%S'
 )
 
 
@@ -13,4 +13,4 @@ logging.basicConfig(
 logging.getLogger('matplotlib.font_manager').setLevel(logging.WARNING)
 logging.getLogger('matplotlib.ticker').setLevel(logging.WARNING)
 logging.getLogger('matplotlib.pyplot').setLevel(logging.WARNING)
-logging.getLogger('matplotlib.colorbar').setLevel(logging.WARNING)
+logging.getLogger('matplotlib.colorbar').setLevel(logging.WARNING)

nbody/utils/model.py

@@ -11,4 +11,4 @@ def model_density_distribution(r_bins: np.ndarray, M: float = 5, a: float = 5) -
     See https://doi.org/10.1086%2F168845 for more information.
     """
     rho = M / (2 * np.pi) * a / (r_bins * (r_bins + a)**3)
-    return rho
+    return rho