0
$\begingroup$

I have been working on a Particle-In-Cell Framework in Python and have noticed an issue where the magnitude of a EM pulse increasing infinitely as the simulation updates.

Currently, I am using the Non-Standard Finite Difference algorithm outlined in this paper: J.L.Vay 2008.

Here is the algorithm used to advance both the electric and magnetic field in vacuum:

@n.njit(parallel=True)
def E_Update(E, B, J, ds, dt, R):
    C = 2.99792458e8
    EPSILON = 8.8541878128e-12
    MU = 1.25663706212e-6
    for i in n.prange(1,E.shape[0]-1):
        for j in n.prange(1,E.shape[1]-1):
            for k in n.prange(1,E.shape[2]-1):
                # x-component update
                # -> Solve DyBz
                Bz_y0 = B[i,j-1,k,2]
                Bz_y1 = B[i,j,k,2]
                DyBz = (Bz_y1-Bz_y0)/ds
                # -> Solve DzBy
                By_z0 = B[i,j,k-1,1]
                By_z1 = B[i,j,k,1]
                DzBy = (By_z1-By_z0)/ds
                # -> Solve Ex
                R[i,j,k,0] = (1/EPSILON)*((DyBz-DzBy) - J[i,j,k,0])*dt + E[i,j,k,0]

                # y-component update
                # -> Solve DxBz
                Bz_x0 = B[i-1,j,k,2]
                Bz_x1 = B[i,j,k,2]
                DxBz = (Bz_x1-Bz_x0)/ds
                # -> Solve DzBx 
                Bx_z0 = B[i,j,k-1,0]
                Bx_z1 = B[i,j,k,0]
                DzBx = (Bx_z1-Bx_z0)/ds
                # -> Solve Ey
                R[i,j,k,1] = (1/EPSILON)*((-1)*(DyBz-DzBy) - J[i,j,k,1])*dt + E[i,j,k,1]
                
                # z-component update
                # -> Solve DxBy
                By_x0 = B[i-1,j,k,1]
                By_x1 = B[i,j,k,1]
                DxBy = (By_x1-By_x0)/ds
                # -> Solve DyBx
                Bx_y0 = B[i,j-1,k,0]
                Bx_y1 = B[i,j,k,0]
                DyBx = (Bx_y1-Bx_y0)/ds
                # -> Solve Ez
                R[i,j,k,2] = (1/EPSILON)*((DxBy-DyBx) - J[i,j,k,2])*dt + E[i,j,k,2]

