add bent waveguide example, implement pml

docs/_toc.yml

@@ -43,6 +43,7 @@ parts:
   - caption: Julia Examples
     chapters:
     - file: julia/waveguide.jl
+    - file: julia/waveguide_bent.jl
     - file: julia/waveguide_anisotropic.jl
     - file: julia/waveguide_overlap_integral.jl
     - file: julia/lithium_niobate_phase_shifter.jl

docs/julia/waveguide bent.jl

@@ -15,7 +15,7 @@
 # ---
 
 # %% [markdown]
-# # Mode solving
+# # Mode solving for bent waveguides
 
 # %% [markdown]
 # ```{caution}
@@ -32,18 +32,21 @@ clip_by_rect = pyimport("shapely.ops").clip_by_rect
 OrderedDict = pyimport("collections").OrderedDict
 mesh_from_OrderedDict = pyimport("femwell.mesh").mesh_from_OrderedDict
 
-radius = 100000
+radius = 1
 wg_width = 0.5
 wg_thickness = 0.22
-sim_width = 3
-sim_height = 4
+sim_left = 0.5
+sim_right = 3
+sim_top = 1
+sim_bottom = 3
+pml_thickness = 3
 core = shapely.geometry.box(radius - wg_width / 2, 0, radius + wg_width / 2, wg_thickness)
 
 env = shapely.geometry.box(
-    radius - sim_width / 2,
-    -sim_height / 2,
-    radius + sim_width / 2,
-    sim_height / 2,
+    radius - wg_width / 2 - sim_left,
+    -sim_bottom - pml_thickness,
+    radius + wg_width / 2 + sim_right + pml_thickness,
+    wg_thickness + sim_top,
 )
 
 polygons = OrderedDict(
@@ -53,14 +56,14 @@ polygons = OrderedDict(
 )
 
 resolutions = Dict(
-    "core" => Dict("resolution" => 0.03, "distance" => 0.5),
-    "slab" => Dict("resolution" => 0.03, "distance" => 0.5),
+    "core" => Dict("resolution" => 0.015, "distance" => 0.5),
+    "slab" => Dict("resolution" => 0.015, "distance" => 0.5),
 )
 
 mesh = mesh_from_OrderedDict(
     polygons,
     resolutions,
-    default_resolution_max = 1,
+    default_resolution_max = 0.2,
     filename = "mesh.msh",
 )
 
@@ -86,18 +89,33 @@ model = GmshDiscreteModel("mesh.msh")
 
 labels = get_face_labeling(model)
 
-epsilons = ["core" => 3.5^2, "box" => 1.444^2, "clad" => 1.44^2]
+epsilons = ["core" => 3.48^2, "box" => 1.46^2, "clad" => 1.0^2]
 ε(tag) = Dict(get_tag_from_name(labels, u) => v for (u, v) in epsilons)[tag]
 
 
 #dΩ = Measure(Ω, 1)
 τ = CellField(get_face_tag(labels, num_cell_dims(model)), Ω)
-
-modes = calculate_modes(model, ε ∘ τ, λ = 1.55, num = 2, order = 1, radius = radius)
-println(n_eff(modes[2]))
+pml_x = x -> 0
+pml_y = x -> 0
+pml_x = x -> 0.1 * max(0, x[1] - (radius + wg_width / 2 + sim_right))^2
+pml_y = x -> 0.1 * max(0, -x[2] - sim_bottom)^2
+
+modes = calculate_modes(model, ε ∘ τ, λ = 1.55, num = 2, order = 1)
+modes = calculate_modes(
+    model,
+    ε ∘ τ,
+    λ = 1.55,
+    num = 2,
+    order = 1,
+    radius = radius,
+    pml = [pml_x, pml_y],
+    k0_guess = modes[1].k,
+)
+println(n_eff(modes[1]))
+println(log10(abs(imag(n_eff(modes[1])))))
 # write_mode_to_vtk("mode", modes[2])
 
-plot_mode(modes[2], absolute = true)
+plot_mode(modes[1], absolute = true)
 #plot_mode(modes[2])
 modes
 
@@ -106,8 +124,8 @@ modes
 
 # %% tags=["remove-stderr"]
 
-modes_p = calculate_modes(model, ε ∘ τ, λ = 1.55, num = 2, order = 1)
-println(n_eff(modes_p[2]))
-plot_mode(modes_p[2], absolute = true)
+#modes_p = calculate_modes(model, ε ∘ τ, λ = 1.55, num = 2, order = 1)
+#println(n_eff(modes_p[1]))
+#plot_mode(modes_p[1], absolute = true)
 
 # %%

src/Maxwell/Waveguide/Waveguide.jl

@@ -74,6 +74,7 @@ function calculate_modes(
     radius::Real = Inf,
     k0_guess::Union{Number,Nothing} = nothing,
     metallic_boundaries = String[],
+    pml::Union{Vector{Function},Nothing} = nothing,
 )
     if count(isnothing, [k0, λ]) != 1
         throw(ArgumentError("Exactly one of k0,λ must be defined"))
@@ -102,7 +103,7 @@ function calculate_modes(
     epsilons = ε(get_cell_points(Measure(Ω, 1)))
     k0_guess =
         isnothing(k0_guess) ? k0^2 * maximum(maximum.(maximum.(real.(epsilons)))) * 1.1 :
-        k0_guess
+        k0_guess^2
 
     lhs, rhs = if radius == Inf
         μ_r = 1
@@ -119,19 +120,35 @@ function calculate_modes(
 
         lhs_straight, rhs_straight
     else
+        V_pml = TestFESpace(model, ReferenceFE(lagrangian, Float64, 2))
+        pml_f = interpolate(pml[1], V_pml)
+        pml_g = interpolate(pml[2], V_pml)
+
+        s_r = 1 - 1im * gradient(pml_f) ⋅ VectorValue(1, 0)
+        s_y = 1 - 1im * gradient(pml_g) ⋅ VectorValue(0, 1)
+        γ_θ = 1 - 1im * pml_f / (x -> x[1])
+
+        d_θ = s_r * s_y / γ_θ
+        D_t =
+            γ_θ *
+            (s_y / s_r * TensorValue(1, 0, 0, 0) + s_r / s_y * TensorValue(0, 0, 0, 1))
+
         radius_function(x) = x[1]
         r = interpolate_everywhere(radius_function, V2)
         lhs_radial((u1, u2), (v1, v2)) =
             ∫(
-                -(r / radius * curl(u1) ⋅ curl(v1)) +
-                (radius / r * gradient(u2) ⋅ v1) +
-                (k0^2 * ε * r / radius * u1 ⋅ v1) -
-                (radius / r * gradient(u2) ⋅ gradient(v2)) +
-                (k0^2 * ε * radius / r * u2 * v2),
+                -(1 / d_θ * r / radius * curl(u1) ⋅ curl(v1)) +
+                (radius / r / (γ_θ * γ_θ) * (D_t ⋅ gradient(u2)) ⋅ v1) +
+                (k0^2 * ε * r / radius * (D_t ⋅ u1) ⋅ v1) -
+                (radius / r / (γ_θ * γ_θ) * (D_t ⋅ gradient(u2)) ⋅ gradient(v2)) +
+                (k0^2 * ε * d_θ * radius / r * u2 * v2),
             )dΩ
 
         rhs_radial((u1, u2), (v1, v2)) =
-            ∫((radius / r * u1 ⋅ v1) - (radius / r * u1 ⋅ gradient(v2)))dΩ
+            ∫(
+                (radius / r * (D_t ⋅ u1) ⋅ v1) -
+                (radius / r / (γ_θ * γ_θ) * (D_t ⋅ u1) ⋅ gradient(v2)),
+            )dΩ
 
         lhs_radial, rhs_radial
     end