diff --git a/.github/workflows/check_docs.yml b/.github/workflows/check_docs.yml index 93673445f..14d8b2618 100644 --- a/.github/workflows/check_docs.yml +++ b/.github/workflows/check_docs.yml @@ -1,9 +1,7 @@ name: Check API Documentation on: - pull_request: - branches: - - devel + workflow_dispatch: concurrency: group: ${{ github.workflow }}-${{ github.ref }} diff --git a/.github/workflows/docs-preview-cleanup.yml b/.github/workflows/docs-preview-cleanup.yml index 742cfc10a..d210afebf 100644 --- a/.github/workflows/docs-preview-cleanup.yml +++ b/.github/workflows/docs-preview-cleanup.yml @@ -1,10 +1,7 @@ name: Remove PR docs preview on: - pull_request: - branches: - - devel - types: [closed] + workflow_dispatch: defaults: run: diff --git a/.github/workflows/docs-preview.yml b/.github/workflows/docs-preview.yml index fd17f1678..3ca5737fc 100644 --- a/.github/workflows/docs-preview.yml +++ b/.github/workflows/docs-preview.yml @@ -1,10 +1,6 @@ name: Deploy PR docs preview on: - pull_request: - branches: - - devel - types: [opened, synchronize, reopened] workflow_dispatch: inputs: pr_number: diff --git a/.github/workflows/pr-feectools-submodule.yml b/.github/workflows/pr-feectools-submodule.yml index 2fa728b8e..7982b12c2 100644 --- a/.github/workflows/pr-feectools-submodule.yml +++ b/.github/workflows/pr-feectools-submodule.yml @@ -1,9 +1,6 @@ name: PR - feectools submodule freshness on: - pull_request: - branches: - - devel workflow_dispatch: concurrency: diff --git a/.github/workflows/profiling-clusters.yml b/.github/workflows/profiling-clusters.yml index 18b3c6b31..0c4e718ad 100644 --- a/.github/workflows/profiling-clusters.yml +++ b/.github/workflows/profiling-clusters.yml @@ -1,12 +1,5 @@ name: Cluster workflow on: - pull_request: - branches: - - devel-unused-for-now - - # schedule: - # # run at 1 a.m. every night - # - cron: "0 1 * * *" workflow_dispatch: jobs: diff --git a/.github/workflows/pypi-release.yml b/.github/workflows/pypi-release.yml index 2afbaa837..3ad572bd9 100644 --- a/.github/workflows/pypi-release.yml +++ b/.github/workflows/pypi-release.yml @@ -26,7 +26,8 @@ jobs: - name: Set up Python uses: actions/setup-python@v4 with: - python-version: "3.10" + # browser/build.py uses tomllib and emits the pure-Python Pyodide wheel. + python-version: "3.11" - name: Install build tools run: | @@ -36,7 +37,13 @@ jobs: - name: Build the package run: python -m build + - name: Build the Pyodide-compatible geometry package + run: python browser/build.py --output dist + + - name: Show release distributions + run: ls -l dist/ + - name: Publish package distributions to PyPI uses: pypa/gh-action-pypi-publish@release/v1 # with: - # repository-url: https://test.pypi.org/legacy/ \ No newline at end of file + # repository-url: https://test.pypi.org/legacy/ diff --git a/.github/workflows/static_analysis.yml b/.github/workflows/static_analysis.yml index e2f031bfc..fb6f71b00 100644 --- a/.github/workflows/static_analysis.yml +++ b/.github/workflows/static_analysis.yml @@ -1,9 +1,7 @@ name: isort & ruff on: - pull_request: - branches: - - devel + workflow_dispatch: concurrency: group: ${{ github.workflow }}-${{ github.ref }} diff --git a/.github/workflows/submod-feectools.yml b/.github/workflows/submod-feectools.yml index 33f8d9e40..f0d692dee 100644 --- a/.github/workflows/submod-feectools.yml +++ b/.github/workflows/submod-feectools.yml @@ -1,9 +1,6 @@ name: Submodule feectools on: - pull_request: - branches: - - devel workflow_dispatch: concurrency: @@ -32,4 +29,4 @@ jobs: - name: Run workflow if submodule changed (all PR commits) if: env.SUBMOD_CHANGED == 'true' run: | - echo "${{ env.SUBMOD_NAME }} has changed, running tests..." \ No newline at end of file + echo "${{ env.SUBMOD_NAME }} has changed, running tests..." diff --git a/.github/workflows/test-PR-examples.yml b/.github/workflows/test-PR-examples.yml index 83cb3f8a5..73e4e2b37 100644 --- a/.github/workflows/test-PR-examples.yml +++ b/.github/workflows/test-PR-examples.yml @@ -1,12 +1,6 @@ name: PR - examples tests in Container on: - pull_request: - branches: - - devel - schedule: - # run at 1 a.m. on Sunday - - cron: "0 1 * * 0" workflow_dispatch: concurrency: diff --git a/.github/workflows/test-PR-generate-params.yml b/.github/workflows/test-PR-generate-params.yml index 7bc4b016b..8da3940eb 100644 --- a/.github/workflows/test-PR-generate-params.yml +++ b/.github/workflows/test-PR-generate-params.yml @@ -1,12 +1,6 @@ name: PR - build, install and generate params on: - pull_request: - branches: - - devel - schedule: - # run at 1 a.m. on Sunday - - cron: "0 1 * * 0" workflow_dispatch: concurrency: diff --git a/.github/workflows/test-PR-models-clones.yml b/.github/workflows/test-PR-models-clones.yml index 7b729f217..a78a05ddd 100644 --- a/.github/workflows/test-PR-models-clones.yml +++ b/.github/workflows/test-PR-models-clones.yml @@ -1,12 +1,6 @@ name: PR - domain cloning in Container on: - pull_request: - branches: - - devel - schedule: - # run at 1 a.m. on Sunday - - cron: "0 1 * * 0" workflow_dispatch: concurrency: @@ -88,4 +82,4 @@ jobs: source /struphy_fortran_/env_fortran_/bin/activate cd /struphy_fortran_/src/struphy mpirun -n 1 pytest -m single --testmon-forceselect -xs --with-mpi --model-name Maxwell - mpirun --oversubscribe -n 4 pytest -x --testmon --with-mpi --nclones 2 models/tests/verification/ \ No newline at end of file + mpirun --oversubscribe -n 4 pytest -x --testmon --with-mpi --nclones 2 models/tests/verification/ diff --git a/.github/workflows/test-PR-models.yml b/.github/workflows/test-PR-models.yml index 4a2e83442..f262e93af 100644 --- a/.github/workflows/test-PR-models.yml +++ b/.github/workflows/test-PR-models.yml @@ -1,12 +1,6 @@ name: PR - model tests in Container on: - pull_request: - branches: - - devel - schedule: - # run at 1 a.m. on Sunday - - cron: "0 1 * * 0" workflow_dispatch: concurrency: @@ -121,4 +115,4 @@ jobs: TESTMON_DATAFILE: ${{ github.workspace }}/.testmondata-model-mpi run: | source /struphy_fortran_/env_fortran_/bin/activate - mpirun --oversubscribe -n 2 pytest -x --testmon --with-mpi /struphy_fortran_/src/struphy/models/tests/verification/ \ No newline at end of file + mpirun --oversubscribe -n 2 pytest -x --testmon --with-mpi /struphy_fortran_/src/struphy/models/tests/verification/ diff --git a/.github/workflows/test-PR-pure-python.yml b/.github/workflows/test-PR-pure-python.yml index 85f52858f..52147c398 100644 --- a/.github/workflows/test-PR-pure-python.yml +++ b/.github/workflows/test-PR-pure-python.yml @@ -1,12 +1,6 @@ name: PR - pure Python tests in Container on: - pull_request: - branches: - - devel - schedule: - # run at 1 a.m. on Sunday - - cron: "0 1 * * 0" workflow_dispatch: concurrency: diff --git a/.github/workflows/test-PR-tutorials.yml b/.github/workflows/test-PR-tutorials.yml index 9acb3754f..fe3d13a79 100644 --- a/.github/workflows/test-PR-tutorials.yml +++ b/.github/workflows/test-PR-tutorials.yml @@ -1,12 +1,6 @@ name: PR - tutorial tests in Container on: - pull_request: - branches: - - devel - schedule: - # run at 1 a.m. on Sunday - - cron: "0 1 * * 0" workflow_dispatch: concurrency: @@ -68,4 +62,4 @@ jobs: cd /struphy_fortran_/tutorials ls jupyter nbconvert --show-config - jupyter nbconvert --execute --embed-images --to html tutorial_*.ipynb \ No newline at end of file + jupyter nbconvert --execute --embed-images --to html tutorial_*.ipynb diff --git a/.github/workflows/test-PR-unit-mpi.yml b/.github/workflows/test-PR-unit-mpi.yml index e3a13dbc5..fbe7cd421 100644 --- a/.github/workflows/test-PR-unit-mpi.yml +++ b/.github/workflows/test-PR-unit-mpi.yml @@ -1,12 +1,6 @@ name: PR - unit tests with MPI in Container on: - pull_request: - branches: - - devel - schedule: - # run at 1 a.m. on Sunday - - cron: "0 1 * * 0" workflow_dispatch: concurrency: diff --git a/.github/workflows/test-PR-unit.yml b/.github/workflows/test-PR-unit.yml index 10cdf62c4..21814d424 100644 --- a/.github/workflows/test-PR-unit.yml +++ b/.github/workflows/test-PR-unit.yml @@ -1,12 +1,6 @@ name: PR - unit tests serial in Container on: - pull_request: - branches: - - devel - schedule: - # run at 1 a.m. on Sunday - - cron: "0 1 * * 0" workflow_dispatch: concurrency: @@ -20,4 +14,4 @@ jobs: os: ubuntu-latest n-procs: 1 secrets: - ghcr-token: ${{ secrets.GHCR_TOKEN }} \ No newline at end of file + ghcr-token: ${{ secrets.GHCR_TOKEN }} diff --git a/.github/workflows/test-browser-geometry.yml b/.github/workflows/test-browser-geometry.yml new file mode 100644 index 000000000..48aed41ff --- /dev/null +++ b/.github/workflows/test-browser-geometry.yml @@ -0,0 +1,45 @@ +name: Browser geometry +on: + pull_request: + paths: + - 'browser/**' + - 'src/struphy/geometry/**' + - 'src/struphy/bsplines/**' + - 'src/struphy/linear_algebra/**' + - 'src/struphy/kernel_arguments/**' + - 'src/struphy/utils/**' + - 'pyproject.toml' + - '.github/workflows/test-browser-geometry.yml' + workflow_dispatch: +permissions: + contents: read +jobs: + geometry: + runs-on: ubuntu-latest + timeout-minutes: 20 + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: '3.13' + - uses: actions/setup-node@v4 + with: + node-version: '22' + cache: npm + cache-dependency-path: browser/package-lock.json + - run: python browser/build.py + - run: python -m pip install browser/dist/*.whl pytest + - run: python -m pytest -q browser/tests + env: + PYTEST_DISABLE_PLUGIN_AUTOLOAD: '1' + - run: npm ci --prefix browser + - run: python browser/prepare.py + - run: npm exec --prefix browser -- playwright install --with-deps chromium + - run: npm test --prefix browser + - uses: actions/upload-artifact@v4 + if: always() + with: + name: browser-geometry + path: | + browser/dist/*.whl + browser/test-results/ diff --git a/browser/.gitignore b/browser/.gitignore new file mode 100644 index 000000000..b310864ce --- /dev/null +++ b/browser/.gitignore @@ -0,0 +1,11 @@ +/dist/ +/.release-dist-check/ +/site/vendor/ +/site/*.whl +/site/wheel.json +/node_modules/ +/.test-env/ + +/site/test_geometry.py +/test-results/ +/playwright-report/ diff --git a/browser/README.md b/browser/README.md new file mode 100644 index 000000000..538920d2d --- /dev/null +++ b/browser/README.md @@ -0,0 +1,149 @@ +# Struphy geometry in the browser + +`struphy-geometry` is a pure-Python wheel built from this repository's shared +geometry sources. Import it as `struphy_geometry`. It can coexist with full +Struphy and retains cuNumPy on its NumPy backend. No Pyccel compilation, MPI, +FEEC installation, GPU, or computation server is required. + +## Run the viewer + +From the repository root (Python >= 3.11 and Node.js are build tools only): + +```sh +npm ci --prefix browser +python3 browser/prepare.py +python3 -m http.server 8765 --directory browser/site +``` + +Open **http://localhost:8765**. Choose a domain, edit its parameters, and click +**Construct domain**. Drag to rotate the sampled surface. The probe displays +physical coordinates, the Jacobian, determinant, and metric. Choose `Tokamak` +to construct flux-aligned geometry from circular flux or an uploaded G-EQDSK +file. **Save geometry** and **Open geometry** exchange `.npz` files. + +The preparation step downloads the pinned Pyodide runtime and dependencies, +verifies downloaded package hashes, and copies them into `site/vendor`. +After preparation, the site uses only same-origin static assets. All geometry +calculations and file parsing run in a browser Web Worker. Files are not +uploaded to a computation server. Any static HTTP(S) host can serve `site/`; +opening `index.html` directly through `file://` is not supported. Vendored +assets are generated and ignored by Git. A fresh load still needs access to +the static host; this is not a service-worker offline-installable application. + +## Use the Python API + +Build just the wheel with `python3 browser/build.py`. Install the resulting +`browser/dist/struphy_geometry--py3-none-any.whl` in Python, or install +its hosted URL with `await micropip.install(...)` in Pyodide. The runtime +requires NumPy, SciPy, and cuNumPy; the viewer pins their tested distribution. + +```python +import cunumpy as xp +from struphy_geometry import domains, CircularFlux, EQDSKFlux +from struphy_geometry import save_geometry, load_geometry, sample_surface + +xp.set_backend("numpy") +domain = domains.IGAPolarTorus(a=0.7, R0=3.0, tor_period=1) +eta = xp.linspace(0.0, 1.0, 24) +xyz = domain(0.8, eta, eta) # a logical coordinate surface +J = domain.jacobian(0.4, 0.2, 0.3, squeeze_out=True) +g = domain.metric(0.4, 0.2, 0.3, squeeze_out=True) + +flux = CircularFlux(R0=3.0, a=1.0) +# For an uploaded file: flux = EQDSKFlux.from_text(uploaded_text) +tokamak = domains.Tokamak( + equilibrium=flux, + num_elements=(8, 24), degree=(2, 3), + xi_param="equal_angle", tor_period=1, +) +archive = save_geometry(tokamak) # bytes, no pickle +restored = load_geometry(archive) # needs no equilibrium object +x, y, z = sample_surface(restored, resolution=24) +``` + +Keep existing evaluation conventions: three coordinate arrays produce tensor +products; a single `(N, 3)` array evaluates a list of points. Use +`squeeze_out=True` for scalar probes. `sample_surface` fixes eta1; it does not +assemble all boundary faces of a general domain. Cartesian coordinates are +in mapping/input-file units; the browser package does not rescale EQDSK data. + +## Included functionality + +- All nine analytic domains in `geometry/domains.py`. +- `IGAPolarCylinder`, `IGAPolarTorus`, generic `Spline`, + `PoloidalSplineStraight`, and `PoloidalSplineTorus`. +- Knot generation, B-spline evaluation, 2D/3D mapping interpolation, and + `spline_interpolation_nd`. +- Jacobians, determinants, inverses, metrics, pullbacks, pushforwards, and + field transformations. +- Tokamak flux tracing with `equal_angle`, `equal_arc_length`, `sfl`, + `equal_area`, and `equal_volume` parametrizations. Custom flux providers + must expose `psi(R, Z, dR=0, dZ=0)`, `psi_axis_RZ`, and `psi_range`. +- EQDSK rectangular poloidal flux interpolation, including first/second + derivatives. The reader accepts scientific E/D notation. It refines the + interpolated magnetic axis for either flux sign and performs no smoothing. + It does not implement the pressure/current profiles of `EQDSKequilibrium`. +- Portable archives of analytic definitions and spline control points. +- Existing Matplotlib `show()` can be used if Matplotlib is installed and + an appropriate display backend is configured; the viewer uses canvas. + +## Deliberate differences from full Struphy + +- The import name is `struphy_geometry`, and the wheel does not depend on or + initialize full Struphy. The full `struphy` distribution remains unchanged + in scope and dependencies. +- Browser `Tokamak()` requires an explicit flux provider. The viewer supplies + `CircularFlux()` when no file is selected; it is a geometry example, not a + force-balanced equilibrium. Full Struphy retains its EQDSK default. +- Browser `Spline()` requires explicit `cx`, `cy`, and `cz` control points. +- `GVECunit` and `DESCunit` are omitted. Create these domains with full + Struphy, then use `struphy.geometry.portable.save_geometry(domain)` to export + them. Browser loading reconstructs the generic spline and preserves mapping + evaluation and derivatives; it does not run those equilibrium solvers. +- PyVista/VTK methods (`create_geometry_mesh`, `show_3d`, `export_geometry`) + are omitted. Mesh arrays and NPZ archives are available instead. +- Analytic `to_dict/from_dict` remains available. Use the archive functions + for spline/equilibrium-derived mappings; these save numerical coefficients + rather than trying to serialize/reconstruct arbitrary equilibrium objects. +- GPU execution and runtime compilation are not browser features. Kernels + run as Python functions with NumPy; cuNumPy's `PyccelKernel` wrapper remains + valid because it accepts ordinary callables. + +Dense sampling and fine flux tracing can be expensive in interpreted Python. +The worker keeps the UI responsive, but does not make the kernels faster. +Start with small grids. Polar axes are coordinate singularities; inspect +inverse metrics away from the axis. Flux tracing retains Struphy's assumptions +about nested surfaces and convergence, so not every EQDSK topology is suitable. + +## How the shared-source build works + +`build.py` has an explicit module allowlist. It generates an isolated namespace, +removes the `pure`/`stack_array` compilation decorators, excludes the native +adapters listed above, and writes a deterministic wheel with metadata, hashes, +and the license. It never copies `.so` files or modifies generated native +kernels. There is no separately maintained copy of the numerical algorithms. +Native compiler imports remain intact in the original sources. + +Small shared-source changes decouple geometry imports from plotting, option +and class-helper infrastructure, allow lightweight flux providers, and fix +the missing exponent factor in `powered_ellipse_df`. Existing compiled +Struphy installations must rebuild their kernels to pick up that fix. + +## Verification + +```sh +python3 browser/build.py +# In an environment with NumPy, SciPy, cuNumPy and pytest: +PYTEST_DISABLE_PLUGIN_AUTOLOAD=1 python3 -m pytest -q browser/tests + +# After preparing the site: +npm exec --prefix browser -- playwright install chromium +npm test --prefix browser +``` + +The same numerical tests run on native Python and in browser Pyodide, covering +all analytic/IGA domains, 3D spline construction, finite-difference Jacobians, +metric identities, transformations, all five flux parametrizations, EQDSK of +both signs, periodic seams, poles, out-of-domain points, and archive roundtrips. +The browser tests block external requests to check that static hosting is +sufficient, and exercise editing, error recovery, and saving/reopening geometry. diff --git a/browser/build.py b/browser/build.py new file mode 100644 index 000000000..695a93bcb --- /dev/null +++ b/browser/build.py @@ -0,0 +1,168 @@ +"""Build a pure-Python geometry wheel from the shared Struphy sources. + +Run with Python >= 3.11: python browser/build.py +No compiler, installed Struphy, or third-party build dependency is needed. +""" + +import argparse +import ast +import base64 +import csv +import hashlib +import io +import zipfile +from pathlib import Path + +import tomllib + +ROOT = Path(__file__).resolve().parents[1] +SOURCES = ( + "geometry/base.py", + "geometry/domains.py", + "geometry/utilities.py", + "geometry/mappings_kernels.py", + "geometry/evaluation_kernels.py", + "geometry/transform_kernels.py", + "geometry/utilities_kernels.py", + "geometry/flux.py", + "geometry/portable.py", + "bsplines/bsplines.py", + "bsplines/bsplines_kernels.py", + "bsplines/evaluation_kernels_2d.py", + "bsplines/evaluation_kernels_3d.py", + "kernel_arguments/pusher_args_kernels.py", + "linear_algebra/linalg_kernels.py", + "linear_algebra/linalg_kron.py", + "utils/class_helpers.py", + "utils/docstring_converter.py", +) + + +class BrowserSource(ast.NodeTransformer): + """Remove explicitly unsupported adapters and compiler-only annotations. + + Numerical function bodies are shared verbatim in meaning with Struphy. + cuNumPy, including its callable wrapper, remains a runtime dependency. + """ + + def visit_If(self, node): + if isinstance(node.test, ast.Name) and node.test.id == "TYPE_CHECKING": + return None + if isinstance(node.test, ast.Compare) and isinstance(node.test.left, ast.Name): + if node.test.left.id == "__name__": + return None + return self.generic_visit(node) + + def visit_ClassDef(self, node): + if node.name in ("GVECunit", "DESCunit"): + return None + if node.name == "Domain": + omitted = {"create_geometry_mesh", "show_3d", "export_geometry"} + node.body = [n for n in node.body if not isinstance(n, ast.FunctionDef) or n.name not in omitted] + return self.generic_visit(node) + + def visit_FunctionDef(self, node): + node.decorator_list = [ + d + for d in node.decorator_list + if not (isinstance(d, ast.Name) and d.id == "pure") + and not (isinstance(d, ast.Call) and isinstance(d.func, ast.Name) and d.func.id == "stack_array") + ] + return self.generic_visit(node) + + def visit_ImportFrom(self, node): + if node.module == "pyccel.decorators": + if any(n.name not in ("pure", "stack_array") for n in node.names): + raise ValueError("Review new Pyccel decorators before adding them to the browser build") + return None + if node.module and node.module.startswith("struphy.fields_background"): + names = {n.name for n in node.names} + if "EQDSKequilibrium" in names: + message = ( + "Tokamak requires equilibrium=CircularFlux() or EQDSKFlux.from_text(text) in the browser package" + ) + elif "GVECequilibrium" in names: + message = ( + "Spline requires explicit cx, cy, cz control points; use load_geometry for exported equilibria" + ) + else: + message = ( + "Equilibrium reconstruction is unavailable; use save_geometry/load_geometry for spline domains" + ) + return ast.Raise( + exc=ast.Call(func=ast.Name(id="ValueError", ctx=ast.Load()), args=[ast.Constant(message)], keywords=[]) + ) + return node + + +def browser_source(source): + tree = BrowserSource().visit(ast.parse(source)) + ast.fix_missing_locations(tree) + # Preserve leading copyright comments as well as the module docstring. + header = [] + for line in source.splitlines(): + if line.startswith("#") or not line.strip(): + header.append(line) + else: + break + text = "\n".join(header) + "\n" + ast.unparse(tree) + "\n" + return text.replace("struphy.", "struphy_geometry.") + + +def build(output): + version = tomllib.loads((ROOT / "pyproject.toml").read_text())["project"]["version"] + package = "struphy_geometry" + files = {} + for path in SOURCES: + files[f"{package}/{path}"] = browser_source((ROOT / "src/struphy" / path).read_text()).encode() + files[f"{package}/{Path(path).parent}/__init__.py"] = b"" + files[f"{package}/__init__.py"] = f'''"""Struphy geometry on NumPy, including Pyodide. No compilation required.""" +from .geometry import domains +from .geometry.base import Domain, Spline, PoloidalSplineStraight, PoloidalSplineTorus, interp_mapping, spline_interpolation_nd +from .geometry.flux import CircularFlux, EQDSKFlux +from .geometry.portable import sample_surface, save_geometry, load_geometry +__version__ = {version!r} +__all__ = ["domains", "Domain", "Spline", "PoloidalSplineStraight", "PoloidalSplineTorus", "interp_mapping", "spline_interpolation_nd", "CircularFlux", "EQDSKFlux", "sample_surface", "save_geometry", "load_geometry"] +'''.encode() + info = f"struphy_geometry-{version}.dist-info" + files[f"{info}/METADATA"] = f"""Metadata-Version: 2.1 +Name: struphy-geometry +Version: {version} +Summary: Portable Struphy domains, mappings, and flux geometry for Python and Pyodide +Requires-Python: >=3.10 +Requires-Dist: numpy +Requires-Dist: scipy +Requires-Dist: cunumpy>=0.1.4,<=0.1.5 +Provides-Extra: plot +Requires-Dist: matplotlib; extra == "plot" +License-File: LICENSE + +Built from the shared Struphy geometry sources. See browser/README.md in the Struphy repository. +""".encode() + files[f"{info}/LICENSE"] = (ROOT / "LICENSE").read_bytes() + files[f"{info}/WHEEL"] = ( + b"Wheel-Version: 1.0\nGenerator: struphy-geometry-build\nRoot-Is-Purelib: true\nTag: py3-none-any\n" + ) + record = io.StringIO(newline="") + writer = csv.writer(record) + for path, data in sorted(files.items()): + digest = base64.urlsafe_b64encode(hashlib.sha256(data).digest()).rstrip(b"=").decode() + writer.writerow([path, "sha256=" + digest, len(data)]) + writer.writerow([f"{info}/RECORD", "", ""]) + files[f"{info}/RECORD"] = record.getvalue().encode() + output.mkdir(parents=True, exist_ok=True) + target = output / f"struphy_geometry-{version}-py3-none-any.whl" + with zipfile.ZipFile(target, "w", compression=zipfile.ZIP_DEFLATED) as wheel: + for path, data in sorted(files.items()): + entry = zipfile.ZipInfo(path, date_time=(2020, 1, 1, 0, 0, 0)) + entry.compress_type = zipfile.ZIP_DEFLATED + entry.external_attr = 0o644 << 16 + wheel.writestr(entry, data) + print(target) + return target + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--output", type=Path, default=Path(__file__).parent / "dist") + build(parser.parse_args().output) diff --git a/browser/package-lock.json b/browser/package-lock.json new file mode 100644 index 000000000..bcef8e70b --- /dev/null +++ b/browser/package-lock.json @@ -0,0 +1,101 @@ +{ + "name": "struphy-geometry-browser", + "lockfileVersion": 3, + "requires": true, + "packages": { + "": { + "name": "struphy-geometry-browser", + "dependencies": { + "pyodide": "314.0.6" + }, + "devDependencies": { + "@playwright/test": "1.63.0" + } + }, + "node_modules/@playwright/test": { + "version": "1.63.0", + "resolved": "https://registry.npmjs.org/@playwright/test/-/test-1.63.0.tgz", + "integrity": "sha512-oxMK4vllB9RK5NQ2l1pq1IfOf2AvnEuj/vYGDj0H2nMtmtZpKtCwt/l00GEO6xjGfpBNAvjovvYdCm50dRQkpQ==", + "dev": true, + "license": "Apache-2.0", + "dependencies": { + "playwright": "1.63.0" + }, + "bin": { + "playwright": "cli.js" + }, + "engines": { + "node": ">=20" + } + }, + "node_modules/@types/emscripten": { + "version": "1.41.6", + "resolved": "https://registry.npmjs.org/@types/emscripten/-/emscripten-1.41.6.tgz", + "integrity": "sha512-uN+9i8bFT5CUcZfyIEYDrSueACEyKGbUs5kC/72DGlZZoinh84sJfVV0i8UOJD1asdzkvLPBRrKs41kZ8MdEXg==", + "license": "MIT" + }, + "node_modules/playwright": { + "version": "1.63.0", + "resolved": "https://registry.npmjs.org/playwright/-/playwright-1.63.0.tgz", + "integrity": "sha512-+7ziBLidS4NaNCdt57SUDT+wYmmd5fmiQejUic/kb+YsYSCPyOOE9sebzMjNmQrsnNpDJqd4WHvV/8lfKfUDUg==", + "dev": true, + "license": "Apache-2.0", + "dependencies": { + "playwright-core": "1.63.0" + }, + "bin": { + "playwright": "cli.js" + }, + "engines": { + "node": ">=20" + } + }, + "node_modules/playwright-core": { + "version": "1.63.0", + "resolved": "https://registry.npmjs.org/playwright-core/-/playwright-core-1.63.0.tgz", + "integrity": "sha512-rYCsBF/M5HjUch52bbtVONEFjv6Xu8sm8h72dNlR5bzIE1fvC/bxgspzkjSfU+MweEMmPM8KJebG6nnyxo5mCg==", + "dev": true, + "license": "Apache-2.0", + "bin": { + "playwright-core": "cli.js" + }, + "engines": { + "node": ">=20" + } + }, + "node_modules/pyodide": { + "version": "314.0.6", + "resolved": "https://registry.npmjs.org/pyodide/-/pyodide-314.0.6.tgz", + "integrity": "sha512-BKDTJyIqFxC4BExLqeRS3f5xvXZIjOt8C3zGLN/Cc7tFxSwvKVhVkchQQ2AGLOtR4YrVOIFxbV8poyDOOmWwxQ==", + "license": "MPL-2.0", + "dependencies": { + "@types/emscripten": "^1.41.4", + "ws": "^8.5.0" + }, + "engines": { + "node": ">=18.0.0" + } + }, + "node_modules/ws": { + "version": "8.21.3", + "resolved": "https://registry.npmjs.org/ws/-/ws-8.21.3.tgz", + "integrity": "sha512-201TZ/kPWxoPr/OKWjquZR1SWKXcvxdH+e1xrx89b3YbmzLMFCLfnaG1HFIgWzJOEWZ7MvpK++odZufgYR50Rw==", + "license": "MIT", + "engines": { + "node": ">=10.0.0" + }, + "peerDependencies": { + "bufferutil": "^4.0.1", + "utf-8-validate": ">=5.0.2" + }, + "peerDependenciesMeta": { + "bufferutil": { + "optional": true + }, + "utf-8-validate": { + "optional": true + } + } + } + } +} diff --git a/browser/package.json b/browser/package.json new file mode 100644 index 000000000..7c6201d5a --- /dev/null +++ b/browser/package.json @@ -0,0 +1,15 @@ +{ + "name": "struphy-geometry-browser", + "private": true, + "type": "module", + "scripts": { + "test": "playwright test", + "serve": "python3 -m http.server 8765 --directory site" + }, + "dependencies": { + "pyodide": "314.0.6" + }, + "devDependencies": { + "@playwright/test": "1.63.0" + } +} diff --git a/browser/playwright.config.js b/browser/playwright.config.js new file mode 100644 index 000000000..010d5bb0f --- /dev/null +++ b/browser/playwright.config.js @@ -0,0 +1,14 @@ +import { defineConfig } from "@playwright/test"; +export default defineConfig({ + testDir: "./tests", + testMatch: "*.spec.js", + timeout: 180000, + workers: 1, + use: { baseURL: "http://127.0.0.1:8765", headless: true }, + webServer: { + command: "python3 -m http.server 8765 --bind 127.0.0.1 --directory site", + url: "http://127.0.0.1:8765", + reuseExistingServer: !process.env.CI, + stderr: "ignore", + }, +}); diff --git a/browser/prepare.py b/browser/prepare.py new file mode 100644 index 000000000..b35b7dd21 --- /dev/null +++ b/browser/prepare.py @@ -0,0 +1,85 @@ +"""Build and vendor a static site; internet is needed only during this step. + +First run npm ci --prefix browser, then python browser/prepare.py. +""" + +import hashlib +import json +import shutil +import time +from concurrent.futures import ThreadPoolExecutor +from pathlib import Path +from urllib.request import urlopen + +from build import build + +HERE = Path(__file__).resolve().parent +SITE = HERE / "site" +VENDOR = SITE / "vendor" + + +def download(url, target, sha256): + if target.exists() and hashlib.sha256(target.read_bytes()).hexdigest() == sha256: + return + print(f"Downloading {target.name}", flush=True) + for attempt in range(3): + try: + with urlopen(url, timeout=60) as response: + data = response.read() + break + except OSError: + if attempt == 2: + raise + time.sleep(1 + attempt) + if hashlib.sha256(data).hexdigest() != sha256: + raise ValueError(f"Checksum mismatch: {url}") + target.write_bytes(data) + + +def main(): + runtime = HERE / "node_modules/pyodide" + if not runtime.exists(): + raise SystemExit("First run npm ci --prefix browser") + VENDOR.mkdir(parents=True, exist_ok=True) + for source in runtime.iterdir(): + if source.suffix in (".js", ".mjs", ".wasm", ".zip", ".json"): + shutil.copy2(source, VENDOR / source.name) + version = json.loads((runtime / "package.json").read_text())["version"] + lock = json.loads((runtime / "pyodide-lock.json").read_text())["packages"] + needed = set() + + def include(name): + if name not in needed: + needed.add(name) + for dependency in lock[name]["depends"]: + include(dependency) + + for name in ("numpy", "scipy", "micropip", "pytest"): + include(name) + downloads = [] + for name in sorted(needed): + package = lock[name] + filename = package["file_name"] + downloads.append( + (f"https://cdn.jsdelivr.net/pyodide/v{version}/full/{filename}", VENDOR / filename, package["sha256"]) + ) + extras = [] + for name, pinned in (("cunumpy", "0.1.4"), ("array-api-compat", "1.15.0")): + with urlopen(f"https://pypi.org/pypi/{name}/{pinned}/json", timeout=60) as response: + metadata = json.load(response) + wheel = next(p for p in metadata["urls"] if p["filename"].endswith("py3-none-any.whl")) + downloads.append((wheel["url"], VENDOR / wheel["filename"], wheel["digests"]["sha256"])) + extras.append("vendor/" + wheel["filename"]) + with ThreadPoolExecutor(max_workers=4) as pool: + list(pool.map(lambda item: download(*item), downloads)) + wheel = build(HERE / "dist") + shutil.copy2(wheel, SITE / wheel.name) + shutil.copy2(HERE / "tests/test_geometry.py", SITE / "test_geometry.py") + (SITE / "wheel.json").write_text(json.dumps({"wheel": wheel.name, "dependencies": extras, "pyodide": version})) + size = sum(p.stat().st_size for p in SITE.rglob("*") if p.is_file()) + print(f"Static site ready: {SITE} ({size / 1024**2:.1f} MiB, including test dependencies)") + print("Serve: python -m http.server 8765 --directory browser/site") + + +if __name__ == "__main__": + main() diff --git a/browser/site/app.js b/browser/site/app.js new file mode 100644 index 000000000..d2ab962a7 --- /dev/null +++ b/browser/site/app.js @@ -0,0 +1,108 @@ +const $ = id => document.getElementById(id); +const worker = new Worker("./worker.js", { type: "module" }); +const pending = new Map(); +let sequence = 0, catalog, eqdskText = "", busy = false, hasGeometry = false; +const status = (text, error = false) => { $("status").textContent = text; $("status").classList.toggle("error", error); }; +function request(action, payload = {}) { + return new Promise((resolve, reject) => { + const id = ++sequence; pending.set(id, { resolve, reject }); worker.postMessage({ id, action, ...payload }); + }); +} +function setBusy(value) { + busy = value; + for (const id of ["construct", "domain", "params", "point", "resolution", "archive", "eqdsk", "clear-flux"]) $(id).disabled = value; + $("save").disabled = value || !hasGeometry; +} +worker.onmessage = async ({ data }) => { + if (data.type === "fatal") { status(data.error, true); setBusy(true); return; } + if (data.type === "ready") { + catalog = data.catalog; + for (const name of Object.keys(catalog).sort()) $("domain").add(new Option(name, name)); + $("domain").value = "IGAPolarTorus"; + chooseDomain(); setBusy(false); await construct(); return; + } + const promise = pending.get(data.id); + if (promise) { pending.delete(data.id); data.error ? promise.reject(new Error(data.error)) : promise.resolve(data.result); } +}; +worker.onerror = error => { status(error.message, true); setBusy(true); }; +function chooseDomain() { + $("params").value = JSON.stringify(catalog[$("domain").value], null, 2); + $("flux-controls").hidden = $("domain").value !== "Tokamak"; +} +$("domain").addEventListener("change", chooseDomain); +$("eqdsk").addEventListener("change", async event => { + const file = event.target.files[0]; + eqdskText = file ? await file.text() : ""; + $("flux-note").textContent = file ? `Using ${file.name}` : "Using analytic circular flux surfaces."; +}); +$("clear-flux").addEventListener("click", () => { eqdskText = ""; $("eqdsk").value = ""; $("flux-note").textContent = "Using analytic circular flux surfaces."; }); +function probeOptions() { return { resolution: Number($("resolution").value), point: $("point").value.split(",").map(Number) }; } +async function construct() { + if (busy) return; + setBusy(true); status("Constructing geometry…"); + try { + const result = await request("render", { name: $("domain").value, params: JSON.parse($("params").value), eqdsk: eqdskText, ...probeOptions() }); + display(result); + } catch (error) { status(error.message, true); } finally { setBusy(false); } +} +$("controls").addEventListener("submit", event => { event.preventDefault(); construct(); }); +$("save").addEventListener("click", async () => { + setBusy(true); + try { + const data = await request("save"); + const bytes = Uint8Array.from(atob(data), c => c.charCodeAt(0)); + const url = URL.createObjectURL(new Blob([bytes], { type: "application/octet-stream" })); + const link = document.createElement("a"); link.href = url; link.download = "geometry.npz"; link.click(); + setTimeout(() => URL.revokeObjectURL(url), 1000); + } catch (error) { status(error.message, true); } finally { setBusy(false); } +}); +$("archive").addEventListener("change", async event => { + const file = event.target.files[0]; if (!file) return; + setBusy(true); status("Loading geometry…"); + try { + const bytes = new Uint8Array(await file.arrayBuffer()); + let binary = ""; for (const byte of bytes) binary += String.fromCharCode(byte); + display(await request("load", { data: btoa(binary), ...probeOptions() })); + } catch (error) { status(error.message, true); } finally { event.target.value = ""; setBusy(false); } +}); +const format = x => Number(x).toPrecision(6); +function display(result) { + hasGeometry = true; surface = result.xyz; draw(); + $("diagnostics").hidden = false; + $("coordinates").textContent = "F(η) = [" + result.mapped.map(format).join(", ") + "]"; + for (const key of ["jacobian", "metric"]) $(key).textContent = result[key].map(row => row.map(x => format(x).padStart(12)).join(" ")).join("\n"); + $("determinant").textContent = "det(DF) = " + format(result.determinant); + status(`${result.name} · surface and probe evaluated in ${result.elapsed.toFixed(3)} s`); +} + +// Local canvas rendering; all mapping construction and evaluation runs in Python. +const canvas = $("view"), ctx = canvas.getContext("2d"); +let surface, yaw = -.6, pitch = .65, zoom = 1, dragging; +function draw() { + const ratio = devicePixelRatio || 1; + canvas.width = canvas.clientWidth * ratio; canvas.height = canvas.clientHeight * ratio; + ctx.scale(ratio, ratio); const width = canvas.clientWidth, height = canvas.clientHeight; + ctx.clearRect(0, 0, width, height); if (!surface) return; + const n = surface[0].length, m = surface[0][0].length; + const bounds = surface.map(a => { const flat = a.flat(); return [Math.min(...flat), Math.max(...flat)]; }); + const center = bounds.map(b => (b[0] + b[1]) / 2); + const extent = Math.max(...bounds.map(b => b[1] - b[0]), 1e-8); + const scale = Math.min(width, height) * .64 * zoom / extent; + const points = Array.from({ length: n }, (_, i) => Array.from({ length: m }, (_, j) => { + const [x,y,z] = surface.map((a,k) => a[i][j] - center[k]); + const u = Math.cos(yaw)*x - Math.sin(yaw)*y, v = Math.sin(yaw)*x + Math.cos(yaw)*y; + return [width/2+scale*u, height/2-scale*(Math.cos(pitch)*z-Math.sin(pitch)*v), Math.sin(pitch)*z+Math.cos(pitch)*v]; + })); + const faces = []; + for (let i=0;i a.reduce((s,p)=>s+p[2],0)-b.reduce((s,p)=>s+p[2],0)); + for (const face of faces) { + ctx.beginPath(); face.forEach((p,i) => i ? ctx.lineTo(p[0],p[1]) : ctx.moveTo(p[0],p[1])); ctx.closePath(); + ctx.fillStyle = "rgba(119,180,195,.86)"; ctx.fill(); ctx.strokeStyle = "rgba(25,80,107,.56)"; ctx.lineWidth = .65; ctx.stroke(); + } +} +canvas.addEventListener("pointerdown", event => { dragging = [event.clientX,event.clientY]; canvas.setPointerCapture(event.pointerId); }); +canvas.addEventListener("pointermove", event => { if (!dragging) return; yaw += (event.clientX-dragging[0])*.008; pitch += (event.clientY-dragging[1])*.008; dragging=[event.clientX,event.clientY]; draw(); }); +for (const event of ["pointerup","pointercancel"]) canvas.addEventListener(event,()=>dragging=null); +canvas.addEventListener("wheel", event => { event.preventDefault(); zoom=Math.max(.3,Math.min(4,zoom*Math.exp(-event.deltaY*.001))); draw(); },{passive:false}); +new ResizeObserver(draw).observe(canvas); diff --git a/browser/site/index.html b/browser/site/index.html new file mode 100644 index 000000000..663403841 --- /dev/null +++ b/browser/site/index.html @@ -0,0 +1,36 @@ + + + + + + Struphy Geometry + + + +
STRUPHY / GEOMETRY

