From f8f96835f941a1d3c8e9cfd2c6d6c77cfce185ea Mon Sep 17 00:00:00 2001 From: dmarek Date: Wed, 3 Dec 2025 16:26:42 -0500 Subject: [PATCH] fix(mode): ensure degenerate modes are orthogonal from mode solver --- CHANGELOG.md | 1 + tests/test_plugins/test_mode_solver.py | 62 +++- tidy3d/components/mode/solver.py | 377 +++++++++++++++++++++++++ 3 files changed, 438 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8f8749ce26..2f8a7fe256 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -26,6 +26,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - 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. - Fixed `Box` gradients to use automatic permittivity detection for inside/outside permittivity. +- 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..cdc412513d 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 @@ -963,7 +963,16 @@ def test_mode_solver_nan_pol_fraction(): md = ms.solve() check_ms_reduction(ms) - + # Inject NaN at mode_index=5 for selected field components + nan_fields = {} + for field_name in ["Ex", "Ez", "Hx", "Hz"]: + field = getattr(md, field_name) + data = field.values.copy() + data[..., 5] = np.nan + nan_fields[field_name] = field.copy(data=data) + + md = md.updated_copy(**nan_fields) + md = ms._filter_polarization(md) assert list(np.where(np.isnan(md.pol_fraction.te))[1]) == [9] @@ -1503,3 +1512,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..e533c3d084 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,7 @@ 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)) + fields = np.stack((E, H), axis=0) neff = neff * np.linalg.norm(kp_to_k) @@ -382,6 +387,13 @@ def conductivity_model_for_good_conductor( mu_offd = np.abs(mu_tensor[off_diagonals]) is_tensorial = np.any(eps_offd > TOL_TENSORIAL) or np.any(mu_offd > TOL_TENSORIAL) + # Determine if ``eps`` and ``mu`` represent reciprocal media + is_reciprocal = True + if is_tensorial: + is_reciprocal = cls._check_reciprocity( + eps_tensor, TOL_TENSORIAL + ) and cls._check_reciprocity(mu_tensor, TOL_TENSORIAL) + # initial vector for eigensolver in correct data type vec_init = cls.set_initial_vec(Nx, Ny, is_tensorial=is_tensorial) @@ -440,6 +452,24 @@ def conductivity_model_for_good_conductor( dmin_pmc=dmin_pmc, ) + # We only can orthogonalize reciprocal modes in this manner. + if is_reciprocal: + dl_f, dl_b = dls + # Normalize all modes so self-dot equals 1.0 and find any self-orthogonal modes + all_modes = list(range(num_modes)) + original_shape = E.shape + E = E.reshape((3, Nx, Ny, num_modes)) + H = H.reshape((3, Nx, Ny, num_modes)) + E, H, self_orthogonal_modes = cls._normalize_modes( + E, H, dl_f, dl_b, all_modes, TOL_NEEDS_ORTHOGONALIZATION + ) + # Identify and orthogonalize degenerate modes, excluding self-orthogonal modes + E, H = cls._orthogonalize_degenerate_modes( + E, H, neff, keff, dl_f, dl_b, self_orthogonal_modes + ) + E = E.reshape(original_shape) + H = H.reshape(original_shape) + return E, H, neff, keff, eps_spec @classmethod @@ -705,6 +735,13 @@ def trim_small_values(cls, mat: sp.csr_matrix, tol: float) -> sp.csr_matrix: mat.eliminate_zeros() return mat + @staticmethod + def _check_reciprocity(material_tensor: np.ndarray, tol: float) -> bool: + """Check if material tensor is symmetric (reciprocal): tensor[i,j] ≈ tensor[j,i].""" + diff = np.abs(material_tensor - material_tensor.transpose(1, 0, 2)) + max_error = np.max(diff) + return max_error <= tol + @classmethod def solver_tensorial( cls, @@ -1089,6 +1126,346 @@ 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 + term2 = Ey[..., mode_1] * Hx[..., mode_2] + Ey[..., mode_2] * Hx[..., mode_1] + term2 *= Ey_Hx_dS + return (1 / 4) * np.sum(term1 - term2) + + @staticmethod + def _cauchy_schwarz_dot_bound( + E: np.ndarray, + H: np.ndarray, + dl_primal: np.ndarray, + dl_dual: np.ndarray, + mode_idx: int, + ) -> complex: + """Calculates the upper bound for the self-dot (overlap) using the Cauchy-Schwarz inequality. + + Hard to predict how the normalization of the eigenvectors from the linear algebra solver + translates to normalization in terms of the EM overlap integral, + so we use this upper bound to determine whether a mode should be considered as self-orthogonal. + + Sketch of derivation for the upper bound: + | S (E x H) . z dA | (Absolute value of the bi-orthogonality overlap integral) + <= S |(E x H) . z| dA (Triangle Inequality) + <= S |E| * |H| dA (Cross Product Magnitude) + <= sqrt(S |E|^2) * sqrt(S |H|^2) (Cauchy-Schwarz) + """ + # Extract all field components + Ex = E[0, ...] + Ey = E[1, ...] + Ez = E[2, ...] + Hx = H[0, ...] + Hy = H[1, ...] + Hz = H[2, ...] + + # 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]) + Ez_dS = np.outer(dl_dual[0], dl_dual[1]) + Hz_dS = np.outer(dl_primal[0], dl_primal[1]) + + E_int = np.sum(Ex[..., mode_idx] * Ex[..., mode_idx].conj() * Ex_Hy_dS) + E_int += np.sum(Ey[..., mode_idx] * Ey[..., mode_idx].conj() * Ey_Hx_dS) + E_int += np.sum(Ez[..., mode_idx] * Ez[..., mode_idx].conj() * Ez_dS) + E_int = E_int / 2 + + H_int = np.sum(Hx[..., mode_idx] * Hx[..., mode_idx].conj() * Ey_Hx_dS) + H_int += np.sum(Hy[..., mode_idx] * Hy[..., mode_idx].conj() * Ex_Hy_dS) + H_int += np.sum(Hz[..., mode_idx] * Hz[..., mode_idx].conj() * Hz_dS) + H_int = H_int / 2 + + return np.sqrt(E_int * H_int) + + @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] # (Nx, Ny, 1, 1) + Ey_i = Ey_sel[..., i : i + 1] + Hy_i = Hy_sel[..., i : i + 1] + Hx_i = Hx_sel[..., i : i + 1] + + # Extract mode j fields (vectorized for j >= i) + Ex_j = Ex_sel[..., j_indices] # (Nx, Ny, 1, n_j) + Ey_j = Ey_sel[..., j_indices] + Hy_j = Hy_sel[..., j_indices] + Hx_j = Hx_sel[..., 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] + + # 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] + + # Sum over spatial dimensions to get S[i, j] for j >= i + S[i, j_indices] = (1 / 4) * np.sum(term1 - term2, axis=(0, 1)) + + # 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], + self_orthogonal_tol: float, + ) -> tuple[np.ndarray, np.ndarray, set[int]]: + """Normalize modes so that their self-dot equals 1.0. + + Returns: + Tuple of (E, H, self_orthogonal_modes) where self_orthogonal_modes is a set of mode indices + that are self-orthogonal and cannot be normalized. + """ + self_orthogonal_modes = set() + for mode_idx in mode_indices: + self_dot = EigSolver._dot(E, H, dl_primal, dl_dual, mode_idx, mode_idx) + self_dot_upper_bound = EigSolver._cauchy_schwarz_dot_bound( + E, H, dl_primal, dl_dual, mode_idx + ) + # Check if mode is self-orthogonal or otherwise not suitable for orthogonalization procedure + if ( + np.abs(self_dot) < self_dot_upper_bound * self_orthogonal_tol + or np.any(np.isnan(E[..., mode_idx])) + or np.any(np.isnan(H[..., mode_idx])) + ): + self_orthogonal_modes.add(mode_idx) + continue + norm_factor = 1.0 / np.sqrt(self_dot) + E[..., mode_idx] *= norm_factor + H[..., mode_idx] *= norm_factor + + return E, H, self_orthogonal_modes + + @classmethod + def _orthogonalize_degenerate_modes( + cls, + E: np.ndarray, + H: np.ndarray, + neff: np.ndarray, + keff: np.ndarray, + dl_primal: np.ndarray, + dl_dual: np.ndarray, + self_orthogonal_modes: set[int], + ) -> tuple[np.ndarray, np.ndarray]: + """Identify and orthogonalize degenerate modes, excluding self-orthogonal modes.""" + # Identify and orthogonalize candidate degenerate modes with a coarser tolerance + degenerate_modes = cls._identify_degenerate_modes( + neff + keff * 1j, TOL_DEGENERATE_CANDIDATE + ) + + # Remove self-orthogonal modes from degenerate mode sets since they cannot be orthogonalized + degenerate_modes = [ + mode_set - self_orthogonal_modes + for mode_set in degenerate_modes + if len(mode_set - self_orthogonal_modes) + > 1 # Only keep sets with at least 2 valid modes + ] + + # 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_primal, dl_dual, TOL_NEEDS_ORTHOGONALIZATION + ) + + E, H = cls._make_orthogonal_basis_for_degenerate_modes( + filtered_degenerate_modes, E, H, dl_primal, dl_dual + ) + + 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`."""