Manual of the SymPIC

1 Introduction

A configuration file of SymPIC is a scheme (a dialect of lisp) source code. Before the PIC iteration, the configuration file will be firstly loaded and evaluated by a built-in scheme interpreter, and configuration parameters are directly read from variables in evaluated results. For more information of the scheme or lisp language, please refer to textbooks such as The Little Schemer , Structure and Interpretation of Computer Programs or ANSI Common Lisp .

2 Environment Setup

The only one environment variable for the SymPIC to run is STDLIB, which should be the location of the stdlib.scm originally in the SymPIC directory.

$ export STDLIB=path-to-sympic-dir/stdlib.scm
$ path-to-SymPIC-dir/comp config.ss

3 Grid Setup and Parallel Subdomain Decomposition

The basic calculation unit in the PIC iteration is decided by the shape of XMAX*YMAX*ZMAX, for example,

(begin
  (define XMAX 4)
  (define YMAX 8)
  (define ZMAX 3))

defines the calculation unit which contains 4x8x3 grids.

Using the Hilbert space filling curve , the computational domain can be expanded. For example,

(begin
  (define NUM_N_HILBERT_DIMENSION 2)
  (define NUM_N_HILBERT 3)
  (define HILBERT_DIR 1))

defines a 2-dimensional 3rd-order Hilbert curve, and its normal vector is pointed to 1 (y) direction. So in this case the Hilbert curve is 8x1x8, and the total calculation domain is 32x8x24.

Note that if you need 2 or 1-dimensional simulation, then shape for both calculation unit and Hilbert curve should also be 2 or 1-dimensional. For example,

(begin
  (define XMAX 6)
  (define YMAX 5)
  (define ZMAX 1)
  (define NUM_N_HILBERT_DIMENSION 2)
  (define NUM_N_HILBERT 2)
  (define HILBERT_DIR 2))

will set the total calculation domain to 24x20x1,

(begin
  (define XMAX 8)
  (define YMAX 1)
  (define ZMAX 1)
  (define NUM_N_HILBERT_DIMENSION 1)
  (define NUM_N_HILBERT 3)
  (define HILBERT_DIR 0))

will set the total calculation domain to 64x1x1. If in one direction has multiple grids, then the shape of calculation unit in this direction should be at least set to 2.

One compute unit (one CPU core or GPU SIMD unit) will be assigned integer multiple of calculation unit task(s), so the total number of compute units should be smaller than (2^NUM_N_HILBERT)^NUM_N_HILBERT_DIMENSION to avoid wasting.

4 Plasma Parameters

4.1 Time Setup

The timestep and total number of timesteps can be set via DELTAT and NUM_TIMESTEP.

4.2 Initialization of Particles

The number of particle species can be specified by the NUM_SPEC variable. For example,

(define NUM_SPEC 2)

will use 2 species of particles.

Mass and charge of each species of sampling points of particles can be set via GET_MASS and GET_CHARGE function, for example,

(begin
  (define NPM 100)
  
(define (GET_MASS i) (case i (0 (* 1 NPM)) (1 (* 1836 NPM)) (else 1)))

  
(define (GET_MASS i) (case i (0 (* -1 NPM E0)) (1 (* NPM E0)) (else 1)))
)

For velocity space, currently Maxwellian and lost cone distribution are supported. The GET_INIT_VT function can specify the thermal velocity of each particle species. For example,

(define GET_INIT_VT
(lambda (i)
  (case i (0 1.00000000000000006e-01) (1 4.00000000000000008e-03) (else 0)))
)

will set the thermal velocities of 0th and 1st particle to 0.1 and 0.004, respectively. To use non-uniform temperature, firstly USE_NON_UNI_TEMPERATURE should be set to 1, and the use GET_INIT_TEMPERATURE_DIST to specify the temperature dist. For example,

