Octree update

CHANGELOG.md

@@ -6,6 +6,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+## [0.5.1] - 2022-09-24
+### Added
+- fixed bugs in octree module
+
 ## [0.5.0] - 2022-09-23
 ### Added
 - writer module

dsmc/octree.py

@@ -1,4 +1,3 @@
-import math
 import numpy as np
 import numpy.typing as npt
 from numba import njit
@@ -65,9 +64,9 @@ def _is_resolved(box : npt.NDArray, N : int, w : float, sigma_T : float, Nmin :
 
 @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]
+    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
 
@@ -99,11 +98,10 @@ def _sort(permutations : npt.NDArray, box : npt.NDArray, positions : npt.NDArray
         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]
+        p = new_permutations[i]
         if _is_inside(positions[p], box):
             _swap(new_permutations, i, runner)
             runner += 1
@@ -114,6 +112,22 @@ def _sort(permutations : npt.NDArray, box : npt.NDArray, positions : npt.NDArray
 @njit
 def get_V(box):
     return (box[0][1] - box[0][0]) * (box[1][1] - box[1][0]) * (box[2][1] - box[2][0])
+
+@njit
+def _create_boxes(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])))
+
+    return [child_geo1, child_geo2, child_geo3, child_geo4, child_geo5, child_geo6, child_geo7, child_geo8]
 
 class Leaf:
     def __init__(self):
@@ -164,49 +178,26 @@ def _create_root(self, positions):
         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
+        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._add_boxes(self.cell_boxes[leaf_id])
+            self.cell_boxes += _create_boxes(self.cell_boxes[leaf_id])
         else:
             pass
-            
-        offset = 0
-            
-        for i in range(leaf.number_children):
+
+        offset = 0 
+
+        for i in range(self.leafs[leaf_id].number_children):
             new_leaf = Leaf()
-            new_leaf.level = leaf.level + 1
+            new_leaf.level = self.leafs[leaf_id].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.elem_offset = self.leafs[leaf_id].elem_offset + offset
+
+            self.permutations, N = _sort(self.permutations, self.cell_boxes[self.leafs[leaf_id].id_first_child + i], positions, new_leaf.elem_offset, self.leafs[leaf_id].number_elements - offset)
 
             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)

examples/temp_relax.py

@@ -1,5 +1,4 @@
 import dsmc
-import dsmc.diagnostics as dia
 import dsmc.particles as prt
 import matplotlib.pyplot as plt
 import numpy as np

setup.py

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

tests/octree/octree_test.py

@@ -22,20 +22,42 @@ def create_particles(N, radius):
 	return positions
 
 if __name__ == "__main__":
-	parser = argparse.ArgumentParser(description='Process some integers.')
-	parser.add_argument('N', type=int)
+    parser = argparse.ArgumentParser(description='Process some integers.')
+    parser.add_argument('N', type=int)
 
-	args = parser.parse_args()
+    args = parser.parse_args()
 
-	radius = 1.0
-	positions = create_particles(args.N, radius)
-	octree = oc.Octree()
-	octree.build(positions)
-	wrt.write_buttom_leafs(octree)
+    radius = 1.0
+    positions = create_particles(args.N, radius)
+    octree = oc.Octree()
+    octree.build(positions)
+    wrt.write_buttom_leafs(octree)
+
+    for i in range(len(octree.leafs)):
+        box = octree.cell_boxes[i]
+        leaf = octree.leafs[i]
+        N = 0
+        for j in range(leaf.elem_offset, leaf.elem_offset + leaf.number_elements):
+            p = octree.permutations[j]
+            pos = positions[p]
+
+            if oc._is_inside(pos, box):
+                N += 1
+            else:
+                raise Exception(pos, box)
+
+        if N != leaf.number_elements:
+            raise Exception(N, leaf.number_elements, box)
 
