Skip to content

[BUG] CppSimulation._write_hdf5 writes p_source_flag as boolean instead of source length #697

@elma16

Description

@elma16

CppSimulation._write_hdf5 writes p_source_flag as a boolean (0 or 1), but the C++ binary expects it to be the number of time steps in the source signal. The legacy kWaveSimulation path (used by kspaceFirstOrder2D/3D) correctly sets p_source_flag = len(source.p[0]). Below is a minimal failing example. The python backend is unaffected.

import numpy as np
import h5py
from kwave.kgrid import kWaveGrid
from kwave.kmedium import kWaveMedium
from kwave.ksource import kSource
from kwave.ksensor import kSensor
from kwave.kspaceFirstOrder import kspaceFirstOrder

# ── Setup ────────────────────────────────────────────────────────────
N = (64, 64)
dx = (1e-4, 1e-4)
kgrid = kWaveGrid(N=N, spacing=dx)
kgrid.makeTime(1500)
Nt = int(kgrid.Nt)

medium = kWaveMedium(sound_speed=1500)
sensor_mask = np.zeros(N)
sensor_mask[:, 0] = 1
Ns = int(np.count_nonzero(sensor_mask))

# ── Create a time-varying pressure source ────────────────────────────
source = kSource()
source.p = np.random.randn(Ns, Nt).astype(np.float32)
source.p_mask = sensor_mask
source.p_mode = "dirichlet"

sensor = kSensor(mask=np.ones(N))

# ── Python backend works ─────────────────────────────────────────────
print("=== Time-varying source with python backend ===")
result_py = kspaceFirstOrder(
    kgrid, medium, source, sensor,
    pml_inside=True, smooth_p0=False, backend="python", quiet=True,
)
print(f"  OK: p shape = {np.array(result_py['p']).shape}\n")

# ── Inspect the HDF5 file written by CppSimulation ──────────────────
print("=== Inspecting CppSimulation HDF5 output (save_only) ===")
result_save = kspaceFirstOrder(
    kgrid, medium, source, sensor,
    pml_inside=True, smooth_p0=False, backend="cpp",
    save_only=True, data_path="/tmp/kwave_mfe",
)
with h5py.File(result_save["input_file"], "r") as f:
    p_flag = int(f["p_source_flag"][()].item())
    nt_val = int(f["Nt"][()].item())
    p_input_shape = f["p_source_input"].shape
    print(f"  p_source_flag  = {p_flag}")
    print(f"  Nt             = {nt_val}")
    print(f"  p_source_input = {p_input_shape}")
    print()
    if p_flag != nt_val:
        print(f"  BUG: p_source_flag={p_flag} but should be Nt={nt_val}")
        print("  The C++ binary uses p_source_flag to validate the dimensions")
        print("  of p_source_input and rejects the file when they don't match.")
        print()
        print("  In the legacy path (kspaceFirstOrder2D), p_source_flag is set to")
        print("  len(source.p[0]) via kWaveSimulation.source_p, not just bool(has_p).")

# ── Verify that the legacy kspaceFirstOrder2D sets p_source_flag correctly ──
print("\n=== Verifying legacy API (kspaceFirstOrder2D) ===")
import warnings
from kwave.kspaceFirstOrder2D import kspaceFirstOrder2D
from kwave.options.simulation_options import SimulationOptions
from kwave.options.simulation_execution_options import SimulationExecutionOptions

source2 = kSource()
source2.p = np.random.randn(Ns, Nt).astype(np.float32)
source2.p_mask = sensor_mask
source2.p_mode = "dirichlet"

sim_opts = SimulationOptions(
    pml_inside=True, smooth_p0=False, save_to_disk=True,
    save_to_disk_exit=True, data_cast="single", data_path="/tmp/kwave_mfe_old",
)
exec_opts = SimulationExecutionOptions(
    is_gpu_simulation=False, delete_data=False, verbose_level=0, show_sim_log=False,
)

import os
import glob

os.makedirs("/tmp/kwave_mfe_old", exist_ok=True)
for f in glob.glob("/tmp/kwave_mfe_old/*.h5"):
    os.remove(f)

kgrid2 = kWaveGrid(N=N, spacing=dx)
kgrid2.makeTime(1500)

with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    try:
        kspaceFirstOrder2D(
            kgrid=kgrid2, medium=medium, source=source2,
            sensor=kSensor(mask=np.ones(N)),
            simulation_options=sim_opts, execution_options=exec_opts,
        )
    except SystemExit:
        pass

for f in glob.glob("/tmp/kwave_mfe_old/*.h5"):
    with h5py.File(f, "r") as hf:
        if "p_source_flag" in hf:
            old_flag = int(hf["p_source_flag"][()].item())
            old_nt = int(hf["Nt"][()].item())
            print(f"  Legacy p_source_flag = {old_flag}")
            print(f"  Legacy Nt            = {old_nt}")
            if old_flag == old_nt:
                print("  ✓ Legacy API correctly sets p_source_flag = Nt")
            break

Metadata

Metadata

Assignees

No one assigned

    Labels

    bugSomething isn't working

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions