No commit message

examples/inverse_design_signal_splitter/inverse_design_signal_splitter.jl

@@ -5,21 +5,21 @@ using Zygote: withgradient, Buffer
 using BSON: @save, @load
 using Optim: Options, minimizer
 using AbbreviatedStackTraces
-using Jello, Luminescent, LuminescentVisualization
+# using Jello, Luminescent, LuminescentVisualization
 Random.seed!(1)
 
-# dir = pwd()
-# include("$dir/src/main.jl")
-# include("$dir/../LuminescentVisualization.jl/src/main.jl")
-# include("$dir/scripts/startup.jl")
+dir = pwd()
+include("$dir/src/main.jl")
+include("$dir/../LuminescentVisualization.jl/src/main.jl")
+include("$dir/scripts/startup.jl")
 
 # loads design layout
 @load "$(@__DIR__)/layout.bson" base signals ports designs
 @load "$(@__DIR__)/modes.bson" modes lb ub λ dx hsub wwg hwg hclad ϵsub ϵclad ϵwg
 
 # training params"
 F = Float32
-dogpu = true
+dogpu = false
 name = "inverse_design_signal_splitter"
 nbasis = 5 # complexity of design region
 contrast = 10.0
@@ -29,6 +29,8 @@ rmin = nothing
 T1 = Δ1 = 1 + norm(ports[1].c - ports[2].c) / λ * sqrt(ϵwg) # simulation duration in [periods] for signal to reach output ports
 Δ2 = 1 # duration to record power at output ports
 T2 = T = T1 + Δ2
+T1 = 0.1
+T2 = 0.2
 ϵmin = ϵclad
 hsub, wwg, hwg, hclad, dx, ub, lb = [hsub, wwg, hwg, hclad, dx, ub, lb] / λ
 
@@ -88,8 +90,8 @@ function make_geometry(model, base, μ, σ, σm)
     place!(base_, model(), round.(Int, designs[1].o / λ / dx) .+ 1)
 
     ϵ = sandwich(copy(base_), round.(Int, [hsub, hwg, hclad] / dx), [ϵsub, ϵwg, ϵclad])
-    ϵ, μ, σ, σm = apply(geometry_padding; ϵ, μ, σ, σm)
-    p = apply(geometry_staggering; ϵ, μ, σ, σm)
+    p = apply(geometry_padding; ϵ, μ, σ, σm)
+    p = apply(geometry_staggering, p)
 end
 
 function metrics(model, ; T1=T1, T2=T2, autodiff=true)
@@ -116,7 +118,7 @@ function loss(v)
 end
 
 p0 = make_geometry(model0, base, μ, σ, σm)
-volume(cpu(p0[1][2]))
+# volume(cpu(p0[1][2]))
 # error()
 
 # gradient free optimization "Optim functions

examples/inverse_design_signal_splitter/scratch.jl

@@ -24,3 +24,4 @@ end
 
 modes = main()
 plot_mode_fields(modes[1])
+

examples/periodic_scattering/3d_periodic_scattering.jl

