Fix alpha mode near unity

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

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:

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)

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

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),
+        )