diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml index f1c2e85..376130e 100644 --- a/.github/workflows/publish.yml +++ b/.github/workflows/publish.yml @@ -25,6 +25,8 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@v4 + with: + persist-credentials: false - name: Install uv uses: astral-sh/setup-uv@v5 @@ -52,6 +54,8 @@ jobs: steps: - uses: actions/checkout@v4 + with: + persist-credentials: false - name: Install uv uses: astral-sh/setup-uv@v5 @@ -92,6 +96,8 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@v4 + with: + persist-credentials: false - name: Download artifacts uses: actions/download-artifact@v4 diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 1e09c59..9f4fd9e 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -10,8 +10,47 @@ permissions: contents: read jobs: + lint: + name: Lint and type check + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v4 + with: + persist-credentials: false + + - name: Install Rust tooling + run: rustup component add rustfmt clippy + + - name: Install uv + uses: astral-sh/setup-uv@v5 + with: + version: "latest" + + - name: Set up Python 3.13 + run: uv python install 3.13 + + - name: Install dependencies + run: uv sync --all-extras + + - name: Check Rust formatting + run: cargo fmt --all -- --check + + - name: Run Rust clippy + run: cargo clippy --all-targets --all-features -- -D warnings + + - name: Check Python formatting + run: uv run ruff format --check . + + - name: Run Python lint + run: uv run ruff check . + + - name: Run Python type check + run: uv run pyright + test-portable: name: Portable tests (Python ${{ matrix.python-version }}) + needs: lint runs-on: ubuntu-latest strategy: fail-fast: false @@ -20,6 +59,8 @@ jobs: steps: - uses: actions/checkout@v4 + with: + persist-credentials: false - name: Install uv uses: astral-sh/setup-uv@v5 @@ -83,6 +124,8 @@ jobs: steps: - uses: actions/checkout@v4 + with: + persist-credentials: false - name: Download local coverage badge uses: actions/download-artifact@v4 @@ -105,6 +148,7 @@ jobs: test-full: name: Full portable tests + needs: lint runs-on: ubuntu-latest steps: @@ -132,6 +176,7 @@ jobs: test-windows: name: Windows portable tests + needs: lint runs-on: windows-latest steps: diff --git a/.gitignore b/.gitignore index 99287e6..e4f0cfb 100644 --- a/.gitignore +++ b/.gitignore @@ -12,3 +12,4 @@ examples/example_outputs_contact_sheet.png examples/material_grid_outputs/ examples/ridge_waveguide_outputs/ examples/soi_hybridization_outputs/ +examples/tidy3d_modal_outputs/ \ No newline at end of file diff --git a/CHANGELOG.md b/CHANGELOG.md index 1f71e15..6aff1c0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,12 @@ # Changelog -## 0.1.0a3 - Unreleased +## 0.1.0a4 - Unreleased + +- Fixed y-normal field mapping so returned global fields use a right-handed + basis and physical +y power overlap is positive. +- Added a regression test for y-normal power-overlap sign. + +## 0.1.0a3 Initial alpha release candidate. diff --git a/Cargo.lock b/Cargo.lock index e87cdf3..4800e51 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -86,7 +86,7 @@ dependencies = [ [[package]] name = "micromode-core" -version = "0.1.0-alpha.3" +version = "0.1.0-alpha.4" dependencies = [ "amd", "nalgebra", diff --git a/Cargo.toml b/Cargo.toml index 5632be9..00edb90 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "micromode-core" -version = "0.1.0-alpha.3" +version = "0.1.0-alpha.4" edition = "2021" description = "Rust core for the MicroMode photonics mode solver." license = "Apache-2.0" @@ -32,3 +32,7 @@ rlu = "0.7" default = [] python = ["pyo3"] extension-module = ["python", "pyo3/extension-module"] + +[lints.clippy] +needless_range_loop = "allow" +too_many_arguments = "allow" diff --git a/README.md b/README.md index d2338db..e1c0297 100644 --- a/README.md +++ b/README.md @@ -57,6 +57,33 @@ data.plot_field("Ex", mode_index=0) data.to_hdf5("modes.h5") ``` +## Examples + + +### Tidy3D Waveguide +![Tidy3D modal monitor example](docs/assets/tidy3d_modal_modes.png) + +The Tidy3D modal monitor example recreates the strip-waveguide setup from +Flexcompute's modal sources and monitors notebook. It solves the first three +x-propagating modes of a silicon waveguide on a silica substrate and plots +`|Ey|` and `|Ez|` on the same y-z mode plane. (See [Tidy3D, "Defining Mode Sources and Monitors"](https://www.flexcompute.com/tidy3d/examples/notebooks/ModalSourcesMonitors/).) + +```bash +uv run --extra dev python examples/tidy3d_modal_sources_monitors.py +``` + +### Hybridization Sweep +![Hybridization sweep example](docs/assets/hybridization_sweep.png) + +The SOI hybridization example sweeps the width of a 220 nm silicon ridge and +solves several modes at each step. It shows how nearby modes exchange character +as the geometry changes by plotting effective index and TE fraction across the +sweep, then rendering representative field profiles. + +```bash +uv run --extra dev python examples/soi_hybridization_sweep.py +``` + ## Physics @@ -96,4 +123,4 @@ around the requested effective index. The Arnoldi stage uses **shift-invert**, adaptive [Ritz-pair](https://en.wikipedia.org/wiki/Ritz_method) checkpointing, early stopping once requested modes are stable, and selective Ritz vector -reconstruction so work is spent on the modes that will actually be returned. \ No newline at end of file +reconstruction so work is spent on the modes that will actually be returned. diff --git a/benchmarks/compare_mode_solver_fixtures.py b/benchmarks/compare_mode_solver_fixtures.py index 9aaa335..70611fa 100644 --- a/benchmarks/compare_mode_solver_fixtures.py +++ b/benchmarks/compare_mode_solver_fixtures.py @@ -114,10 +114,7 @@ def main() -> None: report["cases"].append({"case_id": entry["case_id"], **status}) print(f"Inspected {len(entries)} fixture(s) from {fixture_root}") - print( - "Local validation: " - + ", ".join(f"{key}={value}" for key, value in report["summary"].items() if value) - ) + print("Local validation: " + ", ".join(f"{key}={value}" for key, value in report["summary"].items() if value)) if args.report_json is not None: args.report_json.parent.mkdir(parents=True, exist_ok=True) args.report_json.write_text(json.dumps(report, indent=2, sort_keys=True) + "\n") @@ -325,11 +322,13 @@ def _compare_local_case(root: Path, entry: dict) -> dict: _AXIS_INDEX = {"x": 0, "y": 1, "z": 2} -def _solver_edges_from_field_coords(edges: tuple[np.ndarray, np.ndarray], recipe: dict) -> tuple[np.ndarray, np.ndarray]: +def _solver_edges_from_field_coords( + edges: tuple[np.ndarray, np.ndarray], recipe: dict +) -> tuple[np.ndarray, np.ndarray]: dmin_pmc = tuple(bool(value) for value in recipe.get("dmin_pmc", (False, False))) trim_edges = tuple(recipe.get("trim_edges", ((0, 0), (0, 0)))) out = [] - for axis_edges, has_min_symmetry, (trim_start, trim_end) in zip(edges, dmin_pmc, trim_edges): + for axis_edges, has_min_symmetry, (trim_start, trim_end) in zip(edges, dmin_pmc, trim_edges, strict=True): if not has_min_symmetry: trimmed = axis_edges if trim_start or trim_end: @@ -488,7 +487,7 @@ def _eps_from_recipe( mask = np.ones(eps.shape, dtype=bool) center = box.get("center", (0.0, 0.0, 0.0)) size = box["size"] - for grid, dim in zip(grids, tangent_dims): + for grid, dim in zip(grids, tangent_dims, strict=True): axis = _AXIS_INDEX[dim] mask &= np.abs(grid - center[axis]) <= abs(size[axis]) / 2 eps_value = complex(box["eps"]) @@ -500,7 +499,7 @@ def _eps_from_recipe( center = circle.get("center", (0.0, 0.0, 0.0)) radius = float(circle["radius"]) distance_sq = np.zeros(eps.shape, dtype=float) - for grid, dim in zip(grids, tangent_dims): + for grid, dim in zip(grids, tangent_dims, strict=True): distance_sq += (grid - center[_AXIS_INDEX[dim]]) ** 2 eps[distance_sq <= radius * radius] = complex(circle["eps"]) return eps diff --git a/benchmarks/mode_solver/fixtures.py b/benchmarks/mode_solver/fixtures.py index 4fdec01..501c328 100644 --- a/benchmarks/mode_solver/fixtures.py +++ b/benchmarks/mode_solver/fixtures.py @@ -68,11 +68,7 @@ def _read_xarray_group(group: Any) -> xr.DataArray: if _XARRAY_VALUE_NAME not in group: raise KeyError(f"HDF5 group {group.name!r} is missing {_XARRAY_VALUE_NAME!r}") values = group[_XARRAY_VALUE_NAME][()] - coords = { - name: group[name][()] - for name in group.keys() - if name != _XARRAY_VALUE_NAME and hasattr(group[name], "shape") - } + coords = {name: group[name][()] for name in group if name != _XARRAY_VALUE_NAME and hasattr(group[name], "shape")} dims = _infer_dims(values.shape, coords) xarray_coords = {dim: coords[dim] for dim in dims if dim in coords} return xr.DataArray(values, dims=dims, coords=xarray_coords) @@ -80,7 +76,7 @@ def _read_xarray_group(group: Any) -> xr.DataArray: def _infer_dims(shape: tuple[int, ...], coords: dict[str, np.ndarray]) -> tuple[str, ...]: dims = tuple(dim for dim in _PREFERRED_DIMS if dim in coords) - if len(dims) == len(shape) and all(len(coords[dim]) == size for dim, size in zip(dims, shape)): + if len(dims) == len(shape) and all(len(coords[dim]) == size for dim, size in zip(dims, shape, strict=True)): return dims exact = [dim for dim in _PREFERRED_DIMS if dim in coords and len(coords[dim]) in set(shape)] diff --git a/docs/assets/hybridization_sweep.png b/docs/assets/hybridization_sweep.png new file mode 100644 index 0000000..195b7b2 Binary files /dev/null and b/docs/assets/hybridization_sweep.png differ diff --git a/docs/assets/ridge_modes.png b/docs/assets/ridge_modes.png new file mode 100644 index 0000000..88af95d Binary files /dev/null and b/docs/assets/ridge_modes.png differ diff --git a/docs/assets/tidy3d_modal_modes.png b/docs/assets/tidy3d_modal_modes.png new file mode 100644 index 0000000..2a2c661 Binary files /dev/null and b/docs/assets/tidy3d_modal_modes.png differ diff --git a/examples/README.md b/examples/README.md index 5c423bd..a5d0021 100644 --- a/examples/README.md +++ b/examples/README.md @@ -38,12 +38,22 @@ That example sweeps a 220 nm fully etched SOI ridge width and writes the effective-index and TE-fraction plots to `examples/soi_hybridization_outputs/`. -Run the README ridge-waveguide example: +Recreate the Tidy3D modal sources/monitors mode plot: + +```bash +uv run --extra dev python examples/tidy3d_modal_sources_monitors.py +``` + +That example uses the strip-waveguide geometry from the Tidy3D modal +sources/monitors notebook and writes `|Ey|` and `|Ez|` field plots for the first +three modes to `examples/tidy3d_modal_outputs/`. + +Run the README angled-slab waveguide example: ```bash uv run --extra dev python examples/ridge_waveguide_readme.py ``` -That example rasterizes a 220 nm SOI rib waveguide with a 90 nm slab, 500 nm -top ridge width, inverted angled sidewalls, and subpixel material averaging. +That example rasterizes a 220 nm silicon slab with 500 nm top width and 80 +degree angled sidewalls on an oxide substrate with air around it. It writes publication-style plots to `examples/ridge_waveguide_outputs/`. diff --git a/examples/material_grid_demos.py b/examples/material_grid_demos.py index e082784..9df4431 100644 --- a/examples/material_grid_demos.py +++ b/examples/material_grid_demos.py @@ -1,9 +1,9 @@ from __future__ import annotations import argparse +from collections.abc import Callable, Sequence from dataclasses import dataclass from pathlib import Path -from typing import Callable, Sequence import matplotlib import numpy as np @@ -15,7 +15,6 @@ import matplotlib.pyplot as plt from matplotlib.axes import Axes - FREQ_1550 = sm.C_0 / 1.55 SI_EPS = 3.48**2 SIO2_EPS = 1.44**2 @@ -167,8 +166,8 @@ def make_slot_grid() -> tuple[sm.Materials, np.ndarray]: x_edges, y_edges, xx, yy = demo_grid(nx=42, ny=30) eps = np.full(xx.shape, SIO2_EPS, dtype=np.complex128) rail = np.abs(yy) <= 0.11 - left = (-0.30 <= xx) & (xx <= -0.06) - right = (0.06 <= xx) & (xx <= 0.30) + left = (xx >= -0.30) & (xx <= -0.06) + right = (xx >= 0.06) & (xx <= 0.30) eps[rail & (left | right)] = SI_EPS return sm.Materials.from_diagonal(eps_xx=eps, x_edges=x_edges, y_edges=y_edges), eps @@ -176,8 +175,8 @@ def make_slot_grid() -> tuple[sm.Materials, np.ndarray]: def make_rib_grid() -> tuple[sm.Materials, np.ndarray]: x_edges, y_edges, xx, yy = demo_grid(nx=46, ny=30) eps = np.full(xx.shape, SIO2_EPS, dtype=np.complex128) - slab = (np.abs(xx) <= 0.72) & (-0.16 <= yy) & (yy <= -0.05) - ridge = (np.abs(xx) <= 0.28) & (-0.05 <= yy) & (yy <= 0.18) + slab = (np.abs(xx) <= 0.72) & (yy >= -0.16) & (yy <= -0.05) + ridge = (np.abs(xx) <= 0.28) & (yy >= -0.05) & (yy <= 0.18) eps[slab | ridge] = SI_EPS return sm.Materials.from_diagonal(eps_xx=eps, x_edges=x_edges, y_edges=y_edges), eps @@ -267,7 +266,7 @@ def plot_demo( else: image = draw_image(ax, dims, coords, images[column], cmap="RdBu_r", symmetric=True) plot_eps_contours(ax, coords, eps_background) - ax.set_title(f"mode {mode_index}, {column}\n" f"n_eff={n_eff.real:.5f}{n_eff.imag:+.1e}j") + ax.set_title(f"mode {mode_index}, {column}\nn_eff={n_eff.real:.5f}{n_eff.imag:+.1e}j") fig.colorbar(image, ax=ax, fraction=0.046, pad=0.04) fig.suptitle(f"{demo.title}: {demo.description}", fontsize=14) @@ -296,7 +295,7 @@ def electric_magnitude_image( magnitude_squared = np.abs(ex) ** 2 for component in ("Ey", "Ez"): other_dims, other_coords, values = component_image(data, component, mode_index) - coords_match = all(len(a) == len(b) and np.allclose(a, b) for a, b in zip(coords, other_coords)) + coords_match = all(len(a) == len(b) and np.allclose(a, b) for a, b in zip(coords, other_coords, strict=True)) if other_dims != dims or not coords_match: raise ValueError("field components are not colocated on a common plotting grid") magnitude_squared += np.abs(values) ** 2 @@ -327,7 +326,10 @@ def draw_image( limit = max(max_abs, np.finfo(float).eps) kwargs = {"vmin": -limit, "vmax": limit} else: - kwargs = {"vmin": float(np.nanmin(plot_values)), "vmax": max(float(np.nanmax(plot_values)), np.finfo(float).eps)} + kwargs = { + "vmin": float(np.nanmin(plot_values)), + "vmax": max(float(np.nanmax(plot_values)), np.finfo(float).eps), + } image = ax.imshow( plot_values.T, extent=extent, diff --git a/examples/ridge_waveguide_readme.py b/examples/ridge_waveguide_readme.py index c32e437..0d1f39c 100644 --- a/examples/ridge_waveguide_readme.py +++ b/examples/ridge_waveguide_readme.py @@ -1,8 +1,8 @@ -"""Ridge waveguide example used by the README. +"""Angled slab waveguide example used by the README. -The geometry is a compact SOI rib waveguide. It is MicroMode-native: the -structure is rasterized into a material grid first, then the solver receives -only arrays. +The geometry is a compact silicon slab with 80 degree sidewalls on an oxide +substrate and air around it. It is MicroMode-native: the structure is rasterized +into a material grid first, then the solver receives only arrays. """ from __future__ import annotations @@ -19,7 +19,6 @@ import micromode as mm - WAVELENGTH_UM = 1.55 N_SI = 3.476 N_SIO2 = 1.444 @@ -47,7 +46,7 @@ def publication_style() -> dict[str, object]: def main() -> None: - parser = argparse.ArgumentParser(description="Solve and plot a rasterized ridge waveguide.") + parser = argparse.ArgumentParser(description="Solve and plot a rasterized angled slab waveguide.") parser.add_argument("--step", type=float, default=0.02, help="Grid step in microns.") parser.add_argument("--subpixels", type=int, default=7, help="Subpixel samples per axis for material averaging.") parser.add_argument("--output-dir", type=Path, default=Path("examples/ridge_waveguide_outputs")) @@ -70,31 +69,26 @@ def main() -> None: def ridge_waveguide_materials(step: float = 0.02, subpixels: int = 7) -> tuple[mm.Materials, np.ndarray]: - """Rasterize an SOI rib waveguide with subpixel material averaging.""" + """Rasterize a silicon slab with 80 degree sidewalls and subpixel averaging.""" - film_thickness = 0.22 - slab_height = 0.09 + slab_thickness = 0.22 top_width = 0.50 - sidewall_angle_deg = 82.0 + sidewall_angle_deg = 80.0 substrate_height = 1.0 clad_height = 0.8 domain_width = 3.0 x_edges = centered_edges(width=domain_width, step=step) - y_edges = offset_edges(lower=-substrate_height, upper=film_thickness + clad_height, step=step) + y_edges = offset_edges(lower=-substrate_height, upper=slab_thickness + clad_height, step=step) sample_x, sample_y = subpixel_centers(x_edges, y_edges, subpixels=subpixels) eps_samples = np.full(sample_x.shape, N_AIR**2, dtype=np.complex128) eps_samples[sample_y < 0.0] = N_SIO2**2 - silicon = (sample_y >= 0.0) & (sample_y < slab_height) - - ridge_height = film_thickness - slab_height - bottom_extra = ridge_height / np.tan(np.deg2rad(sidewall_angle_deg)) - vertical_fraction = np.clip((sample_y - slab_height) / ridge_height, 0.0, 1.0) + bottom_extra = slab_thickness / np.tan(np.deg2rad(sidewall_angle_deg)) + vertical_fraction = np.clip(sample_y / slab_thickness, 0.0, 1.0) half_width = 0.5 * top_width + (1.0 - vertical_fraction) * bottom_extra - ridge = (sample_y >= slab_height) & (sample_y < film_thickness) & (np.abs(sample_x) <= half_width) - silicon |= ridge + silicon = (sample_y >= 0.0) & (sample_y < slab_thickness) & (np.abs(sample_x) <= half_width) eps_samples[silicon] = N_SI**2 materials = mm.Materials.from_subpixel_diagonal( @@ -108,12 +102,22 @@ def ridge_waveguide_materials(step: float = 0.02, subpixels: int = 7) -> tuple[m def centered_edges(*, width: float, step: float) -> np.ndarray: - cells = int(round(width / step)) + if step <= 0.0: + raise ValueError("step must be positive") + cells = round(width / step) + if cells < 1: + raise ValueError("width must be at least one step") return np.linspace(-0.5 * width, 0.5 * width, cells + 1) def offset_edges(*, lower: float, upper: float, step: float) -> np.ndarray: - cells = int(round((upper - lower) / step)) + if step <= 0.0: + raise ValueError("step must be positive") + if upper <= lower: + raise ValueError("upper must be greater than lower") + cells = round((upper - lower) / step) + if cells < 1: + raise ValueError("upper - lower must be at least one step") return np.linspace(lower, upper, cells + 1) @@ -153,7 +157,7 @@ def plot_index(materials: mm.Materials, eps: np.ndarray, path: Path) -> None: ) draw_material_outline(ax, x, y, eps, color="black", linewidth=1.3) draw_material_outline(ax, x, y, eps, color="white", linewidth=0.55) - ax.set_title("SOI rib waveguide") + ax.set_title("Angled slab on substrate") ax.set_xlabel("x (um)") ax.set_ylabel("y (um)") ax.set_xlim(-0.9, 0.9) @@ -172,7 +176,7 @@ def plot_modes(materials: mm.Materials, eps: np.ndarray, data: mm.Result, path: with plt.rc_context(publication_style()): fig, axes = plt.subplots(2, 2, figsize=(4.9, 2.65), constrained_layout=False, sharex=True, sharey=True) fig.subplots_adjust(left=0.02, right=0.995, bottom=0.03, top=0.89, wspace=0.06, hspace=-0.04) - for ax, component in zip(axes.ravel(), components): + for ax, component in zip(axes.ravel(), components, strict=True): field = data.field_components[component].isel(f=0, mode_index=0) values = normalize_signed(np.asarray(field.values).squeeze().real) ax.imshow( @@ -199,60 +203,43 @@ def plot_modes(materials: mm.Materials, eps: np.ndarray, data: mm.Result, path: def plot_readme_figure(materials: mm.Materials, eps: np.ndarray, data: mm.Result, path: Path) -> None: - """Create a compact presentation figure for the README.""" + """Create a 5:3 presentation figure showing fundamental field components.""" + style = publication_style() + style["savefig.bbox"] = None x_edges = np.asarray(materials.grid.x_edges, dtype=float) y_edges = np.asarray(materials.grid.y_edges, dtype=float) x = 0.5 * (x_edges[:-1] + x_edges[1:]) y = 0.5 * (y_edges[:-1] + y_edges[1:]) extent = (x_edges[0], x_edges[-1], y_edges[0], y_edges[-1]) - panels = ( - ("permittivity", "eps", 0), - ("mode 0 |E|", None, 0), - ("mode 0 Ex", "Ex", 0), - ("mode 1 |E|", None, 1), - ("mode 1 Ex", "Ex", 1), - ) - - with plt.rc_context(publication_style()): - fig, axes = plt.subplots(1, len(panels), figsize=(7.6, 1.95), constrained_layout=False, sharex=True, sharey=True) - fig.subplots_adjust(left=0.035, right=0.995, bottom=0.12, top=0.82, wspace=0.10) - for ax, (title, component, mode_index) in zip(axes, panels): - values, cmap, vmin, vmax = readme_panel_values(eps, data, component=component, mode_index=mode_index) + components = ("Ex", "Ey", "Ez", "Hx", "Hy", "Hz") + n_eff = complex(np.asarray(data.n_complex.values)[0, 0]) + + with plt.rc_context(style): + fig, axes = plt.subplots(3, 2, figsize=(7.5, 4.5), constrained_layout=False, sharex=True, sharey=True) + fig.subplots_adjust(left=0.025, right=0.995, bottom=0.035, top=0.90, wspace=0.04, hspace=0.20) + fig.suptitle(f"Fundamental mode, n={format_complex(n_eff, precision=4)}", y=0.985) + for ax, component in zip(axes.ravel(), components, strict=True): + values = np.asarray(data.field_components[component].isel(f=0, mode_index=0).values).squeeze() + plot_values = normalize_signed(values.real) ax.imshow( - values.T, + plot_values.T, origin="lower", extent=extent, aspect="equal", - cmap=cmap, - vmin=vmin, - vmax=vmax, + cmap="RdBu_r", + vmin=-1.0, + vmax=1.0, interpolation="bicubic", ) - if component == "eps": - draw_material_outline(ax, x, y, eps, color="black", linewidth=1.2) - draw_material_outline(ax, x, y, eps, color="white", linewidth=0.5) - else: - outline_color = "white" if cmap == "magma" else "#111827" - draw_material_outline(ax, x, y, eps, color=outline_color, linewidth=0.9) - ax.set_title(title) + draw_material_outline(ax, x, y, eps, color="white", linewidth=0.9) + draw_material_outline(ax, x, y, eps, color="#111827", linewidth=0.35) + ax.set_title(f"Re({component})", pad=3) ax.set_xlim(-0.9, 0.9) ax.set_ylim(-0.28, 0.42) ax.set_xticks([]) ax.set_yticks([]) ax.tick_params(bottom=False, left=False, labelbottom=False, labelleft=False) - if component != "eps": - ax.text( - 0.06, - 0.90, - f"n={data.n_eff.values[0, mode_index]:.4f}", - transform=ax.transAxes, - ha="left", - va="top", - color="white" if cmap == "magma" else "black", - alpha=0.5, - fontsize=8, - ) save_figure(fig, path) plt.close(fig) @@ -266,27 +253,9 @@ def draw_material_outline(ax, x: np.ndarray, y: np.ndarray, eps: np.ndarray, **k ax.contour(x, y, eps.real.T, levels=levels, **kwargs) -def readme_panel_values( - eps: np.ndarray, - data: mm.Result, - *, - component: str | None, - mode_index: int, -) -> tuple[np.ndarray, str, float, float]: - if component == "eps": - return eps.real, "Greys", N_AIR**2, N_SI**2 - if component is None: - values = np.sqrt( - sum( - np.abs(np.asarray(data.field_components[name].isel(f=0, mode_index=mode_index).values).squeeze()) ** 2 - for name in ("Ex", "Ey", "Ez") - ) - ) - scale = max(float(np.nanmax(values)), np.finfo(float).eps) - return values / scale, "magma", 0.0, 1.0 - - values = np.asarray(data.field_components[component].isel(f=0, mode_index=mode_index).values).squeeze().real - return normalize_signed(values), "RdBu_r", -1.0, 1.0 +def format_complex(value: complex, *, precision: int) -> str: + sign = "+" if value.imag >= 0.0 else "-" + return f"{value.real:.{precision}f}{sign}{abs(value.imag):.{precision}g}i" def normalize_signed(values: np.ndarray) -> np.ndarray: diff --git a/examples/soi_hybridization_sweep.py b/examples/soi_hybridization_sweep.py index 6a6ecfe..a650785 100644 --- a/examples/soi_hybridization_sweep.py +++ b/examples/soi_hybridization_sweep.py @@ -13,8 +13,6 @@ import matplotlib.pyplot as plt from matplotlib.lines import Line2D -from matplotlib.patches import Patch - WAVELENGTH_UM = 1.55 FREQ_1550 = mm.C_0 / WAVELENGTH_UM @@ -118,7 +116,7 @@ def soi_ridge_materials( ) -> tuple[mm.Materials, np.ndarray]: x = 0.5 * (x_edges[:-1] + x_edges[1:]) y = 0.5 * (y_edges[:-1] + y_edges[1:]) - xx, yy = np.meshgrid(x, y, indexing="ij") + _xx, yy = np.meshgrid(x, y, indexing="ij") eps = np.where(yy < 0.0, N_SIO2**2, N_AIR**2).astype(np.complex128) silicon_fill = rectangle_fill_fraction( x_edges=x_edges, @@ -164,14 +162,8 @@ def write_summary(path: Path, sweep: mm.Sweep, args: argparse.Namespace) -> Path "y_step_um": args.y_step, "width_um": sweep.values.tolist(), "n_eff": sweep.n_eff.tolist(), - "pol_fraction": { - key: value.tolist() - for key, value in sweep.pol_fraction.items() - }, - "pol_fraction_waveguide": { - key: value.tolist() - for key, value in sweep.pol_fraction_waveguide.items() - }, + "pol_fraction": {key: value.tolist() for key, value in sweep.pol_fraction.items()}, + "pol_fraction_waveguide": {key: value.tolist() for key, value in sweep.pol_fraction_waveguide.items()}, } path.write_text(json.dumps(payload, indent=2), encoding="utf-8") return path @@ -260,8 +252,6 @@ def plot_sweep(path: Path, sweep: mm.Sweep) -> None: handlelength=1.8, ) - add_panel_label(axes[0], "a") - add_panel_label(axes[1], "b") save_figure(fig, path) plt.close(fig) @@ -462,7 +452,7 @@ def electric_magnitude_image( magnitude_squared = np.abs(ex) ** 2 for component in ("Ey", "Ez"): other_dims, other_coords, values = component_image(result, component, mode_index) - coords_match = all(len(a) == len(b) and np.allclose(a, b) for a, b in zip(coords, other_coords)) + coords_match = all(len(a) == len(b) and np.allclose(a, b) for a, b in zip(coords, other_coords, strict=True)) if other_dims != dims or not coords_match: raise ValueError("field components are not colocated") magnitude_squared += np.abs(values) ** 2 @@ -526,19 +516,6 @@ def plot_eps_contours(ax, coords: tuple[np.ndarray, np.ndarray], eps: np.ndarray ax.contour(x, y, values.T, levels=[level], colors="white", linewidths=0.9, alpha=0.95) -def add_panel_label(ax, label: str) -> None: - ax.text( - -0.12, - 1.08, - label, - transform=ax.transAxes, - ha="left", - va="top", - fontsize=10, - fontweight="bold", - ) - - def save_figure(fig, path: Path) -> None: fig.savefig(path) fig.savefig(path.with_suffix(".pdf")) diff --git a/examples/tidy3d_modal_sources_monitors.py b/examples/tidy3d_modal_sources_monitors.py new file mode 100644 index 0000000..887a773 --- /dev/null +++ b/examples/tidy3d_modal_sources_monitors.py @@ -0,0 +1,170 @@ +"""Recreate the Tidy3D modal source/monitor mode plot with MicroMode. + +Reference setup: +https://www.flexcompute.com/tidy3d/examples/notebooks/ModalSourcesMonitors/ + +The Tidy3D notebook solves modes of an x-propagating silicon strip waveguide on +a silica substrate. In MicroMode, this is represented by setting normal_axis=0: +the two material-grid coordinates are global y and z, and the returned global +field components are labeled Ey and Ez as in the Tidy3D plot. +""" + +from __future__ import annotations + +import argparse +from pathlib import Path + +import matplotlib +import numpy as np + +import micromode as mm + +matplotlib.use("Agg") + +import matplotlib.pyplot as plt + +WAVELENGTH_UM = 1.55 +WG_HEIGHT_UM = 0.22 +WG_WIDTH_UM = 0.45 +N_SI = 3.48 +N_SIO2 = 1.45 +N_AIR = 1.0 + + +def main() -> None: + args = parse_args() + args.output_dir.mkdir(parents=True, exist_ok=True) + + materials, eps = tidy3d_waveguide_materials(y_step=args.y_step, z_step=args.z_step) + data = mm.solve_modes( + material_grid=materials, + wavelength=WAVELENGTH_UM, + num_modes=3, + target_neff=2.35, + components=("Ey", "Ez"), + krylov_dim=args.krylov_dim, + ) + + print("n_eff:", np.round(data.n_eff.values[0], 6)) + output_path = args.output_dir / "tidy3d_modal_modes.png" + plot_tidy3d_mode_fields(materials, eps, data, output_path) + print(f"Wrote {output_path}") + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description="Recreate the Tidy3D modal sources/monitors mode plot.") + parser.add_argument( + "--output-dir", + type=Path, + default=Path(__file__).resolve().parent / "tidy3d_modal_outputs", + ) + parser.add_argument("--y-step", type=float, default=0.025, help="Mode-plane y grid step in microns.") + parser.add_argument("--z-step", type=float, default=0.025, help="Mode-plane z grid step in microns.") + parser.add_argument("--krylov-dim", type=int, default=64) + return parser.parse_args() + + +def tidy3d_waveguide_materials(*, y_step: float, z_step: float) -> tuple[mm.Materials, np.ndarray]: + """Rasterize the Tidy3D strip waveguide mode plane. + + The Tidy3D mode plane has size [0, 3, 2] in x, y, z. The substrate fills + z < 0, and the silicon strip spans |y| <= 0.225 um and 0 <= z <= 0.22 um. + """ + + y_edges = centered_edges(width=3.0, step=y_step) + z_edges = centered_edges(width=2.0, step=z_step) + y = 0.5 * (y_edges[:-1] + y_edges[1:]) + z = 0.5 * (z_edges[:-1] + z_edges[1:]) + _yy, zz = np.meshgrid(y, z, indexing="ij") + + eps = np.where(zz < 0.0, N_SIO2**2, N_AIR**2).astype(np.complex128) + waveguide = (np.abs(_yy) <= 0.5 * WG_WIDTH_UM) & (zz >= 0.0) & (zz <= WG_HEIGHT_UM) + eps[waveguide] = N_SI**2 + + materials = mm.Materials.from_diagonal( + eps_xx=eps, + x_edges=y_edges, + y_edges=z_edges, + normal_axis=0, + ) + return materials, eps + + +def centered_edges(*, width: float, step: float) -> np.ndarray: + if step <= 0.0: + raise ValueError("step must be positive") + cells = round(width / step) + if cells < 1: + raise ValueError("width must be at least one step") + return np.linspace(-0.5 * width, 0.5 * width, cells + 1) + + +def plot_tidy3d_mode_fields(materials: mm.Materials, eps: np.ndarray, data: mm.Result, path: Path) -> None: + y_edges = np.asarray(materials.grid.x_edges, dtype=float) + z_edges = np.asarray(materials.grid.y_edges, dtype=float) + y = 0.5 * (y_edges[:-1] + y_edges[1:]) + z = 0.5 * (z_edges[:-1] + z_edges[1:]) + extent = (y_edges[0], y_edges[-1], z_edges[0], z_edges[-1]) + + with plt.rc_context(publication_style()): + fig, axes = plt.subplots(3, 2, figsize=(12.0, 12.0), constrained_layout=True, sharex=True, sharey=True) + for mode_index in range(3): + for col, component in enumerate(("Ey", "Ez")): + ax = axes[mode_index, col] + values = field_abs(data, component=component, mode_index=mode_index) + image = ax.imshow( + values.T, + origin="lower", + extent=extent, + aspect="equal", + cmap="magma", + vmin=0.0, + vmax=float(np.nanpercentile(values, 99.5)), + interpolation="nearest", + ) + draw_material_context(ax, y, z, eps) + ax.set_title(f"{component}, mode_index={mode_index}") + ax.set_xlabel("y (um)") + ax.set_ylabel("z (um)") + ax.set_xlim(extent[0], extent[1]) + ax.set_ylim(extent[2], extent[3]) + cbar = fig.colorbar(image, ax=ax, extend="both", fraction=0.046, pad=0.04) + cbar.set_label(f"|{component}|") + save_figure(fig, path) + plt.close(fig) + + +def field_abs(data: mm.Result, *, component: str, mode_index: int) -> np.ndarray: + field = data.field_components[component].isel(f=0, mode_index=mode_index).squeeze(drop=True) + return np.abs(np.asarray(field.transpose("y", "z").values)) + + +def draw_material_context(ax, y: np.ndarray, z: np.ndarray, eps: np.ndarray) -> None: + ax.axhline(0.0, color="#111111", linewidth=0.85, alpha=0.35) + silicon_level = 0.5 * (N_AIR**2 + N_SI**2) + ax.contour(y, z, eps.real.T, levels=[silicon_level], colors="#2f2f2f", linewidths=1.0, alpha=0.75) + + +def publication_style() -> dict[str, object]: + return { + "font.family": "DejaVu Sans", + "font.size": 10, + "axes.titlesize": 12, + "axes.labelsize": 10, + "xtick.labelsize": 10, + "ytick.labelsize": 10, + "figure.dpi": 120, + "savefig.dpi": 220, + "savefig.bbox": "tight", + "savefig.pad_inches": 0.03, + } + + +def save_figure(fig, path: Path) -> None: + fig.savefig(path) + fig.savefig(path.with_suffix(".pdf")) + fig.savefig(path.with_suffix(".svg")) + + +if __name__ == "__main__": + main() diff --git a/pyproject.toml b/pyproject.toml index 75637d0..a0b9d69 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "micromode" -version = "0.1.0a3" +version = "0.1.0a4" description = "Rust-backed photonics mode solver with a small Python API." readme = "README.md" requires-python = ">=3.10,<3.14" @@ -50,10 +50,12 @@ dev = [ "maturin>=1.7,<2", "packaging>=24.2", "pkginfo>=1.12.1.2", + "pyright>=1.1.390,<2", "pytest>=8.1,<10.0.0", "pytest-cov>=5,<8", "pytest-xdist>=3.6,<4.0", - "tomli>=2,<3; python_version < '3.11'", + "ruff>=0.8.0,<1", + "tomli>=2,<3", "twine>=5,<7", ] @@ -79,3 +81,18 @@ markers = [ filterwarnings = ["ignore::DeprecationWarning"] testpaths = ["tests"] python_files = ["test_*.py"] + +[tool.ruff] +line-length = 120 +target-version = "py310" +src = ["python", "tests", "scripts", "benchmarks", "examples"] +extend-exclude = ["dist", "tmp"] + +[tool.ruff.lint] +select = ["E", "F", "I", "B", "UP", "SIM", "RUF"] + +[tool.pyright] +include = ["python/micromode", "scripts"] +exclude = ["**/__pycache__", ".venv", "dist", "tmp"] +pythonVersion = "3.10" +typeCheckingMode = "basic" diff --git a/python/micromode/_rust.py b/python/micromode/_rust.py index d96409f..f713f71 100644 --- a/python/micromode/_rust.py +++ b/python/micromode/_rust.py @@ -9,11 +9,14 @@ C_0 = 2.997_924_58e14 EPSILON_0 = 8.854187812800384e-18 +SparseSolveResult = tuple[np.ndarray, list[np.ndarray], dict[str, object]] try: from ._core import ( solve_diagonal_sparse_py as _solve_diagonal_sparse, + ) + from ._core import ( solve_tensorial_sparse_py as _solve_tensorial_sparse, ) except Exception: # pragma: no cover - exercised when extension is not built locally. @@ -38,7 +41,7 @@ def solve_diagonal_sparse( dmin_pmc: tuple[bool, bool] = (False, False), krylov_dim: int = 32, initial_vector: np.ndarray | None = None, -) -> tuple[np.ndarray, list[np.ndarray]]: +) -> SparseSolveResult: """Run the Rust sparse diagonal mode solver for a prepared 2D Yee grid.""" # This is the normal path for diagonal material grids. The returned fields # are still flattened by mode/component; `raster.py` restores grid axes. @@ -92,14 +95,18 @@ def solve_diagonal_sparse( ) n_complex = _pairs_to_complex(n_pairs) fields = [_pairs_to_complex(component) for component in field_pairs] - return n_complex, fields, _solver_info( - residuals, - power_norms, - lorentz_norms, - lorentz_orthogonality_error, - backend, - operator_size, - operator_nnz, + return ( + n_complex, + fields, + _solver_info( + residuals, + power_norms, + lorentz_norms, + lorentz_orthogonality_error, + backend, + operator_size, + operator_nnz, + ), ) @@ -120,7 +127,7 @@ def solve_tensorial_sparse( dmin_pmc: tuple[bool, bool] = (False, False), krylov_dim: int = 32, initial_vector: np.ndarray | None = None, -) -> tuple[np.ndarray, list[np.ndarray]]: +) -> SparseSolveResult: """Run the Rust sparse full-tensor mode solver for a prepared 2D Yee grid.""" # Full tensor material grids, angled solves, and bent solves use this path # because coordinate transforms can introduce off-diagonal eps/mu terms. @@ -174,14 +181,18 @@ def solve_tensorial_sparse( ) n_complex = _pairs_to_complex(n_pairs) fields = [_pairs_to_complex(component) for component in field_pairs] - return n_complex, fields, _solver_info( - residuals, - power_norms, - lorentz_norms, - lorentz_orthogonality_error, - backend, - operator_size, - operator_nnz, + return ( + n_complex, + fields, + _solver_info( + residuals, + power_norms, + lorentz_norms, + lorentz_orthogonality_error, + backend, + operator_size, + operator_nnz, + ), ) diff --git a/python/micromode/models.py b/python/micromode/models.py index 7c39074..135f8d0 100644 --- a/python/micromode/models.py +++ b/python/micromode/models.py @@ -2,8 +2,9 @@ from __future__ import annotations +from collections.abc import Sequence from dataclasses import dataclass -from typing import Literal, Sequence +from typing import Literal import numpy as np @@ -38,7 +39,7 @@ def __post_init__(self) -> None: object.__setattr__(self, "order", int(self.order)) @classmethod - def from_num_cells(cls, num_cells: tuple[int, int]) -> "PmlSpec": + def from_num_cells(cls, num_cells: tuple[int, int]) -> PmlSpec: return cls(num_cells=num_cells) def as_dict(self) -> dict[str, float | int | tuple[int, int]]: @@ -50,6 +51,14 @@ def as_dict(self) -> dict[str, float | int | tuple[int, int]]: "order": self.order, } + def profile_dict(self) -> dict[str, float | int]: + return { + "sigma_max": self.sigma_max, + "kappa_min": self.kappa_min, + "kappa_max": self.kappa_max, + "order": self.order, + } + @dataclass(frozen=True) class BoundarySpec: @@ -68,11 +77,11 @@ def __post_init__(self) -> None: @property def dmin_pmc(self) -> tuple[bool, bool]: - return tuple(value == "pmc" for value in self.low) + return self.low[0] == "pmc", self.low[1] == "pmc" @property def dmin_pml(self) -> tuple[bool, bool]: - return tuple(value == "pec" for value in self.low) + return self.low[0] == "pec", self.low[1] == "pec" def as_dict(self) -> dict[str, tuple[str, str]]: return {"low": self.low} @@ -155,7 +164,7 @@ def from_diagonal( mu_zz: np.ndarray | None = None, normal_axis: Literal[0, 1, 2] = 2, normal_coordinate: float = 0.0, - ) -> "Materials": + ) -> Materials: grid = Grid( tuple(float(value) for value in x_edges), tuple(float(value) for value in y_edges), @@ -170,7 +179,9 @@ def from_diagonal( np.ones(grid.shape, dtype=np.complex128) if mu_yy is None else mu_yy, np.ones(grid.shape, dtype=np.complex128) if mu_zz is None else mu_zz, ) - return cls(grid=grid, eps_tensor=_diagonal_to_full_tensor(eps_diag), mu_tensor=_diagonal_to_full_tensor(mu_diag)) + return cls( + grid=grid, eps_tensor=_diagonal_to_full_tensor(eps_diag), mu_tensor=_diagonal_to_full_tensor(mu_diag) + ) @classmethod def from_components( @@ -198,7 +209,7 @@ def from_components( mu_zy: np.ndarray | None = None, normal_axis: Literal[0, 1, 2] = 2, normal_coordinate: float = 0.0, - ) -> "Materials": + ) -> Materials: grid = Grid( tuple(float(value) for value in x_edges), tuple(float(value) for value in y_edges), @@ -267,7 +278,7 @@ def from_slice( mu_zy: np.ndarray | None = None, normal_axis: Literal[0, 1, 2] = 2, normal_coordinate: float = 0.0, - ) -> "Materials": + ) -> Materials: """Build a one-dimensional mode slice. ``axis="x"`` makes the supplied material arrays vary along the first @@ -299,10 +310,13 @@ def expand(label: str, values: np.ndarray | None) -> np.ndarray | None: raise ValueError(f"{label} must have shape {(cell_count,)} for a one-dimensional slice") return array[:, None] if axis_index == 0 else array[None, :] + expanded_eps_xx = expand("eps_xx", eps_xx) + if expanded_eps_xx is None: + raise ValueError("eps_xx is required") return cls.from_components( x_edges=x_edges, y_edges=y_edges, - eps_xx=expand("eps_xx", eps_xx), + eps_xx=expanded_eps_xx, eps_yy=expand("eps_yy", eps_yy), eps_zz=expand("eps_zz", eps_zz), eps_xy=expand("eps_xy", eps_xy), @@ -340,7 +354,7 @@ def from_subpixel_diagonal( average: AverageMethod = "arithmetic", normal_axis: Literal[0, 1, 2] = 2, normal_coordinate: float = 0.0, - ) -> "Materials": + ) -> Materials: """Build a grid from higher-resolution samples inside each solver cell. The sample arrays may either have shape ``(nx * sx, ny * sy)`` or @@ -406,10 +420,7 @@ def average_subpixels( elif array.shape == (nx, ny, sx, sy): grouped = array else: - raise ValueError( - "subpixel values must have shape " - f"{(nx * sx, ny * sy)} or {(nx, ny, sx, sy)}" - ) + raise ValueError(f"subpixel values must have shape {(nx * sx, ny * sy)} or {(nx, ny, sx, sy)}") if method == "arithmetic": return grouped.mean(axis=(2, 3)) @@ -435,16 +446,21 @@ def shape(self) -> tuple[int, int]: def is_diagonal(self) -> bool: off_diagonal = np.ones((3, 3), dtype=bool) np.fill_diagonal(off_diagonal, False) + mu_tensor = self._resolved_mu_tensor() return bool( - np.all(np.abs(self.eps_tensor[off_diagonal]) <= 1e-12) - and np.all(np.abs(self.mu_tensor[off_diagonal]) <= 1e-12) + np.all(np.abs(self.eps_tensor[off_diagonal]) <= 1e-12) and np.all(np.abs(mu_tensor[off_diagonal]) <= 1e-12) ) def flat_eps_tensor(self) -> np.ndarray: return self.eps_tensor.reshape(3, 3, -1) def flat_mu_tensor(self) -> np.ndarray: - return self.mu_tensor.reshape(3, 3, -1) + return self._resolved_mu_tensor().reshape(3, 3, -1) + + def _resolved_mu_tensor(self) -> np.ndarray: + if self.mu_tensor is None: + raise RuntimeError("mu_tensor was not initialized") + return self.mu_tensor def diagonal_eps(self) -> np.ndarray: return np.stack([self.eps_tensor[axis, axis] for axis in range(3)], axis=0) diff --git a/python/micromode/raster.py b/python/micromode/raster.py index 559e4b4..deb818d 100644 --- a/python/micromode/raster.py +++ b/python/micromode/raster.py @@ -2,13 +2,14 @@ from __future__ import annotations -from typing import Literal, Sequence +from collections.abc import Sequence +from typing import Literal, cast import numpy as np import xarray as xr from ._rust import C_0, solve_diagonal_sparse, solve_tensorial_sparse -from .models import BoundarySpec, Materials, PmlSpec, SliceAxis, Spec +from .models import BoundaryCondition, BoundarySpec, Materials, PmlSpec, SliceAxis, Spec from .result import Result _COMPONENTS = ("Ex", "Ey", "Ez", "Hx", "Hy", "Hz") @@ -330,7 +331,7 @@ def _solve_one_frequency( # Forward and backward spacings represent the local Yee grid. The derivative # builders need both because E and H components are staggered. dlf = (np.diff(x_edges), np.diff(y_edges)) - dlb = tuple(_dual_steps(steps) for steps in dlf) + dlb = (_dual_steps(dlf[0]), _dual_steps(dlf[1])) eps_tensor = material_grid.flat_eps_tensor() mu_tensor = material_grid.flat_mu_tensor() if target_neff is None: @@ -345,7 +346,7 @@ def _solve_one_frequency( # diagonal grid may become full tensor after this step. eps_tensor, mu_tensor = _transformed_material_tensors( material_grid.eps_tensor, - material_grid.mu_tensor, + material_grid._resolved_mu_tensor(), x_edges=x_edges, y_edges=y_edges, angle_theta=angle_theta, @@ -442,17 +443,21 @@ def _solve_one_frequency_rust_sparse( derivative_scale=C_0 / (2 * np.pi * freq), omega=2 * np.pi * freq, num_pml=pml_spec.num_cells, - pml_profile=pml_spec.as_dict(), + pml_profile=pml_spec.profile_dict(), dmin_pml=boundary_spec.dmin_pml, dmin_pmc=boundary_spec.dmin_pmc, krylov_dim=actual_krylov_dim, initial_vector=_default_initial_vector(2 * nx * ny, shape=(nx, ny)), ) - return n_complex, _fields_to_grid(fields, (nx, ny)), _solver_info_with_context( - solver_info, - backend_kind="diagonal_sparse", - shape=(nx, ny), - krylov_dim=actual_krylov_dim, + return ( + n_complex, + _fields_to_grid(fields, (nx, ny)), + _solver_info_with_context( + solver_info, + backend_kind="diagonal_sparse", + shape=(nx, ny), + krylov_dim=actual_krylov_dim, + ), ) @@ -484,17 +489,21 @@ def _solve_one_frequency_rust_tensorial_sparse( derivative_scale=C_0 / (2 * np.pi * freq), omega=2 * np.pi * freq, num_pml=pml_spec.num_cells, - pml_profile=pml_spec.as_dict(), + pml_profile=pml_spec.profile_dict(), dmin_pml=boundary_spec.dmin_pml, dmin_pmc=boundary_spec.dmin_pmc, krylov_dim=actual_krylov_dim, initial_vector=_default_initial_vector(4 * nx * ny, shape=(nx, ny)), ) - return n_complex, _fields_to_grid(fields, (nx, ny)), _solver_info_with_context( - solver_info, - backend_kind="tensorial_sparse", - shape=(nx, ny), - krylov_dim=actual_krylov_dim, + return ( + n_complex, + _fields_to_grid(fields, (nx, ny)), + _solver_info_with_context( + solver_info, + backend_kind="tensorial_sparse", + shape=(nx, ny), + krylov_dim=actual_krylov_dim, + ), ) @@ -582,7 +591,7 @@ def _resolve_boundary_spec( return BoundarySpec() if isinstance(boundary, BoundarySpec): return boundary - return BoundarySpec(low=boundary) + return BoundarySpec(low=cast(tuple[BoundaryCondition, BoundaryCondition], boundary)) def _solver_info_with_context( @@ -652,7 +661,7 @@ def _fields_to_grid(fields: list[np.ndarray], shape: tuple[int, int]) -> dict[st nx, ny = shape mode_count = fields[0].shape[0] out = {} - for component, values_by_mode in zip(_COMPONENTS, fields): + for component, values_by_mode in zip(_COMPONENTS, fields, strict=True): values = np.asarray(values_by_mode).reshape(mode_count, nx, ny) out[component] = np.moveaxis(values, 0, -1) return out @@ -669,14 +678,17 @@ def _local_fields_to_global(fields: dict[str, np.ndarray], *, normal_axis: int) axis_names = ("x", "y", "z") if normal_axis not in {0, 1, 2}: raise ValueError("normal_axis must be 0, 1, or 2") - local_to_global = tuple(axis for axis in range(3) if axis != normal_axis) + (normal_axis,) + local_to_global = (*(axis for axis in range(3) if axis != normal_axis), normal_axis) out: dict[str, np.ndarray] = {} for prefix in ("E", "H"): for local_axis, global_axis in enumerate(local_to_global): local_name = f"{prefix}{axis_names[local_axis]}" global_name = f"{prefix}{axis_names[global_axis]}" if local_name in fields: - out[global_name] = fields[local_name] + # For y-normal planes the local (x, z, y) order is left-handed; + # flipping the local-y tangential component restores physical +y flux. + sign = -1.0 if normal_axis == 1 and local_axis == 1 else 1.0 + out[global_name] = sign * fields[local_name] return out @@ -736,6 +748,7 @@ def _default_initial_vector(size: int, shape: tuple[int, int] | None = None) -> vector[0, :, :] = 0 if ny > 1: vector[:, 0, :] = 0 - return np.vstack(vector).flatten("F") + stacked = np.concatenate(tuple(vector[ix, :, :] for ix in range(nx)), axis=0) + return stacked.flatten("F") index = np.arange(1, size + 1, dtype=float) return np.sin(0.37 * index) + 1j * np.cos(0.53 * index) diff --git a/python/micromode/result.py b/python/micromode/result.py index e1431a4..7b3c04c 100644 --- a/python/micromode/result.py +++ b/python/micromode/result.py @@ -3,10 +3,11 @@ from __future__ import annotations import json +from collections.abc import Iterable from dataclasses import dataclass from functools import cached_property from pathlib import Path -from typing import Any, Iterable +from typing import Any import numpy as np import xarray as xr @@ -149,13 +150,13 @@ def plot_field( data_array = data_array.transpose(*[dim for dim in _SPATIAL_DIMS if dim in data_array.dims]) values = np.asarray(data_array.values).squeeze() if val == "real": - values = values.real + values = np.real(values) default_cmap = "RdBu_r" elif val in {"abs", "magnitude"}: values = np.abs(values) default_cmap = "magma" elif val == "imag": - values = values.imag + values = np.imag(values) default_cmap = "RdBu_r" else: raise ValueError("val must be one of 'real', 'imag', or 'abs'") @@ -208,7 +209,7 @@ def plot_field_components( nrows = int(np.ceil(len(available) / ncols)) fig, axes = plt.subplots(nrows, ncols, squeeze=False, figsize=(4.0 * ncols, 3.2 * nrows)) flat_axes = axes.ravel() - for ax, component in zip(flat_axes, available): + for ax, component in zip(flat_axes, available, strict=False): self.plot_field(component, f=f, mode_index=mode_index, ax=ax, val=val, **imshow_kwargs) for ax in flat_axes[len(available) :]: ax.set_visible(False) @@ -222,7 +223,7 @@ def plot(self, **kwargs: Any) -> tuple[Any, Any]: def overlap( self, - other: "Result | None" = None, + other: Result | None = None, *, mode_index: int = 0, other_mode_index: int | None = None, @@ -260,7 +261,7 @@ def overlap( def overlap_matrix( self, - other: "Result | None" = None, + other: Result | None = None, *, f: int | float = 0, other_f: int | float | None = None, @@ -271,7 +272,9 @@ def overlap_matrix( other = self if other is None else other other_f = f if other_f is None else other_f - values = np.empty((self.n_complex.sizes["mode_index"], other.n_complex.sizes["mode_index"]), dtype=np.complex128) + values = np.empty( + (self.n_complex.sizes["mode_index"], other.n_complex.sizes["mode_index"]), dtype=np.complex128 + ) for left_index in range(values.shape[0]): for right_index in range(values.shape[1]): values[left_index, right_index] = self.overlap( @@ -317,7 +320,7 @@ def to_hdf5(self, path: str | Path) -> Path: return destination @classmethod - def from_hdf5(cls, path: str | Path) -> "Result": + def from_hdf5(cls, path: str | Path) -> Result: """Load a :class:`Result` saved with :meth:`to_hdf5`.""" try: @@ -329,11 +332,9 @@ def from_hdf5(cls, path: str | Path) -> "Result": n_complex = cls._read_data_array(handle["n_complex"]) n_group = cls._read_data_array(handle["n_group"]) if "n_group" in handle else None dispersion = cls._read_data_array(handle["dispersion"]) if "dispersion" in handle else None - solver_info = json.loads(handle.attrs["solver_info"]) if "solver_info" in handle.attrs else None - field_components = { - component: cls._read_data_array(group) - for component, group in handle["field_components"].items() - } + solver_info = _loads_hdf5_json_attr(handle.attrs["solver_info"]) if "solver_info" in handle.attrs else None + field_group: Any = handle["field_components"] + field_components = {component: cls._read_data_array(group) for component, group in field_group.items()} return cls( n_complex=n_complex, field_components=field_components, @@ -433,16 +434,12 @@ def _add_optional_metric( def _select_field(self, component: str, *, f: int | float, mode_index: int) -> xr.DataArray: data_array = self.field_components[component] - if isinstance(f, int): - data_array = data_array.isel(f=f) - else: - data_array = data_array.sel(f=f, method="nearest") + data_array = data_array.isel(f=f) if isinstance(f, int) else data_array.sel(f=f, method="nearest") return data_array.isel(mode_index=mode_index) def _selected_mode_fields(self, *, f: int | float, mode_index: int) -> dict[str, xr.DataArray]: return { - component: self._select_field(component, f=f, mode_index=mode_index) - for component in self.field_components + component: self._select_field(component, f=f, mode_index=mode_index) for component in self.field_components } def _selected_frequency(self, f: int | float) -> float: @@ -506,6 +503,12 @@ def _json_safe(value: Any) -> Any: return value +def _loads_hdf5_json_attr(value: Any) -> Any: + if isinstance(value, bytes): + value = value.decode("utf-8") + return json.loads(value) + + def _overlap_value( left_result: Result, left: dict[str, xr.DataArray], @@ -531,7 +534,9 @@ def _overlap_value( # uses H* and measures physical flux. Lorentz deliberately does not # conjugate either mode; it is the reciprocal product used to orthogonalize # the mode set in Rust. - missing = [component for component in (*_E_COMPONENTS, *_H_COMPONENTS) if component not in left or component not in right] + missing = [ + component for component in (*_E_COMPONENTS, *_H_COMPONENTS) if component not in left or component not in right + ] if missing: raise ValueError(f"{kind} overlap requires field component(s): {', '.join(missing)}") _validate_overlap_grid(left_result, left, right_result, right) diff --git a/python/micromode/sweep.py b/python/micromode/sweep.py index af0cb32..9609816 100644 --- a/python/micromode/sweep.py +++ b/python/micromode/sweep.py @@ -2,9 +2,9 @@ from __future__ import annotations +from collections.abc import Iterable from dataclasses import dataclass from itertools import permutations -from typing import Iterable import numpy as np diff --git a/scripts/check_dist_artifacts.py b/scripts/check_dist_artifacts.py index 4692c0d..6b879e1 100644 --- a/scripts/check_dist_artifacts.py +++ b/scripts/check_dist_artifacts.py @@ -6,7 +6,6 @@ import re from pathlib import Path - WHEEL_RE = re.compile( r"^micromode-(?P[^-]+)-(?P[^-]+)-(?P[^-]+)-(?P[^.]+(?:\.[^.]+)*)\.whl$" ) @@ -41,7 +40,8 @@ def main() -> None: python_tags: set[str] = set() for wheel in wheels: match = WHEEL_RE.match(wheel.name) - require(match is not None, f"unexpected wheel filename: {wheel.name}") + if match is None: + raise SystemExit(f"unexpected wheel filename: {wheel.name}") platform = match["platform"] python_tag = match["python"] platform_tags.update(platform.split(".")) @@ -56,9 +56,7 @@ def main() -> None: missing = sorted(set(args.require_cpython) - python_tags) require(not missing, f"missing wheels for CPython: {', '.join(missing)}") missing_platforms = sorted( - prefix - for prefix in args.require_platform - if not any(tag.startswith(prefix) for tag in platform_tags) + prefix for prefix in args.require_platform if not any(tag.startswith(prefix) for tag in platform_tags) ) require(not missing_platforms, f"missing wheel platform families: {', '.join(missing_platforms)}") print(f"release artifacts look publishable: {len(sdists)} sdist(s), {len(wheels)} wheel(s)") diff --git a/scripts/check_release_metadata.py b/scripts/check_release_metadata.py index 2cfbe96..4ccd248 100644 --- a/scripts/check_release_metadata.py +++ b/scripts/check_release_metadata.py @@ -7,7 +7,7 @@ from pathlib import Path try: - import tomllib + import tomllib # type: ignore[import-not-found] except ModuleNotFoundError: # pragma: no cover - Python 3.10 compatibility. import tomli as tomllib @@ -56,10 +56,8 @@ def main() -> None: changelog = (ROOT / "CHANGELOG.md").read_text(encoding="utf-8") require(project["version"] in changelog, "Python version is not mentioned in CHANGELOG.md") - require( - re.fullmatch(r"\d+\.\d+\.\d+(a\d+|b\d+|rc\d+)?", project["version"]), - f"Python version {project['version']!r} is not a normal PyPI release version", - ) + version_match = re.fullmatch(r"\d+\.\d+\.\d+(a\d+|b\d+|rc\d+)?", project["version"]) + require(version_match is not None, f"Python version {project['version']!r} is not a normal PyPI release version") require_tag_matches_version(project["version"]) print(f"release metadata looks ready for micromode {project['version']}") @@ -71,7 +69,8 @@ def load_toml(path: str) -> dict: def python_to_cargo_version(version: str) -> str: match = re.fullmatch(r"(\d+\.\d+\.\d+)(?:(a|b|rc)(\d+))?", version) - require(match is not None, f"Python version {version!r} is not a supported release version") + if match is None: + raise ValueError(f"Python version {version!r} is not a supported release version") base, phase, number = match.groups() if phase is None: return base diff --git a/scripts/generate_coverage_badge.py b/scripts/generate_coverage_badge.py index a55b12e..400d0cf 100644 --- a/scripts/generate_coverage_badge.py +++ b/scripts/generate_coverage_badge.py @@ -5,8 +5,8 @@ import argparse import html -from pathlib import Path import xml.etree.ElementTree as ET +from pathlib import Path def main() -> None: @@ -17,7 +17,9 @@ def main() -> None: percent = coverage_percent(args.coverage_xml) args.output_svg.parent.mkdir(parents=True, exist_ok=True) - args.output_svg.write_text(render_badge("coverage", f"{percent:.0f}%", color_for_percent(percent)), encoding="utf-8") + args.output_svg.write_text( + render_badge("coverage", f"{percent:.0f}%", color_for_percent(percent)), encoding="utf-8" + ) def coverage_percent(path: Path) -> float: @@ -53,7 +55,8 @@ def render_badge(label: str, message: str, color: str) -> str: width = label_width + message_width label_center = label_width / 2 message_center = label_width + message_width / 2 - return f""" + return f""" {label}: {message} diff --git a/scripts/smoke_dist.py b/scripts/smoke_dist.py index 0fa99fe..3a249b5 100644 --- a/scripts/smoke_dist.py +++ b/scripts/smoke_dist.py @@ -9,7 +9,6 @@ import tempfile from pathlib import Path - ROOT = Path(__file__).resolve().parents[1] diff --git a/scripts/smoke_wheel.py b/scripts/smoke_wheel.py index c7ae2d8..24beac1 100644 --- a/scripts/smoke_wheel.py +++ b/scripts/smoke_wheel.py @@ -20,8 +20,8 @@ def main() -> None: result = mm.solve_grid( eps_xx=eps, - x_edges=x_edges, - y_edges=y_edges, + x_edges=tuple(float(value) for value in x_edges), + y_edges=tuple(float(value) for value in y_edges), wavelength=1.55, num_modes=1, target_neff=2.4, diff --git a/src/derivatives.rs b/src/derivatives.rs index 70882a2..2d321fa 100644 --- a/src/derivatives.rs +++ b/src/derivatives.rs @@ -349,7 +349,7 @@ pub fn create_sfactor_f( ) -> Vec { let mut sfactor = vec![Complex64::new(1.0, 0.0); n]; for (i, value) in sfactor.iter_mut().enumerate() { - if i <= n_pml - 1 && dmin_pml { + if i < n_pml && dmin_pml { *value = s_value( dls[0], (n_pml as f64 - i as f64 - 0.5) / n_pml as f64, diff --git a/src/eigensolve.rs b/src/eigensolve.rs index d36fa74..b21dbac 100644 --- a/src/eigensolve.rs +++ b/src/eigensolve.rs @@ -151,7 +151,7 @@ fn selected_sparse_shift_invert_native_eigenpairs_impl( initial_vector, options, "native_shift_invert", - profile.as_deref_mut(), + profile, |input, output| factorization.solve_into(input, output, &mut solve_workspace), ) } @@ -171,7 +171,8 @@ where { let n = mat.rows; let krylov_dim = options.krylov_dim.min(n).max(num_modes + 2); - let checkpoint_start = krylov_dim.min(((3 * krylov_dim + 3) / 4).max((num_modes + 8).max(16))); + let checkpoint_start = + krylov_dim.min((3 * krylov_dim).div_ceil(4).max((num_modes + 8).max(16))); let checkpoint_interval = 4; let stability_tolerance = options.tolerance.sqrt(); let mut previous_checkpoint_values: Option> = None; @@ -221,7 +222,7 @@ where let exhausted = norm <= options.tolerance || col + 1 == krylov_dim; let checkpoint_due = actual_dim >= checkpoint_start - && ((actual_dim - checkpoint_start) % checkpoint_interval == 0 || exhausted); + && ((actual_dim - checkpoint_start).is_multiple_of(checkpoint_interval) || exhausted); if checkpoint_due { let candidates = candidate_eigenpairs_from_hessenberg( mat, @@ -263,7 +264,7 @@ where guess_value, backend, options.tolerance, - profile.as_deref_mut(), + profile, ) } diff --git a/src/sparse_matrix.rs b/src/sparse_matrix.rs index 33a78ef..2beb1e6 100644 --- a/src/sparse_matrix.rs +++ b/src/sparse_matrix.rs @@ -101,11 +101,11 @@ impl SparseMatrix { let mut columns = (0..self.cols) .map(|_| BTreeMap::::new()) .collect::>(); - for col in 0..self.cols { + for (col, column) in columns.iter_mut().enumerate().take(self.cols) { for (row, value) in self.column_entries(col) { - columns[col].insert(row, value); + column.insert(row, value); } - *columns[col].entry(col).or_insert(Complex64::new(0.0, 0.0)) -= shift; + *column.entry(col).or_insert(Complex64::new(0.0, 0.0)) -= shift; } Self::from_column_maps(self.rows, self.cols, columns) } @@ -132,12 +132,12 @@ impl SparseMatrix { let mut columns = (0..self.cols) .map(|_| BTreeMap::::new()) .collect::>(); - for col in 0..self.cols { + for (col, column) in columns.iter_mut().enumerate().take(self.cols) { for (row, value) in self.column_entries(col) { - *columns[col].entry(row).or_insert(Complex64::new(0.0, 0.0)) += value; + *column.entry(row).or_insert(Complex64::new(0.0, 0.0)) += value; } for (row, value) in other.column_entries(col) { - *columns[col].entry(row).or_insert(Complex64::new(0.0, 0.0)) += value; + *column.entry(row).or_insert(Complex64::new(0.0, 0.0)) += value; } } Self::from_column_maps(self.rows, self.cols, columns) diff --git a/tests/test_micromode_api.py b/tests/test_micromode_api.py index 2563449..66018c3 100644 --- a/tests/test_micromode_api.py +++ b/tests/test_micromode_api.py @@ -1,5 +1,8 @@ from __future__ import annotations +from collections.abc import Sequence +from typing import Any + import numpy as np import pytest import xarray as xr @@ -7,11 +10,25 @@ import micromode as mm -def _strip_grid(nx: int = 5, ny: int = 4) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - x_edges = np.linspace(-1.0, 1.0, nx + 1) - y_edges = np.linspace(-0.8, 0.8, ny + 1) - x_centers = (x_edges[:-1] + x_edges[1:]) / 2 - y_centers = (y_edges[:-1] + y_edges[1:]) / 2 +def _linspace_edges(start: float, stop: float, count: int) -> tuple[float, ...]: + return tuple(float(value) for value in np.linspace(start, stop, count)) + + +def _edge_centers(edges: Sequence[float]) -> np.ndarray: + edge_array = np.asarray(edges, dtype=float) + return (edge_array[:-1] + edge_array[1:]) / 2 + + +def _solver_info(data: mm.Result) -> dict[str, Any]: + assert data.solver_info is not None + return data.solver_info + + +def _strip_grid(nx: int = 5, ny: int = 4) -> tuple[np.ndarray, tuple[float, ...], tuple[float, ...]]: + x_edges = _linspace_edges(-1.0, 1.0, nx + 1) + y_edges = _linspace_edges(-0.8, 0.8, ny + 1) + x_centers = _edge_centers(x_edges) + y_centers = _edge_centers(y_edges) xx, yy = np.meshgrid(x_centers, y_centers, indexing="ij") eps = np.where((np.abs(xx) <= 0.35) & (np.abs(yy) <= 0.25), 3.4**2, 1.44**2) return eps, x_edges, y_edges @@ -35,10 +52,11 @@ def test_grid_api_solves_with_rust_sparse_backend(): assert data.n_complex.values[0, 0].real >= data.n_complex.values[0, 1].real assert data.field_components["Ex"].shape == (6, 5, 1, 1, 2) assert set(data.field_components) == {"Ex", "Ey", "Ez", "Hx", "Hy", "Hz"} - assert data.solver_info["runs"][0]["phase_convention"] == "dominant_e_real_positive" - assert data.solver_info["runs"][0]["normalization"] == "lorentz_orthogonal_unit_transverse_power" - np.testing.assert_allclose(data.solver_info["runs"][0]["power_norms"], np.ones(2), rtol=1e-10, atol=1e-10) - assert data.solver_info["runs"][0]["lorentz_orthogonality_error"] < 1e-8 + run_info = _solver_info(data)["runs"][0] + assert run_info["phase_convention"] == "dominant_e_real_positive" + assert run_info["normalization"] == "lorentz_orthogonal_unit_transverse_power" + np.testing.assert_allclose(run_info["power_norms"], np.ones(2), rtol=1e-10, atol=1e-10) + assert run_info["lorentz_orthogonality_error"] < 1e-8 assert abs(abs(data.overlap(mode_index=0, kind="power", normalize=False)) - 1.0) < 1e-10 lorentz_matrix = data.overlap_matrix(kind="lorentz").values np.testing.assert_allclose(lorentz_matrix - np.diag(np.diag(lorentz_matrix)), 0.0, atol=1e-8) @@ -77,8 +95,8 @@ def test_materials_api_matches_component_api_for_diagonal_grid(): def test_materials_api_accepts_full_tensor_grid(): - x_edges = np.linspace(-1.0, 1.0, 5) - y_edges = np.linspace(-0.8, 0.8, 4) + x_edges = _linspace_edges(-1.0, 1.0, 5) + y_edges = _linspace_edges(-0.8, 0.8, 4) shape = (len(x_edges) - 1, len(y_edges) - 1) eps_xx = np.full(shape, 2.2**2) eps_yy = np.full(shape, 2.0**2) @@ -126,27 +144,30 @@ def test_solver_specs_report_diagnostics_and_persist_to_hdf5(tmp_path): krylov_dim=18, ) - assert data.solver_info["pml"]["num_cells"] == (1, 1) - assert data.solver_info["pml"]["sigma_max"] == pytest.approx(1.6) - assert data.solver_info["boundary"]["low"] == ("pec", "pmc") - assert data.solver_info["backend"] - assert data.solver_info["runs"][0]["operator_size"] > 0 - assert data.solver_info["runs"][0]["operator_nnz"] > 0 - assert len(data.solver_info["runs"][0]["residuals"]) == 1 - assert data.solver_info["runs"][0]["power_norms"][0] == pytest.approx(1.0) - assert len(data.solver_info["runs"][0]["lorentz_norms"]) == 1 + solver_info = _solver_info(data) + run_info = solver_info["runs"][0] + assert solver_info["pml"]["num_cells"] == (1, 1) + assert solver_info["pml"]["sigma_max"] == pytest.approx(1.6) + assert solver_info["boundary"]["low"] == ("pec", "pmc") + assert solver_info["backend"] + assert run_info["operator_size"] > 0 + assert run_info["operator_nnz"] > 0 + assert len(run_info["residuals"]) == 1 + assert run_info["power_norms"][0] == pytest.approx(1.0) + assert len(run_info["lorentz_norms"]) == 1 loaded = mm.Result.from_hdf5(data.to_hdf5(tmp_path / "diagnostic_result.h5")) - assert loaded.solver_info["pml"]["num_cells"] == [1, 1] - assert loaded.solver_info["boundary"]["low"] == ["pec", "pmc"] - assert loaded.solver_info["runs"][0]["operator_size"] == data.solver_info["runs"][0]["operator_size"] + loaded_solver_info = _solver_info(loaded) + assert loaded_solver_info["pml"]["num_cells"] == [1, 1] + assert loaded_solver_info["boundary"]["low"] == ["pec", "pmc"] + assert loaded_solver_info["runs"][0]["operator_size"] == run_info["operator_size"] def test_slice_api_solves_1d_x_slice_and_matches_single_cell_grid(tmp_path): - x_edges = np.linspace(-1.0, 1.0, 9) - x_centers = (x_edges[:-1] + x_edges[1:]) / 2 + x_edges = _linspace_edges(-1.0, 1.0, 9) + x_centers = _edge_centers(x_edges) eps = np.where(np.abs(x_centers) <= 0.3, 3.4**2, 1.44**2) - y_edges = np.asarray([-0.125, 0.125]) + y_edges = (-0.125, 0.125) sliced = mm.solve_slice( eps_xx=eps, @@ -190,7 +211,7 @@ def test_slice_api_solves_1d_x_slice_and_matches_single_cell_grid(tmp_path): def test_materials_slice_api_supports_y_axis_and_tensor_components(): - y_edges = np.linspace(-0.8, 0.8, 7) + y_edges = _linspace_edges(-0.8, 0.8, 7) shape = (len(y_edges) - 1,) eps_xx = np.full(shape, 2.2**2) eps_yy = np.full(shape, 2.0**2) @@ -223,8 +244,8 @@ def test_materials_slice_api_supports_y_axis_and_tensor_components(): def test_y_normal_solve_returns_global_component_names(): - z_edges = np.linspace(-1.0, 1.0, 9) - z_centers = (z_edges[:-1] + z_edges[1:]) / 2 + z_edges = _linspace_edges(-1.0, 1.0, 9) + z_centers = _edge_centers(z_edges) eps = np.where(np.abs(z_centers) <= 0.3, 3.4**2, 1.44**2) data = mm.solve_slice( @@ -247,9 +268,60 @@ def test_y_normal_solve_returns_global_component_names(): assert np.linalg.norm(data.field_components["Ez"].values) > 0 +def test_x_normal_solve_returns_global_component_names_and_positive_power(): + y_edges = _linspace_edges(-1.0, 1.0, 7) + z_edges = _linspace_edges(-0.8, 0.8, 6) + y_centers = _edge_centers(y_edges) + z_centers = _edge_centers(z_edges) + yy, zz = np.meshgrid(y_centers, z_centers, indexing="ij") + eps = np.where((np.abs(yy) <= 0.35) & (np.abs(zz) <= 0.25), 3.4**2, 1.44**2) + + data = mm.solve_grid( + eps_xx=eps, + x_edges=y_edges, + y_edges=z_edges, + freqs=[mm.C_0 / 1.55], + num_modes=1, + target_neff=2.5, + normal_axis=0, + krylov_dim=16, + ) + + assert data.field_components["Ey"].dims[:3] == ("y", "z", "x") + assert data.field_components["Ey"].attrs["normal_dim"] == "x" + assert np.linalg.norm(data.field_components["Ey"].values) > 0 + assert np.linalg.norm(data.field_components["Ez"].values) > 0 + + power = data.overlap(mode_index=0, kind="power", normalize=False) + assert power.real > 0 + assert abs(abs(power) - 1.0) < 1e-10 + + +def test_y_normal_power_overlap_is_positive(): + z_edges = _linspace_edges(-1.0, 1.0, 9) + z_centers = _edge_centers(z_edges) + eps = np.where(np.abs(z_centers) <= 0.3, 3.4**2, 1.44**2) + + data = mm.solve_slice( + eps_xx=eps, + coord_edges=z_edges, + axis="y", + invariant_width=0.5, + freqs=[mm.C_0 / 1.55], + num_modes=1, + target_neff=2.5, + normal_axis=1, + krylov_dim=16, + ) + + power = data.overlap(mode_index=0, kind="power", normalize=False) + assert power.real > 0 + assert abs(abs(power) - 1.0) < 1e-10 + + def test_materials_subpixel_averaging_helpers(): - x_edges = np.linspace(-1.0, 1.0, 3) - y_edges = np.linspace(-0.5, 0.5, 3) + x_edges = _linspace_edges(-1.0, 1.0, 3) + y_edges = _linspace_edges(-0.5, 0.5, 3) samples = np.asarray( [ [1.0, 3.0, 5.0, 7.0], @@ -303,8 +375,8 @@ def test_angle_and_bend_use_tensorial_rust_path(): def test_full_tensor_grid_supports_angle_and_bend_transform(): - x_edges = np.linspace(-1.0, 1.0, 5) - y_edges = np.linspace(-0.8, 0.8, 4) + x_edges = _linspace_edges(-1.0, 1.0, 5) + y_edges = _linspace_edges(-0.8, 0.8, 4) shape = (len(x_edges) - 1, len(y_edges) - 1) materials = mm.Materials.from_components( eps_xx=np.full(shape, 2.2**2), @@ -328,7 +400,7 @@ def test_full_tensor_grid_supports_angle_and_bend_transform(): krylov_dim=18, ) - assert data.solver_info["runs"][0]["backend_kind"] == "tensorial_sparse" + assert _solver_info(data)["runs"][0]["backend_kind"] == "tensorial_sparse" assert np.isfinite(data.n_complex.values).all() assert np.isfinite(data.field_components["Ex"].values).all() @@ -392,6 +464,27 @@ def test_result_metrics_io_plotting_and_overlap(tmp_path): plt.close("all") +def test_result_from_hdf5_accepts_bytes_solver_info(tmp_path): + eps, x_edges, y_edges = _strip_grid(6, 5) + data = mm.solve_grid( + eps_xx=eps, + x_edges=x_edges, + y_edges=y_edges, + freqs=[mm.C_0 / 1.55], + num_modes=1, + target_neff=2.5, + krylov_dim=16, + ) + path = data.to_hdf5(tmp_path / "bytes_solver_info.h5") + + import h5py + + with h5py.File(path, "a") as handle: + handle.attrs["solver_info"] = b'{"ok": true}' + + assert mm.Result.from_hdf5(path).solver_info == {"ok": True} + + def test_result_overlap_for_synthetic_orthogonal_modes(): dims = ("x", "y", "z", "f", "mode_index") coords = { @@ -473,6 +566,8 @@ def test_spec_validation_for_core_options(): boundary=mm.BoundarySpec(low=("pmc", "pec")), ) + assert isinstance(spec.pml, mm.PmlSpec) + assert isinstance(spec.boundary, mm.BoundarySpec) assert spec.pml.num_cells == (1, 2) assert spec.pml.sigma_max == pytest.approx(1.5) assert spec.boundary.dmin_pmc == (True, False) diff --git a/tests/test_mode_solver_fixtures.py b/tests/test_mode_solver_fixtures.py index d91988d..5baf67e 100644 --- a/tests/test_mode_solver_fixtures.py +++ b/tests/test_mode_solver_fixtures.py @@ -4,7 +4,7 @@ import numpy as np import pytest - +from benchmarks.compare_mode_solver_fixtures import _compare_local_case from benchmarks.mode_solver.fixtures import ( data_path, load_data_array, @@ -14,7 +14,6 @@ sha256_file, summary_path, ) -from benchmarks.compare_mode_solver_fixtures import _compare_local_case ROOT = Path(__file__).resolve().parents[1] SMOKE_FIXTURE_ROOT = ROOT / "fixtures" / "mode_solver" / "smoke" diff --git a/uv.lock b/uv.lock index c4c81e7..e48e271 100644 --- a/uv.lock +++ b/uv.lock @@ -840,7 +840,7 @@ wheels = [ [[package]] name = "micromode" -version = "0.1.0a3" +version = "0.1.0a4" source = { editable = "." } dependencies = [ { name = "h5py" }, @@ -856,10 +856,12 @@ dev = [ { name = "maturin" }, { name = "packaging" }, { name = "pkginfo" }, + { name = "pyright" }, { name = "pytest" }, { name = "pytest-cov" }, { name = "pytest-xdist" }, - { name = "tomli", marker = "python_full_version < '3.11'" }, + { name = "ruff" }, + { name = "tomli" }, { name = "twine" }, ] @@ -871,10 +873,12 @@ requires-dist = [ { name = "numpy", specifier = ">=2.2.6,<2.5.0" }, { name = "packaging", marker = "extra == 'dev'", specifier = ">=24.2" }, { name = "pkginfo", marker = "extra == 'dev'", specifier = ">=1.12.1.2" }, + { name = "pyright", marker = "extra == 'dev'", specifier = ">=1.1.390,<2" }, { name = "pytest", marker = "extra == 'dev'", specifier = ">=8.1,<10.0.0" }, { name = "pytest-cov", marker = "extra == 'dev'", specifier = ">=5,<8" }, { name = "pytest-xdist", marker = "extra == 'dev'", specifier = ">=3.6,<4.0" }, - { name = "tomli", marker = "python_full_version < '3.11' and extra == 'dev'", specifier = ">=2,<3" }, + { name = "ruff", marker = "extra == 'dev'", specifier = ">=0.8.0,<1" }, + { name = "tomli", marker = "extra == 'dev'", specifier = ">=2,<3" }, { name = "twine", marker = "extra == 'dev'", specifier = ">=5,<7" }, { name = "xarray", specifier = ">=2023.8,<2026.5.0" }, ] @@ -913,6 +917,15 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/80/7c/19cd0671d1ba2762fb388fc149697d20d0568ccfeef833b11280a619e526/nh3-0.3.5-cp38-abi3-win_arm64.whl", hash = "sha256:8f85285700a18e9f3fc5bff41fe573fa84f81542ef13b48a89f9fecca0474d3b", size = 611069, upload-time = "2026-04-25T10:44:14.934Z" }, ] +[[package]] +name = "nodeenv" +version = "1.10.0" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/24/bf/d1bda4f6168e0b2e9e5958945e01910052158313224ada5ce1fb2e1113b8/nodeenv-1.10.0.tar.gz", hash = "sha256:996c191ad80897d076bdfba80a41994c2b47c68e224c542b48feba42ba00f8bb", size = 55611, upload-time = "2025-12-20T14:08:54.006Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/88/b2/d0896bdcdc8d28a7fc5717c305f1a861c26e18c05047949fb371034d98bd/nodeenv-1.10.0-py2.py3-none-any.whl", hash = "sha256:5bb13e3eed2923615535339b3c620e76779af4cb4c6a90deccc9e36b274d3827", size = 23438, upload-time = "2025-12-20T14:08:52.782Z" }, +] + [[package]] name = "numpy" version = "2.2.6" @@ -1268,6 +1281,19 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/10/bd/c038d7cc38edc1aa5bf91ab8068b63d4308c66c4c8bb3cbba7dfbc049f9c/pyparsing-3.3.2-py3-none-any.whl", hash = "sha256:850ba148bd908d7e2411587e247a1e4f0327839c40e2e5e6d05a007ecc69911d", size = 122781, upload-time = "2026-01-21T03:57:55.912Z" }, ] +[[package]] +name = "pyright" +version = "1.1.409" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "nodeenv" }, + { name = "typing-extensions" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/51/4e/3aa27f74211522dba7e9cbc3e74de779c6d4b654c54e50a4840623be8014/pyright-1.1.409.tar.gz", hash = "sha256:986ee05beca9e077c165758ad123667c679e050059a2546aa02473930394bc93", size = 4430434, upload-time = "2026-04-23T11:02:03.799Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/16/6b/330d8ebae582b30c2959a1ef4c3bc344ebde48c2ff0c3f113c4710735e11/pyright-1.1.409-py3-none-any.whl", hash = "sha256:aa3ea228cab90c845c7a60d28db7a844c04315356392aa09fafcee98c8c22fb3", size = 6438161, upload-time = "2026-04-23T11:02:01.309Z" }, +] + [[package]] name = "pytest" version = "9.0.3" @@ -1406,6 +1432,31 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/82/3b/64d4899d73f91ba49a8c18a8ff3f0ea8f1c1d75481760df8c68ef5235bf5/rich-15.0.0-py3-none-any.whl", hash = "sha256:33bd4ef74232fb73fe9279a257718407f169c09b78a87ad3d296f548e27de0bb", size = 310654, upload-time = "2026-04-12T08:24:02.83Z" }, ] +[[package]] +name = "ruff" +version = "0.15.13" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/24/21/a7d5c126d5b557715ef81098f3db2fe20f622a039ff2e626af28d674ab80/ruff-0.15.13.tar.gz", hash = "sha256:f9d89f17f7ba7fb2ed42921f0df75da797a9a5d71bc39049e2c687cf2baf44b7", size = 4678180, upload-time = "2026-05-14T13:44:37.869Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/c6/61/11d458dc6ac22504fd8e237b29dfd40504c7fbbcc8930402cfe51a8e63ed/ruff-0.15.13-py3-none-linux_armv6l.whl", hash = "sha256:444b580fc72fd6887e650acd3e575e18cdc79dbcf42fb4030b491057921f61f8", size = 10738279, upload-time = "2026-05-14T13:44:18.7Z" }, + { url = "https://files.pythonhosted.org/packages/86/ca/caa871ee7be718c45256fada4e16a218ee3e33f0c4a46b729a60a24912e6/ruff-0.15.13-py3-none-macosx_10_12_x86_64.whl", hash = "sha256:6590d009e7cb7ebf36f83dbdd44a3fa48a0994ff6f1cdc1b08006abe58f98dc7", size = 11124798, upload-time = "2026-05-14T13:44:06.427Z" }, + { url = "https://files.pythonhosted.org/packages/d3/19/43f5f2e568dddde567fc41f8471f9432c09563e19d3e617a48cfa52f8f0a/ruff-0.15.13-py3-none-macosx_11_0_arm64.whl", hash = "sha256:1c26d2f66163deeb6e08d8b39fbbe983ce3c71cea06a6d7591cfd1421793c629", size = 10460761, upload-time = "2026-05-14T13:44:04.375Z" }, + { url = "https://files.pythonhosted.org/packages/99/df/cf938cd6de3003178f03ad7c1ea2a6c099468c03a35037985070b37e76be/ruff-0.15.13-py3-none-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:9dbd6f94b434f896308e4d57fb7bfde0d02b99f7a64b3bdab0fdfa6a864203a5", size = 10804451, upload-time = "2026-05-14T13:44:25.221Z" }, + { url = "https://files.pythonhosted.org/packages/c7/7d/5d0973129b154ded2225729169d7068f26b467760b146493fde138415f23/ruff-0.15.13-py3-none-manylinux_2_17_armv7l.manylinux2014_armv7l.whl", hash = "sha256:bf3259f3be4d181bda591da5db2571aed6853c6a048157756448020bc6c5cd22", size = 10534285, upload-time = "2026-05-14T13:44:08.888Z" }, + { url = "https://files.pythonhosted.org/packages/1f/e3/6b999bbc66cd51e5f073842bc2a3995e99c5e0e72e16b15e7261f7abf57a/ruff-0.15.13-py3-none-manylinux_2_17_i686.manylinux2014_i686.whl", hash = "sha256:ae9c17e5eb4430c154e76abc25d79a318190f5a997f38fb6b114416c5319ffc9", size = 11312063, upload-time = "2026-05-14T13:44:11.274Z" }, + { url = "https://files.pythonhosted.org/packages/af/5a/642639e9f5db04f1e97fbd6e091c6fd20725bdf072fb114d00eefb9e6eb8/ruff-0.15.13-py3-none-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:2e2e39bff6c341f4b577a21b801326fab0b11847f48fcaa83f00a113c9b3cb55", size = 12183079, upload-time = "2026-05-14T13:44:01.634Z" }, + { url = "https://files.pythonhosted.org/packages/19/4c/7585735f6b53b0f12de13618b2f7d250a844f018822efc899df2e7b8295f/ruff-0.15.13-py3-none-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:e8d9a8e08013542e94d3220bc5b62cc3e5ef87c5f74bff367d3fac14fab013e6", size = 11440833, upload-time = "2026-05-14T13:43:59.043Z" }, + { url = "https://files.pythonhosted.org/packages/e8/31/bf1a0803d077e679cfeee5f2f67290a0fa79c7385b5d9a8c17b9db2c48f0/ruff-0.15.13-py3-none-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:cc411dfebe5eebe55ce041c6ae080eb7668955e866daa2fbb16692a784f1c4ca", size = 11434486, upload-time = "2026-05-14T13:44:27.761Z" }, + { url = "https://files.pythonhosted.org/packages/e1/4e/62c9b999875d4f14db80f277c030578f5e249c9852d65b7ac7ad0b43c041/ruff-0.15.13-py3-none-manylinux_2_31_riscv64.whl", hash = "sha256:768494eb08b9cee54e2fd27969966f74db5a57f6eaa7a90fcb3306af34dfc4bd", size = 11385189, upload-time = "2026-05-14T13:44:13.704Z" }, + { url = "https://files.pythonhosted.org/packages/fc/89/7e959047a104df3eb12863447c110140191fc5b6c4f379ea2e803fcdb0e4/ruff-0.15.13-py3-none-musllinux_1_2_aarch64.whl", hash = "sha256:fb75f9a3a7e42ffe117d734494e6c5e5cb3565d66e12612cb63d0e572a41a5b6", size = 10781380, upload-time = "2026-05-14T13:43:56.734Z" }, + { url = "https://files.pythonhosted.org/packages/ff/52/5fd18f3b88cab63e88aa11516b3b4e1e5f720e5c330f8dbe5c26210f41f8/ruff-0.15.13-py3-none-musllinux_1_2_armv7l.whl", hash = "sha256:8cb74dd33bb2f6613faf7fc03b660053b5ac4f80e706d5788c6335e2a8048d51", size = 10540605, upload-time = "2026-05-14T13:44:20.748Z" }, + { url = "https://files.pythonhosted.org/packages/e8/e0/9e35f338990d3e41a82875ff7053ffe97541dae81c9d02143177f381d572/ruff-0.15.13-py3-none-musllinux_1_2_i686.whl", hash = "sha256:7ef823f817fcd191dc934e984be9cf4094f808effa16f2542ad8e821ba02bbf2", size = 11036554, upload-time = "2026-05-14T13:44:16.256Z" }, + { url = "https://files.pythonhosted.org/packages/c2/13/070fb048c24080fba188f66371e2a92785be257ad02242066dc7255ac6e9/ruff-0.15.13-py3-none-musllinux_1_2_x86_64.whl", hash = "sha256:f345a13937bd7f09f6f5d19fa0721b0c103e00e7f62bc67089a8e5e037719e0b", size = 11528133, upload-time = "2026-05-14T13:44:22.808Z" }, + { url = "https://files.pythonhosted.org/packages/6b/8c/b1e1666aef7fc6555094d73ae6cd981701781ae85b97ceefc0eebd0b4668/ruff-0.15.13-py3-none-win32.whl", hash = "sha256:4044f94208b3b05ba0fc4a4abd0558cf4d6459bd18325eead7fd8cc66f909b41", size = 10721455, upload-time = "2026-05-14T13:44:35.697Z" }, + { url = "https://files.pythonhosted.org/packages/ab/a6/870a3e8a50590bb92be184ad928c2922f088b00d9dc5c5ec7b924ee08c22/ruff-0.15.13-py3-none-win_amd64.whl", hash = "sha256:7064884d442b7d477b4e7473d12da7f08851d2b1982763c5d3f388a19468a1a4", size = 11900409, upload-time = "2026-05-14T13:44:30.389Z" }, + { url = "https://files.pythonhosted.org/packages/9b/36/9c015cd052fca743dae8cb2aeb16b551444787467db42ceab0fc968865af/ruff-0.15.13-py3-none-win_arm64.whl", hash = "sha256:2471da9bd1068c8c064b5fd9c0c4b6dddffd6369cb1cd68b29993b1709ff1b21", size = 11179336, upload-time = "2026-05-14T13:44:33.026Z" }, +] + [[package]] name = "secretstorage" version = "3.5.0"