@@ -3,14 +3,15 @@ simulation of plane wave scattering on periodic array of dielectric spheres
 """
 
 using UnPack, LinearAlgebra, GLMakie
-using Luminescent, LuminescentVisualization
+# using Luminescent, LuminescentVisualization
 
-# dir = pwd()
-# include("$(dir)/src/main.jl")
-# include("$(dir)/scripts/startup.jl")
-# include("$dir/../LuminescentVisualization.jl/src/main.jl")
+dir = pwd()
+include("$(dir)/src/main.jl")
+include("$(dir)/scripts/startup.jl")
+include("$dir/../LuminescentVisualization.jl/src/main.jl")
 
-dogpu = true
+dogpu = false
+F = Float32
 name = "periodic_scattering"
 T = 10 # simulation duration in [periods]
 nx = 20
@@ -42,8 +43,8 @@ monitors = [
 configs = maxwell_setup(boundaries, sources, monitors, dx, sz; ϵmin, F)
 @unpack dt, geometry_padding, geometry_staggering, field_padding, source_instances, monitor_instances, u0, = configs
 
-ϵ, μ, σ, σm = apply(geometry_padding; ϵ, μ, σ, σm)
-p = apply(geometry_staggering; ϵ, μ, σ, σm)
+p = apply(geometry_padding; ϵ, μ, σ, σm)
+p = apply(geometry_staggering; p...)
 
 # move to gpu
 if dogpu
@@ -54,7 +55,8 @@ end
 
 # run simulation
 @showtime u = accumulate(0:dt:T, init=u0) do u, t
-    maxwell_update!(deepcopy(u), p, t, dx, dt, field_padding, source_instances)
+    maxwell_update(u, p, t, dx, dt, field_padding, source_instances)
+    # maxwell_update!(deepcopy(u), p, t, dx, dt, field_padding, source_instances)
 end
 v = [power_flux.((m,), u) for m = monitor_instances]
 
@@ -64,10 +66,8 @@ if dogpu
 end
 
 # make movie, 
-Ez = map(u) do u
-    u[1][3]
-end
-ϵEz = p[1][3]
+Ez = field.(u, :Ez)
+ϵEz = p[:ϵ][3]
 dir = @__DIR__
 recordsim("$dir/$(name).mp4", Ez, v;
     dt,

scripts/setup.jl

@@ -1,6 +1,6 @@
 using Pkg
 # ]add UnPack, BSON,DataStructures, StatsBase, Documenter,ImageTransformations,Interpolations, Zygote, Optim,ArrayPadding, Porcupine, Jello, GLMakie, Functors, Lazy,AbbreviatedStackTraces,ReverseStackTraces,CUDA,Flux
 # Pkg.add(url="https://github.com/paulxshen/ArrayPadding.jl")
-# Pkg.add(url="https://github.com/paulxshen/Porcupine.jl")
+Pkg.add(url="https://github.com/paulxshen/Porcupine.jl")
 Pkg.add(url="https://github.com/paulxshen/Luminescent.jl")
 Pkg.add(url="https://github.com/paulxshen/LuminescentVisualization.jl")

src/boundaries.jl

@@ -92,7 +92,8 @@ function apply(p::AbstractVector{<:InPad}, a::AbstractArray)
     end
 
     a_ = Buffer(a)
-    a_ .= a
+    # a_ .= a
+    a_[axes(a)...] = a
     for p = p
         @unpack l, r, b, m = p
         if isnothing(m)

src/maxwell_setup.jl

@@ -181,8 +181,14 @@ function maxwell_setup(boundaries, sources, monitors, dx, sz, polarization=nothi
     if d == 1
         pf = u -> [u[1] .* u[2]]
     elseif d == 3
-        u0 = [u[1:3], u[4:6]]
+        u0 = MyDict([
+            :E => MyDict([k => zeros(F, Tuple(sizes[k])) for k = (:Ex, :Ey, :Ez)]),
+            :H => MyDict([k => zeros(F, Tuple(sizes[k])) for k = (:Hx, :Hy, :Hz)]),])
+        # u0 = [u[1:3], u[4:6]]
     else
+        u0 = MyDict([
+            :E => MyDict([k => zeros(F, Tuple(sizes[k])) for k = (:Ex, :Ey)]),
+            :H => MyDict([k => zeros(F, Tuple(sizes[k])) for k = (:Hz,)]),])
         if polarization == :TMz
             u -> [-u[1] .* u[3], u[2] .* u[1]]
         else

src/maxwell_update.jl

@@ -1,15 +1,3 @@
-function mark(p; kw...)
-    [mark(p[k], kw[k]) for k = keys(kw)]
-end
-function mark(v, a)
-    l = sum(v) do p
-        p.l
-    end
-    r = sum(v) do p
-        p.r
-    end
-    PaddedArray(a, l, r)
-end
 
 """
     function maxwell_update!(u, p, t, field_padding, source_instances)
