Selecting a mode based on overlap with a structure or a bounding box

docs/_toc.yml

@@ -18,6 +18,7 @@ parts:
     - file: photonics/examples/bent_waveguide.py
     - file: photonics/examples/vary_width.py
     - file: photonics/examples/vary_wavelength.py
+    - file: photonics/examples/selecting_modes.py
     - file: photonics/examples/fiber_overlap.py
     - file: photonics/examples/metal_heater_phase_shifter.py
     - file: photonics/examples/metal_heater_phase_shifter_transient.py

docs/photonics/examples/selecting_modes.py

@@ -0,0 +1,262 @@
+# ---
+# jupyter:
+#   jupytext:
+#     custom_cell_magics: kql
+#     formats: py:percent,ipynb
+#     text_representation:
+#       extension: .py
+#       format_name: percent
+#       format_version: '1.3'
+#       jupytext_version: 1.11.2
+#   kernelspec:
+#     display_name: env_3.11
+#     language: python
+#     name: python3
+# ---
+
+# %% [markdown]
+# # Selecting modes in the mode solver
+
+# Sometimes we have structures where the mode of interest is
+# not the mode with the highest effective index. There are a few
+# ways to select modes of interest in femwell
+
+# %% tags=["hide-input"]
+from collections import OrderedDict
+
+import matplotlib.pyplot as plt
+import numpy as np
+import scipy.constants
+from scipy.constants import epsilon_0, speed_of_light
+from scipy.integrate import solve_ivp
+from shapely.geometry import Polygon
+from skfem import Basis, ElementTriP0, Mesh
+from skfem.io import from_meshio
+
+from femwell.mesh import mesh_from_OrderedDict
+from femwell.mode_solver import (
+    argsort_modes_by_power_in_elements,
+    calculate_coupling_coefficient,
+    calculate_hfield,
+    calculate_overlap,
+    compute_modes,
+    plot_mode,
+)
+from femwell.utils import inside_bbox
+
+# %% [markdown]
+
+# We will use as an example a system with a Si and a SiN sections.
+# This could happen, for example, in a system where we are trying
+# to heat a SiN waveguide with a Si resistor
+
+
+# %% tags=["remove-stderr", "hide-input"]
+
+w_sim = 6
+h_clad = 2
+h_box = 2
+w_sin = 1
+w_si = 0.4
+gap = 1.0
+h_sin = 0.4
+h_si = 0.22
+
+wavelength = 1.55
+k0 = 2 * np.pi / wavelength
+
+polygons = OrderedDict(
+    sin=Polygon(
+        [
+            (-w_sin - gap / 2, 0),
+            (-w_sin - gap / 2, h_sin),
+            (-gap / 2, h_sin),
+            (-gap / 2, 0),
+        ]
+    ),
+    si=Polygon(
+        [
+            (w_si + gap / 2, 0),
+            (w_si + gap / 2, h_si),
+            (gap / 2, h_si),
+            (gap / 2, 0),
+        ]
+    ),
+    clad=Polygon(
+        [
+            (-w_sim / 2, 0),
+            (-w_sim / 2, h_clad),
+            (w_sim / 2, h_clad),
+            (w_sim / 2, 0),
+        ]
+    ),
+    box=Polygon(
+        [
+            (-w_sim / 2, 0),
+            (-w_sim / 2, -h_box),
+            (w_sim / 2, -h_box),
+            (w_sim / 2, 0),
+        ]
+    ),
+)
+
+resolutions = dict(
+    sin={"resolution": 0.03, "distance": 1},
+    si={"resolution": 0.03, "distance": 1},
+)
+
+mesh = from_meshio(
+    mesh_from_OrderedDict(polygons, resolutions, filename="mesh.msh", default_resolution_max=0.2)
+)
+mesh.draw().show()
+
+# %%
+
+basis0 = Basis(mesh, ElementTriP0(), intorder=4)
+
+epsilon = basis0.zeros() + 1.444**2
+epsilon[basis0.get_dofs(elements=("si"))] = 3.4777**2
+epsilon[basis0.get_dofs(elements=("sin"))] = 1.973**2
+
+
+# %% [markdown]
+
+# ## 0. Directly using femwell
+
+# If we use `find_modes`, these are the modes we get:
+
+# %%
+
+# basis0.plot(epsilon, colorbar=True).show()
+lams, basis, xs = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=4)
+
+
+plot_mode(basis, np.real(xs[0]), direction="x")
+plt.show()
+plot_mode(basis, np.real(xs[1]), direction="x")
+plt.show()
+plot_mode(basis, np.real(xs[2]), direction="x")
+plt.show()
+plot_mode(basis, np.real(xs[3]), direction="x")
+plt.show()
+
+print(f"The effective index of the SiN mode is {np.real(lams[2])}")
+
+# %% [markdown]
+
+# We can see how to get the SiN mode (which is the mode of
+# interest for us) we need to go to the third mode found by femwell.
+
+# Are there easier ways to get the SiN modes? Yes!
+
+# %% [markdown]
+
+# ## 1. Hack (not 100% accurate): Erasing the Si waveguide
+
+# One thing we can do to find the SiN mode is to "erase" the Si
+# waveguide, or in other words assign the refractive index of SiO2
+# to the Si waveguide.
+
+# Of course, this is in general not desired, because this way we are
+# missing the effect of the presence of the Si waveguide.
+
+# thi smight not be an issue in this example but there's many
+# examples where this is not an acceptable option.
+
+# %%
+
+epsilon = basis0.zeros() + 1.444**2
+epsilon[basis0.get_dofs(elements=("si"))] = 1.444**2
+epsilon[basis0.get_dofs(elements=("sin"))] = 1.973**2
+
+lams, basis, xs = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=2)
+
+
+plot_mode(basis, np.real(xs[0]), direction="x")
+plt.show()
+plot_mode(basis, np.real(xs[1]), direction="x")
+plt.show()
+
+print(f"The effective index of the SiN mode is {np.real(lams[0])}")
+
+# %% [markdown]
+# ## 2. Giving a guess effective index
+
+# We can use the `n_guess` parameter to `compute_modes` to
+# select modes close to that effective index.
+
+# This is great, but of course we need to know what's that guess
+# effective index. The way to do that would be to use option 1 above
+# and then use that as the n_guess.
+
+# %%
+epsilon = basis0.zeros() + 1.444**2
+epsilon[basis0.get_dofs(elements=("si"))] = 3.4777**2
+epsilon[basis0.get_dofs(elements=("sin"))] = 1.973**2
+
+lams, basis, xs = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=2, n_guess=1.62)
+
+
+plot_mode(basis, np.real(xs[0]), direction="x")
+plt.show()
+plot_mode(basis, np.real(xs[1]), direction="x")
+plt.show()
+
+print(f"The effective index of the SiN mode is {np.real(lams[1])}")
+
+# %% [markdown]
+
+# You can see how using `n_guess` can still give the wrong mode!
+
+# %% [markdown]
+
+# ## 3. Using `argsort_modes_by_power_in_elements`
+
+# This allows to choose a mode that has the biggest overlap with
+# a given structure.
+
+# There are two main ways to specify the structure:
+# 1. Using the name of the polygon of interest
+# 2. Giving a square bounding box of coordinates
+
+# You can also give it directly the selection_basis of the
+# are of interest.
+
+# A requirement for using `argsort_modes_by_power_in_elements` is to
+# calculate the H field of the found modes.
+
+# %%
+
+lams, basis, xs = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=4)
+
+# Calculate H field
+H_modes = list()
+for i, E in enumerate(xs):
+    H_modes.append(
+        calculate_hfield(
+            basis,
+            E,
+            lams[i] * (2 * np.pi / wavelength),
+            omega=2 * np.pi / 1.55 * scipy.constants.speed_of_light,
+        )
+    )
+
+# Option 1: using an element name
+
+ind_mode = argsort_modes_by_power_in_elements(basis, xs, H_modes, elements="sin")[0]
+
+plot_mode(basis, np.real(xs[ind_mode]), direction="x")
+plt.show()
+
+# Option 2: using bounding box
+
+# Format: [xmin, ymin, xmax, ymax]
+bbox = [-2, 0, 0, 0.4]
+
+elements = inside_bbox(bbox)
+ind_mode = argsort_modes_by_power_in_elements(basis, xs, H_modes, elements)[0]
+
+plot_mode(basis, np.real(xs[ind_mode]), direction="x")
+plt.show()
+
+# %%

