From eca3e7b1c2b3b9714c577df6108c0ffaff646b5a Mon Sep 17 00:00:00 2001
From: JiaoLuhuai <103885549+JiaoLuhuai@users.noreply.github.com>
Date: Tue, 27 Sep 2022 16:51:36 +0800
Subject: [PATCH] Add files via upload

---
 pic88.py | 58 ++++++++++++++++++++++++++++++--------------------------
 1 file changed, 31 insertions(+), 27 deletions(-)

diff --git a/pic88.py b/pic88.py
index 32cbe8a..b0b8cd6 100644
--- a/pic88.py
+++ b/pic88.py
@@ -1,26 +1,30 @@
 # Authored by Luhuai Jiao
+# This is a 1D simulation of "two-stream instability" in Plasma Physicis. 
+# Some settings of grids and particles are from "Introduction to Computational Plasma Physics"(《计算等离子体物理导论》，谢华生著)
+
 import taichi as ti
 
-ti.init(arch=ti.gpu)
+ti.init(arch=ti.gpu)  # Try to run on GPU
 PI = 3.141592653589793
-L = 8 * PI
-dt = 0.1
-substepping = 4
-ng = 32
-np = 16384
-vb = 1.0
-vt = 0.3
-wp = 1  #Plasma frequence
-qm = -1
-q = wp * wp / (qm * np / L)
-rho_back = -q * np / L  #background charge
-dx = L / ng
+L = 8 * PI  # simulation domain size
+dt = 0.1  # time step
+substepping = 8
+ng = 32  # number of grids
+np = 16384  # numer of particles
+vb = 1.0  # beam-velocity, one is vb, the other is -vb 
+vt = 0.3  # thermal velocity 
+wp = 1  # Plasma frequence
+qm = -1  # charge-mass ratio
+q = wp * wp / (qm * np / L)  # charge of a particle
+rho_back = -q * np / L  # background charge density
+dx = L / ng  # grid spacing
 inv_dx = 1.0 / dx
-x = ti.Vector.field(1, ti.f32, np)
-v = ti.Vector.field(1, ti.f32, np)
-rho = ti.Vector.field(1, ti.f32, ng)
-e = ti.Vector.field(1, ti.f32, ng)
-v_x_pos1 = ti.Vector.field(2, ti.f32, int(np / 2))  #to plot Vx-X
+x = ti.Vector.field(1, ti.f32, np)  # position
+v = ti.Vector.field(1, ti.f32, np)  # velocity
+rho = ti.Vector.field(1, ti.f32, ng)  # charge density
+e = ti.Vector.field(1, ti.f32, ng)  # electric fields
+# to show x-vx on the screen
+v_x_pos1 = ti.Vector.field(2, ti.f32, int(np / 2))  
 v_x_pos2 = ti.Vector.field(2, ti.f32, int(np / 2))
 
 
@@ -28,41 +32,41 @@
 def initialize():
     for p in x:
         x[p].x = (p + 1) * L / np
-        v[p].x = vt * ti.randn() + (-1)**p * vb
+        v[p].x = vt * ti.randn() + (-1)**p * vb  # two streams
 
 
 @ti.kernel
 def substep():
     for p in x:
         x[p] += v[p] * dt
-        if x[p].x >= L:
+        if x[p].x >= L:  # periodic boundary condition
             x[p] += -L
         if x[p].x < 0:
             x[p] += L
-    rho.fill(rho_back)
-    for p in x:
+    rho.fill(rho_back)  # fill rho with background charge density
+    for p in x:  # Particle state update and scatter to grid (P2G)
         base = (x[p] * inv_dx - 0.5).cast(int)
         fx = x[p] * inv_dx - 0.5 - base.cast(float)
         rho[base] += (1.0 - fx) * q * inv_dx
         rho[base + 1] += fx * q * inv_dx
     e.fill(0.0)
     ti.loop_config(serialize=True)
-    for i in range(ng):
+    for i in range(ng):  # compute electric fields
         e[i] = e[i - 1] + (rho[i - 1] + rho[i]) * dx * 0.5
     s = 0.0
     for i in e:
         s += e[i].x
     for i in e:
         e[i] += -s / ng
-    for p in v:
+    for p in v:  #G2P
         base = (x[p] * inv_dx - 0.5).cast(int)
         fx = x[p] * inv_dx - 0.5 - base.cast(float)
-        a = (e[base] * (1.0 - fx) + e[base + 1] * fx) * qm
+        a = (e[base] * (1.0 - fx) + e[base + 1] * fx) * qm  # compute electric force
         v[p] += a * dt
 
 
 @ti.kernel
-def vx_pos():
+def vx_pos():  # to show x-vx on the screen
     for p in x:
         if p % 2:
             v_x_pos1[int((p - 1) / 2)].x = x[p].x / L
@@ -74,7 +78,7 @@ def vx_pos():
 
 def main():
     initialize()
-    gui = ti.GUI("pic88", (800, 800))
+    gui = ti.GUI("Shortest PIC", (800,800))
     while not gui.get_event(ti.GUI.ESCAPE, ti.GUI.EXIT):
         for s in range(substepping):
             substep()