No commit message

Project.toml

@@ -23,6 +23,10 @@ Optimisers = "3bd65402-5787-11e9-1adc-39752487f4e2"
 Permutations = "2ae35dd2-176d-5d53-8349-f30d82d94d4f"
 Porcupine = "3b53a3d7-3f34-4bba-8df2-4717a8b1e972"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
+
+[compat]
+SparseArrays = "1.11.0"

luminescent/luminescent/pic/solvemodes.py

@@ -11,6 +11,7 @@
 neigs = data["neigs"]
 
 m, n = eps.shape
+# print(m, n)
 m += 1
 n += 1
 x = np.linspace(0.5*dl, (m-.5)*dl, m)
@@ -29,4 +30,5 @@ def ϵfunc(x_, y_):
     "Ex", "Ey", "Ez", "Hx", "Hy", "Hz"]} for m in solver.modes]
 neffs = [np.real(m.neff) for m in solver.modes]
 for i, mode in enumerate(modes):
+    # print(mode["Ex"].shape)
     np.savez(os.path.join(path, f'mode{i}.npz'), **modes[i])

luminescent/tiny.py

@@ -9,5 +9,8 @@
 path = os.path.join("runs", "tiny")
 lumi.make_pic_sim_prob(path, c, wavelengths=wavelengths, keys=[
                        "2,1"],                       nres=15, approx_2D_mode="TE", gpu="CUDA")  # approx_2D_mode="TE")
+path = os.path.join("runs", "tiny3")
+lumi.make_pic_sim_prob(path, c, wavelengths=wavelengths, keys=[
+                       "2,1"],                       nres=15, gpu="CUDA")  # approx_2D_mode="TE")
 # lumi.solve(path)
 # sol = lumi.lumi_solution()

src/core/monitors.jl

@@ -80,18 +80,41 @@ function MonitorInstance(m::Monitor, g, ϵ, temp, mode_solutions=nothing)
     @unpack lb, ub, center, dimsperm, specs, N, center3, lb3, ub3, approx_2D_mode, tags = m
     @unpack deltas, deltas3, field_lims, F, mode_spacing, dl = g
 
+    start = v2i(center + lb - g.lb, deltas)
+    stop = v2i(center + ub - g.lb, deltas)
+    start, stop = min.(start, stop), max.(start, stop)
+
+    sel = abs.(stop - start) .>= 1e-3
+    stop[!sel] .= start[!sel]
+    start += 0.5sel
+    stop -= 0.5sel
+    start = F(start)
+    stop = F(stop)
+    len = int(stop - start + 1)
+
+    roi = dict([k => begin
+        range.(start, stop, len) - lr[:, 1] + 1
+    end for (k, lr) = pairs(field_lims)])
+    _center = round(v2i(center - g.lb, deltas) + 0.5)
+    #  center, g.lb, _center
+
+
     dx = deltas[1][1]
+
     @unpack λmodenums, λmodes = specs
     if !isnothing(λmodenums)
-        start = round((center3 + lb3) / dl)
-        stop = round((center3 + ub3) / dl)
+        start = round((center3 + lb3) / dl + 0.001)
+        stop = round((center3 + ub3) / dl + 0.001)
         start, stop = min.(start, stop), max.(start, stop)
 
         sel = abs.(stop - start) .>= 0.001
-        stop[!sel] .= start[!sel]
         start += 0.5sel
         stop -= 0.5sel
 
+        ratio = int(dx / dl)
+        stop[1:N] = sel[1:N] .* (ratio * len - 1) + start[1:N]
+        start, stop
+
         start += 0.5
         stop += 0.5
         len = int(stop - start + 1)
@@ -104,23 +127,6 @@ function MonitorInstance(m::Monitor, g, ϵ, temp, mode_solutions=nothing)
         end
     end
 
