diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 09e27191..abf338d8 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -6,14 +6,36 @@ on: - main tags: - "[0-9]+.[0-9]+.[0-9]+" - pull_request: + workflow_dispatch: env: GFP_API_KEY: ${{ secrets.GFP_API_KEY }} jobs: - run-notebooks: + # Discover all notebooks and output as JSON array for the matrix + list-notebooks: runs-on: ubuntu-latest + outputs: + notebooks: ${{ steps.find.outputs.notebooks }} + steps: + - name: Checkout repo + uses: actions/checkout@v4 + with: + ref: ${{ github.event.pull_request.head.sha || github.ref }} + - name: Find notebooks + id: find + run: | + nbs=$(find nbs -name "*.ipynb" -not -path "*/.ipynb_checkpoints/*" | sort | jq -Rsc 'split("\n") | map(select(. != ""))') + echo "notebooks=$nbs" >> "$GITHUB_OUTPUT" + + # Run each notebook in parallel + run-notebook: + runs-on: ubuntu-latest + needs: list-notebooks + strategy: + matrix: + notebook: ${{ fromJson(needs.list-notebooks.outputs.notebooks) }} + fail-fast: false steps: - name: Checkout repo uses: actions/checkout@v4 @@ -29,24 +51,22 @@ jobs: run: just dev - name: Create ipykernel run: uv run --no-sync python -m ipykernel install --user --name gsim - - name: Run notebooks + - name: Run notebook run: | - for nb in $(find nbs -name "*.ipynb" -not -path "*/.ipynb_checkpoints/*" | sort); do - dir=$(dirname "$nb") - filename=$(basename "$nb") - (cd "$dir" && xvfb-run -a uv run --no-sync papermill "$filename" "$filename" -k gsim) - done - - name: Upload executed notebooks + dir=$(dirname "${{ matrix.notebook }}") + filename=$(basename "${{ matrix.notebook }}") + cd "$dir" && xvfb-run -a uv run --no-sync papermill "$filename" "$filename" -k gsim + - name: Upload executed notebook uses: actions/upload-artifact@v4 with: - name: nbs - path: nbs + name: nb-${{ hashFiles(matrix.notebook) }} + path: ${{ matrix.notebook }} retention-days: 1 build: name: Build Docs runs-on: ubuntu-latest - needs: run-notebooks + needs: run-notebook steps: - name: Checkout repo uses: actions/checkout@v4 @@ -68,7 +88,8 @@ jobs: - name: Download all notebook artifacts uses: actions/download-artifact@v4 with: - name: nbs + pattern: nb-* + merge-multiple: true path: nbs - name: Build docs run: | diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index e81294d5..6847dafe 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -4,7 +4,7 @@ on: push: tags: - "[0-9]*.[0-9]*.[0-9]*" - pull_request: + workflow_dispatch: jobs: build: diff --git a/mkdocs.yml b/mkdocs.yml index b3289192..b15f3fbd 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -21,6 +21,7 @@ nav: - examples: - Simple Example: nbs/example.md - MEEP Y-Branch: nbs/meep_ybranch.md + - MEEP Ring Coupler: nbs/meep_ring_coupler.md - Palace CPW: nbs/palace_demo_cpw.md - Palace Microstrip: nbs/palace_demo_microstrip.md - api reference: diff --git a/nbs/meep_api_styles.ipynb b/nbs/meep_api_styles.ipynb new file mode 100644 index 00000000..e0bee08d --- /dev/null +++ b/nbs/meep_api_styles.ipynb @@ -0,0 +1,102 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# MEEP API Styles\n", + "\n", + "The `gsim.meep` API supports three equivalent styles for configuring simulations.\n", + "All three produce identical `Simulation` objects — pick whichever reads best for your use case." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from ubcpdk import PDK, cells\n", + "\n", + "PDK.activate()\n", + "c = cells.ebeam_y_1550()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Callable style (recommended)\n", + "\n", + "Updates fields in place — only the fields you pass change, others keep their current values. No extra imports needed." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": "from gsim import meep\n\nsim = meep.Simulation()\n\nsim.geometry(component=c, z_crop=\"auto\")\nsim.materials = {\"si\": 3.47, \"SiO2\": 1.44}\nsim.source(port=\"o1\", wavelength=1.55, wavelength_span=0.01, num_freqs=11)\nsim.monitors = [\"o1\", \"o2\", \"o3\"]\nsim.domain(pml=1.0, margin=0.5)\nsim.solver(resolution=20, simplify_tol=0.01, save_animation=True, verbose_interval=5.0)" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Attribute style\n", + "\n", + "One field per line. Most explicit — good when you need to set fields conditionally." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": "from gsim import meep\n\nsim = meep.Simulation()\n\nsim.geometry.component = c\nsim.geometry.z_crop = \"auto\"\n\nsim.materials = {\"si\": 3.47, \"SiO2\": 1.44}\n\nsim.source.port = \"o1\"\nsim.source.wavelength = 1.55\nsim.source.wavelength_span = 0.01\nsim.source.num_freqs = 11\n\nsim.monitors = [\"o1\", \"o2\", \"o3\"]\n\nsim.domain.pml = 1.0\nsim.domain.margin = 0.5\n\nsim.solver.resolution = 20\nsim.solver.simplify_tol = 0.01\nsim.solver.save_animation = True\nsim.solver.verbose_interval = 5.0" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 3. Constructor style\n", + "\n", + "Replaces the entire sub-object — fields you don't pass reset to defaults. Requires importing model classes, but useful when building configs programmatically or from saved presets." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": "from gsim import meep\nfrom gsim.meep import Domain, FDTD, Geometry, ModeSource\n\nsim = meep.Simulation()\n\nsim.geometry = Geometry(component=c, z_crop=\"auto\")\nsim.materials = {\"si\": 3.47, \"SiO2\": 1.44}\nsim.source = ModeSource(port=\"o1\", wavelength=1.55, wavelength_span=0.01, num_freqs=11)\nsim.monitors = [\"o1\", \"o2\", \"o3\"]\nsim.domain = Domain(pml=1.0, margin=0.5)\nsim.solver = FDTD(\n resolution=20, simplify_tol=0.01, save_animation=True, verbose_interval=5.0\n)" + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Combining styles\n", + "\n", + "All three styles can be freely combined. For example, use callable for bulk setup, then attribute for conditional tweaks:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": "from gsim import meep\n\nsim = meep.Simulation()\n\nsim.geometry(component=c, z_crop=\"auto\")\nsim.materials = {\"si\": 3.47, \"SiO2\": 1.44}\nsim.source(port=\"o1\", wavelength=1.55, wavelength_span=0.01, num_freqs=11)\nsim.monitors = [\"o1\", \"o2\", \"o3\"]\nsim.domain(pml=1.0, margin=0.5)\nsim.solver(resolution=20, simplify_tol=0.01)\n\n# Conditional tweaks with attribute style\ndebug = True\nif debug:\n sim.solver.save_animation = True\n sim.solver.verbose_interval = 5.0" + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} \ No newline at end of file diff --git a/nbs/meep_ring_coupler.ipynb b/nbs/meep_ring_coupler.ipynb new file mode 100644 index 00000000..42b0c632 --- /dev/null +++ b/nbs/meep_ring_coupler.ipynb @@ -0,0 +1,142 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Running MEEP Simulations\n", + "\n", + "[MEEP](https://meep.readthedocs.io/) is an open-source FDTD electromagnetic simulator. This notebook demonstrates using the `gsim.meep` API to run an S-parameter simulation on a photonic Y-branch.\n", + "\n", + "**Requirements:**\n", + "\n", + "- UBC PDK: `uv pip install ubcpdk`\n", + "- [GDSFactory+](https://gdsfactory.com) account for cloud simulation" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "### Load a pcell from UBC PDK" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "from ubcpdk import PDK, cells\n", + "\n", + "PDK.activate()\n", + "\n", + "c = cells.coupler_ring(gap=0.15)\n", + "\n", + "c" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "### Configure and run simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "from gsim import meep\n", + "\n", + "sim = meep.Simulation()\n", + "\n", + "sim.geometry(component=c, z_crop=\"auto\")\n", + "sim.materials = {\"si\": 3.47, \"SiO2\": 1.44}\n", + "sim.source(port=\"o1\", wavelength=1.55, wavelength_span=0.01, num_freqs=11)\n", + "sim.monitors = [\"o1\", \"o2\", \"o3\", \"o4\"]\n", + "sim.domain(pml=1.0, margin=0.5)\n", + "sim.solver(resolution=20, save_animation=True, verbose_interval=5.0)\n", + "sim.solver.stop_after_sources(time=100)\n", + "\n", + "print(sim.validate_config())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "sim.plot_2d(slices=\"xyz\")" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "### Run simulation on cloud" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "# Run on GDSFactory+ cloud\n", + "result = sim.run()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "result.plot(db=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "result.show_animation()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "gsim", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/nbs/meep_sbend.ipynb b/nbs/meep_sbend.ipynb new file mode 100644 index 00000000..462929c8 --- /dev/null +++ b/nbs/meep_sbend.ipynb @@ -0,0 +1,151 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Running MEEP Simulations\n", + "\n", + "[MEEP](https://meep.readthedocs.io/) is an open-source FDTD electromagnetic simulator. This notebook demonstrates using the `gsim.meep` API to run an S-parameter simulation on a photonic Y-branch.\n", + "\n", + "**Requirements:**\n", + "\n", + "- UBC PDK: `uv pip install ubcpdk`\n", + "- [GDSFactory+](https://gdsfactory.com) account for cloud simulation" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "### Load a pcell from UBC PDK" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "from ubcpdk import PDK, cells\n", + "\n", + "PDK.activate()\n", + "\n", + "c = cells.bend_s(size=(20.0, 5.0))\n", + "\n", + "c" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "### Configure and run simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "from gsim import meep\n", + "\n", + "sim = meep.Simulation()\n", + "\n", + "sim.geometry(component=c, z_crop=\"auto\")\n", + "sim.materials = {\"si\": 3.47, \"SiO2\": 1.44}\n", + "sim.source(port=\"o1\", wavelength=1.55, wavelength_span=0.01, num_freqs=11)\n", + "sim.monitors = [\"o1\", \"o2\"]\n", + "sim.domain(pml=1.0, margin=0.5)\n", + "sim.solver(resolution=20, save_animation=True, verbose_interval=5.0)\n", + "sim.solver.stop_after_sources(time=100)\n", + "\n", + "print(sim.validate_config())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "sim.plot_2d(slices=\"xyz\")" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "### Run simulation on cloud" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "# Run on GDSFactory+ cloud\n", + "result = sim.run()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "result.plot(db=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "result.show_animation()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1de30efa", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "gsim", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/nbs/meep_straight.ipynb b/nbs/meep_straight.ipynb new file mode 100644 index 00000000..b26f5694 --- /dev/null +++ b/nbs/meep_straight.ipynb @@ -0,0 +1,151 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Running MEEP Simulations\n", + "\n", + "[MEEP](https://meep.readthedocs.io/) is an open-source FDTD electromagnetic simulator. This notebook demonstrates using the `gsim.meep` API to run an S-parameter simulation on a photonic Y-branch.\n", + "\n", + "**Requirements:**\n", + "\n", + "- UBC PDK: `uv pip install ubcpdk`\n", + "- [GDSFactory+](https://gdsfactory.com) account for cloud simulation" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "### Load a pcell from UBC PDK" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "from ubcpdk import PDK, cells\n", + "\n", + "PDK.activate()\n", + "\n", + "c = cells.straight(length=20.0)\n", + "\n", + "c" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "### Configure and run simulation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "from gsim import meep\n", + "\n", + "sim = meep.Simulation()\n", + "\n", + "sim.geometry(component=c, z_crop=\"auto\")\n", + "sim.materials = {\"si\": 3.47, \"SiO2\": 1.44}\n", + "sim.source(port=\"o1\", wavelength=1.55, wavelength_span=0.01, num_freqs=11)\n", + "sim.monitors = [\"o1\", \"o2\"]\n", + "sim.domain(pml=1.0, margin=0.5)\n", + "sim.solver(resolution=20, save_animation=True, verbose_interval=5.0)\n", + "sim.solver.stop_after_sources(time=100)\n", + "\n", + "print(sim.validate_config())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "sim.plot_2d(slices=\"xyz\")" + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "### Run simulation on cloud" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "# Run on GDSFactory+ cloud\n", + "result = sim.run()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "result.plot(db=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9", + "metadata": {}, + "outputs": [], + "source": [ + "result.show_animation()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4c03f0c9", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "gsim", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/nbs/meep_ybranch.ipynb b/nbs/meep_ybranch.ipynb index 480f15d0..a4e8adbe 100644 --- a/nbs/meep_ybranch.ipynb +++ b/nbs/meep_ybranch.ipynb @@ -52,43 +52,7 @@ "id": "4", "metadata": {}, "outputs": [], - "source": [ - "from gsim import meep\n", - "\n", - "sim = meep.Simulation()\n", - "\n", - "sim.geometry.component = c\n", - "sim.geometry.z_crop = \"auto\"\n", - "\n", - "sim.materials = {\"si\": 3.47, \"SiO2\": 1.44}\n", - "\n", - "sim.source.port = \"o1\"\n", - "sim.source.wavelength = 1.55\n", - "sim.source.bandwidth = 0.01\n", - "sim.source.num_freqs = 11\n", - "\n", - "sim.monitors = [\"o1\", \"o2\", \"o3\"]\n", - "\n", - "sim.domain.pml = 1.0\n", - "sim.domain.margin = 0.5\n", - "sim.domain.margin_z_above = 0.5\n", - "sim.domain.margin_z_below = 0.5\n", - "sim.domain.port_margin = 0.5\n", - "\n", - "sim.solver.resolution = 20\n", - "sim.solver.stopping = \"dft_decay\"\n", - "sim.solver.max_time = 200\n", - "sim.solver.stopping_threshold = 1e-3\n", - "sim.solver.stopping_min_time = 100\n", - "sim.solver.subpixel = False\n", - "sim.solver.simplify_tol = 0.01\n", - "sim.solver.save_geometry = True\n", - "sim.solver.save_fields = True\n", - "sim.solver.save_animation = True\n", - "sim.solver.verbose_interval = 5.0\n", - "\n", - "print(sim.validate_config())" - ] + "source": "from gsim import meep\n\nsim = meep.Simulation()\n\nsim.geometry(component=c, z_crop=\"auto\")\nsim.materials = {\"si\": 3.47, \"SiO2\": 1.44}\nsim.source(port=\"o1\", wavelength=1.55, wavelength_span=0.01, num_freqs=11)\nsim.monitors = [\"o1\", \"o2\", \"o3\"]\nsim.domain(pml=1.0, margin=0.5)\nsim.solver(resolution=20, simplify_tol=0.01, save_animation=True, verbose_interval=5.0)\nsim.solver.stop_after_sources(time=60)\n\nprint(sim.validate_config())" }, { "cell_type": "code", @@ -142,14 +106,23 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "gsim", "language": "python", "name": "python3" }, "language_info": { - "name": "python" + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.10" } }, "nbformat": 4, "nbformat_minor": 5 -} +} \ No newline at end of file diff --git a/nbs/palace_demo_cpw.ipynb b/nbs/palace_demo_cpw.ipynb index 8c5769c9..72b1504e 100644 --- a/nbs/palace_demo_cpw.ipynb +++ b/nbs/palace_demo_cpw.ipynb @@ -188,17 +188,6 @@ "### Run simulation on cloud" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "8", - "metadata": {}, - "outputs": [], - "source": [ - "# Generate Palace config file\n", - "sim.write_config()" - ] - }, { "cell_type": "code", "execution_count": null, @@ -207,7 +196,7 @@ "outputs": [], "source": [ "# Run simulation on GDSFactory+ cloud\n", - "results = sim.simulate()" + "results = sim.run()" ] }, { @@ -263,7 +252,8 @@ "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython3" + "pygments_lexer": "ipython3", + "version": "3.12.10" } }, "nbformat": 4, diff --git a/nbs/palace_demo_microstrip.ipynb b/nbs/palace_demo_microstrip.ipynb index ad4670c4..f796189c 100644 --- a/nbs/palace_demo_microstrip.ipynb +++ b/nbs/palace_demo_microstrip.ipynb @@ -115,17 +115,6 @@ "# sim.plot_mesh(show_groups=[\"metal\", \"P\"], interactive=True)" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "7", - "metadata": {}, - "outputs": [], - "source": [ - "# Generate Palace config file\n", - "sim.write_config()" - ] - }, { "cell_type": "markdown", "id": "8", @@ -142,7 +131,7 @@ "outputs": [], "source": [ "# Run simulation on GDSFactory+ cloud\n", - "results = sim.simulate()" + "results = sim.run()" ] }, { @@ -198,7 +187,8 @@ "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython3" + "pygments_lexer": "ipython3", + "version": "3.12.10" } }, "nbformat": 4, diff --git a/pyproject.toml b/pyproject.toml index d8359e0a..5776dfd0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -54,7 +54,7 @@ dev = [ "towncrier>=24.0.0", "vega-datasets>=0.9.0", "nbstripout>=0.8.1", - "ty>=0.0.13" + "ty>=0.0.15" ] docs = [ "ihp-gdsfactory>=0.1.4", diff --git a/src/gsim/common/__init__.py b/src/gsim/common/__init__.py index e1f52134..fd9c3804 100644 --- a/src/gsim/common/__init__.py +++ b/src/gsim/common/__init__.py @@ -12,50 +12,17 @@ from gsim.common.geometry import Geometry from gsim.common.geometry_model import GeometryModel, Prism, extract_geometry_model -from gsim.common.stack import ( - MATERIALS_DB, - Layer, - LayerStack, - MaterialProperties, - StackLayer, - ValidationResult, - extract_from_pdk, - extract_layer_stack, - get_material_properties, - get_stack, - load_stack_yaml, - material_is_conductor, - material_is_dielectric, - parse_layer_stack, - plot_stack, - print_stack, - print_stack_table, -) +from gsim.common.stack import LayerStack, ValidationResult # Alias for backward compatibility Stack = LayerStack __all__ = [ - "MATERIALS_DB", "Geometry", "GeometryModel", - "Layer", "LayerStack", - "MaterialProperties", "Prism", "Stack", - "StackLayer", "ValidationResult", - "extract_from_pdk", "extract_geometry_model", - "extract_layer_stack", - "get_material_properties", - "get_stack", - "load_stack_yaml", - "material_is_conductor", - "material_is_dielectric", - "parse_layer_stack", - "plot_stack", - "print_stack", - "print_stack_table", ] diff --git a/src/gsim/common/stack/_layer_utils.py b/src/gsim/common/stack/_layer_utils.py new file mode 100644 index 00000000..229a90ce --- /dev/null +++ b/src/gsim/common/stack/_layer_utils.py @@ -0,0 +1,115 @@ +"""Shared layer classification and GDS-layer extraction helpers. + +Used by both ``extractor`` and ``visualization`` modules. +""" + +from __future__ import annotations + +import logging +from typing import Any, Literal + +from gdsfactory.technology import LayerLevel + +logger = logging.getLogger(__name__) + + +def get_gds_layer_tuple(layer_level: LayerLevel) -> tuple[int, int] | None: + """Extract GDS layer tuple (layer, datatype) from a gdsfactory LayerLevel. + + Returns ``None`` when the layer cannot be parsed (callers decide fallback). + """ + layer: Any = layer_level.layer + + if isinstance(layer, tuple): + return (int(layer[0]), int(layer[1])) + + if isinstance(layer, int): + return (int(layer), 0) + + if hasattr(layer, "layer"): + inner = layer.layer + if hasattr(inner, "layer") and hasattr(inner, "datatype"): + return (int(inner.layer), int(inner.datatype)) + if isinstance(inner, int): + datatype = getattr(layer, "datatype", 0) + return (int(inner), int(datatype) if datatype else 0) + if hasattr(inner, "layer"): + innermost = inner.layer + if isinstance(innermost, int): + datatype = getattr(inner, "datatype", 0) + return (int(innermost), int(datatype) if datatype else 0) + + if hasattr(layer, "layer") and hasattr(layer, "datatype"): + return (int(layer.layer), int(layer.datatype)) + + if hasattr(layer, "value"): + if isinstance(layer.value, tuple): + return (int(layer.value[0]), int(layer.value[1])) + if isinstance(layer.value, int): + return (int(layer.value), 0) + + if isinstance(layer, str): + if "/" in layer: + parts = layer.split("/") + return (int(parts[0]), int(parts[1])) + return None + + try: + return (int(layer), 0) + except (TypeError, ValueError): + logger.warning("Could not parse layer %s", layer) + return None + + +def classify_layer_type( + layer_name: str, + material: str | None = None, +) -> Literal["conductor", "via", "dielectric", "substrate"]: + """Classify a layer as conductor, via, dielectric, or substrate. + + When *material* is provided, the materials DB is consulted for additional + classification hints (conductor vs dielectric). When ``None``, only + name-based heuristics are used. + """ + name_lower = layer_name.lower() + + if "via" in name_lower or "cont" in name_lower: + return "via" + + if any( + m in name_lower for m in ["metal", "topmetal", "m1", "m2", "m3", "m4", "m5"] + ): + return "conductor" + + if "substrate" in name_lower or name_lower == "sub": + return "substrate" + + # Material-based lookup (only when caller supplies material) + if material is not None: + from gsim.common.stack.materials import get_material_properties + + props = get_material_properties(material) + if props: + if props.type == "conductor": + return "conductor" + if props.type == "semiconductor" and "substrate" in name_lower: + return "substrate" + if props.type == "dielectric": + return "dielectric" + + material_lower = material.lower() + if material_lower in [ + "aluminum", + "copper", + "tungsten", + "gold", + "al", + "cu", + "w", + ]: + return "conductor" + + if "poly" in name_lower or "active" in name_lower: + return "conductor" + + return "dielectric" diff --git a/src/gsim/common/stack/extractor.py b/src/gsim/common/stack/extractor.py index fd8a6090..63e437a8 100644 --- a/src/gsim/common/stack/extractor.py +++ b/src/gsim/common/stack/extractor.py @@ -8,13 +8,13 @@ import logging from pathlib import Path -from typing import Any, Literal +from typing import Literal import yaml -from gdsfactory.technology import LayerLevel from gdsfactory.technology import LayerStack as GfLayerStack from pydantic import BaseModel, ConfigDict, Field +from gsim.common.stack._layer_utils import classify_layer_type, get_gds_layer_tuple from gsim.common.stack.materials import ( MATERIALS_DB, get_material_properties, @@ -279,87 +279,6 @@ def to_yaml(self, path: Path | None = None) -> str: return yaml_str -def _get_gds_layer_tuple(layer_level: LayerLevel) -> tuple[int, int]: - """Extract GDS layer tuple from LayerLevel.""" - layer: Any = layer_level.layer - - if isinstance(layer, tuple): - return (int(layer[0]), int(layer[1])) - - if isinstance(layer, int): - return (int(layer), 0) - - if hasattr(layer, "layer"): - inner = layer.layer - if hasattr(inner, "layer") and hasattr(inner, "datatype"): - return (int(inner.layer), int(inner.datatype)) # type: ignore[arg-type] - if isinstance(inner, int): - datatype = getattr(layer, "datatype", 0) - return (int(inner), int(datatype) if datatype else 0) - if hasattr(inner, "layer"): - innermost = inner.layer - if isinstance(innermost, int): - datatype = getattr(inner, "datatype", 0) - return (int(innermost), int(datatype) if datatype else 0) - - if hasattr(layer, "layer") and hasattr(layer, "datatype"): - return (int(layer.layer), int(layer.datatype)) # type: ignore[arg-type] - - if hasattr(layer, "value"): - if isinstance(layer.value, tuple): - return (int(layer.value[0]), int(layer.value[1])) # type: ignore[arg-type] - if isinstance(layer.value, int): - return (int(layer.value), 0) - - if isinstance(layer, str): - if "/" in layer: - parts = layer.split("/") - return (int(parts[0]), int(parts[1])) - return (0, 0) - - try: - return (int(layer), 0) # type: ignore[arg-type] - except (TypeError, ValueError): - logger.warning("Could not parse layer %s, using (0, 0)", layer) - return (0, 0) - - -def _classify_layer_type( - layer_name: str, material: str -) -> Literal["conductor", "via", "dielectric", "substrate"]: - """Classify a layer as conductor, via, dielectric, or substrate.""" - name_lower = layer_name.lower() - material_lower = material.lower() - - if "via" in name_lower: - return "via" - - if any( - m in name_lower for m in ["metal", "topmetal", "m1", "m2", "m3", "m4", "m5"] - ): - return "conductor" - - if "substrate" in name_lower or name_lower == "sub": - return "substrate" - - props = get_material_properties(material) - if props: - if props.type == "conductor": - return "conductor" - if props.type == "semiconductor" and "substrate" in name_lower: - return "substrate" - if props.type == "dielectric": - return "dielectric" - - if material_lower in ["aluminum", "copper", "tungsten", "gold", "al", "cu", "w"]: - return "conductor" - - if "poly" in name_lower: - return "conductor" - - return "dielectric" - - def extract_layer_stack( gf_layer_stack: GfLayerStack, pdk_name: str = "unknown", @@ -392,8 +311,8 @@ def extract_layer_stack( thickness = layer_level.thickness if layer_level.thickness is not None else 0.0 zmax = zmin + thickness material = layer_level.material or "unknown" - gds_layer = _get_gds_layer_tuple(layer_level) - layer_type = _classify_layer_type(layer_name, material) + gds_layer = get_gds_layer_tuple(layer_level) or (0, 0) + layer_type = classify_layer_type(layer_name, material) sidewall_angle = getattr(layer_level, "sidewall_angle", 0.0) or 0.0 if layer_type == "substrate" and not include_substrate: diff --git a/src/gsim/common/stack/visualization.py b/src/gsim/common/stack/visualization.py index 18a99731..7bbeb0ef 100644 --- a/src/gsim/common/stack/visualization.py +++ b/src/gsim/common/stack/visualization.py @@ -7,12 +7,12 @@ from __future__ import annotations from dataclasses import dataclass -from typing import Any import plotly.graph_objects as go -from gdsfactory.technology import LayerLevel from gdsfactory.technology import LayerStack as GfLayerStack +from gsim.common.stack._layer_utils import classify_layer_type, get_gds_layer_tuple + @dataclass class StackLayer: @@ -27,52 +27,6 @@ class StackLayer: layer_type: str = "conductor" # conductor, via, dielectric, substrate -def _get_gds_layer_number(layer_level: LayerLevel) -> int | None: - """Extract GDS layer number from LayerLevel.""" - layer: Any = layer_level.layer - - # Handle tuple - if isinstance(layer, tuple): - return int(layer[0]) - - # Handle int - if isinstance(layer, int): - return int(layer) - - # Handle LogicalLayer or enum with nested layer - if hasattr(layer, "layer"): - inner = layer.layer - if hasattr(inner, "layer"): - return int(inner.layer) # type: ignore[arg-type] - if isinstance(inner, int): - return int(inner) - - # Handle enum with value - if hasattr(layer, "value"): - if isinstance(layer.value, tuple): - return int(layer.value[0]) # type: ignore[arg-type] - return int(layer.value) # type: ignore[arg-type] - - return None - - -def _classify_layer(name: str) -> str: - """Classify layer type based on name.""" - name_lower = name.lower() - - if "via" in name_lower or "cont" in name_lower: - return "via" - if "substrate" in name_lower or name_lower == "sub": - return "substrate" - if any( - m in name_lower - for m in ["metal", "topmetal", "m1", "m2", "m3", "m4", "m5", "poly", "active"] - ): - return "conductor" - - return "dielectric" - - def parse_layer_stack(layer_stack: GfLayerStack) -> list[StackLayer]: """Parse a gdsfactory LayerStack into a list of StackLayer objects. @@ -89,8 +43,9 @@ def parse_layer_stack(layer_stack: GfLayerStack) -> list[StackLayer]: thickness = level.thickness if level.thickness is not None else 0.0 zmax = zmin + thickness material = level.material or None - gds_layer = _get_gds_layer_number(level) - layer_type = _classify_layer(name) + tup = get_gds_layer_tuple(level) + gds_layer = tup[0] if tup else None + layer_type = classify_layer_type(name) layers.append( StackLayer( diff --git a/src/gsim/common/viz/__init__.py b/src/gsim/common/viz/__init__.py index e472cc85..030b67d3 100644 --- a/src/gsim/common/viz/__init__.py +++ b/src/gsim/common/viz/__init__.py @@ -6,7 +6,6 @@ 3D backends: - PyVista (desktop, ``plot_prisms_3d``) - Open3D + Plotly (Jupyter, ``plot_prisms_3d_open3d``) - - Three.js / FastAPI (browser, ``serve_threejs_visualization``) 2D backends: - matplotlib (``plot_prism_slices``) @@ -18,7 +17,6 @@ export_3d_mesh, plot_prisms_3d, plot_prisms_3d_open3d, - serve_threejs_visualization, ) __all__ = [ @@ -27,5 +25,4 @@ "plot_prism_slices", "plot_prisms_3d", "plot_prisms_3d_open3d", - "serve_threejs_visualization", ] diff --git a/src/gsim/common/viz/_mesh_helpers.py b/src/gsim/common/viz/_mesh_helpers.py index 690b579c..97feb6fa 100644 --- a/src/gsim/common/viz/_mesh_helpers.py +++ b/src/gsim/common/viz/_mesh_helpers.py @@ -6,6 +6,8 @@ from __future__ import annotations +from typing import Any + import numpy as np from gsim.common.geometry_model import GeometryModel, Prism @@ -72,3 +74,99 @@ def simulation_box_corners( ] return points, lines + + +# --------------------------------------------------------------------------- +# Shared Delaunay triangulation for polygons with holes / concave polygons +# --------------------------------------------------------------------------- + + +def triangulate_polygon_with_holes( + shapely_polygon: Any, + z_base: float, + z_top: float, +) -> tuple[np.ndarray, list[np.ndarray], list[list[int]]] | None: + """Triangulate a Shapely polygon (possibly with holes) via Delaunay. + + Returns: + ``(verts_3d, valid_triangles, boundary_segments)`` where + *verts_3d* is (2*N, 3) — bottom then top copies of the 2D points, + *valid_triangles* is a list of index-triplet arrays (into the first + N points), and *boundary_segments* is a list of ``[i, j]`` pairs. + Returns ``None`` when triangulation fails. + """ + try: + import shapely.geometry as sg + from scipy.spatial import Delaunay + except ImportError: + return None + + all_points: list[tuple[float, float]] = [] + boundary_segments: list[list[int]] = [] + + exterior_coords = list(shapely_polygon.exterior.coords[:-1]) + start_idx = 0 + all_points.extend(exterior_coords) + boundary_segments = [ + [start_idx + i, start_idx + (i + 1) % len(exterior_coords)] + for i in range(len(exterior_coords)) + ] + + for interior in shapely_polygon.interiors: + interior_coords = list(interior.coords[:-1]) + start_idx = len(all_points) + all_points.extend(interior_coords) + boundary_segments.extend( + [start_idx + i, start_idx + (i + 1) % len(interior_coords)] + for i in range(len(interior_coords)) + ) + + if len(all_points) < 3: + return None + + points_2d = np.array(all_points) + tri = Delaunay(points_2d) + + valid_triangles = [ + simplex + for simplex in tri.simplices + if shapely_polygon.contains(sg.Point(*np.mean(points_2d[simplex], axis=0))) + ] + + if not valid_triangles: + return None + + verts_3d = np.array( + [[pt[0], pt[1], z_base] for pt in points_2d] + + [[pt[0], pt[1], z_top] for pt in points_2d] + ) + + return verts_3d, valid_triangles, boundary_segments + + +def collect_triangular_prism_geometry( + prisms: list[Prism], +) -> tuple[np.ndarray, list[int]] | None: + """Collect vertices from many triangular prisms for batch meshing. + + Returns: + ``(combined_vertices, prism_offsets)`` where *combined_vertices* + is (M, 3) and *prism_offsets[i]* is the vertex offset of the + *i*-th prism in that array. Returns ``None`` if *prisms* is empty. + """ + all_vertices: list[np.ndarray] = [] + prism_offsets: list[int] = [] + offset = 0 + + for prism in prisms: + base, top = prism_base_top_vertices(prism) + verts = np.vstack([base, top]) + all_vertices.append(verts) + prism_offsets.append(offset) + offset += len(verts) + + if not all_vertices: + return None + + combined = np.vstack(all_vertices) + return combined, prism_offsets diff --git a/src/gsim/common/viz/render3d.py b/src/gsim/common/viz/render3d.py index ae71fbcc..db9ad48d 100644 --- a/src/gsim/common/viz/render3d.py +++ b/src/gsim/common/viz/render3d.py @@ -12,12 +12,10 @@ export_3d_mesh, plot_prisms_3d, ) -from gsim.common.viz.render3d_threejs import serve_threejs_visualization __all__ = [ "create_web_export", "export_3d_mesh", "plot_prisms_3d", "plot_prisms_3d_open3d", - "serve_threejs_visualization", ] diff --git a/src/gsim/common/viz/render3d_open3d.py b/src/gsim/common/viz/render3d_open3d.py index 528211b3..778fd215 100644 --- a/src/gsim/common/viz/render3d_open3d.py +++ b/src/gsim/common/viz/render3d_open3d.py @@ -13,8 +13,10 @@ from gsim.common.geometry_model import GeometryModel, Prism from gsim.common.viz._colors import generate_layer_colors_with_opacity from gsim.common.viz._mesh_helpers import ( + collect_triangular_prism_geometry, prism_base_top_vertices, simulation_box_corners, + triangulate_polygon_with_holes, ) try: @@ -120,23 +122,58 @@ def plot_prisms_3d_open3d( cx = cy = cz = 0.0 rs = 10.0 + # Standard camera views for view buttons + d = 2.0 # distance multiplier for orthographic views + _views = { + "Iso": dict(eye=dict(x=1.5, y=1.5, z=1.5)), + "Top": dict(eye=dict(x=0, y=0, z=d), up=dict(x=0, y=1, z=0)), + "Front": dict(eye=dict(x=0, y=-d, z=0), up=dict(x=0, y=0, z=1)), + "Right": dict(eye=dict(x=d, y=0, z=0), up=dict(x=0, y=0, z=1)), + } + view_buttons = [ + dict( + label=name, + method="relayout", + args=[ + { + "scene.camera.eye": cam["eye"], + "scene.camera.up": cam.get("up", dict(x=0, y=0, z=1)), + "scene.camera.center": dict(x=0, y=0, z=0), + } + ], + ) + for name, cam in _views.items() + ] + fig.update_layout( scene=dict( xaxis_title="X (um)", yaxis_title="Y (um)", zaxis_title="Z (um)", - aspectmode="data", + aspectmode="cube", camera=dict( eye=dict(x=1.5, y=1.5, z=1.5), center=dict(x=0, y=0, z=0), - projection=dict(type="perspective"), + projection=dict(type="orthographic"), ), xaxis=dict(range=[cx - rs * 0.6, cx + rs * 0.6]), yaxis=dict(range=[cy - rs * 0.6, cy + rs * 0.6]), zaxis=dict(range=[cz - rs * 0.6, cz + rs * 0.6]), ), - title="3D Geometry Visualization", + title="", dragmode="orbit", + updatemenus=[ + dict( + type="buttons", + direction="down", + x=0.01, + y=0.99, + xanchor="left", + yanchor="top", + bgcolor="rgba(255,255,255,0.8)", + buttons=view_buttons, + ), + ], **kwargs, ) @@ -199,16 +236,15 @@ def _convert_prisms_to_open3d( def _merge_triangular_prisms_o3d(prisms: list[Prism]) -> Any: """Merge many triangular prisms into one Open3D mesh for performance.""" - all_vertices: list[np.ndarray] = [] - all_triangles: list[list[int]] = [] - offset = 0 + result = collect_triangular_prism_geometry(prisms) + if result is None: + return None - for prism in prisms: - base, top = prism_base_top_vertices(prism) - verts = np.vstack([base, top]) - all_vertices.append(verts) - n = len(base) + combined, prism_offsets = result + n = 3 # triangular prisms + all_triangles: list[list[int]] = [] + for offset in prism_offsets: # bottom all_triangles.append([offset, offset + 2, offset + 1]) # top @@ -219,16 +255,11 @@ def _merge_triangular_prisms_o3d(prisms: list[Prism]) -> Any: all_triangles.append([offset + i, offset + ni, offset + ni + n]) all_triangles.append([offset + i, offset + ni + n, offset + i + n]) - offset += len(verts) - - if all_vertices: - combined = np.vstack(all_vertices) - mesh = o3d.geometry.TriangleMesh() - mesh.vertices = o3d.utility.Vector3dVector(combined) - mesh.triangles = o3d.utility.Vector3iVector(all_triangles) - mesh.compute_vertex_normals() - return mesh - return None + mesh = o3d.geometry.TriangleMesh() + mesh.vertices = o3d.utility.Vector3dVector(combined) + mesh.triangles = o3d.utility.Vector3iVector(all_triangles) + mesh.compute_vertex_normals() + return mesh def _create_prism_mesh_o3d( @@ -273,62 +304,26 @@ def _create_prism_mesh_with_holes_o3d( top_vertices: np.ndarray, ) -> Any: """Create Open3D mesh from a Shapely polygon with holes via Delaunay.""" - try: - import shapely.geometry as sg - from scipy.spatial import Delaunay - except ImportError: - return _create_prism_mesh_o3d(base_vertices, top_vertices) - - all_points: list[tuple[float, float]] = [] - boundary_segments: list[list[int]] = [] - - exterior_coords = list(shapely_polygon.exterior.coords[:-1]) - start_idx = 0 - all_points.extend(exterior_coords) - boundary_segments = [ - [start_idx + i, start_idx + (i + 1) % len(exterior_coords)] - for i in range(len(exterior_coords)) - ] - - for interior in shapely_polygon.interiors: - interior_coords = list(interior.coords[:-1]) - start_idx = len(all_points) - all_points.extend(interior_coords) - boundary_segments.extend( - [start_idx + i, start_idx + (i + 1) % len(interior_coords)] - for i in range(len(interior_coords)) - ) - - if len(all_points) < 3: - return _create_prism_mesh_o3d(base_vertices, top_vertices) - - points_2d = np.array(all_points) - tri = Delaunay(points_2d) - - valid_triangles = [ - simplex - for simplex in tri.simplices - if shapely_polygon.contains(sg.Point(*np.mean(points_2d[simplex], axis=0))) - ] - - if not valid_triangles: - return _create_prism_mesh_o3d(base_vertices, top_vertices) - z_base = base_vertices[0, 2] if len(base_vertices) > 0 else 0 z_top = top_vertices[0, 2] if len(top_vertices) > 0 else z_base + 1 - verts_3d: list[list[float]] = [[pt[0], pt[1], z_base] for pt in points_2d] + [ - [pt[0], pt[1], z_top] for pt in points_2d - ] + result = triangulate_polygon_with_holes(shapely_polygon, z_base, z_top) + if result is None: + return _create_prism_mesh_o3d(base_vertices, top_vertices) + + verts_3d, valid_triangles, boundary_segments = result + n_pts = len(verts_3d) // 2 - n_pts = len(points_2d) + # bottom triangles faces_3d: list[list[int]] = [ [tri_idx[0], tri_idx[1], tri_idx[2]] for tri_idx in valid_triangles ] + # top triangles (reversed winding) faces_3d.extend( [tri_idx[0] + n_pts, tri_idx[2] + n_pts, tri_idx[1] + n_pts] for tri_idx in valid_triangles ) + # side quads (2 tris each) for seg in boundary_segments: i, j = seg faces_3d.append([i, j, j + n_pts]) diff --git a/src/gsim/common/viz/render3d_pyvista.py b/src/gsim/common/viz/render3d_pyvista.py index f614e5d1..2dd13a28 100644 --- a/src/gsim/common/viz/render3d_pyvista.py +++ b/src/gsim/common/viz/render3d_pyvista.py @@ -12,7 +12,9 @@ from gsim.common.geometry_model import GeometryModel, Prism from gsim.common.viz._colors import generate_layer_colors from gsim.common.viz._mesh_helpers import ( + collect_triangular_prism_geometry, prism_base_top_vertices, + triangulate_polygon_with_holes, ) try: @@ -215,16 +217,15 @@ def _convert_prisms_to_meshes( def _merge_triangular_prisms(prisms: list[Prism]) -> Any: """Merge many triangular prisms into one PyVista mesh for performance.""" - all_vertices: list[np.ndarray] = [] - all_faces: list[int] = [] - offset = 0 + result = collect_triangular_prism_geometry(prisms) + if result is None: + return None - for prism in prisms: - base, top = prism_base_top_vertices(prism) - verts = np.vstack([base, top]) - all_vertices.append(verts) - n = len(base) + combined, prism_offsets = result + n = 3 # triangular prisms + all_faces: list[int] = [] + for offset in prism_offsets: # bottom face all_faces.extend([3, offset + 2, offset + 1, offset + 0]) # top face @@ -233,20 +234,10 @@ def _merge_triangular_prisms(prisms: list[Prism]) -> Any: for i in range(n): ni = (i + 1) % n all_faces.extend( - [ - 4, - offset + i, - offset + ni, - offset + ni + n, - offset + i + n, - ] + [4, offset + i, offset + ni, offset + ni + n, offset + i + n] ) - offset += len(verts) - if all_vertices: - combined = np.vstack(all_vertices) - return pv.PolyData(combined, all_faces) - return None + return pv.PolyData(combined, all_faces) def _create_prism_mesh( @@ -288,68 +279,24 @@ def _create_prism_mesh_with_holes( top_vertices: np.ndarray, ) -> Any: """Create PyVista mesh from a Shapely polygon with holes using Delaunay.""" - try: - import shapely.geometry as sg - from scipy.spatial import Delaunay - except ImportError: - return _create_prism_mesh(base_vertices, top_vertices) - - all_points: list[tuple[float, float]] = [] - boundary_segments: list[list[int]] = [] - - exterior_coords = list(shapely_polygon.exterior.coords[:-1]) - start_idx = 0 - all_points.extend(exterior_coords) - boundary_segments = [ - [start_idx + i, start_idx + (i + 1) % len(exterior_coords)] - for i in range(len(exterior_coords)) - ] - - for interior in shapely_polygon.interiors: - interior_coords = list(interior.coords[:-1]) - start_idx = len(all_points) - all_points.extend(interior_coords) - boundary_segments.extend( - [start_idx + i, start_idx + (i + 1) % len(interior_coords)] - for i in range(len(interior_coords)) - ) - - if len(all_points) < 3: - return _create_prism_mesh(base_vertices, top_vertices) - - points_2d = np.array(all_points) - tri = Delaunay(points_2d) - - valid_triangles = [ - simplex - for simplex in tri.simplices - if shapely_polygon.contains(sg.Point(*np.mean(points_2d[simplex], axis=0))) - ] - - if not valid_triangles: - return _create_prism_mesh(base_vertices, top_vertices) - z_base = base_vertices[0, 2] if len(base_vertices) > 0 else 0 z_top = top_vertices[0, 2] if len(top_vertices) > 0 else z_base + 1 - verts_3d: list[list[float]] = [[pt[0], pt[1], z_base] for pt in points_2d] + [ - [pt[0], pt[1], z_top] for pt in points_2d - ] + result = triangulate_polygon_with_holes(shapely_polygon, z_base, z_top) + if result is None: + return _create_prism_mesh(base_vertices, top_vertices) - n_pts = len(points_2d) + verts_3d, valid_triangles, boundary_segments = result + n_pts = len(verts_3d) // 2 faces_pv: list[int] = [] + # bottom triangles for tri_idx in valid_triangles: faces_pv.extend([3, tri_idx[0], tri_idx[1], tri_idx[2]]) + # top triangles (reversed winding) for tri_idx in valid_triangles: - faces_pv.extend( - [ - 3, - tri_idx[0] + n_pts, - tri_idx[2] + n_pts, - tri_idx[1] + n_pts, - ] - ) + faces_pv.extend([3, tri_idx[0] + n_pts, tri_idx[2] + n_pts, tri_idx[1] + n_pts]) + # side quads for seg in boundary_segments: i, j = seg faces_pv.extend([4, i, j, j + n_pts, i + n_pts]) diff --git a/src/gsim/common/viz/render3d_threejs.py b/src/gsim/common/viz/render3d_threejs.py deleted file mode 100644 index d9daf3c0..00000000 --- a/src/gsim/common/viz/render3d_threejs.py +++ /dev/null @@ -1,233 +0,0 @@ -"""Three.js / FastAPI-based 3D rendering for GeometryModel. - -Provides browser-based interactive 3D visualisation via a local FastAPI server -that serves a Three.js HTML page. -""" - -from __future__ import annotations - -import json -import logging -from pathlib import Path -from typing import Any - -import numpy as np - -from gsim.common.geometry_model import GeometryModel -from gsim.common.viz._colors import generate_layer_colors_with_opacity -from gsim.common.viz._mesh_helpers import simulation_box_corners -from gsim.common.viz.render3d_open3d import _convert_prisms_to_open3d - -logger = logging.getLogger(__name__) - -# --------------------------------------------------------------------------- -# Public API -# --------------------------------------------------------------------------- - - -def serve_threejs_visualization( - geometry_model: GeometryModel, - *, - show_edges: bool = False, - color_by_layer: bool = True, - show_simulation_box: bool = True, - layer_opacity: dict[str, float] | None = None, - port: int = 8000, - auto_open: bool = True, - show_stats: bool = False, - **kwargs: Any, -) -> str: - """Start a FastAPI server with a Three.js 3D viewer. - - Args: - geometry_model: GeometryModel containing prisms and bbox. - show_edges: Show wireframe edges. - color_by_layer: Colour each layer differently. - show_simulation_box: Draw the simulation box. - layer_opacity: Per-layer opacity override. - port: Port to serve on (default 8000). - auto_open: Open the browser automatically. - show_stats: Show FPS counter. - **kwargs: Extra Three.js options. - - Returns: - URL of the running server. - """ - try: - import threading - import time - import webbrowser - - import uvicorn - from fastapi import FastAPI - from fastapi.responses import HTMLResponse - except ImportError as err: - raise ImportError( - "FastAPI and uvicorn required. Install with: pip install fastapi uvicorn" - ) from err - - app = FastAPI(title="3D Geometry Viewer") - - layer_meshes = _convert_prisms_to_open3d(geometry_model) - colors, opacity_dict = generate_layer_colors_with_opacity( - list(layer_meshes.keys()), layer_opacity - ) - - logger.info("Converting geometry: %d layers found", len(layer_meshes)) - for layer_name, meshes in layer_meshes.items(): - logger.info(" Layer '%s': %d meshes", layer_name, len(meshes)) - - threejs_data = _convert_to_threejs_data( - layer_meshes, colors, opacity_dict, color_by_layer - ) - - if show_simulation_box: - threejs_data["simulation_box"] = _create_simulation_box_threejs(geometry_model) - - total_vertices = 0 - total_faces = 0 - for layer in threejs_data.get("layers", []): - for mesh in layer.get("meshes", []): - total_vertices += len(mesh.get("vertices", [])) // 3 - total_faces += len(mesh.get("faces", [])) // 3 - logger.info( - "Three.js data prepared: %d vertices, %d faces", - total_vertices, - total_faces, - ) - - @app.get("/", response_class=HTMLResponse) - def get_visualization(): - """Return the Three.js viewer HTML page.""" - return _generate_threejs_html( - threejs_data, - show_edges=show_edges, - show_stats=show_stats, - **kwargs, - ) - - @app.get("/api/geometry") - def get_geometry_data(): - """Return geometry mesh data as JSON.""" - return threejs_data - - @app.get("/api/info") - def get_info(): - """Return summary info about the geometry.""" - return { - "layers": len(threejs_data.get("layers", [])), - "total_meshes": sum( - len(layer.get("meshes", [])) for layer in threejs_data.get("layers", []) - ), - "has_simulation_box": "simulation_box" in threejs_data, - "layer_names": [ - layer.get("name") for layer in threejs_data.get("layers", []) - ], - } - - server_url = f"http://localhost:{port}" - - def run_server(): - try: - uvicorn.run(app, host="127.0.0.1", port=port, log_level="info") - except Exception as e: - logger.info("Server error: %s", e) - - server_thread = threading.Thread(target=run_server, daemon=True) - server_thread.start() - time.sleep(1) - - if auto_open: - webbrowser.open(server_url) - - logger.info("Server started at: %s", server_url) - logger.info("Geometry API available at: %s/api/geometry", server_url) - logger.info("Server running in background. Keep Python session alive to view.") - - return server_url - - -# --------------------------------------------------------------------------- -# Internal helpers -# --------------------------------------------------------------------------- - - -def _convert_to_threejs_data( - layer_meshes: dict[str, list[Any]], - colors: dict[str, list[float]], - opacity_dict: dict[str, float], - color_by_layer: bool, -) -> dict[str, Any]: - """Convert Open3D meshes to Three.js-compatible JSON data.""" - threejs_data: dict[str, Any] = {"layers": []} - - for layer_name, meshes in layer_meshes.items(): - layer_color = colors[layer_name] if color_by_layer else [0.7, 0.7, 0.7] - layer_opacity = opacity_dict.get(layer_name, 0.8) - - layer_data: dict[str, Any] = { - "name": layer_name, - "color": [int(c * 255) for c in layer_color[:3]], - "opacity": layer_opacity, - "meshes": [], - } - - for i, mesh in enumerate(meshes): - vertices = np.asarray(mesh.vertices).flatten().tolist() - faces = np.asarray(mesh.triangles).flatten().tolist() - - layer_data["meshes"].append( - { - "vertices": vertices, - "faces": faces, - "id": f"{layer_name}_{i}", - } - ) - - threejs_data["layers"].append(layer_data) - - return threejs_data - - -def _create_simulation_box_threejs( - geometry_model: GeometryModel, -) -> dict[str, Any]: - """Create simulation box data for Three.js.""" - points, lines = simulation_box_corners(geometry_model) - return { - "vertices": points.flatten().tolist(), - "lines": [idx for pair in lines for idx in pair], - "color": [0, 0, 0], - } - - -def _generate_threejs_html( - threejs_data: dict[str, Any], - *, - show_edges: bool = False, - show_stats: bool = False, - **kwargs: Any, -) -> str: - """Generate HTML using the viewer.html template.""" - template_path = Path(__file__).parent / "templates" / "viewer.html" - - with open(template_path) as f: - template = f.read() - - width = kwargs.get("width", "100vw") - height = kwargs.get("height", "100vh") - background_color = kwargs.get("background_color", "#f0f0f0") - show_wireframe = str(show_edges).lower() - stats_display = "block" if show_stats else "none" - - data_json = json.dumps(threejs_data, indent=2) - - return template.format( - geometry_data=data_json, - show_wireframe=show_wireframe, - show_stats=str(show_stats).lower(), - width=width, - height=height, - background_color=background_color, - stats_display=stats_display, - ) diff --git a/src/gsim/common/viz/templates/viewer.html b/src/gsim/common/viz/templates/viewer.html deleted file mode 100644 index 23323878..00000000 --- a/src/gsim/common/viz/templates/viewer.html +++ /dev/null @@ -1,482 +0,0 @@ - - - - - - 3D Geometry - Three.js Viewer - - - -
-
-

