From 8f6f7379323d49b9ce540ab8b3715a207009546a Mon Sep 17 00:00:00 2001 From: Kunyuan-LI Date: Sat, 24 Feb 2024 16:43:59 +0100 Subject: [PATCH 1/8] SFWM Aeff calculation --- femwell/examples/Aeff.py | 94 ++++++++++++++++++++ femwell/maxwell/waveguide.py | 160 +++++++++++++---------------------- 2 files changed, 152 insertions(+), 102 deletions(-) create mode 100644 femwell/examples/Aeff.py diff --git a/femwell/examples/Aeff.py b/femwell/examples/Aeff.py new file mode 100644 index 00000000..573d0607 --- /dev/null +++ b/femwell/examples/Aeff.py @@ -0,0 +1,94 @@ +from collections import OrderedDict +import matplotlib.pyplot as plt +import numpy as np +from shapely.geometry import box +from shapely.ops import clip_by_rect +from skfem import Basis, ElementTriP0 +from skfem.io.meshio import from_meshio +from femwell.maxwell.waveguide import compute_modes +from femwell.maxwell.waveguide import calculate_sfwm_Aeff +from femwell.mesh import mesh_from_OrderedDict +from scipy.constants import c, epsilon_0 + + +#Dispersion relations of materials +#Core +def n_X(wavelength): + x = wavelength + return (1 + + 2.19244563 / (1 - (0.20746607 / x) ** 2) + + 0.13435116 / (1 - (0.3985835 / x) ** 2) + + 2.20997784 / (1 - (0.20747044 / x) ** 2) + ) ** 0.5 + +# Box +def n_silicon_dioxide(wavelength): + x = wavelength + return ( + 1 + + 0.6961663 / (1 - (0.0684043 / x) ** 2) + + 0.4079426 / (1 - (0.1162414 / x) ** 2) + + 0.8974794 / (1 - (9.896161 / x) ** 2) + ) ** 0.5 + +Clad=1 + + +core = box(0, 0, 0.5, 0.39) #500x390nm +polygons = OrderedDict( + core=core, + box=clip_by_rect(core.buffer(1.5, resolution=4), -np.inf, -np.inf, np.inf, 0), + clad=clip_by_rect(core.buffer(1.5, resolution=4), -np.inf, 0, np.inf, np.inf), +) + +resolutions = {"core": {"resolution": 0.025, "distance": 2.}} + +mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=0.6)) + +num_modes = 2 + +#For SFWM we have energy conservation and momemtum(k) conservation for 2pumps and signal+idler +lambda_p0 = 1.4 +lambda_s0 = 1.097 +lambda_i0 = 1.686 + +basis0 = Basis(mesh, ElementTriP0()) + +epsilon_p = basis0.zeros() +epsilon_s = basis0.zeros() +epsilon_i = basis0.zeros() + + +for wavelength, epsilon in zip([lambda_p0, lambda_s0, lambda_i0], [epsilon_p, epsilon_s, epsilon_i]): + for subdomain, n_func in { + "core": n_X, + "box": n_silicon_dioxide, + "clad": lambda _: Clad, + }.items(): + n = n_func(wavelength) + epsilon[basis0.get_dofs(elements=subdomain)] = n**2 + + +modes_p = compute_modes(basis0, epsilon_p, wavelength=lambda_p0, num_modes=num_modes, order=1) +modes_s = compute_modes(basis0, epsilon_s, wavelength=lambda_s0, num_modes=num_modes, order=1) +modes_i = compute_modes(basis0, epsilon_i, wavelength=lambda_i0, num_modes=num_modes, order=1) + + +mode_p = max(modes_p, key=lambda mode: mode.te_fraction) +mode_s = max(modes_s, key=lambda mode: mode.te_fraction) +mode_i = max(modes_i, key=lambda mode: mode.te_fraction) + +A_eff = calculate_sfwm_Aeff(basis0, mode_p, mode_s, mode_i) +print(A_eff) + +#Calculation for non-linear coef +chi_3 = 5e-21 # m^2/V^2 #7e-20? +lambda_p0_m = lambda_p0 * 1e-6 #to m +n_p0 = np.real(mode_p.n_eff) +A_eff_m2 = A_eff * 1e-12 #to m^2 + +omega_p0 = 2 * np.pi * c / lambda_p0_m + +gamma = (3 * chi_3 * omega_p0) / (4 * epsilon_0 * c**2 * n_p0**2 * A_eff_m2) + +print("gamma:",gamma) diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index 204b2d74..dd217d3c 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -6,12 +6,11 @@ import matplotlib.pyplot as plt import numpy as np -import scipy.constants import scipy.sparse.linalg from matplotlib.axes import Axes from matplotlib.figure import Figure from numpy.typing import NDArray -from scipy.constants import epsilon_0, speed_of_light +from scipy.constants import epsilon_0, speed_of_light, mu_0 from skfem import ( Basis, BilinearForm, @@ -33,7 +32,7 @@ from skfem.utils import solver_eigen_scipy -@dataclass(frozen=True) +@dataclass(frozen=False)#True class Mode: frequency: float """Frequency of the light""" @@ -70,44 +69,6 @@ def n_eff(self): """Effective refractive index of the mode""" return self.k / self.k0 - @cached_property - def poynting(self): - """Poynting vector of the mode""" - - # Extraction of the fields - (Ex, Ey), Ez = self.basis.interpolate(self.E) - (Hx, Hy), Hz = self.basis.interpolate(self.H) - - # New basis will be the discontinuous variant used by the solved Ez-field - poynting_basis = self.basis.with_element( - ElementVector(ElementDG(self.basis.elem.elems[1]), 3) - ) - - # Calculation of the Poynting vector - Px = Ey * Hz - Ez * Hy - Py = Ez * Hx - Ex * Hz - Pz = Ex * Hy - Ey * Hx - - # Projection of the Poynting vector on the new basis - P_proj = poynting_basis.project(np.stack([Px, Py, Pz], axis=0), dtype=np.complex64) - - return poynting_basis, P_proj - - @property - def Sx(self): - basis, _P = self.poynting - return basis.split_bases()[0], _P[basis.split_indices()[0]] - - @property - def Sy(self): - basis, _P = self.poynting - return basis.split_bases()[1], _P[basis.split_indices()[1]] - - @property - def Sz(self): - basis, _P = self.poynting - return basis.split_bases()[2], _P[basis.split_indices()[2]] - @cached_property def te_fraction(self): """TE-fraction of the mode""" @@ -131,32 +92,6 @@ def tm_fraction(self): return 1 - self.te_fraction - @cached_property - def transversality(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 - - @Functional - def ez(w): - return np.abs( - w.E[0][0] * np.conj(w.E[0][0]) - + w.E[0][1] * np.conj(w.E[0][1]) - + w.E[1] * np.conj(w.E[1]) - ) - - ex_sum = ex.assemble(self.basis, E=self.basis.interpolate(self.E)) - ey_sum = ey.assemble(self.basis, E=self.basis.interpolate(self.E)) - ez_sum = ez.assemble(self.basis, E=self.basis.interpolate(self.E)) - - return (ex_sum + ey_sum) / (ez_sum) - def __repr__(self) -> str: return f"{self.__class__.__name__}(k: {self.k}, n_eff:{self.n_eff})" @@ -177,40 +112,6 @@ def overlap(w): delta_epsilon=self.basis_epsilon_r.interpolate(delta_epsilon), ) - def calculate_effective_area(self, field="xy"): - if field == "xy": - - @Functional(dtype=complex) - def E2(w): - return dot(np.conj(w["E"][0]), w["E"][0]) - - @Functional(dtype=complex) - def E4(w): - return dot(np.conj(w["E"][0]), w["E"][0]) ** 2 - - else: - num = 0 if (field == "x") else 1 - - @Functional(dtype=complex) - def E2(w): - return np.conj(w["E"][0][num]) * w["E"][0][num] - - @Functional(dtype=complex) - def E4(w): - return (np.conj(w["E"][0][num]) * w["E"][0][num]) ** 2 - - return np.real( - E2.assemble( - self.basis, - E=self.basis.interpolate(self.E), - ) - ** 2 - / E4.assemble( - self.basis, - E=self.basis.interpolate(self.E), - ) - ) - def calculate_propagation_loss(self, distance): return -20 / np.log(10) * self.k0 * np.imag(self.n_eff) * distance @@ -320,7 +221,7 @@ def plot_intensity( return fig, ax -@dataclass(frozen=True) +@dataclass(frozen=False)#True class Modes: modes: List @@ -655,6 +556,61 @@ def edge_jump(w): return eta_K + eta_E +def calculate_sfwm_Aeff( + basis: Basis, + mode_p, + mode_s, + mode_i, +) -> np.complex64: + """ + Calculates the overlap integral for SFWM process by considering the interactions + between pump, signal, and idler modes in the xy plane. + + Args: + basis (Basis): Common basis used for all modes in the same geometric structure. + mode_p, mode_s, mode_i: Mode instances for pump, signal, and idler, respectively. + + Returns: + np.complex64: The Aeff result for the SFWM process(1/overlap integral). + """ + def normalize_mode(mode): + @Functional + def E2(w): + return dot(w["E"][0], np.conj(w["E"][0])) + #return dot(w["E"][0], np.conj(w["E"][0])) + #return (w["E"][0]* np.conj(w["E"][0])) + E_xy, _ = mode.basis.interpolate(mode.E) #[0]=xy [1]=z + E_squared_integral = E2.assemble(mode.basis, E=E_xy) + normalization_factor = 1 / np.sqrt(E_squared_integral) + mode.E0 = mode.E * normalization_factor + + normalize_mode(mode_p) + normalize_mode(mode_s) + normalize_mode(mode_i) + + #Normalization needed for uv(x,y)! + E_p, _ = mode_p.basis.interpolate(mode_p.E0) #xy, z + E_s, _ = mode_s.basis.interpolate(mode_s.E0) + E_i, _ = mode_i.basis.interpolate(mode_i.E0) + + @Functional(dtype=np.complex64) + def sfwm_overlap(w): + overlap_SFWM = w["E_p"][0] * w["E_p"][0] * np.conj(w["E_s"][0]) * np.conj(w["E_i"][0]) #up up us* ui*? + return (overlap_SFWM) + + # Assemble the integral over the basis to compute the overlap + overlap_result = sfwm_overlap.assemble( + basis, + E_p=E_p, + E_s=E_s, + E_i=E_i, + ) + + return 1/overlap_result + + + + if __name__ == "__main__": from collections import OrderedDict From e52204c574c44414ab07045289f6e1b16eeb7fff Mon Sep 17 00:00:00 2001 From: Kunyuan-LI Date: Sat, 2 Mar 2024 12:26:31 +0100 Subject: [PATCH 2/8] Aeff.py as the standalone example --- femwell/examples/Aeff.py | 610 ++++++++++++++++++++++++++++++++++- femwell/maxwell/waveguide.py | 160 +++++---- 2 files changed, 710 insertions(+), 60 deletions(-) diff --git a/femwell/examples/Aeff.py b/femwell/examples/Aeff.py index 573d0607..ab3406a4 100644 --- a/femwell/examples/Aeff.py +++ b/femwell/examples/Aeff.py @@ -1,3 +1,7 @@ +"""Aeff analysis based on https://optics.ansys.com/hc/en-us/articles/15100783091731-Spontaneous-Four-wave-Mixing-SWFM-Microring-Resonator-Photon-Source.""" +from dataclasses import dataclass +from functools import cached_property +from typing import List, Tuple from collections import OrderedDict import matplotlib.pyplot as plt import numpy as np @@ -5,11 +9,613 @@ from shapely.ops import clip_by_rect from skfem import Basis, ElementTriP0 from skfem.io.meshio import from_meshio -from femwell.maxwell.waveguide import compute_modes -from femwell.maxwell.waveguide import calculate_sfwm_Aeff from femwell.mesh import mesh_from_OrderedDict from scipy.constants import c, epsilon_0 +import matplotlib.pyplot as plt +import numpy as np +import scipy.sparse.linalg +from matplotlib.axes import Axes +from matplotlib.figure import Figure +from numpy.typing import NDArray +from scipy.constants import epsilon_0, speed_of_light, mu_0 +from skfem import ( + Basis, + BilinearForm, + ElementDG, + ElementTriN1, + ElementTriN2, + ElementTriP0, + ElementTriP1, + ElementTriP2, + ElementVector, + Functional, + InteriorFacetBasis, + LinearForm, + Mesh, + condense, + solve, +) +from skfem.helpers import cross, curl, dot, grad, inner +from skfem.utils import solver_eigen_scipy + + +@dataclass(frozen=False)#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, 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): + @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): + 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 calculate_confinement_factor(self, elements): + @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): + return ( + self.n_eff + + self.calculate_coupling_coefficient(self, delta_epsilon) + * scipy.constants.epsilon_0 + * scipy.constants.speed_of_light + * 0.5 + ) + + def calculate_intensity(self) -> Tuple[NDArray, Basis]: + """Calculates the intensity of a mode. + + The intensity is calculated from the cross-product between the electric and magnetic field, as + described in https://doi.org/10.1364/OE.16.016659. + + The calculation is performed as follows: + 1) The electric and magnetic fields are interpolated on the quadrature points with the simulation basis; + 2) The intensity is calculated directly on the quadrature points; + 3) The intensity is projected on a new discontinuous, piecewise linear triangular basis. + + Returns: + Basis, NDArray: Plot-ready basis and intensity array + """ + (Ex, Ey), _ = self.basis.interpolate(self.E) + (Hx, Hy), _ = self.basis.interpolate(np.conj(self.H)) + intensity = 0.5 * np.real(Ex * Hy - Ey * Hx) + basis2 = self.basis.with_element(ElementDG(ElementTriP1())) + intensity2 = basis2.project(intensity) + + return basis2, intensity2 + + 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 plot_intensity( + self, + ax: Axes = None, + colorbar: bool = True, + normalize: bool = True, + ) -> Tuple[Figure, Axes]: + """Plots the intensity of a mode as outlined in `calculate_intensity`. + + Args: + ax (Axes, optional): Axes onto which the plot is drawn. Defaults to None. + colorbar (bool, optional): Adds a colorbar to the plot. Defaults to True. + normalize (bool, optional): Normalizes the intensity by its maximum value. Defaults to True. + + Returns: + Tuple[Figure, Axes]: Figure and axes of the plot. + """ + intensity_basis, intensity = self.calculate_intensity() + if normalize: + intensity = intensity / intensity.max() + + if ax is None: + fig, ax = plt.subplots() + else: + fig = plt.gcf() + + for subdomain in self.basis.mesh.subdomains.keys() - {"gmsh:bounding_entities"}: + self.basis.mesh.restrict(subdomain).draw(ax=ax, boundaries_only=True, color="w") + intensity_basis.plot(intensity, ax=ax, cmap="inferno") + + if colorbar: + plt.colorbar(ax.collections[-1]) + + return fig, ax + + +@dataclass(frozen=False)#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)) + + @property + def n_effs(self): + return np.array([mode.n_eff for mode in self.modes]) + + +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="scipy", +) -> Modes: + 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 + + @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) * 1.1 + + if metallic_boundaries: + lams, xs = solve( + *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: + lams, xs = solve( + -A, + -B, + solver=solver(k=num_modes, sigma=sigma), + ) + + 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) + ] + ) + + +def calculate_hfield(basis, xs, beta, omega): + @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: Basis, + E_i: np.ndarray, + H_i: np.ndarray, + basis_j: Basis, + E_j: np.ndarray, + H_j: np.ndarray, +) -> np.complex64: + """Calculates the fully vectorial overlap between two modes. + + If the modes do not share the basis, interpolation is performed automatically. + + Args: + basis_i (Basis): Basis of the first mode + E_i (np.ndarray): Electric field of the first mode + H_i (np.ndarray): Magnetic field of the first mode + basis_j (Basis): Basis of the second mode + E_j (np.ndarray): Electric field of the first mode + H_j (np.ndarray): Magnetic field of the first mode + + Returns: + np.complex64: Complex overlap between the two modes + """ + + @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 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_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 eval_error_estimator(basis, u): + @Functional + def interior_residual(w): + h = w.h + x, y = w.x + return h**2 # * load_func(x, y) ** 2 + + eta_K = interior_residual.elemental(basis, w=basis.interpolate(u)) + + # facet jump + fbasis = [InteriorFacetBasis(basis.mesh, basis.elem, side=i) for i in [0, 1]] + w = {"u" + str(i + 1): fbasis[i].interpolate(u) for i in [0, 1]} + + @Functional + def edge_jump(w): + return w.h * ( + np.abs(dot(grad(w["u1"][1]) - grad(w["u2"][1]), w.n)) ** 2 + + np.abs(dot(w["u1"][0] - w["u2"][0], w.n)) ** 2 + ) + + tmp = np.zeros(basis.mesh.facets.shape[1]) + tmp[fbasis[0].find] = edge_jump.elemental(fbasis[0], **w) + eta_E = np.sum(0.5 * tmp[basis.mesh.t2f], axis=0) + + return eta_K + eta_E + + +def calculate_sfwm_Aeff( + basis: Basis, + mode_p, + mode_s, + mode_i, +) -> np.complex64: + """ + Calculates the overlap integral for SFWM process by considering the interactions + between pump, signal, and idler modes in the xy plane. + + Args: + basis (Basis): Common basis used for all modes in the same geometric structure. + mode_p, mode_s, mode_i: Mode instances for pump, signal, and idler, respectively. + + Returns: + np.complex64: The Aeff result for the SFWM process(1/overlap integral). + """ + def normalize_mode(mode): + @Functional + def E2(w): + return dot(w["E"][0], np.conj(w["E"][0])) + #return dot(w["E"][0], np.conj(w["E"][0])) + #return (w["E"][0]* np.conj(w["E"][0])) + E_xy, _ = mode.basis.interpolate(mode.E) #[0]=xy [1]=z + E_squared_integral = E2.assemble(mode.basis, E=E_xy) + normalization_factor = 1 / np.sqrt(E_squared_integral) + mode.E0 = mode.E * normalization_factor + + normalize_mode(mode_p) + normalize_mode(mode_s) + normalize_mode(mode_i) + + #Normalization needed for uv(x,y)! + E_p, _ = mode_p.basis.interpolate(mode_p.E0) #xy, z + E_s, _ = mode_s.basis.interpolate(mode_s.E0) + E_i, _ = mode_i.basis.interpolate(mode_i.E0) + + @Functional(dtype=np.complex64) + def sfwm_overlap(w): + overlap_SFWM = w["E_p"][0] * w["E_p"][0] * np.conj(w["E_s"][0]) * np.conj(w["E_i"][0]) #up up us* ui*? + return (overlap_SFWM) + + # Assemble the integral over the basis to compute the overlap + overlap_result = sfwm_overlap.assemble( + basis, + E_p=E_p, + E_s=E_s, + E_i=E_i, + ) + + return 1/overlap_result + #Dispersion relations of materials #Core diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index dd217d3c..204b2d74 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -6,11 +6,12 @@ import matplotlib.pyplot as plt import numpy as np +import scipy.constants import scipy.sparse.linalg from matplotlib.axes import Axes from matplotlib.figure import Figure from numpy.typing import NDArray -from scipy.constants import epsilon_0, speed_of_light, mu_0 +from scipy.constants import epsilon_0, speed_of_light from skfem import ( Basis, BilinearForm, @@ -32,7 +33,7 @@ from skfem.utils import solver_eigen_scipy -@dataclass(frozen=False)#True +@dataclass(frozen=True) class Mode: frequency: float """Frequency of the light""" @@ -69,6 +70,44 @@ def n_eff(self): """Effective refractive index of the mode""" return self.k / self.k0 + @cached_property + def poynting(self): + """Poynting vector of the mode""" + + # Extraction of the fields + (Ex, Ey), Ez = self.basis.interpolate(self.E) + (Hx, Hy), Hz = self.basis.interpolate(self.H) + + # New basis will be the discontinuous variant used by the solved Ez-field + poynting_basis = self.basis.with_element( + ElementVector(ElementDG(self.basis.elem.elems[1]), 3) + ) + + # Calculation of the Poynting vector + Px = Ey * Hz - Ez * Hy + Py = Ez * Hx - Ex * Hz + Pz = Ex * Hy - Ey * Hx + + # Projection of the Poynting vector on the new basis + P_proj = poynting_basis.project(np.stack([Px, Py, Pz], axis=0), dtype=np.complex64) + + return poynting_basis, P_proj + + @property + def Sx(self): + basis, _P = self.poynting + return basis.split_bases()[0], _P[basis.split_indices()[0]] + + @property + def Sy(self): + basis, _P = self.poynting + return basis.split_bases()[1], _P[basis.split_indices()[1]] + + @property + def Sz(self): + basis, _P = self.poynting + return basis.split_bases()[2], _P[basis.split_indices()[2]] + @cached_property def te_fraction(self): """TE-fraction of the mode""" @@ -92,6 +131,32 @@ def tm_fraction(self): return 1 - self.te_fraction + @cached_property + def transversality(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 + + @Functional + def ez(w): + return np.abs( + w.E[0][0] * np.conj(w.E[0][0]) + + w.E[0][1] * np.conj(w.E[0][1]) + + w.E[1] * np.conj(w.E[1]) + ) + + ex_sum = ex.assemble(self.basis, E=self.basis.interpolate(self.E)) + ey_sum = ey.assemble(self.basis, E=self.basis.interpolate(self.E)) + ez_sum = ez.assemble(self.basis, E=self.basis.interpolate(self.E)) + + return (ex_sum + ey_sum) / (ez_sum) + def __repr__(self) -> str: return f"{self.__class__.__name__}(k: {self.k}, n_eff:{self.n_eff})" @@ -112,6 +177,40 @@ def overlap(w): delta_epsilon=self.basis_epsilon_r.interpolate(delta_epsilon), ) + def calculate_effective_area(self, field="xy"): + if field == "xy": + + @Functional(dtype=complex) + def E2(w): + return dot(np.conj(w["E"][0]), w["E"][0]) + + @Functional(dtype=complex) + def E4(w): + return dot(np.conj(w["E"][0]), w["E"][0]) ** 2 + + else: + num = 0 if (field == "x") else 1 + + @Functional(dtype=complex) + def E2(w): + return np.conj(w["E"][0][num]) * w["E"][0][num] + + @Functional(dtype=complex) + def E4(w): + return (np.conj(w["E"][0][num]) * w["E"][0][num]) ** 2 + + return np.real( + E2.assemble( + self.basis, + E=self.basis.interpolate(self.E), + ) + ** 2 + / E4.assemble( + self.basis, + E=self.basis.interpolate(self.E), + ) + ) + def calculate_propagation_loss(self, distance): return -20 / np.log(10) * self.k0 * np.imag(self.n_eff) * distance @@ -221,7 +320,7 @@ def plot_intensity( return fig, ax -@dataclass(frozen=False)#True +@dataclass(frozen=True) class Modes: modes: List @@ -556,61 +655,6 @@ def edge_jump(w): return eta_K + eta_E -def calculate_sfwm_Aeff( - basis: Basis, - mode_p, - mode_s, - mode_i, -) -> np.complex64: - """ - Calculates the overlap integral for SFWM process by considering the interactions - between pump, signal, and idler modes in the xy plane. - - Args: - basis (Basis): Common basis used for all modes in the same geometric structure. - mode_p, mode_s, mode_i: Mode instances for pump, signal, and idler, respectively. - - Returns: - np.complex64: The Aeff result for the SFWM process(1/overlap integral). - """ - def normalize_mode(mode): - @Functional - def E2(w): - return dot(w["E"][0], np.conj(w["E"][0])) - #return dot(w["E"][0], np.conj(w["E"][0])) - #return (w["E"][0]* np.conj(w["E"][0])) - E_xy, _ = mode.basis.interpolate(mode.E) #[0]=xy [1]=z - E_squared_integral = E2.assemble(mode.basis, E=E_xy) - normalization_factor = 1 / np.sqrt(E_squared_integral) - mode.E0 = mode.E * normalization_factor - - normalize_mode(mode_p) - normalize_mode(mode_s) - normalize_mode(mode_i) - - #Normalization needed for uv(x,y)! - E_p, _ = mode_p.basis.interpolate(mode_p.E0) #xy, z - E_s, _ = mode_s.basis.interpolate(mode_s.E0) - E_i, _ = mode_i.basis.interpolate(mode_i.E0) - - @Functional(dtype=np.complex64) - def sfwm_overlap(w): - overlap_SFWM = w["E_p"][0] * w["E_p"][0] * np.conj(w["E_s"][0]) * np.conj(w["E_i"][0]) #up up us* ui*? - return (overlap_SFWM) - - # Assemble the integral over the basis to compute the overlap - overlap_result = sfwm_overlap.assemble( - basis, - E_p=E_p, - E_s=E_s, - E_i=E_i, - ) - - return 1/overlap_result - - - - if __name__ == "__main__": from collections import OrderedDict From cbda3f16499a14de718988a9bc4ed4c00c82989b Mon Sep 17 00:00:00 2001 From: Kunyuan-LI Date: Sun, 3 Mar 2024 12:34:13 +0100 Subject: [PATCH 3/8] adjust waveguide.py --- femwell/examples/Aeff.py | 615 ++------------------------------------- 1 file changed, 29 insertions(+), 586 deletions(-) diff --git a/femwell/examples/Aeff.py b/femwell/examples/Aeff.py index ab3406a4..9d343aa1 100644 --- a/femwell/examples/Aeff.py +++ b/femwell/examples/Aeff.py @@ -1,575 +1,23 @@ """Aeff analysis based on https://optics.ansys.com/hc/en-us/articles/15100783091731-Spontaneous-Four-wave-Mixing-SWFM-Microring-Resonator-Photon-Source.""" -from dataclasses import dataclass -from functools import cached_property -from typing import List, Tuple from collections import OrderedDict -import matplotlib.pyplot as plt import numpy as np from shapely.geometry import box from shapely.ops import clip_by_rect from skfem import Basis, ElementTriP0 +from skfem import Functional +from skfem.helpers import dot from skfem.io.meshio import from_meshio +from femwell.maxwell.waveguide import compute_modes from femwell.mesh import mesh_from_OrderedDict from scipy.constants import c, epsilon_0 -import matplotlib.pyplot as plt -import numpy as np -import scipy.sparse.linalg -from matplotlib.axes import Axes -from matplotlib.figure import Figure -from numpy.typing import NDArray -from scipy.constants import epsilon_0, speed_of_light, mu_0 -from skfem import ( - Basis, - BilinearForm, - ElementDG, - ElementTriN1, - ElementTriN2, - ElementTriP0, - ElementTriP1, - ElementTriP2, - ElementVector, - Functional, - InteriorFacetBasis, - LinearForm, - Mesh, - condense, - solve, -) -from skfem.helpers import cross, curl, dot, grad, inner -from skfem.utils import solver_eigen_scipy - - -@dataclass(frozen=False)#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, 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): - @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): - 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 calculate_confinement_factor(self, elements): - @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): - return ( - self.n_eff - + self.calculate_coupling_coefficient(self, delta_epsilon) - * scipy.constants.epsilon_0 - * scipy.constants.speed_of_light - * 0.5 - ) - - def calculate_intensity(self) -> Tuple[NDArray, Basis]: - """Calculates the intensity of a mode. - - The intensity is calculated from the cross-product between the electric and magnetic field, as - described in https://doi.org/10.1364/OE.16.016659. - - The calculation is performed as follows: - 1) The electric and magnetic fields are interpolated on the quadrature points with the simulation basis; - 2) The intensity is calculated directly on the quadrature points; - 3) The intensity is projected on a new discontinuous, piecewise linear triangular basis. - - Returns: - Basis, NDArray: Plot-ready basis and intensity array - """ - (Ex, Ey), _ = self.basis.interpolate(self.E) - (Hx, Hy), _ = self.basis.interpolate(np.conj(self.H)) - intensity = 0.5 * np.real(Ex * Hy - Ey * Hx) - basis2 = self.basis.with_element(ElementDG(ElementTriP1())) - intensity2 = basis2.project(intensity) - - return basis2, intensity2 - - 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 plot_intensity( - self, - ax: Axes = None, - colorbar: bool = True, - normalize: bool = True, - ) -> Tuple[Figure, Axes]: - """Plots the intensity of a mode as outlined in `calculate_intensity`. - - Args: - ax (Axes, optional): Axes onto which the plot is drawn. Defaults to None. - colorbar (bool, optional): Adds a colorbar to the plot. Defaults to True. - normalize (bool, optional): Normalizes the intensity by its maximum value. Defaults to True. - - Returns: - Tuple[Figure, Axes]: Figure and axes of the plot. - """ - intensity_basis, intensity = self.calculate_intensity() - if normalize: - intensity = intensity / intensity.max() - - if ax is None: - fig, ax = plt.subplots() - else: - fig = plt.gcf() - - for subdomain in self.basis.mesh.subdomains.keys() - {"gmsh:bounding_entities"}: - self.basis.mesh.restrict(subdomain).draw(ax=ax, boundaries_only=True, color="w") - intensity_basis.plot(intensity, ax=ax, cmap="inferno") - - if colorbar: - plt.colorbar(ax.collections[-1]) - - return fig, ax - - -@dataclass(frozen=False)#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)) - - @property - def n_effs(self): - return np.array([mode.n_eff for mode in self.modes]) - - -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="scipy", -) -> Modes: - 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 - - @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) * 1.1 - - if metallic_boundaries: - lams, xs = solve( - *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: - lams, xs = solve( - -A, - -B, - solver=solver(k=num_modes, sigma=sigma), - ) - - 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) - ] - ) - - -def calculate_hfield(basis, xs, beta, omega): - @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: Basis, - E_i: np.ndarray, - H_i: np.ndarray, - basis_j: Basis, - E_j: np.ndarray, - H_j: np.ndarray, -) -> np.complex64: - """Calculates the fully vectorial overlap between two modes. - - If the modes do not share the basis, interpolation is performed automatically. - - Args: - basis_i (Basis): Basis of the first mode - E_i (np.ndarray): Electric field of the first mode - H_i (np.ndarray): Magnetic field of the first mode - basis_j (Basis): Basis of the second mode - E_j (np.ndarray): Electric field of the first mode - H_j (np.ndarray): Magnetic field of the first mode - - Returns: - np.complex64: Complex overlap between the two modes - """ - - @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 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_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 eval_error_estimator(basis, u): - @Functional - def interior_residual(w): - h = w.h - x, y = w.x - return h**2 # * load_func(x, y) ** 2 - - eta_K = interior_residual.elemental(basis, w=basis.interpolate(u)) - - # facet jump - fbasis = [InteriorFacetBasis(basis.mesh, basis.elem, side=i) for i in [0, 1]] - w = {"u" + str(i + 1): fbasis[i].interpolate(u) for i in [0, 1]} - - @Functional - def edge_jump(w): - return w.h * ( - np.abs(dot(grad(w["u1"][1]) - grad(w["u2"][1]), w.n)) ** 2 - + np.abs(dot(w["u1"][0] - w["u2"][0], w.n)) ** 2 - ) - - tmp = np.zeros(basis.mesh.facets.shape[1]) - tmp[fbasis[0].find] = edge_jump.elemental(fbasis[0], **w) - eta_E = np.sum(0.5 * tmp[basis.mesh.t2f], axis=0) - - return eta_K + eta_E - - def calculate_sfwm_Aeff( - basis: Basis, - mode_p, - mode_s, - mode_i, + basis: Basis, + mode_p, + mode_s, + mode_i ) -> np.complex64: + """ Calculates the overlap integral for SFWM process by considering the interactions between pump, signal, and idler modes in the xy plane. @@ -581,40 +29,35 @@ def calculate_sfwm_Aeff( Returns: np.complex64: The Aeff result for the SFWM process(1/overlap integral). """ + def normalize_mode(mode): @Functional def E2(w): return dot(w["E"][0], np.conj(w["E"][0])) - #return dot(w["E"][0], np.conj(w["E"][0])) - #return (w["E"][0]* np.conj(w["E"][0])) - E_xy, _ = mode.basis.interpolate(mode.E) #[0]=xy [1]=z + + E_xy, _ = mode.basis.interpolate(mode.E) # [0]=xy [1]=z E_squared_integral = E2.assemble(mode.basis, E=E_xy) normalization_factor = 1 / np.sqrt(E_squared_integral) - mode.E0 = mode.E * normalization_factor - - normalize_mode(mode_p) - normalize_mode(mode_s) - normalize_mode(mode_i) - - #Normalization needed for uv(x,y)! - E_p, _ = mode_p.basis.interpolate(mode_p.E0) #xy, z - E_s, _ = mode_s.basis.interpolate(mode_s.E0) - E_i, _ = mode_i.basis.interpolate(mode_i.E0) - + return normalization_factor # Return the normalization factor instead of modifying the mode + + # Calculate normalization factors for each mode + norm_factor_p = normalize_mode(mode_p) + norm_factor_s = normalize_mode(mode_s) + norm_factor_i = normalize_mode(mode_i) + + # Apply normalization factors to the electric fields for the overlap calculation + E_p, _ = mode_p.basis.interpolate(mode_p.E * norm_factor_p) + E_s, _ = mode_s.basis.interpolate(mode_s.E * norm_factor_s) + E_i, _ = mode_i.basis.interpolate(mode_i.E * norm_factor_i) + @Functional(dtype=np.complex64) def sfwm_overlap(w): - overlap_SFWM = w["E_p"][0] * w["E_p"][0] * np.conj(w["E_s"][0]) * np.conj(w["E_i"][0]) #up up us* ui*? - return (overlap_SFWM) + overlap_SFWM = w["E_p"][0] * w["E_p"][0] * np.conj(w["E_s"][0]) * np.conj(w["E_i"][0]) + return overlap_SFWM - # Assemble the integral over the basis to compute the overlap - overlap_result = sfwm_overlap.assemble( - basis, - E_p=E_p, - E_s=E_s, - E_i=E_i, - ) + overlap_result = sfwm_overlap.assemble(basis, E_p=E_p, E_s=E_s, E_i=E_i) - return 1/overlap_result + return 1 / overlap_result #Dispersion relations of materials @@ -685,7 +128,7 @@ def n_silicon_dioxide(wavelength): mode_i = max(modes_i, key=lambda mode: mode.te_fraction) A_eff = calculate_sfwm_Aeff(basis0, mode_p, mode_s, mode_i) -print(A_eff) +print("Aeff in um2:", A_eff) #Calculation for non-linear coef chi_3 = 5e-21 # m^2/V^2 #7e-20? @@ -697,4 +140,4 @@ def n_silicon_dioxide(wavelength): gamma = (3 * chi_3 * omega_p0) / (4 * epsilon_0 * c**2 * n_p0**2 * A_eff_m2) -print("gamma:",gamma) +print("gamma:",gamma) \ No newline at end of file From 2cc8836427d2c22ad1f517f52401e2f7670775cb Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Sun, 3 Mar 2024 18:44:10 -0800 Subject: [PATCH 4/8] fix overlap calculations --- femwell/examples/Aeff.py | 85 +++++++++++++++++++++------------------- 1 file changed, 44 insertions(+), 41 deletions(-) diff --git a/femwell/examples/Aeff.py b/femwell/examples/Aeff.py index 9d343aa1..d9c8096c 100644 --- a/femwell/examples/Aeff.py +++ b/femwell/examples/Aeff.py @@ -1,23 +1,19 @@ """Aeff analysis based on https://optics.ansys.com/hc/en-us/articles/15100783091731-Spontaneous-Four-wave-Mixing-SWFM-Microring-Resonator-Photon-Source.""" from collections import OrderedDict + import numpy as np +from scipy.constants import c, epsilon_0 from shapely.geometry import box from shapely.ops import clip_by_rect -from skfem import Basis, ElementTriP0 -from skfem import Functional +from skfem import Basis, ElementTriP0, Functional from skfem.helpers import dot from skfem.io.meshio import from_meshio + from femwell.maxwell.waveguide import compute_modes from femwell.mesh import mesh_from_OrderedDict -from scipy.constants import c, epsilon_0 -def calculate_sfwm_Aeff( - basis: Basis, - mode_p, - mode_s, - mode_i -) -> np.complex64: - + +def calculate_sfwm_Aeff(basis: Basis, mode_p, mode_s, mode_i) -> np.complex64: """ Calculates the overlap integral for SFWM process by considering the interactions between pump, signal, and idler modes in the xy plane. @@ -30,46 +26,43 @@ def calculate_sfwm_Aeff( np.complex64: The Aeff result for the SFWM process(1/overlap integral). """ - def normalize_mode(mode): + def normalization_factor_mode(mode): @Functional def E2(w): return dot(w["E"][0], np.conj(w["E"][0])) - E_xy, _ = mode.basis.interpolate(mode.E) # [0]=xy [1]=z - E_squared_integral = E2.assemble(mode.basis, E=E_xy) + E = mode.basis.interpolate(mode.E) # [0]=xy [1]=z + E_squared_integral = E2.assemble(mode.basis, E=E) normalization_factor = 1 / np.sqrt(E_squared_integral) return normalization_factor # Return the normalization factor instead of modifying the mode - # Calculate normalization factors for each mode - norm_factor_p = normalize_mode(mode_p) - norm_factor_s = normalize_mode(mode_s) - norm_factor_i = normalize_mode(mode_i) - # Apply normalization factors to the electric fields for the overlap calculation - E_p, _ = mode_p.basis.interpolate(mode_p.E * norm_factor_p) - E_s, _ = mode_s.basis.interpolate(mode_s.E * norm_factor_s) - E_i, _ = mode_i.basis.interpolate(mode_i.E * norm_factor_i) - @Functional(dtype=np.complex64) def sfwm_overlap(w): - overlap_SFWM = w["E_p"][0] * w["E_p"][0] * np.conj(w["E_s"][0]) * np.conj(w["E_i"][0]) - return overlap_SFWM + return dot(w["E_p"][0], w["E_p"][0]) * dot(np.conj(w["E_s"][0]), np.conj(w["E_i"][0])) - overlap_result = sfwm_overlap.assemble(basis, E_p=E_p, E_s=E_s, E_i=E_i) + overlap_result = sfwm_overlap.assemble( + basis, + E_p=mode_p.basis.interpolate(mode_p.E * normalization_factor_mode(mode_p)), + E_s=mode_s.basis.interpolate(mode_s.E * normalization_factor_mode(mode_s)), + E_i=mode_i.basis.interpolate(mode_i.E * normalization_factor_mode(mode_i)), + ) return 1 / overlap_result -#Dispersion relations of materials -#Core +# Dispersion relations of materials +# Core def n_X(wavelength): x = wavelength - return (1 - + 2.19244563 / (1 - (0.20746607 / x) ** 2) - + 0.13435116 / (1 - (0.3985835 / x) ** 2) - + 2.20997784 / (1 - (0.20747044 / x) ** 2) + return ( + 1 + + 2.19244563 / (1 - (0.20746607 / x) ** 2) + + 0.13435116 / (1 - (0.3985835 / x) ** 2) + + 2.20997784 / (1 - (0.20747044 / x) ** 2) ) ** 0.5 + # Box def n_silicon_dioxide(wavelength): x = wavelength @@ -80,27 +73,33 @@ def n_silicon_dioxide(wavelength): + 0.8974794 / (1 - (9.896161 / x) ** 2) ) ** 0.5 -Clad=1 + +Clad = 1 -core = box(0, 0, 0.5, 0.39) #500x390nm +core = box(0, 0, 1.05, 0.5) # 1050x500nm +# core = box(0, 0, .5, 0.39) # 500x390nm polygons = OrderedDict( core=core, box=clip_by_rect(core.buffer(1.5, resolution=4), -np.inf, -np.inf, np.inf, 0), clad=clip_by_rect(core.buffer(1.5, resolution=4), -np.inf, 0, np.inf, np.inf), ) -resolutions = {"core": {"resolution": 0.025, "distance": 2.}} +resolutions = {"core": {"resolution": 0.025, "distance": 2.0}} -mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=0.6)) +mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=0.2)) num_modes = 2 -#For SFWM we have energy conservation and momemtum(k) conservation for 2pumps and signal+idler +# For SFWM we have energy conservation and momemtum(k) conservation for 2pumps and signal+idler lambda_p0 = 1.4 lambda_s0 = 1.097 lambda_i0 = 1.686 +# lambda_p0 = 1.55 +# lambda_s0 = 1.55 +# lambda_i0 = 1.55 + basis0 = Basis(mesh, ElementTriP0()) epsilon_p = basis0.zeros() @@ -108,7 +107,9 @@ def n_silicon_dioxide(wavelength): epsilon_i = basis0.zeros() -for wavelength, epsilon in zip([lambda_p0, lambda_s0, lambda_i0], [epsilon_p, epsilon_s, epsilon_i]): +for wavelength, epsilon in zip( + [lambda_p0, lambda_s0, lambda_i0], [epsilon_p, epsilon_s, epsilon_i] +): for subdomain, n_func in { "core": n_X, "box": n_silicon_dioxide, @@ -122,6 +123,7 @@ def n_silicon_dioxide(wavelength): modes_s = compute_modes(basis0, epsilon_s, wavelength=lambda_s0, num_modes=num_modes, order=1) modes_i = compute_modes(basis0, epsilon_i, wavelength=lambda_i0, num_modes=num_modes, order=1) +# modes_p[0].show(modes_p[0].E.real) mode_p = max(modes_p, key=lambda mode: mode.te_fraction) mode_s = max(modes_s, key=lambda mode: mode.te_fraction) @@ -130,14 +132,15 @@ def n_silicon_dioxide(wavelength): A_eff = calculate_sfwm_Aeff(basis0, mode_p, mode_s, mode_i) print("Aeff in um2:", A_eff) -#Calculation for non-linear coef +# Calculation for non-linear coef chi_3 = 5e-21 # m^2/V^2 #7e-20? -lambda_p0_m = lambda_p0 * 1e-6 #to m +chi_3 = 7e-20 # ? +lambda_p0_m = lambda_p0 * 1e-6 # to m n_p0 = np.real(mode_p.n_eff) -A_eff_m2 = A_eff * 1e-12 #to m^2 +A_eff_m2 = A_eff * 1e-12 # to m^2 omega_p0 = 2 * np.pi * c / lambda_p0_m gamma = (3 * chi_3 * omega_p0) / (4 * epsilon_0 * c**2 * n_p0**2 * A_eff_m2) -print("gamma:",gamma) \ No newline at end of file +print("gamma:", gamma) From 145fbfe05752f44c89c99e268c1a0bd32b0a14a4 Mon Sep 17 00:00:00 2001 From: Kunyuan-LI Date: Tue, 5 Mar 2024 05:54:01 +0100 Subject: [PATCH 5/8] move to docs/photonics/examples --- .../examples/Aeff.py => docs/photonics/examples/SFWM_Aeff.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) rename femwell/examples/Aeff.py => docs/photonics/examples/SFWM_Aeff.py (98%) diff --git a/femwell/examples/Aeff.py b/docs/photonics/examples/SFWM_Aeff.py similarity index 98% rename from femwell/examples/Aeff.py rename to docs/photonics/examples/SFWM_Aeff.py index d9c8096c..aab86eca 100644 --- a/femwell/examples/Aeff.py +++ b/docs/photonics/examples/SFWM_Aeff.py @@ -129,12 +129,12 @@ def n_silicon_dioxide(wavelength): mode_s = max(modes_s, key=lambda mode: mode.te_fraction) mode_i = max(modes_i, key=lambda mode: mode.te_fraction) -A_eff = calculate_sfwm_Aeff(basis0, mode_p, mode_s, mode_i) +A_eff = np.real(calculate_sfwm_Aeff(basis0, mode_p, mode_s, mode_i)) print("Aeff in um2:", A_eff) # Calculation for non-linear coef chi_3 = 5e-21 # m^2/V^2 #7e-20? -chi_3 = 7e-20 # ? +#chi_3 = 7e-20 # ? lambda_p0_m = lambda_p0 * 1e-6 # to m n_p0 = np.real(mode_p.n_eff) A_eff_m2 = A_eff * 1e-12 # to m^2 From 2c9a90be87f97c9a2988bde0667126756f91665e Mon Sep 17 00:00:00 2001 From: Kunyuan-LI Date: Wed, 6 Mar 2024 02:40:20 +0100 Subject: [PATCH 6/8] reevaluation of overlap intengral --- docs/photonics/examples/SFWM_Aeff.py | 10 +-- femwell/tests/test_Aeff_mode.py | 99 ++++++++++++++++++++++++++++ 2 files changed, 104 insertions(+), 5 deletions(-) create mode 100644 femwell/tests/test_Aeff_mode.py diff --git a/docs/photonics/examples/SFWM_Aeff.py b/docs/photonics/examples/SFWM_Aeff.py index aab86eca..f254b375 100644 --- a/docs/photonics/examples/SFWM_Aeff.py +++ b/docs/photonics/examples/SFWM_Aeff.py @@ -39,8 +39,9 @@ def E2(w): # Apply normalization factors to the electric fields for the overlap calculation @Functional(dtype=np.complex64) def sfwm_overlap(w): - return dot(w["E_p"][0], w["E_p"][0]) * dot(np.conj(w["E_s"][0]), np.conj(w["E_i"][0])) - + #return dot(w["E_p"][0], w["E_p"][0]) * dot(np.conj(w["E_s"][0]), np.conj(w["E_i"][0]))#? + return dot(w["E_p"][0], np.conj(w["E_s"][0])) * dot(w["E_p"][0], np.conj(w["E_i"][0]))#?? + #return np.dot(w["E_p"][0][1], np.conj(w["E_s"][0][1]).T) * np.dot(w["E_p"][0][1], np.conj(w["E_i"][0][1]).T)#??? X overlap_result = sfwm_overlap.assemble( basis, E_p=mode_p.basis.interpolate(mode_p.E * normalization_factor_mode(mode_p)), @@ -77,8 +78,8 @@ def n_silicon_dioxide(wavelength): Clad = 1 -core = box(0, 0, 1.05, 0.5) # 1050x500nm -# core = box(0, 0, .5, 0.39) # 500x390nm +#core = box(0, 0, 1.05, 0.5) # 1050x500nm +core = box(0, 0, .5, 0.39) # 500x390nm polygons = OrderedDict( core=core, box=clip_by_rect(core.buffer(1.5, resolution=4), -np.inf, -np.inf, np.inf, 0), @@ -134,7 +135,6 @@ def n_silicon_dioxide(wavelength): # Calculation for non-linear coef chi_3 = 5e-21 # m^2/V^2 #7e-20? -#chi_3 = 7e-20 # ? lambda_p0_m = lambda_p0 * 1e-6 # to m n_p0 = np.real(mode_p.n_eff) A_eff_m2 = A_eff * 1e-12 # to m^2 diff --git a/femwell/tests/test_Aeff_mode.py b/femwell/tests/test_Aeff_mode.py new file mode 100644 index 00000000..477a1343 --- /dev/null +++ b/femwell/tests/test_Aeff_mode.py @@ -0,0 +1,99 @@ +from collections import OrderedDict +import matplotlib.pyplot as plt +import numpy as np +from shapely.ops import clip_by_rect +from skfem import Basis, ElementTriP0 +from skfem.io.meshio import from_meshio +from shapely.geometry import box +from femwell.maxwell.waveguide import compute_modes +from femwell.mesh import mesh_from_OrderedDict + + +def n_X(wavelength): + x = wavelength + return (1 + + 2.19244563 / (1 - (0.20746607 / x) ** 2) + + 0.13435116 / (1 - (0.3985835 / x) ** 2) + + 2.20997784 / (1 - (0.20747044 / x) ** 2) + ) ** 0.5 + +# Box +def n_silicon_dioxide(wavelength): + x = wavelength + return ( + 1 + + 0.6961663 / (1 - (0.0684043 / x) ** 2) + + 0.4079426 / (1 - (0.1162414 / x) ** 2) + + 0.8974794 / (1 - (9.896161 / x) ** 2) + ) ** 0.5 + +Clad=1 + +core = box(0, 0, 0.5, 0.39) +polygons = OrderedDict( + core=core, + box=clip_by_rect(core.buffer(1.5, resolution=4), -np.inf, -np.inf, np.inf, 0), + clad=clip_by_rect(core.buffer(1.5, resolution=4), -np.inf, 0, np.inf, np.inf), +) + +resolutions = {"core": {"resolution": 0.025, "distance": 2.}} + +mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=0.6)) + +num_modes = 2 + +lambda_p0 = 1.4 +lambda_s0 = 1.097 +lambda_i0 = 1.686 + +basis0 = Basis(mesh, ElementTriP0()) + +epsilon_p = basis0.zeros(dtype=complex) +epsilon_s = basis0.zeros(dtype=complex) +epsilon_i = basis0.zeros(dtype=complex) + +for wavelength, epsilon in zip([lambda_p0, lambda_s0, lambda_i0], [epsilon_p, epsilon_s, epsilon_i]): + for subdomain, n_func in { + "core": n_X, + "box": n_silicon_dioxide, + "clad": lambda _: Clad, + }.items(): + n = n_func(wavelength) + epsilon[basis0.get_dofs(elements=subdomain)] = n**2 +modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=2, order=2) + + +modes_p = compute_modes(basis0, epsilon_p, wavelength=lambda_p0, num_modes=num_modes, order=1) +modes_s = compute_modes(basis0, epsilon_s, wavelength=lambda_s0, num_modes=num_modes, order=1) +modes_i = compute_modes(basis0, epsilon_i, wavelength=lambda_i0, num_modes=num_modes, order=1) + +mode_p = max(modes_p, key=lambda mode: mode.te_fraction) +mode_s = max(modes_s, key=lambda mode: mode.te_fraction) +mode_i = max(modes_i, key=lambda mode: mode.te_fraction) + +print(f"Effective refractive index for p mode: {np.real(mode_p.n_eff):.4f}") +mode_p.show(mode_p.E.real, colorbar=True, direction="x") +mode_p.show(mode_p.E.imag, colorbar=True, direction="x") +fig, ax = plt.subplots() +mode_p.plot_intensity(ax=ax) +plt.title("Normalized Intensity for p mode") +plt.tight_layout() +plt.show() + +print(f"Effective refractive index for s mode: {np.real(mode_s.n_eff):.4f}") +mode_s.show(mode_s.E.real, colorbar=True, direction="x") +mode_s.show(mode_s.E.imag, colorbar=True, direction="x") +fig, ax = plt.subplots() +mode_s.plot_intensity(ax=ax) +plt.title("Normalized Intensity for s mode") +plt.tight_layout() +plt.show() + +print(f"Effective refractive index for i mode: {np.real(mode_i.n_eff):.4f}") +mode_i.show(mode_i.E.real, colorbar=True, direction="x") +mode_i.show(mode_i.E.imag, colorbar=True, direction="x") +fig, ax = plt.subplots() +mode_i.plot_intensity(ax=ax) +plt.title("Normalized Intensity for i mode") +plt.tight_layout() +plt.show() From c6aa0871aa749944f3fe492dbcfa8d321d5c5a40 Mon Sep 17 00:00:00 2001 From: Kunyuan-LI Date: Tue, 12 Mar 2024 02:59:59 +0100 Subject: [PATCH 7/8] overlap intergal function modified --- docs/photonics/examples/SFWM_Aeff.py | 34 ++++++++++++++++++++-------- femwell/tests/test_Aeff_mode.py | 6 ++--- 2 files changed, 27 insertions(+), 13 deletions(-) diff --git a/docs/photonics/examples/SFWM_Aeff.py b/docs/photonics/examples/SFWM_Aeff.py index f254b375..26e5e755 100644 --- a/docs/photonics/examples/SFWM_Aeff.py +++ b/docs/photonics/examples/SFWM_Aeff.py @@ -39,9 +39,23 @@ def E2(w): # Apply normalization factors to the electric fields for the overlap calculation @Functional(dtype=np.complex64) def sfwm_overlap(w): + + E_p_xy = w["E_p"][0] # shape: (2, x, 3) + E_s_xy = w["E_s"][0] # shape: (2, x, 3) + E_i_xy = w["E_i"][0] # shape: (2, x, 3) + + overlap_Ex = E_p_xy[0,:,0]*E_p_xy[0,:,0]*np.conj(E_s_xy[0,:,0])*np.conj(E_i_xy[0,:,0]) + \ + E_p_xy[1,:,0]*E_p_xy[1,:,0]*np.conj(E_s_xy[1,:,0])*np.conj(E_i_xy[1,:,0]) + overlap_Ey = E_p_xy[0,:,1]*E_p_xy[0,:,1]*np.conj(E_s_xy[0,:,1])*np.conj(E_i_xy[0,:,1]) + \ + E_p_xy[1,:,1]*E_p_xy[1,:,1]*np.conj(E_s_xy[1,:,1])*np.conj(E_i_xy[1,:,1]) + overlap_Ez = E_p_xy[0,:,2]*E_p_xy[0,:,2]*np.conj(E_s_xy[0,:,2])*np.conj(E_i_xy[0,:,2]) + \ + E_p_xy[1,:,2]*E_p_xy[1,:,2]*np.conj(E_s_xy[1,:,2])*np.conj(E_i_xy[1,:,2]) + + return np.array([overlap_Ex, overlap_Ey, overlap_Ez]).T + #return dot(w["E_p"][0], w["E_p"][0]) * dot(np.conj(w["E_s"][0]), np.conj(w["E_i"][0]))#? - return dot(w["E_p"][0], np.conj(w["E_s"][0])) * dot(w["E_p"][0], np.conj(w["E_i"][0]))#?? - #return np.dot(w["E_p"][0][1], np.conj(w["E_s"][0][1]).T) * np.dot(w["E_p"][0][1], np.conj(w["E_i"][0][1]).T)#??? X + #return dot(w["E_p"][0], np.conj(w["E_s"][0])) * dot(w["E_p"][0], np.conj(w["E_i"][0]))#?? + overlap_result = sfwm_overlap.assemble( basis, E_p=mode_p.basis.interpolate(mode_p.E * normalization_factor_mode(mode_p)), @@ -78,8 +92,8 @@ def n_silicon_dioxide(wavelength): Clad = 1 -#core = box(0, 0, 1.05, 0.5) # 1050x500nm -core = box(0, 0, .5, 0.39) # 500x390nm +core = box(0, 0, 1.05, 0.5) # 1050x500nm +#core = box(0, 0, .5, 0.39) # 500x390nm polygons = OrderedDict( core=core, box=clip_by_rect(core.buffer(1.5, resolution=4), -np.inf, -np.inf, np.inf, 0), @@ -93,13 +107,13 @@ def n_silicon_dioxide(wavelength): num_modes = 2 # For SFWM we have energy conservation and momemtum(k) conservation for 2pumps and signal+idler -lambda_p0 = 1.4 -lambda_s0 = 1.097 -lambda_i0 = 1.686 +# lambda_p0 = 1.4 +# lambda_s0 = 1.097 +# lambda_i0 = 1.686 -# lambda_p0 = 1.55 -# lambda_s0 = 1.55 -# lambda_i0 = 1.55 +lambda_p0 = 1.55 +lambda_s0 = 1.55 +lambda_i0 = 1.55 basis0 = Basis(mesh, ElementTriP0()) diff --git a/femwell/tests/test_Aeff_mode.py b/femwell/tests/test_Aeff_mode.py index 477a1343..3a5e012d 100644 --- a/femwell/tests/test_Aeff_mode.py +++ b/femwell/tests/test_Aeff_mode.py @@ -48,9 +48,9 @@ def n_silicon_dioxide(wavelength): basis0 = Basis(mesh, ElementTriP0()) -epsilon_p = basis0.zeros(dtype=complex) -epsilon_s = basis0.zeros(dtype=complex) -epsilon_i = basis0.zeros(dtype=complex) +epsilon_p = basis0.zeros() +epsilon_s = basis0.zeros() +epsilon_i = basis0.zeros() for wavelength, epsilon in zip([lambda_p0, lambda_s0, lambda_i0], [epsilon_p, epsilon_s, epsilon_i]): for subdomain, n_func in { From 99523ef9557a788b4612df738148f257974ac61c Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Mon, 18 Mar 2024 23:22:09 -0700 Subject: [PATCH 8/8] add four wave mixing example to toc and add markdown --- docs/_toc.yml | 1 + docs/photonics/examples/SFWM_Aeff.py | 109 ++++++++++++++++++++++----- docs/references.bib | 42 +++++++++++ 3 files changed, 134 insertions(+), 18 deletions(-) diff --git a/docs/_toc.yml b/docs/_toc.yml index 3c27db04..eade99a6 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -39,6 +39,7 @@ parts: - file: photonics/examples/coupled_mode_theory.py - file: photonics/examples/depletion_waveguide.py - file: photonics/examples/effective_area.py + - file: photonics/examples/SFWM_Aeff.py - file: photonics/examples/refinement.py - file: photonics/examples/propagation_loss.py - file: electronics/examples/capacitor.py diff --git a/docs/photonics/examples/SFWM_Aeff.py b/docs/photonics/examples/SFWM_Aeff.py index 26e5e755..cd671008 100644 --- a/docs/photonics/examples/SFWM_Aeff.py +++ b/docs/photonics/examples/SFWM_Aeff.py @@ -1,4 +1,23 @@ -"""Aeff analysis based on https://optics.ansys.com/hc/en-us/articles/15100783091731-Spontaneous-Four-wave-Mixing-SWFM-Microring-Resonator-Photon-Source.""" +# --- +# jupyter: +# jupytext: +# formats: py:percent,md:myst +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.15.0 +# kernelspec: +# display_name: Python 3 +# name: python3 +# --- + +# %% [markdown] +# # Optimisation of spontaneous four-wave mixing in a ring microcavity + +# Here we reproduce {cite}`Chuprina2017` + +# %% tags=["remove-stderr", "hide-input", "thebe-init"] from collections import OrderedDict import numpy as np @@ -12,7 +31,28 @@ from femwell.maxwell.waveguide import compute_modes from femwell.mesh import mesh_from_OrderedDict - +# %% [markdown] +# +# $$ +# \int\int\left|u(x,y)\right|^{2}\mathrm{d}x\mathrm{d}y=1 +# $$ +# +# $$ +# A_{\mathrm{eff}} +# = +# \frac{1}{ +# \displaystyle +# \int\int +# \mathrm{d}x \mathrm{d}y +# u_{\mathrm{p}}(x,y) +# u_{\mathrm{p}}(x,y) +# u_{\mathrm{s}}^{*}(x,y) +# u_{\mathrm{i}}^{*}(x,y) +# } +# $$ + + +# %% def calculate_sfwm_Aeff(basis: Basis, mode_p, mode_s, mode_i) -> np.complex64: """ Calculates the overlap integral for SFWM process by considering the interactions @@ -43,19 +83,22 @@ def sfwm_overlap(w): E_p_xy = w["E_p"][0] # shape: (2, x, 3) E_s_xy = w["E_s"][0] # shape: (2, x, 3) E_i_xy = w["E_i"][0] # shape: (2, x, 3) - - overlap_Ex = E_p_xy[0,:,0]*E_p_xy[0,:,0]*np.conj(E_s_xy[0,:,0])*np.conj(E_i_xy[0,:,0]) + \ - E_p_xy[1,:,0]*E_p_xy[1,:,0]*np.conj(E_s_xy[1,:,0])*np.conj(E_i_xy[1,:,0]) - overlap_Ey = E_p_xy[0,:,1]*E_p_xy[0,:,1]*np.conj(E_s_xy[0,:,1])*np.conj(E_i_xy[0,:,1]) + \ - E_p_xy[1,:,1]*E_p_xy[1,:,1]*np.conj(E_s_xy[1,:,1])*np.conj(E_i_xy[1,:,1]) - overlap_Ez = E_p_xy[0,:,2]*E_p_xy[0,:,2]*np.conj(E_s_xy[0,:,2])*np.conj(E_i_xy[0,:,2]) + \ - E_p_xy[1,:,2]*E_p_xy[1,:,2]*np.conj(E_s_xy[1,:,2])*np.conj(E_i_xy[1,:,2]) - + + overlap_Ex = E_p_xy[0, :, 0] * E_p_xy[0, :, 0] * np.conj(E_s_xy[0, :, 0]) * np.conj( + E_i_xy[0, :, 0] + ) + E_p_xy[1, :, 0] * E_p_xy[1, :, 0] * np.conj(E_s_xy[1, :, 0]) * np.conj(E_i_xy[1, :, 0]) + overlap_Ey = E_p_xy[0, :, 1] * E_p_xy[0, :, 1] * np.conj(E_s_xy[0, :, 1]) * np.conj( + E_i_xy[0, :, 1] + ) + E_p_xy[1, :, 1] * E_p_xy[1, :, 1] * np.conj(E_s_xy[1, :, 1]) * np.conj(E_i_xy[1, :, 1]) + overlap_Ez = E_p_xy[0, :, 2] * E_p_xy[0, :, 2] * np.conj(E_s_xy[0, :, 2]) * np.conj( + E_i_xy[0, :, 2] + ) + E_p_xy[1, :, 2] * E_p_xy[1, :, 2] * np.conj(E_s_xy[1, :, 2]) * np.conj(E_i_xy[1, :, 2]) + return np.array([overlap_Ex, overlap_Ey, overlap_Ez]).T - - #return dot(w["E_p"][0], w["E_p"][0]) * dot(np.conj(w["E_s"][0]), np.conj(w["E_i"][0]))#? - #return dot(w["E_p"][0], np.conj(w["E_s"][0])) * dot(w["E_p"][0], np.conj(w["E_i"][0]))#?? - + + # return dot(w["E_p"][0], w["E_p"][0]) * dot(np.conj(w["E_s"][0]), np.conj(w["E_i"][0]))#? + # return dot(w["E_p"][0], np.conj(w["E_s"][0])) * dot(w["E_p"][0], np.conj(w["E_i"][0]))#?? + overlap_result = sfwm_overlap.assemble( basis, E_p=mode_p.basis.interpolate(mode_p.E * normalization_factor_mode(mode_p)), @@ -66,8 +109,13 @@ def sfwm_overlap(w): return 1 / overlap_result +# %% [markdown] +# # Dispersion relations of materials -# Core +# Silicon nitride {cite}`Luke2015` + + +# %% def n_X(wavelength): x = wavelength return ( @@ -78,7 +126,12 @@ def n_X(wavelength): ) ** 0.5 -# Box +# %% [markdown] +# +# Box {cite}`Malitson1965` + + +# %% def n_silicon_dioxide(wavelength): x = wavelength return ( @@ -91,9 +144,14 @@ def n_silicon_dioxide(wavelength): Clad = 1 +# %% [markdown] +# +# Waveguide dimensions 1050x500nm + +# %% -core = box(0, 0, 1.05, 0.5) # 1050x500nm -#core = box(0, 0, .5, 0.39) # 500x390nm +core = box(0, 0, 1.05, 0.5) +# core = box(0, 0, .5, 0.39) # 500x390nm polygons = OrderedDict( core=core, box=clip_by_rect(core.buffer(1.5, resolution=4), -np.inf, -np.inf, np.inf, 0), @@ -155,6 +213,21 @@ def n_silicon_dioxide(wavelength): omega_p0 = 2 * np.pi * c / lambda_p0_m +# %% [markdown] +# $$ +# \gamma=\frac{3\chi^{(3)}\omega_{\mathrm{p}}}{4\varepsilon_{0}c^{2}n_{\mathrm{p0}}^{2}A_{\mathrm{eff}}} +# $$ + +# %% + gamma = (3 * chi_3 * omega_p0) / (4 * epsilon_0 * c**2 * n_p0**2 * A_eff_m2) print("gamma:", gamma) + +# %% [markdown] +# ## Bibliography +# +# ```{bibliography} +# :style: unsrt +# :filter: docname in docnames +# ``` diff --git a/docs/references.bib b/docs/references.bib index db030778..d086b7d0 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -328,3 +328,45 @@ @article{Klenner2016 month = may, pages = {11043} } +@article{Chuprina2017, + title = {Optimisation of spontaneous four-wave mixing in a ring microcavity}, + volume = {47}, + issn = {1468-4799}, + url = {http://dx.doi.org/10.1070/QEL16511}, + doi = {10.1070/qel16511}, + number = {10}, + journal = {Quantum Electronics}, + publisher = {IOP Publishing}, + author = {Chuprina, I N and An, P P and Zubkova, E G and Kovalyuk, V V and Kalachev, A A and Gol'tsman, G N}, + year = {2017}, + month = oct, + pages = {887–891} +} +@article{Luke2015, + title = {Broadband mid-infrared frequency comb generation in a Si\_3N\_4 microresonator}, + volume = {40}, + issn = {1539-4794}, + url = {http://dx.doi.org/10.1364/OL.40.004823}, + doi = {10.1364/ol.40.004823}, + number = {21}, + journal = {Optics Letters}, + publisher = {The Optical Society}, + author = {Luke, Kevin and Okawachi, Yoshitomo and Lamont, Michael R. E. and Gaeta, Alexander L. and Lipson, Michal}, + year = {2015}, + month = oct, + pages = {4823} +} +@article{Malitson1965, + title = {Interspecimen Comparison of the Refractive Index of Fused Silica\textasteriskcentered, \textdagger{}}, + volume = {55}, + issn = {0030-3941}, + url = {http://dx.doi.org/10.1364/JOSA.55.001205}, + doi = {10.1364/josa.55.001205}, + number = {10}, + journal = {Journal of the Optical Society of America}, + publisher = {Optica Publishing Group}, + author = {Malitson, I. H.}, + year = {1965}, + month = oct, + pages = {1205} +}