-    start = v2i(center + lb - g.lb, deltas)
-    stop = v2i(center + ub - g.lb, deltas)
-    start, stop = min.(start, stop), max.(start, stop)
-
-    sel = abs.(stop - start) .>= 1e-3
-    stop[!sel] .= start[!sel]
-    start += 0.5sel
-    stop -= 0.5sel
-    start = F(start)
-    stop = F(stop)
-
-    roi = dict([k => begin
-        range.(start, stop, int(stop - start + 1)) - lr[:, 1] + 1
-    end for (k, lr) = pairs(field_lims)])
-    _center = round(v2i(center - g.lb, deltas) + 0.5)
-    # @show center, g.lb, _center
-
     MonitorInstance(roi, nothing, dimsperm, deltas, _center, fmap(x -> convert.(complex(F), x), λmodes), tags)
 end
 

src/core/sources.jl

@@ -127,8 +127,8 @@ function SourceInstance(s::Source, g, ϵ, temp, mode_solutions=nothing)
     dx = deltas[1][1]
     @unpack λmodenums, λmodes = specs
     if !isnothing(λmodenums)
-        start = round((center3 + lb3) / dl)
-        stop = round((center3 + ub3) / dl)
+        start = round((center3 + lb3) / dl + 0.001)
+        stop = round((center3 + ub3) / dl + 0.001)
         start, stop = min.(start, stop), max.(start, stop)
 
         sel = abs.(stop - start) .>= 0.001

src/del.jl

@@ -10,28 +10,20 @@ function diffpad(a, vl, vr=vl; dims=1, diff=diff, autodiff=true)
     l = !isnothing(vl)
     r = !isnothing(vr)
 
-
-    if autodiff
-        # b = Buffer(a, sz)
-        b = diff(a; dims)
-    else
-        sz = Tuple(size(a) + (l + r - 1) * sel)
-        b = similar(a, sz)
-    end
-
     # @time b[range.(l * sel + 1, sz - r * sel)...] = diff(a; dims)
     # @time pad!(b, vl, l * sel, 0)
     # @time pad!(b, vr, 0, r * sel)
     # println()
 
     # if vl ∈ (0, nothing) && vr ∈ (0, nothing)
     if autodiff
-        b = pad(b, vl, l * sel, 0)
-        b = pad(b, vr, 0, r * sel)
+        b = diff(a; dims)
+        b = pad(b, vl, vr, l * sel, r * sel)
     else
+        sz = Tuple(size(a) + (l + r - 1) * sel)
+        b = similar(a, sz)
         b[range.(l * sel + 1, sz - r * sel)...] = diff(a; dims)
-        pad!(b, vl, l * sel, 0)
-        pad!(b, vr, 0, r * sel)
+        pad!(b, vl, vr, l * sel, r * sel)
     end
     @assert typeof(b) == typeof(a)
     b
@@ -72,7 +64,7 @@ function LinearAlgebra.cross(m::Del, v)
 end
 
 function delcross(diff, Δs, ps, as)
-    autodiff = !haskey(ENV, "autodiff") || ENV["autodiff"] == "1"
+    autodiff = AUTODIFF()
 
 
     _diffpad(a, p, dims) = diffpad(a, p...; dims, diff, autodiff)

src/main.jl

@@ -1,4 +1,4 @@
-using UnPack, LinearAlgebra, Statistics, Random, Jello, DataStructures, FileIO, Porcupine, Dates, NPZ, DataStructures, JSON, Flux, Zygote, CairoMakie, ArrayPadding, Permutations, Functors, ImageBase, Optimisers, BFloat16s, ChainRulesCore
+using UnPack, LinearAlgebra, Statistics, Random, Jello, DataStructures, FileIO, Porcupine, Dates, NPZ, DataStructures, JSON, Flux, Zygote, CairoMakie, ArrayPadding, Permutations, Functors, ImageBase, Optimisers, BFloat16s, ChainRulesCore, SparseArrays
 using Porcupine: keys, values, pairs, fmap, ⊙, trim, round, floor, ceil, invperm, permutedims, dict, cpu, gpu
 using Flux: mae, Adam
 using Zygote: withgradient, Buffer, ignore_derivatives, @ignore_derivatives

