diff --git a/CHANGELOG.md b/CHANGELOG.md index 99aa5e2..6236035 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/dsmc/octree.py b/dsmc/octree.py index 09b3c33..f1dd8ad 100644 --- a/dsmc/octree.py +++ b/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) diff --git a/examples/temp_relax.py b/examples/temp_relax.py index 8961579..bed76ac 100644 --- a/examples/temp_relax.py +++ b/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 diff --git a/setup.py b/setup.py index ea7e25e..005ff89 100644 --- a/setup.py +++ b/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"], diff --git a/tests/octree/octree_test.py b/tests/octree/octree_test.py index 692f661..3856ff1 100644 --- a/tests/octree/octree_test.py +++ b/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") diff --git a/tests/unit/test_data/create_test_particles.py b/tests/unit/test_data/create_test_particles.py new file mode 100644 index 0000000..9025d0b --- /dev/null +++ b/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") \ No newline at end of file diff --git a/tests/unit/test_data/particles.csv b/tests/unit/test_data/particles.csv new file mode 100644 index 0000000..c28b5c6 --- /dev/null +++ b/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 diff --git a/tests/unit/test_dsmc/octree.py b/tests/unit/test_dsmc/octree.py index 1ad85f1..56a1874 100644 --- a/tests/unit/test_dsmc/octree.py +++ b/tests/unit/test_dsmc/octree.py @@ -2,6 +2,7 @@ import numpy as np from dsmc import octree as oc from dsmc import particles as part +import csv class TestOctree(unittest.TestCase): def test__find_bounding_box(self): @@ -85,6 +86,140 @@ def test_sort(self): for i in range(offset + Nnew, len(permutations)): p = permutations[i] self.assertFalse(oc._is_inside(positions[p], box)) + + def test_sort2(self): + box = np.array([[0.0, 1.0], [0.0, 1.0], [0.0, 1.0]]) + N = 20 + Nh = 10 + Nm = 15 + Np = 5 + positions = np.empty((N, 3)) + positions1 = np.empty((Np, 3)) + positions2 = np.empty((Np, 3)) + positions3 = np.empty((Nh, 3)) + permutations = np.array([i for i in range(N)]) + + for i in range(Np): + positions1[i] = np.random.random(3) - np.ones(3) + + for i in range(Np): + positions2[i] = np.random.random(3) + + for i in range(Nh): + positions3[i] = np.random.random(3) - np.ones(3) + + positions = np.concatenate((positions3, positions1, positions2)) + permutations, Nnew = oc._sort(permutations, box, positions, Nh, Nh) + + self.assertEqual(Nnew, Np) + + for i in range(Nh): + p = permutations[i] + pos = positions[p] + a = not (pos[0] <= box[0][1]) + b = not (pos[0] >= box[0][0]) + + c = not (pos[1] <= box[1][1]) + d = not (pos[1] >= box[1][0]) + + e = not (pos[2] <= box[2][1]) + f = not (pos[2] >= box[2][0]) + + self.assertTrue(a or b or c or d or e or f) + + for i in range(Nm, N): + p = permutations[i] + pos = positions[p] + a = not (pos[0] <= box[0][1]) + b = not (pos[0] >= box[0][0]) + + c = not (pos[1] <= box[1][1]) + d = not (pos[1] >= box[1][0]) + + e = not (pos[2] <= box[2][1]) + f = not (pos[2] >= box[2][0]) + + self.assertTrue(a or b or c or d or e or f) + + for i in range(Nh, Nm): + p = permutations[i] + pos = positions[p] + a = (pos[0] <= box[0][1]) + b = (pos[0] >= box[0][0]) + + c = (pos[1] <= box[1][1]) + d = (pos[1] >= box[1][0]) + + e = (pos[2] <= box[2][1]) + f = (pos[2] >= box[2][0]) + + self.assertTrue(a and b and c and d and e and f) + + + def load_particles(self): + positions = [] + + with open("./test_data/particles.csv", "r", newline='') as csvfile: + reader = csv.reader(csvfile, delimiter=',') + for row in reader: + positions.append([float(row[i]) for i in range(3)]) + + return np.array(positions) + + def test_sort3(self): + box = np.array([[-1.0, 1.0], [-1.0, 1.0], [-1.0, 1.0]]) + boxes = [] + oc._create_boxes(box) + N = 100 + permutations = np.array([i for i in range(N)]) + positions = self.load_particles() + + permutations , Np1 = oc._sort(permutations, boxes[0], positions, 0, N) + permutations , Np2 = oc._sort(permutations, boxes[1], positions, Np1, N - Np1) + permutations , Np3 = oc._sort(permutations, boxes[2], positions, Np1 + Np2, N - Np2 - Np1) + #permutations , Np4 = oc._sort(permutations, boxes[3], positions, Np3, N - Np3) + + #permutations , Np5 = oc._sort(permutations, boxes[4], positions, Np4, N - Np4) + #permutations , Np6 = oc._sort(permutations, boxes[5], positions, Np5, N - Np5) + #permutations , Np7 = oc._sort(permutations, boxes[6], positions, Np6, N - Np6) + #permutations , Np8 = oc._sort(permutations, boxes[7], positions, Np7, N - Np7) + + Nnew = np.zeros(8, dtype=int) + + for i in range(Np1): + p = permutations[i] + pos = positions[p] + if oc._is_inside(pos, boxes[0]): + Nnew[0] += 1 + + self.assertEqual(Nnew[0], Np1) + + for i in range(Np1, Np1 + Np2): + p = permutations[i] + pos = positions[p] + if oc._is_inside(pos, boxes[1]): + Nnew[1] += 1 + + self.assertEqual(Nnew[1], Np2) + + for i in range(Np1 + Np2, Np1 + Np2 + Np3): + p = permutations[i] + pos = positions[p] + if oc._is_inside(pos, boxes[2]): + Nnew[2] += 1 + + self.assertEqual(Nnew[2], Np3) + + def test__create_boxes(self): + box_orig = np.array([[-1.0, 1.0], [-1.0, 1.0], [-1.0, 1.0]]) + boxes = oc._create_boxes(box_orig) + + self.assertEqual(8, len(boxes)) + V = 0.0 + + for box in boxes: + V += oc.get_V(box) + + self.assertEqual(oc.get_V(box_orig), V) class TestOctreeOctree(unittest.TestCase): @@ -94,11 +229,3 @@ def test_build(self): octree.w = 1e+18 octree.build(positions) - - def test__add_boxes(self): - octree = oc.Octree() - box = np.array([[-1.0, 1.0], [-1.0, 1.0], [-1.0, 1.0]]) - - octree._add_boxes(box) - - self.assertEqual(8, len(octree.cell_boxes))