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/dsmc/common.py b/dsmc/common.py
new file mode 100644
index 0000000..0990d48
--- /dev/null
+++ b/dsmc/common.py
@@ -0,0 +1,5 @@
+from numba import cfunc
+
+@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/dsmc/dsmc.py b/dsmc/dsmc.py
index 3c13811..4b7d212 100644
--- a/dsmc/dsmc.py
+++ b/dsmc/dsmc.py
@@ -3,26 +3,71 @@
from numba import prange
from . import particles as prt
from . import octree as oc
+from . import common as com
-@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)
-def _boundary(velocities, positions, domain):
+@njit
+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] == 1 or 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
+ 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 or boundary_conds[i][0] == 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)
- return (velocities, 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:
@@ -86,11 +131,41 @@ 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
+@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))
+
class DSMC:
def __init__(self):
self.clear()
@@ -100,6 +175,8 @@ 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 = open, 2 = inflow
+ self.boundary = Boundary()
self.sigma_T = 3.631681e-19
self.mass = None
@@ -107,16 +184,23 @@ 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")
+
+ 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)
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.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)
@@ -127,20 +211,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))
+ N = int(round(com.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))
@@ -153,3 +228,20 @@ def set_mass(self, mass):
def set_weight(self, w):
self.octree.w = w
self.w = w
+
+ 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/dsmc/octree.py b/dsmc/octree.py
index 18afcd4..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
@@ -110,10 +109,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 7e709c3..421fbff 100644
--- a/dsmc/particles.py
+++ b/dsmc/particles.py
@@ -1,9 +1,14 @@
import math
import numpy as np
from numba import njit
+from . import common as com
kb = 1.380649e-23
+@njit
+def calc_vp(T, mass):
+ return np.sqrt(2*kb*T/mass)
+
@njit
def box_muller(T):
"""
@@ -53,15 +58,15 @@ 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
-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
@@ -128,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(com.get_V(box) * n / w))
+
+ self.create_particles(box, mass, T, N, u)
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")
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()
diff --git a/examples/push_bound_inflow_test.py b/examples/push_bound_inflow_test.py
new file mode 100644
index 0000000..d7fdfdb
--- /dev/null
+++ b/examples/push_bound_inflow_test.py
@@ -0,0 +1,57 @@
+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
+ u = np.array([1000.0, 0.0, 0.0])
+ 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_bc_type("xmax", "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_bc_type("xmin", "inflow")
+
+ solver.set_bc_values("xmin", T, n, u)
+
+ 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
diff --git a/examples/push_bound_open_test.py b/examples/push_bound_open_test.py
new file mode 100644
index 0000000..b899097
--- /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 = 1000
+
+ # setup solver
+ solver.set_domain(domain)
+ solver.set_weight(w)
+ solver.set_mass(mass)
+
+ solver.create_particles(domain, T, n)
+
+ 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)
+
+ 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
diff --git a/examples/shock_tube.py b/examples/shock_tube.py
index 37ef606..37d60c0 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
@@ -32,7 +43,7 @@ def write2file(solver, n_file, T_file, p_file, Nbins):
Boxhigh = [(-0.0001, 0.0001), (-0.0001, 0.0001), (0.0, 0.05)]
# low denisty particles
- nlow = 2.41432e+21
+ nlow = nhigh*0.1
Tlow = 300
Boxlow = [(-0.0001, 0.0001), (-0.0001, 0.0001), (0.05, 0.1)]
@@ -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')
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',
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..b6e73ec
--- /dev/null
+++ b/tests/unit/test_dsmc/common.py
@@ -0,0 +1,9 @@
+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)
+ self.assertEqual(8.0, V)
\ No newline at end of file
diff --git a/tests/unit/test_dsmc/octree.py b/tests/unit/test_dsmc/octree.py
index 5eb9a57..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,12 +311,10 @@ 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))
-
- print(boxes_new)
class TestOctreeOctree(unittest.TestCase):
def test_build(self):
diff --git a/tests/unit/test_dsmc/particles.py b/tests/unit/test_dsmc/particles.py
index 1921691..700c46f 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):
@@ -20,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))
@@ -28,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)
@@ -51,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))
@@ -64,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))
@@ -74,3 +78,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)
diff --git a/tools/sod_analytical.py b/tools/sod_analytical.py
new file mode 100644
index 0000000..9e0c1de
--- /dev/null
+++ b/tools/sod_analytical.py
@@ -0,0 +1,141 @@
+# -*- 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
+
+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
+ 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 = p1 * 1e-3
+ 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_val(x, val, name):
+ 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")
+
+ 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]))
+ 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]):
+ 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', 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')
+
+ args = parser.parse_args()
+
+ check_args(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)
+ mass = 6.6422e-26
+
+ # boundary conditions
+ rho = np.zeros(5)
+ p = np.zeros(5)
+ u = np.zeros(5)
+
+ rho[0] = args.rho[0]
+ p[0] = args.p[0]
+ u[0] = 0.0
+
+ rho[-1] = args.rho[1]
+ p[-1] = args.p[1]
+ u[-1] = 0.0
+
+ # calculating states
+ p[1] = p[0]
+ 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])
+
+ 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")
+ 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")