Skip to content

Commit

Permalink
Merge pull request #103 from lucasgrjn/fix_cylindrical
Browse files Browse the repository at this point in the history
Fix cylindrical mode solver (signs + H-field calculations)
  • Loading branch information
HelgeGehring authored Dec 7, 2023
2 parents 7b3486e + 23a9ea5 commit 19fde20
Show file tree
Hide file tree
Showing 4 changed files with 57 additions and 12 deletions.
1 change: 1 addition & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,7 @@ Something missing? Feel free to open an [issue](https://github.com/HelgeGehring/
- Niko Savola (Google, Aalto University)
- Rouven Glauert (Idalab)
- Markus DeMartini (Google)
- Lucas Grosjean (Google, Femto-ST Institute)

Happy about every form of contribution -
pull requests, feature requests, issues, questions, ... :)
11 changes: 9 additions & 2 deletions docs/julia/waveguide_bent.jl
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,7 @@ radiuss = 1:0.5:5
wg_width = 0.5
sim_right = 4
sim_bottom = 4
pml_thickness = 3
neffs = ComplexF64[]
for radius in radiuss
write_mesh(
Expand All @@ -108,8 +109,14 @@ for radius in radiuss
ε(tag) = Dict(get_tag_from_name(labels, u) => v for (u, v) in epsilons)[tag]

τ = CellField(get_face_tag(labels, num_cell_dims(model)), Ω)
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
# pml_x = x -> 0.1 * max(0, x[1] - (radius + wg_width / 2 + sim_right))^2
pml_order = 3
pml_x =
x ->
1 / pml_order *
max(0, (x[1] - (radius + wg_width / 2 + sim_right)) / pml_thickness)^pml_order *
5
pml_y = x -> 0.0

modes = calculate_modes(model, ε τ, λ = 1.55, order = 1)
println(n_eff(modes[1]))
Expand Down
21 changes: 16 additions & 5 deletions docs/julia/waveguide_bent_coupling.jl
Original file line number Diff line number Diff line change
Expand Up @@ -84,13 +84,13 @@ function write_mesh(;
meshwell.polySurface.PolySurface(
model = model,
polygons = clip_by_rect(env, -np.inf, -np.inf, np.inf, 0),
physical_name = "clad",
physical_name = "box",
mesh_order = 2,
),
meshwell.polySurface.PolySurface(
model = model,
polygons = clip_by_rect(env, -np.inf, 0, np.inf, np.inf),
physical_name = "box",
physical_name = "clad",
mesh_order = 2,
),
]
Expand All @@ -115,6 +115,8 @@ radiuss = 10 # [10.,15.,20.]
wg_width = 0.55
sim_right = 1
sim_bottom = 1
pml_thickness = 3
distance = 0.8
overlaps = Float64[]
Δneff_cylindrical = Float64[]
for radius in radiuss
Expand All @@ -123,15 +125,24 @@ for radius in radiuss
wg_width = wg_width,
sim_right = sim_right,
sim_bottom = sim_bottom,
distance = 0.8,
distance = distance,
)
model = GmshDiscreteModel("mesh.msh")
Ω = Triangulation(model)
labels = get_face_labeling(model)

τ = CellField(get_face_tag(labels, num_cell_dims(model)), Ω)
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
# pml_x = x -> 0.1 * max(0, x[1] - (radius + wg_width / 2 + sim_right))^2
pml_order = 3
pml_x =
x ->
1 / pml_order *
max(
0,
(x[1] - (radius + wg_width / 2 + distance / 2 + sim_right)) / pml_thickness,
)^pml_order *
5
pml_y = x -> 0.0


epsilons = ["core1" => 3.48^2, "core2" => 3.48^2, "box" => 1.46^2, "clad" => 1.46^2]
Expand Down
36 changes: 31 additions & 5 deletions src/Maxwell/Waveguide/Waveguide.jl
Original file line number Diff line number Diff line change
Expand Up @@ -14,13 +14,19 @@ export E, H, power, overlap, te_fraction, tm_fraction, coupling_coefficient, per
export frequency, n_eff, n_eff_cylidrical, ω, λ
export plot_mode, plot_field, write_mode_to_vtk

struct Pml_funcs
pml_f::CellField
pml_g::CellField
end

struct Mode
k0::ComplexF64
k::ComplexF64
E::FEFunction
order::Int
ε::CellField
radius::Float64
pml_funcs::Union{Pml_funcs,Nothing}
end

Base.show(io::IO, mode::Mode) = print(io, "Mode(k=$(mode.k), n_eff=$(n_eff(mode)))")
Expand All @@ -31,7 +37,7 @@ n_eff(mode::Mode) =
n_eff_cylidrical(mode::Mode) =
mode.radius != Inf ? mode.k / mode.k0 : error("Not a cylindrical mode")
ω(mode::Mode) = mode.k0 * ustrip(c_0)
frequency(mode::Mode) = 2π * ω(mode)
frequency(mode::Mode) = ω(mode) / (2π)
λ(mode::Mode) = 2π / mode.k0
function E(mode::Mode)
if mode.radius == Inf
Expand All @@ -45,14 +51,32 @@ end
function H(mode::Mode)
if mode.radius == Inf
-1im / ustrip(μ_0) / ω(mode) * (
(1im * mode.k * mode.E[1] - (mode.E[2]))
(1im * mode.k * mode.E[1] + (mode.E[2]))
TensorValue([0.0 1.0 0.0; -1.0 0.0 0.0]) +
curl(mode.E[1]) * VectorValue(0.0, 0.0, 1.0)
)
else
-1im / ustrip(μ_0) / ω(mode) * (
(1im * mode.k * mode.E[1] - mode.radius * (mode.E[2]))
TensorValue([0.0 1.0 0.0; -1.0 0.0 0.0]) * (x -> 1 / x[1]) +
s_r = 1 - 1im * gradient(mode.pml_funcs.pml_f) VectorValue(1, 0)
s_y = 1 - 1im * gradient(mode.pml_funcs.pml_g) VectorValue(0, 1)
s_θ = 1

γ_r = 1
γ_y = 1
γ_θ = 1 - 1im * mode.pml_funcs.pml_f / (x -> x[1])

d_r = γ_r * s_r / (γ_y * s_y * γ_θ * s_θ)
d_y = γ_y * s_y / (γ_r * s_r * γ_θ * s_θ)
d_θ = γ_θ * s_θ / (γ_r * s_r * γ_y * s_y)

pml_operator = (
d_r TensorValue(1, 0, 0, 0, 0, 0, 0, 0, 0) +
d_y TensorValue(0, 0, 0, 0, 1, 0, 0, 0, 0) +
d_θ TensorValue(0, 0, 0, 0, 0, 0, 0, 0, 1)
)

-1im / ustrip(μ_0) / ω(mode) * pml_operator (
(1im * mode.k * mode.E[1] + mode.radius * (mode.E[2]))
TensorValue([0.0 1.0 0.0; -1.0 0.0 0.0]) * (x -> 1 / x[1]) -
curl(mode.E[1]) * VectorValue(0.0, 0.0, 1.0)
)
end
Expand Down Expand Up @@ -204,13 +228,15 @@ function calculate_modes(
order,
ε,
radius,
(radius == Inf ? nothing : Pml_funcs(pml_f, pml_g)),
),
),
),
),
order,
ε,
radius,
(radius == Inf ? nothing : Pml_funcs(pml_f, pml_g)),
) for (k2, E) in zip(vals, eachcol(vecs))
]
end
Expand Down

0 comments on commit 19fde20

Please sign in to comment.