-
Notifications
You must be signed in to change notification settings - Fork 53
Closed
Labels
bugSomething isn't workingSomething isn't working
Description
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
Reactions are currently unavailable
Metadata
Metadata
Assignees
Labels
bugSomething isn't workingSomething isn't working