Currently Prerelease. Expect breaking changes. Report bugs on Github - we usually respond within a day
Generative design meets Maxwell's Equations. Differentiable FDTD package for inverse design & topology optimization in semiconductor photonics, acoustics and RF. GPU and automatic differentiation (AD) compatible. Uses AD by Zygote.jl for adjoint optimization. Integrates with Jello.jl to generate length scale controlled paramaterized geometry . Staggered Yee grid update with fully featured boundary conditions & sources. Customizable physics to potentially incorporate dynamics like heat transfer, charge transport.
Install via
Pkg.add(url="https://github.com/paulxshen/FDTDEngine.jl")
Pkg.add(url="https://github.com/paulxshen/FDTDToolkit.jl")
FDTDToolkit.jl contains visualization utilities
We do a quick 3d simulation of plane wave scattering on periodic array of dielectric spheres (see gallery movie)
"""
simulation of plane wave scattering on periodic array of dielectric spheres
"""
using UnPack, LinearAlgebra, GLMakie
-# using FDTDEngine,FDTDToolkit
-dir = pwd()
-include("$(dir)/src/main.jl")
-include("$dir/../FDTDToolkit.jl/src/main.jl")
+using FDTDEngine,FDTDToolkit
+dogpu = true
+# dogpu = false
+name = "periodic_scattering"
+T = 10 # simulation duration in [periods]
+nx = 20
+dx = 1.0 / nx # pixel resolution in [wavelengths]
-name = "3d_scattering"
-T = 8.0f0 # simulation duration in [periods]
-nx = 16
-dx = 1.0f0 / nx # pixel resolution in [wavelengths]
-
-"geometry"
-l = 2 # domain physical size length
+# geometry
+l = 2 # domain physical size length in [wavelengths]
sz = nx .* (l, l, l) # domain voxel dimensions
ϵ1 = ϵmin = 1 #
-ϵ2 = 2.25f0 #
+ϵ2 = 2.25 #
b = F.([norm(v .- sz ./ 2) < 0.5 / dx for v = Base.product(Base.oneto.(sz)...)]) # sphere
ϵ = ϵ2 * b + ϵ1 * (1 .- b)
-"setup"
+# setup
boundaries = [Periodic(2), Periodic(3)]# unspecified boundaries default to PML
sources = [
- PlaneWave(t -> cos(F(2π) * t), -1; Jz=1) # Jz excited plane wave from -x plane (eg -1)
+ PlaneWave(t -> cos(2π * t), -1; Jz=1) # Jz excited plane wave from -x plane (eg -1)
]
-n = [1, 0, 0] # normal
-δ = 0.2f0 # margin
-# A = (l - δ)^2
+normal = [1, 0, 0] #
+δ = 0.2 # margin
lm = 1 # monitor side length
monitors = [
- Monitor([δ, l / 2, l / 2], [0, lm, lm], n,), # (center, dimensions, normal)
- Monitor([l - δ, l / 2, l / 2], [0, lm, lm], n,),
+ Monitor([δ, l / 2, l / 2], [0, lm, lm]; normal), # (center, dimensions; normal)
+ Monitor([l - δ, l / 2, l / 2], [0, lm, lm]; normal),
]
configs = setup(boundaries, sources, monitors, dx, sz; ϵmin, T)
@unpack μ, σ, σm, dt, geometry_padding, geometry_splits, field_padding, source_instances, monitor_instances, u0, = configs
@@ -43,35 +40,38 @@
ϵ, μ, σ, σm = apply(geometry_padding; ϵ, μ, σ, σm)
p = apply(geometry_splits; ϵ, μ, σ, σm)
+# move to gpu
+if dogpu
+ using CUDA, Flux
+ @assert CUDA.functional()
+ u0, p, field_padding, source_instances = gpu.((u0, p, field_padding, source_instances))
+end
# run simulation
-t = 0:dt:T
-u = [[similar.(a) for a = u0] for t = t]
-u[1] = u0
-@showtime reduce(
- (u, (u1, t)) -> step3!(u1, u, p, t, dx, dt, field_padding, source_instances),
- zip(u[2:end], t[1:end-1]),
- init=u0)
-y = hcat([power.((m,), u) for m = monitor_instances]...)
+@showtime u = accumulate(0:dt:T, init=u0) do u, t
+ step3!(deepcopy(u), p, t, dx, dt, field_padding, source_instances)
+end
+y = [power.((m,), u) for m = monitor_instances]
-# make movie
+# move back to cpu for plotting
+if dogpu
+ u, p, field_padding, source_instances = cpu.((u, p, field_padding, source_instances))
+end
+
+# make movie,
Ez = map(u) do u
u[1][3]
end
-ϵz = p[1][3]
+ϵEz = p[1][3]
dir = @__DIR__
-
-recordsim("$dir/$(name)_nres_$nx.mp4", Ez, y;
+recordsim("$dir/$(name).mp4", Ez, y;
dt,
+ field=:Ez,
monitor_instances,
source_instances,
- geometry=ϵz,
+ geometry=ϵEz,
elevation=30°,
playback=1,
axis1=(; title="$name\nEz"),
axis2=(; title="monitor powers"),
-)
-
<!– Supports 1d (Ez, Hy), 2d TMz (Ez, Hx, Hy), 2d TEz (Hz, Ex, Ey) and 3d. –> Length and time are in units of wavelength and period. This normalization allows usage of relative permitivity and permeability in equations . Fields including electric, magnetic and current density are simply bundled as a vector of vectors of arrays . Boundary conditions pad the field arrays . PML paddings are multilayered, while All other boundaries add single layers. Paddings are stateful and permanent, increasing the size of field and geometry arrays. Finite differencing happens every update step3 and are coordinated to implictly implement a staggered Yee's grid .
If a source has fewer nonzero dimensions than the simulation domain, its signal will get normalized along its singleton dimensions. For example, all planar sources in 3d or line sources in 2d will get scaled up by a factor of 1/dx. This way, discretisation would not affect radiated power.
function PlaneWave(f, dims; fields...)
Constructs plane wave source
Args
- f: time function
- dims: eg -1 for wave coming from -x face
- fields: which fields to excite & their scaling constants (typically a current source, eg Jz=1)
function Source(f, c, lb, ub, label=""; fields...)
-function Source(f, c, L, label=""; fields...)
Constructs custom source. Can be used to specify uniform or modal sources
Args
- f: time function
- c: origin or center of source
- lb: lower bounds wrt to c
- ub: upper bounds wrt to c
- L: source dimensions in [wavelengths]
- fields: which fields to excite & their scaling constants (typically a current source, eg Jz=1)
<!– GaussianBeam –>
Unspecified boundaries default to PML
Periodic(dims)
periodic boundary
function PML(dims, d=0.25f0, σ=20.0f0)
Constructs perfectly matched layers (PML aka ABC, RBC) boundary of depth d wavelengths Doesn't need to be explictly declared as all unspecified boundaries default to PML
PEC(dims)
perfect electrical conductor dims: eg -1 for -x side
PMC(dims)
perfect magnetic conductor
function Monitor(c, L; normal=nothing, label="")
-function Monitor(c, lb, ub; normal=nothing, label="")
Constructs monitor which can span a point, line, surface, or volume monitoring fields or power.
Args
c: origin or center of monitor
L: physical dimensions of monitor
lb: lower bounds wrt to c
ub: upper bounds wrt to c
normal: flux monitor direction (eg normal to flux surface)
function power(m::MonitorInstance, u)
total power (Poynting flux) passing thru monitor surface
function power_density(m::MonitorInstance, u)
power density (avg Poynting flux) passing thru monitor surface
function step3!(u, p, t, field_padding, source_instances)
-function step3!(u1, u, p, t, field_padding, source_instances)
Updates fields for 3d. Please use step3 instead of step3! when doing AD. Mutating step3! Writes new fields either onto old fields or into buffer arrays u1
function step3(u, p, t, field_padding, source_instances)
Updates fields for 3d in a manner amenable to AD. See also Mutating step3!
<!– step1 –> <!– stepTMz –> <!– stepTEz –>
Simply use Flux.gpu to move simulation variables to GPU. This turns Arrays into CUDA arrays which get processed on GPU for both forward and backpropagation passes
Compatible with Zygote.jl and Flux.jl. Please use step3 instead of step3! when doing AD. See inverse design examples
Please contact us for latest scripts. We wrote (Jello.jl)[https://github.com/paulxshen/Jello.jl], an innovative Fourier domain neural model to generate length scale controlled efficiently paramaterized geometry .
Our focus is on inverse design and topology optimization using adjoints from automatic differentiation. We love a streamlined API backed by a clear, concise and extensible codebase. Attention is paid to speed and malloc but that's not our main concern. Meep is the most popular open source FDTD package owing to its maturity and comprehensive features. It supports AD in certain cases using custom adjoints which we avoid in our package in favor of more flexibility . There are numerous commercial software eg Ansys Lumerical, Comsol and various EDA or SI software. TMK none of them supports AD for inverse design
see (examples/)[https://github.com/paulxshen/FDTDEngine.jl/blob/main/examples/]
Discussion & updates at Julia Discourse
Paul Shen <pxshen@alumni.stanford.edu> Consulting and technical support available 2024 (c) Paul Shen