diff --git a/CHANGELOG.md b/CHANGELOG.md index f51376c4f0..24cb5ee21c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,6 +22,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Fixed bug where an extra spatial coordinate could appear in `complex_flux` and `ImpedanceCalculator` results. - Fixed normal for `Box` shape gradient computation to always point outward from boundary which is needed for correct PEC handling. - Fixed `Box` gradients within `GeometryGroup` where the group intersection boundaries were forwarded. +- Improved degenerate mode handling in the mode solver to ensure modes respect the bi-orthogonality condition. ## [2.10.0rc3] - 2025-11-26 diff --git a/tests/test_plugins/test_mode_solver.py b/tests/test_plugins/test_mode_solver.py index 97155fdd56..2d0674c174 100644 --- a/tests/test_plugins/test_mode_solver.py +++ b/tests/test_plugins/test_mode_solver.py @@ -13,7 +13,7 @@ from tidy3d import ScalarFieldDataArray from tidy3d.components.data.monitor_data import ModeSolverData from tidy3d.components.mode.derivatives import create_sfactor_b, create_sfactor_f -from tidy3d.components.mode.solver import compute_modes +from tidy3d.components.mode.solver import TOL_DEGENERATE_CANDIDATE, EigSolver, compute_modes from tidy3d.components.mode_spec import MODE_DATA_KEYS from tidy3d.exceptions import DataError, SetupError from tidy3d.plugins.mode import ModeSolver @@ -1503,3 +1503,52 @@ def test_sort_spec_track_freq(): assert np.allclose(modes_lowest.Ex.abs, modes_lowest_retracked.Ex.abs) assert np.all(modes_lowest.n_eff == modes_lowest_retracked.n_eff) assert np.all(modes_lowest.n_group == modes_lowest_retracked.n_group) + + +def test_degenerate_mode_processing(): + """Ensure degenerate modes returned by mode solver are bi-orthogonal.""" + freq0 = td.C_0 + sim_size = (0, 2, 2) + inf = 10 + W1 = 0.3 + n = 1.5 + num_modes = 4 + mode_spec = td.ModeSpec(num_modes=num_modes) + medium = td.Medium(permittivity=n**2) + geom1 = td.Box.from_bounds((-inf, -W1 / 2, -W1 / 2), (inf, W1 / 2, W1 / 2)) + wg1 = td.Structure(geometry=geom1, medium=medium) + + grid_spec = td.GridSpec.uniform(dl=0.2) + + sim = td.Simulation( + size=sim_size, + structures=[wg1], + grid_spec=grid_spec, + run_time=10 / freq0, + ) + + ms = ModeSolver( + simulation=sim, + plane=sim.geometry, + mode_spec=mode_spec, + freqs=[freq0], + direction="+", + ) + + mode_data = ms.data_raw + + degen_sets = EigSolver._identify_degenerate_modes( + mode_data.n_complex.values[0, :], TOL_DEGENERATE_CANDIDATE + ) + assert len(degen_sets) == 1 + + S = mode_data.outer_dot(mode_data, conjugate=False).isel(f=0).values + threshold = 1e-7 + off_diag_mask = ~np.eye(S.shape[0], dtype=bool) + large_vals = np.abs(S) > threshold + problem_mask = off_diag_mask & large_vals + + indices = np.argwhere(problem_mask) + msg = f"Found {len(indices)} off-diagonal values > {threshold}:\n" + msg += "\n".join(f" |S[{i},{j}]| = {np.abs(S[i, j]):.4e}" for i, j in indices) + assert not np.any(problem_mask), msg diff --git a/tidy3d/components/mode/solver.py b/tidy3d/components/mode/solver.py index bb107d9d16..fccc205781 100644 --- a/tidy3d/components/mode/solver.py +++ b/tidy3d/components/mode/solver.py @@ -18,6 +18,10 @@ TOL_COMPLEX = 1e-10 # Tolerance for eigs TOL_EIGS = fp_eps / 10 +# Relative tolerance of complex effective index for finding candidate degenerate modes +TOL_DEGENERATE_CANDIDATE = 1e3 * TOL_EIGS +# Tolerance used to determine whether possible degenerate modes need manual orthogonalization +TOL_NEEDS_ORTHOGONALIZATION = 1e2 * TOL_EIGS # Tolerance for deciding on the matrix to be diagonal or tensorial TOL_TENSORIAL = 1e-6 # shift target neff by this value, both rel and abs, whichever results in larger shift @@ -285,6 +289,25 @@ def compute_modes( 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)) + + # Normalize all modes so self-dot equals 1.0 + all_modes = list(range(num_modes)) + E, H = cls._normalize_modes(E, H, dl_f, dl_b, all_modes) + + # Identify and orthogonalize candidate degenerate modes with a coarser tolerance + degenerate_modes = cls._identify_degenerate_modes( + neff + keff * 1j, TOL_DEGENERATE_CANDIDATE + ) + # Orthogonalize candidate degenerate modes that result in an overlap greater than + # the required tolerance + filtered_degenerate_modes = cls._identify_degenerate_modes_with_dot( + degenerate_modes, E, H, dl_f, dl_b, TOL_NEEDS_ORTHOGONALIZATION + ) + + E, H = cls._make_orthogonal_basis_for_degenerate_modes( + filtered_degenerate_modes, E, H, dl_f, dl_b + ) + fields = np.stack((E, H), axis=0) neff = neff * np.linalg.norm(kp_to_k) @@ -1089,6 +1112,248 @@ def mode_plane_contain_good_conductor(material_response) -> bool: return False return np.any(np.abs(material_response) > GOOD_CONDUCTOR_THRESHOLD * np.abs(pec_val)) + @staticmethod + def _identify_degenerate_modes( + n_complex: np.ndarray, + tol: float, + ) -> list[set[int]]: + """Inspects the n_complex of modes to find sets of degenerate modes.""" + num_modes = len(n_complex) + ungrouped = set(range(num_modes)) + degenerate_sets = [] + + while ungrouped: + # Start a new group with an ungrouped mode + seed = ungrouped.pop() + current_set = {seed} + + # Find all ungrouped modes similar to the seed + to_remove = set() + for mode_idx in ungrouped: + if np.isclose(n_complex[mode_idx], n_complex[seed], rtol=tol, atol=tol): + current_set.add(mode_idx) + to_remove.add(mode_idx) + + # Remove grouped modes from ungrouped set + ungrouped -= to_remove + + # Only keep groups with more than one mode + if len(current_set) >= 2: + degenerate_sets.append(current_set) + + return degenerate_sets + + @staticmethod + def _identify_degenerate_modes_with_dot( + degenerate_sets: list[set[int]], + E: np.ndarray, + H: np.ndarray, + dl_primal: np.ndarray, + dl_dual: np.ndarray, + tol: float, + ) -> list[set[int]]: + """Inspects the biorthogonality condition between modes to filter sets of degenerate modes. + + For each set of modes identified as degenerate by neff similarity, this function computes + the overlap matrix and uses connected components to find which modes are truly degenerate + (have significant mutual overlap). + + Returns only groups with >= 2 modes that have significant overlap. + """ + import scipy.sparse as sp + + filtered_degenerate_sets = [] + + for degenerate_group in degenerate_sets: + # Convert set to list for indexing + mode_indices = sorted(degenerate_group) + + # S is a symmetric overlap matrix of the degenerate modes + S = EigSolver._outer_dot(E, H, dl_primal, dl_dual, mode_indices) + + # Create adjacency matrix + # True if overlap is significant, False if orthogonal + # Zero out diagonal so a mode doesn't connect to itself + adjacency = (np.abs(S) > tol).astype(int) + np.fill_diagonal(adjacency, 0) + + # Find connected components + # labels: array indicating which group each mode belongs to + _, labels = sp.csgraph.connected_components(sp.csr_matrix(adjacency), directed=False) + + # Group modes by their connected component label + groups = {} + for label_idx, group_id in enumerate(labels): + mode_index = mode_indices[label_idx] + if group_id not in groups: + groups[group_id] = [] + groups[group_id].append(mode_index) + + # Convert groups dict to list of sets, keeping only groups with >= 2 modes + for group_modes in groups.values(): + if len(group_modes) >= 2: + filtered_degenerate_sets.append(set(group_modes)) + + return filtered_degenerate_sets + + @staticmethod + def _dot( + E: np.ndarray, + H: np.ndarray, + dl_primal: np.ndarray, + dl_dual: np.ndarray, + mode_1: int, + mode_2: int, + ) -> complex: + """Dot product based on the bi-orthogonality relationship between E and H.""" + # Extract field components + Ex = E[0, ...] + Ey = E[1, ...] + Hx = H[0, ...] + Hy = H[1, ...] + + # Make the differential area elements + Ex_Hy_dS = np.outer(dl_primal[0], dl_dual[1]) + Ey_Hx_dS = np.outer(dl_dual[0], dl_primal[1]) + + term1 = Ex[..., mode_1] * Hy[..., mode_2] + Ex[..., mode_2] * Hy[..., mode_1] + term1 *= Ex_Hy_dS[..., np.newaxis] + term2 = Ey[..., mode_1] * Hx[..., mode_2] + Ey[..., mode_2] * Hx[..., mode_1] + term2 *= Ey_Hx_dS[..., np.newaxis] + return np.sum(term1 - term2) + + @staticmethod + def _outer_dot( + E: np.ndarray, + H: np.ndarray, + dl_primal: np.ndarray, + dl_dual: np.ndarray, + mode_indices: set[int], + ) -> np.ndarray: + """Vectorized modal overlap matrix calculation for a set of modes. + + This overlap is based on the bi-orthogonality relationship between E and H. + Returns a matrix S where S[i, j] is the overlap between mode i and mode j. + + Note: The overlap matrix is symmetric, so we only compute the upper triangle + (including diagonal) and fill the lower triangle by symmetry. + """ + # Convert set to sorted list for consistent indexing + mode_list = sorted(mode_indices) + n_modes = len(mode_list) + + # Extract field components + Ex_sel = E[0][..., mode_list] # (Nx, Ny, 1, n) + Ey_sel = E[1][..., mode_list] + Hx_sel = H[0][..., mode_list] + Hy_sel = H[1][..., mode_list] + + # Make the differential area elements + Ex_Hy_dS = np.outer(dl_primal[0], dl_dual[1]) + Ey_Hx_dS = np.outer(dl_dual[0], dl_primal[1]) + + # Initialize output matrix + dtype = Ex_sel.dtype + S = np.zeros((n_modes, n_modes), dtype=dtype) + + # Only compute upper triangle (including diagonal) + for i in range(n_modes): + # Vectorize over j >= i + j_indices = np.arange(i, n_modes) + + # Extract mode i fields + Ex_i = Ex_sel[..., i : i + 1, np.newaxis] # (Nx, Ny, 1, 1, 1) + Ey_i = Ey_sel[..., i : i + 1, np.newaxis] + Hy_i = Hy_sel[..., i : i + 1, np.newaxis] + Hx_i = Hx_sel[..., i : i + 1, np.newaxis] + + # Extract mode j fields (vectorized for j >= i) + Ex_j = Ex_sel[..., np.newaxis, j_indices] # (Nx, Ny, 1, 1, n_j) + Ey_j = Ey_sel[..., np.newaxis, j_indices] + Hy_j = Hy_sel[..., np.newaxis, j_indices] + Hx_j = Hx_sel[..., np.newaxis, j_indices] + + # Compute term1: (Ex[i] * Hy[j] + Ex[j] * Hy[i]) * dS + term1 = (Ex_i * Hy_j + Ex_j * Hy_i) * Ex_Hy_dS[..., np.newaxis, np.newaxis, np.newaxis] + + # Compute term2: (Ey[i] * Hx[j] + Ey[j] * Hx[i]) * dS + term2 = (Ey_i * Hx_j + Ey_j * Hx_i) * Ey_Hx_dS[..., np.newaxis, np.newaxis, np.newaxis] + + # Sum over spatial dimensions to get S[i, j] for j >= i + S[i, j_indices] = np.sum(term1 - term2, axis=(0, 1))[0, 0, :] # (1, 1, n_j) -> (n_j,) + + # Fill lower triangle by symmetry + S = S + S.T - np.diag(np.diag(S)) + + return S + + @staticmethod + def _normalize_modes( + E: np.ndarray, + H: np.ndarray, + dl_primal: np.ndarray, + dl_dual: np.ndarray, + mode_indices: list[int], + ) -> tuple[np.ndarray, np.ndarray]: + """Normalize modes so that their self-dot equals 1.0.""" + for mode_idx in mode_indices: + self_dot = EigSolver._dot(E, H, dl_primal, dl_dual, mode_idx, mode_idx) + # Avoid division by zero + dtype = E.dtype + if np.isclose(self_dot, 0, atol=np.finfo(dtype).tiny): + raise ValueError(f"Eigenvector {mode_idx} is self-orthogonal. Cannot normalize.") + norm_factor = 1.0 / np.sqrt(self_dot) + E[..., mode_idx] *= norm_factor + H[..., mode_idx] *= norm_factor + return E, H + + @staticmethod + def _make_orthogonal_basis_for_degenerate_modes( + degenerate_mode_sets: list[set[int]], + E_vec: np.ndarray, + H_vec: np.ndarray, + dl_primal: np.ndarray, + dl_dual: np.ndarray, + ) -> tuple[np.ndarray, np.ndarray]: + """Ensures that groups of degenerate modes are orthogonal, which is not guaranteed for eigenvectors with degenerate eigenvalues.""" + dtype = E_vec.dtype + for degenerate_group in degenerate_mode_sets: + # Convert set to list for indexing + mode_indices = sorted(degenerate_group) + + # S is a symmetric overlap matrix of the degenerate modes. + S = EigSolver._outer_dot(E_vec, H_vec, dl_primal, dl_dual, mode_indices) + + # 1. Diagonalize + # eigenvalues (w) and right eigenvectors (v) + eigvals, eigvecs = np.linalg.eig(S) + # 2. Normalize Eigenvectors using UNCONJUGATED dot product + # Standard np.linalg.eig normalizes so v.conj().T @ v = 1 (Hermitian) + # We need v.T @ v = 1 (Complex Symmetric) + for i in range(eigvecs.shape[1]): + vec = eigvecs[:, i] + # Calculate complex self-dot (no conjugation) + self_dot = np.dot(vec, vec) + + # Avoid division by zero + if np.isclose(self_dot, 0, atol=np.finfo(dtype).tiny): + raise ValueError(f"Eigenvector {i} is self-orthogonal. Cannot orthogonalize.") + + scale_factor = np.sqrt(self_dot) + eigvecs[:, i] = vec / scale_factor + # 3. Calculate Inverse Square Root of Eigenvalues + lambda_inv_sqrt = np.diag(1.0 / np.sqrt(eigvals)) + # 4. Construct W + # Since V is complex orthogonal (V.T @ V = I), V_inv = V.T + # W = V * Lambda^(-1/2) * V.T + W = eigvecs @ lambda_inv_sqrt @ eigvecs.T + # S_new will be identity + # S_new = W.T @ S @ W + E_vec[..., mode_indices] = E_vec[..., mode_indices] @ W + H_vec[..., mode_indices] = H_vec[..., mode_indices] @ W + + return E_vec, H_vec + def compute_modes(*args: Any, **kwargs: Any) -> tuple[Numpy, Numpy, str]: """A wrapper around ``EigSolver.compute_modes``, which is used in :class:`.ModeSolver`."""