From 290972214fe6570aaae952138257b2963f33c11d Mon Sep 17 00:00:00 2001 From: Andrei Danila Date: Mon, 3 Aug 2026 10:22:33 +0100 Subject: [PATCH 01/15] skeleton --- FORWARD_ADJOINT_REFERENCE.md | 603 +++++++++++++++++++++++++++++++++++ stride/problem/base.py | 36 ++- stride/problem/data.py | 11 +- stride/problem/domain.py | 4 +- 4 files changed, 640 insertions(+), 14 deletions(-) create mode 100644 FORWARD_ADJOINT_REFERENCE.md diff --git a/FORWARD_ADJOINT_REFERENCE.md b/FORWARD_ADJOINT_REFERENCE.md new file mode 100644 index 0000000..29e64e2 --- /dev/null +++ b/FORWARD_ADJOINT_REFERENCE.md @@ -0,0 +1,603 @@ +# Stride Forward & Inverse Scripts — Object & Argument Reference + +A practical reference for the objects most used in Stride forward-modelling and +inversion (FWI) scripts, centred on **the arguments accepted by `forward()` and +`adjoint()`** and by the PDE operators they drive. + +Line references point at the public repo +(`/Users/andreidanila/code/sonalis/stride`). A final section documents the +extra functionality and advanced examples in the private fork +(`/Users/andreidanila/code/sonalis/stride-private`). + +--- + +## 1. The two canonical script shapes + +Every script is an `async def main(runtime)` launched with `mosaic.run(main)` +(run on a cluster with `mrun -nw python script.py`). + +**Forward script** ([`breast2D/01_script_forward.py`](stride_examples/examples/breast2D/01_script_forward.py)): + +```python +space = Space(shape=(356, 385), extra=(50, 50), absorbing=(40, 40), spacing=0.5e-3) +time = Time(start=0., step=0.08e-6, num=2500) +problem = Problem(name='anastasio2D', space=space, time=time) + +vp = ScalarField(name='vp', grid=problem.grid); vp.load('...TrueModel.h5') +problem.medium.add(vp) +problem.transducers.default() +problem.geometry.default('elliptical', 128) +problem.acquisitions.default() +for shot in problem.acquisitions.shots: + shot.wavelets.data[0, :] = wavelets.tone_burst(0.5e6, 3, time.num, time.step) + +pde = IsoAcousticDevito.remote(grid=problem.grid, len=runtime.num_workers) +await forward(problem, pde, vp) # <-- entry point 1 +``` + +**Inverse script** ([`breast2D/02_script_inverse.py`](stride_examples/examples/breast2D/02_script_inverse.py)): + +```python +vp = ScalarField.parameter(name='vp', grid=problem.grid, needs_grad=True) +vp.fill(1500.) +problem.medium.add(vp) +problem.acquisitions.load(path=problem.output_folder, project_name=problem.name, version=0) + +pde = IsoAcousticDevito.remote(grid=problem.grid, len=runtime.num_workers) +loss = L2DistanceLoss.remote(len=runtime.num_workers) +optimiser = GradientDescent(vp, step_size=10, + process_grad=ProcessGlobalGradient(), + process_model=ProcessModelIteration(min=1400., max=1700.)) +optimisation_loop = OptimisationLoop() + +for block, freq in optimisation_loop.blocks(num_blocks, max_freqs): + await adjoint(problem, pde, loss, optimisation_loop, optimiser, vp, # <-- entry point 2 + num_iters=8, select_shots=dict(num=16, randomly=True), + f_max=freq, max_freqs=max_freqs) +``` + +Two key rules that explain most of the API: + +- **`.remote(...)`** turns an `Operator`/PDE class into a distributed *tessera* + spread across `len=runtime.num_workers` workers. Construction kwargs + (`grid`, `space`, `time`, `name`, `cached_operator`, …) go here. +- **`needs_grad=True`** on a variable (via `ScalarField.parameter(...)`) is what + switches the whole pipeline from pure forward to gradient-tracking adjoint + mode. `forward()` will automatically save the wavefield when it sees a + `needs_grad` input; `adjoint()` builds the adjoint graph from it. + +--- + +## 2. `forward(problem, pde, *args, **kwargs)` + +Defined in [`stride/__init__.py:53`](stride/__init__.py#L53). Runs `pde` once +per shot, stores the result in `shot.observed`, and (by default) appends it to +the observed HDF5 file. `*args` (e.g. `vp`, `rho`, `alpha`) and all extra +`kwargs` are **forwarded to the PDE**. + +| kwarg | type | default | meaning | +|---|---|---|---| +| `dump` | bool | `True` | Write forward result to disk (`observed` file). When `True` it also pre-loads any existing observed (version 0) so already-run shots are skipped. | +| `shot_ids` | list / int | `None` | Specific shots to run. `None` → all *remaining* shots (those with no observed yet). If none remain it logs a warning and returns. | +| `deallocate` | bool | `False` | Free `shot.observed` right after each shot (memory saving). | +| `safe` | bool | `True` | Discard workers that fail mid-execution instead of aborting. | +| `platform` | str | `'cpu'` | `'cpu'`, `'gpu'`, `'nvidia-acc'`, `'nvidia-cuda'`. Triggers GPU device round-robin across workers. | +| `devices` | list | `None` | Explicit GPU device ids; `None` → auto (`gpu_count()`). | +| `*args` | — | — | Positional PDE inputs (medium fields): `vp`, or `vp, vs, rho` (elastic), etc. Published to all workers. | +| `**kwargs` | — | — | Everything else is passed straight to the PDE `forward` (see §4). | + +Notable behaviour: the loop hands each `shot_id` to a worker, builds a +`sub_problem = problem.sub_problem(shot_id)`, calls +`pde(wavelets, *args, problem=sub_problem, runtime=worker, **kwargs).result()`, +writes `shot.observed.data[:] = traces.data`, and raises if NaN/Inf appear. + +--- + +## 3. `adjoint(problem, pde, loss, optimisation_loop, optimiser, *args, **kwargs)` + +Defined in [`stride/__init__.py:172`](stride/__init__.py#L172). This is the FWI +iteration driver. Per iteration it: selects shots → pre-processes +wavelets/observed → runs the PDE forward (saving the wavefield) → +post-processes modelled+observed traces → evaluates `loss` → runs the adjoint +(`fun.adjoint(...)`) to accumulate the gradient → takes an optimiser step → +dumps the updated variable. + +| kwarg | type | default | meaning | +|---|---|---|---| +| `num_iters` | int | `1` | Iterations to run **within the current block**. | +| `select_shots` | dict | `{}` | Shot-selection rules per iteration, forwarded to `Acquisitions.select_shot_ids` — e.g. `dict(num=16, randomly=True)`, or `dict(start=, end=, every=)`. | +| `lazy_loading` | bool | `False` | Load shot data each iteration and deallocate after, to save memory. | +| `dump` | bool | `True` | Save the updated optimiser variable after each iteration (enables restart). | +| `safe` | bool | `True` | Discard failing workers. | +| `f_min` | float | `None` | High-pass corner for filtering wavelets/traces. `None` → no high-pass. | +| `f_max` | float | `None` | Low-pass corner (frequency continuation). `None` → no low-pass. Usually set per block. | +| `filter_traces` | bool | `True` | Whether the trace-processing pipeline filters modelled/observed. | +| `filter_wavelets` | bool | `= filter_traces` | Whether wavelets/observed are band-limited before the PDE. | +| `filter_wavelets_relaxation` | float | `0.75` | Filter roll-off relaxation for wavelets. | +| `filter_traces_relaxation` | float | `0.75` (or `1.0` if no wavelet filtering) | Filter roll-off relaxation for traces. | +| `step_size` | float / `LineSearch` | `optimiser.step_size` | Step length; a `LineSearch` instance triggers the residual-keeping test-step loop. | +| `platform` / `devices` | str / list | `'cpu'` / `None` | GPU control, same as `forward()`. | +| `*args` | — | — | The medium variable(s) being inverted (also the `wrt` of the gradient), e.g. `vp` (must be a `.parameter(..., needs_grad=True)`). | +| `**kwargs` | — | — | Passed to the four processing pipelines, the PDE, the loss, and `optimiser.step`. Includes everything in §4 (`kernel`, `interpolation_type`, `boundary_type`, `devito_config`, …) and pipeline knobs (§7). | + +The `max_freqs=[...]` kwarg seen in examples is not consumed by `adjoint()` +itself; it is passed through so the pipelines/optimiser see the full schedule. + +**Frequency continuation** is driven outside `adjoint()` by the loop: +```python +for block, freq in optimisation_loop.blocks(num_blocks, max_freqs): + await adjoint(..., f_max=freq, max_freqs=max_freqs) +``` + +--- + +## 4. The PDE operator — `IsoAcousticDevito` + +The single most important object. It is the second-order isotropic **acoustic** +wave equation on Devito. Defined in +[`iso_acoustic/devito.py:24`](stride/physics/iso_acoustic/devito.py#L24). + +### 4.1 Construction (`.remote(...)`) + +```python +pde = IsoAcousticDevito.remote(grid=problem.grid, len=runtime.num_workers) +# or: IsoAcousticDevito.remote(space=space, time=time) +``` + +| construction kwarg | meaning | +|---|---| +| `grid` | Existing `Grid` (preferred). Alternatively pass `space=`/`time=`/`slow_time=`. | +| `len` | Number of worker replicas (`runtime.num_workers`). | +| `name` | Optional PDE name. | +| `cached_operator` | Reuse a compiled operator/grid stored in the worker warehouse across PDE instances. | +| `dev_grid` | Supply an existing `GridDevito` (advanced/shared setups). | + +Class-level: `space_order = 10`, `time_order = 2`. + +### 4.2 Forward-call arguments (the crux) + +Whether via `forward(problem, pde, vp, **kwargs)` or a direct +`await pde(wavelets, vp, problem=sub_problem, **kwargs).result()`, these are the +inputs and options consumed by `before_forward` / `run_forward` / `after_forward` +([`iso_acoustic/devito.py:220`](stride/physics/iso_acoustic/devito.py#L220)): + +**Positional / field inputs** + +| arg | type | default | meaning | +|---|---|---|---| +| `wavelets` | `Traces` | required | Source wavelets (`shot.wavelets`). Supplied automatically by `forward()`/`adjoint()`. | +| `vp` | `ScalarField` | required | Compressional speed of sound, m/s. | +| `rho` | `ScalarField` | `None` | Density, kg/m³. `None` → homogeneous. | +| `alpha` | `ScalarField` | `None` | Attenuation, dB/cm. `None` → lossless. | +| `problem` | `Problem`/`SubProblem` | required | The sub-problem (one shot). Injected by the driver. | + +**Physics / discretisation options** + +| kwarg | type | default | meaning | +|---|---|---|---| +| `kernel` | str | auto (`'OT2'`/`'OT4'`) | Time-stepping order: `'OT2'` (2nd) or `'OT4'` (4th). Auto-chosen from `dt` if unset. | +| `boundary_type` | str | `'sponge_boundary_2'` | Absorbing boundary. Also `'complex_frequency_shift_PML_2'` (lower OT4 stability). | +| `interpolation_type` | str | `'linear'` | Source/receiver interpolation: `'linear'` (bi/tri-linear) or `'hicks'` (sinc). | +| `attenuation_power` | int / None | `0` | Power of the attenuation law when `alpha` given (`0`, `2`, or `None`). | +| `drp` | bool | `False` | Dispersion-relation-preserving coefficients (build-dependent). | +| `diff_source` | bool | `False` | Inject the source as its 1st time-derivative instead of as-is. | +| `adaptive_boxes` | bool | `False` | Adaptive computational boxes (DevitoPRO). | +| `local_prec` | bool | `True` | Local preconditioning (build-dependent). | + +**Wavefield saving / gradient options** + +| kwarg | type | default | meaning | +|---|---|---|---| +| `save_wavefield` | bool | auto | Save forward wavefield for the gradient. Auto-`True` when any input `needs_grad`. | +| `time_bounds` | (int,int) | `(0, time.extended_num)` | Timestep window over which the wavefield is saved. | +| `save_undersampling` | int | auto (bandwidth) | Temporal undersampling factor when saving the wavefield. | +| `save_compression` | str | `None` (2D) / `'bitcomp'` (3D) | Wavefield compression (DevitoPRO/GPU only). | +| `save_interpolation` | bool | build-dependent | Cubic-spline interpolation of the saved wavefield. | +| `stream_wavefield` | bool / str | `True` | Streaming layer strategy: `True`/`False` or `'disk'`, `'host'`, `'device'`, `'disk-host'`, `'host-device'`, `'disk-host-device'`, `'no-layers'`. | +| `spill_wavefield` | bool | `False` | Spill the wavefield to disk/host. | +| `cache_forward` | bool | `False` | Cache the forward wavefield (in memory or `cache_location`) for reuse in the adjoint. | +| `cache_location` | str | `None` | Directory to cache the wavefield to disk. | +| `nbits_compression` | int | `9` | Bits for compressed streaming (maps to `devito_args['nbits']`). | + +**Wavefield dumping (debug / imaging)** + +| kwarg | type | default | meaning | +|---|---|---|---| +| `dump_forward_wavefield` | bool / int | `False` | Dump forward wavefield. `True` → every `save_undersampling`; int → every N steps. | +| `dump_adjoint_wavefield` | bool / int | `False` | Same for the adjoint wavefield. | +| `dump_wavefield_id` | int | shot id | Only dump this shot's wavefields. | + +**Platform / Devito plumbing** + +| kwarg | type | default | meaning | +|---|---|---|---| +| `platform` | str | `None`/`'cpu'` | `None`/`'cpu'` or `'nvidia-acc'` (OpenACC) / `'nvidia-cuda'`. | +| `devito_config` | dict | `{}` | Devito config applied **before** operator generation, e.g. `{'opt': ('advanced', {'index-mode': 'int64'})}`, `{'compiler':'pgcc','language':'openacc','platform':'nvidiaX'}`. `platform='nvidia-acc'` is shorthand for the latter. | +| `devito_args` | dict | `{}` | Args passed when **calling** the compiled operator, e.g. `{'autotune': 'off'}`, `{'deviceid': N}`, `{'nbits': 9}`. | +| `deallocate` | bool | `False` | Free Devito buffers (boundary, `p`, `src`, `rec`, `vp`, …) after the run. | + +### 4.3 Adjoint-call arguments + +The adjoint side (`before_adjoint`/`run_adjoint`/`after_adjoint`, +[`devito.py:689`](stride/physics/iso_acoustic/devito.py#L689)) is invoked by +`fun.adjoint(**kwargs)` inside the `adjoint()` driver — you rarely call it +directly. Its inputs are `(adjoint_source, wavelets, vp, rho=None, alpha=None, +**kwargs)` where `adjoint_source` is produced by the loss. It honours +`dump_adjoint_wavefield`, `dump_wavefield_id`, `cache_forward`, `time_bounds`, +`platform`, `deallocate`, and the same physics kwargs, then returns the +gradients via the `get_grad_*` methods. + +Gradients are produced per variable through the naming convention +(`prepare_grad_vp` / `init_grad_vp` / `get_grad_vp`), so `vp`, `rho`, and +`wavelets` can each be inverted for when flagged `needs_grad` +(see [`problem_type.py:202`](stride/physics/problem_type.py#L202)). + +### 4.4 Calling the PDE directly + +`forward()`/`adjoint()` are conveniences. You can call the tessera yourself, +which is how the homogeneous-medium test sweeps work +([`homogeneous_acoustic/forward_2D.py:109`](stride_examples/examples/homogeneous_acoustic/forward_2D.py#L109)): + +```python +sub_problem = problem.sub_problem(shot.id) +traces = await pde(sub_problem.shot.wavelets, vp, + problem=sub_problem, diff_source=True, + kernel='OT4', interpolation_type='hicks', + boundary_type='complex_frequency_shift_PML_2', + rho=rho, alpha=alpha, attenuation_power=2).result() +data = traces.data # numpy array, shape (num_receivers, time.num) +await pde.clear_operators() # force recompile when config changes +``` + +`forward()` returns nothing (it writes `shot.observed`); a direct call returns a +`Traces` object whose `.data` is the modelled gather. + +### 4.5 Other PDE operators (same call convention) + +- **`IsoElasticDevito`** ([`iso_elastic/devito.py:16`](stride/physics/iso_elastic/devito.py#L16)) — + stress-strain elastic wave equation. Forward inputs are + `(wavelets, vp, vs, rho, problem=...)`; `space_order=10`, `time_order=1`, + `boundary_type='sponge_boundary_1'`. Multi-parameter: pass `vp, vs, rho=rho`. +- **`MarmottantDevito`** ([`marmottant/devito.py`](stride/physics/marmottant/devito.py)) — + microbubble model. + +--- + +## 5. Problem-definition objects + +All in `stride/problem/`. These build the `problem` you hand to +`forward()`/`adjoint()`. + +### `Space` — [`domain.py:9`](stride/problem/domain.py#L9) +Spatial grid = inner domain + padding. +`Space(shape, spacing, extra, absorbing)`. +- `shape` (tuple) inner grid points; `spacing` (tuple or scalar float, m); + `extra` padding points per axis; `absorbing` portion of `extra` used for the + boundary. Exposes `.dim`, `.limit` (physical size), `.extended_shape`, + `.inner`, `.resample(...)`. + +### `Time` — [`domain.py:263`](stride/problem/domain.py#L263) +`Time(start, step, num, stop)` — give any three. `num` must be `int`. Exposes +`.start/.step/.num/.stop`, `.extended_num`, `.extend(...)`, `.resample(...)`. + +### `SlowTime` — [`domain.py:417`](stride/problem/domain.py#L417) +Frame/acquisition sampling for multi-frame data: +`SlowTime(frame_rate|frame_step, acq_rate|acq_step, num_frame, num_acq)`. + +### `Grid` — [`domain.py:528`](stride/problem/domain.py#L528) +`Grid(space, time, slow_time)` — a bundle. `problem.grid` is usually what you +pass to `.remote(grid=...)` and to fields. + +### `Problem` — [`problem.py:13`](stride/problem/problem.py#L13) +`Problem(name, space=, time=)` (or `grid=`). Top-level container. Attributes +auto-created: `.medium`, `.transducers`, `.geometry`, `.acquisitions`, `.grid`. +- `.sub_problem(shot_id)` → a `SubProblem` with `.shot`, `.shot_id` (built per + shot inside the drivers). +- `.plot()`, `.load(...)`, `.dump(...)`, `.output_folder`/`.input_folder` + (default `cwd`). +- `.space_resample(new_spacing)`, `.time_resample(new_step, new_num)` for + multi-resolution FWI. + +### `ScalarField` — [`data.py:790`](stride/problem/data.py#L790) +The medium-field type (`vp`, `rho`, `alpha`). Two construction paths: +- Forward / known model: `ScalarField(name='vp', grid=problem.grid)` then + `.load(path)` or `.fill(value)`. +- Inversion variable: `ScalarField.parameter(name='vp', grid=problem.grid, + needs_grad=True)`. `.parameter()` is injected by `@mosaic.tessera` and returns + the optimisable/distributed variant; `needs_grad=True` enables gradients. +- `time_dependent=`/`slow_time_dependent=` prepend time axes. +- Members: `.data` (inner view), `.extended_data` (with padding), `.fill(v)`, + `.plot()`, `.load()`/`.dump()`, `.needs_grad`, `.grad`, `.clear_grad()`. + +### `Traces` — [`data.py:1383`](stride/problem/data.py#L1383) +Time traces indexed by transducer id — the type of `shot.wavelets` and +`shot.observed`. Write with `shot.wavelets.data[i, :] = ...`. `.data` shape is +`(num_transducers, time.num)`. `.plot(plot_type='gather'|'spectrum')`, +`.get(id)`, `.alike(...)`. `DiskTraces` is the lazy on-disk variant. + +### `Medium` — [`medium.py:8`](stride/problem/medium.py#L8) +Named-field container. `problem.medium.add(vp)`; access `medium.vp` / +`medium['vp']`. `.load()/.dump()/.plot()` iterate all fields. + +### `Transducers` — [`transducers.py:11`](stride/problem/transducers.py#L11) +Registry of transducer devices. `problem.transducers.default()` creates one +`PointTransducer(0)`. + +### `Geometry` — [`geometry.py:106`](stride/problem/geometry.py#L106) +Transducer *locations*. `problem.geometry.default('elliptical', num_locations)` +(2D) auto-computes radius/centre from `space.limit`; +`'ellipsoidal', num_locations, radius, centre, theta=, threshold=` (3D). +`.coordinates`, `.locations`, `.num_locations`, `.plot()`. + +### `Acquisitions` — [`acquisitions.py:733`](stride/problem/acquisitions.py#L733) +The set of shots (the data). +- `.default()` — one single-source shot per location, all locations as receivers. +- `.load(path=, project_name=, version=0, shot_ids=None, fast=False)`. +- `.select_shot_ids(num=, start=, end=, every=1, randomly=False)` — stateful + per-iteration selection (this is what `select_shots` in `adjoint()` feeds). +- `.remaining_shot_ids` — shots with no observed yet (what `forward()` runs). +- `.shots`, `.shot_ids`, `.num_shots`, `.plot()`, `.reset_selection()`. + +### `Shot` — [`acquisitions.py:57`](stride/problem/acquisitions.py#L57) +One acquisition event: `.wavelets` (per source), `.observed` (per receiver), +`.delays`, `.source_coordinates`, `.receiver_coordinates`, +`.num_sources`/`.num_receivers`. + +### Wavelet helpers — [`utils/wavelets.py`](stride/utils/wavelets.py) +- `tone_burst(centre_freq, n_cycles, n_samples, dt, envelope='gaussian')` +- `ricker(centre_freq, n_samples, dt)` +- `continuous_wave(centre_freq, n_samples, dt, ramp_length=4, phase=0)` + +### Asset fetch — [`utils/fetch.py`](stride/utils/fetch.py) +`fetch('anastasio2D', dest='data/...h5')` downloads a known release asset once. + +--- + +## 6. Optimisation objects (inverse scripts) + +### `L2DistanceLoss` — [`loss/l2_distance.py:14`](stride/optimisation/loss/l2_distance.py#L14) +`f = ½‖modelled − observed‖²`. `L2DistanceLoss.remote(len=runtime.num_workers)`. +- Ctor kwarg `d_sample=4` (downsampling of the residual passed to the functional). +- `forward(modelled, observed, **kwargs)` → `FunctionalValue`; consumes + `problem`/`shot_id`, forwards `keep_residual`. +- `adjoint(d_fun, modelled, observed)` → `(grad_modelled, grad_observed)` — the + adjoint source. + +### `GradientDescent` / `LocalOptimiser` — [`optimisers/`](stride/optimisation/optimisers/) +> There is **no `Adam`** in the public repo — only `GradientDescent` (and the +> `LocalOptimiser` base). Adam/NAdam/RAdam/SGD/AdaBelief/AGD live in the private +> fork (§8). + +`GradientDescent(variable, step_size=1., process_grad=..., process_model=...)`. +Base `LocalOptimiser` ctor kwargs: + +| kwarg | default | meaning | +|---|---|---| +| `variable` | — | Variable to optimise (must be `needs_grad`). | +| `step_size` | `1.` | Float or `LineSearch`. | +| `test_step_size` | `1.` | Multiplier on the processed gradient. | +| `force_step` | `False` | Skip step clipping/capping. | +| `max_step` | `None` | Cap on step magnitude. | +| `process_grad` | `ProcessGlobalGradient(**kwargs)` | Gradient pre-processing pipeline. | +| `process_model` | `ProcessModelIteration(**kwargs)` | Model post-processing pipeline. | +| `reset_block` | `False` | Reset optimiser state each block (flag, not method). | +| `reset_iteration` | `False` | Reset optimiser state each iteration. | +| `dump_grad` / `dump_prec` | `False` | Debug dumps. | + +Methods: `await step(step_size=None, grad=None, step_loop=..., **kwargs)`, +`clear_grad()`, `reset()`, `dump()/load()`. + +### `OptimisationLoop` / `Block` / `Iteration` — [`optimisation_loop.py`](stride/optimisation/optimisation_loop.py) +- `OptimisationLoop(name='optimisation_loop')`. Properties `.num_blocks`, + `.current_block`; `.blocks(num, *iters, restart=False, restart_id=-1)` is the + generator you loop over (zips extra sequences like `max_freqs`). +- `Block.iterations(num, *iters, ...)` yields `Iteration`s; `.total_loss`, `.id`. +- `Iteration`: `.add_loss(fun)`, `.add_submitted/.add_completed(shot)`, + `.next_run()` (line-search test steps), `.id`, `.abs_id`, `.total_loss`, + `.prev_run`. + +### Processing pipelines — [`pipelines/default_pipelines.py`](stride/optimisation/pipelines/default_pipelines.py) +`.remote(...)` operators the `adjoint()` driver builds automatically; you tune +them by passing kwargs through `adjoint()`. + +| pipeline | role | default steps (kwarg → default) | +|---|---|---| +| `ProcessWavelets` | pre-process source wavelets | `check_traces`(T), `filter_traces`(T), `shift_traces`, `resonance_filter`(F) | +| `ProcessObserved` | pre-process observed | `check_traces`, `filter_traces` | +| `ProcessWaveletsObserved` | joint step | `differentiate_traces`(T) | +| `ProcessTraces` | modelled+observed before the loss | `check_traces`(T), `filter_offsets`(F), `mute_first_arrival`(T), `mute_traces`(T), `filter_traces`(T), `agc`(F), `norm_per_shot`(T)/`norm_per_trace`(F), `scale_per_*`(F), `time_tweaking`(T), `time_weighting`(T) | +| `ProcessGlobalGradient` | default `process_grad` | `mask_field`(`mask_grad`=T), `smooth_field`(`smooth_grad`=T), `norm_field`(`norm_grad`=T) | +| `ProcessModelIteration` | default `process_model` | `clip` (uses `min=`/`max=`) | + +> Several `ProcessTraces` steps (`resonance_filter`, `differentiate_traces`, +> `filter_offsets`, `mute_first_arrival`, `agc`, `time_tweaking`, +> `time_weighting`) are added *non-raising*: they are no-ops in the public repo +> unless the step is registered — and those step classes ship in the **private +> fork** (§8). + +### Individual steps — [`pipelines/steps/`](stride/optimisation/pipelines/steps/) +Registered (public): `filter_traces`, `norm_per_shot`, `norm_per_trace`, +`scale_per_shot`, `scale_per_trace`, `norm_field`, `smooth_field`, `mask_field`, +`mute_traces`, `clip`, `check_traces`, `dump`, `shift_traces`. Common kwargs: + +| step | key kwargs | +|---|---| +| `FilterTraces` | `f_min`, `f_max`, `filter_type` (`'cos'`/`'butterworth'`/`'fir'`), `filter_relaxation` | +| `MuteTraces` | `f_max`, `filter_relaxation` | +| `Clip` | `min`, `max` | +| `SmoothField` | `smooth_sigma` (default 0.25 = 25% of a cell) | +| `NormField` | `global_norm` (F), `norm_guess_change` (0.5) | +| `MaskField` | `mask`, `mask_rampoff` (10) | +| `ScalePer*` | `scale_to`, `relative_scale` (T) | +| `CheckTraces` | `raise_incorrect` (T), `filter_incorrect` (F) | +| `ShiftTraces` | `f_max`, `filter_relaxation` | + +### `LineSearch` — [`step_length/line_search.py:8`](stride/optimisation/step_length/line_search.py#L8) +Abstract step-length base (`init_search`, `next_step`). Passed as +`step_size=LineSearch(...)`; the optimiser drives it with a `step_loop` callable. +The public repo ships only the abstract base — concrete searches are in the +private fork. + +--- + +## 7. Quick kwarg cheat-sheet + +Passing these through `forward(problem, pde, vp, ...)` / +`adjoint(problem, pde, loss, loop, opt, vp, ...)` reaches the right layer: + +```python +# --- performance / hardware --- +platform='nvidia-cuda' # or 'nvidia-acc', 'gpu', 'cpu' +devices=[0,1] # explicit GPUs +devito_config={'opt': ('advanced', {'index-mode': 'int64'})} +devito_args={'autotune': 'off'} +deallocate=True # free device buffers each shot + +# --- physics --- +kernel='OT4' # or 'OT2' +interpolation_type='hicks' # or 'linear', 'sinc' +boundary_type='complex_frequency_shift_PML_2' +rho=rho_field, alpha=alpha_field, attenuation_power=2 +diff_source=True + +# --- forward wavefield / imaging --- +save_wavefield=True, time_bounds=(0, N), save_undersampling=4 +dump_forward_wavefield=True, dump_adjoint_wavefield=8, dump_wavefield_id=0 +cache_forward=True, cache_location='/scratch' + +# --- inversion control (adjoint only) --- +num_iters=8 +select_shots=dict(num=16, randomly=True) +f_min=0.1e6, f_max=0.5e6, max_freqs=[0.3e6, 0.4e6, 0.5e6, 0.6e6] +lazy_loading=True +step_size=10 # or a LineSearch instance + +# --- pipeline toggles (adjoint only) --- +filter_traces=True, filter_wavelets=True +mask_grad=True, smooth_grad=True, norm_grad=True +smooth_sigma=0.5, norm_guess_change=0.25 +``` + +--- + +## 8. Private fork (`stride-private`) — extensions & advanced examples + +The private repo forks an **older** public baseline and adds specialised +operators, losses, optimisers and pipeline steps. (Its acoustic operator is +older/smaller than the current public one; where the two diverge, the private +defaults differ — e.g. `drp=True` and +`boundary_type='hybrid_interpolating_boundary_2'` by default.) + +### 8.1 Advanced invocation patterns (copy-paste worthy) + +**GPU multi-block acoustic FWI** — `examples/alpha2D/inverse.py:74`: +```python +await adjoint(problem, pde, loss, optimisation_loop, optimiser, vp, + num_iters=num_iters, + select_shots=dict(num=12, randomly=True), + f_max=freq, max_freqs=max_freqs, + kernel='OT4', fw3d_mode=True, + interpolation_type='hicks', platform='nvidia-cuda') +``` +matching forward `alpha2D/forward.py:62`: +```python +await forward(problem, pde, vp, kernel='OT4', fw3d_mode=True, + interpolation_type='hicks', platform='nvidia-cuda') +``` + +**Elastic forward, multi-parameter + Devito opt** — `examples/alpha2D_elastic/forward.py:90`: +```python +await forward(problem, pde, vp, vs, rho=rho, shot_ids=[0], + interpolation_type='sinc', dump=False, deallocate=False, + dump_forward_wavefield=False, + devito_config={'opt': ('advanced', {'index-mode': 'int64'})}) +``` + +**3D elastic, autotune off** — `examples/alpha3D_elastic/forward.py:113`: +```python +await forward(problem, pde, vp, rho=rho, shot_ids=[0], + interpolation_type='sinc', kernel='OT4', + dump=False, deallocate=False, devito_args={'autotune': 'off'}) +``` + +**Optimiser with line search** — `examples/alpha2D/inverse.py:58`: +```python +optimiser = GradientDescent(vp, step_size=LineSearch(), + process_grad=ProcessGlobalGradient(), + process_model=ProcessModelIteration(min=1450., max=3000.)) +``` + +**Custom driver for source/receiver-IR inversion (SRI)** — +`examples/SRI/SRI_transmit_optim_wavelet.py` replaces the library `adjoint()` +with its own `adjoint_finite(problem, pde, loss, optimisation_loop, optimiser, +*args, **kwargs)` that builds `ProcessWavelets.remote(f_min=, f_max=)` / +`ProcessTraces.remote(...)`, convolves the wavelet with each transducer's +`transmit_ir` via a `Convolution.remote(...)` operator, runs the PDE, and calls +`await fun.adjoint(**kwargs)`. Manual per-worker GPU assignment: +```python +devito_args = kwargs.get('devito_args', {}) +devito_args['deviceid'] = devices[worker.indices[1] % num_gpus] +kwargs['devito_args'] = devito_args +``` + +**Multi-resolution FWI** — `examples/alpha2D_resample/` uses +`problem.space_resample()` / `problem.time_resample()` between blocks. + +### 8.2 Added physics operators +- `IsoAcousticAnalyticNumba` (analytic Green's-function acoustic, Numba + CPU/CUDA) — extra kwargs `compute_forward`, `vp_constant`, `gradient_crop`, + `threads_per_block`, `blocks_per_grid`, `diff_source`, `save_undersampling`. +- `TransportDevito` / `SLTransportNumpy` (semi-Lagrangian transport) — kwargs + `interpolating`, `velocity_interpolation`, `taylor_accuracy`, `parallel_time`, + `num_procs`, `fill_holes`; forward inputs `(sigma_0, u, ...)`. +- `FlowInjectionDevito` — kwargs `interpolation_type`, `t_i`. +- New boundaries `HybridHigdonBoundary2`, `HybridInterpolatingBoundary2`. + +### 8.3 Added loss functions (`stride_private/optimisation/loss/`) +| loss | key kwargs | +|---|---| +| `AdaptiveWaveformLoss` (AWI) | `pad`, `augment`, `gamma`, `eta`, `mode='reverse'`, `type='standard'`, `d_sample=4` | +| `OptimalTransportLoss` (GSOT) | `p=2`, `mode='taot'`, `e_start/e_end/e_fac`, `taot_algorithm='dtw'`, `reg`, `max_iter`, `tol`, `d_sample` | +| `FrequencyControllableEnvelopeLoss` (FCEI) | `pad`, `p`, `filter_type='butterworth'`, `filter_order=8`, `num_modify`, `weight_l2`, `weight_fcei` | +| `ReverseTimeMigration` (RTM) | imaging condition | +| `DoubleDifferenceLoss` | double-difference misfit | +| `MultiLoss` | `weight_l2`, `weight_ot`, `weight_fcei`, `fcei_mode`, `mute_traces`, `mute_threshold`, `d_sample` | + +### 8.4 Added optimisers & steps +- Optimisers: `Adam`, `NAdam`, `RAdam` (`betas=(0.9,0.999)`, `eps=1e-8`, `t`), + `SGD`, `AdaBelief`, `AGD` — all on a `LocalOptimiserSaved` base. +- Pipeline steps (these are the ones the public default pipelines reference but + don't ship): `AGC`, `DifferentiateTraces`, `FilterOffsets`, `MuteFirstArrival`, + `ResonanceFilter`, `TimeTweaking`, `TimeWeighting`. +- Constraint `DivergenceFree`; helper operators `Convolution`, `Split`, + `Replicate`, `Unflatten`, `SVDFilter`, `DiffFilter`, `FreqWindowingFilter`. +- Fullwave3D interop (`.ttr/.pgy/.vtr`) in `utils/fullwave.py`. + +--- + +## 9. Mental model (how the pieces connect) + +``` +Space + Time ─► Grid ─► Problem ─┬─ medium (ScalarField vp/rho/alpha) + ├─ transducers (.default) + ├─ geometry (.default 'elliptical'/'ellipsoidal') + └─ acquisitions (Shots: wavelets + observed) + +pde = IsoAcousticDevito.remote(grid=..., len=workers) # the physics + +FORWARD: forward(problem, pde, vp, **pde_kwargs) # writes shot.observed + +INVERSE: loss = L2DistanceLoss.remote() + optimiser = GradientDescent(vp.parameter(needs_grad=True), + process_grad=ProcessGlobalGradient(), + process_model=ProcessModelIteration(min,max)) + loop = OptimisationLoop() + for block, freq in loop.blocks(num_blocks, max_freqs): + adjoint(problem, pde, loss, loop, optimiser, vp, + num_iters=, select_shots=, f_max=freq, **pde_kwargs) + └─ per iter: ProcessWavelets/Observed → pde.forward (saves wavefield) + → ProcessTraces → loss.forward → fun.adjoint (grad) + → optimiser.step → dump updated vp +``` + +*Generated as a code reference for Stride forward/inverse scripting.* diff --git a/stride/problem/base.py b/stride/problem/base.py index cbbbc0e..52c3a09 100644 --- a/stride/problem/base.py +++ b/stride/problem/base.py @@ -4,7 +4,7 @@ from .domain import Space, Time, SlowTime, Grid -__all__ = ['Gridded', 'Saved', 'GriddedSaved', 'ProblemBase'] +__all__ = ['Gridded', 'Meshed', 'Saved', 'GriddedSaved', 'MeshedSaved', 'ProblemBase'] class Gridded: @@ -245,12 +245,16 @@ def load(self, *args, **kwargs): with h5.HDF5(*args, **kwargs, mode='r') as file: description = file.load(filter=kwargs.pop('filter', None), only=kwargs.pop('only', None)) - # TODO If there's already a grid and they don't match, resample instead if 'space' in description and self._grid.space is None: - space = Space(shape=description.space.shape, - spacing=description.space.spacing, - extra=description.space.extra, - absorbing=description.space.absorbing) + if 'shape' in description: + space = Space(shape=description.space.shape, + spacing=description.space.spacing, + extra=description.space.extra, + absorbing=description.space.absorbing) + elif 'nodes' in description: + space = MeshedSpace(nodes=description.space.nodes) + else: + raise Exception self._grid.space = space @@ -289,12 +293,20 @@ def grid_description(self): if self.space is not None: space = self.space - grid_description['space'] = { - 'shape': space.shape, - 'spacing': space.spacing, - 'extra': space.extra, - 'absorbing': space.absorbing, - } + if isinstance(space, Space): + grid_description['space'] = { + 'shape': space.shape, + 'spacing': space.spacing, + 'extra': space.extra, + 'absorbing': space.absorbing, + } + elif isinstance(space, MeshedSpace): + #stand-in attribute + grid_description['space'] = { + 'nodes': space.nodes + } + else: + raise Exception if self.time is not None: time = self.time diff --git a/stride/problem/data.py b/stride/problem/data.py index 8878306..3a158cb 100644 --- a/stride/problem/data.py +++ b/stride/problem/data.py @@ -26,7 +26,8 @@ __all__ = ['Data', 'StructuredData', 'Scalar', 'ScalarField', 'VectorField', 'Traces', - 'DiskTraces', 'ArtifactTraces', 'SparseField', 'SparseCoordinates'] + 'DiskTraces', 'ArtifactTraces', 'SparseField', 'SparseCoordinates', + 'MeshedData', 'MeshedField'] def inv_transform(x): @@ -779,6 +780,14 @@ def __set_desc__(self, description, **kwargs): self._set_data(self.pad_data(data)) +@mosaic.tessera +class MeshedData(StructuredData, GriddedSaved): + pass + +@mosaic.tessera +class MeshedField(MeshedData): + pass + @mosaic.tessera class Scalar(StructuredData): diff --git a/stride/problem/domain.py b/stride/problem/domain.py index 8e7c241..7d11d17 100644 --- a/stride/problem/domain.py +++ b/stride/problem/domain.py @@ -3,7 +3,7 @@ from cached_property import cached_property -__all__ = ['Space', 'Time', 'SlowTime', 'Grid'] +__all__ = ['Space', 'MeshedSpace', 'Time', 'SlowTime', 'Grid'] class Space: @@ -259,6 +259,8 @@ def extended_grid(self): for dim in range(self.dim)] return tuple(axes) +class MeshedSpace: + pass class Time: """ From 135d49663df3b8cc5d91e6085ca6ba034a7ada6c Mon Sep 17 00:00:00 2001 From: Andrei Danila Date: Mon, 3 Aug 2026 10:24:06 +0100 Subject: [PATCH 02/15] minor correction to skeleton --- FORWARD_ADJOINT_REFERENCE.md | 603 ----------------------------------- stride/problem/base.py | 2 +- 2 files changed, 1 insertion(+), 604 deletions(-) delete mode 100644 FORWARD_ADJOINT_REFERENCE.md diff --git a/FORWARD_ADJOINT_REFERENCE.md b/FORWARD_ADJOINT_REFERENCE.md deleted file mode 100644 index 29e64e2..0000000 --- a/FORWARD_ADJOINT_REFERENCE.md +++ /dev/null @@ -1,603 +0,0 @@ -# Stride Forward & Inverse Scripts — Object & Argument Reference - -A practical reference for the objects most used in Stride forward-modelling and -inversion (FWI) scripts, centred on **the arguments accepted by `forward()` and -`adjoint()`** and by the PDE operators they drive. - -Line references point at the public repo -(`/Users/andreidanila/code/sonalis/stride`). A final section documents the -extra functionality and advanced examples in the private fork -(`/Users/andreidanila/code/sonalis/stride-private`). - ---- - -## 1. The two canonical script shapes - -Every script is an `async def main(runtime)` launched with `mosaic.run(main)` -(run on a cluster with `mrun -nw python script.py`). - -**Forward script** ([`breast2D/01_script_forward.py`](stride_examples/examples/breast2D/01_script_forward.py)): - -```python -space = Space(shape=(356, 385), extra=(50, 50), absorbing=(40, 40), spacing=0.5e-3) -time = Time(start=0., step=0.08e-6, num=2500) -problem = Problem(name='anastasio2D', space=space, time=time) - -vp = ScalarField(name='vp', grid=problem.grid); vp.load('...TrueModel.h5') -problem.medium.add(vp) -problem.transducers.default() -problem.geometry.default('elliptical', 128) -problem.acquisitions.default() -for shot in problem.acquisitions.shots: - shot.wavelets.data[0, :] = wavelets.tone_burst(0.5e6, 3, time.num, time.step) - -pde = IsoAcousticDevito.remote(grid=problem.grid, len=runtime.num_workers) -await forward(problem, pde, vp) # <-- entry point 1 -``` - -**Inverse script** ([`breast2D/02_script_inverse.py`](stride_examples/examples/breast2D/02_script_inverse.py)): - -```python -vp = ScalarField.parameter(name='vp', grid=problem.grid, needs_grad=True) -vp.fill(1500.) -problem.medium.add(vp) -problem.acquisitions.load(path=problem.output_folder, project_name=problem.name, version=0) - -pde = IsoAcousticDevito.remote(grid=problem.grid, len=runtime.num_workers) -loss = L2DistanceLoss.remote(len=runtime.num_workers) -optimiser = GradientDescent(vp, step_size=10, - process_grad=ProcessGlobalGradient(), - process_model=ProcessModelIteration(min=1400., max=1700.)) -optimisation_loop = OptimisationLoop() - -for block, freq in optimisation_loop.blocks(num_blocks, max_freqs): - await adjoint(problem, pde, loss, optimisation_loop, optimiser, vp, # <-- entry point 2 - num_iters=8, select_shots=dict(num=16, randomly=True), - f_max=freq, max_freqs=max_freqs) -``` - -Two key rules that explain most of the API: - -- **`.remote(...)`** turns an `Operator`/PDE class into a distributed *tessera* - spread across `len=runtime.num_workers` workers. Construction kwargs - (`grid`, `space`, `time`, `name`, `cached_operator`, …) go here. -- **`needs_grad=True`** on a variable (via `ScalarField.parameter(...)`) is what - switches the whole pipeline from pure forward to gradient-tracking adjoint - mode. `forward()` will automatically save the wavefield when it sees a - `needs_grad` input; `adjoint()` builds the adjoint graph from it. - ---- - -## 2. `forward(problem, pde, *args, **kwargs)` - -Defined in [`stride/__init__.py:53`](stride/__init__.py#L53). Runs `pde` once -per shot, stores the result in `shot.observed`, and (by default) appends it to -the observed HDF5 file. `*args` (e.g. `vp`, `rho`, `alpha`) and all extra -`kwargs` are **forwarded to the PDE**. - -| kwarg | type | default | meaning | -|---|---|---|---| -| `dump` | bool | `True` | Write forward result to disk (`observed` file). When `True` it also pre-loads any existing observed (version 0) so already-run shots are skipped. | -| `shot_ids` | list / int | `None` | Specific shots to run. `None` → all *remaining* shots (those with no observed yet). If none remain it logs a warning and returns. | -| `deallocate` | bool | `False` | Free `shot.observed` right after each shot (memory saving). | -| `safe` | bool | `True` | Discard workers that fail mid-execution instead of aborting. | -| `platform` | str | `'cpu'` | `'cpu'`, `'gpu'`, `'nvidia-acc'`, `'nvidia-cuda'`. Triggers GPU device round-robin across workers. | -| `devices` | list | `None` | Explicit GPU device ids; `None` → auto (`gpu_count()`). | -| `*args` | — | — | Positional PDE inputs (medium fields): `vp`, or `vp, vs, rho` (elastic), etc. Published to all workers. | -| `**kwargs` | — | — | Everything else is passed straight to the PDE `forward` (see §4). | - -Notable behaviour: the loop hands each `shot_id` to a worker, builds a -`sub_problem = problem.sub_problem(shot_id)`, calls -`pde(wavelets, *args, problem=sub_problem, runtime=worker, **kwargs).result()`, -writes `shot.observed.data[:] = traces.data`, and raises if NaN/Inf appear. - ---- - -## 3. `adjoint(problem, pde, loss, optimisation_loop, optimiser, *args, **kwargs)` - -Defined in [`stride/__init__.py:172`](stride/__init__.py#L172). This is the FWI -iteration driver. Per iteration it: selects shots → pre-processes -wavelets/observed → runs the PDE forward (saving the wavefield) → -post-processes modelled+observed traces → evaluates `loss` → runs the adjoint -(`fun.adjoint(...)`) to accumulate the gradient → takes an optimiser step → -dumps the updated variable. - -| kwarg | type | default | meaning | -|---|---|---|---| -| `num_iters` | int | `1` | Iterations to run **within the current block**. | -| `select_shots` | dict | `{}` | Shot-selection rules per iteration, forwarded to `Acquisitions.select_shot_ids` — e.g. `dict(num=16, randomly=True)`, or `dict(start=, end=, every=)`. | -| `lazy_loading` | bool | `False` | Load shot data each iteration and deallocate after, to save memory. | -| `dump` | bool | `True` | Save the updated optimiser variable after each iteration (enables restart). | -| `safe` | bool | `True` | Discard failing workers. | -| `f_min` | float | `None` | High-pass corner for filtering wavelets/traces. `None` → no high-pass. | -| `f_max` | float | `None` | Low-pass corner (frequency continuation). `None` → no low-pass. Usually set per block. | -| `filter_traces` | bool | `True` | Whether the trace-processing pipeline filters modelled/observed. | -| `filter_wavelets` | bool | `= filter_traces` | Whether wavelets/observed are band-limited before the PDE. | -| `filter_wavelets_relaxation` | float | `0.75` | Filter roll-off relaxation for wavelets. | -| `filter_traces_relaxation` | float | `0.75` (or `1.0` if no wavelet filtering) | Filter roll-off relaxation for traces. | -| `step_size` | float / `LineSearch` | `optimiser.step_size` | Step length; a `LineSearch` instance triggers the residual-keeping test-step loop. | -| `platform` / `devices` | str / list | `'cpu'` / `None` | GPU control, same as `forward()`. | -| `*args` | — | — | The medium variable(s) being inverted (also the `wrt` of the gradient), e.g. `vp` (must be a `.parameter(..., needs_grad=True)`). | -| `**kwargs` | — | — | Passed to the four processing pipelines, the PDE, the loss, and `optimiser.step`. Includes everything in §4 (`kernel`, `interpolation_type`, `boundary_type`, `devito_config`, …) and pipeline knobs (§7). | - -The `max_freqs=[...]` kwarg seen in examples is not consumed by `adjoint()` -itself; it is passed through so the pipelines/optimiser see the full schedule. - -**Frequency continuation** is driven outside `adjoint()` by the loop: -```python -for block, freq in optimisation_loop.blocks(num_blocks, max_freqs): - await adjoint(..., f_max=freq, max_freqs=max_freqs) -``` - ---- - -## 4. The PDE operator — `IsoAcousticDevito` - -The single most important object. It is the second-order isotropic **acoustic** -wave equation on Devito. Defined in -[`iso_acoustic/devito.py:24`](stride/physics/iso_acoustic/devito.py#L24). - -### 4.1 Construction (`.remote(...)`) - -```python -pde = IsoAcousticDevito.remote(grid=problem.grid, len=runtime.num_workers) -# or: IsoAcousticDevito.remote(space=space, time=time) -``` - -| construction kwarg | meaning | -|---|---| -| `grid` | Existing `Grid` (preferred). Alternatively pass `space=`/`time=`/`slow_time=`. | -| `len` | Number of worker replicas (`runtime.num_workers`). | -| `name` | Optional PDE name. | -| `cached_operator` | Reuse a compiled operator/grid stored in the worker warehouse across PDE instances. | -| `dev_grid` | Supply an existing `GridDevito` (advanced/shared setups). | - -Class-level: `space_order = 10`, `time_order = 2`. - -### 4.2 Forward-call arguments (the crux) - -Whether via `forward(problem, pde, vp, **kwargs)` or a direct -`await pde(wavelets, vp, problem=sub_problem, **kwargs).result()`, these are the -inputs and options consumed by `before_forward` / `run_forward` / `after_forward` -([`iso_acoustic/devito.py:220`](stride/physics/iso_acoustic/devito.py#L220)): - -**Positional / field inputs** - -| arg | type | default | meaning | -|---|---|---|---| -| `wavelets` | `Traces` | required | Source wavelets (`shot.wavelets`). Supplied automatically by `forward()`/`adjoint()`. | -| `vp` | `ScalarField` | required | Compressional speed of sound, m/s. | -| `rho` | `ScalarField` | `None` | Density, kg/m³. `None` → homogeneous. | -| `alpha` | `ScalarField` | `None` | Attenuation, dB/cm. `None` → lossless. | -| `problem` | `Problem`/`SubProblem` | required | The sub-problem (one shot). Injected by the driver. | - -**Physics / discretisation options** - -| kwarg | type | default | meaning | -|---|---|---|---| -| `kernel` | str | auto (`'OT2'`/`'OT4'`) | Time-stepping order: `'OT2'` (2nd) or `'OT4'` (4th). Auto-chosen from `dt` if unset. | -| `boundary_type` | str | `'sponge_boundary_2'` | Absorbing boundary. Also `'complex_frequency_shift_PML_2'` (lower OT4 stability). | -| `interpolation_type` | str | `'linear'` | Source/receiver interpolation: `'linear'` (bi/tri-linear) or `'hicks'` (sinc). | -| `attenuation_power` | int / None | `0` | Power of the attenuation law when `alpha` given (`0`, `2`, or `None`). | -| `drp` | bool | `False` | Dispersion-relation-preserving coefficients (build-dependent). | -| `diff_source` | bool | `False` | Inject the source as its 1st time-derivative instead of as-is. | -| `adaptive_boxes` | bool | `False` | Adaptive computational boxes (DevitoPRO). | -| `local_prec` | bool | `True` | Local preconditioning (build-dependent). | - -**Wavefield saving / gradient options** - -| kwarg | type | default | meaning | -|---|---|---|---| -| `save_wavefield` | bool | auto | Save forward wavefield for the gradient. Auto-`True` when any input `needs_grad`. | -| `time_bounds` | (int,int) | `(0, time.extended_num)` | Timestep window over which the wavefield is saved. | -| `save_undersampling` | int | auto (bandwidth) | Temporal undersampling factor when saving the wavefield. | -| `save_compression` | str | `None` (2D) / `'bitcomp'` (3D) | Wavefield compression (DevitoPRO/GPU only). | -| `save_interpolation` | bool | build-dependent | Cubic-spline interpolation of the saved wavefield. | -| `stream_wavefield` | bool / str | `True` | Streaming layer strategy: `True`/`False` or `'disk'`, `'host'`, `'device'`, `'disk-host'`, `'host-device'`, `'disk-host-device'`, `'no-layers'`. | -| `spill_wavefield` | bool | `False` | Spill the wavefield to disk/host. | -| `cache_forward` | bool | `False` | Cache the forward wavefield (in memory or `cache_location`) for reuse in the adjoint. | -| `cache_location` | str | `None` | Directory to cache the wavefield to disk. | -| `nbits_compression` | int | `9` | Bits for compressed streaming (maps to `devito_args['nbits']`). | - -**Wavefield dumping (debug / imaging)** - -| kwarg | type | default | meaning | -|---|---|---|---| -| `dump_forward_wavefield` | bool / int | `False` | Dump forward wavefield. `True` → every `save_undersampling`; int → every N steps. | -| `dump_adjoint_wavefield` | bool / int | `False` | Same for the adjoint wavefield. | -| `dump_wavefield_id` | int | shot id | Only dump this shot's wavefields. | - -**Platform / Devito plumbing** - -| kwarg | type | default | meaning | -|---|---|---|---| -| `platform` | str | `None`/`'cpu'` | `None`/`'cpu'` or `'nvidia-acc'` (OpenACC) / `'nvidia-cuda'`. | -| `devito_config` | dict | `{}` | Devito config applied **before** operator generation, e.g. `{'opt': ('advanced', {'index-mode': 'int64'})}`, `{'compiler':'pgcc','language':'openacc','platform':'nvidiaX'}`. `platform='nvidia-acc'` is shorthand for the latter. | -| `devito_args` | dict | `{}` | Args passed when **calling** the compiled operator, e.g. `{'autotune': 'off'}`, `{'deviceid': N}`, `{'nbits': 9}`. | -| `deallocate` | bool | `False` | Free Devito buffers (boundary, `p`, `src`, `rec`, `vp`, …) after the run. | - -### 4.3 Adjoint-call arguments - -The adjoint side (`before_adjoint`/`run_adjoint`/`after_adjoint`, -[`devito.py:689`](stride/physics/iso_acoustic/devito.py#L689)) is invoked by -`fun.adjoint(**kwargs)` inside the `adjoint()` driver — you rarely call it -directly. Its inputs are `(adjoint_source, wavelets, vp, rho=None, alpha=None, -**kwargs)` where `adjoint_source` is produced by the loss. It honours -`dump_adjoint_wavefield`, `dump_wavefield_id`, `cache_forward`, `time_bounds`, -`platform`, `deallocate`, and the same physics kwargs, then returns the -gradients via the `get_grad_*` methods. - -Gradients are produced per variable through the naming convention -(`prepare_grad_vp` / `init_grad_vp` / `get_grad_vp`), so `vp`, `rho`, and -`wavelets` can each be inverted for when flagged `needs_grad` -(see [`problem_type.py:202`](stride/physics/problem_type.py#L202)). - -### 4.4 Calling the PDE directly - -`forward()`/`adjoint()` are conveniences. You can call the tessera yourself, -which is how the homogeneous-medium test sweeps work -([`homogeneous_acoustic/forward_2D.py:109`](stride_examples/examples/homogeneous_acoustic/forward_2D.py#L109)): - -```python -sub_problem = problem.sub_problem(shot.id) -traces = await pde(sub_problem.shot.wavelets, vp, - problem=sub_problem, diff_source=True, - kernel='OT4', interpolation_type='hicks', - boundary_type='complex_frequency_shift_PML_2', - rho=rho, alpha=alpha, attenuation_power=2).result() -data = traces.data # numpy array, shape (num_receivers, time.num) -await pde.clear_operators() # force recompile when config changes -``` - -`forward()` returns nothing (it writes `shot.observed`); a direct call returns a -`Traces` object whose `.data` is the modelled gather. - -### 4.5 Other PDE operators (same call convention) - -- **`IsoElasticDevito`** ([`iso_elastic/devito.py:16`](stride/physics/iso_elastic/devito.py#L16)) — - stress-strain elastic wave equation. Forward inputs are - `(wavelets, vp, vs, rho, problem=...)`; `space_order=10`, `time_order=1`, - `boundary_type='sponge_boundary_1'`. Multi-parameter: pass `vp, vs, rho=rho`. -- **`MarmottantDevito`** ([`marmottant/devito.py`](stride/physics/marmottant/devito.py)) — - microbubble model. - ---- - -## 5. Problem-definition objects - -All in `stride/problem/`. These build the `problem` you hand to -`forward()`/`adjoint()`. - -### `Space` — [`domain.py:9`](stride/problem/domain.py#L9) -Spatial grid = inner domain + padding. -`Space(shape, spacing, extra, absorbing)`. -- `shape` (tuple) inner grid points; `spacing` (tuple or scalar float, m); - `extra` padding points per axis; `absorbing` portion of `extra` used for the - boundary. Exposes `.dim`, `.limit` (physical size), `.extended_shape`, - `.inner`, `.resample(...)`. - -### `Time` — [`domain.py:263`](stride/problem/domain.py#L263) -`Time(start, step, num, stop)` — give any three. `num` must be `int`. Exposes -`.start/.step/.num/.stop`, `.extended_num`, `.extend(...)`, `.resample(...)`. - -### `SlowTime` — [`domain.py:417`](stride/problem/domain.py#L417) -Frame/acquisition sampling for multi-frame data: -`SlowTime(frame_rate|frame_step, acq_rate|acq_step, num_frame, num_acq)`. - -### `Grid` — [`domain.py:528`](stride/problem/domain.py#L528) -`Grid(space, time, slow_time)` — a bundle. `problem.grid` is usually what you -pass to `.remote(grid=...)` and to fields. - -### `Problem` — [`problem.py:13`](stride/problem/problem.py#L13) -`Problem(name, space=, time=)` (or `grid=`). Top-level container. Attributes -auto-created: `.medium`, `.transducers`, `.geometry`, `.acquisitions`, `.grid`. -- `.sub_problem(shot_id)` → a `SubProblem` with `.shot`, `.shot_id` (built per - shot inside the drivers). -- `.plot()`, `.load(...)`, `.dump(...)`, `.output_folder`/`.input_folder` - (default `cwd`). -- `.space_resample(new_spacing)`, `.time_resample(new_step, new_num)` for - multi-resolution FWI. - -### `ScalarField` — [`data.py:790`](stride/problem/data.py#L790) -The medium-field type (`vp`, `rho`, `alpha`). Two construction paths: -- Forward / known model: `ScalarField(name='vp', grid=problem.grid)` then - `.load(path)` or `.fill(value)`. -- Inversion variable: `ScalarField.parameter(name='vp', grid=problem.grid, - needs_grad=True)`. `.parameter()` is injected by `@mosaic.tessera` and returns - the optimisable/distributed variant; `needs_grad=True` enables gradients. -- `time_dependent=`/`slow_time_dependent=` prepend time axes. -- Members: `.data` (inner view), `.extended_data` (with padding), `.fill(v)`, - `.plot()`, `.load()`/`.dump()`, `.needs_grad`, `.grad`, `.clear_grad()`. - -### `Traces` — [`data.py:1383`](stride/problem/data.py#L1383) -Time traces indexed by transducer id — the type of `shot.wavelets` and -`shot.observed`. Write with `shot.wavelets.data[i, :] = ...`. `.data` shape is -`(num_transducers, time.num)`. `.plot(plot_type='gather'|'spectrum')`, -`.get(id)`, `.alike(...)`. `DiskTraces` is the lazy on-disk variant. - -### `Medium` — [`medium.py:8`](stride/problem/medium.py#L8) -Named-field container. `problem.medium.add(vp)`; access `medium.vp` / -`medium['vp']`. `.load()/.dump()/.plot()` iterate all fields. - -### `Transducers` — [`transducers.py:11`](stride/problem/transducers.py#L11) -Registry of transducer devices. `problem.transducers.default()` creates one -`PointTransducer(0)`. - -### `Geometry` — [`geometry.py:106`](stride/problem/geometry.py#L106) -Transducer *locations*. `problem.geometry.default('elliptical', num_locations)` -(2D) auto-computes radius/centre from `space.limit`; -`'ellipsoidal', num_locations, radius, centre, theta=, threshold=` (3D). -`.coordinates`, `.locations`, `.num_locations`, `.plot()`. - -### `Acquisitions` — [`acquisitions.py:733`](stride/problem/acquisitions.py#L733) -The set of shots (the data). -- `.default()` — one single-source shot per location, all locations as receivers. -- `.load(path=, project_name=, version=0, shot_ids=None, fast=False)`. -- `.select_shot_ids(num=, start=, end=, every=1, randomly=False)` — stateful - per-iteration selection (this is what `select_shots` in `adjoint()` feeds). -- `.remaining_shot_ids` — shots with no observed yet (what `forward()` runs). -- `.shots`, `.shot_ids`, `.num_shots`, `.plot()`, `.reset_selection()`. - -### `Shot` — [`acquisitions.py:57`](stride/problem/acquisitions.py#L57) -One acquisition event: `.wavelets` (per source), `.observed` (per receiver), -`.delays`, `.source_coordinates`, `.receiver_coordinates`, -`.num_sources`/`.num_receivers`. - -### Wavelet helpers — [`utils/wavelets.py`](stride/utils/wavelets.py) -- `tone_burst(centre_freq, n_cycles, n_samples, dt, envelope='gaussian')` -- `ricker(centre_freq, n_samples, dt)` -- `continuous_wave(centre_freq, n_samples, dt, ramp_length=4, phase=0)` - -### Asset fetch — [`utils/fetch.py`](stride/utils/fetch.py) -`fetch('anastasio2D', dest='data/...h5')` downloads a known release asset once. - ---- - -## 6. Optimisation objects (inverse scripts) - -### `L2DistanceLoss` — [`loss/l2_distance.py:14`](stride/optimisation/loss/l2_distance.py#L14) -`f = ½‖modelled − observed‖²`. `L2DistanceLoss.remote(len=runtime.num_workers)`. -- Ctor kwarg `d_sample=4` (downsampling of the residual passed to the functional). -- `forward(modelled, observed, **kwargs)` → `FunctionalValue`; consumes - `problem`/`shot_id`, forwards `keep_residual`. -- `adjoint(d_fun, modelled, observed)` → `(grad_modelled, grad_observed)` — the - adjoint source. - -### `GradientDescent` / `LocalOptimiser` — [`optimisers/`](stride/optimisation/optimisers/) -> There is **no `Adam`** in the public repo — only `GradientDescent` (and the -> `LocalOptimiser` base). Adam/NAdam/RAdam/SGD/AdaBelief/AGD live in the private -> fork (§8). - -`GradientDescent(variable, step_size=1., process_grad=..., process_model=...)`. -Base `LocalOptimiser` ctor kwargs: - -| kwarg | default | meaning | -|---|---|---| -| `variable` | — | Variable to optimise (must be `needs_grad`). | -| `step_size` | `1.` | Float or `LineSearch`. | -| `test_step_size` | `1.` | Multiplier on the processed gradient. | -| `force_step` | `False` | Skip step clipping/capping. | -| `max_step` | `None` | Cap on step magnitude. | -| `process_grad` | `ProcessGlobalGradient(**kwargs)` | Gradient pre-processing pipeline. | -| `process_model` | `ProcessModelIteration(**kwargs)` | Model post-processing pipeline. | -| `reset_block` | `False` | Reset optimiser state each block (flag, not method). | -| `reset_iteration` | `False` | Reset optimiser state each iteration. | -| `dump_grad` / `dump_prec` | `False` | Debug dumps. | - -Methods: `await step(step_size=None, grad=None, step_loop=..., **kwargs)`, -`clear_grad()`, `reset()`, `dump()/load()`. - -### `OptimisationLoop` / `Block` / `Iteration` — [`optimisation_loop.py`](stride/optimisation/optimisation_loop.py) -- `OptimisationLoop(name='optimisation_loop')`. Properties `.num_blocks`, - `.current_block`; `.blocks(num, *iters, restart=False, restart_id=-1)` is the - generator you loop over (zips extra sequences like `max_freqs`). -- `Block.iterations(num, *iters, ...)` yields `Iteration`s; `.total_loss`, `.id`. -- `Iteration`: `.add_loss(fun)`, `.add_submitted/.add_completed(shot)`, - `.next_run()` (line-search test steps), `.id`, `.abs_id`, `.total_loss`, - `.prev_run`. - -### Processing pipelines — [`pipelines/default_pipelines.py`](stride/optimisation/pipelines/default_pipelines.py) -`.remote(...)` operators the `adjoint()` driver builds automatically; you tune -them by passing kwargs through `adjoint()`. - -| pipeline | role | default steps (kwarg → default) | -|---|---|---| -| `ProcessWavelets` | pre-process source wavelets | `check_traces`(T), `filter_traces`(T), `shift_traces`, `resonance_filter`(F) | -| `ProcessObserved` | pre-process observed | `check_traces`, `filter_traces` | -| `ProcessWaveletsObserved` | joint step | `differentiate_traces`(T) | -| `ProcessTraces` | modelled+observed before the loss | `check_traces`(T), `filter_offsets`(F), `mute_first_arrival`(T), `mute_traces`(T), `filter_traces`(T), `agc`(F), `norm_per_shot`(T)/`norm_per_trace`(F), `scale_per_*`(F), `time_tweaking`(T), `time_weighting`(T) | -| `ProcessGlobalGradient` | default `process_grad` | `mask_field`(`mask_grad`=T), `smooth_field`(`smooth_grad`=T), `norm_field`(`norm_grad`=T) | -| `ProcessModelIteration` | default `process_model` | `clip` (uses `min=`/`max=`) | - -> Several `ProcessTraces` steps (`resonance_filter`, `differentiate_traces`, -> `filter_offsets`, `mute_first_arrival`, `agc`, `time_tweaking`, -> `time_weighting`) are added *non-raising*: they are no-ops in the public repo -> unless the step is registered — and those step classes ship in the **private -> fork** (§8). - -### Individual steps — [`pipelines/steps/`](stride/optimisation/pipelines/steps/) -Registered (public): `filter_traces`, `norm_per_shot`, `norm_per_trace`, -`scale_per_shot`, `scale_per_trace`, `norm_field`, `smooth_field`, `mask_field`, -`mute_traces`, `clip`, `check_traces`, `dump`, `shift_traces`. Common kwargs: - -| step | key kwargs | -|---|---| -| `FilterTraces` | `f_min`, `f_max`, `filter_type` (`'cos'`/`'butterworth'`/`'fir'`), `filter_relaxation` | -| `MuteTraces` | `f_max`, `filter_relaxation` | -| `Clip` | `min`, `max` | -| `SmoothField` | `smooth_sigma` (default 0.25 = 25% of a cell) | -| `NormField` | `global_norm` (F), `norm_guess_change` (0.5) | -| `MaskField` | `mask`, `mask_rampoff` (10) | -| `ScalePer*` | `scale_to`, `relative_scale` (T) | -| `CheckTraces` | `raise_incorrect` (T), `filter_incorrect` (F) | -| `ShiftTraces` | `f_max`, `filter_relaxation` | - -### `LineSearch` — [`step_length/line_search.py:8`](stride/optimisation/step_length/line_search.py#L8) -Abstract step-length base (`init_search`, `next_step`). Passed as -`step_size=LineSearch(...)`; the optimiser drives it with a `step_loop` callable. -The public repo ships only the abstract base — concrete searches are in the -private fork. - ---- - -## 7. Quick kwarg cheat-sheet - -Passing these through `forward(problem, pde, vp, ...)` / -`adjoint(problem, pde, loss, loop, opt, vp, ...)` reaches the right layer: - -```python -# --- performance / hardware --- -platform='nvidia-cuda' # or 'nvidia-acc', 'gpu', 'cpu' -devices=[0,1] # explicit GPUs -devito_config={'opt': ('advanced', {'index-mode': 'int64'})} -devito_args={'autotune': 'off'} -deallocate=True # free device buffers each shot - -# --- physics --- -kernel='OT4' # or 'OT2' -interpolation_type='hicks' # or 'linear', 'sinc' -boundary_type='complex_frequency_shift_PML_2' -rho=rho_field, alpha=alpha_field, attenuation_power=2 -diff_source=True - -# --- forward wavefield / imaging --- -save_wavefield=True, time_bounds=(0, N), save_undersampling=4 -dump_forward_wavefield=True, dump_adjoint_wavefield=8, dump_wavefield_id=0 -cache_forward=True, cache_location='/scratch' - -# --- inversion control (adjoint only) --- -num_iters=8 -select_shots=dict(num=16, randomly=True) -f_min=0.1e6, f_max=0.5e6, max_freqs=[0.3e6, 0.4e6, 0.5e6, 0.6e6] -lazy_loading=True -step_size=10 # or a LineSearch instance - -# --- pipeline toggles (adjoint only) --- -filter_traces=True, filter_wavelets=True -mask_grad=True, smooth_grad=True, norm_grad=True -smooth_sigma=0.5, norm_guess_change=0.25 -``` - ---- - -## 8. Private fork (`stride-private`) — extensions & advanced examples - -The private repo forks an **older** public baseline and adds specialised -operators, losses, optimisers and pipeline steps. (Its acoustic operator is -older/smaller than the current public one; where the two diverge, the private -defaults differ — e.g. `drp=True` and -`boundary_type='hybrid_interpolating_boundary_2'` by default.) - -### 8.1 Advanced invocation patterns (copy-paste worthy) - -**GPU multi-block acoustic FWI** — `examples/alpha2D/inverse.py:74`: -```python -await adjoint(problem, pde, loss, optimisation_loop, optimiser, vp, - num_iters=num_iters, - select_shots=dict(num=12, randomly=True), - f_max=freq, max_freqs=max_freqs, - kernel='OT4', fw3d_mode=True, - interpolation_type='hicks', platform='nvidia-cuda') -``` -matching forward `alpha2D/forward.py:62`: -```python -await forward(problem, pde, vp, kernel='OT4', fw3d_mode=True, - interpolation_type='hicks', platform='nvidia-cuda') -``` - -**Elastic forward, multi-parameter + Devito opt** — `examples/alpha2D_elastic/forward.py:90`: -```python -await forward(problem, pde, vp, vs, rho=rho, shot_ids=[0], - interpolation_type='sinc', dump=False, deallocate=False, - dump_forward_wavefield=False, - devito_config={'opt': ('advanced', {'index-mode': 'int64'})}) -``` - -**3D elastic, autotune off** — `examples/alpha3D_elastic/forward.py:113`: -```python -await forward(problem, pde, vp, rho=rho, shot_ids=[0], - interpolation_type='sinc', kernel='OT4', - dump=False, deallocate=False, devito_args={'autotune': 'off'}) -``` - -**Optimiser with line search** — `examples/alpha2D/inverse.py:58`: -```python -optimiser = GradientDescent(vp, step_size=LineSearch(), - process_grad=ProcessGlobalGradient(), - process_model=ProcessModelIteration(min=1450., max=3000.)) -``` - -**Custom driver for source/receiver-IR inversion (SRI)** — -`examples/SRI/SRI_transmit_optim_wavelet.py` replaces the library `adjoint()` -with its own `adjoint_finite(problem, pde, loss, optimisation_loop, optimiser, -*args, **kwargs)` that builds `ProcessWavelets.remote(f_min=, f_max=)` / -`ProcessTraces.remote(...)`, convolves the wavelet with each transducer's -`transmit_ir` via a `Convolution.remote(...)` operator, runs the PDE, and calls -`await fun.adjoint(**kwargs)`. Manual per-worker GPU assignment: -```python -devito_args = kwargs.get('devito_args', {}) -devito_args['deviceid'] = devices[worker.indices[1] % num_gpus] -kwargs['devito_args'] = devito_args -``` - -**Multi-resolution FWI** — `examples/alpha2D_resample/` uses -`problem.space_resample()` / `problem.time_resample()` between blocks. - -### 8.2 Added physics operators -- `IsoAcousticAnalyticNumba` (analytic Green's-function acoustic, Numba - CPU/CUDA) — extra kwargs `compute_forward`, `vp_constant`, `gradient_crop`, - `threads_per_block`, `blocks_per_grid`, `diff_source`, `save_undersampling`. -- `TransportDevito` / `SLTransportNumpy` (semi-Lagrangian transport) — kwargs - `interpolating`, `velocity_interpolation`, `taylor_accuracy`, `parallel_time`, - `num_procs`, `fill_holes`; forward inputs `(sigma_0, u, ...)`. -- `FlowInjectionDevito` — kwargs `interpolation_type`, `t_i`. -- New boundaries `HybridHigdonBoundary2`, `HybridInterpolatingBoundary2`. - -### 8.3 Added loss functions (`stride_private/optimisation/loss/`) -| loss | key kwargs | -|---|---| -| `AdaptiveWaveformLoss` (AWI) | `pad`, `augment`, `gamma`, `eta`, `mode='reverse'`, `type='standard'`, `d_sample=4` | -| `OptimalTransportLoss` (GSOT) | `p=2`, `mode='taot'`, `e_start/e_end/e_fac`, `taot_algorithm='dtw'`, `reg`, `max_iter`, `tol`, `d_sample` | -| `FrequencyControllableEnvelopeLoss` (FCEI) | `pad`, `p`, `filter_type='butterworth'`, `filter_order=8`, `num_modify`, `weight_l2`, `weight_fcei` | -| `ReverseTimeMigration` (RTM) | imaging condition | -| `DoubleDifferenceLoss` | double-difference misfit | -| `MultiLoss` | `weight_l2`, `weight_ot`, `weight_fcei`, `fcei_mode`, `mute_traces`, `mute_threshold`, `d_sample` | - -### 8.4 Added optimisers & steps -- Optimisers: `Adam`, `NAdam`, `RAdam` (`betas=(0.9,0.999)`, `eps=1e-8`, `t`), - `SGD`, `AdaBelief`, `AGD` — all on a `LocalOptimiserSaved` base. -- Pipeline steps (these are the ones the public default pipelines reference but - don't ship): `AGC`, `DifferentiateTraces`, `FilterOffsets`, `MuteFirstArrival`, - `ResonanceFilter`, `TimeTweaking`, `TimeWeighting`. -- Constraint `DivergenceFree`; helper operators `Convolution`, `Split`, - `Replicate`, `Unflatten`, `SVDFilter`, `DiffFilter`, `FreqWindowingFilter`. -- Fullwave3D interop (`.ttr/.pgy/.vtr`) in `utils/fullwave.py`. - ---- - -## 9. Mental model (how the pieces connect) - -``` -Space + Time ─► Grid ─► Problem ─┬─ medium (ScalarField vp/rho/alpha) - ├─ transducers (.default) - ├─ geometry (.default 'elliptical'/'ellipsoidal') - └─ acquisitions (Shots: wavelets + observed) - -pde = IsoAcousticDevito.remote(grid=..., len=workers) # the physics - -FORWARD: forward(problem, pde, vp, **pde_kwargs) # writes shot.observed - -INVERSE: loss = L2DistanceLoss.remote() - optimiser = GradientDescent(vp.parameter(needs_grad=True), - process_grad=ProcessGlobalGradient(), - process_model=ProcessModelIteration(min,max)) - loop = OptimisationLoop() - for block, freq in loop.blocks(num_blocks, max_freqs): - adjoint(problem, pde, loss, loop, optimiser, vp, - num_iters=, select_shots=, f_max=freq, **pde_kwargs) - └─ per iter: ProcessWavelets/Observed → pde.forward (saves wavefield) - → ProcessTraces → loss.forward → fun.adjoint (grad) - → optimiser.step → dump updated vp -``` - -*Generated as a code reference for Stride forward/inverse scripting.* diff --git a/stride/problem/base.py b/stride/problem/base.py index 52c3a09..0fa7239 100644 --- a/stride/problem/base.py +++ b/stride/problem/base.py @@ -4,7 +4,7 @@ from .domain import Space, Time, SlowTime, Grid -__all__ = ['Gridded', 'Meshed', 'Saved', 'GriddedSaved', 'MeshedSaved', 'ProblemBase'] +__all__ = ['Gridded', 'Saved', 'GriddedSaved', 'ProblemBase'] class Gridded: From 49a99006ea532da24f5b59e9b45baabab3c4756a Mon Sep 17 00:00:00 2001 From: Andrei Danila Date: Mon, 3 Aug 2026 12:35:06 +0100 Subject: [PATCH 03/15] tests --- stride/tests/__init__.py | 0 stride/tests/conftest.py | 188 ++++++++++++ stride/tests/test_meshed_data.py | 296 +++++++++++++++++++ stride/tests/test_meshed_medium.py | 339 ++++++++++++++++++++++ stride/tests/test_meshed_serialisation.py | 270 +++++++++++++++++ stride/tests/test_meshed_space.py | 236 +++++++++++++++ 6 files changed, 1329 insertions(+) create mode 100644 stride/tests/__init__.py create mode 100644 stride/tests/conftest.py create mode 100644 stride/tests/test_meshed_data.py create mode 100644 stride/tests/test_meshed_medium.py create mode 100644 stride/tests/test_meshed_serialisation.py create mode 100644 stride/tests/test_meshed_space.py diff --git a/stride/tests/__init__.py b/stride/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/stride/tests/conftest.py b/stride/tests/conftest.py new file mode 100644 index 0000000..d4b73ad --- /dev/null +++ b/stride/tests/conftest.py @@ -0,0 +1,188 @@ +""" +Shared fixtures for the FEM/meshed tests. + +The reference mesh is a Kuhn (Freudenthal) tetrahedralisation of a structured +box: every hexahedral cell of an ``(nx, ny, nz)`` node grid is split into six +tetrahedra. This gives a genuine unstructured mesh (flat node list, explicit +connectivity) without needing DOLFINx or gmsh to be installed, which mirrors +what ``ae_modelling.fem.mesh.make_mesh`` produces via +``dolfinx.mesh.create_box``. +""" + +import numpy as np +import pytest + + +# Local corner index within a hexahedron is 4*i + 2*j + k, so the six tets of +# the Kuhn decomposition along the 0 -> 7 diagonal are: +KUHN_TETS = ( + (0, 1, 3, 7), + (0, 1, 5, 7), + (0, 2, 3, 7), + (0, 2, 6, 7), + (0, 4, 5, 7), + (0, 4, 6, 7), +) + + +def box_tetra_mesh(shape=(3, 3, 3), spacing=(1e-3, 1e-3, 1e-3), origin=(0., 0., 0.)): + """ + Build a tetrahedral mesh of an axis-aligned box. + + Parameters + ---------- + shape : tuple + Number of nodes per axis. + spacing : tuple + Node spacing per axis, in metres. + origin : tuple + Lower corner of the box, in metres. + + Returns + ------- + nodes : ndarray + ``(num_nodes, 3)`` float64 node coordinates. + cells : ndarray + ``(num_cells, 4)`` int32 node indices, one row per tetrahedron. + + """ + axes = [np.arange(n) * d + o for n, d, o in zip(shape, spacing, origin)] + mesh = np.meshgrid(*axes, indexing='ij') + nodes = np.stack([each.ravel() for each in mesh], axis=-1).astype(np.float64) + + strides = (shape[1] * shape[2], shape[2], 1) + + cells = [] + for i in range(shape[0] - 1): + for j in range(shape[1] - 1): + for k in range(shape[2] - 1): + corners = [ + (i + ((local >> 2) & 1)) * strides[0] + + (j + ((local >> 1) & 1)) * strides[1] + + (k + (local & 1)) * strides[2] + for local in range(8) + ] + for tet in KUHN_TETS: + cells.append([corners[each] for each in tet]) + + return nodes, np.asarray(cells, dtype=np.int32) + + +def triangle_mesh(shape=(3, 3), spacing=(1e-3, 1e-3), origin=(0., 0.)): + """ + Build a 2D triangular mesh of an axis-aligned rectangle. + + Returns + ------- + nodes : ndarray + ``(num_nodes, 2)`` float64 node coordinates. + cells : ndarray + ``(num_cells, 3)`` int32 node indices, one row per triangle. + + """ + axes = [np.arange(n) * d + o for n, d, o in zip(shape, spacing, origin)] + mesh = np.meshgrid(*axes, indexing='ij') + nodes = np.stack([each.ravel() for each in mesh], axis=-1).astype(np.float64) + + cells = [] + for i in range(shape[0] - 1): + for j in range(shape[1] - 1): + bottom_left = i * shape[1] + j + bottom_right = bottom_left + 1 + top_left = bottom_left + shape[1] + top_right = top_left + 1 + cells.append([bottom_left, bottom_right, top_right]) + cells.append([bottom_left, top_right, top_left]) + + return nodes, np.asarray(cells, dtype=np.int32) + + +@pytest.fixture +def tetra_mesh(): + """Nodes and cells of a 3x3x3-node, 1 mm box (27 nodes, 48 tets).""" + return box_tetra_mesh(shape=(3, 3, 3), spacing=(1e-3, 1e-3, 1e-3)) + + +@pytest.fixture +def tri_mesh(): + """Nodes and cells of a 3x3-node, 1 mm rectangle (9 nodes, 8 triangles).""" + return triangle_mesh(shape=(3, 3), spacing=(1e-3, 1e-3)) + + +@pytest.fixture +def meshed_space(tetra_mesh): + """A 3D MeshedSpace over the reference tetrahedral mesh.""" + from stride.problem.domain import MeshedSpace + + nodes, cells = tetra_mesh + return MeshedSpace(nodes=nodes, cells=cells) + + +@pytest.fixture +def meshed_space_2d(tri_mesh): + """A 2D MeshedSpace over the reference triangular mesh.""" + from stride.problem.domain import MeshedSpace + + nodes, cells = tri_mesh + return MeshedSpace(nodes=nodes, cells=cells) + + +@pytest.fixture +def tagged_meshed_space(tetra_mesh): + """ + A MeshedSpace whose cells carry two material tags. + + Cells in the lower half of the box (in z) are tagged ``1``, the rest ``2``, + mirroring the ``cell_tags`` that ``ae_modelling`` reads out of a gmsh file. + """ + from stride.problem.domain import MeshedSpace + + nodes, cells = tetra_mesh + centroids = nodes[cells].mean(axis=1) + cell_tags = np.where(centroids[:, 2] < 1e-3, 1, 2).astype(np.int32) + + return MeshedSpace(nodes=nodes, cells=cells, cell_tags=cell_tags) + + +@pytest.fixture +def structured_space(): + """A conventional structured Space, used for the serialisation regressions.""" + from stride.problem.domain import Space + + return Space(shape=(6, 8), spacing=(1e-3, 1e-3), extra=(2, 2), absorbing=(1, 1)) + + +@pytest.fixture +def tissue_properties(): + """ + Label -> (conductivity, relative permittivity) map. + + Values are the defaults from ``ae_modelling.tissue.properties``. + + """ + return { + 0: {'name': 'background', 'sigma': 1e-6, 'eps_r': 1e0}, + 1: {'name': 'grey_matter', 'sigma': 1.52e-1, 'eps_r': 2.19e3}, + 2: {'name': 'white_matter', 'sigma': 9.47e-2, 'eps_r': 7.12e2}, + 3: {'name': 'csf', 'sigma': 2e0, 'eps_r': 1.09e2}, + 4: {'name': 'skull', 'sigma': 2.22e-2, 'eps_r': 1.75e2}, + 5: {'name': 'skin', 'sigma': 4.36e-3, 'eps_r': 1.06e3}, + } + + +@pytest.fixture +def label_volume(): + """ + A small voxel label volume plus its world affine. + + Voxels are 1 mm isotropic with the volume origin at the world origin, so + world coordinate ``x`` maps to voxel index ``round(x / 1e-3)``. Labels vary + along z only: the lower half is ``1``, the upper half is ``3``. + """ + volume = np.zeros((4, 4, 4), dtype=np.int64) + volume[:, :, :2] = 1 + volume[:, :, 2:] = 3 + + affine = np.diag([1e-3, 1e-3, 1e-3, 1.]) + + return volume, affine diff --git a/stride/tests/test_meshed_data.py b/stride/tests/test_meshed_data.py new file mode 100644 index 0000000..6cb5c74 --- /dev/null +++ b/stride/tests/test_meshed_data.py @@ -0,0 +1,296 @@ +""" +Tests for MeshedData and MeshedField (stride/problem/data.py). + +MeshedData is to MeshedSpace what StructuredData is to Space: a buffer with an +explicit shape and the arithmetic/gradient machinery. MeshedField is to +MeshedData what ScalarField is to StructuredData: it derives its shape from the +grid, and optionally prepends time and slow-time axes. + +Because a mesh has no padding, the buffer is flat and ``extended_shape == +shape``, exactly as for SparseField. That is the existing class these two +follow most closely. + +Contract under test: + +- ``MeshedData`` takes an explicit ``shape`` or ``data``, like StructuredData +- ``MeshedField`` derives ``(num_nodes,)`` from ``space``, with the same + ``time_dependent`` / ``slow_time_dependent`` / ``dim`` options as SparseField +- padding is a no-op on both +- ``alike`` / ``copy`` / ``detach`` carry the meshed grid across +- inherited arithmetic operates on the nodal buffer +- ``clear_grad`` allocates a nodal gradient +""" + +import numpy as np +import pytest + +from stride.problem.domain import Grid + + +def nodal_grid(space): + return Grid(space, None, None) + + +class TestMeshedDataShape: + + def test_explicit_shape(self, meshed_space): + from stride.problem.data import MeshedData + + data = MeshedData(name='raw', shape=(27,), grid=nodal_grid(meshed_space)) + + assert tuple(data.shape) == (27,) + assert tuple(data.extended_shape) == (27,) + + def test_shape_inferred_from_space(self, meshed_space): + from stride.problem.data import MeshedData + + data = MeshedData(name='raw', grid=nodal_grid(meshed_space)) + + assert tuple(data.shape) == (meshed_space.num_nodes,) + + def test_shape_inferred_from_data(self, meshed_space): + from stride.problem.data import MeshedData + + values = np.arange(27, dtype=np.float32) + data = MeshedData(name='raw', data=values, grid=nodal_grid(meshed_space)) + + assert tuple(data.shape) == (27,) + np.testing.assert_array_equal(data.data, values) + + def test_mismatched_data_length_rejected(self, meshed_space): + from stride.problem.data import MeshedData + + # A nodal field with the wrong number of values cannot be interpolated + # onto the mesh, so this has to fail loudly rather than at solve time. + with pytest.raises(ValueError): + MeshedData(name='raw', data=np.zeros(26, dtype=np.float32), + grid=nodal_grid(meshed_space)) + + def test_default_dtype_is_float32(self, meshed_space): + from stride.problem.data import MeshedData + + data = MeshedData(name='raw', grid=nodal_grid(meshed_space)) + + assert data.dtype == np.float32 + + def test_complex_dtype_is_honoured(self, meshed_space): + from stride.problem.data import MeshedData + + data = MeshedData(name='raw', dtype=np.complex128, grid=nodal_grid(meshed_space)) + data.fill(1 + 2j) + + assert data.data.dtype == np.complex128 + np.testing.assert_allclose(data.data, 1 + 2j) + + +class TestMeshedDataPadding: + + def test_inner_is_the_whole_buffer(self, meshed_space): + from stride.problem.data import MeshedData + + data = MeshedData(name='raw', grid=nodal_grid(meshed_space)) + + assert data.inner == (slice(0, None),) + + def test_pad_data_is_a_noop(self, meshed_space): + from stride.problem.data import MeshedData + + data = MeshedData(name='raw', grid=nodal_grid(meshed_space)) + values = np.arange(27, dtype=np.float32) + + np.testing.assert_array_equal(data.pad_data(values), values) + + def test_data_and_extended_data_agree(self, meshed_space): + from stride.problem.data import MeshedData + + data = MeshedData(name='raw', grid=nodal_grid(meshed_space)) + data.fill(3.) + + np.testing.assert_array_equal(data.data, data.extended_data) + + +class TestMeshedFieldShape: + + def test_scalar_nodal_field(self, meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=nodal_grid(meshed_space)) + + assert tuple(field.shape) == (27,) + assert tuple(field.extended_shape) == (27,) + assert field.inner == (slice(0, None),) + + def test_vector_nodal_field(self, meshed_space): + from stride.problem.data import MeshedField + + # The electric field E = -grad(phi) is the motivating case. + field = MeshedField(name='e_field', dim=3, grid=nodal_grid(meshed_space)) + + assert tuple(field.shape) == (27, 3) + assert field.inner == (slice(0, None), slice(0, None)) + + def test_time_dependent_field(self, meshed_space): + from stride.problem.data import MeshedField + from stride.problem.domain import Time + + time = Time(start=0., step=1e-6, num=11) + field = MeshedField(name='phi', time_dependent=True, + grid=Grid(meshed_space, time, None)) + + assert tuple(field.shape) == (11, 27) + + def test_time_dependent_vector_field(self, meshed_space): + from stride.problem.data import MeshedField + from stride.problem.domain import Time + + time = Time(start=0., step=1e-6, num=11) + field = MeshedField(name='e_field', dim=3, time_dependent=True, + grid=Grid(meshed_space, time, None)) + + assert tuple(field.shape) == (11, 27, 3) + + def test_dim_defaults_to_scalar(self, meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=nodal_grid(meshed_space)) + + assert field.dim == 1 + + def test_num_nodes_is_exposed(self, meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=nodal_grid(meshed_space)) + + assert field.num_nodes == meshed_space.num_nodes + + def test_explicit_shape_overrides_the_grid(self, meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', shape=(5,), grid=nodal_grid(meshed_space)) + + assert tuple(field.shape) == (5,) + + +class TestMeshedFieldCopying: + + def test_alike_keeps_the_meshed_grid(self, meshed_space): + from stride.problem.data import MeshedField + from stride.problem.domain import MeshedSpace + + field = MeshedField(name='sigma', grid=nodal_grid(meshed_space)) + other = field.alike(name='eps') + + assert isinstance(other.space, MeshedSpace) + assert other.space is meshed_space + assert tuple(other.shape) == tuple(field.shape) + assert other.dtype == field.dtype + + def test_copy_duplicates_the_buffer(self, meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=nodal_grid(meshed_space)) + field.fill(2.) + + cpy = field.copy() + cpy.data[:] = 5. + + np.testing.assert_allclose(field.data, 2.) + np.testing.assert_allclose(cpy.data, 5.) + + def test_copy_keeps_vector_shape(self, meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField(name='e_field', dim=3, grid=nodal_grid(meshed_space)) + field.fill(1.) + + assert tuple(field.copy().shape) == (27, 3) + + def test_detach_keeps_shape_and_data(self, meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=nodal_grid(meshed_space)) + field.fill(4.) + + detached = field.detach() + + assert tuple(detached.shape) == (27,) + np.testing.assert_allclose(detached.data, 4.) + + +class TestMeshedFieldArithmetic: + + def test_add(self, meshed_space): + from stride.problem.data import MeshedField + + a = MeshedField(name='a', grid=nodal_grid(meshed_space)) + a.fill(1.) + b = a.copy() + b.fill(2.) + + np.testing.assert_allclose((a + b).data, 3.) + + def test_multiply_by_scalar(self, meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField(name='a', grid=nodal_grid(meshed_space)) + field.fill(3.) + + np.testing.assert_allclose((field * 2).data, 6.) + + def test_in_place_add(self, meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField(name='a', grid=nodal_grid(meshed_space)) + field.fill(1.) + field += 1. + + np.testing.assert_allclose(field.data, 2.) + + def test_operations_preserve_node_count(self, meshed_space): + from stride.problem.data import MeshedField + + a = MeshedField(name='a', grid=nodal_grid(meshed_space)) + a.fill(1.) + + assert tuple((a * 2 + a).shape) == (27,) + + def test_elementwise_over_nodes(self, meshed_space): + from stride.problem.data import MeshedField + + a = MeshedField(name='a', grid=nodal_grid(meshed_space)) + a.allocate() + a.data[:] = np.arange(27) + + np.testing.assert_allclose((a * 2).data, np.arange(27) * 2) + + +class TestMeshedFieldGradient: + + def test_clear_grad_allocates_a_nodal_gradient(self, meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=nodal_grid(meshed_space), + needs_grad=True) + field.clear_grad() + + assert field.grad is not None + assert tuple(field.grad.shape) == (27,) + np.testing.assert_allclose(field.grad.data, 0.) + + def test_gradient_has_a_preconditioner(self, meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=nodal_grid(meshed_space), + needs_grad=True) + field.clear_grad() + + assert field.grad.prec is not None + assert tuple(field.grad.prec.shape) == (27,) + + def test_clear_grad_is_a_noop_without_needs_grad(self, meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=nodal_grid(meshed_space)) + field.clear_grad() + + assert field.grad is None diff --git a/stride/tests/test_meshed_medium.py b/stride/tests/test_meshed_medium.py new file mode 100644 index 0000000..6d25def --- /dev/null +++ b/stride/tests/test_meshed_medium.py @@ -0,0 +1,339 @@ +""" +Tests for loading a medium onto a MeshedSpace. + +This is the stride-side port of the medium path in ``ae-modelling``: + +- ``ae_modelling.tissue.nifti.NIfTILabelSampler.sample_labels`` samples a voxel + segmentation at arbitrary coordinates via the NIfTI affine + -> ``MeshedSpace.sample_labels(volume, affine)``, evaluated at mesh nodes +- ``sample_sigma`` / ``sample_eps`` map labels through a lookup table + -> ``MeshedField.from_labels(labels, lut, ...)`` +- ``ae_modelling.fem.space.DielectricSpace.admittivity`` combines the two into + ``sigma + 1j * omega * eps_r * EPS0`` + -> ordinary MeshedField arithmetic, so nothing new is needed for it +- ``ae_modelling.fem.interpolate.interpolate_medium`` also supports a + label -> value map applied per cell tag + -> ``MeshedField.from_cell_tags(mapping, ...)`` + +The sampling here is nearest-voxel, matching ``world_to_voxel``'s ``np.rint`` +plus clipping. Note that stride must not grow a nibabel dependency: the volume +and its affine are passed in as plain arrays, and reading the ``.nii`` file +stays the caller's job. +""" + +import numpy as np +import pytest + +from stride.problem.domain import Grid + + +# ae_modelling.fem.space.EPS0 +EPS0 = 8.854e-12 + + +def sigma_lut(tissue_properties): + """Conductivity indexed by label, as get_tissue_sigma_array does.""" + lut = np.zeros(max(tissue_properties) + 1) + for label, properties in tissue_properties.items(): + lut[label] = properties['sigma'] + return lut + + +def eps_lut(tissue_properties): + """Relative permittivity indexed by label, as get_tissue_eps_array does.""" + lut = np.zeros(max(tissue_properties) + 1) + for label, properties in tissue_properties.items(): + lut[label] = properties['eps_r'] + return lut + + +class TestTissueLookupTables: + """Guards on the fixture itself, so the sampling tests read unambiguously.""" + + def test_sigma_lut_is_indexed_by_label(self, tissue_properties): + lut = sigma_lut(tissue_properties) + + assert lut.shape == (6,) + assert lut[1] == pytest.approx(1.52e-1) + assert lut[3] == pytest.approx(2e0) + + def test_eps_lut_is_indexed_by_label(self, tissue_properties): + lut = eps_lut(tissue_properties) + + assert lut[1] == pytest.approx(2.19e3) + assert lut[3] == pytest.approx(1.09e2) + + +class TestSampleLabelsAtNodes: + + def test_one_label_per_node(self, meshed_space, label_volume): + volume, affine = label_volume + + labels = meshed_space.sample_labels(volume, affine=affine) + + assert labels.shape == (meshed_space.num_nodes,) + assert labels.dtype.kind == 'i' + + def test_labels_follow_the_volume(self, meshed_space, label_volume): + volume, affine = label_volume + + labels = meshed_space.sample_labels(volume, affine=affine) + + # The fixture volume is 1 for z < 2 mm and 3 above, on a 1 mm voxel + # grid; the mesh nodes sit at z = 0, 1 and 2 mm. + z = meshed_space.nodes[:, 2] + np.testing.assert_array_equal(labels[z < 1.5e-3], 1) + np.testing.assert_array_equal(labels[z > 1.5e-3], 3) + + def test_affine_translation_is_applied(self, tetra_mesh, label_volume): + from stride.problem.domain import MeshedSpace + + volume, _ = label_volume + nodes, cells = tetra_mesh + space = MeshedSpace(nodes=nodes, cells=cells) + + # Shift the volume 2 mm down in z, so every mesh node now lands in the + # volume's upper, label-3 half. + affine = np.diag([1e-3, 1e-3, 1e-3, 1.]) + affine[2, 3] = -2e-3 + + labels = space.sample_labels(volume, affine=affine) + + np.testing.assert_array_equal(labels, 3) + + def test_identity_affine_treats_nodes_as_voxel_indices(self, tetra_mesh, label_volume): + from stride.problem.domain import MeshedSpace + + volume, _ = label_volume + nodes, cells = tetra_mesh + + # Nodes at integer voxel positions, no affine given. + space = MeshedSpace(nodes=nodes * 1e3, cells=cells) + labels = space.sample_labels(volume) + + z = space.nodes[:, 2] + np.testing.assert_array_equal(labels[z < 1.5], 1) + np.testing.assert_array_equal(labels[z > 1.5], 3) + + def test_out_of_volume_nodes_are_clipped(self, tetra_mesh, label_volume): + from stride.problem.domain import MeshedSpace + + volume, affine = label_volume + nodes, cells = tetra_mesh + + # Push the mesh well past the 4x4x4 voxel volume. world_to_voxel clips + # rather than raising, so edge nodes take the nearest in-range label. + space = MeshedSpace(nodes=nodes + 1e-1, cells=cells) + labels = space.sample_labels(volume, affine=affine) + + assert labels.shape == (space.num_nodes,) + np.testing.assert_array_equal(labels, 3) + + def test_2d_volume_sampling(self, meshed_space_2d): + volume = np.zeros((4, 4), dtype=np.int64) + volume[:, 2:] = 5 + affine = np.diag([1e-3, 1e-3, 1.]) + + labels = meshed_space_2d.sample_labels(volume, affine=affine) + + y = meshed_space_2d.nodes[:, 1] + np.testing.assert_array_equal(labels[y < 1.5e-3], 0) + np.testing.assert_array_equal(labels[y > 1.5e-3], 5) + + +class TestFieldFromLabels: + + def test_conductivity_from_labels(self, meshed_space, label_volume, tissue_properties): + from stride.problem.data import MeshedField + + volume, affine = label_volume + labels = meshed_space.sample_labels(volume, affine=affine) + + sigma = MeshedField.from_labels(labels, sigma_lut(tissue_properties), + name='sigma', + grid=Grid(meshed_space, None, None)) + + assert tuple(sigma.shape) == (meshed_space.num_nodes,) + + z = meshed_space.nodes[:, 2] + np.testing.assert_allclose(sigma.data[z < 1.5e-3], 1.52e-1, rtol=1e-6) + np.testing.assert_allclose(sigma.data[z > 1.5e-3], 2e0, rtol=1e-6) + + def test_permittivity_from_labels(self, meshed_space, label_volume, tissue_properties): + from stride.problem.data import MeshedField + + volume, affine = label_volume + labels = meshed_space.sample_labels(volume, affine=affine) + + eps = MeshedField.from_labels(labels, eps_lut(tissue_properties), + name='eps', grid=Grid(meshed_space, None, None)) + + z = meshed_space.nodes[:, 2] + np.testing.assert_allclose(eps.data[z < 1.5e-3], 2.19e3, rtol=1e-6) + np.testing.assert_allclose(eps.data[z > 1.5e-3], 1.09e2, rtol=1e-6) + + def test_lut_may_be_a_mapping(self, meshed_space, label_volume): + from stride.problem.data import MeshedField + + volume, affine = label_volume + labels = meshed_space.sample_labels(volume, affine=affine) + + # Sparse or non-contiguous label sets are common in segmentations, so a + # dict has to work as well as an indexable array. + field = MeshedField.from_labels(labels, {1: 10., 3: 30.}, name='sigma', + grid=Grid(meshed_space, None, None)) + + z = meshed_space.nodes[:, 2] + np.testing.assert_allclose(field.data[z < 1.5e-3], 10.) + np.testing.assert_allclose(field.data[z > 1.5e-3], 30.) + + def test_unmapped_label_raises(self, meshed_space, label_volume): + from stride.problem.data import MeshedField + + volume, affine = label_volume + labels = meshed_space.sample_labels(volume, affine=affine) + + # Silently defaulting an unmapped tissue to zero conductivity would + # produce a plausible-looking but wrong solve. + with pytest.raises(KeyError): + MeshedField.from_labels(labels, {1: 10.}, name='sigma', + grid=Grid(meshed_space, None, None)) + + def test_label_count_must_match_nodes(self, meshed_space, tissue_properties): + from stride.problem.data import MeshedField + + with pytest.raises(ValueError): + MeshedField.from_labels(np.ones(5, dtype=np.int64), + sigma_lut(tissue_properties), + name='sigma', + grid=Grid(meshed_space, None, None)) + + +class TestFieldFromCellTags: + """The dict branch of ae_modelling.fem.interpolate.interpolate_medium.""" + + def test_values_are_assigned_per_cell(self, tagged_meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField.from_cell_tags({1: 0.1, 2: 0.5}, name='sigma', + grid=Grid(tagged_meshed_space, None, None)) + + # One value per cell, as for a DG-0 medium function. + assert tuple(field.shape) == (tagged_meshed_space.num_cells,) + + tags = tagged_meshed_space.cell_tags + np.testing.assert_allclose(field.data[tags == 1], 0.1) + np.testing.assert_allclose(field.data[tags == 2], 0.5) + + def test_missing_tag_in_mapping_raises(self, tagged_meshed_space): + from stride.problem.data import MeshedField + + # Mirrors the assertion in load_mesh that sigma's keys must cover the + # cell tag labels. + with pytest.raises(KeyError): + MeshedField.from_cell_tags({1: 0.1}, name='sigma', + grid=Grid(tagged_meshed_space, None, None)) + + def test_requires_cell_tags_on_the_space(self, meshed_space): + from stride.problem.data import MeshedField + + with pytest.raises(ValueError): + MeshedField.from_cell_tags({1: 0.1}, name='sigma', + grid=Grid(meshed_space, None, None)) + + +class TestAdmittivity: + """DielectricSpace.admittivity, expressed with MeshedField arithmetic.""" + + def _fields(self, meshed_space, label_volume, tissue_properties): + from stride.problem.data import MeshedField + + volume, affine = label_volume + labels = meshed_space.sample_labels(volume, affine=affine) + grid = Grid(meshed_space, None, None) + + sigma = MeshedField.from_labels(labels, sigma_lut(tissue_properties), + name='sigma', dtype=np.complex128, grid=grid) + eps = MeshedField.from_labels(labels, eps_lut(tissue_properties), + name='eps', dtype=np.complex128, grid=grid) + return sigma, eps + + def test_admittivity_is_sigma_plus_j_omega_eps(self, meshed_space, label_volume, + tissue_properties): + omega = 2 * np.pi * 5e5 + sigma, eps = self._fields(meshed_space, label_volume, tissue_properties) + + admittivity = sigma + eps * (1j * omega * EPS0) + + assert admittivity.data.dtype == np.complex128 + + z = meshed_space.nodes[:, 2] + expected = 1.52e-1 + 1j * omega * 2.19e3 * EPS0 + np.testing.assert_allclose(admittivity.data[z < 1.5e-3], expected, rtol=1e-6) + + def test_real_part_is_the_conductivity(self, meshed_space, label_volume, + tissue_properties): + omega = 2 * np.pi * 5e5 + sigma, eps = self._fields(meshed_space, label_volume, tissue_properties) + + admittivity = sigma + eps * (1j * omega * EPS0) + + np.testing.assert_allclose(admittivity.data.real, sigma.data.real, rtol=1e-6) + + def test_zero_frequency_reduces_to_conductivity(self, meshed_space, label_volume, + tissue_properties): + sigma, eps = self._fields(meshed_space, label_volume, tissue_properties) + + admittivity = sigma + eps * (1j * 0. * EPS0) + + np.testing.assert_allclose(admittivity.data, sigma.data, rtol=1e-6) + + +class TestMeshedMedium: + """Medium is space-agnostic; these pin that it stays that way for meshes.""" + + @pytest.fixture + def project(self, tmp_path): + return {'path': str(tmp_path), 'project_name': 'meshed_medium'} + + def _medium(self, meshed_space, label_volume, tissue_properties): + from stride.problem.data import MeshedField + from stride.problem.medium import Medium + + volume, affine = label_volume + labels = meshed_space.sample_labels(volume, affine=affine) + grid = Grid(meshed_space, None, None) + + medium = Medium(grid=grid) + medium.add(MeshedField.from_labels(labels, sigma_lut(tissue_properties), + name='sigma', grid=grid)) + medium.add(MeshedField.from_labels(labels, eps_lut(tissue_properties), + name='eps', grid=grid)) + return medium + + def test_fields_are_accessible_by_name(self, meshed_space, label_volume, + tissue_properties): + medium = self._medium(meshed_space, label_volume, tissue_properties) + + assert set(medium.fields) == {'sigma', 'eps'} + assert tuple(medium.sigma.shape) == (meshed_space.num_nodes,) + assert tuple(medium['eps'].shape) == (meshed_space.num_nodes,) + + def test_medium_round_trip(self, meshed_space, label_volume, tissue_properties, + project): + from stride.problem.data import MeshedField + from stride.problem.domain import MeshedSpace + from stride.problem.medium import Medium + + medium = self._medium(meshed_space, label_volume, tissue_properties) + medium.dump(**project) + + loaded = Medium() + loaded.add(MeshedField(name='sigma')) + loaded.add(MeshedField(name='eps')) + loaded.load(**project) + + assert isinstance(loaded.sigma.space, MeshedSpace) + assert loaded.sigma.space.num_nodes == meshed_space.num_nodes + np.testing.assert_allclose(loaded.sigma.data, medium.sigma.data) + np.testing.assert_allclose(loaded.eps.data, medium.eps.data) diff --git a/stride/tests/test_meshed_serialisation.py b/stride/tests/test_meshed_serialisation.py new file mode 100644 index 0000000..2a7fb81 --- /dev/null +++ b/stride/tests/test_meshed_serialisation.py @@ -0,0 +1,270 @@ +""" +Tests for the meshed branch of GriddedSaved (stride/problem/base.py). + +``GriddedSaved.grid_description`` has to emit a different payload per space +type, and ``GriddedSaved.load`` has to reconstruct the matching space subclass +from what it finds on disk. These tests pin both directions, plus the +structured-Space regressions, since the skeleton edited that code path. + +Contract under test: + +- a structured Space serialises as ``shape``/``spacing``/``extra``/``absorbing`` + and loads back as a Space (unchanged behaviour) +- a MeshedSpace serialises as ``nodes``/``cells`` (and ``cell_tags`` when + present) and loads back as a MeshedSpace +- the branch is chosen from the keys present under ``description.space``, not + from the top-level description +- an unrecognised space payload raises a clear, typed error +- a MeshedField round-trips its data and its mesh through HDF5 +""" + +import numpy as np +import pytest + +from stride.problem.domain import Grid, MeshedSpace, Space + + +@pytest.fixture +def project(tmp_path): + """Path/project_name pair for the HDF5 helpers.""" + return {'path': str(tmp_path), 'project_name': 'meshed'} + + +class TestGridDescriptionStructured: + """The pre-existing Space path must keep working untouched.""" + + def test_structured_space_keys(self, structured_space, project): + from stride.problem.data import ScalarField + + field = ScalarField(name='vp_field', grid=Grid(structured_space, None, None)) + description = field.grid_description() + + assert set(description['space']) == {'shape', 'spacing', 'extra', 'absorbing'} + assert tuple(description['space']['shape']) == (6, 8) + assert tuple(description['space']['extra']) == (2, 2) + + def test_structured_space_has_no_mesh_keys(self, structured_space): + from stride.problem.data import ScalarField + + field = ScalarField(name='vp_field', grid=Grid(structured_space, None, None)) + + assert 'nodes' not in field.grid_description()['space'] + + def test_structured_round_trip(self, structured_space, project): + from stride.problem.data import ScalarField + + field = ScalarField(name='vp_field', grid=Grid(structured_space, None, None)) + field.fill(1500.) + field.dump(**project) + + loaded = ScalarField(name='vp_field') + loaded.load(**project) + + assert isinstance(loaded.space, Space) + assert tuple(loaded.space.shape) == (6, 8) + np.testing.assert_allclose(loaded.space.spacing, (1e-3, 1e-3)) + np.testing.assert_allclose(loaded.space.extra, (2, 2)) + np.testing.assert_allclose(loaded.space.absorbing, (1, 1)) + np.testing.assert_allclose(loaded.data, 1500.) + + +class TestGridDescriptionMeshed: + + def test_meshed_space_keys(self, meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=Grid(meshed_space, None, None)) + space_description = field.grid_description()['space'] + + assert 'nodes' in space_description + assert 'cells' in space_description + + def test_meshed_space_omits_structured_keys(self, meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=Grid(meshed_space, None, None)) + space_description = field.grid_description()['space'] + + # 'shape' under a meshed space would make the load-time branch pick the + # structured Space reconstruction. + assert 'shape' not in space_description + assert 'spacing' not in space_description + + def test_nodes_are_the_node_table(self, meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=Grid(meshed_space, None, None)) + nodes = np.asarray(field.grid_description()['space']['nodes']) + + assert nodes.shape == (27, 3) + np.testing.assert_allclose(nodes, meshed_space.nodes) + + def test_cell_tags_are_included_when_present(self, tagged_meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=Grid(tagged_meshed_space, None, None)) + space_description = field.grid_description()['space'] + + assert 'cell_tags' in space_description + np.testing.assert_array_equal( + np.asarray(space_description['cell_tags']), + tagged_meshed_space.cell_tags, + ) + + def test_unknown_space_type_raises(self, monkeypatch): + from stride.problem.data import MeshedField + + class NotASpace: + pass + + field = MeshedField(name='sigma', shape=(4,)) + field.grid.space = NotASpace() + + # A bare `raise Exception` here would be indistinguishable from a bug + # anywhere else in the dump path. + with pytest.raises((TypeError, ValueError)): + field.grid_description() + + +class TestMeshedRoundTrip: + + def test_space_type_survives(self, meshed_space, project): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=Grid(meshed_space, None, None)) + field.fill(0.152) + field.dump(**project) + + loaded = MeshedField(name='sigma') + loaded.load(**project) + + assert isinstance(loaded.space, MeshedSpace) + assert not isinstance(loaded.space, Space) + + def test_nodes_survive(self, meshed_space, project): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=Grid(meshed_space, None, None)) + field.fill(0.152) + field.dump(**project) + + loaded = MeshedField(name='sigma') + loaded.load(**project) + + assert loaded.space.num_nodes == 27 + assert loaded.space.dim == 3 + np.testing.assert_allclose(loaded.space.nodes, meshed_space.nodes) + + def test_cells_survive(self, meshed_space, project): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=Grid(meshed_space, None, None)) + field.fill(0.152) + field.dump(**project) + + loaded = MeshedField(name='sigma') + loaded.load(**project) + + # Connectivity is what makes the node table a mesh; without it the + # reloaded space cannot be handed back to a FEM solver. + assert loaded.space.num_cells == 48 + np.testing.assert_array_equal(loaded.space.cells, meshed_space.cells) + + def test_cell_tags_survive(self, tagged_meshed_space, project): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=Grid(tagged_meshed_space, None, None)) + field.fill(1.) + field.dump(**project) + + loaded = MeshedField(name='sigma') + loaded.load(**project) + + np.testing.assert_array_equal(loaded.space.cell_tags, + tagged_meshed_space.cell_tags) + + def test_data_survives(self, meshed_space, project): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=Grid(meshed_space, None, None)) + field.allocate() + field.data[:] = np.arange(27) + field.dump(**project) + + loaded = MeshedField(name='sigma') + loaded.load(**project) + + assert tuple(loaded.shape) == (27,) + np.testing.assert_allclose(loaded.data, np.arange(27)) + + def test_bounds_survive(self, meshed_space, project): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=Grid(meshed_space, None, None)) + field.fill(1.) + field.dump(**project) + + loaded = MeshedField(name='sigma') + loaded.load(**project) + + np.testing.assert_allclose(loaded.space.origin, meshed_space.origin) + np.testing.assert_allclose(loaded.space.limit, meshed_space.limit) + + def test_vector_field_round_trip(self, meshed_space, project): + from stride.problem.data import MeshedField + + field = MeshedField(name='e_field', dim=3, grid=Grid(meshed_space, None, None)) + field.allocate() + field.data[:] = np.arange(27 * 3).reshape(27, 3) + field.dump(**project) + + loaded = MeshedField(name='e_field', dim=3) + loaded.load(**project) + + assert tuple(loaded.shape) == (27, 3) + np.testing.assert_allclose(loaded.data, np.arange(27 * 3).reshape(27, 3)) + + def test_two_dimensional_mesh_round_trip(self, meshed_space_2d, project): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=Grid(meshed_space_2d, None, None)) + field.fill(1.) + field.dump(**project) + + loaded = MeshedField(name='sigma') + loaded.load(**project) + + assert loaded.space.dim == 2 + assert loaded.space.num_nodes == 9 + np.testing.assert_allclose(loaded.space.nodes, meshed_space_2d.nodes) + + def test_existing_grid_is_not_overwritten(self, meshed_space, project): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=Grid(meshed_space, None, None)) + field.fill(1.) + field.dump(**project) + + loaded = MeshedField(name='sigma', grid=Grid(meshed_space, None, None)) + loaded.load(**project) + + # GriddedSaved.load only builds a space when the instance has none. + assert loaded.space is meshed_space + + def test_time_dependent_round_trip(self, meshed_space, project): + from stride.problem.data import MeshedField + from stride.problem.domain import Time + + time = Time(start=0., step=1e-6, num=5) + field = MeshedField(name='phi', time_dependent=True, + grid=Grid(meshed_space, time, None)) + field.fill(2.) + field.dump(**project) + + loaded = MeshedField(name='phi', time_dependent=True) + loaded.load(**project) + + assert isinstance(loaded.space, MeshedSpace) + assert loaded.time.num == 5 + assert tuple(loaded.shape) == (5, 27) + np.testing.assert_allclose(loaded.data, 2.) diff --git a/stride/tests/test_meshed_space.py b/stride/tests/test_meshed_space.py new file mode 100644 index 0000000..7e30d60 --- /dev/null +++ b/stride/tests/test_meshed_space.py @@ -0,0 +1,236 @@ +""" +Tests for MeshedSpace (stride/problem/domain.py). + +MeshedSpace is the unstructured counterpart to Space: instead of a shape and a +spacing it is defined by an explicit node list and cell connectivity, ported +from ``ae_modelling.fem.mesh.MeshDomain`` and ``ae_modelling.fem.space``. + +Contract under test: + +- ``MeshedSpace(nodes, cells=None, cell_tags=None, facet_tags=None)`` +- ``dim`` is inferred from ``nodes.shape[1]`` and must be 2 or 3 +- ``num_nodes`` / ``num_cells`` +- ``origin`` / ``limit`` / ``size`` come from the node bounding box, the + analogue of Space's origin/limit/size +- ``shape`` is ``(num_nodes,)`` -- the shape of a scalar nodal field -- with + ``extended_shape == shape``, ``extra`` and ``absorbing`` all-zero and + ``inner == (slice(0, None),)``, so that the StructuredData and GriddedSaved + machinery that reads those attributes keeps working +- ``contains_box(lower, upper)`` is the mesh-covers-the-grid check that + ``ae_modelling.fem.mesh.attach_mesh`` performs with assertions +- ``resample`` raises, because a mesh has no spacing to resample onto +- ``from_dolfinx`` adapts an in-memory DOLFINx mesh (skipped without DOLFINx) +""" + +import numpy as np +import pytest + +from .conftest import box_tetra_mesh + + +try: + import dolfinx # noqa: F401 + HAS_DOLFINX = True +except ImportError: + HAS_DOLFINX = False + + +class TestMeshedSpaceConstruction: + + def test_nodes_and_cells_are_stored(self, tetra_mesh, meshed_space): + nodes, cells = tetra_mesh + + np.testing.assert_array_equal(meshed_space.nodes, nodes) + np.testing.assert_array_equal(meshed_space.cells, cells) + + def test_dim_inferred_from_nodes(self, meshed_space, meshed_space_2d): + assert meshed_space.dim == 3 + assert meshed_space_2d.dim == 2 + + def test_counts(self, meshed_space): + assert meshed_space.num_nodes == 27 + assert meshed_space.num_cells == 48 + + def test_num_cells_is_zero_without_connectivity(self, tetra_mesh): + from stride.problem.domain import MeshedSpace + + nodes, _ = tetra_mesh + space = MeshedSpace(nodes=nodes) + + assert space.cells is None + assert space.num_cells == 0 + assert space.num_nodes == 27 + + def test_nodes_stored_as_float64(self, tetra_mesh): + from stride.problem.domain import MeshedSpace + + nodes, cells = tetra_mesh + space = MeshedSpace(nodes=nodes.astype(np.float32), cells=cells) + + # DOLFINx geometry is double precision; the node table must not be + # silently downcast, or the mesh-vs-grid bounds checks lose precision. + assert space.nodes.dtype == np.float64 + + def test_cell_tags_are_stored(self, tagged_meshed_space): + tags = tagged_meshed_space.cell_tags + + assert tags is not None + assert tags.shape == (tagged_meshed_space.num_cells,) + assert set(np.unique(tags)) == {1, 2} + + def test_cell_tags_default_to_none(self, meshed_space): + assert meshed_space.cell_tags is None + assert meshed_space.facet_tags is None + + def test_ragged_nodes_rejected(self): + from stride.problem.domain import MeshedSpace + + with pytest.raises(ValueError): + MeshedSpace(nodes=np.zeros(10)) + + def test_unsupported_dimensionality_rejected(self): + from stride.problem.domain import MeshedSpace + + # ae_modelling.fem.mesh.make_mesh only handles dim 2 and 3. + with pytest.raises(ValueError): + MeshedSpace(nodes=np.zeros((10, 4))) + + def test_out_of_range_cell_indices_rejected(self, tetra_mesh): + from stride.problem.domain import MeshedSpace + + nodes, cells = tetra_mesh + broken = cells.copy() + broken[0, 0] = len(nodes) + + with pytest.raises(ValueError): + MeshedSpace(nodes=nodes, cells=broken) + + def test_cell_tags_length_must_match_cells(self, tetra_mesh): + from stride.problem.domain import MeshedSpace + + nodes, cells = tetra_mesh + + with pytest.raises(ValueError): + MeshedSpace(nodes=nodes, cells=cells, cell_tags=np.ones(3, dtype=np.int32)) + + +class TestMeshedSpaceGeometry: + + def test_origin_and_limit_from_node_bounds(self, meshed_space): + np.testing.assert_allclose(meshed_space.origin, (0., 0., 0.)) + np.testing.assert_allclose(meshed_space.limit, (2e-3, 2e-3, 2e-3)) + + def test_size_is_the_extent(self, meshed_space): + np.testing.assert_allclose(meshed_space.size, (2e-3, 2e-3, 2e-3)) + + def test_offset_origin_is_respected(self): + from stride.problem.domain import MeshedSpace + + nodes, cells = box_tetra_mesh(shape=(3, 3, 3), spacing=(1e-3, 1e-3, 1e-3), + origin=(-5e-3, 1e-3, 0.)) + space = MeshedSpace(nodes=nodes, cells=cells) + + np.testing.assert_allclose(space.origin, (-5e-3, 1e-3, 0.)) + np.testing.assert_allclose(space.limit, (-3e-3, 3e-3, 2e-3)) + np.testing.assert_allclose(space.size, (2e-3, 2e-3, 2e-3)) + + def test_geometry_is_per_axis(self, meshed_space_2d): + assert len(meshed_space_2d.origin) == 2 + assert len(meshed_space_2d.limit) == 2 + assert len(meshed_space_2d.size) == 2 + + +class TestMeshedSpaceFieldShape: + """A scalar nodal field is flat, one value per node, with no padding.""" + + def test_shape_is_node_count(self, meshed_space): + assert tuple(meshed_space.shape) == (27,) + + def test_extended_shape_matches_shape(self, meshed_space): + assert tuple(meshed_space.extended_shape) == tuple(meshed_space.shape) + + def test_no_extra_or_absorbing_padding(self, meshed_space): + assert tuple(meshed_space.extra) == (0, 0, 0) + assert tuple(meshed_space.absorbing) == (0, 0, 0) + + def test_inner_covers_every_node(self, meshed_space): + assert meshed_space.inner == (slice(0, None),) + + values = np.arange(meshed_space.num_nodes) + np.testing.assert_array_equal(values[meshed_space.inner], values) + + +class TestMeshedSpaceBounds: + """Port of the mesh-covers-the-grid assertions in ae_modelling attach_mesh.""" + + def test_contains_its_own_bounds(self, meshed_space): + assert meshed_space.contains_box(meshed_space.origin, meshed_space.limit) + + def test_contains_an_inset_box(self, meshed_space): + assert meshed_space.contains_box((5e-4, 5e-4, 5e-4), (15e-4, 15e-4, 15e-4)) + + def test_rejects_a_box_that_overhangs(self, meshed_space): + assert not meshed_space.contains_box((0., 0., 0.), (3e-3, 2e-3, 2e-3)) + assert not meshed_space.contains_box((-1e-3, 0., 0.), (2e-3, 2e-3, 2e-3)) + + def test_tolerance_absorbs_round_off(self, meshed_space): + # A box that overhangs by less than atol counts as covered, which is + # what makes the check survive float round-off on gmsh node coordinates. + upper = tuple(each + 1e-12 for each in meshed_space.limit) + + assert meshed_space.contains_box(meshed_space.origin, upper, atol=1e-9) + assert not meshed_space.contains_box(meshed_space.origin, upper, atol=0.) + + +class TestMeshedSpaceResample: + + def test_resample_is_not_supported(self, meshed_space): + # Unlike Space, a mesh has no spacing to resample onto; remeshing is a + # separate operation and must not be silently approximated here. + with pytest.raises(NotImplementedError): + meshed_space.resample(5e-4) + + +@pytest.mark.skipif(not HAS_DOLFINX, reason='DOLFINx not available') +class TestMeshedSpaceFromDolfinx: + + def _dolfinx_box(self): + from mpi4py import MPI + + return dolfinx.mesh.create_box( + MPI.COMM_WORLD, + [np.array([0., 0., 0.]), np.array([2e-3, 2e-3, 2e-3])], + [2, 2, 2], + ) + + def test_nodes_come_from_mesh_geometry(self): + from stride.problem.domain import MeshedSpace + + mesh = self._dolfinx_box() + space = MeshedSpace.from_dolfinx(mesh) + + np.testing.assert_allclose(space.nodes, mesh.geometry.x[:, :3]) + assert space.dim == 3 + assert space.num_nodes == mesh.geometry.x.shape[0] + + def test_bounds_match_the_dolfinx_mesh(self): + from stride.problem.domain import MeshedSpace + + mesh = self._dolfinx_box() + space = MeshedSpace.from_dolfinx(mesh) + + np.testing.assert_allclose(space.origin, mesh.geometry.x.min(axis=0)) + np.testing.assert_allclose(space.limit, mesh.geometry.x.max(axis=0)) + + def test_cells_come_from_topology(self): + from stride.problem.domain import MeshedSpace + + mesh = self._dolfinx_box() + space = MeshedSpace.from_dolfinx(mesh) + + tdim = mesh.topology.dim + num_cells = mesh.topology.index_map(tdim).size_local + + assert space.num_cells == num_cells + assert space.cells.shape[1] == 4 + assert space.cells.max() < space.num_nodes From ce6ea3428ebc5f94e5b969aa74292d034ab6b6ce Mon Sep 17 00:00:00 2001 From: Andrei Danila Date: Tue, 4 Aug 2026 15:27:51 +0100 Subject: [PATCH 04/15] improved tests --- stride/tests/test_meshed_data.py | 64 ++++++++++++++++++++++- stride/tests/test_meshed_medium.py | 22 ++++++++ stride/tests/test_meshed_serialisation.py | 16 ++++++ stride/tests/test_meshed_space.py | 41 ++++++++------- 4 files changed, 123 insertions(+), 20 deletions(-) diff --git a/stride/tests/test_meshed_data.py b/stride/tests/test_meshed_data.py index 6cb5c74..4afa0c3 100644 --- a/stride/tests/test_meshed_data.py +++ b/stride/tests/test_meshed_data.py @@ -57,15 +57,26 @@ def test_shape_inferred_from_data(self, meshed_space): assert tuple(data.shape) == (27,) np.testing.assert_array_equal(data.data, values) - def test_mismatched_data_length_rejected(self, meshed_space): + @pytest.mark.parametrize('num_values', [25, 26, 29]) + def test_mismatched_data_length_rejected(self, meshed_space, num_values): from stride.problem.data import MeshedData # A nodal field with the wrong number of values cannot be interpolated # onto the mesh, so this has to fail loudly rather than at solve time. + # 25 is the case that matters: the inherited StructuredData.pad_data floors its + # pad widths, so an off-by-two would otherwise be silently edge-padded to 27. with pytest.raises(ValueError): - MeshedData(name='raw', data=np.zeros(26, dtype=np.float32), + MeshedData(name='raw', data=np.zeros(num_values, dtype=np.float32), grid=nodal_grid(meshed_space)) + def test_compression_is_rejected(self, meshed_space): + from stride.problem.data import MeshedData + + # maybe_compress crashes outright for buffers of 2.5k-10k float32 elements, and + # above that window the decompressed buffer is read-only. + with pytest.raises(ValueError, match='[Cc]ompression'): + MeshedData(name='raw', compressed=True, grid=nodal_grid(meshed_space)) + def test_default_dtype_is_float32(self, meshed_space): from stride.problem.data import MeshedData @@ -171,6 +182,55 @@ def test_explicit_shape_overrides_the_grid(self, meshed_space): assert tuple(field.shape) == (5,) +class TestMeshedFieldLocation: + """ + A meshed field is either nodal or per-cell, and `shape` alone does not say which, + so `location` is what keeps the two apart. + """ + + def test_defaults_to_nodal(self, meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', grid=nodal_grid(meshed_space)) + + assert field.location == 'nodal' + assert field.num_entities == meshed_space.num_nodes + + def test_cell_location_sizes_from_cells(self, tagged_meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField(name='sigma', location='cell', + grid=nodal_grid(tagged_meshed_space)) + + assert field.location == 'cell' + assert field.num_entities == tagged_meshed_space.num_cells + assert tuple(field.shape) == (48,) + + def test_from_cell_tags_reports_cell_location(self, tagged_meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField.from_cell_tags({1: 0.1, 2: 0.5}, name='sigma', + grid=nodal_grid(tagged_meshed_space)) + + assert field.location == 'cell' + + def test_invalid_location_rejected(self, meshed_space): + from stride.problem.data import MeshedField + + with pytest.raises(ValueError): + MeshedField(name='sigma', location='facet', grid=nodal_grid(meshed_space)) + + def test_alike_preserves_location(self, tagged_meshed_space): + from stride.problem.data import MeshedField + + field = MeshedField.from_cell_tags({1: 0.1, 2: 0.5}, name='sigma', + grid=nodal_grid(tagged_meshed_space)) + other = field.alike(name='eps') + + assert other.location == 'cell' + assert tuple(other.shape) == (48,) + + class TestMeshedFieldCopying: def test_alike_keeps_the_meshed_grid(self, meshed_space): diff --git a/stride/tests/test_meshed_medium.py b/stride/tests/test_meshed_medium.py index 6d25def..de4a017 100644 --- a/stride/tests/test_meshed_medium.py +++ b/stride/tests/test_meshed_medium.py @@ -208,6 +208,28 @@ def test_label_count_must_match_nodes(self, meshed_space, tissue_properties): name='sigma', grid=Grid(meshed_space, None, None)) + def test_label_past_the_end_of_an_array_lut_raises(self, meshed_space, tissue_properties): + from stride.problem.data import MeshedField + + labels = np.full(meshed_space.num_nodes, 9, dtype=np.int64) + + # An array lut would raise IndexError here rather than KeyError; normalise it so both + # lut forms report an unmapped label the same way. + with pytest.raises(KeyError): + MeshedField.from_labels(labels, sigma_lut(tissue_properties), name='sigma', + grid=Grid(meshed_space, None, None)) + + def test_negative_label_raises(self, meshed_space, tissue_properties): + from stride.problem.data import MeshedField + + labels = np.full(meshed_space.num_nodes, -1, dtype=np.int64) + + # Segmentations do use -1 as a sentinel, and a negative index would otherwise + # silently read from the end of the lookup table. + with pytest.raises(KeyError): + MeshedField.from_labels(labels, sigma_lut(tissue_properties), name='sigma', + grid=Grid(meshed_space, None, None)) + class TestFieldFromCellTags: """The dict branch of ae_modelling.fem.interpolate.interpolate_medium.""" diff --git a/stride/tests/test_meshed_serialisation.py b/stride/tests/test_meshed_serialisation.py index 2a7fb81..e74b70d 100644 --- a/stride/tests/test_meshed_serialisation.py +++ b/stride/tests/test_meshed_serialisation.py @@ -251,6 +251,22 @@ def test_existing_grid_is_not_overwritten(self, meshed_space, project): # GriddedSaved.load only builds a space when the instance has none. assert loaded.space is meshed_space + def test_cell_location_round_trip(self, tagged_meshed_space, project): + from stride.problem.data import MeshedField + + field = MeshedField.from_cell_tags({1: 0.1, 2: 0.5}, name='sigma', + grid=Grid(tagged_meshed_space, None, None)) + field.dump(**project) + + loaded = MeshedField(name='sigma') + loaded.load(**project) + + # Without `location` on the description the reloaded field would default to nodal + # and disagree with its own shape. + assert loaded.location == 'cell' + assert tuple(loaded.shape) == (48,) + np.testing.assert_allclose(loaded.data, field.data) + def test_time_dependent_round_trip(self, meshed_space, project): from stride.problem.data import MeshedField from stride.problem.domain import Time diff --git a/stride/tests/test_meshed_space.py b/stride/tests/test_meshed_space.py index 7e30d60..e59b476 100644 --- a/stride/tests/test_meshed_space.py +++ b/stride/tests/test_meshed_space.py @@ -140,24 +140,29 @@ def test_geometry_is_per_axis(self, meshed_space_2d): assert len(meshed_space_2d.size) == 2 -class TestMeshedSpaceFieldShape: - """A scalar nodal field is flat, one value per node, with no padding.""" - - def test_shape_is_node_count(self, meshed_space): - assert tuple(meshed_space.shape) == (27,) - - def test_extended_shape_matches_shape(self, meshed_space): - assert tuple(meshed_space.extended_shape) == tuple(meshed_space.shape) - - def test_no_extra_or_absorbing_padding(self, meshed_space): - assert tuple(meshed_space.extra) == (0, 0, 0) - assert tuple(meshed_space.absorbing) == (0, 0, 0) - - def test_inner_covers_every_node(self, meshed_space): - assert meshed_space.inner == (slice(0, None),) - - values = np.arange(meshed_space.num_nodes) - np.testing.assert_array_equal(values[meshed_space.inner], values) +class TestMeshedSpaceIsNotAGrid: + """ + A MeshedSpace must not masquerade as a structured grid. + + An earlier version of this class asserted the opposite: that MeshedSpace exposed + ``shape``/``extended_shape``/``extra``/``absorbing``/``inner`` as a compatibility surface. That + was wrong. Nothing meshed reads them — MeshedData derives its shape from ``num_nodes``, the way + SparseField uses ``num`` — and their only effect was to let a structured field accept a mesh and + then produce plausible nonsense when plotted or resampled. + """ + + @pytest.mark.parametrize('attribute', ['shape', 'extended_shape', 'extra', + 'absorbing', 'inner', 'spacing', 'grid']) + def test_has_no_grid_attributes(self, meshed_space, attribute): + assert not hasattr(meshed_space, attribute) + + def test_structured_field_rejects_a_mesh(self, meshed_space): + from stride.problem.data import ScalarField + from stride.problem.domain import Grid + + # Fails on the first grid attribute it reaches for, rather than silently constructing. + with pytest.raises(AttributeError): + ScalarField(name='sigma', grid=Grid(meshed_space, None, None)) class TestMeshedSpaceBounds: From 4e0d7536c9d77e64bf52f7cb8b3801b1f32aaaab Mon Sep 17 00:00:00 2001 From: Andrei Danila Date: Tue, 4 Aug 2026 15:32:51 +0100 Subject: [PATCH 05/15] discussed changes to base, data and domain, pending further verification --- stride/problem/base.py | 76 +++++-- stride/problem/data.py | 442 ++++++++++++++++++++++++++++++++++++++- stride/problem/domain.py | 258 ++++++++++++++++++++++- 3 files changed, 758 insertions(+), 18 deletions(-) diff --git a/stride/problem/base.py b/stride/problem/base.py index 0fa7239..521b0f9 100644 --- a/stride/problem/base.py +++ b/stride/problem/base.py @@ -1,7 +1,7 @@ from mosaic import h5 -from .domain import Space, Time, SlowTime, Grid +from .domain import Space, MeshedSpace, Time, SlowTime, Grid __all__ = ['Gridded', 'Saved', 'GriddedSaved', 'ProblemBase'] @@ -246,15 +246,27 @@ def load(self, *args, **kwargs): description = file.load(filter=kwargs.pop('filter', None), only=kwargs.pop('only', None)) if 'space' in description and self._grid.space is None: - if 'shape' in description: - space = Space(shape=description.space.shape, - spacing=description.space.spacing, - extra=description.space.extra, - absorbing=description.space.absorbing) - elif 'nodes' in description: - space = MeshedSpace(nodes=description.space.nodes) + # NOTE the discriminator has to be the keys of description.space, not those of + # description: StructuredData.__get_desc__ also writes a top-level 'shape', so + # testing the outer description would send every field down the structured branch + space_description = description.space + + if 'shape' in space_description: + space = Space(shape=space_description.shape, + spacing=space_description.spacing, + extra=space_description.extra, + absorbing=space_description.absorbing) + + elif 'nodes' in space_description: + space = MeshedSpace( + nodes=self._materialise(space_description.nodes), + cells=self._materialise(space_description.get('cells', None)), + cell_tags=self._materialise(space_description.get('cell_tags', None)), + ) + else: - raise Exception + raise ValueError('Unrecognised space description with keys %s' + % sorted(space_description.keys())) self._grid.space = space @@ -279,6 +291,31 @@ def load(self, *args, **kwargs): kwargs['filename'] = kwargs.pop('filename', file.filename) self.__set_desc__(description, **kwargs) + @staticmethod + def _materialise(value): + """ + Read a lazily-loaded description entry into memory. + + Loading a description defaults to being lazy, which yields the open dataset with a + ``load`` attribute attached rather than an ndarray. Anything that has to outlive the + open file must be materialised explicitly. + + Parameters + ---------- + value : object + Entry of a loaded description. + + Returns + ------- + object + The entry, read into memory if it was lazy. + + """ + if hasattr(value, 'load'): + return value.load() + + return value + def grid_description(self): """ Get a description of the grid of the object. @@ -300,13 +337,24 @@ def grid_description(self): 'extra': space.extra, 'absorbing': space.absorbing, } + elif isinstance(space, MeshedSpace): - #stand-in attribute - grid_description['space'] = { - 'nodes': space.nodes - } + # the connectivity travels with the nodes: a node table on its own is not a mesh + # and cannot be handed back to a solver. Facet tags are deliberately left out, + # being meaningless without the facet connectivity, which is not stored either + space_description = {'nodes': space.nodes} + + if space.cells is not None: + space_description['cells'] = space.cells + + if space.cell_tags is not None: + space_description['cell_tags'] = space.cell_tags + + grid_description['space'] = space_description + else: - raise Exception + raise TypeError('Cannot serialise a grid with space of type %s' + % type(space).__name__) if self.time is not None: time = self.time diff --git a/stride/problem/data.py b/stride/problem/data.py index 3a158cb..6bc430a 100644 --- a/stride/problem/data.py +++ b/stride/problem/data.py @@ -21,6 +21,7 @@ from mosaic.file_manipulation import h5 from .base import GriddedSaved +from .domain import MeshedSpace from ..core import Variable from .. import plotting @@ -781,12 +782,447 @@ def __set_desc__(self, description, **kwargs): @mosaic.tessera -class MeshedData(StructuredData, GriddedSaved): - pass +class MeshedData(StructuredData): + """ + Objects of this type represent data defined over an unstructured mesh. + + This is the mesh counterpart of StructuredData: the buffer is flat, one value per mesh entity, + and there is no inner/extended domain because a mesh has no padding. + + The shape is derived from the :class:`~stride.problem.domain.MeshedSpace` of the grid rather + than from a ``space.shape``, in the same way that :class:`SparseField` derives it from ``num``. + + Parameters + ---------- + name : str + Name of the data. + location : str, optional + Mesh entity the data lives on, ``nodal`` (one value per node) or ``cell`` (one value per + cell, the discontinuous piecewise-constant case), defaults to ``nodal``. + shape : tuple, optional + Shape of the data, derived from the grid if not given. + dtype : data-type, optional + Data type of the data, defaults to float32. + data : ndarray, optional + Data with which to initialise the internal buffer. + grid : Grid or any of MeshedSpace or Time + Grid on which the Problem is defined. + + """ + + def __init__(self, **kwargs): + if kwargs.get('compressed', False): + # maybe_compress crashes for buffers of 2.5k-10k float32 elements (byte_sample uses + # len() where it means nbytes), and above that window the decompressed buffer is a + # read-only np.frombuffer view, so in-place assignment fails + raise ValueError('Compression is not supported for meshed data') + + location = kwargs.pop('location', 'nodal') + if location not in ('nodal', 'cell'): + raise ValueError('Location must be "nodal" or "cell", got %s' % location) + + self._location = location + + data = kwargs.get('data', None) + + super().__init__(**kwargs) + + if self._shape is None and isinstance(self.space, MeshedSpace): + self._init_shape() + + if data is not None and isinstance(self.space, MeshedSpace): + expected = self.num_entities + given = np.asarray(data).shape[0] + + # this has to be explicit: the inherited pad_data floors its pad widths, so an + # off-by-two would otherwise be silently edge-padded to the right length + if given != expected: + raise ValueError('Data has %d values but the mesh has %d %s entities' + % (given, expected, self._location)) + + def _init_shape(self, fill_shape=True): + shape = (self.num_entities,) + + if fill_shape: + self._shape = shape + self._extended_shape = shape + self._inner = (slice(0, None),) + + @property + def location(self): + """ + Mesh entity the data lives on, ``nodal`` or ``cell``. + + """ + return self._location + + @property + def num_entities(self): + """ + Number of mesh entities the data is defined over. + + """ + return self.num_cells if self._location == 'cell' else self.num_nodes + + @property + def num_nodes(self): + """ + Number of nodes in the mesh. + + """ + return self.space.num_nodes + + @property + def num_cells(self): + """ + Number of cells in the mesh. + + """ + return self.space.num_cells + + def alike(self, *args, **kwargs): + """ + Create a data object that shares its characteristics with this object. + + Returns + ------- + MeshedData + Newly created MeshedData. + + """ + kwargs['location'] = kwargs.pop('location', self._location) + + return super().alike(*args, **kwargs) + + def detach(self, *args, **kwargs): + """ + Create a copy of the variable that is detached from the original graph. + + Returns + ------- + MeshedData + Detached variable. + + """ + kwargs['location'] = kwargs.pop('location', self._location) + + return super().detach(*args, **kwargs) + + def as_parameter(self, *args, **kwargs): + """ + Create a copy of the variable, detached and re-initialised as a parameter. + + Returns + ------- + MeshedData + Detached variable. + + """ + kwargs['location'] = kwargs.pop('location', self._location) + + return super().as_parameter(*args, **kwargs) + + def pad_data(self, data, smooth=False): + """ + Padding is a no-op on a mesh, which has no extended domain. + + Parameters + ---------- + data : ndarray + Array to pad. + smooth : bool, optional + Unused. + + Returns + ------- + ndarray + The input, unchanged. + + """ + return data + + def plot(self, **kwargs): + """ + Plotting meshed data is not implemented, so that Medium.plot does not fail on a + medium made of meshed fields. + + Returns + ------- + + """ + pass + + def __get_desc__(self, **kwargs): + description = super().__get_desc__(**kwargs) + description['location'] = self._location + + return description + + def __set_desc__(self, description, **kwargs): + super().__set_desc__(description, **kwargs) + + location = description.get('location', 'nodal') + self._location = location.decode() if isinstance(location, bytes) else location + @mosaic.tessera class MeshedField(MeshedData): - pass + """ + Objects of this type describe a field defined over an unstructured mesh. Meshed fields + can also be time-dependent, and can carry a vector value per entity. + + Parameters + ---------- + name : str + Name of the data. + dim : int, optional + Number of components at every mesh entity, defaults to 1. + time_dependent : bool, optional + Whether or not the field is time-dependent, defaults to False. + slow_time_dependent : bool, optional + Whether or not the field is slow-time dependent, defaults to False. + location : str, optional + Mesh entity the field lives on, ``nodal`` or ``cell``, defaults to ``nodal``. + dtype : data-type, optional + Data type of the data, defaults to float32. + grid : Grid or any of MeshedSpace or Time + Grid on which the Problem is defined. + + """ + + def __init__(self, **kwargs): + # these have to be in place before MeshedData.__init__ calls _init_shape + self._dim = kwargs.pop('dim', 1) + self._time_dependent = kwargs.pop('time_dependent', False) + self._slow_time_dependent = kwargs.pop('slow_time_dependent', False) + + super().__init__(**kwargs) + + def _init_shape(self, fill_shape=True): + shape = () + inner = () + + if self._time_dependent: + shape += (self.time.num,) + inner += (self.time.inner,) + + if self._slow_time_dependent: + shape += (self.slow_time.num,) + inner += (self.slow_time.inner,) + + if self._dim > 1: + shape += (self.num_entities, self._dim) + inner += (slice(0, None), slice(0, None)) + else: + shape += (self.num_entities,) + inner += (slice(0, None),) + + if fill_shape: + self._shape = shape + self._extended_shape = shape + self._inner = inner + + @property + def dim(self): + """ + Number of components at every mesh entity. + + """ + return self._dim + + @property + def time_dependent(self): + """ + Whether or not the field is time dependent. + + """ + return self._time_dependent + + @property + def slow_time_dependent(self): + """ + Whether or not the field is slow-time dependent. + + """ + return self._slow_time_dependent + + def alike(self, *args, **kwargs): + """ + Create a data object that shares its characteristics with this object. + + Returns + ------- + MeshedField + Newly created MeshedField. + + """ + kwargs['dim'] = kwargs.pop('dim', self._dim) + kwargs['time_dependent'] = kwargs.pop('time_dependent', self._time_dependent) + kwargs['slow_time_dependent'] = kwargs.pop('slow_time_dependent', + self._slow_time_dependent) + + return super().alike(*args, **kwargs) + + def detach(self, *args, **kwargs): + """ + Create a copy of the variable that is detached from the original graph. + + Returns + ------- + MeshedField + Detached variable. + + """ + kwargs['dim'] = kwargs.pop('dim', self._dim) + kwargs['time_dependent'] = kwargs.pop('time_dependent', self._time_dependent) + kwargs['slow_time_dependent'] = kwargs.pop('slow_time_dependent', + self._slow_time_dependent) + + return super().detach(*args, **kwargs) + + def as_parameter(self, *args, **kwargs): + """ + Create a copy of the variable, detached and re-initialised as a parameter. + + Returns + ------- + MeshedField + Detached variable. + + """ + kwargs['dim'] = kwargs.pop('dim', self._dim) + kwargs['time_dependent'] = kwargs.pop('time_dependent', self._time_dependent) + kwargs['slow_time_dependent'] = kwargs.pop('slow_time_dependent', + self._slow_time_dependent) + + return super().as_parameter(*args, **kwargs) + + @staticmethod + def values_from_labels(labels, lut): + """ + Map an integer label array through a lookup table. + + This is how a segmentation becomes a material property: ``labels`` comes from sampling a + label volume, and ``lut`` maps each label to a conductivity or permittivity. + + Parameters + ---------- + labels : ndarray + Integer labels, one per mesh entity. + lut : ndarray or dict + Value per label, either indexable by label or a mapping. + + Returns + ------- + ndarray + Value for every entry of ``labels``. + + """ + labels = np.asarray(labels) + unique = np.unique(labels) + + if isinstance(lut, dict): + missing = [int(each) for each in unique if int(each) not in lut] + + if len(missing): + raise KeyError('No value provided for labels %s' % missing) + + return np.asarray([lut[int(each)] for each in labels]) + + lut = np.asarray(lut) + + # an array lut would otherwise raise IndexError past the end and, worse, silently index + # from the end for a negative label, which segmentations do use as a sentinel + out_of_range = [int(each) for each in unique if each < 0 or each >= lut.shape[0]] + + if len(out_of_range): + raise KeyError('Labels %s fall outside a lookup table of size %d' + % (out_of_range, lut.shape[0])) + + return lut[labels] + + @classmethod + def from_labels(cls, labels, lut, **kwargs): + """ + Create a nodal field by mapping sampled labels through a lookup table. + + Note that, being a classmethod, this always builds a local instance. To create a parameter + or a remote instance, use :meth:`values_from_labels` and pass the result as ``data``. + + Parameters + ---------- + labels : ndarray + Integer label per node, as returned by ``MeshedSpace.sample_labels``. + lut : ndarray or dict + Value per label. + + Returns + ------- + MeshedField + Newly created MeshedField. + + """ + labels = np.asarray(labels) + field = cls(**kwargs) + + if labels.shape != (field.num_entities,): + raise ValueError('Expected %d labels, one per %s entity, got shape %s' + % (field.num_entities, field.location, (labels.shape,))) + + field.allocate() + field.data[:] = cls.values_from_labels(labels, lut) + + return field + + @classmethod + def from_cell_tags(cls, mapping, **kwargs): + """ + Create a per-cell field by mapping the cell tags of the mesh through a value mapping. + + This is the discontinuous piecewise-constant medium: one value per cell, taken from the + material tag the mesh generator assigned to it. Unlike the fields built by + :meth:`from_labels`, the resulting shape is ``(num_cells,)``. + + Parameters + ---------- + mapping : dict or ndarray + Value per cell tag. Must cover every tag present in the mesh. + + Returns + ------- + MeshedField + Newly created MeshedField. + + """ + grid = kwargs.get('grid', None) + space = kwargs.get('space', None) if grid is None else grid.space + + if space is None or space.cell_tags is None: + raise ValueError('Creating a field from cell tags needs a space carrying cell tags') + + kwargs['location'] = 'cell' + field = cls(**kwargs) + + field.allocate() + field.data[:] = cls.values_from_labels(space.cell_tags, mapping) + + return field + + def __get_desc__(self, **kwargs): + description = super().__get_desc__(**kwargs) + description['dim'] = self._dim + description['time_dependent'] = self._time_dependent + description['slow_time_dependent'] = self._slow_time_dependent + + return description + + def __set_desc__(self, description, **kwargs): + super().__set_desc__(description, **kwargs) + + self._dim = description.get('dim', 1) + self._time_dependent = description.get('time_dependent', False) + self._slow_time_dependent = description.get('slow_time_dependent', False) + @mosaic.tessera class Scalar(StructuredData): diff --git a/stride/problem/domain.py b/stride/problem/domain.py index 7d11d17..fc23632 100644 --- a/stride/problem/domain.py +++ b/stride/problem/domain.py @@ -1,4 +1,5 @@ +import warnings import numpy as np from cached_property import cached_property @@ -259,8 +260,263 @@ def extended_grid(self): for dim in range(self.dim)] return tuple(axes) + class MeshedSpace: - pass + """ + This defines an unstructured spatial mesh over which the problem is defined. + + Where a :class:`Space` is fully determined by a ``shape`` and a ``spacing``, a MeshedSpace is + determined by an explicit table of ``nodes`` and, optionally, the ``cells`` that connect them. + It is the counterpart used by finite-element physics, and is a sibling of Space rather than a + subclass of it. + + A MeshedSpace deliberately has no ``shape``, ``extended_shape``, ``extra``, ``absorbing``, + ``spacing``, ``inner`` or ``grid``. Those describe a regular grid and a mesh has no analogue of + any of them, so code that requires a structured grid fails immediately when handed a mesh + instead of silently producing a plausible but meaningless result. + + Note also that, unlike Space, ``size`` is not an alias for ``limit``: the mesh origin need not + be at zero, so the extent and the upper bound differ. + + Parameters + ---------- + nodes : ndarray + Node coordinates, of shape ``(num_nodes, dim)``, in metres. + cells : ndarray, optional + Node indices making up each cell, of shape ``(num_cells, nodes_per_cell)``. The number of + nodes per cell is not fixed: linear tetrahedra have 4, quadratic tetrahedra have 10. + cell_tags : ndarray, optional + Material tag for every cell, of shape ``(num_cells,)``. Tags are opaque integers as far as + the MeshedSpace is concerned, and are only meaningful against the label table of whoever + generated the mesh. + facet_tags : optional + Boundary tags, stored as given and not interpreted. These are not serialised, because a + facet tag is meaningless without the facet connectivity, which is not stored either. + + """ + + def __init__(self, nodes=None, cells=None, cell_tags=None, facet_tags=None): + nodes = np.asarray(nodes, dtype=np.float64) + + if nodes.ndim != 2: + raise ValueError('Nodes must be a (num_nodes, dim) array, got shape %s' + % (nodes.shape,)) + + dim = nodes.shape[1] + if dim not in (2, 3): + raise ValueError('Only 2 or 3 dimensions are supported, got %d' % dim) + + if cells is not None: + cells = np.asarray(cells, dtype=np.int32) + + if cells.ndim != 2: + raise ValueError('Cells must be a (num_cells, nodes_per_cell) array, got shape %s' + % (cells.shape,)) + + if cells.size and (cells.min() < 0 or cells.max() >= nodes.shape[0]): + raise ValueError('Cells reference node indices outside [0, %d)' % nodes.shape[0]) + + if cell_tags is not None: + cell_tags = np.asarray(cell_tags) + + if cells is None: + raise ValueError('Cell tags were given without any cells') + + if cell_tags.shape != (cells.shape[0],): + raise ValueError('Cell tags must have one entry per cell, expected %d ' + 'but got shape %s' % (cells.shape[0], (cell_tags.shape,))) + + self.dim = dim + self.nodes = nodes + self.cells = cells + self.cell_tags = cell_tags + self.facet_tags = facet_tags + + # eagerly, as plain tuples: a cached_property would end up in __dict__ and be pickled + # along with the space every time a field travels to a worker + self.origin = tuple(nodes.min(axis=0)) + self.limit = tuple(nodes.max(axis=0)) + + @property + def num_nodes(self): + """ + Number of nodes in the mesh. + + """ + return int(self.nodes.shape[0]) + + @property + def num_cells(self): + """ + Number of cells in the mesh, zero if no connectivity is defined. + + """ + return 0 if self.cells is None else int(self.cells.shape[0]) + + @property + def size(self): + """ + Axis-wise extent of the mesh, as a tuple. + + """ + return tuple(each_limit - each_origin + for each_limit, each_origin in zip(self.limit, self.origin)) + + def contains_box(self, lower, upper, atol=1e-9): + """ + Whether the mesh covers an axis-aligned box. + + This is the check needed before evaluating a finite-element solution onto a structured + grid: every point of the grid has to be owned by some cell of the mesh. + + Parameters + ---------- + lower : tuple or ndarray + Lower corner of the box, in metres. + upper : tuple or ndarray + Upper corner of the box, in metres. + atol : float, optional + Absolute tolerance, to absorb round-off in the node coordinates, defaults to 1e-9. + + Returns + ------- + bool + Whether the mesh covers the box. + + """ + lower = np.asarray(lower, dtype=np.float64) + upper = np.asarray(upper, dtype=np.float64) + + return bool((self.nodes.min(axis=0) <= lower + atol).all() + and (self.nodes.max(axis=0) >= upper - atol).all()) + + def sample_labels(self, volume, affine=None): + """ + Sample a voxelised label volume at the mesh nodes, using nearest-neighbour lookup. + + This is how a segmentation is transferred onto a mesh in order to build a medium. Nodes + that fall outside the volume are clipped to the nearest in-range voxel rather than raising. + + Parameters + ---------- + volume : ndarray + Label volume, with as many dimensions as the mesh. + affine : ndarray, optional + Homogeneous ``(dim+1, dim+1)`` matrix mapping voxel indices to node coordinates. If + not given, the node coordinates are taken to be voxel indices already. + + Returns + ------- + ndarray + Label at every node, of shape ``(num_nodes,)``. + + """ + volume = np.asarray(volume) + + if volume.ndim != self.dim: + raise ValueError('Volume has %d dimensions but the mesh has %d' + % (volume.ndim, self.dim)) + + if affine is None: + indices = self.nodes + else: + affine = np.asarray(affine, dtype=np.float64) + + if affine.shape != (self.dim + 1, self.dim + 1): + raise ValueError('Affine must have shape (%d, %d), got %s' + % (self.dim + 1, self.dim + 1, (affine.shape,))) + + homogeneous = np.hstack([self.nodes, np.ones((self.num_nodes, 1))]) + indices = (np.linalg.inv(affine) @ homogeneous.T).T[:, :self.dim] + + indices = np.rint(indices).astype(np.int64) + for axis in range(self.dim): + indices[:, axis] = np.clip(indices[:, axis], 0, volume.shape[axis] - 1) + + # rint before astype: label volumes read from NIfTI are floats, and truncating turns a + # stored 2.9999 into 2 + return np.rint(volume[tuple(indices.T)]).astype(np.int64) + + def resample(self, *args, **kwargs): + """ + Not available for a MeshedSpace. + + A mesh has no spacing to resample onto, and generating a new mesh is a separate operation + that must not be silently approximated here. + + Returns + ------- + + """ + raise NotImplementedError('A MeshedSpace cannot be resampled, it has no spacing. ' + 'Generate a new mesh instead.') + + @classmethod + def from_dolfinx(cls, mesh, cell_tags=None, facet_tags=None): + """ + Create a MeshedSpace from an in-memory DOLFINx mesh. + + This is only correct for a serial mesh. Under MPI, DOLFINx node coordinates and the cell + dofmap are rank-local and include ghost entities, so the resulting MeshedSpace describes + one partition rather than the whole mesh: node and cell counts under-report, ``origin`` and + ``limit`` bound a sub-box, and nodes shared between ranks appear more than once. A warning + is issued in that case rather than an error, so that experimentation is still possible. + + Parameters + ---------- + mesh : dolfinx.mesh.Mesh + Mesh to adapt. + cell_tags : dolfinx.mesh.MeshTags, optional + Cell tags, which are sparse and get densified to one entry per cell. Cells that carry + no tag are filled with -1. + facet_tags : optional + Facet tags, stored as given. + + Returns + ------- + MeshedSpace + Newly created MeshedSpace. + + """ + if mesh.comm.size > 1: + warnings.warn('MeshedSpace.from_dolfinx is building from rank-local arrays, so the ' + 'resulting space describes this rank\'s partition and not the whole ' + 'mesh. Build the mesh on MPI.COMM_SELF, or gather it, to avoid this.') + + # geometry.x is always padded to three columns, whatever the geometric dimension + dim = mesh.geometry.dim + nodes = np.asarray(mesh.geometry.x)[:, :dim] + + topology_dim = mesh.topology.dim + num_cells = mesh.topology.index_map(topology_dim).size_local + + # the geometry dofmap, not the topology connectivity, is what indexes geometry.x + dofmap = mesh.geometry.dofmap + + if hasattr(dofmap, 'offsets'): + # DOLFINx <= 0.7 exposes a flat array plus offsets. Cells of a single mesh are all + # of the same type, so the first offset step gives the nodes per cell + offsets = np.asarray(dofmap.offsets) + nodes_per_cell = int(offsets[1] - offsets[0]) + cells = np.asarray(dofmap.array).reshape(-1, nodes_per_cell) + + else: + cells = np.asarray(dofmap).reshape(num_cells, -1) + + cells = cells[:num_cells] + + tags = None + if cell_tags is not None: + # MeshTags are a sparse (indices, values) pair, so densify to one entry per cell + tags = np.full(num_cells, -1, dtype=np.int32) + indices = np.asarray(cell_tags.indices) + values = np.asarray(cell_tags.values) + + owned = indices < num_cells + tags[indices[owned]] = values[owned] + + return cls(nodes=nodes, cells=cells, cell_tags=tags, facet_tags=facet_tags) + class Time: """ From b0f8bcaac733918c4677f56d8e006aed1c767635 Mon Sep 17 00:00:00 2001 From: Andrei Danila Date: Tue, 4 Aug 2026 16:11:53 +0100 Subject: [PATCH 06/15] better guards in MeshedData to make sure correct kind of space is used --- stride/problem/data.py | 12 ++++++++++-- stride/tests/test_meshed_space.py | 18 ++++++++++++++++++ 2 files changed, 28 insertions(+), 2 deletions(-) diff --git a/stride/problem/data.py b/stride/problem/data.py index 6bc430a..d9b04eb 100644 --- a/stride/problem/data.py +++ b/stride/problem/data.py @@ -827,10 +827,18 @@ def __init__(self, **kwargs): super().__init__(**kwargs) - if self._shape is None and isinstance(self.space, MeshedSpace): + # a space of the wrong kind has to be caught here. Everything downstream is guarded by + # isinstance checks that would silently skip instead, leaving a field with no shape that + # allocates and fills without complaint. A space of None is valid, and is what an + # instance about to be loaded from file looks like + if self.space is not None and not isinstance(self.space, MeshedSpace): + raise ValueError('Meshed data needs a MeshedSpace, got %s' + % type(self.space).__name__) + + if self._shape is None and self.space is not None: self._init_shape() - if data is not None and isinstance(self.space, MeshedSpace): + if data is not None and self.space is not None: expected = self.num_entities given = np.asarray(data).shape[0] diff --git a/stride/tests/test_meshed_space.py b/stride/tests/test_meshed_space.py index e59b476..9e30d91 100644 --- a/stride/tests/test_meshed_space.py +++ b/stride/tests/test_meshed_space.py @@ -164,6 +164,24 @@ def test_structured_field_rejects_a_mesh(self, meshed_space): with pytest.raises(AttributeError): ScalarField(name='sigma', grid=Grid(meshed_space, None, None)) + def test_meshed_field_rejects_a_structured_space(self, structured_space): + from stride.problem.data import MeshedField + from stride.problem.domain import Grid + + # The reverse direction needs an explicit check: MeshedData sizes itself behind an + # isinstance test, which would otherwise skip and leave a field with no shape that + # still allocates and fills without complaint. + with pytest.raises(ValueError, match='MeshedSpace'): + MeshedField(name='sigma', grid=Grid(structured_space, None, None)) + + def test_no_space_is_still_allowed(self): + from stride.problem.data import MeshedField + + # This is what an instance about to be loaded from file looks like. + field = MeshedField(name='sigma') + + assert field.space is None + class TestMeshedSpaceBounds: """Port of the mesh-covers-the-grid assertions in ae_modelling attach_mesh.""" From e4a7a39b4d7977c33f50dfedf45222c1510847eb Mon Sep 17 00:00:00 2001 From: Andrei Danila Date: Fri, 14 Aug 2026 14:41:49 +0100 Subject: [PATCH 07/15] carlos comments round 1 --- stride/problem/data.py | 22 +- stride/problem/domain.py | 13 +- stride/tests/conftest.py | 37 ++- stride/tests/test_meshed_data.py | 182 +++++--------- stride/tests/test_meshed_medium.py | 275 +++++++++------------- stride/tests/test_meshed_serialisation.py | 107 +++------ stride/tests/test_meshed_space.py | 50 +--- 7 files changed, 247 insertions(+), 439 deletions(-) diff --git a/stride/problem/data.py b/stride/problem/data.py index d9b04eb..8d882fa 100644 --- a/stride/problem/data.py +++ b/stride/problem/data.py @@ -797,8 +797,8 @@ class MeshedData(StructuredData): name : str Name of the data. location : str, optional - Mesh entity the data lives on, ``nodal`` (one value per node) or ``cell`` (one value per - cell, the discontinuous piecewise-constant case), defaults to ``nodal``. + Mesh entity the data lives on, ``node`` (one value per node) or ``cell`` (one value per + cell, the discontinuous piecewise-constant case), defaults to ``node``. shape : tuple, optional Shape of the data, derived from the grid if not given. dtype : data-type, optional @@ -817,9 +817,9 @@ def __init__(self, **kwargs): # read-only np.frombuffer view, so in-place assignment fails raise ValueError('Compression is not supported for meshed data') - location = kwargs.pop('location', 'nodal') - if location not in ('nodal', 'cell'): - raise ValueError('Location must be "nodal" or "cell", got %s' % location) + location = kwargs.pop('location', 'node') + if location not in ('node', 'cell'): + raise ValueError('Location must be "node" or "cell", got %s' % location) self._location = location @@ -859,7 +859,7 @@ def _init_shape(self, fill_shape=True): @property def location(self): """ - Mesh entity the data lives on, ``nodal`` or ``cell``. + Mesh entity the data lives on, ``node`` or ``cell``. """ return self._location @@ -969,7 +969,7 @@ def __get_desc__(self, **kwargs): def __set_desc__(self, description, **kwargs): super().__set_desc__(description, **kwargs) - location = description.get('location', 'nodal') + location = description.get('location', 'node') self._location = location.decode() if isinstance(location, bytes) else location @@ -990,7 +990,7 @@ class MeshedField(MeshedData): slow_time_dependent : bool, optional Whether or not the field is slow-time dependent, defaults to False. location : str, optional - Mesh entity the field lives on, ``nodal`` or ``cell``, defaults to ``nodal``. + Mesh entity the field lives on, ``node`` or ``cell``, defaults to ``node``. dtype : data-type, optional Data type of the data, defaults to float32. grid : Grid or any of MeshedSpace or Time @@ -1111,7 +1111,7 @@ def values_from_labels(labels, lut): Map an integer label array through a lookup table. This is how a segmentation becomes a material property: ``labels`` comes from sampling a - label volume, and ``lut`` maps each label to a conductivity or permittivity. + label volume, and ``lut`` maps each label to a value. Parameters ---------- @@ -1152,7 +1152,7 @@ def values_from_labels(labels, lut): @classmethod def from_labels(cls, labels, lut, **kwargs): """ - Create a nodal field by mapping sampled labels through a lookup table. + Create a node field by mapping sampled labels through a lookup table. Note that, being a classmethod, this always builds a local instance. To create a parameter or a remote instance, use :meth:`values_from_labels` and pass the result as ``data``. @@ -1177,7 +1177,6 @@ def from_labels(cls, labels, lut, **kwargs): raise ValueError('Expected %d labels, one per %s entity, got shape %s' % (field.num_entities, field.location, (labels.shape,))) - field.allocate() field.data[:] = cls.values_from_labels(labels, lut) return field @@ -1211,7 +1210,6 @@ def from_cell_tags(cls, mapping, **kwargs): kwargs['location'] = 'cell' field = cls(**kwargs) - field.allocate() field.data[:] = cls.values_from_labels(space.cell_tags, mapping) return field diff --git a/stride/problem/domain.py b/stride/problem/domain.py index fc23632..acd4ebf 100644 --- a/stride/problem/domain.py +++ b/stride/problem/domain.py @@ -332,11 +332,6 @@ def __init__(self, nodes=None, cells=None, cell_tags=None, facet_tags=None): self.cell_tags = cell_tags self.facet_tags = facet_tags - # eagerly, as plain tuples: a cached_property would end up in __dict__ and be pickled - # along with the space every time a field travels to a worker - self.origin = tuple(nodes.min(axis=0)) - self.limit = tuple(nodes.max(axis=0)) - @property def num_nodes(self): """ @@ -362,6 +357,14 @@ def size(self): return tuple(each_limit - each_origin for each_limit, each_origin in zip(self.limit, self.origin)) + @cached_property + def origin(self): + return tuple(self.nodes.min(axis=0)) + + @cached_property + def limit(self): + return tuple(self.nodes.max(axis=0)) + def contains_box(self, lower, upper, atol=1e-9): """ Whether the mesh covers an axis-aligned box. diff --git a/stride/tests/conftest.py b/stride/tests/conftest.py index d4b73ad..692cf40 100644 --- a/stride/tests/conftest.py +++ b/stride/tests/conftest.py @@ -4,13 +4,13 @@ The reference mesh is a Kuhn (Freudenthal) tetrahedralisation of a structured box: every hexahedral cell of an ``(nx, ny, nz)`` node grid is split into six tetrahedra. This gives a genuine unstructured mesh (flat node list, explicit -connectivity) without needing DOLFINx or gmsh to be installed, which mirrors -what ``ae_modelling.fem.mesh.make_mesh`` produces via -``dolfinx.mesh.create_box``. +connectivity) without needing DOLFINx or gmsh to be installed, mirroring what +``dolfinx.mesh.create_box`` produces. """ import numpy as np import pytest +from stride.problem.domain import MeshedSpace, Space # Local corner index within a hexahedron is 4*i + 2*j + k, so the six tets of @@ -112,8 +112,6 @@ def tri_mesh(): @pytest.fixture def meshed_space(tetra_mesh): """A 3D MeshedSpace over the reference tetrahedral mesh.""" - from stride.problem.domain import MeshedSpace - nodes, cells = tetra_mesh return MeshedSpace(nodes=nodes, cells=cells) @@ -121,8 +119,6 @@ def meshed_space(tetra_mesh): @pytest.fixture def meshed_space_2d(tri_mesh): """A 2D MeshedSpace over the reference triangular mesh.""" - from stride.problem.domain import MeshedSpace - nodes, cells = tri_mesh return MeshedSpace(nodes=nodes, cells=cells) @@ -133,10 +129,8 @@ def tagged_meshed_space(tetra_mesh): A MeshedSpace whose cells carry two material tags. Cells in the lower half of the box (in z) are tagged ``1``, the rest ``2``, - mirroring the ``cell_tags`` that ``ae_modelling`` reads out of a gmsh file. + mirroring the ``cell_tags`` that come out of a gmsh file. """ - from stride.problem.domain import MeshedSpace - nodes, cells = tetra_mesh centroids = nodes[cells].mean(axis=1) cell_tags = np.where(centroids[:, 2] < 1e-3, 1, 2).astype(np.int32) @@ -147,26 +141,27 @@ def tagged_meshed_space(tetra_mesh): @pytest.fixture def structured_space(): """A conventional structured Space, used for the serialisation regressions.""" - from stride.problem.domain import Space - return Space(shape=(6, 8), spacing=(1e-3, 1e-3), extra=(2, 2), absorbing=(1, 1)) @pytest.fixture -def tissue_properties(): +def material_properties(): """ - Label -> (conductivity, relative permittivity) map. + Label -> material property map. - Values are the defaults from ``ae_modelling.tissue.properties``. + Each label carries two independent scalar properties, ``alpha`` and + ``beta``, so that tests can build more than one field from the same set of + labels. The names and values are arbitrary: what matters is that a label + indexes a set of physical values. """ return { - 0: {'name': 'background', 'sigma': 1e-6, 'eps_r': 1e0}, - 1: {'name': 'grey_matter', 'sigma': 1.52e-1, 'eps_r': 2.19e3}, - 2: {'name': 'white_matter', 'sigma': 9.47e-2, 'eps_r': 7.12e2}, - 3: {'name': 'csf', 'sigma': 2e0, 'eps_r': 1.09e2}, - 4: {'name': 'skull', 'sigma': 2.22e-2, 'eps_r': 1.75e2}, - 5: {'name': 'skin', 'sigma': 4.36e-3, 'eps_r': 1.06e3}, + 0: {'name': 'material_0', 'alpha': 1e-6, 'beta': 1e0}, + 1: {'name': 'material_1', 'alpha': 1e-1, 'beta': 2e3}, + 2: {'name': 'material_2', 'alpha': 2e-1, 'beta': 7e2}, + 3: {'name': 'material_3', 'alpha': 5e-1, 'beta': 1e2}, + 4: {'name': 'material_4', 'alpha': 8e-1, 'beta': 3e2}, + 5: {'name': 'material_5', 'alpha': 1e0, 'beta': 5e2}, } diff --git a/stride/tests/test_meshed_data.py b/stride/tests/test_meshed_data.py index 4afa0c3..df53b44 100644 --- a/stride/tests/test_meshed_data.py +++ b/stride/tests/test_meshed_data.py @@ -17,77 +17,59 @@ ``time_dependent`` / ``slow_time_dependent`` / ``dim`` options as SparseField - padding is a no-op on both - ``alike`` / ``copy`` / ``detach`` carry the meshed grid across -- inherited arithmetic operates on the nodal buffer -- ``clear_grad`` allocates a nodal gradient +- inherited arithmetic operates on the node buffer +- ``clear_grad`` allocates a node gradient """ import numpy as np import pytest - -from stride.problem.domain import Grid - - -def nodal_grid(space): - return Grid(space, None, None) +from stride.problem.data import MeshedData, MeshedField +from stride.problem.domain import Grid, MeshedSpace, Time class TestMeshedDataShape: def test_explicit_shape(self, meshed_space): - from stride.problem.data import MeshedData - - data = MeshedData(name='raw', shape=(27,), grid=nodal_grid(meshed_space)) + data = MeshedData(name='raw', shape=(27,), space=meshed_space) assert tuple(data.shape) == (27,) assert tuple(data.extended_shape) == (27,) def test_shape_inferred_from_space(self, meshed_space): - from stride.problem.data import MeshedData - - data = MeshedData(name='raw', grid=nodal_grid(meshed_space)) + data = MeshedData(name='raw', space=meshed_space) assert tuple(data.shape) == (meshed_space.num_nodes,) def test_shape_inferred_from_data(self, meshed_space): - from stride.problem.data import MeshedData - values = np.arange(27, dtype=np.float32) - data = MeshedData(name='raw', data=values, grid=nodal_grid(meshed_space)) + data = MeshedData(name='raw', data=values, space=meshed_space) assert tuple(data.shape) == (27,) np.testing.assert_array_equal(data.data, values) @pytest.mark.parametrize('num_values', [25, 26, 29]) def test_mismatched_data_length_rejected(self, meshed_space, num_values): - from stride.problem.data import MeshedData - - # A nodal field with the wrong number of values cannot be interpolated + # A node field with the wrong number of values cannot be interpolated # onto the mesh, so this has to fail loudly rather than at solve time. # 25 is the case that matters: the inherited StructuredData.pad_data floors its # pad widths, so an off-by-two would otherwise be silently edge-padded to 27. with pytest.raises(ValueError): MeshedData(name='raw', data=np.zeros(num_values, dtype=np.float32), - grid=nodal_grid(meshed_space)) + space=meshed_space) def test_compression_is_rejected(self, meshed_space): - from stride.problem.data import MeshedData - # maybe_compress crashes outright for buffers of 2.5k-10k float32 elements, and # above that window the decompressed buffer is read-only. with pytest.raises(ValueError, match='[Cc]ompression'): - MeshedData(name='raw', compressed=True, grid=nodal_grid(meshed_space)) + MeshedData(name='raw', compressed=True, space=meshed_space) def test_default_dtype_is_float32(self, meshed_space): - from stride.problem.data import MeshedData - - data = MeshedData(name='raw', grid=nodal_grid(meshed_space)) + data = MeshedData(name='raw', space=meshed_space) assert data.dtype == np.float32 def test_complex_dtype_is_honoured(self, meshed_space): - from stride.problem.data import MeshedData - - data = MeshedData(name='raw', dtype=np.complex128, grid=nodal_grid(meshed_space)) + data = MeshedData(name='raw', dtype=np.complex128, space=meshed_space) data.fill(1 + 2j) assert data.data.dtype == np.complex128 @@ -97,24 +79,18 @@ def test_complex_dtype_is_honoured(self, meshed_space): class TestMeshedDataPadding: def test_inner_is_the_whole_buffer(self, meshed_space): - from stride.problem.data import MeshedData - - data = MeshedData(name='raw', grid=nodal_grid(meshed_space)) + data = MeshedData(name='raw', space=meshed_space) assert data.inner == (slice(0, None),) def test_pad_data_is_a_noop(self, meshed_space): - from stride.problem.data import MeshedData - - data = MeshedData(name='raw', grid=nodal_grid(meshed_space)) + data = MeshedData(name='raw', space=meshed_space) values = np.arange(27, dtype=np.float32) np.testing.assert_array_equal(data.pad_data(values), values) def test_data_and_extended_data_agree(self, meshed_space): - from stride.problem.data import MeshedData - - data = MeshedData(name='raw', grid=nodal_grid(meshed_space)) + data = MeshedData(name='raw', space=meshed_space) data.fill(3.) np.testing.assert_array_equal(data.data, data.extended_data) @@ -122,110 +98,90 @@ def test_data_and_extended_data_agree(self, meshed_space): class TestMeshedFieldShape: - def test_scalar_nodal_field(self, meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=nodal_grid(meshed_space)) + def test_scalar_node_field(self, meshed_space): + field = MeshedField(name='alpha', space=meshed_space) assert tuple(field.shape) == (27,) assert tuple(field.extended_shape) == (27,) assert field.inner == (slice(0, None),) - def test_vector_nodal_field(self, meshed_space): - from stride.problem.data import MeshedField - - # The electric field E = -grad(phi) is the motivating case. - field = MeshedField(name='e_field', dim=3, grid=nodal_grid(meshed_space)) + def test_vector_node_field(self, meshed_space): + field = MeshedField(name='vector_field', dim=3, space=meshed_space) assert tuple(field.shape) == (27, 3) assert field.inner == (slice(0, None), slice(0, None)) def test_time_dependent_field(self, meshed_space): - from stride.problem.data import MeshedField - from stride.problem.domain import Time - time = Time(start=0., step=1e-6, num=11) - field = MeshedField(name='phi', time_dependent=True, + field = MeshedField(name='transient', time_dependent=True, grid=Grid(meshed_space, time, None)) assert tuple(field.shape) == (11, 27) def test_time_dependent_vector_field(self, meshed_space): - from stride.problem.data import MeshedField - from stride.problem.domain import Time - time = Time(start=0., step=1e-6, num=11) - field = MeshedField(name='e_field', dim=3, time_dependent=True, + field = MeshedField(name='vector_field', dim=3, time_dependent=True, grid=Grid(meshed_space, time, None)) assert tuple(field.shape) == (11, 27, 3) def test_dim_defaults_to_scalar(self, meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=nodal_grid(meshed_space)) + field = MeshedField(name='alpha', space=meshed_space) assert field.dim == 1 def test_num_nodes_is_exposed(self, meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=nodal_grid(meshed_space)) + field = MeshedField(name='alpha', space=meshed_space) assert field.num_nodes == meshed_space.num_nodes def test_explicit_shape_overrides_the_grid(self, meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', shape=(5,), grid=nodal_grid(meshed_space)) + field = MeshedField(name='alpha', shape=(5,), space=meshed_space) assert tuple(field.shape) == (5,) class TestMeshedFieldLocation: """ - A meshed field is either nodal or per-cell, and `shape` alone does not say which, + A meshed field is either node or per-cell, and `shape` alone does not say which, so `location` is what keeps the two apart. """ - def test_defaults_to_nodal(self, meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=nodal_grid(meshed_space)) + def test_defaults_to_node(self, meshed_space): + field = MeshedField(name='alpha', space=meshed_space) - assert field.location == 'nodal' + assert field.location == 'node' assert field.num_entities == meshed_space.num_nodes def test_cell_location_sizes_from_cells(self, tagged_meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', location='cell', - grid=nodal_grid(tagged_meshed_space)) + field = MeshedField(name='alpha', location='cell', + space=tagged_meshed_space) assert field.location == 'cell' assert field.num_entities == tagged_meshed_space.num_cells assert tuple(field.shape) == (48,) def test_from_cell_tags_reports_cell_location(self, tagged_meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField.from_cell_tags({1: 0.1, 2: 0.5}, name='sigma', - grid=nodal_grid(tagged_meshed_space)) + field = MeshedField.from_cell_tags({1: 0.1, 2: 0.5}, name='alpha', + space=tagged_meshed_space) assert field.location == 'cell' - def test_invalid_location_rejected(self, meshed_space): - from stride.problem.data import MeshedField + def test_from_cell_tags_values_are_consistent(self, tagged_meshed_space): + field = MeshedField.from_cell_tags({1: 0.1, 2: 0.5}, name='alpha', + space=tagged_meshed_space) + np.testing.assert_allclose(field.data, + np.where(tagged_meshed_space.cell_tags == 1, 0.1, 0.5)) + + def test_invalid_location_rejected(self, meshed_space): with pytest.raises(ValueError): - MeshedField(name='sigma', location='facet', grid=nodal_grid(meshed_space)) + MeshedField(name='alpha', location='facet', space=meshed_space) def test_alike_preserves_location(self, tagged_meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField.from_cell_tags({1: 0.1, 2: 0.5}, name='sigma', - grid=nodal_grid(tagged_meshed_space)) - other = field.alike(name='eps') + field = MeshedField.from_cell_tags({1: 0.1, 2: 0.5}, name='alpha', + space=tagged_meshed_space) + other = field.alike(name='beta') assert other.location == 'cell' assert tuple(other.shape) == (48,) @@ -234,11 +190,9 @@ def test_alike_preserves_location(self, tagged_meshed_space): class TestMeshedFieldCopying: def test_alike_keeps_the_meshed_grid(self, meshed_space): - from stride.problem.data import MeshedField - from stride.problem.domain import MeshedSpace - field = MeshedField(name='sigma', grid=nodal_grid(meshed_space)) - other = field.alike(name='eps') + field = MeshedField(name='alpha', space=meshed_space) + other = field.alike(name='beta') assert isinstance(other.space, MeshedSpace) assert other.space is meshed_space @@ -246,9 +200,7 @@ def test_alike_keeps_the_meshed_grid(self, meshed_space): assert other.dtype == field.dtype def test_copy_duplicates_the_buffer(self, meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=nodal_grid(meshed_space)) + field = MeshedField(name='alpha', space=meshed_space) field.fill(2.) cpy = field.copy() @@ -258,17 +210,13 @@ def test_copy_duplicates_the_buffer(self, meshed_space): np.testing.assert_allclose(cpy.data, 5.) def test_copy_keeps_vector_shape(self, meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField(name='e_field', dim=3, grid=nodal_grid(meshed_space)) + field = MeshedField(name='vector_field', dim=3, space=meshed_space) field.fill(1.) assert tuple(field.copy().shape) == (27, 3) def test_detach_keeps_shape_and_data(self, meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=nodal_grid(meshed_space)) + field = MeshedField(name='alpha', space=meshed_space) field.fill(4.) detached = field.detach() @@ -280,9 +228,7 @@ def test_detach_keeps_shape_and_data(self, meshed_space): class TestMeshedFieldArithmetic: def test_add(self, meshed_space): - from stride.problem.data import MeshedField - - a = MeshedField(name='a', grid=nodal_grid(meshed_space)) + a = MeshedField(name='a', space=meshed_space) a.fill(1.) b = a.copy() b.fill(2.) @@ -290,34 +236,26 @@ def test_add(self, meshed_space): np.testing.assert_allclose((a + b).data, 3.) def test_multiply_by_scalar(self, meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField(name='a', grid=nodal_grid(meshed_space)) + field = MeshedField(name='a', space=meshed_space) field.fill(3.) np.testing.assert_allclose((field * 2).data, 6.) def test_in_place_add(self, meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField(name='a', grid=nodal_grid(meshed_space)) + field = MeshedField(name='a', space=meshed_space) field.fill(1.) field += 1. np.testing.assert_allclose(field.data, 2.) def test_operations_preserve_node_count(self, meshed_space): - from stride.problem.data import MeshedField - - a = MeshedField(name='a', grid=nodal_grid(meshed_space)) + a = MeshedField(name='a', space=meshed_space) a.fill(1.) assert tuple((a * 2 + a).shape) == (27,) def test_elementwise_over_nodes(self, meshed_space): - from stride.problem.data import MeshedField - - a = MeshedField(name='a', grid=nodal_grid(meshed_space)) + a = MeshedField(name='a', space=meshed_space) a.allocate() a.data[:] = np.arange(27) @@ -326,10 +264,8 @@ def test_elementwise_over_nodes(self, meshed_space): class TestMeshedFieldGradient: - def test_clear_grad_allocates_a_nodal_gradient(self, meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=nodal_grid(meshed_space), + def test_clear_grad_allocates_a_node_gradient(self, meshed_space): + field = MeshedField(name='alpha', space=meshed_space, needs_grad=True) field.clear_grad() @@ -338,9 +274,7 @@ def test_clear_grad_allocates_a_nodal_gradient(self, meshed_space): np.testing.assert_allclose(field.grad.data, 0.) def test_gradient_has_a_preconditioner(self, meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=nodal_grid(meshed_space), + field = MeshedField(name='alpha', space=meshed_space, needs_grad=True) field.clear_grad() @@ -348,9 +282,7 @@ def test_gradient_has_a_preconditioner(self, meshed_space): assert tuple(field.grad.prec.shape) == (27,) def test_clear_grad_is_a_noop_without_needs_grad(self, meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=nodal_grid(meshed_space)) + field = MeshedField(name='alpha', space=meshed_space) field.clear_grad() assert field.grad is None diff --git a/stride/tests/test_meshed_medium.py b/stride/tests/test_meshed_medium.py index de4a017..97cb3a0 100644 --- a/stride/tests/test_meshed_medium.py +++ b/stride/tests/test_meshed_medium.py @@ -1,67 +1,43 @@ """ Tests for loading a medium onto a MeshedSpace. - -This is the stride-side port of the medium path in ``ae-modelling``: - -- ``ae_modelling.tissue.nifti.NIfTILabelSampler.sample_labels`` samples a voxel - segmentation at arbitrary coordinates via the NIfTI affine - -> ``MeshedSpace.sample_labels(volume, affine)``, evaluated at mesh nodes -- ``sample_sigma`` / ``sample_eps`` map labels through a lookup table - -> ``MeshedField.from_labels(labels, lut, ...)`` -- ``ae_modelling.fem.space.DielectricSpace.admittivity`` combines the two into - ``sigma + 1j * omega * eps_r * EPS0`` - -> ordinary MeshedField arithmetic, so nothing new is needed for it -- ``ae_modelling.fem.interpolate.interpolate_medium`` also supports a - label -> value map applied per cell tag - -> ``MeshedField.from_cell_tags(mapping, ...)`` - -The sampling here is nearest-voxel, matching ``world_to_voxel``'s ``np.rint`` -plus clipping. Note that stride must not grow a nibabel dependency: the volume -and its affine are passed in as plain arrays, and reading the ``.nii`` file -stays the caller's job. """ import numpy as np import pytest -from stride.problem.domain import Grid - - -# ae_modelling.fem.space.EPS0 -EPS0 = 8.854e-12 +from stride.problem.domain import Grid, MeshedSpace +from stride.problem.data import MeshedField +from stride.problem.medium import Medium +def value_lut(material_properties, prop='alpha'): + """ + One property of the material table, as an array indexed by label. -def sigma_lut(tissue_properties): - """Conductivity indexed by label, as get_tissue_sigma_array does.""" - lut = np.zeros(max(tissue_properties) + 1) - for label, properties in tissue_properties.items(): - lut[label] = properties['sigma'] - return lut - + This is the array form of a label -> value lookup: ``lut[label]`` is the + value of ``prop`` for that material. -def eps_lut(tissue_properties): - """Relative permittivity indexed by label, as get_tissue_eps_array does.""" - lut = np.zeros(max(tissue_properties) + 1) - for label, properties in tissue_properties.items(): - lut[label] = properties['eps_r'] + """ + lut = np.zeros(max(material_properties) + 1) + for label, properties in material_properties.items(): + lut[label] = properties[prop] return lut -class TestTissueLookupTables: +class TestValueLookupTables: """Guards on the fixture itself, so the sampling tests read unambiguously.""" - def test_sigma_lut_is_indexed_by_label(self, tissue_properties): - lut = sigma_lut(tissue_properties) + def test_lut_is_indexed_by_label(self, material_properties): + lut = value_lut(material_properties) assert lut.shape == (6,) - assert lut[1] == pytest.approx(1.52e-1) - assert lut[3] == pytest.approx(2e0) + assert lut[1] == pytest.approx(1e-1) + assert lut[3] == pytest.approx(5e-1) - def test_eps_lut_is_indexed_by_label(self, tissue_properties): - lut = eps_lut(tissue_properties) + def test_lut_selects_the_requested_property(self, material_properties): + lut = value_lut(material_properties, prop='beta') - assert lut[1] == pytest.approx(2.19e3) - assert lut[3] == pytest.approx(1.09e2) + assert lut[1] == pytest.approx(2e3) + assert lut[3] == pytest.approx(1e2) class TestSampleLabelsAtNodes: @@ -86,8 +62,6 @@ def test_labels_follow_the_volume(self, meshed_space, label_volume): np.testing.assert_array_equal(labels[z > 1.5e-3], 3) def test_affine_translation_is_applied(self, tetra_mesh, label_volume): - from stride.problem.domain import MeshedSpace - volume, _ = label_volume nodes, cells = tetra_mesh space = MeshedSpace(nodes=nodes, cells=cells) @@ -102,8 +76,6 @@ def test_affine_translation_is_applied(self, tetra_mesh, label_volume): np.testing.assert_array_equal(labels, 3) def test_identity_affine_treats_nodes_as_voxel_indices(self, tetra_mesh, label_volume): - from stride.problem.domain import MeshedSpace - volume, _ = label_volume nodes, cells = tetra_mesh @@ -116,13 +88,10 @@ def test_identity_affine_treats_nodes_as_voxel_indices(self, tetra_mesh, label_v np.testing.assert_array_equal(labels[z > 1.5], 3) def test_out_of_volume_nodes_are_clipped(self, tetra_mesh, label_volume): - from stride.problem.domain import MeshedSpace - volume, affine = label_volume nodes, cells = tetra_mesh - # Push the mesh well past the 4x4x4 voxel volume. world_to_voxel clips - # rather than raising, so edge nodes take the nearest in-range label. + # Push the mesh well past the 4x4x4 voxel volume. space = MeshedSpace(nodes=nodes + 1e-1, cells=cells) labels = space.sample_labels(volume, affine=affine) @@ -143,102 +112,84 @@ def test_2d_volume_sampling(self, meshed_space_2d): class TestFieldFromLabels: - def test_conductivity_from_labels(self, meshed_space, label_volume, tissue_properties): - from stride.problem.data import MeshedField - + def test_field_from_labels(self, meshed_space, label_volume, material_properties): volume, affine = label_volume labels = meshed_space.sample_labels(volume, affine=affine) - sigma = MeshedField.from_labels(labels, sigma_lut(tissue_properties), - name='sigma', - grid=Grid(meshed_space, None, None)) + field = MeshedField.from_labels(labels, value_lut(material_properties), + name='alpha', space=meshed_space) - assert tuple(sigma.shape) == (meshed_space.num_nodes,) + assert tuple(field.shape) == (meshed_space.num_nodes,) z = meshed_space.nodes[:, 2] - np.testing.assert_allclose(sigma.data[z < 1.5e-3], 1.52e-1, rtol=1e-6) - np.testing.assert_allclose(sigma.data[z > 1.5e-3], 2e0, rtol=1e-6) - - def test_permittivity_from_labels(self, meshed_space, label_volume, tissue_properties): - from stride.problem.data import MeshedField + np.testing.assert_allclose(field.data[z < 1.5e-3], 1e-1, rtol=1e-6) + np.testing.assert_allclose(field.data[z > 1.5e-3], 5e-1, rtol=1e-6) + def test_a_second_property_maps_over_the_same_labels(self, meshed_space, label_volume, + material_properties): volume, affine = label_volume labels = meshed_space.sample_labels(volume, affine=affine) - eps = MeshedField.from_labels(labels, eps_lut(tissue_properties), - name='eps', grid=Grid(meshed_space, None, None)) + field = MeshedField.from_labels(labels, value_lut(material_properties, prop='beta'), + name='beta', space=meshed_space) z = meshed_space.nodes[:, 2] - np.testing.assert_allclose(eps.data[z < 1.5e-3], 2.19e3, rtol=1e-6) - np.testing.assert_allclose(eps.data[z > 1.5e-3], 1.09e2, rtol=1e-6) + np.testing.assert_allclose(field.data[z < 1.5e-3], 2e3, rtol=1e-6) + np.testing.assert_allclose(field.data[z > 1.5e-3], 1e2, rtol=1e-6) def test_lut_may_be_a_mapping(self, meshed_space, label_volume): - from stride.problem.data import MeshedField - volume, affine = label_volume labels = meshed_space.sample_labels(volume, affine=affine) - # Sparse or non-contiguous label sets are common in segmentations, so a - # dict has to work as well as an indexable array. - field = MeshedField.from_labels(labels, {1: 10., 3: 30.}, name='sigma', - grid=Grid(meshed_space, None, None)) + # Sparse or non-contiguous label sets are common in labelled volumes, so + # a dict has to work as well as an indexable array. + field = MeshedField.from_labels(labels, {1: 10., 3: 30.}, name='alpha', + space=meshed_space) z = meshed_space.nodes[:, 2] np.testing.assert_allclose(field.data[z < 1.5e-3], 10.) np.testing.assert_allclose(field.data[z > 1.5e-3], 30.) def test_unmapped_label_raises(self, meshed_space, label_volume): - from stride.problem.data import MeshedField - volume, affine = label_volume labels = meshed_space.sample_labels(volume, affine=affine) - # Silently defaulting an unmapped tissue to zero conductivity would - # produce a plausible-looking but wrong solve. + # Silently defaulting an unmapped label to a zero property would produce + # a plausible-looking but wrong solve. with pytest.raises(KeyError): - MeshedField.from_labels(labels, {1: 10.}, name='sigma', - grid=Grid(meshed_space, None, None)) - - def test_label_count_must_match_nodes(self, meshed_space, tissue_properties): - from stride.problem.data import MeshedField + MeshedField.from_labels(labels, {1: 10.}, name='alpha', + space=meshed_space) + def test_label_count_must_match_nodes(self, meshed_space, material_properties): with pytest.raises(ValueError): MeshedField.from_labels(np.ones(5, dtype=np.int64), - sigma_lut(tissue_properties), - name='sigma', - grid=Grid(meshed_space, None, None)) - - def test_label_past_the_end_of_an_array_lut_raises(self, meshed_space, tissue_properties): - from stride.problem.data import MeshedField + value_lut(material_properties), + name='alpha', + space=meshed_space) + def test_label_past_the_end_of_an_array_lut_raises(self, meshed_space, material_properties): labels = np.full(meshed_space.num_nodes, 9, dtype=np.int64) # An array lut would raise IndexError here rather than KeyError; normalise it so both # lut forms report an unmapped label the same way. with pytest.raises(KeyError): - MeshedField.from_labels(labels, sigma_lut(tissue_properties), name='sigma', - grid=Grid(meshed_space, None, None)) - - def test_negative_label_raises(self, meshed_space, tissue_properties): - from stride.problem.data import MeshedField + MeshedField.from_labels(labels, value_lut(material_properties), name='alpha', + space=meshed_space) + def test_negative_label_raises(self, meshed_space, material_properties): labels = np.full(meshed_space.num_nodes, -1, dtype=np.int64) - # Segmentations do use -1 as a sentinel, and a negative index would otherwise + # Labelled volumes do use -1 as a sentinel, and a negative index would otherwise # silently read from the end of the lookup table. with pytest.raises(KeyError): - MeshedField.from_labels(labels, sigma_lut(tissue_properties), name='sigma', - grid=Grid(meshed_space, None, None)) + MeshedField.from_labels(labels, value_lut(material_properties), name='alpha', + space=meshed_space) class TestFieldFromCellTags: - """The dict branch of ae_modelling.fem.interpolate.interpolate_medium.""" - def test_values_are_assigned_per_cell(self, tagged_meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField.from_cell_tags({1: 0.1, 2: 0.5}, name='sigma', - grid=Grid(tagged_meshed_space, None, None)) + field = MeshedField.from_cell_tags({1: 0.1, 2: 0.5}, name='alpha', + space=tagged_meshed_space) # One value per cell, as for a DG-0 medium function. assert tuple(field.shape) == (tagged_meshed_space.num_cells,) @@ -248,67 +199,65 @@ def test_values_are_assigned_per_cell(self, tagged_meshed_space): np.testing.assert_allclose(field.data[tags == 2], 0.5) def test_missing_tag_in_mapping_raises(self, tagged_meshed_space): - from stride.problem.data import MeshedField - - # Mirrors the assertion in load_mesh that sigma's keys must cover the - # cell tag labels. + # A per-material mapping has to cover every cell tag label present on + # the mesh, or some cells would be left without a value. with pytest.raises(KeyError): - MeshedField.from_cell_tags({1: 0.1}, name='sigma', - grid=Grid(tagged_meshed_space, None, None)) + MeshedField.from_cell_tags({1: 0.1}, name='alpha', + space=tagged_meshed_space) def test_requires_cell_tags_on_the_space(self, meshed_space): - from stride.problem.data import MeshedField - with pytest.raises(ValueError): - MeshedField.from_cell_tags({1: 0.1}, name='sigma', - grid=Grid(meshed_space, None, None)) + MeshedField.from_cell_tags({1: 0.1}, name='alpha', space=meshed_space) -class TestAdmittivity: - """DielectricSpace.admittivity, expressed with MeshedField arithmetic.""" +class TestComplexCombination: + """ + Two real node fields combined into a single complex one. - def _fields(self, meshed_space, label_volume, tissue_properties): - from stride.problem.data import MeshedField + A frequency-domain medium coefficient is typically built this way: one field + is the real part, another is scaled onto the imaginary part. + """ + def _fields(self, meshed_space, label_volume, material_properties): volume, affine = label_volume labels = meshed_space.sample_labels(volume, affine=affine) grid = Grid(meshed_space, None, None) - sigma = MeshedField.from_labels(labels, sigma_lut(tissue_properties), - name='sigma', dtype=np.complex128, grid=grid) - eps = MeshedField.from_labels(labels, eps_lut(tissue_properties), - name='eps', dtype=np.complex128, grid=grid) - return sigma, eps + alpha = MeshedField.from_labels(labels, value_lut(material_properties), + name='alpha', dtype=np.complex128, grid=grid) + beta = MeshedField.from_labels(labels, value_lut(material_properties, prop='beta'), + name='beta', dtype=np.complex128, grid=grid) + return alpha, beta - def test_admittivity_is_sigma_plus_j_omega_eps(self, meshed_space, label_volume, - tissue_properties): - omega = 2 * np.pi * 5e5 - sigma, eps = self._fields(meshed_space, label_volume, tissue_properties) + def test_combination_is_alpha_plus_j_scale_beta(self, meshed_space, label_volume, + material_properties): + scale = 1e-4 + alpha, beta = self._fields(meshed_space, label_volume, material_properties) - admittivity = sigma + eps * (1j * omega * EPS0) + combined = alpha + beta * (1j * scale) - assert admittivity.data.dtype == np.complex128 + assert combined.data.dtype == np.complex128 z = meshed_space.nodes[:, 2] - expected = 1.52e-1 + 1j * omega * 2.19e3 * EPS0 - np.testing.assert_allclose(admittivity.data[z < 1.5e-3], expected, rtol=1e-6) + expected = 1e-1 + 1j * scale * 2e3 + np.testing.assert_allclose(combined.data[z < 1.5e-3], expected, rtol=1e-6) - def test_real_part_is_the_conductivity(self, meshed_space, label_volume, - tissue_properties): - omega = 2 * np.pi * 5e5 - sigma, eps = self._fields(meshed_space, label_volume, tissue_properties) + def test_real_part_is_the_first_field(self, meshed_space, label_volume, + material_properties): + scale = 1e-4 + alpha, beta = self._fields(meshed_space, label_volume, material_properties) - admittivity = sigma + eps * (1j * omega * EPS0) + combined = alpha + beta * (1j * scale) - np.testing.assert_allclose(admittivity.data.real, sigma.data.real, rtol=1e-6) + np.testing.assert_allclose(combined.data.real, alpha.data.real, rtol=1e-6) - def test_zero_frequency_reduces_to_conductivity(self, meshed_space, label_volume, - tissue_properties): - sigma, eps = self._fields(meshed_space, label_volume, tissue_properties) + def test_zero_scale_reduces_to_the_first_field(self, meshed_space, label_volume, + material_properties): + alpha, beta = self._fields(meshed_space, label_volume, material_properties) - admittivity = sigma + eps * (1j * 0. * EPS0) + combined = alpha + beta * (1j * 0.) - np.testing.assert_allclose(admittivity.data, sigma.data, rtol=1e-6) + np.testing.assert_allclose(combined.data, alpha.data, rtol=1e-6) class TestMeshedMedium: @@ -318,44 +267,38 @@ class TestMeshedMedium: def project(self, tmp_path): return {'path': str(tmp_path), 'project_name': 'meshed_medium'} - def _medium(self, meshed_space, label_volume, tissue_properties): - from stride.problem.data import MeshedField - from stride.problem.medium import Medium - + def _medium(self, meshed_space, label_volume, material_properties): volume, affine = label_volume labels = meshed_space.sample_labels(volume, affine=affine) grid = Grid(meshed_space, None, None) medium = Medium(grid=grid) - medium.add(MeshedField.from_labels(labels, sigma_lut(tissue_properties), - name='sigma', grid=grid)) - medium.add(MeshedField.from_labels(labels, eps_lut(tissue_properties), - name='eps', grid=grid)) + medium.add(MeshedField.from_labels(labels, value_lut(material_properties), + name='alpha', grid=grid)) + medium.add(MeshedField.from_labels(labels, value_lut(material_properties, prop='beta'), + name='beta', grid=grid)) return medium def test_fields_are_accessible_by_name(self, meshed_space, label_volume, - tissue_properties): - medium = self._medium(meshed_space, label_volume, tissue_properties) + material_properties): + medium = self._medium(meshed_space, label_volume, material_properties) - assert set(medium.fields) == {'sigma', 'eps'} - assert tuple(medium.sigma.shape) == (meshed_space.num_nodes,) - assert tuple(medium['eps'].shape) == (meshed_space.num_nodes,) + assert set(medium.fields) == {'alpha', 'beta'} + assert tuple(medium.alpha.shape) == (meshed_space.num_nodes,) + assert tuple(medium['beta'].shape) == (meshed_space.num_nodes,) - def test_medium_round_trip(self, meshed_space, label_volume, tissue_properties, + def test_medium_round_trip(self, meshed_space, label_volume, material_properties, project): - from stride.problem.data import MeshedField - from stride.problem.domain import MeshedSpace - from stride.problem.medium import Medium - medium = self._medium(meshed_space, label_volume, tissue_properties) + medium = self._medium(meshed_space, label_volume, material_properties) medium.dump(**project) loaded = Medium() - loaded.add(MeshedField(name='sigma')) - loaded.add(MeshedField(name='eps')) + loaded.add(MeshedField(name='alpha')) + loaded.add(MeshedField(name='beta')) loaded.load(**project) - assert isinstance(loaded.sigma.space, MeshedSpace) - assert loaded.sigma.space.num_nodes == meshed_space.num_nodes - np.testing.assert_allclose(loaded.sigma.data, medium.sigma.data) - np.testing.assert_allclose(loaded.eps.data, medium.eps.data) + assert isinstance(loaded.alpha.space, MeshedSpace) + assert loaded.alpha.space.num_nodes == meshed_space.num_nodes + np.testing.assert_allclose(loaded.alpha.data, medium.alpha.data) + np.testing.assert_allclose(loaded.beta.data, medium.beta.data) diff --git a/stride/tests/test_meshed_serialisation.py b/stride/tests/test_meshed_serialisation.py index e74b70d..0f67e9f 100644 --- a/stride/tests/test_meshed_serialisation.py +++ b/stride/tests/test_meshed_serialisation.py @@ -21,8 +21,8 @@ import numpy as np import pytest -from stride.problem.domain import Grid, MeshedSpace, Space - +from stride.problem.domain import Grid, Time, Space, MeshedSpace +from stride.problem.data import ScalarField, MeshedField @pytest.fixture def project(tmp_path): @@ -34,9 +34,7 @@ class TestGridDescriptionStructured: """The pre-existing Space path must keep working untouched.""" def test_structured_space_keys(self, structured_space, project): - from stride.problem.data import ScalarField - - field = ScalarField(name='vp_field', grid=Grid(structured_space, None, None)) + field = ScalarField(name='vp_field', space=structured_space) description = field.grid_description() assert set(description['space']) == {'shape', 'spacing', 'extra', 'absorbing'} @@ -44,16 +42,12 @@ def test_structured_space_keys(self, structured_space, project): assert tuple(description['space']['extra']) == (2, 2) def test_structured_space_has_no_mesh_keys(self, structured_space): - from stride.problem.data import ScalarField - - field = ScalarField(name='vp_field', grid=Grid(structured_space, None, None)) + field = ScalarField(name='vp_field', space=structured_space) assert 'nodes' not in field.grid_description()['space'] def test_structured_round_trip(self, structured_space, project): - from stride.problem.data import ScalarField - - field = ScalarField(name='vp_field', grid=Grid(structured_space, None, None)) + field = ScalarField(name='vp_field', space=structured_space) field.fill(1500.) field.dump(**project) @@ -71,18 +65,14 @@ def test_structured_round_trip(self, structured_space, project): class TestGridDescriptionMeshed: def test_meshed_space_keys(self, meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=Grid(meshed_space, None, None)) + field = MeshedField(name='alpha', space=meshed_space) space_description = field.grid_description()['space'] assert 'nodes' in space_description assert 'cells' in space_description def test_meshed_space_omits_structured_keys(self, meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=Grid(meshed_space, None, None)) + field = MeshedField(name='alpha', space=meshed_space) space_description = field.grid_description()['space'] # 'shape' under a meshed space would make the load-time branch pick the @@ -91,18 +81,14 @@ def test_meshed_space_omits_structured_keys(self, meshed_space): assert 'spacing' not in space_description def test_nodes_are_the_node_table(self, meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=Grid(meshed_space, None, None)) + field = MeshedField(name='alpha', space=meshed_space) nodes = np.asarray(field.grid_description()['space']['nodes']) assert nodes.shape == (27, 3) np.testing.assert_allclose(nodes, meshed_space.nodes) def test_cell_tags_are_included_when_present(self, tagged_meshed_space): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=Grid(tagged_meshed_space, None, None)) + field = MeshedField(name='alpha', space=tagged_meshed_space) space_description = field.grid_description()['space'] assert 'cell_tags' in space_description @@ -112,12 +98,10 @@ def test_cell_tags_are_included_when_present(self, tagged_meshed_space): ) def test_unknown_space_type_raises(self, monkeypatch): - from stride.problem.data import MeshedField - class NotASpace: pass - field = MeshedField(name='sigma', shape=(4,)) + field = MeshedField(name='alpha', shape=(4,)) field.grid.space = NotASpace() # A bare `raise Exception` here would be indistinguishable from a bug @@ -129,26 +113,22 @@ class NotASpace: class TestMeshedRoundTrip: def test_space_type_survives(self, meshed_space, project): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=Grid(meshed_space, None, None)) + field = MeshedField(name='alpha', space=meshed_space) field.fill(0.152) field.dump(**project) - loaded = MeshedField(name='sigma') + loaded = MeshedField(name='alpha') loaded.load(**project) assert isinstance(loaded.space, MeshedSpace) assert not isinstance(loaded.space, Space) def test_nodes_survive(self, meshed_space, project): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=Grid(meshed_space, None, None)) + field = MeshedField(name='alpha', space=meshed_space) field.fill(0.152) field.dump(**project) - loaded = MeshedField(name='sigma') + loaded = MeshedField(name='alpha') loaded.load(**project) assert loaded.space.num_nodes == 27 @@ -156,13 +136,11 @@ def test_nodes_survive(self, meshed_space, project): np.testing.assert_allclose(loaded.space.nodes, meshed_space.nodes) def test_cells_survive(self, meshed_space, project): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=Grid(meshed_space, None, None)) + field = MeshedField(name='alpha', space=meshed_space) field.fill(0.152) field.dump(**project) - loaded = MeshedField(name='sigma') + loaded = MeshedField(name='alpha') loaded.load(**project) # Connectivity is what makes the node table a mesh; without it the @@ -171,67 +149,57 @@ def test_cells_survive(self, meshed_space, project): np.testing.assert_array_equal(loaded.space.cells, meshed_space.cells) def test_cell_tags_survive(self, tagged_meshed_space, project): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=Grid(tagged_meshed_space, None, None)) + field = MeshedField(name='alpha', space=tagged_meshed_space) field.fill(1.) field.dump(**project) - loaded = MeshedField(name='sigma') + loaded = MeshedField(name='alpha') loaded.load(**project) np.testing.assert_array_equal(loaded.space.cell_tags, tagged_meshed_space.cell_tags) def test_data_survives(self, meshed_space, project): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=Grid(meshed_space, None, None)) + field = MeshedField(name='alpha', space=meshed_space) field.allocate() field.data[:] = np.arange(27) field.dump(**project) - loaded = MeshedField(name='sigma') + loaded = MeshedField(name='alpha') loaded.load(**project) assert tuple(loaded.shape) == (27,) np.testing.assert_allclose(loaded.data, np.arange(27)) def test_bounds_survive(self, meshed_space, project): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=Grid(meshed_space, None, None)) + field = MeshedField(name='alpha', space=meshed_space) field.fill(1.) field.dump(**project) - loaded = MeshedField(name='sigma') + loaded = MeshedField(name='alpha') loaded.load(**project) np.testing.assert_allclose(loaded.space.origin, meshed_space.origin) np.testing.assert_allclose(loaded.space.limit, meshed_space.limit) def test_vector_field_round_trip(self, meshed_space, project): - from stride.problem.data import MeshedField - - field = MeshedField(name='e_field', dim=3, grid=Grid(meshed_space, None, None)) + field = MeshedField(name='vector_field', dim=3, space=meshed_space) field.allocate() field.data[:] = np.arange(27 * 3).reshape(27, 3) field.dump(**project) - loaded = MeshedField(name='e_field', dim=3) + loaded = MeshedField(name='vector_field', dim=3) loaded.load(**project) assert tuple(loaded.shape) == (27, 3) np.testing.assert_allclose(loaded.data, np.arange(27 * 3).reshape(27, 3)) def test_two_dimensional_mesh_round_trip(self, meshed_space_2d, project): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=Grid(meshed_space_2d, None, None)) + field = MeshedField(name='alpha', space=meshed_space_2d) field.fill(1.) field.dump(**project) - loaded = MeshedField(name='sigma') + loaded = MeshedField(name='alpha') loaded.load(**project) assert loaded.space.dim == 2 @@ -239,45 +207,38 @@ def test_two_dimensional_mesh_round_trip(self, meshed_space_2d, project): np.testing.assert_allclose(loaded.space.nodes, meshed_space_2d.nodes) def test_existing_grid_is_not_overwritten(self, meshed_space, project): - from stride.problem.data import MeshedField - - field = MeshedField(name='sigma', grid=Grid(meshed_space, None, None)) + field = MeshedField(name='alpha', space=meshed_space) field.fill(1.) field.dump(**project) - loaded = MeshedField(name='sigma', grid=Grid(meshed_space, None, None)) + loaded = MeshedField(name='alpha', space=meshed_space) loaded.load(**project) # GriddedSaved.load only builds a space when the instance has none. assert loaded.space is meshed_space def test_cell_location_round_trip(self, tagged_meshed_space, project): - from stride.problem.data import MeshedField - - field = MeshedField.from_cell_tags({1: 0.1, 2: 0.5}, name='sigma', - grid=Grid(tagged_meshed_space, None, None)) + field = MeshedField.from_cell_tags({1: 0.1, 2: 0.5}, name='alpha', + space=tagged_meshed_space) field.dump(**project) - loaded = MeshedField(name='sigma') + loaded = MeshedField(name='alpha') loaded.load(**project) - # Without `location` on the description the reloaded field would default to nodal + # Without `location` on the description the reloaded field would default to node # and disagree with its own shape. assert loaded.location == 'cell' assert tuple(loaded.shape) == (48,) np.testing.assert_allclose(loaded.data, field.data) def test_time_dependent_round_trip(self, meshed_space, project): - from stride.problem.data import MeshedField - from stride.problem.domain import Time - time = Time(start=0., step=1e-6, num=5) - field = MeshedField(name='phi', time_dependent=True, + field = MeshedField(name='transient', time_dependent=True, grid=Grid(meshed_space, time, None)) field.fill(2.) field.dump(**project) - loaded = MeshedField(name='phi', time_dependent=True) + loaded = MeshedField(name='transient', time_dependent=True) loaded.load(**project) assert isinstance(loaded.space, MeshedSpace) diff --git a/stride/tests/test_meshed_space.py b/stride/tests/test_meshed_space.py index 9e30d91..557870b 100644 --- a/stride/tests/test_meshed_space.py +++ b/stride/tests/test_meshed_space.py @@ -2,8 +2,7 @@ Tests for MeshedSpace (stride/problem/domain.py). MeshedSpace is the unstructured counterpart to Space: instead of a shape and a -spacing it is defined by an explicit node list and cell connectivity, ported -from ``ae_modelling.fem.mesh.MeshDomain`` and ``ae_modelling.fem.space``. +spacing it is defined by an explicit node list and cell connectivity. Contract under test: @@ -12,12 +11,12 @@ - ``num_nodes`` / ``num_cells`` - ``origin`` / ``limit`` / ``size`` come from the node bounding box, the analogue of Space's origin/limit/size -- ``shape`` is ``(num_nodes,)`` -- the shape of a scalar nodal field -- with +- ``shape`` is ``(num_nodes,)`` -- the shape of a scalar node field -- with ``extended_shape == shape``, ``extra`` and ``absorbing`` all-zero and ``inner == (slice(0, None),)``, so that the StructuredData and GriddedSaved machinery that reads those attributes keeps working -- ``contains_box(lower, upper)`` is the mesh-covers-the-grid check that - ``ae_modelling.fem.mesh.attach_mesh`` performs with assertions +- ``contains_box(lower, upper)`` is the mesh-covers-the-grid check performed + when a mesh is attached to a problem grid - ``resample`` raises, because a mesh has no spacing to resample onto - ``from_dolfinx`` adapts an in-memory DOLFINx mesh (skipped without DOLFINx) """ @@ -27,6 +26,10 @@ from .conftest import box_tetra_mesh +from stride.problem.domain import MeshedSpace +from stride.problem.data import MeshedField, ScalarField + + try: import dolfinx # noqa: F401 @@ -52,8 +55,6 @@ def test_counts(self, meshed_space): assert meshed_space.num_cells == 48 def test_num_cells_is_zero_without_connectivity(self, tetra_mesh): - from stride.problem.domain import MeshedSpace - nodes, _ = tetra_mesh space = MeshedSpace(nodes=nodes) @@ -62,8 +63,6 @@ def test_num_cells_is_zero_without_connectivity(self, tetra_mesh): assert space.num_nodes == 27 def test_nodes_stored_as_float64(self, tetra_mesh): - from stride.problem.domain import MeshedSpace - nodes, cells = tetra_mesh space = MeshedSpace(nodes=nodes.astype(np.float32), cells=cells) @@ -83,21 +82,15 @@ def test_cell_tags_default_to_none(self, meshed_space): assert meshed_space.facet_tags is None def test_ragged_nodes_rejected(self): - from stride.problem.domain import MeshedSpace - with pytest.raises(ValueError): MeshedSpace(nodes=np.zeros(10)) def test_unsupported_dimensionality_rejected(self): - from stride.problem.domain import MeshedSpace - - # ae_modelling.fem.mesh.make_mesh only handles dim 2 and 3. + # Only 2D and 3D meshes are supported. with pytest.raises(ValueError): MeshedSpace(nodes=np.zeros((10, 4))) def test_out_of_range_cell_indices_rejected(self, tetra_mesh): - from stride.problem.domain import MeshedSpace - nodes, cells = tetra_mesh broken = cells.copy() broken[0, 0] = len(nodes) @@ -106,8 +99,6 @@ def test_out_of_range_cell_indices_rejected(self, tetra_mesh): MeshedSpace(nodes=nodes, cells=broken) def test_cell_tags_length_must_match_cells(self, tetra_mesh): - from stride.problem.domain import MeshedSpace - nodes, cells = tetra_mesh with pytest.raises(ValueError): @@ -124,8 +115,6 @@ def test_size_is_the_extent(self, meshed_space): np.testing.assert_allclose(meshed_space.size, (2e-3, 2e-3, 2e-3)) def test_offset_origin_is_respected(self): - from stride.problem.domain import MeshedSpace - nodes, cells = box_tetra_mesh(shape=(3, 3, 3), spacing=(1e-3, 1e-3, 1e-3), origin=(-5e-3, 1e-3, 0.)) space = MeshedSpace(nodes=nodes, cells=cells) @@ -157,34 +146,27 @@ def test_has_no_grid_attributes(self, meshed_space, attribute): assert not hasattr(meshed_space, attribute) def test_structured_field_rejects_a_mesh(self, meshed_space): - from stride.problem.data import ScalarField - from stride.problem.domain import Grid # Fails on the first grid attribute it reaches for, rather than silently constructing. with pytest.raises(AttributeError): - ScalarField(name='sigma', grid=Grid(meshed_space, None, None)) + ScalarField(name='alpha', space=meshed_space) def test_meshed_field_rejects_a_structured_space(self, structured_space): - from stride.problem.data import MeshedField - from stride.problem.domain import Grid - # The reverse direction needs an explicit check: MeshedData sizes itself behind an # isinstance test, which would otherwise skip and leave a field with no shape that # still allocates and fills without complaint. with pytest.raises(ValueError, match='MeshedSpace'): - MeshedField(name='sigma', grid=Grid(structured_space, None, None)) + MeshedField(name='alpha', space=structured_space) def test_no_space_is_still_allowed(self): - from stride.problem.data import MeshedField - # This is what an instance about to be loaded from file looks like. - field = MeshedField(name='sigma') + field = MeshedField(name='alpha') assert field.space is None class TestMeshedSpaceBounds: - """Port of the mesh-covers-the-grid assertions in ae_modelling attach_mesh.""" + """The mesh-covers-the-grid assertions made when attaching a mesh to a grid.""" def test_contains_its_own_bounds(self, meshed_space): assert meshed_space.contains_box(meshed_space.origin, meshed_space.limit) @@ -227,8 +209,6 @@ def _dolfinx_box(self): ) def test_nodes_come_from_mesh_geometry(self): - from stride.problem.domain import MeshedSpace - mesh = self._dolfinx_box() space = MeshedSpace.from_dolfinx(mesh) @@ -237,8 +217,6 @@ def test_nodes_come_from_mesh_geometry(self): assert space.num_nodes == mesh.geometry.x.shape[0] def test_bounds_match_the_dolfinx_mesh(self): - from stride.problem.domain import MeshedSpace - mesh = self._dolfinx_box() space = MeshedSpace.from_dolfinx(mesh) @@ -246,8 +224,6 @@ def test_bounds_match_the_dolfinx_mesh(self): np.testing.assert_allclose(space.limit, mesh.geometry.x.max(axis=0)) def test_cells_come_from_topology(self): - from stride.problem.domain import MeshedSpace - mesh = self._dolfinx_box() space = MeshedSpace.from_dolfinx(mesh) From 079bf41b647d127bced826f038663720478746cd Mon Sep 17 00:00:00 2001 From: Andrei Danila Date: Fri, 14 Aug 2026 14:50:11 +0100 Subject: [PATCH 08/15] gh actions added pytest --- .github/workflows/pytest.yml | 62 ++++++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) create mode 100644 .github/workflows/pytest.yml diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml new file mode 100644 index 0000000..a15df71 --- /dev/null +++ b/.github/workflows/pytest.yml @@ -0,0 +1,62 @@ +name: Pytest + +on: + # Trigger the workflow on push or pull request, + # but only for the master branch + push: + branches: + - master + pull_request: + branches: + - master + +jobs: + examples: + name: Pytest + runs-on: ubuntu-latest + + env: + DEVITO_COMPILER: gcc + DEVITO_LANGUAGE: openmp + PYTHON_VERSION: 3.11.12 + + strategy: + # Prevent all build to stop if a single one fails + fail-fast: false + + steps: + - name: Checkout stride + uses: actions/checkout@v3 + with: + path: stride + + - name: Checkout devito + uses: actions/checkout@v3 + with: + repository: devitocodes/devito + path: devito + + - name: Setup conda + uses: mamba-org/setup-micromamba@v2 + with: + environment-file: stride/environment.yml + init-shell: bash + cache-environment: true + post-cleanup: 'all' + + - name: Install dependencies + shell: bash -l {0} + run: | + cd stride + pip install -e . + + - name: Install devito + shell: bash -l {0} + run: | + cd devito + pip install -e . + + - name: Pytest + shell: bash -l {0} + run: | + python -m pytest stride/tests/ -v From 330c0b75328caeea3504c78978146680652da4e6 Mon Sep 17 00:00:00 2001 From: Andrei Danila Date: Fri, 14 Aug 2026 15:25:15 +0100 Subject: [PATCH 09/15] fix actions --- .github/workflows/pytest.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml index a15df71..069f90c 100644 --- a/.github/workflows/pytest.yml +++ b/.github/workflows/pytest.yml @@ -59,4 +59,4 @@ jobs: - name: Pytest shell: bash -l {0} run: | - python -m pytest stride/tests/ -v + python -m pytest tests/ -v From 034b246fd77ff3056023c48eea4a6f9f39b027ad Mon Sep 17 00:00:00 2001 From: Andrei Danila Date: Tue, 25 Aug 2026 17:14:27 +0100 Subject: [PATCH 10/15] added to_dolfinx function in MeshedSpace in order to be able to preserve mesh --- stride/problem/base.py | 20 ++ stride/problem/domain.py | 138 +++++++++++- stride/tests/test_meshed_serialisation.py | 80 +++++++ stride/tests/test_meshed_space.py | 249 ++++++++++++++++++++++ 4 files changed, 485 insertions(+), 2 deletions(-) diff --git a/stride/problem/base.py b/stride/problem/base.py index 521b0f9..d808d88 100644 --- a/stride/problem/base.py +++ b/stride/problem/base.py @@ -258,10 +258,23 @@ def load(self, *args, **kwargs): absorbing=space_description.absorbing) elif 'nodes' in space_description: + # the cell type is stored as a string, so HDF5 hands it back as bytes, and + # the degree comes back as a numpy integer, which is not an int as far as + # isinstance is concerned. Both have to be normalised before construction + cell_type = self._materialise(space_description.get('cell_type', None)) + + if isinstance(cell_type, bytes): + cell_type = cell_type.decode() + + geometry_degree = self._materialise( + space_description.get('geometry_degree', 1)) + space = MeshedSpace( nodes=self._materialise(space_description.nodes), cells=self._materialise(space_description.get('cells', None)), cell_tags=self._materialise(space_description.get('cell_tags', None)), + cell_type=cell_type, + geometry_degree=int(geometry_degree), ) else: @@ -350,6 +363,13 @@ def grid_description(self): if space.cell_tags is not None: space_description['cell_tags'] = space.cell_tags + # the discretisation, without which the mesh cannot be rebuilt: a node and cell + # table alone does not say what shape a cell is or to what degree it is mapped + if space.cell_type is not None: + space_description['cell_type'] = space.cell_type + + space_description['geometry_degree'] = space.geometry_degree + grid_description['space'] = space_description else: diff --git a/stride/problem/domain.py b/stride/problem/domain.py index acd4ebf..c4ba2bd 100644 --- a/stride/problem/domain.py +++ b/stride/problem/domain.py @@ -6,6 +6,24 @@ __all__ = ['Space', 'MeshedSpace', 'Time', 'SlowTime', 'Grid'] +# Topological dimension of every supported cell type. Simplices only for now: adding +# quadrilateral or hexahedron here is most of what supporting them takes. +CELL_TOPOLOGICAL_DIM = { + 'triangle': 2, + 'tetrahedron': 3, +} + +# Cell type implied by (nodes per cell, geometry degree). Within simplices this is unique, so a +# cell type that is not given can be inferred rather than demanded. It stops being unique as soon +# as non-simplices are supported -- six nodes at degree one is a prism, not a triangle -- which is +# why anything not in here raises instead of guessing. +CELL_TYPE_BY_NODES = { + (3, 1): 'triangle', + (6, 2): 'triangle', + (4, 1): 'tetrahedron', + (10, 2): 'tetrahedron', +} + class Space: """ @@ -292,10 +310,22 @@ class MeshedSpace: facet_tags : optional Boundary tags, stored as given and not interpreted. These are not serialised, because a facet tag is meaningless without the facet connectivity, which is not stored either. + cell_type : str, optional + Cell type of the mesh, either ``triangle`` or ``tetrahedron``. Only simplices are + supported. The spelling is the one DOLFINx and basix use, so that it can be handed to + either without translation. This fixes the *topological* dimension, which need not equal + ``dim``: a triangle mesh with three-dimensional nodes is a surface embedded in 3D, so a + cell may live in a space of higher dimension than its own but never a lower one. If not + given, it is inferred from the nodes per cell and the geometry degree, and a combination + that is not unambiguous raises rather than being guessed at. + geometry_degree : int, optional + Degree of the coordinate element: 1 for straight-sided cells, 2 for curved ones. + Defaults to 1. This describes the mesh geometry only, and is unrelated to the degree of + any function space later defined over it. """ - def __init__(self, nodes=None, cells=None, cell_tags=None, facet_tags=None): + def __init__(self, nodes=None, cells=None, cell_tags=None, facet_tags=None, cell_type=None, geometry_degree=1): nodes = np.asarray(nodes, dtype=np.float64) if nodes.ndim != 2: @@ -316,6 +346,14 @@ def __init__(self, nodes=None, cells=None, cell_tags=None, facet_tags=None): if cells.size and (cells.min() < 0 or cells.max() >= nodes.shape[0]): raise ValueError('Cells reference node indices outside [0, %d)' % nodes.shape[0]) + if cell_type is None: + cell_type = CELL_TYPE_BY_NODES.get((cells.shape[1], geometry_degree)) + + if cell_type is None: + raise ValueError('Cannot infer the cell type from %d nodes per cell at ' + 'geometry degree %s. Pass cell_type explicitly.' + % (cells.shape[1], geometry_degree)) + if cell_tags is not None: cell_tags = np.asarray(cell_tags) @@ -326,11 +364,31 @@ def __init__(self, nodes=None, cells=None, cell_tags=None, facet_tags=None): raise ValueError('Cell tags must have one entry per cell, expected %d ' 'but got shape %s' % (cells.shape[0], (cell_tags.shape,))) + if not isinstance(geometry_degree, int): + raise TypeError('geometry_degree must be an int') + if geometry_degree < 1: + raise ValueError('geometry degree must be a positive int') + + if cell_type is not None: + if not isinstance(cell_type, str): + raise TypeError('cell_type must be str') + elif cell_type not in CELL_TOPOLOGICAL_DIM: + raise ValueError('Only simplex cells are supported (%s), got %r' + % (', '.join(sorted(CELL_TOPOLOGICAL_DIM)), cell_type)) + + # a cell can be embedded in a space of higher dimension than its own, but not lower: + # a triangle mesh may be a surface in 3D, a tetrahedron cannot live in 2D + if dim < CELL_TOPOLOGICAL_DIM[cell_type]: + raise ValueError('A %s is %d-dimensional and cannot be embedded in %d dimensions' + % (cell_type, CELL_TOPOLOGICAL_DIM[cell_type], dim)) + self.dim = dim self.nodes = nodes self.cells = cells self.cell_tags = cell_tags self.facet_tags = facet_tags + self.cell_type = cell_type + self.geometry_degree = geometry_degree @property def num_nodes(self): @@ -518,7 +576,83 @@ def from_dolfinx(cls, mesh, cell_tags=None, facet_tags=None): owned = indices < num_cells tags[indices[owned]] = values[owned] - return cls(nodes=nodes, cells=cells, cell_tags=tags, facet_tags=facet_tags) + return cls(nodes=nodes, cells=cells, cell_tags=tags, facet_tags=facet_tags, + cell_type=mesh.topology.cell_name(), + geometry_degree=mesh.geometry.cmap.degree) + + def to_dolfinx(self, comm=None): + """ + Create an in-memory DOLFINx mesh from this MeshedSpace. + + This is the inverse of :meth:`from_dolfinx`, and is only correct in serial. The node and + cell tables stored here describe a whole mesh, so building on a communicator of more + than one rank would partition it into something this space does not describe. Unlike + ``from_dolfinx``, which warns, that case raises. + + DOLFINx is an optional dependency of stride, so it is imported here rather than at + module level. + + Parameters + ---------- + comm : MPI.Intracomm, optional + Communicator on which to build the mesh, defaults to ``MPI.COMM_SELF``. + + Returns + ------- + dolfinx.mesh.Mesh + Newly created mesh. + dolfinx.mesh.MeshTags or None + Cell tags on the new cell ordering, or None if this space carries no tags. Cells + marked -1, which is what ``from_dolfinx`` fills in for untagged cells, are left out + rather than handed back as a label of -1. + + """ + try: + import ufl + import basix.ufl + import dolfinx + from mpi4py import MPI + + except ImportError: + raise ImportError('to_dolfinx needs dolfinx, basix, ufl and mpi4py, which are ' + 'optional dependencies of stride. Install them into the ' + 'environment, for instance with the fenics-dolfinx conda ' + 'package.') from None + + if comm is None: + comm = MPI.COMM_SELF + + if self.cells is None: + raise ValueError('Cannot build a mesh without cells, a node table is not a mesh') + + if comm.size > 1: + raise ValueError('to_dolfinx is serial only, got a communicator of size %d. Build ' + 'on MPI.COMM_SELF and partition afterwards if needed.' % comm.size) + + topology_dim = CELL_TOPOLOGICAL_DIM[self.cell_type] + + element = ufl.Mesh(basix.ufl.element('Lagrange', self.cell_type, self.geometry_degree, + shape=(self.dim,))) + + mesh = dolfinx.mesh.create_mesh(comm, self.cells, element, self.nodes) + + num_cells = mesh.topology.index_map(topology_dim).size_local + if num_cells != self.num_cells: + raise RuntimeError('The rebuilt mesh has %d cells but this space has %d, so the ' + 'cell ordering cannot be recovered' % (num_cells, self.num_cells)) + + cell_tags = None + if self.cell_tags is not None: + original_index = np.asarray(mesh.topology.original_cell_index) + values = np.asarray(self.cell_tags)[original_index] + + tagged = values != -1 + + cell_tags = dolfinx.mesh.meshtags(mesh, topology_dim, + np.arange(num_cells, dtype=np.int32)[tagged], + values[tagged]) + + return mesh, cell_tags class Time: diff --git a/stride/tests/test_meshed_serialisation.py b/stride/tests/test_meshed_serialisation.py index 0f67e9f..8f4c0ad 100644 --- a/stride/tests/test_meshed_serialisation.py +++ b/stride/tests/test_meshed_serialisation.py @@ -16,6 +16,9 @@ from the top-level description - an unrecognised space payload raises a clear, typed error - a MeshedField round-trips its data and its mesh through HDF5 +- ``cell_type``/``geometry_degree`` survive the round trip, so that a space read + back from disk can still be handed to ``to_dolfinx``; a file written before + those fields existed still loads """ import numpy as np @@ -245,3 +248,80 @@ def test_time_dependent_round_trip(self, meshed_space, project): assert loaded.time.num == 5 assert tuple(loaded.shape) == (5, 27) np.testing.assert_allclose(loaded.data, 2.) + + +class TestDiscretisationSurvivesTheRoundTrip: + """ + cell_type and geometry_degree are what make a stored space rebuildable. Without + them a node and cell table does not say what shape a cell is, so to_dolfinx has + nothing to hand basix. + """ + + def test_description_carries_the_discretisation(self, meshed_space): + field = MeshedField(name='alpha', space=meshed_space) + space_description = field.grid_description()['space'] + + assert space_description['cell_type'] == 'tetrahedron' + assert space_description['geometry_degree'] == 1 + + def test_round_trip_keeps_the_cell_type(self, meshed_space, project): + field = MeshedField(name='alpha', space=meshed_space, dtype=np.float64) + field.fill(1.) + field.dump(**project) + + loaded = MeshedField(name='alpha') + loaded.load(**project) + + assert isinstance(loaded.space, MeshedSpace) + assert loaded.space.cell_type == 'tetrahedron' + assert loaded.space.geometry_degree == 1 + + def test_cell_type_survives_as_str_not_bytes(self, meshed_space, project): + """HDF5 hands strings back as bytes, which would break the isinstance check.""" + field = MeshedField(name='alpha', space=meshed_space, dtype=np.float64) + field.dump(**project) + + loaded = MeshedField(name='alpha') + loaded.load(**project) + + assert isinstance(loaded.space.cell_type, str) + assert isinstance(loaded.space.geometry_degree, int) + + def test_a_loaded_space_can_still_be_rebuilt(self, meshed_space, project): + """The point of storing them at all.""" + pytest.importorskip('dolfinx') + + field = MeshedField(name='alpha', space=meshed_space, dtype=np.float64) + field.dump(**project) + + loaded = MeshedField(name='alpha') + loaded.load(**project) + + mesh, _ = loaded.space.to_dolfinx() + + assert mesh.topology.dim == 3 + assert mesh.topology.index_map(3).size_local == meshed_space.num_cells + + def test_a_file_without_the_new_fields_still_loads(self, meshed_space, project): + """ + Backward compatibility: the fields are read with defaults, so a space written + before they existed reconstructs, falling back on inference for the cell type. + """ + field = MeshedField(name='alpha', space=meshed_space, dtype=np.float64) + description = field.grid_description() + + del description['space']['cell_type'] + del description['space']['geometry_degree'] + + assert 'cell_type' not in description['space'] + + # what the load branch does with what it finds, minus the HDF5 layer + space = MeshedSpace( + nodes=description['space']['nodes'], + cells=description['space']['cells'], + cell_type=description['space'].get('cell_type', None), + geometry_degree=int(description['space'].get('geometry_degree', 1)), + ) + + assert space.cell_type == 'tetrahedron' + assert space.geometry_degree == 1 diff --git a/stride/tests/test_meshed_space.py b/stride/tests/test_meshed_space.py index 557870b..23f180a 100644 --- a/stride/tests/test_meshed_space.py +++ b/stride/tests/test_meshed_space.py @@ -18,6 +18,12 @@ - ``contains_box(lower, upper)`` is the mesh-covers-the-grid check performed when a mesh is attached to a problem grid - ``resample`` raises, because a mesh has no spacing to resample onto +- ``cell_type`` / ``geometry_degree`` describe the discretisation, are inferred + from the nodes per cell when not given, and fix the topological dimension, + which may be lower than ``dim`` for a surface mesh +- ``to_dolfinx`` is the inverse of ``from_dolfinx``: it rebuilds a mesh, maps + the cell tags onto whatever ordering ``create_mesh`` chose, and re-sparsifies + the ``-1`` fill that ``from_dolfinx`` introduced - ``from_dolfinx`` adapts an in-memory DOLFINx mesh (skipped without DOLFINx) """ @@ -233,3 +239,246 @@ def test_cells_come_from_topology(self): assert space.num_cells == num_cells assert space.cells.shape[1] == 4 assert space.cells.max() < space.num_nodes + + +class TestMeshedSpaceCellType: + """cell_type / geometry_degree, and the inference that fills them in.""" + + def test_triangle_inferred_from_three_nodes(self, meshed_space_2d): + assert meshed_space_2d.cell_type == 'triangle' + assert meshed_space_2d.geometry_degree == 1 + + def test_tetrahedron_inferred_from_four_nodes(self, meshed_space): + assert meshed_space.cell_type == 'tetrahedron' + assert meshed_space.geometry_degree == 1 + + def test_explicit_cell_type_is_kept(self, tri_mesh): + nodes, cells = tri_mesh + space = MeshedSpace(nodes=nodes, cells=cells, cell_type='triangle') + + assert space.cell_type == 'triangle' + + def test_cell_type_is_none_without_cells(self, tri_mesh): + nodes, _ = tri_mesh + + assert MeshedSpace(nodes=nodes).cell_type is None + + def test_unsupported_cell_type_names_the_limit(self, tri_mesh): + nodes, cells = tri_mesh + + with pytest.raises(ValueError, match='simplex'): + MeshedSpace(nodes=nodes, cells=cells, cell_type='hexahedron') + + def test_uninferrable_cell_asks_for_an_explicit_type(self, tri_mesh): + nodes, _ = tri_mesh + # five nodes per cell is not a simplex at any supported degree + cells = np.zeros((2, 5), dtype=np.int32) + + with pytest.raises(ValueError, match='infer'): + MeshedSpace(nodes=nodes, cells=cells) + + def test_a_triangle_may_be_a_surface_in_three_dimensions(self): + nodes = np.array([[0., 0., 0.], [1e-3, 0., 0.], + [0., 1e-3, 1e-3], [1e-3, 1e-3, 1e-3]]) + cells = np.array([[0, 1, 2], [1, 3, 2]], dtype=np.int32) + space = MeshedSpace(nodes=nodes, cells=cells, cell_type='triangle') + + assert space.cell_type == 'triangle' + assert space.dim == 3 + + def test_a_tetrahedron_cannot_live_in_two_dimensions(self, tri_mesh): + nodes, _ = tri_mesh + cells = np.zeros((2, 4), dtype=np.int32) + + with pytest.raises(ValueError, match='cannot be embedded'): + MeshedSpace(nodes=nodes, cells=cells, cell_type='tetrahedron') + + def test_degree_must_be_a_positive_int(self, tri_mesh): + nodes, cells = tri_mesh + + with pytest.raises(ValueError): + MeshedSpace(nodes=nodes, cells=cells, cell_type='triangle', geometry_degree=0) + + with pytest.raises(TypeError): + MeshedSpace(nodes=nodes, cells=cells, cell_type='triangle', geometry_degree=1.0) + + +@pytest.mark.skipif(not HAS_DOLFINX, reason='DOLFINx not available') +class TestMeshedSpaceToDolfinx: + """ + to_dolfinx has to reproduce a mesh that a solver can use, which is a stronger + claim than the individual arrays matching. The checks below are chosen so that + none of them can be satisfied by a plausible-but-wrong reconstruction. + """ + + def _tagged_square(self, n=4): + """A unit-square mesh tagged 1 left of x=0.5 and 2 right of it.""" + import dolfinx + from mpi4py import MPI + + mesh = dolfinx.mesh.create_unit_square(MPI.COMM_SELF, n, n) + tdim = mesh.topology.dim + cells = np.arange(mesh.topology.index_map(tdim).size_local, dtype=np.int32) + midpoints = dolfinx.mesh.compute_midpoints(mesh, tdim, cells) + tags = self._predicate(midpoints).astype(np.int32) + + return mesh, dolfinx.mesh.meshtags(mesh, tdim, cells, tags) + + @staticmethod + def _predicate(midpoints): + return np.where(midpoints[:, 0] < 0.5, 1, 2) + + @staticmethod + def _total_volume(mesh): + """Sum of triangle areas, straight from the geometry.""" + coordinates = mesh.geometry.x + dofmap = mesh.geometry.dofmap + num_cells = mesh.topology.index_map(mesh.topology.dim).size_local + + total = 0. + for cell in range(num_cells): + points = coordinates[dofmap[cell]][:, :2] + first, second = points[1] - points[0], points[2] - points[0] + total += 0.5*abs(first[0]*second[1] - first[1]*second[0]) + + return total + + def test_rebuilds_a_usable_mesh(self): + mesh, tags = self._tagged_square() + space = MeshedSpace.from_dolfinx(mesh, cell_tags=tags) + + rebuilt, _ = space.to_dolfinx() + + assert rebuilt.topology.dim == 2 + assert rebuilt.geometry.dim == 2 + assert rebuilt.topology.index_map(2).size_local == space.num_cells + + def test_node_coordinates_survive_as_a_set(self): + mesh, tags = self._tagged_square() + space = MeshedSpace.from_dolfinx(mesh, cell_tags=tags) + + rebuilt, _ = space.to_dolfinx() + + # ordering need not survive, the set of positions must + before = set(map(tuple, np.round(space.nodes, 12))) + after = set(map(tuple, np.round(rebuilt.geometry.x[:, :2], 12))) + + assert before == after + assert len(before) == space.num_nodes + + def test_total_volume_is_preserved(self): + """Catches scrambled connectivity, which per-node coordinate checks do not.""" + mesh, tags = self._tagged_square() + space = MeshedSpace.from_dolfinx(mesh, cell_tags=tags) + + rebuilt, _ = space.to_dolfinx() + + np.testing.assert_allclose(self._total_volume(rebuilt), self._total_volume(mesh)) + np.testing.assert_allclose(self._total_volume(rebuilt), 1.) + + def test_tags_land_on_the_same_cells(self): + """ + The load-bearing test. create_mesh may reorder cells, so the tags have to be + permuted. Rather than trusting the permutation, this re-derives what the tag + of every rebuilt cell ought to be from its own midpoint, so it holds + regardless of which way round the permutation was applied. + """ + import dolfinx + + mesh, tags = self._tagged_square() + space = MeshedSpace.from_dolfinx(mesh, cell_tags=tags) + + rebuilt, rebuilt_tags = space.to_dolfinx() + + num_cells = rebuilt.topology.index_map(2).size_local + midpoints = dolfinx.mesh.compute_midpoints( + rebuilt, 2, np.arange(num_cells, dtype=np.int32)) + + expected = self._predicate(midpoints) + actual = np.full(num_cells, -1, dtype=np.int32) + actual[rebuilt_tags.indices] = rebuilt_tags.values + + np.testing.assert_array_equal(actual, expected) + + def test_untagged_cells_do_not_come_back_as_minus_one(self): + """ + from_dolfinx densifies sparse tags with -1. Handing that back would turn + 'untagged' into a material label of -1, which a lookup table would then + either fail on or silently honour. + """ + import dolfinx + from mpi4py import MPI + + mesh = dolfinx.mesh.create_unit_square(MPI.COMM_SELF, 4, 4) + tagged = np.array([0, 1, 2], dtype=np.int32) + tags = dolfinx.mesh.meshtags(mesh, 2, tagged, np.full(3, 7, dtype=np.int32)) + + space = MeshedSpace.from_dolfinx(mesh, cell_tags=tags) + assert (space.cell_tags == -1).any(), 'fixture should leave most cells untagged' + + _, rebuilt_tags = space.to_dolfinx() + + assert len(rebuilt_tags.indices) == 3 + assert -1 not in rebuilt_tags.values + np.testing.assert_array_equal(np.unique(rebuilt_tags.values), [7]) + + def test_no_tags_gives_no_tags(self): + mesh, _ = self._tagged_square() + space = MeshedSpace.from_dolfinx(mesh) + + _, rebuilt_tags = space.to_dolfinx() + + assert rebuilt_tags is None + + def test_round_trips_twice(self): + mesh, tags = self._tagged_square() + space = MeshedSpace.from_dolfinx(mesh, cell_tags=tags) + + rebuilt, rebuilt_tags = space.to_dolfinx() + again = MeshedSpace.from_dolfinx(rebuilt, cell_tags=rebuilt_tags) + + assert again.cell_type == space.cell_type + assert again.geometry_degree == space.geometry_degree + assert again.num_nodes == space.num_nodes + assert again.num_cells == space.num_cells + np.testing.assert_array_equal(np.sort(again.cell_tags), np.sort(space.cell_tags)) + + def test_a_surface_mesh_keeps_its_dimensions(self): + """tdim 2 with gdim 3 has to survive, which is why cell_type is stored.""" + nodes = np.array([[0., 0., 0.], [1., 0., 0.], [0., 1., 1.], [1., 1., 1.]]) + cells = np.array([[0, 1, 2], [1, 3, 2]], dtype=np.int32) + space = MeshedSpace(nodes=nodes, cells=cells, cell_type='triangle') + + rebuilt, _ = space.to_dolfinx() + + assert rebuilt.topology.dim == 2 + assert rebuilt.geometry.dim == 3 + + def test_without_cells_it_refuses(self, tri_mesh): + nodes, _ = tri_mesh + + with pytest.raises(ValueError, match='not a mesh'): + MeshedSpace(nodes=nodes).to_dolfinx() + + def test_a_parallel_communicator_is_refused(self): + from mpi4py import MPI + + mesh, _ = self._tagged_square() + space = MeshedSpace.from_dolfinx(mesh) + + class _Parallel: + size = 4 + + with pytest.raises(ValueError, match='serial only'): + space.to_dolfinx(comm=_Parallel()) + + def test_an_explicit_communicator_is_honoured(self): + """comm=None means 'pick a default', not 'ignore what you were given'.""" + from mpi4py import MPI + + mesh, _ = self._tagged_square() + space = MeshedSpace.from_dolfinx(mesh) + + rebuilt, _ = space.to_dolfinx(comm=MPI.COMM_SELF) + + assert rebuilt.comm.size == 1 From 5c7cacfb2bb523a255bed1408838db023b32394c Mon Sep 17 00:00:00 2001 From: Andrei Danila Date: Wed, 26 Aug 2026 10:04:28 +0100 Subject: [PATCH 11/15] fixed pytest action --- .github/workflows/pytest.yml | 1 + stride/problem/domain.py | 7 +------ 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml index 069f90c..e141298 100644 --- a/.github/workflows/pytest.yml +++ b/.github/workflows/pytest.yml @@ -59,4 +59,5 @@ jobs: - name: Pytest shell: bash -l {0} run: | + cd .. python -m pytest tests/ -v diff --git a/stride/problem/domain.py b/stride/problem/domain.py index c4ba2bd..1bda0bb 100644 --- a/stride/problem/domain.py +++ b/stride/problem/domain.py @@ -6,17 +6,12 @@ __all__ = ['Space', 'MeshedSpace', 'Time', 'SlowTime', 'Grid'] -# Topological dimension of every supported cell type. Simplices only for now: adding -# quadrilateral or hexahedron here is most of what supporting them takes. + CELL_TOPOLOGICAL_DIM = { 'triangle': 2, 'tetrahedron': 3, } -# Cell type implied by (nodes per cell, geometry degree). Within simplices this is unique, so a -# cell type that is not given can be inferred rather than demanded. It stops being unique as soon -# as non-simplices are supported -- six nodes at degree one is a prism, not a triangle -- which is -# why anything not in here raises instead of guessing. CELL_TYPE_BY_NODES = { (3, 1): 'triangle', (6, 2): 'triangle', From 4861703c519ecb366873a1443816faaa1de839cd Mon Sep 17 00:00:00 2001 From: Andrei Danila Date: Wed, 26 Aug 2026 10:22:47 +0100 Subject: [PATCH 12/15] cleaned imports and flake8 --- stride/tests/test_meshed_medium.py | 1 + stride/tests/test_meshed_serialisation.py | 1 + stride/tests/test_meshed_space.py | 11 +---------- 3 files changed, 3 insertions(+), 10 deletions(-) diff --git a/stride/tests/test_meshed_medium.py b/stride/tests/test_meshed_medium.py index 97cb3a0..6368085 100644 --- a/stride/tests/test_meshed_medium.py +++ b/stride/tests/test_meshed_medium.py @@ -9,6 +9,7 @@ from stride.problem.data import MeshedField from stride.problem.medium import Medium + def value_lut(material_properties, prop='alpha'): """ One property of the material table, as an array indexed by label. diff --git a/stride/tests/test_meshed_serialisation.py b/stride/tests/test_meshed_serialisation.py index 8f4c0ad..aa44a7e 100644 --- a/stride/tests/test_meshed_serialisation.py +++ b/stride/tests/test_meshed_serialisation.py @@ -27,6 +27,7 @@ from stride.problem.domain import Grid, Time, Space, MeshedSpace from stride.problem.data import ScalarField, MeshedField + @pytest.fixture def project(tmp_path): """Path/project_name pair for the HDF5 helpers.""" diff --git a/stride/tests/test_meshed_space.py b/stride/tests/test_meshed_space.py index 23f180a..7a60714 100644 --- a/stride/tests/test_meshed_space.py +++ b/stride/tests/test_meshed_space.py @@ -35,10 +35,9 @@ from stride.problem.domain import MeshedSpace from stride.problem.data import MeshedField, ScalarField - - try: import dolfinx # noqa: F401 + from mpi4py import MPI HAS_DOLFINX = True except ImportError: HAS_DOLFINX = False @@ -206,7 +205,6 @@ def test_resample_is_not_supported(self, meshed_space): class TestMeshedSpaceFromDolfinx: def _dolfinx_box(self): - from mpi4py import MPI return dolfinx.mesh.create_box( MPI.COMM_WORLD, @@ -313,8 +311,6 @@ class TestMeshedSpaceToDolfinx: def _tagged_square(self, n=4): """A unit-square mesh tagged 1 left of x=0.5 and 2 right of it.""" - import dolfinx - from mpi4py import MPI mesh = dolfinx.mesh.create_unit_square(MPI.COMM_SELF, n, n) tdim = mesh.topology.dim @@ -383,7 +379,6 @@ def test_tags_land_on_the_same_cells(self): of every rebuilt cell ought to be from its own midpoint, so it holds regardless of which way round the permutation was applied. """ - import dolfinx mesh, tags = self._tagged_square() space = MeshedSpace.from_dolfinx(mesh, cell_tags=tags) @@ -406,8 +401,6 @@ def test_untagged_cells_do_not_come_back_as_minus_one(self): 'untagged' into a material label of -1, which a lookup table would then either fail on or silently honour. """ - import dolfinx - from mpi4py import MPI mesh = dolfinx.mesh.create_unit_square(MPI.COMM_SELF, 4, 4) tagged = np.array([0, 1, 2], dtype=np.int32) @@ -461,7 +454,6 @@ def test_without_cells_it_refuses(self, tri_mesh): MeshedSpace(nodes=nodes).to_dolfinx() def test_a_parallel_communicator_is_refused(self): - from mpi4py import MPI mesh, _ = self._tagged_square() space = MeshedSpace.from_dolfinx(mesh) @@ -474,7 +466,6 @@ class _Parallel: def test_an_explicit_communicator_is_honoured(self): """comm=None means 'pick a default', not 'ignore what you were given'.""" - from mpi4py import MPI mesh, _ = self._tagged_square() space = MeshedSpace.from_dolfinx(mesh) From a92bdab7b88716601eaafd8568dd0e2d3f1b80a2 Mon Sep 17 00:00:00 2001 From: Andrei Danila Date: Wed, 26 Aug 2026 10:57:19 +0100 Subject: [PATCH 13/15] fix pytest path in CI --- .github/workflows/pytest.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml index e141298..97229da 100644 --- a/.github/workflows/pytest.yml +++ b/.github/workflows/pytest.yml @@ -59,5 +59,5 @@ jobs: - name: Pytest shell: bash -l {0} run: | - cd .. - python -m pytest tests/ -v + cd stride + python -m pytest stride/tests/ -v From 9431ceeba3f39bb36666fd4b7741e763cf2b38cc Mon Sep 17 00:00:00 2001 From: Andrei Danila Date: Wed, 2 Sep 2026 13:14:41 +0100 Subject: [PATCH 14/15] minor change for in create_mesh signature for dolfinx version compatibility and env.yaml files for dolfinx builds --- environment-fem-complex.yml | 36 ++++++++++++++++++++++++++++++++++++ environment-fem-real.yml | 36 ++++++++++++++++++++++++++++++++++++ stride/problem/domain.py | 2 +- 3 files changed, 73 insertions(+), 1 deletion(-) create mode 100644 environment-fem-complex.yml create mode 100644 environment-fem-real.yml diff --git a/environment-fem-complex.yml b/environment-fem-complex.yml new file mode 100644 index 0000000..a1e9196 --- /dev/null +++ b/environment-fem-complex.yml @@ -0,0 +1,36 @@ +# Finite-element extras for stride, complex-scalar PETSc build. +# +# This file is ADDITIVE: it lists only what the FEM problem types need on top of the base +# environment, and is applied over an existing environment rather than creating one: +# +# conda env create -n stride-dolfinx-complex -f environment.yml +# conda env update -n stride-dolfinx-complex -f environment-fem-complex.yml +# +# Copying the base dependency list into this file instead is what lets the two drift apart, so +# please keep it additive. +# +# DOLFINx cannot be installed with pip -- there are no wheels -- so these dependencies are +# conda-only and deliberately not in requirements-optional.txt. Everything in stride works +# without them; only the FEM problem types need them. +# +# PETSc is built for either real or complex scalars and the two cannot coexist. Use this file +# for the complex build, environment-fem-real.yml for the real one. A capacitive medium +# needs the complex build, because a real one discards the imaginary part of the admittivity +# with only a warning. +# +# The MPI implementation is pinned here rather than in the base environment, which needs no MPI +# at all. Mixing implementations between the base and this overlay causes hangs at exit inside +# MPI finalisation, so if the base environment ever grows an MPI dependency, the two pins must +# agree. +name: stride-dolfinx-complex +channels: + - conda-forge +dependencies: + - fenics-dolfinx==0.9.0 + - petsc=*=*complex* + - mpi=*=*openmpi* + - openmpi + - petsc4py + - mpi4py + - gmsh # mesh generation, plus its Python bindings + - python-gmsh diff --git a/environment-fem-real.yml b/environment-fem-real.yml new file mode 100644 index 0000000..daa606f --- /dev/null +++ b/environment-fem-real.yml @@ -0,0 +1,36 @@ +# Finite-element extras for stride, real-scalar PETSc build. +# +# This file is ADDITIVE: it lists only what the FEM problem types need on top of the base +# environment, and is applied over an existing environment rather than creating one: +# +# conda env create -n stride-dolfinx-real -f environment.yml +# conda env update -n stride-dolfinx-real -f environment-fem-real.yml +# +# Copying the base dependency list into this file instead is what lets the two drift apart, so +# please keep it additive. +# +# DOLFINx cannot be installed with pip -- there are no wheels -- so these dependencies are +# conda-only and deliberately not in requirements-optional.txt. Everything in stride works +# without them; only the FEM problem types need them. +# +# PETSc is built for either real or complex scalars and the two cannot coexist. Use this file +# for the real build, environment-fem-complex.yml for the complex one. A capacitive medium +# needs the complex build, because a real one discards the imaginary part of the admittivity +# with only a warning. +# +# The MPI implementation is pinned here rather than in the base environment, which needs no MPI +# at all. Mixing implementations between the base and this overlay causes hangs at exit inside +# MPI finalisation, so if the base environment ever grows an MPI dependency, the two pins must +# agree. +name: stride-dolfinx-real +channels: + - conda-forge +dependencies: + - fenics-dolfinx==0.9.0 + - petsc=*=*real* + - mpi=*=*openmpi* + - openmpi + - petsc4py + - mpi4py + - gmsh # mesh generation, plus its Python bindings + - python-gmsh diff --git a/stride/problem/domain.py b/stride/problem/domain.py index 1bda0bb..cb7057f 100644 --- a/stride/problem/domain.py +++ b/stride/problem/domain.py @@ -629,7 +629,7 @@ def to_dolfinx(self, comm=None): element = ufl.Mesh(basix.ufl.element('Lagrange', self.cell_type, self.geometry_degree, shape=(self.dim,))) - mesh = dolfinx.mesh.create_mesh(comm, self.cells, element, self.nodes) + mesh = dolfinx.mesh.create_mesh(comm, cells=self.cells, x=self.nodes, e=element) num_cells = mesh.topology.index_map(topology_dim).size_local if num_cells != self.num_cells: From c863ac4465fb5c90d511488fc821fc917cc306a2 Mon Sep 17 00:00:00 2001 From: Andrei Danila Date: Mon, 7 Sep 2026 12:14:33 +0100 Subject: [PATCH 15/15] cicd included dolfinx --- .github/workflows/pytest.yml | 66 ++++++++++++++++++++++++++++++++++++ environment-fem-complex.yml | 4 +-- environment-fem-real.yml | 4 +-- 3 files changed, 68 insertions(+), 6 deletions(-) diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml index 97229da..61f3d12 100644 --- a/.github/workflows/pytest.yml +++ b/.github/workflows/pytest.yml @@ -61,3 +61,69 @@ jobs: run: | cd stride python -m pytest stride/tests/ -v + + pytest-dolfinx: + name: Pytest DOLFINx + runs-on: ubuntu-latest + + env: + DEVITO_COMPILER: gcc + DEVITO_LANGUAGE: openmp + PYTHON_VERSION: 3.11.12 + + strategy: + # Prevent all build to stop if a single one fails + fail-fast: false + + steps: + - name: Checkout stride + uses: actions/checkout@v3 + with: + path: stride + + - name: Checkout devito + uses: actions/checkout@v3 + with: + repository: devitocodes/devito + path: devito + + - name: Setup conda + uses: mamba-org/setup-micromamba@v2 + with: + environment-file: stride/environment.yml + init-shell: bash + cache-environment: true + post-cleanup: 'all' + + # applied over the base environment rather than replacing it, the same way the file is meant + # to be used locally. The complex build is the one that covers both cases: it runs the + # resistive tests as well, whereas a real build has to skip the capacitive ones + - name: Install FEM dependencies + shell: bash -l {0} + run: | + micromamba install -y -n stride -f stride/environment-fem-complex.yml + + # a missing or real-scalar build would turn the FEM tests into skips rather than failures, + # so the suite would go green having run none of them + - name: Check the build is complex + shell: bash -l {0} + run: | + python -c "import dolfinx, numpy; assert numpy.dtype(dolfinx.default_scalar_type).kind == 'c'" + + - name: Install dependencies + shell: bash -l {0} + run: | + cd stride + pip install -e . + + - name: Install devito + shell: bash -l {0} + run: | + cd devito + pip install -e . + + - name: Pytest + shell: bash -l {0} + run: | + cd stride + python -m pytest stride/tests/ -v diff --git a/environment-fem-complex.yml b/environment-fem-complex.yml index a1e9196..838cdba 100644 --- a/environment-fem-complex.yml +++ b/environment-fem-complex.yml @@ -14,9 +14,7 @@ # without them; only the FEM problem types need them. # # PETSc is built for either real or complex scalars and the two cannot coexist. Use this file -# for the complex build, environment-fem-real.yml for the real one. A capacitive medium -# needs the complex build, because a real one discards the imaginary part of the admittivity -# with only a warning. +# for the complex build, environment-fem-real.yml for the real one. # # The MPI implementation is pinned here rather than in the base environment, which needs no MPI # at all. Mixing implementations between the base and this overlay causes hangs at exit inside diff --git a/environment-fem-real.yml b/environment-fem-real.yml index daa606f..7ea6ab5 100644 --- a/environment-fem-real.yml +++ b/environment-fem-real.yml @@ -14,9 +14,7 @@ # without them; only the FEM problem types need them. # # PETSc is built for either real or complex scalars and the two cannot coexist. Use this file -# for the real build, environment-fem-complex.yml for the complex one. A capacitive medium -# needs the complex build, because a real one discards the imaginary part of the admittivity -# with only a warning. +# for the real build, environment-fem-complex.yml for the complex one. # # The MPI implementation is pinned here rather than in the base environment, which needs no MPI # at all. Mixing implementations between the base and this overlay causes hangs at exit inside