cylindric mode solver seems working

docs/julia/waveguide bent.jl

@@ -0,0 +1,113 @@
+# ---
+# jupyter:
+#   jupytext:
+#     custom_cell_magics: kql
+#     formats: jl:percent,ipynb
+#     text_representation:
+#       extension: .jl
+#       format_name: percent
+#       format_version: '1.3'
+#       jupytext_version: 1.11.2
+#   kernelspec:
+#     display_name: base
+#     language: julia
+#     name: julia-1.9
+# ---
+
+# %% [markdown]
+# # Mode solving
+
+# %% [markdown]
+# ```{caution}
+# **This example is under construction**
+# As Julia-Dicts are not ordered, the mesh might become incorrect when adjusted (for now, better do the meshing in python)
+# ```
+
+# %% tags=["hide-output", "thebe-init"]
+using PyCall
+np = pyimport("numpy")
+shapely = pyimport("shapely")
+shapely.affinity = pyimport("shapely.affinity")
+clip_by_rect = pyimport("shapely.ops").clip_by_rect
+OrderedDict = pyimport("collections").OrderedDict
+mesh_from_OrderedDict = pyimport("femwell.mesh").mesh_from_OrderedDict
+
+radius = 100000
+wg_width = 0.5
+wg_thickness = 0.22
+sim_width = 3
+sim_height = 4
+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,
+)
+
+polygons = OrderedDict(
+    core = core,
+    box = clip_by_rect(env, -np.inf, -np.inf, np.inf, 0),
+    clad = clip_by_rect(env, -np.inf, 0, np.inf, np.inf),
+)
+
+resolutions = Dict(
+    "core" => Dict("resolution" => 0.03, "distance" => 0.5),
+    "slab" => Dict("resolution" => 0.03, "distance" => 0.5),
+)
+
+mesh = mesh_from_OrderedDict(
+    polygons,
+    resolutions,
+    default_resolution_max = 1,
+    filename = "mesh.msh",
+)
+
+# %% tags=["remove-stderr", "hide-output", "thebe-init"]
+using Gridap
+using Gridap.Geometry
+using Gridap.Visualization
+using Gridap.ReferenceFEs
+using GridapGmsh
+using GridapMakie, CairoMakie
+
+using Femwell.Maxwell.Waveguide
+
+CairoMakie.inline!(true)
+
+# %% tags=["remove-stderr"]
+model = GmshDiscreteModel("mesh.msh")
+Ω = Triangulation(model)
+#fig = plot(Ω)
+#fig.axis.aspect=DataAspect()
+#wireframe!(Ω, color=:black, linewidth=1)
+#display(fig)
+
+labels = get_face_labeling(model)
+
+epsilons = ["core" => 3.5^2, "box" => 1.444^2, "clad" => 1.44^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]))
+# write_mode_to_vtk("mode", modes[2])
+
+plot_mode(modes[2], absolute = true)
+#plot_mode(modes[2])
+modes
+
+# %% [markdown]
+# For comparison, we the effective refractive index of a straight waveguide:
+
+# %% 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)
+
+# %%

docs/julia/waveguide.jl

@@ -32,18 +32,12 @@ clip_by_rect = pyimport("shapely.ops").clip_by_rect
 OrderedDict = pyimport("collections").OrderedDict
 mesh_from_OrderedDict = pyimport("femwell.mesh").mesh_from_OrderedDict
 
-radius = 5
 wg_width = 0.5
 wg_thickness = 0.22
 slab_width = 3
 slab_thickness = 0.11
-core = shapely.geometry.box(radius - wg_width / 2, 0, radius + wg_width / 2, wg_thickness)
-slab = shapely.geometry.box(
-    radius - slab_width / 2,
-    0,
-    radius + slab_width / 2,
-    slab_thickness,
-)
+core = shapely.geometry.box(-wg_width / 2, 0, wg_width / 2, wg_thickness)
+slab = shapely.geometry.box(slab_width / 2, 0, slab_width / 2, slab_thickness)
 env = shapely.affinity.scale(core.buffer(3, resolution = 8), xfact = 0.5)
 
 polygons = OrderedDict(
@@ -87,16 +81,16 @@ model = GmshDiscreteModel("mesh.msh")
 
 labels = get_face_labeling(model)
 
-epsilons = ["core" => 3.5^2 + 0.001im, "slab" => 3.5^2, "box" => 1.444^2, "clad" => 1.0^2]
+epsilons = ["core" => 3.5^2, "slab" => 1.44^2, "box" => 1.444^2, "clad" => 1.44^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 = 1, order = 1, radius = 5)
+modes = calculate_modes(model, ε ∘ τ, λ = 1.55, num = 1, order = 1)
 println(n_eff(modes[1]))
-write_mode_to_vtk("mode", modes[1])
+# write_mode_to_vtk("mode", modes[2])
 
 plot_mode(modes[1])
 #plot_mode(modes[2])
@@ -117,7 +111,8 @@ println(perturbed_neff(modes[1], ε_p ∘ τ))
 
 # %% tags=["remove-stderr"]
 
-#modes_p = calculate_modes(model, ε ∘ τ + ε_p ∘ τ, λ = 1.55, num = 2, order = 1)
-#println(n_eff(modes_p[1]))
+modes_p = calculate_modes(model, ε ∘ τ + ε_p ∘ τ, λ = 1.55, num = 2, order = 1)
+println(n_eff(modes_p[1]))
+plot_mode(modes_p[1])
 
 # %%

src/Maxwell/Waveguide/Waveguide.jl

@@ -99,46 +99,58 @@ function calculate_modes(
     Ω = Triangulation(model)
     dΩ = Measure(Ω, 2 * order)
 
-    μ_r = 1
-    radius_factor(x) = (1 + x[1] / radius)^2
-    twothree = TensorValue([1 0 0; 0 1 0])
-    lhs((u1, u2), (v1, v2)) =
-        ∫(
-            1 / μ_r * (curl(u1) ⊙ curl(v1) / k0^2 + ∇(u2) ⊙ v1) +
-            radius_factor * (
-                -u1 ⋅ twothree ⋅ ε ⋅ transpose(twothree) ⋅ v1 +
-                u1 ⋅ twothree ⋅ ε ⋅ transpose(twothree) ⋅ ∇(v2) -
-                u2 * VectorValue(0, 0, 1) ⋅ ε ⋅ VectorValue(0, 0, 1) * v2 * k0^2
-            ),
-        )dΩ
-    rhs((u1, u2), (v1, v2)) = ∫(-1 / μ_r * u1 ⊙ v1 / k0^2)dΩ
-
-    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),
-        )dΩ
-
-    rhs_radial((u1, u2), (v1, v2)) =
-        ∫((radius / r * u1 ⋅ v1) - (radius / r * u1 ⋅ gradient(v2)))dΩ
-
     epsilons = ε(get_cell_points(Measure(Ω, 1)))
     k0_guess =
         isnothing(k0_guess) ? k0^2 * maximum(maximum.(maximum.(real.(epsilons)))) * 1.1 :
         k0_guess
 
+    lhs, rhs = if radius == Inf
+        μ_r = 1
+        twothree = TensorValue([1 0 0; 0 1 0])
+        lhs_straight((u1, u2), (v1, v2)) =
+            ∫(
+                1 / μ_r * (curl(u1) ⊙ curl(v1) / k0^2 + ∇(u2) ⊙ v1) + (
+                    -u1 ⋅ twothree ⋅ ε ⋅ transpose(twothree) ⋅ v1 +
+                    u1 ⋅ twothree ⋅ ε ⋅ transpose(twothree) ⋅ ∇(v2) -
+                    u2 * VectorValue(0, 0, 1) ⋅ ε ⋅ VectorValue(0, 0, 1) * v2 * k0^2
+                ),
+            )dΩ
+        rhs_straight((u1, u2), (v1, v2)) = ∫(-1 / μ_r * u1 ⊙ v1 / k0^2)dΩ
+
+        lhs_straight, rhs_straight
+    else
+        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),
+            )dΩ
+
+        rhs_radial((u1, u2), (v1, v2)) =
+            ∫((radius / r * u1 ⋅ v1) - (radius / r * u1 ⋅ gradient(v2)))dΩ
+
+        lhs_radial, rhs_radial
+    end
+
+
     assem = Gridap.FESpaces.SparseMatrixAssembler(U, V)
-    A = assemble_matrix(lhs_radial, assem, U, V)
-    B = assemble_matrix(rhs_radial, assem, U, V)
+    A = assemble_matrix(lhs, assem, U, V)
+    B = assemble_matrix(rhs, assem, U, V)
     if all(imag(A.nzval) .== 0) && all(imag(B.nzval) .== 0)
         A, B = real(A), real(B)
     end
     vals, vecs = eigs(A, B, sigma = k0_guess, nev = num)
-    vecs[num_free_dofs(V1)+1:end, :] ./= 1im * sqrt.(vals)' / k0^2
+
+    if radius == Inf
+        vecs[num_free_dofs(V1)+1:end, :] ./= 1im * sqrt.(vals)' / k0^2
+    else
+        vecs[num_free_dofs(V1)+1:end, :] ./= 1im * sqrt.(vals)'
+    end
+
 
     return [
         Mode(
@@ -157,7 +169,13 @@ function plot_field(field)
     display(fig)
 end
 
-function plot_mode(mode::Mode; vertical = false, vectors = false, same_range = false)
+function plot_mode(
+    mode::Mode;
+    vertical = false,
+    vectors = false,
+    same_range = false,
+    absolute = false,
+)
     Ω = get_triangulation(mode.E)
     model = get_active_model(Ω)
     labels = get_face_labeling(model)
@@ -190,8 +208,8 @@ function plot_mode(mode::Mode; vertical = false, vectors = false, same_range = f
             plt = plot!(
                 ax,
                 Ω,
-                real((E(mode) ⋅ vector)),
-                colorrange = (-colorrange, colorrange),
+                (absolute ? abs : real)((E(mode) ⋅ vector)),
+                #colorrange = (-colorrange, colorrange),
             )
             wireframe!(fig[x, y], ∂Ω, color = :black)
             Colorbar(fig[x+!vertical, y+vertical], plt, vertical = vertical)