From 0f36ab0bc0ee26df99300768863497b127ac5586 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Thu, 22 Sep 2022 15:04:09 +0200 Subject: [PATCH 01/19] extended create_particles --- dsmc/particles.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/dsmc/particles.py b/dsmc/particles.py index 852ff0c..31f17fa 100644 --- a/dsmc/particles.py +++ b/dsmc/particles.py @@ -127,8 +127,13 @@ def VelPos(self, vel_pos): self._velocities = vel_pos[0] self._positions = vel_pos[1] self._N = len(self._positions) - + def create_particles(self, X, mass, T, N): - self._velocities = get_velocities(T, mass, N) - self._positions = calc_positions(X[0], X[1], X[2], N) - self._N = N + if self._N == 0: + self._velocities = get_velocities(T, mass, N) + 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._positions = np.concatenate(self._positions, calc_positions(X[0], X[1], X[2], N)) + self._N += N From 6e2093b36af97894908d05af73a72fcbc4e968f2 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Thu, 22 Sep 2022 15:40:40 +0200 Subject: [PATCH 02/19] started solver implementation --- dsmc/__init__.py | 1 + dsmc/dsmc.py | 38 ++++++++++++++++++++++++++++++++---- dsmc/particles.py | 2 +- examples/shock_tube.py | 32 ++++++++++++++++++++++++++++++ tests/unit/test_dsmc/dsmc.py | 6 +++--- 5 files changed, 71 insertions(+), 8 deletions(-) create mode 100644 examples/shock_tube.py diff --git a/dsmc/__init__.py b/dsmc/__init__.py index e69de29..9ebec8a 100644 --- a/dsmc/__init__.py +++ b/dsmc/__init__.py @@ -0,0 +1 @@ +from .dsmc import DSMC diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index b30dac7..b587c9a 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -1,7 +1,37 @@ import math -import numpy +import numpy as np from numba import njit +from . import particles as prt +from . import octree as oc -@njit -def test(): - return 1.0 +class DSMC: + def __init__(self): + self.clear() + + def clear(self): + self.particles = prt.Particles() + self.octree = oc.Octree() + self.w = None + self.domain = None + + def advance(self, dt): + if self.domain == None: + raise Exception("simulation domain not defined") + if particles.N == 0: + raise Exception("no particles created") + if self.w == None: + raise Exception("particle weight not set") + + def create_particles(self, box, mass, T, n): + N = int(round(n / self.w)) + print("creating {} particles".format(N)) + self.particles.create_particles(box, mass, T, N) + + print("now containing {} particles, {} total".format(N, self.particles.N)) + + def set_domain(self, domain): + self.domain = np.array(domain) + + def set_weight(self, w): + self.octree.w = w + self.w = w diff --git a/dsmc/particles.py b/dsmc/particles.py index 31f17fa..a978c34 100644 --- a/dsmc/particles.py +++ b/dsmc/particles.py @@ -128,7 +128,7 @@ 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): if self._N == 0: self._velocities = get_velocities(T, mass, N) self._positions = calc_positions(X[0], X[1], X[2], N) diff --git a/examples/shock_tube.py b/examples/shock_tube.py new file mode 100644 index 0000000..2a6e61b --- /dev/null +++ b/examples/shock_tube.py @@ -0,0 +1,32 @@ +import dsmc + +if __name__ == '__main__': + # general parameters + solver = dsmc.DSMC() + domain = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (0, 50e-3)) + dt = 1e-7 + w = 2.4134e+10 + mass = 6.6422e-26 + niter = 300 + + # low denisty particles + nhigh = 2.5e+14 + Thigh = 300 + Boxhigh = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (25e-3, 50e-3)) + + # high denisty particles + nlow = 2.5e+13 + Tlow = 300 + Boxlow = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (0, 25e-3)) + + solver.set_domain(domain) + solver.set_weight(w) + + solver.create_particles(Boxlow, mass, Tlow, nlow) + solver.create_particles(Boxhigh, mass, Thigh, nhigh) + + for it in range(niter): + print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) + + print("") + print('done') diff --git a/tests/unit/test_dsmc/dsmc.py b/tests/unit/test_dsmc/dsmc.py index 128cccb..db710a8 100644 --- a/tests/unit/test_dsmc/dsmc.py +++ b/tests/unit/test_dsmc/dsmc.py @@ -1,6 +1,6 @@ import unittest -from dsmc import dsmc +from dsmc import dsmc as ds class TestDSMC(unittest.TestCase): - def test_test(self): - self.assertEqual(1.0, dsmc.test()) + def test_Constructor(self): + ds.DSMC() From d1b16f466b9414b32815dd468070f9ccb67e903e Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Thu, 22 Sep 2022 15:47:48 +0200 Subject: [PATCH 03/19] fixed bug --- dsmc/particles.py | 4 ++-- tests/unit/test_dsmc/particles.py | 11 +++++++++++ 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/dsmc/particles.py b/dsmc/particles.py index a978c34..5f588f4 100644 --- a/dsmc/particles.py +++ b/dsmc/particles.py @@ -134,6 +134,6 @@ def create_particles(self, X, mass, T, N): 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._positions = np.concatenate(self._positions, calc_positions(X[0], X[1], X[2], N)) + self._velocities = np.concatenate((self._velocities, get_velocities(T, mass, N))) + self._positions = np.concatenate((self._positions, calc_positions(X[0], X[1], X[2], N))) self._N += N diff --git a/tests/unit/test_dsmc/particles.py b/tests/unit/test_dsmc/particles.py index e27497f..1921691 100644 --- a/tests/unit/test_dsmc/particles.py +++ b/tests/unit/test_dsmc/particles.py @@ -63,3 +63,14 @@ 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]) + + particles.create_particles(X, mass, T, N) + + self.assertEqual(2*N, len(particles.Pos)) + self.assertEqual(2*N, len(particles.Vel)) + self.assertEqual(2*N, particles.N) + + for i in range(2*N): + 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]) From 96b65b223152c7854af2e61d057734df813dd65f Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Thu, 22 Sep 2022 15:50:35 +0200 Subject: [PATCH 04/19] fixed bug --- dsmc/dsmc.py | 4 ++-- examples/shock_tube.py | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index b587c9a..09ee7df 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -15,9 +15,9 @@ def clear(self): self.domain = None def advance(self, dt): - if self.domain == None: + if self.domain is None: raise Exception("simulation domain not defined") - if particles.N == 0: + if self.particles.N == 0: raise Exception("no particles created") if self.w == None: raise Exception("particle weight not set") diff --git a/examples/shock_tube.py b/examples/shock_tube.py index 2a6e61b..9642967 100644 --- a/examples/shock_tube.py +++ b/examples/shock_tube.py @@ -27,6 +27,7 @@ for it in range(niter): print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) + solver.advance(dt) print("") print('done') From 56c6a81d66c2987409d4f5564e538b70449604a5 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Thu, 22 Sep 2022 16:04:59 +0200 Subject: [PATCH 05/19] implemented velocity update --- dsmc/dsmc.py | 9 +++++++++ dsmc/particles.py | 5 ++--- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index 09ee7df..7c84ba4 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -4,6 +4,10 @@ from . import particles as prt from . import octree as oc +@njit +def _push(velocities, positions, dt): + return positions + velocities*dt + class DSMC: def __init__(self): self.clear() @@ -21,6 +25,11 @@ def advance(self, dt): raise Exception("no particles created") if self.w == None: raise Exception("particle weight not set") + + self.octree.build(self.particles.Pos) + # update velocities + self.particles.VelPos = (self.particles.Vel, _push(self.particles.Vel, self.particles.Pos, dt)) + # do boundary def create_particles(self, box, mass, T, n): N = int(round(n / self.w)) diff --git a/dsmc/particles.py b/dsmc/particles.py index 5f588f4..9ca72eb 100644 --- a/dsmc/particles.py +++ b/dsmc/particles.py @@ -121,9 +121,8 @@ def VelPos(self): @VelPos.setter def VelPos(self, vel_pos): - assert vel_pos[0] == vel_pos[1] - assert isinstance(vel_pos[0], np.array) - assert isinstance(vel_pos[1], np.array) + assert len(vel_pos[0]) == len(vel_pos[1]) + self._velocities = vel_pos[0] self._positions = vel_pos[1] self._N = len(self._positions) From 66dd8ae1b73c8bd088bdc8cb9decbe4da61e0637 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Thu, 22 Sep 2022 16:28:32 +0200 Subject: [PATCH 06/19] implemented boundary --- dsmc/dsmc.py | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index 7c84ba4..86b769f 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -5,9 +5,23 @@ from . import octree as oc @njit -def _push(velocities, positions, dt): +def _push(velocities, positions, dt): return positions + velocities*dt +@njit +def _boundary(velocities, positions, domain): + for p in range(len(positions)): + while not oc._is_inside(positions[p], domain): + 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 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) + class DSMC: def __init__(self): self.clear() @@ -28,8 +42,8 @@ def advance(self, dt): self.octree.build(self.particles.Pos) # update velocities - self.particles.VelPos = (self.particles.Vel, _push(self.particles.Vel, self.particles.Pos, dt)) - # do boundary + positions = _push(self.particles.Vel, self.particles.Pos, dt) + self.particles.VelPos = _boundary(self.particles.Vel, positions, self.domain) def create_particles(self, box, mass, T, n): N = int(round(n / self.w)) From dde585ce5bb4a9f3028d8001f4f0eb1a41f761f9 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Thu, 22 Sep 2022 17:08:08 +0200 Subject: [PATCH 07/19] added shocktube example --- examples/plot.py | 42 ++++++++++++++++++++++++++++++++++++++++++ examples/shock_tube.py | 14 ++++++++++---- 2 files changed, 52 insertions(+), 4 deletions(-) create mode 100644 examples/plot.py diff --git a/examples/plot.py b/examples/plot.py new file mode 100644 index 0000000..a2231bf --- /dev/null +++ b/examples/plot.py @@ -0,0 +1,42 @@ +import matplotlib.pyplot as plt +import csv +import numpy as np + +if __name__ == '__main__': + res = [] + + with open('test.csv') as file: + reader = csv.reader(file, delimiter=',') + + for line in reader: + l = [m for m in line if m] + data = (np.array(l)).astype(float) + data = np.sort(data) + + N = 100 + sor = np.zeros((N, )) + x = np.zeros((N, )) + step = data[-1] / N + q = 0 + print(step, data.shape) + + for i in range(len(data)): + if data[i] < step: + x[q] = data[i] + sor[q] += 1 + else: + q += 1 + step += step + sor[q] += 1 + x[q] = data[i] + + + x.resize(q) + sor.resize(q) + + plt.plot(x, sor) + plt.show() + plt.hist(data) + plt.show() + break + diff --git a/examples/shock_tube.py b/examples/shock_tube.py index 9642967..9d7c5f6 100644 --- a/examples/shock_tube.py +++ b/examples/shock_tube.py @@ -4,7 +4,7 @@ # general parameters solver = dsmc.DSMC() domain = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (0, 50e-3)) - dt = 1e-7 + dt = 1e-6 w = 2.4134e+10 mass = 6.6422e-26 niter = 300 @@ -25,9 +25,15 @@ solver.create_particles(Boxlow, mass, Tlow, nlow) solver.create_particles(Boxhigh, mass, Thigh, nhigh) - for it in range(niter): - print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) - solver.advance(dt) + with open("test.csv", "w") as file: + for it in range(niter): + print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True) + solver.advance(dt) + + for pos in solver.particles.Pos: + file.write("{:.4e},".format(pos[2])) + + file.write("\n") print("") print('done') From e4242a00053761d5c2ce0e9834f1708e6e110594 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 08:53:56 +0200 Subject: [PATCH 08/19] modified plot script --- examples/plot.py | 26 ++++++++++++-------------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/examples/plot.py b/examples/plot.py index a2231bf..8e1d08c 100644 --- a/examples/plot.py +++ b/examples/plot.py @@ -7,6 +7,7 @@ with open('test.csv') as file: reader = csv.reader(file, delimiter=',') + res = [] for line in reader: l = [m for m in line if m] @@ -16,27 +17,24 @@ N = 100 sor = np.zeros((N, )) x = np.zeros((N, )) - step = data[-1] / N + dx = 0.05/N + x[0] = dx q = 0 - print(step, data.shape) for i in range(len(data)): - if data[i] < step: - x[q] = data[i] - sor[q] += 1 - else: + while data[i] > x[q]: q += 1 - step += step - sor[q] += 1 - x[q] = data[i] + x[q] = x[q - 1] + dx + + sor[q] += 1 x.resize(q) sor.resize(q) - plt.plot(x, sor) - plt.show() - plt.hist(data) + res.append((x, sor)) + + for i in range(len(res)): + if i%100 == 0: + plt.plot(res[i][0], res[i][1]) plt.show() - break - From 40410b1ab589957964eb9d7984d9b40bce7cfed8 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 08:55:46 +0200 Subject: [PATCH 09/19] added sigma T to dsmc --- dsmc/dsmc.py | 1 + 1 file changed, 1 insertion(+) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index 86b769f..7994fb6 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -31,6 +31,7 @@ def clear(self): self.octree = oc.Octree() self.w = None self.domain = None + self.sigma_T = 3.631681e-19 def advance(self, dt): if self.domain is None: From 2dd13d6f5681cbd9f4bf3d99711e1b0269e384b9 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 09:20:46 +0200 Subject: [PATCH 10/19] implemneted _calc_prob --- dsmc/dsmc.py | 23 +++++++++++++++++++++++ tests/unit/test_dsmc/dsmc.py | 14 ++++++++++++++ 2 files changed, 37 insertions(+) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index 7994fb6..cd3b764 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -21,6 +21,29 @@ def _boundary(velocities, positions, domain): velocities[p][i] *= -1.0 return (velocities, positions) + +@njit +def _calc_prob(vel1 : np.ndarray, vel2 : np.ndarray, sigma_T : float, Vc : float, dt : float, w : float, N : int) -> np.single: + """ + Parameters + ---------- + vel1 : velocity + vel2 : velocity + sigma_T : float + total cross section + Vc : float + cell volume + w : float + weight + N : int + number of particles + + Returns + ------- + collision proability : float + """ + return np.linalg.norm(vel1 - vel2) * sigma_T * dt * w * N / Vc; + class DSMC: def __init__(self): diff --git a/tests/unit/test_dsmc/dsmc.py b/tests/unit/test_dsmc/dsmc.py index db710a8..37be91a 100644 --- a/tests/unit/test_dsmc/dsmc.py +++ b/tests/unit/test_dsmc/dsmc.py @@ -1,6 +1,20 @@ import unittest +import numpy as np from dsmc import dsmc as ds class TestDSMC(unittest.TestCase): def test_Constructor(self): ds.DSMC() + + def test__calc_prob(self): + vel1 : np.ndarray = np.array([1.0, 2.0, 3.0]) + vel2 : np.ndarray = np.array([1.0, 2.0, 3.0]) + sigma_T : float = 1.0 + Vc : float = 1.0 + dt : float = 1.0 + w : float = 1.0 + N : int = 1 + + res = ds._calc_prob(vel1, vel2, sigma_T, Vc, dt, w, N) + + self.assertEqual(np.linalg.norm(vel1 - vel2), res) From 7b5c62e54f7af7082a2c8567e151adfb497db797 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 09:35:47 +0200 Subject: [PATCH 11/19] implemented _calc_post_col_vels --- dsmc/dsmc.py | 19 +++++++++++++++++++ tests/unit/test_dsmc/dsmc.py | 19 +++++++++++++++++++ 2 files changed, 38 insertions(+) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index cd3b764..589bfdc 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -43,7 +43,26 @@ def _calc_prob(vel1 : np.ndarray, vel2 : np.ndarray, sigma_T : float, Vc : float collision proability : float """ return np.linalg.norm(vel1 - vel2) * sigma_T * dt * w * N / Vc; + +@njit +def _calc_post_col_vels(velocity1 : np.ndarray, velocity2 : np.ndarray, mass1 : float, mass2 : float, rel_vel_module : float, rand_number1 : float, rand_number2 : float) -> tuple[np.ndarray, np.ndarray]: + mass12 = (mass1 + mass2) + mass1_12 = (mass1 / mass12) + mass2_12 = (mass2 / mass12) + cos_xi = (2.0 * rand_number1 - 1.0) + sin_xi = (np.sqrt(1.0 - cos_xi * cos_xi)) + epsilon = (2.0 * np.pi * rand_number2) + + centre_of_mass_velocity = (velocity1 * mass1 + velocity2 * mass2) * (1.0 / mass12) + + rel_velocity_new = np.empty((3, )) + + rel_velocity_new[0] = rel_vel_module * cos_xi + rel_velocity_new[1] = rel_vel_module * sin_xi * np.cos(epsilon) + rel_velocity_new[2] = rel_vel_module * sin_xi * np.sin(epsilon) + + return (centre_of_mass_velocity + rel_velocity_new * mass2_12 , centre_of_mass_velocity - rel_velocity_new * mass1_12) class DSMC: def __init__(self): diff --git a/tests/unit/test_dsmc/dsmc.py b/tests/unit/test_dsmc/dsmc.py index 37be91a..90cdf10 100644 --- a/tests/unit/test_dsmc/dsmc.py +++ b/tests/unit/test_dsmc/dsmc.py @@ -18,3 +18,22 @@ def test__calc_prob(self): res = ds._calc_prob(vel1, vel2, sigma_T, Vc, dt, w, N) self.assertEqual(np.linalg.norm(vel1 - vel2), res) + + def test__calc_post_col_vels(self): + velocity1 : np.ndarray = np.array([1.0, 2.0, 3.0]) + velocity2 : np.ndarray = np.array([1.0, 2.0, 3.0]) + mass1 : float = 1.0 + mass2 : float = 1.0 + rel_vel_module : float = 1.0 + rand_number1 : float = 1.0 + rand_number2 : float = 1.0 + + res = ds._calc_post_col_vels(velocity1, velocity2, mass1, mass2, rel_vel_module, rand_number1, rand_number2) + + self.assertEqual(1.5, res[0][0]) + self.assertEqual(2.0, res[0][1]) + self.assertEqual(3.0, res[0][2]) + + self.assertEqual(0.5, res[1][0]) + self.assertEqual(2.0, res[1][1]) + self.assertEqual(3.0, res[1][2]) From da6965c7e52f4e2e1fd76c67c7c01dbaab98cb06 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 10:05:07 +0200 Subject: [PATCH 12/19] implemented function get_V --- dsmc/octree.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/dsmc/octree.py b/dsmc/octree.py index 81b18ac..20cf959 100644 --- a/dsmc/octree.py +++ b/dsmc/octree.py @@ -110,6 +110,10 @@ def _sort(permutations : npt.NDArray, box : npt.NDArray, positions : npt.NDArray Nnew += 1 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]) class Leaf: def __init__(self): From 614a3a5ce1fd788477efde0591b647a6b6507728 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 10:06:11 +0200 Subject: [PATCH 13/19] implemented velocity update in dsmc --- dsmc/dsmc.py | 34 +++++++++++++++++++++++++++++++--- 1 file changed, 31 insertions(+), 3 deletions(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index 589bfdc..b399ed7 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -63,6 +63,24 @@ def _calc_post_col_vels(velocity1 : np.ndarray, velocity2 : np.ndarray, mass1 : rel_velocity_new[2] = rel_vel_module * sin_xi * np.sin(epsilon) return (centre_of_mass_velocity + rel_velocity_new * mass2_12 , centre_of_mass_velocity - rel_velocity_new * mass1_12) + +@njit +def _update_velocities(permutations : np.ndarray, velocities : np.ndarray, mass : float, sigma_T : float, Vc : float, dt : float, w : float, offset : int, N : int): + i = 1 + while i < N: + p1 = permutations[offset + i - 1] + p2 = permutations[offset + i] + P = _calc_prob(velocities[p1], velocities[p2], sigma_T, Vc, dt, w, N) + R = np.random.random(3) + + if R[0] < P: + new_vels = _calc_post_col_vels(velocities[p1], velocities[p2], mass, mass, np.linalg.norm(velocities[p1] - velocities[p2]), R[1], R[2]) + velocities[p1] = new_vels[0] + velocities[p2] = new_vels[1] + + i += 2 + + return velocities class DSMC: def __init__(self): @@ -74,6 +92,7 @@ def clear(self): self.w = None self.domain = None self.sigma_T = 3.631681e-19 + self.mass = None def advance(self, dt): if self.domain is None: @@ -84,20 +103,29 @@ def advance(self, dt): raise Exception("particle weight not set") self.octree.build(self.particles.Pos) - # update velocities + + for i in range(len(self.octree.leafs)): + leaf = self.octree.leafs[i] + if not leaf.number_children: + Vc = oc.get_V(self.octree.cell_boxes[i]) + self.particles.VelPos = (_update_velocities(self.octree.permutations, self.particles.Vel, self.mass, self.sigma_T, Vc, dt, self.w, leaf.elem_offset, leaf.number_elements), self.particles.Pos) + positions = _push(self.particles.Vel, self.particles.Pos, dt) self.particles.VelPos = _boundary(self.particles.Vel, positions, self.domain) - def create_particles(self, box, mass, T, n): + def create_particles(self, box, T, n): N = int(round(n / self.w)) print("creating {} particles".format(N)) - self.particles.create_particles(box, mass, T, N) + self.particles.create_particles(box, self.mass, T, N) print("now containing {} particles, {} total".format(N, self.particles.N)) def set_domain(self, domain): self.domain = np.array(domain) + def set_mass(self, mass): + self.mass = mass + def set_weight(self, w): self.octree.w = w self.w = w From 71569cfe5c611d6e568345ad956c4a4e6c02c7a7 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 10:07:02 +0200 Subject: [PATCH 14/19] implemented new shock tube interface --- examples/shock_tube.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/examples/shock_tube.py b/examples/shock_tube.py index 9d7c5f6..2952a5b 100644 --- a/examples/shock_tube.py +++ b/examples/shock_tube.py @@ -5,7 +5,7 @@ solver = dsmc.DSMC() domain = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (0, 50e-3)) dt = 1e-6 - w = 2.4134e+10 + w = 2.4134e+11 mass = 6.6422e-26 niter = 300 @@ -21,9 +21,10 @@ solver.set_domain(domain) solver.set_weight(w) + solver.set_mass(mass) - solver.create_particles(Boxlow, mass, Tlow, nlow) - solver.create_particles(Boxhigh, mass, Thigh, nhigh) + solver.create_particles(Boxlow, Tlow, nlow) + solver.create_particles(Boxhigh, Thigh, nhigh) with open("test.csv", "w") as file: for it in range(niter): From 3bc3519aab415940854f1e1314a206f8d0505c4b Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 10:14:06 +0200 Subject: [PATCH 15/19] added _update_velocities --- dsmc/dsmc.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index b399ed7..9879ab3 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -65,7 +65,7 @@ def _calc_post_col_vels(velocity1 : np.ndarray, velocity2 : np.ndarray, mass1 : return (centre_of_mass_velocity + rel_velocity_new * mass2_12 , centre_of_mass_velocity - rel_velocity_new * mass1_12) @njit -def _update_velocities(permutations : np.ndarray, velocities : np.ndarray, mass : float, sigma_T : float, Vc : float, dt : float, w : float, offset : int, N : int): +def _update_velocities(permutations : np.ndarray, velocities : np.ndarray, mass : float, sigma_T : float, Vc : float, dt : float, w : float, offset : int, N : int) -> np.ndarray: i = 1 while i < N: p1 = permutations[offset + i - 1] @@ -81,6 +81,14 @@ def _update_velocities(permutations : np.ndarray, velocities : np.ndarray, mass i += 2 return velocities + +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]: + Vc = oc.get_V(cell_boxes[i]) + velocities = _update_velocities(permutations, velocities, mass, sigma_T, Vc, dt, w, elem_offsets[i], number_elements[i]) + + return velocities class DSMC: def __init__(self): From 27880a14b0da85e61865429dac8d7a8f23c72e32 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 10:33:38 +0200 Subject: [PATCH 16/19] fixed bug when building octree --- dsmc/octree.py | 1 + 1 file changed, 1 insertion(+) diff --git a/dsmc/octree.py b/dsmc/octree.py index 20cf959..09b3c33 100644 --- a/dsmc/octree.py +++ b/dsmc/octree.py @@ -141,6 +141,7 @@ def clear(self): self.level = 0 def build(self, positions): + self.clear() self._create_root(positions) self.permutations = np.array([i for i in range(len(positions))]) From 619f69d647cc8c7f9acb70384591562910c78775 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 10:34:01 +0200 Subject: [PATCH 17/19] refactored velocity update function in dsmc --- dsmc/dsmc.py | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/dsmc/dsmc.py b/dsmc/dsmc.py index 9879ab3..5af06e9 100644 --- a/dsmc/dsmc.py +++ b/dsmc/dsmc.py @@ -81,7 +81,8 @@ def _update_velocities(permutations : np.ndarray, velocities : np.ndarray, mass i += 2 return velocities - + +@njit 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]: @@ -112,15 +113,19 @@ def advance(self, dt): self.octree.build(self.particles.Pos) - for i in range(len(self.octree.leafs)): - leaf = self.octree.leafs[i] - if not leaf.number_children: - Vc = oc.get_V(self.octree.cell_boxes[i]) - self.particles.VelPos = (_update_velocities(self.octree.permutations, self.particles.Vel, self.mass, self.sigma_T, Vc, dt, self.w, leaf.elem_offset, leaf.number_elements), self.particles.Pos) - + 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) + def _update_velocities(self, dt): + 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) + def create_particles(self, box, T, n): N = int(round(n / self.w)) print("creating {} particles".format(N)) From 17405f9d03e7899b575b5b491f097f23c1f375e9 Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 10:39:30 +0200 Subject: [PATCH 18/19] extended shocktube test case --- examples/plot.py | 6 ++---- examples/shock_tube.py | 8 ++++---- 2 files changed, 6 insertions(+), 8 deletions(-) diff --git a/examples/plot.py b/examples/plot.py index 8e1d08c..b4e91ae 100644 --- a/examples/plot.py +++ b/examples/plot.py @@ -34,7 +34,5 @@ res.append((x, sor)) - for i in range(len(res)): - if i%100 == 0: - plt.plot(res[i][0], res[i][1]) - plt.show() + plt.plot(res[-1][0], res[-1][1]) + plt.show() diff --git a/examples/shock_tube.py b/examples/shock_tube.py index 2952a5b..9b61c08 100644 --- a/examples/shock_tube.py +++ b/examples/shock_tube.py @@ -4,18 +4,18 @@ # general parameters solver = dsmc.DSMC() domain = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (0, 50e-3)) - dt = 1e-6 - w = 2.4134e+11 + dt = 1e-7 + w = 2.4134e+15 mass = 6.6422e-26 niter = 300 # low denisty particles - nhigh = 2.5e+14 + nhigh = 2.5e+20 Thigh = 300 Boxhigh = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (25e-3, 50e-3)) # high denisty particles - nlow = 2.5e+13 + nlow = 2.5e+19 Tlow = 300 Boxlow = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (0, 25e-3)) From 99b827965d13c5367f1af837b7b4f0de1294c7ee Mon Sep 17 00:00:00 2001 From: Leo Basov Date: Fri, 23 Sep 2022 10:41:03 +0200 Subject: [PATCH 19/19] updated changelog and setup.py --- CHANGELOG.md | 4 ++++ setup.py | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0af8d15..74f67f4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.3.0] - 2022-09-23 +### Added +- dsmc module + ## [0.3.0] - 2022-09-22 ### Added - octree module diff --git a/setup.py b/setup.py index 13e9750..7963a05 100644 --- a/setup.py +++ b/setup.py @@ -1,7 +1,7 @@ from setuptools import setup setup( name='dsmc', - version='0.3.0', + version='0.4.0', author='Leo Basov', python_requires='>=3.6, <4', packages=["dsmc"],