implemented cuboid domain object

CHANGELOG.md

@@ -6,6 +6,19 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+## [0.9.0] - 2022-10-09
+## Added
+- mesh2d as for diagnostic purposes
+- particle writer
+- docs
+- hypersonic test case
+- backup boundary module
+
+## [0.8.0] - 2022-10-03
+### Added
+- planar writer
+- planar diagnostic
+
 ## [0.7.0] - 2022-10-01
 ### Added
 - inflow boundary condition

doc/Boundary.ipynb

@@ -0,0 +1,167 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "d49edbe1-fe92-4a67-82e9-78a877139d5d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import sympy as sym"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7f598afc-282c-4bb5-8d47-80c77fcc1178",
+   "metadata": {},
+   "source": [
+    "**Boundary collision method**\n",
+    "\n",
+    "Given a plane $\\Pi$ and and a line $L$ as\n",
+    "\n",
+    "$$\n",
+    "\\Pi \\rightarrow (p - p_0) \\cdot \\vec{n}_p = 0 \\\\\n",
+    "L \\rightarrow l = l_0 + t \\cdot \\vec{n}_l\n",
+    "$$\n",
+    "\n",
+    "where $\\vec{n}_1$ is the normal vector of a plane given as\n",
+    "\n",
+    "$$\n",
+    "\\vec{n}_p = (p_1 - p_0) \\times (p_2 - p_0)\n",
+    "$$\n",
+    "\n",
+    "is a vector in the direction of the line as\n",
+    "\n",
+    "$$\n",
+    "\\vec{n}_l = l_1 - l_0.\n",
+    "$$\n",
+    "\n",
+    "The positons of the intersection point found by setting the point $p$ in the equation for the plane as being a point in the equation of the line:\n",
+    "\n",
+    "$$\n",
+    "(l - p_0) \\cdot \\vec{n}_p = 0 = (l_0 + t \\cdot \\vec{n}_l - p_0) \\cdot \\vec{n}_p \\implies t = - \\frac{(l_0 - p_0) \\cdot \\vec{n}_p}{\\vec{n}_p \\cdot \\vec{n}_l}\n",
+    "$$\n",
+    "\n",
+    "inserting the line equation the intersection point can be found as\n",
+    "\n",
+    "$$\n",
+    "l^* = l_0 + t \\cdot \\vec{n}_l.\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "4abb0e41-5e15-4326-b663-df0b6b5cfe0b",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[0. 0. 0.]\n"
+     ]
+    }
+   ],
+   "source": [
+    "l1 = np.array([1.0, 1.0, 0.0])\n",
+    "l0 = np.array([-1.0, -1.0, 0.0])\n",
+    "\n",
+    "p0 = np.array([0.0, -1.0, -1.0])\n",
+    "p1 = np.array([0.0, -1.0, 1.0])\n",
+    "p2 = np.array([0.0, 1.0, 1.0])\n",
+    "\n",
+    "nl = l1 - l0\n",
+    "n_p = np.cross((p1 - p0), (p2 - p0))\n",
+    "\n",
+    "t = - ((l0 - p0).dot(n_p) / n_p.dot(nl))\n",
+    "\n",
+    "ls = l0 +t*nl\n",
+    "\n",
+    "print(ls)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d43f6e71-7324-4b21-9c0d-371fb76d01f6",
+   "metadata": {},
+   "source": [
+    "**Reflections**\n",
+    "\n",
+    "The reflected line can be found as\n",
+    "\n",
+    "$$\n",
+    "L_R \\rightarrow l_R = l^* + \\lambda \\cdot \\vec{n}_R\n",
+    "$$\n",
+    "\n",
+    "with \n",
+    "\n",
+    "$$\n",
+    "\\vec{n}_R = \\vec{n}_l - 2 \\frac{\\vec{n}_p \\cdot \\vec{n}_l}{|| \\vec{n}_p ||^2} \\cdot \\vec{n}_p.\n",
+    "$$\n",
+    "\n",
+    "The reflected point can now the found as \n",
+    "\n",
+    "$$\n",
+    "l_R = l^* + (1 - t) \\cdot \\vec{n}_R\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "91de06b7-7511-471d-9b64-f9c2a784d682",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([ 1., -1.,  0.])"
+      ]
+     },
+     "execution_count": 3,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "nR = nl - 2*(n_p.dot(nl) / n_p.dot(n_p))*n_p\n",
+    "lR = ls - (1  -t )*nR\n",
+    "\n",
+    "lR"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a81e9e3a-d8d0-463f-a5e4-b8700e696875",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

