From 69509c85477ea12b73ddb9664c089933fa872e1b Mon Sep 17 00:00:00 2001 From: mdmaas Date: Wed, 19 Aug 2026 21:02:54 -0300 Subject: [PATCH 1/4] incorporate palace binary downloader --- README.md | 8 +- pyproject.toml | 2 + src/gsim/palace/__init__.py | 9 +- src/gsim/palace/base.py | 20 ++- src/gsim/palace/runtime.py | 289 ++++++++++++++++++++++++++++--- tests/palace/test_sim_classes.py | 242 ++++++++++++++++++++++---- 6 files changed, 504 insertions(+), 66 deletions(-) diff --git a/README.md b/README.md index ed4edd5c..6bbc1697 100644 --- a/README.md +++ b/README.md @@ -32,11 +32,9 @@ with minimal boilerplate. pip install gsim ``` -For local Palace execution on Linux x86_64, install the optional prebuilt binary separately: - -```bash -pip install "palacetoolkit-palace-cpu @ https://github.com/EpsilonForge/PalaceToolkit/releases/download/palace-cpu-v0.1.2/palacetoolkit_palace_cpu-0.1.0-py3-none-linux_x86_64.whl" -``` +For local Palace execution on Linux x86_64, gsim can auto-download and cache a prebuilt Palace CPU binary on first use. +No additional package installation is required. If you prefer to supply a binary yourself, set `PALACE_BIN` (path to an +executable) or put `palace` on your `PATH`; a Palace SIF image may be given via `PALACE_SIF` for Apptainer-based runs. For development (requires [uv](https://docs.astral.sh/uv/)): diff --git a/pyproject.toml b/pyproject.toml index 9d18d06c..5e16c5a9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -275,6 +275,8 @@ select = ["ALL"] ] "tests/**/*.py" = [ "ANN", # flake8-annotations + "ARG001", # unused-function-argument + "ARG005", # unused-lambda-argument (test stubs) "D", # pydocstyle "INP001", # implicit-namespace-package "PLC0415", # allow imports inside tests diff --git a/src/gsim/palace/__init__.py b/src/gsim/palace/__init__.py index 836c67cd..b155a643 100644 --- a/src/gsim/palace/__init__.py +++ b/src/gsim/palace/__init__.py @@ -136,8 +136,12 @@ load_sparams, ) -# Runtime / binary resolution (optional palace-toolkit-cpu dependency) -from gsim.palace.runtime import resolve_palace_binary, resolve_palace_library_dir +# Runtime / binary resolution (self-contained; can auto-download a Palace CPU runtime) +from gsim.palace.runtime import ( + install_palace_runtime, + resolve_palace_binary, + resolve_palace_library_dir, +) from gsim.viz import ( close_interactive_view, close_interactive_views, @@ -206,6 +210,7 @@ "get_material_properties", "get_port_map", "get_stack", + "install_palace_runtime", "interactive_mode", "load_boundary_field_data", "load_field_context", diff --git a/src/gsim/palace/base.py b/src/gsim/palace/base.py index 0da686a0..b7151138 100644 --- a/src/gsim/palace/base.py +++ b/src/gsim/palace/base.py @@ -2115,13 +2115,23 @@ def run_local( resolved_exe = bundled lib_dir = resolve_palace_library_dir() if verbose: + from gsim.palace.runtime import ( + _cached_binary as _cached, + ) from gsim.palace.runtime import ( _palace_cpu_available as _cpu_avail, ) + from gsim.palace.runtime import ( + _palace_toolkit_available as _toolkit_avail, + ) source = ( "palace-toolkit-cpu" if _cpu_avail() + else "gsim cached runtime" + if _cached() is not None + else "palace-toolkit" + if _toolkit_avail() else "PALACE_BIN / PATH" ) logger.info( @@ -2149,9 +2159,9 @@ def run_local( if resolved_exe is None: raise FileNotFoundError( - "Palace executable not found. Set PALACE_BIN, " - "PALACE_EXECUTABLE, or install the optional " - "palacetoolkit-palace-cpu wheel documented in the gsim README." + "Palace executable not found. Set PALACE_BIN or " + "PALACE_EXECUTABLE, install a Palace binary, or let gsim " + "auto-download the prebuilt CPU runtime (Linux x86_64)." ) exe_path = Path(resolved_exe) @@ -2163,8 +2173,8 @@ def run_local( raise FileNotFoundError( f"Palace executable not found: {exe_path}. " "Install Palace directly or provide correct path via " - "palace_executable, or install the optional " - "palacetoolkit-palace-cpu wheel documented in the gsim README." + "palace_executable, or let gsim auto-download the " + "prebuilt CPU runtime (Linux x86_64)." ) exe_path = Path(resolved) diff --git a/src/gsim/palace/runtime.py b/src/gsim/palace/runtime.py index c3f51ef4..47195eab 100644 --- a/src/gsim/palace/runtime.py +++ b/src/gsim/palace/runtime.py @@ -1,38 +1,218 @@ """Palace runtime/binary resolution. -Provides a unified resolver for locating a Palace executable, with -optional delegation to the ``palacetoolkit_palace_cpu`` package (the -``palace-toolkit-cpu`` distribution) when installed. +Provides a unified resolver for locating a Palace executable, plus a +self-contained installer that downloads and caches a prebuilt Palace CPU +binary. gsim absorbs this functionality so users do **not** need to install +any direct-URL wheel or third-party runtime package to run Palace locally on +Linux x86_64. Resolution order ----------------- 1. ``PALACE_BIN`` environment variable. 2. ``PALACE_EXECUTABLE`` environment variable, or ``"palace"`` in ``PATH``. -3. ``palacetoolkit_palace_cpu`` packaged binary (when the optional - ``palace-toolkit-cpu`` extra is installed). -4. ``None`` if nothing was found. +3. ``palacetoolkit_palace_cpu`` packaged binary (if the legacy + ``palace-toolkit-cpu`` wheel happens to be installed). +4. gsim's own cached/downloaded Palace CPU runtime (Linux x86_64). +5. Delegation to ``palacetoolkit`` (if the ``palace-toolkit`` distribution + happens to be installed). +6. ``None`` if nothing was found. + +The auto-download is only attempted when ``PALACETOOLKIT_AUTO_DOWNLOAD_BINARY`` +is not disabled, and only on Linux x86_64 (the platform the prebuilt Palace CPU +wheel is provided for). """ from __future__ import annotations import importlib.util +import json import logging import os +import platform import shutil +import stat import subprocess +import tempfile +from contextlib import suppress from pathlib import Path +from urllib.request import Request, urlopen +from zipfile import ZipFile logger = logging.getLogger(__name__) +_DEFAULT_BINARY_TAG = "0.17.0" +_AUTO_DOWNLOAD_ENV = "PALACETOOLKIT_AUTO_DOWNLOAD_BINARY" +_TAG_ENV = "PALACETOOLKIT_PALACE_CPU_TAG" +_CACHE_ENV = "PALACETOOLKIT_RUNTIME_DIR" + + +def _is_linux_x86_64() -> bool: + """Return whether the current platform is Linux on x86_64. + + The prebuilt Palace CPU runtime is only provided for this platform. + """ + return platform.system() == "Linux" and platform.machine() == "x86_64" + + +def _runtime_cache_dir() -> Path: + """Return the directory used to cache downloaded Palace runtimes.""" + root = os.environ.get(_CACHE_ENV, "").strip() + if root: + return Path(root).expanduser().resolve() + return (Path.home() / ".cache" / "palacetoolkit" / "runtime").resolve() + + +def _binary_tag() -> str: + """Return the Palace CPU runtime version tag to download.""" + return os.environ.get(_TAG_ENV, _DEFAULT_BINARY_TAG).strip() or _DEFAULT_BINARY_TAG + + +def _binary_wheel_url(tag: str) -> str: + """Return the GitHub release URL for the given Palace CPU runtime tag.""" + return ( + "https://github.com/EpsilonForge/PalaceToolkit/releases/download/" + f"palace-cpu-v{tag}/" + f"palacetoolkit_palace_cpu-{tag}-py3-none-linux_x86_64.whl" + ) + + +def _binary_wheel_url_from_release(tag: str, timeout: float) -> str | None: + """Discover the current wheel URL from the GitHub release API (best-effort).""" + api_url = ( + "https://api.github.com/repos/EpsilonForge/PalaceToolkit/releases/tags/" + f"palace-cpu-v{tag}" + ) + request = Request( # noqa: S310 + api_url, headers={"Accept": "application/vnd.github+json"} + ) + with urlopen(request, timeout=timeout) as response: # noqa: S310 + payload = json.loads(response.read().decode("utf-8")) + + for asset in payload.get("assets", []): + name = str(asset.get("name", "")) + if name.endswith("linux_x86_64.whl") and "palacetoolkit_palace_cpu-" in name: + url = str(asset.get("browser_download_url", "")) + if url: + return url + return None + + +def _cached_runtime_prefix(tag: str | None = None) -> Path: + """Return the cache directory for a specific runtime tag.""" + resolved_tag = tag or _binary_tag() + return _runtime_cache_dir() / f"palace-cpu-v{resolved_tag}" + + +def _set_executable(path: Path) -> None: + """Make the given path executable for all users.""" + mode = path.stat().st_mode + path.chmod(mode | stat.S_IXUSR | stat.S_IXGRP | stat.S_IXOTH) + + +def install_palace_runtime(force: bool = False, timeout: float = 180.0) -> Path: + """Download and cache the prebuilt Palace CPU runtime. + + Returns: + Path to the cached ``palace`` launcher executable. + + Raises: + RuntimeError: If the platform is unsupported or the download/install + fails. + """ + if not _is_linux_x86_64(): + raise RuntimeError( + "Prebuilt runtime download is only supported on Linux x86_64" + ) + + tag = _binary_tag() + prefix = _cached_runtime_prefix(tag) + bin_palace = prefix / "bin" / "palace" + lib_dir = prefix / "lib" + if not force and bin_palace.is_file() and lib_dir.is_dir(): + return bin_palace + + prefix.mkdir(parents=True, exist_ok=True) + downloads = _runtime_cache_dir() / "downloads" + downloads.mkdir(parents=True, exist_ok=True) + wheel_name = f"palacetoolkit_palace_cpu-{tag}-py3-none-linux_x86_64.whl" + wheel_path = downloads / wheel_name + + if force or not wheel_path.is_file(): + url = _binary_wheel_url(tag) + with suppress(Exception): + discovered = _binary_wheel_url_from_release(tag, timeout=timeout) + if discovered: + url = discovered + with urlopen(url, timeout=timeout) as response: # noqa: S310 + wheel_path.write_bytes(response.read()) + + with tempfile.TemporaryDirectory( + prefix="palace-runtime-", dir=_runtime_cache_dir() + ) as tmp: + tmp_path = Path(tmp) + with ZipFile(wheel_path, "r") as wheel_zip: + wheel_zip.extractall(tmp_path) + + payload_root = tmp_path / "palacetoolkit_palace_cpu" + if not payload_root.is_dir(): + raise RuntimeError( + "Downloaded wheel does not contain palacetoolkit_palace_cpu payload" + ) + + bin_src = payload_root / "bin" + lib_src = payload_root / "lib" + if not bin_src.is_dir() or not lib_src.is_dir(): + raise RuntimeError( + "Downloaded wheel is missing expected bin/lib runtime directories" + ) + + if prefix.exists(): + shutil.rmtree(prefix) + prefix.mkdir(parents=True, exist_ok=True) + shutil.copytree(bin_src, prefix / "bin") + shutil.copytree(lib_src, prefix / "lib") + + if not bin_palace.is_file(): + raise RuntimeError("Cached runtime install did not produce bin/palace") + _set_executable(bin_palace) + bin_native = prefix / "bin" / "palace-x86_64.bin" + if bin_native.is_file(): + _set_executable(bin_native) + return bin_palace + + +def _cached_binary() -> Path | None: + """Return the cached ``palace`` launcher path, or ``None`` if not present.""" + candidate = _cached_runtime_prefix() / "bin" / "palace" + return candidate if candidate.is_file() else None + + +def _cached_library_dir() -> Path | None: + """Return the cached runtime ``lib`` directory, or ``None`` if not present.""" + candidate = _cached_runtime_prefix() / "lib" + return candidate if candidate.is_dir() else None + + +def _auto_download_enabled() -> bool: + """Return whether auto-download of the Palace runtime is enabled.""" + raw = os.environ.get(_AUTO_DOWNLOAD_ENV, "1").strip().lower() + return raw not in {"0", "false", "no", "off"} + def _palace_cpu_available() -> bool: - """Check whether the optional ``palacetoolkit_palace_cpu`` package is installed.""" + """Check whether the legacy ``palacetoolkit_palace_cpu`` package is installed.""" return importlib.util.find_spec("palacetoolkit_palace_cpu") is not None +def _palace_toolkit_available() -> bool: + """Check whether the ``palacetoolkit`` package (``palace-toolkit``) is installed.""" + return importlib.util.find_spec("palacetoolkit") is not None + + def resolve_palace_binary( *, prefer_bundled: bool = False, + download_if_missing: bool = True, ) -> Path | None: """Return a path to a runnable Palace executable, or ``None``. @@ -40,9 +220,10 @@ def resolve_palace_binary( ---------- prefer_bundled: If ``True``, skip the ``PALACE_BIN`` / ``PALACE_EXECUTABLE`` / - ``PATH`` checks and go straight to the palace-toolkit-cpu bundled - binary (useful when the caller explicitly wants the bundled - runtime). + ``PATH`` checks and go straight to gsim's cached/bundled runtime. + download_if_missing: + If ``True`` (default), auto-download and cache a prebuilt Palace CPU + runtime on Linux x86_64 when no binary is found elsewhere. Returns: ------- @@ -74,7 +255,7 @@ def resolve_palace_binary( ) return Path(resolved).resolve() - # 3. Optional palace-toolkit-cpu bundled binary + # 3. Legacy palace-toolkit-cpu packaged binary if _palace_cpu_available(): from palacetoolkit_palace_cpu import palace_binary_path @@ -94,34 +275,101 @@ def resolve_palace_binary( "resolve_palace_binary: palace-toolkit-cpu not installed — skipping" ) + # 4. gsim's own cached runtime + cached = _cached_binary() + if cached is not None and _binary_is_runnable(cached, _cached_library_dir()): + logger.info( + "resolve_palace_binary: using gsim cached runtime %s", + cached, + ) + return cached.resolve() + + # 5. Auto-download a prebuilt Palace CPU runtime (Linux x86_64) + if download_if_missing and _is_linux_x86_64() and _auto_download_enabled(): + with suppress(Exception): + downloaded = install_palace_runtime(force=False) + if _binary_is_runnable(downloaded, _cached_library_dir()): + logger.info( + "resolve_palace_binary: using gsim downloaded runtime %s", + downloaded, + ) + return downloaded.resolve() + + # 6. Delegation to the palace-toolkit package (if installed) as a fallback + if _palace_toolkit_available(): + try: + from palacetoolkit.palace_runtime import ( + resolve_palace_binary as _ptk_resolve_binary, + ) + + candidate = _ptk_resolve_binary() + except Exception as exc: + logger.debug( + "resolve_palace_binary: palacetoolkit resolver failed: %s", exc + ) + candidate = None + if candidate is not None: + candidate = Path(candidate) + if candidate.is_file() and _binary_is_runnable(candidate): + logger.info( + "resolve_palace_binary: using palace-toolkit runtime %s", + candidate, + ) + return candidate.resolve() + return None def resolve_palace_library_dir() -> Path | None: """Return the Palace library directory (for ``LD_LIBRARY_PATH``). - Only available when ``palace-toolkit-cpu`` is installed and provides a - bundled ``lib/`` directory alongside its binary. - Returns: ------- Path | None """ - if not _palace_cpu_available(): - return None + if _palace_cpu_available(): + from palacetoolkit_palace_cpu import palace_library_path - from palacetoolkit_palace_cpu import palace_library_path + lib_dir = palace_library_path() + if lib_dir.is_dir(): + return lib_dir.resolve() - lib_dir = palace_library_path() - return lib_dir.resolve() if lib_dir.is_dir() else None + cached = _cached_library_dir() + if cached is not None: + return cached.resolve() + if _palace_toolkit_available(): + try: + from palacetoolkit.palace_runtime import ( + resolve_palace_library_dir as _ptk_resolve_lib, + ) + + lib_dir = _ptk_resolve_lib() + except Exception as exc: + logger.debug( + "resolve_palace_library_dir: palacetoolkit resolver failed: %s", + exc, + ) + return None + if lib_dir is not None and lib_dir.is_dir(): + return lib_dir.resolve() -def _binary_is_runnable(binary: Path, timeout: float = 15.0) -> bool: + return None + + +def _binary_is_runnable( + binary: Path, lib_dir: Path | None = None, timeout: float = 15.0 +) -> bool: """Smoke test: file exists, is executable, and responds to --version or --help.""" bin_str = str(binary) if not binary.is_file() or not os.access(bin_str, os.X_OK): return False + run_env = os.environ.copy() + if lib_dir is not None and lib_dir.is_dir(): + prior = run_env.get("LD_LIBRARY_PATH", "") + run_env["LD_LIBRARY_PATH"] = f"{lib_dir}:{prior}" if prior else str(lib_dir) + for flag in ("--version", "--help"): try: result = subprocess.run( # noqa: S603 @@ -130,6 +378,7 @@ def _binary_is_runnable(binary: Path, timeout: float = 15.0) -> bool: text=True, timeout=timeout, check=False, + env=run_env, ) if result.returncode == 0: return True diff --git a/tests/palace/test_sim_classes.py b/tests/palace/test_sim_classes.py index bddc6b05..6961c292 100644 --- a/tests/palace/test_sim_classes.py +++ b/tests/palace/test_sim_classes.py @@ -6,13 +6,10 @@ from __future__ import annotations +import os import sys from pathlib import Path from types import SimpleNamespace -from typing import TYPE_CHECKING - -if TYPE_CHECKING: - from collections.abc import Generator import pytest @@ -583,35 +580,34 @@ def _mock_gcloud(monkeypatch: pytest.MonkeyPatch) -> None: "print_job_summary", "run_simulation", ): - setattr(gcloud, name, lambda *a, **kw: None) # noqa: ARG005 + setattr(gcloud, name, lambda *a, **kw: None) gcloud.RunResult = type("RunResult", (), {}) # ty: ignore[unresolved-attribute] monkeypatch.setitem(sys.modules, "gsim.gcloud", gcloud) @pytest.fixture -def _no_palacetoolkit(monkeypatch: pytest.MonkeyPatch) -> Generator[None]: - """Ensure palacetoolkit_palace_cpu is unavailable during the test.""" - if "palacetoolkit_palace_cpu" in sys.modules: - old = sys.modules["palacetoolkit_palace_cpu"] - monkeypatch.delitem(sys.modules, "palacetoolkit_palace_cpu", raising=False) - yield - sys.modules["palacetoolkit_palace_cpu"] = old - else: - yield +def _no_local_runtime(monkeypatch: pytest.MonkeyPatch) -> None: + """Isolate gsim's own cached/downloaded runtime during resolver tests.""" + import gsim.palace.runtime as rt + + monkeypatch.setattr(rt, "_cached_binary", lambda: None) + monkeypatch.setattr(rt, "_cached_library_dir", lambda: None) + monkeypatch.setattr(rt, "_is_linux_x86_64", lambda: False) + monkeypatch.setattr(rt, "_auto_download_enabled", lambda: False) + monkeypatch.setattr(rt, "_palace_cpu_available", lambda: False) + monkeypatch.setattr(rt, "_palace_toolkit_available", lambda: False) class TestResolvePalaceBinary: - @pytest.mark.usefixtures("_mock_gcloud", "_no_palacetoolkit") + @pytest.mark.usefixtures("_mock_gcloud", "_no_local_runtime") def test_returns_none_when_nothing_found(self) -> None: from gsim.palace.runtime import resolve_palace_binary with pytest.MonkeyPatch().context() as mp: mp.delenv("PALACE_BIN", raising=False) mp.delenv("PALACE_EXECUTABLE", raising=False) - with mp.context() as mp2: - mp2.setattr("gsim.palace.runtime._palace_cpu_available", lambda: False) - result = resolve_palace_binary() - assert result is None + result = resolve_palace_binary() + assert result is None @pytest.mark.usefixtures("_mock_gcloud") def test_uses_palace_bin_env(self) -> None: @@ -621,7 +617,7 @@ def test_uses_palace_bin_env(self) -> None: with pytest.MonkeyPatch().context() as mp: mp.setenv("PALACE_BIN", str(fake_bin)) - mp.setattr("gsim.palace.runtime._binary_is_runnable", lambda _: True) + mp.setattr("gsim.palace.runtime._binary_is_runnable", lambda *a, **k: True) mp.setattr("pathlib.Path.is_file", lambda _: True) result = resolve_palace_binary() assert result is not None @@ -640,11 +636,76 @@ def palace_binary_path() -> Path: with pytest.MonkeyPatch().context() as mp: mp.setitem(sys.modules, "palacetoolkit_palace_cpu", _FakePalaceCPU()) mp.setattr("gsim.palace.runtime._palace_cpu_available", lambda: True) - mp.setattr("gsim.palace.runtime._binary_is_runnable", lambda _: True) + mp.setattr("gsim.palace.runtime._palace_toolkit_available", lambda: False) + mp.setattr("gsim.palace.runtime._cached_binary", lambda: None) + mp.setattr("gsim.palace.runtime._binary_is_runnable", lambda *a, **k: True) mp.setattr("pathlib.Path.is_file", lambda _: True) result = resolve_palace_binary() assert result is not None + @pytest.mark.usefixtures("_mock_gcloud") + def test_delegates_to_palacetoolkit_package(self) -> None: + from gsim.palace.runtime import resolve_palace_binary + + fake_ptk_bin = Path("/opt/palacetoolkit/runtime/bin/palace") + + import types + + ptk = types.ModuleType("palacetoolkit") + ptk.__path__ = [] # type: ignore[attr-defined] + ptk_runtime = types.ModuleType("palacetoolkit.palace_runtime") + setattr(ptk_runtime, "resolve_palace_binary", lambda: fake_ptk_bin) # noqa: B010 + + with pytest.MonkeyPatch().context() as mp: + mp.setattr("gsim.palace.runtime._palace_cpu_available", lambda: False) + mp.setattr("gsim.palace.runtime._palace_toolkit_available", lambda: True) + mp.setattr("gsim.palace.runtime._cached_binary", lambda: None) + mp.setattr("gsim.palace.runtime._binary_is_runnable", lambda *a, **k: True) + mp.setattr("pathlib.Path.is_file", lambda _: True) + mp.setitem(sys.modules, "palacetoolkit", ptk) + mp.setitem(sys.modules, "palacetoolkit.palace_runtime", ptk_runtime) + result = resolve_palace_binary() + assert result is not None + + @pytest.mark.usefixtures("_mock_gcloud") + def test_uses_gsim_cached_runtime(self) -> None: + from gsim.palace.runtime import resolve_palace_binary + + fake_bin = Path("/home/user/.cache/palacetoolkit/runtime/bin/palace") + + with pytest.MonkeyPatch().context() as mp: + mp.delenv("PALACE_BIN", raising=False) + mp.delenv("PALACE_EXECUTABLE", raising=False) + mp.setattr("gsim.palace.runtime._palace_cpu_available", lambda: False) + mp.setattr("gsim.palace.runtime._palace_toolkit_available", lambda: False) + mp.setattr("gsim.palace.runtime._cached_binary", lambda: fake_bin) + mp.setattr("gsim.palace.runtime._binary_is_runnable", lambda *a, **k: True) + result = resolve_palace_binary() + assert result is not None + + @pytest.mark.usefixtures("_mock_gcloud") + def test_downloads_runtime_when_missing(self) -> None: + from gsim.palace.runtime import resolve_palace_binary + + fake_downloaded = Path("/home/user/.cache/palacetoolkit/runtime/bin/palace") + + with pytest.MonkeyPatch().context() as mp: + mp.delenv("PALACE_BIN", raising=False) + mp.delenv("PALACE_EXECUTABLE", raising=False) + mp.setattr("gsim.palace.runtime._palace_cpu_available", lambda: False) + mp.setattr("gsim.palace.runtime._palace_toolkit_available", lambda: False) + mp.setattr("gsim.palace.runtime._cached_binary", lambda: None) + mp.setattr("gsim.palace.runtime._is_linux_x86_64", lambda: True) + mp.setattr("gsim.palace.runtime._auto_download_enabled", lambda: True) + mp.setattr( + "gsim.palace.runtime.install_palace_runtime", + lambda **k: fake_downloaded, + ) + mp.setattr("gsim.palace.runtime._cached_library_dir", lambda: None) + mp.setattr("gsim.palace.runtime._binary_is_runnable", lambda *a, **k: True) + result = resolve_palace_binary() + assert result == fake_downloaded.resolve() + @pytest.mark.usefixtures("_mock_gcloud") def test_prefer_bundled_skips_env(self) -> None: from gsim.palace.runtime import resolve_palace_binary @@ -652,18 +713,93 @@ def test_prefer_bundled_skips_env(self) -> None: with pytest.MonkeyPatch().context() as mp: mp.setenv("PALACE_BIN", "/usr/bin/palace") mp.setattr("gsim.palace.runtime._palace_cpu_available", lambda: False) + mp.setattr("gsim.palace.runtime._palace_toolkit_available", lambda: False) + mp.setattr("gsim.palace.runtime._cached_binary", lambda: None) + mp.setattr("gsim.palace.runtime._auto_download_enabled", lambda: False) result = resolve_palace_binary(prefer_bundled=True) assert result is None +class TestInstallPalaceRuntime: + @pytest.mark.usefixtures("_mock_gcloud") + def test_returns_cached_binary_when_present(self, tmp_path: Path) -> None: + import gsim.palace.runtime as rt + + tag = "0.17.0" + with pytest.MonkeyPatch().context() as mp: + mp.setattr(rt, "_runtime_cache_dir", lambda: tmp_path) + mp.setattr(rt, "_binary_tag", lambda: tag) + prefix = tmp_path / f"palace-cpu-v{tag}" + (prefix / "bin").mkdir(parents=True) + (prefix / "lib").mkdir(parents=True) + bin_palace = prefix / "bin" / "palace" + bin_palace.write_text("#!/bin/sh\nexit 0\n") + bin_palace.chmod(0o755) + result = rt.install_palace_runtime(force=False) + assert result == bin_palace + + @pytest.mark.usefixtures("_mock_gcloud") + def test_raises_on_non_linux_x86_64(self) -> None: + import gsim.palace.runtime as rt + + with pytest.MonkeyPatch().context() as mp: + mp.setattr(rt, "_is_linux_x86_64", lambda: False) + with pytest.raises(RuntimeError): + rt.install_palace_runtime() + + @pytest.mark.usefixtures("_mock_gcloud") + def test_downloads_and_extracts_runtime(self, tmp_path: Path) -> None: + import io + import zipfile + + import gsim.palace.runtime as rt + + tag = "0.9.9" + cache_dir = tmp_path / "cache" + + # Build a fake wheel in memory: payload with bin/palace and lib/libfoo.so + wheel_buf = io.BytesIO() + with zipfile.ZipFile(wheel_buf, "w") as zf: + zf.writestr("palacetoolkit_palace_cpu/bin/palace", "#!/bin/sh\nexit 0\n") + zf.writestr("palacetoolkit_palace_cpu/bin/palace-x86_64.bin", "x") + zf.writestr("palacetoolkit_palace_cpu/lib/libfoo.so", "libdata") + wheel_buf.seek(0) + + class _FakeResponse: + def __enter__(self): + return self + + def __exit__(self, *exc): + return False + + def read(self): + return wheel_buf.getvalue() + + with pytest.MonkeyPatch().context() as mp: + mp.setattr(rt, "_runtime_cache_dir", lambda: cache_dir) + mp.setattr(rt, "_binary_tag", lambda: tag) + mp.setattr(rt, "_is_linux_x86_64", lambda: True) + mp.setattr( + rt, "_binary_wheel_url", lambda t: "https://example.invalid/x.whl" + ) + mp.setattr(rt, "_binary_wheel_url_from_release", lambda t, timeout: None) + mp.setattr(rt, "urlopen", lambda *a, **k: _FakeResponse()) + + result = rt.install_palace_runtime(force=False) + + prefix = cache_dir / f"palace-cpu-v{tag}" + assert result == prefix / "bin" / "palace" + assert (prefix / "bin" / "palace").is_file() + assert (prefix / "lib" / "libfoo.so").is_file() + assert os.access(result, os.X_OK) + + class TestResolvePalaceLibraryDir: - @pytest.mark.usefixtures("_mock_gcloud", "_no_palacetoolkit") - def test_returns_none_without_palacetoolkit(self) -> None: + @pytest.mark.usefixtures("_mock_gcloud", "_no_local_runtime") + def test_returns_none_without_runtime(self) -> None: from gsim.palace.runtime import resolve_palace_library_dir - with pytest.MonkeyPatch().context() as mp: - mp.setattr("gsim.palace.runtime._palace_cpu_available", lambda: False) - assert resolve_palace_library_dir() is None + assert resolve_palace_library_dir() is None @pytest.mark.usefixtures("_mock_gcloud") def test_delegates_to_palacetoolkit(self) -> None: @@ -679,7 +815,43 @@ def palace_library_path() -> Path: with pytest.MonkeyPatch().context() as mp: mp.setitem(sys.modules, "palacetoolkit_palace_cpu", _FakePalaceCPU()) mp.setattr("gsim.palace.runtime._palace_cpu_available", lambda: True) + mp.setattr("gsim.palace.runtime._palace_toolkit_available", lambda: False) + mp.setattr("pathlib.Path.is_dir", lambda _: True) + result = resolve_palace_library_dir() + assert result is not None + + @pytest.mark.usefixtures("_mock_gcloud") + def test_uses_gsim_cached_library_dir(self) -> None: + from gsim.palace.runtime import resolve_palace_library_dir + + fake_lib = Path("/home/user/.cache/palacetoolkit/runtime/lib") + + with pytest.MonkeyPatch().context() as mp: + mp.setattr("gsim.palace.runtime._palace_cpu_available", lambda: False) + mp.setattr("gsim.palace.runtime._palace_toolkit_available", lambda: False) + mp.setattr("gsim.palace.runtime._cached_library_dir", lambda: fake_lib) + result = resolve_palace_library_dir() + assert result == fake_lib.resolve() + + @pytest.mark.usefixtures("_mock_gcloud") + def test_delegates_to_palacetoolkit_package(self) -> None: + from gsim.palace.runtime import resolve_palace_library_dir + + fake_lib = Path("/opt/palacetoolkit/runtime/lib") + + import types + + ptk = types.ModuleType("palacetoolkit") + ptk.__path__ = [] # type: ignore[attr-defined] + ptk_runtime = types.ModuleType("palacetoolkit.palace_runtime") + setattr(ptk_runtime, "resolve_palace_library_dir", lambda: fake_lib) # noqa: B010 + + with pytest.MonkeyPatch().context() as mp: + mp.setattr("gsim.palace.runtime._palace_cpu_available", lambda: False) + mp.setattr("gsim.palace.runtime._palace_toolkit_available", lambda: True) mp.setattr("pathlib.Path.is_dir", lambda _: True) + mp.setitem(sys.modules, "palacetoolkit", ptk) + mp.setitem(sys.modules, "palacetoolkit.palace_runtime", ptk_runtime) result = resolve_palace_library_dir() assert result is not None @@ -687,11 +859,13 @@ def palace_library_path() -> Path: class TestPalacetoolkitAvailable: @pytest.mark.usefixtures("_mock_gcloud") def test_true_when_installed(self) -> None: - from gsim.palace.runtime import _palace_cpu_available - - # Since we mocked gcloud but not palacetoolkit_palace_cpu, if it's - # actually installed on the system, this will be True. We can't force - # it to be True via mock here without patching importlib, which is - # fragile. Instead we just verify the function runs. - result = _palace_cpu_available() - assert isinstance(result, bool) + from gsim.palace.runtime import ( + _palace_cpu_available, + _palace_toolkit_available, + ) + + # If either package is actually installed on the system, this will be + # True. We can't force it via mock here without patching importlib, + # which is fragile. Instead we just verify the functions run. + assert isinstance(_palace_cpu_available(), bool) + assert isinstance(_palace_toolkit_available(), bool) From 26d1d33c99e0cce9b7ed9b6b22e2a6893ab01149 Mon Sep 17 00:00:00 2001 From: mdmaas Date: Wed, 19 Aug 2026 21:35:56 -0300 Subject: [PATCH 2/4] fix pre-commit --- pyproject.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 5e16c5a9..069e1e94 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -275,7 +275,6 @@ select = ["ALL"] ] "tests/**/*.py" = [ "ANN", # flake8-annotations - "ARG001", # unused-function-argument "ARG005", # unused-lambda-argument (test stubs) "D", # pydocstyle "INP001", # implicit-namespace-package From a9068341eeeaef8a55327ee4a7708b809073ff51 Mon Sep 17 00:00:00 2001 From: mdmaas Date: Wed, 19 Aug 2026 22:25:53 -0300 Subject: [PATCH 3/4] fixed macOS and Windows tests: no auto-download available --- src/gsim/palace/runtime.py | 20 ++++++++++++-------- tests/palace/test_sim_classes.py | 5 ++++- 2 files changed, 16 insertions(+), 9 deletions(-) diff --git a/src/gsim/palace/runtime.py b/src/gsim/palace/runtime.py index 47195eab..b75e46ff 100644 --- a/src/gsim/palace/runtime.py +++ b/src/gsim/palace/runtime.py @@ -110,20 +110,19 @@ def _set_executable(path: Path) -> None: def install_palace_runtime(force: bool = False, timeout: float = 180.0) -> Path: - """Download and cache the prebuilt Palace CPU runtime. + """Return a Palace runtime, downloading and caching it if needed. + + An already-cached runtime is returned on any platform; only the download + itself is restricted to Linux x86_64 (the only platform the prebuilt + Palace CPU wheel is provided for). Returns: Path to the cached ``palace`` launcher executable. Raises: - RuntimeError: If the platform is unsupported or the download/install - fails. + RuntimeError: If no runtime is cached, the platform is unsupported, + or the download/install fails. """ - if not _is_linux_x86_64(): - raise RuntimeError( - "Prebuilt runtime download is only supported on Linux x86_64" - ) - tag = _binary_tag() prefix = _cached_runtime_prefix(tag) bin_palace = prefix / "bin" / "palace" @@ -131,6 +130,11 @@ def install_palace_runtime(force: bool = False, timeout: float = 180.0) -> Path: if not force and bin_palace.is_file() and lib_dir.is_dir(): return bin_palace + if not _is_linux_x86_64(): + raise RuntimeError( + "Prebuilt runtime download is only supported on Linux x86_64" + ) + prefix.mkdir(parents=True, exist_ok=True) downloads = _runtime_cache_dir() / "downloads" downloads.mkdir(parents=True, exist_ok=True) diff --git a/tests/palace/test_sim_classes.py b/tests/palace/test_sim_classes.py index 6961c292..a11f7bb6 100644 --- a/tests/palace/test_sim_classes.py +++ b/tests/palace/test_sim_classes.py @@ -739,10 +739,13 @@ def test_returns_cached_binary_when_present(self, tmp_path: Path) -> None: assert result == bin_palace @pytest.mark.usefixtures("_mock_gcloud") - def test_raises_on_non_linux_x86_64(self) -> None: + def test_raises_on_non_linux_x86_64(self, tmp_path: Path) -> None: import gsim.palace.runtime as rt with pytest.MonkeyPatch().context() as mp: + # Use an empty cache dir so the fallthrough to the platform guard + # is deterministic regardless of what is cached on the host. + mp.setattr(rt, "_runtime_cache_dir", lambda: tmp_path) mp.setattr(rt, "_is_linux_x86_64", lambda: False) with pytest.raises(RuntimeError): rt.install_palace_runtime() From 22fee6048b89a68b357a2c2ca27562badddbc563 Mon Sep 17 00:00:00 2001 From: mdmaas Date: Tue, 25 Aug 2026 23:45:31 -0300 Subject: [PATCH 4/4] feat: PN junction depletion model with auto capacitance/high-res modes --- CHANGELOG.md | 10 + docs/api/common.md | 32 ++ nbs/palace_2d_twmzm.ipynb | 330 ++++++++++++------- pyproject.toml | 2 +- src/gsim/common/cross_section.py | 6 + src/gsim/common/stack/__init__.py | 17 +- src/gsim/common/stack/doping.py | 224 ++++++++++++- src/gsim/common/stack/junction.py | 439 +++++++++++++++++++++++++ src/gsim/palace/base.py | 62 ++++ tests/common/test_cross_section.py | 38 +++ tests/common/test_junction_physics.py | 215 ++++++++++++ tests/common/test_junction_profile.py | 180 ++++++++++ tests/palace/test_pn_junction_modes.py | 159 +++++++++ 13 files changed, 1598 insertions(+), 116 deletions(-) create mode 100644 src/gsim/common/stack/junction.py create mode 100644 tests/common/test_junction_physics.py create mode 100644 tests/common/test_junction_profile.py create mode 100644 tests/palace/test_pn_junction_modes.py diff --git a/CHANGELOG.md b/CHANGELOG.md index d15c3ff3..2574bc8d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,15 @@ # Changelog +## Unreleased + +- PN-junction depletion model from Sze *Physics of Semiconductor Devices* (`PNJunctionConfig`, + `make_pn_junction_profile`): computes built-in voltage, depletion width `W` (abrupt or linearly graded), asymmetric + P/N split `x_p`/`x_n`, and capacitance `C_j = eps_s A / W`. The depletion region is represented automatically — meshed + as a dielectric strip in high-res mode when `W >= ~1/5` of the flanking doped sections, otherwise applied as a lumped + Impedance boundary via `sim.set_pn_junction()`. The 2D TWMZM demo now illustrates both modes. +- Fix: `build_doped_cross_section()` now registers doping/rib materials on `stack.materials`; previously doped domains + silently resolved to eps=1.0 without conductivity in generated Palace configs. + ## 0.1.0 - Electrostatic simulation end-to-end for Palace ([#146](https://github.com/gdsfactory/gsim/pull/146)) diff --git a/docs/api/common.md b/docs/api/common.md index 332cd738..514e65d1 100644 --- a/docs/api/common.md +++ b/docs/api/common.md @@ -38,6 +38,38 @@ inherited_members: false members: false +## PN Junction + +Depletion model after Sze & Ng, *Physics of Semiconductor Devices*, ch. 2. + +::: gsim.common.stack.PNJunctionConfig + options: + show_source: false + +::: gsim.common.stack.make_pn_junction_profile + options: + show_source: false + +::: gsim.common.stack.built_in_voltage + options: + show_source: false + +::: gsim.common.stack.depletion_width + options: + show_source: false + +::: gsim.common.stack.depletion_extents + options: + show_source: false + +::: gsim.common.stack.junction_capacitance_per_area + options: + show_source: false + +::: gsim.common.stack.select_junction_mode + options: + show_source: false + ## Visualization ::: gsim.common.viz.plot_prisms_3d diff --git a/nbs/palace_2d_twmzm.ipynb b/nbs/palace_2d_twmzm.ipynb index bf3bf6a9..0806ecc9 100644 --- a/nbs/palace_2d_twmzm.ipynb +++ b/nbs/palace_2d_twmzm.ipynb @@ -30,7 +30,7 @@ "id": "1", "metadata": {}, "source": [ - "## Geometry & materials parameters\n", + "## Geometry & materials parameters\n", "\n", "Define the layout dimensions, the substrate/metal stack, and the doping\n", "(geometry + material) for the rib and graded slab regions. These are consumed\n", @@ -44,6 +44,7 @@ "- `WG` (1,0): Waveguide core (220 nm Si, 400 nm wide)\n", "- `SLAB90` (3,0): 90 nm slab regions\n", "- `N` (20,0) / `P` (21,0): PN junction doping\n", + "- `(22,0)`: Depletion/junction strip (drawn only in high-res mode)\n", "- `NPP` (24,0) / `PP` (23,0): N+/P+ graded contact doping (via `make_doping_profile`)\n", "- `M1` (41,0): CPW electrodes (Al, 1 um thick)\n" ] @@ -80,10 +81,19 @@ "METAL1_ZMIN = 1.1 # metal1 bottom (top of the oxide stack)\n", "METAL1_THICKNESS = 1.0 # CPW electrode thickness on metal1\n", "\n", - "# --- PN junction / doping material model -------------------------------------\n", + "# --- PN junction depletion model (Sze ch. 2) ---------------------------------\n", + "# W = sqrt(2 eps_s (V_bi + V_R)/q * (Na+Nd)/(Na Nd)); x_p/x_n split the\n", + "# depletion into the P/N sides; C_j = eps_s A / W. Doping in cm^-3.\n", "SI_PERMITTIVITY = 11.9\n", "FMAX_RF_MATERIAL = 200e9 # validity range of the constant-eps doping models (Hz)\n", - "RIB_DOPING_SIGMA = 1.6e3 # p_rib / n_rib junction conductivity (S/m)\n", + "PN_RIB_SIGMA = 1.6e3 # Drude conductivity of the P/N rib regions (S/m)\n", + "PN_JUNCTION = {\n", + " \"na_cm3\": 1e19,\n", + " \"nd_cm3\": 1e19,\n", + " \"v_reverse\": 0.0,\n", + " \"permittivity\": SI_PERMITTIVITY,\n", + "}\n", + "JUNCTION_GDS_LAYER = (22, 0) # depletion-strip GDS layer (drawn only in high-res)\n", "\n", "# Graded slab doping {side: [(width_um, sigma_S_per_m), ...]}, from the rib edge.\n", "DOPING_PROFILE = {\n", @@ -104,6 +114,68 @@ "cell_type": "markdown", "id": "3", "metadata": {}, + "source": [ + "## PN-junction width & automatic representation mode\n", + "\n", + "`make_pn_junction_profile()` splits the rib into P / depletion / N regions\n", + "using the depletion approximation (Sze & Ng, *Physics of Semiconductor\n", + "Devices*, 3rd ed., ch. 2):\n", + "\n", + "$$V_{bi} = \\frac{k_B T}{q}\\ln\\frac{N_A N_D}{n_i^2}, \\qquad\n", + "W = \\sqrt{\\frac{2 \\varepsilon_s (V_{bi}+V_R)}{q}\\frac{N_A+N_D}{N_A N_D}}, \\qquad\n", + "C_j = \\frac{\\varepsilon_s}{W}$$\n", + "\n", + "Two representations, selected automatically from $W$:\n", + "\n", + "- **capacitance**: $W$ is far thinner than the neighbouring doped sections\n", + " ($W < \\tfrac{1}{5}$ of a flank). The P/N geometry stays unchanged and the\n", + " computed $C_j$ is applied as a lumped Impedance boundary.\n", + "- **high_res**: $W$ is comparable to the flanks. The depletion strip is drawn\n", + " as a contiguous dielectric rectangle ($\\varepsilon_s$, no carriers) that is\n", + " resolved on the actual mesh.\n", + "\n", + "Both modes return the doped regions, the P/N regions, and the junction\n", + "metadata (widths, capacitance, chosen mode).\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "from gsim.common.stack.junction import PNJunctionConfig\n", + "\n", + "junc = PNJunctionConfig.model_validate(PN_JUNCTION)\n", + "flank = RIB_WIDTH / 2\n", + "print(\n", + " f\"Bias point: V_bi = {junc.v_bi:.3f} V, W = {junc.w_um * 1e3:.1f} nm \"\n", + " f\"(x_p = {junc.xp_um * 1e3:.1f} nm, x_n = {junc.xn_um * 1e3:.1f} nm)\"\n", + ")\n", + "print(\n", + " f\"C_j = eps_s A / W = {junc.capacitance(LENGTH, RIB_HEIGHT) * 1e15:.2f} fF \"\n", + " f\"(A = {LENGTH} x {RIB_HEIGHT} um)\"\n", + ")\n", + "print(\n", + " f\"Flank size = {flank * 1e3:.0f} nm -> auto mode selects \"\n", + " f\"'{junc.select_mode(flank, flank)}'\\n\"\n", + ")\n", + "\n", + "print(\"How doping moves W across the auto-selection threshold:\")\n", + "for n_cm3 in (1e19, 5e18, 2e18, 1e18):\n", + " j = PNJunctionConfig(na_cm3=n_cm3, nd_cm3=n_cm3)\n", + " mode = j.select_mode(flank, flank)\n", + " print(\n", + " f\" Na = Nd = {n_cm3:.1e} cm^-3 : W = {j.w_um * 1e3:6.1f} nm, \"\n", + " f\"C = {j.capacitance(LENGTH, RIB_HEIGHT) * 1e15:6.2f} fF -> {mode}\"\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, "source": [ "## Build TW-MZM cross-section geometry\n", "\n", @@ -114,14 +186,14 @@ { "cell_type": "code", "execution_count": null, - "id": "4", + "id": "6", "metadata": {}, "outputs": [], "source": [ "import gdsfactory as gf\n", "\n", "from gsim.common.cross_section import build_optical_cross_section\n", - "from gsim.common.stack.doping import make_doping_profile\n", + "from gsim.common.stack.doping import make_doping_profile, make_pn_junction_profile\n", "\n", "gf.gpdk.PDK.activate()\n", "\n", @@ -134,12 +206,14 @@ " return r\n", "\n", "\n", - "def _add_device_core(comp: gf.Component) -> None:\n", + "def _add_device_core(comp: gf.Component) -> dict:\n", " \"\"\"Rib + slab + PN junction — shared by the RF and optical components.\n", "\n", - " The P/N rectangles are the same \"doping profile\" polygons in both, so the\n", - " optical cross-section still shows the junction shape. The optical stack\n", - " maps all four regions to plain silicon.\n", + " The PN-junction regions come from ``make_pn_junction_profile()``: the\n", + " depletion width W (and its x_p/x_n split) follows from the configured\n", + " doping/bias, and the representation mode is auto-selected. In\n", + " capacitance mode the P/N rectangles stay adjacent; in high-res mode a\n", + " contiguous depleted-Si strip of width W is drawn between them.\n", " \"\"\"\n", " # 1. Rib waveguide core (PN junction sits inside it)\n", " wg = comp << centered_rect(LENGTH, RIB_WIDTH, LAYER.WG)\n", @@ -149,17 +223,26 @@ " slab = comp << centered_rect(LENGTH, 2 * SLAB_HALF + RIB_WIDTH, LAYER.SLAB90)\n", " slab.y = 0.0\n", "\n", - " # 3. PN junction (P above / N below the rib centre)\n", - " p_half = comp << gf.c.rectangle((LENGTH, RIB_WIDTH / 2), layer=LAYER.P)\n", - " p_half.y = RIB_CENTER_Y + RIB_WIDTH / 4\n", - " n_half = comp << gf.c.rectangle((LENGTH, RIB_WIDTH / 2), layer=LAYER.N)\n", - " n_half.y = RIB_CENTER_Y - RIB_WIDTH / 4\n", + " # 3. PN junction: P / depletion / N regions around the rib centre\n", + " return make_pn_junction_profile(\n", + " comp,\n", + " length=LENGTH,\n", + " center_y=RIB_CENTER_Y,\n", + " rib_width=RIB_WIDTH,\n", + " junction=PN_JUNCTION,\n", + " p_region=(\"p_rib\", tuple(LAYER.P), PN_RIB_SIGMA),\n", + " n_region=(\"n_rib\", tuple(LAYER.N), PN_RIB_SIGMA),\n", + " junction_region=(\"junction\", JUNCTION_GDS_LAYER),\n", + " zmin=0.0,\n", + " zmax=RIB_HEIGHT,\n", + " fmax=FMAX_RF_MATERIAL,\n", + " )\n", "\n", "\n", - "def _build_rf_component() -> tuple[gf.Component, dict]:\n", + "def _build_rf_component() -> tuple[gf.Component, dict, dict]:\n", " \"\"\"Full TW-MZM cross-section: device core + graded doping + CPW + vias.\"\"\"\n", " comp = gf.Component()\n", - " _add_device_core(comp)\n", + " pn_result = _add_device_core(comp)\n", "\n", " # 4. Graded N+/P+ slab doping (contiguous, no gaps)\n", " doping_result = make_doping_profile(\n", @@ -195,35 +278,30 @@ " via_g_to_n.x = 0.0\n", " via_g_to_n.y = VIA_G_TO_N_Y\n", "\n", - " return comp, doping_result\n", + " return comp, doping_result, pn_result\n", "\n", "\n", - "def _build_optical_component() -> gf.Component:\n", + "def _build_optical_component() -> tuple[gf.Component, dict]:\n", " \"\"\"Optical-only cross-section: rib + slab + PN junction (all silicon).\n", "\n", " No electrodes, vias, or graded doping — the optical mode sees a single\n", " homogeneous Si body embedded in the uniform SiO2 cladding stack.\n", " \"\"\"\n", " comp = gf.Component()\n", - " _add_device_core(comp)\n", - " return comp\n", - "\n", + " pn_result = _add_device_core(comp)\n", + " return comp, pn_result\n", "\n", - "# --- RF component (electrodes, vias, graded doping) ------------------------\n", - "comp, doping_result = _build_rf_component()\n", "\n", - "# --- Optical component (rib + slab + PN junction only) ---------------------\n", - "comp_optical = _build_optical_component()\n", + "# --- RF component (electrodes, vias, graded doping) ---------------------------\n", + "comp, doping_result, pn_result = _build_rf_component()\n", "\n", - "# -- Plot ----------------------------------------------------------------------\n", - "_cc = comp.copy()\n", - "_cc.draw_ports()\n", - "_cc.plot()" + "# --- Optical-only component ----------------------------------------------------\n", + "comp_optical, pn_result_optical = _build_optical_component()" ] }, { "cell_type": "markdown", - "id": "5", + "id": "7", "metadata": {}, "source": [ "## Inspect 2D cross-section\n", @@ -236,7 +314,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6", + "id": "8", "metadata": {}, "outputs": [], "source": [ @@ -247,29 +325,36 @@ "\n", "from gsim.common.cross_section import build_doped_cross_section\n", "\n", + "# The reusable `gsim.common.cross_section.build_doped_cross_section()` helper\n", + "# assembles the base PDK stack and overrides `metal1`. The graded-slab doping\n", + "# specs and the PN-junction specs (P/N plus the depletion strip when meshed)\n", + "# are merged into a single `doping=` input.\n", + "\n", + "doping_input = {\n", + " \"layer_specs\": {**doping_result[\"layer_specs\"], **pn_result[\"layer_specs\"]},\n", + " \"materials\": {**doping_result[\"materials\"], **pn_result[\"materials\"]},\n", + "}\n", + "\n", "stack, section = build_doped_cross_section(\n", " comp,\n", " axis=CROSS_SECTION_AXIS,\n", " value=CROSS_SECTION_VALUE,\n", " substrate_thickness=BOX_THICKNESS,\n", " metal1=(METAL1_ZMIN, METAL1_THICKNESS),\n", - " doping=doping_result,\n", - " rib_layers=[\n", - " (\"p_rib\", LAYER.P, RIB_DOPING_SIGMA),\n", - " (\"n_rib\", LAYER.N, RIB_DOPING_SIGMA),\n", - " ],\n", - " rib_height=RIB_HEIGHT,\n", + " doping=doping_input,\n", " permittivity=SI_PERMITTIVITY,\n", " fmax=FMAX_RF_MATERIAL,\n", - ")" + ")\n", + "\n", + "print(\"PN-junction representation:\", pn_result[\"junction\"][\"mode\"])" ] }, { "cell_type": "markdown", - "id": "7", + "id": "9", "metadata": {}, "source": [ - "## Optical-only cross-section\n", + "### Optical-only cross-section\n", "\n", "The optical analysis uses a simplified component: the same rib + slab + PN\n", "junction (identical \"doping profile\" polygons) but **no** electrodes, vias, or\n", @@ -277,29 +362,36 @@ "sees one homogeneous Si body embedded in a uniform SiO2 cladding.\n", "\n", "`gsim.common.cross_section.build_optical_cross_section()` assembles the\n", - "minimal all-dielectric `LayerStack` and extracts the 2D cross-section at $x=0$." + "minimal all-dielectric `LayerStack` and extracts the 2D cross-section at $x=0$.\n", + "When the auto-selected PN representation is high-res, the depletion strip is\n", + "included as an extra silicon region.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "8", + "id": "10", "metadata": {}, "outputs": [], "source": [ "# Uniform SiO2 cladding height above z=0 (um)\n", "OPT_CLAD_TOP = 3.0\n", "\n", + "device_layers = {\n", + " \"core\": (LAYER.WG, 0.0, RIB_HEIGHT),\n", + " \"slab\": (LAYER.SLAB90, 0.0, SLAB_THICKNESS),\n", + " \"p_rib\": (LAYER.P, 0.0, RIB_HEIGHT),\n", + " \"n_rib\": (LAYER.N, 0.0, RIB_HEIGHT),\n", + "}\n", + "if \"junction\" in pn_result_optical[\"layer_specs\"]:\n", + " # High-res mode drew the depletion strip; it is plain silicon optically.\n", + " device_layers[\"junction\"] = (JUNCTION_GDS_LAYER, 0.0, RIB_HEIGHT)\n", + "\n", "stack_opt, section_opt = build_optical_cross_section(\n", " comp_optical,\n", " axis=CROSS_SECTION_AXIS,\n", " value=CROSS_SECTION_VALUE,\n", - " device_layers={\n", - " \"core\": (LAYER.WG, 0.0, RIB_HEIGHT),\n", - " \"slab\": (LAYER.SLAB90, 0.0, SLAB_THICKNESS),\n", - " \"p_rib\": (LAYER.P, 0.0, RIB_HEIGHT),\n", - " \"n_rib\": (LAYER.N, 0.0, RIB_HEIGHT),\n", - " },\n", + " device_layers=device_layers,\n", " substrate_thickness=BOX_THICKNESS,\n", " cladding_top=OPT_CLAD_TOP,\n", ")\n", @@ -310,7 +402,7 @@ }, { "cell_type": "markdown", - "id": "9", + "id": "11", "metadata": {}, "source": [ "### Optical cross-section plot\n", @@ -322,7 +414,7 @@ { "cell_type": "code", "execution_count": null, - "id": "10", + "id": "12", "metadata": {}, "outputs": [], "source": [ @@ -354,7 +446,7 @@ }, { "cell_type": "markdown", - "id": "11", + "id": "13", "metadata": {}, "source": [ "### Optical material properties (1550 nm)\n", @@ -367,7 +459,7 @@ { "cell_type": "code", "execution_count": null, - "id": "12", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -389,7 +481,7 @@ }, { "cell_type": "markdown", - "id": "13", + "id": "15", "metadata": {}, "source": [ "## Plot the 2D cross-section\n", @@ -402,7 +494,7 @@ { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "16", "metadata": {}, "outputs": [], "source": [ @@ -430,7 +522,7 @@ { "cell_type": "code", "execution_count": null, - "id": "15", + "id": "17", "metadata": {}, "outputs": [], "source": [ @@ -451,7 +543,7 @@ }, { "cell_type": "markdown", - "id": "16", + "id": "18", "metadata": {}, "source": [ "## RF simulation (50 GHz)\n", @@ -463,7 +555,7 @@ { "cell_type": "code", "execution_count": null, - "id": "17", + "id": "19", "metadata": {}, "outputs": [], "source": [ @@ -499,14 +591,13 @@ " \"pp_slab_1\",\n", "]\n", "\n", - "# --- PN junction lumped model -------------------------------------------------\n", - "PN_JUNCTION_CAPACITANCE = 1e-15 # F, on the p_rib / n_rib interface" + "# --- PN junction lumped model -------------------------------------------------" ] }, { "cell_type": "code", "execution_count": null, - "id": "18", + "id": "20", "metadata": {}, "outputs": [], "source": [ @@ -541,7 +632,7 @@ { "cell_type": "code", "execution_count": null, - "id": "19", + "id": "21", "metadata": {}, "outputs": [], "source": [ @@ -556,12 +647,31 @@ { "cell_type": "code", "execution_count": null, - "id": "20", + "id": "22", "metadata": {}, "outputs": [], "source": [ - "# Generate Palace config file (mesh must be present)\n", - "sim.add_impedance_boundary(\"p_rib\", \"n_rib\", capacitance=PN_JUNCTION_CAPACITANCE)\n", + "# Generate Palace config file (mesh must be present).\n", + "#\n", + "# Capacitance mode: apply C_j = eps_s A / W as a lumped Impedance boundary\n", + "# on the p_rib/n_rib interface via `set_pn_junction()`.\n", + "# High-res mode: the depletion strip already exists as dielectric geometry on\n", + "# the mesh — no lumped boundary is needed (and adding one would double-count).\n", + "if pn_result[\"junction\"][\"mode\"] == \"capacitance\":\n", + " applied_c = sim.set_pn_junction(\n", + " PN_JUNCTION,\n", + " layer_p=\"p_rib\",\n", + " layer_n=\"n_rib\",\n", + " length_um=LENGTH,\n", + " height_um=RIB_HEIGHT,\n", + " )\n", + " print(f\"Applied lumped junction capacitance: {applied_c * 1e15:.2f} fF\")\n", + "else:\n", + " print(\n", + " f\"High-res mode: depletion strip W = \"\n", + " f\"{pn_result['junction']['w_um'] * 1e3:.1f} nm meshed as dielectric.\"\n", + " )\n", + "\n", "sim.write_config()\n", "print(\"Config written to:\", sim.output_dir)" ] @@ -569,7 +679,7 @@ { "cell_type": "code", "execution_count": null, - "id": "21", + "id": "23", "metadata": {}, "outputs": [], "source": [ @@ -581,7 +691,7 @@ { "cell_type": "code", "execution_count": null, - "id": "22", + "id": "24", "metadata": {}, "outputs": [], "source": [ @@ -611,7 +721,7 @@ { "cell_type": "code", "execution_count": null, - "id": "23", + "id": "25", "metadata": {}, "outputs": [], "source": [ @@ -629,7 +739,7 @@ }, { "cell_type": "markdown", - "id": "24", + "id": "26", "metadata": {}, "source": [ "## Optical simulation (1550 nm)\n", @@ -649,7 +759,7 @@ { "cell_type": "code", "execution_count": null, - "id": "25", + "id": "27", "metadata": { "lines_to_next_cell": 2 }, @@ -680,7 +790,7 @@ { "cell_type": "code", "execution_count": null, - "id": "26", + "id": "28", "metadata": {}, "outputs": [], "source": [ @@ -735,7 +845,7 @@ { "cell_type": "code", "execution_count": null, - "id": "27", + "id": "29", "metadata": {}, "outputs": [], "source": [ @@ -750,7 +860,7 @@ { "cell_type": "code", "execution_count": null, - "id": "28", + "id": "30", "metadata": {}, "outputs": [], "source": [ @@ -761,7 +871,7 @@ { "cell_type": "code", "execution_count": null, - "id": "29", + "id": "31", "metadata": {}, "outputs": [], "source": [ @@ -779,7 +889,7 @@ { "cell_type": "code", "execution_count": null, - "id": "30", + "id": "32", "metadata": {}, "outputs": [], "source": [ @@ -795,7 +905,7 @@ }, { "cell_type": "markdown", - "id": "31", + "id": "33", "metadata": {}, "source": [ "## Summary\n", @@ -803,53 +913,51 @@ "The geometry has been built and meshed for both RF (50 GHz) and optical (193 THz / 1550 nm) analysis.\n", "\n", "**Cross-section elements:**\n", - "| Component | Layer | y-range (um) | z-range (um) | Material | sigma (S/m) |\n", - "|---|---|---|---|---|---|\n", - "| Rib core | WG (1,0) | [-20.2, -19.8] | [0, 0.22] | Si (intrinsic) | 2 |\n", - "| Slab (90 nm) | SLAB90 (3,0) | [-40.2, +0.2] | [0, 0.09] | Si (intrinsic) | 2 |\n", - "| PN junction (P) | P (21,0) | [-20.0, -19.8] | [0, 0.22] | doped Si (p_rib) | 1.6x10^3 |\n", - "| PN junction (N) | N (20,0) | [-20.2, -20.0] | [0, 0.22] | doped Si (n_rib) | 1.6x10^3 |\n", - "| P+ graded inner | PP (23,0) | [-18.3, -16.3] | [0, 0.09] | doped Si (pp_slab_0) | 2x10^4 |\n", - "| P+ graded outer | PP (23,1) | [-15.3, -13.3] | [0, 0.09] | doped Si (pp_slab_1) | 8x10^4 |\n", - "| N+ graded inner | NPP (24,0) | [-21.7, -23.7] | [0, 0.09] | doped Si (npp_slab_0) | 2x10^4 |\n", - "| N+ graded outer | NPP (24,1) | [-24.7, -26.7] | [0, 0.09] | doped Si (npp_slab_1) | 8x10^4 |\n", - "| Vias (S to P+) | VIAC/VIA1/VIA2 | [-15.0, -9.0] | [0.09, 3.2] | W/Al | 3.5x10^7 |\n", - "| Vias (G to N+) | VIAC/VIA1/VIA2 | [-31.0, -25.0] | [0.09, 3.2] | W/Al | 3.5x10^7 |\n", - "| CPW signal | M1 (41,0) | [-10, +10] | [3.2, 4.2] | Al (1 um) | 3.5x10^7 |\n", - "| CPW ground (top) | M1 (41,0) | [+30, +70] | [3.2, 4.2] | Al (1 um) | 3.5x10^7 |\n", - "| CPW ground (bot) | M1 (41,0) | [-70, -30] | [3.2, 4.2] | Al (1 um) | 3.5x10^7 |\n", + "| Component | Layer | Material | sigma (S/m) |\n", + "|---|---|---|---|\n", + "| Rib core | WG (1,0) | Si (intrinsic) | 2 |\n", + "| Slab (90 nm) | SLAB90 (3,0) | Si (intrinsic) | 2 |\n", + "| PN junction (P) | P (21,0) | doped Si (p_rib) | 1.6x10^3 |\n", + "| Depletion strip | (22,0) | Si (eps 11.9, undoped) | — |\n", + "| PN junction (N) | N (20,0) | doped Si (n_rib) | 1.6x10^3 |\n", + "| P+ graded inner/outer | PP (23,0)/(23,1) | doped Si | 2e4 / 8e4 |\n", + "| N+ graded inner/outer | NPP (24,0)/(24,1) | doped Si | 2e4 / 8e4 |\n", + "| Vias (S->P+, G->N+) | VIAC/VIA1/VIA2 | W/Al | 3.5x10^7 |\n", + "| CPW signal / grounds | M1 (41,0) | Al (1 um) | 3.5x10^7 |\n", + "\n", + "**PN-junction model** (Sze & Ng, *Physics of Semiconductor Devices*, 3rd ed., ch. 2):\n", + "- Built-in voltage $V_{bi} = \\frac{k_BT}{q}\\ln(N_AN_D/n_i^2)$.\n", + "- Abrupt-junction depletion width $W = \\sqrt{\\frac{2\\varepsilon_s(V_{bi}+V_R)}{q}\\frac{N_A+N_D}{N_A N_D}}$,\n", + " split asymmetrically $x_p = W N_D/(N_A+N_D)$ into P and $x_n = W N_A/(N_A+N_D)$ into N;\n", + " linearly graded junctions use $W = [12\\varepsilon_s(V_{bi}+V_R)/(qa)]^{1/3}$.\n", + "- Junction capacitance $C_j = \\varepsilon_s A / W$.\n", + "\n", + "**Representation modes (auto-selected from $W$ vs the flanking doped sections):**\n", + "- **capacitance**: $W$ below ~1/5 of a flank -> P/N geometry unchanged; $C_j$\n", + " applied as a lumped Impedance boundary (`sim.set_pn_junction()`).\n", + "- **high_res**: $W$ comparable to the flanks -> contiguous depleted-Si strip of\n", + " width $W$ drawn between P and N ($\\varepsilon_s$, no carriers) and resolved\n", + " on the actual mesh; no lumped boundary.\n", "\n", "**Material modelling notes:**\n", - "- Doping regions are modelled as **semiconductors** (finite sigma from the Drude free-carrier model), not metals. This avoids short-circuiting the PN junction.\n", - "- Conductivities are derived from $\\sigma = q\\mu N$ with typical dopant concentrations ($N \\sim 10^{19}\\ \\text{cm}^{-3}$ for the junction, $\\sim 10^{20}\\ \\text{cm}^{-3}$ for the contacts).\n", - "- Doping on each side of the rib uses a **configurable piecewise gradient** via `make_doping_profile()`, and the whole cross-section assembly is wrapped by `build_doped_cross_section()`.\n", - "- The **depletion region** and voltage-dependent capacitance are NOT modelled here — this is a linear small-signal analysis at a fixed bias point.\n", - "- The **plasma-dispersion effect** is not applied to the optical simulation; the rib is treated as intrinsic Si at 1550 nm.\n", + "- Doped regions are modelled as **semiconductors** (finite Drude $\\sigma$),\n", + " not metals, so they do not short-circuit the junction.\n", + "- The depletion region is now represented either lumped or geometrically\n", + " (auto-selected); earlier revisions omitted it entirely.\n", + "- The **plasma-dispersion effect** is not applied to the optical simulation;\n", + " the rib is treated as intrinsic Si at 1550 nm.\n", "\n", "**Optical-only component:**\n", "- The optical run (`sim_optical`) uses a **simplified component** — the same\n", - " rib + slab + PN junction (\"doping profile\") polygons, but no electrodes,\n", - " vias, or graded doping.\n", + " rib + slab + PN junction polygons, but no electrodes, vias, or graded doping.\n", "- Every device region maps to **plain silicon** in `build_optical_cross_section()`,\n", - " so the optical mode sees one homogeneous Si body in a **uniform SiO2 cladding**.\n", - "- The background medium is SiO2, not air: `set_airbox(material=\"sio2\", ...)`\n", - " fills the padded 2D domain with the cladding material (default is air).\n", - "- Material dispersion is evaluated at `F_OPT` by the boundary-mode config\n", - " generator: Si -> eps~12.09 (n~3.478), SiO2 -> eps~2.09 (n~1.444) at 1550 nm.\n", - "\n", - "**Next steps (user action):**\n", - "1. Verify the zoomed cross-section plot shows the rib (centred at y=-20), PN junction, graded doping, and vias.\n", - "2. To run locally, provide a Palace CPU runner via `PALACE_BIN` or as `palace` on PATH. 2D\n", - " mode analysis defaults to a single MPI rank + OpenMP threads; pass `num_processes=1` explicitly if you want\n", - " to be explicit about it.\n", - "3. Run `sim_optical.run_local(verbose=True)` for the optical mode.\n", - "4. Use `gsim.palace.plot_fields_2d()` to visualise mode profiles, and `gsim.palace.plot_plane_section()` for cross-section physical groups.\n" + " so the optical mode sees one homogeneous Si body in uniform SiO2.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "32", + "id": "34", "metadata": {}, "outputs": [], "source": [ diff --git a/pyproject.toml b/pyproject.toml index 069e1e94..2b27102d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -89,7 +89,7 @@ requires = ["build", "setuptools>=61", "uv", "wheel"] build-backend = "setuptools.build_meta" [tool.codespell] -ignore-words-list = "doubleclick,euclidian,te" +ignore-words-list = "doubleclick,euclidian,te,nd" [tool.interrogate] docstring-style = "google" diff --git a/src/gsim/common/cross_section.py b/src/gsim/common/cross_section.py index 170d38e5..61b094ef 100644 --- a/src/gsim/common/cross_section.py +++ b/src/gsim/common/cross_section.py @@ -305,6 +305,12 @@ def build_doped_cross_section( for name, layer in layer_specs.items(): stack.layers[name] = layer + # Register the doping/rib materials on the stack so downstream consumers + # (Palace config generator, Meep, ...) resolve their eps/sigma instead of + # silently falling back to vacuum. + for name, mat in materials.items(): + stack.materials[name] = mat.to_dict() if hasattr(mat, "to_dict") else mat + section = extract_plane_section( component.copy(), stack, diff --git a/src/gsim/common/stack/__init__.py b/src/gsim/common/stack/__init__.py index f947263a..db2b7aa4 100644 --- a/src/gsim/common/stack/__init__.py +++ b/src/gsim/common/stack/__init__.py @@ -23,7 +23,7 @@ import gdsfactory as gf import yaml -from gsim.common.stack.doping import make_doping_profile +from gsim.common.stack.doping import make_doping_profile, make_pn_junction_profile from gsim.common.stack.extractor import ( Layer, LayerStack, @@ -31,6 +31,14 @@ extract_from_pdk, extract_layer_stack, ) +from gsim.common.stack.junction import ( + PNJunctionConfig, + built_in_voltage, + depletion_extents, + depletion_width, + junction_capacitance_per_area, + select_junction_mode, +) from gsim.common.stack.materials import ( MATERIALS_DB, DispersionModel, @@ -170,25 +178,32 @@ def load_stack_yaml(yaml_path: str | Path) -> LayerStack: "LayerStack", "LorentzianTerm", "MaterialProperties", + "PNJunctionConfig", "ResolvedMaterial", "SellmeierTerm", "StackLayer", "ValidationResult", "ValidityRange", + "built_in_voltage", + "depletion_extents", + "depletion_width", "extract_from_pdk", "extract_layer_stack", "get_material_properties", "get_stack", + "junction_capacitance_per_area", "load_overlay", "load_stack_yaml", "make_doped_material", "make_doped_materials", "make_doping_profile", + "make_pn_junction_profile", "merge_overlay", "parse_layer_stack", "plot_stack", "print_stack", "print_stack_table", "resolve_material_at_wavelength", + "select_junction_mode", "should_enable_dispersion", ] diff --git a/src/gsim/common/stack/doping.py b/src/gsim/common/stack/doping.py index 51de63c8..2d8d50ec 100644 --- a/src/gsim/common/stack/doping.py +++ b/src/gsim/common/stack/doping.py @@ -41,15 +41,22 @@ from __future__ import annotations -from typing import TYPE_CHECKING, Any, cast +import logging +from typing import TYPE_CHECKING, Any, Literal, cast import gdsfactory as gf -from gsim.common.stack.materials import make_doped_materials +from gsim.common.stack.junction import ( + JUNCTION_MODE_FRACTION, + PNJunctionConfig, +) +from gsim.common.stack.materials import MaterialProperties, make_doped_materials if TYPE_CHECKING: from gsim.common.stack.extractor import Layer +logger = logging.getLogger(__name__) + _SideConfig = dict[str, dict[str, Any]] @@ -163,4 +170,215 @@ def make_doping_profile( return result -__all__ = ["make_doping_profile"] +def _as_junction_config( + junction: PNJunctionConfig | dict[str, Any], +) -> PNJunctionConfig: + """Accept a config object or plain dict for the junction parameters.""" + if isinstance(junction, PNJunctionConfig): + return junction + return PNJunctionConfig.model_validate(junction) + + +def _add_rect( + comp: gf.Component, + *, + length: float, + y0: float, + y1: float, + gds_layer: tuple[int, int], +) -> float: + """Draw a rectangle spanning ``[y0, y1]`` and return its y-centre.""" + rect = comp << gf.c.rectangle((length, y1 - y0), layer=gds_layer) + rect.y = (y0 + y1) / 2 + return (y0 + y1) / 2 + + +def make_pn_junction_profile( + comp: gf.Component, + *, + length: float, + center_y: float, + rib_width: float, + junction: PNJunctionConfig | dict[str, Any], + p_region: tuple[str, tuple[int, int], float], + n_region: tuple[str, tuple[int, int], float], + junction_region: tuple[str, tuple[int, int]] | None = None, + zmin: float = 0.0, + zmax: float | None = None, + fmax: float = 200e9, + mode: Literal["auto", "capacitance", "high_res"] = "auto", + mode_fraction: float = JUNCTION_MODE_FRACTION, + mesh_resolution: str | float = "fine", +) -> dict[str, dict[str, Any]]: + """Build P / depletion-junction / N rib regions around ``center_y``. + + The depletion width ``W`` (and its asymmetric split ``xp``/``xn`` into + the P and N halves) comes from :class:`PNJunctionConfig`, which + implements the textbook abrupt/linearly-graded junction formulas + (Sze, *Physics of Semiconductor Devices*, ch. 2). + + Two representation modes are supported: + + - ``"high_res"``: three contiguous rectangles are drawn — N + ``[cy - rib_width/2, cy - xn]``, depleted-junction dielectric strip + ``[cy - xn, cy + xp]``, P ``[cy + xp, cy + rib_width/2]``. The + junction strip is registered as a patterned dielectric with a real + GDS layer so it appears on the simulation mesh. + - ``"capacitance"``: geometry is unchanged from a plain P/N split + (adjacent half-rectangles); no junction polygon is drawn and callers + apply the computed capacitance as a lumped impedance boundary instead + (see ``PalaceSimMixin.set_pn_junction``). + + With ``mode="auto"`` the choice falls out of + :func:`gsim.common.stack.junction.select_junction_mode`: the strip is + meshed only when ``W >= mode_fraction * min(P flank, N flank)``, where + each flank is ``rib_width / 2``. + + Args: + comp: gdsfactory component the rectangles are added to. + length: Rectangle length along the propagation direction (um). + center_y: Y coordinate of the metallurgical junction / rib centre. + rib_width: Full rib width (um); P occupies the upper half, N the + lower half. + junction: Depletion-model parameters + (:class:`PNJunctionConfig` or its dict form). + p_region: ``(name, gds_layer, sigma_S_per_m)`` for the P region. + n_region: ``(name, gds_layer, sigma_S_per_m)`` for the N region. + junction_region: ``(name, gds_layer)`` used to register the + depletion strip in high-res mode. Required when the selected + mode is ``"high_res"``; ignored in capacitance mode. + zmin: Bottom z of the regions (um). + zmax: Top z of the regions (um); defaults to ``zmin + 0.22``. + fmax: Upper frequency of the Drude-model validity range (Hz). + mode: ``"auto"``, ``"capacitance"`` or ``"high_res"``. + mode_fraction: Auto-mode threshold fraction (~1/5 default). + mesh_resolution: Mesh resolution assigned to the generated layers. + + Returns: + Dict with keys: + + - ``layer_specs``: ``{name: Layer}`` for every drawn region. + - ``materials``: ``{name: MaterialProperties}`` (Drude models for + P/N, plain dielectric for the junction strip). + - ``centres``: ``{role: y_centre}`` for drawn regions. + - ``junction``: computed quantities (widths, capacitance, chosen + mode and selection reason). + """ + from gsim.common.stack.extractor import Layer + from gsim.common.stack.junction import select_junction_mode + + cfg = _as_junction_config(junction) + p_name, p_layer, p_sigma = p_region + n_name, n_layer, n_sigma = n_region + + ztop = 0.22 if zmax is None else zmax + if ztop <= zmin: + raise ValueError("zmax must exceed zmin.") + if length <= 0: + raise ValueError("length must be positive.") + if cfg.xp_um + cfg.xn_um > rib_width: + raise ValueError( + f"Depletion width W={cfg.w_um:.4g} um does not fit in the " + f"{rib_width:.4g} um rib." + ) + + flank_um = rib_width / 2 + if mode == "auto": + mode = select_junction_mode( + cfg.w_um, flank_um, flank_um, fraction=mode_fraction + ) + reason = ( + f"W={cfg.w_um:.4g} um vs threshold " + f"{mode_fraction * flank_um:.4g} um (= {mode_fraction} * flank)" + ) + else: + reason = f"forced by caller (mode={mode!r})" + logger.info("PN junction mode: %s (%s)", mode, reason) + + result: dict[str, dict[str, Any]] = { + "layer_specs": {}, + "materials": {}, + "centres": {}, + } + layer_specs = cast("dict[str, Layer]", result["layer_specs"]) + materials: dict[str, Any] = result["materials"] + centres: dict[str, float] = result["centres"] + + def _doped_spec(name: str, gds_layer: tuple[int, int], _sigma: float) -> Layer: + return Layer( + name=name, + gds_layer=gds_layer, + zmin=zmin, + zmax=ztop, + thickness=ztop - zmin, + material=name, + layer_type="dielectric", + mesh_resolution=mesh_resolution, + ) + + xp, xn = cfg.xp_um, cfg.xn_um + + # N region: lower half, trimmed by xn when the strip is meshed. + n_y0 = center_y - flank_um + n_y1 = center_y if mode == "capacitance" else center_y - xn + centres["n"] = _add_rect( + comp, length=length, y0=n_y0, y1=n_y1, gds_layer=tuple(n_layer) + ) + layer_specs[n_name] = _doped_spec(n_name, tuple(n_layer), n_sigma) + + # P region: upper half, trimmed by xp when the strip is meshed. + p_y0 = center_y if mode == "capacitance" else center_y + xp + p_y1 = center_y + flank_um + centres["p"] = _add_rect( + comp, length=length, y0=p_y0, y1=p_y1, gds_layer=tuple(p_layer) + ) + layer_specs[p_name] = _doped_spec(p_name, tuple(p_layer), p_sigma) + + materials.update( + make_doped_materials( + [(p_name, p_sigma), (n_name, n_sigma)], + permittivity=cfg.permittivity, + fmax=fmax, + source_prefix="doped Si", + ) + ) + + if mode == "high_res": + if junction_region is None: + raise ValueError( + "mode='high_res' requires junction_region=(name, gds_layer)." + ) + j_name, j_layer = junction_region + centres["junction"] = _add_rect( + comp, + length=length, + y0=center_y - xn, + y1=center_y + xp, + gds_layer=tuple(j_layer), + ) + layer_specs[j_name] = Layer( + name=j_name, + gds_layer=tuple(j_layer), + zmin=zmin, + zmax=ztop, + thickness=ztop - zmin, + material=j_name, + layer_type="dielectric", + mesh_resolution=mesh_resolution, + ) + # Depleted silicon has no free carriers: pure real permittivity. + materials[j_name] = MaterialProperties( + permittivity=cfg.permittivity, + dispersion_models=[], + ) + + result["junction"] = { + **cfg.to_metadata(), + "c_f": cfg.capacitance(length, ztop - zmin), + "mode": mode, + "selection_reason": reason, + } + return result + + +__all__ = ["make_doping_profile", "make_pn_junction_profile"] diff --git a/src/gsim/common/stack/junction.py b/src/gsim/common/stack/junction.py new file mode 100644 index 00000000..3e796aca --- /dev/null +++ b/src/gsim/common/stack/junction.py @@ -0,0 +1,439 @@ +"""PN-junction depletion model (Sze, *Physics of Semiconductor Devices*). + +This module implements the textbook depletion approximation for an abrupt or +linearly graded PN junction: + +- S. M. Sze and K. K. Ng, *Physics of Semiconductor Devices*, 3rd ed., + Wiley (2007), chapter 2 ("p-n Junction Diodes"). + +Provided quantities (all concentrations in ``cm^-3``, lengths in ``um``): + +1. Built-in potential:: + + V_bi = (k_B T / q) ln(Na Nd / ni^2) (Sze eq. 2.60) + +2. Depletion width under reverse bias VR (abrupt junction):: + + W = sqrt( 2 eps_s (V_bi + VR) / q * (Na + Nd)/(Na Nd) ) (eq. 2.66) + x_p = W Nd / (Na + Nd) (spilled into the P side) + x_n = W Na / (Na + Nd) (spilled into the N side) + +3. Depletion width for a linearly graded junction with grade constant + ``a = |dN/dx|`` near the metallurgical junction:: + + W = [ 12 eps_s (V_bi + VR) / (q a) ]^(1/3) (eq. 2.72) + +4. Junction capacitance per unit area (parallel-plate form of the depletion + charge, valid for W much smaller than the device lateral dimensions):: + + C_j = eps_s / W + +The same module also provides :func:`select_junction_mode`, which decides +whether the depletion strip can be resolved on the simulation mesh +(``"high_res"``) or should be collapsed into a lumped capacitance boundary +(``"capacitance"``). + +Example: +------- + >>> from gsim.common.stack.junction import PNJunctionConfig + >>> junc = PNJunctionConfig(na_cm3=1e19, nd_cm3=1e19, v_reverse=0.0) + >>> junc.v_bi # built-in potential [V] + >>> junc.w_um # total depletion width [um] + >>> junc.xp_um # depletion extent into the P side [um] + >>> junc.xn_um # depletion extent into the N side [um] + >>> junc.capacitance(length_um=10.0, height_um=0.22) # absolute C [F] +""" + +from __future__ import annotations + +import math +from typing import Any, Literal, Self + +from pydantic import BaseModel, ConfigDict, Field, model_validator +from scipy.constants import Boltzmann as KB # noqa: N814 +from scipy.constants import elementary_charge as Q # noqa: N812 +from scipy.constants import epsilon_0 as EPS0 # noqa: N812 + +__all__ = [ + "DEFAULT_SI_PERMITTIVITY", + "JUNCTION_MODE_FRACTION", + "NI_SI_300K_CM3", + "PNJunctionConfig", + "built_in_voltage", + "depletion_extents", + "depletion_width", + "junction_capacitance_per_area", + "select_junction_mode", +] + +#: Intrinsic carrier concentration of silicon at 300 K in cm^-3. +#: Classic textbook value used by Sze; override for other materials/T. +NI_SI_300K_CM3: float = 1.5e10 + +#: Default relative permittivity of depleted (intrinsic) silicon. +DEFAULT_SI_PERMITTIVITY: float = 11.9 + +#: A depletion width is considered mesh-resolvable when it reaches this +#: fraction of the smallest doped section flanking the junction. +JUNCTION_MODE_FRACTION: float = 0.2 + +JunctionMode = Literal["capacitance", "high_res"] + + +def built_in_voltage( + na_cm3: float, + nd_cm3: float, + *, + temperature_k: float = 300.0, + ni_cm3: float = NI_SI_300K_CM3, +) -> float: + """Compute the built-in potential ``V_bi`` of a PN junction in volts. + + Implements ``V_bi = (k_B T / q) ln(Na Nd / ni^2)`` (Sze ch. 2). + + Args: + na_cm3: Acceptor concentration on the P side in cm^-3 (> 0). + nd_cm3: Donor concentration on the N side in cm^-3 (> 0). + temperature_k: Lattice temperature in kelvin (> 0). + ni_cm3: Intrinsic carrier concentration in cm^-3 (> 0). + + Returns: + Built-in potential in volts. + + Raises: + ValueError: If any input is non-positive or ``Na*Nd <= ni**2``. + """ + if na_cm3 <= 0 or nd_cm3 <= 0: + raise ValueError("Doping concentrations must be positive (cm^-3).") + if temperature_k <= 0: + raise ValueError("temperature_k must be positive.") + if ni_cm3 <= 0: + raise ValueError("ni_cm3 must be positive.") + product = na_cm3 * nd_cm3 + if product <= ni_cm3**2: + raise ValueError( + f"Na*Nd ({product:.3g} cm^-6) must exceed ni^2 " + f"({ni_cm3**2:.3g} cm^-6); degenerate case has no junction." + ) + vt = KB * temperature_k / Q + return float(vt * math.log(product / ni_cm3**2)) + + +def _validate_bias(v_reverse: float, v_bi: float) -> None: + """Reject bias points beyond flat-band (no physical solution).""" + if v_bi + v_reverse <= 0: + raise ValueError( + f"V_bi + v_reverse = {v_bi + v_reverse:.4g} V must be > 0 " + "(applied forward bias beyond flat-band has no solution)." + ) + + +def _eps_si(permittivity: float) -> float: + """Return absolute permittivity in F/m from a relative value.""" + if permittivity < 1.0: + raise ValueError("permittivity must be >= 1.") + return permittivity * EPS0 + + +def depletion_width( + na_cm3: float, + nd_cm3: float, + *, + v_reverse: float = 0.0, + temperature_k: float = 300.0, + ni_cm3: float = NI_SI_300K_CM3, + permittivity: float = DEFAULT_SI_PERMITTIVITY, + grading: Literal["abrupt", "linear"] = "abrupt", + grade_const_cm4: float | None = None, +) -> float: + """Compute the total depletion width ``W`` in micrometers. + + Args: + na_cm3: Acceptor concentration in cm^-3 (> 0). + nd_cm3: Donor concentration in cm^-3 (> 0). + v_reverse: Applied reverse-bias voltage in volts (positive = reverse). + Negative values model forward bias down to (but excluding) + flat-band. + temperature_k: Lattice temperature in kelvin. + ni_cm3: Intrinsic carrier concentration in cm^-3. + permittivity: Relative permittivity of the semiconductor. + grading: ``"abrupt"`` (step junction) or ``"linear"`` (linearly + graded). + grade_const_cm4: Grade constant ``a = |dN/dx|`` in cm^-4 for + ``grading="linear"``. + + Returns: + Total depletion width in micrometers. + + Raises: + ValueError: On non-positive inputs, missing grade constant, or bias + beyond flat-band. + """ + v_bi = built_in_voltage(na_cm3, nd_cm3, temperature_k=temperature_k, ni_cm3=ni_cm3) + _validate_bias(v_reverse, v_bi) + eps_s = _eps_si(permittivity) + + if grading == "linear": + if grade_const_cm4 is None or grade_const_cm4 <= 0: + raise ValueError("grading='linear' requires grade_const_cm4 > 0.") + # a in m^-4 (1 cm^-4 = 1e6 m^-4); W comes out in meters. + a_m4 = grade_const_cm4 * 1e6 + w_m = (12.0 * eps_s * (v_bi + v_reverse) / (Q * a_m4)) ** (1.0 / 3.0) + return float(w_m * 1e6) + + if grading != "abrupt": + raise ValueError(f"Unknown grading type: {grading!r}") + + # Abrupt junction: W = sqrt(2 eps_s (V_bi+VR)/q * (Na+Nd)/(NaNd)). + na_m3 = na_cm3 * 1e6 + nd_m3 = nd_cm3 * 1e6 + w_m = math.sqrt( + 2.0 * eps_s * (v_bi + v_reverse) / Q * (na_m3 + nd_m3) / (na_m3 * nd_m3) + ) + return float(w_m * 1e6) + + +def depletion_extents( + na_cm3: float, + nd_cm3: float, + *, + w_um: float, + grading: Literal["abrupt", "linear"] = "abrupt", +) -> tuple[float, float]: + """Split a total depletion width into P-side/N-side extents in micrometers. + + For an abrupt junction the depletion spills asymmetrically:: + + x_p = W Nd / (Na + Nd), x_n = W Na / (Na + Nd) + + A linearly graded junction is symmetric around the metallurgical + junction, so ``x_p = x_n = W/2``. + + Args: + na_cm3: Acceptor concentration in cm^-3 (> 0). + nd_cm3: Donor concentration in cm^-3 (> 0). + w_um: Total depletion width in micrometers (from + :func:`depletion_width`). + grading: Junction grading type. + + Returns: + ``(xp_um, xn_um)`` — extents spilled into the P and N sides. + """ + if na_cm3 <= 0 or nd_cm3 <= 0: + raise ValueError("Doping concentrations must be positive (cm^-3).") + if w_um < 0: + raise ValueError("w_um must be non-negative.") + if grading == "linear": + return w_um / 2.0, w_um / 2.0 + total = na_cm3 + nd_cm3 + return w_um * nd_cm3 / total, w_um * na_cm3 / total + + +def junction_capacitance_per_area( + permittivity: float, + w_um: float, +) -> float: + """Depletion capacitance per unit area ``C_j = eps_s / W`` in F/m^2. + + Args: + permittivity: Relative permittivity of the semiconductor. + w_um: Total depletion width in micrometers (> 0). + + Returns: + Capacitance per unit area in F/m^2. + """ + if w_um <= 0: + raise ValueError("w_um must be positive.") + return _eps_si(permittivity) / (w_um * 1e-6) + + +def select_junction_mode( + w_um: float, + p_extent_um: float, + n_extent_um: float, + *, + fraction: float = JUNCTION_MODE_FRACTION, +) -> JunctionMode: + """Choose how to represent the depletion region in a simulation. + + The depletion strip is meshed explicitly (``"high_res"``) when its width + is comparable to the doped sections flanking it — specifically when + ``w_um >= fraction * min(p_extent, n_extent)``. Otherwise the region is + far thinner than its neighbours and meshing it would only bloat the + model, so a lumped capacitance boundary is used instead + (``"capacitance"``). + + Args: + w_um: Total depletion width in micrometers (> 0). + p_extent_um: Size of the doped section flanking the junction on the + P side (micrometers, > 0). + n_extent_um: Size of the doped section flanking the junction on the + N side (micrometers, > 0). + fraction: Resolvability threshold as a fraction of the smaller flank + (default ~1/5). + + Returns: + ``"high_res"`` when the geometry should carry the depletion strip, + ``"capacitance"`` otherwise. + """ + if w_um <= 0: + raise ValueError("w_um must be positive.") + if p_extent_um <= 0 or n_extent_um <= 0: + raise ValueError("Flank extents must be positive.") + if not 0 < fraction <= 1: + raise ValueError("fraction must lie in (0, 1].") + threshold_um = fraction * min(p_extent_um, n_extent_um) + return "high_res" if w_um >= threshold_um else "capacitance" + + +class PNJunctionConfig(BaseModel): + """Parameters of a PN-junction depletion model (depletion approximation). + + Concentrations use the semiconductor-industry convention (cm^-3); + derived lengths are exposed in micrometers and capacitances in farads. + See module docstring for the underlying formulas (Sze ch. 2). + + Attributes: + na_cm3: Acceptor concentration on the P side (cm^-3). + nd_cm3: Donor concentration on the N side (cm^-3). + v_reverse: Applied reverse bias in volts (positive = reverse; + negative values model forward bias below flat-band). + temperature_k: Lattice temperature in kelvin. + ni_cm3: Intrinsic carrier concentration (cm^-3). + permittivity: Relative permittivity of the depleted semiconductor. + grading: ``"abrupt"`` or ``"linear"`` junction profile. + grade_const_cm4: Grade constant ``a = |dN/dx|`` in cm^-4, required + when ``grading="linear"``. + """ + + model_config = ConfigDict(validate_assignment=True) + + na_cm3: float = Field(gt=0, description="Acceptor concentration (cm^-3)") + nd_cm3: float = Field(gt=0, description="Donor concentration (cm^-3)") + v_reverse: float = Field( + default=0.0, description="Applied reverse bias [V] (positive = reverse)" + ) + temperature_k: float = Field(default=300.0, gt=0, description="Temperature [K]") + ni_cm3: float = Field( + default=NI_SI_300K_CM3, gt=0, description="Intrinsic carriers (cm^-3)" + ) + permittivity: float = Field( + default=DEFAULT_SI_PERMITTIVITY, + ge=1.0, + description="Relative permittivity of the semiconductor", + ) + grading: Literal["abrupt", "linear"] = Field(default="abrupt") + grade_const_cm4: float | None = Field( + default=None, gt=0, description="Grade constant a = |dN/dx| (cm^-4)" + ) + + @model_validator(mode="after") + def _validate_physics(self) -> Self: + """Check grading configuration and bias range.""" + if self.grading == "linear" and self.grade_const_cm4 is None: + raise ValueError("grading='linear' requires grade_const_cm4.") + _validate_bias(self.v_reverse, self.v_bi) + return self + + @property + def v_bi(self) -> float: + """Built-in potential in volts.""" + return built_in_voltage( + self.na_cm3, + self.nd_cm3, + temperature_k=self.temperature_k, + ni_cm3=self.ni_cm3, + ) + + @property + def w_um(self) -> float: + """Total depletion width in micrometers at the configured bias.""" + return depletion_width( + self.na_cm3, + self.nd_cm3, + v_reverse=self.v_reverse, + temperature_k=self.temperature_k, + ni_cm3=self.ni_cm3, + permittivity=self.permittivity, + grading=self.grading, + grade_const_cm4=self.grade_const_cm4, + ) + + @property + def xp_um(self) -> float: + """Depletion extent spilled into the P side (micrometers).""" + xp, _xn = depletion_extents( + self.na_cm3, self.nd_cm3, w_um=self.w_um, grading=self.grading + ) + return xp + + @property + def xn_um(self) -> float: + """Depletion extent spilled into the N side (micrometers).""" + _xp, xn = depletion_extents( + self.na_cm3, self.nd_cm3, w_um=self.w_um, grading=self.grading + ) + return xn + + @property + def c_per_area(self) -> float: + """Junction capacitance per unit area in F/m^2 (``eps_s / W``).""" + return junction_capacitance_per_area(self.permittivity, self.w_um) + + def capacitance(self, length_um: float, height_um: float) -> float: + """Absolute junction capacitance for a rectangular junction face. + + Treats the depletion strip as a parallel-plate capacitor of area + ``length x height`` filled with the depleted semiconductor: + ``C = eps_s * A / W``. + + Args: + length_um: Device length along the propagation direction (um). + height_um: Junction z-extent (um), e.g. the rib height. + + Returns: + Absolute capacitance in farads. + """ + if length_um <= 0 or height_um <= 0: + raise ValueError("length_um and height_um must be positive.") + area_m2 = length_um * height_um * 1e-12 + return float(self.c_per_area * area_m2) + + def select_mode( + self, + p_extent_um: float, + n_extent_um: float, + *, + fraction: float = JUNCTION_MODE_FRACTION, + ) -> JunctionMode: + """Auto-select the representation mode for this junction. + + Thin wrapper around :func:`select_junction_mode` using this config's + computed depletion width. + + Args: + p_extent_um: Size of the doped flank on the P side (um). + n_extent_um: Size of the doped flank on the N side (um). + fraction: Resolvability threshold fraction (~1/5 default). + + Returns: + ``"high_res"`` or ``"capacitance"``. + """ + return select_junction_mode( + self.w_um, p_extent_um, n_extent_um, fraction=fraction + ) + + def to_metadata(self) -> dict[str, Any]: + """Return a plain-dict summary of the computed junction quantities.""" + return { + "na_cm3": self.na_cm3, + "nd_cm3": self.nd_cm3, + "v_reverse": self.v_reverse, + "temperature_k": self.temperature_k, + "v_bi": self.v_bi, + "w_um": self.w_um, + "xp_um": self.xp_um, + "xn_um": self.xn_um, + "c_per_area_f_m2": self.c_per_area, + "grading": self.grading, + } diff --git a/src/gsim/palace/base.py b/src/gsim/palace/base.py index b7151138..f1e0d52a 100644 --- a/src/gsim/palace/base.py +++ b/src/gsim/palace/base.py @@ -427,6 +427,68 @@ def add_impedance_boundary( ) ) + def set_pn_junction( + self, + junction: Any, + *, + layer_p: str, + layer_n: str, + length_um: float, + height_um: float, + name: str | None = None, + ) -> float: + """Apply the depletion capacitance of a PN junction between two layers. + + Capacitance-mode modelling of a PN junction: the depletion width + ``W`` is computed from the doping concentrations and bias point via + :class:`gsim.common.stack.junction.PNJunctionConfig` (Sze, + *Physics of Semiconductor Devices*, ch. 2), converted to an absolute + parallel-plate capacitance ``C = eps_s * A / W``, and applied as a + lumped Impedance boundary on the shared P/N interface. + + Use this when the depletion strip is too thin to resolve on the mesh + (the auto-selection in + :func:`gsim.common.stack.doping.make_pn_junction_profile` picks this + regime); for well-resolved depletion regions prefer drawing them as + dielectric geometry (``mode="high_res"``) instead. + + Args: + junction: ``PNJunctionConfig`` or its dict form (doping + concentrations, bias, temperature, permittivity). + layer_p: Name of the P-doped layer. + layer_n: Name of the N-doped layer. + length_um: Device length along the propagation direction (um). + height_um: Junction z-extent (um), e.g. the rib height. + name: Optional display name for the boundary. + + Returns: + The absolute capacitance applied [F]. + + Example: + >>> sim.set_pn_junction( + ... {"na_cm3": 1e19, "nd_cm3": 1e19}, + ... layer_p="p_rib", + ... layer_n="n_rib", + ... length_um=10.0, + ... height_um=0.22, + ... ) + """ + from gsim.common.stack.junction import PNJunctionConfig + + cfg = ( + junction + if isinstance(junction, PNJunctionConfig) + else PNJunctionConfig.model_validate(junction) + ) + capacitance = cfg.capacitance(length_um=length_um, height_um=height_um) + self.add_impedance_boundary( + layer_p, + layer_n, + capacitance=capacitance, + name=name, + ) + return capacitance + # ------------------------------------------------------------------------- # Material methods # ------------------------------------------------------------------------- diff --git a/tests/common/test_cross_section.py b/tests/common/test_cross_section.py index cfd39260..2bc04e56 100644 --- a/tests/common/test_cross_section.py +++ b/tests/common/test_cross_section.py @@ -387,6 +387,44 @@ def test_builds_stack_with_doping_and_rib_layers(self): layers = {r.layer_name for r in section} assert {"core", "pp_slab_0", "npp_slab_0"} <= layers + def test_doping_materials_registered_on_stack(self): + """Doping/rib materials must land on stack.materials for solver config. + + Regression: the merged materials dict used to be computed but never + attached, so doped domains silently resolved to eps=1.0 without + conductivity in the generated Palace config. + """ + comp, LAYER = self._build_component() + doping = self._doping_result(comp) + + stack, _section = build_doped_cross_section( + comp, + axis="x", + value=0.0, + substrate_thickness=2.0, + include_substrate=False, + doping=doping, + metal1=(1.1, 1.0), + rib_layers=[ + ("p_rib", LAYER.P, 1.6e3), + ("n_rib", LAYER.N, 1.6e3), + ], + permittivity=11.9, + fmax=200e9, + verbose=False, + ) + + for name, sigma in ( + ("pp_slab_0", 2e4), + ("p_rib", 1.6e3), + ("n_rib", 1.6e3), + ): + assert name in stack.materials, f"{name} missing from stack.materials" + props = stack.materials[name] + assert isinstance(props, dict) + assert props["permittivity"] == pytest.approx(11.9) + assert props["conductivity"] == pytest.approx(sigma) + def test_metal1_override_applied(self): comp, _LAYER = self._build_component() stack, _ = build_doped_cross_section( diff --git a/tests/common/test_junction_physics.py b/tests/common/test_junction_physics.py new file mode 100644 index 00000000..1b830b31 --- /dev/null +++ b/tests/common/test_junction_physics.py @@ -0,0 +1,215 @@ +"""Tests for the PN-junction depletion physics (Sze ch. 2 formulas). + +The expected values are recomputed here from the textbook expressions with +scipy.constants so the tests validate the wiring independently of the +implementation internals. +""" + +from __future__ import annotations + +import math + +import pytest +from pydantic import ValidationError +from scipy.constants import Boltzmann as KB # noqa: N814 +from scipy.constants import elementary_charge as Q # noqa: N812 +from scipy.constants import epsilon_0 as EPS0 # noqa: N812 + +from gsim.common.stack.junction import ( + PNJunctionConfig, + built_in_voltage, + depletion_extents, + depletion_width, + junction_capacitance_per_area, + select_junction_mode, +) + +VT_300 = KB * 300.0 / Q + + +class TestBuiltInVoltage: + def test_symmetric_silicon_value(self): + v_bi = built_in_voltage(1e19, 1e19) + expected = VT_300 * math.log(1e38 / (1.5e10) ** 2) + assert v_bi == pytest.approx(expected, rel=1e-12) + assert v_bi == pytest.approx(1.05, abs=0.03) + + def test_temperature_dependence(self): + cold = built_in_voltage(1e18, 1e18, temperature_k=250.0) + hot = built_in_voltage(1e18, 1e18, temperature_k=350.0) + expected_cold = KB * 250.0 / Q * math.log(1e36 / (1.5e10) ** 2) + expected_hot = KB * 350.0 / Q * math.log(1e36 / (1.5e10) ** 2) + assert cold == pytest.approx(expected_cold, rel=1e-12) + assert hot == pytest.approx(expected_hot, rel=1e-12) + + def test_rejects_nonphysical_inputs(self): + with pytest.raises(ValueError): + built_in_voltage(-1e18, 1e18) + with pytest.raises(ValueError): + built_in_voltage(1e18, 0.0) + with pytest.raises(ValueError): + built_in_voltage(1e18, 1e18, temperature_k=0.0) + + def test_rejects_degenerate_doping(self): + with pytest.raises(ValueError, match="ni"): + built_in_voltage(1e9, 1e9) + + +class TestDepletionWidthAbrupt: + def test_symmetric_hand_check(self): + w_um = depletion_width(1e18, 1e18) + na_m3 = nd_m3 = 1e18 * 1e6 + eps_s = 11.9 * EPS0 + expected_m = math.sqrt( + 2 + * eps_s + * VT_300 + * math.log(1e36 / (1.5e10) ** 2) + / Q + * (na_m3 + nd_m3) + / (na_m3 * nd_m3) + ) + assert w_um == pytest.approx(expected_m * 1e6, rel=1e-12) + + def test_reverse_bias_sqrt_scaling(self): + w0 = depletion_width(1e19, 5e17) + vbi = built_in_voltage(1e19, 5e17) + w_r = depletion_width(1e19, 5e17, v_reverse=2.0) + assert w_r / w0 == pytest.approx(math.sqrt((vbi + 2.0) / vbi), rel=1e-12) + + def test_one_sided_limit(self): + # NA >> ND: nearly all the depletion spills into the lightly doped side. + w = depletion_width(1e20, 1e17) + xp, xn = depletion_extents(1e20, 1e17, w_um=w) + assert xn == pytest.approx(w, rel=1e-3) + assert xp == pytest.approx(w * 1e-3, rel=1e-2) + + def test_forward_bias_below_flatband(self): + vbi = built_in_voltage(1e18, 1e18) + w_eq = depletion_width(1e18, 1e18) + w_fw = depletion_width(1e18, 1e18, v_reverse=-vbi / 2) + assert w_fw < w_eq + with pytest.raises(ValueError, match="flat-band"): + depletion_width(1e18, 1e18, v_reverse=-(vbi + 0.01)) + + +class TestDepletionWidthGraded: + def test_cubic_root_law(self): + a_cm4 = 1e21 + vbi = built_in_voltage(1e18, 1e18) + w = depletion_width(1e18, 1e18, grading="linear", grade_const_cm4=a_cm4) + eps_s = 11.9 * EPS0 + expected_m = (12 * eps_s * vbi / (Q * a_cm4 * 1e6)) ** (1 / 3) + assert w == pytest.approx(expected_m * 1e6, rel=1e-12) + + def test_graded_bias_scaling(self): + kwargs = dict(grading="linear", grade_const_cm4=1e21) + w0 = depletion_width(1e18, 1e18, **kwargs) + w_r = depletion_width(1e18, 1e18, v_reverse=1.0, **kwargs) + vbi = built_in_voltage(1e18, 1e18) + assert w_r / w0 == pytest.approx(((vbi + 1.0) / vbi) ** (1 / 3), rel=1e-12) + + def test_graded_is_symmetric(self): + w = depletion_width(1e18, 1e19, grading="linear", grade_const_cm4=1e20) + xp, xn = depletion_extents(1e18, 1e19, w_um=w, grading="linear") + assert xp == pytest.approx(w / 2) + assert xn == pytest.approx(w / 2) + + def test_requires_grade_constant(self): + with pytest.raises(ValueError, match="grade_const"): + depletion_width(1e18, 1e18, grading="linear") + + def test_unknown_grading(self): + with pytest.raises(ValueError, match="grading"): + depletion_width(1e18, 1e18, grading="exponential") # type: ignore[arg-type] + + +class TestCapacitance: + def test_per_area_inverse_w(self): + eps_r = 11.9 + for w_um in (0.01, 0.05, 0.2): + c = junction_capacitance_per_area(eps_r, w_um) + assert c == pytest.approx(eps_r * EPS0 / (w_um * 1e-6), rel=1e-12) + + def test_absolute_capacitance(self): + junc = PNJunctionConfig(na_cm3=1e19, nd_cm3=1e19) + c = junc.capacitance(length_um=10.0, height_um=0.22) + area_m2 = 10.0 * 0.22 * 1e-12 + assert c == pytest.approx(junc.c_per_area * area_m2, rel=1e-12) + # Same order as typical TW-MZM junction caps (~fF per 10 um). + assert 1e-15 < c < 1e-13 + + def test_capacitance_scales_with_bias(self): + junc0 = PNJunctionConfig(na_cm3=1e19, nd_cm3=1e19) + junc_r = PNJunctionConfig(na_cm3=1e19, nd_cm3=1e19, v_reverse=3.0) + assert junc_r.capacitance(10.0, 0.22) < junc0.capacitance(10.0, 0.22) + + +class TestSelectJunctionMode: + def test_comparable_width_selects_high_res(self): + assert select_junction_mode(0.05, 0.2, 0.2) == "high_res" + assert select_junction_mode(0.0401, 0.2, 0.2) == "high_res" + + def test_too_thin_selects_capacitance(self): + assert select_junction_mode(0.0166, 0.2, 0.2) == "capacitance" + assert select_junction_mode(0.0399, 0.2, 0.2) == "capacitance" + + def test_threshold_is_fraction_of_smaller_flank(self): + assert select_junction_mode(0.0099, 0.05, 0.4, fraction=0.2) == "capacitance" + assert select_junction_mode(0.0101, 0.05, 0.4, fraction=0.2) == "high_res" + + def test_custom_fraction(self): + assert select_junction_mode(0.09, 0.2, 0.2, fraction=0.5) == "capacitance" + assert select_junction_mode(0.11, 0.2, 0.2, fraction=0.5) == "high_res" + + def test_invalid_inputs(self): + with pytest.raises(ValueError): + select_junction_mode(0.0, 0.2, 0.2) + with pytest.raises(ValueError): + select_junction_mode(0.1, 0.0, 0.2) + with pytest.raises(ValueError): + select_junction_mode(0.1, 0.2, 0.2, fraction=1.5) + + +class TestPNJunctionConfig: + def test_derived_quantities_consistent(self): + cfg = PNJunctionConfig(na_cm3=2e18, nd_cm3=8e18, v_reverse=0.5) + assert cfg.v_bi == pytest.approx(built_in_voltage(2e18, 8e18)) + assert cfg.w_um == pytest.approx( + depletion_width(2e18, 8e18, v_reverse=0.5), rel=1e-12 + ) + total = cfg.xp_um + cfg.xn_um + assert total == pytest.approx(cfg.w_um, rel=1e-12) + # Asymmetric split: more depletion on the lighter-doped side. + assert cfg.xp_um > cfg.xn_um + + def test_dict_construction(self): + cfg = PNJunctionConfig.model_validate({"na_cm3": 1e19, "nd_cm3": 1e19}) + assert cfg.na_cm3 == 1e19 + + def test_linear_requires_grade_const(self): + with pytest.raises(ValidationError, match="grade_const"): + PNJunctionConfig(na_cm3=1e18, nd_cm3=1e18, grading="linear") + + def test_rejects_beyond_flatband(self): + vbi = built_in_voltage(1e18, 1e18) + with pytest.raises(ValidationError): + PNJunctionConfig(na_cm3=1e18, nd_cm3=1e18, v_reverse=-vbi - 0.05) + + def test_rejects_bad_concentrations(self): + with pytest.raises(ValidationError): + PNJunctionConfig(na_cm3=0.0, nd_cm3=1e18) + + def test_to_metadata_keys(self): + meta = PNJunctionConfig(na_cm3=1e18, nd_cm3=1e18).to_metadata() + for key in ( + "na_cm3", + "nd_cm3", + "v_bi", + "w_um", + "xp_um", + "xn_um", + "c_per_area_f_m2", + "grading", + ): + assert key in meta diff --git a/tests/common/test_junction_profile.py b/tests/common/test_junction_profile.py new file mode 100644 index 00000000..29061113 --- /dev/null +++ b/tests/common/test_junction_profile.py @@ -0,0 +1,180 @@ +"""Tests for ``make_pn_junction_profile`` geometry, materials and mode selection.""" + +from __future__ import annotations + +import logging + +import gdsfactory as gf +import pytest + +from gsim.common.cross_section import extract_plane_section +from gsim.common.stack.doping import make_pn_junction_profile +from gsim.common.stack.extractor import LayerStack +from gsim.common.stack.junction import PNJunctionConfig + +CY = -20.0 +RIB_WIDTH = 0.4 +LENGTH = 10.0 + +P_REGION = ("p_rib", (21, 0), 1.6e3) +N_REGION = ("n_rib", (20, 0), 1.6e3) +JUNCTION_REGION = ("junction", (22, 0)) + + +def _thin_junction() -> PNJunctionConfig: + """Na=Nd=1e19 cm^-3 at zero bias -> W ~ 17 nm < threshold.""" + return PNJunctionConfig(na_cm3=1e19, nd_cm3=1e19) + + +def _wide_junction() -> PNJunctionConfig: + """Light doping + reverse bias -> W ~ 71 nm > threshold (40 nm).""" + return PNJunctionConfig(na_cm3=1e18, nd_cm3=1e18, v_reverse=1.0) + + +def _build(junction, **kwargs): + comp = gf.Component() + kwargs.setdefault("p_region", P_REGION) + kwargs.setdefault("n_region", N_REGION) + kwargs.setdefault("zmin", 0.0) + kwargs.setdefault("zmax", 0.22) + result = make_pn_junction_profile( + comp, + length=LENGTH, + center_y=CY, + rib_width=RIB_WIDTH, + junction=junction, + **kwargs, + ) + return comp, result + + +def _section_rects(comp, result): + """Extract the x=0 plane section from a profile-built component.""" + stack = LayerStack(pdk_name="test") + stack.layers.update(result["layer_specs"]) + for name, mat in result["materials"].items(): + stack.materials[name] = mat.to_dict() + rects = extract_plane_section(comp.copy(), stack, axis="x", value=0.0) + return sorted(rects, key=lambda r: r.y0) + + +class TestAutoModeSelection: + def test_thin_junction_selects_capacitance(self): + _comp, res = _build(_thin_junction()) + assert res["junction"]["mode"] == "capacitance" + assert "threshold" in res["junction"]["selection_reason"] + + def test_wide_junction_selects_high_res(self): + _comp, res = _build(_wide_junction(), junction_region=JUNCTION_REGION) + assert res["junction"]["mode"] == "high_res" + + def test_auto_logs_selection_reason(self, caplog): + with caplog.at_level(logging.INFO, logger="gsim.common.stack.doping"): + _comp, _res = _build(_thin_junction()) + assert any("capacitance" in rec.message for rec in caplog.records) + + def test_forced_mode_overrides_auto(self): + _comp, res = _build( + _thin_junction(), mode="high_res", junction_region=JUNCTION_REGION + ) + assert res["junction"]["mode"] == "high_res" + assert "forced" in res["junction"]["selection_reason"] + _comp, res = _build(_wide_junction(), mode="capacitance") + assert res["junction"]["mode"] == "capacitance" + + +class TestCapacitanceModeGeometry: + def test_no_junction_polygon_or_spec(self): + comp, res = _build(_thin_junction()) + assert "junction" not in res["layer_specs"] + assert "junction" not in res["materials"] + # No polygon may exist on the junction GDS layer. + polys = comp.get_polygons(layers=(JUNCTION_REGION[1],)) + assert not any(v for v in polys.values()) + + def test_p_n_adjacent_halves(self): + comp, res = _build(_thin_junction()) + rects = _section_rects(comp, res) + names = [r.layer_name for r in rects] + assert set(names) == {"p_rib", "n_rib"} + by_name = {r.layer_name: r for r in rects} + assert by_name["p_rib"].y0 == pytest.approx(CY) + assert by_name["n_rib"].y1 == pytest.approx(CY) + + def test_junction_metadata_present(self): + junc = _thin_junction() + _comp, res = _build(junc) + meta = res["junction"] + assert meta["w_um"] == pytest.approx(junc.w_um) + assert meta["c_f"] == pytest.approx(junc.capacitance(LENGTH, 0.22)) + assert meta["xp_um"] + meta["xn_um"] == pytest.approx(meta["w_um"]) + + +class TestHighResModeGeometry: + def test_three_contiguous_regions(self): + comp, res = _build(_wide_junction(), junction_region=JUNCTION_REGION) + rects = _section_rects(comp, res) + names = [r.layer_name for r in rects] + assert names == ["n_rib", "junction", "p_rib"] + + n_r, j_r, p_r = rects + # Contiguity with no gaps or overlaps. + assert n_r.y1 == pytest.approx(j_r.y0) + assert j_r.y1 == pytest.approx(p_r.y0) + + junc = _wide_junction() + # Depletion strip spans [cy - xn, cy + xp] (within layout DBU rounding). + assert j_r.y0 == pytest.approx(CY - junc.xn_um, abs=2e-3) + assert j_r.y1 == pytest.approx(CY + junc.xp_um, abs=2e-3) + assert (j_r.y1 - j_r.y0) == pytest.approx(junc.w_um, abs=4e-3) + # Flanks fill the rest of the rib. + assert (p_r.y1 - p_r.y0) == pytest.approx(RIB_WIDTH / 2 - junc.xp_um, abs=4e-3) + assert (n_r.y1 - n_r.y0) == pytest.approx(RIB_WIDTH / 2 - junc.xn_um, abs=4e-3) + # Full rib span is covered exactly once. + assert p_r.y1 - n_r.y0 == pytest.approx(RIB_WIDTH) + + def test_material_models(self): + _comp, res = _build(_wide_junction(), junction_region=JUNCTION_REGION) + # Doped regions carry Drude conductivity. + for name in ("p_rib", "n_rib"): + mat = res["materials"][name] + assert mat.conductivity == pytest.approx(1.6e3) + assert mat.permittivity == pytest.approx(11.9) + # Junction strip: depleted silicon -> pure real permittivity, no carriers. + jmat = res["materials"]["junction"] + assert jmat.permittivity == pytest.approx(11.9) + assert jmat.conductivity is None + assert jmat.dispersion_models == [] + + def test_high_res_requires_junction_region(self): + with pytest.raises(ValueError, match="junction_region"): + _build(_wide_junction()) + + def test_layer_specs_reference_materials(self): + _comp, res = _build(_wide_junction(), junction_region=JUNCTION_REGION) + for name in ("p_rib", "n_rib", "junction"): + spec = res["layer_specs"][name] + assert spec.material == name + assert spec.zmin == 0.0 + assert spec.zmax == 0.22 + + +class TestValidation: + def test_depletion_wider_than_rib_rejected(self): + big = PNJunctionConfig(na_cm3=1e16, nd_cm3=1e16, v_reverse=5.0) + if big.w_um <= RIB_WIDTH: + pytest.skip("picked parameters do not exceed the rib width") + with pytest.raises(ValueError, match="fit"): + _build(big) + + def test_invalid_zmax_rejected(self): + with pytest.raises(ValueError): + _build(_thin_junction(), zmax=-1.0) + + def test_accepts_dict_junction_config(self): + _comp, res = _build( + {"na_cm3": 1e18, "nd_cm3": 1e18, "v_reverse": 1.0}, + junction_region=JUNCTION_REGION, + ) + assert res["junction"]["w_um"] == pytest.approx(_wide_junction().w_um) + assert res["junction"]["mode"] == "high_res" diff --git a/tests/palace/test_pn_junction_modes.py b/tests/palace/test_pn_junction_modes.py new file mode 100644 index 00000000..6c691946 --- /dev/null +++ b/tests/palace/test_pn_junction_modes.py @@ -0,0 +1,159 @@ +"""End-to-end tests: PN-junction capacitance vs high-res mesh representation. + +Capacitance mode must produce a Palace ``Boundaries.Impedance`` entry with +``Cs = C / interface_length`` and no junction domain. High-res mode must +produce a ``junction`` dielectric domain (pure real permittivity) and no +Impedance boundary. +""" + +from __future__ import annotations + +import json +from pathlib import Path + +import gdsfactory as gf +import pytest + +from gsim.common.cross_section import build_doped_cross_section +from gsim.common.stack.doping import make_pn_junction_profile +from gsim.common.stack.junction import PNJunctionConfig +from gsim.palace import BoundaryModeSim + +F_RF = 50e9 + + +def _thin_junction() -> PNJunctionConfig: + """W ~ 17 nm -> below the auto threshold -> capacitance mode.""" + return PNJunctionConfig(na_cm3=1e19, nd_cm3=1e19) + + +def _wide_junction() -> PNJunctionConfig: + """W ~ 71 nm -> above the auto threshold -> high-res mode.""" + return PNJunctionConfig(na_cm3=1e18, nd_cm3=1e18, v_reverse=1.0) + + +def _build_device(junction: PNJunctionConfig, **profile_kwargs): + """Build the rib+slab+doping device and return (comp, stack).""" + gf.gpdk.PDK.activate() + comp = gf.Component() + wg = comp << gf.c.rectangle((10.0, 0.4), centered=True, layer=(1, 0)) + wg.y = -20.0 + slab = comp << gf.c.rectangle((10.0, 100.0), centered=True, layer=(3, 0)) + slab.y = -5.0 + + pn = make_pn_junction_profile( + comp, + length=10.0, + center_y=-20.0, + rib_width=0.4, + junction=junction, + p_region=("p_rib", (21, 0), 1.6e3), + n_region=("n_rib", (20, 0), 1.6e3), + junction_region=("junction", (22, 0)), + zmin=0.0, + zmax=0.22, + **profile_kwargs, + ) + stack, _section = build_doped_cross_section( + comp, + axis="x", + value=0.0, + substrate_thickness=2.0, + doping=pn, + verbose=False, + ) + return comp, stack, pn + + +def _make_sim(junction: PNJunctionConfig, tmp_path: Path, apply_capacitance: bool): + comp, stack, pn = _build_device(junction) + sim = BoundaryModeSim() + sim.set_output_dir(str(tmp_path / "palace-sim-pn")) + sim.set_stack(stack) + sim.set_airbox(margin_x=3.0, margin_y=3.0, z_above=2.0, z_below=2.0) + sim.set_geometry(comp) + sim.set_cross_section("x=0") + sim.set_boundary_mode(freq=F_RF, num_modes=1, save=0) + sim.mesh(preset="coarse", refined_mesh_size=0.05, max_mesh_size=40.0) + if apply_capacitance: + applied = sim.set_pn_junction( + junction, + layer_p="p_rib", + layer_n="n_rib", + length_um=10.0, + height_um=0.22, + ) + assert applied == pytest.approx(junction.capacitance(10.0, 0.22)) + sim.write_config() + config_path = Path(sim.output_dir) / "config.json" + return sim, json.loads(config_path.read_text()), pn + + +@pytest.fixture(scope="module") +def cap_mode(tmp_path_factory): + """Thin depletion: auto-selected capacitance mode with lumped C.""" + return _make_sim(_thin_junction(), tmp_path_factory.mktemp("cap"), True) + + +@pytest.fixture(scope="module") +def hires_mode(tmp_path_factory): + """Wide depletion: auto-selected high-res mode, no lumped C.""" + return _make_sim(_wide_junction(), tmp_path_factory.mktemp("hires"), False) + + +class TestCapacitanceMode: + def test_no_junction_domain_on_mesh(self, cap_mode): + sim, _config, _pn = cap_mode + groups = sim._last_mesh_result.groups + assert "junction" not in groups["volumes"] + + def test_impedance_boundary_in_config(self, cap_mode): + _sim, config, _pn = cap_mode + impedance = config.get("Boundaries", {}).get("Impedance", []) + assert len(impedance) == 1 + assert "Cs" in impedance[0] + assert impedance[0]["Cs"] > 0 + + def test_cs_value_matches_computed_capacitance(self, cap_mode): + _sim, config, pn = cap_mode + # Interface p_rib|n_rib is the vertical rib edge; its curve length is + # the 0.22 um rib height, so Cs = C / 0.22um. + expected_cs = pn["junction"]["c_f"] / (0.22 * 1e-6) + cs = config["Boundaries"]["Impedance"][0]["Cs"] + assert cs == pytest.approx(expected_cs, rel=1e-9) + + def test_doped_domains_present(self, cap_mode): + sim, _config, _pn = cap_mode + groups = sim._last_mesh_result.groups + assert {"p_rib", "n_rib"} <= set(groups["volumes"]) + + +class TestHighResMode: + def test_junction_dielectric_domain_on_mesh(self, hires_mode): + sim, _config, _pn = hires_mode + groups = sim._last_mesh_result.groups + assert "junction" in groups["volumes"] + assert groups["volumes"]["junction"].get("is_shaped_dielectric") is True + + def test_no_impedance_boundary(self, hires_mode): + _sim, config, _pn = hires_mode + assert not config.get("Boundaries", {}).get("Impedance") + + def test_junction_material_is_pure_dielectric(self, hires_mode): + sim, config, _pn = hires_mode + groups = sim._last_mesh_result.groups + junc_attr = groups["volumes"]["junction"]["phys_group"] + materials = config["Domains"]["Materials"] + entries = [m for m in materials if junc_attr in m.get("Attributes", [])] + assert len(entries) == 1, f"Expected one material for attr {junc_attr}" + entry = entries[0] + assert abs(float(entry["Permittivity"]) - 11.9) < 1e-6 + assert not entry.get("Conductivity"), ( + "Depleted silicon must have zero conductivity" + ) + + def test_p_n_junction_strip_contiguous(self, hires_mode): + """All three regions survive as separate domains.""" + sim, _config, _pn = hires_mode + volumes = sim._last_mesh_result.groups["volumes"] + assert {"p_rib", "n_rib", "junction"} <= set(volumes)