diff --git a/kwave/kWaveSimulation_helper/save_to_disk_func.py b/kwave/kWaveSimulation_helper/save_to_disk_func.py index 178b7a71..bebdec07 100644 --- a/kwave/kWaveSimulation_helper/save_to_disk_func.py +++ b/kwave/kWaveSimulation_helper/save_to_disk_func.py @@ -229,6 +229,9 @@ def grab_medium_props(integer_variables, float_variables, medium, is_elastic_cod integer_variables.absorbing_flag = 0 if medium.absorbing: + from kwave.solvers.validation import warn_cpp_alpha_mode_unsupported + + warn_cpp_alpha_mode_unsupported(medium.alpha_mode, stacklevel=4) if is_elastic_code: # pragma: no cover # add to the variable list float_variables["chi"] = None diff --git a/kwave/solvers/cpp_simulation.py b/kwave/solvers/cpp_simulation.py index 545898b4..fe8d7f72 100644 --- a/kwave/solvers/cpp_simulation.py +++ b/kwave/solvers/cpp_simulation.py @@ -9,6 +9,7 @@ import stat import subprocess import tempfile +import warnings import numpy as np @@ -36,6 +37,10 @@ def __init__(self, kgrid, medium, source, sensor, *, pml_size, pml_alpha, use_sg self.use_sg = use_sg self.ndim = kgrid.dim + from kwave.solvers.validation import warn_cpp_alpha_mode_unsupported + + warn_cpp_alpha_mode_unsupported(medium.alpha_mode, stacklevel=3) + def prepare(self, data_path=None): """Write HDF5 input file. Returns (input_file, output_file).""" if data_path is None: diff --git a/kwave/solvers/kspace_solver.py b/kwave/solvers/kspace_solver.py index fd822816..baa17a3f 100644 --- a/kwave/solvers/kspace_solver.py +++ b/kwave/solvers/kspace_solver.py @@ -411,6 +411,7 @@ def _setup_physics_operators(self): # Absorption/dispersion alpha_coeff_raw = getattr(self.medium, "alpha_coeff", 0) + alpha_mode = self.medium.alpha_mode if not _is_enabled(alpha_coeff_raw): self._absorption = lambda div_u: 0 self._dispersion = lambda rho: 0 @@ -419,17 +420,23 @@ def _setup_physics_operators(self): alpha_power = float(xp.array(getattr(self.medium, "alpha_power", 1.5)).flatten()[0]) alpha_np = 100 * alpha_coeff * (1e-6 / (2 * np.pi)) ** alpha_power / (20 * np.log10(np.e)) + no_absorption = alpha_mode == "no_absorption" + no_dispersion = alpha_mode == "no_dispersion" + if abs(alpha_power - 2.0) < 1e-10: # Stokes - self._absorption = lambda div_u: -2 * alpha_np * self.c0 * self.rho0 * div_u + self._absorption = (lambda div_u: 0) if no_absorption else (lambda div_u: -2 * alpha_np * self.c0 * self.rho0 * div_u) self._dispersion = lambda rho: 0 else: # Power-law with fractional Laplacian tau = -2 * alpha_np * self.c0 ** (alpha_power - 1) - eta = 2 * alpha_np * self.c0**alpha_power * xp.tan(np.pi * alpha_power / 2) nabla1 = self._fractional_laplacian(alpha_power - 2) - nabla2 = self._fractional_laplacian(alpha_power - 1) - # Fractional Laplacian already includes full k-space structure; no kappa needed - self._absorption = lambda div_u: tau * self._diff(self.rho0 * div_u, nabla1) - self._dispersion = lambda rho: eta * self._diff(rho, nabla2) + self._absorption = (lambda div_u: 0) if no_absorption else (lambda div_u: tau * self._diff(self.rho0 * div_u, nabla1)) + + if no_dispersion: + self._dispersion = lambda rho: 0 + else: + eta = 2 * alpha_np * self.c0**alpha_power * xp.tan(np.pi * alpha_power / 2) + nabla2 = self._fractional_laplacian(alpha_power - 1) + self._dispersion = lambda rho: eta * self._diff(rho, nabla2) # Nonlinearity BonA_raw = getattr(self.medium, "BonA", 0) diff --git a/kwave/solvers/validation.py b/kwave/solvers/validation.py index a628b1a1..afd464ad 100644 --- a/kwave/solvers/validation.py +++ b/kwave/solvers/validation.py @@ -47,6 +47,15 @@ def validate_medium(medium, kgrid): power = float(np.asarray(medium.alpha_power).flat[0]) if power < 0 or power > 3: warnings.warn(f"medium.alpha_power={power} is outside typical range [0, 3].", stacklevel=3) + alpha_mode = medium.alpha_mode + if abs(power - 1.0) < 0.05 and alpha_mode != "no_dispersion": + raise ValueError( + f"medium.alpha_power={power} is too close to 1.0. The dispersion term " + f"contains tan(pi*alpha_power/2), which diverges at alpha_power=1. " + f"For backend='python': set medium.alpha_mode='no_dispersion' to disable " + f"the dispersion term. For backend='cpp' or to use both absorption and " + f"dispersion, choose an alpha_power value further from 1.0." + ) def validate_pml(pml_size, kgrid): @@ -97,6 +106,18 @@ def validate_source(source, kgrid): raise ValueError("All velocity source components must either be single (1D) or many (2D) time series.") +def warn_cpp_alpha_mode_unsupported(alpha_mode, stacklevel=3): + """Warn that the C++ backend cannot honor medium.alpha_mode.""" + if alpha_mode in ("no_absorption", "no_dispersion"): + warnings.warn( + f"medium.alpha_mode='{alpha_mode}' is not supported by the C++ backend " + f"and will be silently ignored. The C++ binary always computes both " + f"absorption and dispersion when absorbing_flag > 0. Use backend='python' " + f"to honor alpha_mode.", + stacklevel=stacklevel, + ) + + def validate_sensor(sensor, kgrid): if sensor is None or not hasattr(sensor, "mask") or sensor.mask is None: return diff --git a/tests/test_alpha_power_near_unity.py b/tests/test_alpha_power_near_unity.py new file mode 100644 index 00000000..66dff615 --- /dev/null +++ b/tests/test_alpha_power_near_unity.py @@ -0,0 +1,132 @@ +"""Tests for alpha_power near 1.0 (tan(pi*y/2) singularity).""" +import numpy as np +import pytest + +from kwave.kgrid import kWaveGrid +from kwave.kmedium import kWaveMedium +from kwave.ksensor import kSensor +from kwave.ksource import kSource +from kwave.kspaceFirstOrder import kspaceFirstOrder + + +def _make_sim(): + """Create a minimal 2D simulation with a point source.""" + N = (64, 64) + dx = (1e-4, 1e-4) + kgrid = kWaveGrid(N, dx) + kgrid.makeTime(1500) + + source = kSource() + source.p0 = np.zeros(N) + source.p0[32, 32] = 1.0 + + sensor = kSensor(mask=np.ones(N)) + return kgrid, source, sensor + + +def test_alpha_power_near_unity_raises_without_no_dispersion(): + """alpha_power close to 1.0 must raise when dispersion is enabled.""" + kgrid, source, sensor = _make_sim() + medium = kWaveMedium(sound_speed=1500, density=1000, alpha_coeff=0.5, alpha_power=1.01) + + with pytest.raises(ValueError, match="too close to 1.0"): + kspaceFirstOrder(kgrid, medium, source, sensor, pml_inside=True, quiet=True) + + +def test_alpha_power_near_unity_works_with_no_dispersion(): + """alpha_power close to 1.0 must produce valid output when dispersion is disabled.""" + kgrid, source, sensor = _make_sim() + medium = kWaveMedium( + sound_speed=1500, + density=1000, + alpha_coeff=0.5, + alpha_power=1.01, + alpha_mode="no_dispersion", + ) + + result = kspaceFirstOrder(kgrid, medium, source, sensor, pml_inside=True, quiet=True) + p = np.asarray(result["p"]) + assert not np.any(np.isnan(p)), "Output contains NaN with alpha_power=1.01 and no_dispersion" + assert not np.all(p == 0), "Output is all zeros — absorption had no effect" + + +def test_alpha_mode_no_absorption(): + """alpha_mode='no_absorption' should disable absorption but keep dispersion.""" + kgrid, source, sensor = _make_sim() + medium_lossy = kWaveMedium(sound_speed=1500, density=1000, alpha_coeff=0.5, alpha_power=1.5) + medium_no_abs = kWaveMedium( + sound_speed=1500, + density=1000, + alpha_coeff=0.5, + alpha_power=1.5, + alpha_mode="no_absorption", + ) + + result_lossy = kspaceFirstOrder(kgrid, medium_lossy, source, sensor, pml_inside=True, quiet=True) + result_no_abs = kspaceFirstOrder(kgrid, medium_no_abs, source, sensor, pml_inside=True, quiet=True) + + p_lossy = np.asarray(result_lossy["p"]) + p_no_abs = np.asarray(result_no_abs["p"]) + + assert not np.any(np.isnan(p_no_abs)) + # Disabling absorption should change the output + assert not np.allclose(p_lossy, p_no_abs), "Disabling absorption should change the output" + + +def test_alpha_mode_no_dispersion(): + """alpha_mode='no_dispersion' should disable dispersion but keep absorption.""" + kgrid, source, sensor = _make_sim() + medium_lossy = kWaveMedium(sound_speed=1500, density=1000, alpha_coeff=0.5, alpha_power=1.5) + medium_no_disp = kWaveMedium( + sound_speed=1500, + density=1000, + alpha_coeff=0.5, + alpha_power=1.5, + alpha_mode="no_dispersion", + ) + + result_lossy = kspaceFirstOrder(kgrid, medium_lossy, source, sensor, pml_inside=True, quiet=True) + result_no_disp = kspaceFirstOrder(kgrid, medium_no_disp, source, sensor, pml_inside=True, quiet=True) + + p_lossy = np.asarray(result_lossy["p"]) + p_no_disp = np.asarray(result_no_disp["p"]) + + assert not np.any(np.isnan(p_no_disp)) + # Both should have absorption so similar amplitude, but waveforms differ + assert not np.allclose(p_lossy, p_no_disp), "Disabling dispersion should change the output" + + +def test_alpha_power_normal_range_unaffected(): + """alpha_power=1.5 (well away from singularity) should work as before.""" + kgrid, source, sensor = _make_sim() + medium = kWaveMedium(sound_speed=1500, density=1000, alpha_coeff=0.5, alpha_power=1.5) + + result = kspaceFirstOrder(kgrid, medium, source, sensor, pml_inside=True, quiet=True) + p = np.asarray(result["p"]) + assert not np.any(np.isnan(p)) + assert p.shape[0] > 0 + + +def test_cpp_backend_warns_on_alpha_mode(tmp_path): + """C++ backend should warn when alpha_mode is set (it cannot honor it).""" + kgrid, source, sensor = _make_sim() + medium = kWaveMedium( + sound_speed=1500, + density=1000, + alpha_coeff=0.5, + alpha_power=1.5, + alpha_mode="no_dispersion", + ) + + with pytest.warns(UserWarning, match="not supported by the C\\+\\+ backend"): + kspaceFirstOrder( + kgrid, + medium, + source, + sensor, + pml_inside=True, + quiet=True, + backend="cpp", + save_only=True, + data_path=str(tmp_path), + )