@n.njit(parallel=True)
def B_Update(E, B, ds, dt, a, b, g, R):
    C = 2.99792458e8
    EPSILON = 8.8541878128e-12
    MU = 1.25663706212e-6
    C1 = ((1-(dt/(2*MU)))/(1+(dt/(2*MU))))
    C2 = (1/(1+(dt/(2*MU))))
    for i in n.prange(1,E.shape[0]-1):
        for j in n.prange(1,E.shape[1]-1):
            for k in n.prange(1,E.shape[2]-1):
                
                # x-component update
                # -> Solve DzEy 
                Ey_z0 = E[i,j,k,1]
                Ey_z1 = E[i,j,k+1,1] 
                DzEy = (Ey_z1-Ey_z0)/ds 
                # -> Solve S1zEy
                S1zEyQ1 = E[i+1,j,k,1] # Ey(i+1,j,k)
                S1zEyQ2 = E[i-1,j,k,1] # Ey(i-1,j,k)
                S1zEyQ3 = E[i,j+1,k,1] # Ey(i,j+1,k)
                S1zEyQ4 = E[i,j-1,k,1] # Ey(i,j-1,k)
                S1zEy = S1zEyQ1 + S1zEyQ2 + S1zEyQ3 + S1zEyQ4
                # -> Solve S2zEy
                S2zEyQ1 = E[i+1,j+1,k,1] # Ey[i+1,j+1,k]
                S2zEyQ2 = E[i+1,j-1,k,1] # Ey[i+1,j-1,k]
                S2zEyQ3 = E[i-1,j+1,k,1] # Ey[i-1,j+1,k]
                S2zEyQ4 = E[i-1,j-1,k,1] # Ey[i-1,j-1,k]
                S2zEy = S2zEyQ1 + S2zEyQ2 + S2zEyQ3 + S2zEyQ4
                # -> Solve DSzEy
                DSzEy = (a+b*S1zEy+g*S2zEy)*(DzEy)
                # -> Solve DyEz
                Ez_y0 = E[i,j,k,2]
                Ez_y1 = E[i,j+1,k,2]
                DyEz = (Ez_y1-Ez_y0)/ds
                # -> Solve S1yEz
                S1yEzQ1 = E[i+1,j,k,2] # Ez(i+1,j,k)
                S1yEzQ2 = E[i-1,j,k,2] # Ez(i-1,j,k)
                S1yEzQ3 = E[i,j,k+1,2] # Ez(i,j,k+1)
                S1yEzQ4 = E[i,j,k-1,2] # Ez(i,j,k-1)
                S1yEz = S1yEzQ1 + S1yEzQ2 + S1yEzQ3 + S1yEzQ4
                # -> Solve S2yEz
                S2yEzQ1 = E[i+1,j,k+1,2] # Ez(i+1,j,k+1)
                S2yEzQ2 = E[i+1,j,k-1,2] # Ez(i+1,j,k-1)
                S2yEzQ3 = E[i-1,j,k+1,2] # Ez(i-1,j,k+1)
                S2yEzQ4 = E[i-1,j,k-1,2] # Ez(i-1,j,k-1)
                S2yEz = S2yEzQ1 + S2yEzQ2 + S2yEzQ3 + S2yEzQ4
                # -> Solve DSyEz
                DSyEz = (a+b*S1yEz+g*S2yEz)*(DyEz)
                # -> Solve DtBx
                R[i,j,k,0] = (1/MU)*((DSzEy-DSyEz))*dt + B[i,j,k,0]
                
                # y-component update
                # -> Solve DzEx
                Ex_z0 = E[i,j,k,0]
                Ex_z1 = E[i,j,k+1,0]
                DzEx = (Ex_z1-Ex_z0)/ds
                # -> Solve S1zEx
                S1zExQ1 = E[i+1,j,k,0] # Ex(i+1,j,k)
                S1zExQ2 = E[i-1,j,k,0] # Ex(i-1,j,k)
                S1zExQ3 = E[i,j+1,k,0] # Ex(i,j+1,k)
                S1zExQ4 = E[i,j-1,k,0] # Ex(i,j-1,k)
                S1zEx = S1zExQ1 + S1zExQ2 + S1zExQ3 + S1zExQ4
                # -> Solve S2zEx
                S2zExQ1 = E[i+1,j+1,k,0] # Ex(i+1,j+1,k)
                S2zExQ2 = E[i+1,j-1,k,0] # Ex(i+1,j-1,k)
                S2zExQ3 = E[i-1,j+1,k,0] # Ex(i-1,j+1,k)
                S2zExQ4 = E[i-1,j-1,k,0] # Ex(i-1,j-1,k)
                S2zEx = S2zExQ1 + S2zExQ2 + S2zExQ3 + S2zExQ4
                # -> Solve DSzEx
                DSzEx = (a+b*S1zEx+g*S2zEx)*(DzEx)
                # -> Solve DxEz
                Ez_x0 = E[i,j,k,2]
                Ez_x1 = E[i+1,j,k,2]
                DxEz = (Ez_x1-Ez_x0)/ds
                # -> Solve S1xEz
                S1xEzQ1 = E[i,j+1,k,2] # Ez(i,j+1,k)
                S1xEzQ2 = E[i,j-1,k,2] # Ez(i,j-1,k)
                S1xEzQ3 = E[i,j,k+1,2] # Ez(i,j,k+1)
                S1xEzQ4 = E[i,j,k-1,2] # Ez(i,j,k-1)
                S1xEz = S1xEzQ1 + S1xEzQ2 + S1xEzQ3 + S1xEzQ4
                # -> Solve S2xEz
                S2xEzQ1 = E[i,j+1,k+1,2] # Ez(i,j+1,k+1)
                S2xEzQ2 = E[i,j+1,k-1,2] # Ez(i,j+1,k-1)
                S2xEzQ3 = E[i,j-1,k+1,2] # Ez(i,j-1,k+1)
                S2xEzQ4 = E[i,j-1,k-1,2] # Ez(i,j-1,k-1)
                S2xEz = S2xEzQ1+S2xEzQ2+S2xEzQ3+S2xEzQ4
                # -> Solve DSxEz
                DSxEz = (a+b*S1xEz+g*S2xEz)*(DxEz)
                # -> Solve DtBy
                R[i,j,k,1] = (1/MU)*((DSxEz-DSzEx))*dt + B[i,j,k,1]
                
                # z-component update
                # -> Solve DyEx
                Ex_y0 = E[i,j,k,0]
                Ex_y1 = E[i,j+1,k,0]
                DyEx = (Ex_y1-Ex_y0)/ds
                # -> Solve S1yEx
                S1yExQ1 = E[i+1,j,k,0] # Ex(i+1,j,k)
                S1yExQ2 = E[i-1,j,k,0] # Ex(i-1,j,k)
                S1yExQ3 = E[i,j,k+1,0] # Ex(i,j,k+1)
                S1yExQ4 = E[i,j,k-1,0] # Ex(i,j,k-1)
                S1yEx = S1yExQ1 + S1yExQ2 + S1yExQ3 + S1yExQ4
                # -> Solve S2yEx
                S2yExQ1 = E[i+1,j,k+1,0] # Ex(i+1,j,k+1)
                S2yExQ2 = E[i+1,j,k-1,0] # Ex(i+1,j,k-1)
                S2yExQ3 = E[i-1,j,k+1,0] # Ex(i-1,j,k+1)
                S2yExQ4 = E[i-1,j,k-1,0] # Ex(i-1,j,k-1)
                S2yEx = S2yExQ1 + S2yExQ2 + S2yExQ3 + S2yExQ4
                # -> Solve DSyEx
                DSyEx = (a+b*S1yEx+g*S2yEx)*(DyEx)
                # -> Solve DxEy
                Ey_x0 = E[i,j,k,1]
                Ey_x1 = E[i+1,j,k,1]
                DxEy = (Ey_x1-Ey_x0)/ds
                # -> Solve S1xEy
                S1xEyQ1 = E[i,j+1,k,1] # Ey(i,j+1,k)
                S1xEyQ2 = E[i,j-1,k,1] # Ey(i,j-1,k)
                S1xEyQ3 = E[i,j,k+1,1] # Ey(i,j,k+1)
                S1xEyQ4 = E[i,j,k-1,1] # Ey(i,j,k-1)
                S1xEy = S1xEyQ1 + S1xEyQ2 + S1xEyQ3 + S1xEyQ4
                # -> Solve S2xEy
                S2xEyQ1 = E[i,j+1,k+1,1] # Ey(i,j+1,k+1)
                S2xEyQ2 = E[i,j+1,k-1,1] # Ey(i,j+1,k-1)
                S2xEyQ3 = E[i,j-1,k+1,1] # Ey(i,j-1,k+1)
                S2xEyQ4 = E[i,j-1,k-1,1] # Ey(i,j-1,k-1)
                S2xEy = S2xEyQ1 + S2xEyQ2 + S2xEyQ3 + S2xEyQ4
                # -> Solve DSxEy
                DSxEy = (a+b*S1xEy+g*S2xEy)*(DxEy)
                # -> Solve DtBz
                R[i,j,k,2] = (1/MU)*((DSyEx-DSxEy))*dt + B[i,j,k,2]

Each update step involves advancing the B-Field by $\Delta t / 2$, then the E-Field by $\Delta t$, and finally the B-Field once more by $\Delta t / 2$.

So far, I've tried to resolve the issue by:

  1. Reverting to standard FDTD (i.e. setting a,b, and g in the B_Update method to 1,0, and 0 respectively)
  2. Ensuring the CFL Condition is met by keeping the time-step several orders of magnitude lower than traversal time of light across the length of a Yee-Cell.
  3. Using natural units for the simulation ($c = \mu = \epsilon = 1)$

The only resource I've found that touches on the issue of FDTD simulations "exploding" is a YouTube lecture about the basics of the algorithm. The stability of the FDTD method depends on the existence of certain values used in the update equations at the proper time and spatial coordinates. I could be wrong, but I believe that the update functions I've written meet this requirement for stability.

Any suggested changes or resource recommendations would be greatly appreciated.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.