Plotter

CHANGELOG.md

@@ -6,7 +6,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
-## [0.3.0] - 2022-09-23
+## [0.5.0] - 2022-09-23
+### Added
+- writer module
+
+## [0.4.0] - 2022-09-23
 ### Added
 - dsmc module
 

dsmc/diagnostics.py

@@ -0,0 +1,60 @@
+import numpy as np
+from numba import njit
+from . import particles as prt
+from . import octree as oc
+
+def sort_bin(positions, axis, Nbin):
+    bins = [[] for _ in range(Nbin)]
+    b = 0
+    box = oc._find_bounding_box(positions)
+    dx = (box[axis][1] - box[axis][0]) / (Nbin - 1)
+    xx = [dx]
+    x = dx
+    sub_pos = np.array([pos[axis] for pos in positions])
+    sorted_pos = np.argsort(sub_pos)
+
+    for i in range(len(sorted_pos)):
+        p = sorted_pos[i]
+        while positions[p][axis] > x:
+            x += dx
+            b += 1
+            xx.append(x)
+
+        bins[b].append(p)
+
+    return bins, box, xx
+
+def calc_n(bins, box, axis, w):
+    Nbins = len(bins)
+    dx = (box[axis][1] - box[axis][0]) / Nbins
+    V = 1
+    n = np.empty((Nbins, ))
+    for i in range(3):
+        if i == axis:
+            V *= (box[i][1] - box[i][0]) / Nbins
+        else:
+            V *= (box[i][1] - box[i][0])
+
+    for i in range(Nbins):
+        n[i] = len(bins[i]) * w / V
+
+    return n
+
+def calc_T(bins, velocities, mass):
+    Nbins = len(bins)
+    T = np.empty((Nbins, ))
+
+    for i in range(Nbins):
+        vels = np.array([velocities[p] for p in bins[i]])
+        T[i] = prt.calc_temperature(vels, mass)
+
+    return T
+
+def calc_p(n, T):
+    Nbins = len(n)
+    p = np.empty((Nbins, ))
+
+    for i in range(Nbins):
+        p[i] = n[i]*T[i]*prt.kb
+
+    return p

dsmc/dsmc.py

@@ -126,11 +126,21 @@ 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):
-        N = int(round(n / self.w))
+    def create_particles(self, box, T, n, u = None):
+        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)
+
         print("now containing {} particles, {} total".format(N, self.particles.N))
 
     def set_domain(self, domain):

dsmc/particles.py

@@ -57,7 +57,7 @@ def get_vel(T, mass):
 
 @njit
 def get_velocities(T, mass, N):
-    velocities = np.empty((N, 3), dtype=float)
+    velocities = np.empty((N, 3))
 
     for i in range(N):
         velocities[i] = get_vel(T, mass)
