FXC-1617: switch tensorial mode solver to tidy3d-extras

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).

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

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

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 = "<<<MODE_SOLVER_MONITOR>>>"
 
@@ -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,