-	with open("particles.csv", "w") as file:
-		file.write("x, y, z\n")
-		for pos in positions:
-			file.write("{}, {}, {}\n".format(pos[0], pos[1], pos[2]))
+    with open("particles.csv", "w") as file:
+        file.write("x, y, z\n")
+
+        for i in range(len(octree.leafs)):
+            leaf = octree.leafs[i]
+            if leaf.number_children == 0 and leaf.number_elements > 0:
+                for j in range(leaf.elem_offset, leaf.elem_offset + leaf.number_elements):
+                    p = octree.permutations[j]
+                    pos = positions[p]
+                    file.write("{}, {}, {}\n".format(pos[0], pos[1], pos[2]))
 
-	print("done")
+    print("done")

tests/unit/test_data/create_test_particles.py

@@ -0,0 +1,11 @@
+import numpy as np
+
+if __name__ == "__main__":
+    N = 100
+    print("writing {} particles".format(N))
+    with open("particles.csv", "w") as file:
+        for i in range(N):
+            pos = np.random.random(3)*2.0 - np.ones(3)
+            file.write("{},{},{}\n".format(pos[0], pos[1], pos[2]))
+
+    print("done")

tests/unit/test_data/particles.csv

@@ -0,0 +1,100 @@
+0.16720244270111206,0.8208120902385547,-0.15821775757746326
+-0.3539408610094026,0.6094479586781756,-0.2241998042084452
+-0.253602773408937,0.37842838770619047,0.5080805107290469
+0.901801119661219,-0.3547753241315563,-0.7994223817818384
+-0.014540968173746727,-0.6057811852068014,-0.7747150654889452
+0.6652045024170465,0.7496016129090863,0.09269932097207323
+-0.08371757238708022,-0.49569628427213286,-0.5230501995932113
+-0.45628210172061734,-0.7032547855919364,-0.01052801065881459
+-0.4050878274555494,0.16220023223033286,-0.13943218567361249
+-0.44015418512841675,0.4399168343737543,-0.13664372635262478
+0.3822321359090435,-0.03867223327253977,-0.04984480804603453
+-0.5467265015048175,0.7694212897999899,0.5107554092663027
+0.39267358552702203,0.8855795814797593,-0.3874196750436716
+-0.8339217114875397,-0.7277069297889218,0.6932782778950106
+-0.06677305356215668,0.08155865919002991,-0.22499099071269968
+0.7461692085104521,-0.3569676906172021,0.19992850640433013
+-0.3183568330760138,0.8889477004291035,-0.22695163370221572
+0.14226452585971483,0.6749945724336466,-0.10977565633862918
+0.7567553000958556,0.6860849556496309,0.7019398659703175
+-0.45544407366595374,0.8415794987089842,0.4617571900758932
+0.008921032851599175,0.09329133267050382,0.12341660395588416
+-0.11332760346288917,0.4042049964883072,-0.6363589932576024
+0.45801429375353897,-0.9119723962473707,0.4512774411458276
+0.46255470634607354,0.8505216562953395,-0.9625208705796366
+-0.1634911149788587,0.4954169015383527,-0.7514390159518116
+0.9652111916953872,-0.326652317518793,-0.46448051286137493
+-0.012190671492502414,0.3513603108139549,0.21120932266533754
+-0.3219199072922947,-0.6707629224151808,-0.06547558541371457
+0.0003677083723194752,-0.28631549200531037,-0.9589132043625856
+-0.48982911826707265,0.33709944881539533,-0.5705459675253599
+0.24553297910626437,0.8153935258740543,0.7754269872933011
+-0.3597841927312453,-0.30153996145772344,-0.5545584807322563
+0.4514246571306877,0.4022956382825962,-0.49522635923734093
+0.28698107560336505,-0.2816866625438861,0.8615966614249093
+0.09405268228661168,0.3419520966030001,0.7822256998740897
+-0.6883526842679561,-0.716974021117027,0.8038683522503012
+0.8511390374742094,0.898837413751163,0.11732018920271692
+0.10212127918229652,-0.6901138895567802,0.044572452623179215
+0.9742771687163521,0.3900116089700414,0.2849645916426635
+-0.41446601047533616,-0.48501116366070063,-0.21171170708052967
+-0.8985239501161144,0.3231647945259839,0.20482889661142267
+0.2691336467016665,0.019006468602156934,-0.5833077026738216
+0.7753963192748963,0.4078932463220293,-0.17231148296520438
+-0.19696942316303256,0.10696738046497845,0.9229353551086716
+0.5610982192389629,-0.36377857925636947,-0.39267739735571716
+-0.3199431452231085,-0.23776421953278826,-0.7917537970082731
+0.7121864877364648,-0.8806372932737705,0.44221940740157173
+-0.4773400242683048,-0.4846890064033458,0.32609824265601484
+-0.25066582819700156,0.17875327390040763,0.6792664170237515
+-0.34955904993833786,-0.5269472436510836,0.05666029913257997
+-0.9841821357736698,0.4621791724011153,-0.38399251190936035
+-0.7925590286957733,0.8971899467301814,0.9465523534014599
+0.7791035904527936,-0.3065439183432015,0.9668490966958605
+0.7550138781918301,0.9954600710856267,0.15445789110262798
+0.015132235590737064,-0.1279001510759381,-0.7784382777493579
+0.7831998282225763,-0.45035646564466214,0.31822998226315025
+0.25092976269744693,0.9085270898010307,-0.31735361664304573
+0.3805923268337683,0.03787861289930383,-0.6266312154338713
+-0.6040348148625538,0.4530362126423204,0.498118925527959
+0.4136176171187369,-0.9050493502417782,-0.29965017480855916
+-0.42774070010191867,-0.06461591748042528,-0.7976491333955078
+0.09766097044471311,-0.16598582888531488,-0.7608659826062489
+0.33798463823175595,-0.42150973385006063,-0.3713352945819577
+-0.58932187049589,0.3813552375040625,-0.26284015727708554
+-0.4907965754459813,0.20341783229024646,0.6622776161420185
+0.8098477083222237,-0.8719199252159624,0.7488605548345317
+-0.8193929114745129,0.840626931453603,-0.6315912117747844
+-0.24239201418824297,0.47098841467093844,-0.669183593651901
+0.5041947220259142,-0.29544662686416756,-0.5914541788966325
+-0.26262096117978784,0.9801755149383851,0.527269048439474
+0.4275209506977513,-0.667659249339922,0.1347518319999672
+0.9855519426068913,-0.6135765984848576,-0.4287141280807911
+0.20982860758802357,0.07462089616107126,0.5943311631472485
+0.32596254556079907,-0.5523299607926777,0.5000136968687108
+0.6298276916469507,0.7866892509990404,-0.9079439336773338
+0.78142640981854,0.5399877759560965,0.6973205654874475
+-0.3485450193568893,-0.6303033064361121,0.17527538428193168
+-0.8260677718113378,0.40584791574484824,-0.9942842646615264
+-0.20591928046068153,-0.5923212325551002,-0.968679333547988
+0.45672042754779163,0.2149643271371091,0.21071030029661353
+0.057075371525048046,-0.6650594086767447,-0.5251178367633449
+0.1809775156746154,0.4818178472834924,0.3972083715175747
+-0.6207155824174482,-0.4160776878844272,0.6917477102880314
+-0.904590420565192,-0.193648165470109,0.5189271677890341
+0.7951185512115042,-0.6569231762630592,0.3204489978755769
+0.7584068620500117,-0.3176223099304243,-0.870445024417934
+-0.4922227934260852,-0.5001060615332293,-0.5348416435824386
+0.8605684066110679,0.16239313316117032,-0.7627141519458742
+0.6590758735534021,0.9222932137960997,-0.021059658756133137
+-0.5501989000546204,0.26214029221325275,-0.6910176960391732
+0.967198486330149,-0.41117570121974034,0.6945615910511642
+0.9468060139718262,0.5879482989935261,0.08334965744250988
+0.4824650338598835,-0.21178248688885692,-0.9756564237429781
+0.253591146580193,0.08009088613267878,0.4274284592732791
+-0.4389966392317204,0.22920620294282146,-0.0002775036237094852
+-0.7146674493249148,-0.538980218854364,0.879679600158005
+0.8244476665630003,-0.7994155142294497,-0.6359195332795282
+0.6824501498828124,-0.8031118692347363,0.9876378371063732
+0.718260570290431,0.3877012248400389,0.5868224029728766
+-0.7941746749515597,-0.49028054610731187,0.9579007583397083