🔬 3D Geometry Visualization

-

Mouse: Rotate view

-

Wheel: Zoom in/out (enhanced sensitivity)

-

Right-click: Pan camera

-

Double-click: Reset view

-
- -
- -
- -
- -
- -
-

📊 Layers

-
-
- -
-
FPS: --
-
Vertices: --
-
Faces: --
-
-
- - - - - - - - - diff --git a/src/gsim/gcloud.py b/src/gsim/gcloud.py index 2b2cb508..a6bd8a02 100644 --- a/src/gsim/gcloud.py +++ b/src/gsim/gcloud.py @@ -7,15 +7,19 @@ from gsim import gcloud # Run simulation (uploads, starts, waits, downloads) - results = gcloud.run_simulation("./sim", job_type="palace") + result = gcloud.run_simulation("./sim", job_type="palace") + print(result.sim_dir) # palace_abc123/ + print(result.files) # {"port-S.csv": Path(...), ...} # Or use solver-specific wrappers: from gsim import palace as pa - results = pa.run_simulation("./sim") + result = pa.run_simulation("./sim") """ from __future__ import annotations +import shutil +from dataclasses import dataclass, field from pathlib import Path from typing import TYPE_CHECKING @@ -26,6 +30,38 @@ from typing import Literal +@dataclass +class RunResult: + """Result of a cloud simulation run. + + Attributes: + sim_dir: Root directory (``{job_type}_{job_name}/``). + files: Flat mapping of filename → Path inside ``output/``. + job_name: Cloud job identifier. + """ + + sim_dir: Path + files: dict[str, Path] = field(default_factory=dict) + job_name: str = "" + + +def _flatten_results(raw_results: dict) -> dict[str, Path]: + """Flatten gdsfactoryplus download results to a filename → Path dict. + + The SDK may return directories (extracted archives) or individual files. + This walks everything and returns a flat mapping. + """ + flat: dict[str, Path] = {} + for result_path in raw_results.values(): + if result_path.is_dir(): + for file_path in result_path.rglob("*"): + if file_path.is_file() and not file_path.name.startswith("."): + flat[file_path.name] = file_path + else: + flat[result_path.name] = result_path + return flat + + def _handle_failed_job(job, output_dir: Path, verbose: bool) -> None: """Handle a failed simulation job by downloading logs and raising informative error. @@ -54,16 +90,7 @@ def _handle_failed_job(job, output_dir: Path, verbose: bool) -> None: print("Downloading logs from failed job...") # noqa: T201 raw_results = sim.download_results(job, output_dir=output_dir) - - # Flatten results to find all files - all_files: dict[str, Path] = {} - for result_path in raw_results.values(): - if result_path.is_dir(): - for file_path in result_path.rglob("*"): - if file_path.is_file() and not file_path.name.startswith("."): - all_files[file_path.name] = file_path - else: - all_files[result_path.name] = result_path + all_files = _flatten_results(raw_results) # Look for log files and display them log_files = ["palace.log", "stdout.log", "stderr.log", "output.log"] @@ -109,48 +136,55 @@ def upload_simulation_dir(input_dir: str | Path, job_type: str): def run_simulation( - output_dir: str | Path, + config_dir: str | Path, job_type: Literal["palace", "meep"] = "palace", verbose: bool = True, on_started: Callable | None = None, -) -> dict[str, Path]: + parent_dir: str | Path | None = None, +) -> RunResult: """Run a simulation on GDSFactory+ cloud. This function handles the complete workflow: - 1. Uploads simulation files + 1. Uploads simulation files from *config_dir* 2. Starts the simulation job - 3. Waits for completion - 4. Downloads results + 3. Creates a structured directory ``sim-data-{job_name}/`` + with ``input/`` (config files) and ``output/`` (results) sub-dirs + 4. Waits for completion + 5. Downloads results into ``output/`` Args: - output_dir: Directory containing the simulation files - job_type: Type of simulation (default: "palace") - verbose: Print progress messages (default True) - on_started: Optional callback called with job object when simulation starts + config_dir: Directory containing the simulation config files. + job_type: Type of simulation (default: "palace"). + verbose: Print progress messages (default True). + on_started: Optional callback called with job object when simulation starts. + parent_dir: Where to create the sim directory. + Defaults to the current working directory. Returns: - Dict mapping result filename to local Path. + RunResult with sim_dir, files dict, and job_name. Raises: RuntimeError: If simulation fails Example: - >>> results = gcloud.run_simulation("./sim", job_type="palace") + >>> result = gcloud.run_simulation("./sim", job_type="palace") Uploading simulation... done Job started: palace-abc123 Waiting for completion... done (2m 34s) Downloading results... done + >>> print(result.sim_dir) + sim-data-palace-abc123/ """ - output_dir = Path(output_dir) + config_dir = Path(config_dir) - if not output_dir.exists(): - raise FileNotFoundError(f"Output directory not found: {output_dir}") + if not config_dir.exists(): + raise FileNotFoundError(f"Config directory not found: {config_dir}") # Upload if verbose: print("Uploading simulation... ", end="", flush=True) # noqa: T201 - pre_job = upload_simulation_dir(output_dir, job_type) + pre_job = upload_simulation_dir(config_dir, job_type) if verbose: print("done") # noqa: T201 @@ -164,6 +198,20 @@ def run_simulation( if on_started: on_started(job) + # Create structured directory + root = Path(parent_dir) if parent_dir else Path.cwd() + sim_dir = root / f"sim-data-{job.job_name}" + input_dir = sim_dir / "input" + output_dir = sim_dir / "output" + input_dir.mkdir(parents=True, exist_ok=True) + output_dir.mkdir(parents=True, exist_ok=True) + + # Move config files into input/ + for item in list(config_dir.iterdir()): + shutil.move(str(item), str(input_dir / item.name)) + # Remove now-empty config_dir (may fail if it was CWD, etc.) + shutil.rmtree(config_dir, ignore_errors=True) + # Wait finished_job = sim.wait_for_simulation(job) @@ -171,29 +219,14 @@ def run_simulation( if finished_job.exit_code != 0: _handle_failed_job(finished_job, output_dir, verbose) - # Download - raw_results = sim.download_results( - finished_job, output_dir=f"sim-data-{finished_job.job_name}" - ) - - # Flatten results: gdsfactoryplus returns extracted directories, - # but we want a dict of filename -> Path for individual files - results: dict[str, Path] = {} - for result_path in raw_results.values(): - if result_path.is_dir(): - # Recursively find all files in the extracted directory - for file_path in result_path.rglob("*"): - if file_path.is_file() and not file_path.name.startswith("."): - results[file_path.name] = file_path - else: - results[result_path.name] = result_path + # Download into output/ + raw_results = sim.download_results(finished_job, output_dir=output_dir) + files = _flatten_results(raw_results) - if verbose and results: - # Find common parent directory for display - first_path = next(iter(results.values())) - print(f"Downloaded {len(results)} files to {first_path.parent}") # noqa: T201 + if verbose and files: + print(f"Downloaded {len(files)} files to {output_dir}") # noqa: T201 - return results + return RunResult(sim_dir=sim_dir, files=files, job_name=job.job_name) def print_job_summary(job) -> None: diff --git a/src/gsim/meep/__init__.py b/src/gsim/meep/__init__.py index dec415b4..3a841889 100644 --- a/src/gsim/meep/__init__.py +++ b/src/gsim/meep/__init__.py @@ -9,12 +9,13 @@ from gsim import meep sim = meep.Simulation() - sim.geometry.component = ybranch - sim.source.port = "o1" + sim.geometry(component=ybranch, z_crop="auto") + sim.materials = {"si": 3.47, "SiO2": 1.44} + sim.source(port="o1", wavelength=1.55, wavelength_span=0.01, num_freqs=11) sim.monitors = ["o1", "o2"] - sim.solver.stopping = "dft_decay" - sim.solver.max_time = 200 - result = sim.run("./meep-sim") + sim.domain(pml=1.0, margin=0.5) + sim.solver(resolution=32, simplify_tol=0.01) + result = sim.run() """ from gsim.meep.models import ( diff --git a/src/gsim/meep/meep-overview.md b/src/gsim/meep/meep-overview.md deleted file mode 100644 index b15cbabd..00000000 --- a/src/gsim/meep/meep-overview.md +++ /dev/null @@ -1,135 +0,0 @@ -# gsim.meep — Overview & Roadmap - -## Architecture - -``` -GDS + PDK ──► 3D Geometry ──► Client Viz ──► SimConfig JSON ──► Cloud Runner -(common/) (GeometryModel) (no meep) (no meep) (meep in Docker) -``` - -**Key constraint:** gsim is the client SDK — no meep dependency. MEEP runs only in Docker/cloud. - -**Upload package:** `layout.gds` + `sim_config.json` + `run_meep.py` -> Docker -> `s_parameters.csv` + `meep_debug.json` + diagnostic PNGs. - ---- - -## API: Declarative `Simulation` - -```python -from gsim import meep - -sim = meep.Simulation() -sim.geometry.component = ybranch -sim.geometry.z_crop = "auto" -sim.materials = {"si": 3.47, "SiO2": 1.44} # float shorthand -sim.source.port = "o1" -sim.source.wavelength = 1.55 -sim.source.bandwidth = 0.1 -sim.source.num_freqs = 11 -sim.monitors = ["o1", "o2"] -sim.domain.pml = 1.0 -sim.domain.margin = 0.5 -sim.solver.resolution = 32 -sim.solver.stopping = "dft_decay" -sim.solver.max_time = 200 -sim.solver.simplify_tol = 0.01 -sim.solver.save_geometry = True -sim.solver.save_fields = True -sim.plot_2d(slices="xyz") -result = sim.run("./meep-sim") -``` - -### Design principles - -- **6 typed physics objects** -- `Geometry`, `Material`, `ModeSource`, `Domain`, `FDTD` + `monitors: list[str]` -- assigned to a `Simulation` container. No ordering dependencies. -- **Attribute style or constructor style** -- `sim.source.port = "o1"` or `sim.source = ModeSource(port="o1")`. Both work via Pydantic `validate_assignment=True`. -- **Float shorthand for materials** -- `{"si": 3.47}` auto-normalizes to `Material(n=3.47)` via validator. -- **Flat stopping fields** -- `stopping` mode string + flat fields (`max_time`, `stopping_threshold`, etc.) on `FDTD`. -- **Source defines spectral window** -- `WavelengthConfig` derived from `ModeSource.wavelength/bandwidth/num_freqs`. Monitors are just port name strings. -- **JSON contract unchanged** -- `write_config()` translates new API -> existing `SimConfig` -> JSON. Runner template untouched. - ---- - -## Module Structure - -### `gsim.meep` -- Public API - -| Module | Purpose | -| ------------------- | ----------------------------------------------------------------------------------------------------------------------- | -| `__init__.py` | Exports: `Simulation`, all model classes | -| `models/api.py` | Declarative models: `Geometry`, `Material`, `ModeSource`, `Domain`, `FDTD`, `Symmetry` | -| `simulation.py` | `Simulation` container -- `write_config()`, `run()`, `validate_config()`, `plot_2d()`/`plot_3d()`, `estimate_meep_np()` | -| `viz.py` | Standalone viz helpers -- `build_geometry_model()`, `plot_2d()`, `plot_3d()` (calls `common/viz`) | -| `models/config.py` | `SimConfig` + all sub-configs (JSON serialization layer) | -| `models/results.py` | `SParameterResult` -- CSV + debug JSON + diagnostic PNGs | -| `ports.py` | `extract_port_info()` -- port center/direction/normal from gdsfactory | -| `materials.py` | `resolve_materials()` -- material names -> (n, k) via common DB | -| `script.py` | `generate_meep_script()` -- cloud runner template (string in Python) | -| `overlay.py` | `SimOverlay` + `PortOverlay` + `DielectricOverlay` -- viz metadata | - -### `gsim.common` -- Shared infrastructure - -Solver-agnostic: `LayeredComponentBase`, `GeometryModel`/`Prism`, `LayerStack`, `MaterialProperties`, `viz/` (Matplotlib 2D, PyVista/Open3D/Three.js 3D). - -### Visualization pipeline - -``` -Simulation.plot_2d() - -> viz.build_geometry_model(component, stack, domain_config) - -> LayeredComponentBase + extract_geometry_model() - -> crop_geometry_model() if stack is z-cropped - -> viz.build_overlay(gm, component, stack, domain_config) - -> overlay.build_sim_overlay() - -> common/viz/render2d.plot_prism_slices(gm, overlay=overlay) -``` - -- **Z slices** -- draw actual prism polygon outlines at the slice plane (XY view) -- **X/Y slices** -- compute Shapely line-polygon intersections for accurate cross-sections (XZ/YZ views); each intersection segment becomes a rectangle from z_base to z_top -- **Overlay** -- sim cell boundary, PML regions, port markers, dielectric backgrounds - ---- - -## Key Design Decisions - -- **GDS-file approach** -- Send raw GDS to cloud (not polygon coords in JSON). Complex geometries have 1000s of nodes; DerivedLayers need gdsfactory to resolve. -- **Port z-center = highest refractive index** -- Photonic ports center on waveguide core (not conductor layer like RF). -- **Z-crop** -- Auto-crops stack around core layer. Full UBC stack is 15um; cropped to ~1.5um (massive compute savings). -- **Port extension into PML** -- `gf.components.extend_ports()` at `write_config()` time. Original bbox stored in `SimConfig.component_bbox` for correct cell sizing. -- **Symmetries disabled for S-params** -- `add_mode_monitor` uses `use_symmetry=false` internally; `get_eigenmode_coefficients` doesn't apply `S.multiplicity()`. Source port coefficients underestimated ~2x. gplugins also never uses `mp.Mirror`. -- **`dft_min_run_time` default 100** -- Prevents false convergence before pulse traverses device. Must exceed `device_length * n_group`. -- **Stack resolved lazily** -- `_ensure_stack()` falls back to `get_stack()` (active PDK defaults) when no explicit stack kwargs are set. Same path for `write_config()` and viz. - ---- - -## TODO - -### Near-term - -- [ ] **Cloud end-to-end test** -- Full round-trip (upload -> run meep -> download results) hasn't been tested with real cloud instance. -- [ ] **Monitor z-span touches PML** -- After z-crop, layer stack z-extent exactly matches PML inner boundary. At coarse resolution, outermost monitor pixels may straddle PML. Fix: shrink monitor z-span by ~`2/resolution` inset. -- [ ] **Empty `core2` layer in config** -- UBC PDK's ebeam_y_1550 includes a `core2` entry (GDS layer 31) with zero geometry. Harmless but clutters config. Consider filtering empty layers. - -### Medium-term - -- [ ] **Custom monitors** -- Monitors are port-name strings only. Add `FieldMonitor(center, size)` and `FluxMonitor` for custom measurement locations. -- [ ] **Per-port margin/mode_index** -- `port_margin` is global in `Domain`. Could add per-port control for margin and higher-order modes. -- [ ] **Typed boundaries** -- Only PML today. Add `Periodic`, `Bloch`, `PEC`, `PMC` per-axis boundary specs. -- [ ] **Mesh refinement regions** -- Single global `resolution`. Add local refinement for thin features. -- [ ] **Dispersive materials** -- `Material(n, k)` is non-dispersive. Add Sellmeier/Lorentz/Drude support. -- [ ] **`updated_copy()` for parameter sweeps** -- Pydantic `model_copy(update={...})` makes this nearly free. - -### Deferred - -- [ ] **JSON schema v2** -- Restructure JSON groups (`"solver.wavelength"`, per-port margins, etc.) with version field for runner backward compat. Requires atomic client + runner update. -- [ ] **Port symmetries** -- gplugins-style: run fewer source ports, copy S-params between symmetric port pairs (e.g. S31=S21 for Y-branch). -- [ ] **`add_flux` + `eig_parity` path** -- Alternative to `add_mode_monitor` that works correctly with `mp.Mirror` symmetry. Needs auto-detection of correct parity from symmetry config and waveguide polarization. - ---- - -## Docker / Cloud - -- **Base:** `continuumio/miniconda3` -> conda env with `pymeep=*=mpi_mpich_*`, `gdsfactory`, `nlopt` -- **Entrypoint:** downloads input -> `mpirun -np $NP python run_meep.py` -> uploads outputs -- **`meep_np` hardcoded to 2** -- auto-compute logic exists (`total_voxels // 200k`) but `lscpu` reports host cores inside containers, not Batch-allocated vCPUs. TODO: pass Batch vCPU allocation as `MEEP_NP` env var from job submission, then restore auto-compute with proper clamping -- **Instance recommendation:** `c7a.16xlarge` (AMD EPYC Genoa, DDR5, ~460 GB/s mem BW) for typical photonics; `hpc7a.96xlarge` for large problems -- **Local testing:** `./nbs/test_meep_local.sh` (regenerate -> Docker build -> run -> results) diff --git a/src/gsim/meep/models/api.py b/src/gsim/meep/models/api.py index 2cee14cc..25861b8f 100644 --- a/src/gsim/meep/models/api.py +++ b/src/gsim/meep/models/api.py @@ -31,6 +31,12 @@ class Geometry(BaseModel): description='Z-crop mode: "auto" | layer_name | None (no crop)', ) + def __call__(self, **kwargs: Any) -> Geometry: + """Update fields in place. Returns self for chaining.""" + for k, v in kwargs.items(): + setattr(self, k, v) + return self + # --------------------------------------------------------------------------- # Material @@ -65,10 +71,11 @@ class ModeSource(BaseModel): gt=0, description="Center wavelength in um", ) - bandwidth: float = Field( + wavelength_span: float = Field( default=0.1, ge=0, - description="Measurement wavelength bandwidth in um", + description="Wavelength span of the measurement frequency grid in um. " + "Together with num_freqs, sets the spacing between monitor frequency points.", ) num_freqs: int = Field( default=11, @@ -76,6 +83,12 @@ class ModeSource(BaseModel): description="Number of frequency points", ) + def __call__(self, **kwargs: Any) -> ModeSource: + """Update fields in place. Returns self for chaining.""" + for k, v in kwargs.items(): + setattr(self, k, v) + return self + # --------------------------------------------------------------------------- # Domain @@ -118,11 +131,30 @@ class Domain(BaseModel): ge=0, description="Extend ports into PML (um). 0 = auto (margin + pml).", ) + source_port_offset: float = Field( + default=0.1, + ge=0, + description="Distance to offset source from port center into device (um). " + "Matches gplugins port_source_offset.", + ) + distance_source_to_monitors: float = Field( + default=0.2, + ge=0, + description="Distance between source and its port monitor (um). " + "The source-port monitor is placed this far past the source, " + "deeper into the device. Matches gplugins distance_source_to_monitors.", + ) symmetries: list[Symmetry] = Field( default_factory=list, description="Mirror symmetry planes. Not yet used in production runs.", ) + def __call__(self, **kwargs: Any) -> Domain: + """Update fields in place. Returns self for chaining.""" + for k, v in kwargs.items(): + setattr(self, k, v) + return self + # --------------------------------------------------------------------------- # FDTD solver @@ -137,27 +169,45 @@ class FDTD(BaseModel): resolution: int = Field(default=32, ge=4, description="Pixels per micrometer") # Stopping criteria (flat fields instead of variant classes) - stopping: Literal["fixed", "decay", "dft_decay"] = Field( - default="dft_decay", - description="Stopping mode: fixed time, field decay, or DFT convergence", + stopping: Literal["fixed", "field_decay", "dft_decay", "energy_decay"] = Field( + default="field_decay", + description=( + "Stopping mode: 'field_decay' (recommended, matches MEEP tutorials) " + "monitors a field component at a point; 'energy_decay' monitors " + "total field energy; 'dft_decay' waits for DFT convergence; " + "'fixed' runs for max_time." + ), ) max_time: float = Field( - default=200.0, gt=0, description="Max run time after sources (um/c)" + default=2000.0, gt=0, description="Max run time after sources (um/c)" ) stopping_threshold: float = Field( - default=1e-3, gt=0, lt=1, description="Decay/convergence threshold" + default=0.05, gt=0, lt=1, description="Decay/convergence threshold" ) stopping_min_time: float = Field( - default=100.0, ge=0, description="Min run time for dft_decay mode" + default=100.0, + ge=0, + description=( + "Minimum absolute sim time for dft_decay mode (not time-after-sources). " + "Must exceed pulse transit time to avoid false convergence." + ), ) stopping_component: str = Field( - default="Ey", description="Field component for decay mode" + default="Ey", description="Field component for field_decay mode" ) stopping_dt: float = Field( - default=50.0, gt=0, description="Decay measurement window for decay mode" + default=50.0, + gt=0, + description="Decay measurement window for field_decay/energy_decay modes", ) stopping_monitor_port: str | None = Field( - default=None, description="Port to monitor for decay mode" + default=None, description="Port to monitor for field_decay mode" + ) + wall_time_max: float = Field( + default=0.0, + ge=0, + description="Wall-clock time limit in seconds (0=unlimited). " + "Orthogonal safety net — stops the sim if real time exceeds this.", ) subpixel: bool = Field(default=False, description="Toggle subpixel averaging") @@ -173,6 +223,108 @@ class FDTD(BaseModel): description="Shapely simplification tolerance in um (0=off)", ) + # -- Convenience methods for stopping configuration -- + + def stop_when_energy_decayed( + self, dt: float = 50.0, decay_by: float = 0.05 + ) -> FDTD: + """Stop when total field energy in the cell decays (recommended). + + Monitors total electromagnetic energy and stops when it has decayed + by ``decay_by`` from its peak value. More robust than ``dft_decay`` + for devices where DFTs can falsely converge on near-zero fields. + + Args: + dt: Time window between energy checks (MEEP time units). + decay_by: Fractional energy decay threshold (e.g. 0.05 = 5%). + + Returns: + self (for fluent chaining). + """ + self.stopping = "energy_decay" + self.stopping_dt = dt + self.stopping_threshold = decay_by + return self + + def stop_when_dft_decayed(self, tol: float = 1e-3, min_time: float = 100.0) -> FDTD: + """Stop when all DFT monitors converge. + + Args: + tol: DFT convergence tolerance. + min_time: Minimum absolute sim time before checking convergence. + + Returns: + self (for fluent chaining). + """ + self.stopping = "dft_decay" + self.stopping_threshold = tol + self.stopping_min_time = min_time + return self + + def stop_when_fields_decayed( + self, + dt: float = 50.0, + component: str = "Ey", + decay_by: float = 0.05, + monitor_port: str | None = None, + ) -> FDTD: + """Stop when a field component decays at a point (recommended). + + Matches the standard MEEP tutorial stopping condition. Monitors + |component|² at a point and stops when it decays by ``decay_by`` + from its peak value. + + Args: + dt: Decay measurement time window. + component: Field component name (e.g. "Ey", "Hz"). + decay_by: Fractional decay threshold (e.g. 0.05 = 5%). + monitor_port: Port to monitor (None = first non-source port). + + Returns: + self (for fluent chaining). + """ + self.stopping = "field_decay" + self.stopping_dt = dt + self.stopping_component = component + self.stopping_threshold = decay_by + self.stopping_monitor_port = monitor_port + return self + + def stop_after_sources(self, time: float) -> FDTD: + """Run for a fixed sim-time after sources turn off. + + Args: + time: Run time after sources in MEEP time units (um/c). + + Returns: + self (for fluent chaining). + """ + self.stopping = "fixed" + self.max_time = time + return self + + def stop_after_walltime(self, seconds: float) -> FDTD: + """Set a wall-clock time limit (safety net). + + This is orthogonal to the sim-time stopping mode — it caps + how long the FDTD run is allowed to take in real (wall) seconds. + Combine with any other stopping method. + + Args: + seconds: Maximum wall-clock seconds for the FDTD run. + + Returns: + self (for fluent chaining). + """ + self.wall_time_max = seconds + return self + + def __call__(self, **kwargs: Any) -> FDTD: + """Update fields in place. Returns self for chaining.""" + for k, v in kwargs.items(): + setattr(self, k, v) + return self + # Diagnostics — output control for plots, fields, animations save_geometry: bool = Field(default=True) save_fields: bool = Field(default=True) diff --git a/src/gsim/meep/models/config.py b/src/gsim/meep/models/config.py index 52abcfa7..0b05576b 100644 --- a/src/gsim/meep/models/config.py +++ b/src/gsim/meep/models/config.py @@ -19,7 +19,7 @@ class SymmetryEntry(BaseModel): model_config = ConfigDict(validate_assignment=True) direction: Literal["X", "Y", "Z"] - phase: Literal[1, -1] = Field(default=1) + phase: Literal[1, -1] class DomainConfig(BaseModel): @@ -33,60 +33,80 @@ class DomainConfig(BaseModel): Cell size formula: cell_x = bbox_width + 2*(margin_xy + dpml) cell_y = bbox_height + 2*(margin_xy + dpml) - cell_z = z_extent + 2*dpml (z-margins baked into z_extent via set_z_crop) + cell_z = z_extent + 2*dpml (z-margins baked into z_extent) """ model_config = ConfigDict(validate_assignment=True) - dpml: float = Field(default=1.0, ge=0, description="PML thickness in um") + dpml: float = Field(ge=0, description="PML thickness in um") margin_xy: float = Field( - default=0.5, ge=0, description="XY margin between geometry and PML in um" + ge=0, description="XY margin between geometry and PML in um" ) margin_z_above: float = Field( - default=0.5, ge=0, description="Z margin above core kept by set_z_crop in um" + ge=0, description="Z margin above core kept by set_z_crop in um" ) margin_z_below: float = Field( - default=0.5, ge=0, description="Z margin below core kept by set_z_crop in um" + ge=0, description="Z margin below core kept by set_z_crop in um" ) port_margin: float = Field( - default=0.5, ge=0, - description="Margin on each side of port waveguide width for mode monitors in um", + description="Margin on each side of port width for monitors (um)", ) extend_ports: float = Field( - default=0.0, ge=0, description="Length to extend waveguide ports into PML in um. " "0 = auto (margin_xy + dpml).", ) + source_port_offset: float = Field( + ge=0, + description="Distance to offset source from port center into device (um).", + ) + distance_source_to_monitors: float = Field( + ge=0, + description="Distance between source and its port monitor (um). " + "Source-port monitor is placed this far past the source into the device.", + ) class StoppingConfig(BaseModel): """Controls when the MEEP simulation stops. - ``fixed`` mode runs for a fixed time after sources turn off. - ``decay`` mode monitors field decay at a point and stops when the - fields have decayed by ``threshold``, with ``max_time`` as - a numeric time cap (whichever fires first). + ``field_decay`` mode (recommended, matches MEEP tutorials) monitors + a field component at a point and stops when |component|² decays by + ``threshold`` from its peak, with ``max_time`` as a numeric time cap + (whichever fires first). + + ``energy_decay`` mode monitors total electromagnetic energy in the + cell and stops when it decays by ``threshold`` from its peak. + ``dft_decay`` mode monitors convergence of all DFT monitors and - stops when they stabilize, with built-in min/max time bounds. - Best for S-parameter extraction. + stops when they stabilize. ``dft_min_run_time`` is an *absolute* + sim time (not time-after-sources) — with a broadband source turning + off at ~t=78, a min_run_time=100 starts checking at t=100, only ~22 + time units after the source ends. + + ``fixed`` mode runs for a fixed time after sources turn off. """ model_config = ConfigDict(validate_assignment=True) - mode: Literal["fixed", "decay", "dft_decay"] = Field(default="fixed") - max_time: float = Field(default=100.0, gt=0, serialization_alias="run_after_sources") - decay_dt: float = Field(default=50.0, gt=0) - decay_component: str = Field(default="Ey") - threshold: float = Field(default=1e-3, gt=0, lt=1, serialization_alias="decay_by") + mode: Literal["fixed", "field_decay", "dft_decay", "energy_decay"] + max_time: float = Field(gt=0, serialization_alias="run_after_sources") + decay_dt: float = Field(gt=0) + decay_component: str + threshold: float = Field(gt=0, lt=1, serialization_alias="decay_by") decay_monitor_port: str | None = Field(default=None) dft_min_run_time: float = Field( - default=100, ge=0, - description="Minimum run time after sources for dft_decay mode. " - "Must exceed pulse transit time through the device to avoid " - "false convergence on near-zero fields at output ports.", + description="Minimum absolute sim time for dft_decay mode (not " + "time-after-sources). Must exceed pulse transit time through the " + "device to avoid false convergence on near-zero fields.", + ) + wall_time_max: float = Field( + default=0.0, + ge=0, + description="Wall-clock time limit in seconds (0=unlimited). " + "Orthogonal safety net for all stopping modes.", ) @@ -104,7 +124,9 @@ class SourceConfig(BaseModel): bandwidth: float | None = Field( default=None, - description="Source Gaussian bandwidth in wavelength um. None = auto (~3x monitor bw).", + description=( + "Source Gaussian bandwidth in wavelength um. None = auto (~3x monitor bw)." + ), ) port: str | None = Field( default=None, @@ -148,11 +170,9 @@ class WavelengthConfig(BaseModel): model_config = ConfigDict(validate_assignment=True) - wavelength: float = Field(default=1.55, gt=0, description="Center wavelength in um") - bandwidth: float = Field( - default=0.1, ge=0, description="Wavelength bandwidth in um" - ) - num_freqs: int = Field(default=11, ge=1, description="Number of frequency points") + wavelength: float = Field(gt=0, description="Center wavelength in um") + bandwidth: float = Field(ge=0, description="Wavelength bandwidth in um") + num_freqs: int = Field(ge=1, description="Number of frequency points") @computed_field # type: ignore[prop-decorator] @property @@ -176,9 +196,7 @@ class ResolutionConfig(BaseModel): model_config = ConfigDict(validate_assignment=True) - pixels_per_um: int = Field( - default=32, ge=4, description="Grid pixels per micrometer" - ) + pixels_per_um: int = Field(ge=4, description="Grid pixels per micrometer") @classmethod def coarse(cls) -> ResolutionConfig: @@ -208,15 +226,13 @@ class AccuracyConfig(BaseModel): model_config = ConfigDict(validate_assignment=True) - eps_averaging: bool = Field(default=False, description="Toggle subpixel averaging") + eps_averaging: bool = Field(description="Toggle subpixel averaging") subpixel_maxeval: int = Field( - default=0, ge=0, description="Cap on integration evaluations (0=unlimited)" - ) - subpixel_tol: float = Field( - default=1e-4, gt=0, description="Subpixel integration tolerance" + ge=0, description="Cap on integration evaluations (0=unlimited)" ) + subpixel_tol: float = Field(gt=0, description="Subpixel integration tolerance") simplify_tol: float = Field( - default=0.0, ge=0, description="Shapely simplification tolerance in um (0=off)" + ge=0, description="Shapely simplification tolerance in um (0=off)" ) @@ -225,24 +241,21 @@ class DiagnosticsConfig(BaseModel): model_config = ConfigDict(validate_assignment=True) - save_geometry: bool = Field(default=True, description="Pre-run geometry cross-section plots") - save_fields: bool = Field(default=True, description="Post-run field snapshot") - save_epsilon_raw: bool = Field( - default=False, description="Raw epsilon .npy (advanced)" - ) + save_geometry: bool = Field(description="Pre-run geometry cross-section plots") + save_fields: bool = Field(description="Post-run field snapshot") + save_epsilon_raw: bool = Field(description="Raw epsilon .npy (advanced)") save_animation: bool = Field( - default=False, description="Field animation MP4 (heavy, needs ffmpeg)" + description="Field animation MP4 (heavy, needs ffmpeg)" ) animation_interval: float = Field( - default=0.5, gt=0, + gt=0, description="MEEP time units between animation frames", ) preview_only: bool = Field( - default=False, description="Init sim and save geometry diagnostics, skip FDTD run", ) verbose_interval: float = Field( - default=0, ge=0, + ge=0, description="MEEP time units between progress prints (0=off)", ) @@ -297,34 +310,30 @@ class SimConfig(BaseModel): model_config = ConfigDict(validate_assignment=True) - gds_filename: str = Field( - default="layout.gds", description="GDS file with 2D layout" - ) + gds_filename: str = Field(description="GDS file with 2D layout") component_bbox: list[float] | None = Field( default=None, - description="Original component bbox [xmin, ymin, xmax, ymax] before port extension.", + description=( + "Original component bbox [xmin, ymin, xmax, ymax] before port extension." + ), ) - layer_stack: list[LayerStackEntry] = Field(default_factory=list) - dielectrics: list[dict[str, Any]] = Field(default_factory=list) - ports: list[PortData] = Field(default_factory=list) - materials: dict[str, MaterialData] = Field(default_factory=dict) - wavelength: WavelengthConfig = Field( - default_factory=WavelengthConfig, serialization_alias="fdtd" - ) - source: SourceConfig = Field(default_factory=SourceConfig) - stopping: StoppingConfig = Field(default_factory=StoppingConfig) - resolution: ResolutionConfig = Field(default_factory=ResolutionConfig) - domain: DomainConfig = Field(default_factory=DomainConfig) - accuracy: AccuracyConfig = Field(default_factory=AccuracyConfig) + layer_stack: list[LayerStackEntry] + dielectrics: list[dict[str, Any]] + ports: list[PortData] + materials: dict[str, MaterialData] + wavelength: WavelengthConfig = Field(serialization_alias="fdtd") + source: SourceConfig + stopping: StoppingConfig + resolution: ResolutionConfig + domain: DomainConfig + accuracy: AccuracyConfig verbose_interval: float = Field( - default=0, ge=0, description="MEEP time units between progress prints (0=off)" + ge=0, description="MEEP time units between progress prints (0=off)" ) - diagnostics: DiagnosticsConfig = Field(default_factory=DiagnosticsConfig) - symmetries: list[SymmetryEntry] = Field(default_factory=list) + diagnostics: DiagnosticsConfig + symmetries: list[SymmetryEntry] split_chunks_evenly: bool = Field(default=False) - meep_np: int = Field( - default=1, ge=1, description="Recommended MPI process count" - ) + meep_np: int = Field(default=1, ge=1, description="Recommended MPI process count") def to_json(self, path: str | Path) -> Path: """Write config to JSON file. diff --git a/src/gsim/meep/script.py b/src/gsim/meep/script.py index 3d08b362..0a9b99b8 100644 --- a/src/gsim/meep/script.py +++ b/src/gsim/meep/script.py @@ -59,7 +59,7 @@ def build_materials(config): materials = {} for name, props in config["materials"].items(): n = props["refractive_index"] - k = props.get("extinction_coeff", 0.0) + k = props["extinction_coeff"] if k > 0: materials[name] = mp.Medium(index=n, D_conductivity=2 * cmath.pi * k / n) else: @@ -202,7 +202,7 @@ def build_background_slabs(config, materials): geometry list so that patterned prisms (added later) take precedence. """ slabs = [] - for diel in sorted(config.get("dielectrics", []), key=lambda d: d["zmin"]): + for diel in sorted(config["dielectrics"], key=lambda d: d["zmin"]): mat = materials.get(diel["material"]) if mat is None: continue @@ -224,9 +224,9 @@ def build_symmetries(config): """Build MEEP symmetry objects from config.""" direction_map = {"X": mp.X, "Y": mp.Y, "Z": mp.Z} symmetries = [] - for sym in config.get("symmetries", []): + for sym in config["symmetries"]: symmetries.append( - mp.Mirror(direction_map[sym["direction"]], phase=sym.get("phase", 1)) + mp.Mirror(direction_map[sym["direction"]], phase=sym["phase"]) ) return symmetries @@ -238,6 +238,37 @@ def build_symmetries(config): } +def _make_time_cap(cap): + """Wrap a numeric time cap as a callable for until_after_sources lists. + + When until_after_sources receives a list, every element must be callable. + This wraps a float (time units) into a function that returns True once + ``cap`` time units have elapsed since the first call (i.e. after sources + turn off). + """ + _t0 = [None] + def _check(sim_obj): + if _t0[0] is None: + _t0[0] = sim_obj.meep_time() + return (sim_obj.meep_time() - _t0[0]) >= cap + return _check + + +def _make_wall_time_cap(wall_seconds): + """Wrap a wall-clock time limit as a callable for until_after_sources lists. + + Returns True once ``wall_seconds`` of real (wall) time have elapsed + since the first call. This provides a safety net orthogonal to any + sim-time stopping condition. + """ + _deadline = [None] + def _check(sim_obj): + if _deadline[0] is None: + _deadline[0] = time.time() + wall_seconds + return time.time() >= _deadline[0] + return _check + + def resolve_decay_monitor_point(config): """Return center mp.Vector3 for decay monitoring. @@ -245,9 +276,9 @@ def resolve_decay_monitor_point(config): otherwise picks the first non-source port. Falls back to the first port if all are sources. """ - stopping = config.get("stopping", {}) + stopping = config["stopping"] target_name = stopping.get("decay_monitor_port") - ports = config.get("ports", []) + ports = config["ports"] if target_name: for p in ports: @@ -256,7 +287,7 @@ def resolve_decay_monitor_point(config): # Fallback: first non-source port, or first port for p in ports: - if not p.get("is_source", False): + if not p["is_source"]: return mp.Vector3(*p["center"]) if ports: return mp.Vector3(*ports[0]["center"]) @@ -272,11 +303,11 @@ def build_geometry(config, materials): 2. Extrude to 3D as mp.Prism with correct z-range and material 3. Handle polygon holes via Delaunay triangulation """ - gds_filename = config.get("gds_filename", "layout.gds") + gds_filename = config["gds_filename"] component = load_gds_component(gds_filename) - accuracy = config.get("accuracy", {}) - simplify_tol = accuracy.get("simplify_tol", 0.0) + accuracy = config["accuracy"] + simplify_tol = accuracy["simplify_tol"] geometry = [] total_vertices = 0 @@ -288,7 +319,7 @@ def build_geometry(config, materials): zmax = layer_entry["zmax"] height = zmax - zmin gds_layer = layer_entry["gds_layer"] - sidewall_angle_deg = layer_entry.get("sidewall_angle", 0.0) + sidewall_angle_deg = layer_entry["sidewall_angle"] sw_rad = math.radians(sidewall_angle_deg) if sidewall_angle_deg else 0 if height <= 0: @@ -344,22 +375,35 @@ def get_port_z_span(config): def build_sources(config): - """Build MEEP source from config port data.""" + """Build MEEP source from config port data. + + The source is offset from the port center by ``source_port_offset`` + along the propagation direction (into the device). This separates + the soft source from the port monitor so eigenmode coefficients + measure the true incident amplitude rather than half of it. + """ fdtd = config["fdtd"] fcen = fdtd["fcen"] df = fdtd["df"] fwidth = config["source"]["fwidth"] z_span = get_port_z_span(config) - port_margin = config.get("domain", {}).get("port_margin", 0.5) + port_margin = config["domain"]["port_margin"] + source_port_offset = config["domain"].get("source_port_offset", 0.1) sources = [] for port in config["ports"]: if not port["is_source"]: continue - center = mp.Vector3(*port["center"]) + # Offset source center into the device along propagation direction + center_list = list(port["center"]) normal_axis = port["normal_axis"] direction = port["direction"] + if direction == "+": + center_list[normal_axis] += source_port_offset + else: + center_list[normal_axis] -= source_port_offset + center = mp.Vector3(*center_list) size = [0, 0, 0] transverse_axis = 1 - normal_axis @@ -390,18 +434,45 @@ def build_sources(config): def build_monitors(config, sim): - """Build mode monitors at all ports and return flux regions.""" + """Build mode monitors at all ports and return flux regions. + + The source-port monitor is offset further into the device (past + the source) by ``source_port_offset + distance_source_to_monitors`` + so the forward-going mode from the source passes through it at + full amplitude. This matches the gplugins approach. + """ fdtd = config["fdtd"] fcen = fdtd["fcen"] df = fdtd["df"] nfreq = fdtd["num_freqs"] z_span = get_port_z_span(config) - port_margin = config.get("domain", {}).get("port_margin", 0.5) + port_margin = config["domain"]["port_margin"] + source_port_offset = config["domain"].get("source_port_offset", 0.1) + distance_source_to_monitors = config["domain"].get( + "distance_source_to_monitors", 0.2 + ) monitors = {} for port in config["ports"]: - center = mp.Vector3(*port["center"]) + center_list = list(port["center"]) normal_axis = port["normal_axis"] + direction = port["direction"] + + # All monitors shift inward from port center (matches gplugins): + # source-port monitor: source_port_offset + distance_source_to_monitors + # non-source monitors: source_port_offset + if port["is_source"]: + offset = source_port_offset + distance_source_to_monitors + else: + offset = source_port_offset + + if offset > 0: + if direction == "+": + center_list[normal_axis] += offset + else: + center_list[normal_axis] -= offset + + center = mp.Vector3(*center_list) size = [0, 0, 0] transverse_axis = 1 - normal_axis @@ -602,9 +673,9 @@ def save_debug_log(config, s_params, debug_data, wall_seconds=0.0, if not mp.am_master(): return fdtd = config["fdtd"] - domain = config.get("domain", {}) - resolution = config.get("resolution", {}).get("pixels_per_um", 0) - stopping = config.get("stopping", {}) + domain = config["domain"] + resolution = config["resolution"]["pixels_per_um"] + stopping = config["stopping"] meep_time = 0.0 timesteps = 0 @@ -622,7 +693,7 @@ def save_debug_log(config, s_params, debug_data, wall_seconds=0.0, "meep_time": meep_time, "timesteps": timesteps, "wall_seconds": wall_seconds, - "stopping_mode": stopping.get("mode", "fixed"), + "stopping_mode": stopping["mode"], }, "eigenmode_info": debug_data.get("eigenmode_info", {}), "incident_coefficients": debug_data.get("incident_coefficients", {}), @@ -648,8 +719,8 @@ def save_geometry_diagnostics(sim, config, cell_center): # plot2D is collective — all ranks call it, only master saves pass - domain = config.get("domain", {}) - dpml = domain.get("dpml", 1.0) + domain = config["domain"] + dpml = domain["dpml"] z_min = min(l["zmin"] for l in config["layer_stack"]) z_max = max(l["zmax"] for l in config["layer_stack"]) z_core = (z_min + z_max) / 2 @@ -929,7 +1000,7 @@ def main(): # Compute simulation cell from component bounds + layer z-range # Use original component bbox if available (port extension changes GDS bbox) - component_bbox = config.get("component_bbox") + component_bbox = config["component_bbox"] if component_bbox is not None: bbox_left, bbox_bottom, bbox_right, bbox_top = component_bbox else: @@ -940,9 +1011,9 @@ def main(): z_min = min(l["zmin"] for l in config["layer_stack"]) z_max = max(l["zmax"] for l in config["layer_stack"]) - domain = config.get("domain", {}) - dpml = domain.get("dpml", 1.0) - margin_xy = domain.get("margin_xy", 0.5) + domain = config["domain"] + dpml = domain["dpml"] + margin_xy = domain["margin_xy"] # XY: margin_xy is gap between geometry bbox and PML cell_x = (bbox_right - bbox_left) + 2 * (margin_xy + dpml) @@ -960,12 +1031,12 @@ def main(): logger.info("PML: %.2f um, margin_xy: %.2f", dpml, margin_xy) logger.info("Resolution: %s pixels/um", resolution) - accuracy = config.get("accuracy", {}) - diagnostics = config.get("diagnostics", {}) - preview_only = diagnostics.get("preview_only", False) + accuracy = config["accuracy"] + diagnostics = config["diagnostics"] + preview_only = diagnostics["preview_only"] # In preview mode, skip expensive subpixel averaging - eps_avg = False if preview_only else accuracy.get("eps_averaging", True) + eps_avg = False if preview_only else accuracy["eps_averaging"] # Symmetries are NOT used for S-parameter extraction runs. # MEEP's get_eigenmode_coefficients with add_mode_monitor produces @@ -974,7 +1045,7 @@ def main(): # consistent with gplugins which also never uses mp.Mirror for # S-parameter extraction. Symmetries are only applied in preview-only # mode (geometry validation, no FDTD/S-params). - cfg_symmetries = config.get("symmetries", []) + cfg_symmetries = config["symmetries"] if cfg_symmetries and not preview_only: logger.info("Symmetries present in config but IGNORED for S-parameter " "extraction (causes incorrect eigenmode normalization). " @@ -989,21 +1060,21 @@ def main(): resolution=resolution, boundary_layers=[mp.PML(dpml)], symmetries=use_symmetries, - split_chunks_evenly=config.get("split_chunks_evenly", False), + split_chunks_evenly=config["split_chunks_evenly"], eps_averaging=eps_avg, ) - spx_maxeval = accuracy.get("subpixel_maxeval", 0) + spx_maxeval = accuracy["subpixel_maxeval"] if spx_maxeval > 0: sim_kwargs["subpixel_maxeval"] = spx_maxeval - spx_tol = accuracy.get("subpixel_tol", 1e-4) + spx_tol = accuracy["subpixel_tol"] if spx_tol != 1e-4: sim_kwargs["subpixel_tol"] = spx_tol sim = mp.Simulation(**sim_kwargs) # --- Diagnostics & preview mode --- - diag_geometry = diagnostics.get("save_geometry", True) - diag_fields = diagnostics.get("save_fields", True) - diag_epsilon = diagnostics.get("save_epsilon_raw", False) + diag_geometry = diagnostics["save_geometry"] + diag_fields = diagnostics["save_fields"] + diag_epsilon = diagnostics["save_epsilon_raw"] if diag_geometry or diag_epsilon or preview_only: logger.info("Initializing simulation for diagnostics...") @@ -1023,12 +1094,12 @@ def main(): logger.info("Building monitors...") monitors = build_monitors(config, sim) - stopping = config.get("stopping", {}) - run_after = stopping.get("run_after_sources", 100) + stopping = config["stopping"] + run_after = stopping["run_after_sources"] # Build verbose step functions step_funcs = [] - verbose_interval = config.get("verbose_interval", 0) + verbose_interval = config["verbose_interval"] if verbose_interval > 0: _wall_start = time.time() def _verbose_print(sim_obj): @@ -1037,8 +1108,8 @@ def _verbose_print(sim_obj): step_funcs.append(mp.at_every(verbose_interval, _verbose_print)) # Animation field capture step function (raw data, no plotting yet) - diag_animation = diagnostics.get("save_animation", False) - animation_interval = diagnostics.get("animation_interval", 0.5) + diag_animation = diagnostics["save_animation"] + animation_interval = diagnostics["animation_interval"] _frame_counter = [0] # mutable container for closure _anim_plane = None @@ -1062,43 +1133,78 @@ def _capture_frame(sim_obj): animation_interval, ) - stop_mode = stopping.get("mode", "fixed") + stop_mode = stopping["mode"] + wall_time_max = stopping.get("wall_time_max", 0) + if wall_time_max > 0: + logger.info("Wall-clock safety net: %.0f seconds", wall_time_max) wall_start = time.time() if stop_mode == "dft_decay": - decay_by = stopping.get("decay_by", 1e-3) - min_time = stopping.get("dft_min_run_time", 100) + decay_by = stopping["decay_by"] + min_time = stopping["dft_min_run_time"] logger.info( "Running simulation (dft_decay mode: tol=%s, " "min=%.1f, max=%.1f)...", decay_by, min_time, run_after, ) - sim.run( - *step_funcs, - until_after_sources=mp.stop_when_dft_decayed( - tol=decay_by, - minimum_run_time=min_time, - maximum_run_time=run_after, - ), + dft_fn = mp.stop_when_dft_decayed( + tol=decay_by, + minimum_run_time=min_time, + maximum_run_time=run_after, ) - elif stop_mode == "decay": - dt = stopping.get("decay_dt", 50.0) - comp_name = stopping.get("decay_component", "Ey") + if wall_time_max > 0: + conds = [dft_fn, _make_wall_time_cap(wall_time_max)] + sim.run(*step_funcs, until_after_sources=conds) + else: + sim.run(*step_funcs, until_after_sources=dft_fn) + elif stop_mode == "energy_decay": + dt = stopping["decay_dt"] + decay_by = stopping["decay_by"] + logger.info( + "Running simulation (energy_decay mode: dt=%s, " + "decay_by=%s, cap=%.1f)...", + dt, decay_by, run_after, + ) + energy_fn = mp.stop_when_energy_decayed(dt=dt, decay_by=decay_by) + conds = [energy_fn, _make_time_cap(run_after)] + if wall_time_max > 0: + conds.append(_make_wall_time_cap(wall_time_max)) + sim.run(*step_funcs, until_after_sources=conds) + elif stop_mode == "field_decay": + dt = stopping["decay_dt"] + comp_name = stopping["decay_component"] comp = _COMPONENT_MAP.get(comp_name, mp.Ey) - decay_by = stopping.get("decay_by", 1e-3) + decay_by = stopping["decay_by"] monitor_pt = resolve_decay_monitor_point(config) - logger.info("Running simulation (decay mode: component=%s, dt=%s, " + logger.info("Running simulation (field_decay mode: component=%s, dt=%s, " "decay_by=%s, cap=%.1f)...", comp_name, dt, decay_by, run_after) # Decay condition + numeric time cap (list = OR logic, first wins) decay_fn = mp.stop_when_fields_decayed(dt, comp, monitor_pt, decay_by) - sim.run(*step_funcs, until_after_sources=[decay_fn, run_after]) + conds = [decay_fn, _make_time_cap(run_after)] + if wall_time_max > 0: + conds.append(_make_wall_time_cap(wall_time_max)) + sim.run(*step_funcs, until_after_sources=conds) else: logger.info("Running simulation (until_after_sources=%.1f)...", run_after) - sim.run(*step_funcs, until_after_sources=run_after) + if wall_time_max > 0: + conds = [ + _make_time_cap(run_after), + _make_wall_time_cap(wall_time_max), + ] + sim.run(*step_funcs, until_after_sources=conds) + else: + sim.run(*step_funcs, until_after_sources=run_after) wall_seconds = time.time() - wall_start + if wall_time_max > 0 and wall_seconds >= wall_time_max * 0.95: + logger.warning( + "Wall-clock limit likely triggered (%.1fs elapsed, limit %.0fs). " + "Results may be incomplete — consider increasing wall_time_max or " + "using a coarser resolution.", + wall_seconds, wall_time_max, + ) if diag_fields: save_field_snapshot(sim, config, cell_center) diff --git a/src/gsim/meep/simulation.py b/src/gsim/meep/simulation.py index ea938ae8..afb84126 100644 --- a/src/gsim/meep/simulation.py +++ b/src/gsim/meep/simulation.py @@ -98,7 +98,7 @@ class Simulation(BaseModel): sim.monitors = ["o1", "o2"] sim.solver.stopping = "dft_decay" sim.solver.max_time = 200 - result = sim.run("./meep-sim") + result = sim.run() # creates sim-data-{job_name}/ in CWD """ model_config = ConfigDict( @@ -308,7 +308,7 @@ def _wavelength_config(self) -> Any: return WavelengthConfig( wavelength=self.source.wavelength, - bandwidth=self.source.bandwidth, + bandwidth=self.source.wavelength_span, num_freqs=self.source.num_freqs, ) @@ -334,6 +334,7 @@ def _stopping_config(self) -> Any: decay_component=s.stopping_component, decay_dt=s.stopping_dt, decay_monitor_port=s.stopping_monitor_port, + wall_time_max=s.wall_time_max, ) def _domain_config(self) -> Any: @@ -347,6 +348,8 @@ def _domain_config(self) -> Any: margin_z_below=self.domain.margin_z_below, port_margin=self.domain.port_margin, extend_ports=self.domain.extend_ports, + source_port_offset=self.domain.source_port_offset, + distance_source_to_monitors=self.domain.distance_source_to_monitors, ) def _resolution_config(self) -> Any: @@ -617,12 +620,21 @@ def write_config(self, output_dir: str | Path) -> Path: # run # ------------------------------------------------------------------------- - def run(self, output_dir: str | Path | None = None, *, verbose: bool = True) -> Any: + def run( + self, + parent_dir: str | Path | None = None, + *, + verbose: bool = True, + ) -> Any: """Run MEEP simulation on the cloud. + Config files are written to a temporary directory, uploaded, then + moved into a structured ``sim-data-{job_name}/input/`` directory. + Results are downloaded to ``sim-data-{job_name}/output/``. + Args: - output_dir: Directory to write config and download results. - If None, a temporary directory is created automatically. + parent_dir: Where to create the sim directory. + Defaults to the current working directory. verbose: Print progress info. Returns: @@ -633,18 +645,23 @@ def run(self, output_dir: str | Path | None = None, *, verbose: bool = True) -> from gsim.gcloud import run_simulation from gsim.meep.models.results import SParameterResult - if output_dir is None: - output_dir = Path(tempfile.mkdtemp(prefix="meep_")) - - output_dir = self.write_config(output_dir) + tmp = Path(tempfile.mkdtemp(prefix="meep_")) + try: + self.write_config(tmp) + result = run_simulation( + config_dir=tmp, + job_type="meep", + verbose=verbose, + parent_dir=parent_dir, + ) + except Exception: + # Clean up temp dir on failure (run_simulation moves it on success) + import shutil - results = run_simulation( - output_dir, - job_type="meep", - verbose=verbose, - ) + shutil.rmtree(tmp, ignore_errors=True) + raise - csv_path = results.get("s_parameters.csv") + csv_path = result.files.get("s_parameters.csv") if csv_path is not None: return SParameterResult.from_csv(csv_path) diff --git a/src/gsim/palace/__init__.py b/src/gsim/palace/__init__.py index 6b242127..8de7ce5f 100644 --- a/src/gsim/palace/__init__.py +++ b/src/gsim/palace/__init__.py @@ -22,7 +22,7 @@ # Generate mesh and run sim.mesh("./sim", preset="fine") - results = sim.simulate() + results = sim.run() """ from __future__ import annotations diff --git a/src/gsim/palace/driven.py b/src/gsim/palace/driven.py index 52c5860d..583f7e86 100644 --- a/src/gsim/palace/driven.py +++ b/src/gsim/palace/driven.py @@ -46,7 +46,7 @@ class DrivenSim(PalaceSimMixin, BaseModel): >>> sim.add_cpw_port("P3", "P4", layer="topmetal2", length=5.0) >>> sim.set_driven(fmin=1e9, fmax=100e9, num_points=40) >>> sim.mesh("./sim", preset="default") - >>> results = sim.simulate() + >>> results = sim.run() Attributes: geometry: Wrapped gdsfactory Component (from common) @@ -663,7 +663,7 @@ def write_config(self) -> Path: # Simulation # ------------------------------------------------------------------------- - def simulate( + def run( self, *, verbose: bool = True, @@ -672,6 +672,10 @@ def simulate( Requires mesh() and write_config() to be called first. + Config files are uploaded, then moved into a structured + ``sim-data-{job_name}/input/`` directory. Results are downloaded + to ``sim-data-{job_name}/output/``. + Args: verbose: Print progress messages @@ -684,7 +688,7 @@ def simulate( RuntimeError: If simulation fails Example: - >>> results = sim.simulate() + >>> results = sim.run() >>> print(f"S-params saved to: {results['port-S.csv']}") """ from gsim.gcloud import run_simulation @@ -692,15 +696,14 @@ def simulate( if self._output_dir is None: raise ValueError("Output directory not set. Call set_output_dir() first.") - output_dir = self._output_dir - - return run_simulation( - output_dir, + result = run_simulation( + config_dir=self._output_dir, job_type="palace", verbose=verbose, ) + return result.files - def simulate_local( + def run_local( self, *, verbose: bool = True, diff --git a/src/gsim/palace/eigenmode.py b/src/gsim/palace/eigenmode.py index aeb123fa..e75394dd 100644 --- a/src/gsim/palace/eigenmode.py +++ b/src/gsim/palace/eigenmode.py @@ -45,7 +45,7 @@ class EigenmodeSim(PalaceSimMixin, BaseModel): >>> sim.set_eigenmode(num_modes=10, target=50e9) >>> sim.set_output_dir("./sim") >>> sim.mesh(preset="default") - >>> results = sim.simulate() + >>> results = sim.run() Attributes: geometry: Wrapped gdsfactory Component (from common) @@ -530,7 +530,7 @@ def mesh( # Simulation # ------------------------------------------------------------------------- - def simulate( + def run( self, output_dir: str | Path | None = None, *, diff --git a/src/gsim/palace/electrostatic.py b/src/gsim/palace/electrostatic.py index 4e4f8ab0..0f1d0626 100644 --- a/src/gsim/palace/electrostatic.py +++ b/src/gsim/palace/electrostatic.py @@ -46,7 +46,7 @@ class ElectrostaticSim(PalaceSimMixin, BaseModel): >>> sim.set_electrostatic() >>> sim.set_output_dir("./sim") >>> sim.mesh(preset="default") - >>> results = sim.simulate() + >>> results = sim.run() Attributes: geometry: Wrapped gdsfactory Component (from common) @@ -371,7 +371,7 @@ def mesh( # Simulation # ------------------------------------------------------------------------- - def simulate( + def run( self, output_dir: str | Path | None = None, *, diff --git a/tests/meep/test_meep_models.py b/tests/meep/test_meep_models.py index a3ea3e98..61a5c8a6 100644 --- a/tests/meep/test_meep_models.py +++ b/tests/meep/test_meep_models.py @@ -32,25 +32,23 @@ class TestWavelengthConfig: """Test WavelengthConfig frequency/wavelength conversion.""" - def test_defaults(self): - cfg = WavelengthConfig() - assert cfg.wavelength == 1.55 - assert cfg.bandwidth == 0.1 - assert cfg.num_freqs == 11 + def test_requires_all_fields(self): + with pytest.raises(ValidationError): + WavelengthConfig() def test_fcen(self): - cfg = WavelengthConfig(wavelength=1.55) + cfg = WavelengthConfig(wavelength=1.55, bandwidth=0.1, num_freqs=11) assert abs(cfg.fcen - 1.0 / 1.55) < 1e-10 def test_df(self): - cfg = WavelengthConfig(wavelength=1.55, bandwidth=0.1) + cfg = WavelengthConfig(wavelength=1.55, bandwidth=0.1, num_freqs=11) wl_min = 1.55 - 0.05 wl_max = 1.55 + 0.05 expected_df = 1.0 / wl_min - 1.0 / wl_max assert abs(cfg.df - expected_df) < 1e-10 def test_model_dump(self): - cfg = WavelengthConfig() + cfg = WavelengthConfig(wavelength=1.55, bandwidth=0.1, num_freqs=11) d = cfg.model_dump() assert "wavelength" in d assert "fcen" in d @@ -81,7 +79,7 @@ def test_custom(self): assert cfg.pixels_per_um == 48 def test_model_dump(self): - cfg = ResolutionConfig() + cfg = ResolutionConfig(pixels_per_um=32) d = cfg.model_dump() assert d["pixels_per_um"] == 32 @@ -90,12 +88,25 @@ class TestSimConfig: """Test SimConfig serialization.""" def test_to_json(self, tmp_path): - wl_cfg = WavelengthConfig() + wl_cfg = WavelengthConfig(wavelength=1.55, bandwidth=0.1, num_freqs=11) fwidth = SourceConfig().compute_fwidth(wl_cfg.fcen, wl_cfg.df) source_cfg = SourceConfig(fwidth=fwidth) - stopping_cfg = StoppingConfig(mode="dft_decay", max_time=200.0) + stopping_cfg = StoppingConfig( + mode="dft_decay", + max_time=200.0, + decay_dt=50.0, + decay_component="Ey", + threshold=0.05, + dft_min_run_time=100, + ) + from gsim.meep.models.config import ( + AccuracyConfig, + DiagnosticsConfig, + ) + cfg = SimConfig( gds_filename="layout.gds", + verbose_interval=0, layer_stack=[ # ty: ignore[invalid-argument-type] { "layer_name": "core", @@ -121,7 +132,34 @@ def test_to_json(self, tmp_path): wavelength=wl_cfg, source=source_cfg, stopping=stopping_cfg, - resolution=ResolutionConfig(), + resolution=ResolutionConfig(pixels_per_um=32), + domain=DomainConfig( + dpml=1.0, + margin_xy=0.5, + margin_z_above=0.5, + margin_z_below=0.5, + port_margin=0.5, + extend_ports=0.0, + source_port_offset=0.1, + distance_source_to_monitors=0.2, + ), + accuracy=AccuracyConfig( + eps_averaging=False, + subpixel_maxeval=0, + subpixel_tol=1e-4, + simplify_tol=0.0, + ), + diagnostics=DiagnosticsConfig( + save_geometry=True, + save_fields=True, + save_epsilon_raw=False, + save_animation=False, + animation_interval=0.5, + preview_only=False, + verbose_interval=0, + ), + dielectrics=[], + symmetries=[], ) path = tmp_path / "config.json" cfg.to_json(path) @@ -519,18 +557,20 @@ def test_import_all_public_api(self): class TestDomainConfig: """Test DomainConfig model.""" - def test_defaults(self): - cfg = DomainConfig() - assert cfg.dpml == 1.0 - assert cfg.margin_xy == 0.5 - assert cfg.margin_z_above == 0.5 - assert cfg.margin_z_below == 0.5 - assert cfg.port_margin == 0.5 - assert cfg.extend_ports == 0.0 + def test_requires_all_fields(self): + with pytest.raises(ValidationError): + DomainConfig() def test_custom(self): cfg = DomainConfig( - dpml=0.5, margin_xy=0.2, margin_z_above=0.3, margin_z_below=0.4 + dpml=0.5, + margin_xy=0.2, + margin_z_above=0.3, + margin_z_below=0.4, + port_margin=0.5, + extend_ports=0.0, + source_port_offset=0.1, + distance_source_to_monitors=0.2, ) assert cfg.dpml == 0.5 assert cfg.margin_xy == 0.2 @@ -538,17 +578,42 @@ def test_custom(self): assert cfg.margin_z_below == 0.4 def test_extend_ports_custom(self): - cfg = DomainConfig(extend_ports=2.5) + cfg = DomainConfig( + dpml=1.0, + margin_xy=0.5, + margin_z_above=0.5, + margin_z_below=0.5, + port_margin=0.5, + extend_ports=2.5, + source_port_offset=0.1, + distance_source_to_monitors=0.2, + ) assert cfg.extend_ports == 2.5 def test_extend_ports_serialization(self): - cfg = DomainConfig(extend_ports=3.0) + cfg = DomainConfig( + dpml=1.0, + margin_xy=0.5, + margin_z_above=0.5, + margin_z_below=0.5, + port_margin=0.5, + extend_ports=3.0, + source_port_offset=0.1, + distance_source_to_monitors=0.2, + ) d = cfg.model_dump() assert d["extend_ports"] == 3.0 def test_model_dump(self): cfg = DomainConfig( - dpml=2.0, margin_xy=0.5, margin_z_above=1.0, margin_z_below=1.5 + dpml=2.0, + margin_xy=0.5, + margin_z_above=1.0, + margin_z_below=1.5, + port_margin=0.5, + extend_ports=0.0, + source_port_offset=0.1, + distance_source_to_monitors=0.2, ) d = cfg.model_dump() assert d["dpml"] == 2.0 @@ -560,23 +625,73 @@ def test_model_dump(self): class TestSimConfigComponentBbox: """Test SimConfig.component_bbox field.""" - def test_default_none(self): - cfg = SimConfig() + @pytest.fixture + def _sim_kwargs(self): + from gsim.meep.models.config import AccuracyConfig, DiagnosticsConfig + + return dict( + gds_filename="layout.gds", + verbose_interval=0, + layer_stack=[], + dielectrics=[], + ports=[], + materials={}, + wavelength=WavelengthConfig(wavelength=1.55, bandwidth=0.1, num_freqs=11), + source=SourceConfig(), + stopping=StoppingConfig( + mode="fixed", + max_time=100.0, + decay_dt=50.0, + decay_component="Ey", + threshold=0.05, + dft_min_run_time=100, + ), + resolution=ResolutionConfig(pixels_per_um=32), + domain=DomainConfig( + dpml=1.0, + margin_xy=0.5, + margin_z_above=0.5, + margin_z_below=0.5, + port_margin=0.5, + extend_ports=0.0, + source_port_offset=0.1, + distance_source_to_monitors=0.2, + ), + accuracy=AccuracyConfig( + eps_averaging=False, + subpixel_maxeval=0, + subpixel_tol=1e-4, + simplify_tol=0.0, + ), + diagnostics=DiagnosticsConfig( + save_geometry=True, + save_fields=True, + save_epsilon_raw=False, + save_animation=False, + animation_interval=0.5, + preview_only=False, + verbose_interval=0, + ), + symmetries=[], + ) + + def test_default_none(self, _sim_kwargs): + cfg = SimConfig(**_sim_kwargs) assert cfg.component_bbox is None - def test_with_bbox(self): - cfg = SimConfig(component_bbox=[-5.0, -2.0, 5.0, 2.0]) + def test_with_bbox(self, _sim_kwargs): + cfg = SimConfig(**_sim_kwargs, component_bbox=[-5.0, -2.0, 5.0, 2.0]) assert cfg.component_bbox == [-5.0, -2.0, 5.0, 2.0] - def test_json_roundtrip(self, tmp_path): - cfg = SimConfig(component_bbox=[-1.0, -0.5, 1.0, 0.5]) + def test_json_roundtrip(self, tmp_path, _sim_kwargs): + cfg = SimConfig(**_sim_kwargs, component_bbox=[-1.0, -0.5, 1.0, 0.5]) path = tmp_path / "config.json" cfg.to_json(path) data = json.loads(path.read_text()) assert data["component_bbox"] == [-1.0, -0.5, 1.0, 0.5] - def test_json_roundtrip_none(self, tmp_path): - cfg = SimConfig() + def test_json_roundtrip_none(self, tmp_path, _sim_kwargs): + cfg = SimConfig(**_sim_kwargs) path = tmp_path / "config.json" cfg.to_json(path) data = json.loads(path.read_text()) @@ -587,14 +702,52 @@ class TestDomainInSimConfig: """Test domain serialization in SimConfig JSON.""" def test_domain_in_sim_config_json(self, tmp_path): + from gsim.meep.models.config import AccuracyConfig, DiagnosticsConfig + cfg = SimConfig( gds_filename="layout.gds", + verbose_interval=0, layer_stack=[], + dielectrics=[], ports=[], materials={}, - wavelength=WavelengthConfig(), - resolution=ResolutionConfig(), - domain=DomainConfig(dpml=0.5, margin_xy=0.2), + wavelength=WavelengthConfig(wavelength=1.55, bandwidth=0.1, num_freqs=11), + source=SourceConfig(), + stopping=StoppingConfig( + mode="fixed", + max_time=100.0, + decay_dt=50.0, + decay_component="Ey", + threshold=0.05, + dft_min_run_time=100, + ), + resolution=ResolutionConfig(pixels_per_um=32), + domain=DomainConfig( + dpml=0.5, + margin_xy=0.2, + margin_z_above=0.5, + margin_z_below=0.5, + port_margin=0.5, + extend_ports=0.0, + source_port_offset=0.1, + distance_source_to_monitors=0.2, + ), + accuracy=AccuracyConfig( + eps_averaging=False, + subpixel_maxeval=0, + subpixel_tol=1e-4, + simplify_tol=0.0, + ), + diagnostics=DiagnosticsConfig( + save_geometry=True, + save_fields=True, + save_epsilon_raw=False, + save_animation=False, + animation_interval=0.5, + preview_only=False, + verbose_interval=0, + ), + symmetries=[], ) path = tmp_path / "config.json" cfg.to_json(path) @@ -618,9 +771,9 @@ def test_creation(self): assert s.direction == "X" assert s.phase == -1 - def test_default_phase(self): - s = SymmetryEntry(direction="Y") - assert s.phase == 1 + def test_requires_phase(self): + with pytest.raises(ValidationError): + SymmetryEntry(direction="Y") def test_model_dump(self): s = SymmetryEntry(direction="Z", phase=-1) @@ -644,37 +797,121 @@ def test_invalid_phase(self): class TestStoppingConfig: """Test StoppingConfig model.""" - def test_default_is_fixed(self): - cfg = StoppingConfig() - assert cfg.mode == "fixed" - assert cfg.max_time == 100.0 - assert cfg.decay_dt == 50.0 - assert cfg.decay_component == "Ey" - assert cfg.threshold == 1e-3 - assert cfg.decay_monitor_port is None - - def test_decay_mode(self): - cfg = StoppingConfig(mode="decay", threshold=1e-4, decay_dt=25.0) - assert cfg.mode == "decay" + def test_requires_all_fields(self): + with pytest.raises(ValidationError): + StoppingConfig() + + def test_field_decay_mode(self): + cfg = StoppingConfig( + mode="field_decay", + threshold=1e-4, + decay_dt=25.0, + max_time=100.0, + decay_component="Ey", + dft_min_run_time=100, + ) + assert cfg.mode == "field_decay" assert cfg.threshold == 1e-4 assert cfg.decay_dt == 25.0 + def test_energy_decay_mode(self): + cfg = StoppingConfig( + mode="energy_decay", + decay_dt=100, + threshold=1e-4, + max_time=100.0, + decay_component="Ey", + dft_min_run_time=100, + ) + assert cfg.mode == "energy_decay" + assert cfg.decay_dt == 100 + assert cfg.threshold == 1e-4 + def test_invalid_mode(self): with pytest.raises(ValidationError): - StoppingConfig(mode="invalid") # ty: ignore[invalid-argument-type] + StoppingConfig( # ty: ignore[invalid-argument-type] + mode="invalid", + max_time=100.0, + decay_dt=50.0, + decay_component="Ey", + threshold=0.05, + dft_min_run_time=100, + ) def test_threshold_bounds(self): with pytest.raises(ValidationError): - StoppingConfig(threshold=0) + StoppingConfig( + mode="fixed", + max_time=100.0, + decay_dt=50.0, + decay_component="Ey", + threshold=0, + dft_min_run_time=100, + ) with pytest.raises(ValidationError): - StoppingConfig(threshold=1.0) + StoppingConfig( + mode="fixed", + max_time=100.0, + decay_dt=50.0, + decay_component="Ey", + threshold=1.0, + dft_min_run_time=100, + ) + + def test_wall_time_max_default(self): + cfg = StoppingConfig( + mode="fixed", + max_time=100.0, + decay_dt=50.0, + decay_component="Ey", + threshold=0.05, + dft_min_run_time=100, + ) + assert cfg.wall_time_max == 0.0 + + def test_wall_time_max_set(self): + cfg = StoppingConfig( + mode="field_decay", + max_time=200.0, + decay_dt=50.0, + decay_component="Ey", + threshold=0.05, + dft_min_run_time=100, + wall_time_max=3600, + ) + assert cfg.wall_time_max == 3600 + + def test_wall_time_max_serialized(self): + cfg = StoppingConfig( + mode="fixed", + max_time=100.0, + decay_dt=50.0, + decay_component="Ey", + threshold=0.05, + dft_min_run_time=100, + wall_time_max=1800, + ) + data = cfg.model_dump() + assert data["wall_time_max"] == 1800 + + def test_wall_time_max_negative_rejected(self): + with pytest.raises(ValidationError): + StoppingConfig( + mode="fixed", + max_time=100.0, + decay_dt=50.0, + decay_component="Ey", + threshold=0.05, + dft_min_run_time=100, + wall_time_max=-1, + ) class TestWavelengthConfigStopping: """Test that WavelengthConfig no longer embeds StoppingConfig.""" def test_model_dump_excludes_stopping(self): - cfg = WavelengthConfig() + cfg = WavelengthConfig(wavelength=1.55, bandwidth=0.1, num_freqs=11) d = cfg.model_dump() assert "stopping" not in d assert "run_after_sources" not in d @@ -729,7 +966,7 @@ def test_auto_fwidth_wider_than_monitor(self): """Auto source fwidth should always be wider than monitor df.""" cfg = SourceConfig() fcen = 1.0 / 1.55 - wl = WavelengthConfig(wavelength=1.55, bandwidth=0.1) + wl = WavelengthConfig(wavelength=1.55, bandwidth=0.1, num_freqs=11) fwidth = cfg.compute_fwidth(fcen, wl.df) assert fwidth > wl.df @@ -783,6 +1020,13 @@ def test_script_has_component_map(self): assert "mp.Ez" in script assert "mp.Ey" in script + def test_script_has_stop_when_energy_decayed(self): + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + assert "stop_when_energy_decayed" in script + assert "energy_decay" in script + def test_script_valid_python(self): from gsim.meep.script import generate_meep_script @@ -860,7 +1104,14 @@ def test_build_sim_overlay(self): ) domain_cfg = DomainConfig( - dpml=1.0, margin_xy=0.5, margin_z_above=0.0, margin_z_below=0.0 + dpml=1.0, + margin_xy=0.5, + margin_z_above=0.0, + margin_z_below=0.0, + port_margin=0.5, + extend_ports=0.0, + source_port_offset=0.1, + distance_source_to_monitors=0.2, ) port_data = [ @@ -919,7 +1170,16 @@ def test_build_sim_overlay_with_dielectrics(self): bbox=((-2.0, -1.0, -1.0), (2.0, 1.0, 1.22)), ) - domain_cfg = DomainConfig(dpml=1.0, margin_xy=0.5) + domain_cfg = DomainConfig( + dpml=1.0, + margin_xy=0.5, + margin_z_above=0.5, + margin_z_below=0.5, + port_margin=0.5, + extend_ports=0.0, + source_port_offset=0.1, + distance_source_to_monitors=0.2, + ) dielectrics = [ {"name": "substrate", "material": "silicon", "zmin": -1.0, "zmax": 0.0}, {"name": "oxide", "material": "SiO2", "zmin": 0.0, "zmax": 1.22}, @@ -977,6 +1237,22 @@ def test_script_handles_component_bbox(self): assert "bbox_top" in script assert "bbox_bottom" in script + def test_script_contains_wall_time_cap(self): + """Verify runner script has wall-clock time cap logic.""" + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + assert "_make_wall_time_cap" in script + assert "wall_time_max" in script + + def test_script_wall_time_all_modes(self): + """Wall-clock cap should be wired into all 4 stopping modes.""" + from gsim.meep.script import generate_meep_script + + script = generate_meep_script() + # Each mode branch should reference _make_wall_time_cap + assert script.count("_make_wall_time_cap") >= 4 + # --------------------------------------------------------------------------- # Render2d overlay tests diff --git a/tests/meep/test_simulation.py b/tests/meep/test_simulation.py index b84009cc..58741a95 100644 --- a/tests/meep/test_simulation.py +++ b/tests/meep/test_simulation.py @@ -44,14 +44,14 @@ def test_defaults(self): s = ModeSource() assert s.port is None assert s.wavelength == 1.55 - assert s.bandwidth == 0.1 + assert s.wavelength_span == 0.1 assert s.num_freqs == 11 def test_custom(self): - s = ModeSource(port="o1", wavelength=1.31, bandwidth=0.05, num_freqs=21) + s = ModeSource(port="o1", wavelength=1.31, wavelength_span=0.05, num_freqs=21) assert s.port == "o1" assert s.wavelength == 1.31 - assert s.bandwidth == 0.05 + assert s.wavelength_span == 0.05 assert s.num_freqs == 21 @@ -64,6 +64,8 @@ def test_defaults(self): assert d.margin_z_below == 0.5 assert d.port_margin == 0.5 assert d.extend_ports == 0.0 + assert d.source_port_offset == 0.1 + assert d.distance_source_to_monitors == 0.2 assert d.symmetries == [] def test_custom(self): @@ -82,9 +84,9 @@ def test_symmetries(self): class TestStoppingFields: def test_defaults(self): f = FDTD() - assert f.stopping == "dft_decay" - assert f.max_time == 200.0 - assert f.stopping_threshold == 1e-3 + assert f.stopping == "field_decay" + assert f.max_time == 2000.0 + assert f.stopping_threshold == 0.05 assert f.stopping_min_time == 100.0 assert f.stopping_component == "Ey" assert f.stopping_dt == 50.0 @@ -95,19 +97,30 @@ def test_fixed_mode(self): assert f.stopping == "fixed" assert f.max_time == 100 - def test_decay_mode(self): - f = FDTD(stopping="decay", stopping_threshold=1e-4, stopping_dt=25) - assert f.stopping == "decay" + def test_field_decay_mode(self): + f = FDTD(stopping="field_decay", stopping_threshold=1e-4, stopping_dt=25) + assert f.stopping == "field_decay" assert f.stopping_threshold == 1e-4 assert f.stopping_dt == 25 def test_dft_decay_custom(self): - f = FDTD(max_time=300, stopping_threshold=1e-4, stopping_min_time=80) + f = FDTD( + stopping="dft_decay", + max_time=300, + stopping_threshold=1e-4, + stopping_min_time=80, + ) assert f.stopping == "dft_decay" assert f.max_time == 300 assert f.stopping_threshold == 1e-4 assert f.stopping_min_time == 80 + def test_energy_decay_mode(self): + f = FDTD(stopping="energy_decay", stopping_dt=100, stopping_threshold=1e-4) + assert f.stopping == "energy_decay" + assert f.stopping_dt == 100 + assert f.stopping_threshold == 1e-4 + def test_invalid_mode(self): with pytest.raises(ValidationError): FDTD(stopping="invalid") # ty: ignore[invalid-argument-type] @@ -117,14 +130,14 @@ class TestFDTD: def test_defaults(self): f = FDTD() assert f.resolution == 32 - assert f.stopping == "dft_decay" - assert f.max_time == 200.0 + assert f.stopping == "field_decay" + assert f.max_time == 2000.0 assert f.subpixel is False assert f.simplify_tol == 0.0 def test_with_custom_stopping(self): - f = FDTD(stopping="decay", stopping_threshold=1e-4) - assert f.stopping == "decay" + f = FDTD(stopping="field_decay", stopping_threshold=1e-4) + assert f.stopping == "field_decay" assert f.stopping_threshold == 1e-4 def test_resolution_custom(self): @@ -206,17 +219,108 @@ def test_solver_resolution(self): def test_solver_stopping_replace(self): sim = Simulation() - sim.solver.stopping = "decay" + sim.solver.stopping = "field_decay" sim.solver.stopping_threshold = 1e-4 - assert sim.solver.stopping == "decay" + assert sim.solver.stopping == "field_decay" assert sim.solver.stopping_threshold == 1e-4 + def test_stop_when_energy_decayed(self): + f = FDTD() + result = f.stop_when_energy_decayed(dt=100, decay_by=1e-4) + assert result is f + assert f.stopping == "energy_decay" + assert f.stopping_dt == 100 + assert f.stopping_threshold == 1e-4 + + def test_stop_when_dft_decayed(self): + f = FDTD() + result = f.stop_when_dft_decayed(tol=1e-4, min_time=200) + assert result is f + assert f.stopping == "dft_decay" + assert f.stopping_threshold == 1e-4 + assert f.stopping_min_time == 200 + + def test_stop_when_fields_decayed(self): + f = FDTD() + result = f.stop_when_fields_decayed( + dt=25, component="Hz", decay_by=1e-4, monitor_port="o2" + ) + assert result is f + assert f.stopping == "field_decay" + assert f.stopping_dt == 25 + assert f.stopping_component == "Hz" + assert f.stopping_threshold == 1e-4 + assert f.stopping_monitor_port == "o2" + + def test_stop_after_sources(self): + f = FDTD() + result = f.stop_after_sources(time=500) + assert result is f + assert f.stopping == "fixed" + assert f.max_time == 500 + + def test_wall_time_max_default(self): + f = FDTD() + assert f.wall_time_max == 0.0 + + def test_stop_after_walltime(self): + f = FDTD() + result = f.stop_after_walltime(seconds=3600) + assert result is f + assert f.wall_time_max == 3600 + + def test_wall_time_max_orthogonal(self): + """wall_time_max can be combined with any stopping mode.""" + f = FDTD() + f.stop_when_fields_decayed(dt=50, decay_by=0.05) + f.stop_after_walltime(seconds=1800) + assert f.stopping == "field_decay" + assert f.wall_time_max == 1800 + def test_geometry_component(self): sim = Simulation() sim.geometry.component = "placeholder" assert sim.geometry.component == "placeholder" +class TestCallableAPI: + def test_source_callable(self): + sim = Simulation() + result = sim.source(port="o1", wavelength=1.31, wavelength_span=0.05) + assert result is sim.source + assert sim.source.port == "o1" + assert sim.source.wavelength == 1.31 + assert sim.source.wavelength_span == 0.05 + + def test_domain_callable(self): + sim = Simulation() + sim.domain(pml=0.5, margin=0.2) + assert sim.domain.pml == 0.5 + assert sim.domain.margin == 0.2 + + def test_geometry_callable(self): + sim = Simulation() + sim.geometry(component="placeholder", z_crop="auto") + assert sim.geometry.component == "placeholder" + assert sim.geometry.z_crop == "auto" + + def test_solver_callable(self): + sim = Simulation() + sim.solver(resolution=64, simplify_tol=0.01) + assert sim.solver.resolution == 64 + assert sim.solver.simplify_tol == 0.01 + + def test_callable_validation(self): + sim = Simulation() + with pytest.raises(ValidationError): + sim.solver(resolution=1) # ge=4 + + def test_callable_invalid_field(self): + sim = Simulation() + with pytest.raises(ValidationError): + sim.source(nonexistent_field="value") + + class TestWholeObjectAssignment: def test_source(self): sim = Simulation() @@ -280,7 +384,7 @@ def test_from_source_defaults(self): def test_from_source_custom(self): sim = Simulation() - sim.source = ModeSource(wavelength=1.31, bandwidth=0.05, num_freqs=21) + sim.source = ModeSource(wavelength=1.31, wavelength_span=0.05, num_freqs=21) wl = sim._wavelength_config() assert wl.wavelength == 1.31 assert wl.bandwidth == 0.05 @@ -301,14 +405,14 @@ def test_stopping_fixed(self): assert cfg.mode == "fixed" assert cfg.max_time == 200 - def test_stopping_decay(self): + def test_stopping_field_decay(self): sim = Simulation() - sim.solver.stopping = "decay" + sim.solver.stopping = "field_decay" sim.solver.stopping_threshold = 1e-4 sim.solver.stopping_dt = 25 sim.solver.stopping_monitor_port = "o2" cfg = sim._stopping_config() - assert cfg.mode == "decay" + assert cfg.mode == "field_decay" assert cfg.threshold == 1e-4 assert cfg.decay_dt == 25 assert cfg.decay_monitor_port == "o2" @@ -325,6 +429,14 @@ def test_stopping_dft_decay(self): assert cfg.threshold == 1e-4 assert cfg.dft_min_run_time == 80 + def test_stopping_energy_decay(self): + sim = Simulation() + sim.solver.stop_when_energy_decayed(dt=100, decay_by=1e-4) + cfg = sim._stopping_config() + assert cfg.mode == "energy_decay" + assert cfg.decay_dt == 100 + assert cfg.threshold == 1e-4 + def test_domain_translation(self): sim = Simulation() sim.domain = Domain(pml=0.5, margin=0.3, margin_z_above=1.0, margin_z_below=0.8) @@ -333,6 +445,8 @@ def test_domain_translation(self): assert cfg.margin_xy == 0.3 assert cfg.margin_z_above == 1.0 assert cfg.margin_z_below == 0.8 + assert cfg.source_port_offset == 0.1 + assert cfg.distance_source_to_monitors == 0.2 def test_resolution_translation(self): sim = Simulation() @@ -350,7 +464,7 @@ def test_accuracy_translation(self): def test_source_translation(self): sim = Simulation() - sim.source = ModeSource(port="o1", bandwidth=0.05) + sim.source = ModeSource(port="o1") cfg = sim._source_config() assert cfg.port == "o1" assert cfg.bandwidth is None # always auto fwidth