Dsmc update

CHANGELOG.md

@@ -6,9 +6,17 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
-## [0.5.1] - 2022-09-24
+## [0.6.0] - 2022-09-28
 ### Added
-- fixed bugs in octree module
+- functionality for aspect ratio check in octree
+
+### Fixed
+- bug in temperature initialization
+- several bugs in octree
+
+## [0.5.1] - 2022-09-24
+### Fixed
+- bugs in octree module
 
 ## [0.5.0] - 2022-09-23
 ### Added

dsmc/diagnostics.py

@@ -1,5 +1,4 @@
 import numpy as np
-from numba import njit
 from . import particles as prt
 from . import octree as oc
 
@@ -26,7 +25,6 @@ def sort_bin(positions, axis, Nbin):
 
 def calc_n(bins, box, axis, w):
     Nbins = len(bins)
-    dx = (box[axis][1] - box[axis][0]) / Nbins
     V = 1
     n = np.empty((Nbins, ))
     for i in range(3):

dsmc/dsmc.py

@@ -1,16 +1,18 @@
-import math
 import numpy as np
 from numba import njit
+from numba import prange
 from . import particles as prt
 from . import octree as oc
 
-@njit
-def _push(velocities, positions, dt):  
-    return positions + velocities*dt
+@njit(parallel=True)
+def _push(velocities, positions, dt):
+    for p in prange(len(positions)):
+        positions[p] = positions[p] + velocities[p]*dt
+    return positions
 
-@njit 
+@njit(parallel=True)
 def _boundary(velocities, positions, domain):
-    for p in range(len(positions)):
+    for p in prange(len(positions)):
         while not oc._is_inside(positions[p], domain):
             for i in range(3):
                 if positions[p][i] < domain[i][0]:
@@ -19,11 +21,11 @@ def _boundary(velocities, positions, domain):
                 if positions[p][i] > domain[i][1]:
                     positions[p][i] = 2.0 * domain[i][1] - positions[p][i]
                     velocities[p][i] *= -1.0
-                    
+
     return (velocities, positions)
- 
-@njit   
-def _calc_prob(vel1 : np.ndarray, vel2 : np.ndarray, sigma_T : float, Vc : float, dt : float, w : float, N : int) -> np.single:
+
+@njit
+def _calc_prob(rel_vel : float, sigma_T : float, Vc : float, dt : float, w : float, N : int) -> np.single:
     """
     Parameters
     ----------
@@ -37,14 +39,14 @@ def _calc_prob(vel1 : np.ndarray, vel2 : np.ndarray, sigma_T : float, Vc : float
         weight
     N : int
         number of particles
-        
+
     Returns
     -------
     collision proability : float
     """
-    return np.linalg.norm(vel1 - vel2) * sigma_T * dt * w * N / Vc;
-  
-@njit  
+    return rel_vel * sigma_T * dt * w * N / Vc;
+
+@njit
 def _calc_post_col_vels(velocity1 : np.ndarray, velocity2 : np.ndarray, mass1 : float, mass2 : float, rel_vel_module : float, rand_number1 : float, rand_number2 : float) -> tuple[np.ndarray, np.ndarray]:
     mass12 = (mass1 + mass2)
     mass1_12 = (mass1 / mass12)