src/pic/run.jl

@@ -48,34 +48,52 @@ function picrun(path; array=Array, kw...)
     sz = round.(L / dl)
 
     GC.gc(true)
-    Nf = study == "inverse_design" ? F : Float16
+    if AUTODIFF()
+        Nf = F
+    else
+        # Nf = [Bool, UInt8][findfirst(>=(length(enum), 2 .^ [1, 8]))]
+        Nf = Float16
+    end
+    enum = [matprops(v.material).ϵ |> F for v = values(layer_stack)] |> unique |> sort
+    ϵmin = minimum(enum) |> Nf
+    # if AUTODIFF()
     ϵ3 = zeros(Nf, Tuple(sz))
+    # else
+    #     ϵ3 = SparseArrays.zeros(Nf, Tuple(sz))
+    # end
+    # @show Nf
     layer_stack = sort(collect(pairs(layer_stack)), by=kv -> -kv[2].mesh_order) |> OrderedDict
-    ϵmin = 100
     for (k, v) = pairs(layer_stack)
+        @unpack material, thickness = v
+        ϵ = matprops(material).ϵ |> Nf
+        # if AUTODIFF() || ϵ != ϵmin
         a = stack(map(sort(collect(readdir(joinpath(temp, string(k)), join=true)))) do file
-            a = Nf.(Gray.(FileIO.load(file)))
+            a = Nf.(Bool.(FileIO.load(file)))
             a = downsample(a, 2)
-            reverse(a', dims=2)
+            a = reverse(a', dims=2)
+            # if !AUTODIFF()
+            #     a = sparse(a)
+            # end
+            a
         end)
         a = downsample(a, (1, 1, 2))
+        # @show typeof(a)
         # a = a[Base.OneTo.(min.(size(a), sz))...]
-        @unpack material, thickness = v
 
         start = 1 + round.([(v.origin / dl)..., (v.zmin - zmin) / dl])
         I = range.(start, start + size(a) - 1)
         overhang = max.(last.(I) .- sz, 0)
         a = a[Base.oneto.(size(a) - overhang)...]
         I = range.(start, start + size(a) - 1)
 
-        ϵ = matprops(material).ϵ |> Nf
-        ϵmin = min(ϵ, ϵmin) |> Nf
         ϵ3[I...] .*= 1 - a
         ϵ3[I...] .+= a .* ϵ
+        # end
     end
     ϵ3 = max.(ϵmin, ϵ3)
     # @show eltype(ϵ3)
     @assert eltype(ϵ3) == Nf
+    # @show typeof(ϵ3)
     GC.gc(true)
 
     ϵ2 = ϵ3[:, :, 1+round((zcenter - zmin) / dl)]

src/pictest.jl

@@ -8,11 +8,13 @@ using CUDA
 # picrun(joinpath("runs", "straight");)# array=cu)
 # picrun(joinpath("runs", "bend_R5"), array=cu)
 # picrun(joinpath("runs", "mode_converter"))
-picrun(joinpath("runs", "demux"))
 
 # picrun(joinpath("runs", "splitter"); array=cu)
 # picrun(joinpath("runs", "splitter"))
 
 # picrun(joinpath("runs", "tiny"); array=cu)
+# picrun(joinpath("runs", "tiny3"); array=cu)
 # picrun(joinpath("runs", "back");)# array=cu)
 # models[1]()
+
+picrun(joinpath("runs", "demux"))

src/sim/setup.jl

@@ -312,7 +312,7 @@ function setup(dl, boundaries, sources, monitors, deltas, mode_deltas;
             Tss *= ceil(Tssmin / Tss)
         end
     end
-    # Tss *= 4
+    Tss = max(Tss, Tssmin)
     @show Ttrans, Tss
     Ttrans, Tss = convert.(F, (Ttrans, Tss))
     prob = (;