Octree

CHANGELOG.md

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

doc/octree_resolution.ipynb

@@ -0,0 +1,68 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "9db716f6-c301-4d2e-b06f-19cd98de2e5f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import math"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "772d42fa-dc37-48e9-a599-392077582e98",
+   "metadata": {},
+   "source": [
+    "The mean free path can be written as\n",
+    "$$\n",
+    "\\lambda = \\frac{1}{\\sqrt{2} \\cdot n_{c} \\cdot \\sigma_T}\n",
+    "$$\n",
+    "where $n_c$ is the number density in the cell and $\\sigma_T$ is the total cross section.\n",
+    "The number of particles in cell $N_{c}$ can be calculated as (asuming cubic cells)\n",
+    "$$\n",
+    "N_{c} = \\frac{1}{w} \\cdot L_c^3 \\cdot n_c\n",
+    "$$\n",
+    "where $L_C$ is the side length of the cell and $w$ being the particle weight.\n",
+    "Assuming that the cell side fullfills the criteria\n",
+    "$$\n",
+    "\\frac{L_c}{\\lambda} \\leq \\frac{1}{2}.\n",
+    "$$\n",
+    "This leads to the criterium\n",
+    "$$\n",
+    "N \\leq \\frac{\\sqrt{2}}{32 \\cdot w \\cdot \\sigma^2_T \\cdot n^3_c}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4f03b0f7-64a9-48e5-989c-2e253ce09c1d",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "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.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

dsmc/octree.py

@@ -0,0 +1,207 @@
+import math
+import numpy as np
+import numpy.typing as npt
+from numba import njit
+
+fmin = np.finfo(float).min
+fmax = np.finfo(float).max
+
+@njit
+def _find_bounding_box(positions : npt.NDArray) -> npt.NDArray:
+    box = np.array([[fmax, fmin], [fmax, fmin], [fmax, fmin]])
+
+    for pos in positions:
+        for i in range(3):
+            if pos[i] < box[i][0]:
+                box[i][0] = pos[i]
+            if pos[i] > box[i][1]:
+                box[i][1] = pos[i]
+
+    return box
+
+@njit    
+def _calc_N_res(w : float, sigma_T : float, n : float) -> int:
+    """
+    Parameters
+    ----------
+    w : float
+        particle weight
+    sigma_t : float
+              total cross section [m^2]
+    n : float
+        number density [1/m^3]
+    """
+
+    return int(round(np.sqrt(2.0) / (32.0 * w * sigma_T**3 * n**2)))
+
+@njit  
+def _calc_n(box : npt.NDArray, N : float, w : float) -> float:
+    """Calculates number density in cell
+
+    Parameters
+    ----------
+    box : np.array(3, 3)
+          cell
+    N : int
+        number of particles in cell
+    w : float
+        particle weight
+
+    Returns
+    -------
+    number density : float
+    """
+    return np.prod(np.array([box[i][1] - box[i][0] for i in range(3)])) * N / w
+
+@njit  
+def _is_resolved(box : npt.NDArray, N : int, w : float, sigma_T : float, Nmin : int, Nmax : int) -> bool:
+    if N == 0:
+        return False
+
+    n = _calc_n(box, N, w)
+    Nres = _calc_N_res(w, sigma_T, n)
+
+    return N > 2 * min(Nmin, max(Nmin, Nres))
+
+@njit
+def _is_inside(position : npt.NDArray, box : npt.NDArray) -> bool:
+    a : bool = position[0] > box[0][0] and position[0] <= box[0][1]
+    b : bool = position[1] > box[1][0] and position[1] <= box[1][1]
+    c : bool = position[2] > box[2][0] and position[2] <= box[2][1]
+
+    return a and b and c
+
+@njit    
+def _swap(arr, pos1, pos2):
+    temp = arr[pos2]
+    arr[pos2] = arr[pos1]
+    arr[pos1] = temp
+
+@njit
+def _sort(permutations : npt.NDArray, box : npt.NDArray, positions : npt.NDArray, offset : int, N : int) -> tuple[npt.NDArray, int]:
+    '''sort particles in cell
+
+    Parameters
+    ----------
+    box : np.ndarray
+          cell
+    positions : ndarray((3, ))
+        particle positions
+    offset : int
+        number offset
+    N : int
+        number of particles to be considered
+
+    Returns
+    -------
+    new_permutations : ndarray[int]
+    N : int
+        number of found positions
+    '''
+    new_permutations = np.copy(permutations)
+    temp = np.empty((3,))
+    runner = offset
+    Nnew = 0
+    for i in range(offset, offset + N):
+        p = permutations[i]
+        if _is_inside(positions[p], box):
+            _swap(new_permutations, i, runner)
+            runner += 1
+            Nnew += 1
+
+    return new_permutations, Nnew
+
+class Leaf:
+    def __init__(self):
+        self.level = 0
+        self.elem_offset = 0
+        self.number_elements = 0
+        self.id_parent = None
+        self.id_first_child = None
+        self.number_children = 0
+
+class Octree:
+    def __init__(self):
+        self.clear()
+
+    def clear(self):
+        self.cell_boxes = []
+        self.leafs = []
+        self.sigma_T = 3.631681e-19
+        self.w = 1.0
+        self.Nmin = 8
+        self.Nmax = 64
+        self.max_level = 10
+        self.permutations = []
+        self.cell_offsets = []
+        self.level = 0
+
+    def build(self, positions):
+        self._create_root(positions)
+        self.permutations = np.array([i for i in range(len(positions))])
+
+        for level in range(self.max_level):
+            self.level += 1
+            self.cell_offsets.append(self.cell_offsets[-1])
+            for i in range(self.cell_offsets[level], self.cell_offsets[level + 1]):
+                self._progress(i, positions)
+
+            if self.cell_offsets[level + 1] == self.cell_offsets[level + 2]:
+                break            
+
+    def _create_root(self, positions):
+        box = _find_bounding_box(positions)
+        leaf = Leaf()
+        leaf.number_elements = len(positions)
+
+        self.cell_offsets += [0, 1]
+        self.leafs.append(leaf)
+        self.cell_boxes.append(box)
+
+    def _progress(self, leaf_id, positions):
+        leaf = self.leafs[leaf_id]
+        if _is_resolved(self.cell_boxes[leaf_id], leaf.number_elements, self.w, self.sigma_T, self.Nmin, self.Nmax):
+            leaf.number_children = 8
+            leaf.id_first_child = leaf_id + 1
+            self.cell_offsets[-1] += 8
+            self._add_boxes(self.cell_boxes[leaf_id])
+        else:
+            pass
+
+        offset = 0
+
+        for i in range(leaf.number_children):
+            new_leaf = Leaf()
+            new_leaf.level = leaf.level + 1
+            new_leaf.id_parent = leaf_id
+
+            self.permutations, N = _sort(self.permutations, self.cell_boxes[leaf_id + 1 + i], positions, leaf.elem_offset, leaf.number_elements)
+
+            new_leaf.number_elements = N
+            new_leaf.elem_offset = leaf.elem_offset + offset
+            offset += N
+
+            self.leafs.append(new_leaf)
+
+    def _add_boxes(self, box):
+        half = np.array([0.5*(box[i][0] + box[i][1]) for i in range(3)])
+
+        child_geo1 = np.array(((half[0], box[0][1]), (half[1], box[1][1]), (half[2], box[2][1])))
+        child_geo2 = np.array(((box[0][0], half[0]), (half[1], box[1][1]), (half[2], box[2][1])))
+        child_geo3 = np.array(((box[0][0], half[0]), (box[1][0], half[1]), (half[2], box[2][1])))
+        child_geo4 = np.array(((half[0], box[0][1]), (box[1][0], half[1]), (half[2], box[2][1])))
+
+        child_geo5 = np.array(((half[0], box[0][1]), (half[1], box[1][1]), (box[2][0], half[2])))
+        child_geo6 = np.array(((box[0][0], half[0]), (half[1], box[1][1]), (box[2][0], half[2])))
+        child_geo7 = np.array(((box[0][0], half[0]), (box[1][0], half[1]), (box[2][0], half[2])))
+        child_geo8 = np.array(((half[0], box[0][1]), (box[1][0], half[1]), (box[2][0], half[2])))
+
+        self.cell_boxes.append(child_geo1)
+        self.cell_boxes.append(child_geo2)
+        self.cell_boxes.append(child_geo3)
+        self.cell_boxes.append(child_geo4)
+
+        self.cell_boxes.append(child_geo5)
+        self.cell_boxes.append(child_geo6)
+        self.cell_boxes.append(child_geo7)
+        self.cell_boxes.append(child_geo8)

setup.py

@@ -1,7 +1,7 @@
 from setuptools import setup
 setup(
     name='dsmc',
-    version='0.1.0',
+    version='0.3.0',
     author='Leo Basov',
     python_requires='>=3.6, <4',
     packages=["dsmc"],

tests/unit/main.py

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