From 8a56591339af683649cae9a41ffcf0532f29c307 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Sat, 1 Apr 2023 18:12:36 -0700 Subject: [PATCH 01/46] start implementing object oriented mode solver --- femwell/mode_solver.py | 144 +++++++++++++++++++++++------------------ 1 file changed, 81 insertions(+), 63 deletions(-) diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index 4e2e74ed..c03acbe3 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -1,6 +1,10 @@ """Waveguide analysis based on https://doi.org/10.1080/02726340290084012.""" +from dataclasses import dataclass import matplotlib.pyplot as plt import numpy as np +from typing import List +from numpy.typing import NDArray +from pyparsing import col import scipy.constants import scipy.sparse.linalg from scipy.constants import epsilon_0, speed_of_light @@ -23,6 +27,55 @@ from skfem.helpers import cross, curl, dot, grad, inner from skfem.utils import solver_eigen_scipy +@dataclass +class Modes: + modes: List + + def __getitem__(self, idx): + return self.modes[idx] + + def __len__(self): + return len(self.modes) + + def __repr__(self) -> str: + modes = '\n\t' + '\n\t'.join(repr(mode) for mode in self.modes) + '\n' + return f"{self.__class__.__name__}(modes=({modes}))" + +@dataclass(frozen=True) +class Mode: + frequency: float + k: float + basis: Basis + epsilon_r: NDArray + E: NDArray + H: NDArray + + @property + def omega(self): + return 2 * np.pi * self.frequency + + @property + def k0(self): + return self.omega / speed_of_light + + @property + def n_eff(self): + return self.k / self.k0 + + def __repr__(self) -> str: + return f"{self.__class__.__name__}(k: {self.k}, n_eff:{self.n_eff})" + + def calculate_overlap(self, mode): + return calculate_overlap(self.basis, self.E, self.H, mode.basis, mode.E, mode.H) + + def plot(self, field, plot_vectors=False, colorbar=True, direction='y', title='E'): + plot_mode(self.basis, field, plot_vectors=plot_vectors, colorbar=colorbar, title=title, direction=direction) + + def show(self, field, **kwargs): + self.plot(field=field, **kwargs) + plt.show() + + def compute_modes( basis_epsilon_r, @@ -37,6 +90,7 @@ def compute_modes( solver="slepc", normalize=True, cache_path=None, + return_objects=False ): if solver == "scipy": solver = solver_eigen_scipy @@ -109,14 +163,21 @@ def bform(e_t, e_z, v_t, v_z, w): lams[:, np.newaxis] ) # undo the scaling E_3,new = beta * E_3 + hs = [] if normalize: for i, lam in enumerate(lams): H = calculate_hfield( basis, xs[i], np.sqrt(lam), omega=k0 * scipy.constants.speed_of_light ) - xs[i] /= np.sqrt(calculate_overlap(basis, xs[i], H, basis, xs[i], H)) + power = calculate_overlap(basis, xs[i], H, basis, xs[i], H) + xs[i] /= np.sqrt(power) + H /= np.sqrt(power) + hs.append(H) - return np.sqrt(lams)[:num_modes] / k0, basis, xs[:num_modes] + if return_objects: + return Modes(modes=[Mode(frequency=speed_of_light/wavelength,k=np.sqrt(lams[i]), basis=basis, epsilon_r=epsilon_r, E=xs[i], H=hs[i]) for i in range(num_modes)]) + else: + return np.sqrt(lams)[:num_modes] / k0, basis, xs[:num_modes] def calculate_hfield(basis, xs, beta, omega=1): @@ -326,7 +387,7 @@ 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): +def argsort_modes_by_power_in_elements(modes, elements): """Sorts the modes in the "modes" list by the power contained within the given elements. @@ -334,11 +395,11 @@ def argsort_modes_by_power_in_elements(mode_basis, E_modes, H_modes, elements): the indices sorted from highest to lowest power. """ - selection_basis = Basis(mode_basis.mesh, mode_basis.elem, elements=elements) + selection_basis = Basis(modes[0].basis.mesh, modes[0].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) + calculate_overlap(selection_basis, mode.E, mode.H, selection_basis, mode.E, mode.H) + for mode in modes ] return np.argsort(np.abs(overlaps))[::-1] @@ -401,46 +462,22 @@ def argsort_modes_by_power_in_elements(mode_basis, E_modes, H_modes, elements): epsilon[basis0.get_dofs(elements="box")] = 1.444**2 # 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, solver="scipy" + modes = compute_modes( + basis0, epsilon, wavelength=1.55, mu_r=1, num_modes=6, order=2, radius=3 ) - print(lams) + print(modes) - plot_mode(basis, np.real(xs[0])) - plt.show() - plot_mode(basis, np.real(xs[1])) - plt.show() - plot_mode(basis, np.imag(xs[0])) - plt.show() + modes[0].show(np.real(modes[0].E)) + modes[0].show(np.imag(modes[0].E)) - xbs = calculate_hfield(basis, xs[0], lams[0] * (2 * np.pi / 1.55)) + modes[0].show(np.real(modes[0].H)) + modes[0].show(np.imag(modes[0].H)) - plot_mode(basis, np.real(xbs)) - plt.show() - plot_mode(basis, np.imag(xbs)) - plt.show() + integrals = np.zeros((len(modes),) * 2, dtype=complex) - 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] - E_j = xs[j] - H_i = calculate_hfield( - basis, - E_i, - lams[i] * (2 * np.pi / 1.55), - omega=2 * np.pi / 1.55 * scipy.constants.speed_of_light, - ) - H_j = calculate_hfield( - basis, - E_j, - lams[j] * (2 * np.pi / 1.55), - 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) + for i in range(len(modes)): + for j in range(len(modes)): + integrals[i, j] = modes[i].calculate_overlap(modes[j]) plt.imshow(np.real(integrals)) plt.colorbar() @@ -451,29 +488,10 @@ 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, + argsort_modes_by_power_in_elements( + modes = modes, + elements=sel_fun, ) ) From 031e8abefb8f7c8d41d708e0b5a4fc30105003e1 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Sat, 1 Apr 2023 18:15:36 -0700 Subject: [PATCH 02/46] run pre-commit --- femwell/mode_solver.py | 48 +++++++++++++++++++++++++++++------------- 1 file changed, 33 insertions(+), 15 deletions(-) diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index c03acbe3..e59d2617 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -1,12 +1,13 @@ """Waveguide analysis based on https://doi.org/10.1080/02726340290084012.""" from dataclasses import dataclass +from typing import List + import matplotlib.pyplot as plt import numpy as np -from typing import List -from numpy.typing import NDArray -from pyparsing import col import scipy.constants import scipy.sparse.linalg +from numpy.typing import NDArray +from pyparsing import col from scipy.constants import epsilon_0, speed_of_light from skfem import ( Basis, @@ -27,6 +28,7 @@ from skfem.helpers import cross, curl, dot, grad, inner from skfem.utils import solver_eigen_scipy + @dataclass class Modes: modes: List @@ -38,9 +40,10 @@ def __len__(self): return len(self.modes) def __repr__(self) -> str: - modes = '\n\t' + '\n\t'.join(repr(mode) for mode in self.modes) + '\n' + modes = "\n\t" + "\n\t".join(repr(mode) for mode in self.modes) + "\n" return f"{self.__class__.__name__}(modes=({modes}))" + @dataclass(frozen=True) class Mode: frequency: float @@ -68,13 +71,19 @@ def __repr__(self) -> str: def calculate_overlap(self, mode): return calculate_overlap(self.basis, self.E, self.H, mode.basis, mode.E, mode.H) - def plot(self, field, plot_vectors=False, colorbar=True, direction='y', title='E'): - plot_mode(self.basis, field, plot_vectors=plot_vectors, colorbar=colorbar, title=title, direction=direction) + def plot(self, field, plot_vectors=False, colorbar=True, direction="y", title="E"): + plot_mode( + self.basis, + field, + plot_vectors=plot_vectors, + colorbar=colorbar, + title=title, + direction=direction, + ) def show(self, field, **kwargs): self.plot(field=field, **kwargs) plt.show() - def compute_modes( @@ -90,7 +99,7 @@ def compute_modes( solver="slepc", normalize=True, cache_path=None, - return_objects=False + return_objects=False, ): if solver == "scipy": solver = solver_eigen_scipy @@ -171,11 +180,23 @@ def bform(e_t, e_z, v_t, v_z, w): ) power = calculate_overlap(basis, xs[i], H, basis, xs[i], H) xs[i] /= np.sqrt(power) - H /= np.sqrt(power) + H /= np.sqrt(power) hs.append(H) if return_objects: - return Modes(modes=[Mode(frequency=speed_of_light/wavelength,k=np.sqrt(lams[i]), basis=basis, epsilon_r=epsilon_r, E=xs[i], H=hs[i]) for i in range(num_modes)]) + return Modes( + modes=[ + Mode( + frequency=speed_of_light / wavelength, + k=np.sqrt(lams[i]), + basis=basis, + epsilon_r=epsilon_r, + E=xs[i], + H=hs[i], + ) + for i in range(num_modes) + ] + ) else: return np.sqrt(lams)[:num_modes] / k0, basis, xs[:num_modes] @@ -462,9 +483,7 @@ def argsort_modes_by_power_in_elements(modes, elements): epsilon[basis0.get_dofs(elements="box")] = 1.444**2 # basis0.plot(epsilon, colorbar=True).show() - modes = compute_modes( - basis0, epsilon, wavelength=1.55, mu_r=1, num_modes=6, order=2, radius=3 - ) + modes = compute_modes(basis0, epsilon, wavelength=1.55, mu_r=1, num_modes=6, order=2, radius=3) print(modes) modes[0].show(np.real(modes[0].E)) @@ -488,10 +507,9 @@ def sel_fun(x): print(x) return (x[0] < 0) * (x[0] > -1) * (x[1] > 0) * (x[1] < 0.5) - print( argsort_modes_by_power_in_elements( - modes = modes, + modes=modes, elements=sel_fun, ) ) From 84eac9597a998ac5771a5ae7677d76088120fea9 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 3 Apr 2023 13:21:05 -0700 Subject: [PATCH 03/46] add te and tm-fraction --- femwell/mode_solver.py | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index e59d2617..db090159 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -1,5 +1,7 @@ """Waveguide analysis based on https://doi.org/10.1080/02726340290084012.""" from dataclasses import dataclass +from functools import cache, cached_property +from time import time from typing import List import matplotlib.pyplot as plt @@ -65,6 +67,14 @@ def k0(self): def n_eff(self): return self.k / self.k0 + @cached_property + def te_fraction(self): + return calculate_te_frac(self.basis, self.E) + + @cached_property + def tm_fraction(self): + return 1 - calculate_te_frac(self.basis, self.E) + def __repr__(self) -> str: return f"{self.__class__.__name__}(k: {self.k}, n_eff:{self.n_eff})" @@ -483,7 +493,16 @@ def argsort_modes_by_power_in_elements(modes, elements): epsilon[basis0.get_dofs(elements="box")] = 1.444**2 # basis0.plot(epsilon, colorbar=True).show() - modes = compute_modes(basis0, epsilon, wavelength=1.55, mu_r=1, num_modes=6, order=2, radius=3) + modes = compute_modes( + basis0, + epsilon, + wavelength=1.55, + mu_r=1, + num_modes=6, + order=2, + radius=3, + return_objects=True, + ) print(modes) modes[0].show(np.real(modes[0].E)) From 0c1020e76a5def42b082508afbb19b5d910c4eb0 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 3 Apr 2023 15:06:18 -0700 Subject: [PATCH 04/46] simplify Mode class --- femwell/mode_solver.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index db090159..4668623b 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -73,7 +73,7 @@ def te_fraction(self): @cached_property def tm_fraction(self): - return 1 - calculate_te_frac(self.basis, self.E) + return 1 - self.te_fraction def __repr__(self) -> str: return f"{self.__class__.__name__}(k: {self.k}, n_eff:{self.n_eff})" From 4b49f23cfb4fc8797f6a23a578911639a5353412 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 3 Apr 2023 16:06:15 -0700 Subject: [PATCH 05/46] add wavelength to Mode class --- femwell/mode_solver.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index 4668623b..2cf6d48d 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -63,6 +63,10 @@ def omega(self): def k0(self): return self.omega / speed_of_light + @property + def wavelength(self): + return speed_of_light / self.frequency + @property def n_eff(self): return self.k / self.k0 From 3ec9bdc3e747d97b792b78adc4ff178e6545609c Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 3 Apr 2023 16:18:16 -0700 Subject: [PATCH 06/46] Mode class: add basis_epsilon_r --- femwell/mode_solver.py | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index 2cf6d48d..79637587 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -50,8 +50,9 @@ def __repr__(self) -> str: class Mode: frequency: float k: float - basis: Basis + basis_epsilon_r: Basis epsilon_r: NDArray + basis: Basis E: NDArray H: NDArray @@ -73,7 +74,18 @@ def n_eff(self): @cached_property def te_fraction(self): - return calculate_te_frac(self.basis, self.E) + @Functional + def ex(w): + return np.abs(w.E[0][0]) ** 2 + + @Functional + def ey(w): + return np.abs(w.E[0][1]) ** 2 + + ex_sum = ex.assemble(self.basis, E=self.basis.interpolate(self.E)) + ey_sum = ey.assemble(self.basis, E=self.basis.interpolate(self.E)) + + return ex_sum / (ex_sum + ey_sum) @cached_property def tm_fraction(self): @@ -203,8 +215,9 @@ def bform(e_t, e_z, v_t, v_z, w): Mode( frequency=speed_of_light / wavelength, k=np.sqrt(lams[i]), - basis=basis, + basis_epsilon_r=basis_epsilon_r, epsilon_r=epsilon_r, + basis=basis, E=xs[i], H=hs[i], ) @@ -508,6 +521,7 @@ def argsort_modes_by_power_in_elements(modes, elements): return_objects=True, ) print(modes) + print(modes[0].te_fraction) modes[0].show(np.real(modes[0].E)) modes[0].show(np.imag(modes[0].E)) From f7091d14c5f2be319b6d4d9b74c2e5ed24dbc20e Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 3 Apr 2023 16:42:10 -0700 Subject: [PATCH 07/46] simplify example --- docs/photonics/examples/bent_waveguide.py | 43 +++++++---------------- 1 file changed, 12 insertions(+), 31 deletions(-) diff --git a/docs/photonics/examples/bent_waveguide.py b/docs/photonics/examples/bent_waveguide.py index d20eaf76..23f46c0a 100644 --- a/docs/photonics/examples/bent_waveguide.py +++ b/docs/photonics/examples/bent_waveguide.py @@ -19,7 +19,6 @@ import matplotlib.pyplot as plt import numpy as np -import scipy.constants from shapely import box from shapely.ops import clip_by_rect from skfem import Basis, ElementDG, ElementTriP1 @@ -27,12 +26,7 @@ from tqdm import tqdm from femwell.mesh import mesh_from_OrderedDict -from femwell.mode_solver import ( - calculate_hfield, - calculate_overlap, - compute_modes, - plot_mode, -) +from femwell.mode_solver import compute_modes # - @@ -93,14 +87,8 @@ # and for mode overlap calculations between straight and bent waveguides. # + -lams_straight, basis_straight, xs_straight = compute_modes( - basis0, epsilon, wavelength=wavelength, num_modes=1, order=2, radius=np.inf -) -H_straight = calculate_hfield( - basis_straight, - xs_straight[0], - 2 * np.pi / wavelength * lams_straight[0], - omega=2 * np.pi / wavelength * scipy.constants.speed_of_light, +modes_straight = compute_modes( + basis0, epsilon, wavelength=wavelength, num_modes=1, order=2, radius=np.inf, return_objects=True ) # - @@ -112,9 +100,9 @@ radiuss = np.linspace(40, 5, 21) radiuss_lams = [] overlaps = [] -lam_guess = lams_straight[0] +lam_guess = modes_straight[0].n_eff for radius in tqdm(radiuss): - lams, basis, xs = compute_modes( + modes = compute_modes( basis0, epsilon, wavelength=wavelength, @@ -123,19 +111,12 @@ radius=radius, n_guess=lam_guess, solver="scipy", + return_objects=True, ) - lam_guess = lams[0] - H_bent = calculate_hfield( - basis_straight, - xs[0], - 2 * np.pi / wavelength * lams[0], - omega=2 * np.pi / wavelength * scipy.constants.speed_of_light, - ) - radiuss_lams.append(lams[0]) + lam_guess = modes[0].n_eff + radiuss_lams.append(modes[0].n_eff) - overlaps.append( - calculate_overlap(basis_straight, xs[0], H_bent, basis_straight, xs_straight[0], H_straight) - ) + overlaps.append(modes_straight[0].calculate_overlap(modes[0])) # - # And now we plot it! @@ -164,7 +145,7 @@ # Here we show only the real part of the mode. # + tags=["hide-input"] -for i, lam in enumerate(lams): - print(f"Effective refractive index: {lam:.14f}") - plot_mode(basis, xs[i].real, colorbar=True, direction="x") +for mode in modes: + print(f"Effective refractive index: {mode.n_eff:.14f}") + mode.plot(mode.E.real, colorbar=True, direction="x") plt.show() From 0e357b928e4d2940a74fe5f945eaf536d2799d6d Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 3 Apr 2023 17:20:26 -0700 Subject: [PATCH 08/46] add calculate_coupling_coefficient to Mode class --- femwell/mode_solver.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index 79637587..ee767974 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -97,6 +97,11 @@ def __repr__(self) -> str: def calculate_overlap(self, mode): return calculate_overlap(self.basis, self.E, self.H, mode.basis, mode.E, mode.H) + def calculate_coupling_coefficient(self, mode, delta_epsilon): + return calculate_coupling_coefficient( + self.basis_epsilon_r, delta_epsilon, self.basis, self.E, mode.E + ) + def plot(self, field, plot_vectors=False, colorbar=True, direction="y", title="E"): plot_mode( self.basis, From 499e5742b8601261bd1f77c984f2b78d7c6429be Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 3 Apr 2023 17:20:56 -0700 Subject: [PATCH 09/46] adjust coupled_mode_theory example to object-oriented mode solver --- .../photonics/examples/coupled_mode_theory.py | 120 ++++++++---------- 1 file changed, 51 insertions(+), 69 deletions(-) diff --git a/docs/photonics/examples/coupled_mode_theory.py b/docs/photonics/examples/coupled_mode_theory.py index b6b25b0b..589b1b5a 100644 --- a/docs/photonics/examples/coupled_mode_theory.py +++ b/docs/photonics/examples/coupled_mode_theory.py @@ -28,6 +28,7 @@ # %% tags=["hide-input"] from collections import OrderedDict +from itertools import chain import matplotlib.pyplot as plt import numpy as np @@ -38,13 +39,7 @@ from skfem.io import from_meshio from femwell.mesh import mesh_from_OrderedDict -from femwell.mode_solver import ( - calculate_coupling_coefficient, - calculate_hfield, - calculate_overlap, - compute_modes, - plot_mode, -) +from femwell.mode_solver import compute_modes # %% [markdown] # Let's set up the geometry! @@ -119,17 +114,20 @@ epsilon = basis0.zeros() + 1.444**2 epsilon[basis0.get_dofs(elements=("core_1", "core_2"))] = 3.4777**2 # basis0.plot(epsilon, colorbar=True).show() -lams_both, basis, xs_both = compute_modes( - basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=2 +modes_both = compute_modes( + basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=2, return_objects=True +) +modes_both[0].show(modes_both[0].E.real, direction="x") +modes_both[1].show(modes_both[1].E.real, direction="x") +print( + "Refractive index of symmetric and assymetric mode:", + modes_both[0].n_eff, + ", ", + modes_both[1].n_eff, ) -plot_mode(basis, np.real(xs_both[0]), direction="x") -plt.show() -plot_mode(basis, np.real(xs_both[1]), direction="x") -plt.show() -print("Refractive index of symmetric and assymetric mode:", lams_both) # https://www.fiberoptics4sale.com/blogs/wave-optics/directional-couplers print( - f"Maximum power transfer after {np.pi / (2 * np.pi / wavelength * np.real(lams_both[0] - lams_both[1]))} um prop length" + f"Maximum power transfer after {np.pi / (2 * np.pi / wavelength * np.real(modes_both[0].n_eff - modes_both[1].n_eff))} um prop length" ) # %% [markdown] @@ -139,18 +137,20 @@ epsilon = basis0.zeros() + 1.444**2 epsilon[basis0.get_dofs(elements="core_1")] = 3.4777**2 # basis0.plot(epsilon, colorbar=True).show() -lams_1, basis, xs_1 = compute_modes(basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=1) -print("Effective refractive index of the mode of the first waveguide", lams_1) -plot_mode(basis, np.real(xs_1[0]), direction="x") -plt.show() +modes_1 = compute_modes( + basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=1, return_objects=True +) +print("Effective refractive index of the mode of the first waveguide", modes_1[0].n_eff) +modes_1[0].show(modes_1[0].E.real, direction="x") epsilon_2 = basis0.zeros() + 1.444**2 epsilon_2[basis0.get_dofs(elements="core_2")] = 3.4777**2 # basis0.plot(epsilon_2, colorbar=True).show() -lams_2, basis, xs_2 = compute_modes(basis0, epsilon_2, wavelength=wavelength, mu_r=1, num_modes=1) -print("Effective refractive index of the mode of the second waveguide", lams_2) -plot_mode(basis, np.real(xs_2[0]), direction="x") -plt.show() +modes_2 = compute_modes( + basis0, epsilon_2, wavelength=wavelength, mu_r=1, num_modes=1, return_objects=True +) +print("Effective refractive index of the mode of the second waveguide", modes_2[0].n_eff) +modes_2[0].show(modes_2[0].E.real, direction="x") # %% length = 200 @@ -158,41 +158,27 @@ # %% epsilons = [epsilon, epsilon_2] -modes = [(lam, x, 0) for lam, x in zip(lams_1, xs_1)] + [ - (lam, x, 1) for lam, x in zip(lams_2, xs_2) -] - -overlap_integrals = np.zeros((len(modes), len(modes)), dtype=complex) -for i, (lam_i, E_i, epsilon_i) in enumerate(modes): - for j, (lam_j, E_j, epsilon_j) in enumerate(modes): - H_i = calculate_hfield( - basis, - E_i, - lam_i * (2 * np.pi / 1.55), - omega=2 * np.pi / wavelength * speed_of_light, - ) - H_j = calculate_hfield( - basis, - E_j, - lam_j * (2 * np.pi / 1.55), - omega=2 * np.pi / wavelength * speed_of_light, - ) - overlap_integrals[i, j] = calculate_overlap(basis, E_i, H_i, basis, E_j, H_j) + +num_modes = len(modes_1) + len(modes_2) +overlap_integrals = np.zeros((num_modes, num_modes), dtype=complex) +for i, mode_i in enumerate(chain(modes_1, modes_2)): + for j, mode_j in enumerate(chain(modes_1, modes_2)): + overlap_integrals[i, j] = mode_i.calculate_overlap(mode_j) print("overlap", overlap_integrals) # plt.imshow(np.abs(overlap_integrals)) # plt.colorbar() # plt.show() -coupling_coefficients = np.zeros((len(modes), len(modes)), dtype=complex) -for i, (lam_i, E_i, epsilon_i) in enumerate(modes): - for j, (lam_j, E_j, epsilon_j) in enumerate(modes): +coupling_coefficients = np.zeros((num_modes, num_modes), dtype=complex) +for i, mode_i in enumerate(chain(modes_1, modes_2)): + for j, mode_j in enumerate(chain(modes_1, modes_2)): coupling_coefficients[i, j] = ( k0 * speed_of_light * epsilon_0 - * calculate_coupling_coefficient( - basis0, epsilons[(epsilon_j + 1) % 2] - 1.444**2, basis, E_i, E_j + * mode_i.calculate_coupling_coefficient( + mode_j, epsilons[(j // len(modes_1) + 1) % 2] - 1.444**2 ) * 0.5 ) @@ -225,8 +211,10 @@ ) print(kappas) -delta = 0.5 * (np.real(lams_1[0]) * k0 + kappas[1, 1] - (np.real(lams_2[0]) * k0 + kappas[0, 0])) -print(delta, np.real(lams_1[0]) * k0, kappas[1, 1]) +delta = 0.5 * ( + np.real(modes_1[0].n_eff) * k0 + kappas[1, 1] - (np.real(modes_2[0].n_eff) * k0 + kappas[0, 0]) +) +print(delta, np.real(modes_1[0].n_eff) * k0, kappas[1, 1]) beta_c = (kappas[0, 1] * kappas[1, 0] + delta**2) ** 0.5 @@ -246,8 +234,11 @@ # %% def fun(t, y): phase_matrix = [ - [np.exp(2j * np.pi / wavelength * (lam_i - lam_j) * t) for lam_j, E_j, epsilon_j in modes] - for lam_i, E_i, epsilon_i in modes + [ + np.exp(2j * np.pi / wavelength * (mode_i.n_eff - mode_j.n_eff) * t) + for mode_j in chain(modes_1, modes_2) + ] + for mode_i in chain(modes_1, modes_2) ] matrix = ( np.linalg.inv(overlap_integrals * phase_matrix) @@ -271,28 +262,19 @@ def fun(t, y): # %% R = [] -lam_i = lams_1[0] -E_i = xs_1[0] +lam_i = modes_1[0].n_eff +E_i = modes_1[0].E -for lam_j, E_j in zip(lams_both, xs_both): - H_i = calculate_hfield( - basis, - E_i, - lam_i * (2 * np.pi / 1.55), - omega=2 * np.pi / wavelength * speed_of_light, - ) - H_j = calculate_hfield( - basis, - E_j, - lam_j * (2 * np.pi / 1.55), - omega=2 * np.pi / wavelength * speed_of_light, - ) - R.append(np.abs(calculate_overlap(basis, E_i, H_i, basis, E_j, H_j) ** 2)) +for mode_j in modes_both: + R.append(np.abs(mode_i.calculate_overlap(mode_j) ** 2)) print(R) P = ( R[0] ** 2 + R[1] ** 2 - + 2 * R[0] * R[1] * np.cos(2 * np.pi / wavelength * (lams_both[0] - lams_both[1]) * ts) + + 2 + * R[0] + * R[1] + * np.cos(2 * np.pi / wavelength * (modes_both[0].n_eff - modes_both[1].n_eff) * ts) ) plt.plot(ts, P) From 2603d5e8af300a8d2776728b4087602f128a2764 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 3 Apr 2023 17:25:07 -0700 Subject: [PATCH 10/46] Modes class: return fix and axs in plot --- femwell/mode_solver.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index ee767974..97ef4177 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -103,7 +103,7 @@ def calculate_coupling_coefficient(self, mode, delta_epsilon): ) def plot(self, field, plot_vectors=False, colorbar=True, direction="y", title="E"): - plot_mode( + return plot_mode( self.basis, field, plot_vectors=plot_vectors, From f718bbc4c250ab15e864825cfc05560fbf533c89 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 3 Apr 2023 17:25:23 -0700 Subject: [PATCH 11/46] fiber_overlap example: use object oriented mode solver --- docs/photonics/examples/fiber_overlap.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/photonics/examples/fiber_overlap.py b/docs/photonics/examples/fiber_overlap.py index d32b00ce..973fde6e 100644 --- a/docs/photonics/examples/fiber_overlap.py +++ b/docs/photonics/examples/fiber_overlap.py @@ -61,9 +61,9 @@ # Thus, we know, that we chose the cladding thick enough if the field vanishes at the outer boundaries. # + -lams, basis, xs = compute_modes(basis0, epsilon, wavelength=1.55, mu_r=1, num_modes=1) +modes = compute_modes(basis0, epsilon, wavelength=1.55, mu_r=1, num_modes=1, return_objects=True) -fig, axs = plot_mode(basis, np.real(xs[0]), direction="x") +fig, axs = modes[0].plot(modes[0].E.real, direction="x") plt.tight_layout() plt.show() # - @@ -84,7 +84,7 @@ ) efficiency = overlap( - basis_fiber, basis.interpolate(xs[0])[0][1], basis_fiber.interpolate(x_fiber) + basis_fiber, modes[0].basis.interpolate(modes[0].E)[0][1], basis_fiber.interpolate(x_fiber) ) efficiencies.append(efficiency) From dae9d70bd80cdf7fdec20113d0bf1a14c266723c Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 3 Apr 2023 17:34:31 -0700 Subject: [PATCH 12/46] Mode class: add calculate_propagation_loss --- femwell/mode_solver.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index 97ef4177..b184fa08 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -102,6 +102,9 @@ def calculate_coupling_coefficient(self, mode, delta_epsilon): self.basis_epsilon_r, delta_epsilon, self.basis, self.E, mode.E ) + def calculate_propagation_loss(self, distance): + return -20 / np.log(10) * self.k0 * np.imag(self.n_eff) * distance + def plot(self, field, plot_vectors=False, colorbar=True, direction="y", title="E"): return plot_mode( self.basis, From feb8d90840d7c5ef9c5205e6382f0c4ca87d5e9f Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 3 Apr 2023 17:35:16 -0700 Subject: [PATCH 13/46] adjust leaky waveguide example to object oriented mode solver --- docs/photonics/examples/leaky_waveguide.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/docs/photonics/examples/leaky_waveguide.py b/docs/photonics/examples/leaky_waveguide.py index e6e1120f..64e947ba 100644 --- a/docs/photonics/examples/leaky_waveguide.py +++ b/docs/photonics/examples/leaky_waveguide.py @@ -29,7 +29,7 @@ from skfem.io.meshio import from_meshio from femwell.mesh import mesh_from_OrderedDict -from femwell.mode_solver import compute_modes, plot_mode +from femwell.mode_solver import compute_modes from femwell.visualization import plot_domains # %% [markdown] @@ -105,13 +105,15 @@ # %% wavelength = 1.55 -lams, basis, xs = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=1, order=1) -for i, lam in enumerate(lams): +modes = compute_modes( + basis0, epsilon, wavelength=wavelength, num_modes=1, order=1, return_objects=True +) +for mode in modes: print( - f"Effective refractive index: {lam:.12f}, " - f"Loss: {-20/np.log(10)*2*np.pi/wavelength*np.imag(lam):4f} / dB/um" + f"Effective refractive index: {mode.n_eff:.12f}, " + f"Loss: {mode.calculate_propagation_loss(distance=1):4f} / dB/um" ) - plot_mode(basis, xs[i].real, colorbar=True, direction="x") + mode.plot(mode.E.real, colorbar=True, direction="x") plt.show() # %% From c1b3ee1718d3898de826a8ddac15e2b1ada6e491 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 3 Apr 2023 17:54:33 -0700 Subject: [PATCH 14/46] Mode class: add calcualte_power, Modes calss: add sorted --- femwell/mode_solver.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index b184fa08..6a435311 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -31,7 +31,7 @@ from skfem.utils import solver_eigen_scipy -@dataclass +@dataclass(frozen=True) class Modes: modes: List @@ -45,6 +45,9 @@ def __repr__(self) -> str: modes = "\n\t" + "\n\t".join(repr(mode) for mode in self.modes) + "\n" return f"{self.__class__.__name__}(modes=({modes}))" + def sorted(self, key): + return Modes(modes=sorted(self.modes, key=key)) + @dataclass(frozen=True) class Mode: @@ -105,6 +108,13 @@ def calculate_coupling_coefficient(self, mode, delta_epsilon): def calculate_propagation_loss(self, distance): return -20 / np.log(10) * self.k0 * np.imag(self.n_eff) * distance + def calculate_power(self, elements=None): + if not elements: + basis = self.basis + else: + basis = self.basis.with_elements(elements) + return calculate_overlap(basis, self.E, self.H, basis, self.E, self.H) + def plot(self, field, plot_vectors=False, colorbar=True, direction="y", title="E"): return plot_mode( self.basis, From bfc327aa1527b5254f8200f8bc81d7affb1a4d28 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 3 Apr 2023 17:54:48 -0700 Subject: [PATCH 15/46] selecting_modes object oriented --- docs/photonics/examples/selecting_modes.py | 63 +++++++--------------- 1 file changed, 19 insertions(+), 44 deletions(-) diff --git a/docs/photonics/examples/selecting_modes.py b/docs/photonics/examples/selecting_modes.py index c9a9ef53..9041c107 100644 --- a/docs/photonics/examples/selecting_modes.py +++ b/docs/photonics/examples/selecting_modes.py @@ -128,19 +128,12 @@ # %% # basis0.plot(epsilon, colorbar=True).show() -lams, basis, xs = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=4) +modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=4, return_objects=True) +for mode in modes: + mode.show(mode.E.real, direction="x") -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])}") +print(f"The effective index of the SiN mode is {np.real(modes[2].n_eff)}") # %% [markdown] @@ -169,15 +162,12 @@ 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) - +modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=2, return_objects=True) -plot_mode(basis, np.real(xs[0]), direction="x") -plt.show() -plot_mode(basis, np.real(xs[1]), direction="x") -plt.show() +for mode in modes: + mode.show(mode.E.real, direction="x") -print(f"The effective index of the SiN mode is {np.real(lams[0])}") +print(f"The effective index of the SiN mode is {np.real(modes[0].n_eff)}") # %% [markdown] # ## 2. Giving a guess effective index @@ -194,15 +184,14 @@ 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) - +modes = compute_modes( + basis0, epsilon, wavelength=wavelength, num_modes=2, n_guess=1.62, return_objects=True +) -plot_mode(basis, np.real(xs[0]), direction="x") -plt.show() -plot_mode(basis, np.real(xs[1]), direction="x") -plt.show() +for mode in modes: + mode.show(mode.E.real, direction="x") -print(f"The effective index of the SiN mode is {np.real(lams[1])}") +print(f"The effective index of the SiN mode is {np.real(modes[1].n_eff)}") # %% [markdown] @@ -227,26 +216,13 @@ # %% -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, - ) - ) +modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=4, return_objects=True) # Option 1: using an element name -ind_mode = argsort_modes_by_power_in_elements(basis, xs, H_modes, elements="sin")[0] +modes_sorted = modes.sorted(key=lambda mode: mode.calculate_power(elements="sin")) -plot_mode(basis, np.real(xs[ind_mode]), direction="x") -plt.show() +modes_sorted[0].show(modes_sorted[0].E.real, direction="x") # Option 2: using bounding box @@ -254,9 +230,8 @@ bbox = [-2, 0, 0, 0.4] elements = inside_bbox(bbox) -ind_mode = argsort_modes_by_power_in_elements(basis, xs, H_modes, elements)[0] +modes_sorted = modes.sorted(key=lambda mode: mode.calculate_power(elements=elements)) -plot_mode(basis, np.real(xs[ind_mode]), direction="x") -plt.show() +modes_sorted[0].show(modes_sorted[0].E.real, direction="x") # %% From b48d9554e6a3a8c6b78c34d4659cbdf9af29d108 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 3 Apr 2023 18:00:44 -0700 Subject: [PATCH 16/46] remove unneeded imports --- femwell/mode_solver.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index 6a435311..6ab38079 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -1,7 +1,6 @@ """Waveguide analysis based on https://doi.org/10.1080/02726340290084012.""" from dataclasses import dataclass -from functools import cache, cached_property -from time import time +from functools import cached_property from typing import List import matplotlib.pyplot as plt @@ -9,7 +8,6 @@ import scipy.constants import scipy.sparse.linalg from numpy.typing import NDArray -from pyparsing import col from scipy.constants import epsilon_0, speed_of_light from skfem import ( Basis, From 3b517614d2fc1abc85a992502814fadab145bdc2 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 3 Apr 2023 18:04:50 -0700 Subject: [PATCH 17/46] simplify coupled mode theory example --- docs/photonics/examples/coupled_mode_theory.py | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/docs/photonics/examples/coupled_mode_theory.py b/docs/photonics/examples/coupled_mode_theory.py index 589b1b5a..c25f0e98 100644 --- a/docs/photonics/examples/coupled_mode_theory.py +++ b/docs/photonics/examples/coupled_mode_theory.py @@ -260,13 +260,7 @@ def fun(t, y): # ## two modes # %% -R = [] - -lam_i = modes_1[0].n_eff -E_i = modes_1[0].E - -for mode_j in modes_both: - R.append(np.abs(mode_i.calculate_overlap(mode_j) ** 2)) +R = [np.abs(modes_1[0].calculate_overlap(mode_j) ** 2) for mode_j in modes_both] print(R) P = ( R[0] ** 2 From 1e6a5d6a8c2a3cf093a54fc7e94e2f4ea92cb6af Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 3 Apr 2023 18:07:30 -0700 Subject: [PATCH 18/46] simplify mode_solver --- femwell/mode_solver.py | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index 6ab38079..7d8473f4 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -557,12 +557,6 @@ def argsort_modes_by_power_in_elements(modes, elements): # 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) - print( - argsort_modes_by_power_in_elements( - modes=modes, - elements=sel_fun, - ) - ) + print(modes.sorted(lambda mode: mode.calculate_power(elements=sel_fun))) From 3bef0acf68dd51cd19d89b1753a2e21628253f7f Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 3 Apr 2023 18:12:47 -0700 Subject: [PATCH 19/46] add deprecation warning --- femwell/mode_solver.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index 7d8473f4..3d0e24b3 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -241,6 +241,13 @@ def bform(e_t, e_z, v_t, v_z, w): ] ) else: + import warnings + + warnings.warn( + "return_objects will become the default behaviour with version 0.1.0." + "The current behaviour with return_objects=False is deprecated and will be removed.", + DeprecationWarning, + ) return np.sqrt(lams)[:num_modes] / k0, basis, xs[:num_modes] From 10f51df7275f8482cce9bf2dfe2e14802bf6b15e Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 3 Apr 2023 18:36:32 -0700 Subject: [PATCH 20/46] adjust vary_width and vary_wavelength --- docs/photonics/examples/vary_wavelength.py | 25 +++++++++++++--------- docs/photonics/examples/vary_width.py | 21 +++++++++++------- femwell/mode_solver.py | 2 +- 3 files changed, 29 insertions(+), 19 deletions(-) diff --git a/docs/photonics/examples/vary_wavelength.py b/docs/photonics/examples/vary_wavelength.py index 57e3d065..115aeed3 100644 --- a/docs/photonics/examples/vary_wavelength.py +++ b/docs/photonics/examples/vary_wavelength.py @@ -79,7 +79,7 @@ def n_silicon_dioxide(wavelength): wavelengths = np.linspace(1.2, 1.9, 20) num_modes = 2 -all_lams = np.zeros((wavelengths.shape[0], num_modes)) +all_neffs = np.zeros((wavelengths.shape[0], num_modes)) all_te_fracs = np.zeros((wavelengths.shape[0], num_modes)) for i, wavelength in enumerate(tqdm(wavelengths)): basis0 = Basis(mesh, ElementTriP0()) @@ -91,11 +91,16 @@ def n_silicon_dioxide(wavelength): }.items(): epsilon[basis0.get_dofs(elements=subdomain)] = n**2 - lams, basis, xs = compute_modes( - basis0, epsilon, wavelength=wavelength, num_modes=num_modes, normalize=False + modes = compute_modes( + basis0, + epsilon, + wavelength=wavelength, + num_modes=num_modes, + normalize=False, + return_objects=True, ) - all_lams[i] = np.real(lams) - all_te_fracs[i, :] = [calculate_te_frac(basis, xs[idx]) for idx in range(num_modes)] + all_neffs[i] = np.real([mode.n_eff for mode in modes]) + all_te_fracs[i, :] = [mode.te_fraction for mode in modes] # - # Now we only look at the real part of the effective refractive indices and plot all the dispersion relation, @@ -103,21 +108,21 @@ def n_silicon_dioxide(wavelength): # So convenient! # + tags=["hide-input"] -all_lams = np.real(all_lams) +all_neffs = np.real(all_neffs) fig, axs = plt.subplots(1, 3) axs[0].set_xlabel("Wavelength / µm") axs[0].set_ylabel("Effective refractive index") -axs[0].set_ylim(1.444, np.max(all_lams) + 0.1 * (np.max(all_lams) - 1.444)) -for lams, te_fracs in zip(all_lams.T, all_te_fracs.T): +axs[0].set_ylim(1.444, np.max(all_neffs) + 0.1 * (np.max(all_neffs) - 1.444)) +for lams, te_fracs in zip(all_neffs.T, all_te_fracs.T): axs[0].plot(wavelengths, lams) sc = axs[0].scatter(wavelengths, lams, c=te_fracs, cmap="cool", vmin=0, vmax=1) axs[1].set_xlabel("Wavelength / µm") axs[1].set_ylabel("Group velocity $v_g$") -for lams, te_fracs in zip(all_lams.T, all_te_fracs.T): +for lams, te_fracs in zip(all_neffs.T, all_te_fracs.T): fit = Polynomial.fit(wavelengths, lams, deg=4) y = speed_of_light / (lams - wavelengths * fit.deriv(1)(wavelengths)) axs[1].plot(wavelengths, y) @@ -125,7 +130,7 @@ def n_silicon_dioxide(wavelength): axs[2].set_xlabel("Wavelength / µm") axs[2].set_ylabel("Group velocity dispersion coefficient $D$") -for lams, te_fracs in zip(all_lams.T, all_te_fracs.T): +for lams, te_fracs in zip(all_neffs.T, all_te_fracs.T): fit = Polynomial.fit(wavelengths, lams, deg=4) y = wavelengths**2 * fit.deriv(2)(wavelengths) axs[2].plot(wavelengths, y) diff --git a/docs/photonics/examples/vary_width.py b/docs/photonics/examples/vary_width.py index 0c498add..3ea6ab4f 100644 --- a/docs/photonics/examples/vary_width.py +++ b/docs/photonics/examples/vary_width.py @@ -37,7 +37,7 @@ num_modes = 8 widths = np.linspace(0.5, 3.5, 100) -all_lams = np.zeros((widths.shape[0], num_modes)) +all_neffs = np.zeros((widths.shape[0], num_modes)) all_te_fracs = np.zeros((widths.shape[0], num_modes)) for i, width in enumerate(tqdm(widths)): core = box(0, 0, width, 0.5) @@ -56,19 +56,24 @@ for subdomain, n in {"core": 1.9963, "box": 1.444, "clad": 1}.items(): epsilon[basis0.get_dofs(elements=subdomain)] = n**2 - lams, basis, xs = compute_modes( - basis0, epsilon, wavelength=wavelength, num_modes=num_modes, normalize=False + modes = compute_modes( + basis0, + epsilon, + wavelength=wavelength, + num_modes=num_modes, + normalize=False, + return_objects=True, ) - all_lams[i] = np.real(lams) - all_te_fracs[i, :] = [calculate_te_frac(basis, xs[idx]) for idx in range(num_modes)] + all_neffs[i] = np.real([mode.n_eff for mode in modes]) + all_te_fracs[i, :] = [mode.te_fraction for mode in modes] # - # + tags=["hide-input"] -all_lams = np.real(all_lams) +all_neffs = np.real(all_neffs) plt.xlabel("Width of waveguide / µm") plt.ylabel("Effective refractive index") -plt.ylim(1.444, np.max(all_lams) + 0.1 * (np.max(all_lams) - 1.444)) -for lams, te_fracs in zip(all_lams.T, all_te_fracs.T): +plt.ylim(1.444, np.max(all_neffs) + 0.1 * (np.max(all_neffs) - 1.444)) +for lams, te_fracs in zip(all_neffs.T, all_te_fracs.T): plt.plot(widths, lams) plt.scatter(widths, lams, c=te_fracs, cmap="cool") plt.colorbar().set_label("TE fraction") diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index 3d0e24b3..70f8abe0 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -215,7 +215,7 @@ def bform(e_t, e_z, v_t, v_z, w): ) # undo the scaling E_3,new = beta * E_3 hs = [] - if normalize: + if normalize or return_objects: for i, lam in enumerate(lams): H = calculate_hfield( basis, xs[i], np.sqrt(lam), omega=k0 * scipy.constants.speed_of_light From 16d7be0b9820c83ab6a5bb214b0383269f192811 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 3 Apr 2023 18:39:35 -0700 Subject: [PATCH 21/46] remove unneeded imports --- docs/photonics/examples/vary_wavelength.py | 2 +- docs/photonics/examples/vary_width.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/photonics/examples/vary_wavelength.py b/docs/photonics/examples/vary_wavelength.py index 115aeed3..1448a507 100644 --- a/docs/photonics/examples/vary_wavelength.py +++ b/docs/photonics/examples/vary_wavelength.py @@ -29,7 +29,7 @@ from tqdm import tqdm from femwell.mesh import mesh_from_OrderedDict -from femwell.mode_solver import calculate_te_frac, compute_modes +from femwell.mode_solver import compute_modes # - diff --git a/docs/photonics/examples/vary_width.py b/docs/photonics/examples/vary_width.py index 3ea6ab4f..1c5259f4 100644 --- a/docs/photonics/examples/vary_width.py +++ b/docs/photonics/examples/vary_width.py @@ -28,7 +28,7 @@ from tqdm import tqdm from femwell.mesh import mesh_from_OrderedDict -from femwell.mode_solver import calculate_te_frac, compute_modes +from femwell.mode_solver import compute_modes # - From 2a2c44b0bbab17c5d211407a048d577b30e8087e Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 3 Apr 2023 19:02:15 -0700 Subject: [PATCH 22/46] add comments --- femwell/mode_solver.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index 70f8abe0..bb6e1ab8 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -50,31 +50,44 @@ def sorted(self, key): @dataclass(frozen=True) class Mode: frequency: float + """Frequency of the light""" k: float + """Propagation constant of the mode""" basis_epsilon_r: Basis + """Basis used for epsilon_r""" epsilon_r: NDArray + """Epsilon_r with which the mode was calculated""" basis: Basis + """Basis on which the mode was calculated and E/H are defined""" E: NDArray + """Electric field of the mode""" H: NDArray + """Magnetic field of the mode""" @property def omega(self): + """Angular frequency of the light""" return 2 * np.pi * self.frequency @property def k0(self): + """Vacuum propagation constant of the light""" return self.omega / speed_of_light @property def wavelength(self): + """Vacuum wavelength of the light""" return speed_of_light / self.frequency @property def n_eff(self): + """Effective refractive index of the mode""" return self.k / self.k0 @cached_property def te_fraction(self): + """TE-fraction of the mode""" + @Functional def ex(w): return np.abs(w.E[0][0]) ** 2 @@ -90,6 +103,8 @@ def ey(w): @cached_property def tm_fraction(self): + """TM-fraction of the mode""" + return 1 - self.te_fraction def __repr__(self) -> str: From d7e656a999e10968f7e9d937a0daf881ca80c745 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Fri, 7 Apr 2023 19:43:30 -0700 Subject: [PATCH 23/46] move mode solver to new file, now scale invariant --- docs/photonics/examples/bent_waveguide.py | 5 +- .../photonics/examples/coupled_mode_theory.py | 14 +- docs/photonics/examples/fiber_overlap.py | 4 +- docs/photonics/examples/leaky_waveguide.py | 6 +- docs/photonics/examples/selecting_modes.py | 19 +- docs/photonics/examples/vary_wavelength.py | 7 +- docs/photonics/examples/vary_width.py | 9 +- femwell/maxwell/waveguide.py | 350 ++++++++++++++++++ femwell/mode_solver.py | 244 +----------- 9 files changed, 375 insertions(+), 283 deletions(-) create mode 100644 femwell/maxwell/waveguide.py diff --git a/docs/photonics/examples/bent_waveguide.py b/docs/photonics/examples/bent_waveguide.py index 23f46c0a..7ac6765f 100644 --- a/docs/photonics/examples/bent_waveguide.py +++ b/docs/photonics/examples/bent_waveguide.py @@ -25,8 +25,8 @@ from skfem.io.meshio import from_meshio from tqdm import tqdm +from femwell.maxwell.waveguide import compute_modes from femwell.mesh import mesh_from_OrderedDict -from femwell.mode_solver import compute_modes # - @@ -88,7 +88,7 @@ # + modes_straight = compute_modes( - basis0, epsilon, wavelength=wavelength, num_modes=1, order=2, radius=np.inf, return_objects=True + basis0, epsilon, wavelength=wavelength, num_modes=1, order=2, radius=np.inf ) # - @@ -111,7 +111,6 @@ radius=radius, n_guess=lam_guess, solver="scipy", - return_objects=True, ) lam_guess = modes[0].n_eff radiuss_lams.append(modes[0].n_eff) diff --git a/docs/photonics/examples/coupled_mode_theory.py b/docs/photonics/examples/coupled_mode_theory.py index c25f0e98..4e57938c 100644 --- a/docs/photonics/examples/coupled_mode_theory.py +++ b/docs/photonics/examples/coupled_mode_theory.py @@ -38,8 +38,8 @@ from skfem import Basis, ElementTriP0, Mesh from skfem.io import from_meshio +from femwell.maxwell.waveguide import compute_modes from femwell.mesh import mesh_from_OrderedDict -from femwell.mode_solver import compute_modes # %% [markdown] # Let's set up the geometry! @@ -114,9 +114,7 @@ epsilon = basis0.zeros() + 1.444**2 epsilon[basis0.get_dofs(elements=("core_1", "core_2"))] = 3.4777**2 # basis0.plot(epsilon, colorbar=True).show() -modes_both = compute_modes( - basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=2, return_objects=True -) +modes_both = compute_modes(basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=2) modes_both[0].show(modes_both[0].E.real, direction="x") modes_both[1].show(modes_both[1].E.real, direction="x") print( @@ -137,18 +135,14 @@ epsilon = basis0.zeros() + 1.444**2 epsilon[basis0.get_dofs(elements="core_1")] = 3.4777**2 # basis0.plot(epsilon, colorbar=True).show() -modes_1 = compute_modes( - basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=1, return_objects=True -) +modes_1 = compute_modes(basis0, epsilon, wavelength=wavelength, mu_r=1, num_modes=1) print("Effective refractive index of the mode of the first waveguide", modes_1[0].n_eff) modes_1[0].show(modes_1[0].E.real, direction="x") epsilon_2 = basis0.zeros() + 1.444**2 epsilon_2[basis0.get_dofs(elements="core_2")] = 3.4777**2 # basis0.plot(epsilon_2, colorbar=True).show() -modes_2 = compute_modes( - basis0, epsilon_2, wavelength=wavelength, mu_r=1, num_modes=1, return_objects=True -) +modes_2 = compute_modes(basis0, epsilon_2, wavelength=wavelength, mu_r=1, num_modes=1) print("Effective refractive index of the mode of the second waveguide", modes_2[0].n_eff) modes_2[0].show(modes_2[0].E.real, direction="x") diff --git a/docs/photonics/examples/fiber_overlap.py b/docs/photonics/examples/fiber_overlap.py index 973fde6e..b5df9dbf 100644 --- a/docs/photonics/examples/fiber_overlap.py +++ b/docs/photonics/examples/fiber_overlap.py @@ -25,8 +25,8 @@ from tqdm import tqdm from femwell.fiber import e_field_gaussian, overlap +from femwell.maxwell.waveguide import compute_modes from femwell.mesh import mesh_from_OrderedDict -from femwell.mode_solver import compute_modes, plot_mode # - @@ -61,7 +61,7 @@ # Thus, we know, that we chose the cladding thick enough if the field vanishes at the outer boundaries. # + -modes = compute_modes(basis0, epsilon, wavelength=1.55, mu_r=1, num_modes=1, return_objects=True) +modes = compute_modes(basis0, epsilon, wavelength=1.55, mu_r=1, num_modes=1) fig, axs = modes[0].plot(modes[0].E.real, direction="x") plt.tight_layout() diff --git a/docs/photonics/examples/leaky_waveguide.py b/docs/photonics/examples/leaky_waveguide.py index 64e947ba..b60ac1ed 100644 --- a/docs/photonics/examples/leaky_waveguide.py +++ b/docs/photonics/examples/leaky_waveguide.py @@ -28,8 +28,8 @@ from skfem import Basis, ElementDG, ElementTriP1 from skfem.io.meshio import from_meshio +from femwell.maxwell.waveguide import compute_modes from femwell.mesh import mesh_from_OrderedDict -from femwell.mode_solver import compute_modes from femwell.visualization import plot_domains # %% [markdown] @@ -105,9 +105,7 @@ # %% wavelength = 1.55 -modes = compute_modes( - basis0, epsilon, wavelength=wavelength, num_modes=1, order=1, return_objects=True -) +modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=1, order=1) for mode in modes: print( f"Effective refractive index: {mode.n_eff:.12f}, " diff --git a/docs/photonics/examples/selecting_modes.py b/docs/photonics/examples/selecting_modes.py index 9041c107..5364434a 100644 --- a/docs/photonics/examples/selecting_modes.py +++ b/docs/photonics/examples/selecting_modes.py @@ -33,15 +33,8 @@ from skfem import Basis, ElementTriP0, Mesh from skfem.io import from_meshio +from femwell.maxwell.waveguide import compute_modes 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] @@ -128,7 +121,7 @@ # %% # basis0.plot(epsilon, colorbar=True).show() -modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=4, return_objects=True) +modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=4) for mode in modes: mode.show(mode.E.real, direction="x") @@ -162,7 +155,7 @@ epsilon[basis0.get_dofs(elements=("si"))] = 1.444**2 epsilon[basis0.get_dofs(elements=("sin"))] = 1.973**2 -modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=2, return_objects=True) +modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=2) for mode in modes: mode.show(mode.E.real, direction="x") @@ -184,9 +177,7 @@ epsilon[basis0.get_dofs(elements=("si"))] = 3.4777**2 epsilon[basis0.get_dofs(elements=("sin"))] = 1.973**2 -modes = compute_modes( - basis0, epsilon, wavelength=wavelength, num_modes=2, n_guess=1.62, return_objects=True -) +modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=2, n_guess=1.62) for mode in modes: mode.show(mode.E.real, direction="x") @@ -216,7 +207,7 @@ # %% -modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=4, return_objects=True) +modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=4) # Option 1: using an element name diff --git a/docs/photonics/examples/vary_wavelength.py b/docs/photonics/examples/vary_wavelength.py index 1448a507..1490bf8e 100644 --- a/docs/photonics/examples/vary_wavelength.py +++ b/docs/photonics/examples/vary_wavelength.py @@ -92,12 +92,7 @@ def n_silicon_dioxide(wavelength): epsilon[basis0.get_dofs(elements=subdomain)] = n**2 modes = compute_modes( - basis0, - epsilon, - wavelength=wavelength, - num_modes=num_modes, - normalize=False, - return_objects=True, + basis0, epsilon, wavelength=wavelength, num_modes=num_modes, normalize=False ) all_neffs[i] = np.real([mode.n_eff for mode in modes]) all_te_fracs[i, :] = [mode.te_fraction for mode in modes] diff --git a/docs/photonics/examples/vary_width.py b/docs/photonics/examples/vary_width.py index 1c5259f4..516e3a0d 100644 --- a/docs/photonics/examples/vary_width.py +++ b/docs/photonics/examples/vary_width.py @@ -27,8 +27,8 @@ from skfem.io.meshio import from_meshio from tqdm import tqdm +from femwell.maxwell.waveguide import compute_modes from femwell.mesh import mesh_from_OrderedDict -from femwell.mode_solver import compute_modes # - @@ -57,12 +57,7 @@ epsilon[basis0.get_dofs(elements=subdomain)] = n**2 modes = compute_modes( - basis0, - epsilon, - wavelength=wavelength, - num_modes=num_modes, - normalize=False, - return_objects=True, + basis0, epsilon, wavelength=wavelength, num_modes=num_modes, normalize=False ) all_neffs[i] = np.real([mode.n_eff for mode in modes]) all_te_fracs[i, :] = [mode.te_fraction for mode in modes] diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py new file mode 100644 index 00000000..1db4ad58 --- /dev/null +++ b/femwell/maxwell/waveguide.py @@ -0,0 +1,350 @@ +"""Waveguide analysis based on https://doi.org/10.1080/02726340290084012.""" +from dataclasses import dataclass +from functools import cached_property +from typing import List + +import matplotlib.pyplot as plt +import numpy as np +import scipy.constants +import scipy.sparse.linalg +from numpy.typing import NDArray +from scipy.constants import epsilon_0, speed_of_light +from skfem import ( + Basis, + BilinearForm, + ElementDG, + ElementTriN1, + ElementTriN2, + ElementTriP0, + ElementTriP1, + ElementTriP2, + ElementVector, + Functional, + LinearForm, + Mesh, + condense, + solve, +) +from skfem.helpers import cross, curl, dot, grad, inner +from skfem.utils import solver_eigen_scipy + +from femwell.mode_solver import ( + calculate_coupling_coefficient, + calculate_energy_current_density, + calculate_hfield, + calculate_overlap, + calculate_scalar_product, + plot_mode, +) + + +@dataclass(frozen=True) +class Modes: + modes: List + + def __getitem__(self, idx): + return self.modes[idx] + + def __len__(self): + return len(self.modes) + + def __repr__(self) -> str: + modes = "\n\t" + "\n\t".join(repr(mode) for mode in self.modes) + "\n" + return f"{self.__class__.__name__}(modes=({modes}))" + + def sorted(self, key): + return Modes(modes=sorted(self.modes, key=key)) + + +@dataclass(frozen=True) +class Mode: + frequency: float + """Frequency of the light""" + k: float + """Propagation constant of the mode""" + basis_epsilon_r: Basis + """Basis used for epsilon_r""" + epsilon_r: NDArray + """Epsilon_r with which the mode was calculated""" + basis: Basis + """Basis on which the mode was calculated and E/H are defined""" + E: NDArray + """Electric field of the mode""" + H: NDArray + """Magnetic field of the mode""" + + @property + def omega(self): + """Angular frequency of the light""" + return 2 * np.pi * self.frequency + + @property + def k0(self): + """Vacuum propagation constant of the light""" + return self.omega / speed_of_light + + @property + def wavelength(self): + """Vacuum wavelength of the light""" + return speed_of_light / self.frequency + + @property + def n_eff(self): + """Effective refractive index of the mode""" + return self.k / self.k0 + + @cached_property + def te_fraction(self): + """TE-fraction of the mode""" + + @Functional + def ex(w): + return np.abs(w.E[0][0]) ** 2 + + @Functional + def ey(w): + return np.abs(w.E[0][1]) ** 2 + + ex_sum = ex.assemble(self.basis, E=self.basis.interpolate(self.E)) + ey_sum = ey.assemble(self.basis, E=self.basis.interpolate(self.E)) + + return ex_sum / (ex_sum + ey_sum) + + @cached_property + def tm_fraction(self): + """TM-fraction of the mode""" + + return 1 - self.te_fraction + + def __repr__(self) -> str: + return f"{self.__class__.__name__}(k: {self.k}, n_eff:{self.n_eff})" + + def calculate_overlap(self, mode): + return calculate_overlap(self.basis, self.E, self.H, mode.basis, mode.E, mode.H) + + def calculate_coupling_coefficient(self, mode, delta_epsilon): + return calculate_coupling_coefficient( + self.basis_epsilon_r, delta_epsilon, self.basis, self.E, mode.E + ) + + def calculate_propagation_loss(self, distance): + return -20 / np.log(10) * self.k0 * np.imag(self.n_eff) * distance + + def calculate_power(self, elements=None): + if not elements: + basis = self.basis + else: + basis = self.basis.with_elements(elements) + return calculate_overlap(basis, self.E, self.H, basis, self.E, self.H) + + def plot(self, field, plot_vectors=False, colorbar=True, direction="y", title="E"): + return plot_mode( + self.basis, + field, + plot_vectors=plot_vectors, + colorbar=colorbar, + title=title, + direction=direction, + ) + + def show(self, field, **kwargs): + self.plot(field=field, **kwargs) + plt.show() + + +def compute_modes( + basis_epsilon_r, + epsilon_r, + wavelength, + mu_r=1, + num_modes=1, + order=1, + metallic_boundaries=False, + radius=np.inf, + n_guess=None, + solver="slepc", +): + if solver == "scipy": + solver = solver_eigen_scipy + elif solver == "slepc": + from femwell.solver import solver_eigen_slepc + + solver = solver_eigen_slepc + else: + raise ValueError("`solver` must either be `scipy` or `slepc`") + + k0 = 2 * np.pi / wavelength + + if order == 1: + element = ElementTriN1() * ElementTriP1() + elif order == 2: + element = ElementTriN2() * ElementTriP2() + else: + raise AssertionError("Only order 1 and 2 implemented by now.") + + basis = basis_epsilon_r.with_element(element) + basis_epsilon_r = basis.with_element(basis_epsilon_r.elem) # adjust quadrature + + print(k0) + + @BilinearForm(dtype=epsilon_r.dtype) + def aform(e_t, e_z, v_t, v_z, w): + epsilon = w.epsilon * (1 + w.x[0] / radius) ** 2 + + return ( + 1 / mu_r * curl(e_t) * curl(v_t) / k0**2 + - epsilon * dot(e_t, v_t) + + 1 / mu_r * dot(grad(e_z), v_t) + + epsilon * inner(e_t, grad(v_z)) + - epsilon * e_z * v_z * k0**2 + ) + + @BilinearForm(dtype=epsilon_r.dtype) + def bform(e_t, e_z, v_t, v_z, w): + return -1 / mu_r * dot(e_t, v_t) / k0**2 + + A = aform.assemble(basis, epsilon=basis_epsilon_r.interpolate(epsilon_r)) + B = bform.assemble(basis, epsilon=basis_epsilon_r.interpolate(epsilon_r)) + + if n_guess: + sigma = sigma = k0**2 * n_guess**2 + else: + sigma = sigma = k0**2 * np.max(epsilon_r) ** 2 + + if metallic_boundaries: + lams, xs = solve( + *condense(-A, -B, D=basis.get_dofs(), x=basis.zeros(dtype=complex)), + solver=solver(k=num_modes, sigma=sigma), + ) + else: + lams, xs = solve( + -A, + -B, + solver=solver(k=num_modes, sigma=sigma), + ) + + idx = np.abs(np.real(lams)).argsort()[::-1] + lams = lams[idx] + xs = xs[:, idx] + + xs = xs.T + print([np.sum(np.abs(xs[0, basis.split_indices()[i]])) for i in range(2)]) + xs[:, basis.split_indices()[1]] /= 1j * np.sqrt( + lams[:, np.newaxis] / k0**4 + ) # undo the scaling E_3,new = beta * E_3 + + hs = [] + for i, lam in enumerate(lams): + H = calculate_hfield(basis, xs[i], np.sqrt(lam), omega=k0 * scipy.constants.speed_of_light) + power = calculate_overlap(basis, xs[i], H, basis, xs[i], H) + xs[i] /= np.sqrt(power) + H /= np.sqrt(power) + hs.append(H) + + return Modes( + modes=[ + Mode( + frequency=speed_of_light / wavelength, + k=np.sqrt(lams[i]), + basis_epsilon_r=basis_epsilon_r, + epsilon_r=epsilon_r, + basis=basis, + E=xs[i], + H=hs[i], + ) + for i in range(num_modes) + ] + ) + + +if __name__ == "__main__": + from collections import OrderedDict + + from shapely.geometry import Polygon + + from femwell.mesh import mesh_from_OrderedDict + + x_min = 0 + w_sim = 3 + h_clad = 0.7 + h_box = 0.5 + w_core = 1 + h_core = 0.22 + offset_heater = 2.2 + h_heater = 0.14 + w_heater = 2 + + polygons = OrderedDict( + core=Polygon( + [ + (x_min - w_core / 2, 0), + (x_min - w_core / 2, h_core), + (x_min + w_core / 2, h_core), + (x_min + w_core / 2, 0), + ] + ), + clad=Polygon( + [ + (x_min - w_sim / 2, 0), + (x_min - w_sim / 2, h_clad), + (x_min + w_sim / 2, h_clad), + (x_min + w_sim / 2, 0), + ] + ), + box=Polygon( + [ + (x_min - w_sim / 2, 0), + (x_min - w_sim / 2, -h_box), + (x_min + w_sim / 2, -h_box), + (x_min + w_sim / 2, 0), + ] + ), + ) + + resolutions = dict(core={"resolution": 0.05, "distance": 1}) + + mesh_from_OrderedDict(polygons, resolutions, filename="mesh.msh", default_resolution_max=0.2) + + scale = 1e0 + mesh = Mesh.load("mesh.msh").scaled(scale) + basis = Basis(mesh, ElementTriN2() * ElementTriP2()) + basis0 = basis.with_element(ElementTriP0()) + epsilon = basis0.zeros(dtype=complex) + epsilon[basis0.get_dofs(elements="core")] = 3.4777**2 + epsilon[basis0.get_dofs(elements="clad")] = 1.444**2 + epsilon[basis0.get_dofs(elements="box")] = 1.444**2 + # basis0.plot(epsilon, colorbar=True).show() + + modes = compute_modes( + basis0, + epsilon, + wavelength=1.55 * scale, + mu_r=1, + num_modes=6, + order=2, + # radius=3 * scale, + ) + print(modes) + print(modes[0].te_fraction) + + modes[0].show(np.real(modes[0].E)) + modes[0].show(np.imag(modes[0].E)) + + modes[0].show(np.real(modes[0].H)) + modes[0].show(np.imag(modes[0].H)) + + integrals = np.zeros((len(modes),) * 2, dtype=complex) + + for i in range(len(modes)): + for j in range(len(modes)): + integrals[i, j] = modes[i].calculate_overlap(modes[j]) + + plt.imshow(np.real(integrals)) + plt.colorbar() + plt.show() + + # Create basis to select a certain simulation extent + def sel_fun(x): + return (x[0] < 0) * (x[0] > -1) * (x[1] > 0) * (x[1] < 0.5) + + print(modes.sorted(lambda mode: mode.calculate_power(elements=sel_fun))) diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index bb6e1ab8..4215595c 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -1,13 +1,9 @@ """Waveguide analysis based on https://doi.org/10.1080/02726340290084012.""" -from dataclasses import dataclass -from functools import cached_property -from typing import List import matplotlib.pyplot as plt import numpy as np import scipy.constants import scipy.sparse.linalg -from numpy.typing import NDArray from scipy.constants import epsilon_0, speed_of_light from skfem import ( Basis, @@ -21,7 +17,6 @@ ElementVector, Functional, LinearForm, - Mesh, condense, solve, ) @@ -29,120 +24,6 @@ from skfem.utils import solver_eigen_scipy -@dataclass(frozen=True) -class Modes: - modes: List - - def __getitem__(self, idx): - return self.modes[idx] - - def __len__(self): - return len(self.modes) - - def __repr__(self) -> str: - modes = "\n\t" + "\n\t".join(repr(mode) for mode in self.modes) + "\n" - return f"{self.__class__.__name__}(modes=({modes}))" - - def sorted(self, key): - return Modes(modes=sorted(self.modes, key=key)) - - -@dataclass(frozen=True) -class Mode: - frequency: float - """Frequency of the light""" - k: float - """Propagation constant of the mode""" - basis_epsilon_r: Basis - """Basis used for epsilon_r""" - epsilon_r: NDArray - """Epsilon_r with which the mode was calculated""" - basis: Basis - """Basis on which the mode was calculated and E/H are defined""" - E: NDArray - """Electric field of the mode""" - H: NDArray - """Magnetic field of the mode""" - - @property - def omega(self): - """Angular frequency of the light""" - return 2 * np.pi * self.frequency - - @property - def k0(self): - """Vacuum propagation constant of the light""" - return self.omega / speed_of_light - - @property - def wavelength(self): - """Vacuum wavelength of the light""" - return speed_of_light / self.frequency - - @property - def n_eff(self): - """Effective refractive index of the mode""" - return self.k / self.k0 - - @cached_property - def te_fraction(self): - """TE-fraction of the mode""" - - @Functional - def ex(w): - return np.abs(w.E[0][0]) ** 2 - - @Functional - def ey(w): - return np.abs(w.E[0][1]) ** 2 - - ex_sum = ex.assemble(self.basis, E=self.basis.interpolate(self.E)) - ey_sum = ey.assemble(self.basis, E=self.basis.interpolate(self.E)) - - return ex_sum / (ex_sum + ey_sum) - - @cached_property - def tm_fraction(self): - """TM-fraction of the mode""" - - return 1 - self.te_fraction - - def __repr__(self) -> str: - return f"{self.__class__.__name__}(k: {self.k}, n_eff:{self.n_eff})" - - def calculate_overlap(self, mode): - return calculate_overlap(self.basis, self.E, self.H, mode.basis, mode.E, mode.H) - - def calculate_coupling_coefficient(self, mode, delta_epsilon): - return calculate_coupling_coefficient( - self.basis_epsilon_r, delta_epsilon, self.basis, self.E, mode.E - ) - - def calculate_propagation_loss(self, distance): - return -20 / np.log(10) * self.k0 * np.imag(self.n_eff) * distance - - def calculate_power(self, elements=None): - if not elements: - basis = self.basis - else: - basis = self.basis.with_elements(elements) - return calculate_overlap(basis, self.E, self.H, basis, self.E, self.H) - - def plot(self, field, plot_vectors=False, colorbar=True, direction="y", title="E"): - return plot_mode( - self.basis, - field, - plot_vectors=plot_vectors, - colorbar=colorbar, - title=title, - direction=direction, - ) - - def show(self, field, **kwargs): - self.plot(field=field, **kwargs) - plt.show() - - def compute_modes( basis_epsilon_r, epsilon_r, @@ -156,7 +37,6 @@ def compute_modes( solver="slepc", normalize=True, cache_path=None, - return_objects=False, ): if solver == "scipy": solver = solver_eigen_scipy @@ -230,7 +110,7 @@ def bform(e_t, e_z, v_t, v_z, w): ) # undo the scaling E_3,new = beta * E_3 hs = [] - if normalize or return_objects: + if normalize: for i, lam in enumerate(lams): H = calculate_hfield( basis, xs[i], np.sqrt(lam), omega=k0 * scipy.constants.speed_of_light @@ -240,30 +120,13 @@ def bform(e_t, e_z, v_t, v_z, w): H /= np.sqrt(power) hs.append(H) - if return_objects: - return Modes( - modes=[ - Mode( - frequency=speed_of_light / wavelength, - k=np.sqrt(lams[i]), - basis_epsilon_r=basis_epsilon_r, - epsilon_r=epsilon_r, - basis=basis, - E=xs[i], - H=hs[i], - ) - for i in range(num_modes) - ] - ) - else: - import warnings + import warnings - warnings.warn( - "return_objects will become the default behaviour with version 0.1.0." - "The current behaviour with return_objects=False is deprecated and will be removed.", - DeprecationWarning, - ) - return np.sqrt(lams)[:num_modes] / k0, basis, xs[:num_modes] + warnings.warn( + "femwell.maxwell.waveguide will replace this module with version 0.1.0.", + DeprecationWarning, + ) + return np.sqrt(lams)[:num_modes] / k0, basis, xs[:num_modes] def calculate_hfield(basis, xs, beta, omega=1): @@ -489,96 +352,3 @@ def argsort_modes_by_power_in_elements(modes, elements): ] return np.argsort(np.abs(overlaps))[::-1] - - -if __name__ == "__main__": - from collections import OrderedDict - - from shapely.geometry import Polygon - - from femwell.mesh import mesh_from_OrderedDict - - x_min = 0 - w_sim = 3 - h_clad = 0.7 - h_box = 0.5 - w_core = 1 - h_core = 0.22 - offset_heater = 2.2 - h_heater = 0.14 - w_heater = 2 - - polygons = OrderedDict( - core=Polygon( - [ - (x_min - w_core / 2, 0), - (x_min - w_core / 2, h_core), - (x_min + w_core / 2, h_core), - (x_min + w_core / 2, 0), - ] - ), - clad=Polygon( - [ - (x_min - w_sim / 2, 0), - (x_min - w_sim / 2, h_clad), - (x_min + w_sim / 2, h_clad), - (x_min + w_sim / 2, 0), - ] - ), - box=Polygon( - [ - (x_min - w_sim / 2, 0), - (x_min - w_sim / 2, -h_box), - (x_min + w_sim / 2, -h_box), - (x_min + w_sim / 2, 0), - ] - ), - ) - - resolutions = dict(core={"resolution": 0.05, "distance": 1}) - - mesh_from_OrderedDict(polygons, resolutions, filename="mesh.msh", default_resolution_max=0.2) - - mesh = Mesh.load("mesh.msh") - basis = Basis(mesh, ElementTriN2() * ElementTriP2()) - basis0 = basis.with_element(ElementTriP0()) - epsilon = basis0.zeros(dtype=complex) - epsilon[basis0.get_dofs(elements="core")] = 3.4777**2 - epsilon[basis0.get_dofs(elements="clad")] = 1.444**2 - epsilon[basis0.get_dofs(elements="box")] = 1.444**2 - # basis0.plot(epsilon, colorbar=True).show() - - modes = compute_modes( - basis0, - epsilon, - wavelength=1.55, - mu_r=1, - num_modes=6, - order=2, - radius=3, - return_objects=True, - ) - print(modes) - print(modes[0].te_fraction) - - modes[0].show(np.real(modes[0].E)) - modes[0].show(np.imag(modes[0].E)) - - modes[0].show(np.real(modes[0].H)) - modes[0].show(np.imag(modes[0].H)) - - integrals = np.zeros((len(modes),) * 2, dtype=complex) - - for i in range(len(modes)): - for j in range(len(modes)): - integrals[i, j] = modes[i].calculate_overlap(modes[j]) - - plt.imshow(np.real(integrals)) - plt.colorbar() - plt.show() - - # Create basis to select a certain simulation extent - def sel_fun(x): - return (x[0] < 0) * (x[0] > -1) * (x[1] > 0) * (x[1] < 0.5) - - print(modes.sorted(lambda mode: mode.calculate_power(elements=sel_fun))) From 754889788d7b8b40edb06d2623436b7459acc901 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Fri, 7 Apr 2023 19:47:31 -0700 Subject: [PATCH 24/46] simplify compute_modes --- femwell/maxwell/waveguide.py | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index 1db4ad58..c2af814b 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -223,21 +223,17 @@ def bform(e_t, e_z, v_t, v_z, w): solver=solver(k=num_modes, sigma=sigma), ) - idx = np.abs(np.real(lams)).argsort()[::-1] - lams = lams[idx] - xs = xs[:, idx] - - xs = xs.T - print([np.sum(np.abs(xs[0, basis.split_indices()[i]])) for i in range(2)]) - xs[:, basis.split_indices()[1]] /= 1j * np.sqrt( - lams[:, np.newaxis] / k0**4 + xs[basis.split_indices()[1], :] /= 1j * np.sqrt( + lams[np.newaxis, :] / k0**4 ) # undo the scaling E_3,new = beta * E_3 hs = [] for i, lam in enumerate(lams): - H = calculate_hfield(basis, xs[i], np.sqrt(lam), omega=k0 * scipy.constants.speed_of_light) - power = calculate_overlap(basis, xs[i], H, basis, xs[i], H) - xs[i] /= np.sqrt(power) + H = calculate_hfield( + basis, xs[:, i], np.sqrt(lam), omega=k0 * scipy.constants.speed_of_light + ) + power = calculate_overlap(basis, xs[:, i], H, basis, xs[:, i], H) + xs[:, i] /= np.sqrt(power) H /= np.sqrt(power) hs.append(H) @@ -249,7 +245,7 @@ def bform(e_t, e_z, v_t, v_z, w): basis_epsilon_r=basis_epsilon_r, epsilon_r=epsilon_r, basis=basis, - E=xs[i], + E=xs[:, i], H=hs[i], ) for i in range(num_modes) From 625826ab0cfbefe6d51e44004e5b16f5b221a94d Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Fri, 7 Apr 2023 19:47:56 -0700 Subject: [PATCH 25/46] remove debug output --- femwell/maxwell/waveguide.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index c2af814b..a783ffd6 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -185,8 +185,6 @@ def compute_modes( basis = basis_epsilon_r.with_element(element) basis_epsilon_r = basis.with_element(basis_epsilon_r.elem) # adjust quadrature - print(k0) - @BilinearForm(dtype=epsilon_r.dtype) def aform(e_t, e_z, v_t, v_z, w): epsilon = w.epsilon * (1 + w.x[0] / radius) ** 2 From 8d4417f56322e8e7093ed0c858d28518f475894f Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Fri, 7 Apr 2023 19:56:11 -0700 Subject: [PATCH 26/46] add radius again --- femwell/maxwell/waveguide.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index a783ffd6..368bb647 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -316,7 +316,7 @@ def bform(e_t, e_z, v_t, v_z, w): mu_r=1, num_modes=6, order=2, - # radius=3 * scale, + radius=3 * scale, ) print(modes) print(modes[0].te_fraction) From 5673eb63b74c7c223c9592f7115f2ac5ca4ec441 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Fri, 7 Apr 2023 20:02:24 -0700 Subject: [PATCH 27/46] coulomb -> use object oriented mode solver --- docs/photonics/examples/culomb.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/photonics/examples/culomb.py b/docs/photonics/examples/culomb.py index 75bf4318..2b9c1b72 100644 --- a/docs/photonics/examples/culomb.py +++ b/docs/photonics/examples/culomb.py @@ -27,8 +27,8 @@ from tqdm import tqdm from femwell.culomb import solve_coulomb +from femwell.maxwell.waveguide import compute_modes from femwell.mesh import mesh_from_OrderedDict -from femwell.mode_solver import compute_modes # - @@ -108,8 +108,8 @@ epsilon **= 2 # basis_epsilon_r.plot(epsilon, colorbar=True).show() - neffs, basis_modes, modes = compute_modes(basis_epsilon_r, epsilon, 1.55, 1, 1, order=1) - voltages_neffs.append(neffs[0]) + modes = compute_modes(basis_epsilon_r, epsilon, wavelength=1.55, order=1) + voltages_neffs.append(modes[0].n_eff) # from mode_solver import plot_mode # plot_mode(basis_modes, modes[0]) From 6b007af0d1e2889089b1f73c62ce582f3617942d Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Fri, 7 Apr 2023 20:07:23 -0700 Subject: [PATCH 28/46] adjust examples to object oriented mode solver --- docs/photonics/examples/vary_wavelength.py | 6 ++---- docs/photonics/examples/vary_width.py | 4 +--- femwell/thermal.py | 12 ++++-------- 3 files changed, 7 insertions(+), 15 deletions(-) diff --git a/docs/photonics/examples/vary_wavelength.py b/docs/photonics/examples/vary_wavelength.py index 1490bf8e..b50073b0 100644 --- a/docs/photonics/examples/vary_wavelength.py +++ b/docs/photonics/examples/vary_wavelength.py @@ -28,8 +28,8 @@ from skfem.io.meshio import from_meshio from tqdm import tqdm +from femwell.maxwell.waveguide import compute_modes from femwell.mesh import mesh_from_OrderedDict -from femwell.mode_solver import compute_modes # - @@ -91,9 +91,7 @@ def n_silicon_dioxide(wavelength): }.items(): epsilon[basis0.get_dofs(elements=subdomain)] = n**2 - modes = compute_modes( - basis0, epsilon, wavelength=wavelength, num_modes=num_modes, normalize=False - ) + modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=num_modes) all_neffs[i] = np.real([mode.n_eff for mode in modes]) all_te_fracs[i, :] = [mode.te_fraction for mode in modes] # - diff --git a/docs/photonics/examples/vary_width.py b/docs/photonics/examples/vary_width.py index 516e3a0d..d0df6e08 100644 --- a/docs/photonics/examples/vary_width.py +++ b/docs/photonics/examples/vary_width.py @@ -56,9 +56,7 @@ for subdomain, n in {"core": 1.9963, "box": 1.444, "clad": 1}.items(): epsilon[basis0.get_dofs(elements=subdomain)] = n**2 - modes = compute_modes( - basis0, epsilon, wavelength=wavelength, num_modes=num_modes, normalize=False - ) + modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=num_modes) all_neffs[i] = np.real([mode.n_eff for mode in modes]) all_te_fracs[i, :] = [mode.te_fraction for mode in modes] # - diff --git a/femwell/thermal.py b/femwell/thermal.py index 67a348cf..e00cfd51 100644 --- a/femwell/thermal.py +++ b/femwell/thermal.py @@ -156,7 +156,6 @@ def unit_load(v, _): mesh_from_OrderedDict(polygons, resolutions, filename="mesh.msh", default_resolution_max=0.4) mesh = Mesh.load("mesh.msh") - print(mesh.boundaries) currents = np.linspace(0.007, 10e-3, 10) / polygons["heater"].area neffs = [] @@ -185,7 +184,7 @@ def unit_load(v, _): # basis.plot(temperature, colorbar=True) # plt.show() - from femwell.mode_solver import compute_modes, plot_mode + from femwell.maxwell.waveguide import compute_modes temperature0 = basis0.project(basis.interpolate(temperature)) epsilon = basis0.zeros() + (1.444 + 1.00e-5 * temperature0) ** 2 @@ -194,14 +193,11 @@ def unit_load(v, _): ) ** 2 # basis0.plot(epsilon, colorbar=True).show() - lams, basis, xs = compute_modes(basis0, epsilon, wavelength=1.55, mu_r=1, num_modes=5) + modes = compute_modes(basis0, epsilon, wavelength=1.55, num_modes=5) - print(lams) + # modes[0].show(modes[0].E) - # plot_mode(basis, xs[0]) - # plt.show() - - neffs.append(np.real(lams[0])) + neffs.append(np.real(modes[0].n_eff)) print(f"Phase shift: {2 * np.pi / 1.55 * (neffs[-1] - neffs[0]) * 320}") plt.xlabel("Power") From 41d882d85f7d0ca1a07c719ea23be2e3821019a1 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Fri, 7 Apr 2023 20:10:45 -0700 Subject: [PATCH 29/46] remove old example --- femwell/examples/overlap_integrals.py | 81 --------------------------- 1 file changed, 81 deletions(-) delete mode 100644 femwell/examples/overlap_integrals.py diff --git a/femwell/examples/overlap_integrals.py b/femwell/examples/overlap_integrals.py deleted file mode 100644 index 2a97d3bd..00000000 --- a/femwell/examples/overlap_integrals.py +++ /dev/null @@ -1,81 +0,0 @@ -import tempfile - -import matplotlib.pyplot as plt -import numpy as np -from skfem import Basis, ElementTriP0, Mesh -from tqdm.auto import tqdm - -from femwell.mode_solver import ( - calculate_hfield, - calculate_overlap, - compute_modes, - plot_mode, -) -from femwell.waveguide import mesh_waveguide - -widths = [0.5, 0.8] - -modes = [] - -for width in widths: - with tempfile.TemporaryDirectory() as tmpdirname: - mesh = mesh_waveguide( - wsim=2, - hclad=0.7, - hbox=0.5, - wcore=width, - hcore=0.22, - filename=f"{tmpdirname}/mesh.msh", - ) - mesh = Mesh.load(f"{tmpdirname}/mesh.msh") - - basis0 = Basis(mesh, ElementTriP0(), intorder=4) - epsilon = basis0.zeros() - epsilon[basis0.get_dofs(elements="core")] = 3.4777**2 - epsilon[basis0.get_dofs(elements="clad")] = 1.444**2 - epsilon[basis0.get_dofs(elements="box")] = 1.444**2 - # basis0.plot(epsilon, colorbar=True).show() - - lams, basis, xs = compute_modes(basis0, epsilon, wavelength=1.55, mu_r=1, num_modes=10) - - fig, axs = plot_mode(basis, np.real(xs[0]), colorbar=False) - plt.show() - - modes.append([lams, basis, xs]) - -for mode_a, mode_b in zip(modes[:-1], modes[1:]): - lams_a, basis_a, xs_a = mode_a - lams_b, basis_b, xs_b = mode_b - - integrals = np.zeros((len(lams_a), len(lams_b)), dtype=complex) - for i in range(len(lams_a)): - for j in range(len(lams_b)): - E_i = xs_a[i] - E_j = xs_b[j] - H_i = calculate_hfield(basis_a, E_i, lams_a[i] * (2 * np.pi / 1.55)) - H_j = calculate_hfield(basis_b, E_j, lams_b[j] * (2 * np.pi / 1.55)) - integrals[i, j] = calculate_overlap(basis_a, E_i, H_i, basis_b, E_j, H_j) - - plt.imshow(np.abs(integrals)) - plt.colorbar() - plt.show() - - z = np.abs(integrals) ** 2 - - fig, axCenter = plt.subplots(figsize=(8, 8)) - fig.subplots_adjust(0.05, 0.1, 0.95, 0.95) - - from mpl_toolkits.axes_grid1 import make_axes_locatable - - divider = make_axes_locatable(axCenter) - axvert = divider.append_axes("right", size="30%", pad=0.5) - axhoriz = divider.append_axes("top", size="20%", pad=0.25) - - axCenter.imshow(z, origin="lower", cmap="jet") - axhoriz.plot(range(np.shape(z)[1]), np.sum(z, 0)) - axvert.plot(np.sum(z, 1), range(np.shape(z)[0])) - - axhoriz.margins(x=0) - axvert.margins(y=0) - - plt.show() From 9d626e5d031849feae7bc5cfa8ece4976d6c7604 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Fri, 7 Apr 2023 20:12:07 -0700 Subject: [PATCH 30/46] adjust example to object oriented mode solver --- docs/photonics/examples/depletion_waveguide.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/photonics/examples/depletion_waveguide.py b/docs/photonics/examples/depletion_waveguide.py index 82341b50..14e531d8 100644 --- a/docs/photonics/examples/depletion_waveguide.py +++ b/docs/photonics/examples/depletion_waveguide.py @@ -23,8 +23,8 @@ from skfem import Basis, ElementTriP0 from skfem.io.meshio import from_meshio +from femwell.maxwell.waveguide import compute_modes from femwell.mesh import mesh_from_OrderedDict -from femwell.mode_solver import compute_modes from femwell.pn_analytical import index_pn_junction, k_to_alpha_dB # - @@ -106,8 +106,8 @@ def define_index(mesh, V, xpn, NA, ND, wavelength): neff_vs_V = [] for V in voltages: basis0, epsilon = define_index(mesh, V, xpn, NA, ND, wavelength) - lams, basis, xs = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=1, order=2) - neff_vs_V.append(lams) + modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=1, order=2) + neff_vs_V.append(modes[0].n_eff) # - # + tags=["hide-input"] From cb9b777d3a7f4d4f9452780e6615e37f9d6758a6 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Fri, 7 Apr 2023 20:14:36 -0700 Subject: [PATCH 31/46] update examples --- docs/photonics/examples/metal_heater_phase_shifter.py | 9 ++++----- femwell/waveguide.py | 7 +++---- 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/docs/photonics/examples/metal_heater_phase_shifter.py b/docs/photonics/examples/metal_heater_phase_shifter.py index 6247fea3..7629d38f 100644 --- a/docs/photonics/examples/metal_heater_phase_shifter.py +++ b/docs/photonics/examples/metal_heater_phase_shifter.py @@ -24,8 +24,8 @@ from skfem.io import from_meshio from tqdm import tqdm +from femwell.maxwell.waveguide import compute_modes from femwell.mesh import mesh_from_OrderedDict -from femwell.mode_solver import compute_modes, plot_mode from femwell.thermal import solve_thermal # - @@ -144,13 +144,12 @@ ) ** 2 # basis0.plot(epsilon, colorbar=True).show() - lams, basis, xs = compute_modes(basis0, epsilon, wavelength=1.55, num_modes=1) + modes = compute_modes(basis0, epsilon, wavelength=1.55, num_modes=1) if current_density == current_densities[-1]: - plot_mode(basis, xs[0].real) - plt.show() + modes[0].show(modes[0].E.real) - neffs.append(np.real(lams[0])) + neffs.append(np.real(modes[0].n_eff)) print(f"Phase shift: {2 * np.pi / 1.55 * (neffs[-1] - neffs[0]) * 320}") plt.xlabel("Current / mA") diff --git a/femwell/waveguide.py b/femwell/waveguide.py index d9558528..0d7bd0fc 100644 --- a/femwell/waveguide.py +++ b/femwell/waveguide.py @@ -8,8 +8,8 @@ from skfem import Basis, ElementTriP0, ElementVector, Mesh from skfem.io.meshio import from_meshio +from femwell.maxwell.waveguide import compute_modes from femwell.mesh import mesh_from_OrderedDict -from femwell.mode_solver import compute_modes, plot_mode def mesh_waveguide(filename, wsim, hclad, hbox, wcore, hcore): @@ -66,7 +66,6 @@ def mesh_waveguide(filename, wsim, hclad, hbox, wcore, hcore): epsilon[basis0.get_dofs(elements="box")] = 1.444**2 basis0.plot(epsilon, colorbar=True).show() - lams, basis, xs = compute_modes(basis0, epsilon, wavelength=1.55, mu_r=1, num_modes=5) + modes = compute_modes(basis0, epsilon, wavelength=1.55, mu_r=1, num_modes=5) - fig, axs = plot_mode(basis, xs[0], colorbar=False) - plt.show() + modes[0].show(modes[0].E.real) From 2ef9a56329c4aba22e92ebaba1181d0da3562644 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Sat, 8 Apr 2023 12:40:55 -0700 Subject: [PATCH 32/46] update benchmark --- docs/benchmarks/mode_solver.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/benchmarks/mode_solver.py b/docs/benchmarks/mode_solver.py index e9e25db1..ee2cbfd2 100644 --- a/docs/benchmarks/mode_solver.py +++ b/docs/benchmarks/mode_solver.py @@ -30,8 +30,8 @@ from skfem import Basis, ElementTriP0 from skfem.io.meshio import from_meshio +from femwell.maxwell.waveguide import compute_modes from femwell.mesh import mesh_from_OrderedDict -from femwell.mode_solver import compute_modes core_width = 3 core_thickness = 1 @@ -86,9 +86,9 @@ for subdomain, n in {"core": 3.44, "box": 3.40, "clad": 1}.items(): epsilon[basis0.get_dofs(elements=subdomain)] = n**2 - lams, basis, xs = compute_modes(basis0, epsilon, wavelength=1.15, num_modes=1, order=2) + modes = compute_modes(basis0, epsilon, wavelength=1.15, num_modes=1, order=2) - neff_values_femwell.append(np.real(lams[0])) + neff_values_femwell.append(np.real(modes[0].n_eff)) pd.DataFrame( { From 61629afcdb2275411fbca74596f3b0083bf800bf Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Sat, 8 Apr 2023 12:41:15 -0700 Subject: [PATCH 33/46] update waveguide_modes example --- docs/photonics/examples/waveguide_modes.py | 32 ++++++---------------- femwell/maxwell/waveguide.py | 11 +++++++- 2 files changed, 18 insertions(+), 25 deletions(-) diff --git a/docs/photonics/examples/waveguide_modes.py b/docs/photonics/examples/waveguide_modes.py index c4c65f43..b02b3fe5 100644 --- a/docs/photonics/examples/waveguide_modes.py +++ b/docs/photonics/examples/waveguide_modes.py @@ -26,14 +26,8 @@ from skfem import Basis, ElementTriP0 from skfem.io.meshio import from_meshio +from femwell.maxwell.waveguide import compute_modes from femwell.mesh import mesh_from_OrderedDict -from femwell.mode_solver import ( - calculate_hfield, - calculate_overlap, - compute_modes, - confinement_factor, - plot_mode, -) from femwell.visualization import plot_domains # - @@ -86,11 +80,10 @@ # + wavelength = 1.55 -lams, basis, xs = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=2, order=2) -for i, lam in enumerate(lams): - print(f"Effective refractive index: {lam:.4f}") - plot_mode(basis, xs[i].real, colorbar=True, direction="x") - plt.show() +modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=2, order=2) +for mode in modes: + print(f"Effective refractive index: {mode.n_eff:.4f}") + mode.show(mode.E.real, colorbar=True, direction="x") powers_in_waveguide = [] confinement_factors_waveguide = [] @@ -100,17 +93,8 @@ # What percentage of the mode is within the core for the calculated modes? # + -basis_waveguide = Basis(basis.mesh, basis.elem, elements="core") -basis0_waveguide = basis_waveguide.with_element(ElementTriP0()) -for i, lam in enumerate(lams): - H = calculate_hfield( - basis, xs[i], 2 * np.pi / wavelength * lam, omega=2 * np.pi / wavelength * speed_of_light - ) - powers_in_waveguide.append( - calculate_overlap(basis_waveguide, xs[i], H, basis_waveguide, xs[i], H) - ) - confinement_factors_waveguide.append( - confinement_factor(basis0_waveguide, epsilon, basis_waveguide, xs[i]) - ) +for mode in modes: + powers_in_waveguide.append(mode.calculate_power(elements="core")) + confinement_factors_waveguide.append(mode.calculate_confinement_factor(elements="core")) print(powers_in_waveguide) print(confinement_factors_waveguide) diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index 368bb647..0d803e57 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -34,6 +34,7 @@ calculate_hfield, calculate_overlap, calculate_scalar_product, + confinement_factor, plot_mode, ) @@ -119,7 +120,7 @@ def tm_fraction(self): def __repr__(self) -> str: return f"{self.__class__.__name__}(k: {self.k}, n_eff:{self.n_eff})" - def calculate_overlap(self, mode): + def calculate_overlap(self, mode, elements=None): return calculate_overlap(self.basis, self.E, self.H, mode.basis, mode.E, mode.H) def calculate_coupling_coefficient(self, mode, delta_epsilon): @@ -137,6 +138,14 @@ def calculate_power(self, elements=None): basis = self.basis.with_elements(elements) return calculate_overlap(basis, self.E, self.H, basis, self.E, self.H) + def calculate_confinement_factor(self, elements): + return confinement_factor( + self.basis_epsilon_r.with_elements(elements), + self.epsilon_r, + self.basis.with_elements(elements), + self.E, + ) + def plot(self, field, plot_vectors=False, colorbar=True, direction="y", title="E"): return plot_mode( self.basis, From 63f5446f75fa4601b0400911e71072d8f5622768 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Sun, 9 Apr 2023 16:50:43 -0700 Subject: [PATCH 34/46] update example --- .../metal_heater_phase_shifter_transient.py | 31 ++++------- femwell/maxwell/waveguide.py | 52 ++++++++++++------- 2 files changed, 42 insertions(+), 41 deletions(-) diff --git a/docs/photonics/examples/metal_heater_phase_shifter_transient.py b/docs/photonics/examples/metal_heater_phase_shifter_transient.py index 8965969a..276f8723 100644 --- a/docs/photonics/examples/metal_heater_phase_shifter_transient.py +++ b/docs/photonics/examples/metal_heater_phase_shifter_transient.py @@ -26,8 +26,8 @@ from skfem.io import from_meshio from tqdm import tqdm +from femwell.maxwell.waveguide import compute_modes from femwell.mesh import mesh_from_OrderedDict -from femwell.mode_solver import calculate_coupling_coefficient, compute_modes, plot_mode from femwell.thermal_transient import solve_thermal_transient # Simulating the TiN TOPS heater in https://doi.org/10.1364/OE.27.010456 @@ -163,7 +163,7 @@ def unit_load(v, w): epsilon_0 = basis0.zeros() + 1.444**2 epsilon_0[basis0.get_dofs(elements="core")] = 3.4777**2 -lams_0, basis_modes_0, xs_0 = compute_modes( +modes_0 = compute_modes( basis0, epsilon_0, wavelength=wavelength, mu_r=1, num_modes=1, solver="slepc" ) @@ -173,8 +173,6 @@ def unit_load(v, w): # basis.plot(temperature, vmin=0, vmax=np.max(temperatures)) # plt.show() - from femwell.mode_solver import compute_modes - temperature0 = basis0.project(basis.interpolate(temperature)) epsilon = basis0.zeros() + (1.444 + 1.00e-5 * temperature0) ** 2 epsilon[basis0.get_dofs(elements="core")] = ( @@ -182,27 +180,14 @@ def unit_load(v, w): ) ** 2 # basis0.plot(epsilon, colorbar=True).show() - lams, basis_modes, xs = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=1) + modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=1) # from femwell.mode_solver import plot_mode # plot_mode(basis_modes, xs[0]) # plt.show() - neffs.append(np.real(lams[0])) - neffs_approximated.append( - np.real( - lams_0[0] - + calculate_coupling_coefficient( - basis0, - (epsilon - epsilon_0) * scipy.constants.epsilon_0, - basis_modes_0, - xs_0[0], - xs_0[0], - ) - * scipy.constants.speed_of_light - * 0.5 - ) - ) + neffs.append(np.real(modes[0].n_eff)) + neffs_approximated.append(np.real(modes_0[0].calculate_pertubated_neff(epsilon - epsilon_0))) fig = plt.figure() ax = fig.add_subplot(111) @@ -211,10 +196,12 @@ def unit_load(v, w): ax.plot(times * 1e6, current(times) * 1000, "b-o") ax2 = ax.twinx() ax2.set_ylabel("Phase shift") -ax2.plot(times * 1e6, 2 * np.pi / wavelength * (neffs - lams_0[0]) * 320, "r-o", label="Exact") +ax2.plot( + times * 1e6, 2 * np.pi / wavelength * (neffs - modes_0[0].n_eff) * 320, "r-o", label="Exact" +) ax2.plot( times * 1e6, - 2 * np.pi / wavelength * (neffs_approximated - lams_0[0]) * 320, + 2 * np.pi / wavelength * (neffs_approximated - modes_0[0].n_eff) * 320, "g-o", label="Approximation", ) diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index 0d803e57..4190cf4a 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -39,24 +39,6 @@ ) -@dataclass(frozen=True) -class Modes: - modes: List - - def __getitem__(self, idx): - return self.modes[idx] - - def __len__(self): - return len(self.modes) - - def __repr__(self) -> str: - modes = "\n\t" + "\n\t".join(repr(mode) for mode in self.modes) + "\n" - return f"{self.__class__.__name__}(modes=({modes}))" - - def sorted(self, key): - return Modes(modes=sorted(self.modes, key=key)) - - @dataclass(frozen=True) class Mode: frequency: float @@ -146,6 +128,20 @@ def calculate_confinement_factor(self, elements): self.E, ) + def calculate_pertubated_neff(self, delta_epsilon): + return ( + self.n_eff + + calculate_coupling_coefficient( + self.basis_epsilon_r, + delta_epsilon * scipy.constants.epsilon_0, + self.basis, + self.E, + self.E, + ) + * scipy.constants.speed_of_light + * 0.5 + ) + def plot(self, field, plot_vectors=False, colorbar=True, direction="y", title="E"): return plot_mode( self.basis, @@ -161,6 +157,24 @@ def show(self, field, **kwargs): plt.show() +@dataclass(frozen=True) +class Modes: + modes: List + + def __getitem__(self, idx) -> Mode: + return self.modes[idx] + + def __len__(self) -> int: + return len(self.modes) + + def __repr__(self) -> str: + modes = "\n\t" + "\n\t".join(repr(mode) for mode in self.modes) + "\n" + return f"{self.__class__.__name__}(modes=({modes}))" + + def sorted(self, key): + return Modes(modes=sorted(self.modes, key=key)) + + def compute_modes( basis_epsilon_r, epsilon_r, @@ -172,7 +186,7 @@ def compute_modes( radius=np.inf, n_guess=None, solver="slepc", -): +) -> Modes: if solver == "scipy": solver = solver_eigen_scipy elif solver == "slepc": From dc47d6f9990cc54f3af9c40cc68e92fac5ffce56 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Sun, 9 Apr 2023 16:53:55 -0700 Subject: [PATCH 35/46] remove unnecessary parameter --- .../examples/metal_heater_phase_shifter_transient.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/docs/photonics/examples/metal_heater_phase_shifter_transient.py b/docs/photonics/examples/metal_heater_phase_shifter_transient.py index 276f8723..dbf2f78f 100644 --- a/docs/photonics/examples/metal_heater_phase_shifter_transient.py +++ b/docs/photonics/examples/metal_heater_phase_shifter_transient.py @@ -163,9 +163,7 @@ def unit_load(v, w): epsilon_0 = basis0.zeros() + 1.444**2 epsilon_0[basis0.get_dofs(elements="core")] = 3.4777**2 -modes_0 = compute_modes( - basis0, epsilon_0, wavelength=wavelength, mu_r=1, num_modes=1, solver="slepc" -) +modes_0 = compute_modes(basis0, epsilon_0, wavelength=wavelength, num_modes=1, solver="slepc") neffs = [] neffs_approximated = [] From 26a4c3882b7ce5f9b33706a0a99b976af6147514 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Sun, 9 Apr 2023 16:57:19 -0700 Subject: [PATCH 36/46] remove file which is already an example --- femwell/thermal_transient.py | 288 ----------------------------------- 1 file changed, 288 deletions(-) delete mode 100644 femwell/thermal_transient.py diff --git a/femwell/thermal_transient.py b/femwell/thermal_transient.py deleted file mode 100644 index 87fd539c..00000000 --- a/femwell/thermal_transient.py +++ /dev/null @@ -1,288 +0,0 @@ -from typing import Dict - -import numpy as np -from scipy.sparse.linalg import splu -from skfem import Basis, BilinearForm, ElementTriP0, LinearForm, Mesh, asm, enforce -from skfem.helpers import dot - -from femwell.thermal import solve_thermal - -""" -implemented like in https://www-users.cse.umn.edu/~arnold/8445.f11/notes.pdf page 81 in the middle of the page -""" - - -def solve_thermal_transient( - basis0, - thermal_conductivity_p0, - thermal_diffusivity_p0, - specific_conductivity: Dict[str, float], - current_densities_0, - current_densities, - fixed_boundaries, - dt, - steps, -): - basis, temperature = solve_thermal( - basis0, - thermal_conductivity_p0, - specific_conductivity, - {domain: current for domain, current in current_densities_0.items()}, - fixed_boundaries=fixed_boundaries, - ) - - @BilinearForm - def diffusivity_laplace(u, v, w): - return w["thermal_conductivity"] * dot(u.grad, v.grad) - - @BilinearForm - def mass(u, v, w): - return w["thermal_conductivity"] / w["thermal_diffusivity"] * u * v - - L = diffusivity_laplace.assemble( - basis, - thermal_conductivity=basis0.interpolate(thermal_conductivity_p0), - ) - M = mass.assemble( - basis, - thermal_diffusivity=basis0.interpolate(thermal_diffusivity_p0), - thermal_conductivity=basis0.interpolate(thermal_conductivity_p0), - ) - - theta = 0.5 # Crank–Nicolson - - x = basis.zeros() - for key, value in fixed_boundaries.items(): - x[basis.get_dofs(key)] = value - M0, L0 = enforce(M, L, D=basis.get_dofs(set(fixed_boundaries.keys()))) - A = M0 + theta * L0 * dt - B = M0 - (1 - theta) * L0 * dt - - backsolve = splu(A).solve # .T as splu prefers CSC - - t = 0 - temperatures = [temperature] - for i in range(steps): - joule_heating_rhs = basis.zeros() - for ( - domain, - current_density, - ) in current_densities.items(): # sum up the sources for the heating - current_density_p0 = basis0.zeros() - current_density_p0[basis0.get_dofs(elements=domain)] = current_density(t) - - @LinearForm - def joule_heating(v, w): - return w["current_density"] ** 2 / specific_conductivity[domain] * v - - joule_heating_rhs += asm( - joule_heating, - basis, - current_density=basis0.interpolate(current_density_p0), - ) - - t, temperature = t + dt, backsolve(B @ temperature + joule_heating_rhs * dt) - temperatures.append(temperature) - - return basis, temperatures - - -if __name__ == "__main__": - from collections import OrderedDict - - import matplotlib.pyplot as plt - import scipy.constants - from shapely.geometry import LineString, Polygon - from skfem.io import from_meshio - - from femwell.mesh import mesh_from_OrderedDict - - # Simulating the TiN TOPS heater in https://doi.org/10.1364/OE.27.010456 - - w_sim = 8 * 2 - h_clad = 2.8 - h_box = 1 - w_core = 0.5 - h_core = 0.22 - offset_heater = 2.2 - h_heater = 0.14 - w_heater = 2 - h_silicon = 3 - - wavelength = 1.55 - - polygons = OrderedDict( - bottom=LineString([(-w_sim / 2, -h_box), (w_sim / 2, -h_box)]), - core=Polygon( - [ - (-w_core / 2, 0), - (-w_core / 2, h_core), - (w_core / 2, h_core), - (w_core / 2, 0), - ] - ), - heater=Polygon( - [ - (-w_heater / 2, offset_heater), - (-w_heater / 2, offset_heater + h_heater), - (w_heater / 2, offset_heater + h_heater), - (w_heater / 2, offset_heater), - ] - ), - 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), - ] - ), - # silicon=Polygon([ - # (-w_sim / 2, - h_box - h_silicon), - # (-w_sim / 2, - h_box), - # (w_sim / 2, - h_box), - # (w_sim / 2, - h_box - h_silicon), - # ]), - ) - - resolutions = dict( - core={"resolution": 0.05, "distance": 1}, - clad={"resolution": 1, "distance": 1}, - box={"resolution": 1, "distance": 1}, - silicon={"resolution": 1, "distance": 1}, - heater={"resolution": 0.05, "distance": 1}, - ) - - mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=0.3)) - - basis0 = Basis(mesh, ElementTriP0(), intorder=4) - thermal_conductivity_p0 = basis0.zeros() - for domain, value in { - "core": 148, - "box": 1.38, - "clad": 1.38, - "heater": 28, - }.items(): # , 'silicon': 28 - thermal_conductivity_p0[basis0.get_dofs(elements=domain)] = value - thermal_conductivity_p0 *= 1e-12 # 1e-12 -> conversion from 1/m^2 -> 1/um^2 - - thermal_diffusivity_p0 = basis0.zeros() - for domain, value in { - "heater": 28 / 598 / 5240, - "box": 1.38 / 709 / 2203, - "clad": 1.38 / 709 / 2203, - "core": 148 / 711 / 2330, - # "silicon": 148 / 711 / 2330, - }.items(): - thermal_diffusivity_p0[basis0.get_dofs(elements=domain)] = value - thermal_diffusivity_p0 *= 1e12 # 1e-12 -> conversion from m^2 -> um^2 - - dt = 0.1e-5 - steps = 100 - current = ( - lambda t: 0.007 / polygons["heater"].area * ((t < dt * steps / 10) + (t > dt * steps / 2)) - ) - basis, temperatures = solve_thermal_transient( - basis0, - thermal_conductivity_p0, - thermal_diffusivity_p0, - specific_conductivity={"heater": 2.3e6}, - current_densities_0={"heater": current(0)}, - current_densities={"heater": current}, - fixed_boundaries={"bottom": 0}, - dt=dt, - steps=steps, - ) - - @LinearForm - def unit_load(v, w): - return v - - M = unit_load.assemble(basis) - - times = np.array([dt * i for i in range(steps + 1)]) - plt.xlabel("Time [us]") - plt.ylabel("Average temperature") - plt.plot(times * 1e6, M @ np.array(temperatures).T / np.sum(M)) - plt.show() - - # for i in range(0, steps, 10): - # fig, ax = plt.subplots(subplot_kw=dict(aspect=1)) - # for subdomain in mesh.subdomains.keys() - {'gmsh:bounding_entities'}: - # mesh.restrict(subdomain).draw(ax=ax, boundaries_only=True) - # basis.plot(temperatures[i], ax=ax, vmin=0, vmax=np.max(temperatures), shading='gouraud').show() - - # Calculate modes - - from tqdm.auto import tqdm - - from femwell.mode_solver import ( - calculate_coupling_coefficient, - compute_modes, - plot_mode, - ) - - epsilon_0 = basis0.zeros() + 1.444**2 - epsilon_0[basis0.get_dofs(elements="core")] = 3.4777**2 - lams_0, basis_modes_0, xs_0 = compute_modes( - basis0, epsilon_0, wavelength=wavelength, mu_r=1, num_modes=1 - ) - - neffs = [] - neffs_approximated = [] - for i, temperature in enumerate(tqdm(temperatures)): - # basis.plot(temperature, vmin=0, vmax=np.max(temperatures)) - # plt.show() - - temperature0 = basis0.project(basis.interpolate(temperature)) - epsilon = basis0.zeros() + (1.444 + 1.00e-5 * temperature0) ** 2 - epsilon[basis0.get_dofs(elements="core")] = ( - 3.4777 + 1.86e-4 * temperature0[basis0.get_dofs(elements="core")] - ) ** 2 - # basis0.plot(epsilon, colorbar=True).show() - - lams, basis_modes, xs = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=1) - - # plot_mode(basis_modes, xs[0]) - # plt.show() - - neffs.append(np.real(lams[0])) - neffs_approximated.append( - np.real( - lams_0[0] - + calculate_coupling_coefficient( - basis0, - (epsilon - epsilon_0) * scipy.constants.epsilon_0, - basis_modes_0, - xs_0[0], - xs_0[0], - ) - * scipy.constants.speed_of_light - * 2e-3 - ) - ) - - fig = plt.figure() - ax = fig.add_subplot(111) - ax.set_xlabel("Time [us]") - ax.set_ylabel("Current [mA]") - ax.plot(times * 1e6, current(times), "b-o") - ax2 = ax.twinx() - ax2.set_ylabel("Phase shift") - ax2.plot(times * 1e6, 2 * np.pi / 1.55 * (neffs - lams_0[0]) * 320, "r-o", label="Exact") - ax2.plot( - times * 1e6, - 2 * np.pi / 1.55 * (neffs_approximated - lams_0[0]) * 320, - "g-o", - label="Approximation", - ) - plt.legend() - plt.show() From 749b6974d8bb23c86c6f3049c56738e4b6c66bd2 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Sun, 9 Apr 2023 17:04:05 -0700 Subject: [PATCH 37/46] update coplanar_waveguide.py --- femwell/examples/coplanar_waveguide.py | 25 ++++++++----------------- femwell/maxwell/waveguide.py | 4 ++++ 2 files changed, 12 insertions(+), 17 deletions(-) diff --git a/femwell/examples/coplanar_waveguide.py b/femwell/examples/coplanar_waveguide.py index 6dec63b8..543142ff 100644 --- a/femwell/examples/coplanar_waveguide.py +++ b/femwell/examples/coplanar_waveguide.py @@ -9,8 +9,8 @@ from shapely.geometry import LineString, box from skfem import Basis, ElementTriP0, Mesh +from femwell.maxwell.waveguide import compute_modes from femwell.mesh import mesh_from_OrderedDict -from femwell.mode_solver import calculate_hfield, compute_modes, plot_mode def mesh_waveguide(filename, wsim, hclad, hsi, wcore, hcore): @@ -109,7 +109,7 @@ def mesh_coax(filename, radius_inner, radius_outer): basis0.plot(np.real(epsilon), colorbar=True).show() conductors = ["isolator2___isolator"] - lams, basis, xs = compute_modes( + modes = compute_modes( basis0, epsilon, wavelength=scipy.constants.speed_of_light / frequency * 1e3, @@ -117,10 +117,9 @@ def mesh_coax(filename, radius_inner, radius_outer): num_modes=len(conductors), metallic_boundaries=True, ) - print("propagation constants", 1 / lams) + print("propagation constants", 1 / modes.n_effs) - fig, axs = plot_mode(basis, np.real(xs[0]), plot_vectors=True) - plt.show() + modes[0].show(modes[0].E.real, plot_vectors=True) from skfem import * from skfem.helpers import * @@ -129,20 +128,12 @@ def mesh_coax(filename, radius_inner, radius_outer): def current_form(w): return inner(np.array([w.n[1], -w.n[0]]), w.H) - currents = np.zeros((len(conductors), len(lams))) + currents = np.zeros((len(conductors), len(modes))) - for mode_i in range(len(lams)): - xbs = calculate_hfield( - basis, - xs[mode_i], - lams[mode_i] * (2 * np.pi / (scipy.constants.speed_of_light / frequency * 1e3)), - omega=2 * np.pi * frequency * 1e-3, - ) + for mode_i, mode in enumerate(modes): + mode.show(mode.H.real, plot_vectors=True) - fig, axs = plot_mode(basis, np.real(xbs), plot_vectors=True) - plt.show() - - (ht, ht_basis), (hz, hz_basis) = basis.split(xbs) + (ht, ht_basis), (hz, hz_basis) = mode.basis.split(mode.H) for conductors_i, conductor in enumerate(conductors): facet_basis = ht_basis.boundary(facets=mesh.boundaries[conductor]) current = abs(current_form.assemble(facet_basis, H=facet_basis.interpolate(ht))) diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index 4190cf4a..9990a9fe 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -174,6 +174,10 @@ def __repr__(self) -> str: def sorted(self, key): return Modes(modes=sorted(self.modes, key=key)) + @property + def n_effs(self): + return np.array([mode.n_eff for mode in self.modes]) + def compute_modes( basis_epsilon_r, From b03291b4bdf075378a9945f68c4f9eb27a935bc5 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Sun, 9 Apr 2023 17:05:33 -0700 Subject: [PATCH 38/46] module-level warning in mode_solver --- femwell/mode_solver.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index 4215595c..d53e7b37 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -1,5 +1,7 @@ """Waveguide analysis based on https://doi.org/10.1080/02726340290084012.""" +import warnings + import matplotlib.pyplot as plt import numpy as np import scipy.constants @@ -23,6 +25,11 @@ from skfem.helpers import cross, curl, dot, grad, inner from skfem.utils import solver_eigen_scipy +warnings.warn( + "femwell.maxwell.waveguide will replace this module with version 0.1.0.", + DeprecationWarning, +) + def compute_modes( basis_epsilon_r, @@ -120,12 +127,6 @@ def bform(e_t, e_z, v_t, v_z, w): H /= np.sqrt(power) hs.append(H) - import warnings - - warnings.warn( - "femwell.maxwell.waveguide will replace this module with version 0.1.0.", - DeprecationWarning, - ) return np.sqrt(lams)[:num_modes] / k0, basis, xs[:num_modes] From d11eac14b0c9e4e8c634c0122005a345fc4a934c Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Sun, 9 Apr 2023 17:10:43 -0700 Subject: [PATCH 39/46] keyword only parameters in compute_modes --- femwell/maxwell/waveguide.py | 1 + 1 file changed, 1 insertion(+) diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index 9990a9fe..e09ecbca 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -183,6 +183,7 @@ def compute_modes( basis_epsilon_r, epsilon_r, wavelength, + *, mu_r=1, num_modes=1, order=1, From 648a55bf31f071c079d713bf5ccd135006cd6adc Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Sun, 9 Apr 2023 17:31:50 -0700 Subject: [PATCH 40/46] move warining back into compute_modes as the functions are used in maxwell.waveguide --- femwell/mode_solver.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py index d53e7b37..ea48ef17 100644 --- a/femwell/mode_solver.py +++ b/femwell/mode_solver.py @@ -1,7 +1,5 @@ """Waveguide analysis based on https://doi.org/10.1080/02726340290084012.""" -import warnings - import matplotlib.pyplot as plt import numpy as np import scipy.constants @@ -25,11 +23,6 @@ from skfem.helpers import cross, curl, dot, grad, inner from skfem.utils import solver_eigen_scipy -warnings.warn( - "femwell.maxwell.waveguide will replace this module with version 0.1.0.", - DeprecationWarning, -) - def compute_modes( basis_epsilon_r, @@ -45,6 +38,12 @@ def compute_modes( normalize=True, cache_path=None, ): + import warnings + + warnings.warn( + "femwell.maxwell.waveguide will replace this module with version 0.1.0.", + DeprecationWarning, + ) if solver == "scipy": solver = solver_eigen_scipy elif solver == "slepc": From 2027eddfa0b727e0ff81a810f3af45227291e69c Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Sun, 9 Apr 2023 21:47:01 -0700 Subject: [PATCH 41/46] allow specifing metallic bounaries (probably a temoorary solution) --- femwell/maxwell/waveguide.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index e09ecbca..d8dc61d0 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -239,7 +239,12 @@ def bform(e_t, e_z, v_t, v_z, w): if metallic_boundaries: lams, xs = solve( - *condense(-A, -B, D=basis.get_dofs(), x=basis.zeros(dtype=complex)), + *condense( + -A, + -B, + D=basis.get_dofs(None if metallic_boundaries is True else metallic_boundaries), + x=basis.zeros(dtype=complex), + ), solver=solver(k=num_modes, sigma=sigma), ) else: From 28aa4aefaad84946b86f8e4c3ebf2af5d6aeb47c Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Sun, 9 Apr 2023 21:54:48 -0700 Subject: [PATCH 42/46] remove non-object oriented mode solver --- femwell/{eme.py => eme.py.bak} | 0 femwell/maxwell/waveguide.py | 202 ++++++++++++++++++- femwell/mode_solver.py | 354 --------------------------------- 3 files changed, 192 insertions(+), 364 deletions(-) rename femwell/{eme.py => eme.py.bak} (100%) delete mode 100644 femwell/mode_solver.py diff --git a/femwell/eme.py b/femwell/eme.py.bak similarity index 100% rename from femwell/eme.py rename to femwell/eme.py.bak diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index d8dc61d0..9c18586a 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -28,16 +28,6 @@ from skfem.helpers import cross, curl, dot, grad, inner from skfem.utils import solver_eigen_scipy -from femwell.mode_solver import ( - calculate_coupling_coefficient, - calculate_energy_current_density, - calculate_hfield, - calculate_overlap, - calculate_scalar_product, - confinement_factor, - plot_mode, -) - @dataclass(frozen=True) class Mode: @@ -284,6 +274,198 @@ def bform(e_t, e_z, v_t, v_z, w): ) +def calculate_hfield(basis, xs, beta, omega=1): + @BilinearForm(dtype=np.complex64) + def aform(e_t, e_z, v_t, v_z, w): + return ( + (-1j * beta * e_t[1] + e_z.grad[1]) * v_t[0] + + (1j * beta * e_t[0] - e_z.grad[0]) * v_t[1] + + e_t.curl * v_z + ) + + @BilinearForm(dtype=np.complex64) + def bform(e_t, e_z, v_t, v_z, w): + return dot(e_t, v_t) + e_z * v_z + + return ( + scipy.sparse.linalg.spsolve( + bform.assemble(basis), aform.assemble(basis) @ xs.astype(complex) + ) + * -1j + / scipy.constants.mu_0 + / omega + ) + + +def calculate_energy_current_density(basis, xs): + basis_energy = basis.with_element(ElementTriP0()) + + @LinearForm(dtype=complex) + def aform(v, w): + e_t, e_z = w["e"] + return abs(e_t[0]) ** 2 * v + abs(e_t[1]) ** 2 * v + abs(e_z) * v + + a_operator = aform.assemble(basis_energy, e=basis.interpolate(xs)) + + @BilinearForm(dtype=complex) + def bform(e, v, w): + return e * v + + b_operator = bform.assemble(basis_energy) + + return basis_energy, scipy.sparse.linalg.spsolve(b_operator, a_operator) + + +def calculate_overlap(basis_i, E_i, H_i, basis_j, E_j, H_j): + @Functional(dtype=np.complex64) + def overlap(w): + return cross(np.conj(w["E_i"][0]), w["H_j"][0]) + cross(w["E_j"][0], np.conj(w["H_i"][0])) + + if basis_i == basis_j: + return 0.5 * overlap.assemble( + basis_i, + E_i=basis_i.interpolate(E_i), + H_i=basis_i.interpolate(H_i), + E_j=basis_j.interpolate(E_j), + H_j=basis_j.interpolate(H_j), + ) + basis_j_fix = basis_j.with_element(ElementVector(ElementTriP1())) + + (et, et_basis), (ez, ez_basis) = basis_j.split(E_j) + E_j = basis_j_fix.project(et_basis.interpolate(et), dtype=np.cfloat) + (et_x, et_x_basis), (et_y, et_y_basis) = basis_j_fix.split(E_j) + + (et, et_basis), (ez, ez_basis) = basis_j.split(H_j) + H_j = basis_j_fix.project(et_basis.interpolate(et), dtype=np.cfloat) + (ht_x, ht_x_basis), (ht_y, ht_y_basis) = basis_j_fix.split(H_j) + + @Functional(dtype=np.complex64) + def overlap(w): + return cross( + np.conj(w["E_i"][0]), + np.array((ht_x_basis.interpolator(ht_x)(w.x), ht_y_basis.interpolator(ht_y)(w.x))), + ) + cross( + np.array((et_x_basis.interpolator(et_x)(w.x), et_y_basis.interpolator(et_y)(w.x))), + np.conj(w["H_i"][0]), + ) + + return 0.5 * overlap.assemble( + basis_i, E_i=basis_i.interpolate(E_i), H_i=basis_i.interpolate(H_i) + ) + + +def calculate_scalar_product(basis_i, E_i, basis_j, H_j): + @Functional + def overlap(w): + return cross(np.conj(w["E_i"][0]), w["H_j"][0]) + + if basis_i == basis_j: + return overlap.assemble(basis_i, E_i=basis_i.interpolate(E_i), H_j=basis_j.interpolate(H_j)) + basis_j_fix = basis_j.with_element(ElementVector(ElementTriP1())) + + (et, et_basis), (ez, ez_basis) = basis_j.split(H_j) + H_j = basis_j_fix.project(et_basis.interpolate(et), dtype=np.cfloat) + (ht_x, ht_x_basis), (ht_y, ht_y_basis) = basis_j_fix.split(H_j) + + @Functional(dtype=np.complex64) + def overlap(w): + return cross( + np.conj(w["E_i"][0]), + np.array((ht_x_basis.interpolator(ht_x)(w.x), ht_y_basis.interpolator(ht_y)(w.x))), + ) + + return overlap.assemble(basis_i, E_i=basis_i.interpolate(E_i)) + + +def calculate_coupling_coefficient(basis_epsilon, delta_epsilon, basis, E_i, E_j): + @Functional(dtype=complex) + def overlap(w): + return w["delta_epsilon"] * ( + dot(np.conj(w["E_i"][0]), w["E_j"][0]) + np.conj(w["E_i"][1]) * w["E_j"][1] + ) + + return overlap.assemble( + basis, + E_i=basis.interpolate(E_i), + E_j=basis.interpolate(E_j), + delta_epsilon=basis_epsilon.interpolate(delta_epsilon), + ) + + +def confinement_factor(basis_epsilon, epsilon, basis, E): + @Functional + def factor(w): + return ( + ( + np.sqrt(w["epsilon"]) + * (dot(np.conj(w["E"][0]), w["E"][0]) + np.conj(w["E"][1]) * w["E"][1]) + ) + * speed_of_light + * epsilon_0 + ) + + return factor.assemble( + basis, + E=basis.interpolate(E), + epsilon=basis_epsilon.interpolate(epsilon), + ) + + +def plot_mode(basis, mode, plot_vectors=False, colorbar=True, title="E", direction="y"): + from mpl_toolkits.axes_grid1 import make_axes_locatable + + (et, et_basis), (ez, ez_basis) = basis.split(mode) + + if plot_vectors: + rc = (2, 1) if direction != "x" else (1, 2) + fig, axs = plt.subplots(*rc, subplot_kw=dict(aspect=1)) + for ax in axs: + basis.mesh.draw(ax=ax, boundaries=True, boundaries_only=True) + for subdomain in basis.mesh.subdomains.keys() - {"gmsh:bounding_entities"}: + basis.mesh.restrict(subdomain).draw(ax=ax, boundaries_only=True) + et_basis.plot(et, ax=axs[0]) + ez_basis.plot(ez, ax=axs[1]) + + divider = make_axes_locatable(axs[1]) + cax = divider.append_axes("right", size="5%", pad=0.05) + plt.colorbar(axs[1].collections[0], cax=cax) + return fig, axs + + plot_basis = et_basis.with_element(ElementVector(ElementDG(ElementTriP1()))) + et_xy = plot_basis.project(et_basis.interpolate(et)) + (et_x, et_x_basis), (et_y, et_y_basis) = plot_basis.split(et_xy) + + rc = (3, 1) if direction != "x" else (1, 3) + fig, axs = plt.subplots(*rc, subplot_kw=dict(aspect=1)) + for ax in axs: + basis.mesh.draw(ax=ax, boundaries=True, boundaries_only=True) + for subdomain in basis.mesh.subdomains.keys() - {"gmsh:bounding_entities"}: + basis.mesh.restrict(subdomain).draw(ax=ax, boundaries_only=True) + + for ax, component in zip(axs, "xyz"): + ax.set_title(f"${title}_{component}$") + + maxabs = max(np.max(np.abs(data.value)) for data in basis.interpolate(mode)) + vmin = -maxabs if colorbar == "same" else None + vmax = maxabs if colorbar == "same" else None + + et_x_basis.plot(et_x, shading="gouraud", ax=axs[0], vmin=vmin, vmax=vmax) + et_y_basis.plot(et_y, shading="gouraud", ax=axs[1], vmin=vmin, vmax=vmax) + ez_basis.plot(ez, shading="gouraud", ax=axs[2], vmin=vmin, vmax=vmax) + + if colorbar: + if colorbar == "same": + plt.colorbar(axs[0].collections[-1], ax=axs.ravel().tolist()) + else: + for ax in axs: + divider = make_axes_locatable(ax) + cax = divider.append_axes("right", size="5%", pad=0.05) + plt.colorbar(ax.collections[-1], cax=cax) + plt.tight_layout() + + return fig, axs + + if __name__ == "__main__": from collections import OrderedDict diff --git a/femwell/mode_solver.py b/femwell/mode_solver.py deleted file mode 100644 index ea48ef17..00000000 --- a/femwell/mode_solver.py +++ /dev/null @@ -1,354 +0,0 @@ -"""Waveguide analysis based on https://doi.org/10.1080/02726340290084012.""" - -import matplotlib.pyplot as plt -import numpy as np -import scipy.constants -import scipy.sparse.linalg -from scipy.constants import epsilon_0, speed_of_light -from skfem import ( - Basis, - BilinearForm, - ElementDG, - ElementTriN1, - ElementTriN2, - ElementTriP0, - ElementTriP1, - ElementTriP2, - ElementVector, - Functional, - LinearForm, - condense, - solve, -) -from skfem.helpers import cross, curl, dot, grad, inner -from skfem.utils import solver_eigen_scipy - - -def compute_modes( - basis_epsilon_r, - epsilon_r, - wavelength, - mu_r=1, - num_modes=1, - order=1, - metallic_boundaries=False, - radius=np.inf, - n_guess=None, - solver="slepc", - normalize=True, - cache_path=None, -): - import warnings - - warnings.warn( - "femwell.maxwell.waveguide will replace this module with version 0.1.0.", - DeprecationWarning, - ) - if solver == "scipy": - solver = solver_eigen_scipy - elif solver == "slepc": - from femwell.solver import solver_eigen_slepc - - solver = solver_eigen_slepc - else: - raise ValueError("`solver` must either be `scipy` or `slepc`") - - if cache_path: - from femwell.solver import solver_cached - - solver = solver_cached(solver, cache_path) - - k0 = 2 * np.pi / wavelength - - if order == 1: - element = ElementTriN1() * ElementTriP1() - elif order == 2: - element = ElementTriN2() * ElementTriP2() - else: - raise AssertionError("Only order 1 and 2 implemented by now.") - - basis = basis_epsilon_r.with_element(element) - basis_epsilon_r = basis.with_element(basis_epsilon_r.elem) # adjust quadrature - - @BilinearForm(dtype=epsilon_r.dtype) - def aform(e_t, e_z, v_t, v_z, w): - epsilon = w.epsilon * (1 + w.x[0] / radius) ** 2 - - return ( - 1 / mu_r * curl(e_t) * curl(v_t) - - k0**2 * epsilon * dot(e_t, v_t) - + 1 / mu_r * dot(grad(e_z), v_t) - + epsilon * inner(e_t, grad(v_z)) - - epsilon * e_z * v_z - ) - - @BilinearForm(dtype=epsilon_r.dtype) - def bform(e_t, e_z, v_t, v_z, w): - return -1 / mu_r * dot(e_t, v_t) - - A = aform.assemble(basis, epsilon=basis_epsilon_r.interpolate(epsilon_r)) - B = bform.assemble(basis, epsilon=basis_epsilon_r.interpolate(epsilon_r)) - - if n_guess: - sigma = sigma = k0**2 * n_guess**2 - else: - sigma = sigma = k0**2 * np.max(epsilon_r) ** 2 - - if metallic_boundaries: - lams, xs = solve( - *condense(-A, -B, D=basis.get_dofs(), x=basis.zeros(dtype=complex)), - solver=solver(k=num_modes, sigma=sigma), - ) - else: - lams, xs = solve( - -A, - -B, - solver=solver(k=num_modes, sigma=sigma), - ) - - idx = np.abs(np.real(lams)).argsort()[::-1] - lams = lams[idx] - xs = xs[:, idx] - - xs = xs.T - xs[:, basis.split_indices()[1]] /= 1j * np.sqrt( - lams[:, np.newaxis] - ) # undo the scaling E_3,new = beta * E_3 - - hs = [] - if normalize: - for i, lam in enumerate(lams): - H = calculate_hfield( - basis, xs[i], np.sqrt(lam), omega=k0 * scipy.constants.speed_of_light - ) - power = calculate_overlap(basis, xs[i], H, basis, xs[i], H) - xs[i] /= np.sqrt(power) - H /= np.sqrt(power) - hs.append(H) - - return np.sqrt(lams)[:num_modes] / k0, basis, xs[:num_modes] - - -def calculate_hfield(basis, xs, beta, omega=1): - @BilinearForm(dtype=np.complex64) - def aform(e_t, e_z, v_t, v_z, w): - return ( - (-1j * beta * e_t[1] + e_z.grad[1]) * v_t[0] - + (1j * beta * e_t[0] - e_z.grad[0]) * v_t[1] - + e_t.curl * v_z - ) - - @BilinearForm(dtype=np.complex64) - def bform(e_t, e_z, v_t, v_z, w): - return dot(e_t, v_t) + e_z * v_z - - return ( - scipy.sparse.linalg.spsolve( - bform.assemble(basis), aform.assemble(basis) @ xs.astype(complex) - ) - * -1j - / scipy.constants.mu_0 - / omega - ) - - -def calculate_energy_current_density(basis, xs): - basis_energy = basis.with_element(ElementTriP0()) - - @LinearForm(dtype=complex) - def aform(v, w): - e_t, e_z = w["e"] - return abs(e_t[0]) ** 2 * v + abs(e_t[1]) ** 2 * v + abs(e_z) * v - - a_operator = aform.assemble(basis_energy, e=basis.interpolate(xs)) - - @BilinearForm(dtype=complex) - def bform(e, v, w): - return e * v - - b_operator = bform.assemble(basis_energy) - - return basis_energy, scipy.sparse.linalg.spsolve(b_operator, a_operator) - - -def calculate_overlap(basis_i, E_i, H_i, basis_j, E_j, H_j): - @Functional(dtype=np.complex64) - def overlap(w): - return cross(np.conj(w["E_i"][0]), w["H_j"][0]) + cross(w["E_j"][0], np.conj(w["H_i"][0])) - - if basis_i == basis_j: - return 0.5 * overlap.assemble( - basis_i, - E_i=basis_i.interpolate(E_i), - H_i=basis_i.interpolate(H_i), - E_j=basis_j.interpolate(E_j), - H_j=basis_j.interpolate(H_j), - ) - basis_j_fix = basis_j.with_element(ElementVector(ElementTriP1())) - - (et, et_basis), (ez, ez_basis) = basis_j.split(E_j) - E_j = basis_j_fix.project(et_basis.interpolate(et), dtype=np.cfloat) - (et_x, et_x_basis), (et_y, et_y_basis) = basis_j_fix.split(E_j) - - (et, et_basis), (ez, ez_basis) = basis_j.split(H_j) - H_j = basis_j_fix.project(et_basis.interpolate(et), dtype=np.cfloat) - (ht_x, ht_x_basis), (ht_y, ht_y_basis) = basis_j_fix.split(H_j) - - @Functional(dtype=np.complex64) - def overlap(w): - return cross( - np.conj(w["E_i"][0]), - np.array((ht_x_basis.interpolator(ht_x)(w.x), ht_y_basis.interpolator(ht_y)(w.x))), - ) + cross( - np.array((et_x_basis.interpolator(et_x)(w.x), et_y_basis.interpolator(et_y)(w.x))), - np.conj(w["H_i"][0]), - ) - - return 0.5 * overlap.assemble( - basis_i, E_i=basis_i.interpolate(E_i), H_i=basis_i.interpolate(H_i) - ) - - -def calculate_scalar_product(basis_i, E_i, basis_j, H_j): - @Functional - def overlap(w): - return cross(np.conj(w["E_i"][0]), w["H_j"][0]) - - if basis_i == basis_j: - return overlap.assemble(basis_i, E_i=basis_i.interpolate(E_i), H_j=basis_j.interpolate(H_j)) - basis_j_fix = basis_j.with_element(ElementVector(ElementTriP1())) - - (et, et_basis), (ez, ez_basis) = basis_j.split(H_j) - H_j = basis_j_fix.project(et_basis.interpolate(et), dtype=np.cfloat) - (ht_x, ht_x_basis), (ht_y, ht_y_basis) = basis_j_fix.split(H_j) - - @Functional(dtype=np.complex64) - def overlap(w): - return cross( - np.conj(w["E_i"][0]), - np.array((ht_x_basis.interpolator(ht_x)(w.x), ht_y_basis.interpolator(ht_y)(w.x))), - ) - - return overlap.assemble(basis_i, E_i=basis_i.interpolate(E_i)) - - -def calculate_coupling_coefficient(basis_epsilon, delta_epsilon, basis, E_i, E_j): - @Functional(dtype=complex) - def overlap(w): - return w["delta_epsilon"] * ( - dot(np.conj(w["E_i"][0]), w["E_j"][0]) + np.conj(w["E_i"][1]) * w["E_j"][1] - ) - - return overlap.assemble( - basis, - E_i=basis.interpolate(E_i), - E_j=basis.interpolate(E_j), - delta_epsilon=basis_epsilon.interpolate(delta_epsilon), - ) - - -def confinement_factor(basis_epsilon, epsilon, basis, E): - @Functional - def factor(w): - return ( - ( - np.sqrt(w["epsilon"]) - * (dot(np.conj(w["E"][0]), w["E"][0]) + np.conj(w["E"][1]) * w["E"][1]) - ) - * speed_of_light - * epsilon_0 - ) - - return factor.assemble( - basis, - E=basis.interpolate(E), - epsilon=basis_epsilon.interpolate(epsilon), - ) - - -def calculate_te_frac(basis, x): - @Functional - def ex(w): - return np.abs(w.E[0][0]) ** 2 - - @Functional - def ey(w): - return np.abs(w.E[0][1]) ** 2 - - ex_sum = ex.assemble(basis, E=basis.interpolate(x)) - ey_sum = ey.assemble(basis, E=basis.interpolate(x)) - - return ex_sum / (ex_sum + ey_sum) - - -def plot_mode(basis, mode, plot_vectors=False, colorbar=True, title="E", direction="y"): - from mpl_toolkits.axes_grid1 import make_axes_locatable - - (et, et_basis), (ez, ez_basis) = basis.split(mode) - - if plot_vectors: - rc = (2, 1) if direction != "x" else (1, 2) - fig, axs = plt.subplots(*rc, subplot_kw=dict(aspect=1)) - for ax in axs: - basis.mesh.draw(ax=ax, boundaries=True, boundaries_only=True) - for subdomain in basis.mesh.subdomains.keys() - {"gmsh:bounding_entities"}: - basis.mesh.restrict(subdomain).draw(ax=ax, boundaries_only=True) - et_basis.plot(et, ax=axs[0]) - ez_basis.plot(ez, ax=axs[1]) - - divider = make_axes_locatable(axs[1]) - cax = divider.append_axes("right", size="5%", pad=0.05) - plt.colorbar(axs[1].collections[0], cax=cax) - return fig, axs - - plot_basis = et_basis.with_element(ElementVector(ElementDG(ElementTriP1()))) - et_xy = plot_basis.project(et_basis.interpolate(et)) - (et_x, et_x_basis), (et_y, et_y_basis) = plot_basis.split(et_xy) - - rc = (3, 1) if direction != "x" else (1, 3) - fig, axs = plt.subplots(*rc, subplot_kw=dict(aspect=1)) - for ax in axs: - basis.mesh.draw(ax=ax, boundaries=True, boundaries_only=True) - for subdomain in basis.mesh.subdomains.keys() - {"gmsh:bounding_entities"}: - basis.mesh.restrict(subdomain).draw(ax=ax, boundaries_only=True) - - for ax, component in zip(axs, "xyz"): - ax.set_title(f"${title}_{component}$") - - maxabs = max(np.max(np.abs(data.value)) for data in basis.interpolate(mode)) - vmin = -maxabs if colorbar == "same" else None - vmax = maxabs if colorbar == "same" else None - - et_x_basis.plot(et_x, shading="gouraud", ax=axs[0], vmin=vmin, vmax=vmax) - et_y_basis.plot(et_y, shading="gouraud", ax=axs[1], vmin=vmin, vmax=vmax) - ez_basis.plot(ez, shading="gouraud", ax=axs[2], vmin=vmin, vmax=vmax) - - if colorbar: - if colorbar == "same": - plt.colorbar(axs[0].collections[-1], ax=axs.ravel().tolist()) - else: - for ax in axs: - divider = make_axes_locatable(ax) - cax = divider.append_axes("right", size="5%", pad=0.05) - plt.colorbar(ax.collections[-1], cax=cax) - plt.tight_layout() - - return fig, axs - - -def argsort_modes_by_power_in_elements(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(modes[0].basis.mesh, modes[0].basis.elem, elements=elements) - - overlaps = [ - calculate_overlap(selection_basis, mode.E, mode.H, selection_basis, mode.E, mode.H) - for mode in modes - ] - - return np.argsort(np.abs(overlaps))[::-1] From 5769cd2bd67c2a3999a0a52516a2dcce1f1a3adb Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Sun, 9 Apr 2023 21:58:32 -0700 Subject: [PATCH 43/46] move confinement factor --- femwell/maxwell/waveguide.py | 40 +++++++++++++++--------------------- 1 file changed, 16 insertions(+), 24 deletions(-) diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index 9c18586a..0c5b3491 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -111,11 +111,22 @@ def calculate_power(self, elements=None): return calculate_overlap(basis, self.E, self.H, basis, self.E, self.H) def calculate_confinement_factor(self, elements): - return confinement_factor( - self.basis_epsilon_r.with_elements(elements), - self.epsilon_r, - self.basis.with_elements(elements), - self.E, + @Functional + def factor(w): + return np.sqrt(w["epsilon"]) * ( + dot(np.conj(w["E"][0]), w["E"][0]) + np.conj(w["E"][1]) * w["E"][1] + ) + + basis = self.basis.with_elements(elements) + basis_epsilon_r = self.basis_epsilon_r.with_elements(elements) + return ( + speed_of_light + * epsilon_0 + * factor.assemble( + basis, + E=basis.interpolate(self.E), + epsilon=basis_epsilon_r.interpolate(self.epsilon_r), + ) ) def calculate_pertubated_neff(self, delta_epsilon): @@ -392,25 +403,6 @@ def overlap(w): ) -def confinement_factor(basis_epsilon, epsilon, basis, E): - @Functional - def factor(w): - return ( - ( - np.sqrt(w["epsilon"]) - * (dot(np.conj(w["E"][0]), w["E"][0]) + np.conj(w["E"][1]) * w["E"][1]) - ) - * speed_of_light - * epsilon_0 - ) - - return factor.assemble( - basis, - E=basis.interpolate(E), - epsilon=basis_epsilon.interpolate(epsilon), - ) - - def plot_mode(basis, mode, plot_vectors=False, colorbar=True, title="E", direction="y"): from mpl_toolkits.axes_grid1 import make_axes_locatable From e675a8d1ec60c9d488bf89c19a4a02ba2d7cab86 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Sun, 9 Apr 2023 22:06:37 -0700 Subject: [PATCH 44/46] Revert "remove file which is already an example" This reverts commit 26a4c3882b7ce5f9b33706a0a99b976af6147514. --- femwell/thermal_transient.py | 288 +++++++++++++++++++++++++++++++++++ 1 file changed, 288 insertions(+) create mode 100644 femwell/thermal_transient.py diff --git a/femwell/thermal_transient.py b/femwell/thermal_transient.py new file mode 100644 index 00000000..87fd539c --- /dev/null +++ b/femwell/thermal_transient.py @@ -0,0 +1,288 @@ +from typing import Dict + +import numpy as np +from scipy.sparse.linalg import splu +from skfem import Basis, BilinearForm, ElementTriP0, LinearForm, Mesh, asm, enforce +from skfem.helpers import dot + +from femwell.thermal import solve_thermal + +""" +implemented like in https://www-users.cse.umn.edu/~arnold/8445.f11/notes.pdf page 81 in the middle of the page +""" + + +def solve_thermal_transient( + basis0, + thermal_conductivity_p0, + thermal_diffusivity_p0, + specific_conductivity: Dict[str, float], + current_densities_0, + current_densities, + fixed_boundaries, + dt, + steps, +): + basis, temperature = solve_thermal( + basis0, + thermal_conductivity_p0, + specific_conductivity, + {domain: current for domain, current in current_densities_0.items()}, + fixed_boundaries=fixed_boundaries, + ) + + @BilinearForm + def diffusivity_laplace(u, v, w): + return w["thermal_conductivity"] * dot(u.grad, v.grad) + + @BilinearForm + def mass(u, v, w): + return w["thermal_conductivity"] / w["thermal_diffusivity"] * u * v + + L = diffusivity_laplace.assemble( + basis, + thermal_conductivity=basis0.interpolate(thermal_conductivity_p0), + ) + M = mass.assemble( + basis, + thermal_diffusivity=basis0.interpolate(thermal_diffusivity_p0), + thermal_conductivity=basis0.interpolate(thermal_conductivity_p0), + ) + + theta = 0.5 # Crank–Nicolson + + x = basis.zeros() + for key, value in fixed_boundaries.items(): + x[basis.get_dofs(key)] = value + M0, L0 = enforce(M, L, D=basis.get_dofs(set(fixed_boundaries.keys()))) + A = M0 + theta * L0 * dt + B = M0 - (1 - theta) * L0 * dt + + backsolve = splu(A).solve # .T as splu prefers CSC + + t = 0 + temperatures = [temperature] + for i in range(steps): + joule_heating_rhs = basis.zeros() + for ( + domain, + current_density, + ) in current_densities.items(): # sum up the sources for the heating + current_density_p0 = basis0.zeros() + current_density_p0[basis0.get_dofs(elements=domain)] = current_density(t) + + @LinearForm + def joule_heating(v, w): + return w["current_density"] ** 2 / specific_conductivity[domain] * v + + joule_heating_rhs += asm( + joule_heating, + basis, + current_density=basis0.interpolate(current_density_p0), + ) + + t, temperature = t + dt, backsolve(B @ temperature + joule_heating_rhs * dt) + temperatures.append(temperature) + + return basis, temperatures + + +if __name__ == "__main__": + from collections import OrderedDict + + import matplotlib.pyplot as plt + import scipy.constants + from shapely.geometry import LineString, Polygon + from skfem.io import from_meshio + + from femwell.mesh import mesh_from_OrderedDict + + # Simulating the TiN TOPS heater in https://doi.org/10.1364/OE.27.010456 + + w_sim = 8 * 2 + h_clad = 2.8 + h_box = 1 + w_core = 0.5 + h_core = 0.22 + offset_heater = 2.2 + h_heater = 0.14 + w_heater = 2 + h_silicon = 3 + + wavelength = 1.55 + + polygons = OrderedDict( + bottom=LineString([(-w_sim / 2, -h_box), (w_sim / 2, -h_box)]), + core=Polygon( + [ + (-w_core / 2, 0), + (-w_core / 2, h_core), + (w_core / 2, h_core), + (w_core / 2, 0), + ] + ), + heater=Polygon( + [ + (-w_heater / 2, offset_heater), + (-w_heater / 2, offset_heater + h_heater), + (w_heater / 2, offset_heater + h_heater), + (w_heater / 2, offset_heater), + ] + ), + 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), + ] + ), + # silicon=Polygon([ + # (-w_sim / 2, - h_box - h_silicon), + # (-w_sim / 2, - h_box), + # (w_sim / 2, - h_box), + # (w_sim / 2, - h_box - h_silicon), + # ]), + ) + + resolutions = dict( + core={"resolution": 0.05, "distance": 1}, + clad={"resolution": 1, "distance": 1}, + box={"resolution": 1, "distance": 1}, + silicon={"resolution": 1, "distance": 1}, + heater={"resolution": 0.05, "distance": 1}, + ) + + mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=0.3)) + + basis0 = Basis(mesh, ElementTriP0(), intorder=4) + thermal_conductivity_p0 = basis0.zeros() + for domain, value in { + "core": 148, + "box": 1.38, + "clad": 1.38, + "heater": 28, + }.items(): # , 'silicon': 28 + thermal_conductivity_p0[basis0.get_dofs(elements=domain)] = value + thermal_conductivity_p0 *= 1e-12 # 1e-12 -> conversion from 1/m^2 -> 1/um^2 + + thermal_diffusivity_p0 = basis0.zeros() + for domain, value in { + "heater": 28 / 598 / 5240, + "box": 1.38 / 709 / 2203, + "clad": 1.38 / 709 / 2203, + "core": 148 / 711 / 2330, + # "silicon": 148 / 711 / 2330, + }.items(): + thermal_diffusivity_p0[basis0.get_dofs(elements=domain)] = value + thermal_diffusivity_p0 *= 1e12 # 1e-12 -> conversion from m^2 -> um^2 + + dt = 0.1e-5 + steps = 100 + current = ( + lambda t: 0.007 / polygons["heater"].area * ((t < dt * steps / 10) + (t > dt * steps / 2)) + ) + basis, temperatures = solve_thermal_transient( + basis0, + thermal_conductivity_p0, + thermal_diffusivity_p0, + specific_conductivity={"heater": 2.3e6}, + current_densities_0={"heater": current(0)}, + current_densities={"heater": current}, + fixed_boundaries={"bottom": 0}, + dt=dt, + steps=steps, + ) + + @LinearForm + def unit_load(v, w): + return v + + M = unit_load.assemble(basis) + + times = np.array([dt * i for i in range(steps + 1)]) + plt.xlabel("Time [us]") + plt.ylabel("Average temperature") + plt.plot(times * 1e6, M @ np.array(temperatures).T / np.sum(M)) + plt.show() + + # for i in range(0, steps, 10): + # fig, ax = plt.subplots(subplot_kw=dict(aspect=1)) + # for subdomain in mesh.subdomains.keys() - {'gmsh:bounding_entities'}: + # mesh.restrict(subdomain).draw(ax=ax, boundaries_only=True) + # basis.plot(temperatures[i], ax=ax, vmin=0, vmax=np.max(temperatures), shading='gouraud').show() + + # Calculate modes + + from tqdm.auto import tqdm + + from femwell.mode_solver import ( + calculate_coupling_coefficient, + compute_modes, + plot_mode, + ) + + epsilon_0 = basis0.zeros() + 1.444**2 + epsilon_0[basis0.get_dofs(elements="core")] = 3.4777**2 + lams_0, basis_modes_0, xs_0 = compute_modes( + basis0, epsilon_0, wavelength=wavelength, mu_r=1, num_modes=1 + ) + + neffs = [] + neffs_approximated = [] + for i, temperature in enumerate(tqdm(temperatures)): + # basis.plot(temperature, vmin=0, vmax=np.max(temperatures)) + # plt.show() + + temperature0 = basis0.project(basis.interpolate(temperature)) + epsilon = basis0.zeros() + (1.444 + 1.00e-5 * temperature0) ** 2 + epsilon[basis0.get_dofs(elements="core")] = ( + 3.4777 + 1.86e-4 * temperature0[basis0.get_dofs(elements="core")] + ) ** 2 + # basis0.plot(epsilon, colorbar=True).show() + + lams, basis_modes, xs = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=1) + + # plot_mode(basis_modes, xs[0]) + # plt.show() + + neffs.append(np.real(lams[0])) + neffs_approximated.append( + np.real( + lams_0[0] + + calculate_coupling_coefficient( + basis0, + (epsilon - epsilon_0) * scipy.constants.epsilon_0, + basis_modes_0, + xs_0[0], + xs_0[0], + ) + * scipy.constants.speed_of_light + * 2e-3 + ) + ) + + fig = plt.figure() + ax = fig.add_subplot(111) + ax.set_xlabel("Time [us]") + ax.set_ylabel("Current [mA]") + ax.plot(times * 1e6, current(times), "b-o") + ax2 = ax.twinx() + ax2.set_ylabel("Phase shift") + ax2.plot(times * 1e6, 2 * np.pi / 1.55 * (neffs - lams_0[0]) * 320, "r-o", label="Exact") + ax2.plot( + times * 1e6, + 2 * np.pi / 1.55 * (neffs_approximated - lams_0[0]) * 320, + "g-o", + label="Approximation", + ) + plt.legend() + plt.show() From 231451730a5f838f4eac7e73a653a5121d41f01f Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Sun, 9 Apr 2023 22:09:23 -0700 Subject: [PATCH 45/46] move functions in class --- femwell/maxwell/waveguide.py | 37 +++++++++++++----------------------- 1 file changed, 13 insertions(+), 24 deletions(-) diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index 0c5b3491..1fcf176b 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -96,8 +96,17 @@ def calculate_overlap(self, mode, elements=None): return calculate_overlap(self.basis, self.E, self.H, mode.basis, mode.E, mode.H) def calculate_coupling_coefficient(self, mode, delta_epsilon): - return calculate_coupling_coefficient( - self.basis_epsilon_r, delta_epsilon, self.basis, self.E, mode.E + @Functional(dtype=complex) + def overlap(w): + return w["delta_epsilon"] * ( + dot(np.conj(w["E_i"][0]), w["E_j"][0]) + np.conj(w["E_i"][1]) * w["E_j"][1] + ) + + return overlap.assemble( + self.basis, + E_i=self.basis.interpolate(self.E), + E_j=self.basis.interpolate(mode.E), + delta_epsilon=self.basis_epsilon_r.interpolate(delta_epsilon), ) def calculate_propagation_loss(self, distance): @@ -132,13 +141,8 @@ def factor(w): def calculate_pertubated_neff(self, delta_epsilon): return ( self.n_eff - + calculate_coupling_coefficient( - self.basis_epsilon_r, - delta_epsilon * scipy.constants.epsilon_0, - self.basis, - self.E, - self.E, - ) + + self.calculate_coupling_coefficient(self, delta_epsilon) + * scipy.constants.epsilon_0 * scipy.constants.speed_of_light * 0.5 ) @@ -388,21 +392,6 @@ def overlap(w): return overlap.assemble(basis_i, E_i=basis_i.interpolate(E_i)) -def calculate_coupling_coefficient(basis_epsilon, delta_epsilon, basis, E_i, E_j): - @Functional(dtype=complex) - def overlap(w): - return w["delta_epsilon"] * ( - dot(np.conj(w["E_i"][0]), w["E_j"][0]) + np.conj(w["E_i"][1]) * w["E_j"][1] - ) - - return overlap.assemble( - basis, - E_i=basis.interpolate(E_i), - E_j=basis.interpolate(E_j), - delta_epsilon=basis_epsilon.interpolate(delta_epsilon), - ) - - def plot_mode(basis, mode, plot_vectors=False, colorbar=True, title="E", direction="y"): from mpl_toolkits.axes_grid1 import make_axes_locatable From eed92f72a715d8353059ce96b24f3c9e0d30b9ef Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Sun, 9 Apr 2023 22:14:59 -0700 Subject: [PATCH 46/46] remove wrong default value --- femwell/maxwell/waveguide.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index 1fcf176b..89edc0b4 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -289,7 +289,7 @@ def bform(e_t, e_z, v_t, v_z, w): ) -def calculate_hfield(basis, xs, beta, omega=1): +def calculate_hfield(basis, xs, beta, omega): @BilinearForm(dtype=np.complex64) def aform(e_t, e_z, v_t, v_z, w): return (