Boundary conditions

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

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])

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,23 +175,32 @@ 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
 
     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))

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)])

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)

dsmc/writer.py

@@ -19,8 +19,9 @@ def _wrtie_body(f, octree):
 	f.write("<UnstructuredGrid>\n")
 	f.write("<Piece NumberOfPoints=\"{}\" NumberOfCells=\"{}\">\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("</Piece>\n")
 	f.write("</UnstructuredGrid>\n")
@@ -93,6 +94,16 @@ def _write_cells(f, octree, leaf_ids):
 
 	f.write("</DataArray>\n")
 	f.write("</Cells>\n")
+
+def _write_cell_data(f, octree, leaf_ids):
+    f.write("<CellData Scalars=\"number_density\">\n")
+    f.write("<DataArray type=\"Float32\" Name=\"particle_numbers\" format=\"ascii\">\n")
+
+    for i in leaf_ids:
+        f.write("{} ".format(octree.leafs[i].number_elements))
+
+    f.write("</DataArray>\n")
+    f.write("</CellData>\n")
 
 def _write_footer(f):
 	f.write("</VTKFile>\n")

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()