dsmc/boundary.py

@@ -0,0 +1,180 @@
+import numpy as np
+from numba import njit
+from . import common as co
+
+@njit
+def _check_if_parallel(v1, v2, diff=1e-6):
+    n1 = np.linalg.norm(v1)
+    n2 = np.linalg.norm(v2)
+
+    if n1 < 1.0e-13 or n2 < 1.0e-13:
+        return False
+
+    V1 = np.copy(v1) / n1
+    V2 = np.copy(v2) / n2
+
+    return V1.dot(V2) > diff
+
+@njit
+def _intersect(l0, l1, p0, p1, p2):
+    """
+    Args:
+        l0 : first position on line
+        l1 : second position on line
+        p0 : first position on plane
+        p1 : first position on plane
+        p2 : first position on plane
+
+    Returns:
+        (intersected, n_l, n_r, t)
+    """
+    n_l = l1 - l0
+    n_p = np.cross((p1 - p0), (p2 - p1))
+
+    if _check_if_parallel(n_l, n_p):
+        return (True, n_l, n_p, - ((l0 - p0).dot(n_p) / n_p.dot(n_l)))
+    else:
+        return (False, n_l, n_p, 0.0)
+
+@njit    
+def _calc_nr(n_l, n_p):
+    return n_l - 2.0 * (n_p.dot(n_l) / n_p.dot(n_p))*n_p
+
+@njit     
+def _reflect(vel, pos, pos_old, p0, p1, p2, domain):
+    intersected, n_l, n_p, t = _intersect(pos_old, pos, p0, p1, p2)
+
+    if intersected and t < 1.0 and t > 0.0:
+        k = 1.0
+        p = pos_old + n_l*(t*k)
+        while not co.is_inside(p, domain):
+            k *= 0.9
+            p = pos_old + n_l*(t*k)
+        pos_old = p
+        n_r = _calc_nr(n_l, n_p)
+        pos = pos_old  + (1.0 - t)*n_r
+        vel = (np.linalg.norm(vel) / np.linalg.norm(n_r))*n_r
+
+    return (vel, pos, pos_old)
+
+@njit
+def _get_plane(domain, i, j):
+    if i == 0:
+        if j == 0:
+            p0 = np.array([domain[i][j], domain[1][0], domain[2][0]])
+            p1 = np.array([domain[i][j], domain[1][0], domain[2][1]])
+            p2 = np.array([domain[i][j], domain[1][1], domain[2][0]])
+        elif j == 1:
+            p0 = np.array([domain[i][j], domain[1][0], domain[2][0]])
+            p1 = np.array([domain[i][j], domain[1][1], domain[2][0]])
+            p2 = np.array([domain[i][j], domain[1][0], domain[2][1]])
+    elif i == 1:
+        if j == 0:
+            p0 = np.array([domain[0][0], domain[i][j], domain[2][0]])
+            p1 = np.array([domain[0][1], domain[i][j], domain[2][0]])
+            p2 = np.array([domain[0][0], domain[i][j], domain[2][1]])
+        if j == 1:
+            p0 = np.array([domain[0][0], domain[i][j], domain[2][0]])
+            p1 = np.array([domain[0][0], domain[i][j], domain[2][1]])
+            p2 = np.array([domain[0][1], domain[i][j], domain[2][0]])
+    elif i == 2:
+        if j == 0:
+            p0 = np.array([domain[0][0], domain[1][0], domain[i][j]])
+            p1 = np.array([domain[0][0], domain[1][1], domain[i][j]])
+            p2 = np.array([domain[0][1], domain[1][0], domain[i][j]])
+        if j == 1:
+            p0 = np.array([domain[0][0], domain[1][0], domain[i][j]])
+            p1 = np.array([domain[0][1], domain[1][0], domain[i][j]])
+            p2 = np.array([domain[0][0], domain[1][1], domain[i][j]])
+
+    return (p0, p1, p2)
+
+@njit
+def _boundary(velocities, positions, old_positions, domain, boundary_conds):
+    kept_parts = np.ones(positions.shape[0], dtype=np.uint)
+
+    for p in range(len(positions)):
+        counter = 0
+        while not co.is_inside(positions[p], domain) and kept_parts[p]:
+            if counter > 10:
+                kept_parts[p] = 0
+                break
+
+            for i in range(3):
+                for j in range(2):
+                    p0, p1, p2 = _get_plane(domain, i, j)
+                    if boundary_conds[i][j] == 0:
+                        velocities[p], positions[p], old_positions[p] = _reflect(velocities[p], positions[p], old_positions[p], p0, p1, p2, domain)
+                        counter += 1
+                    elif boundary_conds[i][j] == 1 or boundary_conds[i][j] == 2:
+                        if _intersect(old_positions[p], positions[p], p0, p1, p2)[0]:
+                            kept_parts[p] = 0
+
+    N = int(sum(kept_parts))
+    p = 0
+    new_velocities = np.empty((N, 3))
+    new_positions = np.empty((N, 3))
+    new_old_positions = np.empty((N, 3))
+
+    for i in range(positions.shape[0]):
+        if kept_parts[i] == 1:
+            new_velocities[p] = velocities[i]
+            new_positions[p] = positions[i]
+            new_old_positions[p] = old_positions[p]
+            p += 1
+        else:
+            continue
+
+    return (new_velocities, new_positions, new_old_positions)
+
+@njit
+def _get_boundary(boundary):
+    if boundary == "xmin":
+        return (0, 0)
+    elif boundary == "xmax":
+        return (0, 1)
+    elif boundary == "ymin":
+        return (1, 0)
+    elif boundary == "ymax":
+        return (1, 1)
+    elif boundary == "zmin":
+        return (2, 0)
+    elif boundary == "zmax":
+        return (2, 1)
+
+@njit
+def _get_bc_type(bc_type):
+    if bc_type == "ela":
+        return 0
+    elif bc_type == "open":
+        return 1
+    elif bc_type == "inflow":
+        return 2
+
+class Boundary:
+    def __init__(self):
+        self.T = np.ones((3, 2))*300.0
+        self.n = np.ones((3, 2))*1e+18
+        self.u = np.zeros((3, 2, 3))
+        self.boundary_conds = np.array([[0, 0], [0, 0], [0, 0]], dtype=np.uint) # 0 = ela, 1 = open, 2 = inflow 
+        self.domain = None
+
+    def boundary(self, velocities, positions, old_positions):        
+        return _boundary(velocities, positions, old_positions, self.domain, self.boundary_conds)
+
+    def set_bc_type(self, boundary, bc_type):
+        bound = _get_boundary(boundary)
+        bc = _get_bc_type(bc_type)
+
+        self.boundary_conds[bound[0]][bound[1]] = bc
+
+        print("boundary [" + boundary + "] set to [" + bc_type + "]")
+
+    def set_bc_values(self, boundary, T, n, u):
+        i, j = _get_boundary(boundary)
+
+        self.T[i][j] = T
+        self.n[i][j] = n
+        self.u[i][j] = u
+
+        print("boundary [" + boundary + "] set to values T : {}, n : {}, u : {}".format(T, n, u))

dsmc/common.py

@@ -1,4 +1,5 @@
 from numba import cfunc, njit
+import numpy.typing as npt
 
 @cfunc("double(double[::3, ::2])")
 def get_V(a):
@@ -8,4 +9,12 @@ def get_V(a):
 def swap(arr, pos1, pos2):
     temp = arr[pos2]
     arr[pos2] = arr[pos1]
-    arr[pos1] = temp
+    arr[pos1] = temp
+
+@njit
+def is_inside(position : npt.NDArray, box : npt.NDArray) -> bool:
+    a : bool = position[0] >= box[0][0] and position[0] <= box[0][1]
+    b : bool = position[1] >= box[1][0] and position[1] <= box[1][1]
+    c : bool = position[2] >= box[2][0] and position[2] <= box[2][1]
+
+    return a and b and c