Particle Module

CHANGELOG.md

@@ -6,6 +6,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+## [0.2.0] - 2022-09-15
+### Added
+- particle module
+
 ## [0.1.0] - 2022-09-15
 ### Added
 - unit test environment

dsmc/particles.py

@@ -0,0 +1,134 @@
+import math
+import numpy as np
+from numba import njit
+
+kb = 1.380649e-23
+
+@njit
+def box_muller(T):
+    """
+    Parameters
+    ----------
+    T : float
+        temperature [K]
+
+    Returns
+    -------
+    x : float
+        m * v^2 / (2 * kb)
+    """
+    r1 = np.random.random()
+    r2 = np.random.random()
+    r3 = np.random.random()
+
+    return T*(-np.log(r1) - np.log(r2) * math.pow(np.cos((np.pi/2.0)*r3), 2))
+
+@njit
+def x2velocity(x, mass):
+    """
+    Parameters
+    ----------
+    x : float
+        m * v^2 / (2 * kb)
+    mass : float
+           particle mass [kg]
+
+    Returns
+    -------
+    speed of particle : float
+    """
+    return math.sqrt((2.0/3.0) * x * kb /mass)
+
+@njit
+def get_vel(T, mass):
+    """
+    Parameters
+    ----------
+    T : float
+        temperature [K]
+    mass : float
+           particle mass [kg]
+
+    Returns
+    -------
+    velocity : np.array, shape = (3, 1)
+    """
+    return np.array([(-1)**(int(2*np.random.random())) * x2velocity(box_muller(T), mass) for _ in range(3)])
+
+@njit
+def get_velocities(T, mass, N):
+    velocities = np.empty((N, 3), dtype=float)
+
+    for i in range(N):
+        velocities[i] = get_vel(T, mass)
+
+    return velocities
+
+@njit
+def calc_temperature(velocities, mass):
+    tot_e = 0.0
+
+    if not len(velocities):
+        return 0.0
+
+    for i in range(len(velocities)):
+        tot_e += 0.5 * mass * np.dot(velocities[i], velocities[i])
+
+    return tot_e / ((3.0/2.0) * len(velocities) * kb)
+
+@njit
+def calc_positions(x, y, z, N):
+    """
+    Parameters
+    ----------
+    x : tuple(2), dtype = float
+        xmin, xmax
+    y : tuple(2), dtype = float
+        ymin, ymax
+    z : tuple(2), dtype = float
+        zmin, zmax
+    """
+    positions = np.empty((N, 3), dtype=float)
+
+    for i in range(N):
+        positions[i][0] = x[0] + np.random.random() * (x[1] - x[0])
+        positions[i][1] = y[0] + np.random.random() * (y[1] - y[0])
+        positions[i][2] = z[0] + np.random.random() * (z[1] - z[0])
+
+    return positions
+
+class Particles:
+    def __init__(self):
+        self._velocities = None
+        self._positions = None
+        self._N = 0 # number of particles
+
+    @property
+    def Vel(self):
+        return self._velocities
+
+    @property
+    def Pos(self):
+        return self._positions
+
+    @property
+    def N(self):
+        return self._N
+
+    @property
+    def VelPos(self):
+        return (self._velocities, self._positions)
+
+    @VelPos.setter
+    def VelPos(self, vel_pos):
+        assert vel_pos[0] == vel_pos[1]
+        assert isinstance(vel_pos[0], np.array)
+        assert isinstance(vel_pos[1], np.array)
+        self._velocities = vel_pos[0]
+        self._positions = vel_pos[1]
+        self._N = len(self._positions)
+
+    def create_particles(self, X, mass, T, N):
+        self._velocities = get_velocities(T, mass, N)
+        self._positions = calc_positions(X[0], X[1], X[2], N)
+        self._N  = N

tests/unit/main.py

@@ -3,6 +3,7 @@
 import unittest
 
 from test_dsmc.dsmc import *
+from test_dsmc.particles import *
 
 if __name__ == '__main__':
 	unittest.main()

tests/unit/test_dsmc/particles.py

@@ -0,0 +1,65 @@
+import unittest
+from dsmc import particles as pa
+
+class TestParticles(unittest.TestCase):
+    def test_box_muller(self):
+        T = 300
+        result = pa.box_muller(T)
+
+    def test_x2velocity(self):
+        x = 1.0
+        mass = 1.0e-26
+        result = pa.x2velocity(x, mass)
+
+    def test_get_vel(self):
+        T = 300
+        mass = 1.0e-26
+        result = pa.get_vel(T, mass)
+
+    def test_get_velocities(self):
+        T = 300
+        mass = 1.0e-26
+        N = 1000
+        velocities = pa.get_velocities(T, mass, N)
+
+        self.assertEqual(N, len(velocities))
+
+    def test_calc_temperature(self):
+        T = 300
+        mass = 1.0e-26
+        N = 10000
+        velocities = pa.get_velocities(T, mass, N)
+        T_new = pa.calc_temperature(velocities, mass)
+        diff = abs((T_new - T)/T)
+
+        self.assertTrue(diff < 0.1)
+
+    def test_calc_positions(self):
+        x = (-1.0, 1.0)
+        y = (2.0, 3.0)
+        z = (-2.0, -1.0)
+        N = 1000
+        positions = pa.calc_positions(x, y, z, N)
+
+        self.assertEqual(N, len(positions))
+
+    def test_create_particles(self):
+        x = (-1.0, 1.0)
+        y = (2.0, 3.0)
+        z = (-2.0, -1.0)
+        X = (x, y, z)
+        N = 1000
+        mass = 1.0e-26
+        T = 300
+        particles = pa.Particles()
+
+        particles.create_particles(X, mass, T, N)
+
+        self.assertEqual(N, len(particles.Pos))
+        self.assertEqual(N, len(particles.Vel))
+        self.assertEqual(N, particles.N)
+
+        for i in range(N):
+            self.assertTrue(particles.Pos[i][0] >= x[0] and particles.Pos[i][0] <= x[1])
+            self.assertTrue(particles.Pos[i][1] >= y[0] and particles.Pos[i][1] <= y[1])
+            self.assertTrue(particles.Pos[i][2] >= z[0] and particles.Pos[i][2] <= z[1])