k-wave-python/tests/test_alpha_power_near_unity.py at ccce9194fee8e935199c867ec7efa5dba26d63fd · waltsims/k-wave-python · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
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),
        )