(begin
  (define USE_NON_UNI_TEMPERATURE 1)
  (define GET_INIT_TEMPERATURE_DIST
(lambda (i z y x l)
  (case i (0 (* (+ 1 (cos x)) (case l (0 5.00000000000000000e-01) (1 1.19999999999999996e+00) (2 1.00000000000000000e+00) (else 0)))) (1 (- 1 (sin (/ y 2)))) (else 0)))
))

sets the thermal speed in x, y, z directions of the 0th particle to (1+cos(x))*0.5*GET_INIT_VT(0), (1+cos(x))*1.2*GET_INIT_VT(0), (1+cos(x))*1.0*GET_INIT_VT(0), respectively, sets the thermal speed of the 1st particle to (1-sin(y/2))*GET_INIT_VT(1). The GET_INIT_LOAD function will set number of sampling points per grid. For example,

(define GET_INIT_LOAD
(lambda (i)
  (case i (0 100) (1 50)))
)

sets the number of sampling points per grid of 0th and 1st particle to 100 and 50, respectively.

Non-uniform density can be set via USE_NON_UNI_DENSITY and GET_INIT_DENSITY_DIST. For example,

(begin
  (define USE_NON_UNI_DENSITY 1)
  (define GET_INIT_DENSITY_DIST
(lambda (i z y x)
  (case i (0 (exp (- (+ (* x x) (* y y))))) (1 (exp (- (+ (* z z) (* y y))))) (else 0)))
))
the number of sampling points per grid of 0th and 1st particle to GET_INIT_LOAD(0)*exp(-(x*x+y*y)) and GET_INIT_LOAD(1)*exp(-(y*y+z*z)), respectively.

4.3 Initialization of Fields

Following parameters are used for setting initial and external fields. For example for initial E and B,

(begin
  (define USE_INIT_EB0 1)
  (define GET_INIT_E0
(lambda (z y x l)
  (case l (0 1) (else 0)))
)
  (define GET_INIT_B0
(lambda (z y x l)
  (case l (1 1) (else 0)))
))

will set the initial E and B to [1,0,0] and [0,1,0], respectively. For external E and B, the following code

(begin
  (define USE_INIT_EXT_EB 1)
  (define GET_INIT_E
(lambda (z y x l)
  (case l (1 (log z)) (else 0)))
)
  (define GET_INIT_B
(lambda (z y x l)
  (case l (2 2) (else 0)))
))

will set the external E and B to [0,log(z),0] and [0,0,2]

4.4 Unit Normalization

We may now introduce the unit normalization. The unit constraint in the SymPIC is the speed of light and the grid size, which are set to 1. So if the real grid size is \(\Delta x\), and the speed of light is \(\mathrm{c}\), then the normalized time is \(\Delta x/\mathrm{c}\). Usually the time step DELTAT is set to 0.5, which means the real time step is \(0.5\Delta x/\mathrm{c}\), to comply the CFL constraint, however if the gyro-period or the plasma frequency of electrons are too large then the DELTAT should be smaller enough to resolve the fastest physics.

Suppose the charge and mass of the \(s\)-th species are \(q_s=\)GET_CHARGE(\(s\)) and \(m_s=\)GET_MASS(\(s\)), then in the SymPIC code, the motion equation for the \(p\)-th particle of the \(s\)-th species is $$\ddot{\mathbf{x}}_{sp}=\frac{q_s}{m_s}\left(\mathbf{E}(\mathbf{x}_{sp})+\dot{\mathbf{x}}_{sp}\times\mathbf{B}(\mathbf{x}_{sp})\right)$$ and the Maxwell's equations are $$\begin{eqnarray*} \dot{\mathbf{B}}&=&-\nabla\times\mathbf{E}\\ \dot{\mathbf{E}}&=&\nabla\times\mathbf{B}-\sum_\mathrm{s,p}q_s\dot {\mathbf{x}}_{sp}\delta (\mathbf{x}-\mathbf{x}_{sp}) \end{eqnarray*}$$