@@ -55,100 +57,99 @@ def _calc_post_col_vels(velocity1 : np.ndarray, velocity2 : np.ndarray, mass1 :
     epsilon = (2.0 * np.pi * rand_number2)
 
     centre_of_mass_velocity = (velocity1 * mass1 + velocity2 * mass2) * (1.0 / mass12)
-    
+
     rel_velocity_new = np.empty((3, ))
 
     rel_velocity_new[0] = rel_vel_module * cos_xi
     rel_velocity_new[1] = rel_vel_module * sin_xi * np.cos(epsilon)
     rel_velocity_new[2] = rel_vel_module * sin_xi * np.sin(epsilon)
 
     return (centre_of_mass_velocity + rel_velocity_new * mass2_12 , centre_of_mass_velocity - rel_velocity_new * mass1_12)
-   
-@njit 
+
+@njit
 def _update_velocities(permutations : np.ndarray, velocities : np.ndarray, mass : float, sigma_T : float, Vc : float, dt : float, w : float, offset : int, N : int) -> np.ndarray:
-    i = 1
-    while i < N:
+    for i in range(1, N, 2):
         p1 = permutations[offset + i - 1]
         p2 = permutations[offset + i]
-        P = _calc_prob(velocities[p1], velocities[p2], sigma_T, Vc, dt, w, N)
+        rel_vel = np.linalg.norm(velocities[p1] - velocities[p2])
+        P = _calc_prob(rel_vel, sigma_T, Vc, dt, w, N)
         R = np.random.random(3)
-        
+
         if R[0] < P:
-            new_vels = _calc_post_col_vels(velocities[p1], velocities[p2], mass, mass, np.linalg.norm(velocities[p1] - velocities[p2]), R[1], R[2])
+            new_vels = _calc_post_col_vels(velocities[p1], velocities[p2], mass, mass, rel_vel, R[1], R[2])
             velocities[p1] = new_vels[0]
             velocities[p2] = new_vels[1]
-
-        i += 2
-
+
     return velocities
- 
-@njit   
+
+@njit
 def _update_vels(permutations : np.ndarray, velocities : np.ndarray, mass : float, sigma_T : float, dt : float, w : float, elem_offsets : np.ndarray, number_elements : np.ndarray, number_children : np.ndarray, cell_boxes : np.ndarray, Nleafs : int) -> np.ndarray:
     for i in range(Nleafs):
-        if not number_children[i]:
+        if not number_children[i] and number_elements[i]:
             Vc = oc.get_V(cell_boxes[i])
             velocities = _update_velocities(permutations, velocities, mass, sigma_T, Vc, dt, w, elem_offsets[i], number_elements[i])
-            
+
     return velocities
 
 class DSMC:
     def __init__(self):
         self.clear()
-        
+
     def clear(self):
         self.particles = prt.Particles()
         self.octree = oc.Octree()
         self.w = None
         self.domain = None
         self.sigma_T = 3.631681e-19
         self.mass = None
-        
-    def advance(self, dt):
+
+    def advance(self, dt, collisions=True, octree=True):
         if self.domain is None:
             raise Exception("simulation domain not defined")
         if self.particles.N == 0:
             raise Exception("no particles created")
         if self.w == None:
             raise Exception("particle weight not set")
-
-        self.octree.build(self.particles.Pos)
-
-        self.particles.VelPos = (self._update_velocities(dt), self.particles.Pos)
+
+        if octree:
+            self.octree.build(self.particles.Pos)
+        if collisions and octree:
+            self.particles.VelPos = (self._update_velocities(dt), self.particles.Pos)
         positions = _push(self.particles.Vel, self.particles.Pos, dt)
         self.particles.VelPos = _boundary(self.particles.Vel, positions, self.domain)
-        
+
     def _update_velocities(self, dt):
         Nleafs : int = len(self.octree.leafs)
         elem_offsets : np.ndarray = np.array([leaf.elem_offset for leaf in self.octree.leafs], dtype=int)
         number_elements : np.ndarray = np.array([leaf.number_elements for leaf in self.octree.leafs], dtype=int)
         number_children : np.ndarray = np.array([leaf.number_children for leaf in self.octree.leafs], dtype=int)
         cell_boxes : np.ndarray = np.array([box for box in self.octree.cell_boxes])
-        
+
         return _update_vels(self.octree.permutations, self.particles.Vel, self.mass, self.sigma_T, dt, self.w, elem_offsets, number_elements, number_children, cell_boxes, Nleafs)
-        
+
     def create_particles(self, box, T, n, u = None):
         box = np.array(box)
         N = int(round(oc.get_V(box) * n / self.w))
         print("creating {} particles".format(N))
         self.particles.create_particles(box, self.mass, T, N)
-        
+
         if u is not None:
             velocities = self.particles.Vel
             u = np.array(u)
-            
+
             for i in range(len(velocities)):
                 velocities[i] += u
-            
+
             self.particles.VelPos = (velocities, self.particles.Pos)
-        
+
         print("now containing {} particles, {} total".format(N, self.particles.N))
-        
+
     def set_domain(self, domain):
         self.domain = np.array(domain)
-        
+
     def set_mass(self, mass):
         self.mass = mass
-        
+
     def set_weight(self, w):
         self.octree.w = w
         self.w = w

dsmc/octree.py

@@ -1,6 +1,7 @@
 import numpy as np
 import numpy.typing as npt
 from numba import njit
+import numba as nb
 
 fmin = np.finfo(float).min
 fmax = np.finfo(float).max
@@ -128,6 +129,58 @@ def _create_boxes(box):
     child_geo8 = np.array(((half[0], box[0][1]), (box[1][0], half[1]), (box[2][0], half[2])))
 
     return [child_geo1, child_geo2, child_geo3, child_geo4, child_geo5, child_geo6, child_geo7, child_geo8]
+
+@njit
+def _get_min_aspect_ratio(box, axis):
+    half = np.array([0.5*(box[i][1] - box[i][0]) for i in range(3)])
+
+    match axis:
+        case 0:
+            return min(half[0] / half[1], half[0] / half[2]);
+        case 1:
+            return min(half[1] / half[0], half[1] / half[2]);
+        case 2:
+            return min(half[2] / half[1], half[2] / half[0]);
+
+@njit
+def _devide(box, axis):
+    half = 0.5*(box[axis][0] + box[axis][1])
+    box1 = np.copy(box)
+    box2 = np.copy(box)
+
+    box1[axis][0] = box[axis][0]
+    box1[axis][1] = half
+
+    box2[axis][0] = half
+    box2[axis][1] = box[axis][1]
+
+    return (box1, box2)
+
+@njit
+def _create_combined_boxes(box, min_aspect_ratio):
+    boxes = np.empty((15, 3, 2))
+    boxes[0] = box
+    N = 0
+    Nold = 0
+    q = 1
+
+    for i in range(3):
+        if _get_min_aspect_ratio(box, i) > min_aspect_ratio:
+            for b in range(Nold, Nold + 2**N):
+                new_boxes = _devide(boxes[b], i)
+                boxes[q] = new_boxes[0]
+                boxes[q + 1] = new_boxes[1]
+                q += 2
+            Nold += 2**N
+            N += 1
+
+    N = 2**N
+    new_boxes = np.empty((N, 3, 2))
+
+    for b in range(N):
+        new_boxes[b] = boxes[Nold + b]
+
+    return new_boxes
 
 class Leaf:
     def __init__(self):
@@ -141,6 +194,7 @@ def __init__(self):
 class Octree:
     def __init__(self):
         self.clear()
+        self.min_aspect_ratio = 2.0/3.0
 
     def clear(self):
         self.cell_boxes = []
@@ -179,10 +233,17 @@ def _create_root(self, positions):
 
     def _progress(self, leaf_id, positions):
         if _is_resolved(self.cell_boxes[leaf_id], self.leafs[leaf_id].number_elements, self.w, self.sigma_T, self.Nmin, self.Nmax):
-            self.leafs[leaf_id].number_children = 8
+
             self.leafs[leaf_id].id_first_child = self.cell_offsets[-1]
-            self.cell_offsets[-1] += 8
-            self.cell_boxes += _create_boxes(self.cell_boxes[leaf_id])
+
+            new_boxes = _create_combined_boxes(self.cell_boxes[leaf_id], self.min_aspect_ratio)
+            self.cell_offsets[-1] += len(new_boxes)
+            self.leafs[leaf_id].number_children = len(new_boxes)
+
+            for box in new_boxes:
+                self.cell_boxes.append(box)
+
+            #raise Exception()
         else:
             pass
 

dsmc/particles.py

@@ -37,7 +37,7 @@ def x2velocity(x, mass):
     -------
     speed of particle : float
     """
-    return math.sqrt((2.0/3.0) * x * kb /mass)
+    return math.sqrt(2.0 * x * kb /mass)
 
 @njit
 def get_vel(T, mass):
@@ -53,7 +53,8 @@ def get_vel(T, mass):
     -------
     velocity : np.array, shape = (3, 1)
     """
-    return np.array([(-1)**(int(2*np.random.random())) * x2velocity(box_muller(T), mass) for _ in range(3)])
+    v = np.random.random(3)
+    return v * x2velocity(box_muller(T), mass) / np.linalg.norm(v)
 
 @njit
 def get_velocities(T, mass, N):

dsmc/writer.py

@@ -1,7 +1,5 @@
-from . import octree as oc
-
-def write_buttom_leafs(octree):
-	f = open("octree.vtu", "w")
+def write_buttom_leafs(octree, file_name="octree.vtu"):
+	f = open(file_name, "w")
 
 	_write_header(f)
 	_wrtie_body(f, octree)

examples/T_sample.py

@@ -0,0 +1,37 @@
+import dsmc
+import dsmc.particles as prt
+import matplotlib.pyplot as plt
+import numpy as np
+import math
+
+def maxwell(x, T):
+	return 2.0 * np.sqrt(x) * np.exp(-x/T) / (math.pow(T, 3.0/2.0) * np.sqrt(math.pi))
+
+def calc_x(velocities, mass):
+	return np.array([mass*vel.dot(vel)/(2.0*prt.kb) for vel in velocities])
+
+if __name__ == '__main__':
+	solver = dsmc.DSMC()
+	domain = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (0, 50e-3))
+	w = 1e6
+	mass = 6.6422e-26
+	T = 300
+	n = 1e+20
+	Nbins = 100
+
+	solver.set_domain(domain)
+	solver.set_weight(w)
+	solver.set_mass(mass)
+
+	solver.create_particles(domain, T, n)
+
+	print(prt.calc_temperature(solver.particles.Vel, mass))
+
+	x = calc_x(solver.particles.Vel, solver.mass)
+
+	xm = np.linspace(0, 3500, 1000)
+	dist = [maxwell(xi, 300) for xi in xm]
+
+	plt.plot(xm, dist)
+	plt.hist(x, Nbins, density=True)
+	plt.show()