Explore mapped domains.

Construct, inspect, and save geometry entirely in your browser.

+
+
+ + + + + + +
+

Files stay in this browser. Saved splines include their control points, so the original equilibrium is not needed to reopen them.

+
+
+
Drag to rotate · scroll to zoom
+

Surface at η₁ = 1 · the outer radial surface for cylinders and tori

+

Loading Python and numerical packages…

+ +
+
+
Python / cuNumPy / SciPy / Pyodide · no computation server
+ + + diff --git a/browser/site/style.css b/browser/site/style.css new file mode 100644 index 000000000..735b9569f --- /dev/null +++ b/browser/site/style.css @@ -0,0 +1,7 @@ +:root { font-family: Inter, system-ui, sans-serif; color: #163047; background: #f3f6f8; font-synthesis: none; } +* { box-sizing: border-box; } body { margin: 0; } header, main, footer { max-width: 1320px; margin: auto; padding: 28px 36px; } +header { padding-bottom: 18px; } .eyebrow { color: #35728c; font-size: 12px; letter-spacing: .15em; font-weight: 700; } h1 { font-size: clamp(28px, 4vw, 42px); letter-spacing: -.04em; margin: 12px 0; } header p { color: #526878; margin: 0; } +main { display: grid; grid-template-columns: 310px minmax(0, 1fr); gap: 30px; padding-top: 8px; } form { background: white; padding: 22px; border: 1px solid #dce5eb; border-radius: 12px; align-self: start; } label { display: block; font-size: 12px; font-weight: 650; margin: 0 0 7px; } input, select, textarea { width: 100%; padding: 10px; border: 1px solid #c6d5df; border-radius: 6px; color: #173047; background: #fff; margin-bottom: 18px; } textarea { font: 12px/1.6 ui-monospace, monospace; resize: vertical; } input:focus, select:focus, textarea:focus { outline: 2px solid #68acc4; outline-offset: 1px; } +button { padding: 11px 13px; border: 1px solid #c6d5df; border-radius: 6px; font: inherit; font-size: 12px; background: #fff; color: #173047; cursor: pointer; } button:disabled { opacity: .5; cursor: default; } #construct { width: 100%; background: #17546a; color: white; border-color: #17546a; font-weight: 600; } .archive { display: flex; gap: 8px; margin: 13px 0; } .upload { margin: 0; padding: 11px 10px; border: 1px solid #c6d5df; border-radius: 6px; cursor: pointer; } .upload input { width: 1px; height: 1px; position: absolute; opacity: 0; } .upload:focus-within { outline: 2px solid #68acc4; } .note, #flux-note { font-size: 11px; color: #617582; line-height: 1.6; } #clear-flux { margin-bottom: 15px; } +.viewer { height: 450px; min-height: 280px; background: #e5edf2; border: 1px solid #d3dfe7; border-radius: 12px; position: relative; overflow: hidden; } canvas { width: 100%; height: 100%; display: block; touch-action: none; cursor: grab; } .hint { position: absolute; bottom: 16px; left: 18px; color: #527084; font-size: 11px; pointer-events: none; } .caption { font-size: 11px; color: #607789; } #status { padding: 12px 15px; background: #e0ecec; border-radius: 6px; font-size: 13px; overflow-wrap: anywhere; } #status.error { background: #fbe8e4; color: #923a28; } h2 { font-size: 16px; } h3 { font-size: 12px; margin-bottom: 8px; color: #577184; } .matrices { display: flex; gap: 35px; flex-wrap: wrap; } pre { font-size: 12px; line-height: 1.8; margin: 0; } #coordinates, #determinant { font: 12px/1.8 ui-monospace, monospace; } footer { padding-top: 0; font-size: 11px; color: #6e8493; } +@media (max-width: 760px) { header, main, footer { padding-left: 18px; padding-right: 18px; } main { grid-template-columns: 1fr; } .viewer { height: 340px; } textarea { height: 150px; } } diff --git a/browser/site/worker.js b/browser/site/worker.js new file mode 100644 index 000000000..30196aca7 --- /dev/null +++ b/browser/site/worker.js @@ -0,0 +1,98 @@ +import { loadPyodide } from "./vendor/pyodide.mjs"; + +let pyodide; +const ready = (async () => { + const config = await (await fetch("./wheel.json")).json(); + pyodide = await loadPyodide({ indexURL: new URL("./vendor/", import.meta.url).href }); + await pyodide.loadPackage(["numpy", "scipy", "micropip"]); + const micropip = pyodide.pyimport("micropip"); + try { + for (const wheel of [...config.dependencies].reverse()) { + await micropip.install(new URL(wheel, import.meta.url).href); + } + await micropip.install(new URL(config.wheel, import.meta.url).href); + } finally { micropip.destroy(); } + pyodide.runPython(` +import base64, inspect, json, time +import numpy as np +from struphy_geometry import domains, CircularFlux, EQDSKFlux, sample_surface, save_geometry, load_geometry +current_domain = None + +def catalog(): + result = {} + for name, cls in vars(domains).items(): + if inspect.isclass(cls) and cls.__module__ == domains.__name__: + result[name] = {k: v.default for k, v in inspect.signature(cls).parameters.items() if k != "equilibrium"} + return result + +def snapshot(domain, request): + resolution = int(request.get("resolution", 24)) + if not 4 <= resolution <= 64: + raise ValueError("Resolution must be between 4 and 64") + point = request.get("point", [.37, .23, .41]) + if np.shape(point) != (3,) or not np.isfinite(point).all() or not np.all((np.array(point) > 0) & (np.array(point) < 1)): + raise ValueError("Probe coordinates must be three finite values strictly between 0 and 1") + started = time.perf_counter() + xyz = sample_surface(domain, resolution) + mapped = domain(*point, squeeze_out=True) + jacobian = domain.jacobian(*point, squeeze_out=True) + metric = domain.metric(*point, squeeze_out=True) + determinant = float(domain.jacobian_det(*point, squeeze_out=True)) + if not all(np.isfinite(x).all() for x in (*xyz, mapped, jacobian, metric, determinant)): + raise ValueError("The mapping is singular or non-finite at the requested points") + return dict(name=type(domain).__name__, xyz=[a.tolist() for a in xyz], mapped=mapped.tolist(), + jacobian=jacobian.tolist(), metric=metric.tolist(), determinant=determinant, + elapsed=time.perf_counter()-started) + +def handle(request): + global current_domain + action = request["action"] + if action == "render": + name = request["name"] + if name not in catalog(): + raise ValueError("Unknown domain") + params = request["params"] + if name == "Tokamak": + params["equilibrium"] = EQDSKFlux.from_text(request["eqdsk"]) if request.get("eqdsk") else CircularFlux() + candidate = getattr(domains, name)(**params) + result = snapshot(candidate, request) + current_domain = candidate + return result + if action == "load": + candidate = load_geometry(base64.b64decode(request["data"])) + result = snapshot(candidate, request) + current_domain = candidate + return result + if action == "save": + if current_domain is None: + raise ValueError("Construct a domain first") + return base64.b64encode(save_geometry(current_domain)).decode() + raise ValueError("Unknown action") +`); + postMessage({ type: "ready", catalog: JSON.parse(pyodide.runPython("json.dumps(catalog())")) }); +})(); +ready.catch(error => postMessage({ type: "fatal", error: String(error) })); + +// Queue requests: each operation owns the Python globals until it completes. +let queue = ready; +self.onmessage = ({ data }) => { + queue = queue.catch(() => {}).then(async () => { + try { + if (data.action === "tests") { + await ready; + await pyodide.loadPackage("pytest"); + pyodide.FS.writeFile("/test_geometry.py", await (await fetch("./test_geometry.py")).text()); + const output = []; + pyodide.setStdout({ batched: line => output.push(line) }); + pyodide.setStderr({ batched: line => output.push(line) }); + const code = pyodide.runPython(`import os, pytest\nos.environ['PYTEST_DISABLE_PLUGIN_AUTOLOAD'] = '1'\nint(pytest.main(['-q', '/test_geometry.py', '-p', 'no:cacheprovider']))`); + postMessage({ id: data.id, result: { code, output: output.join("\n") } }); + } else { + await ready; + pyodide.globals.set("request_json", JSON.stringify(data)); + const result = JSON.parse(pyodide.runPython("json.dumps(handle(json.loads(request_json)), allow_nan=False)")); + postMessage({ id: data.id, result }); + } + } catch (error) { postMessage({ id: data.id, error: String(error) }); } + }); +}; diff --git a/browser/tests/browser.spec.js b/browser/tests/browser.spec.js new file mode 100644 index 000000000..a140b0667 --- /dev/null +++ b/browser/tests/browser.spec.js @@ -0,0 +1,62 @@ +import { test, expect } from "@playwright/test"; +import { fileURLToPath } from "node:url"; + +// Block external requests: the prepared site must work with only static assets. +test.beforeEach(async ({ context }) => { + await context.route("**/*", route => { + const url = new URL(route.request().url()); + return url.hostname === "127.0.0.1" || url.protocol === "blob:" ? route.continue() : route.abort(); + }); +}); + +test("numerical suite in real browser Pyodide", async ({ page }) => { + await page.goto("/"); + const result = await page.evaluate(() => new Promise((resolve, reject) => { + const worker = new Worker("./worker.js", { type: "module" }); + worker.onerror = error => { worker.terminate(); reject(new Error(error.message)); }; + worker.onmessage = ({ data }) => { + if (data.type === "fatal" || data.error) { worker.terminate(); reject(new Error(data.error)); } + else if (data.type === "ready") worker.postMessage({ id: 1, action: "tests" }); + else if (data.id === 1) { worker.terminate(); resolve(data.result); } + }; + })); + console.log(result.output); + expect(result.code, result.output).toBe(0); +}); + +test("construct, change parameters, report errors, save and reopen", async ({ page }) => { + const errors = []; + page.on("pageerror", error => errors.push(error.message)); + await page.goto("/"); + await expect(page.getByRole("status")).toContainText("IGAPolarTorus ·", { timeout: 120000 }); + await expect(page.locator("#diagnostics")).toBeVisible(); + await page.locator("#domain").selectOption("PoweredEllipticCylinder"); + await page.getByRole("button", { name: "Construct domain" }).click(); + await expect(page.getByRole("status")).toContainText("PoweredEllipticCylinder ·"); + await page.locator("#domain").selectOption("Tokamak"); + await page.locator("#params").fill(JSON.stringify({ num_elements: [4, 12], degree: [2, 3], xi_param: "equal_angle" })); + await page.getByRole("button", { name: "Construct domain" }).click(); + await expect(page.getByRole("status")).toContainText("Tokamak ·"); + await page.locator("#eqdsk").setInputFiles(fileURLToPath(new URL( + "../../src/struphy/fields_background/mhd_equil/eqdsk/data/AUGNLED_g031213.00830.high", import.meta.url, + ))); + await expect(page.locator("#flux-note")).toContainText("AUGNLED"); + await page.getByRole("button", { name: "Construct domain" }).click(); + await expect(page.getByRole("status")).toContainText("Tokamak ·", { timeout: 90000 }); + const before = await page.locator("#coordinates").textContent(); + const downloaded = page.waitForEvent("download"); + await page.getByRole("button", { name: "Save geometry" }).click(); + const download = await downloaded; + await page.locator("#archive").setInputFiles(await download.path()); + await expect(page.getByRole("status")).toContainText("PoloidalSplineTorus ·"); + await expect(page.locator("#coordinates")).toHaveText(before); + await page.locator("#params").fill("invalid json"); + await page.getByRole("button", { name: "Construct domain" }).click(); + await expect(page.locator("#status")).toHaveClass("error"); + await expect(page.getByRole("button", { name: "Construct domain" })).toBeEnabled(); + await page.locator("#domain").selectOption("Cuboid"); + await page.getByRole("button", { name: "Construct domain" }).click(); + await expect(page.getByRole("status")).toContainText("Cuboid ·"); + await page.screenshot({ path: "test-results/geometry.png", fullPage: true }); + expect(errors).toEqual([]); +}); diff --git a/browser/tests/conftest.py b/browser/tests/conftest.py new file mode 100644 index 000000000..77bbee578 --- /dev/null +++ b/browser/tests/conftest.py @@ -0,0 +1,10 @@ +"""Test the built wheel, not the full Struphy installation.""" + +import sys +from pathlib import Path + +if sys.platform != "emscripten": + wheels = sorted((Path(__file__).resolve().parents[1] / "dist").glob("*.whl")) + if len(wheels) != 1: + raise RuntimeError("Run python browser/build.py and keep exactly one wheel in browser/dist") + sys.path.insert(0, str(wheels[0])) diff --git a/browser/tests/test_geometry.py b/browser/tests/test_geometry.py new file mode 100644 index 000000000..5054f2fb7 --- /dev/null +++ b/browser/tests/test_geometry.py @@ -0,0 +1,187 @@ +"""Numerical regression tests shared by native Python and browser Pyodide.""" + +import io +import json +import sys + +import numpy as np +import pytest +from struphy_geometry import ( + CircularFlux, + EQDSKFlux, + Spline, + domains, + interp_mapping, + load_geometry, + sample_surface, + save_geometry, +) + +ANALYTIC = [ + "Cuboid", + "Orthogonal", + "Colella", + "HollowCylinder", + "PoweredEllipticCylinder", + "HollowTorus", + "ShafranovShiftCylinder", + "ShafranovSqrtCylinder", + "ShafranovDshapedCylinder", +] +MAPPINGS = ANALYTIC + ["IGAPolarCylinder", "IGAPolarTorus"] + + +@pytest.mark.parametrize("name", MAPPINGS) +def test_derivatives_and_shapes(name): + d = getattr(domains, name)() + q = np.array([0.37, 0.23, 0.41]) + J = d.jacobian(*q, squeeze_out=True) + h = 1e-6 + fd = np.column_stack( + [ + (d(*(q + np.eye(3)[i] * h), squeeze_out=True) - d(*(q - np.eye(3)[i] * h), squeeze_out=True)) / (2 * h) + for i in range(3) + ] + ) + np.testing.assert_allclose(J, fd, rtol=1e-5, atol=1e-7) + np.testing.assert_allclose(d.metric(*q, squeeze_out=True), J.T @ J, atol=1e-10) + np.testing.assert_allclose(d.jacobian_inv(*q, squeeze_out=True) @ J, np.eye(3), atol=1e-10) + np.testing.assert_allclose(d.jacobian_det(*q, squeeze_out=True), np.linalg.det(J), rtol=1e-10) + np.testing.assert_allclose( + d.metric_inv(*q, squeeze_out=True) @ d.metric(*q, squeeze_out=True), np.eye(3), atol=1e-10 + ) + assert d(np.array([q, q])).shape == (3, 2) + assert d(np.array([0.2, 0.4]), np.array([0.2, 0.4]), np.array([0.2, 0.4])).shape == (3, 2, 2, 2) + + +@pytest.mark.parametrize("name", MAPPINGS) +def test_archive_roundtrip(name): + d = getattr(domains, name)() + restored = load_geometry(save_geometry(d)) + q = np.random.default_rng(17).uniform(0.1, 0.9, (8, 3)) + np.testing.assert_allclose(restored(q), d(q), atol=1e-12) + np.testing.assert_allclose(restored.jacobian(q), d.jacobian(q), atol=1e-12) + if name in ANALYTIC: + np.testing.assert_allclose(d.from_dict(d.to_dict())(q), d(q)) + + +@pytest.mark.parametrize("kind", ["0", "1", "2", "3", "v"]) +def test_pull_push(kind): + d = domains.HollowTorus(tor_period=1) + q = (0.31, 0.22, 0.43) + a = 2.0 if kind in ("0", "3") else [1.0, 2.0, 3.0] + pulled = d.pull(a, *q, kind=kind, squeeze_out=True) + pushed = d.push(pulled, *q, kind=kind, squeeze_out=True) + np.testing.assert_allclose(pushed, a, atol=1e-10) + + +@pytest.mark.parametrize( + "forward,backward,a", + [ + ("0_to_3", "3_to_0", 2.0), + ("1_to_2", "2_to_1", [1.0, 2.0, 3.0]), + ("v_to_1", "1_to_v", [1.0, 2.0, 3.0]), + ("v_to_2", "2_to_v", [1.0, 2.0, 3.0]), + ], +) +def test_transform(forward, backward, a): + d = domains.IGAPolarTorus() + q = (0.31, 0.22, 0.43) + result = d.transform(a, *q, kind=forward, squeeze_out=True) + np.testing.assert_allclose(d.transform(result, *q, kind=backward, squeeze_out=True), a, atol=1e-10) + + +def test_spline_3d_construction(): + counts, degree, periodic = (3, 4, 3), (2, 2, 2), (False, False, False) + + def component(i): + def f(*etas): + return np.meshgrid(*etas, indexing="ij")[i] + + return f + + coefficients = interp_mapping(counts, degree, periodic, *[component(i) for i in range(3)]) + d = Spline(counts, degree, periodic, *coefficients) + q = np.array([[0.2, 0.3, 0.4], [0.6, 0.7, 0.8]]) + np.testing.assert_allclose(d(q), q.T, atol=1e-12) + np.testing.assert_allclose(load_geometry(save_geometry(d))(q), q.T, atol=1e-12) + + +@pytest.mark.parametrize("mode", ["equal_angle", "equal_arc_length", "sfl", "equal_area", "equal_volume"]) +def test_flux_tracing(mode): + flux = CircularFlux() + d = domains.Tokamak( + equilibrium=flux, + num_elements=(4, 12), + degree=(2, 3), + xi_param=mode, + num_elements_pre=(6, 24), + p_pre=(3, 3), + psi_power=0.5, + psi_shifts=(4.0, 10.0), + tor_period=1, + ) + q = np.array([[0.2, 0.13, 0.2], [0.7, 0.73, 0.4]]) + xyz = d(q) + expected = 0.04 + q[:, 0] ** 2 * (0.9 - 0.04) + np.testing.assert_allclose(flux.psi(np.hypot(xyz[0], xyz[1]), xyz[2]), expected, rtol=0.04, atol=0.002) + restored = load_geometry(save_geometry(d)) + np.testing.assert_allclose(restored(q), xyz, atol=1e-12) + np.testing.assert_allclose(restored.jacobian(q), d.jacobian(q), atol=1e-12) + + +def eqdsk_text(sign=1): + n = 9 + R, Z = np.linspace(1.0, 5.0, n), np.linspace(-2.0, 2.0, n) + flux = sign * ((R[None, :] - 3.0) ** 2 + Z[:, None] ** 2) + header = [4.0, 4.0, 3.0, 1.0, 0.0, 3.0, 0.0, 0.0, sign, 1.0, *([0.0] * 10)] + values = [*header, *([0.0] * (4 * n)), *flux.ravel(), *([0.0] * n)] + return f"Synthetic circular flux 0 {n} {n}\n" + "\n".join( + "".join(f"{x:16.9E}".replace("E", "D") for x in values[i : i + 5]) for i in range(0, len(values), 5) + ) + + +@pytest.mark.parametrize("sign", [1, -1]) +def test_eqdsk(sign): + flux = EQDSKFlux.from_text(eqdsk_text(sign)) + np.testing.assert_allclose(flux.psi_axis_RZ, [3.0, 0.0], atol=1e-10) + np.testing.assert_allclose(flux.psi([3.2, 3.4], 0.3), sign * np.array([0.13, 0.25]), atol=1e-10) + np.testing.assert_allclose(flux.psi(3.2, 0.3, dR=1), sign * 0.4, atol=1e-10) + d = domains.Tokamak(equilibrium=flux, num_elements=(4, 12)) + assert np.isfinite(d(0.5, 0.3, 0.1)).all() + + +def test_seams_poles_and_outside(): + d = domains.IGAPolarTorus(tor_period=1) + np.testing.assert_allclose(d(0.5, 0.0, 0.3), d(0.5, 1.0, 0.3), atol=1e-12) + np.testing.assert_allclose(d(0.5, 0.3, 0.0), d(0.5, 0.3, 1.0), atol=1e-12) + np.testing.assert_allclose(d(0.0, 0.2, 0.3), d(0.0, 0.7, 0.3), atol=1e-12) + q = np.array([[0.2, 0.3, 0.4], [-0.1, 0.3, 0.4], [1.1, 0.3, 0.4]]) + assert d(q.copy(), remove_outside=True).shape == (3, 1) + outside = d(q.copy(), remove_outside=False) + np.testing.assert_array_equal(outside[:, 1:], -1.0) + assert sample_surface(d, 8)[0].shape == (8, 8) + + +def test_invalid_inputs_and_excluded_adapters(): + with pytest.raises(ValueError, match="header"): + EQDSKFlux.from_text("bad") + with pytest.raises(ValueError, match="Incomplete"): + EQDSKFlux.from_text("bad 0 9 9\n") + with pytest.raises(ValueError): + CircularFlux(a=-1) + with pytest.raises(ValueError, match="explicit"): + Spline() + with pytest.raises(ValueError, match="equilibrium"): + domains.Tokamak() + assert not hasattr(domains, "GVECunit") + assert not hasattr(domains, "DESCunit") + stream = io.BytesIO() + np.savez(stream, metadata=json.dumps({"format": "struphy-geometry", "version": 99})) + with pytest.raises(ValueError, match="version"): + load_geometry(stream.getvalue()) + + +def test_no_native_runtime_dependencies(): + forbidden = {"struphy", "pyccel", "feectools", "mpi4py", "cupy", "vtk", "pyvista", "h5py", "scope_profiler"} + assert not forbidden.intersection(name.split(".")[0] for name in sys.modules) diff --git a/doc/sections/domains.rst b/doc/sections/domains.rst index 47fd3bc95..2e878d2ad 100644 --- a/doc/sections/domains.rst +++ b/doc/sections/domains.rst @@ -7,6 +7,11 @@ Struphy models are implemented in curvilinear coordinates and can be run on a variaty of mapped domains. Besides analytical mappings, there are discrete spline mappings available (IGA approach). +A standalone ``struphy-geometry`` wheel and a static browser viewer are available +for running domains, spline interpolation, and flux tracing entirely in Pyodide. +See the `browser geometry guide `_ +for building the package, loading EQDSK files, and exchanging portable geometry archives. + The (physical) domain :math:`\Omega \subset \mathbb R^3` is an open subset of :math:`\mathbb R^3`, defined by a diffeomorphism diff --git a/pyproject.toml b/pyproject.toml index 3b0b039ff..7d47bc984 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,7 +24,7 @@ classifiers = [ ] dependencies = [ "numpy<=2.5.0", - "cunumpy>=0.1.4, <=0.1.5", + "cunumpy>=0.1.5, <=0.1.5", "pyccel>=2.2.0, <=2.2.3", "feectools >= 0.1.10, <=0.1.10", "scipy<=1.18.0", diff --git a/src/struphy/geometry/base.py b/src/struphy/geometry/base.py index ed3c8177c..ab95c4caa 100644 --- a/src/struphy/geometry/base.py +++ b/src/struphy/geometry/base.py @@ -7,10 +7,8 @@ from abc import ABCMeta, abstractmethod import cunumpy as xp -import h5py import numpy as np from cunumpy import PyccelKernel -from pyvista import Plotter, StructuredGrid from scipy.sparse import csc_matrix, kron from scipy.sparse.linalg import splu, spsolve @@ -29,8 +27,12 @@ def display(*objects, **kwargs): from struphy.geometry import evaluation_kernels, transform_kernels from struphy.kernel_arguments.pusher_args_kernels import DomainArguments from struphy.linear_algebra import linalg_kron +from struphy.utils.class_helpers import ( + __class_with_params_repr_no_defaults__, + all_class_params_are_default, + all_subclasses, +) from struphy.utils.docstring_converter import rst_to_html, rst_to_latex, rst_to_markdown -from struphy.utils.utils import __class_with_params_repr_no_defaults__, all_class_params_are_default, all_subclasses logger = logging.getLogger("struphy") @@ -1706,6 +1708,8 @@ def create_geometry_mesh( pyvista.StructuredGrid """ + from pyvista import StructuredGrid + grids_log = [ xp.linspace(1e-6, 1.0, nx), xp.linspace(0.0, 1.0, ny), @@ -1727,6 +1731,8 @@ def show_3d( nz: int = 32, ): """Show the 3D geometry using PyVista.""" + from pyvista import Plotter + mesh = self.create_geometry_mesh(nx, ny, nz) plotter = Plotter() plotter.add_mesh(mesh, show_edges=True) @@ -2161,7 +2167,6 @@ def to_dict(self) -> dict: @classmethod def from_dict(cls, dct): - from struphy.fields_background.base import FluidEquilibrium from struphy.geometry.utilities import get_domain_by_name name = dct["type"] @@ -2169,6 +2174,8 @@ def from_dict(cls, dct): params = dict(dct["params"]) for key, value in params.items(): if isinstance(value, dict) and "type" in value and "params" in value: + from struphy.fields_background.base import FluidEquilibrium + params[key] = FluidEquilibrium.from_dict(value) return domain_cls(**params) diff --git a/src/struphy/geometry/domains.py b/src/struphy/geometry/domains.py index 3b5c50296..3da7a406d 100644 --- a/src/struphy/geometry/domains.py +++ b/src/struphy/geometry/domains.py @@ -9,11 +9,10 @@ """ import copy +from typing import TYPE_CHECKING import cunumpy as xp -from struphy.fields_background.base import AxisymmMHDequilibrium -from struphy.fields_background.equils import EQDSKequilibrium from struphy.geometry.base import ( Domain, PoloidalSplineStraight, @@ -23,6 +22,9 @@ ) from struphy.geometry.utilities import field_line_tracing +if TYPE_CHECKING: + from struphy.fields_background.base import AxisymmMHDequilibrium + class Tokamak(PoloidalSplineTorus): r"""Mappings for Tokamak MHD equilibria constructed via :ref:`field-line tracing ` of a poloidal flux function :math:`\psi`. @@ -69,7 +71,7 @@ def doc_mapping(cls): def __init__( self, - equilibrium: AxisymmMHDequilibrium = None, + equilibrium: "AxisymmMHDequilibrium" = None, num_elements: tuple = (8, 32), degree: tuple = (2, 3), psi_power: float = 0.75, @@ -84,9 +86,16 @@ def __init__( if r_min != 0.0: r0 = r_min if equilibrium is None: + from struphy.fields_background.equils import EQDSKequilibrium + equilibrium = EQDSKequilibrium() else: - assert isinstance(equilibrium, AxisymmMHDequilibrium) + # Geometry only needs a flux function and its axis/boundary metadata. + # This also accepts lightweight flux providers without the MHD framework. + if not callable(getattr(equilibrium, "psi", None)) or not all( + hasattr(equilibrium, name) for name in ("psi_range", "psi_axis_RZ") + ): + raise TypeError("equilibrium must provide psi, psi_range, and psi_axis_RZ") # use the params setter self.params = copy.deepcopy(locals()) diff --git a/src/struphy/geometry/flux.py b/src/struphy/geometry/flux.py new file mode 100644 index 000000000..20cc71176 --- /dev/null +++ b/src/struphy/geometry/flux.py @@ -0,0 +1,104 @@ +"""Small axisymmetric flux providers for constructing geometry without an MHD solver.""" + +import re +from pathlib import Path + +import numpy as np +from scipy.interpolate import RectBivariateSpline +from scipy.optimize import root + + +class CircularFlux: + """Analytic circular flux surfaces: psi = ((R - R0)**2 + (Z - Z0)**2) / a**2. + + Implements the flux interface accepted by ``Tokamak``. This is a geometry + example, not a force-balanced MHD equilibrium. + """ + + def __init__(self, R0=3.0, Z0=0.0, a=1.0): + if not np.isfinite([R0, Z0, a]).all() or a <= 0 or R0 <= a: + raise ValueError("Require finite parameters and R0 > a > 0") + self.R0, self.Z0, self.a = float(R0), float(Z0), float(a) + self.psi_axis_RZ = (self.R0, self.Z0) + self.psi_range = (0.0, 1.0) + + def psi(self, R, Z, dR=0, dZ=0): + R, Z = np.broadcast_arrays(np.asarray(R), np.asarray(Z)) + if dR < 0 or dZ < 0 or int(dR) != dR or int(dZ) != dZ: + raise ValueError("Derivative orders must be non-negative integers") + if dR == dZ == 0: + out = ((R - self.R0) ** 2 + (Z - self.Z0) ** 2) / self.a**2 + elif (dR, dZ) == (1, 0): + out = 2 * (R - self.R0) / self.a**2 + elif (dR, dZ) == (0, 1): + out = 2 * (Z - self.Z0) / self.a**2 + elif (dR, dZ) in ((2, 0), (0, 2)): + out = np.full(R.shape, 2 / self.a**2) + else: + out = np.zeros(R.shape) + return out.item() if out.ndim == 0 else out + + +class EQDSKFlux: + """Poloidal flux from a G-EQDSK file, with no unit rescaling or smoothing. + + Use ``from_text`` with an uploaded browser file. The reader accepts E/D + scientific notation and reads the rectangular flux grid and axis metadata; + pressure, current, and boundary contour data are not interpreted. The + interpolated magnetic axis is refined by default, for either flux sign. + """ + + def __init__(self, filename, *, refine_axis=True): + self._read(Path(filename).read_text(), refine_axis=refine_axis) + + @classmethod + def from_text(cls, text, *, refine_axis=True): + obj = cls.__new__(cls) + obj._read(text, refine_axis=refine_axis) + return obj + + def _read(self, text, *, refine_axis): + lines = text.splitlines() + try: + nR, nZ = map(int, lines[0].split()[-2:]) + except (IndexError, ValueError) as exc: + raise ValueError("Invalid G-EQDSK header") from exc + if nR < 4 or nZ < 4: + raise ValueError("G-EQDSK flux grid must have at least four points per axis") + values = np.array( + [ + float(s.replace("D", "E").replace("d", "e")) + for s in re.findall(r"[+-]?(?:\d+\.\d*|\.\d+|\d+)[EeDd][+-]?\d+", "\n".join(lines[1:])) + ] + ) + required = 20 + 4 * nR + nR * nZ + if values.size < required or not np.isfinite(values[:required]).all(): + raise ValueError("Incomplete or non-finite G-EQDSK flux data") + rdim, zdim, _, rleft, zmid = values[:5] + Raxis, Zaxis, _, psi_edge = values[5:9] + if rdim <= 0 or zdim <= 0: + raise ValueError("G-EQDSK grid dimensions must be positive") + self.R = np.linspace(rleft, rleft + rdim, nR) + self.Z = np.linspace(zmid - zdim / 2, zmid + zdim / 2, nZ) + grid = values[20 + 4 * nR : required].reshape(nZ, nR).T + self._spline = RectBivariateSpline(self.R, self.Z, grid, kx=3, ky=3, s=0) + axis = np.array([Raxis, Zaxis]) + if refine_axis: + solution = root(lambda q: [self.psi(*q, dR=1), self.psi(*q, dZ=1)], axis) + if not solution.success: + raise ValueError("Could not locate the interpolated magnetic axis") + axis = solution.x + if not (self.R[0] < axis[0] < self.R[-1] and self.Z[0] < axis[1] < self.Z[-1]): + raise ValueError("Magnetic axis is outside the flux grid") + self.psi_axis_RZ = tuple(float(x) for x in axis) + self.psi_range = (float(self.psi(*axis)), float(psi_edge)) + if self.psi_range[0] == self.psi_range[1]: + raise ValueError("Axis and edge flux must differ") + + def psi(self, R, Z, dR=0, dZ=0): + """Evaluate flux or derivatives at broadcast-compatible R/Z coordinates.""" + R, Z = np.broadcast_arrays(np.asarray(R), np.asarray(Z)) + # FITPACK wrappers differ in scalar output shape across SciPy versions. + # Evaluate flat point lists and restore our explicit broadcasting contract. + out = np.asarray(self._spline.ev(R.ravel(), Z.ravel(), dx=dR, dy=dZ)).reshape(R.shape) + return out.item() if out.ndim == 0 else out diff --git a/src/struphy/geometry/mappings_kernels.py b/src/struphy/geometry/mappings_kernels.py index 62a99d1f0..8826241c6 100644 --- a/src/struphy/geometry/mappings_kernels.py +++ b/src/struphy/geometry/mappings_kernels.py @@ -764,10 +764,10 @@ def powered_ellipse_df( ): """Jacobian matrix for :meth:`struphy.geometry.mappings_kernels.powered_ellipse`.""" - df_out[0, 0] = (eta1 ** (s - 1)) * rx * cos(2 * pi * eta2) + df_out[0, 0] = s * (eta1 ** (s - 1)) * rx * cos(2 * pi * eta2) df_out[0, 1] = -2 * pi * (eta1**s) * rx * sin(2 * pi * eta2) df_out[0, 2] = 0.0 - df_out[1, 0] = (eta1 ** (s - 1)) * ry * sin(2 * pi * eta2) + df_out[1, 0] = s * (eta1 ** (s - 1)) * ry * sin(2 * pi * eta2) df_out[1, 1] = 2 * pi * (eta1**s) * ry * cos(2 * pi * eta2) df_out[1, 2] = 0.0 df_out[2, 0] = 0.0 diff --git a/src/struphy/geometry/portable.py b/src/struphy/geometry/portable.py new file mode 100644 index 000000000..86d71d59e --- /dev/null +++ b/src/struphy/geometry/portable.py @@ -0,0 +1,90 @@ +"""Array-only mesh sampling and portable, non-pickle geometry archives.""" + +import io +import json + +import numpy as np + +from struphy.geometry import domains +from struphy.geometry.base import PoloidalSplineStraight, PoloidalSplineTorus, Spline + +_ANALYTIC = ( + "Cuboid", + "Orthogonal", + "Colella", + "HollowCylinder", + "PoweredEllipticCylinder", + "HollowTorus", + "ShafranovShiftCylinder", + "ShafranovSqrtCylinder", + "ShafranovDshapedCylinder", +) + + +def sample_surface(domain, resolution=32, radial_coordinate=1.0): + """Return x/y/z arrays of shape (resolution, resolution) at fixed eta1. + + This samples a logical coordinate surface, not necessarily the entire domain + boundary. It is the outer radial surface for cylindrical/toroidal mappings. + Use the domain directly to sample other faces, slices, or volume grids. + """ + if not isinstance(resolution, (int, np.integer)) or resolution < 2: + raise ValueError("resolution must be an integer >= 2") + if not 0 <= radial_coordinate <= 1: + raise ValueError("radial_coordinate must lie in [0, 1]") + eta = np.linspace(0.0, 1.0, resolution) + xyz = domain(float(radial_coordinate), eta, eta, squeeze_out=True) + return tuple(np.asarray(component) for component in xyz) + + +def save_geometry(domain): + """Return NPZ bytes containing an evaluated mapping's definition. + + Spline archives contain control points, never equilibrium objects. Loading + GVEC/DESC/Tokamak-derived splines consequently needs no original solver or + input file. The reconstructed object is a generic spline of the same kind. + """ + metadata = {"format": "struphy-geometry", "version": 1} + arrays = {} + if domain.kind_map in (0, 1, 2): + metadata.update( + kind_map=int(domain.kind_map), + num_elements=[int(x) for x in domain.num_elements], + degree=[int(x) for x in domain.degree], + spl_kind=[bool(x) for x in domain.spl_kind], + ) + arrays = {key: np.asarray(getattr(domain, key)) for key in ("cx", "cy", "cz", "params_numpy")} + elif type(domain).__name__ in _ANALYTIC: + metadata.update(type=type(domain).__name__, params=domain.params) + else: + raise ValueError(f"Unsupported mapping type: {type(domain).__name__}") + stream = io.BytesIO() + np.savez_compressed(stream, metadata=np.array(json.dumps(metadata)), **arrays) + return stream.getvalue() + + +def load_geometry(data): + """Load a mapping from NPZ bytes produced by ``save_geometry`` (no pickle).""" + with np.load(io.BytesIO(data), allow_pickle=False) as archive: + meta = json.loads(str(archive["metadata"])) + if meta.get("format") != "struphy-geometry" or meta.get("version") != 1: + raise ValueError("Unsupported geometry archive format/version") + if "kind_map" not in meta: + if meta.get("type") not in _ANALYTIC: + raise ValueError("Unknown analytic mapping") + return getattr(domains, meta["type"])(**meta["params"]) + kind = meta["kind_map"] + if kind not in (0, 1, 2): + raise ValueError("Unknown spline mapping kind") + ndim = 3 if kind == 0 else 2 + kwargs = {key: tuple(meta[key][:ndim]) for key in ("num_elements", "degree", "spl_kind")} + kwargs.update(cx=archive["cx"].copy(), cy=archive["cy"].copy()) + if kind == 0: + return Spline(**kwargs, cz=archive["cz"].copy()) + # Domain extends poloidal control points with a singleton third axis. + kwargs["cx"] = kwargs["cx"].reshape(kwargs["cx"].shape[:2]) + kwargs["cy"] = kwargs["cy"].reshape(kwargs["cy"].shape[:2]) + param = float(archive["params_numpy"][0]) + if kind == 1: + return PoloidalSplineStraight(**kwargs, Lz=param) + return PoloidalSplineTorus(**kwargs, tor_period=param) diff --git a/src/struphy/geometry/utilities.py b/src/struphy/geometry/utilities.py index 625b384d2..50c052384 100644 --- a/src/struphy/geometry/utilities.py +++ b/src/struphy/geometry/utilities.py @@ -2,7 +2,7 @@ "Domain-related utility functions." import logging -from typing import Callable +from typing import Callable, Literal import cunumpy as xp import numpy as np @@ -17,7 +17,6 @@ # if TYPE_CHECKING: from struphy.geometry.base import Domain, PoloidalSplineTorus from struphy.geometry.utilities_kernels import weighted_arc_lengths_flux_surface -from struphy.io.options import LiteralOptions from struphy.linear_algebra.linalg_kron import kron_lusolve_2d logger = logging.getLogger("struphy") @@ -386,7 +385,7 @@ class TransformedPformComponent: def __init__( self, fun: Callable | list, - given_in_basis: LiteralOptions.GivenInBasis, + given_in_basis: Literal["0", "1", "2", "3", "v", "physical", "physical_at_eta", "norm", None], out_form: str, comp: int = 0, domain: Domain = None, diff --git a/src/struphy/utils/class_helpers.py b/src/struphy/utils/class_helpers.py new file mode 100644 index 000000000..bcfe4cd6a --- /dev/null +++ b/src/struphy/utils/class_helpers.py @@ -0,0 +1,24 @@ +"""Class introspection helpers without simulation or MPI dependencies.""" + +import inspect + + +def __class_with_params_repr_no_defaults__(cls_instance): + sig = inspect.signature(cls_instance.__class__.__init__) + defaults = {k: v.default for k, v in sig.parameters.items() if k != "self"} + out = f"{cls_instance.__class__.__name__}(" + for k, v in cls_instance.params.items(): + if k in defaults and v != defaults[k]: + out += f"{k}={v}, " + out += ")" + return out + + +def all_class_params_are_default(cls_instance): + return cls_instance.__repr_no_defaults__() == cls_instance.__class__.__name__ + "()" + + +def all_subclasses(cls): + subclasses = cls.__subclasses__() + subclasses = subclasses + [g for s in subclasses for g in all_subclasses(s)] + return subclasses diff --git a/src/struphy/utils/utils.py b/src/struphy/utils/utils.py index 1d7757043..9f5dbe429 100644 --- a/src/struphy/utils/utils.py +++ b/src/struphy/utils/utils.py @@ -13,6 +13,11 @@ from feectools.ddm.mpi import mpi as MPI import struphy +from struphy.utils.class_helpers import ( + __class_with_params_repr_no_defaults__, + all_class_params_are_default, + all_subclasses, +) logger = logging.getLogger("struphy") @@ -167,21 +172,6 @@ def __dataclass_repr_all_stacked__(obj): return out -def __class_with_params_repr_no_defaults__(cls_instance): - sig = inspect.signature(cls_instance.__class__.__init__) - defaults = {k: v.default for k, v in sig.parameters.items() if k != "self"} - out = f"{cls_instance.__class__.__name__}(" - for k, v in cls_instance.params.items(): - if k in defaults and v != defaults[k]: - out += f"{k}={v}, " - out += ")" - return out - - -def all_class_params_are_default(cls_instance): - return cls_instance.__repr_no_defaults__() == cls_instance.__class__.__name__ + "()" - - def ruff_autofix_and_format(code: str) -> str: with tempfile.NamedTemporaryFile(suffix=".py", mode="w+") as tmp: tmp.write(code) @@ -205,12 +195,6 @@ def ruff_autofix_and_format(code: str) -> str: return result -def all_subclasses(cls): - subclasses = cls.__subclasses__() - subclasses = subclasses + [g for s in subclasses for g in all_subclasses(s)] - return subclasses - - if __name__ == "__main__": from struphy import set_logging_level