In practice, we may set the rest mass of the electrons to 1, and in this case all units are fixed. Some units of physical quantities are shown in the following table, where \(-q_e\) and \(m_e\) are the charge and rest mass of electrons, \(\mu_0\) is the permeability of the vacuum.

Physical quantity Unit
Number Density \(\Delta x^{-3}\)
Charge \(\sqrt{\frac{m_e\Delta x}{\mu_0}}\)
Magnetic Field \(\mathrm{c}\sqrt{\frac{\mu_0m_e}{\Delta x^3}}\)
Electric Field \(\mathrm{c}^2\sqrt{\frac{\mu_0m_e}{\Delta x^3}}\)
Current Density \(\mathrm{c}\sqrt{\frac{m_e}{\mu_0\Delta x^5}}\)
Current \(\mathrm{c}\sqrt{\frac{m_e}{\mu_0\Delta x}}\)

4.5 Boundary Conditions

By default, periodic boundary in all x, y, z directions are used. To use PEC and ABC boundaries, USE_NP_BOUNDARY should be set to 1. Use USE_PEC_DIR and USE_ABC_DIR to specify the directions of PEC and ABC boundary needed. Their effects are shown as follows.

Value of USE_PEC_DIR/USE_ABC_DIR 0 1 2 3 4 5 6 7
Applied direction(s) None x y x and y z x and z y and z x, y and z

5 Source Setup

To be added.

6 Parallel Runtime Setup

Use the NUM_RUNTIME parameter to specify the number of runtimes in one MPI process that are used, and GET_DEV_TYPE and GET_DEV_ID are used to specify the parallel device type and id. For example, if the compiled code contains only 1 kind of parallel runtimes, No. 0 runtime is OpenMP, then

(begin
  (define NUM_RUNTIME 1)
  
(define (GET_DEV_TYPE x_runtime rank) 0)

  
(define (GET_DEV_ID x_runtime rank) 0)
)

will set each MPI process only use 1 runtime which is OpenMP.

If the compiled code contains 2 types of parallel runtimes, No. 0 runtime is CUDA, and No. 1 is OpenMP, then

(begin
  (define NUM_RUNTIME 3)
  (define GET_DEV_TYPE
(lambda (x_runtime rank)
  (case x_runtime (0 0) (1 0) (2 1) (else 0)))
)
  (define GET_DEV_ID
(lambda (x_runtime rank)
  (define rk_rem (* 2 (remainder rank 2)))
  
(if (< x_runtime 2)
  (+ rk_rem x_runtime)
  0)
)
))

will set the even number of MPI process using No. 0, 1 CUDA GPU and OpenMP, odd number of MPI process using No. 2, 3 CUDA GPU and OpenMP. Note that for multi-GPU system, if we assign more than one CUDA GPU in one process, due to the low efficiency of CUDA runtime this could harm the performance, however using more MPI processes will cause more communications that will also harm the performance. So some trade-off will need to be considered.

7 Output Setup

Electromagnetic fields as well as some density profiles are outputed every NUM_DUMP_TIMESTEP time steps.

8 Checkpoint Setup

To use checkpoint, USE_CHECKPOINT should be set to 1. In this case checkpoint will be outputed every NUM_CHECKPOINT_TIMESTEP time steps. The program will automatically search available checkpoint when starting.

9 Algorithm Selection

The SymPIC is designed to support multiple types of PIC schemes. Currently it support the 2nd-order non-relativistic and the 1st-order relativistic electromagnetic symplectic schemes, and the 2nd-order non-relativistic adiabatic electron-kinetic ion model. The default algorithm is the 2nd-order explicit charge-conservative electromagnetic PIC scheme. To use the 1st-order relativistic electromagnetic PIC scheme, the USE_REL should be set to 1. To use the 2nd-order low frequency adiabatic electron-kinetic ion model, USE_ITG should be set to 1.

10 Examples

1. Landau Damping (CN)
2. Dispersion Relation of Ion Bernstein Waves