From e292379788fb0afbd487977bbb96825fa9b97a71 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 11:56:35 +0200 Subject: [PATCH 01/18] added diagnostics --- dsmc/diagnostics.py | 39 +++++++++++++++++++++++++++++ tests/unit/main.py | 1 + tests/unit/test_dsmc/diagnostics.py | 20 +++++++++++++++ 3 files changed, 60 insertions(+) create mode 100644 dsmc/diagnostics.py create mode 100644 tests/unit/test_dsmc/diagnostics.py diff --git a/dsmc/diagnostics.py b/dsmc/diagnostics.py new file mode 100644 index 0000000..5ff8ca5 --- /dev/null +++ b/dsmc/diagnostics.py @@ -0,0 +1,39 @@ +import numpy as np +from numba import njit +from . import particles as prt +from . import octree as oc + +def sort_bin(positions, axis, Nbin): + bins = [[] for _ in range(Nbin)] + b = 0 + box = oc._find_bounding_box(positions) + dx = (box[axis][1] - box[axis][0]) / (Nbin - 1) + x = dx + sub_pos = np.array([pos[axis] for pos in positions]) + sorted_pos = np.argsort(sub_pos) + + for i in range(len(sorted_pos)): + p = sorted_pos[i] + while positions[p][axis] > x: + x += dx + b += 1 + + bins[b].append(p) + + return bins, box + +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): + if i == axis: + V *= (box[i][1] - box[i][0]) / Nbins + else: + V *= (box[i][1] - box[i][0]) + + for i in range(Nbins): + n[i] = len(bins[i]) * w / V + + return n diff --git a/tests/unit/main.py b/tests/unit/main.py index da2edd5..f9c8d93 100644 --- a/tests/unit/main.py +++ b/tests/unit/main.py @@ -5,6 +5,7 @@ from test_dsmc.dsmc import * from test_dsmc.particles import * from test_dsmc.octree import * +from test_dsmc.diagnostics import * if __name__ == '__main__': unittest.main() diff --git a/tests/unit/test_dsmc/diagnostics.py b/tests/unit/test_dsmc/diagnostics.py new file mode 100644 index 0000000..e4587c1 --- /dev/null +++ b/tests/unit/test_dsmc/diagnostics.py @@ -0,0 +1,20 @@ +import numpy as np +import dsmc.diagnostics as dia +import unittest + +class TestDiagnostics(unittest.TestCase): + def test_sort_bin(self): + positions = np.array(([1, 2, 3], [9, 8, 7], [10, 11, 12], [4, 5, 6])) + Nbins = 4 + + bins1, box = dia.sort_bin(positions, 0, Nbins) + + self.assertEqual(Nbins, len(bins1)) + + for b in bins1: + self.assertEqual(1, len(b)) + + self.assertEqual(0, bins1[0][0]) + self.assertEqual(3, bins1[1][0]) + self.assertEqual(1, bins1[2][0]) + self.assertEqual(2, bins1[3][0]) From 11481da8f38c6e38cb7620f8de2b88398e36351c Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 12:02:33 +0200 Subject: [PATCH 02/18] fixed bug in particle creation --- dsmc/dsmc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index 5af06e9..033b4e9 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -127,7 +127,7 @@ def _update_velocities(self, dt): 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): - N = int(round(n / self.w)) + N = int(round(oc.get_V(box) * n / self.w)) print("creating {} particles".format(N)) self.particles.create_particles(box, self.mass, T, N) From 0affd51c82e3c92b421a729a85b0a0f7796a885a Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 12:11:03 +0200 Subject: [PATCH 03/18] updated shocktube example --- examples/plot.py | 28 +++------------------------- examples/shock_tube.py | 14 +++++++++++--- 2 files changed, 14 insertions(+), 28 deletions(-) diff --git a/examples/plot.py b/examples/plot.py index b4e91ae..34c097a 100644 --- a/examples/plot.py +++ b/examples/plot.py @@ -7,32 +7,10 @@ with open('test.csv') as file: reader = csv.reader(file, delimiter=',') - res = [] for line in reader: l = [m for m in line if m] - data = (np.array(l)).astype(float) - data = np.sort(data) + n = [float(l[i]) for i in range(1, len(l))] - N = 100 - sor = np.zeros((N, )) - x = np.zeros((N, )) - dx = 0.05/N - x[0] = dx - q = 0 - - for i in range(len(data)): - while data[i] > x[q]: - q += 1 - x[q] = x[q - 1] + dx - - sor[q] += 1 - - - x.resize(q) - sor.resize(q) - - res.append((x, sor)) - - plt.plot(res[-1][0], res[-1][1]) - plt.show() + plt.plot(n) + plt.show() diff --git a/examples/shock_tube.py b/examples/shock_tube.py index 9b61c08..e8ff68c 100644 --- a/examples/shock_tube.py +++ b/examples/shock_tube.py @@ -1,11 +1,13 @@ import dsmc +import dsmc.diagnostics as dia +import matplotlib.pyplot as plt if __name__ == '__main__': # general parameters solver = dsmc.DSMC() domain = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (0, 50e-3)) dt = 1e-7 - w = 2.4134e+15 + w = 2.4134e+7 mass = 6.6422e-26 niter = 300 @@ -31,8 +33,14 @@ print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) solver.advance(dt) - for pos in solver.particles.Pos: - file.write("{:.4e},".format(pos[2])) + bins, box = dia.sort_bin(solver.particles.Pos, 2, 100) + N = [len(b) for b in bins] + n = dia.calc_n(bins, box, 2, solver.w) + + file.write("{}".format(it*dt)) + + for ni in n: + file.write(",{}".format(ni)) file.write("\n") From bc36b8ac92cff731e6d6d3a9d7291ea48387f095 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 12:30:48 +0200 Subject: [PATCH 04/18] implemented calc T --- dsmc/diagnostics.py | 10 ++++++++++ examples/plot.py | 10 +++++++--- examples/shock_tube.py | 35 ++++++++++++++++++++++------------- 3 files changed, 39 insertions(+), 16 deletions(-) diff --git a/dsmc/diagnostics.py b/dsmc/diagnostics.py index 5ff8ca5..89cc454 100644 --- a/dsmc/diagnostics.py +++ b/dsmc/diagnostics.py @@ -37,3 +37,13 @@ def calc_n(bins, box, axis, w): n[i] = len(bins[i]) * w / V return n + +def calc_T(bins, velocities, mass): + Nbins = len(bins) + T = np.empty((Nbins, )) + + for i in range(Nbins): + vels = np.array([velocities[p] for p in bins[i]]) + T[i] = prt.calc_temperature(vels, mass) + + return T diff --git a/examples/plot.py b/examples/plot.py index 34c097a..18ac59c 100644 --- a/examples/plot.py +++ b/examples/plot.py @@ -1,16 +1,20 @@ import matplotlib.pyplot as plt import csv import numpy as np +import argparse if __name__ == '__main__': - res = [] + parser = argparse.ArgumentParser(description='Process some integers.') + parser.add_argument('file_name', type=str) - with open('test.csv') as file: + args = parser.parse_args() + + with open(args.file_name) as file: reader = csv.reader(file, delimiter=',') for line in reader: l = [m for m in line if m] - n = [float(l[i]) for i in range(1, len(l))] + n = [float(l[i]) for i in range(len(l))] plt.plot(n) plt.show() diff --git a/examples/shock_tube.py b/examples/shock_tube.py index e8ff68c..7ccdc28 100644 --- a/examples/shock_tube.py +++ b/examples/shock_tube.py @@ -10,6 +10,7 @@ w = 2.4134e+7 mass = 6.6422e-26 niter = 300 + Nbins = 100 # low denisty particles nhigh = 2.5e+20 @@ -28,21 +29,29 @@ solver.create_particles(Boxlow, Tlow, nlow) solver.create_particles(Boxhigh, Thigh, nhigh) - with open("test.csv", "w") as file: - for it in range(niter): - print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) - solver.advance(dt) - - bins, box = dia.sort_bin(solver.particles.Pos, 2, 100) - N = [len(b) for b in bins] - n = dia.calc_n(bins, box, 2, solver.w) - - file.write("{}".format(it*dt)) + n_file = open("n.csv", "w") + T_file = open("T.csv", "w") + p_file = open("p.csv", "w") + + for it in range(niter): + print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) + solver.advance(dt) - for ni in n: - file.write(",{}".format(ni)) + bins, box = dia.sort_bin(solver.particles.Pos, 2, Nbins) + N = [len(b) for b in bins] + n = dia.calc_n(bins, box, 2, solver.w) + T = dia.calc_T(bins, solver.particles.Vel, mass) + + for i in range(Nbins): + n_file.write("{},".format(n[i])) + T_file.write("{},".format(T[i])) - file.write("\n") + n_file.write("\n") + T_file.write("\n") + n_file.close() + T_file.close() + p_file.close() + print("") print('done') From fb446213a348eadeef04bb49e4b14eaa82ef676a Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 12:34:10 +0200 Subject: [PATCH 05/18] implemented write p --- dsmc/diagnostics.py | 9 +++++++++ examples/shock_tube.py | 3 +++ 2 files changed, 12 insertions(+) diff --git a/dsmc/diagnostics.py b/dsmc/diagnostics.py index 89cc454..3b113a3 100644 --- a/dsmc/diagnostics.py +++ b/dsmc/diagnostics.py @@ -47,3 +47,12 @@ def calc_T(bins, velocities, mass): T[i] = prt.calc_temperature(vels, mass) return T + +def calc_p(n, T): + Nbins = len(n) + p = np.empty((Nbins, )) + + for i in range(Nbins): + p[i] = n[i]*T[i]*prt.kb + + return p diff --git a/examples/shock_tube.py b/examples/shock_tube.py index 7ccdc28..45eb2cf 100644 --- a/examples/shock_tube.py +++ b/examples/shock_tube.py @@ -41,13 +41,16 @@ N = [len(b) for b in bins] n = dia.calc_n(bins, box, 2, solver.w) T = dia.calc_T(bins, solver.particles.Vel, mass) + p = dia.calc_p(n, T) for i in range(Nbins): n_file.write("{},".format(n[i])) T_file.write("{},".format(T[i])) + p_file.write("{},".format(p[i])) n_file.write("\n") T_file.write("\n") + p_file.write("\n") n_file.close() T_file.close() From 63c30aeafe3c19d739bb72e5f60a5bb6a058767b Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 13:34:11 +0200 Subject: [PATCH 06/18] implemented new velocity creation --- dsmc/dsmc.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index 033b4e9..ae69dd4 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -126,10 +126,19 @@ def _update_velocities(self, dt): 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): + def create_particles(self, box, T, n, u = None): N = int(round(oc.get_V(box) * n / self.w)) print("creating {} particles".format(N)) - self.particles.create_particles(box, self.mass, T, N) + self.particles.create_particles(np.array(box), self.mass, T, N) + + if u is not None: + velocities = self.particles.Vel + u = np.array(box) + + 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)) From 5eeeee6024e5f2a6067df34a865845eb5cdf0958 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 14:05:13 +0200 Subject: [PATCH 07/18] fixed bug --- dsmc/dsmc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index ae69dd4..6512fcf 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -133,7 +133,7 @@ def create_particles(self, box, T, n, u = None): if u is not None: velocities = self.particles.Vel - u = np.array(box) + u = np.array(u) for i in range(len(velocities)): velocities[i] += u From cbddc8f8724d727b6994bdb7fa8f0019ca2dbe63 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 14:16:57 +0200 Subject: [PATCH 08/18] added temp relax test case --- examples/temp_relax.py | 50 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) create mode 100644 examples/temp_relax.py diff --git a/examples/temp_relax.py b/examples/temp_relax.py new file mode 100644 index 0000000..8961579 --- /dev/null +++ b/examples/temp_relax.py @@ -0,0 +1,50 @@ +import dsmc +import dsmc.diagnostics as dia +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__': + # general parameters + solver = dsmc.DSMC() + domain = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (0, 50e-3)) + dt = 1e-5 + w = 2.4134e+7 + mass = 6.6422e-26 + niter = 500 + Nbins = 100 + + # particles + n = 2.5e+20 + T = 300 + Box = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (25e-3, 50e-3)) + u = np.array((1000.0, 0.0, 0.0)) + + + solver.set_domain(domain) + solver.set_weight(w) + solver.set_mass(mass) + + solver.create_particles(Box, T, n, u) + + for it in range(niter): + print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) + solver.advance(dt) + 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() + + print("") + print('done') From 8a641b33e14b9e3fe9fffe53545a1218002c7083 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 14:39:33 +0200 Subject: [PATCH 09/18] removed unsuported data type --- dsmc/particles.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dsmc/particles.py b/dsmc/particles.py index 9ca72eb..54723fc 100644 --- a/dsmc/particles.py +++ b/dsmc/particles.py @@ -57,7 +57,7 @@ def get_vel(T, mass): @njit def get_velocities(T, mass, N): - velocities = np.empty((N, 3), dtype=float) + velocities = np.empty((N, 3)) for i in range(N): velocities[i] = get_vel(T, mass) @@ -88,7 +88,7 @@ def calc_positions(x, y, z, N): z : tuple(2), dtype = float zmin, zmax """ - positions = np.empty((N, 3), dtype=float) + positions = np.empty((N, 3)) for i in range(N): positions[i][0] = x[0] + np.random.random() * (x[1] - x[0]) From 07ef757036f5c571d6cfe09dbaf242d682e70715 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 14:52:00 +0200 Subject: [PATCH 10/18] fixed bug in dsmc --- dsmc/dsmc.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index 6512fcf..ad25fa4 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -127,9 +127,10 @@ def _update_velocities(self, dt): 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(np.array(box), self.mass, T, N) + self.particles.create_particles(box, self.mass, T, N) if u is not None: velocities = self.particles.Vel From f490a3fd2a567ed11aa65c0cbb0739688d667443 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 15:03:11 +0200 Subject: [PATCH 11/18] updated shock tube test case --- examples/shock_tube.py | 49 +++++++++++++++++++++++------------------- 1 file changed, 27 insertions(+), 22 deletions(-) diff --git a/examples/shock_tube.py b/examples/shock_tube.py index 45eb2cf..5d49de4 100644 --- a/examples/shock_tube.py +++ b/examples/shock_tube.py @@ -2,26 +2,43 @@ import dsmc.diagnostics as dia import matplotlib.pyplot as plt +def write2file(solver, n_file, T_file, p_file): + bins, box = dia.sort_bin(solver.particles.Pos, 2, Nbins) + N = [len(b) for b in bins] + n = dia.calc_n(bins, box, 2, solver.w) + T = dia.calc_T(bins, solver.particles.Vel, mass) + p = dia.calc_p(n, T) + + for i in range(Nbins): + n_file.write("{},".format(n[i])) + T_file.write("{},".format(T[i])) + p_file.write("{},".format(p[i])) + + n_file.write("\n") + T_file.write("\n") + p_file.write("\n") + if __name__ == '__main__': # general parameters solver = dsmc.DSMC() - domain = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (0, 50e-3)) + domain = [(-0.0001, 0.0001), (-0.0001, 0.0001), (0.0, 0.1)] dt = 1e-7 - w = 2.4134e+7 + w = 1e+9 mass = 6.6422e-26 niter = 300 Nbins = 100 # low denisty particles - nhigh = 2.5e+20 + nhigh = 2.41432e+22 Thigh = 300 - Boxhigh = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (25e-3, 50e-3)) + Boxhigh = [(-0.0001, 0.0001), (-0.0001, 0.0001), (0.05, 0.1)] # high denisty particles - nlow = 2.5e+19 + nlow = 2.41432e+21 Tlow = 300 - Boxlow = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (0, 25e-3)) + Boxlow = [(-0.0001, 0.0001), (-0.0001, 0.0001), (0.0, 0.05)] + # setup solver solver.set_domain(domain) solver.set_weight(w) solver.set_mass(mass) @@ -29,29 +46,17 @@ solver.create_particles(Boxlow, Tlow, nlow) solver.create_particles(Boxhigh, Thigh, nhigh) + # open files n_file = open("n.csv", "w") T_file = open("T.csv", "w") p_file = open("p.csv", "w") for it in range(niter): print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) - solver.advance(dt) - - bins, box = dia.sort_bin(solver.particles.Pos, 2, Nbins) - N = [len(b) for b in bins] - n = dia.calc_n(bins, box, 2, solver.w) - T = dia.calc_T(bins, solver.particles.Vel, mass) - p = dia.calc_p(n, T) - - for i in range(Nbins): - n_file.write("{},".format(n[i])) - T_file.write("{},".format(T[i])) - p_file.write("{},".format(p[i])) - - n_file.write("\n") - T_file.write("\n") - p_file.write("\n") + solver.advance(dt) + write2file(solver, n_file, T_file, p_file) + # close files n_file.close() T_file.close() p_file.close() From e8b5efb9c50eca563fc286b8e01b408db9ec52ef Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 15:11:03 +0200 Subject: [PATCH 12/18] updated test case --- dsmc/diagnostics.py | 4 +++- examples/shock_tube.py | 10 +++++----- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/dsmc/diagnostics.py b/dsmc/diagnostics.py index 3b113a3..dd6aeb4 100644 --- a/dsmc/diagnostics.py +++ b/dsmc/diagnostics.py @@ -8,6 +8,7 @@ def sort_bin(positions, axis, Nbin): b = 0 box = oc._find_bounding_box(positions) dx = (box[axis][1] - box[axis][0]) / (Nbin - 1) + xx = [dx] x = dx sub_pos = np.array([pos[axis] for pos in positions]) sorted_pos = np.argsort(sub_pos) @@ -17,10 +18,11 @@ def sort_bin(positions, axis, Nbin): while positions[p][axis] > x: x += dx b += 1 + xx.append(x) bins[b].append(p) - return bins, box + return bins, box, xx def calc_n(bins, box, axis, w): Nbins = len(bins) diff --git a/examples/shock_tube.py b/examples/shock_tube.py index 5d49de4..b10df19 100644 --- a/examples/shock_tube.py +++ b/examples/shock_tube.py @@ -3,7 +3,7 @@ import matplotlib.pyplot as plt def write2file(solver, n_file, T_file, p_file): - bins, box = dia.sort_bin(solver.particles.Pos, 2, Nbins) + bins, box, x = dia.sort_bin(solver.particles.Pos, 2, Nbins) N = [len(b) for b in bins] n = dia.calc_n(bins, box, 2, solver.w) T = dia.calc_T(bins, solver.particles.Vel, mass) @@ -28,15 +28,15 @@ def write2file(solver, n_file, T_file, p_file): niter = 300 Nbins = 100 - # low denisty particles + # high denisty particles nhigh = 2.41432e+22 Thigh = 300 - Boxhigh = [(-0.0001, 0.0001), (-0.0001, 0.0001), (0.05, 0.1)] + Boxhigh = [(-0.0001, 0.0001), (-0.0001, 0.0001), (0.0, 0.05)] - # high denisty particles + # low denisty particles nlow = 2.41432e+21 Tlow = 300 - Boxlow = [(-0.0001, 0.0001), (-0.0001, 0.0001), (0.0, 0.05)] + Boxlow = [(-0.0001, 0.0001), (-0.0001, 0.0001), (0.05, 0.1)] # setup solver solver.set_domain(domain) From 4184879bd9472a9141053cb42ce5f7e3b409e644 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 16:01:52 +0200 Subject: [PATCH 13/18] implemented writer for octree --- dsmc/writer.py | 100 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 100 insertions(+) create mode 100644 dsmc/writer.py diff --git a/dsmc/writer.py b/dsmc/writer.py new file mode 100644 index 0000000..19a2a45 --- /dev/null +++ b/dsmc/writer.py @@ -0,0 +1,100 @@ +from . import octree as oc + +def write_buttom_leafs(octree): + f = open("octree.vtu", "w") + + _write_header(f) + _wrtie_body(f, octree) + _write_footer(f) + + f.close() + +def _write_header(f): + f.write("\n") + +def _wrtie_body(f, octree): + leaf_ids = [] + for i in range(len(octree.leafs)): + if octree.leafs[i].number_children == 0: + leaf_ids.append(i) + + f.write("\n") + f.write("\n".format(len(leaf_ids) * 8, len(leaf_ids))) + + _write_points(f, octree, leaf_ids); + _write_cells(f, octree, leaf_ids); + + f.write("\n") + f.write("\n") + +def _write_points(f, octree, leaf_ids): + f.write("\n") + f.write("\n") + + for i in leaf_ids: + box = octree.cell_boxes[i] + + f.write("{} ".format(box[0][0])) + f.write("{} ".format(box[1][0])) + f.write("{} ".format(box[2][0])) + + f.write("{} ".format(box[0][1])) + f.write("{} ".format(box[1][0])) + f.write("{} ".format(box[2][0])) + + f.write("{} ".format(box[0][1])) + f.write("{} ".format(box[1][1])) + f.write("{} ".format(box[2][0])) + + f.write("{} ".format(box[0][0])) + f.write("{} ".format(box[1][1])) + f.write("{} ".format(box[2][0])) + + f.write("{} ".format(box[0][0])) + f.write("{} ".format(box[1][0])) + f.write("{} ".format(box[2][1])) + + f.write("{} ".format(box[0][1])) + f.write("{} ".format(box[1][0])) + f.write("{} ".format(box[2][1])) + + f.write("{} ".format(box[0][1])) + f.write("{} ".format(box[1][1])) + f.write("{} ".format(box[2][1])) + + f.write("{} ".format(box[0][0])) + f.write("{} ".format(box[1][1])) + f.write("{} ".format(box[2][1])) + + f.write("\n") + f.write("\n") + + +def _write_cells(f, octree, leaf_ids): + 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_footer(f): + f.write("\n") From 34689ac2269200d81334dcf9155ce7b357ee93b0 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 16:05:51 +0200 Subject: [PATCH 14/18] added octree test --- examples/octree_test.py | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) create mode 100644 examples/octree_test.py diff --git a/examples/octree_test.py b/examples/octree_test.py new file mode 100644 index 0000000..349d391 --- /dev/null +++ b/examples/octree_test.py @@ -0,0 +1,36 @@ +import dsmc.writer as wrt +import dsmc.octree as oc +import numpy as np +import argparse + +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__": + parser = argparse.ArgumentParser(description='Process some integers.') + parser.add_argument('N', type=int) + + args = parser.parse_args() + + radius = 1.0 + positions = create_particles(args.N, radius) + octree = oc.Octree() + octree.build(positions) + wrt.write_buttom_leafs(octree) + + print("done") From f4ecb48d0be323e576a837110a9e4b0cb4f80b63 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 16:18:32 +0200 Subject: [PATCH 15/18] extenden octree test --- examples/octree_test.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/examples/octree_test.py b/examples/octree_test.py index 349d391..692f661 100644 --- a/examples/octree_test.py +++ b/examples/octree_test.py @@ -33,4 +33,9 @@ def create_particles(N, radius): octree.build(positions) wrt.write_buttom_leafs(octree) + 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])) + print("done") From e05cd40cf1037f97b209088adfe1fed0a8a35d75 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 16:19:33 +0200 Subject: [PATCH 16/18] moved octree test to tests --- {examples => tests/octree}/octree_test.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename {examples => tests/octree}/octree_test.py (100%) diff --git a/examples/octree_test.py b/tests/octree/octree_test.py similarity index 100% rename from examples/octree_test.py rename to tests/octree/octree_test.py From da25aea794e4d24437792cfd989d39daae084200 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 16:20:46 +0200 Subject: [PATCH 17/18] updated changelog and setup.py --- CHANGELOG.md | 6 +++++- setup.py | 2 +- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 74f67f4..99aa5e2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,7 +6,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] -## [0.3.0] - 2022-09-23 +## [0.5.0] - 2022-09-23 +### Added +- writer module + +## [0.4.0] - 2022-09-23 ### Added - dsmc module diff --git a/setup.py b/setup.py index 7963a05..ea7e25e 100644 --- a/setup.py +++ b/setup.py @@ -1,7 +1,7 @@ from setuptools import setup setup( name='dsmc', - version='0.4.0', + version='0.5.0', author='Leo Basov', python_requires='>=3.6, <4', packages=["dsmc"], From b79ab52ef06b5a57e9d8f1e28e5f9aa5c8b8dedb Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 16:22:52 +0200 Subject: [PATCH 18/18] updated unit tests --- tests/unit/test_dsmc/diagnostics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit/test_dsmc/diagnostics.py b/tests/unit/test_dsmc/diagnostics.py index e4587c1..ade3bfd 100644 --- a/tests/unit/test_dsmc/diagnostics.py +++ b/tests/unit/test_dsmc/diagnostics.py @@ -7,7 +7,7 @@ def test_sort_bin(self): positions = np.array(([1, 2, 3], [9, 8, 7], [10, 11, 12], [4, 5, 6])) Nbins = 4 - bins1, box = dia.sort_bin(positions, 0, Nbins) + bins1, box, x = dia.sort_bin(positions, 0, Nbins) self.assertEqual(Nbins, len(bins1))