An example code for axisymmetric FDTD

unstable-fdtd.jl

ε_0    = 8.85418781e-12  # vacuum permittivity [F/m]
μ_0    = 1.256637062e-6  # vacuum permeability [H/m]
c      = 299_792_458.    # speed of light [m/s]
NR     = 101             # number of grid points along radial direction [1]
NZ     = 101             # number of grid points along axial direction [1]
RADIUS = 0.010           # radius along r-axis [m]
LENGTH = 0.010           # lenght along z-axis [m]

DZ     = LENGTH / (NZ-1) # cell-size along axial direction [m]
DR     = RADIUS / (NR-1) # cell-size along radial direction [m]
DT     = DZ / √2c        # time-step at Courant limit (assumes DR = DZ) [s] 

B_φ    =  zeros(Float64, NR-1, NZ-1) # magnetic field [A/m]
E_z    =  zeros(Float64, NR,   NZ-1) # electric field [V/m]
E_r    =  zeros(Float64, NR-1, NZ)   # electric field [V/m]
J_z    =  zeros(Float64, NR,   NZ-1) # current density [A/m^2]
J_r    =  zeros(Float64, NR-1, NZ)   # current density [A/m^2]
R      =  zeros(Float64, NR)         # radial coordinate of nodes [m]
S      =  zeros(Float64, NR)         # surface area in axial direction [m^2]

R .= [(j-1) * DR for j=1:NR]
S .= [π * ((j-0.5) * DR)^2 - π * ((j-1.5) * DR)^2 for j=1:NR]
S[1] = π * (0.5*DR)^2

gauss_pulse(x, σ, μ) = (1/σ/sqrt(2π)) * exp(-0.5(x - μ)^2 / σ^2)
frequency = 300e9
period   = 1.0 / frequency

function boundary_conditions!()
    for k in 1:NZ-1
        E_z[NR, k] = 0.0 # perfect electric conductor at r=RADIUS
    end
    
    for j in 1:NR-1
        E_r[j, 1]   = E_r[j,   2] # perfect magnetic conductor at z=0
        E_r[j,NZ]   = E_r[j,NZ-1] # perfect magnetic conductor at z=LENGTH
    end
    
    return nothing
end

function update_magnetic!()
    for k = 1:NZ-1
        for j = 1:NR-1
            B_φ[j,k] += DT * (E_z[j+1,k] - E_z[j,k]) / DR
            B_φ[j,k] -= DT * (E_r[j,k+1] - E_r[j,k]) / DZ
        end
    end
    
    return nothing
end

function update_electric!()
    for k = 1:NZ-1
        for j = 2:NR-1
            d1 = 1.0 + 0.5 / (j-1)
            d2 = 1.0 - 0.5 / (j-1)
            E_z[j,k] = E_z[j,k] -
                        c^2 * DT * μ_0 * J_z[j,k] +
                        c^2 * DT * (d1 * B_φ[j,k] - d2 * B_φ[j-1,k]) / DR
        end
    end

    for k = 1:NZ-1
        d0 = 4.0 / DR # the algorithm is stable for d0 = 2.0 / DR
        E_z[1,k] = E_z[1,k] -
                    c^2 * DT * μ_0 * J_z[1,k] +
                    c^2 * DT * (d0 * B_φ[1,k])
    end

    for k in 2:NZ-1
        for j in 1:NR-1
            E_r[j,k] = E_r[j,k] - 
                        c^2 * DT * μ_0 * J_r[j,k] -
                        c^2 * DT * (B_φ[j,k] - B_φ[j,k-1]) / DZ
        end
    end
    
    return nothing
end

for t=0:DT:20period
    J_z[10, 10:90] .= gauss_pulse(t + 0.5DT, 0.5period, 2.0period)
    J_z ./= S

    update_magnetic!()
    update_electric!()
    boundary_conditions!()
end

@show extrema(E_z)