Caputo-Wismer imaging

Yonderdrake includes a small application layer for power-law acoustic attenuation on ordinary Firedrake meshes:

from yonderdrake.applications import (
    CaputoWismerMaterial,
    SensorArray,
    record_sensor_data,
    reconstruct_initial_pressure,
)

materials = (
    CaputoWismerMaterial(tissue_mask, wave_speed, damping, alpha),
)
sensors = SensorArray.ring(
    V,
    num_sensors=64,
    radius=0.9,
    width=0.04,
)
data = record_sensor_data(
    initial_pressure,
    sensors,
    materials=materials,
    dt=dt,
    num_steps=num_steps,
    num_modes=32,
    absorbing_speed=outer_wave_speed,
)
image = reconstruct_initial_pressure(
    V,
    data,
    sensors,
    materials=materials,
    dt=dt,
    num_modes=32,
    absorbing_speed=outer_wave_speed,
    method="kaltenbacher",
    attenuation="dissipative",
)

The user controls the sensor count, placement, and averaging width. Ring sensors are equally spaced in angle:

ring = SensorArray.ring(V2, 64, radius=0.9, width=0.04)
sphere = SensorArray.sphere(V3, 128, radius=0.9, width=0.06)
custom = SensorArray(V2, locations_xy, width=0.04)
Ring, sphere, and custom sensor arrays shown in a row

The constructors place sensor centres. width controls the local averaging footprint around each centre. The small 3D panel is an orthographic view of the spherical arrangement.

The sphere uses an approximately uniform golden-angle arrangement. Custom locations are an array with one (x, y) or (x, y, z) row per sensor. sensors.locations contains the positions used by the model.

The field and sensors share one scalar continuous Lagrange space. The wave and imaging demos use CG2 by default. Forward stepping and iterative reconstruction both use every rank in a distributed mesh.

The helpers use a centred wave update, so choose dt to satisfy the usual CFL restriction. The default is the regularized adjoint reconstruction adapted from Kaltenbacher and Schlintl (2022). It differentiates Yonderdrake’s actual Birk-Song wave update, including heterogeneous attenuation and the absorbing boundary.

quick_backprojection = reconstruct_initial_pressure(
    V,
    data,
    sensors,
    materials=materials,
    dt=dt,
    absorbing_speed=outer_wave_speed,
    method="adjoint",
    attenuation="dissipative",
)

lossless_backprojection = reconstruct_initial_pressure(
    V,
    data,
    sensors,
    materials=materials,
    dt=dt,
    absorbing_speed=outer_wave_speed,
    method="adjoint",
    attenuation="none",
)

method="kaltenbacher" is recommended and selected by default. It iterates to the requested tolerance unless the user supplies an iteration cap. method="adjoint" is the cheaper one-pass backprojection. The attenuation option is independent: "dissipative" uses the supplied Caputo-Wismer loss model and is the default, while "none" deliberately tests a lossless reconstruction model.