So I just quickly went through and wrapped the code in a let block, removed all the const annotations and put @inbounds in front of all the loops and the first time I run the code I get
elapsed CPU time: 0.668674 seconds
elapsed real time: 0.799 seconds
and the second time (now that all the functions are compiled) I get
elapsed CPU time: 0.543445 seconds
elapsed real time: 0.551 seconds
Here’s my modified code:
using CPUTime
let
CPUtic()
start = time()
#Set constants - these are physical domain sizes for the FDTD problem (nearly a cubic space)
SizeX = 48
SizeY = 47
SizeZ = 47
MaxTime = 150 #this is the time domain size for the FDTD problem
Cdtds = 1.0/(sqrt(3.0)) #Courant number
imp0 = 377.0 #characteristic impedance of free space
abccoef = (Cdtds - 1.0)/(Cdtds + 1.0) #coefficient associated with boundary conditions
ppw = 15 #points per wave
#initialize arrays to store electric and magnetic field values in 3D space
Ex = zeros(Float64,(SizeX,SizeY,SizeZ))
Ey = zeros(Float64,(SizeX,SizeY,SizeZ))
Ez = zeros(Float64,(SizeX,SizeY,SizeZ))
Hx = zeros(Float64,(SizeX,SizeY,SizeZ))
Hy = zeros(Float64,(SizeX,SizeY,SizeZ))
Hz = zeros(Float64,(SizeX,SizeY,SizeZ))
#Define arrays for plotting results.
Eout = zeros(Float64,(SizeY,SizeZ,MaxTime)) #Ouput will be electric field on a 2D cross section
tvar = collect(1:MaxTime) #time array
#Set electric field update coefficients
Cexe = fill(1.0,(SizeX,SizeY,SizeZ))
Cexh = fill(Cdtds * imp0,(SizeX,SizeY,SizeZ))
Ceye = fill(1.0,(SizeX,SizeY,SizeZ))
Ceyh = fill(Cdtds * imp0,(SizeX,SizeY,SizeZ))
Ceze = fill(1.0,(SizeX,SizeY,SizeZ))
Cezh = fill(Cdtds * imp0,(SizeX,SizeY,SizeZ))
#Set magnetic field update coefficients
Chxh = fill(1.0,(SizeX,SizeY,SizeZ))
Chxe = fill(Cdtds / imp0,(SizeX,SizeY,SizeZ))
Chyh = fill(1.0,(SizeX,SizeY,SizeZ))
Chye = fill(Cdtds / imp0,(SizeX,SizeY,SizeZ))
Chzh = fill(1.0,(SizeX,SizeY,SizeZ))
Chze = fill(Cdtds / imp0,(SizeX,SizeY,SizeZ))
#Initialize temporary variables for absorbing boundary condition (ABC) function
Eyx0 = zeros(Float64, (SizeY-1,SizeZ))
Ezx0 = zeros(Float64, (SizeY,SizeZ-1))
Eyx1 = zeros(Float64, (SizeY-1,SizeZ))
Ezx1 = zeros(Float64, (SizeY,SizeZ-1))
Exy0 = zeros(Float64, (SizeX-1,SizeZ))
Ezy0 = zeros(Float64, (SizeX,SizeZ-1))
Exy1 = zeros(Float64, (SizeX-1,SizeZ))
Ezy1 = zeros(Float64, (SizeX,SizeZ-1))
Exz0 = zeros(Float64, (SizeX-1,SizeY))
Eyz0 = zeros(Float64, (SizeX,SizeY-1))
Exz1 = zeros(Float64, (SizeX-1,SizeY))
Eyz1 = zeros(Float64, (SizeX,SizeY-1))
#update magnetic field (H) over all space at a given time step given electric and magnetic fields from the previous time step
#achieved by nested for loops
function updateH(Ex,Ey,Ez,Hx,Hy,Hz)
@inbounds for pp = 1:SizeZ-1
for nn = 1:SizeY-1
for mm = 1:SizeX
Hx[mm,nn,pp] = (Chxh[mm,nn,pp] * Hx[mm,nn,pp] +
Chxe[mm,nn,pp] * ((Ey[mm,nn,pp+1] - Ey[mm,nn,pp]) - (Ez[mm,nn+1,pp] - Ez[mm,nn,pp])))
end;end;end
@inbounds for pp = 1:SizeZ-1
for nn = 1:SizeY
for mm = 1:SizeX-1
Hy[mm,nn,pp] = (Chyh[mm,nn,pp] * Hy[mm,nn,pp] +
Chye[mm,nn,pp] * ((Ez[mm+1,nn,pp] - Ez[mm,nn,pp]) - (Ex[mm,nn,pp+1] - Ex[mm,nn,pp])))
end;end;end
@inbounds for pp = 1:SizeZ
for nn = 1:SizeY-1
for mm = 1:SizeX-1
Hz[mm,nn,pp] = (Chzh[mm,nn,pp] * Hz[mm,nn,pp] +
Chze[mm,nn,pp] * ((Ex[mm,nn+1,pp] - Ex[mm,nn,pp]) - (Ey[mm+1,nn,pp] - Ey[mm,nn,pp])))
end;end;end
return Hx,Hy,Hz
end
#update electric field (E) over all space at a given time step given electric and magnetic fields from the previous time step
function updateE(Ex,Ey,Ez,Hx,Hy,Hz) #update electric field
@inbounds for pp = 2:SizeZ-1
for nn = 2:SizeY-1
for mm = 1:SizeX-1
Ex[mm,nn,pp] = (Cexe[mm,nn,pp] * Ex[mm,nn,pp] +
Cexh[mm,nn,pp] * ((Hz[mm,nn,pp] - Hz[mm,nn-1,pp]) - (Hy[mm,nn,pp] - Hy[mm,nn,pp-1])))
end;end;end
@inbounds for pp = 2:SizeZ-1
for nn = 1:SizeY-1
for mm = 2:SizeX-1
Ey[mm,nn,pp] = (Ceye[mm,nn,pp] * Ey[mm,nn,pp] +
Ceyh[mm,nn,pp] * ((Hx[mm,nn,pp] - Hx[mm,nn,pp-1]) - (Hz[mm,nn,pp] - Hz[mm-1,nn,pp])))
end;end;end
@inbounds for pp = 1:SizeZ-1
for nn = 2:SizeY-1
for mm = 2:SizeX-1
Ez[mm,nn,pp] = (Ceze[mm,nn,pp] * Ez[mm,nn,pp] +
Cezh[mm,nn,pp] * ((Hy[mm,nn,pp] - Hy[mm-1,nn,pp]) - (Hx[mm,nn,pp] - Hx[mm,nn-1,pp])))
end;end;end
return Ex,Ey,Ez
end
#apply absorbing boundary conditions (ABCs) at the boundary planes of all x,y,z dimensions
function abc3dfirst(Ex,Ey,Ez)
# global Eyx0, Ezx0, Eyx1, Ezx1, Exy0, Ezy0, Exy1, Ezy1, Exz0, Eyz0, Exz1, Eyz1
# =============================================================================
#ABC at x boundary plane 1
mm = 1
@inbounds for pp = 1:SizeZ
for nn = 1:SizeY-1
Ey[mm,nn,pp] = Eyx0[nn,pp] + abccoef * (Ey[mm+1,nn,pp] - Ey[mm,nn,pp])
Eyx0[nn,pp] = Ey[mm+1,nn,pp]
end;end
for pp = 1:SizeZ-1
for nn = 1:SizeY
Ez[mm,nn,pp] = Ezx0[nn,pp] + abccoef * (Ez[mm+1,nn,pp] - Ez[mm,nn,pp])
Ezx0[nn,pp] = Ez[mm+1,nn,pp]
end;end
#ABC at x boundary plane 2
mm = SizeX
@inbounds for pp = 1:SizeZ
for nn = 1:SizeY-1
Ey[mm,nn,pp] = Eyx1[nn,pp] + abccoef * (Ey[mm-1,nn,pp] - Ey[mm,nn,pp])
Eyx1[nn,pp] = Ey[mm-1,nn,pp]
end;end
@inbounds for pp = 1:SizeZ-1
for nn = 1:SizeY
Ez[mm,nn,pp] = Ezx1[nn,pp] + abccoef * (Ez[mm-1,nn,pp] - Ez[mm,nn,pp])
Ezx1[nn,pp] = Ez[mm-1,nn,pp]
end;end
# =============================================================================
#ABC at y boundary plane 1
nn = 1
@inbounds for pp = 1:SizeZ
for mm = 1:SizeX-1
Ex[mm,nn,pp] = Exy0[mm,pp] + abccoef * (Ex[mm,nn+1,pp] - Ex[mm,nn,pp])
Exy0[mm,pp] = Ex[mm,nn+1,pp]
end;end
@inbounds for pp = 1:SizeZ-1
for mm = 1:SizeX
Ez[mm,nn,pp] = Ezy0[mm,pp] + abccoef * (Ez[mm,nn+1,pp] - Ez[mm,nn,pp])
Ezy0[mm,pp] = Ez[mm,nn+1,pp]
end;end
#ABC at x boundary plane 2
nn = SizeY
@inbounds for pp = 1:SizeZ
for mm = 1:SizeX-1
Ex[mm,nn,pp] = Exy1[mm,pp] + abccoef * (Ex[mm,nn-1,pp] - Ex[mm,nn,pp])
Exy1[mm,pp] = Ex[mm,nn-1,pp]
end;end
@inbounds for pp = 1:SizeZ-1
for mm = 1:SizeX
Ez[mm,nn,pp] = Ezy1[mm,pp] + abccoef * (Ez[mm,nn-1,pp] - Ez[mm,nn,pp])
Ezy1[mm,pp] = Ez[mm,nn-1,pp]
end;end
# =============================================================================
#ABC at z boundary plane 1
pp = 1
@inbounds for nn = 1:SizeY
for mm = 1:SizeX-1
Ex[mm,nn,pp] = Exz0[mm,nn] + abccoef * (Ex[mm,nn,pp+1] - Ex[mm,nn,pp])
Exz0[mm,nn] = Ex[mm,nn,pp+1]
end;end
@inbounds for nn = 1:SizeY-1
for mm = 1:SizeX
Ey[mm,nn,pp] = Eyz0[mm,nn] + abccoef * (Ey[mm,nn,pp+1] - Ey[mm,nn,pp])
Eyz0[mm,nn] = Ey[mm,nn,pp+1]
end;end
#ABC at z boundary plane 2
pp = SizeZ
@inbounds for nn = 1:SizeY
for mm = 1:SizeX-1
Ex[mm,nn,pp] = Exz1[mm,nn] + abccoef * (Ex[mm,nn,pp-1] - Ex[mm,nn,pp])
Exz1[mm,nn] = Ex[mm,nn,pp-1]
end;end
@inbounds for nn = 1:SizeY-1
for mm = 1:SizeX
Ey[mm,nn,pp] = Eyz1[mm,nn] + abccoef * (Ey[mm,nn,pp-1] - Ey[mm,nn,pp])
Eyz1[mm,nn] = Ey[mm,nn,pp-1]
end;end
return Ex,Ey,Ez
end
#define the driving source - a Ricker wavelet: https://en.wikipedia.org/wiki/Mexican_hat_wavelet
function ricker(t,location)
src = 0
delay = 5
if t > delay
arg = π * ((Cdtds * (t - (delay+1)) - location) / ppw - 1)
src = (1.0 - 2.0 * arg^2) * exp(-(arg^2))
end
return src
end
#Main control code block. Step through applicable functions
for qtime = 1:MaxTime #time stepping loop
Hx,Hy,Hz = updateH(Ex,Ey,Ez,Hx,Hy,Hz) #magnetic field update equation
Ex,Ey,Ez = updateE(Ex,Ey,Ez,Hx,Hy,Hz) #electric field update equation
Ex[24,23,23] += ricker(qtime,0) #update Ricker wavelet source
Ex,Ey,Ez = abc3dfirst(Ex,Ey,Ez) #apply absorbing boundary conditions
Eout[:,:,qtime] = Ex[24,:,:] #Save constant-x plane time domain data to a 3D matrix
#source[qtime] = ricker(qtime,0) #source output plot for troubleshooting - normally commented out
end
CPUtoc()
println("elapsed real time: ", round(time() - start; digits=3)," seconds")
println("Computation Complete")
end
There’s probably some remaining low-hanging fruit to improve the speed but that’s a good start. I should also mention that this is not written in a recommended julia style, but that’s less important.