diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index b61c142..d8c256f 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -8,6 +8,10 @@ name: encomp CI/release # plus source-only doc/static tests # docs-rtd -> rebuild the docs in a ReadTheDocs-equivalent environment (no project # install, no compiled plugin) so an RTD-only break fails the PR +# oracle-parity -> install Python CoolProp only as a test oracle and compare it with +# the bundled implementation +# free-threaded-smoke -> visible, non-publishing Python 3.14t readiness probe; allowed +# to fail until the released dependency stack supports it completely # dependency-floors -> install at the LOWEST versions the declared bounds allow and run the # suite, so a too-low bound fails here rather than on a user's machine # rust -> cargo fmt, clippy, and in-crate tests @@ -96,6 +100,24 @@ jobs: READTHEDOCS: "True" run: uv run --no-sync python -m sphinx -T -W --keep-going -b html docs docs/_build/html + oracle-parity: + name: bundled CoolProp parity oracle + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@9c091bb21b7c1c1d1991bb908d89e4e9dddfe3e0 # v7.0.0 + - uses: astral-sh/setup-uv@fac544c07dec837d0ccb6301d7b5580bf5edae39 # v8.2.0 + with: + python-version: "3.13" + - name: Install development tools and the test-only Python CoolProp oracle + run: uv sync --locked --no-install-project --group oracle + - name: Build the one bundled CoolProp implementation and native artifact + run: | + python -m pip install cmake + uv run --no-sync python scripts/build_libcoolprop.py + uv run --no-sync maturin develop --release + - name: Compare native scalar/plugin results with Python CoolProp 8.0.0 + run: uv run --no-sync pytest encomp/tests/test_coolprop_plugin.py encomp/tests/test_rust_parity.py -q + dependency-floors: # The `test` job below installs the wheel's dependencies from PyPI at their LATEST # versions, so it can never notice that a declared lower bound is too low. Resolve the @@ -157,12 +179,16 @@ jobs: import importlib.metadata as metadata from pathlib import Path - names = ("encomp", "pint", "pydantic", "typeguard", "typing_extensions", "numpy", "polars", "sympy", "uncertainties", "coolprop") + names = ("encomp", "pint", "pydantic", "typeguard", "typing_extensions", "numpy", "polars", "sympy", "uncertainties") for name in names: print(f"{name:20s} {metadata.version(name)}") + import importlib.util import encomp.coolprop as cp + assert importlib.util.find_spec("CoolProp") is None, "Python CoolProp must not be installed" + assert not any(r.lower().startswith("coolprop") for r in (metadata.requires("encomp") or [])) + package = Path(cp.__file__).parent print(f"\nencomp.coolprop imported from {package}") for path in sorted(package.iterdir()): @@ -195,6 +221,44 @@ jobs: - run: cargo clippy --manifest-path rust/Cargo.toml --all-targets -- -D warnings - run: PYO3_PYTHON="$(which python)" cargo test --manifest-path rust/Cargo.toml --no-default-features + free-threaded-smoke: + name: Python 3.14t smoke (experimental) + runs-on: ubuntu-latest + # This preparatory signal must not block normal-GIL releases. Python 3.15t wheels + # remain gated on a released Polars/PyO3 stack and the full verification in PLAN.md. + continue-on-error: true + steps: + - uses: actions/checkout@9c091bb21b7c1c1d1991bb908d89e4e9dddfe3e0 # v7.0.0 + - uses: astral-sh/setup-uv@fac544c07dec837d0ccb6301d7b5580bf5edae39 # v8.2.0 + with: + python-version: "3.14t" + - name: Install and build without the Python CoolProp package + run: | + python -m pip install cmake + uv sync --locked --no-install-project + uv run --no-sync python scripts/build_libcoolprop.py + uv run --no-sync maturin develop --release + - name: Keep the GIL disabled while importing and exercising both native interfaces + env: + PYTHON_GIL: "0" + run: | + uv run --no-sync python -W error - <<'PY' + import importlib.util + import sys + + import polars as pl + + from encomp import coolprop as cp + from encomp.fluids import Water + from encomp.units import Quantity as Q + + assert not sys._is_gil_enabled() + assert importlib.util.find_spec("CoolProp") is None + assert Water(P=Q(5e6, "Pa"), T=Q(400.0, "K")).D.m > 900.0 + assert pl.DataFrame({"P": [5e6], "T": [400.0]}).select(cp.water("DMASS", "P", "T"))[0, 0] > 900.0 + assert not sys._is_gil_enabled() + PY + build-wheels: name: build ${{ matrix.os }} runs-on: ${{ matrix.os }} @@ -240,16 +304,26 @@ jobs: shell: bash run: | python -m pip install --upgrade pip - # install the PREBUILT encomp wheel (abi3 cp313 installs on 3.14 too) by - # name so it is not rebuilt from source; deps + test deps come from PyPI - version=$(python -c 'from pathlib import Path; wheels = sorted(Path("wheelhouse").glob("encomp-*.whl")); assert len(wheels) == 1, f"expected exactly one encomp wheel, found {len(wheels)}"; print(wheels[0].name.split("-")[1])') - pip install --find-links wheelhouse "encomp==$version" pytest hypothesis - python -c "from encomp import coolprop as cp; assert cp.self_check(), 'plugin self_check failed'" + # Install the downloaded artifact by its exact path. PR builds deliberately + # retain the current released version, so asking pip for `encomp==$version` + # can select that same-version wheel from PyPI instead of this candidate. + wheel=$(python -c 'from pathlib import Path; wheels = sorted(Path("wheelhouse").glob("encomp-*.whl")); assert len(wheels) == 1, f"expected exactly one encomp wheel, found {len(wheels)}"; print(wheels[0])') + pip install "$wheel" pytest hypothesis + python - <<'PY' + import importlib.metadata as metadata + import importlib.util + + from encomp import coolprop as cp + + assert cp.self_check(), "plugin self_check failed" + assert importlib.util.find_spec("CoolProp") is None, "Python CoolProp was unexpectedly installed" + assert not any(r.lower().startswith("coolprop") for r in (metadata.requires("encomp") or [])) + PY python -m pytest --pyargs encomp.tests -q publish: name: publish to PyPI - needs: [typecheck, docs-rtd, dependency-floors, rust, build-wheels, test] + needs: [typecheck, docs-rtd, oracle-parity, dependency-floors, rust, build-wheels, test] if: github.event_name == 'push' && startsWith(github.ref, 'refs/tags/v') runs-on: ubuntu-latest environment: pypi # configure the Trusted Publisher against this env diff --git a/AGENTS.md b/AGENTS.md index b99857a..5c2026f 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -73,9 +73,9 @@ assert pressure.to("kPa").m == 100.0 comparison and pickling overrides on `Quantity` intentionally narrow pint's signatures — that's the point of the library, not an LSP bug to fix. - CoolProp evaluation is **one property per DAG node** by design; don't batch multiple - output properties into one plugin call. All CoolProp paths hold the GIL — the Rust - plugin (`rust/`) exists precisely to evaluate without it; don't add Python-side - parallelism around CoolProp. + output properties into one plugin call. Scalars detach from the GIL in the direct + PyO3 bridge and own thread-local native states; arrays run through the Rust/Polars + plugin. Don't add Python-side parallelism around CoolProp. ## Tooling diff --git a/README.md b/README.md index 9a13626..ac193fb 100644 --- a/README.md +++ b/README.md @@ -44,7 +44,9 @@ The remaining modules (`encomp.gases`, `encomp.conversion`, `encomp.constants`, pip install encomp ``` -`encomp` ships as a single per-platform wheel that bundles the compiled Rust plugin and the CoolProp shared library, so supported platforms have nothing to build. Wheels are provided for Windows (x86_64), Linux (x86_64 and arm64), and macOS (Apple Silicon only). PyPI does not publish an sdist; unsupported platforms, including Intel Macs, need a build from the git repository; see [Tests](#tests). +`encomp` ships as a single per-platform wheel containing its dual-purpose native Rust artifact and one bundled CoolProp shared library, so supported platforms have nothing to build and do not install the Python `CoolProp` package. Wheels are provided for Windows (x86_64), Linux (x86_64 and arm64), and macOS (Apple Silicon only). PyPI does not publish an sdist; unsupported platforms, including Intel Macs, need a build from the git repository; see [Tests](#tests). + +On the macOS arm64 wheel measured for the native-only change, the resolver's combined wheel download falls from 19.61 MB (`encomp` plus the separate Python CoolProp wheel) to 9.08 MB for `encomp` alone, a 53.7% reduction. The `encomp` artifact itself grows by about 35 KB for the scalar bridge and bundled license; removing the separate 10.57 MB wheel is where the installation-size saving comes from. ## The `Quantity` class @@ -378,7 +380,7 @@ w: Water[pl.Expr] = Water(P=Q(pl.col("P"), "Pa"), T=Q(pl.col("T"), "K")) df.select(w.D.m.alias("rho"), w.H.m.alias("h"), w.S.m.alias("s")) ``` -`pl.Expr` (lazy) inputs are evaluated exclusively through the plugin (there is no `map_batches` fallback). Eager `float` / NumPy / `pl.Series` inputs use the Python CoolProp path, except arrays of at least `EAGER_PLUGIN_MIN_SIZE` (1000) elements, which also route through the plugin. The two paths are verified to agree on value, `NaN`/null handling, and dtype; results are bit-identical when the installed `coolprop` matches the bundled build (8.0.0). +`encomp` contains one thermodynamic implementation: the bundled CoolProp shared library. Scalars call it through a direct PyO3 bridge with a bounded per-thread state cache; NumPy arrays, `pl.Series`, and `pl.Expr` inputs use the batched Polars plugin path at every size. Both private interfaces live in the same `_internal` artifact and share the same loaded library. The Python `CoolProp` package is not installed or imported. `EAGER_PLUGIN_MIN_SIZE` remains temporarily as a deprecated no-op for compatibility. The plugin is also usable directly on any Polars expression, independent of the `Fluid` class (the `encomp.coolprop` package): @@ -402,15 +404,9 @@ df.select( ### Implementation -The GIL is not the only serialization point: the CoolProp C-API takes a global handle-table lock on every call, so per-row calls serialize even in Rust. The plugin uses the batched C-API (`AbstractState_update_and_1_out`): one call per chunk, the handle lock taken once at construction, then the flash loop runs lock-free in C++. Independent chunks and independent property expressions parallelize. - -Benchmarks (CoolProp 8.0, 14-thread pool): +Scalar calls reuse a bounded thread-local cache of native `AbstractState` handles; a handle is never shared between threads. Arrays use the batched C API (`AbstractState_update_and_1_out`), with handle-table locking confined to construction and destruction. The flash loop itself runs lock-free in C++, so independent chunks and property expressions can execute in parallel. -| workload | vs `map_batches` | notes | -| --- | --- | --- | -| single property `D`, 1M rows | ~2.1x | also ~2x faster than vectorized `PropsSI` | -| 4 independent properties, one `collect()`, 1M rows | ~4.6x | `map_batches` is serial on the GIL; the plugin runs ~4 cores | -| 8 enthalpy evaluations, 1M rows | ~4.9x | ~6 cores vs 1, roughly half the peak memory | +On the normal-GIL macOS benchmark used for this change, public scalar calls improved by 4-33%. The existing 100k-row HEOS array path changed from 470.49 ms to 469.92 ms, while three independent properties changed from 534.06 ms to 537.51 ms (+0.6%). See `encomp/coolprop/README.md` and `scripts/benchmark_coolprop.py` for the full reproducible measurements. Each `fluid(...)` / `humid_air(...)` is an independent plugin node, so selecting *K* properties of one state runs *K* flashes of it — Polars cannot reuse the shared flash across opaque plugin nodes. Independent properties still parallelize, so this is total work, not wall-clock. See `encomp/coolprop/README.md` for the design, thread-safety model, and caveats. @@ -439,7 +435,7 @@ For array magnitudes, convert the expression to a NumPy-aware function with `enc ## Tests -Development checkouts build the native CoolProp plugin locally. Install Rust, CMake, git, and a C++ compiler, then from the repository root run: +Development checkouts build the dual-purpose native CoolProp artifact locally. Install Rust, CMake, git, and a C++ compiler, then from the repository root run: ```bash python scripts/build_libcoolprop.py @@ -447,7 +443,7 @@ uv sync --all-extras --all-groups uv run pytest ``` -See `encomp/coolprop/README.md` for the plugin build details. +See `encomp/coolprop/README.md` for native build details. The test suite ships inside the wheel, so an installed `encomp` doubles as its own post-install smoke test (it needs only `pytest` and `hypothesis`): diff --git a/docs/usage.md b/docs/usage.md index fd86874..cf290f1 100644 --- a/docs/usage.md +++ b/docs/usage.md @@ -716,7 +716,7 @@ density = Fluid("HEOS::CO2[0.7]&O2[0.3]", P=Q(10, "bar"), T=Q(300, "K")).assume_ ### Using vector inputs -CoolProp evaluates vector inputs in a single backend call. +CoolProp evaluates vector inputs through encomp's batched native plugin path. The inputs are {py:class}`encomp.units.Quantity` instances with vector magnitudes: one-dimensional NumPy arrays or `pl.Series`, all of the same length (or a single scalar, which is repeated). ```python @@ -757,7 +757,7 @@ Missing or out-of-range results surface as `NaN` (for a numpy magnitude) or `nul ### Parallel evaluation with Polars {py:class}`encomp.fluids.Fluid` properties also accept `Quantity`-wrapped Polars expressions (`pl.Expr`) and return a `pl.Expr`. Independent property nodes in one `select` / `with_columns` / `collect()` (eager or lazy) are evaluated in parallel by the `encomp.coolprop` plugin -- a native Rust extension over the CoolProp C-API that runs without holding the GIL. -`pl.Expr` (lazy) inputs are evaluated exclusively through this plugin (there is no `map_batches` fallback). Eager `float` / NumPy / `pl.Series` inputs use the Python CoolProp path, except vector magnitudes of at least `EAGER_PLUGIN_MIN_SIZE` (1000) elements, which also route through the plugin (results are bit-identical when the installed `coolprop` matches the bundled build, 8.0.0). +`encomp` bundles CoolProp once. Scalar inputs use a direct PyO3 bridge with per-thread cached native states; NumPy arrays, `pl.Series`, and `pl.Expr` inputs use the batched plugin at every size (there is no Python `CoolProp` or `map_batches` fallback). `EAGER_PLUGIN_MIN_SIZE` is retained only as a deprecated no-op. ```python import polars as pl diff --git a/encomp/coolprop/LICENSE.CoolProp b/encomp/coolprop/LICENSE.CoolProp new file mode 100644 index 0000000..2bcd117 --- /dev/null +++ b/encomp/coolprop/LICENSE.CoolProp @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2012-2018 Ian H. Bell and other CoolProp developers + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/encomp/coolprop/README.md b/encomp/coolprop/README.md index 2981d10..939e6c4 100644 --- a/encomp/coolprop/README.md +++ b/encomp/coolprop/README.md @@ -1,11 +1,11 @@ # encomp.coolprop -Parallel CoolProp property evaluation as **native Polars expression plugins** -(Rust, `pyo3-polars`). Independent property nodes in one `collect()` run in -parallel on the Polars thread pool — GIL-free — instead of a `map_batches` Python -UDF that holds the GIL and serializes. Usable directly, or as the evaluation path -of `encomp.fluids` for `pl.Expr` inputs and large eager arrays (≥1000 elements); -small eager inputs (scalars, short arrays) use the Python CoolProp path. +One bundled CoolProp implementation behind a dual-purpose native Rust artifact. +Its private PyO3 interface handles scalar evaluation, metadata, parsing, and +validation; its Polars plugin ABI handles every array and expression size. +Independent property nodes in one `collect()` run in parallel on the Polars +thread pool without holding the GIL. The Python `CoolProp` package is neither a +runtime dependency nor imported by `encomp`. ## Usage @@ -39,43 +39,54 @@ The API mirrors `encomp.fluids.Fluid`: - `water(output, in1, in2)` — the IF97 water/steam shorthand, mirroring `encomp.fluids.Water`. Equivalent to `fluid(..., name="IF97::Water")`. It takes no `composition` (water is pure) and no `assume_phase` (IF97 would ignore it). -- `humid_air(output, in1, in2, in3)` — same naming rule for the three HAPropsSI inputs. +- `humid_air(output, in1, in2, in3)` — same naming rule for the three humid-air inputs. - `FluidInput` / `HumidAirInput` (state inputs), `FluidParam` / `HumidAirParam` / `Backend` / `Phase` / `AssumedPhase` are `Literal`s; `CName` / `Composition` mirror `encomp.fluids`. The matching `frozenset`s and `is_*` `TypeIs` predicates are exported too, along with `resolve_fluid_spec` (name → backend/fluids/fractions). -## Performance (CoolProp 8.0, 14-thread pool) +## Performance regression gate -``` -CORRECTNESS plugin vs raw PropsSI (IF97 water): 0.0e+00 exact (D/H/S/Cp) +Apple M4 Pro, CPython 3.14.3, Polars 1.42.1, CoolProp 8.0.0; medians of seven +warmed runs. Lower is better. The baseline is the unmodified `v1.8.0` wheel. -SPEED (property D) raw PropsSI | encomp map_batches | rust plugin | vs m_b - N=1 0.003 ms | 0.114 ms | 0.029 ms | 3.98x - N=1M 177.6 ms | 184.1 ms | 86.6 ms | 2.12x +| workload | v1.8.0 | native-only | change | +| --- | ---: | ---: | ---: | +| public IF97 scalar | 62.49 us | 59.68 us | -4.5% | +| public HEOS scalar | 97.36 us | 65.14 us | -33.1% | +| public humid-air scalar | 57.35 us | 55.30 us | -3.6% | +| HEOS array, 100k rows | 470.49 ms | 469.92 ms | -0.1% | +| three independent HEOS properties, 100k rows | 534.06 ms | 537.51 ms | +0.6% | -PARALLELISM 4 independent props, ONE collect(), 1M rows - encomp map_batches : combined 823 ms, 1.00x serial (GIL) - rust plugin : combined 178 ms, 2.57x PARALLEL (4.6x faster than encomp) +The direct bridge itself takes 0.20 us for a warm IF97 state, 4.88 us for a +varying-temperature cached HEOS state, and 2.05 us for humid air. The matching +high-level bundled C calls take 1.06 us, 31.05 us, and 2.33 us respectively; +the state-owning bridge is faster where it can reuse `AbstractState`. A cold +first IF97 call is 28.54 us, and a fill-plus-eviction pass across 17 HEOS +configurations averages 29.37 us/configuration. Reproduce with +`scripts/benchmark_coolprop.py`. -8 enthalpy calcs, 1M rows: 4.9x faster, ~6 cores vs 1, ~half the peak memory. -``` +The matching macOS arm64 package check reduces the combined wheel download from +19.61 MB (`encomp` plus Python CoolProp) to 9.08 MB for `encomp` alone (-53.7%). +The encomp wheel itself is about 35 KB larger because the old Python binding was +a separate distribution; removing that 10.57 MB wheel produces the net saving. ## Design -The GIL is not the only serialization point: CoolProp's C-API takes a global -handle-table lock on every call, so per-row calls serialize even in Rust. The -plugin uses the batched C-API (`AbstractState_update_and_1_out`): one -call per chunk, the handle lock taken once at construction, then the flash loop -runs lock-free in C++ — so independent chunks/expressions parallelize. +The direct scalar interface keeps a bounded LRU of `AbstractState` handles in +thread-local storage. A handle is mutated only by its owning thread and is freed +on eviction or thread shutdown. The plugin uses the batched C-API +(`AbstractState_update_and_1_out`): one call per chunk, with the handle-table lock +taken only for construction/destruction and the flash loop lock-free in C++. Thread-safety (is it safe to evaluate concurrently?): the math is pure — a flash on a per-handle `AbstractState` mutates no global state — safe to run concurrently **iff** (a) each thread owns its handle (never shared; a state caches its last flash), (b) handle create/destroy is synchronized (the global handle table is the one shared structure — `coolprop.rs` guards `factory`/`free` with a narrow mutex, -never the hot path), and (c) global config isn't mutated during evaluation. Results -are bit-identical across thread-pool sizes (1/4/8 are covered by `test_thread_count_parity`). +never the hot path), and (c) global config isn't mutated during evaluation. The +PyO3 module stores no Python objects in global state and explicitly declares that +it does not require the GIL. All `unsafe` is confined to `rust/src/coolprop.rs` (the FFI boundary); `lib.rs` has none. Every `unsafe` block carries a `// SAFETY:` comment (rationale + error modes), @@ -84,9 +95,8 @@ enforced by `clippy::undocumented_unsafe_blocks`. The safe wrapper uses an RAII buffers. libCoolProp is loaded at runtime via `libloading`, so the same plugin works on macOS/Linux/Windows by shipping that platform's `.dylib`/`.so`/`.dll`. -HumidAir caveat: `humid_air` is correct and GIL-free, but `HAPropsSI` is internally -more serialized than the Fluid path — it parallelizes only ~1.25x (vs ~2.5x) and is -much slower per call (iterative solver). Still better than the Python path. +Humid-air caveat: `humid_air` is correct and GIL-free, but CoolProp's humid-air +solver is internally iterative and more serialized than the fluid path. ## Caveats @@ -100,12 +110,9 @@ much slower per call (iterative solver). Still better than the Python path. (bumping `pyo3-polars` to the new Rust polars) is needed only if polars bumps the ffi ABI **major** — rare and announced. `self_check()` evaluates a known value at first use, so a genuinely incompatible polars surfaces as a clear error, never a silent wrong result. -- **Version match**: CoolProp enum integers can differ across major versions. `encomp` - requires Python `coolprop>=8.0.0,<9`, while the bundled Rust lib is built at 8.0.0. - The Rust side resolves the *output* parameter index via CoolProp at runtime, never - hardcoded. The one integer that crosses from Python into the bundled lib is the - `input_pairs` enum index, computed by the Python `coolprop` (`generate_update_pair`); - that enum is stable within a CoolProp major, which is why the requirement caps it. +- **One version source**: `encomp/coolprop/_build_info.py` pins the bundled version and + upstream commit. Parameter indexes and input-pair enum values are resolved at runtime + from that bundled library; no numeric CoolProp enum is hardcoded in Python or Rust. - **One property per node (no output batching).** Each `fluid(...)` / `humid_air(...)` is an independent plugin node, so selecting K properties of one state runs K flashes of it — Polars cannot reuse (CSE) the shared flash across opaque plugin nodes. Independent @@ -113,9 +120,10 @@ much slower per call (iterative solver). Still better than the Python path. ## Build / install -This is part of `encomp`: the whole package (Python + this compiled plugin + -bundled `libCoolProp`) ships in ONE per-platform wheel. `libCoolProp` is built from -source once and bundled, so every build starts with `scripts/build_libcoolprop.py`. +This is part of `encomp`: the whole package (Python + the dual-purpose `_internal` +artifact + bundled `libCoolProp`) ships in one per-platform wheel. `libCoolProp` is +built from source once and bundled, so every build starts with +`scripts/build_libcoolprop.py`. Dev (from the repo root, editable): @@ -148,6 +156,6 @@ PYO3_PYTHON=$(pwd)/.venv/bin/python cargo test --manifest-path rust/Cargo.toml - ## Files - `encomp/coolprop/__init__.py` — public Python API (`fluid`, `water`, `humid_air`). -- `rust/src/lib.rs` — the `cp_evaluate` / `ha_evaluate` plugin expressions (crate at repo root). +- `rust/src/lib.rs` — the PyO3 scalar API plus `cp_evaluate` / `ha_evaluate` plugin expressions. - `rust/src/coolprop.rs` — the CoolProp C-API bindings + thread-safety model (all `unsafe`). - `scripts/build_libcoolprop.py` — builds + bundles CoolProp's shared library. diff --git a/encomp/coolprop/__init__.py b/encomp/coolprop/__init__.py index bed1646..5c8853e 100644 --- a/encomp/coolprop/__init__.py +++ b/encomp/coolprop/__init__.py @@ -25,6 +25,7 @@ from __future__ import annotations +import importlib import importlib.util import logging import math @@ -32,19 +33,48 @@ import warnings from functools import cache, lru_cache from pathlib import Path -from typing import Any, Literal, TypeIs, cast, get_args +from typing import Literal, Protocol, TypeIs, cast, get_args -import CoolProp.CoolProp as _CoolProp import polars as pl from polars.plugins import register_plugin_function from .._polars_dtype import canonical_unit_string as _canonical_unit_string -# CoolProp ships incomplete type stubs; treat the module as Any (mirrors encomp). -_cp: Any = _CoolProp - _LOGGER = logging.getLogger(__name__) + +class _NativeModule(Protocol): + """Typed view of the private PyO3 API in ``_internal``.""" + + def initialize(self, lib_path: str) -> None: ... + def parameter_index(self, name: str) -> int: ... + def parameter_information(self, name: str, field: str) -> str: ... + def resolve_input_pair(self, name1: str, name2: str) -> tuple[int, bool]: ... + def resolve_fluid_name(self, name: str) -> tuple[str, str, list[float] | None]: ... + def validate_fluid(self, backend: str, fluid: str, fractions: list[float] | None = None) -> None: ... + def prepare_fluid( + self, + backend: str, + fluid: str, + fractions: list[float] | None = None, + phase: str | None = None, + ) -> int: ... + def prepare_humid_air(self, output: str, name1: str, name2: str, name3: str) -> int: ... + def fluid_scalar( + self, + config_id: int, + input_pair: int, + value1: float, + value2: float, + output: int, + ) -> float: ... + def humid_air_scalar(self, config_id: int, value1: float, value2: float, value3: float) -> float: ... + def lib_version(self) -> str: ... + def clear_scalar_cache(self) -> None: ... + def scalar_cache_info(self) -> tuple[int, int, int, int, int]: ... + def handle_counts(self) -> tuple[int, int]: ... + + # mole fractions may sum this far from 1 (float rounding) before it is an error COMPOSITION_SUM_TOLERANCE = 0.01 @@ -453,6 +483,20 @@ def lib_path() -> str: raise RuntimeError(f"bundled CoolProp library ({', '.join(_LIB_NAMES)}) not found in {_HERE}") +@lru_cache(maxsize=1) +def _native() -> _NativeModule: + """Import and initialize the private PyO3 side of the dual-purpose artifact. + + The import stays lazy so ReadTheDocs can import the source tree without a locally + built native artifact. Released wheels always contain both this module and the + bundled library. + """ + + module = cast(_NativeModule, importlib.import_module("encomp.coolprop._internal")) + module.initialize(lib_path()) + return module + + @lru_cache(maxsize=1) def _plugin_path() -> str: """Absolute path to the compiled plugin extension (``_internal.*``). @@ -479,12 +523,14 @@ def _plugin_path() -> str: @cache def _resolve_pair(name1: str, name2: str) -> tuple[int, bool]: - # CoolProp's generate_update_pair gives the canonical input_pairs index and - # value order; the swap decision is value-independent, so resolve it once. - pair, a, _ = _cp.generate_update_pair(_cp.get_parameter_index(name1), 1.0, _cp.get_parameter_index(name2), 2.0) - if int(pair) == 0: - raise ValueError(f"unsupported CoolProp input pair: {name1!r}, {name2!r}") - return int(pair), (a == 2.0) + # Rust normalizes aliases with the bundled library's parameter indexes, matches + # a reviewed pair table, and obtains the pair enum from that same library. + return _native().resolve_input_pair(name1, name2) + + +@cache +def _parameter_index(name: str) -> int: + return _native().parameter_index(name) def _as_expr(x: str | pl.Expr) -> pl.Expr: @@ -513,7 +559,7 @@ def _expected_si_unit(name: str) -> str: metadata, so the plugin can accept a unit-typed input by comparing the two strings for equality. """ - return _canonical_unit_string(_cp.get_parameter_information(_cp.get_parameter_index(name), "units")) + return _canonical_unit_string(_native().parameter_information(name, "units")) @cache @@ -610,14 +656,10 @@ def resolve_fluid_spec( fractions: list[float] | None = mole_fractions is_concentration = False else: - backend_raw, fluid_str = _cp.extract_backend(name) - backend = "HEOS" if backend_raw == "?" else backend_raw - species, fracs = _cp.extract_fractions(fluid_str) - fluids = "&".join(species) - # extract_fractions validates + returns fractions only when the name carries them; a - # single species with a fraction is an incompressible concentration (absolute basis). - fractions = [float(x) for x in fracs] if fracs else None - is_concentration = fractions is not None and len(species) == 1 + backend, fluids, fractions = _native().resolve_fluid_name(name) + # The native parser validates + returns fractions only when the name carries + # them; one remaining species is an incompressible concentration (absolute basis). + is_concentration = fractions is not None and "&" not in fluids # mixture mole fractions must sum to 1; a concentration is absolute and passes through if fractions is not None and not is_concentration: total = sum(fractions) @@ -630,6 +672,58 @@ def resolve_fluid_spec( return backend, fluids, fractions +@cache +def _prepared_fluid_config( + backend: str, + fluids: str, + fractions: tuple[float, ...] | None, + phase: Phase | None, +) -> int: + return _native().prepare_fluid(backend, fluids, None if fractions is None else list(fractions), phase) + + +@cache +def _prepared_humid_air_config(output: str, name1: str, name2: str, name3: str) -> int: + return _native().prepare_humid_air(output, name1, name2, name3) + + +def _fluid_scalar( # pyright: ignore[reportUnusedFunction] - private sibling-module bridge + output: str, + name1: str, + value1: float, + name2: str, + value2: float, + *, + name: CName, + assume_phase: AssumedPhase | None = None, + composition: Composition | None = None, +) -> float: + """Private direct native scalar adapter used by :mod:`encomp.fluids`.""" + + pair, swap = _resolve_pair(name1, name2) + a, b = (value2, value1) if swap else (value1, value2) + backend, fluids, fractions = resolve_fluid_spec(name, composition) + phase = _phase_from_assumed(assume_phase) if assume_phase is not None else None + fraction_key = None if fractions is None else tuple(fractions) + config_id = _prepared_fluid_config(backend, fluids, fraction_key, phase) + return _native().fluid_scalar(config_id, pair, a, b, _parameter_index(output)) + + +def _humid_air_scalar( # pyright: ignore[reportUnusedFunction] - private sibling-module bridge + output: str, + name1: str, + value1: float, + name2: str, + value2: float, + name3: str, + value3: float, +) -> float: + """Private direct native humid-air scalar adapter used by :mod:`encomp.fluids`.""" + + config_id = _prepared_humid_air_config(output, name1, name2, name3) + return _native().humid_air_scalar(config_id, value1, value2, value3) + + def fluid( output: FluidParam, input1: FluidInput | pl.Expr, @@ -678,7 +772,7 @@ def fluid( ) _reject_duplicate_input_keys( "fluid", - ((name1, _cp.get_parameter_index(name1)), (name2, _cp.get_parameter_index(name2))), + ((name1, _parameter_index(name1)), (name2, _parameter_index(name2))), ) pair_idx, swap = _resolve_pair(name1, name2) a, b = (input2, input1) if swap else (input1, input2) # canonical order @@ -795,12 +889,7 @@ def humid_air( @lru_cache(maxsize=1) def lib_version() -> str: """CoolProp version of the bundled library (cached; the bundled lib never changes).""" - import ctypes - - lib = ctypes.CDLL(lib_path()) - buf = ctypes.create_string_buffer(256) - lib.get_global_param_string(b"version", buf, 256) - return buf.value.decode() + return _native().lib_version() @lru_cache(maxsize=1) diff --git a/encomp/coolprop/_build_info.py b/encomp/coolprop/_build_info.py new file mode 100644 index 0000000..892d142 --- /dev/null +++ b/encomp/coolprop/_build_info.py @@ -0,0 +1,5 @@ +"""Committed identity of the one CoolProp implementation bundled in encomp wheels.""" + +BUNDLED_COOLPROP_VERSION = "8.0.0" +BUNDLED_COOLPROP_TAG = f"v{BUNDLED_COOLPROP_VERSION}" +BUNDLED_COOLPROP_COMMIT = "ae81610e7d23efc57f9d051c8e70a4d66e87537f" diff --git a/encomp/fluids.py b/encomp/fluids.py index 1f6f9de..4728ed5 100644 --- a/encomp/fluids.py +++ b/encomp/fluids.py @@ -1,39 +1,31 @@ -""" -Classes and functions relating to fluid properties. -Uses CoolProp as backend. -""" +"""Fluid-property classes backed by encomp's bundled CoolProp library.""" import hashlib import logging from abc import ABC, abstractmethod from collections import OrderedDict -from collections.abc import Callable, Mapping +from collections.abc import Mapping from functools import cache from threading import Lock from typing import Annotated, Any, ClassVar, Generic, Literal, Self, TypedDict, Unpack, cast -# CoolProp.CoolProp is a compiled extension module exporting both PropsSI and -# HAPropsSI. Importing the module (rather than the untyped pure-Python -# CoolProp.HumidAirProp) avoids a missing-stub warning. The functions are -# untyped, so they are exposed as typed aliases below. -import CoolProp.CoolProp as _CoolProp import numpy as np import polars as pl from .coolprop import ( + ASSUMED_PHASES, FLUID_INPUTS, HUMID_AIR_INPUTS, AssumedPhase, - Backend, CName, Composition, FluidParam, HumidAirParam, - Phase, - is_backend, + _fluid_scalar, # pyright: ignore[reportPrivateUsage] + _humid_air_scalar, # pyright: ignore[reportPrivateUsage] + _native, # pyright: ignore[reportPrivateUsage] is_fluid_param, is_humid_air_param, - is_phase, resolve_fluid_spec, ) from .settings import SETTINGS @@ -86,10 +78,6 @@ "clear_expr_evaluation_cache", ] -_cp: Any = _CoolProp -PropsSI: Callable[..., float | Numpy1DArray] = _cp.PropsSI -HAPropsSI: Callable[..., float | Numpy1DArray] = _cp.HAPropsSI - # Strict Literals for CoolProp property names, single source of truth in the # encomp.coolprop plugin (fluid + humid-air namespaces). The fluid-name / composition / # assumed-phase types (CName, CommonFluidName, Composition, FractionValue, AssumedPhase) @@ -109,6 +97,7 @@ "Variable", "N/A", ] +BackendKind = Literal["fluid", "humid_air"] class FluidState(TypedDict, Generic[MT], total=False): # noqa: UP046 @@ -187,29 +176,16 @@ class HumidAirState(TypedDict, Generic[MT], total=False): # noqa: UP046 psi_w: Quantity[Any, MT] | Quantity[Any, float] -# user-facing phase name (AssumedPhase, imported from encomp.coolprop) -> CoolProp phase -# index for the low-level AbstractState.specify_phase path. (The rust plugin path instead -# uses the ``phase_*`` string, resolved by encomp.coolprop._phase_from_assumed.) -_ASSUMED_PHASE_MAP: dict[str, Any] = { - "gas": _cp.iphase_gas, - "liquid": _cp.iphase_liquid, - "supercritical": _cp.iphase_supercritical, - "supercritical_gas": _cp.iphase_supercritical_gas, - "supercritical_liquid": _cp.iphase_supercritical_liquid, - "twophase": _cp.iphase_twophase, -} - # backends that ignore AbstractState.specify_phase (region-explicit), so assuming a -# phase has no effect and should not switch evaluation to the slower low-level loop +# phase has no effect _PHASE_IGNORING_BACKENDS = frozenset({"IF97"}) EAGER_PLUGIN_MIN_SIZE = 1000 -"""Eager numpy / ``pl.Series`` inputs with at least this many elements are evaluated -through the GIL-free rust plugin (faster and lower peak memory than CoolProp's -vectorized ``PropsSI``, and much faster than the low-level per-row Python loop for an -assumed phase / composition). Scalars and smaller arrays stay on the Python CoolProp -path, where the plugin's fixed dispatch overhead would dominate. The two paths are -verified to agree on dtype, value, and NaN/inf handling (``test_fluids``).""" +"""Deprecated compatibility constant; now a no-op. + +All eager arrays use the native Rust/Polars batch path, regardless of size. Scalars use +the direct PyO3 bridge, so there is no longer a Python CoolProp cutoff to configure. +""" # Caches the CONSTRUCTED pl.Expr per (fluid config, output, input-expr digests), so # repeated evaluations of the same property return the IDENTICAL expression object. @@ -239,7 +215,7 @@ def clear_expr_evaluation_cache() -> None: @cache -def _resolve_fluid_name(backend: str, fluids: str) -> None: +def _resolve_fluid_name(backend: str, fluids: str, fractions: tuple[float, ...] | None) -> None: # Constructing the AbstractState is the only reliable way to ask CoolProp whether a # name resolves (it covers pure fluids, INCOMP, and mixtures alike). It costs tens of # microseconds for the usual backends -- and seconds the very first time a tabular @@ -250,14 +226,14 @@ def _resolve_fluid_name(backend: str, fluids: str) -> None: # functools.cache stores return values, never exceptions, so only a resolved name is # remembered: a failure (an invalid name, or a transient CoolProp error) is re-checked # on the next call instead of being cached as "this fluid does not exist". - _cp.AbstractState(backend, fluids) + _native().validate_fluid(backend, fluids, None if fractions is None else list(fractions)) def _validate_fluid_name(name: CName, composition: Composition | None = None) -> None: - backend, fluids, _ = resolve_fluid_spec(name, composition) + backend, fluids, fractions = resolve_fluid_spec(name, composition) try: - _resolve_fluid_name(backend, fluids) + _resolve_fluid_name(backend, fluids, None if fractions is None else tuple(fractions)) except Exception as e: raise ValueError( f"Fluid '{name}' could not be initialized, ensure that the name is a valid CoolProp fluid name" @@ -287,9 +263,7 @@ def _get_expr_evaluation_cache_key( fluid_any._assumed_phase, comp_key, output, - bool(fluid_any._append_name_to_cp_inputs), - bool(fluid_any._evaluate_invalid_separately), - id(type(fluid).BACKEND["backend"]), + type(fluid).BACKEND_KIND, tuple((prop, _expr_cache_digest(expr)) for prop, expr in points), ) @@ -315,32 +289,16 @@ class CoolPropFluid(ABC, Generic[MT]): # noqa: UP046 name: CName points: list[tuple[CProperty, Quantity[Any, MT] | Quantity[Any, float]]] - BACKEND: ClassVar[dict[Literal["backend"], Callable[..., float | Numpy1DArray]]] = {"backend": PropsSI} - - # PropsSI expects the fluid name as the first input, but not HAPropsSI - _append_name_to_cp_inputs: bool = True - - # HAPropsSI fails if one or more inputs are incorrect, - # PropsSI returns NaN for invalid inputs in case valid inputs are also present - _evaluate_invalid_separately: bool = False + BACKEND_KIND: ClassVar[BackendKind] = "fluid" - # assumed phase (set via Fluid.assume_phase). When not None, property - # evaluation switches from the high-level PropsSI backend to the low-level - # AbstractState API with specify_phase, which skips the (for mixtures very - # expensive) phase-stability search. None means CoolProp determines the phase. + # Assumed phase passed to both native interfaces. None means CoolProp determines + # the phase; a value skips the expensive mixture phase-stability search. _assumed_phase: str | None = None - # fixed mixture composition (set via Fluid(..., composition=...)) as MOLE - # fractions, one float per species. When not None, evaluation uses the - # low-level AbstractState API (PropsSI cannot take an explicit composition). - # None means the composition is fixed in the fluid name, or the fluid is pure. + # Fixed mixture composition as mole fractions, one float per species. None means + # the composition is fixed in the fluid name, or the fluid is pure. _composition: Composition | None = None - # True when an assumed phase or a composition is configured, i.e. evaluation - # must use the low-level AbstractState path. A single flag so the normal - # PropsSI path pays one attribute read, not several is-not-None checks. - _lowlevel: bool = False - # display only head of vector inputs in the __repr__ method _repr_cutoff: int = 3 @@ -373,7 +331,7 @@ class CoolPropFluid(ABC, Generic[MT]): # noqa: UP046 8.0: "Not imposed", } - # unit and description for properties in function PropsSI + # CoolProp fluid-property names, units, and descriptions # (name1, name2, ...): (unit, description) # names are case-sensitive PROPERTY_MAP: dict[tuple[CProperty, ...], tuple[UnitString, str]] = { @@ -571,18 +529,8 @@ def __init__(self, name: CName, **kwargs: Quantity[Any, MT] | Quantity[Any, floa the ``INCOMP::`` prefix, optionally with a concentration (``"INCOMP::MEG[0.5]"``). An unknown name raises ``ValueError`` at construction. - The names recognized by the installed CoolProp build are listed by - - .. code:: python - - import CoolProp.CoolProp as CP - - CP.get_global_param_string("FluidsList") # pure and pseudo-pure - CP.get_global_param_string("incompressible_list_pure") - CP.get_global_param_string("incompressible_list_solution") - CP.get_global_param_string("predefined_mixtures") - - Examples: ``Water``, ``Air``, ``Nitrogen``, ``CarbonDioxide``, ``Ammonia``, + Examples recognized by the bundled build include ``Water``, ``Air``, + ``Nitrogen``, ``CarbonDioxide``, ``Ammonia``, ``R134a``, ``Toluene``, ``INCOMP::T66``, ``INCOMP::MPG[0.5]``, ``R410A.mix``. Refer to the CoolProp documentation for more information: @@ -770,156 +718,27 @@ def check_exception(self, prop: CProperty, e: ValueError) -> None: self._warn_coolprop_nan(prop, msg) - # ------------------------------------------------------------------ # - # Low-level AbstractState path. Used whenever an assumed phase and/or a - # ``composition=`` dict is set (PropsSI cannot do either). The composition is - # fixed (baked into the name, or passed as floats), so the AbstractState is - # built once and the loop only flashes each (P, T) row. - # ------------------------------------------------------------------ # - - def _constant_state(self) -> Any: # noqa: ANN401 - CoolProp AbstractState is untyped - # name/composition parsing is delegated to resolve_fluid_spec (the single source); - # here we just build the AbstractState and pin the assumed phase. - backend, fluids, fractions = resolve_fluid_spec(self.name, self._composition) - state = _cp.AbstractState(backend, fluids) - if fractions is not None: - state.set_mole_fractions(fractions) - if self._assumed_phase is not None: - state.specify_phase(_ASSUMED_PHASE_MAP[self._assumed_phase]) - return state - - def _lowlevel_loop_constant( - self, output: CProperty, p1_name: CProperty, p1_arr: Numpy1DArray, p2_name: CProperty, p2_arr: Numpy1DArray - ) -> Numpy1DArray: - state = self._constant_state() - k1 = _cp.get_parameter_index(p1_name) - k2 = _cp.get_parameter_index(p2_name) - out_idx = _cp.get_parameter_index(output) - generate_update_pair = _cp.generate_update_pair - update = state.update - keyed_output = state.keyed_output - - out = np.full(p1_arr.size, np.nan) - for i in range(p1_arr.size): - v1 = p1_arr[i] - v2 = p2_arr[i] - if not (np.isfinite(v1) and np.isfinite(v2)): - continue - try: - pair, a, b = generate_update_pair(k1, v1, k2, v2) - update(pair, a, b) - out[i] = keyed_output(out_idx) - except ValueError as e: - self.check_exception(output, e) - - out[~np.isfinite(out)] = np.nan - return out - - def _evaluate_lowlevel( - self, output: CProperty, points: tuple[tuple[CProperty, float | Numpy1DArray | pl.Expr], ...] - ) -> float | Numpy1DArray | pl.Expr: - p1_name = points[0][0] - p2_name = points[1][0] - - def compute(p1: Numpy1DArray, p2: Numpy1DArray) -> Numpy1DArray: - return self._lowlevel_loop_constant(output, p1_name, p1, p2_name, p2) - - return self._lowlevel_dispatch(output, points, compute) - - def _lowlevel_dispatch( - self, - output: CProperty, - points: tuple[tuple[CProperty, float | Numpy1DArray | pl.Expr], ...], - compute: Callable[[Numpy1DArray, Numpy1DArray], Numpy1DArray], - ) -> float | Numpy1DArray | pl.Expr: - points = tuple((name, self._reduce_single_element(mag)) for name, mag in points) - p1 = points[0][1] - p2 = points[1][1] - all_vals: list[float | Numpy1DArray | pl.Expr] = [p1, p2] - - if any(isinstance(m, pl.Expr) for m in all_vals): - if any(isinstance(m, np.ndarray) for m in all_vals): - raise TypeError("cannot mix numpy array and pl.Expr inputs") - return self._lowlevel_expr(output, points) - - if all(isinstance(m, float) for m in all_vals): - out = compute(np.array([cast("float", p1)]), np.array([cast("float", p2)])) - return float(out[0]) - - arrays = [m for m in all_vals if isinstance(m, np.ndarray)] - shape = arrays[0].shape - n = arrays[0].size - if any(a.shape != shape for a in arrays): - raise ValueError("all array inputs must share the same shape") - - def to_arr(m: float | Numpy1DArray | pl.Expr) -> Numpy1DArray: - if isinstance(m, np.ndarray): - return m.flatten().astype(float) # ty: ignore[invalid-argument-type, unresolved-attribute] - return np.full(n, cast("float", m)) - - p1_arr, p2_arr = to_arr(p1), to_arr(p2) - # large eager arrays: skip the per-row Python loop, use the rust plugin (it - # honors the assumed phase / composition too). - if n >= EAGER_PLUGIN_MIN_SIZE and self._rust_representable(): - named = ((points[0][0], p1_arr), (points[1][0], p2_arr)) - return self._rust_eager(output, named).reshape(shape) - - out = compute(p1_arr, p2_arr) - return out.reshape(shape) - - def _rust_spec(self) -> tuple[Backend, str, list[float] | None, Phase | None] | None: - """Plugin spec ``(backend, fluids, fractions, phase)``, or None if the rust plugin - cannot represent it (an unsupported backend). Name/composition resolution is - delegated to :func:`encomp.coolprop.resolve_fluid_spec` -- the same function - :func:`~encomp.coolprop.fluid` uses -- so both agree.""" - backend, fluids, fractions = resolve_fluid_spec(self.name, self._composition) - if not is_backend(backend): - return None - phase: Phase | None - if self._assumed_phase is None: - phase = None - else: - phase_candidate = f"phase_{self._assumed_phase}" - if not is_phase(phase_candidate): - return None - phase = phase_candidate - return backend, fluids, fractions, phase - - def _rust_representable(self) -> bool: - """Whether the rust plugin can evaluate this fluid's config (the output name is - resolved by CoolProp at runtime, so it never affects representability). - - Humid air is always representable (per-row HAPropsSI in the plugin); a Fluid - is representable unless its backend is one the plugin's spec rejects. - """ - if self.BACKEND["backend"] is HAPropsSI: - return True - return self._rust_spec() is not None - - def _rust_plugin_expr(self, output: CProperty, points: tuple[tuple[str, pl.Expr], ...]) -> pl.Expr: - """The raw encomp.coolprop plugin expr (aliased to ``output``): a + def _plugin_expr(self, output: CProperty, points: tuple[tuple[str, pl.Expr], ...]) -> pl.Expr: + """Build an encomp.coolprop plugin expression aliased to ``output``: a precision-preserving dtype (Float32 in -> Float32 out) with null (not NaN) for failed/out-of-range rows. - Shared by the lazy ``pl.Expr`` path (:meth:`_rust_expr`) and the eager - large-array path (:meth:`_rust_eager`); the latter reads it via ``to_numpy()``, - which maps null back to NaN, so its numpy masking is unchanged. Output property - names are resolved by - CoolProp at runtime, so any valid name works regardless of the static + Shared by lazy expressions and eager arrays; the latter reads it via + ``to_numpy()``, which maps null back to NaN. Output property names are resolved + by CoolProp at runtime, so any valid name works regardless of the static ``FluidParam`` / ``HumidAirParam`` sets. Raises if the plugin is unavailable - (there is no map_batches fallback) or cannot represent the request. + or cannot represent the request. """ from . import coolprop as _cprust if not _cprust.self_check(): raise RuntimeError( "the encomp.coolprop rust plugin failed to load, but it is required to " - "evaluate pl.Expr (lazy) inputs and large eager arrays (there is no " - "fallback). Reinstall encomp with its compiled plugin." + "evaluate arrays and pl.Expr inputs. Reinstall encomp with its compiled plugin." ) names = [p[0] for p in points] exprs = [p[1] for p in points] - if self.BACKEND["backend"] is HAPropsSI: + if self.BACKEND_KIND == "humid_air": n1, n2, n3 = names # the plugin reads each input's property from its (aliased) name expr = _cprust.humid_air( @@ -929,14 +748,9 @@ def _rust_plugin_expr(self, output: CProperty, points: tuple[tuple[str, pl.Expr] exprs[2].alias(n3), ) else: - if self._rust_spec() is None: - raise TypeError( - f"the encomp.coolprop rust plugin cannot represent this evaluation for " - f"fluid {self.name!r} (an unsupported backend)." - ) n1, n2 = names # pass the high-level spec straight through; coolprop.fluid resolves the - # name + composition the same way _rust_spec does (both via resolve_fluid_spec) + # name + composition through the bundled native implementation expr = _cprust.fluid( cast("FluidParam", output), exprs[0].alias(n1), @@ -947,53 +761,25 @@ def _rust_plugin_expr(self, output: CProperty, points: tuple[tuple[str, pl.Expr] ) return expr.alias(output) - def _rust_expr(self, output: CProperty, points: tuple[tuple[CProperty, pl.Expr], ...]) -> pl.Expr: - """The plugin expr for lazy ``pl.Expr`` inputs. The plugin's output dtype already - preserves the input precision (Float32 in -> Float32 out) and emits null (not - NaN) for failed/out-of-range rows directly, so no ``fill_nan(None)`` wrapper is - needed: that wrapper lowers to ``when(is_not_nan).then(x).otherwise(null)``, which - re-evaluates the whole plugin subtree (no common-subexpression elimination) and - runs the CoolProp flash 2-3x.""" - return self._rust_plugin_expr(output, points) - - def _rust_eager(self, output: CProperty, points: tuple[tuple[str, Numpy1DArray], ...]) -> Numpy1DArray: - """Evaluate eager numpy ``points`` through the rust plugin, returning a 1-D - Float64 array with the *same* invalid handling as :meth:`evaluate_multiple` - (non-finite inputs and non-finite results -> NaN). ``points`` are equal-length - arrays; the caller reshapes the result. + def _evaluate_array(self, output: CProperty, points: tuple[tuple[str, Numpy1DArray], ...]) -> Numpy1DArray: + """Evaluate eager NumPy ``points`` through the native plugin, returning a 1-D + Float64 array with non-finite inputs and results mapped to NaN. ``points`` are + equal-length arrays; the caller reshapes the result. """ names = [name for name, _ in points] arrs = [np.ascontiguousarray(arr, dtype=float).ravel() for _, arr in points] # generic column names avoid any collision between the input property names frame = pl.DataFrame({f"__in{i}": a for i, a in enumerate(arrs)}) plugin_points = tuple((names[i], pl.col(f"__in{i}")) for i in range(len(arrs))) - expr = self._rust_plugin_expr(output, plugin_points) + expr = self._plugin_expr(output, plugin_points) out = np.array(frame.select(expr)[output].to_numpy(), dtype=float) - # match evaluate_multiple exactly: rows with any non-finite input are NaN, and - # any inf/_HUGE that slipped through becomes NaN. + # Rows with any non-finite input are NaN, and any inf/_HUGE that slipped + # through becomes NaN. finite_inputs = np.logical_and.reduce([np.isfinite(a) for a in arrs]) out[~finite_inputs] = np.nan out[~np.isfinite(out)] = np.nan return out - def _lowlevel_expr( - self, - output: CProperty, - points: tuple[tuple[CProperty, float | Numpy1DArray | pl.Expr], ...], - ) -> pl.Expr: - # fixed composition and/or assumed phase with pl.Expr inputs: rust-plugin only. - expr_points: tuple[tuple[CProperty, pl.Expr], ...] = tuple( - (name, v if isinstance(v, pl.Expr) else pl.lit(v)) for name, v in points - ) - key = _get_expr_evaluation_cache_key(self, output, expr_points) - cached = _expr_evaluation_cache_get(key) - if cached is not None: - return cached - - expr = self._rust_expr(output, expr_points) - _expr_evaluation_cache_set(key, expr) - return expr - @staticmethod def _reduce_single_element(x: float | Numpy1DArray | pl.Expr) -> float | Numpy1DArray | pl.Expr: if isinstance(x, np.ndarray) and x.size == 1: @@ -1001,32 +787,34 @@ def _reduce_single_element(x: float | Numpy1DArray | pl.Expr) -> float | Numpy1D return x def evaluate_single(self, output: CProperty, *points: tuple[CProperty, float]) -> float: - """Evaluate ``output`` for scalar inputs through the CoolProp backend, returning ``NaN`` on an invalid state.""" + """Evaluate a scalar through the direct native bridge, returning ``NaN`` on an invalid state.""" # A non-finite input cannot fix a state, and CoolProp does not reliably say so: it - # returns 0.0 ("Liquid") for PropsSI("PHASE", "T", nan, ...), the -1 single-phase - # sentinel for "Q", a finite constant for state-independent outputs like TCRIT, and - # HAPropsSI echoes an input straight back. Mask the input here, exactly as - # evaluate_multiple, _lowlevel_loop_constant and _rust_eager already do, so the - # scalar path cannot disagree with the vector/plugin paths. + # can otherwise return a phase sentinel, a state-independent constant, or echo an + # input. Mask here so scalar and batch/plugin paths agree. if not all(np.isfinite(value) for _, value in points): return np.nan - inputs = list(flatten(points)) - - if self._append_name_to_cp_inputs: - inputs.append(self.name) - try: - val = self.BACKEND["backend"](output, *inputs) - - if not isinstance(val, float): - raise TypeError(f"Unexpected value type: {type(val)}, expected float") - - if val == np.inf or val == -np.inf: - val = np.nan - - return val + if self.BACKEND_KIND == "humid_air": + (name1, value1), (name2, value2), (name3, value3) = points + value = _humid_air_scalar(output, name1, value1, name2, value2, name3, value3) + else: + (name1, value1), (name2, value2) = points + value = _fluid_scalar( + output, + name1, + value1, + name2, + value2, + name=self.name, + assume_phase=cast("AssumedPhase | None", self._assumed_phase), + composition=self._composition, + ) + if not np.isfinite(value): + self._warn_coolprop_nan(output, "native evaluation returned no finite result") + return np.nan + return value except ValueError as e: self.check_exception(output, e) @@ -1035,113 +823,24 @@ def evaluate_single(self, output: CProperty, *points: tuple[CProperty, float]) - def evaluate_expression(self, output: CProperty, *points: tuple[CProperty, pl.Expr]) -> pl.Expr: """Build (or reuse) the cached ``pl.Expr`` plugin node that evaluates ``output``.""" - # pl.Expr (lazy) inputs are evaluated only through the GIL-free rust plugin; - # _rust_expr raises if it is unavailable (no map_batches fallback). + # Expressions are evaluated only through the GIL-free native plugin. key = _get_expr_evaluation_cache_key(self, output, points) cached_expr = _expr_evaluation_cache_get(key) if cached_expr is not None: return cached_expr - expr = self._rust_expr(output, points) + # A fill_nan wrapper would duplicate the opaque plugin subtree because Polars + # cannot apply common-subexpression elimination across plugin nodes. + expr = self._plugin_expr(output, points) _expr_evaluation_cache_set(key, expr) return expr - def evaluate_multiple_separately( - self, - output: CProperty, - props: list[CProperty], - arrs_flat_masked: list[Numpy1DArray], - N: int, - ) -> Numpy1DArray: - """Evaluate ``output`` row by row, so a single invalid row yields ``NaN`` rather than failing the batch.""" - - vals: list[float] = [] - - for i in range(N): - arrs_flat_masked_i = [n[i] for n in arrs_flat_masked] - - inputs_i = list(flatten(list(zip(props, arrs_flat_masked_i, strict=False)))) - - if self._append_name_to_cp_inputs: - inputs_i.append(self.name) - - try: - val_i = self.BACKEND["backend"](output, *inputs_i) - - if not isinstance(val_i, float): - raise TypeError(f"Unexpected value type: {type(val_i)}, expected float") - - except ValueError as e: - self.check_exception(output, e) - val_i = np.nan - - vals.append(val_i) - - return np.array(vals) - - def evaluate_multiple(self, output: CProperty, *points: tuple[CProperty, Numpy1DArray]) -> Numpy1DArray: - """Evaluate ``output`` for array inputs in one vectorized backend call; non-finite rows are ``NaN``.""" - - props: list[CProperty] = [pt[0] for pt in points] - arrs = [pt[1] for pt in points] - shape = arrs[0].shape - - arrs_flat = [n.flatten() for n in arrs] - - mask: np.ndarray = np.logical_and.reduce([np.isfinite(n) for n in arrs_flat]) - - def get_empty_like(x: np.ndarray) -> np.ndarray: - empty = np.empty_like(x).astype(float) - empty[:] = np.nan - return empty - - val = get_empty_like(arrs_flat[0]) - - # number of finite (not nan, inf, ...) values - N = mask.astype(int).sum() - - if N > 0: - arrs_flat_masked = [n[mask] for n in arrs_flat] - - inputs = list(flatten(list(zip(props, arrs_flat_masked, strict=False)))) - - if self._append_name_to_cp_inputs: - inputs.append(self.name) - - # this can fail if the numeric values - # are *all* incorrect, for example negative pressure - try: - val_masked = self.BACKEND["backend"](output, *inputs) - - if isinstance(val_masked, float): - raise TypeError(f"Unexpected value type: {type(val_masked)}, expected np.ndarray") - - except ValueError as e: - self.check_exception(output, e) - - if self._evaluate_invalid_separately: - val_masked = self.evaluate_multiple_separately(output, props, arrs_flat_masked, N) - else: - val_masked = get_empty_like(arrs_flat_masked[0]) - - val[mask] = val_masked - - def validate_output(x: Numpy1DArray) -> Numpy1DArray: - x[x == np.inf] = np.nan - x[x == -np.inf] = np.nan - return x.reshape(shape) - - return validate_output(val) - def evaluate( self, output: CProperty, *points: tuple[CProperty, float | Numpy1DArray | pl.Expr] ) -> float | Numpy1DArray | pl.Expr: """Evaluate ``output`` in CoolProp's own unit, dispatching on the magnitude type of ``points``.""" - if self._lowlevel: - return self._evaluate_lowlevel(output, points) - if all(isinstance(pt[1], float) for pt in points): scalar_points = [cast("tuple[CProperty, float]", n) for n in points] return self.evaluate_single(output, *scalar_points) @@ -1190,13 +889,7 @@ def expand_scalars(x: float | np.ndarray) -> Numpy1DArray: points_arr: tuple[tuple[CProperty, Numpy1DArray], ...] = tuple( (p, expand_scalars(cast(Any, v))) for p, v in reduced_points ) - - # large eager arrays go through the GIL-free rust plugin (faster + leaner); - # smaller arrays stay on the vectorized PropsSI path (lower fixed overhead). - if n >= EAGER_PLUGIN_MIN_SIZE and self._rust_representable(): - return self._rust_eager(output, points_arr).reshape(shape) - - return self.evaluate_multiple(output, *points_arr) + return self._evaluate_array(output, points_arr).reshape(shape) def _eager_series_output_dtype(self) -> pl.DataType: # eager pl.Series output preserves polars precision: Float32 only when every @@ -1265,8 +958,7 @@ def get( convert_magnitude: bool = True, ) -> Quantity[Any, MT]: """ - Wraps the CoolProp backend function (``PropsSI``, or ``HAPropsSI`` - for :py:class:`encomp.fluids.HumidAir`), handles input + Evaluates through encomp's bundled native CoolProp artifact and handles input and output with :py:class:`encomp.units.Quantity` objects. Parameters @@ -1276,7 +968,7 @@ def get( points : list[tuple[CProperty, Quantity[Any, MT] | Quantity[Any, float]]] | None Fixed state variables: name and value of the property. The number of points must match the number expected - by the CoolProp backend function. + by the fluid implementation. If None, the points from the ``__init__`` method are used convert_magnitude : bool Whether to convert the output to the same magnitude type as the input, @@ -1413,12 +1105,11 @@ def __init__( def _init_composition(self, name: CName, composition: Composition) -> None: # all name + composition parsing/validation lives in resolve_fluid_spec (the single - # source); store the canonical name + resolved fractions for the low-level/rust paths + # source); store the canonical name + resolved fractions for both native paths backend, fluids, fractions = resolve_fluid_spec(name, composition) assert fractions is not None # a composition= dict always resolves to a mixture self.name = f"{backend}::{fluids}" self._composition = dict(zip(fluids.split("&"), fractions, strict=True)) - self._lowlevel = True def assume_phase(self, phase: AssumedPhase | None) -> Self: """Force the equation of state to assume ``phase``, skipping CoolProp's @@ -1427,15 +1118,13 @@ def assume_phase(self, phase: AssumedPhase | None) -> Self: With ``P, T`` inputs CoolProp normally runs a phase-stability search to decide whether the state is single- or two-phase. For HEOS/GERG mixtures that search dominates the cost (~5 ms/point); assuming a phase you - already know skips it and switches evaluation to the low-level - ``AbstractState`` API -- on the order of 100-1000x faster. + already know skips it while retaining the same native evaluation API. Important caveats: * **Only the HEOS/GERG backends honour it.** ``IF97`` (the default for :class:`Water`) is region-explicit and ignores an assumed phase, so - this call is a no-op there: it emits a warning and keeps the fast - vectorized path (rather than switching to the slower per-point loop). + this call is a no-op there and emits a warning. Use ``Fluid("HEOS::Water", ...)`` if you need an assumed phase for water. * **The assumed phase must be correct for the operating domain.** Forcing a phase the fluid is not actually in does NOT raise -- on HEOS you get @@ -1456,20 +1145,19 @@ def assume_phase(self, phase: AssumedPhase | None) -> Self: ``"supercritical_gas"``, ``"supercritical_liquid"``, ``"twophase"``, or ``None`` to clear. """ - if phase is not None and phase not in _ASSUMED_PHASE_MAP: - raise ValueError(f"unknown phase {phase!r}, expected one of {sorted(_ASSUMED_PHASE_MAP)} or None") + if phase is not None and phase not in ASSUMED_PHASES: + raise ValueError(f"unknown phase {phase!r}, expected one of {sorted(ASSUMED_PHASES)} or None") backend = self.name.split("::", 1)[0] if "::" in self.name else self.name if phase is not None and backend.upper() in _PHASE_IGNORING_BACKENDS: _LOGGER.warning( f"the {backend} backend is region-explicit and ignores an assumed phase; " - f"assume_phase({phase!r}) is a no-op here (the fast vectorized path is kept). " + f"assume_phase({phase!r}) is a no-op here. " "Use a HEOS-backed fluid (e.g. Fluid('HEOS::Water', ...)) to assume a phase." ) return self self._assumed_phase = phase - self._lowlevel = phase is not None or self._composition is not None return self @property @@ -1724,18 +1412,16 @@ def __repr__(self) -> str: class HumidAir(CoolPropFluid[MT]): - """Humid-air wrapper around CoolProp's ``HAPropsSI`` function. + """Humid-air properties from encomp's bundled CoolProp implementation. ``HumidAir.M`` follows CoolProp's humid-air naming and returns dynamic viscosity. This differs from ``Fluid.M``, which returns molar mass. """ - BACKEND = {"backend": HAPropsSI} - _append_name_to_cp_inputs = False - _evaluate_invalid_separately = True + BACKEND_KIND: ClassVar[BackendKind] = "humid_air" STATE_INPUTS: ClassVar[frozenset[str]] = frozenset(HUMID_AIR_INPUTS) - # unit and description for properties in function HAPropsSI + # CoolProp humid-air property names, units, and descriptions PROPERTY_MAP: dict[tuple[CProperty, ...], tuple[str, str]] = { ("B", "Twb", "T_wb", "WetBulb"): ("K", "Wet-Bulb Temperature"), ("C", "cp"): ("J/kg/K", "Mixture specific heat per unit dry air"), @@ -1771,7 +1457,7 @@ class HumidAir(CoolPropFluid[MT]): ALL_PROPERTIES: set[CProperty] = set(flatten(list(PROPERTY_MAP))) - # HAPropsSI has different parameter names + # Humid air has a property namespace distinct from AbstractState fluids. # density is not defined, need to use either Vda (volume per dry air) # or Vha (per humid air) RETURN_UNITS: dict[CProperty, str] = { @@ -1794,8 +1480,7 @@ class HumidAir(CoolPropFluid[MT]): def __init__(self, **kwargs: Unpack[HumidAirState[MT]]) -> None: """ - Interface to the CoolProp function for humid air, - ``CoolProp.CoolProp.HAPropsSI``. + Interface to the bundled CoolProp implementation for humid air. Needs three fixed points instead of two. Parameters diff --git a/encomp/tests/_coolprop_golden.py b/encomp/tests/_coolprop_golden.py new file mode 100644 index 0000000..5382150 --- /dev/null +++ b/encomp/tests/_coolprop_golden.py @@ -0,0 +1,88 @@ +"""Compact CoolProp 8.0.0 references generated before removing its Python binding.""" + +from dataclasses import dataclass +from typing import Literal + +from encomp.coolprop import AssumedPhase, CName, Composition +from encomp.fluids import CProperty + + +@dataclass(frozen=True) +class GoldenCase: + id: str + kind: Literal["fluid", "humid_air", "water"] + inputs: tuple[tuple[CProperty, float, str], ...] + output: CProperty + expected: float + name: CName | None = None + composition: Composition | None = None + phase: AssumedPhase | None = None + + +GOLDEN_CASES = ( + GoldenCase( + id="water-if97-pt", + kind="water", + inputs=(("P", 5e6, "Pa"), ("T", 400.0, "K")), + output="DMASS", + expected=939.9062482796301, + ), + GoldenCase( + id="heos-pure-pt", + kind="fluid", + name="HEOS::CarbonDioxide", + inputs=(("P", 5e6, "Pa"), ("T", 310.0, "K")), + output="HMASS", + # Fluid converts CoolProp's J/kg result to its preferred kJ/kg unit. + expected=462.41253710627754, + ), + GoldenCase( + id="incompressible-concentration", + kind="fluid", + name="INCOMP::MEG[0.5]", + inputs=(("P", 2e5, "Pa"), ("T", 300.0, "K")), + output="DMASS", + expected=1061.1793077204613, + ), + GoldenCase( + id="embedded-mixture", + kind="fluid", + name="HEOS::CarbonDioxide[0.7]&Oxygen[0.3]", + inputs=(("P", 5e6, "Pa"), ("T", 300.0, "K")), + output="DMASS", + expected=97.49264109222536, + ), + GoldenCase( + id="explicit-composition", + kind="fluid", + name="HEOS", + composition={"CarbonDioxide": 0.6, "Oxygen": 0.4}, + inputs=(("P", 5e6, "Pa"), ("T", 300.0, "K")), + output="DMASS", + expected=91.2386351314185, + ), + GoldenCase( + id="assumed-gas-phase", + kind="fluid", + name="HEOS::Water", + phase="gas", + inputs=(("P", 1e5, "Pa"), ("T", 500.0, "K")), + output="DMASS", + expected=0.435140075089223, + ), + GoldenCase( + id="non-pt-input-pair", + kind="fluid", + name="HEOS::Water", + inputs=(("P", 2e6, "Pa"), ("H", 749714.8321521956, "J/kg")), + output="DMASS", + expected=891.0411792241698, + ), + GoldenCase( + id="humid-air", + kind="humid_air", + inputs=(("P", 101325.0, "Pa"), ("T", 300.0, "K"), ("R", 0.5, "")), + output="W", + expected=0.01109552970536823, + ), +) diff --git a/encomp/tests/test_coolprop_plugin.py b/encomp/tests/test_coolprop_plugin.py index 6631f78..0518e7f 100644 --- a/encomp/tests/test_coolprop_plugin.py +++ b/encomp/tests/test_coolprop_plugin.py @@ -12,14 +12,13 @@ from pathlib import Path from typing import Any, cast -import CoolProp.CoolProp as _CP import numpy as np import polars as pl import pytest from encomp import coolprop as cp -CP: Any = _CP +CP: Any = pytest.importorskip("CoolProp.CoolProp", reason="Python CoolProp oracle group is not installed") RTOL = 1e-5 @@ -45,17 +44,9 @@ def _plugin_built() -> bool: def _bundled_lib_expected_version() -> str: - # the bundled libCoolProp is built at the LOWER BOUND of the coolprop requirement - # (scripts/build_libcoolprop.py derives its git tag from the same bound); read it - # from the installed package metadata so this also holds for installed wheels - import importlib.metadata - import re - - requires = importlib.metadata.requires("encomp") or [] - requirement = next(r for r in requires if r.startswith("coolprop")) - match = re.search(r">=\s*([\w.]+)", requirement) - assert match, f"no lower bound in the coolprop requirement: {requirement!r}" - return match.group(1) + from encomp.coolprop._build_info import BUNDLED_COOLPROP_VERSION + + return BUNDLED_COOLPROP_VERSION def test_self_check_version_and_lib_path() -> None: diff --git a/encomp/tests/test_fluids.py b/encomp/tests/test_fluids.py index ae9e26b..453e749 100644 --- a/encomp/tests/test_fluids.py +++ b/encomp/tests/test_fluids.py @@ -2,7 +2,6 @@ # pyright: reportConstantRedefinition=false, reportPrivateUsage=false import logging -from collections.abc import Callable from typing import Any, assert_type, cast import numpy as np @@ -239,7 +238,10 @@ def test_humid_air_out_of_range_logs_warning(caplog: pytest.LogCaptureFixture, m assert np.isnan(float(value.m)) assert 'CoolProp could not calculate "W" for fluid "Humid air"' in caplog.text - assert "outside the range of validity" in caplog.text + # The C HAPropsSI interface exposes details only through a process-global error + # outbox. The native bridge deliberately does not read it, so concurrent failures + # cannot attach another thread's message to this warning. + assert "native evaluation returned no finite result" in caplog.text caplog.clear() monkeypatch.setattr(SETTINGS, "ignore_coolprop_warnings", True) @@ -617,21 +619,21 @@ def test_polars_fluids() -> None: def _count_water_h_evaluations(monkeypatch: pytest.MonkeyPatch) -> list[int]: - # pl.Expr inputs are built into a rust-plugin expr by _rust_expr; the expr cache - # dedups identical requests *before* that call, so counting _rust_expr invocations + # pl.Expr inputs are built into a native plugin expression; the expression cache + # deduplicates identical requests before that call, so counting invocations # counts the distinct (cache-missing) expressions actually built. calls = [0] - original_rust_expr = CoolPropFluid._rust_expr + original_plugin_expr = CoolPropFluid._plugin_expr - def counted_rust_expr( + def counted_plugin_expr( self: CoolPropFluid[pl.Expr], output: str, points: tuple[tuple[str, pl.Expr], ...] ) -> pl.Expr: if getattr(self, "name", "") == "IF97::Water" and output == "H": calls[0] += 1 - return cast(pl.Expr, cast(Any, original_rust_expr)(self, output, points)) + return cast(pl.Expr, cast(Any, original_plugin_expr)(self, output, points)) - monkeypatch.setattr(CoolPropFluid, "_rust_expr", counted_rust_expr) + monkeypatch.setattr(CoolPropFluid, "_plugin_expr", counted_plugin_expr) return calls @@ -693,8 +695,7 @@ def test_polars_fluids_rust_backend() -> None: T = np.linspace(300.0, 500.0, 7) df = pl.DataFrame({"P": P, "T": T}) - # pl.Expr inputs go through the rust plugin; eager numpy inputs go through the - # Python CoolProp path. The two must match to float precision. + # pl.Expr and eager NumPy inputs both use the native batch path and must match. clear_expr_evaluation_cache() w = Water[pl.Expr](P=Q(pl.col("P"), "Pa"), T=Q(pl.col("T"), "K")) rust = df.select(w.D.m.alias("D"), w.H.m.alias("H"), w.S.m.alias("S"), w.C.m.alias("C")) @@ -710,86 +711,74 @@ def test_polars_fluids_rust_backend() -> None: assert bad.is_null().all() -def _eager_both_paths( - monkeypatch: pytest.MonkeyPatch, build: Callable[[], Q[Any, Any]] -) -> tuple[np.ndarray, np.ndarray]: - # run build() once forcing the Python path and once forcing the rust plugin, by - # moving the eager size threshold above / below the input length. - monkeypatch.setattr("encomp.fluids.EAGER_PLUGIN_MIN_SIZE", 10**18) - py = build().m - monkeypatch.setattr("encomp.fluids.EAGER_PLUGIN_MIN_SIZE", 1) - ru = build().m - return py, ru - - -def _assert_eager_parity(label: str, py: np.ndarray, ru: np.ndarray) -> None: - assert py.dtype == np.float64 and ru.dtype == np.float64, f"{label}: dtype {py.dtype}/{ru.dtype}" - assert not np.isinf(py).any() and not np.isinf(ru).any(), f"{label}: inf present" - assert np.array_equal(np.isnan(py), np.isnan(ru)), f"{label}: NaN positions differ" - finite = ~np.isnan(py) - assert np.allclose(py[finite], ru[finite], rtol=1e-9, atol=1e-12), f"{label}: values differ" - - -def test_eager_plugin_parity(monkeypatch: pytest.MonkeyPatch) -> None: - # eager numpy arrays >= EAGER_PLUGIN_MIN_SIZE route through the rust plugin; the - # result must match the Python CoolProp path exactly -- dtype (float64), value, and - # NaN/inf/invalid handling -- including injected non-finite / out-of-range rows. +def test_eager_plugin_matches_direct_scalar_bridge() -> None: + # Every eager array size now uses the plugin. Compare it with the independent + # direct-PyO3 scalar interface, including invalid rows and special configs. rng = np.random.default_rng(0) - n = 1500 + n = 24 P = np.full(n, 50e5) T = 300.0 + 250.0 * rng.random(n) T[0] = np.nan # non-finite input -> NaN P[1] = -1.0 # invalid (negative pressure) -> NaN T[2] = 50.0 # below the IF97 range -> NaN - fluid_cases = { - "IF97": lambda: Water(P=Q(P, "Pa"), T=Q(T, "K")).D, - "HEOS": lambda: Fluid("HEOS::Water", P=Q(P, "Pa"), T=Q(T, "K")).D, - "assumed-phase": lambda: Fluid("HEOS::Water", P=Q(P, "Pa"), T=Q(T, "K")).assume_phase("gas").D, - "composition": lambda: ( - Fluid("HEOS", composition={"CO2": 0.5, "O2": 0.5}, P=Q(P, "Pa"), T=Q(T, "K")) + def compare_array_and_scalars(array: np.ndarray, scalars: np.ndarray) -> None: + assert array.dtype == np.float64 + assert not np.isinf(array).any() + assert np.array_equal(np.isnan(array), np.isnan(scalars)) + finite = np.isfinite(array) + assert np.allclose(array[finite], scalars[finite], rtol=1e-9, atol=1e-12) + + array = cast(np.ndarray, Water(P=Q(P, "Pa"), T=Q(T, "K")).D.m) + scalars = np.array([Water(P=Q(float(p), "Pa"), T=Q(float(t), "K")).D.m for p, t in zip(P, T, strict=True)]) + compare_array_and_scalars(array, scalars) + + array = cast( + np.ndarray, + Fluid("HEOS", composition={"CO2": 0.5, "O2": 0.5}, P=Q(P, "Pa"), T=Q(T, "K")) + .assume_phase("supercritical_gas") + .D.m, + ) + scalars = np.array( + [ + Fluid("HEOS", composition={"CO2": 0.5, "O2": 0.5}, P=Q(float(p), "Pa"), T=Q(float(t), "K")) .assume_phase("supercritical_gas") - .D - ), - } - for label, build in fluid_cases.items(): - py, ru = _eager_both_paths(monkeypatch, build) - _assert_eager_parity(label, py, ru) + .D.m + for p, t in zip(P, T, strict=True) + ] + ) + compare_array_and_scalars(array, scalars) - # humid air: the plugin loops HAPropsSI per row, matching evaluate_multiple_separately Pa = np.full(n, 101325.0) Tdb = 290.0 + 20.0 * rng.random(n) R = 0.2 + 0.6 * rng.random(n) Tdb[0] = np.nan R[1] = 5.0 # relative humidity > 1 -> invalid - py, ru = _eager_both_paths(monkeypatch, lambda: HumidAir(P=Q(Pa, "Pa"), T=Q(Tdb, "K"), R=Q(R, "")).W) - _assert_eager_parity("humid-air", py, ru) + array = cast(np.ndarray, HumidAir(P=Q(Pa, "Pa"), T=Q(Tdb, "K"), R=Q(R, "")).W.m) + scalars = np.array( + [ + HumidAir(P=Q(float(p), "Pa"), T=Q(float(t), "K"), R=Q(float(r), "")).W.m + for p, t, r in zip(Pa, Tdb, R, strict=True) + ] + ) + compare_array_and_scalars(array, scalars) -def test_eager_plugin_threshold_dispatch(monkeypatch: pytest.MonkeyPatch) -> None: - # below EAGER_PLUGIN_MIN_SIZE -> Python (evaluate_multiple); at/above -> rust plugin +def test_eager_plugin_threshold_is_deprecated_noop(monkeypatch: pytest.MonkeyPatch) -> None: monkeypatch.setattr("encomp.fluids.EAGER_PLUGIN_MIN_SIZE", 1000) - counts = {"rust": 0, "multiple": 0} - orig_rust = CoolPropFluid._rust_eager - orig_mult = CoolPropFluid.evaluate_multiple + calls = 0 + original_array = CoolPropFluid._evaluate_array def spy_rust(self: CoolPropFluid[Any], output: str, points: tuple[tuple[str, np.ndarray], ...]) -> np.ndarray: - counts["rust"] += 1 - return cast(np.ndarray, cast(Any, orig_rust)(self, output, points)) - - def spy_mult(self: CoolPropFluid[Any], output: str, *points: tuple[str, np.ndarray]) -> np.ndarray: - counts["multiple"] += 1 - return cast(np.ndarray, cast(Any, orig_mult)(self, output, *points)) - - monkeypatch.setattr(CoolPropFluid, "_rust_eager", spy_rust) - monkeypatch.setattr(CoolPropFluid, "evaluate_multiple", spy_mult) + nonlocal calls + calls += 1 + return cast(np.ndarray, cast(Any, original_array)(self, output, points)) - Water(P=Q(np.full(999, 50e5), "Pa"), T=Q(np.full(999, 400.0), "K")).D.m - assert counts == {"rust": 0, "multiple": 1} - - counts["rust"] = counts["multiple"] = 0 - Water(P=Q(np.full(1000, 50e5), "Pa"), T=Q(np.full(1000, 400.0), "K")).D.m - assert counts == {"rust": 1, "multiple": 0} + monkeypatch.setattr(CoolPropFluid, "_evaluate_array", spy_rust) + Water(P=Q(np.full(2, 50e5), "Pa"), T=Q(np.full(2, 400.0), "K")).D.m + monkeypatch.setattr("encomp.fluids.EAGER_PLUGIN_MIN_SIZE", 10**18) + Water(P=Q(np.full(3, 50e5), "Pa"), T=Q(np.full(3, 400.0), "K")).D.m + assert calls == 2 def test_polars_dtype_preservation() -> None: @@ -894,13 +883,13 @@ def test_assume_phase() -> None: Fluid(mix, P=Q(50.0, "bar"), T=Q(350.0, "degC")).assume_phase(cast(Any, "plasma")) -def test_assume_phase_changes_value(monkeypatch: pytest.MonkeyPatch) -> None: +def test_assume_phase_changes_value() -> None: # test_assume_phase only checks single-phase regions where assumed == auto, so a # regression that silently dropped assume_phase would pass it. Here a DISCRIMINATING # state (HEOS::Water at 1 bar, 110 C) is auto-phase GAS (steam, ~0.57 kg/m3); forcing # "liquid" must return the metastable subcooled-liquid root (~950 kg/m3). This proves - # assume_phase actually takes effect, and that the scalar low-level path and the - # large-eager rust path agree on the forced value. (An independent CoolProp + # assume_phase actually takes effect, and that the direct scalar and batch paths + # agree on the forced value. (An independent CoolProp # specify_phase reference for the same state lives in test_coolprop_plugin.) auto = Water(P=Q(1.0, "bar"), T=Q(110.0, "degC")).D.to("kg/m3") assert float(auto.m) < 10.0 # auto-phase water here is steam @@ -909,8 +898,7 @@ def test_assume_phase_changes_value(monkeypatch: pytest.MonkeyPatch) -> None: assert float(scalar.m) > 900.0 # forcing liquid gives the subcooled-liquid root assert float(scalar.m) != approx(float(auto.m), rel=0.5) # unmistakably not the auto (steam) value - # same config, large-eager array -> rust plugin path; must match the scalar low-level path - monkeypatch.setattr("encomp.fluids.EAGER_PLUGIN_MIN_SIZE", 1000) + # same config, eager array -> plugin path; must match the direct scalar path n = 1200 eager = ( Fluid("HEOS::Water", P=Q(np.full(n, 1e5), "Pa"), T=Q(np.full(n, 383.15), "K")) @@ -977,8 +965,8 @@ def test_composition() -> None: def test_assume_phase_if97_noop(caplog: pytest.LogCaptureFixture) -> None: # Water uses IF97, which is region-explicit and ignores an assumed phase: - # assume_phase must be a no-op (warn + keep the fast vectorized path), not a - # pessimisation. The returned value must equal the auto-phase value. + # assume_phase must be a no-op with a warning. The returned value must equal + # the automatic-phase value. auto = Water(P=Q(50.0, "bar"), T=Q(150.0, "degC")).D w = Water(P=Q(50.0, "bar"), T=Q(150.0, "degC")) with caplog.at_level(logging.WARNING): @@ -1046,13 +1034,12 @@ def test_composition_rejects_non_float_fractions() -> None: def test_incompressible_mixture_all_paths_agree() -> None: # regression: an INCOMP concentration mixture (INCOMP::MEG[0.5], ...) used to lose its - # concentration on the rust path (expr / large-eager), silently diverging from the - # PropsSI scalar/small-eager path. All paths must agree, mass- AND volume-based. + # concentration identically in scalar, eager-array, and expression forms. assert encomp_coolprop.self_check() for name, temp in [("INCOMP::MEG[0.5]", 300.0), ("INCOMP::MITSW[0.035]", 290.0), ("INCOMP::AEG[0.4]", 300.0)]: pressure = np.full(6, 3e5) temperature = np.linspace(temp - 12.0, temp + 12.0, 6) - ref = Fluid(name, P=Q(pressure, "Pa"), T=Q(temperature, "K")).D.m # small-eager -> PropsSI + ref = Fluid(name, P=Q(pressure, "Pa"), T=Q(temperature, "K")).D.m clear_expr_evaluation_cache() expr = ( @@ -1062,10 +1049,10 @@ def test_incompressible_mixture_all_paths_agree() -> None: ) assert np.asarray(expr, dtype=float) == approx(ref, rel=1e-4), f"expr {name}" - clear_expr_evaluation_cache() # large-eager (>= EAGER_PLUGIN_MIN_SIZE) -> plugin + clear_expr_evaluation_cache() big = Fluid(name, P=Q(np.full(1200, 3e5), "Pa"), T=Q(np.full(1200, temp), "K")).D.m small = float(Fluid(name, P=Q(3e5, "Pa"), T=Q(temp, "K")).D.m) - assert float(np.asarray(big, dtype=float)[0]) == approx(small, rel=1e-4), f"large-eager {name}" + assert float(np.asarray(big, dtype=float)[0]) == approx(small, rel=1e-4), f"array {name}" def test_dimensional_property_missing_is_nan_never_zero() -> None: @@ -1093,14 +1080,14 @@ def check_polars(series: pl.Series, n_missing: int) -> None: present = series.drop_nulls().to_numpy() assert np.all(np.isfinite(present)) and not np.any(present == 0.0), "present rows: non-zero, never 0.0" - reps = 300 # 4 * 300 = 1200 >= EAGER_PLUGIN_MIN_SIZE -> plugin + reps = 300 # exercise the same native array path with a larger vector big_p, big_t, big_missing = np.tile(pressure, reps), np.tile(temperature, reps), int(missing.sum()) * reps - # numpy magnitude: eager (PropsSI) and large-eager (plugin) + # NumPy magnitudes of different lengths use the same native array path. check_numpy(Water(P=Q(pressure, "Pa"), T=Q(temperature, "K")).D.m, missing) check_numpy(Water(P=Q(big_p, "Pa"), T=Q(big_t, "K")).D.m, np.tile(missing, reps)) - # pl.Series magnitude: small (PropsSI) and large (plugin) -- missing surfaces as null + # pl.Series magnitudes of different lengths use the same native array path. check_polars( Water[pl.Series](P=Q(pl.Series(pressure), "Pa"), T=Q(pl.Series(temperature), "K")).D.m, int(missing.sum()) ) diff --git a/encomp/tests/test_native_coolprop.py b/encomp/tests/test_native_coolprop.py new file mode 100644 index 0000000..cbf4796 --- /dev/null +++ b/encomp/tests/test_native_coolprop.py @@ -0,0 +1,198 @@ +"""Tests for the direct PyO3 side of encomp's single native CoolProp artifact.""" + +# pyright: reportPrivateUsage=false + +from __future__ import annotations + +import importlib.metadata +import math +import subprocess +import sys +from concurrent.futures import ThreadPoolExecutor +from typing import Any, Literal, cast + +import numpy as np +import polars as pl +import pytest +from pytest import approx, raises + +from encomp import coolprop as cp +from encomp.coolprop._build_info import BUNDLED_COOLPROP_VERSION +from encomp.fluids import Fluid, HumidAir, Water +from encomp.units import Quantity as Q +from encomp.utypes import Numpy1DArray + +from ._coolprop_golden import GOLDEN_CASES, GoldenCase + + +def _golden_result( + case: GoldenCase, mode: Literal["expr", "numpy", "scalar", "series"] +) -> float | Numpy1DArray | pl.Series: + frame_data: dict[str, list[float]] = {} + points: dict[str, Any] = {} + for name, value, unit in case.inputs: + if mode == "scalar": + magnitude: Any = value + elif mode == "numpy": + magnitude = np.array([value, value]) + elif mode == "series": + magnitude = pl.Series(name, [value, value]) + else: + frame_data[name] = [value, value] + magnitude = pl.col(name) + points[name] = Q(magnitude, unit) + + if case.kind == "water": + state: Any = cast(Any, Water)(**points) + elif case.kind == "humid_air": + state = cast(Any, HumidAir)(**points) + else: + assert case.name is not None + state = cast(Any, Fluid)(case.name, composition=case.composition, **points) + if case.phase is not None: + state.assume_phase(case.phase) + + result = state.get(case.output).m + if mode == "expr": + return pl.DataFrame(frame_data).select(result.alias("result"))["result"] + return cast("float | Numpy1DArray | pl.Series", result) + + +@pytest.mark.parametrize("case", GOLDEN_CASES, ids=lambda case: case.id) +@pytest.mark.parametrize("mode", ("scalar", "numpy", "series", "expr")) +def test_committed_golden_cases(case: GoldenCase, mode: Literal["expr", "numpy", "scalar", "series"]) -> None: + result = _golden_result(case, mode) + if mode == "scalar": + assert isinstance(result, float) + assert result == approx(case.expected, rel=5e-10) + return + + if isinstance(result, pl.Series): + assert result.dtype == pl.Float64 + values = result.to_numpy() + else: + assert isinstance(result, np.ndarray) + assert result.dtype == np.float64 + values = result + assert values.shape == (2,) + assert values == approx(np.full(2, case.expected), rel=5e-10) + + +def test_native_metadata_and_pair_resolution() -> None: + native = cp._native() + assert native.lib_version() == BUNDLED_COOLPROP_VERSION + assert native.parameter_information("P", "units") == "Pa" + assert native.parameter_index("D") == native.parameter_index("DMASS") + + for first, second in (("P", "T"), ("P", "H"), ("P", "Q"), ("D", "T"), ("H", "S")): + pair, swap = native.resolve_input_pair(first, second) + reverse_pair, reverse_swap = native.resolve_input_pair(second, first) + assert pair > 0 and reverse_pair == pair + assert swap is not reverse_swap + + with raises(ValueError, match="unsupported CoolProp input pair"): + native.resolve_input_pair("P", "P") + with raises(ValueError, match="unknown parameter"): + native.parameter_index("NOT_A_PROPERTY") + + +def test_native_name_parsing_and_validation() -> None: + native = cp._native() + assert native.resolve_fluid_name("Water") == ("HEOS", "Water", None) + assert native.resolve_fluid_name("IF97::Water") == ("IF97", "Water", None) + assert native.resolve_fluid_name("HEOS::CO2[0.5]&O2[0.5]") == ("HEOS", "CO2&O2", [0.5, 0.5]) + assert native.resolve_fluid_name("INCOMP::MEG[0.5]") == ("INCOMP", "MEG", [0.5]) + assert native.resolve_fluid_name("INCOMP::EG-20%") == ("INCOMP", "EG", [0.2]) + + native.validate_fluid("IF97", "Water") + native.validate_fluid("HEOS", "CO2&O2", [0.5, 0.5]) + with raises(ValueError, match="factory"): + native.validate_fluid("HEOS", "DefinitelyNotAFluid") + + +def test_direct_scalar_and_batch_paths_agree() -> None: + pressure = 5e6 + temperature = 400.0 + scalar = Water(P=Q(pressure, "Pa"), T=Q(temperature, "K")).D.m + eager = Water(P=Q(np.array([pressure, pressure]), "Pa"), T=Q(np.array([temperature, temperature]), "K")).D.m + assert scalar == approx(float(eager[0]), rel=1e-12) + + scalar_heos = Fluid("HEOS::CarbonDioxide", P=Q(pressure, "Pa"), T=Q(temperature, "K")).D.m + eager_heos = Fluid( + "HEOS::CarbonDioxide", + P=Q(np.array([pressure, pressure]), "Pa"), + T=Q(np.array([temperature, temperature]), "K"), + ).D.m + assert scalar_heos == approx(float(eager_heos[0]), rel=1e-12) + + scalar_ha = HumidAir(P=Q(101325.0, "Pa"), T=Q(300.0, "K"), R=Q(0.5, "")).W.m + eager_ha = HumidAir( + P=Q(np.array([101325.0, 101325.0]), "Pa"), + T=Q(np.array([300.0, 300.0]), "K"), + R=Q(np.array([0.5, 0.5]), ""), + ).W.m + assert scalar_ha == approx(float(eager_ha[0]), rel=1e-12) + + +def test_thread_local_cache_eviction_and_destruction() -> None: + native = cp._native() + native.clear_scalar_cache() + _, freed_before = native.handle_counts() + + fluids = ( + "Water", + "CarbonDioxide", + "Nitrogen", + "Oxygen", + "Methane", + "Ammonia", + "R134a", + "n-Propane", + "Hydrogen", + "Argon", + "Ethane", + "CarbonMonoxide", + "Helium", + "Neon", + "Krypton", + "Xenon", + "SulfurHexafluoride", + ) + for fluid in fluids: + value = Fluid(f"HEOS::{fluid}", P=Q(101325.0, "Pa"), T=Q(300.0, "K")).D.m + assert math.isfinite(value) + + size, capacity, _hits, misses, evictions = native.scalar_cache_info() + assert size == capacity == 16 + assert misses == len(fluids) + assert evictions == 1 + + native.clear_scalar_cache() + assert native.scalar_cache_info()[0] == 0 + _, freed_after = native.handle_counts() + assert freed_after - freed_before >= len(fluids) + + +def test_concurrent_scalar_success_and_failure_are_isolated() -> None: + def evaluate(index: int) -> float: + pressure = 5e6 if index % 5 else -1.0 + value = Fluid("HEOS::CarbonDioxide", P=Q(pressure, "Pa"), T=Q(300.0, "K")).D.m + assert isinstance(value, float) + return value + + with ThreadPoolExecutor(max_workers=8) as pool: + values = list(pool.map(evaluate, range(200))) + + for index, value in enumerate(values): + assert math.isnan(value) if index % 5 == 0 else math.isfinite(value) + + +def test_runtime_import_does_not_load_python_coolprop() -> None: + code = ( + "import sys; import encomp.coolprop, encomp.fluids; " + "assert not any(n == 'CoolProp' or n.startswith('CoolProp.') for n in sys.modules)" + ) + subprocess.run([sys.executable, "-c", code], check=True) + + requirements = importlib.metadata.requires("encomp") or [] + assert not any(requirement.lower().startswith("coolprop") for requirement in requirements) diff --git a/encomp/tests/test_registry_conformance.py b/encomp/tests/test_registry_conformance.py index d6dbcf4..439a4b2 100644 --- a/encomp/tests/test_registry_conformance.py +++ b/encomp/tests/test_registry_conformance.py @@ -11,8 +11,6 @@ import typing from collections.abc import Iterable -import CoolProp.CoolProp as _CoolProp - from .. import coolprop as cp from ..fluids import CoolPropFluid, Fluid, FluidState, HumidAir, HumidAirState from ..units import Quantity @@ -80,10 +78,10 @@ def test_typed_properties_are_known_names() -> None: assert name in cls.ALL_PROPERTIES, f"{cls.__name__}.{name} is not in its PROPERTY_MAP" -def test_fluid_params_are_valid_per_coolprop() -> None: - # CoolProp's parameter index is the runtime authority for the fluid namespace; +def test_fluid_params_are_valid_per_bundled_coolprop() -> None: + # The bundled library's parameter index is the runtime authority for the fluid namespace; # every FluidParam Literal must resolve (guards against typos in the Literal). # (HAPropsSI has no equivalent lookup, so the humid-air names are locked only # against fluids.HumidAir.PROPERTY_MAP above.) for name in sorted(cp.FLUID_PARAMS): - _CoolProp.get_parameter_index(name) # raises ValueError for unknown names + cp._native().parameter_index(name) # pyright: ignore[reportPrivateUsage] diff --git a/encomp/tests/test_rust_parity.py b/encomp/tests/test_rust_parity.py index b48080c..f610985 100644 --- a/encomp/tests/test_rust_parity.py +++ b/encomp/tests/test_rust_parity.py @@ -1,19 +1,19 @@ -"""Comprehensive parity between the rust plugin and the Python CoolProp path. +"""Comprehensive parity between encomp's native paths and Python CoolProp oracle. -``pl.Expr`` (lazy) inputs are evaluated by the encomp.coolprop rust plugin; eager -numpy inputs go through the Python CoolProp path. Both evaluate CoolProp 8.0 and -cast to Float32, so for the same fluid / property / inputs the results should match -to float precision. We compare across many fluids (IF97 water, a broad set of HEOS -pure fluids, and a mixture), several properties, and (P,T) grids spanning liquid / +``pl.Expr`` inputs are evaluated by the encomp.coolprop Rust plugin; eager arrays +use that same batch path. The optional Python package is only a numerical oracle. +Both evaluate CoolProp 8.0, so results should match closely across IF97 and HEOS +pure fluids, a mixture, several properties, and (P,T) grids spanning liquid / gas / supercritical (reduced coordinates from each fluid's critical point), plus a few non-PT input pairs. The finite overlap is compared with a relative tolerance. """ +# pyright: reportPrivateUsage=false + from __future__ import annotations -from typing import Any +from typing import Any, cast -import CoolProp.CoolProp as _CP import numpy as np import polars as pl import pytest @@ -22,9 +22,8 @@ from encomp.fluids import Fluid, clear_expr_evaluation_cache from encomp.units import Quantity as Q -# CoolProp.CoolProp is a compiled, untyped extension module; alias it as Any so -# its dynamic functions (PropsSI, ...) do not surface unknown-type errors. -CP: Any = _CP +# Optional test-only oracle; absent from runtime and installed-wheel environments. +CP: Any = pytest.importorskip("CoolProp.CoolProp", reason="Python CoolProp oracle group is not installed") # skip (not fail) when the plugin binaries simply have not been built in this @@ -76,44 +75,130 @@ def _grid(name: str, kwargs: dict[str, Any]) -> pl.DataFrame: return pl.DataFrame({"P": pp.ravel(), "T": tt.ravel()}) -def _evaluate(name: str, kwargs: dict[str, Any], prop: str, mode: str, df: pl.DataFrame) -> np.ndarray: +def _evaluate_native(name: str, kwargs: dict[str, Any], prop: str, df: pl.DataFrame) -> np.ndarray: clear_expr_evaluation_cache() - if mode == "rust": # pl.Expr inputs -> rust plugin - fluid = Fluid(name, P=Q(pl.col("P"), "Pa"), T=Q(pl.col("T"), "K"), **kwargs) - return df.select(getattr(fluid, prop).m.alias("x"))["x"].to_numpy().astype(float) - # python reference: eager numpy inputs -> Python CoolProp path - fluid = Fluid(name, P=Q(df["P"].to_numpy(), "Pa"), T=Q(df["T"].to_numpy(), "K"), **kwargs) - return np.asarray(getattr(fluid, prop).m, dtype=float) + expression = cast(Any, encomp_coolprop.fluid)(prop, "P", "T", name=name, **kwargs) + return df.select(expression.alias("x"))["x"].to_numpy().astype(float) + + +def _evaluate_oracle(name: str, kwargs: dict[str, Any], prop: str, df: pl.DataFrame) -> np.ndarray: + # Independent test-only Python CoolProp numerical reference. Explicit + # composition is rendered into the equivalent name syntax for PropsSI. + composition = kwargs.get("composition") + if composition is not None: + components = "&".join(f"{species}[{fraction}]" for species, fraction in composition.items()) + name = f"{name}::{components}" + return np.asarray(CP.PropsSI(prop, "P", df["P"].to_numpy(), "T", df["T"].to_numpy(), name), dtype=float) @pytest.mark.parametrize(("label", "name", "kwargs"), _CONFIGS, ids=[c[0] for c in _CONFIGS]) -def test_rust_python_parity(label: str, name: str, kwargs: dict[str, Any]) -> None: +def test_native_python_oracle_parity(label: str, name: str, kwargs: dict[str, Any]) -> None: assert encomp_coolprop.self_check() df = _grid(name, kwargs) for prop in PROPS: - py = _evaluate(name, kwargs, prop, "python", df) - ru = _evaluate(name, kwargs, prop, "rust", df) - assert py.shape == ru.shape + oracle = _evaluate_oracle(name, kwargs, prop, df) + native = _evaluate_native(name, kwargs, prop, df) + assert oracle.shape == native.shape # the two paths must agree on WHICH points are finite -- otherwise a path that # returns NaN where the other returns a value would be silently excluded below - finite_py, finite_ru = np.isfinite(py), np.isfinite(ru) - assert np.array_equal(finite_py, finite_ru), ( - f"{label}.{prop}: finite/NaN masks differ ({int((finite_py != finite_ru).sum())} points)" + finite_oracle, finite_native = np.isfinite(oracle), np.isfinite(native) + assert np.array_equal(finite_oracle, finite_native), ( + f"{label}.{prop}: finite/NaN masks differ ({int((finite_oracle != finite_native).sum())} points)" ) - both = finite_py & finite_ru + both = finite_oracle & finite_native assert both.sum() >= 8, f"{label}.{prop}: only {int(both.sum())} finite-overlap points" - rel = np.abs(ru[both] - py[both]) / np.maximum(np.abs(py[both]), 1e-9) + rel = np.abs(native[both] - oracle[both]) / np.maximum(np.abs(oracle[both]), 1e-9) worst = int(np.argmax(rel)) assert rel.max() <= RTOL, ( - f"{label}.{prop}: max rel {rel.max():.2e} (rust {ru[both][worst]:.6g} vs python {py[both][worst]:.6g})" + f"{label}.{prop}: max rel {rel.max():.2e} " + f"(native {native[both][worst]:.6g} vs oracle {oracle[both][worst]:.6g})" ) -def test_rust_python_parity_input_pairs() -> None: - """Parity for non-PT input pairs (P,H and P,S) -- the rust path resolves the - canonical pair + column order, the python path uses PropsSI.""" +_INPUT_PAIR_PARAMETERS = ( + ("Q", "T"), + ("Qmass", "T"), + ("P", "Q"), + ("P", "Qmass"), + ("P", "T"), + ("Dmolar", "T"), + ("Dmass", "T"), + ("Hmolar", "T"), + ("Hmass", "T"), + ("Smolar", "T"), + ("Smass", "T"), + ("T", "Umolar"), + ("T", "Umass"), + ("Dmass", "Hmass"), + ("Dmolar", "Hmolar"), + ("Dmass", "Smass"), + ("Dmolar", "Smolar"), + ("Dmass", "Umass"), + ("Dmolar", "Umolar"), + ("Dmass", "P"), + ("Dmolar", "P"), + ("Dmass", "Q"), + ("Dmass", "Qmass"), + ("Dmolar", "Q"), + ("Dmolar", "Qmass"), + ("Hmass", "P"), + ("Hmolar", "P"), + ("P", "Smass"), + ("P", "Smolar"), + ("P", "Umass"), + ("P", "Umolar"), + ("Hmass", "Smass"), + ("Hmolar", "Smolar"), + ("Smass", "Umass"), + ("Smolar", "Umolar"), +) + + +@pytest.mark.parametrize(("name1", "name2"), _INPUT_PAIR_PARAMETERS) +def test_native_input_pair_table_matches_python_oracle(name1: str, name2: str) -> None: + """Every reviewed native pair matches CoolProp's generate_update_pair table.""" + + first, second = 1.25, 2.5 + oracle_pair, oracle_first, oracle_second = CP.generate_update_pair( + CP.get_parameter_index(name1), + first, + CP.get_parameter_index(name2), + second, + ) + native_pair, swap = encomp_coolprop._native().resolve_input_pair(name1, name2) + native_values = (second, first) if swap else (first, second) + assert native_pair == int(oracle_pair) + assert native_values == (oracle_first, oracle_second) + + +@pytest.mark.parametrize( + "name", + ( + "Water", + "IF97::Water", + "HEOS::CO2[0.5]&O2[0.5]", + "HEOS::CO2[0]&O2[1]", + "INCOMP::MEG[0.5]", + "INCOMP::EG-20%", + ), +) +def test_native_fraction_parser_matches_python_oracle(name: str) -> None: + """The native fallback parser matches CoolProp's documented Python parser.""" + + oracle_backend, fraction_spec = CP.extract_backend(name) + oracle_fluids, oracle_fractions = CP.extract_fractions(fraction_spec) + expected = ( + "HEOS" if oracle_backend == "?" else oracle_backend, + "&".join(oracle_fluids), + list(oracle_fractions) or None, + ) + assert encomp_coolprop._native().resolve_fluid_name(name) == expected + + +def test_native_python_oracle_parity_input_pairs() -> None: + """Parity for non-PT pairs where native resolves order and the oracle uses PropsSI.""" assert encomp_coolprop.self_check() p = np.geomspace(2e5, 200e5, 25) @@ -123,45 +208,36 @@ def test_rust_python_parity_input_pairs() -> None: for second, unit, vals in [("HMASS", "J/kg", h), ("SMASS", "J/kg/K", s)]: df = pl.DataFrame({"P": p, second: vals}) - - def density(mode: str, second: str = second, unit: str = unit, df: pl.DataFrame = df) -> np.ndarray: - clear_expr_evaluation_cache() - second_point: dict[str, Any] - if mode == "rust": - second_point = {second: Q(pl.col(second), unit)} - fluid = Fluid("HEOS::Water", P=Q(pl.col("P"), "Pa"), **second_point) - return df.select(fluid.D.m.alias("d"))["d"].to_numpy().astype(float) - second_point = {second: Q(df[second].to_numpy(), unit)} - fluid = Fluid("HEOS::Water", P=Q(df["P"].to_numpy(), "Pa"), **second_point) - return np.asarray(fluid.D.m, dtype=float) - - py = density("python") - ru = density("rust") - finite_py, finite_ru = np.isfinite(py), np.isfinite(ru) - assert np.array_equal(finite_py, finite_ru), ( - f"P,{second}: finite/NaN masks differ ({int((finite_py != finite_ru).sum())} points)" + clear_expr_evaluation_cache() + second_point: dict[str, Any] = {second: Q(pl.col(second), unit)} + fluid = Fluid("HEOS::Water", P=Q(pl.col("P"), "Pa"), **second_point) + native = df.select(fluid.D.m.alias("d"))["d"].to_numpy().astype(float) + oracle = np.asarray( + CP.PropsSI("DMASS", "P", df["P"].to_numpy(), second, df[second].to_numpy(), "HEOS::Water"), + dtype=float, + ) + finite_oracle, finite_native = np.isfinite(oracle), np.isfinite(native) + assert np.array_equal(finite_oracle, finite_native), ( + f"P,{second}: finite/NaN masks differ ({int((finite_oracle != finite_native).sum())} points)" ) - both = finite_py & finite_ru + both = finite_oracle & finite_native assert both.sum() >= 15, f"P,{second}: only {int(both.sum())} finite points" - rel = np.abs(ru[both] - py[both]) / np.maximum(np.abs(py[both]), 1e-9) + rel = np.abs(native[both] - oracle[both]) / np.maximum(np.abs(oracle[both]), 1e-9) assert rel.max() <= RTOL, f"P,{second}: max rel {rel.max():.2e}" -def test_eager_rust_non_pt_parity(monkeypatch: pytest.MonkeyPatch) -> None: - """The large-eager numpy path (_rust_eager, n >= EAGER_PLUGIN_MIN_SIZE) with a NON-PT - pair must match the Python PropsSI path -- the eager-plugin parity test is PT-only, so - this covers the eager canonical-pair/column-order resolution for P,H.""" +def test_eager_native_non_pt_parity() -> None: + """The eager plugin's non-PT canonical order matches the oracle PropsSI path.""" assert encomp_coolprop.self_check() - n = 1200 # >= EAGER_PLUGIN_MIN_SIZE (1000) -> eager rust plugin + n = 1200 p = np.geomspace(2e5, 200e5, n) t = np.linspace(300.0, 620.0, n) h = CP.PropsSI("HMASS", "P", p, "T", t, "HEOS::Water") eager = np.asarray(Fluid("HEOS::Water", P=Q(p, "Pa"), H=Q(h, "J/kg")).D.m, dtype=float) - monkeypatch.setattr("encomp.fluids.EAGER_PLUGIN_MIN_SIZE", 10**18) # force Python PropsSI - reference = np.asarray(Fluid("HEOS::Water", P=Q(p, "Pa"), H=Q(h, "J/kg")).D.m, dtype=float) + reference = np.asarray(CP.PropsSI("DMASS", "P", p, "HMASS", h, "HEOS::Water"), dtype=float) assert np.array_equal(np.isfinite(eager), np.isfinite(reference)) both = np.isfinite(eager) & np.isfinite(reference) diff --git a/notebooks/README.md b/notebooks/README.md index 6ff7f25..04e9688 100644 --- a/notebooks/README.md +++ b/notebooks/README.md @@ -5,5 +5,5 @@ CI examples. The published tutorial notebook lives in `docs/notebooks/getting-started.ipynb`. The notebooks here may require development dependencies such as `hvplot` and -`altair`, plus the local CoolProp plugin build described in +`altair`, plus the local native CoolProp build described in `encomp/coolprop/README.md`. diff --git a/pyproject.toml b/pyproject.toml index fcd06f4..a8ce097 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -48,16 +48,6 @@ dependencies = [ "polars>=1.42", "sympy>=1.13", "uncertainties>=3.2", - # encomp uses TWO copies of CoolProp: the Python `coolprop` package (scalar / - # small-eager PropsSI path) and the bundled libCoolProp the Rust plugin dlopens - # (lazy / large-eager path). The only integer crossing from Python into the bundled - # lib is the `input_pairs` enum index (via generate_update_pair); that enum is stable - # within a CoolProp major, so allow any 8.x but cap the major (a 9.x could renumber - # it -> silently wrong results). The bundled lib is built at the floor 8.0.0 - # (scripts/build_libcoolprop.py); exact value-parity between the two paths holds when - # the installed coolprop equals that floor (a newer 8.x may differ slightly for a - # fluid whose reference EOS was updated). - "coolprop>=8.0.0,<9", ] [dependency-groups] @@ -89,6 +79,9 @@ dev = [ "microsoft-python-type-stubs", "numpydoc>=1.10.0", ] +# Test-only numerical oracle. Runtime and installed-wheel tests deliberately omit +# this group: encomp evaluates exclusively through its bundled native library. +oracle = ["coolprop==8.0.0"] [project.urls] Homepage = "https://github.com/wlaur/encomp" @@ -112,6 +105,7 @@ include = [ # the bundled CoolProp lib: libCoolProp.{dylib,so} on macOS/Linux, but MSVC names # the Windows DLL CoolProp.dll (no "lib" prefix) -- match both spellings. { path = "encomp/coolprop/*CoolProp.*", format = "wheel" }, + { path = "encomp/coolprop/LICENSE.CoolProp", format = "wheel" }, ] # Per-platform binary wheels (the plugin + bundled libCoolProp are native). abi3 -> @@ -204,6 +198,7 @@ target-version = "py313" max-complexity = 10 [tool.pytest.ini_options] +testpaths = ["encomp/tests"] filterwarnings = [ "ignore::pint.UnitStrippedWarning", ] diff --git a/rust/src/coolprop.rs b/rust/src/coolprop.rs index 6113719..a3facb3 100644 --- a/rust/src/coolprop.rs +++ b/rust/src/coolprop.rs @@ -32,7 +32,8 @@ use libloading::Library; use std::collections::HashSet; -use std::ffi::{CString, c_char, c_double, c_long}; +use std::ffi::{CString, c_char, c_double, c_int, c_long}; +use std::sync::atomic::{AtomicU64, Ordering}; use std::sync::{Mutex, MutexGuard, RwLock, RwLockReadGuard, RwLockWriteGuard}; const BUFLEN: c_long = 1024; @@ -53,7 +54,11 @@ type Batch1Fn = unsafe extern "C" fn( ); type SetFracFn = unsafe extern "C" fn(c_long, *const c_double, c_long, *mut c_long, *mut c_char, c_long); type SpecPhaseFn = unsafe extern "C" fn(c_long, *const c_char, *mut c_long, *mut c_char, c_long); +type UpdateFn = unsafe extern "C" fn(c_long, c_long, c_double, c_double, *mut c_long, *mut c_char, c_long); +type KeyedOutputFn = unsafe extern "C" fn(c_long, c_long, *mut c_long, *mut c_char, c_long) -> c_double; type IndexFn = unsafe extern "C" fn(*const c_char) -> c_long; +type StringFn = unsafe extern "C" fn(*const c_char, *mut c_char, c_int) -> c_long; +type ExtractBackendFn = unsafe extern "C" fn(*const c_char, *mut c_char, c_long, *mut c_char, c_long) -> c_int; type HAPropsSIFn = unsafe extern "C" fn( *const c_char, *const c_char, @@ -100,8 +105,13 @@ pub struct CoolProp { batch1: Batch1Fn, set_fractions: SetFracFn, specify_phase: SpecPhaseFn, + update: UpdateFn, + keyed_output: KeyedOutputFn, get_param_index: IndexFn, get_input_pair_index: IndexFn, + get_global_param_string: StringFn, + get_parameter_information_string: StringFn, + extract_backend: ExtractBackendFn, ha_props_si: HAPropsSIFn, /// NARROW lock: guards ONLY handle create/destroy (the shared handle table), /// never the evaluation path -- so concurrent flashing stays parallel. @@ -109,6 +119,8 @@ pub struct CoolProp { /// keys ("F\0\0", or "HA") already warmed up. Read-locked on the hot /// path (concurrent, no contention once warm); write-locked only for a first-time warmup. warmed: RwLock>, + handles_created: AtomicU64, + handles_freed: AtomicU64, } impl CoolProp { @@ -128,15 +140,34 @@ impl CoolProp { // yield plain `extern "C"` fn pointers that do not borrow `lib`, so `lib` // can be moved into `_lib` afterwards and the pointers stay valid for its // lifetime. - let (factory, free, batch1, set_fractions, specify_phase, get_param_index, get_input_pair_index, ha_props_si) = unsafe { + let ( + factory, + free, + batch1, + set_fractions, + specify_phase, + update, + keyed_output, + get_param_index, + get_input_pair_index, + get_global_param_string, + get_parameter_information_string, + extract_backend, + ha_props_si, + ) = unsafe { ( sym!(FactoryFn, b"AbstractState_factory\0"), sym!(FreeFn, b"AbstractState_free\0"), sym!(Batch1Fn, b"AbstractState_update_and_1_out\0"), sym!(SetFracFn, b"AbstractState_set_fractions\0"), sym!(SpecPhaseFn, b"AbstractState_specify_phase\0"), + sym!(UpdateFn, b"AbstractState_update\0"), + sym!(KeyedOutputFn, b"AbstractState_keyed_output\0"), sym!(IndexFn, b"get_param_index\0"), sym!(IndexFn, b"get_input_pair_index\0"), + sym!(StringFn, b"get_global_param_string\0"), + sym!(StringFn, b"get_parameter_information_string\0"), + sym!(ExtractBackendFn, b"C_extract_backend\0"), sym!(HAPropsSIFn, b"HAPropsSI\0"), ) }; @@ -147,11 +178,18 @@ impl CoolProp { batch1, set_fractions, specify_phase, + update, + keyed_output, get_param_index, get_input_pair_index, + get_global_param_string, + get_parameter_information_string, + extract_backend, ha_props_si, handle_lock: Mutex::new(()), warmed: RwLock::new(HashSet::new()), + handles_created: AtomicU64::new(0), + handles_freed: AtomicU64::new(0), }) } @@ -180,6 +218,62 @@ impl CoolProp { Ok(index) } + /// Read a process-global metadata string such as the bundled library version. + /// Successful calls use only caller-owned buffers. Errors deliberately return a + /// narrow Rust message rather than consulting CoolProp's process-global error + /// outbox, whose text is not safe to associate with one concurrent caller. + pub fn global_param_string(&self, name: &str) -> Result { + let name = cstr(name)?; + let mut out = [0 as c_char; BUFLEN as usize]; + // SAFETY: `name` is NUL terminated and alive for the call; `out` is a writable + // BUFLEN-byte caller-owned buffer, and BUFLEN fits the C `int` parameter. + let ok = unsafe { (self.get_global_param_string)(name.as_ptr(), out.as_mut_ptr(), BUFLEN as c_int) }; + if ok != 1 { + return Err(CpError("CoolProp global metadata lookup failed".into())); + } + Ok(read_msg(&out)) + } + + /// Return one parameter-information field (currently used for SI units). + pub fn parameter_information(&self, name: &str, field: &str) -> Result { + let name = cstr(name)?; + let field = cstr(field)?; + let mut out = [0 as c_char; BUFLEN as usize]; + let field_bytes = field.as_bytes_with_nul(); + if field_bytes.len() > out.len() { + return Err(CpError("CoolProp parameter-information selector is too long".into())); + } + for (dest, source) in out.iter_mut().zip(field_bytes) { + *dest = *source as c_char; + } + // CoolProp uses this buffer as both input (the information selector) and + // output (the resulting string). Both pointers remain valid for the call. + // SAFETY: `name` is NUL terminated; `out` is initialized with a NUL-terminated + // selector and provides BUFLEN writable bytes. + let ok = unsafe { (self.get_parameter_information_string)(name.as_ptr(), out.as_mut_ptr(), BUFLEN as c_int) }; + if ok != 1 { + return Err(CpError(format!( + "unknown parameter {name:?} or information field {field:?}" + ))); + } + Ok(read_msg(&out)) + } + + /// Split an optional ``BACKEND::`` prefix using CoolProp's own parser. + pub fn extract_backend(&self, name: &str) -> Result<(String, String), CpError> { + let name = cstr(name)?; + let mut backend = [0 as c_char; BUFLEN as usize]; + let mut fluid = [0 as c_char; BUFLEN as usize]; + // SAFETY: `name` is NUL terminated; both output buffers are caller-owned, + // writable, and their exact lengths are supplied. + let code = + unsafe { (self.extract_backend)(name.as_ptr(), backend.as_mut_ptr(), BUFLEN, fluid.as_mut_ptr(), BUFLEN) }; + if code != 0 { + return Err(CpError("CoolProp backend/fluid name exceeds the native buffer".into())); + } + Ok((read_msg(&backend), read_msg(&fluid))) + } + /// Warm up ONE (backend, fluid) config, lazily. The FIRST call (any config) also runs /// CoolProp's process-global lazy inits (fluid library, index maps, flash machinery); /// every distinct config additionally pays one warmup flash so its backend-specific init @@ -260,7 +354,7 @@ impl CoolProp { // state we touch (8.0 per-thread backends), so the loop is lock-free. let val = unsafe { (self.ha_props_si)(o.as_ptr(), a.as_ptr(), v1[i], b.as_ptr(), v2[i], c.as_ptr(), v3[i]) }; - out[i] = if val.is_finite() { val } else { f64::NAN }; // CoolProp returns _HUGE on error + out[i] = if valid_result(val) { val } else { f64::NAN }; // CoolProp returns _HUGE on error } Ok(()) } @@ -281,8 +375,22 @@ impl CoolProp { if err != 0 { return Err(CpError(format!("factory({backend},{fluid}): {}", read_msg(&msg)))); } + self.handles_created.fetch_add(1, Ordering::Relaxed); Ok(State { cp: self, handle }) } + + pub fn handle_counts(&self) -> (u64, u64) { + ( + self.handles_created.load(Ordering::Relaxed), + self.handles_freed.load(Ordering::Relaxed), + ) + } +} + +/// CoolProp's C wrapper returns `_HUGE` (currently 1e300) for several failures. +/// No physical property encomp supports approaches this sentinel. +fn valid_result(value: f64) -> bool { + value.is_finite() && value.abs() < 1e290 } /// An owned AbstractState handle. RAII-frees on drop (under the narrow lock). @@ -380,6 +488,52 @@ impl State<'_> { } Ok(()) } + + /// Fast scalar evaluation through the same batched C entry point as the array + /// plugin. The successful path is one FFI call. If that call leaves its output + /// untouched, retry through the individually error-reporting update/keyed-output + /// functions so Python receives the existing narrow `ValueError` behavior without + /// putting the global PropsSI error outbox back on the hot path. + pub fn update_and_1_out_scalar( + &mut self, + input_pair: i64, + value1: f64, + value2: f64, + out_key: c_long, + ) -> Result { + let mut out = [f64::NAN]; + self.update_and_1_out(input_pair, &[value1], &[value2], out_key, &mut out)?; + if valid_result(out[0]) { + return Ok(out[0]); + } + + let pair = c_long::try_from(input_pair) + .map_err(|_| CpError(format!("input pair {input_pair} exceeds the C API integer range")))?; + let mut err: c_long = 0; + let mut msg = [0 as c_char; BUFLEN as usize]; + // SAFETY: `self.handle` is live and thread-owned; scalar values are passed by + // value; error buffers are caller-owned and valid for BUFLEN bytes. + unsafe { + (self.cp.update)(self.handle, pair, value1, value2, &mut err, msg.as_mut_ptr(), BUFLEN); + } + if err != 0 { + return Err(CpError(read_msg(&msg))); + } + + err = 0; + msg.fill(0); + // SAFETY: `self.handle` remains live after the successful update; the output + // key came from this same library; error buffers are caller-owned. + let value = unsafe { (self.cp.keyed_output)(self.handle, out_key, &mut err, msg.as_mut_ptr(), BUFLEN) }; + if err != 0 { + return Err(CpError(read_msg(&msg))); + } + if valid_result(value) { + Ok(value) + } else { + Err(CpError("CoolProp returned a non-finite result".into())) + } + } } impl Drop for State<'_> { @@ -390,6 +544,7 @@ impl Drop for State<'_> { // SAFETY: frees a handle we own exactly once (Drop runs once); the // handle-table write is serialised by `handle_lock` held here. unsafe { (self.cp.free)(self.handle, &mut err, msg.as_mut_ptr(), BUFLEN) }; + self.cp.handles_freed.fetch_add(1, Ordering::Relaxed); } } diff --git a/rust/src/lib.rs b/rust/src/lib.rs index 1330531..bf9a0cb 100644 --- a/rust/src/lib.rs +++ b/rust/src/lib.rs @@ -8,17 +8,78 @@ mod coolprop; -use coolprop::{CoolProp, CpError}; +use coolprop::{CoolProp, CpError, State}; use polars::prelude::*; +use pyo3::exceptions::{PyRuntimeError, PyValueError}; +use pyo3::prelude::*; use pyo3_polars::derive::polars_expr; use serde::Deserialize; -use std::sync::{Mutex, OnceLock}; +use std::cell::RefCell; +use std::collections::VecDeque; +use std::sync::{Arc, Mutex, OnceLock}; // One loaded libCoolProp per process (lib_path is constant for a session). static CP: OnceLock = OnceLock::new(); // serialises the one-time load + warmup (see `coolprop`); untouched on the hot path static INIT: Mutex<()> = Mutex::new(()); +/// Prepared scalar configurations contain only owned Rust values. Python caches +/// their integer IDs, so a hot scalar call crosses PyO3 with four integers and two +/// floats rather than repeatedly allocating backend/fluid/property strings. +#[derive(Debug, PartialEq)] +struct FluidConfig { + backend: String, + fluid: String, + fractions: Option>, + phase: Option, +} + +#[derive(Debug, PartialEq)] +struct HumidAirConfig { + output: String, + name1: String, + name2: String, + name3: String, +} + +static FLUID_CONFIGS: Mutex>> = Mutex::new(Vec::new()); +static HUMID_AIR_CONFIGS: Mutex>> = Mutex::new(Vec::new()); + +const SCALAR_STATE_CACHE_CAPACITY: usize = 16; + +struct ScalarStateCache { + /// Least recently used at the front, most recently used at the back. + entries: VecDeque<(usize, State<'static>)>, + hits: u64, + misses: u64, + evictions: u64, +} + +impl ScalarStateCache { + const fn new() -> Self { + Self { + entries: VecDeque::new(), + hits: 0, + misses: 0, + evictions: 0, + } + } + + fn clear(&mut self) { + // Dropping each State frees its native handle under CoolProp's narrow + // handle-table lock. This also makes destruction directly testable. + self.entries.clear(); + self.hits = 0; + self.misses = 0; + self.evictions = 0; + } +} + +thread_local! { + /// Mutable AbstractState handles never leave their owning OS thread. + static SCALAR_STATES: RefCell = const { RefCell::new(ScalarStateCache::new()) }; +} + fn perr(e: CpError) -> PolarsError { PolarsError::ComputeError(e.0.into()) } @@ -27,6 +88,14 @@ fn compute_error(msg: impl Into) -> PolarsError { PolarsError::ComputeError(msg.into().into()) } +fn cp_error(e: CpError) -> PyErr { + PyValueError::new_err(e.0) +} + +fn runtime_error(message: impl Into) -> PyErr { + PyRuntimeError::new_err(message.into()) +} + fn validate_input_count(kind: &str, got: usize, expected: usize) -> PolarsResult<()> { if got != expected { return Err(compute_error(format!("{kind}: expected {expected} inputs, got {got}"))); @@ -129,7 +198,7 @@ fn to_f64(s: &Series) -> PolarsResult> { Ok(s.f64()?.iter().map(|o| o.unwrap_or(f64::NAN)).collect()) } -fn coolprop(lib_path: &str) -> PolarsResult<&'static CoolProp> { +fn load_coolprop(lib_path: &str) -> Result<&'static CoolProp, CpError> { if let Some(c) = CP.get() { return Ok(c); // fast path: already initialised, no lock } @@ -139,14 +208,415 @@ fn coolprop(lib_path: &str) -> PolarsResult<&'static CoolProp> { // (CoolProp::ensure_warmed_*), so an unused backend is never warmed. let _guard = INIT .lock() - .map_err(|_| compute_error("CoolProp init lock was poisoned"))?; + .map_err(|_| CpError("CoolProp init lock was poisoned".into()))?; if let Some(c) = CP.get() { return Ok(c); } - let c = CoolProp::load(lib_path).map_err(perr)?; + let c = CoolProp::load(lib_path)?; CP.set(c) - .map_err(|_| compute_error("CoolProp was initialized concurrently"))?; - CP.get().ok_or_else(|| compute_error("CoolProp failed to initialize")) + .map_err(|_| CpError("CoolProp was initialized concurrently".into()))?; + CP.get().ok_or_else(|| CpError("CoolProp failed to initialize".into())) +} + +fn coolprop(lib_path: &str) -> PolarsResult<&'static CoolProp> { + load_coolprop(lib_path).map_err(perr) +} + +fn initialized_coolprop() -> Result<&'static CoolProp, CpError> { + CP.get() + .ok_or_else(|| CpError("native CoolProp is not initialized; call initialize(lib_path) first".into())) +} + +fn get_fluid_config(config_id: usize) -> Result, CpError> { + FLUID_CONFIGS + .lock() + .map_err(|_| CpError("fluid configuration registry lock was poisoned".into()))? + .get(config_id) + .cloned() + .ok_or_else(|| CpError(format!("unknown fluid configuration id {config_id}"))) +} + +fn get_humid_air_config(config_id: usize) -> Result, CpError> { + HUMID_AIR_CONFIGS + .lock() + .map_err(|_| CpError("humid-air configuration registry lock was poisoned".into()))? + .get(config_id) + .cloned() + .ok_or_else(|| CpError(format!("unknown humid-air configuration id {config_id}"))) +} + +fn configured_state(cp: &'static CoolProp, config: &FluidConfig) -> Result, CpError> { + cp.ensure_warmed_fluid(&config.backend, &config.fluid)?; + let mut state = cp.state(&config.backend, &config.fluid)?; + if let Some(fractions) = &config.fractions { + state.set_fractions(fractions)?; + } + if let Some(phase) = &config.phase { + state.specify_phase(phase)?; + } + Ok(state) +} + +fn evaluate_scalar_fluid( + cp: &'static CoolProp, + config_id: usize, + input_pair: i64, + value1: f64, + value2: f64, + output: i64, +) -> Result { + SCALAR_STATES.with(|cache_cell| { + let mut cache = cache_cell.borrow_mut(); + let position = cache.entries.iter().position(|(id, _)| *id == config_id); + let mut entry = if let Some(position) = position { + cache.hits += 1; + cache + .entries + .remove(position) + .ok_or_else(|| CpError("scalar state cache entry disappeared".into()))? + } else { + cache.misses += 1; + let config = get_fluid_config(config_id)?; + let state = configured_state(cp, &config)?; + if cache.entries.len() == SCALAR_STATE_CACHE_CAPACITY { + cache.entries.pop_front(); + cache.evictions += 1; + } + (config_id, state) + }; + let output = std::ffi::c_long::try_from(output) + .map_err(|_| CpError(format!("output parameter {output} exceeds the C API integer range")))?; + let value = entry.1.update_and_1_out_scalar(input_pair, value1, value2, output)?; + cache.entries.push_back(entry); + Ok(value) + }) +} + +#[derive(Clone, Copy)] +struct ResolvedInputPair { + first: std::ffi::c_long, + second: std::ffi::c_long, + pair: i64, +} + +/// Mirrors CoolProp's reviewed `generate_update_pair` table, but asks the loaded +/// library for every parameter and pair number at runtime. Enum integers therefore +/// never cross a version boundary or live in encomp source. +const INPUT_PAIR_NAMES: &[(&str, &str, &str)] = &[ + ("Q", "T", "QT_INPUTS"), + ("Qmass", "T", "QmassT_INPUTS"), + ("P", "Q", "PQ_INPUTS"), + ("P", "Qmass", "PQmass_INPUTS"), + ("P", "T", "PT_INPUTS"), + ("Dmolar", "T", "DmolarT_INPUTS"), + ("Dmass", "T", "DmassT_INPUTS"), + ("Hmolar", "T", "HmolarT_INPUTS"), + ("Hmass", "T", "HmassT_INPUTS"), + ("Smolar", "T", "SmolarT_INPUTS"), + ("Smass", "T", "SmassT_INPUTS"), + ("T", "Umolar", "TUmolar_INPUTS"), + ("T", "Umass", "TUmass_INPUTS"), + ("Dmass", "Hmass", "DmassHmass_INPUTS"), + ("Dmolar", "Hmolar", "DmolarHmolar_INPUTS"), + ("Dmass", "Smass", "DmassSmass_INPUTS"), + ("Dmolar", "Smolar", "DmolarSmolar_INPUTS"), + ("Dmass", "Umass", "DmassUmass_INPUTS"), + ("Dmolar", "Umolar", "DmolarUmolar_INPUTS"), + ("Dmass", "P", "DmassP_INPUTS"), + ("Dmolar", "P", "DmolarP_INPUTS"), + ("Dmass", "Q", "DmassQ_INPUTS"), + ("Dmass", "Qmass", "DmassQmass_INPUTS"), + ("Dmolar", "Q", "DmolarQ_INPUTS"), + ("Dmolar", "Qmass", "DmolarQmass_INPUTS"), + ("Hmass", "P", "HmassP_INPUTS"), + ("Hmolar", "P", "HmolarP_INPUTS"), + ("P", "Smass", "PSmass_INPUTS"), + ("P", "Smolar", "PSmolar_INPUTS"), + ("P", "Umass", "PUmass_INPUTS"), + ("P", "Umolar", "PUmolar_INPUTS"), + ("Hmass", "Smass", "HmassSmass_INPUTS"), + ("Hmolar", "Smolar", "HmolarSmolar_INPUTS"), + ("Smass", "Umass", "SmassUmass_INPUTS"), + ("Smolar", "Umolar", "SmolarUmolar_INPUTS"), +]; + +static RESOLVED_INPUT_PAIRS: OnceLock, String>> = OnceLock::new(); + +fn resolved_input_pairs(cp: &CoolProp) -> Result<&[ResolvedInputPair], CpError> { + RESOLVED_INPUT_PAIRS + .get_or_init(|| { + INPUT_PAIR_NAMES + .iter() + .map(|(first, second, pair)| { + Ok(ResolvedInputPair { + first: cp.param_index(first).map_err(|e| e.0)?, + second: cp.param_index(second).map_err(|e| e.0)?, + pair: cp.input_pair_index(pair).map_err(|e| e.0)?, + }) + }) + .collect() + }) + .as_deref() + .map_err(|message| CpError(message.clone())) +} + +fn resolve_pair_native(cp: &CoolProp, name1: &str, name2: &str) -> Result<(i64, bool), CpError> { + let key1 = cp.param_index(name1)?; + let key2 = cp.param_index(name2)?; + for resolved in resolved_input_pairs(cp)? { + if key1 == resolved.first && key2 == resolved.second { + return Ok((resolved.pair, false)); + } + if key1 == resolved.second && key2 == resolved.first { + return Ok((resolved.pair, true)); + } + } + Err(CpError(format!( + "unsupported CoolProp input pair: {name1:?}, {name2:?}" + ))) +} + +/// Fallback implementation of CoolProp's documented fraction grammar. Upstream +/// exposes backend splitting through C but not `extract_fractions`, so this parser +/// is parity-tested against the Python binding in the oracle CI job. +fn extract_fractions(fluid: &str) -> Result<(String, Option>), CpError> { + if fluid.contains('[') && fluid.contains(']') { + let parts: Vec<&str> = fluid.split('&').collect(); + let mut names = Vec::with_capacity(parts.len()); + let mut fractions = Vec::with_capacity(parts.len()); + for part in &parts { + let without_close = part + .strip_suffix(']') + .ok_or_else(|| CpError(format!("Fluid entry [{part}] must end with ']' character")))?; + let mut pieces = without_close.split('['); + let name = pieces.next().unwrap_or_default(); + let fraction_text = pieces + .next() + .ok_or_else(|| CpError(format!("Could not break [{without_close}] into name/fraction")))?; + if name.is_empty() || pieces.next().is_some() { + return Err(CpError(format!("Could not break [{without_close}] into name/fraction"))); + } + let fraction = fraction_text + .parse::() + .map_err(|_| CpError(format!("fraction [{fraction_text}] was not converted fully")))?; + if !fraction.is_finite() || !(0.0..=1.0).contains(&fraction) { + return Err(CpError(format!( + "fraction [{fraction_text}] was not converted to a value between 0 and 1 inclusive" + ))); + } + // CoolProp removes zero-fraction components from multi-fluid names, but + // retains a zero single-fluid INCOMP concentration. + if fraction > 10.0 * f64::EPSILON || parts.len() == 1 { + names.push(name); + fractions.push(fraction); + } + } + return Ok((names.join("&"), Some(fractions))); + } + + if fluid.contains('-') && fluid.contains('%') { + let parts: Vec<&str> = fluid.split('-').collect(); + if parts.len() != 2 || parts[0].is_empty() { + return Err(CpError(format!( + "format of incompressible solution {fluid:?} is invalid; expected EG-20%" + ))); + } + let fraction_text = parts[1] + .strip_suffix('%') + .ok_or_else(|| CpError(format!("invalid incompressible concentration {fluid:?}")))?; + let fraction = fraction_text + .parse::() + .map_err(|_| CpError(format!("invalid incompressible concentration {fluid:?}")))? + * 0.01; + if !fraction.is_finite() { + return Err(CpError(format!("invalid incompressible concentration {fluid:?}"))); + } + return Ok((parts[0].to_string(), Some(vec![fraction]))); + } + + Ok((fluid.to_string(), None)) +} + +#[pyfunction(name = "initialize")] +fn py_initialize(lib_path: &str) -> PyResult<()> { + load_coolprop(lib_path).map(|_| ()).map_err(cp_error) +} + +#[pyfunction(name = "parameter_index")] +fn py_parameter_index(name: &str) -> PyResult { + let index = initialized_coolprop() + .and_then(|cp| cp.param_index(name)) + .map_err(cp_error)?; + #[allow(clippy::useless_conversion)] + Ok(i64::from(index)) +} + +#[pyfunction(name = "parameter_information")] +fn py_parameter_information(name: &str, field: &str) -> PyResult { + initialized_coolprop() + .and_then(|cp| cp.parameter_information(name, field)) + .map_err(cp_error) +} + +#[pyfunction(name = "resolve_input_pair")] +fn py_resolve_input_pair(name1: &str, name2: &str) -> PyResult<(i64, bool)> { + resolve_pair_native(initialized_coolprop().map_err(cp_error)?, name1, name2).map_err(cp_error) +} + +#[pyfunction(name = "resolve_fluid_name")] +fn py_resolve_fluid_name(name: &str) -> PyResult<(String, String, Option>)> { + let cp = initialized_coolprop().map_err(cp_error)?; + let (backend, fluid) = cp.extract_backend(name).map_err(cp_error)?; + let (fluid, fractions) = extract_fractions(&fluid).map_err(cp_error)?; + let backend = if backend == "?" { "HEOS".to_string() } else { backend }; + Ok((backend, fluid, fractions)) +} + +#[pyfunction(name = "validate_fluid", signature = (backend, fluid, fractions=None))] +fn py_validate_fluid(py: Python<'_>, backend: String, fluid: String, fractions: Option>) -> PyResult<()> { + let cp = initialized_coolprop().map_err(cp_error)?; + py.detach(move || { + cp.ensure_warmed_fluid(&backend, &fluid)?; + let mut state = cp.state(&backend, &fluid)?; + if let Some(fractions) = fractions { + state.set_fractions(&fractions)?; + } + Ok(()) + }) + .map_err(cp_error) +} + +#[pyfunction(name = "prepare_fluid", signature = (backend, fluid, fractions=None, phase=None))] +fn py_prepare_fluid( + backend: String, + fluid: String, + fractions: Option>, + phase: Option, +) -> PyResult { + let candidate = FluidConfig { + backend, + fluid, + fractions, + phase, + }; + let mut configs = FLUID_CONFIGS + .lock() + .map_err(|_| runtime_error("fluid configuration registry lock was poisoned"))?; + if let Some(id) = configs.iter().position(|config| config.as_ref() == &candidate) { + return Ok(id); + } + let id = configs.len(); + configs.push(Arc::new(candidate)); + Ok(id) +} + +#[pyfunction(name = "prepare_humid_air")] +fn py_prepare_humid_air(output: String, name1: String, name2: String, name3: String) -> PyResult { + let candidate = HumidAirConfig { + output, + name1, + name2, + name3, + }; + let mut configs = HUMID_AIR_CONFIGS + .lock() + .map_err(|_| runtime_error("humid-air configuration registry lock was poisoned"))?; + if let Some(id) = configs.iter().position(|config| config.as_ref() == &candidate) { + return Ok(id); + } + let id = configs.len(); + configs.push(Arc::new(candidate)); + Ok(id) +} + +#[pyfunction(name = "fluid_scalar")] +fn py_fluid_scalar( + py: Python<'_>, + config_id: usize, + input_pair: i64, + value1: f64, + value2: f64, + output: i64, +) -> PyResult { + let cp = initialized_coolprop().map_err(cp_error)?; + py.detach(move || evaluate_scalar_fluid(cp, config_id, input_pair, value1, value2, output)) + .map_err(cp_error) +} + +#[pyfunction(name = "humid_air_scalar")] +fn py_humid_air_scalar(py: Python<'_>, config_id: usize, value1: f64, value2: f64, value3: f64) -> PyResult { + let cp = initialized_coolprop().map_err(cp_error)?; + let config = get_humid_air_config(config_id).map_err(cp_error)?; + py.detach(move || { + cp.ensure_warmed_humid_air()?; + let mut out = [f64::NAN]; + cp.ha_props_si_batch( + &config.output, + &config.name1, + &[value1], + &config.name2, + &[value2], + &config.name3, + &[value3], + &mut out, + )?; + Ok(out[0]) + }) + .map_err(cp_error) +} + +#[pyfunction(name = "lib_version")] +fn py_lib_version() -> PyResult { + initialized_coolprop() + .and_then(|cp| cp.global_param_string("version")) + .map_err(cp_error) +} + +#[pyfunction(name = "clear_scalar_cache")] +fn py_clear_scalar_cache() { + SCALAR_STATES.with(|cache| cache.borrow_mut().clear()); +} + +#[pyfunction(name = "scalar_cache_info")] +fn py_scalar_cache_info() -> (usize, usize, u64, u64, u64) { + SCALAR_STATES.with(|cache| { + let cache = cache.borrow(); + ( + cache.entries.len(), + SCALAR_STATE_CACHE_CAPACITY, + cache.hits, + cache.misses, + cache.evictions, + ) + }) +} + +#[pyfunction(name = "handle_counts")] +fn py_handle_counts() -> PyResult<(u64, u64)> { + initialized_coolprop().map(CoolProp::handle_counts).map_err(cp_error) +} + +/// The PyO3 API and Polars plugin C ABI intentionally coexist in this cdylib. +/// No Python object is stored in native global state; scalar handles are bounded +/// and thread-local, so the module is safe to import without enabling the GIL. +#[pymodule] +#[pyo3(gil_used = false)] +fn _internal(module: &Bound<'_, PyModule>) -> PyResult<()> { + module.add_function(wrap_pyfunction!(py_initialize, module)?)?; + module.add_function(wrap_pyfunction!(py_parameter_index, module)?)?; + module.add_function(wrap_pyfunction!(py_parameter_information, module)?)?; + module.add_function(wrap_pyfunction!(py_resolve_input_pair, module)?)?; + module.add_function(wrap_pyfunction!(py_resolve_fluid_name, module)?)?; + module.add_function(wrap_pyfunction!(py_validate_fluid, module)?)?; + module.add_function(wrap_pyfunction!(py_prepare_fluid, module)?)?; + module.add_function(wrap_pyfunction!(py_prepare_humid_air, module)?)?; + module.add_function(wrap_pyfunction!(py_fluid_scalar, module)?)?; + module.add_function(wrap_pyfunction!(py_humid_air_scalar, module)?)?; + module.add_function(wrap_pyfunction!(py_lib_version, module)?)?; + module.add_function(wrap_pyfunction!(py_clear_scalar_cache, module)?)?; + module.add_function(wrap_pyfunction!(py_scalar_cache_info, module)?)?; + module.add_function(wrap_pyfunction!(py_handle_counts, module)?)?; + Ok(()) } /// Output dtype preserves the input precision: Float32 only when every non-scalar @@ -343,6 +813,23 @@ mod tests { assert_eq!(ca.get(2), None); assert_eq!(ca.get(3), Some(-2.5)); } + + #[test] + fn fraction_parser_covers_mixtures_and_incompressibles() { + assert_eq!(extract_fractions("Water").unwrap(), ("Water".into(), None)); + assert_eq!( + extract_fractions("CO2[0.5]&O2[0.5]").unwrap(), + ("CO2&O2".into(), Some(vec![0.5, 0.5])) + ); + assert_eq!( + extract_fractions("CO2[0]&O2[1]").unwrap(), + ("O2".into(), Some(vec![1.0])) + ); + assert_eq!(extract_fractions("MEG[0.5]").unwrap(), ("MEG".into(), Some(vec![0.5]))); + assert_eq!(extract_fractions("EG-20%").unwrap(), ("EG".into(), Some(vec![0.2]))); + assert!(extract_fractions("CO2[1.2]&O2[-0.2]").is_err()); + assert!(extract_fractions("CO2[abc]&O2[1]").is_err()); + } } #[derive(Deserialize)] diff --git a/scripts/benchmark_coolprop.py b/scripts/benchmark_coolprop.py new file mode 100644 index 0000000..0f88a20 --- /dev/null +++ b/scripts/benchmark_coolprop.py @@ -0,0 +1,248 @@ +#!/usr/bin/env python +"""Benchmark the direct scalar bridge and existing Polars batch path. + +Run after a release build (``uv run maturin develop --release``). The optional +Python CoolProp oracle is reported only when the ``oracle`` dependency group is +installed; it is never needed by encomp itself. +""" + +from __future__ import annotations + +import argparse +import ctypes +import importlib +import json +import statistics +import time +import timeit +from collections.abc import Callable +from typing import Any + +import numpy as np +import polars as pl + +from encomp import coolprop as cp +from encomp.fluids import Fluid, HumidAir +from encomp.units import Quantity as Q + + +def median_seconds(function: Callable[[], object], *, number: int, repeat: int) -> float: + return statistics.median(timeit.repeat(function, number=number, repeat=repeat)) / number + + +def benchmark(*, array_size: int, repeat: int, quick: bool) -> dict[str, float | int | list[int]]: + native = cp._native() # pyright: ignore[reportPrivateUsage] + pair = native.resolve_input_pair("P", "T")[0] + density = native.parameter_index("DMASS") + if97_config = native.prepare_fluid("IF97", "Water") + heos_config = native.prepare_fluid("HEOS", "CarbonDioxide") + humid_config = native.prepare_humid_air("W", "P", "T", "R") + + native.clear_scalar_cache() + started_first_call = time.perf_counter_ns() + native.fluid_scalar(if97_config, pair, 101325.0, 300.0, density) + first_call_us = (time.perf_counter_ns() - started_first_call) / 1e3 + + library = ctypes.CDLL(cp.lib_path()) + props_si: Any = library.PropsSI + props_si.argtypes = [ + ctypes.c_char_p, + ctypes.c_char_p, + ctypes.c_double, + ctypes.c_char_p, + ctypes.c_double, + ctypes.c_char_p, + ] + props_si.restype = ctypes.c_double + ha_props_si: Any = library.HAPropsSI + ha_props_si.argtypes = [ + ctypes.c_char_p, + ctypes.c_char_p, + ctypes.c_double, + ctypes.c_char_p, + ctypes.c_double, + ctypes.c_char_p, + ctypes.c_double, + ] + ha_props_si.restype = ctypes.c_double + + if97 = Fluid("IF97::Water", P=Q(101325.0, "Pa"), T=Q(300.0, "K")) + heos = Fluid("HEOS::CarbonDioxide", P=Q(5e6, "Pa"), T=Q(300.0, "K")) + humid = HumidAir(P=Q(101325.0, "Pa"), T=Q(300.0, "K"), R=Q(0.5, "")) + _ = (if97.D.m, heos.D.m, humid.W.m) + + scalar_number = 2_000 if quick else 20_000 + fast_number = 10_000 if quick else 100_000 + results: dict[str, float | int | list[int]] = { + "array_size": array_size, + "native_if97_first_call_us": first_call_us, + "native_if97_us": median_seconds( + lambda: native.fluid_scalar(if97_config, pair, 101325.0, 300.0, density), + number=fast_number, + repeat=repeat, + ) + * 1e6, + "public_if97_us": median_seconds(lambda: if97.D.m, number=scalar_number, repeat=repeat) * 1e6, + "native_heos_us": median_seconds( + lambda: native.fluid_scalar(heos_config, pair, 5e6, 300.0, density), + number=scalar_number, + repeat=repeat, + ) + * 1e6, + "public_heos_us": median_seconds(lambda: heos.D.m, number=scalar_number, repeat=repeat) * 1e6, + "native_humid_air_us": median_seconds( + lambda: native.humid_air_scalar(humid_config, 101325.0, 300.0, 0.5), + number=fast_number, + repeat=repeat, + ) + * 1e6, + "public_humid_air_us": median_seconds(lambda: humid.W.m, number=scalar_number, repeat=repeat) * 1e6, + "bundled_cabi_if97_us": median_seconds( + lambda: props_si(b"DMASS", b"P", 101325.0, b"T", 300.0, b"IF97::Water"), + number=fast_number, + repeat=repeat, + ) + * 1e6, + "bundled_cabi_heos_us": median_seconds( + lambda: props_si(b"DMASS", b"P", 5e6, b"T", 300.0, b"HEOS::CarbonDioxide"), + number=scalar_number, + repeat=repeat, + ) + * 1e6, + "bundled_cabi_humid_air_us": median_seconds( + lambda: ha_props_si(b"W", b"P", 101325.0, b"T", 300.0, b"R", 0.5), + number=fast_number, + repeat=repeat, + ) + * 1e6, + } + + temperatures = tuple(float(value) for value in np.linspace(290.0, 310.0, 127)) + varying_index = 0 + + def varying_heos() -> float: + nonlocal varying_index + value = native.fluid_scalar(heos_config, pair, 5e6, temperatures[varying_index % len(temperatures)], density) + varying_index += 1 + return value + + results["native_heos_varying_state_us"] = ( + median_seconds( + varying_heos, + number=scalar_number, + repeat=repeat, + ) + * 1e6 + ) + + def cold_if97() -> float: + native.clear_scalar_cache() + return native.fluid_scalar(if97_config, pair, 101325.0, 300.0, density) + + results["native_if97_cache_miss_us"] = ( + median_seconds( + cold_if97, + number=1, + repeat=max(repeat, 11), + ) + * 1e6 + ) + + eviction_configs = [ + native.prepare_fluid("HEOS", fluid) + for fluid in ( + "Water", + "CarbonDioxide", + "Nitrogen", + "Oxygen", + "Methane", + "Ammonia", + "R134a", + "n-Propane", + "Hydrogen", + "Argon", + "Ethane", + "CarbonMonoxide", + "Helium", + "Neon", + "Krypton", + "Xenon", + "SulfurHexafluoride", + ) + ] + for config in eviction_configs: + native.fluid_scalar(config, pair, 101325.0, 300.0, density) + native.clear_scalar_cache() + + def fill_and_evict() -> None: + native.clear_scalar_cache() + for config in eviction_configs: + native.fluid_scalar(config, pair, 101325.0, 300.0, density) + + results["native_cache_fill_and_eviction_us_per_config"] = ( + median_seconds(fill_and_evict, number=1, repeat=repeat) / len(eviction_configs) * 1e6 + ) + + pressure = np.full(array_size, 5e6) + temperature = np.linspace(290.0, 310.0, array_size) + eager = Fluid("HEOS::CarbonDioxide", P=Q(pressure, "Pa"), T=Q(temperature, "K")) + _ = eager.D.m + array_seconds = median_seconds(lambda: eager.D.m, number=1, repeat=repeat) + results["heos_array_ms"] = array_seconds * 1e3 + results["heos_array_mrows_s"] = array_size / array_seconds / 1e6 + + frame = pl.DataFrame({"P": pressure, "T": temperature}) + expressions = [cp.fluid(prop, "P", "T", name="HEOS::CarbonDioxide") for prop in ("DMASS", "HMASS", "SMASS")] + frame.select(expressions) + multi_seconds = median_seconds(lambda: frame.select(expressions), number=1, repeat=repeat) + results["three_property_polars_ms"] = multi_seconds * 1e3 + results["three_property_polars_mprops_s"] = 3 * array_size / multi_seconds / 1e6 + results["scalar_cache_info"] = list(native.scalar_cache_info()) + + try: + oracle: Any = importlib.import_module("CoolProp.CoolProp") + except ModuleNotFoundError: + pass + else: + results["oracle_if97_us"] = ( + median_seconds( + lambda: oracle.PropsSI("DMASS", "P", 101325.0, "T", 300.0, "IF97::Water"), + number=fast_number, + repeat=repeat, + ) + * 1e6 + ) + results["oracle_heos_us"] = ( + median_seconds( + lambda: oracle.PropsSI("DMASS", "P", 5e6, "T", 300.0, "HEOS::CarbonDioxide"), + number=scalar_number, + repeat=repeat, + ) + * 1e6 + ) + results["oracle_humid_air_us"] = ( + median_seconds( + lambda: oracle.HAPropsSI("W", "P", 101325.0, "T", 300.0, "R", 0.5), + number=fast_number, + repeat=repeat, + ) + * 1e6 + ) + + return results + + +def main() -> None: + parser = argparse.ArgumentParser() + parser.add_argument("--array-size", type=int, default=100_000) + parser.add_argument("--repeat", type=int, default=7) + parser.add_argument("--quick", action="store_true") + args = parser.parse_args() + started = time.perf_counter() + results = benchmark(array_size=args.array_size, repeat=args.repeat, quick=args.quick) + results["elapsed_s"] = time.perf_counter() - started + print(json.dumps(results, indent=2, sort_keys=True)) + + +if __name__ == "__main__": + main() diff --git a/scripts/build_libcoolprop.py b/scripts/build_libcoolprop.py index d6cd6a4..59b6b05 100644 --- a/scripts/build_libcoolprop.py +++ b/scripts/build_libcoolprop.py @@ -10,7 +10,7 @@ from __future__ import annotations import os -import re +import runpy import shutil import subprocess import sys @@ -18,32 +18,9 @@ PROJECT = Path(__file__).resolve().parent.parent # repo root PKG = PROJECT / "encomp" / "coolprop" - - -def _coolprop_version() -> str: - """The CoolProp git tag to build, derived from the LOWER BOUND of the ``coolprop`` - requirement in pyproject.toml (single source of truth). The bundled C++ library is - built at this floor; the Python ``coolprop`` package may be any version the - requirement allows (>= the floor, same major). The two must share the CoolProp major - so the ``input_pairs`` enum index encomp computes on the Python side stays valid for - the bundled library. - - Parsed with a small regex because only one dependency spec is needed here.""" - text = (PROJECT / "pyproject.toml").read_text() - match = re.search(r"""["']coolprop\s*(?:==|>=)\s*([\w.]+)""", text) - if match: - return f"v{match.group(1)}" - raise RuntimeError("no 'coolprop==' or 'coolprop>=' requirement found in pyproject.toml [project.dependencies]") - - -COOLPROP_VERSION = _coolprop_version() # e.g. "v8.0.0" (the floor of the coolprop requirement) - -# A git tag is mutable, and this clone produces the C++ library that ships inside the wheel. -# Pin the commit each supported tag must resolve to, and verify it after cloning. Submodules -# are pinned by the superproject commit, so this covers them too. -COOLPROP_COMMITS = { - "v8.0.0": "ae81610e7d23efc57f9d051c8e70a4d66e87537f", -} +BUILD_INFO = runpy.run_path(str(PKG / "_build_info.py")) +COOLPROP_VERSION = str(BUILD_INFO["BUNDLED_COOLPROP_TAG"]) +COOLPROP_COMMIT = str(BUILD_INFO["BUNDLED_COOLPROP_COMMIT"]) def run(*cmd: str, cwd: Path | None = None) -> None: @@ -52,14 +29,6 @@ def run(*cmd: str, cwd: Path | None = None) -> None: def verify_commit(src: Path) -> None: - expected = COOLPROP_COMMITS.get(COOLPROP_VERSION) - if expected is None: - raise RuntimeError( - f"no pinned commit for CoolProp {COOLPROP_VERSION}. Resolve the tag " - f"(git ls-remote --tags https://github.com/CoolProp/CoolProp.git {COOLPROP_VERSION}) " - "and add the commit it peels to under COOLPROP_COMMITS." - ) - head = subprocess.run( ["git", "rev-parse", "HEAD"], cwd=src, @@ -68,9 +37,9 @@ def verify_commit(src: Path) -> None: check=True, ).stdout.strip() - if head != expected: + if head != COOLPROP_COMMIT: raise RuntimeError( - f"CoolProp {COOLPROP_VERSION} resolved to commit {head}, expected {expected}. " + f"CoolProp {COOLPROP_VERSION} resolved to commit {head}, expected {COOLPROP_COMMIT}. " "The upstream tag moved; do not build from it until this is understood." ) diff --git a/uv.lock b/uv.lock index d56e55f..30c1deb 100644 --- a/uv.lock +++ b/uv.lock @@ -517,7 +517,6 @@ name = "encomp" version = "1.8.0" source = { editable = "." } dependencies = [ - { name = "coolprop" }, { name = "numpy" }, { name = "pint" }, { name = "polars" }, @@ -558,10 +557,12 @@ dev = [ { name = "sphinx-inline-tabs" }, { name = "ty" }, ] +oracle = [ + { name = "coolprop" }, +] [package.metadata] requires-dist = [ - { name = "coolprop", specifier = ">=8.0.0,<9" }, { name = "numpy", specifier = ">=2.1" }, { name = "pint", specifier = ">=0.25.3" }, { name = "polars", specifier = ">=1.42" }, @@ -602,6 +603,7 @@ dev = [ { name = "sphinx-inline-tabs" }, { name = "ty" }, ] +oracle = [{ name = "coolprop", specifier = "==8.0.0" }] [[package]] name = "executing"