@@ -88,7 +88,7 @@ def calc_positions(x, y, z, N):
     z : tuple(2), dtype = float
         zmin, zmax
     """
-    positions = np.empty((N, 3), dtype=float)
+    positions = np.empty((N, 3))
 
     for i in range(N):
         positions[i][0] = x[0] + np.random.random() * (x[1] - x[0])

dsmc/writer.py

@@ -0,0 +1,100 @@
+from . import octree as oc
+
+def write_buttom_leafs(octree):
+	f = open("octree.vtu", "w")
+
+	_write_header(f)
+	_wrtie_body(f, octree)
+	_write_footer(f)
+
+	f.close()
+
+def _write_header(f):
+	f.write("<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n")
+
+def _wrtie_body(f, octree):
+	leaf_ids = []
+	for i in range(len(octree.leafs)):
+		if octree.leafs[i].number_children == 0:
+			leaf_ids.append(i)
+
+	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);
+
+	f.write("</Piece>\n")
+	f.write("</UnstructuredGrid>\n")
+
+def _write_points(f, octree, leaf_ids):
+	f.write("<Points>\n")
+	f.write("<DataArray type=\"Float32\" NumberOfComponents=\"3\" format=\"ascii\">\n")
+
+	for i in leaf_ids:
+		box = octree.cell_boxes[i]
+
+		f.write("{} ".format(box[0][0]))
+		f.write("{} ".format(box[1][0]))
+		f.write("{} ".format(box[2][0]))
+
+		f.write("{} ".format(box[0][1]))
+		f.write("{} ".format(box[1][0]))
+		f.write("{} ".format(box[2][0]))
+
+		f.write("{} ".format(box[0][1]))
+		f.write("{} ".format(box[1][1]))
+		f.write("{} ".format(box[2][0]))
+
+		f.write("{} ".format(box[0][0]))
+		f.write("{} ".format(box[1][1]))
+		f.write("{} ".format(box[2][0]))
+
+		f.write("{} ".format(box[0][0]))
+		f.write("{} ".format(box[1][0]))
+		f.write("{} ".format(box[2][1]))
+
+		f.write("{} ".format(box[0][1]))
+		f.write("{} ".format(box[1][0]))
+		f.write("{} ".format(box[2][1]))
+
+		f.write("{} ".format(box[0][1]))
+		f.write("{} ".format(box[1][1]))
+		f.write("{} ".format(box[2][1]))
+
+		f.write("{} ".format(box[0][0]))
+		f.write("{} ".format(box[1][1]))
+		f.write("{} ".format(box[2][1]))
+
+	f.write("</DataArray>\n")
+	f.write("</Points>\n")
+
+
+def _write_cells(f, octree, leaf_ids):
+	k = 0
+
+	f.write("<Cells>\n")
+	f.write("<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">\n")
+
+	for i in range(len(leaf_ids)):
+		for _ in range(8):
+			f.write("{} ".format(k))
+			k += 1
+
+	f.write("</DataArray>\n")
+	f.write("<DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">\n")
+
+	for i in range(len(leaf_ids)):
+		f.write("{} ".format((i + 1) * 8))
+
+	f.write("</DataArray>\n")
+	f.write("<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n")
+
+	for _ in range(len(leaf_ids)):
+		f.write("12 ")
+
+	f.write("</DataArray>\n")
+	f.write("</Cells>\n")
+
+def _write_footer(f):
+	f.write("</VTKFile>\n")

examples/plot.py

@@ -1,38 +1,20 @@
 import matplotlib.pyplot as plt
 import csv
 import numpy as np
+import argparse
 
 if __name__ == '__main__':
-    res = []
+    parser = argparse.ArgumentParser(description='Process some integers.')
+    parser.add_argument('file_name', type=str)
 
-    with open('test.csv') as file:
+    args = parser.parse_args()
+
+    with open(args.file_name) as file:
         reader = csv.reader(file, delimiter=',')
-        res = []
 
         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, ))
-            dx = 0.05/N
-            x[0] = dx
-            q = 0
-
-            for i in range(len(data)):
-                while data[i] > x[q]:
-                    q += 1
-                    x[q] = x[q - 1] + dx
-
-                sor[q] += 1
-
-
-            x.resize(q)
-            sor.resize(q)
-
-            res.append((x, sor))
+            n = [float(l[i]) for i in range(len(l))]
 
-    plt.plot(res[-1][0], res[-1][1])
-    plt.show()
+        plt.plot(n)
+        plt.show()

examples/shock_tube.py

@@ -1,40 +1,65 @@
 import dsmc
+import dsmc.diagnostics as dia
+import matplotlib.pyplot as plt
+
+def write2file(solver, n_file, T_file, p_file):
+    bins, box, x = dia.sort_bin(solver.particles.Pos, 2, Nbins)
+    N = [len(b) for b in bins]
+    n = dia.calc_n(bins, box, 2, solver.w)
+    T = dia.calc_T(bins, solver.particles.Vel, mass)
+    p = dia.calc_p(n, T)
+
+    for i in range(Nbins):
+        n_file.write("{},".format(n[i]))
+        T_file.write("{},".format(T[i]))
+        p_file.write("{},".format(p[i]))
+
+    n_file.write("\n")
+    T_file.write("\n")
+    p_file.write("\n")
 
 if __name__ == '__main__':
     # general parameters
     solver = dsmc.DSMC()
-    domain = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (0, 50e-3))
+    domain = [(-0.0001, 0.0001), (-0.0001, 0.0001), (0.0, 0.1)]
     dt = 1e-7
-    w = 2.4134e+15
+    w = 1e+9
     mass = 6.6422e-26
     niter = 300
+    Nbins = 100
 
-    # low denisty particles
-    nhigh = 2.5e+20
+    # high denisty particles
+    nhigh = 2.41432e+22 
     Thigh = 300
-    Boxhigh = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (25e-3, 50e-3))
+    Boxhigh = [(-0.0001, 0.0001), (-0.0001, 0.0001), (0.0, 0.05)]
 
-    # high denisty particles
-    nlow = 2.5e+19
+    # low denisty particles
+    nlow = 2.41432e+21
     Tlow = 300
-    Boxlow = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (0, 25e-3))
+    Boxlow = [(-0.0001, 0.0001), (-0.0001, 0.0001), (0.05, 0.1)]
 
+    # setup solver
     solver.set_domain(domain)
     solver.set_weight(w)
     solver.set_mass(mass)
 
     solver.create_particles(Boxlow, Tlow, nlow)
     solver.create_particles(Boxhigh, Thigh, nhigh)
 
-    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")
+    # open files
+    n_file = open("n.csv", "w")
+    T_file = open("T.csv", "w")
+    p_file = open("p.csv", "w")
+
+    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)
 
+	# close files
+    n_file.close()
+    T_file.close()
+    p_file.close()        
+
     print("")
     print('done')

examples/temp_relax.py

@@ -0,0 +1,50 @@
+import dsmc
+import dsmc.diagnostics as dia
+import dsmc.particles as prt
+import matplotlib.pyplot as plt
+import numpy as np
+import math
+
+def maxwell(x, T):
+    return 2.0 * np.sqrt(x) * np.exp(-x/T) / (math.pow(T, 3.0/2.0) * np.sqrt(math.pi))
+
+def calc_x(velocities, mass):
+    return np.array([mass*vel.dot(vel)/(2.0*prt.kb) for vel in velocities])
+
+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-5
+    w = 2.4134e+7
+    mass = 6.6422e-26
+    niter = 500
+    Nbins = 100
+
+    # particles
+    n = 2.5e+20
+    T = 300
+    Box = ((-0.1e-3, 0.1e-3), (-0.1e-3, 0.1e-3), (25e-3, 50e-3))
+    u = np.array((1000.0, 0.0, 0.0))
+
+
+    solver.set_domain(domain)
+    solver.set_weight(w)
+    solver.set_mass(mass)
+
+    solver.create_particles(Box, T, n, u)
+
+    for it in range(niter):
+        print("iteration {:4}/{}".format(it + 1, niter), end="\r", flush=True)
+        solver.advance(dt)
+        x = calc_x(solver.particles.Vel, solver.mass)
+
+    xm = np.linspace(0, 3500, 1000)
+    dist = [maxwell(xi, 300) for xi in xm]
+
+    plt.plot(xm, dist)
+    plt.hist(x, Nbins, density=True)
+    plt.show()  
+
+    print("")
+    print('done')