From 4601d6e4d7caf0e51b8795a0ebd288eda11260ca Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Wed, 28 Sep 2022 11:25:08 +0200 Subject: [PATCH 01/33] implemented writing of particle numbers to ocetree --- dsmc/writer.py | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/dsmc/writer.py b/dsmc/writer.py index c181d62..03c6d7f 100644 --- a/dsmc/writer.py +++ b/dsmc/writer.py @@ -19,8 +19,9 @@ def _wrtie_body(f, octree): f.write("\n") f.write("\n".format(len(leaf_ids) * 8, len(leaf_ids))) - _write_points(f, octree, leaf_ids); - _write_cells(f, octree, leaf_ids); + _write_points(f, octree, leaf_ids) + _write_cells(f, octree, leaf_ids) + _write_cell_data(f, octree, leaf_ids) f.write("\n") f.write("\n") @@ -93,6 +94,16 @@ def _write_cells(f, octree, leaf_ids): f.write("\n") f.write("\n") + +def _write_cell_data(f, octree, leaf_ids): + f.write("\n") + f.write("\n") + + for i in leaf_ids: + f.write("{} ".format(octree.leafs[i].number_elements)) + + f.write("\n") + f.write("\n") def _write_footer(f): f.write("\n") From a964fe3d50f8f8f9fa0721f46dba959227c12e8c Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Wed, 28 Sep 2022 11:36:53 +0200 Subject: [PATCH 02/33] removed unecessary print statement --- tests/unit/test_dsmc/octree.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/unit/test_dsmc/octree.py b/tests/unit/test_dsmc/octree.py index 5eb9a57..f04fd22 100644 --- a/tests/unit/test_dsmc/octree.py +++ b/tests/unit/test_dsmc/octree.py @@ -315,8 +315,6 @@ def test__create_combined_boxes_4(self): self.assertEqual(oc.get_V(box), V) self.assertEqual(8, len(boxes_new)) - - print(boxes_new) class TestOctreeOctree(unittest.TestCase): def test_build(self): From 80381fc142d5d109292f881c214753e6f0ab556f Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Wed, 28 Sep 2022 11:37:26 +0200 Subject: [PATCH 03/33] implemented outflow boundary --- dsmc/dsmc.py | 40 +++++++++++++++++++++++++++++++--------- 1 file changed, 31 insertions(+), 9 deletions(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index 3c13811..1fe95b0 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -11,18 +11,39 @@ def _push(velocities, positions, dt): return positions @njit(parallel=True) -def _boundary(velocities, positions, domain): +def _boundary(velocities, positions, domain, boundary_conds): + kept_parts = np.ones(positions.shape[0], dtype=np.uint) + for p in prange(len(positions)): - while not oc._is_inside(positions[p], domain): + while not oc._is_inside(positions[p], domain) and kept_parts[p]: for i in range(3): if positions[p][i] < domain[i][0]: - positions[p][i] = 2.0 * domain[i][0] - positions[p][i] - velocities[p][i] *= -1.0 + if boundary_conds[i][0] == 0: + positions[p][i] = 2.0 * domain[i][0] - positions[p][i] + velocities[p][i] *= -1.0 + elif boundary_conds[i][0] == 2: + kept_parts[p] = 0 if positions[p][i] > domain[i][1]: - positions[p][i] = 2.0 * domain[i][1] - positions[p][i] - velocities[p][i] *= -1.0 - - return (velocities, positions) + if boundary_conds[i][1] == 0: + positions[p][i] = 2.0 * domain[i][1] - positions[p][i] + velocities[p][i] *= -1.0 + elif boundary_conds[i][1] == 2: + kept_parts[p] = 0 + + N = sum(kept_parts) + p = 0 + new_velocities = np.empty((N, 3)) + new_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] + p += 1 + else: + continue + + return (new_velocities, new_positions) @njit def _calc_prob(rel_vel : float, sigma_T : float, Vc : float, dt : float, w : float, N : int) -> np.single: @@ -100,6 +121,7 @@ def clear(self): self.octree = oc.Octree() self.w = None self.domain = None + self.boundary_conds = np.array([[0, 0], [0, 0], [0, 0]], dtype=np.uint) # 0 = ela, 1 = periodic, 2 = open, 3 = inflow self.sigma_T = 3.631681e-19 self.mass = None @@ -116,7 +138,7 @@ def advance(self, dt, collisions=True, octree=True): if collisions and octree: self.particles.VelPos = (self._update_velocities(dt), self.particles.Pos) positions = _push(self.particles.Vel, self.particles.Pos, dt) - self.particles.VelPos = _boundary(self.particles.Vel, positions, self.domain) + self.particles.VelPos = _boundary(self.particles.Vel, positions, self.domain, self.boundary_conds) def _update_velocities(self, dt): Nleafs : int = len(self.octree.leafs) From 0dc9713112cf1f989de33e5368189dfb90107d4c Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Wed, 28 Sep 2022 11:50:19 +0200 Subject: [PATCH 04/33] added setting of boundary conditions --- dsmc/dsmc.py | 36 +++++++++++++++++++++++++++++++++++- 1 file changed, 35 insertions(+), 1 deletion(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index 1fe95b0..1b43464 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -112,6 +112,32 @@ def _update_vels(permutations : np.ndarray, velocities : np.ndarray, mass : floa return velocities +@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 == "period": + return 1 + elif bc_type == "open": + return 2 + elif bc_type == "inflow": + return 3 + class DSMC: def __init__(self): self.clear() @@ -129,7 +155,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 created") + raise Exception("no particles in domain") if self.w == None: raise Exception("particle weight not set") @@ -175,3 +201,11 @@ def set_mass(self, mass): def set_weight(self, w): self.octree.w = w self.w = w + + def set_boundary(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 + "]") From 37277b5c3b8602fec9e383c27c6d75d5e07156e9 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Wed, 28 Sep 2022 15:51:58 +0200 Subject: [PATCH 05/33] added open boundary exapmle --- examples/push_bound_open_test.py | 52 ++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) create mode 100644 examples/push_bound_open_test.py diff --git a/examples/push_bound_open_test.py b/examples/push_bound_open_test.py new file mode 100644 index 0000000..8313b44 --- /dev/null +++ b/examples/push_bound_open_test.py @@ -0,0 +1,52 @@ +import dsmc +import dsmc.octree as oc +import dsmc.writer as wrt +import numpy as np + +if __name__ == '__main__': + # general parameters + solver = dsmc.DSMC() + tree = oc.Octree() + domain = [(-0.5, 0.5), (-0.5, 0.5), (-0.5, 0.5)] + positions = np.array([(-0.5, -0.5, -0.5), (0.5, 0.5, 0.5)]) + dt = 1e-5 + w = 1e+17 + n = 1e18 + mass = 6.6422e-26 + T = 300 + niter = 300 + + # setup solver + solver.set_domain(domain) + solver.set_weight(w) + solver.set_mass(mass) + + solver.create_particles(domain, T, n) + + solver.set_boundary("xmax", "open") + solver.set_boundary("xmin", "open") + solver.set_boundary("ymax", "open") + solver.set_boundary("ymin", "open") + solver.set_boundary("zmax", "open") + solver.set_boundary("zmin", "open") + + tree.build(positions) + + for it in range(niter): + print("iteration {:4}/{:4}".format(it + 1, niter), end="\r", flush=True) + + try: + solver.advance(dt, collisions=False, octree=False) + except Exception as e: + print(e) + break + + with open("particles_{}.csv".format(it), "w") as file: + file.write("x, y, z\n") + for pos in solver.particles.Pos: + file.write("{},{},{}\n".format(pos[0], pos[1], pos[2])) + + wrt.write_buttom_leafs(tree, "octree.vtu") + + + print("done") \ No newline at end of file From 18274f741ca0380c525e54e6432f7c3ce325ed2f Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Wed, 28 Sep 2022 15:59:18 +0200 Subject: [PATCH 06/33] changed open boundary index --- examples/push_bound_open_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/push_bound_open_test.py b/examples/push_bound_open_test.py index 8313b44..362817d 100644 --- a/examples/push_bound_open_test.py +++ b/examples/push_bound_open_test.py @@ -14,7 +14,7 @@ n = 1e18 mass = 6.6422e-26 T = 300 - niter = 300 + niter = 1000 # setup solver solver.set_domain(domain) From 8f904d3f66f6a50448ec5c1dd924f2bfa8a731dd Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Wed, 28 Sep 2022 15:59:31 +0200 Subject: [PATCH 07/33] changed open boundary index --- dsmc/dsmc.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index 1b43464..a4b92d8 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -21,13 +21,13 @@ def _boundary(velocities, positions, domain, boundary_conds): if boundary_conds[i][0] == 0: positions[p][i] = 2.0 * domain[i][0] - positions[p][i] velocities[p][i] *= -1.0 - elif boundary_conds[i][0] == 2: + elif boundary_conds[i][0] == 1: kept_parts[p] = 0 if positions[p][i] > domain[i][1]: if boundary_conds[i][1] == 0: positions[p][i] = 2.0 * domain[i][1] - positions[p][i] velocities[p][i] *= -1.0 - elif boundary_conds[i][1] == 2: + elif boundary_conds[i][1] == 1: kept_parts[p] = 0 N = sum(kept_parts) @@ -131,12 +131,10 @@ def _get_boundary(boundary): def _get_bc_type(bc_type): if bc_type == "ela": return 0 - elif bc_type == "period": - return 1 elif bc_type == "open": - return 2 + return 1 elif bc_type == "inflow": - return 3 + return 2 class DSMC: def __init__(self): @@ -147,7 +145,7 @@ def clear(self): self.octree = oc.Octree() self.w = None self.domain = None - self.boundary_conds = np.array([[0, 0], [0, 0], [0, 0]], dtype=np.uint) # 0 = ela, 1 = periodic, 2 = open, 3 = inflow + self.boundary_conds = np.array([[0, 0], [0, 0], [0, 0]], dtype=np.uint) # 0 = ela, 1 = open, 2 = inflow self.sigma_T = 3.631681e-19 self.mass = None From b3f922f08e1ace901c11d18abdde9835b4c1a02a Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Wed, 28 Sep 2022 16:15:15 +0200 Subject: [PATCH 08/33] implemented calc_vp --- dsmc/particles.py | 4 ++++ tests/unit/test_dsmc/particles.py | 8 ++++++++ 2 files changed, 12 insertions(+) diff --git a/dsmc/particles.py b/dsmc/particles.py index 7e709c3..467f119 100644 --- a/dsmc/particles.py +++ b/dsmc/particles.py @@ -4,6 +4,10 @@ kb = 1.380649e-23 +@njit +def calc_vp(T, mass): + return np.sqrt(2*kb*T/mass) + @njit def box_muller(T): """ diff --git a/tests/unit/test_dsmc/particles.py b/tests/unit/test_dsmc/particles.py index 1921691..0adf835 100644 --- a/tests/unit/test_dsmc/particles.py +++ b/tests/unit/test_dsmc/particles.py @@ -1,5 +1,6 @@ import unittest from dsmc import particles as pa +import numpy as np class TestParticles(unittest.TestCase): def test_box_muller(self): @@ -74,3 +75,10 @@ def test_create_particles(self): self.assertTrue(particles.Pos[i][0] >= x[0] and particles.Pos[i][0] <= x[1]) self.assertTrue(particles.Pos[i][1] >= y[0] and particles.Pos[i][1] <= y[1]) self.assertTrue(particles.Pos[i][2] >= z[0] and particles.Pos[i][2] <= z[1]) + + def test_calc_vp(self): + T = 300 + mass = 1e-26 + vp = pa.calc_vp(T, mass) + + self.assertEqual(np.sqrt(2.0*pa.kb*T/mass), vp) From 77c6f78be509d32cd01f428b048458a64162aea7 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Wed, 28 Sep 2022 16:53:08 +0200 Subject: [PATCH 09/33] added inflow boundary --- dsmc/dsmc.py | 29 ++++++++++++++++------------- dsmc/particles.py | 27 ++++++++++++++++++++++----- tests/unit/test_dsmc/particles.py | 11 +++++++---- 3 files changed, 45 insertions(+), 22 deletions(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index a4b92d8..e79e2e0 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -21,13 +21,13 @@ def _boundary(velocities, positions, domain, boundary_conds): if boundary_conds[i][0] == 0: positions[p][i] = 2.0 * domain[i][0] - positions[p][i] velocities[p][i] *= -1.0 - elif boundary_conds[i][0] == 1: + 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: positions[p][i] = 2.0 * domain[i][1] - positions[p][i] velocities[p][i] *= -1.0 - elif boundary_conds[i][1] == 1: + elif boundary_conds[i][1] == 1 or boundary_conds[i][0] == 2: kept_parts[p] = 0 N = sum(kept_parts) @@ -135,6 +135,12 @@ def _get_bc_type(bc_type): 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)) class DSMC: def __init__(self): @@ -146,6 +152,7 @@ def clear(self): self.w = None self.domain = None self.boundary_conds = np.array([[0, 0], [0, 0], [0, 0]], dtype=np.uint) # 0 = ela, 1 = open, 2 = inflow + self.boundary = Boundary() self.sigma_T = 3.631681e-19 self.mass = None @@ -156,6 +163,11 @@ def advance(self, dt, collisions=True, octree=True): raise Exception("no particles in domain") if self.w == None: raise Exception("particle weight not set") + + for i in range(3): + for j in range(2): + 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) @@ -173,20 +185,11 @@ def _update_velocities(self, dt): return _update_vels(self.octree.permutations, self.particles.Vel, self.mass, self.sigma_T, dt, self.w, elem_offsets, number_elements, number_children, cell_boxes, Nleafs) - def create_particles(self, box, T, n, u = None): + def create_particles(self, box, T, n, u = np.zeros(3)): box = np.array(box) N = int(round(oc.get_V(box) * n / self.w)) print("creating {} particles".format(N)) - self.particles.create_particles(box, self.mass, T, N) - - if u is not None: - velocities = self.particles.Vel - u = np.array(u) - - for i in range(len(velocities)): - velocities[i] += u - - self.particles.VelPos = (velocities, self.particles.Pos) + self.particles.create_particles(box, self.mass, T, N, u) print("now containing {} particles, {} total".format(N, self.particles.N)) diff --git a/dsmc/particles.py b/dsmc/particles.py index 467f119..3f554e5 100644 --- a/dsmc/particles.py +++ b/dsmc/particles.py @@ -1,6 +1,7 @@ import math import numpy as np from numba import njit +from . import octree as oc kb = 1.380649e-23 @@ -61,11 +62,11 @@ def get_vel(T, mass): return v * x2velocity(box_muller(T), mass) / np.linalg.norm(v) @njit -def get_velocities(T, mass, N): +def get_velocities(T, mass, N, u): velocities = np.empty((N, 3)) for i in range(N): - velocities[i] = get_vel(T, mass) + velocities[i] = get_vel(T, mass) + u return velocities @@ -132,12 +133,28 @@ def VelPos(self, vel_pos): self._positions = vel_pos[1] self._N = len(self._positions) - def create_particles(self, X, mass, T, N): + def create_particles(self, X, mass, T, N, u = np.zeros(3)): if self._N == 0: - self._velocities = get_velocities(T, mass, N) + self._velocities = get_velocities(T, mass, N, u) self._positions = calc_positions(X[0], X[1], X[2], N) self._N = N else: - self._velocities = np.concatenate((self._velocities, get_velocities(T, mass, N))) + 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._N += N + + + def inflow(self, mass, T, u, n, w, dt, domain, axis, minmax): + L = max(calc_vp(T, mass) * dt * 10, np.linalg.norm(u) * dt) + box = np.copy(domain) + + if minmax == 0: + box[axis][1] = box[axis][0] + box[axis][0] = box[axis][1] - L + elif minmax == 1: + box[axis][0] = box[axis][1] + box[axis][1] = box[axis][0] + L + + N = int(round(oc.get_V(box) * n / w)) + + self.create_particles(box, mass, T, N, u) diff --git a/tests/unit/test_dsmc/particles.py b/tests/unit/test_dsmc/particles.py index 0adf835..700c46f 100644 --- a/tests/unit/test_dsmc/particles.py +++ b/tests/unit/test_dsmc/particles.py @@ -21,7 +21,8 @@ def test_get_velocities(self): T = 300 mass = 1.0e-26 N = 1000 - velocities = pa.get_velocities(T, mass, N) + u = np.zeros(3) + velocities = pa.get_velocities(T, mass, N, u) self.assertEqual(N, len(velocities)) @@ -29,7 +30,8 @@ def test_calc_temperature(self): T = 300 mass = 1.0e-26 N = 10000 - velocities = pa.get_velocities(T, mass, N) + u = np.zeros(3) + velocities = pa.get_velocities(T, mass, N, u) T_new = pa.calc_temperature(velocities, mass) diff = abs((T_new - T)/T) @@ -52,9 +54,10 @@ def test_create_particles(self): N = 1000 mass = 1.0e-26 T = 300 + u = np.zeros(3) particles = pa.Particles() - particles.create_particles(X, mass, T, N) + particles.create_particles(X, mass, T, N, u) self.assertEqual(N, len(particles.Pos)) self.assertEqual(N, len(particles.Vel)) @@ -65,7 +68,7 @@ def test_create_particles(self): self.assertTrue(particles.Pos[i][1] >= y[0] and particles.Pos[i][1] <= y[1]) self.assertTrue(particles.Pos[i][2] >= z[0] and particles.Pos[i][2] <= z[1]) - particles.create_particles(X, mass, T, N) + particles.create_particles(X, mass, T, N, u) self.assertEqual(2*N, len(particles.Pos)) self.assertEqual(2*N, len(particles.Vel)) From a1919d7fdafc491eb66cbe095fa9592550b118f0 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Wed, 28 Sep 2022 16:53:52 +0200 Subject: [PATCH 10/33] added test for inflow boundary --- examples/push_bound_inflow_test.py | 56 ++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) create mode 100644 examples/push_bound_inflow_test.py diff --git a/examples/push_bound_inflow_test.py b/examples/push_bound_inflow_test.py new file mode 100644 index 0000000..37abe50 --- /dev/null +++ b/examples/push_bound_inflow_test.py @@ -0,0 +1,56 @@ +import dsmc +import dsmc.octree as oc +import dsmc.writer as wrt +import numpy as np + +if __name__ == '__main__': + # general parameters + solver = dsmc.DSMC() + tree = oc.Octree() + domain = [(-0.5, 0.5), (-0.5, 0.5), (-0.5, 0.5)] + positions = np.array([(-0.5, -0.5, -0.5), (0.5, 0.5, 0.5)]) + dt = 1e-5 + w = 1e+16 + n = 1e18 + mass = 6.6422e-26 + T = 300 + niter = 1000 + + # setup solver + solver.set_domain(domain) + solver.set_weight(w) + solver.set_mass(mass) + + solver.create_particles(domain, T, n) + + solver.set_boundary("xmax", "open") + + solver.set_boundary("ymax", "open") + solver.set_boundary("ymin", "open") + solver.set_boundary("zmax", "open") + solver.set_boundary("zmin", "open") + + solver.set_boundary("xmin", "inflow") + + solver.boundary.u[0][0] = np.array([1000, 0, 0]) + + tree.build(positions) + + for it in range(niter): + print("iteration {:4}/{:4}".format(it + 1, niter), end="\r", flush=True) + + try: + solver.advance(dt, collisions=False, octree=False) + except Exception as e: + print(e) + break + + with open("particles_{}.csv".format(it), "w") as file: + file.write("x, y, z\n") + for pos in solver.particles.Pos: + file.write("{},{},{}\n".format(pos[0], pos[1], pos[2])) + + wrt.write_buttom_leafs(tree, "octree.vtu") + + + print("done") \ No newline at end of file From 5b1b892d07d9a52ba4cdaf363f86965cb7a9de01 Mon Sep 17 00:00:00 2001 From: Basov Date: Thu, 29 Sep 2022 09:42:18 +0200 Subject: [PATCH 11/33] implemented settning of boundary values --- dsmc/dsmc.py | 11 ++++++++++- examples/push_bound_inflow_test.py | 15 ++++++++------- examples/push_bound_open_test.py | 12 ++++++------ 3 files changed, 24 insertions(+), 14 deletions(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index e79e2e0..3b59096 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -203,10 +203,19 @@ def set_weight(self, w): self.octree.w = w self.w = w - def set_boundary(self, boundary, bc_type): + 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.boundary.T[i][j] = T + self.boundary.n[i][j] = n + self.boundary.u[i][j] = u + + print("boundary [" + boundary + "] set to values T : {}, n : {}, u : {}".format(T, n, u)) diff --git a/examples/push_bound_inflow_test.py b/examples/push_bound_inflow_test.py index 37abe50..d7fdfdb 100644 --- a/examples/push_bound_inflow_test.py +++ b/examples/push_bound_inflow_test.py @@ -12,6 +12,7 @@ dt = 1e-5 w = 1e+16 n = 1e18 + u = np.array([1000.0, 0.0, 0.0]) mass = 6.6422e-26 T = 300 niter = 1000 @@ -23,16 +24,16 @@ solver.create_particles(domain, T, n) - solver.set_boundary("xmax", "open") + solver.set_bc_type("xmax", "open") - solver.set_boundary("ymax", "open") - solver.set_boundary("ymin", "open") - solver.set_boundary("zmax", "open") - solver.set_boundary("zmin", "open") + solver.set_bc_type("ymax", "open") + solver.set_bc_type("ymin", "open") + solver.set_bc_type("zmax", "open") + solver.set_bc_type("zmin", "open") - solver.set_boundary("xmin", "inflow") + solver.set_bc_type("xmin", "inflow") - solver.boundary.u[0][0] = np.array([1000, 0, 0]) + solver.set_bc_values("xmin", T, n, u) tree.build(positions) diff --git a/examples/push_bound_open_test.py b/examples/push_bound_open_test.py index 362817d..b899097 100644 --- a/examples/push_bound_open_test.py +++ b/examples/push_bound_open_test.py @@ -23,12 +23,12 @@ solver.create_particles(domain, T, n) - solver.set_boundary("xmax", "open") - solver.set_boundary("xmin", "open") - solver.set_boundary("ymax", "open") - solver.set_boundary("ymin", "open") - solver.set_boundary("zmax", "open") - solver.set_boundary("zmin", "open") + solver.set_bc_type("xmax", "open") + solver.set_bc_type("xmin", "open") + solver.set_bc_type("ymax", "open") + solver.set_bc_type("ymin", "open") + solver.set_bc_type("zmax", "open") + solver.set_bc_type("zmin", "open") tree.build(positions) From 9f2eaba84936333ff558c114bcf52157b9521d99 Mon Sep 17 00:00:00 2001 From: Basov Date: Thu, 29 Sep 2022 10:19:13 +0200 Subject: [PATCH 12/33] added boundary check --- dsmc/dsmc.py | 31 ++++++++++++++++++++++++++++--- 1 file changed, 28 insertions(+), 3 deletions(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index 3b59096..bfc404f 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -4,13 +4,13 @@ from . import particles as prt from . import octree as oc -@njit(parallel=True) +@njit def _push(velocities, positions, dt): for p in prange(len(positions)): positions[p] = positions[p] + velocities[p]*dt return positions -@njit(parallel=True) +@njit def _boundary(velocities, positions, domain, boundary_conds): kept_parts = np.ones(positions.shape[0], dtype=np.uint) @@ -45,6 +45,29 @@ def _boundary(velocities, positions, domain, boundary_conds): return (new_velocities, new_positions) +@njit +def _check_positions(velocities, positions, old_positions, domain): + kept_parts = np.ones(positions.shape[0], dtype=np.uint) + + for i in prange(positions.shape[0]): + if (not oc._is_inside(positions[i], domain)) and (not oc._is_inside(old_positions[i], domain)): + 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 _calc_prob(rel_vel : float, sigma_T : float, Vc : float, dt : float, w : float, N : int) -> np.single: """ @@ -173,8 +196,10 @@ 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 = _boundary(self.particles.Vel, positions, self.domain, self.boundary_conds) + 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) def _update_velocities(self, dt): Nleafs : int = len(self.octree.leafs) From 6209448c907e488e9c4ef5a032d53dbb39c6df44 Mon Sep 17 00:00:00 2001 From: Basov Date: Thu, 29 Sep 2022 11:01:01 +0200 Subject: [PATCH 13/33] fixed bug in particle creation --- dsmc/particles.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dsmc/particles.py b/dsmc/particles.py index 3f554e5..f018b92 100644 --- a/dsmc/particles.py +++ b/dsmc/particles.py @@ -58,7 +58,7 @@ def get_vel(T, mass): ------- velocity : np.array, shape = (3, 1) """ - v = np.random.random(3) + v = np.random.random(3)*2.0 - np.ones(3) return v * x2velocity(box_muller(T), mass) / np.linalg.norm(v) @njit From 404d9a320da3245e91704579b03430f25b650850 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Thu, 29 Sep 2022 11:54:23 +0200 Subject: [PATCH 14/33] updated shock tube test case --- examples/shock_tube.py | 35 +++++++++++++++++++++-------------- 1 file changed, 21 insertions(+), 14 deletions(-) diff --git a/examples/shock_tube.py b/examples/shock_tube.py index 37ef606..23e468a 100644 --- a/examples/shock_tube.py +++ b/examples/shock_tube.py @@ -1,7 +1,13 @@ import dsmc import dsmc.diagnostics as dia +import time + +def write2file(solver, Nbins): + # open files + n_file = open("n.csv", "w") + T_file = open("T.csv", "w") + p_file = open("p.csv", "w") -def write2file(solver, n_file, T_file, p_file, Nbins): bins, box, x = dia.sort_bin(solver.particles.Pos, 2, Nbins) n = dia.calc_n(bins, box, 2, solver.w) T = dia.calc_T(bins, solver.particles.Vel, mass) @@ -15,16 +21,21 @@ def write2file(solver, n_file, T_file, p_file, Nbins): n_file.write("\n") T_file.write("\n") p_file.write("\n") + + # close files + n_file.close() + T_file.close() + p_file.close() if __name__ == '__main__': # general parameters solver = dsmc.DSMC() domain = [(-0.0001, 0.0001), (-0.0001, 0.0001), (0.0, 0.1)] dt = 1e-7 - w = 1e+9 + w = 1e+8 mass = 6.6422e-26 niter = 300 - Nbins = 100 + Nbins = 1000 # high denisty particles nhigh = 2.41432e+22 @@ -44,20 +55,16 @@ def write2file(solver, n_file, T_file, p_file, Nbins): solver.create_particles(Boxlow, Tlow, nlow) solver.create_particles(Boxhigh, Thigh, nhigh) - # open files - n_file = open("n.csv", "w") - T_file = open("T.csv", "w") - p_file = open("p.csv", "w") + # 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) - write2file(solver, n_file, T_file, p_file, Nbins) - - # close files - n_file.close() - T_file.close() - p_file.close() - + print("") + print("--- %s seconds ---" % (time.time() - start_time)) + + write2file(solver, Nbins) + print('done') From a01ce944efc4112b2050365e7389a8ff1e8e5330 Mon Sep 17 00:00:00 2001 From: Basov Date: Thu, 29 Sep 2022 12:57:55 +0200 Subject: [PATCH 15/33] added tools --- tools/sod_analytical.py | 96 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 96 insertions(+) create mode 100644 tools/sod_analytical.py diff --git a/tools/sod_analytical.py b/tools/sod_analytical.py new file mode 100644 index 0000000..00dbe61 --- /dev/null +++ b/tools/sod_analytical.py @@ -0,0 +1,96 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu Sep 29 12:03:03 2022 + +@author: Leo Basov +""" + +import numpy as np +import matplotlib.pyplot as plt +import argparse + +def _calc_u3(p1, p3, rhoL, Gamma, gamma, beta): + num = (1 - Gamma**2) * p1**(1.0/gamma) + denum = Gamma**2 * rhoL + return (p1**beta - p3**beta) * np.sqrt(num / denum) + +def _calc_u4(p3, p5, rhoR, Gamma): + return (p3 - p5) * np.sqrt((1 - Gamma) / (rhoR*(p3 - Gamma*p5))) + +def _calc_p3(p, rho, Gamma, gamma, beta): + p1 = p[0] + p3 = p[0] + p5 = p[4] + dp = 1e-5 + rhoL = rho[0] + rhoR = rho[-1] + err = 1e+5 + err_last = 1e+6 + + while err < err_last: + err_last = err + err = abs(_calc_u3(p1, p3, rhoL, Gamma, gamma, beta) - _calc_u4(p3, p5, rhoR, Gamma)) + p3 -= dp + + return p3 + +if __name__ == '__main__': + parser = argparse.ArgumentParser() + parser.add_argument('t', type=float, help='time of results') + + args = parser.parse_args() + + # number points per segment + N = 10 + + # gas = Ar + gamma = 1.67 + Gamma = (gamma - 1.0) / (gamma + 1.0) + beta = (gamma - 1.0) / (2.0 * gamma) + + # boundary conditions + rho = np.zeros(5) + p = np.zeros(5) + u = np.zeros(5) + + rho[0] = 1.0 + p[0] = 1.0 + u[0] = 0.0 + + rho[-1] = 0.125 + p[-1] = 0.1 + u[-1] = 0.0 + + # calculating states + p[2] = _calc_p3(p, rho, Gamma, gamma, beta) + p[3] = p[2] + + # speed of characterisics + u[1] = np.sqrt(gamma * p[0] / rho[0]) # speed of rarefication going left + u[2] = 0 + u[3] = (p[2] - p[4]) / np.sqrt((rho[4]/2) * ((gamma + 1) * p[2] + (gamma - 1) * p[4])) + u[4] = np.sqrt(gamma * p[-1] / rho[-1]) # speed of pressure increase going right + + rho[2] = rho[0]*(p[2] / p[0])**(1.0/gamma) + rho[3] = rho[4] * (p[3] + Gamma*p[4]) / (p[4] + Gamma*p[3]) + + # calc x + x = np.array([0.0, 0.5 - args.t*u[1], 0.5, 0.5 + args.t*u[3], 0.5 + args.t*u[4], 1,0]) + + # calc rho + rho1 = np.ones(N) * rho[0] + rho2 = np.linspace(rho[0], rho[2], N) + rho3 = np.ones(N) * rho[2] + rho4 = np.ones(N) * rho[3] + rho5 = np.ones(N) * rho[4] + + + plt.plot(np.linspace(x[0], x[1], N), rho1) + plt.plot(np.linspace(x[1], x[2], N), rho2) + plt.plot(np.linspace(x[2], x[3], N), rho3) + plt.plot(np.linspace(x[3], x[4], N), rho4) + plt.plot(np.linspace(x[4], x[5], N), rho5) + + plt.show() + + print("done") \ No newline at end of file From 5c58fec859c3e23c9a11720bf2344201d3b11d62 Mon Sep 17 00:00:00 2001 From: Basov Date: Thu, 29 Sep 2022 13:12:44 +0200 Subject: [PATCH 16/33] added plotting to script --- tools/sod_analytical.py | 27 ++++++++++++--------------- 1 file changed, 12 insertions(+), 15 deletions(-) diff --git a/tools/sod_analytical.py b/tools/sod_analytical.py index 00dbe61..c389a99 100644 --- a/tools/sod_analytical.py +++ b/tools/sod_analytical.py @@ -34,6 +34,15 @@ def _calc_p3(p, rho, Gamma, gamma, beta): return p3 +def plot_rho(x, rho): + plt.plot((x[0], x[1]), (rho[0], rho[1])) + plt.plot((x[1], x[2]), (rho[1], rho[2])) + plt.plot((x[2], x[3]), (rho[2], rho[2])) + plt.plot((x[3], x[3]), (rho[2], rho[3])) + plt.plot((x[3], x[4]), (rho[3], rho[3])) + plt.plot((x[4], x[4]), (rho[3], rho[4])) + plt.plot((x[4], x[5]), (rho[4], rho[4])) + if __name__ == '__main__': parser = argparse.ArgumentParser() parser.add_argument('t', type=float, help='time of results') @@ -71,25 +80,13 @@ def _calc_p3(p, rho, Gamma, gamma, beta): u[3] = (p[2] - p[4]) / np.sqrt((rho[4]/2) * ((gamma + 1) * p[2] + (gamma - 1) * p[4])) u[4] = np.sqrt(gamma * p[-1] / rho[-1]) # speed of pressure increase going right + rho[1] = rho[0] rho[2] = rho[0]*(p[2] / p[0])**(1.0/gamma) rho[3] = rho[4] * (p[3] + Gamma*p[4]) / (p[4] + Gamma*p[3]) # calc x - x = np.array([0.0, 0.5 - args.t*u[1], 0.5, 0.5 + args.t*u[3], 0.5 + args.t*u[4], 1,0]) - - # calc rho - rho1 = np.ones(N) * rho[0] - rho2 = np.linspace(rho[0], rho[2], N) - rho3 = np.ones(N) * rho[2] - rho4 = np.ones(N) * rho[3] - rho5 = np.ones(N) * rho[4] - - - plt.plot(np.linspace(x[0], x[1], N), rho1) - plt.plot(np.linspace(x[1], x[2], N), rho2) - plt.plot(np.linspace(x[2], x[3], N), rho3) - plt.plot(np.linspace(x[3], x[4], N), rho4) - plt.plot(np.linspace(x[4], x[5], N), rho5) + x = np.array([0.0, 0.5 - args.t*u[1], 0.5, 0.5 + args.t*u[3], 0.5 + args.t*u[3] + args.t*u[4], 1,0]) + plot_rho(x, rho) plt.show() From 77fc4ac0d6c20672090b9d94a1a63b6ad4fb3b07 Mon Sep 17 00:00:00 2001 From: Basov Date: Thu, 29 Sep 2022 13:15:53 +0200 Subject: [PATCH 17/33] added plot_t --- examples/plot_n.py | 104 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 104 insertions(+) create mode 100644 examples/plot_n.py diff --git a/examples/plot_n.py b/examples/plot_n.py new file mode 100644 index 0000000..6deaf4d --- /dev/null +++ b/examples/plot_n.py @@ -0,0 +1,104 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu Sep 29 13:13:54 2022 + +@author: baso_le +""" + +import matplotlib.pyplot as plt +import csv +import numpy as np +import argparse + +def _calc_u3(p1, p3, rhoL, Gamma, gamma, beta): + num = (1 - Gamma**2) * p1**(1.0/gamma) + denum = Gamma**2 * rhoL + return (p1**beta - p3**beta) * np.sqrt(num / denum) + +def _calc_u4(p3, p5, rhoR, Gamma): + return (p3 - p5) * np.sqrt((1 - Gamma) / (rhoR*(p3 - Gamma*p5))) + +def _calc_p3(p, rho, Gamma, gamma, beta): + p1 = p[0] + p3 = p[0] + p5 = p[4] + dp = 1e-5 + rhoL = rho[0] + rhoR = rho[-1] + err = 1e+5 + err_last = 1e+6 + + while err < err_last: + err_last = err + err = abs(_calc_u3(p1, p3, rhoL, Gamma, gamma, beta) - _calc_u4(p3, p5, rhoR, Gamma)) + p3 -= dp + + return p3 + +def plot_rho(x, rho): + plt.plot((x[0], x[1]), (rho[0], rho[1])) + plt.plot((x[1], x[2]), (rho[1], rho[2])) + plt.plot((x[2], x[3]), (rho[2], rho[2])) + plt.plot((x[3], x[3]), (rho[2], rho[3])) + plt.plot((x[3], x[4]), (rho[3], rho[3])) + plt.plot((x[4], x[4]), (rho[3], rho[4])) + plt.plot((x[4], x[5]), (rho[4], rho[4])) + +if __name__ == '__main__': + parser = argparse.ArgumentParser(description='Process some integers.') + parser.add_argument('file_name', type=str) + parser.add_argument('t', type=float) + + args = parser.parse_args() + + with open(args.file_name) as file: + reader = csv.reader(file, delimiter=',') + + for line in reader: + l = [m for m in line if m] + n = [float(l[i]) for i in range(len(l))] + + plt.plot(n) + + # number points per segment + N = 10 + + # gas = Ar + gamma = 1.67 + Gamma = (gamma - 1.0) / (gamma + 1.0) + beta = (gamma - 1.0) / (2.0 * gamma) + + # boundary conditions + rho = np.zeros(5) + p = np.zeros(5) + u = np.zeros(5) + + rho[0] = 1.0 + p[0] = 1.0 + u[0] = 0.0 + + rho[-1] = 0.125 + p[-1] = 0.1 + u[-1] = 0.0 + + # calculating states + p[2] = _calc_p3(p, rho, Gamma, gamma, beta) + p[3] = p[2] + + # speed of characterisics + u[1] = np.sqrt(gamma * p[0] / rho[0]) # speed of rarefication going left + u[2] = 0 + u[3] = (p[2] - p[4]) / np.sqrt((rho[4]/2) * ((gamma + 1) * p[2] + (gamma - 1) * p[4])) + u[4] = np.sqrt(gamma * p[-1] / rho[-1]) # speed of pressure increase going right + + rho[1] = rho[0] + rho[2] = rho[0]*(p[2] / p[0])**(1.0/gamma) + rho[3] = rho[4] * (p[3] + Gamma*p[4]) / (p[4] + Gamma*p[3]) + + # calc x + x = np.array([0.0, 0.5 - args.t*u[1], 0.5, 0.5 + args.t*u[3], 0.5 + args.t*u[3] + args.t*u[4], 1,0]) + plot_rho(x, rho) + + plt.show() + + print("done") From 68ee0ece64a64c66737906954593b1079ef17f1d Mon Sep 17 00:00:00 2001 From: Basov Date: Thu, 29 Sep 2022 14:27:33 +0200 Subject: [PATCH 18/33] extended tools --- tools/sod_analytical.py | 35 ++++++++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/tools/sod_analytical.py b/tools/sod_analytical.py index c389a99..01d9209 100644 --- a/tools/sod_analytical.py +++ b/tools/sod_analytical.py @@ -21,7 +21,7 @@ def _calc_p3(p, rho, Gamma, gamma, beta): p1 = p[0] p3 = p[0] p5 = p[4] - dp = 1e-5 + dp = p1 * 1e-3 rhoL = rho[0] rhoR = rho[-1] err = 1e+5 @@ -34,7 +34,7 @@ def _calc_p3(p, rho, Gamma, gamma, beta): return p3 -def plot_rho(x, rho): +def plot_val(x, val, name): plt.plot((x[0], x[1]), (rho[0], rho[1])) plt.plot((x[1], x[2]), (rho[1], rho[2])) plt.plot((x[2], x[3]), (rho[2], rho[2])) @@ -42,13 +42,31 @@ def plot_rho(x, rho): plt.plot((x[3], x[4]), (rho[3], rho[3])) plt.plot((x[4], x[4]), (rho[3], rho[4])) plt.plot((x[4], x[5]), (rho[4], rho[4])) + + plt.ylabel(name) + plt.xlabel("x") + + plt.show() + +def check_args(args): + if args.p is not None and (args.p[0] <= args.p[1]): + raise Exception("p1 must be > p2") + elif args.rho is not None and (args.rho[0] <= args.rho[1]): + raise Exception("rho1 must be > rho2") + elif args.L is not None and (args.L <= 0.0): + raise Exception("L must be > 0") if __name__ == '__main__': parser = argparse.ArgumentParser() parser.add_argument('t', type=float, help='time of results') + parser.add_argument('-p', type=float, help='pressure', nargs = 2, default=(1.0, 0.1)) + parser.add_argument('-rho', type=float, help='density', nargs = 2, default=(1.0, 0.125)) + parser.add_argument('-L', type=float, help='tube length', default=1.0) args = parser.parse_args() + check_args(args) + # number points per segment N = 10 @@ -62,12 +80,12 @@ def plot_rho(x, rho): p = np.zeros(5) u = np.zeros(5) - rho[0] = 1.0 - p[0] = 1.0 + rho[0] = args.rho[0] + p[0] = args.p[0] u[0] = 0.0 - rho[-1] = 0.125 - p[-1] = 0.1 + rho[-1] = args.rho[1] + p[-1] = args.p[1] u[-1] = 0.0 # calculating states @@ -85,9 +103,8 @@ def plot_rho(x, rho): rho[3] = rho[4] * (p[3] + Gamma*p[4]) / (p[4] + Gamma*p[3]) # calc x - x = np.array([0.0, 0.5 - args.t*u[1], 0.5, 0.5 + args.t*u[3], 0.5 + args.t*u[3] + args.t*u[4], 1,0]) - plot_rho(x, rho) + x = np.array([0.0, args.L*0.5 - args.t*u[1], args.L*0.5, args.L*0.5 + args.t*u[3], args.L*0.5 + args.t*u[3] + args.t*u[4], args.L]) - plt.show() + plot_val(x, rho, "rho") print("done") \ No newline at end of file From 2f66fe9a6e9652a942177e9e8c54f27e19b3f257 Mon Sep 17 00:00:00 2001 From: Basov Date: Thu, 29 Sep 2022 14:30:24 +0200 Subject: [PATCH 19/33] fixed bug in tools --- tools/sod_analytical.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/tools/sod_analytical.py b/tools/sod_analytical.py index 01d9209..27c9061 100644 --- a/tools/sod_analytical.py +++ b/tools/sod_analytical.py @@ -35,13 +35,13 @@ def _calc_p3(p, rho, Gamma, gamma, beta): return p3 def plot_val(x, val, name): - plt.plot((x[0], x[1]), (rho[0], rho[1])) - plt.plot((x[1], x[2]), (rho[1], rho[2])) - plt.plot((x[2], x[3]), (rho[2], rho[2])) - plt.plot((x[3], x[3]), (rho[2], rho[3])) - plt.plot((x[3], x[4]), (rho[3], rho[3])) - plt.plot((x[4], x[4]), (rho[3], rho[4])) - plt.plot((x[4], x[5]), (rho[4], rho[4])) + plt.plot((x[0], x[1]), (val[0], val[1])) + plt.plot((x[1], x[2]), (val[1], val[2])) + plt.plot((x[2], x[3]), (val[2], val[2])) + plt.plot((x[3], x[3]), (val[2], val[3])) + plt.plot((x[3], x[4]), (val[3], val[3])) + plt.plot((x[4], x[4]), (val[3], val[4])) + plt.plot((x[4], x[5]), (val[4], val[4])) plt.ylabel(name) plt.xlabel("x") @@ -62,6 +62,7 @@ def check_args(args): parser.add_argument('-p', type=float, help='pressure', nargs = 2, default=(1.0, 0.1)) parser.add_argument('-rho', type=float, help='density', nargs = 2, default=(1.0, 0.125)) parser.add_argument('-L', type=float, help='tube length', default=1.0) + parser.add_argument('-plt', type=bool, help='plot values', default=True) args = parser.parse_args() @@ -89,6 +90,7 @@ def check_args(args): u[-1] = 0.0 # calculating states + p[1] = p[0] p[2] = _calc_p3(p, rho, Gamma, gamma, beta) p[3] = p[2] @@ -105,6 +107,8 @@ def check_args(args): # calc x x = np.array([0.0, args.L*0.5 - args.t*u[1], args.L*0.5, args.L*0.5 + args.t*u[3], args.L*0.5 + args.t*u[3] + args.t*u[4], args.L]) - plot_val(x, rho, "rho") + if args.plt: + plot_val(x, rho, "rho") + plot_val(x, p, "p") print("done") \ No newline at end of file From 9d07a1eb54a7783c8e8b4a02b86936334f97630b Mon Sep 17 00:00:00 2001 From: Basov Date: Thu, 29 Sep 2022 14:40:05 +0200 Subject: [PATCH 20/33] implemented wrtiting to file --- tools/sod_analytical.py | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/tools/sod_analytical.py b/tools/sod_analytical.py index 27c9061..792e22f 100644 --- a/tools/sod_analytical.py +++ b/tools/sod_analytical.py @@ -48,6 +48,16 @@ def plot_val(x, val, name): plt.show() +def write_values(x, val, name, file): + with open(file + "_" + name + ".csv", "w") as file: + file.write("{}, {}\n".format(x[0], val[0])) + file.write("{}, {}\n".format(x[1], val[1])) + file.write("{}, {}\n".format(x[2], val[2])) + file.write("{}, {}\n".format(x[3], val[2])) + file.write("{}, {}\n".format(x[3], val[3])) + file.write("{}, {}\n".format(x[4], val[3])) + file.write("{}, {}\n".format(x[4], val[4])) + def check_args(args): if args.p is not None and (args.p[0] <= args.p[1]): raise Exception("p1 must be > p2") @@ -62,7 +72,8 @@ def check_args(args): parser.add_argument('-p', type=float, help='pressure', nargs = 2, default=(1.0, 0.1)) parser.add_argument('-rho', type=float, help='density', nargs = 2, default=(1.0, 0.125)) parser.add_argument('-L', type=float, help='tube length', default=1.0) - parser.add_argument('-plt', type=bool, help='plot values', default=True) + parser.add_argument('-plt', type=bool, help='plot values', const=True, nargs='?') + parser.add_argument('-w', type=str, help='write to file') args = parser.parse_args() @@ -110,5 +121,10 @@ def check_args(args): if args.plt: plot_val(x, rho, "rho") plot_val(x, p, "p") + + if args.w: + print("writing to file " + args.w + "_X.csv") + write_values(x, rho, "rho", args.w) + write_values(x, p, "p", args.w) print("done") \ No newline at end of file From 7ce1e8eb175071b01eb258f07d396766e4770ef0 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Thu, 29 Sep 2022 14:42:04 +0200 Subject: [PATCH 21/33] removed unecesary file --- examples/plot_n.py | 104 --------------------------------------------- 1 file changed, 104 deletions(-) delete mode 100644 examples/plot_n.py diff --git a/examples/plot_n.py b/examples/plot_n.py deleted file mode 100644 index 6deaf4d..0000000 --- a/examples/plot_n.py +++ /dev/null @@ -1,104 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Thu Sep 29 13:13:54 2022 - -@author: baso_le -""" - -import matplotlib.pyplot as plt -import csv -import numpy as np -import argparse - -def _calc_u3(p1, p3, rhoL, Gamma, gamma, beta): - num = (1 - Gamma**2) * p1**(1.0/gamma) - denum = Gamma**2 * rhoL - return (p1**beta - p3**beta) * np.sqrt(num / denum) - -def _calc_u4(p3, p5, rhoR, Gamma): - return (p3 - p5) * np.sqrt((1 - Gamma) / (rhoR*(p3 - Gamma*p5))) - -def _calc_p3(p, rho, Gamma, gamma, beta): - p1 = p[0] - p3 = p[0] - p5 = p[4] - dp = 1e-5 - rhoL = rho[0] - rhoR = rho[-1] - err = 1e+5 - err_last = 1e+6 - - while err < err_last: - err_last = err - err = abs(_calc_u3(p1, p3, rhoL, Gamma, gamma, beta) - _calc_u4(p3, p5, rhoR, Gamma)) - p3 -= dp - - return p3 - -def plot_rho(x, rho): - plt.plot((x[0], x[1]), (rho[0], rho[1])) - plt.plot((x[1], x[2]), (rho[1], rho[2])) - plt.plot((x[2], x[3]), (rho[2], rho[2])) - plt.plot((x[3], x[3]), (rho[2], rho[3])) - plt.plot((x[3], x[4]), (rho[3], rho[3])) - plt.plot((x[4], x[4]), (rho[3], rho[4])) - plt.plot((x[4], x[5]), (rho[4], rho[4])) - -if __name__ == '__main__': - parser = argparse.ArgumentParser(description='Process some integers.') - parser.add_argument('file_name', type=str) - parser.add_argument('t', type=float) - - args = parser.parse_args() - - with open(args.file_name) as file: - reader = csv.reader(file, delimiter=',') - - for line in reader: - l = [m for m in line if m] - n = [float(l[i]) for i in range(len(l))] - - plt.plot(n) - - # number points per segment - N = 10 - - # gas = Ar - gamma = 1.67 - Gamma = (gamma - 1.0) / (gamma + 1.0) - beta = (gamma - 1.0) / (2.0 * gamma) - - # boundary conditions - rho = np.zeros(5) - p = np.zeros(5) - u = np.zeros(5) - - rho[0] = 1.0 - p[0] = 1.0 - u[0] = 0.0 - - rho[-1] = 0.125 - p[-1] = 0.1 - u[-1] = 0.0 - - # calculating states - p[2] = _calc_p3(p, rho, Gamma, gamma, beta) - p[3] = p[2] - - # speed of characterisics - u[1] = np.sqrt(gamma * p[0] / rho[0]) # speed of rarefication going left - u[2] = 0 - u[3] = (p[2] - p[4]) / np.sqrt((rho[4]/2) * ((gamma + 1) * p[2] + (gamma - 1) * p[4])) - u[4] = np.sqrt(gamma * p[-1] / rho[-1]) # speed of pressure increase going right - - rho[1] = rho[0] - rho[2] = rho[0]*(p[2] / p[0])**(1.0/gamma) - rho[3] = rho[4] * (p[3] + Gamma*p[4]) / (p[4] + Gamma*p[3]) - - # calc x - x = np.array([0.0, 0.5 - args.t*u[1], 0.5, 0.5 + args.t*u[3], 0.5 + args.t*u[3] + args.t*u[4], 1,0]) - plot_rho(x, rho) - - plt.show() - - print("done") From b5368a9e516d7b963c6a4c43d419ae5c0f458fe7 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Thu, 29 Sep 2022 14:44:36 +0200 Subject: [PATCH 22/33] modifed shocktube input --- examples/shock_tube.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/shock_tube.py b/examples/shock_tube.py index 23e468a..80d9583 100644 --- a/examples/shock_tube.py +++ b/examples/shock_tube.py @@ -43,7 +43,7 @@ def write2file(solver, Nbins): Boxhigh = [(-0.0001, 0.0001), (-0.0001, 0.0001), (0.0, 0.05)] # low denisty particles - nlow = 2.41432e+21 + nlow = nhigh*0.125 Tlow = 300 Boxlow = [(-0.0001, 0.0001), (-0.0001, 0.0001), (0.05, 0.1)] From ea08607f18f44ff130a59c7853566304c2bf7e72 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Thu, 29 Sep 2022 15:08:13 +0200 Subject: [PATCH 23/33] updated sod tool --- tools/sod_analytical.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tools/sod_analytical.py b/tools/sod_analytical.py index 792e22f..162c147 100644 --- a/tools/sod_analytical.py +++ b/tools/sod_analytical.py @@ -57,6 +57,7 @@ def write_values(x, val, name, file): file.write("{}, {}\n".format(x[3], val[3])) file.write("{}, {}\n".format(x[4], val[3])) file.write("{}, {}\n".format(x[4], val[4])) + file.write("{}, {}\n".format(x[5], val[4])) def check_args(args): if args.p is not None and (args.p[0] <= args.p[1]): @@ -127,4 +128,4 @@ def check_args(args): write_values(x, rho, "rho", args.w) write_values(x, p, "p", args.w) - print("done") \ No newline at end of file + print("done") From 94a3456cb293536b9d5e53103e4eca13f7ecab26 Mon Sep 17 00:00:00 2001 From: Basov Date: Thu, 29 Sep 2022 16:39:58 +0200 Subject: [PATCH 24/33] extended tools --- tools/sod_analytical.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/tools/sod_analytical.py b/tools/sod_analytical.py index 162c147..d9660df 100644 --- a/tools/sod_analytical.py +++ b/tools/sod_analytical.py @@ -9,6 +9,8 @@ import matplotlib.pyplot as plt import argparse +kb = 1.380649e-23 + def _calc_u3(p1, p3, rhoL, Gamma, gamma, beta): num = (1 - Gamma**2) * p1**(1.0/gamma) denum = Gamma**2 * rhoL @@ -69,10 +71,10 @@ def check_args(args): if __name__ == '__main__': parser = argparse.ArgumentParser() - parser.add_argument('t', type=float, help='time of results') - parser.add_argument('-p', type=float, help='pressure', nargs = 2, default=(1.0, 0.1)) - parser.add_argument('-rho', type=float, help='density', nargs = 2, default=(1.0, 0.125)) - parser.add_argument('-L', type=float, help='tube length', default=1.0) + parser.add_argument('-t', type=float, help='time of results', default=3.0e-5) + parser.add_argument('-p', type=float, help='pressure', nargs = 2, default=(100, 10.0)) + parser.add_argument('-rho', type=float, help='density', nargs = 2, default=(0.0016036396304, 0.0016036396304*0.1)) + parser.add_argument('-L', type=float, help='tube length', default=0.1) parser.add_argument('-plt', type=bool, help='plot values', const=True, nargs='?') parser.add_argument('-w', type=str, help='write to file') @@ -87,6 +89,7 @@ def check_args(args): gamma = 1.67 Gamma = (gamma - 1.0) / (gamma + 1.0) beta = (gamma - 1.0) / (2.0 * gamma) + mass = 6.6422e-26 # boundary conditions rho = np.zeros(5) @@ -116,12 +119,17 @@ def check_args(args): rho[2] = rho[0]*(p[2] / p[0])**(1.0/gamma) rho[3] = rho[4] * (p[3] + Gamma*p[4]) / (p[4] + Gamma*p[3]) + n = [r/mass for r in rho] + T = [p[i] / (n[i] * kb) for i in range(5)] + # calc x x = np.array([0.0, args.L*0.5 - args.t*u[1], args.L*0.5, args.L*0.5 + args.t*u[3], args.L*0.5 + args.t*u[3] + args.t*u[4], args.L]) if args.plt: plot_val(x, rho, "rho") + plot_val(x, n, "n") plot_val(x, p, "p") + plot_val(x, T, "T") if args.w: print("writing to file " + args.w + "_X.csv") From c3924701eef78865670056b8b74c1b8a6e1a24ac Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Thu, 29 Sep 2022 17:18:19 +0200 Subject: [PATCH 25/33] fixed bug in script --- tools/sod_analytical.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tools/sod_analytical.py b/tools/sod_analytical.py index d9660df..9e0c1de 100644 --- a/tools/sod_analytical.py +++ b/tools/sod_analytical.py @@ -135,5 +135,7 @@ def check_args(args): print("writing to file " + args.w + "_X.csv") write_values(x, rho, "rho", args.w) write_values(x, p, "p", args.w) + write_values(x, n, "n", args.w) + write_values(x, T, "T", args.w) print("done") From f85278064cd2b7b035ffa53ce69f3e116ee04065 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Thu, 29 Sep 2022 17:19:00 +0200 Subject: [PATCH 26/33] extended plot script --- examples/plot.py | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/examples/plot.py b/examples/plot.py index 18ac59c..e0e2afe 100644 --- a/examples/plot.py +++ b/examples/plot.py @@ -6,6 +6,7 @@ if __name__ == '__main__': parser = argparse.ArgumentParser(description='Process some integers.') parser.add_argument('file_name', type=str) + parser.add_argument('-ref', type=str) args = parser.parse_args() @@ -15,6 +16,19 @@ for line in reader: l = [m for m in line if m] n = [float(l[i]) for i in range(len(l))] - - plt.plot(n) - plt.show() + + plt.plot(np.linspace(0, 0.1, len(n)), n) + + if args.ref: + x = [] + val = [] + with open(args.ref) as file: + reader = csv.reader(file, delimiter=',') + + for line in reader: + x.append(float(line[0])) + val.append(float(line[1])) + + plt.plot(x, val) + + plt.show() From 9912974aca2eff4cf97347ac4d4efa94b7cb2773 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Thu, 29 Sep 2022 17:19:40 +0200 Subject: [PATCH 27/33] updated shock tube test case --- examples/shock_tube.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/shock_tube.py b/examples/shock_tube.py index 80d9583..37d60c0 100644 --- a/examples/shock_tube.py +++ b/examples/shock_tube.py @@ -43,7 +43,7 @@ def write2file(solver, Nbins): Boxhigh = [(-0.0001, 0.0001), (-0.0001, 0.0001), (0.0, 0.05)] # low denisty particles - nlow = nhigh*0.125 + nlow = nhigh*0.1 Tlow = 300 Boxlow = [(-0.0001, 0.0001), (-0.0001, 0.0001), (0.05, 0.1)] From 1d687f397a53b1f44fff3513c641c6e50f15b759 Mon Sep 17 00:00:00 2001 From: LeoBasov <40352679+LeoBasov@users.noreply.github.com> Date: Fri, 30 Sep 2022 21:58:59 +0200 Subject: [PATCH 28/33] added common module --- dsmc/common.py | 5 +++++ tests/unit/main.py | 1 + tests/unit/test_dsmc/common.py | 8 ++++++++ 3 files changed, 14 insertions(+) create mode 100644 dsmc/common.py create mode 100644 tests/unit/test_dsmc/common.py diff --git a/dsmc/common.py b/dsmc/common.py new file mode 100644 index 0000000..018d46e --- /dev/null +++ b/dsmc/common.py @@ -0,0 +1,5 @@ +from numba import njit + +@njit(float([[:] [:], [:]])) +def get_V(box): + return (box[0][1] - box[0][0]) * (box[1][1] - box[1][0]) * (box[2][1] - box[2][0]) \ No newline at end of file diff --git a/tests/unit/main.py b/tests/unit/main.py index f9c8d93..df9783a 100644 --- a/tests/unit/main.py +++ b/tests/unit/main.py @@ -6,6 +6,7 @@ from test_dsmc.particles import * from test_dsmc.octree import * from test_dsmc.diagnostics import * +from test_dsmc.common import * if __name__ == '__main__': unittest.main() diff --git a/tests/unit/test_dsmc/common.py b/tests/unit/test_dsmc/common.py new file mode 100644 index 0000000..bdf2eb7 --- /dev/null +++ b/tests/unit/test_dsmc/common.py @@ -0,0 +1,8 @@ +import numpy as np +import dsmc.common as com +import unittest + +class TestCommon(unittest.TestCase): + def test_pass(self): + box = np.array([[-1.0, 1.0], [-1.0, 1.0], [-1.0, 1.0]]) + V= com.get_V(box) \ No newline at end of file From 16f8dde62d05747642da387bae647f5679803034 Mon Sep 17 00:00:00 2001 From: LeoBasov <40352679+LeoBasov@users.noreply.github.com> Date: Fri, 30 Sep 2022 23:10:58 +0200 Subject: [PATCH 29/33] fixed new get_V function --- dsmc/common.py | 8 ++++---- tests/unit/test_dsmc/common.py | 3 ++- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/dsmc/common.py b/dsmc/common.py index 018d46e..0990d48 100644 --- a/dsmc/common.py +++ b/dsmc/common.py @@ -1,5 +1,5 @@ -from numba import njit +from numba import cfunc -@njit(float([[:] [:], [:]])) -def get_V(box): - return (box[0][1] - box[0][0]) * (box[1][1] - box[1][0]) * (box[2][1] - box[2][0]) \ No newline at end of file +@cfunc("double(double[::3, ::2])") +def get_V(a): + return (a[0][1] - a[0][0]) * (a[1][1] - a[1][0]) * (a[2][1] - a[2][0]) \ No newline at end of file diff --git a/tests/unit/test_dsmc/common.py b/tests/unit/test_dsmc/common.py index bdf2eb7..b6e73ec 100644 --- a/tests/unit/test_dsmc/common.py +++ b/tests/unit/test_dsmc/common.py @@ -5,4 +5,5 @@ class TestCommon(unittest.TestCase): def test_pass(self): box = np.array([[-1.0, 1.0], [-1.0, 1.0], [-1.0, 1.0]]) - V= com.get_V(box) \ No newline at end of file + V = com.get_V(box) + self.assertEqual(8.0, V) \ No newline at end of file From 5e1a40dc73d96f7f8b818606fb8375ee5313959b Mon Sep 17 00:00:00 2001 From: LeoBasov <40352679+LeoBasov@users.noreply.github.com> Date: Fri, 30 Sep 2022 23:20:43 +0200 Subject: [PATCH 30/33] moved function get_V to common --- dsmc/dsmc.py | 5 +++-- dsmc/octree.py | 4 ---- dsmc/particles.py | 4 ++-- 3 files changed, 5 insertions(+), 8 deletions(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index bfc404f..4b7d212 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -3,6 +3,7 @@ from numba import prange from . import particles as prt from . import octree as oc +from . import common as com @njit def _push(velocities, positions, dt): @@ -130,7 +131,7 @@ def _update_velocities(permutations : np.ndarray, velocities : np.ndarray, mass def _update_vels(permutations : np.ndarray, velocities : np.ndarray, mass : float, sigma_T : float, dt : float, w : float, elem_offsets : np.ndarray, number_elements : np.ndarray, number_children : np.ndarray, cell_boxes : np.ndarray, Nleafs : int) -> np.ndarray: for i in range(Nleafs): if not number_children[i] and number_elements[i]: - Vc = oc.get_V(cell_boxes[i]) + Vc = com.get_V(cell_boxes[i]) velocities = _update_velocities(permutations, velocities, mass, sigma_T, Vc, dt, w, elem_offsets[i], number_elements[i]) return velocities @@ -212,7 +213,7 @@ def _update_velocities(self, dt): def create_particles(self, box, T, n, u = np.zeros(3)): box = np.array(box) - N = int(round(oc.get_V(box) * n / self.w)) + 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) diff --git a/dsmc/octree.py b/dsmc/octree.py index 18afcd4..0dc634b 100644 --- a/dsmc/octree.py +++ b/dsmc/octree.py @@ -110,10 +110,6 @@ def _sort(permutations : npt.NDArray, box : npt.NDArray, positions : npt.NDArray return new_permutations, Nnew -@njit -def get_V(box): - return (box[0][1] - box[0][0]) * (box[1][1] - box[1][0]) * (box[2][1] - box[2][0]) - @njit def _create_boxes(box): half = np.array([0.5*(box[i][0] + box[i][1]) for i in range(3)]) diff --git a/dsmc/particles.py b/dsmc/particles.py index f018b92..421fbff 100644 --- a/dsmc/particles.py +++ b/dsmc/particles.py @@ -1,7 +1,7 @@ import math import numpy as np from numba import njit -from . import octree as oc +from . import common as com kb = 1.380649e-23 @@ -155,6 +155,6 @@ def inflow(self, mass, T, u, n, w, dt, domain, axis, minmax): box[axis][0] = box[axis][1] box[axis][1] = box[axis][0] + L - N = int(round(oc.get_V(box) * n / w)) + N = int(round(com.get_V(box) * n / w)) self.create_particles(box, mass, T, N, u) From 08719a932933d7ffb3a1f9a2bcd3166ecf3483db Mon Sep 17 00:00:00 2001 From: LeoBasov <40352679+LeoBasov@users.noreply.github.com> Date: Fri, 30 Sep 2022 23:21:44 +0200 Subject: [PATCH 31/33] removed unused package --- dsmc/octree.py | 1 - 1 file changed, 1 deletion(-) diff --git a/dsmc/octree.py b/dsmc/octree.py index 0dc634b..cd18815 100644 --- a/dsmc/octree.py +++ b/dsmc/octree.py @@ -1,7 +1,6 @@ import numpy as np import numpy.typing as npt from numba import njit -import numba as nb fmin = np.finfo(float).min fmax = np.finfo(float).max From fccc30bb9d4d4baa04fd83d565318dadd82a71c9 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Sat, 1 Oct 2022 10:01:16 +0200 Subject: [PATCH 32/33] updated unit tests --- tests/unit/test_dsmc/octree.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/tests/unit/test_dsmc/octree.py b/tests/unit/test_dsmc/octree.py index f04fd22..0fb249d 100644 --- a/tests/unit/test_dsmc/octree.py +++ b/tests/unit/test_dsmc/octree.py @@ -1,7 +1,7 @@ import unittest import numpy as np from dsmc import octree as oc -from dsmc import particles as part +from dsmc import common as com import csv class TestOctree(unittest.TestCase): @@ -217,9 +217,9 @@ def test__create_boxes(self): V = 0.0 for box in boxes: - V += oc.get_V(box) + V += com.get_V(box) - self.assertEqual(oc.get_V(box_orig), V) + self.assertEqual(com.get_V(box_orig), V) def test__get_min_aspect_ratio_1(self): box = np.array([(0.0, 1.0), (0.0, 10.0), (0.0, 100.0)]) @@ -275,9 +275,9 @@ def test__create_combined_boxes_1(self): V = 0.0 for b in boxes_new: - V += oc.get_V(b) + V += com.get_V(b) - self.assertEqual(oc.get_V(box), V) + self.assertEqual(com.get_V(box), V) self.assertEqual(len(boxes_old), len(boxes_new)) def test__create_combined_boxes_2(self): @@ -287,9 +287,9 @@ def test__create_combined_boxes_2(self): V = 0.0 for b in boxes_new: - V += oc.get_V(b) + V += com.get_V(b) - self.assertEqual(oc.get_V(box), V) + self.assertEqual(com.get_V(box), V) self.assertEqual(4, len(boxes_new)) def test__create_combined_boxes_3(self): @@ -299,9 +299,9 @@ def test__create_combined_boxes_3(self): V = 0.0 for b in boxes_new: - V += oc.get_V(b) + V += com.get_V(b) - self.assertEqual(oc.get_V(box), V) + self.assertEqual(com.get_V(box), V) self.assertEqual(4, len(boxes_new)) def test__create_combined_boxes_4(self): @@ -311,9 +311,9 @@ def test__create_combined_boxes_4(self): V = 0.0 for b in boxes_new: - V += oc.get_V(b) + V += com.get_V(b) - self.assertEqual(oc.get_V(box), V) + self.assertEqual(com.get_V(box), V) self.assertEqual(8, len(boxes_new)) class TestOctreeOctree(unittest.TestCase): From 9d2e54cd36b162da1c311e61453cfd7df44db349 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Sat, 1 Oct 2022 10:08:39 +0200 Subject: [PATCH 33/33] updated changelog and setup.py --- CHANGELOG.md | 5 +++++ setup.py | 4 ++-- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c85a38a..5a0257d 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.7.0] - 2022-10-01 +### Added +- inflow boundary condition +- open boundary condition + ## [0.6.0] - 2022-09-28 ### Added - functionality for aspect ratio check in octree diff --git a/setup.py b/setup.py index 689ad62..9a176ef 100644 --- a/setup.py +++ b/setup.py @@ -1,9 +1,9 @@ from setuptools import setup setup( name='dsmc', - version='0.6.0', + version='0.7.0', author='Leo Basov', - python_requires='>=3.6, <4', + python_requires='>=3.10, <4', packages=["dsmc"], install_requires=[ 'numpy',