From 5c13499c40ae1afab235a7d25969b30f4a97f527 Mon Sep 17 00:00:00 2001 From: LeoBasov <40352679+LeoBasov@users.noreply.github.com> Date: Sun, 2 Oct 2022 17:49:26 +0200 Subject: [PATCH 01/54] set up hyptersonic test case --- examples/hypersonic_flow.py | 46 +++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) create mode 100644 examples/hypersonic_flow.py diff --git a/examples/hypersonic_flow.py b/examples/hypersonic_flow.py new file mode 100644 index 0000000..df50453 --- /dev/null +++ b/examples/hypersonic_flow.py @@ -0,0 +1,46 @@ +import dsmc +import time +import numpy as np + +if __name__ == '__main__': + # general parameters + solver = dsmc.DSMC() + domain = [(-3.0, 3.0), (-1.5, 1.5), (-0.025, 0.025)] + obj = [(-0.25, 0.25), (-0.25, 0.25), (-0.5, 0.5)] + dt = 1e-6 + w = 2.3e+13 + mass = 6.6422e-26 + T = 273.0 + n = 2.6e+19 + u = np.array([0.0, -3043.0, 0.0]) + niter = 1 + + # setup solver + solver.set_domain(domain) + solver.set_weight(w) + solver.set_mass(mass) + + # set up boundary + solver.set_bc_type("xmin", "inflow") + solver.set_bc_type("xmax", "inflow") + solver.set_bc_type("ymin", "inflow") + + solver.set_bc_type("ymax", "open") + + solver.set_bc_values("xmin", T, n, u) + solver.set_bc_values("xmax", T, n, u) + solver.set_bc_values("ymin", T, n, u) + + # set object + solver.add_object(obj) + + # start timing + start_time = time.time() + + for it in range(niter): + print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) + solver.advance(dt) + + print("") + print("--- %s seconds ---" % (time.time() - start_time)) + print('done') \ No newline at end of file From 01745314260ded6c3e30697a378625a35c7abc0d Mon Sep 17 00:00:00 2001 From: LeoBasov <40352679+LeoBasov@users.noreply.github.com> Date: Sun, 2 Oct 2022 17:55:44 +0200 Subject: [PATCH 02/54] implemented particle and coll object check --- dsmc/dsmc.py | 61 +++++++++++++++++++++++++++++++++++-- examples/hypersonic_flow.py | 3 ++ 2 files changed, 62 insertions(+), 2 deletions(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index 4b7d212..9d92198 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -14,18 +14,21 @@ def _push(velocities, positions, dt): @njit def _boundary(velocities, positions, domain, boundary_conds): kept_parts = np.ones(positions.shape[0], dtype=np.uint) + old_positions = np.copy(positions) for p in prange(len(positions)): while not oc._is_inside(positions[p], domain) and kept_parts[p]: for i in range(3): if positions[p][i] < domain[i][0]: if boundary_conds[i][0] == 0: + old_positions[p][i] = positions[p][i] positions[p][i] = 2.0 * domain[i][0] - positions[p][i] velocities[p][i] *= -1.0 elif boundary_conds[i][0] == 1 or boundary_conds[i][0] == 2: kept_parts[p] = 0 if positions[p][i] > domain[i][1]: if boundary_conds[i][1] == 0: + old_positions[p][i] = positions[p][i] positions[p][i] = 2.0 * domain[i][1] - positions[p][i] velocities[p][i] *= -1.0 elif boundary_conds[i][1] == 1 or boundary_conds[i][0] == 2: @@ -35,16 +38,18 @@ def _boundary(velocities, positions, domain, boundary_conds): p = 0 new_velocities = np.empty((N, 3)) new_positions = np.empty((N, 3)) + new_old_positions = np.empty((N, 3)) for i in range(positions.shape[0]): if kept_parts[i] == 1: new_velocities[p] = velocities[i] new_positions[p] = positions[i] + new_old_positions[p] = old_positions[p] p += 1 else: continue - return (new_velocities, new_positions) + return (new_velocities, new_positions, new_old_positions) @njit def _check_positions(velocities, positions, old_positions, domain): @@ -69,6 +74,46 @@ def _check_positions(velocities, positions, old_positions, domain): return (new_velocities, new_positions) +@njit +def _check_created_particles(velocities, positions, obj): + kept_parts = np.ones(positions.shape[0], dtype=np.uint) + + for i in prange(positions.shape[0]): + if (not oc._is_inside(positions[i], obj)): + kept_parts[i] = 0 + + N = sum(kept_parts) + p = 0 + new_velocities = np.empty((N, 3)) + new_positions = np.empty((N, 3)) + + for i in prange(positions.shape[0]): + if kept_parts[i] == 1: + new_velocities[p] = velocities[i] + new_positions[p] = positions[i] + p += 1 + else: + continue + + return (new_velocities, new_positions) + +@njit +def _object(velocities, positions, old_positions, coll_obj): + for p in range(positions.shape[0]): + if oc._is_inside(positions[p], coll_obj): + for i in range(3): + if (old_positions[p][i] < coll_obj[i][0]): + old_positions[p][i] = positions[p][i] + positions[p][i] = 2.0 * coll_obj[i][0] - positions[p][i] + velocities[p][i] *= -1.0 + + if (old_positions[p][i] > coll_obj[i][1]): + old_positions[p][i] = positions[p][i] + positions[p][i] = 2.0 * coll_obj[i][1] - positions[p][i] + velocities[p][i] *= -1.0 + + return (velocities, positions, old_positions) + @njit def _calc_prob(rel_vel : float, sigma_T : float, Vc : float, dt : float, w : float, N : int) -> np.single: """ @@ -179,6 +224,7 @@ def clear(self): self.boundary = Boundary() self.sigma_T = 3.631681e-19 self.mass = None + self.objects = [] def advance(self, dt, collisions=True, octree=True): if self.domain is None: @@ -200,7 +246,12 @@ def advance(self, dt, collisions=True, octree=True): old_positions = np.copy(self.particles.Pos) positions = _push(self.particles.Vel, self.particles.Pos, dt) self.particles.VelPos = _check_positions(self.particles.Vel, positions, old_positions, self.domain) - self.particles.VelPos = _boundary(self.particles.Vel, self.particles.Pos, self.domain, self.boundary_conds) + velocities, positions, old_positions = _boundary(self.particles.Vel, self.particles.Pos, self.domain, self.boundary_conds) + + for obj in self.objects: + velocities, positions, old_positions = _object(velocities, positions, old_positions, obj) + + self.particles.VelPos = (velocities, positions) def _update_velocities(self, dt): Nleafs : int = len(self.octree.leafs) @@ -216,6 +267,9 @@ def create_particles(self, box, T, n, u = np.zeros(3)): N = int(round(com.get_V(box) * n / self.w)) print("creating {} particles".format(N)) self.particles.create_particles(box, self.mass, T, N, u) + + for obj in self.objects: + self.particles.VelPos = _check_created_particles(self.particles.Vel, self.particles.Pos, obj) print("now containing {} particles, {} total".format(N, self.particles.N)) @@ -245,3 +299,6 @@ def set_bc_values(self, boundary, T, n, u): self.boundary.u[i][j] = u print("boundary [" + boundary + "] set to values T : {}, n : {}, u : {}".format(T, n, u)) + + def add_object(self, coll_object): + self.objects.append(np.array(coll_object)) \ No newline at end of file diff --git a/examples/hypersonic_flow.py b/examples/hypersonic_flow.py index df50453..fef920d 100644 --- a/examples/hypersonic_flow.py +++ b/examples/hypersonic_flow.py @@ -34,6 +34,9 @@ # set object solver.add_object(obj) + # create particles + solver.create_particles(domain, T, n) + # start timing start_time = time.time() From 236ef0365062b4ebb932ccc72aa25a433bb89c8e Mon Sep 17 00:00:00 2001 From: LeoBasov <40352679+LeoBasov@users.noreply.github.com> Date: Sun, 2 Oct 2022 17:58:19 +0200 Subject: [PATCH 03/54] fixed bug --- dsmc/dsmc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index 9d92198..0636efd 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -79,7 +79,7 @@ def _check_created_particles(velocities, positions, obj): kept_parts = np.ones(positions.shape[0], dtype=np.uint) for i in prange(positions.shape[0]): - if (not oc._is_inside(positions[i], obj)): + if oc._is_inside(positions[i], obj): kept_parts[i] = 0 N = sum(kept_parts) From b12a56a8443c27d40cdd546ee0aeb5c1dec6ef45 Mon Sep 17 00:00:00 2001 From: LeoBasov <40352679+LeoBasov@users.noreply.github.com> Date: Sun, 2 Oct 2022 18:08:46 +0200 Subject: [PATCH 04/54] extended test case --- examples/hypersonic_flow.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/hypersonic_flow.py b/examples/hypersonic_flow.py index fef920d..4a751ff 100644 --- a/examples/hypersonic_flow.py +++ b/examples/hypersonic_flow.py @@ -1,4 +1,5 @@ import dsmc +import dsmc.writer as wrt import time import numpy as np @@ -8,7 +9,7 @@ domain = [(-3.0, 3.0), (-1.5, 1.5), (-0.025, 0.025)] obj = [(-0.25, 0.25), (-0.25, 0.25), (-0.5, 0.5)] dt = 1e-6 - w = 2.3e+13 + w = 2.3e+14 mass = 6.6422e-26 T = 273.0 n = 2.6e+19 @@ -43,6 +44,7 @@ for it in range(niter): print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) solver.advance(dt) + wrt.write_buttom_leafs(solver.octree) print("") print("--- %s seconds ---" % (time.time() - start_time)) From 3764b1b831f473987faac4f60dc502a5bd64b419 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Sun, 2 Oct 2022 19:55:01 +0200 Subject: [PATCH 05/54] implemented diagnostics sorting --- dsmc/diagnostics.py | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/dsmc/diagnostics.py b/dsmc/diagnostics.py index 56b0b7f..ef741f4 100644 --- a/dsmc/diagnostics.py +++ b/dsmc/diagnostics.py @@ -1,7 +1,33 @@ +from numba import njit import numpy as np from . import particles as prt from . import octree as oc +@njit +def sort_flat(positions, box, Nx, Ny): + vals = np.zeros((Nx, Ny)) + dx = (box[0][1] - box[0][0]) / Nx + dy = (box[1][1] - box[1][0]) / Ny + + for p in range(positions.shape[0]): + pos = positions[p] + X = box[0][0] + Y = box[1][0] + + for x in range(Nx): + for y in range(Ny): + a = (pos[0] < (X + dx)) and (pos[0] >= X) + b = (pos[1] < (Y + dy)) and (pos[1] >= Y) + + if a and b: + + vals[x][y] += 1 + + Y += dy + X += dx + + return vals.astype(np.int_) + def sort_bin(positions, axis, Nbin): bins = [[] for _ in range(Nbin)] b = 0 From d47162fa630abd2ce27137069c16211cf9696919 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Mon, 3 Oct 2022 20:21:29 +0200 Subject: [PATCH 06/54] updated version number --- CHANGELOG.md | 5 +++++ setup.py | 2 +- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5a0257d..49a2f7d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.8.0] - 2022-10-03 +### Added +- planar writer +- planar diagnostic + ## [0.7.0] - 2022-10-01 ### Added - inflow boundary condition diff --git a/setup.py b/setup.py index 0b8444c..1064716 100644 --- a/setup.py +++ b/setup.py @@ -1,7 +1,7 @@ from setuptools import setup setup( name='dsmc', - version='0.7.0', + version='0.8.0', author='Leo Basov', python_requires='>=3.10, <4', packages=["dsmc", "dsmc.writer"], From c4034744c56c7f69c37ac1b2e42d70f7bb289986 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Mon, 3 Oct 2022 20:40:42 +0200 Subject: [PATCH 07/54] updated hypersonic test case --- examples/hypersonic_flow.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/examples/hypersonic_flow.py b/examples/hypersonic_flow.py index 4a751ff..c19a89a 100644 --- a/examples/hypersonic_flow.py +++ b/examples/hypersonic_flow.py @@ -1,5 +1,6 @@ import dsmc import dsmc.writer as wrt +import dsmc.diagnostics as dia import time import numpy as np @@ -16,6 +17,14 @@ u = np.array([0.0, -3043.0, 0.0]) niter = 1 + # set up diagnostics values + x0 = domain[0][0] + x1 = domain[0][1] + y0 = domain[1][0] + y1 = domain[1][1] + Nx = 100 + Ny = 50 + # setup solver solver.set_domain(domain) solver.set_weight(w) @@ -44,7 +53,8 @@ for it in range(niter): print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) solver.advance(dt) - wrt.write_buttom_leafs(solver.octree) + boxes, values = dia.analyse_2d(solver.particles.Pos, x0, x1, y0, y1, Nx, Ny) + wrt.write_planar(boxes, values, "hypersonic_{}.vtu".format(it)) print("") print("--- %s seconds ---" % (time.time() - start_time)) From 4ce987c24dcc698c9b7acc7438db4d8e30991894 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Tue, 4 Oct 2022 09:38:48 +0200 Subject: [PATCH 08/54] set up domain obj test --- tests/domain/domain_obj.py | 71 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) create mode 100644 tests/domain/domain_obj.py diff --git a/tests/domain/domain_obj.py b/tests/domain/domain_obj.py new file mode 100644 index 0000000..ff52525 --- /dev/null +++ b/tests/domain/domain_obj.py @@ -0,0 +1,71 @@ +import dsmc +import time +import numpy as np + +def create_pos_and_vels(): + positions = np.zeros((6, 3)) + velocitiies = np.zeros((6, 3)) + + # x + positions[0][0] = -0.75 + velocitiies[0][0] = 1.0 + + positions[1][0] = 0.75 + velocitiies[1][0] = -1.0 + + # y + positions[2][1] = -0.75 + velocitiies[2][1] = 1.0 + + positions[3][1] = 0.75 + velocitiies[3][1] = -1.0 + + # z + positions[4][2] = -0.75 + velocitiies[4][2] = 1.0 + + positions[5][2] = 0.75 + velocitiies[5][2] = -1.0 + + return (velocitiies, positions) + +def write_partices(positions, it): + with open("particles_{}.csv".format(it), "w") as file: + file.write("x, y, z\n") + for pos in positions: + file.write("{}, {}, {}\n".format(pos[0], pos[1], pos[2])) + +if __name__ == '__main__': + # general parameters + solver = dsmc.DSMC() + domain = [(-1.0, 1.0), (-1.0, 1.0), (-1.0, 1.0)] + obj = [(-0.5, 0.5), (-0.5, 0.5), (-0.5, 0.5)] + dt = 0.01 + w = 2.3e+14 + mass = 6.6422e-26 + T = 273.0 + n = 2.6e+19 + niter = 100 + + # setup solver + solver.set_domain(domain) + solver.set_weight(w) + solver.set_mass(mass) + + # set object + solver.add_object(obj) + + # create particles + solver.particles.VelPos = create_pos_and_vels() + + # start timing + start_time = time.time() + + for it in range(niter): + print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) + solver.advance(dt) + write_partices(solver.particles.Pos, it) + + print("") + print("--- %s seconds ---" % (time.time() - start_time)) + print('done') \ No newline at end of file From f0f9400c5807c88416676f17c282f5183e2b1b61 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Tue, 4 Oct 2022 09:45:32 +0200 Subject: [PATCH 09/54] added writing of boxes for domain obj test --- tests/domain/domain_obj.py | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/tests/domain/domain_obj.py b/tests/domain/domain_obj.py index ff52525..03b9d3e 100644 --- a/tests/domain/domain_obj.py +++ b/tests/domain/domain_obj.py @@ -1,4 +1,6 @@ import dsmc +import dsmc.octree as oc +import dsmc.writer as wrt import time import numpy as np @@ -45,7 +47,17 @@ def write_partices(positions, it): mass = 6.6422e-26 T = 273.0 n = 2.6e+19 - niter = 100 + niter = 1 + + # trees + tree_inner = oc.Octree() + tree_outer = oc.Octree() + + inner_pos = np.array([[-0.5, -0.5, -0.5], [0.5, 0.5, 0.5]]) + outer_pos = np.array([[-1.0, -1.0, -1.0], [1.0, 1.0, 1.0]]) + + tree_inner.build(inner_pos) + tree_outer.build(outer_pos) # setup solver solver.set_domain(domain) @@ -65,6 +77,9 @@ def write_partices(positions, it): print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) solver.advance(dt) write_partices(solver.particles.Pos, it) + + wrt.write_buttom_leafs(tree_inner, "inner_box.vtu") + wrt.write_buttom_leafs(tree_outer, "outer_box.vtu") print("") print("--- %s seconds ---" % (time.time() - start_time)) From 4578a629d3a2342f073d42aff065e99c8d05db83 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Tue, 4 Oct 2022 09:47:52 +0200 Subject: [PATCH 10/54] fixed bug in obj test --- tests/domain/domain_obj.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/domain/domain_obj.py b/tests/domain/domain_obj.py index 03b9d3e..22c8ad1 100644 --- a/tests/domain/domain_obj.py +++ b/tests/domain/domain_obj.py @@ -47,7 +47,7 @@ def write_partices(positions, it): mass = 6.6422e-26 T = 273.0 n = 2.6e+19 - niter = 1 + niter = 200 # trees tree_inner = oc.Octree() @@ -75,7 +75,7 @@ def write_partices(positions, it): for it in range(niter): print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) - solver.advance(dt) + solver.advance(dt, octree=False) write_partices(solver.particles.Pos, it) wrt.write_buttom_leafs(tree_inner, "inner_box.vtu") From 5ec20f6a2b7cff515a1760258d523f3017e2ccb9 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Tue, 4 Oct 2022 10:08:08 +0200 Subject: [PATCH 11/54] updated boundary --- dsmc/dsmc.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index 0636efd..7a3ccd3 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -7,14 +7,14 @@ @njit def _push(velocities, positions, dt): + old_positions = np.copy(positions) for p in prange(len(positions)): positions[p] = positions[p] + velocities[p]*dt - return positions + return (velocities, positions, old_positions) @njit -def _boundary(velocities, positions, domain, boundary_conds): +def _boundary(velocities, positions, old_positions, domain, boundary_conds): kept_parts = np.ones(positions.shape[0], dtype=np.uint) - old_positions = np.copy(positions) for p in prange(len(positions)): while not oc._is_inside(positions[p], domain) and kept_parts[p]: @@ -63,16 +63,18 @@ def _check_positions(velocities, positions, old_positions, domain): p = 0 new_velocities = np.empty((N, 3)) new_positions = np.empty((N, 3)) + new_old_positions = np.empty((N, 3)) for i in prange(positions.shape[0]): if kept_parts[i] == 1: new_velocities[p] = velocities[i] new_positions[p] = positions[i] + new_old_positions[p] = old_positions[i] p += 1 else: continue - return (new_velocities, new_positions) + return (new_velocities, new_positions, new_old_positions) @njit def _check_created_particles(velocities, positions, obj): @@ -243,10 +245,9 @@ def advance(self, dt, collisions=True, octree=True): self.octree.build(self.particles.Pos) if collisions and octree: self.particles.VelPos = (self._update_velocities(dt), self.particles.Pos) - old_positions = np.copy(self.particles.Pos) - positions = _push(self.particles.Vel, self.particles.Pos, dt) - self.particles.VelPos = _check_positions(self.particles.Vel, positions, old_positions, self.domain) - velocities, positions, old_positions = _boundary(self.particles.Vel, self.particles.Pos, self.domain, self.boundary_conds) + velocities, positions, old_positions = _push(self.particles.Vel, self.particles.Pos, dt) + velocities, positions, old_positions = _check_positions(velocities, positions, old_positions, self.domain) + velocities, positions, old_positions = _boundary(velocities, positions, old_positions, self.domain, self.boundary_conds) for obj in self.objects: velocities, positions, old_positions = _object(velocities, positions, old_positions, obj) From e1f2bb9dbd8341cb697d68127604d4706f9e726b Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Tue, 4 Oct 2022 10:08:52 +0200 Subject: [PATCH 12/54] fixed bug in hypersonic test case --- examples/hypersonic_flow.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/hypersonic_flow.py b/examples/hypersonic_flow.py index c19a89a..c6f9264 100644 --- a/examples/hypersonic_flow.py +++ b/examples/hypersonic_flow.py @@ -45,7 +45,7 @@ solver.add_object(obj) # create particles - solver.create_particles(domain, T, n) + solver.create_particles(domain, T, n, u) # start timing start_time = time.time() From e66ae2ddf5ce001013c5fc62df8ff6bc28782522 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Tue, 4 Oct 2022 10:47:20 +0200 Subject: [PATCH 13/54] fixed bug --- dsmc/dsmc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index 7a3ccd3..0148ce1 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -34,7 +34,7 @@ def _boundary(velocities, positions, old_positions, domain, boundary_conds): elif boundary_conds[i][1] == 1 or boundary_conds[i][0] == 2: kept_parts[p] = 0 - N = sum(kept_parts) + N = int(sum(kept_parts)) p = 0 new_velocities = np.empty((N, 3)) new_positions = np.empty((N, 3)) From c71a39ac77fecd26d7b4065ffd1d686684a56ef9 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Tue, 4 Oct 2022 10:47:50 +0200 Subject: [PATCH 14/54] fixed bug in test case --- examples/hypersonic_flow.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/examples/hypersonic_flow.py b/examples/hypersonic_flow.py index c6f9264..f498391 100644 --- a/examples/hypersonic_flow.py +++ b/examples/hypersonic_flow.py @@ -10,12 +10,13 @@ domain = [(-3.0, 3.0), (-1.5, 1.5), (-0.025, 0.025)] obj = [(-0.25, 0.25), (-0.25, 0.25), (-0.5, 0.5)] dt = 1e-6 - w = 2.3e+14 + w = 0.2 * 2.3e+15 mass = 6.6422e-26 T = 273.0 n = 2.6e+19 u = np.array([0.0, -3043.0, 0.0]) - niter = 1 + niter = 2001 + niter_sample = 10 # set up diagnostics values x0 = domain[0][0] @@ -33,13 +34,13 @@ # set up boundary solver.set_bc_type("xmin", "inflow") solver.set_bc_type("xmax", "inflow") - solver.set_bc_type("ymin", "inflow") + solver.set_bc_type("ymax", "inflow") - solver.set_bc_type("ymax", "open") + solver.set_bc_type("ymin", "open") solver.set_bc_values("xmin", T, n, u) solver.set_bc_values("xmax", T, n, u) - solver.set_bc_values("ymin", T, n, u) + solver.set_bc_values("ymax", T, n, u) # set object solver.add_object(obj) @@ -53,9 +54,15 @@ for it in range(niter): print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) solver.advance(dt) + + for it in range(niter_sample): + print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) + solver.advance(dt) boxes, values = dia.analyse_2d(solver.particles.Pos, x0, x1, y0, y1, Nx, Ny) wrt.write_planar(boxes, values, "hypersonic_{}.vtu".format(it)) + wrt.write_buttom_leafs(solver.octree) + print("") print("--- %s seconds ---" % (time.time() - start_time)) print('done') \ No newline at end of file From 73f31f2a125aadcd52bc1281854e800f71733ee5 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Tue, 4 Oct 2022 11:03:44 +0200 Subject: [PATCH 15/54] removed duplicate import --- dsmc/diagnostics.py | 1 - 1 file changed, 1 deletion(-) diff --git a/dsmc/diagnostics.py b/dsmc/diagnostics.py index 7a4a14b..9782cc7 100644 --- a/dsmc/diagnostics.py +++ b/dsmc/diagnostics.py @@ -1,6 +1,5 @@ from numba import njit import numpy as np -from numba import njit from . import particles as prt from . import octree as oc from . import common as com From b60c88552630baae636e25d9f8f772f2c05d10c8 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Tue, 4 Oct 2022 11:18:54 +0200 Subject: [PATCH 16/54] added mesh2d --- dsmc/mesh/__init__.py | 2 ++ dsmc/mesh/mesh2d.py | 39 +++++++++++++++++++++++++++ setup.py | 2 +- tests/unit/main.py | 1 + tests/unit/test_dsmc/mesh/__init__.py | 2 ++ tests/unit/test_dsmc/mesh/mesh2d.py | 7 +++++ 6 files changed, 52 insertions(+), 1 deletion(-) create mode 100644 dsmc/mesh/__init__.py create mode 100644 dsmc/mesh/mesh2d.py create mode 100644 tests/unit/test_dsmc/mesh/__init__.py create mode 100644 tests/unit/test_dsmc/mesh/mesh2d.py diff --git a/dsmc/mesh/__init__.py b/dsmc/mesh/__init__.py new file mode 100644 index 0000000..633f866 --- /dev/null +++ b/dsmc/mesh/__init__.py @@ -0,0 +1,2 @@ +# -*- coding: utf-8 -*- + diff --git a/dsmc/mesh/mesh2d.py b/dsmc/mesh/mesh2d.py new file mode 100644 index 0000000..187498d --- /dev/null +++ b/dsmc/mesh/mesh2d.py @@ -0,0 +1,39 @@ +from enum import Enum +import numpy as np + +def _get_cell_id(val1, val2, n_cells1, n_cells2, min1, min2, cell_size): + #uint cell_id1, cell_id2; + + if (val1 < min1): + return (False, 0) + elif (val1 > (min1 + n_cells1 * cell_size)): + return (False, 0) + else: + cell_id1 = np.floor((val1 - min1) / cell_size) + + if (val2 < min2): + return (False, 0) + elif (val2 > (min2 + n_cells2 * cell_size)): + return (False, 0) + else: + cell_id2 = np.floor((val2 - min2) / cell_size) + + return (True, cell_id2 * n_cells1 + cell_id1) + + +class Plane(Enum): + XY = 0 + YZ = 1 + XZ = 2 + +class Mesh2d: + def __init__(self): + self.clear() + + def clear(self): + cell_size = 1.0 + min1 = 0.0 + min2 = 0.0 + n_cells1 = 1 + n_cells2 = 1 + plane = Plane.XY \ No newline at end of file diff --git a/setup.py b/setup.py index 1064716..2b71d95 100644 --- a/setup.py +++ b/setup.py @@ -4,7 +4,7 @@ version='0.8.0', author='Leo Basov', python_requires='>=3.10, <4', - packages=["dsmc", "dsmc.writer"], + packages=["dsmc", "dsmc.writer", "dsmc.mesh"], install_requires=[ 'numpy', 'llvmlite', diff --git a/tests/unit/main.py b/tests/unit/main.py index df9783a..3a94703 100644 --- a/tests/unit/main.py +++ b/tests/unit/main.py @@ -7,6 +7,7 @@ from test_dsmc.octree import * from test_dsmc.diagnostics import * from test_dsmc.common import * +from test_dsmc.mesh.mesh2d import * if __name__ == '__main__': unittest.main() diff --git a/tests/unit/test_dsmc/mesh/__init__.py b/tests/unit/test_dsmc/mesh/__init__.py new file mode 100644 index 0000000..633f866 --- /dev/null +++ b/tests/unit/test_dsmc/mesh/__init__.py @@ -0,0 +1,2 @@ +# -*- coding: utf-8 -*- + diff --git a/tests/unit/test_dsmc/mesh/mesh2d.py b/tests/unit/test_dsmc/mesh/mesh2d.py new file mode 100644 index 0000000..674fba8 --- /dev/null +++ b/tests/unit/test_dsmc/mesh/mesh2d.py @@ -0,0 +1,7 @@ +import numpy as np +import dsmc.mesh.mesh2d as msh +import unittest + +class TestMesh2(unittest.TestCase): + def test_pass(self): + self.assertTrue(False) \ No newline at end of file From 77226cbc8e429019abdfc2bb3c5a80902f297a31 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Tue, 4 Oct 2022 11:21:51 +0200 Subject: [PATCH 17/54] fixed bug in mesh constructor --- dsmc/mesh/mesh2d.py | 12 ++++++------ tests/unit/test_dsmc/mesh/mesh2d.py | 11 +++++++++-- 2 files changed, 15 insertions(+), 8 deletions(-) diff --git a/dsmc/mesh/mesh2d.py b/dsmc/mesh/mesh2d.py index 187498d..a26935b 100644 --- a/dsmc/mesh/mesh2d.py +++ b/dsmc/mesh/mesh2d.py @@ -31,9 +31,9 @@ def __init__(self): self.clear() def clear(self): - cell_size = 1.0 - min1 = 0.0 - min2 = 0.0 - n_cells1 = 1 - n_cells2 = 1 - plane = Plane.XY \ No newline at end of file + self.cell_size = 1.0 + self.min1 = 0.0 + self.min2 = 0.0 + self.n_cells1 = 1 + self.n_cells2 = 1 + self.plane = Plane.XY \ No newline at end of file diff --git a/tests/unit/test_dsmc/mesh/mesh2d.py b/tests/unit/test_dsmc/mesh/mesh2d.py index 674fba8..bf85c24 100644 --- a/tests/unit/test_dsmc/mesh/mesh2d.py +++ b/tests/unit/test_dsmc/mesh/mesh2d.py @@ -3,5 +3,12 @@ import unittest class TestMesh2(unittest.TestCase): - def test_pass(self): - self.assertTrue(False) \ No newline at end of file + def test_constructor(self): + mesh = msh.Mesh2d() + + self.assertEqual(1.0, mesh.cell_size) + self.assertEqual(0.0, mesh.min1) + self.assertEqual(0.0, mesh.min2) + self.assertEqual(1, mesh.n_cells1) + self.assertEqual(1, mesh.n_cells2) + self.assertEqual(msh.Plane.XY, mesh.plane) \ No newline at end of file From 97ea11de5a9892b31f4c3cc0ef37759159c346ce Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Tue, 4 Oct 2022 11:26:28 +0200 Subject: [PATCH 18/54] implemented test__get_cell_id1 --- dsmc/mesh/mesh2d.py | 2 -- tests/unit/test_dsmc/mesh/mesh2d.py | 22 ++++++++++++++++++++-- 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/dsmc/mesh/mesh2d.py b/dsmc/mesh/mesh2d.py index a26935b..9898526 100644 --- a/dsmc/mesh/mesh2d.py +++ b/dsmc/mesh/mesh2d.py @@ -2,8 +2,6 @@ import numpy as np def _get_cell_id(val1, val2, n_cells1, n_cells2, min1, min2, cell_size): - #uint cell_id1, cell_id2; - if (val1 < min1): return (False, 0) elif (val1 > (min1 + n_cells1 * cell_size)): diff --git a/tests/unit/test_dsmc/mesh/mesh2d.py b/tests/unit/test_dsmc/mesh/mesh2d.py index bf85c24..b3a2f5e 100644 --- a/tests/unit/test_dsmc/mesh/mesh2d.py +++ b/tests/unit/test_dsmc/mesh/mesh2d.py @@ -1,4 +1,3 @@ -import numpy as np import dsmc.mesh.mesh2d as msh import unittest @@ -11,4 +10,23 @@ def test_constructor(self): self.assertEqual(0.0, mesh.min2) self.assertEqual(1, mesh.n_cells1) self.assertEqual(1, mesh.n_cells2) - self.assertEqual(msh.Plane.XY, mesh.plane) \ No newline at end of file + self.assertEqual(msh.Plane.XY, mesh.plane) + + def test__get_cell_id1(self): + val1_f = -1.0 + val1_t = 0.1 + val2_f = -1.0 + val2_t = 0.1 + n_cells1 = 10 + n_cells2 = 10 + min1 = -0.5 + min2 = -0.5 + cell_size = 0.1 + + res1_f = msh._get_cell_id(val1_t, val2_f, n_cells1, n_cells2, min1, min2, cell_size) + res2_f = msh._get_cell_id(val1_f, val2_t, n_cells1, n_cells2, min1, min2, cell_size) + res3_t = msh._get_cell_id(val1_t, val2_t, n_cells1, n_cells2, min1, min2, cell_size) + + self.assertFalse(res1_f[0]) + self.assertFalse(res2_f[0]) + self.assertTrue(res3_t[0]) \ No newline at end of file From d2dc6c36ad50255551ec7db28d5b9e251f40dc5c Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Tue, 4 Oct 2022 11:41:20 +0200 Subject: [PATCH 19/54] added unit test --- dsmc/mesh/mesh2d.py | 6 +++--- tests/unit/test_dsmc/mesh/mesh2d.py | 20 +++++++++++++++++++- 2 files changed, 22 insertions(+), 4 deletions(-) diff --git a/dsmc/mesh/mesh2d.py b/dsmc/mesh/mesh2d.py index 9898526..102c196 100644 --- a/dsmc/mesh/mesh2d.py +++ b/dsmc/mesh/mesh2d.py @@ -1,5 +1,5 @@ from enum import Enum -import numpy as np +import math def _get_cell_id(val1, val2, n_cells1, n_cells2, min1, min2, cell_size): if (val1 < min1): @@ -7,14 +7,14 @@ def _get_cell_id(val1, val2, n_cells1, n_cells2, min1, min2, cell_size): elif (val1 > (min1 + n_cells1 * cell_size)): return (False, 0) else: - cell_id1 = np.floor((val1 - min1) / cell_size) + cell_id1 = math.floor((val1 - min1) / cell_size) if (val2 < min2): return (False, 0) elif (val2 > (min2 + n_cells2 * cell_size)): return (False, 0) else: - cell_id2 = np.floor((val2 - min2) / cell_size) + cell_id2 = math.floor((val2 - min2) / cell_size) return (True, cell_id2 * n_cells1 + cell_id1) diff --git a/tests/unit/test_dsmc/mesh/mesh2d.py b/tests/unit/test_dsmc/mesh/mesh2d.py index b3a2f5e..1337e82 100644 --- a/tests/unit/test_dsmc/mesh/mesh2d.py +++ b/tests/unit/test_dsmc/mesh/mesh2d.py @@ -29,4 +29,22 @@ def test__get_cell_id1(self): self.assertFalse(res1_f[0]) self.assertFalse(res2_f[0]) - self.assertTrue(res3_t[0]) \ No newline at end of file + self.assertTrue(res3_t[0]) + + def test__get_cell_id2(self): + val1 = (-0.45, -0.5) + val2 = (-0.35, -0.35) + n_cells1 = 10 + n_cells2 = 10 + min1 = -0.5 + min2 = -0.5 + cell_size = 0.1 + + res1 = msh._get_cell_id(val1[0], val1[1], n_cells1, n_cells2, min1, min2, cell_size) + res2 = msh._get_cell_id(val2[0], val2[1], n_cells1, n_cells2, min1, min2, cell_size) + + self.assertTrue(res1[0]) + self.assertTrue(res2[0]) + + self.assertEqual(0, res1[1]) + self.assertEqual(11, res2[1]) \ No newline at end of file From bd1367869333f356b65560d71db5c5f32a518bc8 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Tue, 4 Oct 2022 11:58:49 +0200 Subject: [PATCH 20/54] implemented get_id --- dsmc/mesh/mesh2d.py | 13 ++++++++++- tests/unit/test_dsmc/mesh/mesh2d.py | 35 ++++++++++++++++++++++++++++- 2 files changed, 46 insertions(+), 2 deletions(-) diff --git a/dsmc/mesh/mesh2d.py b/dsmc/mesh/mesh2d.py index 102c196..dfaef8c 100644 --- a/dsmc/mesh/mesh2d.py +++ b/dsmc/mesh/mesh2d.py @@ -1,6 +1,8 @@ from enum import Enum import math +from numba import njit +@njit def _get_cell_id(val1, val2, n_cells1, n_cells2, min1, min2, cell_size): if (val1 < min1): return (False, 0) @@ -34,4 +36,13 @@ def clear(self): self.min2 = 0.0 self.n_cells1 = 1 self.n_cells2 = 1 - self.plane = Plane.XY \ No newline at end of file + self.plane = Plane.XY + + def get_cell_id(self, position): + match self.plane: + case Plane.XY: + return _get_cell_id(position[0], position[1], self.n_cells1, self.n_cells2, self.min1, self.min2, self.cell_size) + case Plane.YZ: + return _get_cell_id(position[1], position[2], self.n_cells1, self.n_cells2, self.min1, self.min2, self.cell_size) + case Plane.XZ: + return _get_cell_id(position[0], position[2], self.n_cells1, self.n_cells2, self.min1, self.min2, self.cell_size) \ No newline at end of file diff --git a/tests/unit/test_dsmc/mesh/mesh2d.py b/tests/unit/test_dsmc/mesh/mesh2d.py index 1337e82..0076a6e 100644 --- a/tests/unit/test_dsmc/mesh/mesh2d.py +++ b/tests/unit/test_dsmc/mesh/mesh2d.py @@ -1,5 +1,6 @@ import dsmc.mesh.mesh2d as msh import unittest +import numpy as np class TestMesh2(unittest.TestCase): def test_constructor(self): @@ -47,4 +48,36 @@ def test__get_cell_id2(self): self.assertTrue(res2[0]) self.assertEqual(0, res1[1]) - self.assertEqual(11, res2[1]) \ No newline at end of file + self.assertEqual(11, res2[1]) + + def test_get_cell_id(self): + mesh = msh.Mesh2d() + + mesh.n_cells1 = 10 + mesh.n_cells2 = 10 + mesh.min1 = -0.5 + mesh.min2 = -0.5 + mesh.cell_size = 0.1 + + pos = np.array([-0.45, -0.35, -0.25]) + + # XY + mesh.plane = msh.Plane.XY + res = mesh.get_cell_id(pos) + + self.assertTrue(res[0]) + self.assertEqual(10, res[1]) + + # YZ + mesh.plane = msh.Plane.YZ + res = mesh.get_cell_id(pos) + + self.assertTrue(res[0]) + self.assertEqual(21, res[1]) + + # XZ + mesh.plane = msh.Plane.XZ + res = mesh.get_cell_id(pos) + + self.assertTrue(res[0]) + self.assertEqual(20, res[1]) \ No newline at end of file From 5d2da5bb727099c7d7f6440498c2b6dbede1f2d4 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Tue, 4 Oct 2022 13:19:52 +0200 Subject: [PATCH 21/54] implemented _sort --- dsmc/mesh/mesh2d.py | 11 +++++++++++ tests/unit/test_dsmc/mesh/mesh2d.py | 26 +++++++++++++++++++++++++- 2 files changed, 36 insertions(+), 1 deletion(-) diff --git a/dsmc/mesh/mesh2d.py b/dsmc/mesh/mesh2d.py index dfaef8c..cefd02d 100644 --- a/dsmc/mesh/mesh2d.py +++ b/dsmc/mesh/mesh2d.py @@ -1,6 +1,7 @@ from enum import Enum import math from numba import njit +import numpy as np @njit def _get_cell_id(val1, val2, n_cells1, n_cells2, min1, min2, cell_size): @@ -20,6 +21,16 @@ def _get_cell_id(val1, val2, n_cells1, n_cells2, min1, min2, cell_size): return (True, cell_id2 * n_cells1 + cell_id1) +@njit +def _sort(values1, values2, n_cells1, n_cells2, min1, min2, cell_size): + inside = np.empty((len(values1)), np.bool_) + ids = np.empty((len(values1)), np.int_) + + for i in range(len(values1)): + inside[i], ids[i] = _get_cell_id(values1[i], values2[i], n_cells1, n_cells2, min1, min2, cell_size) + + return (inside, ids) + class Plane(Enum): XY = 0 diff --git a/tests/unit/test_dsmc/mesh/mesh2d.py b/tests/unit/test_dsmc/mesh/mesh2d.py index 0076a6e..d7ae39d 100644 --- a/tests/unit/test_dsmc/mesh/mesh2d.py +++ b/tests/unit/test_dsmc/mesh/mesh2d.py @@ -80,4 +80,28 @@ def test_get_cell_id(self): res = mesh.get_cell_id(pos) self.assertTrue(res[0]) - self.assertEqual(20, res[1]) \ No newline at end of file + self.assertEqual(20, res[1]) + + def test__sort(self): + mesh = msh.Mesh2d() + + mesh.n_cells1 = 10 + mesh.n_cells2 = 10 + mesh.min1 = -0.5 + mesh.min2 = -0.5 + mesh.cell_size = 0.1 + + pos1 = np.array([-0.45, -0.35, 0.0]) + pos2 = np.array([-0.45, -3.50, 0.0]) + positions = np.array([pos1, pos2]) + + values1 = np.array([positions[i][0] for i in range(len(positions))]) + values2 = np.array([positions[i][1] for i in range(len(positions))]) + + inside, ids = msh._sort(values1, values2, mesh.n_cells1, mesh.n_cells2, mesh.min1, mesh.min2, mesh.cell_size) + + self.assertTrue(inside[0]) + self.assertFalse(inside[1]) + + self.assertEqual(10, ids[0]) + self.assertEqual(0, ids[1]) \ No newline at end of file From 92e0deda2a31c2f8eafd4b9d6344162344411240 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Wed, 5 Oct 2022 09:27:22 +0200 Subject: [PATCH 22/54] added multiple cell size for 2d mesh --- dsmc/mesh/mesh2d.py | 23 ++++++++++++----------- tests/unit/test_dsmc/mesh/mesh2d.py | 27 ++++++++++++++++----------- 2 files changed, 28 insertions(+), 22 deletions(-) diff --git a/dsmc/mesh/mesh2d.py b/dsmc/mesh/mesh2d.py index cefd02d..b07650c 100644 --- a/dsmc/mesh/mesh2d.py +++ b/dsmc/mesh/mesh2d.py @@ -4,30 +4,30 @@ import numpy as np @njit -def _get_cell_id(val1, val2, n_cells1, n_cells2, min1, min2, cell_size): +def _get_cell_id(val1, val2, n_cells1, n_cells2, min1, min2, cell_size1, cell_size2): if (val1 < min1): return (False, 0) - elif (val1 > (min1 + n_cells1 * cell_size)): + elif (val1 > (min1 + n_cells1 * cell_size1)): return (False, 0) else: - cell_id1 = math.floor((val1 - min1) / cell_size) + cell_id1 = math.floor((val1 - min1) / cell_size1) if (val2 < min2): return (False, 0) - elif (val2 > (min2 + n_cells2 * cell_size)): + elif (val2 > (min2 + n_cells2 * cell_size2)): return (False, 0) else: - cell_id2 = math.floor((val2 - min2) / cell_size) + cell_id2 = math.floor((val2 - min2) / cell_size2) return (True, cell_id2 * n_cells1 + cell_id1) @njit -def _sort(values1, values2, n_cells1, n_cells2, min1, min2, cell_size): +def _sort(values1, values2, n_cells1, n_cells2, min1, min2, cell_size1, cell_size2): inside = np.empty((len(values1)), np.bool_) ids = np.empty((len(values1)), np.int_) for i in range(len(values1)): - inside[i], ids[i] = _get_cell_id(values1[i], values2[i], n_cells1, n_cells2, min1, min2, cell_size) + inside[i], ids[i] = _get_cell_id(values1[i], values2[i], n_cells1, n_cells2, min1, min2, cell_size1, cell_size2) return (inside, ids) @@ -42,7 +42,8 @@ def __init__(self): self.clear() def clear(self): - self.cell_size = 1.0 + self.cell_size1 = 1.0 + self.cell_size2 = 1.0 self.min1 = 0.0 self.min2 = 0.0 self.n_cells1 = 1 @@ -52,8 +53,8 @@ def clear(self): def get_cell_id(self, position): match self.plane: case Plane.XY: - return _get_cell_id(position[0], position[1], self.n_cells1, self.n_cells2, self.min1, self.min2, self.cell_size) + return _get_cell_id(position[0], position[1], self.n_cells1, self.n_cells2, self.min1, self.min2, self.cell_size1, self.cell_size2) case Plane.YZ: - return _get_cell_id(position[1], position[2], self.n_cells1, self.n_cells2, self.min1, self.min2, self.cell_size) + return _get_cell_id(position[1], position[2], self.n_cells1, self.n_cells2, self.min1, self.min2, self.cell_size1, self.cell_size2) case Plane.XZ: - return _get_cell_id(position[0], position[2], self.n_cells1, self.n_cells2, self.min1, self.min2, self.cell_size) \ No newline at end of file + return _get_cell_id(position[0], position[2], self.n_cells1, self.n_cells2, self.min1, self.min2, self.cell_size1, self.cell_size2) \ No newline at end of file diff --git a/tests/unit/test_dsmc/mesh/mesh2d.py b/tests/unit/test_dsmc/mesh/mesh2d.py index d7ae39d..13a5c04 100644 --- a/tests/unit/test_dsmc/mesh/mesh2d.py +++ b/tests/unit/test_dsmc/mesh/mesh2d.py @@ -6,7 +6,8 @@ class TestMesh2(unittest.TestCase): def test_constructor(self): mesh = msh.Mesh2d() - self.assertEqual(1.0, mesh.cell_size) + self.assertEqual(1.0, mesh.cell_size1) + self.assertEqual(1.0, mesh.cell_size2) self.assertEqual(0.0, mesh.min1) self.assertEqual(0.0, mesh.min2) self.assertEqual(1, mesh.n_cells1) @@ -22,11 +23,12 @@ def test__get_cell_id1(self): n_cells2 = 10 min1 = -0.5 min2 = -0.5 - cell_size = 0.1 + cell_size1 = 0.1 + cell_size2 = 0.1 - res1_f = msh._get_cell_id(val1_t, val2_f, n_cells1, n_cells2, min1, min2, cell_size) - res2_f = msh._get_cell_id(val1_f, val2_t, n_cells1, n_cells2, min1, min2, cell_size) - res3_t = msh._get_cell_id(val1_t, val2_t, n_cells1, n_cells2, min1, min2, cell_size) + res1_f = msh._get_cell_id(val1_t, val2_f, n_cells1, n_cells2, min1, min2, cell_size1, cell_size2) + res2_f = msh._get_cell_id(val1_f, val2_t, n_cells1, n_cells2, min1, min2, cell_size1, cell_size2) + res3_t = msh._get_cell_id(val1_t, val2_t, n_cells1, n_cells2, min1, min2, cell_size1, cell_size2) self.assertFalse(res1_f[0]) self.assertFalse(res2_f[0]) @@ -39,10 +41,11 @@ def test__get_cell_id2(self): n_cells2 = 10 min1 = -0.5 min2 = -0.5 - cell_size = 0.1 + cell_size1 = 0.1 + cell_size2 = 0.1 - res1 = msh._get_cell_id(val1[0], val1[1], n_cells1, n_cells2, min1, min2, cell_size) - res2 = msh._get_cell_id(val2[0], val2[1], n_cells1, n_cells2, min1, min2, cell_size) + res1 = msh._get_cell_id(val1[0], val1[1], n_cells1, n_cells2, min1, min2, cell_size1, cell_size2) + res2 = msh._get_cell_id(val2[0], val2[1], n_cells1, n_cells2, min1, min2, cell_size1, cell_size2) self.assertTrue(res1[0]) self.assertTrue(res2[0]) @@ -57,7 +60,8 @@ def test_get_cell_id(self): mesh.n_cells2 = 10 mesh.min1 = -0.5 mesh.min2 = -0.5 - mesh.cell_size = 0.1 + mesh.cell_size1 = 0.1 + mesh.cell_size2 = 0.1 pos = np.array([-0.45, -0.35, -0.25]) @@ -89,7 +93,8 @@ def test__sort(self): mesh.n_cells2 = 10 mesh.min1 = -0.5 mesh.min2 = -0.5 - mesh.cell_size = 0.1 + mesh.cell_size1 = 0.1 + mesh.cell_size2 = 0.1 pos1 = np.array([-0.45, -0.35, 0.0]) pos2 = np.array([-0.45, -3.50, 0.0]) @@ -98,7 +103,7 @@ def test__sort(self): values1 = np.array([positions[i][0] for i in range(len(positions))]) values2 = np.array([positions[i][1] for i in range(len(positions))]) - inside, ids = msh._sort(values1, values2, mesh.n_cells1, mesh.n_cells2, mesh.min1, mesh.min2, mesh.cell_size) + inside, ids = msh._sort(values1, values2, mesh.n_cells1, mesh.n_cells2, mesh.min1, mesh.min2, mesh.cell_size1, mesh.cell_size2) self.assertTrue(inside[0]) self.assertFalse(inside[1]) From 99ec53410305f2118246553c8797d1ecaf743f22 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Wed, 5 Oct 2022 09:46:18 +0200 Subject: [PATCH 23/54] implemented sort --- dsmc/mesh/mesh2d.py | 22 +++++++++++++++++++ tests/unit/test_dsmc/mesh/mesh2d.py | 33 ++++++++++++++++++++++++++++- 2 files changed, 54 insertions(+), 1 deletion(-) diff --git a/dsmc/mesh/mesh2d.py b/dsmc/mesh/mesh2d.py index b07650c..a409084 100644 --- a/dsmc/mesh/mesh2d.py +++ b/dsmc/mesh/mesh2d.py @@ -49,6 +49,28 @@ def clear(self): self.n_cells1 = 1 self.n_cells2 = 1 self.plane = Plane.XY + self.cells = None # this holds an 2d array, cells[] : cells, cells[][] position ids in cell + + def sort(self, position): + self.cells = [[] for _ in range(self.n_cells1 * self.n_cells2)] + + match self.plane: + case Plane.XY: + positions1 = np.array([position[i][0] for i in range(len(position))]) + positions2 = np.array([position[i][1] for i in range(len(position))]) + case Plane.YZ: + positions1 = np.array([position[i][1] for i in range(len(position))]) + positions2 = np.array([position[i][2] for i in range(len(position))]) + case Plane.XZ: + positions1 = np.array([position[i][0] for i in range(len(position))]) + positions2 = np.array([position[i][2] for i in range(len(position))]) + + inside, cell_ids = _sort(positions1, positions2, self.n_cells1, self.n_cells2, self.min1, self.min2, self.cell_size1, self.cell_size2) + sorted_ids = np.argsort(cell_ids) + + for idx in sorted_ids: + if inside[idx]: + self.cells[cell_ids[idx]].append(idx) def get_cell_id(self, position): match self.plane: diff --git a/tests/unit/test_dsmc/mesh/mesh2d.py b/tests/unit/test_dsmc/mesh/mesh2d.py index 13a5c04..e993b28 100644 --- a/tests/unit/test_dsmc/mesh/mesh2d.py +++ b/tests/unit/test_dsmc/mesh/mesh2d.py @@ -109,4 +109,35 @@ def test__sort(self): self.assertFalse(inside[1]) self.assertEqual(10, ids[0]) - self.assertEqual(0, ids[1]) \ No newline at end of file + self.assertEqual(0, ids[1]) + + def test_sort(self): + mesh = msh.Mesh2d() + + mesh.n_cells1 = 10 + mesh.n_cells2 = 10 + mesh.min1 = -0.5 + mesh.min2 = -0.5 + mesh.cell_size1 = 0.1 + mesh.cell_size2 = 0.1 + + positions = [] + + positions.append([-0.45, -0.45, 0]) + positions.append([-0.45, -0.45, 0]) + + positions.append([0.45, 0.45, 0]) + + mesh.sort(positions) + + self.assertEqual(len(mesh.cells), mesh.n_cells1 * mesh.n_cells2) + self.assertEqual(len(mesh.cells[0]), 2) + self.assertEqual(len(mesh.cells[99]), 1) + + self.assertEqual(mesh.cells[0][0], 0) + self.assertEqual(mesh.cells[0][1], 1) + self.assertEqual(mesh.cells[99][0], 2) + + for i in range(len(mesh.cells)): + if i != 0 and i != 99: + self.assertEqual(0, len(mesh.cells[i])) \ No newline at end of file From f7f89de5b42acdeef147e6e7ca2cd2a7ce0eb0a4 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Wed, 5 Oct 2022 10:10:49 +0200 Subject: [PATCH 24/54] updated diagnostics --- dsmc/diagnostics.py | 98 ++++++++----------------------------- examples/hypersonic_flow.py | 30 +++++++----- 2 files changed, 39 insertions(+), 89 deletions(-) diff --git a/dsmc/diagnostics.py b/dsmc/diagnostics.py index 9782cc7..6f5a8ac 100644 --- a/dsmc/diagnostics.py +++ b/dsmc/diagnostics.py @@ -1,92 +1,36 @@ -from numba import njit import numpy as np from . import particles as prt from . import octree as oc -from . import common as com - -@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 - boxes = np.empty((Nx*Ny, 2, 2)) - y = y0 - - for i in range(Ny): - x = x0 - 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 (permutations, leafs, boxes) +#from . import common as com 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 calc_vals(self, positions, velocities, ids, box): + #V = com.get_V(box) + self.number_elements = len(ids) + -@njit -def sort_flat(positions, box, Nx, Ny): - vals = np.zeros((Nx, Ny)) - dx = (box[0][1] - box[0][0]) / Nx - dy = (box[1][1] - box[1][0]) / Ny +def analyse_2d(positions, velocities, mesh2d, h): + boxes = np.empty((len(mesh2d.cells), 3, 2)) + values = [Values() for _ in range(len(mesh2d.cells))] - for p in range(positions.shape[0]): - pos = positions[p] - X = box[0][0] - Y = box[1][0] + for y in range(mesh2d.n_cells2): + for x in range(mesh2d.n_cells1): + N = x + y * mesh2d.n_cells1 + + boxes[N][0][0] = mesh2d.min1 + mesh2d.cell_size1*x + boxes[N][0][1] = mesh2d.min1 + mesh2d.cell_size1*(x + 1) + boxes[N][1][0] = mesh2d.min2 + mesh2d.cell_size2*y + boxes[N][1][1] = mesh2d.min2 + mesh2d.cell_size2*(y + 1) + boxes[N][2][0] = 0.0 + boxes[N][2][1] = h + + values[N].calc_vals(positions, velocities, mesh2d.cells[N], boxes[N]) - for x in range(Nx): - for y in range(Ny): - a = (pos[0] < (X + dx)) and (pos[0] >= X) - b = (pos[1] < (Y + dy)) and (pos[1] >= Y) - - if a and b: + return (boxes, values) - vals[x][y] += 1 - - Y += dy - X += dx - - return vals.astype(np.int_) def sort_bin(positions, axis, Nbin): bins = [[] for _ in range(Nbin)] diff --git a/examples/hypersonic_flow.py b/examples/hypersonic_flow.py index f498391..b1f974c 100644 --- a/examples/hypersonic_flow.py +++ b/examples/hypersonic_flow.py @@ -1,6 +1,7 @@ import dsmc import dsmc.writer as wrt import dsmc.diagnostics as dia +import dsmc.mesh.mesh2d as msh2d import time import numpy as np @@ -15,16 +16,20 @@ T = 273.0 n = 2.6e+19 u = np.array([0.0, -3043.0, 0.0]) - niter = 2001 + niter = 1#2001 niter_sample = 10 - # set up diagnostics values - x0 = domain[0][0] - x1 = domain[0][1] - y0 = domain[1][0] - y1 = domain[1][1] - Nx = 100 - Ny = 50 + # set up mesh2 + mesh = msh2d.Mesh2d() + + mesh.n_cells1 = 100 + mesh.n_cells2 = 50 + mesh.min1 = domain[0][0] + mesh.min2 = domain[1][0] + mesh.cell_size1 = 0.06 + mesh.cell_size2 = 0.06 + + h = domain[2][1] - domain[2][0] # setup solver solver.set_domain(domain) @@ -51,14 +56,15 @@ # start timing start_time = time.time() - for it in range(niter): + """for it in range(niter): print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) - solver.advance(dt) + solver.advance(dt)""" for it in range(niter_sample): - print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) + print("iteration {:4}/{}".format(it + 1, niter_sample), end="\r", flush=True) solver.advance(dt) - boxes, values = dia.analyse_2d(solver.particles.Pos, x0, x1, y0, y1, Nx, Ny) + mesh.sort(solver.particles.Pos) + boxes, values = dia.analyse_2d(solver.particles.Pos, solver.particles.Vel, mesh, h) wrt.write_planar(boxes, values, "hypersonic_{}.vtu".format(it)) wrt.write_buttom_leafs(solver.octree) From e8162ca28755a3f3faf2a0555ca2b36044183274 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Wed, 5 Oct 2022 10:24:43 +0200 Subject: [PATCH 25/54] rempaced exception by warning in dsmc --- dsmc/dsmc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index 0148ce1..d99a694 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -232,7 +232,7 @@ def advance(self, dt, collisions=True, octree=True): if self.domain is None: raise Exception("simulation domain not defined") if self.particles.N == 0: - raise Exception("no particles in domain") + print("warning: no particles in domain") if self.w == None: raise Exception("particle weight not set") From 32bc1560dd55941a778de2f6a77132e436273261 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Wed, 5 Oct 2022 12:09:32 +0200 Subject: [PATCH 26/54] updated hypersonic test case --- examples/hypersonic_flow.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/examples/hypersonic_flow.py b/examples/hypersonic_flow.py index b1f974c..aae1028 100644 --- a/examples/hypersonic_flow.py +++ b/examples/hypersonic_flow.py @@ -16,8 +16,8 @@ T = 273.0 n = 2.6e+19 u = np.array([0.0, -3043.0, 0.0]) - niter = 1#2001 - niter_sample = 10 + niter = 10000 + niter_sample = 100 # set up mesh2 mesh = msh2d.Mesh2d() @@ -56,16 +56,18 @@ # start timing start_time = time.time() - """for it in range(niter): + for it in range(niter): print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) - solver.advance(dt)""" + solver.advance(dt) for it in range(niter_sample): print("iteration {:4}/{}".format(it + 1, niter_sample), end="\r", flush=True) solver.advance(dt) - mesh.sort(solver.particles.Pos) - boxes, values = dia.analyse_2d(solver.particles.Pos, solver.particles.Vel, mesh, h) - wrt.write_planar(boxes, values, "hypersonic_{}.vtu".format(it)) + + if (it + 1)%10 == 0: + mesh.sort(solver.particles.Pos) + boxes, values = dia.analyse_2d(solver.particles.Pos, solver.particles.Vel, mesh, h) + wrt.write_planar(boxes, values, "hypersonic_{}.vtu".format(it)) wrt.write_buttom_leafs(solver.octree) From 46713411065e40a57b0c47612fd9a09283cc063e Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Wed, 5 Oct 2022 13:14:57 +0200 Subject: [PATCH 27/54] updated dsmc order and hyper sonic test case --- dsmc/dsmc.py | 13 +++++++------ examples/hypersonic_flow.py | 2 +- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index d99a694..5cca160 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -241,27 +241,28 @@ def advance(self, dt, collisions=True, octree=True): if self.boundary_conds[i][j] == 2: self.particles.inflow(self.mass, self.boundary.T[i][j], self.boundary.u[i][j], self.boundary.n[i][j], self.w, dt, self.domain, i, j) - if octree: - self.octree.build(self.particles.Pos) - if collisions and octree: - self.particles.VelPos = (self._update_velocities(dt), self.particles.Pos) velocities, positions, old_positions = _push(self.particles.Vel, self.particles.Pos, dt) velocities, positions, old_positions = _check_positions(velocities, positions, old_positions, self.domain) velocities, positions, old_positions = _boundary(velocities, positions, old_positions, self.domain, self.boundary_conds) for obj in self.objects: velocities, positions, old_positions = _object(velocities, positions, old_positions, obj) + + if octree: + self.octree.build(positions) + if collisions and octree: + velocities = self._update_velocities(dt, velocities) self.particles.VelPos = (velocities, positions) - def _update_velocities(self, dt): + def _update_velocities(self, dt, velocities): Nleafs : int = len(self.octree.leafs) elem_offsets : np.ndarray = np.array([leaf.elem_offset for leaf in self.octree.leafs], dtype=int) number_elements : np.ndarray = np.array([leaf.number_elements for leaf in self.octree.leafs], dtype=int) number_children : np.ndarray = np.array([leaf.number_children for leaf in self.octree.leafs], dtype=int) cell_boxes : np.ndarray = np.array([box for box in self.octree.cell_boxes]) - 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) + return _update_vels(self.octree.permutations, velocities, 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 = np.zeros(3)): box = np.array(box) diff --git a/examples/hypersonic_flow.py b/examples/hypersonic_flow.py index aae1028..4bf0a79 100644 --- a/examples/hypersonic_flow.py +++ b/examples/hypersonic_flow.py @@ -16,7 +16,7 @@ T = 273.0 n = 2.6e+19 u = np.array([0.0, -3043.0, 0.0]) - niter = 10000 + niter = 1000 niter_sample = 100 # set up mesh2 From c511f8f4d7bcd4b2eddc8659b88a0ed4bba5d4f4 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Thu, 6 Oct 2022 13:13:35 +0200 Subject: [PATCH 28/54] modified particle creation --- dsmc/particles.py | 11 +++++------ tests/unit/test_dsmc/particles.py | 3 ++- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/dsmc/particles.py b/dsmc/particles.py index 421fbff..4516065 100644 --- a/dsmc/particles.py +++ b/dsmc/particles.py @@ -83,7 +83,7 @@ def calc_temperature(velocities, mass): return tot_e / ((3.0/2.0) * len(velocities) * kb) @njit -def calc_positions(x, y, z, N): +def calc_positions(X, N): """ Parameters ---------- @@ -97,9 +97,8 @@ def calc_positions(x, y, z, N): positions = np.empty((N, 3)) for i in range(N): - positions[i][0] = x[0] + np.random.random() * (x[1] - x[0]) - positions[i][1] = y[0] + np.random.random() * (y[1] - y[0]) - positions[i][2] = z[0] + np.random.random() * (z[1] - z[0]) + for j in range(3): + positions[i][j] = X[j][0] + np.random.random() * (X[j][1] - X[j][0]) return positions @@ -136,11 +135,11 @@ def VelPos(self, vel_pos): def create_particles(self, X, mass, T, N, u = np.zeros(3)): if self._N == 0: self._velocities = get_velocities(T, mass, N, u) - self._positions = calc_positions(X[0], X[1], X[2], N) + self._positions = calc_positions(X, N) self._N = N else: self._velocities = np.concatenate((self._velocities, get_velocities(T, mass, N, u))) - self._positions = np.concatenate((self._positions, calc_positions(X[0], X[1], X[2], N))) + self._positions = np.concatenate((self._positions, calc_positions(X, N))) self._N += N diff --git a/tests/unit/test_dsmc/particles.py b/tests/unit/test_dsmc/particles.py index 700c46f..0bed9aa 100644 --- a/tests/unit/test_dsmc/particles.py +++ b/tests/unit/test_dsmc/particles.py @@ -41,8 +41,9 @@ def test_calc_positions(self): x = (-1.0, 1.0) y = (2.0, 3.0) z = (-2.0, -1.0) + X = np.array([x, y, z]) N = 1000 - positions = pa.calc_positions(x, y, z, N) + positions = pa.calc_positions(X, N) self.assertEqual(N, len(positions)) From ca4732cc046130d9385867e9796bf934f364a5d1 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Thu, 6 Oct 2022 13:39:21 +0200 Subject: [PATCH 29/54] added particle writer --- dsmc/writer/__init__.py | 1 + dsmc/writer/particles.py | 5 +++++ 2 files changed, 6 insertions(+) create mode 100644 dsmc/writer/particles.py diff --git a/dsmc/writer/__init__.py b/dsmc/writer/__init__.py index 6847e5f..a53d8ad 100644 --- a/dsmc/writer/__init__.py +++ b/dsmc/writer/__init__.py @@ -1,2 +1,3 @@ from .octree import write_buttom_leafs from .planar import write as write_planar +from .particles import write_positions diff --git a/dsmc/writer/particles.py b/dsmc/writer/particles.py new file mode 100644 index 0000000..b6f2c5d --- /dev/null +++ b/dsmc/writer/particles.py @@ -0,0 +1,5 @@ +def write_positions(particles, file_name="positions.csv"): + with open(file_name, "w") as file: + file.write("x, y, z\n") + for pos in particles.Pos: + file.write("{}, {}, {}\n".format(pos[0], pos[1], pos[2])) \ No newline at end of file From dc582f0e4d81f4d99537a473d0fa432029416689 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Thu, 6 Oct 2022 13:39:37 +0200 Subject: [PATCH 30/54] added hyper sonic mini test case --- examples/hypersonic_flow_mini.py | 70 ++++++++++++++++++++++++++++++++ 1 file changed, 70 insertions(+) create mode 100644 examples/hypersonic_flow_mini.py diff --git a/examples/hypersonic_flow_mini.py b/examples/hypersonic_flow_mini.py new file mode 100644 index 0000000..109eb1c --- /dev/null +++ b/examples/hypersonic_flow_mini.py @@ -0,0 +1,70 @@ +import dsmc +import dsmc.writer as wrt +import dsmc.octree as oc +import time +import numpy as np + +if __name__ == '__main__': + # general parameters + solver = dsmc.DSMC() + domain = [(-3.0, 3.0), (-1.5, 1.5), (-0.025, 0.025)] + obj = [(-0.25, 0.25), (-0.25, 0.25), (-0.5, 0.5)] + dt = 1e-6 + w = 0.2 * 2.3e+15 + mass = 6.6422e-26 + T = 273.0 + n = 2.6e+19 + u = np.array([0.0, -3043.0, 0.0]) + niter = 500 + + h = domain[2][1] - domain[2][0] + + # setup solver + solver.set_domain(domain) + solver.set_weight(w) + solver.set_mass(mass) + + # set object + solver.add_object(obj) + + # trees + tree_inner = oc.Octree() + tree_outer = oc.Octree() + + inner_pos = np.array([[-0.25, -0.25, -0.025], [0.25, 0.25, 0.025]]) + outer_pos = np.array([[-3.0, -1.5, -0.025], [3.0, 1.5, 0.025]]) + + tree_inner.build(inner_pos) + tree_outer.build(outer_pos) + + # create particles + positions = np.zeros((2, 3)) + velocities = np.zeros((2, 3)) + + positions[0][0] = -0.01 + positions[0][1] = 1.4 + positions[0][2] = -0.01 + + positions[1][0] = 0.01 + positions[1][1] = 1.3 + positions[1][2] = 0.01 + + velocities[0][1] = -3043.0 + velocities[1][1] = -3043.0 + + solver.particles.VelPos = (velocities, positions) + + # start timing + start_time = time.time() + + for it in range(niter): + print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) + solver.advance(dt) + wrt.write_positions(solver.particles, "pos_{}.csv".format(it)) + + wrt.write_buttom_leafs(tree_inner, "innter.vtu") + wrt.write_buttom_leafs(tree_outer, "outer.vtu") + + print("") + print("--- %s seconds ---" % (time.time() - start_time)) + print('done') \ No newline at end of file From 5ae58e7c7a3402a9c1b6d861ba8fca0969e08ce7 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 09:00:32 +0200 Subject: [PATCH 31/54] added writing of boundary geometry for hypersonic test case --- examples/hypersonic_flow.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/examples/hypersonic_flow.py b/examples/hypersonic_flow.py index 4bf0a79..d084f39 100644 --- a/examples/hypersonic_flow.py +++ b/examples/hypersonic_flow.py @@ -2,6 +2,7 @@ import dsmc.writer as wrt import dsmc.diagnostics as dia import dsmc.mesh.mesh2d as msh2d +import dsmc.octree as oc import time import numpy as np @@ -53,6 +54,19 @@ # create particles solver.create_particles(domain, T, n, u) + # trees + tree_inner = oc.Octree() + tree_outer = oc.Octree() + + inner_pos = np.array([[-0.25, -0.25, -0.025], [0.25, 0.25, 0.025]]) + outer_pos = np.array([[-3.0, -1.5, -0.025], [3.0, 1.5, 0.025]]) + + tree_inner.build(inner_pos) + tree_outer.build(outer_pos) + + wrt.write_buttom_leafs(tree_inner, "innter.vtu") + wrt.write_buttom_leafs(tree_outer, "outer.vtu") + # start timing start_time = time.time() From 56be2684c3102e0426c8c6456d38d5f3eca1aff6 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 09:31:59 +0200 Subject: [PATCH 32/54] added to docs --- doc/Boundary.ipynb | 133 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 133 insertions(+) create mode 100644 doc/Boundary.ipynb diff --git a/doc/Boundary.ipynb b/doc/Boundary.ipynb new file mode 100644 index 0000000..51dbb40 --- /dev/null +++ b/doc/Boundary.ipynb @@ -0,0 +1,133 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "d49edbe1-fe92-4a67-82e9-78a877139d5d", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import sympy as sym" + ] + }, + { + "cell_type": "markdown", + "id": "7f598afc-282c-4bb5-8d47-80c77fcc1178", + "metadata": {}, + "source": [ + "**Boundary collision method**\n", + "\n", + "For a given particle position $\\vec{P}$ and an old particle positions $\\vec{P}_{\\mathrm{old}}$ and for the plane there are 3 points: $\\vec{A}$, $\\vec{B}$ and $\\vec{C}$.\n", + "The equation for the line representing the particle movement can be writen as\n", + "\n", + "$$\n", + "\\vec{l} = \\vec{P}_{\\mathrm{old}} + t \\cdot \\vec{D}\n", + "$$\n", + "\n", + "where $\\vec{D} = \\vec{P}_{\\mathrm{old}} - \\vec{P}$ and the equation representing the plane can be written as\n", + "\n", + "$$\n", + "\\vec{p} = \\vec{A} + k \\cdot \\vec{Q}_1 + s \\cdot \\vec{Q}_2\n", + "$$\n", + "\n", + "with $\\vec{Q}_1 = \\vec{B} - \\vec{A}$ and $\\vec{Q}_2 = \\vec{C} - \\vec{A}$.\n", + "For the line to intersect the plane both need to be equal $\\vec{l} = \\vec{p}$ which leads to the equation\n", + "\n", + "$$\n", + "k \\cdot \\vec{Q}_1 + s \\cdot \\vec{Q}_2 - t \\cdot \\vec{D} = \\vec{R}\n", + "$$\n", + "\n", + "where $\\vec{R} = \\vec{P}_{\\mathrm{old}} - \\vec{A}$.\n", + "The equation can be written in matrix form as\n", + "\n", + "$$\n", + "\\begin{bmatrix}\n", + "\\vec{Q}_{1, x} & \\vec{Q}_{2, x} & -\\vec{D}_x \\\\\n", + "\\vec{Q}_{1, y} & \\vec{Q}_{2, y} & -\\vec{D}_y \\\\ \n", + "\\vec{Q}_{1, z} & \\vec{Q}_{2, z} & -\\vec{D}_z\n", + "\\end{bmatrix}\n", + "\\cdot\n", + "\\begin{pmatrix}\n", + "k \\\\\n", + "s \\\\\n", + "t\n", + "\\end{pmatrix}\n", + "=\n", + "\\begin{pmatrix}\n", + "\\vec{R}_x\\\\\n", + "\\vec{R}_y \\\\\n", + "\\vec{R}_z\n", + "\\end{pmatrix}.\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "4abb0e41-5e15-4326-b663-df0b6b5cfe0b", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0. 0. 0.] [0. 0. 0.]\n" + ] + } + ], + "source": [ + "P = np.array([1.0, 1.0, 0.0])\n", + "Pold = np.array([-1.0, -1.0, 0.0])\n", + "\n", + "A = np.array([0.0, -1.0, -1.0])\n", + "B = np.array([0.0, -1.0, 1.0])\n", + "C = np.array([0.0, 1.0, 1.0])\n", + "\n", + "D = Pold - P\n", + "Q1 = B - A\n", + "Q2 = C - A\n", + "R = Pold - A\n", + "M = np.array([[Q1[i], Q2[i], -D[i]] for i in range(3)])\n", + "\n", + "k, s, t = np.linalg.solve(M, R)\n", + "\n", + "p = A + k*Q1 + s*Q2\n", + "l = Pold + t *D\n", + "\n", + "print(p, l)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d204319a-eb91-47f0-8c5d-8cd05e2f3c2c", + "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.8.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From e8f97e660f508fa2d96a262419f7424e7757d23b Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 10:14:44 +0200 Subject: [PATCH 33/54] updated docs --- doc/Boundary.ipynb | 93 ++++++++++++++++++++++++---------------------- 1 file changed, 49 insertions(+), 44 deletions(-) diff --git a/doc/Boundary.ipynb b/doc/Boundary.ipynb index 51dbb40..58eae27 100644 --- a/doc/Boundary.ipynb +++ b/doc/Boundary.ipynb @@ -18,53 +18,41 @@ "source": [ "**Boundary collision method**\n", "\n", - "For a given particle position $\\vec{P}$ and an old particle positions $\\vec{P}_{\\mathrm{old}}$ and for the plane there are 3 points: $\\vec{A}$, $\\vec{B}$ and $\\vec{C}$.\n", - "The equation for the line representing the particle movement can be writen as\n", + "Given a plane $\\Pi$ and and a line $L$ as\n", "\n", "$$\n", - "\\vec{l} = \\vec{P}_{\\mathrm{old}} + t \\cdot \\vec{D}\n", + "\\Pi \\rightarrow (p - p_0) \\cdot \\vec{n}_p = 0 \\\\\n", + "L \\rightarrow l = l_0 + t \\cdot \\vec{n}_l\n", "$$\n", "\n", - "where $\\vec{D} = \\vec{P}_{\\mathrm{old}} - \\vec{P}$ and the equation representing the plane can be written as\n", + "where $\\vec{n}_1$ is the normal vector of a plane given as\n", "\n", "$$\n", - "\\vec{p} = \\vec{A} + k \\cdot \\vec{Q}_1 + s \\cdot \\vec{Q}_2\n", + "\\vec{n}_p = \\frac{(p_1 - p_0) \\times (p_2 - p_0)}{|| (p_1 - p_0) \\times (p_2 - p_0) ||}\n", "$$\n", "\n", - "with $\\vec{Q}_1 = \\vec{B} - \\vec{A}$ and $\\vec{Q}_2 = \\vec{C} - \\vec{A}$.\n", - "For the line to intersect the plane both need to be equal $\\vec{l} = \\vec{p}$ which leads to the equation\n", + "is a unit vector in the direction of the line as\n", "\n", "$$\n", - "k \\cdot \\vec{Q}_1 + s \\cdot \\vec{Q}_2 - t \\cdot \\vec{D} = \\vec{R}\n", + "\\vec{n}_l = \\frac{l_1 - l_0}{|| l_1 - l_0 ||}.\n", "$$\n", "\n", - "where $\\vec{R} = \\vec{P}_{\\mathrm{old}} - \\vec{A}$.\n", - "The equation can be written in matrix form as\n", + "The positons of the intersection point found by setting the point $p$ in the equation for the plane as being a point in the equation of the line:\n", "\n", "$$\n", - "\\begin{bmatrix}\n", - "\\vec{Q}_{1, x} & \\vec{Q}_{2, x} & -\\vec{D}_x \\\\\n", - "\\vec{Q}_{1, y} & \\vec{Q}_{2, y} & -\\vec{D}_y \\\\ \n", - "\\vec{Q}_{1, z} & \\vec{Q}_{2, z} & -\\vec{D}_z\n", - "\\end{bmatrix}\n", - "\\cdot\n", - "\\begin{pmatrix}\n", - "k \\\\\n", - "s \\\\\n", - "t\n", - "\\end{pmatrix}\n", - "=\n", - "\\begin{pmatrix}\n", - "\\vec{R}_x\\\\\n", - "\\vec{R}_y \\\\\n", - "\\vec{R}_z\n", - "\\end{pmatrix}.\n", + "(l - p_0) \\cdot \\vec{n}_p = 0 = (l_0 + t \\cdot \\vec{n}_l - p_0) \\cdot \\vec{n}_p \\implies t = - \\frac{(l_0 - p_0) \\cdot \\vec{n}_p}{\\vec{n}_p \\cdot \\vec{n}_l}\n", + "$$\n", + "\n", + "inserting the line equation the intersection point can be found as\n", + "\n", + "$$\n", + "l^* = l_0 - \\frac{(l_0 - p_0) \\cdot \\vec{n}_p}{\\vec{n}_p \\cdot \\vec{n}_l} \\cdot \\vec{n}_l.\n", "$$" ] }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 2, "id": "4abb0e41-5e15-4326-b663-df0b6b5cfe0b", "metadata": { "tags": [] @@ -74,36 +62,53 @@ "name": "stdout", "output_type": "stream", "text": [ - "[0. 0. 0.] [0. 0. 0.]\n" + "[0. 0. 0.]\n" ] } ], "source": [ - "P = np.array([1.0, 1.0, 0.0])\n", - "Pold = np.array([-1.0, -1.0, 0.0])\n", + "l1 = np.array([1.0, 1.0, 0.0])\n", + "l0 = np.array([-1.0, -1.0, 0.0])\n", + "\n", + "p0 = np.array([0.0, -1.0, -1.0])\n", + "p1 = np.array([0.0, -1.0, 1.0])\n", + "p2 = np.array([0.0, 1.0, 1.0])\n", "\n", - "A = np.array([0.0, -1.0, -1.0])\n", - "B = np.array([0.0, -1.0, 1.0])\n", - "C = np.array([0.0, 1.0, 1.0])\n", + "nl = l1 - l0\n", + "nl = nl / np.linalg.norm(nl)\n", + "\n", + "n_p = np.cross((p1 - p0), (p2 - p0))\n", + "n_p = n_p / np.linalg.norm(n_p)\n", + "\n", + "ls = l0 - ((l0 - p0).dot(n_p) / n_p.dot(nl))*nl\n", + "\n", + "print(ls)" + ] + }, + { + "cell_type": "markdown", + "id": "d43f6e71-7324-4b21-9c0d-371fb76d01f6", + "metadata": {}, + "source": [ + "**Reflections**\n", "\n", - "D = Pold - P\n", - "Q1 = B - A\n", - "Q2 = C - A\n", - "R = Pold - A\n", - "M = np.array([[Q1[i], Q2[i], -D[i]] for i in range(3)])\n", + "The reflected line can be found as\n", "\n", - "k, s, t = np.linalg.solve(M, R)\n", + "$$\n", + "L_R \\rightarrow l_R = l^* + \\lambda \\cdot \\vec{n}_R\n", + "$$\n", "\n", - "p = A + k*Q1 + s*Q2\n", - "l = Pold + t *D\n", + "with \n", "\n", - "print(p, l)" + "$$\n", + "\\vec{n}_R = \\vec{n}_l - 2 \\frac{\\vec{n}_p \\cdot \\vec{n}_l}{|| \\vec{n}_p ||^2} \\cdot \\vec{n}_p.\n", + "$$" ] }, { "cell_type": "code", "execution_count": null, - "id": "d204319a-eb91-47f0-8c5d-8cd05e2f3c2c", + "id": "91de06b7-7511-471d-9b64-f9c2a784d682", "metadata": {}, "outputs": [], "source": [] From 80041dbc5aba48f6f7a9f50e5561af64dc310388 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 10:42:33 +0200 Subject: [PATCH 34/54] updated docs --- doc/Boundary.ipynb | 47 +++++++++++++++++++++++++++++++++++++--------- 1 file changed, 38 insertions(+), 9 deletions(-) diff --git a/doc/Boundary.ipynb b/doc/Boundary.ipynb index 58eae27..f2bc8b3 100644 --- a/doc/Boundary.ipynb +++ b/doc/Boundary.ipynb @@ -28,13 +28,13 @@ "where $\\vec{n}_1$ is the normal vector of a plane given as\n", "\n", "$$\n", - "\\vec{n}_p = \\frac{(p_1 - p_0) \\times (p_2 - p_0)}{|| (p_1 - p_0) \\times (p_2 - p_0) ||}\n", + "\\vec{n}_p = (p_1 - p_0) \\times (p_2 - p_0)\n", "$$\n", "\n", - "is a unit vector in the direction of the line as\n", + "is a vector in the direction of the line as\n", "\n", "$$\n", - "\\vec{n}_l = \\frac{l_1 - l_0}{|| l_1 - l_0 ||}.\n", + "\\vec{n}_l = l_1 - l_0.\n", "$$\n", "\n", "The positons of the intersection point found by setting the point $p$ in the equation for the plane as being a point in the equation of the line:\n", @@ -46,7 +46,7 @@ "inserting the line equation the intersection point can be found as\n", "\n", "$$\n", - "l^* = l_0 - \\frac{(l_0 - p_0) \\cdot \\vec{n}_p}{\\vec{n}_p \\cdot \\vec{n}_l} \\cdot \\vec{n}_l.\n", + "l^* = l_0 + t \\cdot \\vec{n}_l.\n", "$$" ] }, @@ -75,12 +75,11 @@ "p2 = np.array([0.0, 1.0, 1.0])\n", "\n", "nl = l1 - l0\n", - "nl = nl / np.linalg.norm(nl)\n", - "\n", "n_p = np.cross((p1 - p0), (p2 - p0))\n", - "n_p = n_p / np.linalg.norm(n_p)\n", "\n", - "ls = l0 - ((l0 - p0).dot(n_p) / n_p.dot(nl))*nl\n", + "t = - ((l0 - p0).dot(n_p) / n_p.dot(nl))\n", + "\n", + "ls = l0 +t*nl\n", "\n", "print(ls)" ] @@ -102,14 +101,44 @@ "\n", "$$\n", "\\vec{n}_R = \\vec{n}_l - 2 \\frac{\\vec{n}_p \\cdot \\vec{n}_l}{|| \\vec{n}_p ||^2} \\cdot \\vec{n}_p.\n", + "$$\n", + "\n", + "The reflected point can now the found as \n", + "\n", + "$$\n", + "l_R = l^* + (1 - t) \\cdot \\vec{n}_R\n", "$$" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "id": "91de06b7-7511-471d-9b64-f9c2a784d682", "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 1., -1., 0.])" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "nR = nl - 2*(n_p.dot(nl) / n_p.dot(n_p))*n_p\n", + "lR = ls - (1 -t )*nR\n", + "\n", + "lR" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a81e9e3a-d8d0-463f-a5e4-b8700e696875", + "metadata": {}, "outputs": [], "source": [] } From 353c704bae326961a95dfe9121fe9c98e532ba41 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 11:02:11 +0200 Subject: [PATCH 35/54] added boundary --- dsmc/boundary.py | 21 +++++++++++++++++++++ tests/unit/main.py | 1 + tests/unit/test_dsmc/boundary.py | 14 ++++++++++++++ 3 files changed, 36 insertions(+) create mode 100644 dsmc/boundary.py create mode 100644 tests/unit/test_dsmc/boundary.py diff --git a/dsmc/boundary.py b/dsmc/boundary.py new file mode 100644 index 0000000..a3f543d --- /dev/null +++ b/dsmc/boundary.py @@ -0,0 +1,21 @@ +import numpy as np +from numba import njit + +@njit +def _check_if_parallel(v1, v2, diff=1e-6): + V1 = np.copy(v1) / np.linalg.norm(v1) + V2 = np.copy(v2) / np.linalg.norm(v2) + + return V1.dot(V2) > diff + +def _intersect(l0, l1, p0, p1, p2): + """ + Args: + l0 : first position on line + l1 : second position on line + p0 : first position on plane + p1 : first position on plane + p2 : first position on plane + """ + n_l = l1 - l0 + n_p = (p1 - p0).cross(p2 - p1) \ No newline at end of file diff --git a/tests/unit/main.py b/tests/unit/main.py index 3a94703..4cae687 100644 --- a/tests/unit/main.py +++ b/tests/unit/main.py @@ -8,6 +8,7 @@ from test_dsmc.diagnostics import * from test_dsmc.common import * from test_dsmc.mesh.mesh2d import * +from test_dsmc.boundary import * if __name__ == '__main__': unittest.main() diff --git a/tests/unit/test_dsmc/boundary.py b/tests/unit/test_dsmc/boundary.py new file mode 100644 index 0000000..15e7000 --- /dev/null +++ b/tests/unit/test_dsmc/boundary.py @@ -0,0 +1,14 @@ +import numpy as np +import dsmc.boundary as bo +import unittest + +class TestCommon(unittest.TestCase): + def test__check_if_parallel(self): + v1 = np.array([1.0, 0.0, 0.0]) + v2 = np.array([1.0, 0.0, 0.0]) + v3 = np.array([0.0, 1.0, 1.0]) + v4 = np.array([1.0e-6, 10.0, 10.0]) + + self.assertTrue(bo._check_if_parallel(v1, v2)) + self.assertFalse(bo._check_if_parallel(v1, v3)) + self.assertFalse(bo._check_if_parallel(v1, v4)) \ No newline at end of file From 470a0767141637a456812be6dbe7bfd4d7052c42 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 11:10:32 +0200 Subject: [PATCH 36/54] implemented _intersect --- dsmc/boundary.py | 8 +++++++- tests/unit/test_dsmc/boundary.py | 24 +++++++++++++++++++++++- 2 files changed, 30 insertions(+), 2 deletions(-) diff --git a/dsmc/boundary.py b/dsmc/boundary.py index a3f543d..ace9733 100644 --- a/dsmc/boundary.py +++ b/dsmc/boundary.py @@ -8,6 +8,7 @@ def _check_if_parallel(v1, v2, diff=1e-6): return V1.dot(V2) > diff +@njit def _intersect(l0, l1, p0, p1, p2): """ Args: @@ -18,4 +19,9 @@ def _intersect(l0, l1, p0, p1, p2): p2 : first position on plane """ n_l = l1 - l0 - n_p = (p1 - p0).cross(p2 - p1) \ No newline at end of file + n_p = np.cross((p1 - p0), (p2 - p1)) + + if _check_if_parallel(n_l, n_p): + return (False, n_l, n_p, 0.0) + else: + return (True, n_l, n_p, - ((l0 - p0).dot(n_p) / n_p.dot(n_l))) \ No newline at end of file diff --git a/tests/unit/test_dsmc/boundary.py b/tests/unit/test_dsmc/boundary.py index 15e7000..fe38a58 100644 --- a/tests/unit/test_dsmc/boundary.py +++ b/tests/unit/test_dsmc/boundary.py @@ -11,4 +11,26 @@ def test__check_if_parallel(self): self.assertTrue(bo._check_if_parallel(v1, v2)) self.assertFalse(bo._check_if_parallel(v1, v3)) - self.assertFalse(bo._check_if_parallel(v1, v4)) \ No newline at end of file + self.assertFalse(bo._check_if_parallel(v1, v4)) + + def test__intersect(self): + l0 = np.array([-1.0, -1.0, 0.0]) + l1 = np.array([1.0, 1.0, 0.0]) + + p0 = np.array([0.0, -1.0, -1.0]) + p1 = np.array([0.0, -1.0, 1.0]) + p2 = np.array([0.0, 1.0, 1.0]) + + intersected, n_l, n_r, t = bo._intersect(l0, l1, p0, p1, p2) + + self.assertTrue(intersected) + + self.assertEqual(2.0, n_l[0]) + self.assertEqual(2.0, n_l[1]) + self.assertEqual(0.0, n_l[2]) + + self.assertEqual(-4.0, n_r[0]) + self.assertEqual(0.0, n_r[1]) + self.assertEqual(0.0, n_r[2]) + + self.assertEqual(0.5, t) \ No newline at end of file From a8379e4aed8334d8f05291886830f98644559adf Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 11:30:39 +0200 Subject: [PATCH 37/54] implemented _calc_nr --- dsmc/boundary.py | 21 ++++++++++++++++++++- tests/unit/test_dsmc/boundary.py | 25 ++++++++++++++++++++----- 2 files changed, 40 insertions(+), 6 deletions(-) diff --git a/dsmc/boundary.py b/dsmc/boundary.py index ace9733..1020b4b 100644 --- a/dsmc/boundary.py +++ b/dsmc/boundary.py @@ -17,6 +17,9 @@ def _intersect(l0, l1, p0, p1, p2): p0 : first position on plane p1 : first position on plane p2 : first position on plane + + Returns: + (intersected, n_l, n_r, t) """ n_l = l1 - l0 n_p = np.cross((p1 - p0), (p2 - p1)) @@ -24,4 +27,20 @@ def _intersect(l0, l1, p0, p1, p2): if _check_if_parallel(n_l, n_p): return (False, n_l, n_p, 0.0) else: - return (True, n_l, n_p, - ((l0 - p0).dot(n_p) / n_p.dot(n_l))) \ No newline at end of file + return (True, n_l, n_p, - ((l0 - p0).dot(n_p) / n_p.dot(n_l))) + +@njit +def _calc_nr(n_l, n_p): + return n_l - 2.0 * (n_p.dot(n_l) / n_p.dot(n_p))*n_p + +def _reflect(vel, pos, pos_old, p0, p1, p2): + intersected, n_l, n_p, t = _intersect(pos_old, pos, p0, p1, p2) + + if intersected: + pos_old = pos_old + n_l*t + n_r = _calc_nr(n_l, n_p) + pos = pos_old + (1.0 - t)*n_r + vel = pos_old + (np.linalg.norm(vel) / np.linalg.norm(n_r))*n_r + + return (vel, pos, pos_old) + \ No newline at end of file diff --git a/tests/unit/test_dsmc/boundary.py b/tests/unit/test_dsmc/boundary.py index fe38a58..99f0cbb 100644 --- a/tests/unit/test_dsmc/boundary.py +++ b/tests/unit/test_dsmc/boundary.py @@ -21,7 +21,7 @@ def test__intersect(self): p1 = np.array([0.0, -1.0, 1.0]) p2 = np.array([0.0, 1.0, 1.0]) - intersected, n_l, n_r, t = bo._intersect(l0, l1, p0, p1, p2) + intersected, n_l, n_p, t = bo._intersect(l0, l1, p0, p1, p2) self.assertTrue(intersected) @@ -29,8 +29,23 @@ def test__intersect(self): self.assertEqual(2.0, n_l[1]) self.assertEqual(0.0, n_l[2]) - self.assertEqual(-4.0, n_r[0]) - self.assertEqual(0.0, n_r[1]) - self.assertEqual(0.0, n_r[2]) + self.assertEqual(-4.0, n_p[0]) + self.assertEqual(0.0, n_p[1]) + self.assertEqual(0.0, n_p[2]) - self.assertEqual(0.5, t) \ No newline at end of file + self.assertEqual(0.5, t) + + def test__calc_nr(self): + l0 = np.array([-1.0, -1.0, 0.0]) + l1 = np.array([1.0, 1.0, 0.0]) + + p0 = np.array([0.0, -1.0, -1.0]) + p1 = np.array([0.0, -1.0, 1.0]) + p2 = np.array([0.0, 1.0, 1.0]) + + intersected, n_l, n_p, t = bo._intersect(l0, l1, p0, p1, p2) + n_r = bo._calc_nr(n_l, n_p) + + self.assertEqual(-2.0, n_r[0]) + self.assertEqual(2.0, n_r[1]) + self.assertEqual(0.0, n_r[2]) \ No newline at end of file From 3de87ce853148ba0f8ae2d07b0126b505d885750 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 11:34:24 +0200 Subject: [PATCH 38/54] implemented reflect --- dsmc/boundary.py | 3 ++- tests/unit/test_dsmc/boundary.py | 25 ++++++++++++++++++++++++- 2 files changed, 26 insertions(+), 2 deletions(-) diff --git a/dsmc/boundary.py b/dsmc/boundary.py index 1020b4b..b75ce44 100644 --- a/dsmc/boundary.py +++ b/dsmc/boundary.py @@ -32,7 +32,8 @@ def _intersect(l0, l1, p0, p1, p2): @njit def _calc_nr(n_l, n_p): return n_l - 2.0 * (n_p.dot(n_l) / n_p.dot(n_p))*n_p - + +@njit def _reflect(vel, pos, pos_old, p0, p1, p2): intersected, n_l, n_p, t = _intersect(pos_old, pos, p0, p1, p2) diff --git a/tests/unit/test_dsmc/boundary.py b/tests/unit/test_dsmc/boundary.py index 99f0cbb..de1290f 100644 --- a/tests/unit/test_dsmc/boundary.py +++ b/tests/unit/test_dsmc/boundary.py @@ -48,4 +48,27 @@ def test__calc_nr(self): self.assertEqual(-2.0, n_r[0]) self.assertEqual(2.0, n_r[1]) - self.assertEqual(0.0, n_r[2]) \ No newline at end of file + self.assertEqual(0.0, n_r[2]) + + def test__reflect(self): + vel = np.array([123.0, 123.0, 0.0]) + pos_old = np.array([-1.0, -1.0, 0.0]) + pos = np.array([1.0, 1.0, 0.0]) + + p0 = np.array([0.0, -1.0, -1.0]) + p1 = np.array([0.0, -1.0, 1.0]) + p2 = np.array([0.0, 1.0, 1.0]) + + new_vel, new_pos, new_pos_old = bo._reflect(vel, pos, pos_old, p0, p1, p2) + + self.assertEqual(-123.0, new_vel[0]) + self.assertEqual(123.0, new_vel[1]) + self.assertEqual(0.0, new_vel[2]) + + self.assertEqual(-1.0, new_pos[0]) + self.assertEqual(1.0, new_pos[1]) + self.assertEqual(0.0, new_pos[2]) + + self.assertEqual(0.0, new_pos_old[0]) + self.assertEqual(0.0, new_pos_old[1]) + self.assertEqual(0.0, new_pos_old[2]) \ No newline at end of file From ac900830b853ea3abba67d7744f34423d4c9eacd Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 12:08:33 +0200 Subject: [PATCH 39/54] implemented _get_plane --- dsmc/boundary.py | 80 +++++++++++++++++++- tests/unit/test_dsmc/boundary.py | 123 ++++++++++++++++++++++++++++++- 2 files changed, 201 insertions(+), 2 deletions(-) diff --git a/dsmc/boundary.py b/dsmc/boundary.py index b75ce44..4c68b63 100644 --- a/dsmc/boundary.py +++ b/dsmc/boundary.py @@ -44,4 +44,82 @@ def _reflect(vel, pos, pos_old, p0, p1, p2): vel = pos_old + (np.linalg.norm(vel) / np.linalg.norm(n_r))*n_r return (vel, pos, pos_old) - \ No newline at end of file + +@njit +def _get_plane(domain, i, j): + if i == 0: + p0 = np.array([domain[i][j], domain[1][0], domain[2][0]]) + p1 = np.array([domain[i][j], domain[1][1], domain[2][0]]) + p2 = np.array([domain[i][j], domain[1][0], domain[2][1]]) + elif i == 1: + p0 = np.array([domain[0][0], domain[i][j], domain[2][0]]) + p1 = np.array([domain[0][1], domain[i][j], domain[2][0]]) + p2 = np.array([domain[0][0], domain[i][j], domain[2][1]]) + elif i == 2: + p0 = np.array([domain[0][0], domain[1][0], domain[i][j]]) + p1 = np.array([domain[0][1], domain[1][0], domain[i][j]]) + p2 = np.array([domain[0][0], domain[1][1], domain[i][j]]) + + return (p0, p1, p2) + +def _boundary(velocities, positions, old_positions, domain, boundary_conds): + kept_parts = np.ones(positions.shape[0], dtype=np.uint) + + for p in range(len(positions)): + for i in range(3): + for j in range(2): + p0, p1, p2 = _get_plane(domain, i, j) + if boundary_conds[i][j] == 0: + velocities[p], positions[p], old_positions[p] = _reflect(velocities[p], positions[p], old_positions[p], p0, p1, p2) + elif boundary_conds[i][j] == 1 or boundary_conds[i][j] == 2: + if _intersect(old_positions[p], positions[p], p0, p1, p2)[0]: + kept_parts[p] = 0 + + N = int(sum(kept_parts)) + p = 0 + new_velocities = np.empty((N, 3)) + new_positions = np.empty((N, 3)) + new_old_positions = np.empty((N, 3)) + + for i in range(positions.shape[0]): + if kept_parts[i] == 1: + new_velocities[p] = velocities[i] + new_positions[p] = positions[i] + new_old_positions[p] = old_positions[p] + p += 1 + else: + continue + + return (new_velocities, new_positions, new_old_positions) + +@njit +def _get_boundary(boundary): + if boundary == "xmin": + return (0, 0) + elif boundary == "xmax": + return (0, 1) + elif boundary == "ymin": + return (1, 0) + elif boundary == "ymax": + return (1, 1) + elif boundary == "zmin": + return (2, 0) + elif boundary == "zmax": + return (2, 1) + +@njit +def _get_bc_type(bc_type): + if bc_type == "ela": + return 0 + elif bc_type == "open": + return 1 + elif bc_type == "inflow": + return 2 + +class Boundary: + def __init__(self): + self.T = np.ones((3, 2))*300.0 + self.n = np.ones((3, 2))*1e+18 + self.u = np.zeros((3, 2, 3)) + self.boundary_conds = np.array([[0, 0], [0, 0], [0, 0]], dtype=np.uint) # 0 = ela, 1 = open, 2 = inflow + self.domain = None \ No newline at end of file diff --git a/tests/unit/test_dsmc/boundary.py b/tests/unit/test_dsmc/boundary.py index de1290f..2fd4b13 100644 --- a/tests/unit/test_dsmc/boundary.py +++ b/tests/unit/test_dsmc/boundary.py @@ -71,4 +71,125 @@ def test__reflect(self): self.assertEqual(0.0, new_pos_old[0]) self.assertEqual(0.0, new_pos_old[1]) - self.assertEqual(0.0, new_pos_old[2]) \ No newline at end of file + self.assertEqual(0.0, new_pos_old[2]) + + def test__get_plane1(self): + domain = np.array([(0, 1), (2, 4), (4, 7)]) + axis = 0 + + p0, p1, p2 = bo._get_plane(domain, axis, 0) + + self.assertEqual((3,), p0.shape) + self.assertEqual((3,), p1.shape) + self.assertEqual((3,), p2.shape) + + self.assertEqual(0.0, p0[0]) + self.assertEqual(2.0, p0[1]) + self.assertEqual(4.0, p0[2]) + + self.assertEqual(0.0, p1[0]) + self.assertEqual(4.0, p1[1]) + self.assertEqual(4.0, p1[2]) + + self.assertEqual(0.0, p2[0]) + self.assertEqual(2.0, p2[1]) + self.assertEqual(7.0, p2[2]) + + p0, p1, p2 = bo._get_plane(domain, axis, 1) + + self.assertEqual((3,), p0.shape) + self.assertEqual((3,), p1.shape) + self.assertEqual((3,), p2.shape) + + self.assertEqual(1.0, p0[0]) + self.assertEqual(2.0, p0[1]) + self.assertEqual(4.0, p0[2]) + + self.assertEqual(1.0, p1[0]) + self.assertEqual(4.0, p1[1]) + self.assertEqual(4.0, p1[2]) + + self.assertEqual(1.0, p2[0]) + self.assertEqual(2.0, p2[1]) + self.assertEqual(7.0, p2[2]) + + def test__get_plane2(self): + domain = np.array([(0, 1), (2, 4), (4, 7)]) + axis = 1 + + p0, p1, p2 = bo._get_plane(domain, axis, 0) + + self.assertEqual((3,), p0.shape) + self.assertEqual((3,), p1.shape) + self.assertEqual((3,), p2.shape) + + self.assertEqual(0.0, p0[0]) + self.assertEqual(2.0, p0[1]) + self.assertEqual(4.0, p0[2]) + + self.assertEqual(1.0, p1[0]) + self.assertEqual(2.0, p1[1]) + self.assertEqual(4.0, p1[2]) + + self.assertEqual(0.0, p2[0]) + self.assertEqual(2.0, p2[1]) + self.assertEqual(7.0, p2[2]) + + p0, p1, p2 = bo._get_plane(domain, axis, 1) + + self.assertEqual((3,), p0.shape) + self.assertEqual((3,), p1.shape) + self.assertEqual((3,), p2.shape) + + self.assertEqual(0.0, p0[0]) + self.assertEqual(4.0, p0[1]) + self.assertEqual(4.0, p0[2]) + + self.assertEqual(1.0, p1[0]) + self.assertEqual(4.0, p1[1]) + self.assertEqual(4.0, p1[2]) + + self.assertEqual(0.0, p2[0]) + self.assertEqual(4.0, p2[1]) + self.assertEqual(7.0, p2[2]) + + def test__get_plane3(self): + domain = np.array([(0, 1), (2, 4), (4, 7)]) + axis = 2 + + p0, p1, p2 = bo._get_plane(domain, axis, 0) + + self.assertEqual((3,), p0.shape) + self.assertEqual((3,), p1.shape) + self.assertEqual((3,), p2.shape) + + self.assertEqual(0.0, p0[0]) + self.assertEqual(2.0, p0[1]) + self.assertEqual(4.0, p0[2]) + + self.assertEqual(1.0, p1[0]) + self.assertEqual(2.0, p1[1]) + self.assertEqual(4.0, p1[2]) + + self.assertEqual(0.0, p2[0]) + self.assertEqual(4.0, p2[1]) + self.assertEqual(4.0, p2[2]) + + p0, p1, p2 = bo._get_plane(domain, axis, 1) + + self.assertEqual((3,), p0.shape) + self.assertEqual((3,), p1.shape) + self.assertEqual((3,), p2.shape) + + self.assertEqual(0.0, p0[0]) + self.assertEqual(2.0, p0[1]) + self.assertEqual(7.0, p0[2]) + + self.assertEqual(1.0, p1[0]) + self.assertEqual(2.0, p1[1]) + self.assertEqual(7.0, p1[2]) + + self.assertEqual(0.0, p2[0]) + self.assertEqual(4.0, p2[1]) + self.assertEqual(7.0, p2[2]) + \ No newline at end of file From 9a89d114f60454b8d6e5267efae8e4c96d67b98a Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 12:12:01 +0200 Subject: [PATCH 40/54] updated boundary --- dsmc/boundary.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/dsmc/boundary.py b/dsmc/boundary.py index 4c68b63..96df51c 100644 --- a/dsmc/boundary.py +++ b/dsmc/boundary.py @@ -62,6 +62,7 @@ def _get_plane(domain, i, j): return (p0, p1, p2) +@njit def _boundary(velocities, positions, old_positions, domain, boundary_conds): kept_parts = np.ones(positions.shape[0], dtype=np.uint) @@ -122,4 +123,7 @@ def __init__(self): self.n = np.ones((3, 2))*1e+18 self.u = np.zeros((3, 2, 3)) self.boundary_conds = np.array([[0, 0], [0, 0], [0, 0]], dtype=np.uint) # 0 = ela, 1 = open, 2 = inflow - self.domain = None \ No newline at end of file + self.domain = None + + def boundary(self, velocities, positions, old_positions): + return self._boundary(velocities, positions, old_positions, self.domain, self.boundary_conds) \ No newline at end of file From 2f429c730f6b3efd56b1f5a4c730818ddb4e9f5e Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 13:33:16 +0200 Subject: [PATCH 41/54] added boundary test --- tests/domain/boundary.py | 51 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) create mode 100644 tests/domain/boundary.py diff --git a/tests/domain/boundary.py b/tests/domain/boundary.py new file mode 100644 index 0000000..aaff66a --- /dev/null +++ b/tests/domain/boundary.py @@ -0,0 +1,51 @@ +import dsmc.dsmc as ds +import dsmc.particles as prt +import dsmc.octree as oc +import dsmc.writer as wrt +import dsmc.boundary as bo +import dsmc.common as co +import time +import numpy as np + +if __name__ == '__main__': + # general parameters + boundary = bo.Boundary() + particles = prt.Particles() + domain = np.array([(-1.0, 1.0), (-1.0, 1.0), (-1.0, 1.0)]) + dt = 0.0001 + w = 2.3e+18 + mass = 6.6422e-26 + T = 273.0 + n = 2.6e+19 + N = int(co.get_V(domain)*n/w) + niter = 1000 + + # trees + tree_outer = oc.Octree() + outer_pos = np.array([[-1.0, -1.0, -1.0], [1.0, 1.0, 1.0]]) + tree_outer.build(outer_pos) + wrt.write_buttom_leafs(tree_outer, "outer_box.vtu") + + # setup solver + boundary.domain = domain + + # create particles + particles.create_particles(domain, mass, T, N) + + # start timing + start_time = time.time() + + wrt.write_positions(particles, "pos_{}.csv".format(0)) + + for it in range(niter): + print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) + velocities, positions, old_positions = ds._push(particles.Vel, particles.Pos, dt) + velocities, positions, old_positions = boundary.boundary(velocities, positions, old_positions) + + particles.VelPos = (velocities, positions) + wrt.write_positions(particles, "pos_{}.csv".format(it + 1)) + + + print("") + print("--- %s seconds ---" % (time.time() - start_time)) + print('done') \ No newline at end of file From b4bee5d9c720a2346c774997632ee68275f1d7cf Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 13:35:06 +0200 Subject: [PATCH 42/54] fixed new boundary --- dsmc/boundary.py | 37 ++++++++++++++++++++++++++----------- 1 file changed, 26 insertions(+), 11 deletions(-) diff --git a/dsmc/boundary.py b/dsmc/boundary.py index 96df51c..cfaa701 100644 --- a/dsmc/boundary.py +++ b/dsmc/boundary.py @@ -37,7 +37,7 @@ def _calc_nr(n_l, n_p): def _reflect(vel, pos, pos_old, p0, p1, p2): intersected, n_l, n_p, t = _intersect(pos_old, pos, p0, p1, p2) - if intersected: + if intersected and t < 1: pos_old = pos_old + n_l*t n_r = _calc_nr(n_l, n_p) pos = pos_old + (1.0 - t)*n_r @@ -48,17 +48,32 @@ def _reflect(vel, pos, pos_old, p0, p1, p2): @njit def _get_plane(domain, i, j): if i == 0: - p0 = np.array([domain[i][j], domain[1][0], domain[2][0]]) - p1 = np.array([domain[i][j], domain[1][1], domain[2][0]]) - p2 = np.array([domain[i][j], domain[1][0], domain[2][1]]) + if j == 0: + p0 = np.array([domain[i][j], domain[1][0], domain[2][0]]) + p1 = np.array([domain[i][j], domain[1][1], domain[2][0]]) + p2 = np.array([domain[i][j], domain[1][0], domain[2][1]]) + elif j == 1: + p0 = np.array([domain[i][j], domain[1][0], domain[2][0]]) + p1 = np.array([domain[i][j], domain[1][0], domain[2][1]]) + p2 = np.array([domain[i][j], domain[1][1], domain[2][0]]) elif i == 1: - p0 = np.array([domain[0][0], domain[i][j], domain[2][0]]) - p1 = np.array([domain[0][1], domain[i][j], domain[2][0]]) - p2 = np.array([domain[0][0], domain[i][j], domain[2][1]]) + if j == 0: + p0 = np.array([domain[0][0], domain[i][j], domain[2][0]]) + p1 = np.array([domain[0][0], domain[i][j], domain[2][1]]) + p2 = np.array([domain[0][1], domain[i][j], domain[2][0]]) + if j == 1: + p0 = np.array([domain[0][0], domain[i][j], domain[2][0]]) + p1 = np.array([domain[0][1], domain[i][j], domain[2][0]]) + p2 = np.array([domain[0][0], domain[i][j], domain[2][1]]) elif i == 2: - p0 = np.array([domain[0][0], domain[1][0], domain[i][j]]) - p1 = np.array([domain[0][1], domain[1][0], domain[i][j]]) - p2 = np.array([domain[0][0], domain[1][1], domain[i][j]]) + if j == 0: + p0 = np.array([domain[0][0], domain[1][0], domain[i][j]]) + p1 = np.array([domain[0][1], domain[1][0], domain[i][j]]) + p2 = np.array([domain[0][0], domain[1][1], domain[i][j]]) + if j == 1: + p0 = np.array([domain[0][0], domain[1][0], domain[i][j]]) + p1 = np.array([domain[0][0], domain[1][1], domain[i][j]]) + p2 = np.array([domain[0][1], domain[1][0], domain[i][j]]) return (p0, p1, p2) @@ -126,4 +141,4 @@ def __init__(self): self.domain = None def boundary(self, velocities, positions, old_positions): - return self._boundary(velocities, positions, old_positions, self.domain, self.boundary_conds) \ No newline at end of file + return _boundary(velocities, positions, old_positions, self.domain, self.boundary_conds) \ No newline at end of file From 5f896a6a4b6ea161167776c4dac0461c6d240e25 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 14:28:35 +0200 Subject: [PATCH 43/54] updated unit test --- tests/unit/test_dsmc/boundary.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/tests/unit/test_dsmc/boundary.py b/tests/unit/test_dsmc/boundary.py index 2fd4b13..a294398 100644 --- a/tests/unit/test_dsmc/boundary.py +++ b/tests/unit/test_dsmc/boundary.py @@ -106,12 +106,12 @@ def test__get_plane1(self): self.assertEqual(4.0, p0[2]) self.assertEqual(1.0, p1[0]) - self.assertEqual(4.0, p1[1]) - self.assertEqual(4.0, p1[2]) + self.assertEqual(2.0, p1[1]) + self.assertEqual(7.0, p1[2]) self.assertEqual(1.0, p2[0]) - self.assertEqual(2.0, p2[1]) - self.assertEqual(7.0, p2[2]) + self.assertEqual(4.0, p2[1]) + self.assertEqual(4.0, p2[2]) def test__get_plane2(self): domain = np.array([(0, 1), (2, 4), (4, 7)]) @@ -127,13 +127,13 @@ def test__get_plane2(self): self.assertEqual(2.0, p0[1]) self.assertEqual(4.0, p0[2]) - self.assertEqual(1.0, p1[0]) + self.assertEqual(0.0, p1[0]) self.assertEqual(2.0, p1[1]) - self.assertEqual(4.0, p1[2]) + self.assertEqual(7.0, p1[2]) - self.assertEqual(0.0, p2[0]) + self.assertEqual(1.0, p2[0]) self.assertEqual(2.0, p2[1]) - self.assertEqual(7.0, p2[2]) + self.assertEqual(4.0, p2[2]) p0, p1, p2 = bo._get_plane(domain, axis, 1) @@ -185,11 +185,11 @@ def test__get_plane3(self): self.assertEqual(2.0, p0[1]) self.assertEqual(7.0, p0[2]) - self.assertEqual(1.0, p1[0]) - self.assertEqual(2.0, p1[1]) + self.assertEqual(0.0, p1[0]) + self.assertEqual(4.0, p1[1]) self.assertEqual(7.0, p1[2]) - self.assertEqual(0.0, p2[0]) - self.assertEqual(4.0, p2[1]) + self.assertEqual(1.0, p2[0]) + self.assertEqual(2.0, p2[1]) self.assertEqual(7.0, p2[2]) \ No newline at end of file From b91b1d895590d6c730c94c820ff3ab3e4cd75806 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 15:24:34 +0200 Subject: [PATCH 44/54] updated boundary --- dsmc/boundary.py | 40 +++++++++++++++++++++++----------------- tests/domain/boundary.py | 8 ++++++++ 2 files changed, 31 insertions(+), 17 deletions(-) diff --git a/dsmc/boundary.py b/dsmc/boundary.py index cfaa701..f9f988d 100644 --- a/dsmc/boundary.py +++ b/dsmc/boundary.py @@ -3,8 +3,14 @@ @njit def _check_if_parallel(v1, v2, diff=1e-6): - V1 = np.copy(v1) / np.linalg.norm(v1) - V2 = np.copy(v2) / np.linalg.norm(v2) + n1 = np.linalg.norm(v1) + n2 = np.linalg.norm(v2) + + if n1 < 1.0e-13 or n2 < 1.0e-13: + return False + + V1 = np.copy(v1) / n1 + V2 = np.copy(v2) / n2 return V1.dot(V2) > diff @@ -25,9 +31,9 @@ def _intersect(l0, l1, p0, p1, p2): n_p = np.cross((p1 - p0), (p2 - p1)) if _check_if_parallel(n_l, n_p): - return (False, n_l, n_p, 0.0) - else: return (True, n_l, n_p, - ((l0 - p0).dot(n_p) / n_p.dot(n_l))) + else: + return (False, n_l, n_p, 0.0) @njit def _calc_nr(n_l, n_p): @@ -37,7 +43,7 @@ def _calc_nr(n_l, n_p): def _reflect(vel, pos, pos_old, p0, p1, p2): intersected, n_l, n_p, t = _intersect(pos_old, pos, p0, p1, p2) - if intersected and t < 1: + if intersected and t < 1 and t > 0: pos_old = pos_old + n_l*t n_r = _calc_nr(n_l, n_p) pos = pos_old + (1.0 - t)*n_r @@ -49,31 +55,31 @@ def _reflect(vel, pos, pos_old, p0, p1, p2): def _get_plane(domain, i, j): if i == 0: if j == 0: - p0 = np.array([domain[i][j], domain[1][0], domain[2][0]]) - p1 = np.array([domain[i][j], domain[1][1], domain[2][0]]) - p2 = np.array([domain[i][j], domain[1][0], domain[2][1]]) - elif j == 1: p0 = np.array([domain[i][j], domain[1][0], domain[2][0]]) p1 = np.array([domain[i][j], domain[1][0], domain[2][1]]) p2 = np.array([domain[i][j], domain[1][1], domain[2][0]]) + elif j == 1: + p0 = np.array([domain[i][j], domain[1][0], domain[2][0]]) + p1 = np.array([domain[i][j], domain[1][1], domain[2][0]]) + p2 = np.array([domain[i][j], domain[1][0], domain[2][1]]) elif i == 1: if j == 0: - p0 = np.array([domain[0][0], domain[i][j], domain[2][0]]) - p1 = np.array([domain[0][0], domain[i][j], domain[2][1]]) - p2 = np.array([domain[0][1], domain[i][j], domain[2][0]]) - if j == 1: p0 = np.array([domain[0][0], domain[i][j], domain[2][0]]) p1 = np.array([domain[0][1], domain[i][j], domain[2][0]]) p2 = np.array([domain[0][0], domain[i][j], domain[2][1]]) + if j == 1: + p0 = np.array([domain[0][0], domain[i][j], domain[2][0]]) + p1 = np.array([domain[0][0], domain[i][j], domain[2][1]]) + p2 = np.array([domain[0][1], domain[i][j], domain[2][0]]) elif i == 2: if j == 0: - p0 = np.array([domain[0][0], domain[1][0], domain[i][j]]) - p1 = np.array([domain[0][1], domain[1][0], domain[i][j]]) - p2 = np.array([domain[0][0], domain[1][1], domain[i][j]]) - if j == 1: p0 = np.array([domain[0][0], domain[1][0], domain[i][j]]) p1 = np.array([domain[0][0], domain[1][1], domain[i][j]]) p2 = np.array([domain[0][1], domain[1][0], domain[i][j]]) + if j == 1: + p0 = np.array([domain[0][0], domain[1][0], domain[i][j]]) + p1 = np.array([domain[0][1], domain[1][0], domain[i][j]]) + p2 = np.array([domain[0][0], domain[1][1], domain[i][j]]) return (p0, p1, p2) diff --git a/tests/domain/boundary.py b/tests/domain/boundary.py index aaff66a..fd22aef 100644 --- a/tests/domain/boundary.py +++ b/tests/domain/boundary.py @@ -32,6 +32,14 @@ # create particles particles.create_particles(domain, mass, T, N) + velocities, positions = particles.VelPos + + """for vel in velocities: + vel[0] = 0 + vel[1] = 0""" + + particles.VelPos = velocities, positions + # start timing start_time = time.time() From 6f85d68fba0b35b58007c8afef416ba2aa0de2fc Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 16:09:26 +0200 Subject: [PATCH 45/54] moved is_inside to common --- dsmc/common.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/dsmc/common.py b/dsmc/common.py index 17bc81d..24e3fba 100644 --- a/dsmc/common.py +++ b/dsmc/common.py @@ -1,4 +1,5 @@ from numba import cfunc, njit +import numpy.typing as npt @cfunc("double(double[::3, ::2])") def get_V(a): @@ -8,4 +9,12 @@ def get_V(a): def swap(arr, pos1, pos2): temp = arr[pos2] arr[pos2] = arr[pos1] - arr[pos1] = temp \ No newline at end of file + arr[pos1] = temp + +@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 \ No newline at end of file From 298d86c4fad09d1d1a366a9a6d4463bf32e3a527 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 16:11:55 +0200 Subject: [PATCH 46/54] fixed new boundary --- tests/domain/boundary.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/tests/domain/boundary.py b/tests/domain/boundary.py index fd22aef..ae88535 100644 --- a/tests/domain/boundary.py +++ b/tests/domain/boundary.py @@ -13,7 +13,7 @@ particles = prt.Particles() domain = np.array([(-1.0, 1.0), (-1.0, 1.0), (-1.0, 1.0)]) dt = 0.0001 - w = 2.3e+18 + w = 2.3e+16 mass = 6.6422e-26 T = 273.0 n = 2.6e+19 @@ -34,9 +34,8 @@ velocities, positions = particles.VelPos - """for vel in velocities: - vel[0] = 0 - vel[1] = 0""" + #velocities = np.array([[150.0, 150.0, 0.0]]) + #positions = np.array([[0.0, 0.0, 0.0]]) particles.VelPos = velocities, positions @@ -46,7 +45,7 @@ wrt.write_positions(particles, "pos_{}.csv".format(0)) for it in range(niter): - print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) + print("iteration {:4}/{}, N particles {}/{}".format(it + 1, niter, particles.N, N), end="\r", flush=True) velocities, positions, old_positions = ds._push(particles.Vel, particles.Pos, dt) velocities, positions, old_positions = boundary.boundary(velocities, positions, old_positions) From e9acab4b87980e48fdf500c12aca22e91d2a9ae2 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 16:12:03 +0200 Subject: [PATCH 47/54] fixed new boundary --- dsmc/boundary.py | 35 ++++++++++++++++++++++++----------- 1 file changed, 24 insertions(+), 11 deletions(-) diff --git a/dsmc/boundary.py b/dsmc/boundary.py index f9f988d..ff4cd87 100644 --- a/dsmc/boundary.py +++ b/dsmc/boundary.py @@ -1,5 +1,6 @@ import numpy as np from numba import njit +from . import common as co @njit def _check_if_parallel(v1, v2, diff=1e-6): @@ -40,11 +41,16 @@ def _calc_nr(n_l, n_p): return n_l - 2.0 * (n_p.dot(n_l) / n_p.dot(n_p))*n_p @njit -def _reflect(vel, pos, pos_old, p0, p1, p2): +def _reflect(vel, pos, pos_old, p0, p1, p2, domain): intersected, n_l, n_p, t = _intersect(pos_old, pos, p0, p1, p2) - if intersected and t < 1 and t > 0: - pos_old = pos_old + n_l*t + if intersected and t < 1.0 and t > 0.0: + k = 1.0 + p = pos_old + n_l*(t*k) + while not co.is_inside(p, domain): + k *= 0.9 + p = pos_old + n_l*(t*k) + pos_old = p n_r = _calc_nr(n_l, n_p) pos = pos_old + (1.0 - t)*n_r vel = pos_old + (np.linalg.norm(vel) / np.linalg.norm(n_r))*n_r @@ -88,14 +94,21 @@ def _boundary(velocities, positions, old_positions, domain, boundary_conds): kept_parts = np.ones(positions.shape[0], dtype=np.uint) for p in range(len(positions)): - for i in range(3): - for j in range(2): - p0, p1, p2 = _get_plane(domain, i, j) - if boundary_conds[i][j] == 0: - velocities[p], positions[p], old_positions[p] = _reflect(velocities[p], positions[p], old_positions[p], p0, p1, p2) - elif boundary_conds[i][j] == 1 or boundary_conds[i][j] == 2: - if _intersect(old_positions[p], positions[p], p0, p1, p2)[0]: - kept_parts[p] = 0 + counter = 0 + while not co.is_inside(positions[p], domain) and kept_parts[p]: + if counter > 10: + kept_parts[p] = 0 + break + + for i in range(3): + for j in range(2): + p0, p1, p2 = _get_plane(domain, i, j) + if boundary_conds[i][j] == 0: + velocities[p], positions[p], old_positions[p] = _reflect(velocities[p], positions[p], old_positions[p], p0, p1, p2, domain) + counter += 1 + elif boundary_conds[i][j] == 1 or boundary_conds[i][j] == 2: + if _intersect(old_positions[p], positions[p], p0, p1, p2)[0]: + kept_parts[p] = 0 N = int(sum(kept_parts)) p = 0 From 91163220b75288e1b49e2a51a854fc59ce8d771b Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 16:24:09 +0200 Subject: [PATCH 48/54] updated unit tests --- tests/unit/test_dsmc/boundary.py | 62 ++++++++++++++++---------------- 1 file changed, 31 insertions(+), 31 deletions(-) diff --git a/tests/unit/test_dsmc/boundary.py b/tests/unit/test_dsmc/boundary.py index a294398..ed983fb 100644 --- a/tests/unit/test_dsmc/boundary.py +++ b/tests/unit/test_dsmc/boundary.py @@ -18,8 +18,8 @@ def test__intersect(self): l1 = np.array([1.0, 1.0, 0.0]) p0 = np.array([0.0, -1.0, -1.0]) - p1 = np.array([0.0, -1.0, 1.0]) - p2 = np.array([0.0, 1.0, 1.0]) + p1 = np.array([0.0, 1.0, -1.0]) + p2 = np.array([0.0, -1.0, 1.0]) intersected, n_l, n_p, t = bo._intersect(l0, l1, p0, p1, p2) @@ -29,7 +29,7 @@ def test__intersect(self): self.assertEqual(2.0, n_l[1]) self.assertEqual(0.0, n_l[2]) - self.assertEqual(-4.0, n_p[0]) + self.assertEqual(4.0, n_p[0]) self.assertEqual(0.0, n_p[1]) self.assertEqual(0.0, n_p[2]) @@ -55,11 +55,11 @@ def test__reflect(self): pos_old = np.array([-1.0, -1.0, 0.0]) pos = np.array([1.0, 1.0, 0.0]) - p0 = np.array([0.0, -1.0, -1.0]) - p1 = np.array([0.0, -1.0, 1.0]) - p2 = np.array([0.0, 1.0, 1.0]) + domain = np.array([(-2.0, 0.0), (-1.0, 1.0), (-1.0, 1.0)]) + + p0, p1, p2 = bo._get_plane(domain, 0, 1) - new_vel, new_pos, new_pos_old = bo._reflect(vel, pos, pos_old, p0, p1, p2) + new_vel, new_pos, new_pos_old = bo._reflect(vel, pos, pos_old, p0, p1, p2, domain) self.assertEqual(-123.0, new_vel[0]) self.assertEqual(123.0, new_vel[1]) @@ -88,12 +88,12 @@ def test__get_plane1(self): self.assertEqual(4.0, p0[2]) self.assertEqual(0.0, p1[0]) - self.assertEqual(4.0, p1[1]) - self.assertEqual(4.0, p1[2]) + self.assertEqual(2.0, p1[1]) + self.assertEqual(7.0, p1[2]) self.assertEqual(0.0, p2[0]) - self.assertEqual(2.0, p2[1]) - self.assertEqual(7.0, p2[2]) + self.assertEqual(4.0, p2[1]) + self.assertEqual(4.0, p2[2]) p0, p1, p2 = bo._get_plane(domain, axis, 1) @@ -106,12 +106,12 @@ def test__get_plane1(self): self.assertEqual(4.0, p0[2]) self.assertEqual(1.0, p1[0]) - self.assertEqual(2.0, p1[1]) - self.assertEqual(7.0, p1[2]) + self.assertEqual(4.0, p1[1]) + self.assertEqual(4.0, p1[2]) self.assertEqual(1.0, p2[0]) - self.assertEqual(4.0, p2[1]) - self.assertEqual(4.0, p2[2]) + self.assertEqual(2.0, p2[1]) + self.assertEqual(7.0, p2[2]) def test__get_plane2(self): domain = np.array([(0, 1), (2, 4), (4, 7)]) @@ -127,13 +127,13 @@ def test__get_plane2(self): self.assertEqual(2.0, p0[1]) self.assertEqual(4.0, p0[2]) - self.assertEqual(0.0, p1[0]) + self.assertEqual(1.0, p1[0]) self.assertEqual(2.0, p1[1]) - self.assertEqual(7.0, p1[2]) + self.assertEqual(4.0, p1[2]) - self.assertEqual(1.0, p2[0]) + self.assertEqual(0.0, p2[0]) self.assertEqual(2.0, p2[1]) - self.assertEqual(4.0, p2[2]) + self.assertEqual(7.0, p2[2]) p0, p1, p2 = bo._get_plane(domain, axis, 1) @@ -145,13 +145,13 @@ def test__get_plane2(self): self.assertEqual(4.0, p0[1]) self.assertEqual(4.0, p0[2]) - self.assertEqual(1.0, p1[0]) + self.assertEqual(0.0, p1[0]) self.assertEqual(4.0, p1[1]) - self.assertEqual(4.0, p1[2]) + self.assertEqual(7.0, p1[2]) - self.assertEqual(0.0, p2[0]) + self.assertEqual(1.0, p2[0]) self.assertEqual(4.0, p2[1]) - self.assertEqual(7.0, p2[2]) + self.assertEqual(4.0, p2[2]) def test__get_plane3(self): domain = np.array([(0, 1), (2, 4), (4, 7)]) @@ -167,12 +167,12 @@ def test__get_plane3(self): self.assertEqual(2.0, p0[1]) self.assertEqual(4.0, p0[2]) - self.assertEqual(1.0, p1[0]) - self.assertEqual(2.0, p1[1]) + self.assertEqual(0.0, p1[0]) + self.assertEqual(4.0, p1[1]) self.assertEqual(4.0, p1[2]) - self.assertEqual(0.0, p2[0]) - self.assertEqual(4.0, p2[1]) + self.assertEqual(1.0, p2[0]) + self.assertEqual(2.0, p2[1]) self.assertEqual(4.0, p2[2]) p0, p1, p2 = bo._get_plane(domain, axis, 1) @@ -185,11 +185,11 @@ def test__get_plane3(self): self.assertEqual(2.0, p0[1]) self.assertEqual(7.0, p0[2]) - self.assertEqual(0.0, p1[0]) - self.assertEqual(4.0, p1[1]) + self.assertEqual(1.0, p1[0]) + self.assertEqual(2.0, p1[1]) self.assertEqual(7.0, p1[2]) - self.assertEqual(1.0, p2[0]) - self.assertEqual(2.0, p2[1]) + self.assertEqual(0.0, p2[0]) + self.assertEqual(4.0, p2[1]) self.assertEqual(7.0, p2[2]) \ No newline at end of file From 10e7a106efbee1898b929816fc2999df61d790fb Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 16:59:46 +0200 Subject: [PATCH 49/54] adde values setting for boundary --- dsmc/boundary.py | 21 +++++++++++++++++++-- 1 file changed, 19 insertions(+), 2 deletions(-) diff --git a/dsmc/boundary.py b/dsmc/boundary.py index ff4cd87..6d1453a 100644 --- a/dsmc/boundary.py +++ b/dsmc/boundary.py @@ -159,5 +159,22 @@ def __init__(self): self.boundary_conds = np.array([[0, 0], [0, 0], [0, 0]], dtype=np.uint) # 0 = ela, 1 = open, 2 = inflow self.domain = None - def boundary(self, velocities, positions, old_positions): - return _boundary(velocities, positions, old_positions, self.domain, self.boundary_conds) \ No newline at end of file + def boundary(self, velocities, positions, old_positions): + return _boundary(velocities, positions, old_positions, self.domain, self.boundary_conds) + + def set_bc_type(self, boundary, bc_type): + bound = _get_boundary(boundary) + bc = _get_bc_type(bc_type) + + self.boundary_conds[bound[0]][bound[1]] = bc + + print("boundary [" + boundary + "] set to [" + bc_type + "]") + + def set_bc_values(self, boundary, T, n, u): + i, j = _get_boundary(boundary) + + self.T[i][j] = T + self.n[i][j] = n + self.u[i][j] = u + + print("boundary [" + boundary + "] set to values T : {}, n : {}, u : {}".format(T, n, u)) \ No newline at end of file From 0ea2de06a203d2482ce2dec21b7bec6782d9792b Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 17:03:15 +0200 Subject: [PATCH 50/54] updated boudndary test case --- tests/domain/boundary.py | 32 ++++++++++++++++++++++++++++---- 1 file changed, 28 insertions(+), 4 deletions(-) diff --git a/tests/domain/boundary.py b/tests/domain/boundary.py index ae88535..477af40 100644 --- a/tests/domain/boundary.py +++ b/tests/domain/boundary.py @@ -7,6 +7,33 @@ import time import numpy as np +def create_pos_and_vels(): + positions = np.zeros((6, 3)) + velocitiies = np.zeros((6, 3)) + + # x + positions[0][0] = -0.75 + velocitiies[0][0] = 100.0 + + positions[1][0] = 0.75 + velocitiies[1][0] = -100.0 + + # y + positions[2][1] = -0.75 + velocitiies[2][1] = 100.0 + + positions[3][1] = 0.75 + velocitiies[3][1] = -100.0 + + # z + positions[4][2] = -0.75 + velocitiies[4][2] = 100.0 + + positions[5][2] = 0.75 + velocitiies[5][2] = -100.0 + + return (velocitiies, positions) + if __name__ == '__main__': # general parameters boundary = bo.Boundary() @@ -33,11 +60,8 @@ particles.create_particles(domain, mass, T, N) velocities, positions = particles.VelPos - - #velocities = np.array([[150.0, 150.0, 0.0]]) - #positions = np.array([[0.0, 0.0, 0.0]]) - particles.VelPos = velocities, positions + particles.VelPos = create_pos_and_vels() # start timing start_time = time.time() From b5bfa3166af5912ed8343bab219d254407664816 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 17:41:01 +0200 Subject: [PATCH 51/54] updated test case --- tests/domain/boundary.py | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/tests/domain/boundary.py b/tests/domain/boundary.py index 477af40..5d70ba9 100644 --- a/tests/domain/boundary.py +++ b/tests/domain/boundary.py @@ -51,7 +51,7 @@ def create_pos_and_vels(): tree_outer = oc.Octree() outer_pos = np.array([[-1.0, -1.0, -1.0], [1.0, 1.0, 1.0]]) tree_outer.build(outer_pos) - wrt.write_buttom_leafs(tree_outer, "outer_box.vtu") + #wrt.write_buttom_leafs(tree_outer, "outer_box.vtu") # setup solver boundary.domain = domain @@ -61,20 +61,31 @@ def create_pos_and_vels(): velocities, positions = particles.VelPos - particles.VelPos = create_pos_and_vels() + #particles.VelPos = create_pos_and_vels() # start timing start_time = time.time() wrt.write_positions(particles, "pos_{}.csv".format(0)) + E0 = 0.0 + + for vel in particles.Vel: + E0 += vel.dot(vel) + for it in range(niter): - print("iteration {:4}/{}, N particles {}/{}".format(it + 1, niter, particles.N, N), end="\r", flush=True) + E = 0.0 + velocities, positions, old_positions = ds._push(particles.Vel, particles.Pos, dt) velocities, positions, old_positions = boundary.boundary(velocities, positions, old_positions) particles.VelPos = (velocities, positions) wrt.write_positions(particles, "pos_{}.csv".format(it + 1)) + + for vel in particles.Vel: + E += vel.dot(vel) + + print("iteration {:4}/{}, N particles {}/{}, Efrac {}".format(it + 1, niter, particles.N, N, E/E0), end="\r", flush=True) print("") From 0ca8b968fceea131138c043befcb05daebfaf57b Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 7 Oct 2022 17:41:28 +0200 Subject: [PATCH 52/54] fixes velocity calculation in new boundary --- dsmc/boundary.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dsmc/boundary.py b/dsmc/boundary.py index 6d1453a..b2f979a 100644 --- a/dsmc/boundary.py +++ b/dsmc/boundary.py @@ -53,7 +53,7 @@ def _reflect(vel, pos, pos_old, p0, p1, p2, domain): pos_old = p n_r = _calc_nr(n_l, n_p) pos = pos_old + (1.0 - t)*n_r - vel = pos_old + (np.linalg.norm(vel) / np.linalg.norm(n_r))*n_r + vel = (np.linalg.norm(vel) / np.linalg.norm(n_r))*n_r return (vel, pos, pos_old) From 58455b2f5d397523a79c22f30cc6a5ca582c8f08 Mon Sep 17 00:00:00 2001 From: LeoBasov <40352679+LeoBasov@users.noreply.github.com> Date: Sun, 9 Oct 2022 20:09:49 +0200 Subject: [PATCH 53/54] updated changelog and setup.py version number --- CHANGELOG.md | 8 ++++++++ setup.py | 2 +- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 49a2f7d..761b90c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.9.0] - 2022-10-09 +## Added +- mesh2d as for diagnostic purposes +- particle writer +- docs +- hypersonic test case +- backup boundary module + ## [0.8.0] - 2022-10-03 ### Added - planar writer diff --git a/setup.py b/setup.py index 2b71d95..800507d 100644 --- a/setup.py +++ b/setup.py @@ -1,7 +1,7 @@ from setuptools import setup setup( name='dsmc', - version='0.8.0', + version='0.9.0', author='Leo Basov', python_requires='>=3.10, <4', packages=["dsmc", "dsmc.writer", "dsmc.mesh"], From 3193f83a0b605266735ce1f233c45fc727cc1f6e Mon Sep 17 00:00:00 2001 From: LeoBasov <40352679+LeoBasov@users.noreply.github.com> Date: Sun, 9 Oct 2022 20:30:48 +0200 Subject: [PATCH 54/54] moved boundary interaction --- dsmc/dsmc.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index 5cca160..01209ae 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -243,11 +243,13 @@ def advance(self, dt, collisions=True, octree=True): velocities, positions, old_positions = _push(self.particles.Vel, self.particles.Pos, dt) velocities, positions, old_positions = _check_positions(velocities, positions, old_positions, self.domain) - velocities, positions, old_positions = _boundary(velocities, positions, old_positions, self.domain, self.boundary_conds) + for obj in self.objects: velocities, positions, old_positions = _object(velocities, positions, old_positions, obj) + velocities, positions, old_positions = _boundary(velocities, positions, old_positions, self.domain, self.boundary_conds) + if octree: self.octree.build(positions) if collisions and octree: