From 88a3b9b22950ce607dde83768f9378fcd3bd8a36 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Sat, 4 Nov 2023 17:02:36 -0700 Subject: [PATCH] add bent waveguide example, implement pml --- docs/_toc.yml | 1 + docs/julia/waveguide bent.jl | 56 ++++++++++++++++++++---------- src/Maxwell/Waveguide/Waveguide.jl | 31 +++++++++++++---- 3 files changed, 62 insertions(+), 26 deletions(-) diff --git a/docs/_toc.yml b/docs/_toc.yml index 72606832..d7a87001 100644 --- a/docs/_toc.yml +++ b/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 diff --git a/docs/julia/waveguide bent.jl b/docs/julia/waveguide bent.jl index 2a38bf04..e921565f 100644 --- a/docs/julia/waveguide bent.jl +++ b/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) # %% diff --git a/src/Maxwell/Waveguide/Waveguide.jl b/src/Maxwell/Waveguide/Waveguide.jl index 70bffb84..946eb4b0 100644 --- a/src/Maxwell/Waveguide/Waveguide.jl +++ b/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