diff --git a/CHANGELOG.md b/CHANGELOG.md index 3cf63806cd..26f41183fa 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,6 +24,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - 1D lumped elements (with zero lateral extent) are no longer allowed. Use a small finite lateral extent (e.g., `1e-6`) instead. - `ModeSortSpec.sort_key` is now required with a default of `"n_eff"` (previously optional with `None` default). `ModeSortSpec.sort_order` is now optional with a default of `None`, which automatically selects the natural order based on `sort_key` and `sort_reference`: ascending when a reference is provided (closest first), otherwise descending for `n_eff` and polarization fractions (higher values first), ascending for `k_eff` and `mode_area` (lower values first). - Changed the interpretation of `waist_distance` and `waist_distances` for backward-propagating Gaussian beams (`GaussianBeam`, `AstigmaticGaussianBeam`, `GaussianBeamProfile`, `AstigmaticGaussianBeamProfile`). Previously, the waist position was interpreted relative to the directed propagation axis, meaning switching `direction` from `+` to `-` would also flip the waist position in the global reference frame. Now, the waist position is defined consistently for both directions: a positive `waist_distance` always places the beam waist behind the source/monitor plane (toward the negative normal axis), regardless of propagation direction. This ensures reciprocity between Gaussian sources and overlap monitors used in port-based S-matrix calculations. Users with existing simulations using backward-propagating Gaussian beams with non-zero waist distances may need to adjust their values. +- The fully tensorial mode solver (required for fully anisotropic media or non-zero `angle_theta` in `ModeSpec`) is now only available through `tidy3d-extras` or by running through the Tidy3D server. Users attempting to run locally without `tidy3d-extras` will receive a clear error message with instructions. ### Changed - `ModeSortSpec.sort_key` is now required with a default of `"n_eff"` (previously optional with `None` default). `ModeSortSpec.sort_order` is now optional with a default of `None`, which automatically selects the natural order based on `sort_key` and `sort_reference`: ascending when a reference is provided (closest first), otherwise descending for `n_eff` and polarization fractions (higher values first), ascending for `k_eff` and `mode_area` (lower values first). diff --git a/docs/faq b/docs/faq index 2ae3f1f737..9bc073c3fa 160000 --- a/docs/faq +++ b/docs/faq @@ -1 +1 @@ -Subproject commit 2ae3f1f7375c337a45c95ea1307320b1fad654f6 +Subproject commit 9bc073c3fa5ba8e8a56785296fe059ffc10683bd diff --git a/tests/test_components/test_mode.py b/tests/test_components/test_mode.py index 9d13bf4315..e47fea4a35 100644 --- a/tests/test_components/test_mode.py +++ b/tests/test_components/test_mode.py @@ -510,3 +510,102 @@ def test_filter_pol_with_default_sort_spec(): filter_pol="te", sort_spec=td.ModeSortSpec(sort_reference=1.5), ) + + +def _make_tensorial_mode_sim(angle_theta=0, fully_anisotropic=False): + """Helper to create a mode simulation that requires a tensorial solver.""" + structures = [] + if fully_anisotropic: + # Use a FullyAnisotropicMedium + structures.append( + td.Structure( + geometry=td.Box(size=(1, 1, td.inf)), + medium=td.FullyAnisotropicMedium( + permittivity=np.eye(3) * 4.0 + np.array([[0, 0.1, 0], [0.1, 0, 0], [0, 0, 0]]) + ), + ) + ) + else: + structures.append( + td.Structure( + geometry=td.Box(size=(1, 1, td.inf)), + medium=td.Medium(permittivity=4.0), + ) + ) + + return td.ModeSimulation( + size=(2, 2, 0), + freqs=[td.C_0], + mode_spec=td.ModeSpec(num_modes=1, angle_theta=angle_theta), + grid_spec=td.GridSpec.uniform(dl=0.2), + structures=structures, + ) + + +def test_tensorial_mode_solver_error_without_extras(monkeypatch): + """Test that attempting to run a tensorial mode solver locally without tidy3d-extras raises an error.""" + from tidy3d.components.mode import mode_solver as mode_solver_module + from tidy3d.components.mode.solver import compute_modes + from tidy3d.packaging import tidy3d_extras + + # Mock _get_solver_func to always return the base solver (simulating no tidy3d-extras) + monkeypatch.setattr(mode_solver_module, "_get_solver_func", lambda: compute_modes) + + # Also disable local subpixel to prevent tidy3d-extras from bypassing our patched code path + # when tidy3d-extras is installed and licensed + monkeypatch.setitem(tidy3d_extras, "use_local_subpixel", False) + + # Test with angle_theta (angled mode) - should raise NotImplementedError from base solver + sim = _make_tensorial_mode_sim(angle_theta=np.pi / 6) + with pytest.raises(NotImplementedError, match="tensorial mode solver"): + sim.run_local() + + # Test with fully anisotropic medium + sim_aniso = _make_tensorial_mode_sim(fully_anisotropic=True) + with pytest.raises(NotImplementedError, match="tensorial mode solver"): + sim_aniso.run_local() + + +def test_tensorial_mode_solver_with_extras(): + """Test that tensorial mode solver works when tidy3d-extras is available and usable.""" + from tidy3d.components.mode.mode_solver import _get_solver_func + from tidy3d.components.mode.solver import compute_modes + + # Check if _get_solver_func returns something other than the base solver + solver_func = _get_solver_func() + if solver_func is compute_modes: + pytest.skip("tidy3d-extras not available or not properly initialized/licensed") + + # Test with angle_theta (angled mode) + sim = _make_tensorial_mode_sim(angle_theta=np.pi / 6) + result = sim.run_local() + assert result is not None + # ModeSimulationData has n_eff through modes_raw + assert result.modes_raw.n_eff is not None + + # Test with fully anisotropic medium + sim_aniso = _make_tensorial_mode_sim(fully_anisotropic=True) + result_aniso = sim_aniso.run_local() + assert result_aniso is not None + assert result_aniso.modes_raw.n_eff is not None + + +def test_diagonal_mode_solver_still_works(): + """Test that the diagonal (non-tensorial) mode solver still works without tidy3d-extras.""" + # Simple waveguide simulation that doesn't require tensorial solver + sim = td.ModeSimulation( + size=(2, 2, 0), + freqs=[td.C_0], + mode_spec=td.ModeSpec(num_modes=1), # No angle, no anisotropy + grid_spec=td.GridSpec.uniform(dl=0.2), + structures=[ + td.Structure( + geometry=td.Box(size=(1, 1, td.inf)), + medium=td.Medium(permittivity=4.0), + ) + ], + ) + result = sim.run_local() + assert result is not None + # ModeSimulationData has n_eff through modes_raw + assert result.modes_raw.n_eff is not None diff --git a/tests/test_plugins/test_mode_solver.py b/tests/test_plugins/test_mode_solver.py index 02adf555db..c19e649f70 100644 --- a/tests/test_plugins/test_mode_solver.py +++ b/tests/test_plugins/test_mode_solver.py @@ -13,7 +13,7 @@ from tidy3d import Coords, Grid, ModeIndexDataArray, ScalarFieldDataArray, ScalarModeFieldDataArray 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 TOL_DEGENERATE_CANDIDATE, EigSolver, compute_modes +from tidy3d.components.mode.solver import TOL_DEGENERATE_CANDIDATE, EigSolver from tidy3d.components.mode_spec import MODE_DATA_KEYS from tidy3d.exceptions import DataError, SetupError, ValidationError from tidy3d.plugins.mode import ModeSolver @@ -236,21 +236,6 @@ def mock_download(resource_id, remote_filename, to_file, *args, **kwargs): ) -def test_compute_modes(): - """Test direct call to `compute_modes`.""" - eps_cross = np.random.rand(10, 10) - coords = np.arange(11) - mode_spec = td.ModeSpec(num_modes=3, target_neff=2.0) - _ = compute_modes( - eps_cross=[eps_cross] * 9, - coords=[coords, coords], - freq=td.C_0 / 1.0, - mode_spec=mode_spec, - direction="-", - precision="single", - ) - - def compare_colocation(ms): """Compare mode-solver fields with colocation applied during run or post-run.""" data_col = ms.solve() @@ -660,125 +645,6 @@ def test_mode_solver_unstructured_custom_medium(nx, cond_factor, interp, tol, tm assert error_up < tol -@td.packaging.disable_local_subpixel -def test_mode_solver_straight_vs_angled(): - """Compare results for a straight and angled nominally identical waveguides. - Note: results do not match perfectly because of the numerical grid. - """ - simulation = td.Simulation( - size=SIM_SIZE, - grid_spec=td.GridSpec.auto(wavelength=1.0, min_steps_per_wvl=16), - structures=[WAVEGUIDE], - run_time=1e-12, - symmetry=(0, 0, 1), - boundary_spec=td.BoundarySpec.all_sides(boundary=td.Periodic()), - sources=[SRC], - ) - mode_spec = td.ModeSpec(num_modes=5, group_index_step=True) - freqs = [td.C_0 / 0.9, td.C_0 / 1.0, td.C_0 / 1.1] - ms = ModeSolver( - simulation=simulation, - plane=PLANE, - mode_spec=mode_spec, - freqs=freqs, - direction="-", - ) - - angle = np.pi / 6 - width, height = WAVEGUIDE.geometry.size[0], WAVEGUIDE.geometry.size[2] - vertices = np.array( - [[-width / 2, -100, 0], [width / 2, -100, 0], [width / 2, 100, 0], [-width / 2, 100, 0]] - ) - vertices = PLANE.rotate_points(vertices.T, axis=[0, 0, 1], angle=-angle).T - vertices = [verts[:2] for verts in vertices] - wg_angled = td.Structure( - geometry=td.PolySlab(vertices=vertices, slab_bounds=(-height / 2, height / 2)), - medium=WG_MEDIUM, - ) - mode_spec_angled = mode_spec.updated_copy(angle_theta=angle) - src_angled = td.ModeSource( - source_time=td.GaussianPulse(freq0=2e14, fwidth=1e13), - center=PLANE.center, - size=PLANE.size, - mode_spec=mode_spec_angled, - direction="-", - mode_index=0, - ) - sim_angled = simulation.updated_copy(structures=[wg_angled], sources=[src_angled]) - # sim_angled.plot(z=0) - # plt.show() - - ms_angled = ModeSolver( - simulation=sim_angled, - plane=PLANE, - mode_spec=mode_spec_angled, - freqs=freqs, - direction="-", - ) - - check_ms_reduction(ms) - check_ms_reduction(ms_angled) - - for key, val in ms.data.modes_info.items(): - tol = 1e-2 - if key == "TE (Ex) fraction": - tol = 0.1 - elif key == "wg TE fraction": - tol = 1.3e-2 - elif key == "mode area": - tol = 2.1e-2 - elif key == "dispersion (ps/(nm km))": - tol = 0.7 - # print( - # key, - # (np.abs(val - ms_angled.data.modes_info[key]) / np.abs(val)).values.max(), - # (np.abs(val - ms_angled.data.modes_info[key]) / np.abs(ms_angled.data.modes_info[key])).values.max(), - # ) - assert np.allclose(val, ms_angled.data.modes_info[key], rtol=tol) - - -def test_mode_solver_angle_bend(): - """Run mode solver with angle and bend and symmetry""" - simulation = td.Simulation( - size=SIM_SIZE, - grid_spec=td.GridSpec(wavelength=1.0), - structures=[WAVEGUIDE], - run_time=1e-12, - symmetry=(-1, 0, 1), - boundary_spec=td.BoundarySpec.all_sides(boundary=td.Periodic()), - sources=[SRC], - ) - mode_spec = td.ModeSpec( - num_modes=3, - target_neff=2.0, - bend_radius=3, - bend_axis=0, - angle_theta=np.pi / 3, - angle_phi=np.pi, - sort_spec=td.ModeSortSpec(track_freq="highest"), - ) - # put plane entirely in the symmetry quadrant rather than sitting on its center - plane = td.Box(center=(0, 0.5, 0), size=(1, 0, 1)) - ms = ModeSolver( - simulation=simulation, plane=plane, mode_spec=mode_spec, freqs=[td.C_0 / 1.0], direction="-" - ) - compare_colocation(ms) - verify_pol_fraction(ms) - verify_dtype(ms) - _ = ms.data.to_dataframe() - check_ms_reduction(ms) - - # Plot field - _, ax = plt.subplots(1) - ms.plot_field("Ex", ax=ax, mode_index=1) - plt.close() - - # Create source and monitor - st = td.GaussianPulse(freq0=1.0e12, fwidth=1.0e12) - _ = ms.to_source(source_time=st, direction="-") - _ = ms.to_monitor(freqs=np.array([1.0, 2.0]) * 1e12, name="mode_mnt") - - def test_mode_bend_radius(): """Test that the bend radius is correctly applied to the center of the mode plane in the case of an auto-grid that is not symmetric w.r.t. that center, and that nominally identical diff --git a/tidy3d/components/mode/mode_solver.py b/tidy3d/components/mode/mode_solver.py index 5acc5f4e64..f8bdaa4861 100644 --- a/tidy3d/components/mode/mode_solver.py +++ b/tidy3d/components/mode/mode_solver.py @@ -85,7 +85,12 @@ if TYPE_CHECKING: from matplotlib.colors import Colormap -from tidy3d.packaging import supports_local_subpixel, tidy3d_extras +from tidy3d.packaging import ( + Tidy3dImportError, + _check_tidy3d_extras_available, + supports_local_subpixel, + tidy3d_extras, +) # Importing the local solver may not work if e.g. scipy is not installed IMPORT_ERROR_MSG = """Could not import local solver, 'ModeSolver' objects can still be constructed @@ -99,6 +104,42 @@ log.warning(IMPORT_ERROR_MSG) LOCAL_SOLVER_IMPORTED = False + +def _get_solver_func(): + """Get the best available mode solver function. + + Returns the tidy3d-extras compute_modes if available (handles all cases including + fully tensorial), otherwise falls back to the base compute_modes which will raise + an informative error if a tensorial solve is attempted. + + Returns + ------- + callable + The compute_modes function to use for solving. + + Raises + ------ + ImportError + If the local solver could not be imported (e.g., scipy not installed). + """ + if not LOCAL_SOLVER_IMPORTED: + raise ImportError(IMPORT_ERROR_MSG) + + # Try to get tidy3d-extras solver (handles all cases including tensorial) + try: + _check_tidy3d_extras_available(quiet=True) + if tidy3d_extras["mod"] is not None: + from tidy3d_extras.mode import EigSolver as ExtrasEigSolver + + return ExtrasEigSolver.compute_modes + except (Tidy3dImportError, ImportError, AttributeError): + # tidy3d-extras is optional; if unavailable or incompatible, fall back to base solver. + pass + + # Fall back to base solver (will raise error if tensorial solve is attempted) + return compute_modes + + FIELD = tuple[ArrayComplex3D, ArrayComplex3D, ArrayComplex3D] MODE_MONITOR_NAME = "<<>>" @@ -1644,10 +1685,9 @@ def _solve_single_freq( The fields are rotated from propagation coordinates back to global coordinates. """ - if not LOCAL_SOLVER_IMPORTED: - raise ImportError(IMPORT_ERROR_MSG) + solver_func = _get_solver_func() - solver_fields, n_complex, eps_spec = compute_modes( + solver_fields, n_complex, eps_spec = solver_func( eps_cross=self._solver_eps(freq), coords=coords, freq=freq, @@ -1701,14 +1741,13 @@ def _solve_single_freq_relative( Modes are computed as linear combinations of ``basis_fields``. """ - if not LOCAL_SOLVER_IMPORTED: - raise ImportError(IMPORT_ERROR_MSG) + solver_func = _get_solver_func() solver_basis_fields = self._postprocess_solver_fields_inverse( fields=basis_fields, normal_axis=self.normal_axis, plane=self.plane ) - solver_fields, n_complex, eps_spec = compute_modes( + solver_fields, n_complex, eps_spec = solver_func( eps_cross=self._solver_eps(freq), coords=coords, freq=freq, diff --git a/tidy3d/components/mode/solver.py b/tidy3d/components/mode/solver.py index e533c3d084..f8b419e460 100644 --- a/tidy3d/components/mode/solver.py +++ b/tidy3d/components/mode/solver.py @@ -38,6 +38,15 @@ # Consider a material to be good conductor if |ep| (or |mu|) > GOOD_CONDUCTOR_THRESHOLD * |pec_val| GOOD_CONDUCTOR_THRESHOLD = 0.9 +# Error message for tensorial mode solver without tidy3d-extras +TENSORIAL_SOLVER_ERROR_MSG = ( + "The fully tensorial mode solver is required for fully anisotropic media " + "or non-zero 'angle_theta' in the 'ModeSpec', but it is not available in the base " + "'tidy3d' package. To run locally, please install 'tidy3d-extras' using " + r"'pip install tidy3d\[extras]'. Alternatively, you can run the mode solver " + "through the Tidy3D server using 'web.run(...)'." +) + class EigSolver(Tidy3dBaseModel): """Interface for computing eigenvalues given permittivity and mode spec. @@ -757,132 +766,13 @@ def solver_tensorial( Nxy=None, dmin_pmc=None, ): - """EM eigenmode solver assuming ``eps`` or ``mu`` have off-diagonal elements.""" - import scipy.sparse as sp - - mode_solver_type = "tensorial" - N = eps.shape[-1] - dxf, dxb, dyf, dyb = der_mats - - # 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 - ) - - mat = sp.bmat( - [ - [axax, axay, axbx, axby], - [ayax, ayay, aybx, ayby], - [bxax, bxay, bxbx, bxby], - [byax, byay, bybx, byby], - ] - ) - - # The eigenvalues for the matrix above are 1j * (neff + 1j * keff) - # Multiply the matrix by -1j, so that eigenvalues are (neff + 1j * keff) - mat *= -1j - - # change matrix sign for backward direction - if direction == "-": - mat *= -1 - - # Cast matrix to target data type - mat_dtype = cls.matrix_data_type(eps, mu, der_mats, mat_precision, is_tensorial=True) - mat = cls.type_conversion(mat, mat_dtype) - - # Trim small values in single precision case - if mat_precision == "single": - cls.trim_small_values(mat, tol=fp_eps) + """EM eigenmode solver assuming ``eps`` or ``mu`` have off-diagonal elements. - # Casting starting vector to target data type - vec_init = cls.type_conversion(vec_init, mat_dtype) - - # Starting eigenvalue guess in target data type - eig_guess = cls.type_conversion(np.array([neff_guess]), mat_dtype)[0] - - # Call the eigensolver. - vals, vecs = cls.solver_eigs( - mat, - num_modes, - vec_init, - guess_value=eig_guess, - mode_solver_type=mode_solver_type, - ) - neff, keff = cls.eigs_to_effective_index(vals, mode_solver_type) - # 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) - - # 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 E, H, neff, keff + Note: The fully tensorial mode solver is only available through ``tidy3d-extras`` + or by running the mode solver through the Tidy3D server. This method raises an error + directing users to those options. + """ + raise NotImplementedError(TENSORIAL_SOLVER_ERROR_MSG) @classmethod def solver_eigs(