femwell/mode_solver.py

@@ -326,6 +326,24 @@ def plot_mode(basis, mode, plot_vectors=False, colorbar=True, title="E", directi
     return fig, axs
 
 
+def argsort_modes_by_power_in_elements(mode_basis, E_modes, H_modes, elements):
+    """Sorts the modes in the "modes" list by the power contained
+    within the given elements.
+
+    Returns:
+        the indices sorted from highest to lowest power.
+    """
+
+    selection_basis = Basis(mode_basis.mesh, mode_basis.elem, elements=elements)
+
+    overlaps = [
+        calculate_overlap(selection_basis, E_f, H_f, selection_basis, E_f, H_f)
+        for E_f, H_f in zip(E_modes, H_modes)
+    ]
+
+    return np.argsort(np.abs(overlaps))[::-1]
+
+
 if __name__ == "__main__":
     from collections import OrderedDict
 
@@ -384,7 +402,7 @@ def plot_mode(basis, mode, plot_vectors=False, colorbar=True, title="E", directi
     # basis0.plot(epsilon, colorbar=True).show()
 
     lams, basis, xs = compute_modes(
-        basis0, epsilon, wavelength=1.55, mu_r=1, num_modes=6, order=2, radius=3
+        basis0, epsilon, wavelength=1.55, mu_r=1, num_modes=6, order=2, radius=3, solver="scipy"
     )
     print(lams)
 
@@ -403,6 +421,8 @@ def plot_mode(basis, mode, plot_vectors=False, colorbar=True, title="E", directi
     plt.show()
 
     integrals = np.zeros((len(lams),) * 2, dtype=complex)
+    H_modes = list()
+
     for i in range(len(lams)):
         for j in range(len(lams)):
             E_i = xs[i]
@@ -420,7 +440,40 @@ def plot_mode(basis, mode, plot_vectors=False, colorbar=True, title="E", directi
                 omega=2 * np.pi / 1.55 * scipy.constants.speed_of_light,
             )
             integrals[i, j] = calculate_overlap(basis, E_i, H_i, basis, E_j, H_j)
+        H_modes.append(H_i)
 
     plt.imshow(np.real(integrals))
     plt.colorbar()
     plt.show()
+
+    # Create basis to select a certain simulation extent
+    def sel_fun(x):
+        print(x)
+        return (x[0] < 0) * (x[0] > -1) * (x[1] > 0) * (x[1] < 0.5)
+
+    selection_basis = Basis(
+        basis.mesh,
+        basis.elem,
+        # elements = lambda x: x[0] < 0 and x[0] > -1 and x[1] > 0 and x[1] < 0.5
+        elements=lambda x: sel_fun(x),
+    )
+
+    print(
+        select_mode_by_overlap(
+            mode_basis=basis,
+            E_modes=xs,
+            H_modes=H_modes,
+            elements_list=["core"],
+            selection_basis=None,
+        )
+    )
+
+    print(
+        select_mode_by_overlap(
+            mode_basis=basis,
+            E_modes=xs,
+            H_modes=H_modes,
+            elements_list=None,
+            selection_basis=selection_basis,
+        )
+    )

femwell/utils.py

@@ -10,6 +10,7 @@
 else:
     from scipy.sparse import spmatrix
 
+from skfem import Basis
 from skfem.utils import CondensedSystem, bmat
 
 
@@ -92,3 +93,18 @@ def mpc_symmetric(
             lambda x: np.concatenate((x, T @ x[len(U) :] + g)),
         ),
     )
+
+
+def inside_bbox(bbox):
+    """
+    Creates a selection function that is True for elements in the
+    given bounding box coordinates.
+
+    Args:
+        bbox: 4 element list with [xmin, ymin, xmax, ymax]
+    """
+
+    def sel_fun(x):
+        return (x[0] < bbox[2]) * (x[0] > bbox[0]) * (x[1] > bbox[1]) * (x[1] < bbox[3])
+
+    return sel_fun