Fix NaN when alpha_power is near 1.0, respect alpha_mode in all paths

Three issues caused NaN output with attenuation exponents near 1.0:

1. The Python solver ignored medium.alpha_mode entirely, always computing
   the dispersion term containing tan(pi*alpha_power/2) which diverges
   at alpha_power=1. Now respects 'no_dispersion' and 'no_absorption'.

2. No validation caught alpha_power values *near* 1.0 (only == 1 exactly).
   Now raises ValueError when |alpha_power - 1| < 0.05 unless
   alpha_mode='no_dispersion'.

3. The C++ backend silently ignored alpha_mode (never written to HDF5).
   Now warns that alpha_mode is unsupported by the C++ binary, directing
   users to backend='python'.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: waltsims

kwave/kWaveSimulation_helper/save_to_disk_func.py

@@ -1,5 +1,6 @@
 import logging
 import os
+import warnings
 
 import numpy as np
 from scipy.io import savemat
@@ -229,6 +230,15 @@ def grab_medium_props(integer_variables, float_variables, medium, is_elastic_cod
         integer_variables.absorbing_flag = 0
 
     if medium.absorbing:
+        alpha_mode = getattr(medium, "alpha_mode", None)
+        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=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,16 @@ def __init__(self, kgrid, medium, source, sensor, *, pml_size, pml_alpha, use_sg
         self.use_sg = use_sg
         self.ndim = kgrid.dim
 
+        alpha_mode = getattr(medium, "alpha_mode", None)
+        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=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 = getattr(self.medium, "alpha_mode", None)
         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 -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,14 @@ 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 = getattr(medium, "alpha_mode", None)
+            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"Set medium.alpha_mode='no_dispersion' to disable the dispersion term, "
+                    f"or use an alpha_power value further from 1.0."
+                )
 
 
 def validate_pml(pml_size, kgrid):