From 06ab7b2a9139c1a096b739d16c37620b42fb59fc Mon Sep 17 00:00:00 2001 From: LeoBasov <40352679+LeoBasov@users.noreply.github.com> Date: Mon, 3 Oct 2022 15:55:32 +0200 Subject: [PATCH 01/11] implemented _sort_2d' --- dsmc/diagnostics.py | 29 ++++++++++++++++++++++++++ tests/diagnostics/dia_test.py | 39 +++++++++++++++++++++++++++++++++++ 2 files changed, 68 insertions(+) create mode 100644 tests/diagnostics/dia_test.py diff --git a/dsmc/diagnostics.py b/dsmc/diagnostics.py index 56b0b7f..73f0c86 100644 --- a/dsmc/diagnostics.py +++ b/dsmc/diagnostics.py @@ -2,6 +2,35 @@ from . import particles as prt from . import octree as oc +def _sort_2d(positions, x0, x1, y0, y1, Nx, Ny): + dx = (x1 - x0) / Nx + dy = (y1 - y0) / Ny + sorted_vals = np.zeros((Nx, Ny)) + + for pos in positions: + x = x0 + y = y0 + brk = False + + for i in range(Nx): + if brk: + break + y = y0 + + for j in range(Ny): + print(x, y) + a = pos[0] >= x and pos[0] <= (x + dx) + b = pos[1] >= y and pos[1] <= (y + dy) + if a and b: + sorted_vals[i][j] += 1 + brk = True + break + y += dy + x += dx + + return sorted_vals + + def sort_bin(positions, axis, Nbin): bins = [[] for _ in range(Nbin)] b = 0 diff --git a/tests/diagnostics/dia_test.py b/tests/diagnostics/dia_test.py new file mode 100644 index 0000000..f4ce5ae --- /dev/null +++ b/tests/diagnostics/dia_test.py @@ -0,0 +1,39 @@ +import numpy as np +import dsmc.diagnostics as dia + +def create_particles(N, radius): + positions = np.empty((N + 1, 3)) + + for i in range(N): + phi = np.random.random() * 2.0 * np.pi + theta = np.random.random() * np.pi + r = np.random.normal(0.0, 0.01) + theta1 = np.random.random() * np.pi - 0.5 * np.pi + dis = np.array((r * np.sin(theta) * np.cos(phi), r * np.sin(theta) * np.sin(phi), r * np.cos(theta))) + offset = np.array((radius * np.sin(theta1), radius * np.cos(theta1), 0.0)) + + dis += offset; + positions[i] = dis + + positions[N] = np.array((0.0, -1.0, 0.0)) + + return positions + +if __name__ == "__main__": + N = 1000 + radius = 1.0 + positions = create_particles(N, radius) + + # set up sort + Nx = 10 + Ny = 10 + x0 = -1.5 + x1 = 1.5 + y0 = -1.5 + y1 = 1.5 + + sorted_vals = dia._sort_2d(positions, x0, x1, y0, y1, Nx, Ny) + + print(sorted_vals) + + print("done") \ No newline at end of file From ea7923ba63362ac134ed9460f074465996805e53 Mon Sep 17 00:00:00 2001 From: LeoBasov <40352679+LeoBasov@users.noreply.github.com> Date: Mon, 3 Oct 2022 16:16:10 +0200 Subject: [PATCH 02/11] updated sort --- dsmc/diagnostics.py | 8 ++++---- tests/diagnostics/dia_test.py | 7 ++++--- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/dsmc/diagnostics.py b/dsmc/diagnostics.py index 73f0c86..83fce44 100644 --- a/dsmc/diagnostics.py +++ b/dsmc/diagnostics.py @@ -5,9 +5,10 @@ def _sort_2d(positions, x0, x1, y0, y1, Nx, Ny): dx = (x1 - x0) / Nx dy = (y1 - y0) / Ny - sorted_vals = np.zeros((Nx, Ny)) + sorted_vals = [[[] for _ in range(Ny)] for _ in range(Nx)] - for pos in positions: + for p in range(len(positions)): + pos = positions[p] x = x0 y = y0 brk = False @@ -18,11 +19,10 @@ def _sort_2d(positions, x0, x1, y0, y1, Nx, Ny): y = y0 for j in range(Ny): - print(x, y) a = pos[0] >= x and pos[0] <= (x + dx) b = pos[1] >= y and pos[1] <= (y + dy) if a and b: - sorted_vals[i][j] += 1 + sorted_vals[i][j].append(p) brk = True break y += dy diff --git a/tests/diagnostics/dia_test.py b/tests/diagnostics/dia_test.py index f4ce5ae..4836fbc 100644 --- a/tests/diagnostics/dia_test.py +++ b/tests/diagnostics/dia_test.py @@ -1,5 +1,6 @@ import numpy as np import dsmc.diagnostics as dia +import time def create_particles(N, radius): positions = np.empty((N + 1, 3)) @@ -20,7 +21,7 @@ def create_particles(N, radius): return positions if __name__ == "__main__": - N = 1000 + N = 100000 radius = 1.0 positions = create_particles(N, radius) @@ -32,8 +33,8 @@ def create_particles(N, radius): y0 = -1.5 y1 = 1.5 + start_time = time.time() sorted_vals = dia._sort_2d(positions, x0, x1, y0, y1, Nx, Ny) - - print(sorted_vals) + print("--- %s seconds ---" % (time.time() - start_time)) print("done") \ No newline at end of file From cb55622bf57f357ed651ed1db99801bba96f2aeb Mon Sep 17 00:00:00 2001 From: LeoBasov <40352679+LeoBasov@users.noreply.github.com> Date: Mon, 3 Oct 2022 16:36:50 +0200 Subject: [PATCH 03/11] moved swap operation to common --- dsmc/common.py | 10 ++++++++-- dsmc/octree.py | 9 ++------- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/dsmc/common.py b/dsmc/common.py index 0990d48..17bc81d 100644 --- a/dsmc/common.py +++ b/dsmc/common.py @@ -1,5 +1,11 @@ -from numba import cfunc +from numba import cfunc, njit @cfunc("double(double[::3, ::2])") def get_V(a): - return (a[0][1] - a[0][0]) * (a[1][1] - a[1][0]) * (a[2][1] - a[2][0]) \ No newline at end of file + return (a[0][1] - a[0][0]) * (a[1][1] - a[1][0]) * (a[2][1] - a[2][0]) + +@njit +def swap(arr, pos1, pos2): + temp = arr[pos2] + arr[pos2] = arr[pos1] + arr[pos1] = temp \ No newline at end of file diff --git a/dsmc/octree.py b/dsmc/octree.py index cd18815..dd40ed0 100644 --- a/dsmc/octree.py +++ b/dsmc/octree.py @@ -1,6 +1,7 @@ import numpy as np import numpy.typing as npt from numba import njit +from . import common as com fmin = np.finfo(float).min fmax = np.finfo(float).max @@ -69,12 +70,6 @@ def _is_inside(position : npt.NDArray, box : npt.NDArray) -> bool: 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]: @@ -103,7 +98,7 @@ def _sort(permutations : npt.NDArray, box : npt.NDArray, positions : npt.NDArray for i in range(offset, offset + N): p = new_permutations[i] if _is_inside(positions[p], box): - _swap(new_permutations, i, runner) + com.swap(new_permutations, i, runner) runner += 1 Nnew += 1 From 26f74543ad1c7fff7e74edfd1d289aeff31a7125 Mon Sep 17 00:00:00 2001 From: LeoBasov <40352679+LeoBasov@users.noreply.github.com> Date: Mon, 3 Oct 2022 17:11:12 +0200 Subject: [PATCH 04/11] implemented new sorting --- dsmc/diagnostics.py | 61 +++++++++++++++++++++++------------ tests/diagnostics/dia_test.py | 5 ++- 2 files changed, 42 insertions(+), 24 deletions(-) diff --git a/dsmc/diagnostics.py b/dsmc/diagnostics.py index 83fce44..003f1b1 100644 --- a/dsmc/diagnostics.py +++ b/dsmc/diagnostics.py @@ -1,34 +1,53 @@ import numpy as np +from numba import njit from . import particles as prt from . import octree as oc +from . import common as com -def _sort_2d(positions, x0, x1, y0, y1, Nx, Ny): +@njit +def _in_2d_box(pos, box): + a = pos[0] >= box[0][0] and pos[0] <= box[0][1] + b = pos[1] >= box[1][0] and pos[1] <= box[1][1] + return a and b + +@njit +def _set_up_2d_boxes(x0, x1, y0, y1, Nx, Ny): dx = (x1 - x0) / Nx dy = (y1 - y0) / Ny - sorted_vals = [[[] for _ in range(Ny)] for _ in range(Nx)] + boxes = np.empty((Nx*Ny, 2, 2)) + y = y0 - for p in range(len(positions)): - pos = positions[p] + for i in range(Ny): x = x0 - y = y0 - brk = False - - for i in range(Nx): - if brk: - break - y = y0 - - for j in range(Ny): - a = pos[0] >= x and pos[0] <= (x + dx) - b = pos[1] >= y and pos[1] <= (y + dy) - if a and b: - sorted_vals[i][j].append(p) - brk = True - break - y += dy + for j in range(Nx): + k = j + i*Nx + boxes[k][0][0] = x + boxes[k][0][1] = x + dx + boxes[k][1][0] = y + boxes[k][1][1] = y + dy x += dx + y += dy + + return boxes + +@njit +def _sort_2d(positions, x0, x1, y0, y1, Nx, Ny): + boxes = _set_up_2d_boxes(x0, x1, y0, y1, Nx, Ny) # [] : bix id, [][] : axis, [][][] : 0 min , 1 max + permutations = [i for i in range(len(positions))] + leafs = np.zeros((Nx * Ny, 2), dtype=np.int_) # leaf[0] : ofseat, leaf[1] : n-parts + + for i in range(len(leafs)): + leafs[i][0] = leafs[i - 1][0] + leafs[i - 1][1] if i > 0 else 0 + leafs[i][1] = 0 + runner = leafs[i][0] + for j in range(leafs[i][0], len(positions)): + p = permutations[j] + if _in_2d_box(positions[p], boxes[i]): + com.swap(permutations, j, runner) + runner += 1 + leafs[i][1] += 1 - return sorted_vals + return (permutations, leafs, boxes) def sort_bin(positions, axis, Nbin): diff --git a/tests/diagnostics/dia_test.py b/tests/diagnostics/dia_test.py index 4836fbc..2dac231 100644 --- a/tests/diagnostics/dia_test.py +++ b/tests/diagnostics/dia_test.py @@ -34,7 +34,6 @@ def create_particles(N, radius): y1 = 1.5 start_time = time.time() - sorted_vals = dia._sort_2d(positions, x0, x1, y0, y1, Nx, Ny) - print("--- %s seconds ---" % (time.time() - start_time)) - + permutations, leafs, boxes = dia._sort_2d(positions, x0, x1, y0, y1, Nx, Ny) + print("--- %s seconds ---" % (time.time() - start_time)) print("done") \ No newline at end of file From e2328148badf5834874e6a2e5bb61a526494229a Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Mon, 3 Oct 2022 19:03:25 +0200 Subject: [PATCH 05/11] moved writer to sub module --- dsmc/writer/__init__.py | 1 + dsmc/{writer.py => writer/octree.py} | 0 2 files changed, 1 insertion(+) create mode 100644 dsmc/writer/__init__.py rename dsmc/{writer.py => writer/octree.py} (100%) diff --git a/dsmc/writer/__init__.py b/dsmc/writer/__init__.py new file mode 100644 index 0000000..e170db2 --- /dev/null +++ b/dsmc/writer/__init__.py @@ -0,0 +1 @@ +from .octree import write_buttom_leafs \ No newline at end of file diff --git a/dsmc/writer.py b/dsmc/writer/octree.py similarity index 100% rename from dsmc/writer.py rename to dsmc/writer/octree.py From 586faf26021baa4dc625e1f3d70f6c62ef20465b Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Mon, 3 Oct 2022 19:10:42 +0200 Subject: [PATCH 06/11] started implementation of planar writer --- dsmc/writer/planar.py | 86 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 86 insertions(+) create mode 100644 dsmc/writer/planar.py diff --git a/dsmc/writer/planar.py b/dsmc/writer/planar.py new file mode 100644 index 0000000..f218a32 --- /dev/null +++ b/dsmc/writer/planar.py @@ -0,0 +1,86 @@ +def write(boxes, values, file_name="planar.vtu"): + f = open(file_name, "w") + + _write_header(f) + _wrtie_body(f, boxes, values) + _write_footer(f) + + f.close() + +def _write_header(f): + f.write("\n") + +def _wrtie_body(f, boxes, values): + f.write("\n") + f.write("\n".format(len(boxes) * 4, len(boxes))) + + _write_points(f, boxes) + _write_cells(f, boxes) + #_write_cell_data(f, boxes) + + f.write("\n") + f.write("\n") + +def _write_points(f, boxes): + f.write("\n") + f.write("\n") + + for box in boxes: + f.write("{} ".format(box[0][0])) + f.write("{} ".format(box[1][0])) + f.write("{} ".format(0.0)) + + f.write("{} ".format(box[0][1])) + f.write("{} ".format(box[1][0])) + f.write("{} ".format(0.0)) + + f.write("{} ".format(box[0][1])) + f.write("{} ".format(box[1][1])) + f.write("{} ".format(0.0)) + + f.write("{} ".format(box[0][0])) + f.write("{} ".format(box[1][1])) + f.write("{} ".format(0.0)) + + f.write("\n") + f.write("\n") + + +def _write_cells(f, boxes): + k = 0 + + f.write("\n") + f.write("\n") + + for i in range(len(leaf_ids)): + for _ in range(8): + f.write("{} ".format(k)) + k += 1 + + f.write("\n") + f.write("\n") + + for i in range(len(leaf_ids)): + f.write("{} ".format((i + 1) * 8)) + + f.write("\n") + f.write("\n") + + for _ in range(len(leaf_ids)): + f.write("12 ") + + f.write("\n") + f.write("\n") + +def _write_cell_data(f, octree, leaf_ids): + f.write("\n") + f.write("\n") + + for i in leaf_ids: + f.write("{} ".format(octree.leafs[i].number_elements)) + + f.write("\n") + f.write("\n") + +def _write_footer(f): + f.write("\n") From d7b890f5e9cf3fbe099d31f7a84c326c1ff0071d Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Mon, 3 Oct 2022 19:37:02 +0200 Subject: [PATCH 07/11] updated planar writer --- dsmc/writer/planar.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/dsmc/writer/planar.py b/dsmc/writer/planar.py index f218a32..bc4c8d3 100644 --- a/dsmc/writer/planar.py +++ b/dsmc/writer/planar.py @@ -52,22 +52,22 @@ def _write_cells(f, boxes): f.write("\n") f.write("\n") - for i in range(len(leaf_ids)): - for _ in range(8): + for i in range(len(boxes)): + for _ in range(4): f.write("{} ".format(k)) k += 1 f.write("\n") f.write("\n") - for i in range(len(leaf_ids)): - f.write("{} ".format((i + 1) * 8)) + for i in range(len(boxes)): + f.write("{} ".format((i + 1) * 4)) f.write("\n") f.write("\n") - for _ in range(len(leaf_ids)): - f.write("12 ") + for _ in range(len(boxes)): + f.write("9 ") f.write("\n") f.write("\n") From d6c984d29ee7d2e499a96c27910f205a7819eaab Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Mon, 3 Oct 2022 19:37:28 +0200 Subject: [PATCH 08/11] updated init --- dsmc/writer/__init__.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/dsmc/writer/__init__.py b/dsmc/writer/__init__.py index e170db2..6847e5f 100644 --- a/dsmc/writer/__init__.py +++ b/dsmc/writer/__init__.py @@ -1 +1,2 @@ -from .octree import write_buttom_leafs \ No newline at end of file +from .octree import write_buttom_leafs +from .planar import write as write_planar From 7ebeb3b37dc0a0fb9843ab38382109458e8958f1 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Mon, 3 Oct 2022 19:59:58 +0200 Subject: [PATCH 09/11] updated setup file --- setup.py | 2 +- tests/diagnostics/dia_test.py | 7 ++++++- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/setup.py b/setup.py index 9a176ef..0b8444c 100644 --- a/setup.py +++ b/setup.py @@ -4,7 +4,7 @@ version='0.7.0', author='Leo Basov', python_requires='>=3.10, <4', - packages=["dsmc"], + packages=["dsmc", "dsmc.writer"], install_requires=[ 'numpy', 'llvmlite', diff --git a/tests/diagnostics/dia_test.py b/tests/diagnostics/dia_test.py index 2dac231..52dbc5b 100644 --- a/tests/diagnostics/dia_test.py +++ b/tests/diagnostics/dia_test.py @@ -1,5 +1,6 @@ import numpy as np import dsmc.diagnostics as dia +import dsmc.writer as wrt import time def create_particles(N, radius): @@ -21,7 +22,7 @@ def create_particles(N, radius): return positions if __name__ == "__main__": - N = 100000 + N = 1000 radius = 1.0 positions = create_particles(N, radius) @@ -36,4 +37,8 @@ def create_particles(N, radius): start_time = time.time() permutations, leafs, boxes = dia._sort_2d(positions, x0, x1, y0, y1, Nx, Ny) print("--- %s seconds ---" % (time.time() - start_time)) + + print("writing to file") + wrt.write_planar(boxes, N) + print("done") \ No newline at end of file From 333d9c7d07c3cd818c27936787438611dd292abd Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Mon, 3 Oct 2022 20:09:23 +0200 Subject: [PATCH 10/11] updated test case --- dsmc/writer/planar.py | 10 +++++----- tests/diagnostics/dia_test.py | 17 +++++++++++++---- 2 files changed, 18 insertions(+), 9 deletions(-) diff --git a/dsmc/writer/planar.py b/dsmc/writer/planar.py index bc4c8d3..e6101d2 100644 --- a/dsmc/writer/planar.py +++ b/dsmc/writer/planar.py @@ -16,7 +16,7 @@ def _wrtie_body(f, boxes, values): _write_points(f, boxes) _write_cells(f, boxes) - #_write_cell_data(f, boxes) + _write_cell_data(f, boxes, values) f.write("\n") f.write("\n") @@ -72,12 +72,12 @@ def _write_cells(f, boxes): f.write("\n") f.write("\n") -def _write_cell_data(f, octree, leaf_ids): - f.write("\n") +def _write_cell_data(f, boxes, values): + f.write("\n") f.write("\n") - for i in leaf_ids: - f.write("{} ".format(octree.leafs[i].number_elements)) + for vals in values: + f.write("{} ".format(vals.number_elements)) f.write("\n") f.write("\n") diff --git a/tests/diagnostics/dia_test.py b/tests/diagnostics/dia_test.py index 52dbc5b..17338d6 100644 --- a/tests/diagnostics/dia_test.py +++ b/tests/diagnostics/dia_test.py @@ -21,14 +21,18 @@ def create_particles(N, radius): return positions +class Values: + def __init__(self): + self.number_elements = 0 + if __name__ == "__main__": - N = 1000 + N = 10000 radius = 1.0 positions = create_particles(N, radius) # set up sort - Nx = 10 - Ny = 10 + Nx = 100 + Ny = 100 x0 = -1.5 x1 = 1.5 y0 = -1.5 @@ -39,6 +43,11 @@ def create_particles(N, radius): print("--- %s seconds ---" % (time.time() - start_time)) print("writing to file") - wrt.write_planar(boxes, N) + values = [Values() for _ in range(len(leafs))] + + for i in range(len(leafs)): + values[i].number_elements = leafs[i][1] + + wrt.write_planar(boxes, values) print("done") \ No newline at end of file From 7d9c3ba8e4deba571a89986bec2491cb453fd072 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Mon, 3 Oct 2022 20:13:52 +0200 Subject: [PATCH 11/11] updated diagnnostics --- dsmc/diagnostics.py | 13 +++++++++++++ tests/diagnostics/dia_test.py | 10 +--------- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/dsmc/diagnostics.py b/dsmc/diagnostics.py index 003f1b1..97a7e58 100644 --- a/dsmc/diagnostics.py +++ b/dsmc/diagnostics.py @@ -48,6 +48,19 @@ def _sort_2d(positions, x0, x1, y0, y1, Nx, Ny): leafs[i][1] += 1 return (permutations, leafs, boxes) + +class Values: + def __init__(self): + self.number_elements = 0 + +def analyse_2d(positions, x0, x1, y0, y1, Nx, Ny): + permutations, leafs, boxes = _sort_2d(positions, x0, x1, y0, y1, Nx, Ny) + values = [Values() for _ in range(len(leafs))] + + for i in range(len(leafs)): + values[i].number_elements = leafs[i][1] + + return (boxes, values) def sort_bin(positions, axis, Nbin): diff --git a/tests/diagnostics/dia_test.py b/tests/diagnostics/dia_test.py index 17338d6..768b85d 100644 --- a/tests/diagnostics/dia_test.py +++ b/tests/diagnostics/dia_test.py @@ -21,10 +21,6 @@ def create_particles(N, radius): return positions -class Values: - def __init__(self): - self.number_elements = 0 - if __name__ == "__main__": N = 10000 radius = 1.0 @@ -39,14 +35,10 @@ def __init__(self): y1 = 1.5 start_time = time.time() - permutations, leafs, boxes = dia._sort_2d(positions, x0, x1, y0, y1, Nx, Ny) + boxes, values = dia.analyse_2d(positions, x0, x1, y0, y1, Nx, Ny) print("--- %s seconds ---" % (time.time() - start_time)) print("writing to file") - values = [Values() for _ in range(len(leafs))] - - for i in range(len(leafs)): - values[i].number_elements = leafs[i][1] wrt.write_planar(boxes, values)