From aff72ddd65ae2a65f949d26ff147753c5863049a Mon Sep 17 00:00:00 2001 From: simbilod Date: Wed, 9 Nov 2022 22:19:32 -0800 Subject: [PATCH 1/3] add simple empirical waveguide loss model to tidy3d --- gdsfactory/simulation/gtidy3d/modes.py | 100 ++++++++++++++++++++----- 1 file changed, 81 insertions(+), 19 deletions(-) diff --git a/gdsfactory/simulation/gtidy3d/modes.py b/gdsfactory/simulation/gtidy3d/modes.py index eae1ed6dc3..a4bc8155bf 100644 --- a/gdsfactory/simulation/gtidy3d/modes.py +++ b/gdsfactory/simulation/gtidy3d/modes.py @@ -157,7 +157,11 @@ class Waveguide(BaseModel): cache: filepath for caching modes. If None does not use file cache. precision: single or double. filter_pol: te, tm or None. - + loss_model: whether to include a scattering loss region at the interfaces. Default is False + sidewall_sigma: size of the region to append to the sidewalls in the loss model. Default is 10 nm. + sidewall_k: imaginary index addition for the sidewall loss region. Default is 0 (no extra loss). + top_sigma: size of the loss region to append to the top surfaces in the loss model. Default is 10 nm. + top_k: imaginary index addition for the top surfaces loss region. Default is 0 (no extra loss). :: __________________________ @@ -196,6 +200,12 @@ class Waveguide(BaseModel): precision: Precision = "single" filter_pol: Optional[FilterPol] = None + loss_model: Optional[bool] = False + sidewall_sigma: Optional[float] = 10 * nm + sidewall_k: Optional[float] = 0.1 + top_sigma: Optional[float] = 10 * nm + top_k: Optional[float] = 0.1 + class Config: """Config for Waveguide.""" @@ -257,6 +267,8 @@ def get_n(self, Y, Z): complex_solver = True elif self.dn_dict is not None: complex_solver = True + elif self.loss_model: + complex_solver = True if complex_solver: mat_dtype = np.complex128 if self.precision == "double" else np.complex64 elif self.precision == "double": @@ -273,6 +285,28 @@ def get_n(self, Y, Z): n[inds_core] = ncore n[inds_slab] = ncore if slab_thickness else nclad + if self.loss_model: + inds_top = ( + (Z >= t_box + wg_thickness - self.top_sigma/2) & (Z <= t_box + wg_thickness + self.top_sigma/2) & (-w / 2 <= Y) & (Y <= w / 2) + ) + inds_top_slab_left = ( + (Z >= t_box + slab_thickness - self.top_sigma/2) & (Z <= t_box + slab_thickness + self.top_sigma/2) & (-w / 2 >= Y) + ) + inds_top_slab_right = ( + (Z >= t_box + slab_thickness - self.top_sigma/2) & (Z <= t_box + slab_thickness + self.top_sigma/2) & (Y >= w / 2) + ) + inds_sidewall_left = ( + (Z >= t_box + slab_thickness) & (Z <= t_box + wg_thickness) & (Y >= -w / 2 - self.sidewall_sigma/2) & (Y <= -w / 2 + self.sidewall_sigma/2) + ) + inds_sidewall_right = ( + (Z >= t_box + slab_thickness) & (Z <= t_box + wg_thickness) & (Y >= w / 2 - self.sidewall_sigma/2) & (Y <= w / 2 + self.sidewall_sigma/2) + ) + n[inds_top] += 1j*self.top_k + n[inds_top_slab_left] += 1j*self.top_k + n[inds_top_slab_right] += 1j*self.top_k + n[inds_sidewall_left] += 1j*self.sidewall_k + n[inds_sidewall_right] += 1j*self.sidewall_k + if self.dn_dict is not None: dn = griddata( (self.dn_dict["x"], self.dn_dict["y"]), @@ -286,7 +320,8 @@ def get_n(self, Y, Z): return n - def plot_index(self) -> None: + + def plot_index(self, func=None) -> None: x, y, Xx, Yx, Xy, Yy, Xz, Yz = create_mesh( -self.w_sim / 2, 0.0, @@ -304,7 +339,10 @@ def plot_index(self) -> None: Xx, Yx, ) - plot(Xx, Yx, nx) + if func == None: + plot(Xx, Yx, nx) + else: + plot(Xx, Yx, func(nx)) plt.show() def compute_modes( @@ -620,10 +658,16 @@ def get_overlap( - wg2.Ey[..., mode_index2] * np.conj(wg1.Hx[..., mode_index1]) ) + def get_loss(self): + """Returns loss for computed modes in dB/cm""" + if not hasattr(self, "neffs"): + self.compute_modes() + wavelength = self.wavelength * 1e-6 # convert to m + alphas = 4 * np.pi * np.imag(self.neffs) / wavelength # lin/m loss + return 10*np.log10(np.exp(1)) * alphas * 1e-2 # dB/cm loss class WaveguideCoupler(Waveguide): """Waveguide coupler Model. - Parameters: wavelength: (um). wg_width1: left waveguide width in um. @@ -1039,19 +1083,37 @@ def plot_sweep_width( "group_index", ) + if __name__ == "__main__": - widths = np.arange(400, 601, 50) * 1e-3 - widths = np.array([500]) * nm - thicknesses = np.array([210, 220, 230]) * nm - widths = np.array([490, 500, 510]) * nm - widths = np.array([495, 500, 505]) * nm - thicknesses = np.array([220]) * nm - - df = plot_sweep_width( - steps=3, - wavelength=1.55, - wg_thickness=220 * nm, - slab_thickness=0 * nm, - ncore=si, - nclad=sio2, - ) + + wg = Waveguide(nmodes=2, + wg_width=500*nm, + wavelength=1.55, + wg_thickness=220 * nm, + slab_thickness=90 * nm, + ncore=si, + nclad=sio2, + loss_model=True, + sidewall_k = 1E-4, + top_k = 0, + sidewall_sigma = 20*nm, + top_sigma = 0, + resolution = 400, + cache=None, + precision='double' + ) + wg.plot_index() + wg.plot_index(func=np.imag) + wg.compute_mode_properties() + print(wg.neffs) + print(wg.get_loss()) + # wg.plot_Ex() + # wg.compute_mode_properties() + # for mode_number in range(nmodes): + # n = np.real(wg.neffs[mode_number]) + # fraction_te = wg.fraction_te[mode_number] + # plt.scatter(wg_width, n, c=fraction_te, vmin=0, vmax=1, cmap=cmap) + # n[inds_top_slab_left] += self.top_k + # n[inds_top_slab_right] += self.top_k + # n[inds_sidewall_left] += self.sidewall_k + # n[inds_sidewall_right] += self.sidewall_k \ No newline at end of file From e97ee3bedc6299938e9e12728a0d6ae169586166 Mon Sep 17 00:00:00 2001 From: simbilod Date: Wed, 9 Nov 2022 22:45:02 -0800 Subject: [PATCH 2/3] pre-commit --- gdsfactory/simulation/gtidy3d/modes.py | 34 ++------------------------ 1 file changed, 2 insertions(+), 32 deletions(-) diff --git a/gdsfactory/simulation/gtidy3d/modes.py b/gdsfactory/simulation/gtidy3d/modes.py index a4bc8155bf..485c051914 100644 --- a/gdsfactory/simulation/gtidy3d/modes.py +++ b/gdsfactory/simulation/gtidy3d/modes.py @@ -640,10 +640,8 @@ def get_overlap( self, wg: "Waveguide", mode_index1: int = 0, mode_index2: int = 0 ) -> float: """Returns mode overlap integral. - Args: wg: other waveguide. - """ wg1 = self wg2 = wg @@ -683,11 +681,8 @@ class WaveguideCoupler(Waveguide): nmodes: number of modes to compute. bend_radius: optional bend radius (um). cache: filepath for caching modes. If None does not use file cache. - :: - -w1-gap/2 - wg_width1 wg_width2 <-------> <-------> _______ | _______ __ @@ -700,7 +695,6 @@ class WaveguideCoupler(Waveguide): gap <---------------------> w_sim - """ wg_width: Optional[float] = None @@ -714,11 +708,9 @@ def w_sim(self): def get_n(self, Y, Z): """Return index matrix for a waveguide coupler. - Args: Y: 2D array. Z: 2D array. - """ w1 = self.wg_width1 w2 = self.wg_width2 @@ -771,16 +763,13 @@ def sweep_bend_loss( **kwargs, ) -> Tuple[np.ndarray, np.ndarray]: """Returns overlap integral squared for the bend mode mismatch loss. - The loss is squared because you hit the bend loss twice (from bend to straight and from straight to bend). - Args: bend_radius_min: min bend radius (um). bend_radius_max: max bend radius (um). steps: number of steps. mode_index: where 0 is the fundamental mode. - Keyword Args: wavelength: (um). wg_width: waveguide width in um. @@ -793,7 +782,6 @@ def sweep_bend_loss( xmargin: margin from waveguide edge to each side (um). resolution: pixels/um. nmodes: number of modes to compute. - """ r = np.linspace(bend_radius_min, bend_radius_max, steps) integral = np.zeros_like(r) @@ -822,13 +810,11 @@ def sweep_neff( **kwargs, ) -> pd.DataFrame: """Sweep waveguide width and compute effective index. - Args: wavelength: (um). thicknesses: in um. widths: in um. mode_index: integer, where 0 is the fundamental mode. - Keyword Args: mode_index: integer. ncore: core refractive index. @@ -840,7 +826,6 @@ def sweep_neff( resolution: pixels/um. nmodes: number of modes to compute. bend_radius: optional bend radius (um). - """ widths_thicknesses = list(it.product(widths, thicknesses)) @@ -868,12 +853,10 @@ def group_index( wavelength: float, wavelength_step: float = 0.01, mode_index: int = 0, **kwargs ) -> float: """Returns group_index. - Args: wavelength: (um). wavelength_step: in um. mode_index: integer, where 0 is the fundamental mode. - Keyword Args: wg_width: waveguide width. wg_thickness: thickness waveguide (um). @@ -886,7 +869,6 @@ def group_index( resolution: pixels/um. nmodes: number of modes to compute. bend_radius: optional bend radius (um). - """ wc = Waveguide(wavelength=wavelength, **kwargs) wf = Waveguide( @@ -914,12 +896,10 @@ def sweep_group_index( **kwargs, ) -> pd.DataFrame: """Sweep waveguide width and compute group index. - Args: wavelength: (um). thicknesses: in um. widths: in um. - Keyword Args: mode_index: integer. ncore: core refractive index. @@ -931,7 +911,6 @@ def sweep_group_index( resolution: pixels/um. nmodes: number of modes to compute. bend_radius: optional bend radius (um). - """ widths_thicknesses = list(it.product(widths, thicknesses)) @@ -960,16 +939,12 @@ def sweep_width( **kwargs, ) -> pd.DataFrame: """Sweep waveguide width and compute effective index. - Returns pandas dataframe with effective index (neff) and fraction_te. - Args: width1: starting waveguide width in um. width2: end waveguide width in um. steps: number of points. nmodes: number of modes to compute. - - Keyword Args: wavelength: (um). wg_width: waveguide width in um. @@ -983,7 +958,6 @@ def sweep_width( resolution: pixels/um. nmodes: number of modes to compute. bend_radius: optional bend radius (um). - """ width = np.linspace(width1, width2, steps) neff = {} @@ -1016,16 +990,13 @@ def plot_sweep_width( **kwargs, ) -> pd.DataFrame: """Sweep waveguide width and compute effective index. - Returns pandas dataframe with effective index (neff) and fraction_te. - Args: width1: starting waveguide width in um. width2: end waveguide width in um. steps: number of points. nmodes: number of modes to compute. cmap: colormap for the TE fraction. - Keyword Args: wavelength: (um). wg_width: waveguide width in um. @@ -1039,7 +1010,6 @@ def plot_sweep_width( resolution: pixels/um. Can be a single number or tuple (x, y). nmodes: number of modes to compute. bend_radius: optional bend radius (um). - """ width = np.linspace(width1, width2, steps) neff = {} @@ -1095,9 +1065,9 @@ def plot_sweep_width( nclad=sio2, loss_model=True, sidewall_k = 1E-4, - top_k = 0, + top_k = 1E-4, sidewall_sigma = 20*nm, - top_sigma = 0, + top_sigma = 20*nm, resolution = 400, cache=None, precision='double' From fc2280808ef40c5d6e71cb583bc3ad2310b3f78e Mon Sep 17 00:00:00 2001 From: simbilod Date: Wed, 9 Nov 2022 22:49:27 -0800 Subject: [PATCH 3/3] pre-commit all clear --- gdsfactory/simulation/gtidy3d/modes.py | 102 ++++++++++++++++--------- 1 file changed, 67 insertions(+), 35 deletions(-) diff --git a/gdsfactory/simulation/gtidy3d/modes.py b/gdsfactory/simulation/gtidy3d/modes.py index 485c051914..ec4de4e907 100644 --- a/gdsfactory/simulation/gtidy3d/modes.py +++ b/gdsfactory/simulation/gtidy3d/modes.py @@ -287,25 +287,38 @@ def get_n(self, Y, Z): if self.loss_model: inds_top = ( - (Z >= t_box + wg_thickness - self.top_sigma/2) & (Z <= t_box + wg_thickness + self.top_sigma/2) & (-w / 2 <= Y) & (Y <= w / 2) + (Z >= t_box + wg_thickness - self.top_sigma / 2) + & (Z <= t_box + wg_thickness + self.top_sigma / 2) + & (-w / 2 <= Y) + & (Y <= w / 2) ) inds_top_slab_left = ( - (Z >= t_box + slab_thickness - self.top_sigma/2) & (Z <= t_box + slab_thickness + self.top_sigma/2) & (-w / 2 >= Y) + (Z >= t_box + slab_thickness - self.top_sigma / 2) + & (Z <= t_box + slab_thickness + self.top_sigma / 2) + & (-w / 2 >= Y) ) inds_top_slab_right = ( - (Z >= t_box + slab_thickness - self.top_sigma/2) & (Z <= t_box + slab_thickness + self.top_sigma/2) & (Y >= w / 2) + (Z >= t_box + slab_thickness - self.top_sigma / 2) + & (Z <= t_box + slab_thickness + self.top_sigma / 2) + & (Y >= w / 2) ) inds_sidewall_left = ( - (Z >= t_box + slab_thickness) & (Z <= t_box + wg_thickness) & (Y >= -w / 2 - self.sidewall_sigma/2) & (Y <= -w / 2 + self.sidewall_sigma/2) + (Z >= t_box + slab_thickness) + & (Z <= t_box + wg_thickness) + & (Y >= -w / 2 - self.sidewall_sigma / 2) + & (Y <= -w / 2 + self.sidewall_sigma / 2) ) inds_sidewall_right = ( - (Z >= t_box + slab_thickness) & (Z <= t_box + wg_thickness) & (Y >= w / 2 - self.sidewall_sigma/2) & (Y <= w / 2 + self.sidewall_sigma/2) + (Z >= t_box + slab_thickness) + & (Z <= t_box + wg_thickness) + & (Y >= w / 2 - self.sidewall_sigma / 2) + & (Y <= w / 2 + self.sidewall_sigma / 2) ) - n[inds_top] += 1j*self.top_k - n[inds_top_slab_left] += 1j*self.top_k - n[inds_top_slab_right] += 1j*self.top_k - n[inds_sidewall_left] += 1j*self.sidewall_k - n[inds_sidewall_right] += 1j*self.sidewall_k + n[inds_top] += 1j * self.top_k + n[inds_top_slab_left] += 1j * self.top_k + n[inds_top_slab_right] += 1j * self.top_k + n[inds_sidewall_left] += 1j * self.sidewall_k + n[inds_sidewall_right] += 1j * self.sidewall_k if self.dn_dict is not None: dn = griddata( @@ -320,7 +333,6 @@ def get_n(self, Y, Z): return n - def plot_index(self, func=None) -> None: x, y, Xx, Yx, Xy, Yy, Xz, Yz = create_mesh( -self.w_sim / 2, @@ -339,7 +351,7 @@ def plot_index(self, func=None) -> None: Xx, Yx, ) - if func == None: + if func is None: plot(Xx, Yx, nx) else: plot(Xx, Yx, func(nx)) @@ -640,6 +652,7 @@ def get_overlap( self, wg: "Waveguide", mode_index1: int = 0, mode_index2: int = 0 ) -> float: """Returns mode overlap integral. + Args: wg: other waveguide. """ @@ -657,15 +670,17 @@ def get_overlap( ) def get_loss(self): - """Returns loss for computed modes in dB/cm""" + """Returns loss for computed modes in dB/cm.""" if not hasattr(self, "neffs"): self.compute_modes() wavelength = self.wavelength * 1e-6 # convert to m - alphas = 4 * np.pi * np.imag(self.neffs) / wavelength # lin/m loss - return 10*np.log10(np.exp(1)) * alphas * 1e-2 # dB/cm loss + alphas = 4 * np.pi * np.imag(self.neffs) / wavelength # lin/m loss + return 10 * np.log10(np.exp(1)) * alphas * 1e-2 # dB/cm loss + class WaveguideCoupler(Waveguide): """Waveguide coupler Model. + Parameters: wavelength: (um). wg_width1: left waveguide width in um. @@ -708,6 +723,7 @@ def w_sim(self): def get_n(self, Y, Z): """Return index matrix for a waveguide coupler. + Args: Y: 2D array. Z: 2D array. @@ -763,13 +779,16 @@ def sweep_bend_loss( **kwargs, ) -> Tuple[np.ndarray, np.ndarray]: """Returns overlap integral squared for the bend mode mismatch loss. + The loss is squared because you hit the bend loss twice (from bend to straight and from straight to bend). + Args: bend_radius_min: min bend radius (um). bend_radius_max: max bend radius (um). steps: number of steps. mode_index: where 0 is the fundamental mode. + Keyword Args: wavelength: (um). wg_width: waveguide width in um. @@ -810,11 +829,13 @@ def sweep_neff( **kwargs, ) -> pd.DataFrame: """Sweep waveguide width and compute effective index. + Args: wavelength: (um). thicknesses: in um. widths: in um. mode_index: integer, where 0 is the fundamental mode. + Keyword Args: mode_index: integer. ncore: core refractive index. @@ -853,10 +874,12 @@ def group_index( wavelength: float, wavelength_step: float = 0.01, mode_index: int = 0, **kwargs ) -> float: """Returns group_index. + Args: wavelength: (um). wavelength_step: in um. mode_index: integer, where 0 is the fundamental mode. + Keyword Args: wg_width: waveguide width. wg_thickness: thickness waveguide (um). @@ -896,10 +919,12 @@ def sweep_group_index( **kwargs, ) -> pd.DataFrame: """Sweep waveguide width and compute group index. + Args: wavelength: (um). thicknesses: in um. widths: in um. + Keyword Args: mode_index: integer. ncore: core refractive index. @@ -939,12 +964,15 @@ def sweep_width( **kwargs, ) -> pd.DataFrame: """Sweep waveguide width and compute effective index. + Returns pandas dataframe with effective index (neff) and fraction_te. + Args: width1: starting waveguide width in um. width2: end waveguide width in um. steps: number of points. nmodes: number of modes to compute. + Keyword Args: wavelength: (um). wg_width: waveguide width in um. @@ -990,13 +1018,16 @@ def plot_sweep_width( **kwargs, ) -> pd.DataFrame: """Sweep waveguide width and compute effective index. + Returns pandas dataframe with effective index (neff) and fraction_te. + Args: width1: starting waveguide width in um. width2: end waveguide width in um. steps: number of points. nmodes: number of modes to compute. cmap: colormap for the TE fraction. + Keyword Args: wavelength: (um). wg_width: waveguide width in um. @@ -1056,22 +1087,23 @@ def plot_sweep_width( if __name__ == "__main__": - wg = Waveguide(nmodes=2, - wg_width=500*nm, - wavelength=1.55, - wg_thickness=220 * nm, - slab_thickness=90 * nm, - ncore=si, - nclad=sio2, - loss_model=True, - sidewall_k = 1E-4, - top_k = 1E-4, - sidewall_sigma = 20*nm, - top_sigma = 20*nm, - resolution = 400, - cache=None, - precision='double' - ) + wg = Waveguide( + nmodes=2, + wg_width=500 * nm, + wavelength=1.55, + wg_thickness=220 * nm, + slab_thickness=90 * nm, + ncore=si, + nclad=sio2, + loss_model=True, + sidewall_k=1e-4, + top_k=1e-4, + sidewall_sigma=20 * nm, + top_sigma=20 * nm, + resolution=400, + cache=None, + precision="double", + ) wg.plot_index() wg.plot_index(func=np.imag) wg.compute_mode_properties() @@ -1083,7 +1115,7 @@ def plot_sweep_width( # n = np.real(wg.neffs[mode_number]) # fraction_te = wg.fraction_te[mode_number] # plt.scatter(wg_width, n, c=fraction_te, vmin=0, vmax=1, cmap=cmap) - # n[inds_top_slab_left] += self.top_k - # n[inds_top_slab_right] += self.top_k - # n[inds_sidewall_left] += self.sidewall_k - # n[inds_sidewall_right] += self.sidewall_k \ No newline at end of file + # n[inds_top_slab_left] += self.top_k + # n[inds_top_slab_right] += self.top_k + # n[inds_sidewall_left] += self.sidewall_k + # n[inds_sidewall_right] += self.sidewall_k