@@ -19,9 +7,11 @@ Updates fields for 3d. Please use `maxwell_update` instead of `maxwell_update!`
 """
 function maxwell_update!(u1, u, p, t, dx, dt, field_padding, source_instances)
     ∇ = StaggeredDel([dx, dx, dx])
-    ϵ, μ, σ, σm = p
-    E, H = u
-    E1, H1 = u1
+    @unpack ϵ, μ, σ, σm = p
+    # ϵ, μ, σ, σm = p |> values
+    @unpack E, H = u
+    E1 = u1[:E]
+    H1 = u1[:H]
     # J = apply(source_instances, t; Jx=zeros(F, size(E[1])), Jy=zeros(F, size(E[2])), Jz=zeros(F, size(E[3])))
     J = apply(source_instances, t; Jx=0, Jy=0, Jz=0)
 
@@ -35,7 +25,7 @@ function maxwell_update!(u1, u, p, t, dx, dt, field_padding, source_instances)
     end
     Ex, Ey, Ez = E1
     apply!(field_padding; Ex, Ey, Ez)
-    H = collect.(H)
+    H = unmark(; H...)
 
     # then update H
     E = mark(field_padding; Ex, Ey, Ez)
@@ -46,93 +36,62 @@ function maxwell_update!(u1, u, p, t, dx, dt, field_padding, source_instances)
     end
     Hx, Hy, Hz = H1
     apply!(field_padding; Hx, Hy, Hz)
-    E = collect.(E)
+    E = unmark(; E...)
 
     [E, H]
     # u1
 end
 maxwell_update!(u, p, t, dx, dt, field_padding, source_instances) = maxwell_update!(u, u, p, t, dx, dt, field_padding, source_instances)
 
-# function maxwell_update!(u, p, t, dx, dt, field_padding, source_instances)
-#     ∇ = StaggeredDel([dx, dx, dx])
-#     ϵ, μ, σ, σm = p
-#     E, H = u
-#     J = apply(source_instances, t; Jx=0, Jy=0, Jz=0)
-
-#     # first update E
-#     Hx, Hy, Hz = H
-#     H = mark(field_padding; Hx, Hy, Hz)
-#     dEdt = (∇ × H - σ * E - J) / ϵ
-
-#     E_ = bufferfrom.(E)
-#     for i = 1:3
-#         E_[i][:, :, :] = E[i] + dEdt[i] * dt
-#     end
-#     Ex, Ey, Ez = [copy(E_[1]), copy(E_[2]), copy(E_[3])]
-#     # E .= E + dEdt * dt
-
-#     apply!(field_padding; Ex, Ey, Ez)
-#     H = collect.(H)
-
-#     # then update H
-#     E = mark(field_padding; Ex, Ey, Ez)
-#     dHdt = -(∇ × E .+ σm * H) / μ
-
-#     H_ = bufferfrom.(H)
-#     for i = 1:3
-#         H_[i][:, :, :] = H[i] + dHdt[i] * dt
-#     end
-#     Hx, Hy, Hz = [copy(H_[1]), copy(H_[2]), copy(H_[3])]
-#     # H .= H + dHdt * dt
-
-#     apply!(field_padding; Hx, Hy, Hz)
-#     H = [Hx, Hy, Hz]
-#     E = collect.(E)
-
-#     [E, H]
-# end
 """
     function maxwell_update(u, p, t, field_padding, source_instances)
 
 Updates fields for 3d in a manner amenable to AD. See also Mutating `maxwell_update!`
 """
 function maxwell_update(u, p, t, dx, dt, field_padding, source_instances; ignore_boundary_autodiff=false)
     ∇ = StaggeredDel([dx, dx, dx])
-    ϵ, μ, σ, σm = p
-    E, H = u
+    @unpack ϵ, μ, σ, σm = p
+    @unpack E, H = u
     J = apply(source_instances, t; Jx=0, Jy=0, Jz=0)
 
     # first update E
-    Hx, Hy, Hz = H
-    H = mark(field_padding; Hx, Hy, Hz)
+    # H = mark(field_padding; H...)
+    H = mark(field_padding, H)
+    # @info typeof(E)
     dEdt = (∇ × H - σ * E - J) / ϵ
-    Ex, Ey, Ez = E + dEdt * dt
-
-    if ignore_boundary_autodiff
-        ignore_derivatives() do
-            apply!(field_padding; Ex, Ey, Ez)
-        end
-    else
-        Ex, Ey, Ez = apply(field_padding; Ex, Ey, Ez)
-    end
-    H = collect.(H)
+    # @info typeof(dEdt)
+    E = E + dEdt * dt
+
+    # Ex, Ey, Ez = E
+    # E = apply(field_padding; Ex, Ey, Ez)
+    # Hx, Hy, Hz = H
+    # H = unmark(; Hx, Hy, Hz)
+
+    # E = apply(field_padding; E...)
+    # H = unmark(; H...)
+    E = apply(field_padding, E)
+    H = unmark(H)
 
     # then update H
-    E = mark(field_padding; Ex, Ey, Ez)
+    # E = mark(field_padding; Ex, Ey, Ez)
+    E = mark(field_padding, E)
     dHdt = -(∇ × E .+ σm * H) / μ
-    Hx, Hy, Hz = H + dHdt * dt
-
-    if ignore_boundary_autodiff
-        ignore_derivatives() do
-            apply!(field_padding; Hx, Hy, Hz)
-        end
-    else
-        Hx, Hy, Hz = apply(field_padding; Hx, Hy, Hz)
-    end
-    H = [Hx, Hy, Hz]
-    E = collect.(E)
-
-    [E, H]
+    H = H + dHdt * dt
+
+    # Ex, Ey, Ez = E
+    # Ex, Ey, Ez = unmark(; Ex, Ey, Ez)
+    # E = (; Ex, Ey, Ez)
+    # Hx, Hy, Hz = H
+    # H = apply(field_padding; Hx, Hy, Hz)
+
+    # H = apply(field_padding; H...)
+    # E = unmark(; E...)
+    H = apply(field_padding, H)
+    E = unmark(E)
+
+    # [E, H]
+    MyDict([:E => E, :H => H])
+    # (; E, H)
 end
 maxwell_update! = maxwell_update!
 # Flux.trainable(m::PaddedArray) = (; a=m.a)

src/monitors.jl

@@ -113,8 +113,9 @@ end
 
 function field(u, k)
     if k in keys(fij)
-        i, j = fij[k]
-        u[i][j]
+        u[k[1]][k]
+        # i, j = fij[k]
+        # u[i][j]
     elseif k == "|E|2"
         sum(u[1]) do a
             a .^ 2
@@ -181,12 +182,9 @@ end
  power_flux density (avg Poynting flux) passing thru monitor 
 """
 function power_flux_density(m::OrthogonalMonitorInstance, u)
-    E = [u[1][i][m.fi[k]...] for (i, k) = enumerate([:Ex, :Ey, :Ez])]
-    H = [u[2][i][m.fi[k]...] for (i, k) = enumerate([:Hx, :Hy, :Hz])]
-
-    # E = getindex.(u[1], Ref.(m.i)...)
-    # H = getindex.(u[2], Ref.(m.i)...)
-    # E, H = u
+    d = ndims(m)
+    E = [u[:E][k][m.fi[k]...] for k = keys(u[:E])]
+    H = [u[:H][k][m.fi[k]...] for k = keys(u[:H])]
 
     r = mean(sum((E × H) .* m.n))
     r

src/sources.jl

@@ -156,11 +156,11 @@ function apply(v::AbstractVector{<:SourceInstance}, t::Real; kw...)
                 kw[k] .+ sum(a)
             end
         end
-        for k = keys(kw)
+        for k = _keys(kw)
         # end for (k, a) = pairs(kw)
     ]
 end
 
 # function apply(d,t; kw...)
-#     [apply(d[k], kw[k]) for k = keys(kw)]
+#     [apply(d[k], kw[k]) for k = _keys(kw)]
 # end