From a9381115ec05eb357a6e572aea7283bb61a415f1 Mon Sep 17 00:00:00 2001 From: momchil Date: Thu, 27 Jan 2022 17:37:13 -0800 Subject: [PATCH] Bringing in upgraded mode solver Latest solver updates from old tidy3d bend modes in ModeSpec Mode solver always takes grid boundaries instead of grid_size so that it works with nonuniform coords Changed structure close to PML warning a bit --- .pylintrc | 2 +- lint.py | 5 +- test_static.sh | 11 +- tests/test_plugins.py | 9 +- tidy3d/components/mode.py | 27 +- tidy3d/components/monitor.py | 8 +- tidy3d/components/simulation.py | 18 +- tidy3d/components/types.py | 1 + tidy3d/plugins/mode/mode_solver.py | 29 +-- tidy3d/plugins/mode/solver.py | 391 +++++++++++++++++++++-------- tidy3d/plugins/mode/transforms.py | 114 +++++++++ 11 files changed, 469 insertions(+), 146 deletions(-) create mode 100644 tidy3d/plugins/mode/transforms.py diff --git a/.pylintrc b/.pylintrc index 5d99d0c2b..9d7b17120 100644 --- a/.pylintrc +++ b/.pylintrc @@ -1,7 +1,7 @@ [MASTER] extension-pkg-allow-list=pydantic ignore=material_library.py, plugins -good-names=ax, im, Lx, Ly, Lz, x0, y0, z0, x, y, z, f, t, y1, y2, x1, x2, xs, ys, zs, Ax, Nx, Ny, Nz, dl, rr, E, H, xx, yy, zz, dx, dy, Jx, Jy, Jz, Ex, Ey, Ez, Mx, My, Mz, Hx, Hy, Hz, dz, e, fp, dt, a, c, k0, Jx_k, Jy_k, My_k, Mx_k, N_theta, N_phi, L_theta, L_phi +good-names=ax, im, Lx, Ly, Lz, x0, y0, z0, x, y, z, u, v, w, f, t, y1, y2, x1, x2, xs, ys, zs, Ax, Nx, Ny, Nz, N, dl, rr, E, H, xx, yy, zz, dx, dy, Jx, Jy, Jz, Ex, Ey, Ez, Mx, My, Mz, Hx, Hy, Hz, dz, e, fp, dt, a, c, kx, ky, kz, k0, Jx_k, Jy_k, My_k, Mx_k, N_theta, N_phi, L_theta, L_phi [BASIC] diff --git a/lint.py b/lint.py index 4e7aa628c..5df0fb93e 100644 --- a/lint.py +++ b/lint.py @@ -35,7 +35,10 @@ def main(): results = Run([path], do_exit=False) - final_score = results.linter.stats.global_note + try: + final_score = results.linter.stats.global_note + except AttributeError: + final_score = results.linter.stats["global_note"] if final_score < threshold: diff --git a/test_static.sh b/test_static.sh index 157f71122..e0749f0da 100644 --- a/test_static.sh +++ b/test_static.sh @@ -2,11 +2,10 @@ black . python lint.py -pytest -rA tests/test_components.py -pytest -rA tests/test_grid.py -pytest -rA tests/test_IO.py -pytest -rA tests/test_material_library.py -# pytest -rA tests/test_core.py -pytest -rA tests/test_plugins.py +pytest -ra tests/test_components.py +pytest -ra tests/test_grid.py +pytest -ra tests/test_IO.py +pytest -ra tests/test_material_library.py +pytest -ra tests/test_plugins.py pytest --doctest-modules tidy3d/components --ignore=tidy3d/components/base.py \ No newline at end of file diff --git a/tests/test_plugins.py b/tests/test_plugins.py index b166f0f11..0b55c58ad 100644 --- a/tests/test_plugins.py +++ b/tests/test_plugins.py @@ -46,7 +46,14 @@ def test_mode_solver(): simulation = td.Simulation(size=(2, 2, 2), grid_size=(0.1, 0.1, 0.1), structures=[waveguide]) plane = td.Box(center=(0, 0, 0), size=(0, 1, 1)) ms = ModeSolver(simulation=simulation, plane=plane, freq=td.constants.C_0 / 1.5) - modes = ms.solve(mode_spec=td.ModeSpec(num_modes=2)) + mode_spec = td.ModeSpec( + num_modes=3, + target_neff=2.0, + bend_radius=3.0, + bend_axis=0, + num_pml=(10, 10), + ) + modes = ms.solve(mode_spec=mode_spec) def _test_coeffs(): diff --git a/tidy3d/components/mode.py b/tidy3d/components/mode.py index c9176e6ad..2412abc6a 100644 --- a/tidy3d/components/mode.py +++ b/tidy3d/components/mode.py @@ -4,8 +4,10 @@ import pydantic as pd +from ..constants import MICROMETER from .base import Tidy3dBaseModel -from .types import Symmetry +from .types import Symmetry, Axis2D +from ..log import SetupError class ModeSpec(Tidy3dBaseModel): @@ -38,3 +40,26 @@ class ModeSpec(Tidy3dBaseModel): title="Number of PML layers", description="Number of standard pml layers to add in the first two non-propagation axes.", ) + + bend_radius: float = pd.Field( + None, + title="Bend radius", + description="A curvature radius for simulation of waveguide bends. Can be negative, in " + "which case the mode plane center has a smaller value than the curvature center along the " + "axis that is perpendicular to both the normal axis and the bend axis.", + units=MICROMETER, + ) + + bend_axis: Axis2D = pd.Field( + None, + title="Bend axis", + description="Index into the first two non-propagating axes defining the normal to the " + "plane in which the bend lies. This must be provided if ``bend_radius`` is not ``None``.", + ) + + @pd.validator("bend_axis", always=True) + def bend_axis_given(cls, val, values): + """check that ``bend_axis`` is provided if ``bend_radius`` is not ``None``""" + if val is None and values.get("bend_radius") is not None: + raise SetupError("bend_axis must also be defined if bend_radius is defined.") + return val diff --git a/tidy3d/components/monitor.py b/tidy3d/components/monitor.py index 7955aa366..f318fa0af 100644 --- a/tidy3d/components/monitor.py +++ b/tidy3d/components/monitor.py @@ -4,7 +4,7 @@ import pydantic -from .types import Literal, Ax, Direction, EMField, ArrayLike +from .types import Literal, Ax, EMField, ArrayLike from .geometry import Box from .validators import assert_plane from .mode import ModeSpec @@ -265,12 +265,6 @@ class ModeMonitor(PlanarMonitor, FreqMonitor): ... name='mode_monitor') """ - direction: List[Direction] = pydantic.Field( - ["+", "-"], - title="Direction", - description="Specifies which direction(s) of mode propagation for the monitor to measure.", - ) - mode_spec: ModeSpec = pydantic.Field( ..., title="Mode Specification", diff --git a/tidy3d/components/simulation.py b/tidy3d/components/simulation.py index bfbbaadb9..174e3c5d9 100644 --- a/tidy3d/components/simulation.py +++ b/tidy3d/components/simulation.py @@ -265,16 +265,16 @@ def _structures_not_close_pml(cls, val, values): # pylint:disable=too-many-loca if (not structures) or (not sources): return val - def warn(structure): + def warn(istruct, side): """Warning message for a structure too close to PML.""" log.warning( - f"a structure\n\n{structure}\n\nwas detected as being less " - "than half of a central wavelength from a PML on - side" + f"Structure at structures[{istruct}] was detected as being less " + f"than half of a central wavelength from a PML on side {side}. " "To avoid inaccurate results, please increase gap between " "any structures and PML or fully extend structure through the pml." ) - for structure in structures: + for istruct, structure in enumerate(structures): struct_bound_min, struct_bound_max = structure.geometry.bounds for source in sources: @@ -282,15 +282,17 @@ def warn(structure): f_average = (fmin_src + fmax_src) / 2.0 lambda0 = C_0 / f_average - for sim_val, struct_val, pml in zip(sim_bound_min, struct_bound_min, val): + zipped = zip(["x", "y", "z"], sim_bound_min, struct_bound_min, val) + for axis, sim_val, struct_val, pml in zipped: if pml.num_layers > 0 and struct_val > sim_val: if abs(sim_val - struct_val) < lambda0 / 2: - warn(structure) + warn(istruct, axis + "-min") - for sim_val, struct_val, pml in zip(sim_bound_max, struct_bound_max, val): + zipped = zip(["x", "y", "z"], sim_bound_max, struct_bound_max, val) + for axis, sim_val, struct_val, pml in zipped: if pml.num_layers > 0 and struct_val < sim_val: if abs(sim_val - struct_val) < lambda0 / 2: - warn(structure) + warn(istruct, axis + "-max") return val diff --git a/tidy3d/components/types.py b/tidy3d/components/types.py index 95784439c..e3441a4eb 100644 --- a/tidy3d/components/types.py +++ b/tidy3d/components/types.py @@ -97,6 +97,7 @@ def __modify_schema__(cls, field_schema): Bound = Tuple[Coordinate, Coordinate] GridSize = Union[pydantic.PositiveFloat, List[pydantic.PositiveFloat]] Axis = Literal[0, 1, 2] +Axis2D = Literal[0, 1] Vertices = List[Coordinate2D] Shapely = BaseGeometry diff --git a/tidy3d/plugins/mode/mode_solver.py b/tidy3d/plugins/mode/mode_solver.py index d85507755..f0f03ba86 100644 --- a/tidy3d/plugins/mode/mode_solver.py +++ b/tidy3d/plugins/mode/mode_solver.py @@ -14,9 +14,7 @@ from ...components import ModeMonitor from ...components import ModeSource, GaussianPulse from ...components.types import Direction -from ...components.data import ScalarFieldData, FieldData, Tidy3dData -from ...log import SetupError -from ...constants import C_0 +from ...components.data import ScalarFieldData, FieldData from .solver import compute_modes @@ -56,9 +54,9 @@ class ModeInfo(Tidy3dBaseModel): """stores information about a (solved) mode. Attributes ---------- - field_data: xr.Dataset + field_data: FieldData Contains information about the fields of the modal profile. - mode: Mode + mode_spec: ModeSpec Specifications of the mode. n_eff: float Real part of the effective refractive index of mode. @@ -100,13 +98,13 @@ def solve(self, mode_spec: ModeSpec) -> List[ModeInfo]: Parameters ---------- - mode : Mode - ``Mode`` object containing specifications of mode. + mode_spec : ModeSpec + ``ModeSpec`` object containing specifications of the mode solver. Returns ------- - ModeInfo - Object containing mode profile and effective index data. + List[ModeInfo] + A list of ``ModeInfo`` objects for each mode. """ normal_axis = self.plane.size.index(0.0) @@ -124,6 +122,10 @@ def solve(self, mode_spec: ModeSpec) -> List[ModeInfo]: (eps_xx, eps_yy, eps_zz), axis=normal_axis ) + # get the in-plane grid coordinates on which eps and the mode fields live + plane_grid = self.simulation.discretize(self.plane) + _, plane_coords = self.plane.pop_axis(plane_grid.boundaries.to_list, axis=normal_axis) + # note: from this point on, in waveguide coordinates (propagating in z) # construct eps_cross section to feed to mode solver @@ -135,16 +137,11 @@ def solve(self, mode_spec: ModeSpec) -> List[ModeInfo]: # if mode_spec.symmetries[1] != 0: # eps_cross = np.stack(tuple(e[:, Ny // 2] for e in eps_cross)) - # note, internally discretizing, need to make consistent. mode_fields, n_eff_complex = compute_modes( eps_cross=eps_cross, + coords=plane_coords, freq=self.freq, - grid_size=self.simulation.grid_size, - pml_layers=mode_spec.num_pml, - num_modes=mode_spec.num_modes, - target_neff=mode_spec.target_neff, - symmetries=mode_spec.symmetries, - coords=None, + mode_spec=mode_spec, ) def rotate_field_coords(e_field, h_field): diff --git a/tidy3d/plugins/mode/solver.py b/tidy3d/plugins/mode/solver.py index cf6f5cd5c..958182576 100644 --- a/tidy3d/plugins/mode/solver.py +++ b/tidy3d/plugins/mode/solver.py @@ -1,3 +1,4 @@ +"""Mode solver for propagating EM modes.""" from typing import Tuple import numpy as np @@ -5,71 +6,45 @@ import scipy.sparse.linalg as spl from ...components.types import Numpy -from ...constants import EPSILON_0, ETA_0, C_0, MU_0, fp_eps, pec_val +from ...constants import ETA_0, C_0, fp_eps, pec_val from .derivatives import create_D_matrices as D_mats from .derivatives import create_S_matrices as S_mats - -# from .Mode import Mode +from .transforms import radial_transform, angled_transform def compute_modes( eps_cross, + coords, freq, - grid_size, - pml_layers, - num_modes=1, - target_neff=None, - symmetries=(0, 0), - coords=None, + mode_spec, + angle_theta=0.0, + angle_phi=0.0, ) -> Tuple[Numpy, Numpy]: """Solve for the modes of a waveguide cross section. Parameters ---------- eps_cross : array_like or tuple of array_like - Either a single 2D array defining the relative permittivity in the - cross-section, or three 2D arrays defining the permittivity at the Ex, - Ey, and Ez locations of the Yee cell, respectively. + Either a single 2D array defining the relative permittivity in the cross-section, or three + 2D arrays defining the permittivity at the Ex, Ey, and Ez locations of the Yee grid. + coords : List[Numpy] + Two 1D arrays with each with size one larger than the corresponding axis of ``eps_cross``. + Defines a (potentially non-uniform) Cartesian grid on which the modes are computed. freq : float (Hertz) Frequency at which the eigenmodes are computed. - grid_size : list or tuple of float - (micron) Step size in x, y and z. The mesh step in z is currently - unused, but it could be needed if numerical dispersion is to be taken - into account. - pml_layers : list or tuple of int - Number of pml layers in x and y. - num_modes : int, optional - Number of modes to be computed. - target_neff : None or float, optional - Look for modes closest to target_neff. If ``None``, modes with the - largest effective index are returned. - symmetries : array_like, optional - Array of two integers defining reflection symmetry to be applied - at the xmin and the ymin locations. Note then that this assumes that - ``eps_cross`` is only supplied in the quadrants in which it is needed - and *not* in their symmetric counterparts. Each element can be ``0`` - (no symmetry), ``1`` (even, i.e. 'PMC' symmetry) or ``-1`` (odd, i.e. - 'PEC' symmetry). - coords : List of array_like or None, optional - If provided, overrides ``grid_size``, and must be a list of two arrays - with size one larger than the corresponding axis of ``eps_cross`. - Defines a non-uniform Cartesian grid on which the modes are computed. + mode_spec : ModeSpec + ``ModeSpec`` object containing specifications of the mode solver. Returns ------- - list of dict - A list of all the computed modes. Each entry is a dictionary with the - real and imaginary part of the effective index of the waveguide and - all the E and H field components. - - Raises - ------ - RuntimeError - Description - ValueError - Description + Tuple[Numpy, Numpy] + The first array gives the E and H fields for all modes, the second one give sthe complex + effective index. """ + num_modes = mode_spec.num_modes + bend_radius = mode_spec.bend_radius + bend_axis = mode_spec.bend_axis omega = 2 * np.pi * freq k0 = omega / C_0 @@ -81,45 +56,86 @@ def compute_modes( else: raise ValueError except Exception as e: - printf("Wrong input to mode solver pemittivity!") - raise (e) + print("Wrong input to mode solver pemittivity!") + raise e Nx, Ny = eps_xx.shape N = eps_xx.size - if coords is None: - coords_x = grid_size[0] * np.arange(Nx + 1) - coords_y = grid_size[1] * np.arange(Ny + 1) - coords = [coords_x, coords_y] - else: - if coords[0].size != Nx + 1 or coords[1].size != Ny + 1: - raise ValueError("Mismatch between 'coords' and 'esp_cross' shapes.") - - """ The forward derivative matrices already impose PEC boundary at the - xmax and ymax interfaces. Here, we also impose PEC boundaries on the - xmin and ymin interfaces through the permittivity at those positions, - unless a PMC symmetry is specifically requested. The PMC symmetry is + if coords[0].size != Nx + 1 or coords[1].size != Ny + 1: + raise ValueError("Mismatch between 'coords' and 'esp_cross' shapes.") + new_coords = [np.copy(c) for c in coords] + + """We work with full tensorial epsilon in mu to handle the most general cases that can + be introduced by coordinate transformations. In the solver, we distinguish the case when + these tensors are still diagonal, in which case the matrix for diagonalization has shape + (2N, 2N), and the full tensorial case, in which case it has shape (4N, 4N).""" + eps_tensor = np.zeros((3, 3, N), dtype=np.complex128) + mu_tensor = np.zeros((3, 3, N), dtype=np.complex128) + for dim, eps in enumerate([eps_xx, eps_yy, eps_zz]): + eps_tensor[dim, dim, :] = eps.ravel() + mu_tensor[dim, dim, :] = 1.0 + + # Get Jacobian of all coordinate transformations. Initialize as identity (same as mu so far) + jac_e = np.copy(mu_tensor) + jac_h = np.copy(mu_tensor) + + if bend_radius is not None: + new_coords, jac_e, jac_h = radial_transform(new_coords, bend_radius, bend_axis) + + if angle_theta > 0: + new_coords, jac_e_tmp, jac_h_tmp = angled_transform(new_coords, angle_theta, angle_phi) + jac_e = np.einsum("ij...,jp...->ip...", jac_e_tmp, jac_e) + jac_h = np.einsum("ij...,jp...->ip...", jac_h_tmp, jac_h) + + """We also need to keep track of the transformation of the k-vector. This is + the eigenvalue of the momentum operator assuming some sort of translational invariance and is + different from just the transformation of the derivative operator. For example, in a bent + waveguide, there is strictly speaking no k-vector in the original coordinates as the system + is not translationally invariant there. However, if we define kz = R k_phi, then the + effective index approaches that for a straight-waveguide in the limit of infinite radius. + Since we use w = R phi in the radial_transform, there is nothing else neede in the k transform. + For the angled_transform, the transformation between k-vectors follows from writing the field as + E' exp(i k_p w) in transformed coordinates, and identifying this with + E exp(i k_x x + i k_y y + i k_z z) in the original ones.""" + kxy = np.cos(angle_theta) ** 2 + kz = np.cos(angle_theta) * np.sin(angle_theta) + kp_to_k = np.array([kxy * np.sin(angle_phi), kxy * np.cos(angle_phi), kz]) + + # Transform epsilon and mu + jac_e_det = np.linalg.det(np.moveaxis(jac_e, [0, 1], [-2, -1])) + jac_h_det = np.linalg.det(np.moveaxis(jac_h, [0, 1], [-2, -1])) + eps_tensor = np.einsum("ij...,jp...->ip...", jac_e, eps_tensor) # J.dot(eps) + eps_tensor = np.einsum("ij...,pj...->ip...", eps_tensor, jac_e) # (J.dot(eps)).dot(J.T) + eps_tensor /= jac_e_det + mu_tensor = np.einsum("ij...,jp...->ip...", jac_h, mu_tensor) + mu_tensor = np.einsum("ij...,pj...->ip...", mu_tensor, jac_h) + mu_tensor /= jac_h_det + + """ The forward derivative matrices already impose PEC boundary at the xmax and ymax interfaces. + Here, we also impose PEC boundaries on the xmin and ymin interfaces through the permittivity at + those positions, unless a PMC symmetry is specifically requested. The PMC symmetry is imposed by modifying the backward derivative matrices.""" dmin_pmc = [False, False] - if symmetries[0] != 1: + if mode_spec.symmetries[0] != 1: # PEC at the xmin edge - eps_yy[0, :] = pec_val - eps_zz[0, :] = pec_val + eps_tensor[1, 1, :Ny] = pec_val + eps_tensor[2, 2, :Ny] = pec_val else: # Modify the backwards x derivative dmin_pmc[0] = True if Ny > 1: - if symmetries[1] != 1: + if mode_spec.symmetries[1] != 1: # PEC at the ymin edge - eps_xx[:, 0] = pec_val - eps_zz[:, 0] = pec_val + eps_tensor[0, 0, ::Ny] = pec_val + eps_tensor[2, 2, ::Ny] = pec_val else: # Modify the backwards y derivative dmin_pmc[1] = True # Primal grid steps for E-field derivatives - dLf = [c[1:] - c[:-1] for c in coords] + dLf = [c[1:] - c[:-1] for c in new_coords] # Dual grid steps for H-field derivatives dLtmp = [(dL[:-1] + dL[1:]) / 2 for dL in dLf] dLb = [np.hstack((d1[0], d2)) for d1, d2 in zip(dLf, dLtmp)] @@ -128,37 +144,109 @@ def compute_modes( Dmats = D_mats((Nx, Ny), dLf, dLb, dmin_pmc) # PML matrices; do not impose PML on the bottom when symmetry present - dmin_pml = np.array(symmetries) == 0 - Smats = S_mats(omega, (Nx, Ny), pml_layers, dLf, dLb, dmin_pml) + dmin_pml = np.array(mode_spec.symmetries) == 0 + Smats = S_mats(omega, (Nx, Ny), mode_spec.num_pml, dLf, dLb, dmin_pml) + + # Add the PML on top of the derivatives; normalize by k0 to match the EM-possible notation + SDmats = [Smat.dot(Dmat) / k0 for Smat, Dmat in zip(Smats, Dmats)] + + # Determine initial guess value for the solver in transformed coordinates + if mode_spec.target_neff is None: + eps_physical = np.array(eps_cross) + eps_physical = eps_physical[np.abs(eps_physical) < np.abs(pec_val)] + n_max = np.sqrt(np.max(np.abs(eps_physical))) + target = n_max + else: + target = mode_spec.target_neff + target_neff_p = target / np.linalg.norm(kp_to_k) + + # Solve for the modes + E, H, neff, keff = solver_em(eps_tensor, mu_tensor, SDmats, num_modes, target_neff_p) + + # Transform back to original axes, E = J^T E' + E = np.sum(jac_e[..., None] * E[:, None, ...], axis=0) + E = E.reshape(3, Nx, Ny, 1, num_modes) + H = np.sum(jac_h[..., None] * H[:, None, ...], axis=0) + H = H.reshape(3, Nx, Ny, 1, num_modes) + neff = neff * np.linalg.norm(kp_to_k) + + F = np.stack((E, H), axis=0) + + return F, neff + 1j * keff - # Add the PML on top of the derivatives - SDmats = [Smat.dot(Dmat) for Smat, Dmat in zip(Smats, Dmats)] - # Normalize by k0 to match the EM-possible notation - Dxf, Dxb, Dyf, Dyb = [mat / k0 for mat in SDmats] +def solver_em(eps_tensor, mu_tensor, SDmats, num_modes, neff_guess): + """Solve for the electromagnetic modes of a system defined by in-plane permittivity and + permeability and assuming translational invariance in the normal direction. + + Parameters + ---------- + eps_tensor : np.ndarray + Shape (3, 3, N), the permittivity tensor at every point in the plane. + mu_tensor : np.ndarray + Shape (3, 3, N), the permittivity tensor at every point in the plane. + SDmats : List[scipy.sparse.csr_matrix] + The sparce derivative matrices Dxf, Dxb, Dyf, Dyb, including the PML. + num_modes : int + Number of modes to solve for. + neff_guess : float + Initial guess for the effective index. + + Returns + ------- + E : np.ndarray + Electric field of the eigenmodes, shape (3, N, num_modes). + H : np.ndarray + Magnetic field of the eigenmodes, shape (3, N, num_modes). + neff : np.ndarray + Real part of the effective index, shape (num_modes, ). + keff : np.ndarray + Imaginary part of the effective index, shape (num_modes, ). + """ + + off_diagonals = (np.ones((3, 3)) - np.eye(3)).astype(bool) + eps_offd = np.abs(eps_tensor[off_diagonals]) + mu_offd = np.abs(mu_tensor[off_diagonals]) + if np.any(eps_offd > 1e-6) or np.any(mu_offd > 1e-6): + return solver_tensorial(eps_tensor, mu_tensor, SDmats, num_modes, neff_guess) + else: + return solver_diagonal(eps_tensor, mu_tensor, SDmats, num_modes, neff_guess) - # Compute matrix for diagonalization - inv_eps_zz = sp.spdiags(1 / eps_zz.flatten(), [0], N, N) + +def solver_diagonal(eps, mu, SDmats, num_modes, neff_guess): + """EM eigenmode solver assuming ``eps`` and ``mu`` are diagonal everywhere.""" + + N = eps.shape[-1] + + # Unpack eps, mu and derivatives + eps_xx = eps[0, 0, :] + eps_yy = eps[1, 1, :] + eps_zz = eps[2, 2, :] + mu_xx = mu[0, 0, :] + mu_yy = mu[1, 1, :] + mu_zz = mu[2, 2, :] + Dxf, Dxb, Dyf, Dyb = SDmats + + # Compute the matrix for diagonalization + inv_eps_zz = sp.spdiags(1 / eps_zz, [0], N, N) + inv_mu_zz = sp.spdiags(1 / mu_zz, [0], N, N) P11 = -Dxf.dot(inv_eps_zz).dot(Dyb) - P12 = Dxf.dot(inv_eps_zz).dot(Dxb) + sp.eye(N) - P21 = -Dyf.dot(inv_eps_zz).dot(Dyb) - sp.eye(N) + P12 = Dxf.dot(inv_eps_zz).dot(Dxb) + sp.spdiags(mu_yy, [0], N, N) + P21 = -Dyf.dot(inv_eps_zz).dot(Dyb) - sp.spdiags(mu_xx, [0], N, N) P22 = Dyf.dot(inv_eps_zz).dot(Dxb) - Q11 = -Dxb.dot(Dyf) - Q12 = Dxb.dot(Dxf) + sp.spdiags(eps_yy.flatten(), [0], N, N) - Q21 = -Dyb.dot(Dyf) - sp.spdiags(eps_xx.flatten(), [0], N, N) - Q22 = Dyb.dot(Dxf) + Q11 = -Dxb.dot(inv_mu_zz).dot(Dyf) + Q12 = Dxb.dot(inv_mu_zz).dot(Dxf) + sp.spdiags(eps_yy, [0], N, N) + Q21 = -Dyb.dot(inv_mu_zz).dot(Dyf) - sp.spdiags(eps_xx, [0], N, N) + Q22 = Dyb.dot(inv_mu_zz).dot(Dxf) Pmat = sp.bmat([[P11, P12], [P21, P22]]) Qmat = sp.bmat([[Q11, Q12], [Q21, Q22]]) A = Pmat.dot(Qmat) - if target_neff is None: - n_max = np.sqrt(np.max(eps_cross)) - guess_value = -(n_max ** 2) - else: - guess_value = -(target_neff ** 2) - - vals, vecs = solver_eigs(A, num_modes, guess_value=guess_value) + # Call the eigensolver. The eigenvalues are -(neff + 1j * keff)**2 + vals, vecs = solver_eigs(A, num_modes, guess_value=-(neff_guess ** 2)) + if vals.size == 0: + raise RuntimeError("Could not find any eigenmodes for this waveguide") vre, vim = -np.real(vals), -np.imag(vals) # Sort by descending real part @@ -168,38 +256,131 @@ def compute_modes( vecs = vecs[:, sort_inds] # Real and imaginary part of the effective index - neff_tmp = np.sqrt(vre / 2 + np.sqrt(vre ** 2 + vim ** 2) / 2) - keff = vim / 2 / (neff_tmp + 1e-10) - - # Correct formula taking numerical dispersion into account - # neff = 2/grid_size[2]*np.arcsin((neff_tmp + 1j*keff)*grid_size[2]/2) - neff = neff_tmp + neff = np.sqrt(vre / 2 + np.sqrt(vre ** 2 + vim ** 2) / 2) + keff = vim / 2 / (neff + 1e-10) # Field components from eigenvectors Ex = vecs[:N, :] Ey = vecs[N:, :] - # Get the other field components; normalize according to CEM - Hs = -Qmat.dot(vecs) / (neff + 1j * keff)[np.newaxis, :] / ETA_0 - Hx = Hs[:N, :] - Hy = Hs[N:, :] - - Hz = Dxf.dot(Ey) - Dyf.dot(Ex) + # Get the other field components + Hs = Qmat.dot(vecs) + Hx = Hs[:N, :] / (1j * neff - keff) + Hy = Hs[N:, :] / (1j * neff - keff) + Hz = inv_mu_zz.dot((Dxf.dot(Ey) - Dyf.dot(Ex))) Ez = inv_eps_zz.dot((Dxb.dot(Hy) - Dyb.dot(Hx))) - # Store all the information about the modes. + # Bundle up + E = np.stack((Ex, Ey, Ez), axis=0) + H = np.stack((Hx, Hy, Hz), axis=0) - def reshape(array): - return array.reshape(Nx, Ny, 1, num_modes) + # Return to standard H field units (see CEM notes for H normalization used in solver) + H *= -1j / ETA_0 + + return E, H, neff, keff + + +def solver_tensorial(eps, mu, SDmats, num_modes, neff_guess): + """EM eigenmode solver assuming ``eps`` or ``mu`` have off-diagonal elements.""" + + N = eps.shape[-1] + Dxf, Dxb, Dyf, Dyb = SDmats + + # Compute all blocks of the matrix for diagonalization + inv_eps_zz = sp.spdiags(1 / eps[2, 2, :], [0], N, N) + inv_mu_zz = sp.spdiags(1 / mu[2, 2, :], [0], N, N) + axax = -Dxf.dot(sp.spdiags(eps[2, 0, :] / eps[2, 2, :], [0], N, N)) - sp.spdiags( + mu[1, 2, :] / mu[2, 2, :], [0], N, N + ).dot(Dyf) + axay = -Dxf.dot(sp.spdiags(eps[2, 1, :] / eps[2, 2, :], [0], N, N)) + sp.spdiags( + mu[1, 2, :] / mu[2, 2, :], [0], N, N + ).dot(Dxf) + axbx = -Dxf.dot(inv_eps_zz).dot(Dyb) + sp.spdiags( + mu[1, 0, :] - mu[1, 2, :] * mu[2, 0, :] / mu[2, 2, :], [0], N, N + ) + axby = Dxf.dot(inv_eps_zz).dot(Dxb) + sp.spdiags( + mu[1, 1, :] - mu[1, 2, :] * mu[2, 1, :] / mu[2, 2, :], [0], N, N + ) + ayax = -Dyf.dot(sp.spdiags(eps[2, 0, :] / eps[2, 2, :], [0], N, N)) + sp.spdiags( + mu[0, 2, :] / mu[2, 2, :], [0], N, N + ).dot(Dyf) + ayay = -Dyf.dot(sp.spdiags(eps[2, 1, :] / eps[2, 2, :], [0], N, N)) - sp.spdiags( + mu[0, 2, :] / mu[2, 2, :], [0], N, N + ).dot(Dxf) + aybx = -Dyf.dot(inv_eps_zz).dot(Dyb) + sp.spdiags( + -mu[0, 0, :] + mu[0, 2, :] * mu[2, 0, :] / mu[2, 2, :], [0], N, N + ) + ayby = Dyf.dot(inv_eps_zz).dot(Dxb) + sp.spdiags( + -mu[0, 1, :] + mu[0, 2, :] * mu[2, 1, :] / mu[2, 2, :], [0], N, N + ) + bxbx = -Dxb.dot(sp.spdiags(mu[2, 0, :] / mu[2, 2, :], [0], N, N)) - sp.spdiags( + eps[1, 2, :] / eps[2, 2, :], [0], N, N + ).dot(Dyb) + bxby = -Dxb.dot(sp.spdiags(mu[2, 1, :] / mu[2, 2, :], [0], N, N)) + sp.spdiags( + eps[1, 2, :] / eps[2, 2, :], [0], N, N + ).dot(Dxb) + bxax = -Dxb.dot(inv_mu_zz).dot(Dyf) + sp.spdiags( + eps[1, 0, :] - eps[1, 2, :] * eps[2, 0, :] / eps[2, 2, :], [0], N, N + ) + bxay = Dxb.dot(inv_mu_zz).dot(Dxf) + sp.spdiags( + eps[1, 1, :] - eps[1, 2, :] * eps[2, 1, :] / eps[2, 2, :], [0], N, N + ) + bybx = -Dyb.dot(sp.spdiags(mu[2, 0, :] / mu[2, 2, :], [0], N, N)) + sp.spdiags( + eps[0, 2, :] / eps[2, 2, :], [0], N, N + ).dot(Dyb) + byby = -Dyb.dot(sp.spdiags(mu[2, 1, :] / mu[2, 2, :], [0], N, N)) - sp.spdiags( + eps[0, 2, :] / eps[2, 2, :], [0], N, N + ).dot(Dxb) + byax = -Dyb.dot(inv_mu_zz).dot(Dyf) + sp.spdiags( + -eps[0, 0, :] + eps[0, 2, :] * eps[2, 0, :] / eps[2, 2, :], [0], N, N + ) + byay = Dyb.dot(inv_mu_zz).dot(Dxf) + sp.spdiags( + -eps[0, 1, :] + eps[0, 2, :] * eps[2, 1, :] / eps[2, 2, :], [0], N, N + ) + + A = sp.bmat( + [ + [axax, axay, axbx, axby], + [ayax, ayay, aybx, ayby], + [bxax, bxay, bxbx, bxby], + [byax, byay, bybx, byby], + ] + ) + + # Call the eigensolver. The eigenvalues are 1j * (neff + 1j * keff) + vals, vecs = solver_eigs(A, num_modes, guess_value=1j * neff_guess) + if vals.size == 0: + raise RuntimeError("Could not find any eigenmodes for this waveguide") + # Real and imaginary part of the effective index + neff, keff = np.imag(vals), -np.real(vals) - (Ex, Ey, Ez, Hx, Hy, Hz) = tuple(map(reshape, (Ex, Ey, Ez, Hx, Hy, Hz))) + # Sort by descending real part + sort_inds = np.argsort(neff)[::-1] + neff = neff[sort_inds] + keff = keff[sort_inds] + vecs = vecs[:, sort_inds] + # Field components from eigenvectors + Ex = vecs[:N, :] + Ey = vecs[N : 2 * N, :] + Hx = vecs[2 * N : 3 * N, :] + Hy = vecs[3 * N :, :] + + # Get the other field components + hxy_term = (-mu[2, 0, :] * Hx.T - mu[2, 1, :] * Hy.T).T + Hz = inv_mu_zz.dot(Dxf.dot(Ey) - Dyf.dot(Ex) + hxy_term) + exy_term = (-eps[2, 0, :] * Ex.T - eps[2, 1, :] * Ey.T).T + Ez = inv_eps_zz.dot(Dxb.dot(Hy) - Dyb.dot(Hx) + exy_term) + + # Bundle up E = np.stack((Ex, Ey, Ez), axis=0) H = np.stack((Hx, Hy, Hz), axis=0) - F = np.stack((E, H), axis=0) + # Return to standard H field units (see CEM notes for H normalization used in solver) + # The minus sign here is suspicious, need to check how modes are used in Mode objects + H *= -1j / ETA_0 - return F, neff + 1j * keff + return E, H, neff, keff def solver_eigs(A, num_modes, guess_value=1.0): diff --git a/tidy3d/plugins/mode/transforms.py b/tidy3d/plugins/mode/transforms.py new file mode 100644 index 000000000..c3403e7c6 --- /dev/null +++ b/tidy3d/plugins/mode/transforms.py @@ -0,0 +1,114 @@ +""" Coordinate transformations. + +The Jacobian of a transformation from coordinates r = (x, y, z) into coordinates +r' = (u, v, w) is defined as J_ij = dr'_i/dr_j. Here, z and w are the propagation axes in the +original and transformed planes, respectively, and the coords are only provided in (x, y) and +transformed to (u, v). The Yee grid positions also have to be taken into account. The Jacobian +for the transformation of eps and E is evaluated at the r' positions of E-field components. +Similarly, the jacobian for mu and H is evaluated at the r' positions of H-field components. +Currently, the half-step offset in w is ignored, which should be a pretty good approximation.""" + +import numpy as np + + +def radial_transform(coords, radius, bend_axis): + """Compute the new coordinates and the Jacobian of a polar coordinate transformation. After + offsetting the plane such that its center is a distance of ``radius`` away from the center of + curvature, we have, e.g. for ``bend_axis=='y'``: + + u = (x**2 + z**2) + v = y + w = R acos(x / u) + + These are all evaluated at z = 0 below. + + Parameters + ---------- + coords : tuple + A tuple of two arrays of size Nx + 1, Ny + 1, respectively. + radius : float + Radius of the bend. + bend_axis : 0 or 1 + Axis normal to the bend plane. + + Returns + ------- + new_coords: tuple + Transformed coordinates, same shape as ``coords``. + jac_e: np.ndarrray + Jacobian of the transformation at the E-field positions, shape ``(3, 3, Nx * Ny)``. + jac_h: np.ndarrray + Jacobian of the transformation at the H-field positions, shape ``(3, 3, Nx * Ny)``. + k_to_kp: np.ndarray + A matrix of shape (3, 3) that transforms the k-vector from the original coordinates to the + transformed ones. + """ + + Nx, Ny = coords[0].size - 1, coords[1].size - 1 + norm_axis = 0 if bend_axis == 1 else 1 + + # Center the new coordinates such that the radius is at the center of the plane + u = coords[0] + (norm_axis == 0) * (radius - coords[0][Nx // 2]) + v = coords[1] + (norm_axis == 1) * (radius - coords[1][Ny // 2]) + new_coords = (u, v) + + """The only nontrivial derivative is dwdz and it only depends on the coordinate in the + norm_axis direction (orthogonal to both bend_axis and z). We need to compute that derivative + at the En and Hn positions. + """ + dwdz_e = radius / new_coords[norm_axis][:-1] + dwdz_h = radius / (new_coords[norm_axis][:-1] + new_coords[norm_axis][1:]) * 2 + + jac_e = np.zeros((3, 3, Nx, Ny)) + jac_e[0, 0, :, :] = 1 + jac_e[1, 1, :, :] = 1 + jac_e[2, 2, :, :] = np.expand_dims(dwdz_e, axis=bend_axis) + jac_e = jac_e.reshape((3, 3, -1)) + + jac_h = np.zeros((3, 3, Nx, Ny)) + jac_h[0, 0, :, :] = 1 + jac_h[1, 1, :, :] = 1 + jac_h[2, 2, :, :] = np.expand_dims(dwdz_h, axis=bend_axis) + jac_h = jac_h.reshape((3, 3, -1)) + + return new_coords, jac_e, jac_h + + +def angled_transform(coords, angle_theta, angle_phi): + """Compute the new coordinates and the Jacobian for a transformation that "straightens" + an angled waveguide such that it is translationally invariant in w. The transformation is + u = x - tan(angle) * z + + Parameters + ---------- + coords : tuple + A tuple of two arrays of size Nx + 1, Ny + 1, respectively. + angle_theta : float, optional + (radian) Polar angle from the normal axis. + angle_phi : float, optional + (radian) Azimuth angle in the plane orthogonal to the normal axis. + + Returns + ------- + new_coords: tuple + Transformed coordinates, same shape as ``coords``. + jac_e: np.ndarrray + Jacobian of the transformation at the E-field positions, shape ``(3, 3, Nx * Ny)``. + jac_h: np.ndarrray + Jacobian of the transformation at the H-field positions, shape ``(3, 3, Nx * Ny)``. + """ + + Nx, Ny = coords[0].size - 1, coords[1].size - 1 + + # The new coordinates are exactly the same at z = 0 + new_coords = (np.copy(c) for c in coords) + + # The only nontrivial derivatives are dudz, dvdz and they are constant everywhere + jac = np.zeros((3, 3, Nx * Ny)) + jac[0, 0, :] = 1 + jac[1, 1, :] = 1 + jac[2, 2, :] = 1 + jac[0, 2, :] = -np.tan(angle_theta) * np.cos(angle_phi) + jac[1, 2, :] = -np.tan(angle_theta) * np.sin(angle_phi) + + return new_coords, jac, jac