Skip to content

Commit

Permalink
passing tests and fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
@djps
djps committed Sep 10, 2024
1 parent 696b5d9 commit 936dfb3
Show file tree
Hide file tree
Showing 3 changed files with 136 additions and 93 deletions.
41 changes: 24 additions & 17 deletions kwave/kWaveSimulation_helper/create_storage_variables.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -217,8 +217,17 @@ def create_shift_operators(record: Recorder, record_old: Recorder, kgrid: kWaveG
record.y_shift_neg = ifftshift(np.exp(-1j * kgrid.k_vec.y * kgrid.dy / 2)).T
elif kgrid.dim == 3:
record.x_shift_neg = ifftshift(np.exp(-1j * kgrid.k_vec.x * kgrid.dx / 2))
record.y_shift_neg = ifftshift(np.exp(-1j * kgrid.k_vec.y * kgrid.dy / 2)).T
record.z_shift_neg = np.transpose(ifftshift(np.exp(-1j * kgrid.k_vec.z * kgrid.dz / 2)), [1, 2, 0])
record.y_shift_neg = ifftshift(np.exp(-1j * kgrid.k_vec.y * kgrid.dy / 2))
record.z_shift_neg = ifftshift(np.exp(-1j * kgrid.k_vec.z * kgrid.dz / 2))

record.x_shift_neg = np.expand_dims(record.x_shift_neg, axis=-1)

record.y_shift_neg = np.expand_dims(record.y_shift_neg, axis=0)

record.z_shift_neg = np.squeeze(record.z_shift_neg)
record.z_shift_neg = np.expand_dims(record.z_shift_neg, axis=0)
record.z_shift_neg = np.expand_dims(record.z_shift_neg, axis=0)

else:
if kgrid.dim == 1:
record.x_shift_neg = 1
Expand Down Expand Up @@ -446,29 +455,27 @@ def compute_triangulation_points(flags, kgrid, record, mask):
Barycentric coordinates)
"""

if kgrid.dim == 1:
# align sensor data as a column vector to be the same as kgrid.x_vec
# so that calls to interp return data in the correct dimension
sensor_x = np.reshape((mask, (-1, 1)))

elif kgrid.dim == 2:
sensor_x = mask[0, :]
sensor_y = mask[1, :]

elif kgrid.dim == 3:
sensor_x = mask[0, :]
sensor_y = mask[1, :]
sensor_z = mask[2, :]


if not flags.binary_sensor_mask:

if kgrid.dim == 1:

# assign pseudonym for Cartesain grid points in 1D (this is later used for data casting)
record.grid_x = kgrid.x_vec

else:

if kgrid.dim == 1:
# align sensor data as a column vector to be the same as kgrid.x_vec
# so that calls to interp return data in the correct dimension
sensor_x = np.reshape((mask, (-1, 1)))
elif kgrid.dim == 2:
sensor_x = mask[0, :]
sensor_y = mask[1, :]
elif kgrid.dim == 3:
sensor_x = mask[0, :]
sensor_y = mask[1, :]
sensor_z = mask[2, :]

# update command line status
print(' calculating Delaunay triangulation...')

Expand Down
49 changes: 34 additions & 15 deletions kwave/kWaveSimulation_helper/scale_source_terms_func.py
View file
Original file line number Diff line number Diff line change
@@ -1,14 +1,18 @@
import logging

from copy import deepcopy

import numpy as np

from kwave.kgrid import kWaveGrid
from kwave.ksource import kSource
from kwave.utils.dotdictionary import dotdict


def scale_source_terms_func(c0, dt, kgrid: kWaveGrid, source, p_source_pos_index,
s_source_pos_index, u_source_pos_index,
def scale_source_terms_func(c0, dt, kgrid: kWaveGrid, source,
p_source_pos_index,
s_source_pos_index,
u_source_pos_index,
transducer_input_signal, flags: dotdict):
"""
Subscript for the first-order k-Wave simulation functions to scale source terms to the correct units.
Expand Down Expand Up @@ -211,31 +215,46 @@ def scale_stress_sources(source, c0, flags, dt, dx, N, s_source_pos_index):
"""

source.sxx = scale_stress_source(source, c0, flags.source_sxx, flags.source_p0, source.sxx, dt, N, dx, s_source_pos_index)
source.syy = scale_stress_source(source, c0, flags.source_syy, flags.source_p0, source.syy, dt, N, dx, s_source_pos_index)
source.szz = scale_stress_source(source, c0, flags.source_szz, flags.source_p0, source.szz, dt, N, dx, s_source_pos_index)
# source.sxy = scale_stress_source(source, c0, flags.source_sxy, True, source.sxy, dt, N, dx, s_source_pos_index)
# source.sxz = scale_stress_source(source, c0, flags.source_sxz, True, source.sxz, dt, N, dx, s_source_pos_index)
# source.syz = scale_stress_source(source, c0, flags.source_syz, True, source.syz, dt, N, dx, s_source_pos_index)
if flags.source_sxx:
print("scale source.sxx")
source.sxx = scale_stress_source(source, c0, flags.source_sxx, flags.source_p0, deepcopy(source.sxx), dt, N, dx, s_source_pos_index)
if flags.source_syy:
print("scale source.syy")
source.syy = scale_stress_source(source, c0, flags.source_syy, flags.source_p0, deepcopy(source.syy), dt, N, dx, s_source_pos_index)
if flags.source_szz:
print("scale source.szz")
source.szz = scale_stress_source(source, c0, flags.source_szz, flags.source_p0, deepcopy(source.szz), dt, N, dx, s_source_pos_index)
if flags.source_sxy:
source.sxy = scale_stress_source(source, c0, flags.source_sxy, flags.source_p0, source.sxy, dt, N, dx, s_source_pos_index)
if flags.source_sxz:
source.sxz = scale_stress_source(source, c0, flags.source_sxz, flags.source_p0, source.sxz, dt, N, dx, s_source_pos_index)
if flags.source_syz:
source.syz = scale_stress_source(source, c0, flags.source_syz, flags.source_p0, source.syz, dt, N, dx, s_source_pos_index)


def scale_stress_source(source, c0, is_source_exists, is_p0_exists, source_s, dt, N, dx, s_source_pos_index):
if is_source_exists:
if source.s_mode == "dirichlet" or is_p0_exists:
source_s = source_s / N
print("source.s_mode == dirichlet or is_p0_exists")
else:
if c0.size == 1:
print("if c0.size == 1")
# compute the scale parameter based on the homogeneous sound
# speed
source_s = source_s * (2 * dt * c0 / (N * dx))

else:
# compute the scale parameter seperately for each source
# position based on the sound speed at that position
ind = range(source_s[:, 0].size)
mask = s_source_pos_index.flatten("F")[ind]
scale = (2.0 * dt * np.expand_dims(c0.ravel(order="F")[mask.ravel(order="F")], axis=-1) ) / (N * dx)
source_s[ind, :] *= scale
s_index = range(0, np.shape(source_s)[0])
print("s_index:", s_index)
mask = s_source_pos_index.flatten("F")[s_index]
print("mask:", mask)
scale = (2.0 * dt * np.expand_dims(c0.ravel(order="F")[mask.ravel(order="F")], axis=-1)) / (N * dx)
print("scale:", scale, "from:", dt, N, dx, np.expand_dims(c0.ravel(order="F")[mask.ravel(order="F")], axis=-1))
print(np.shape(source_s[s_index, :]), np.shape(scale))
source_s[s_index, :] = source_s[s_index, :] * scale

return source_s

Expand Down Expand Up @@ -278,9 +297,9 @@ def scale_velocity_sources(flags, source, kgrid, c0, dt, dx, dy, dz, u_source_po
# flags.source_ux, source.u_mode, source.ux, kgrid, c0, dt, dx, u_source_pos_index, flags.nonuniform_grid
# )

source.ux = scale_velocity_source(flags.source_ux, source.u_mode, source.ux, c0, dt, u_source_pos_index, dx)
source.uy = scale_velocity_source(flags.source_uy, source.u_mode, source.uy, c0, dt, u_source_pos_index, dy)
source.uz = scale_velocity_source(flags.source_uz, source.u_mode, source.uz, c0, dt, u_source_pos_index, dz)
source.ux = scale_velocity_source(flags.source_ux, source.u_mode, deepcopy(source.ux), c0, dt, u_source_pos_index, dx)
source.uy = scale_velocity_source(flags.source_uy, source.u_mode, deepcopy(source.uy), c0, dt, u_source_pos_index, dy)
source.uz = scale_velocity_source(flags.source_uz, source.u_mode, deepcopy(source.uz), c0, dt, u_source_pos_index, dz)


# def scale_velocity_source_x(is_source_ux, source_u_mode, source_val, kgrid, c0, dt, dx, u_source_pos_index, is_nonuniform_grid):
Expand Down
139 changes: 78 additions & 61 deletions tests/test_pstd_elastic_3d_check_mpml_stability.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,12 +15,16 @@
from kwave.utils.signals import tone_burst
from kwave.utils.mapgen import make_spherical_section

import scipy.io as sio
# import scipy.io as sio

def test_pstd_elastic_3d_check_mpml_stability():

test_pass: bool = True

# mat_contents = sio.loadmat('C:/Users/dsinden/dev/octave/k-Wave/testing/unit/mpml_stability.mat')

# mpml_smask = mat_contents['s_mask']

# create the computational grid
PML_SIZE: int = 10
Nx: int = 80 - 2 * PML_SIZE
Expand All @@ -35,43 +39,52 @@ def test_pstd_elastic_3d_check_mpml_stability():
sound_speed_compression = 1500.0 * np.ones((Nx, Ny, Nz)) # [m/s]
sound_speed_shear = np.zeros((Nx, Ny, Nz)) # [m/s]
density = 1000.0 * np.ones((Nx, Ny, Nz)) # [kg/m^3]
# define the properties of the lower layer of the propagation medium
sound_speed_compression[Nx // 2 - 1:, :, :] = 2000.0 # [m/s]
sound_speed_shear[Nx // 2 - 1:, :, :] = 1000.0 # [m/s]
density[Nx // 2 - 1:, :, :] = 1200.0 # [kg/m^3]
medium = kWaveMedium(sound_speed_compression,
sound_speed_compression=sound_speed_compression,
sound_speed_shear=sound_speed_shear,
density=density)

# define the properties of the lower layer of the propagation medium
medium.sound_speed_compression[Nx // 2 - 1:, :, :] = 2000.0 # [m/s]
medium.sound_speed_shear[Nx // 2 - 1:, :, :] = 1000.0 # [m/s]
medium.density[Nx // 2 - 1:, :, :] = 1200.0 # [kg/m^3]

# create the time array
cfl = 0.3 # Courant-Friedrichs-Lewy number
cfl = 0.3 # Courant-Friedrichs-Lewy number
t_end = 8e-6 # [s]
kgrid.makeTime(medium.sound_speed_compression.max(), cfl, t_end)

# define the source mask
s_rad: int = 15
s_height: int = 8
offset: int = 14
offset: int = 15
ss, _ = make_spherical_section(s_rad, s_height)

source = kSource()
ss_width = np.shape(ss)[1]
ss_half_width: int = np.floor(ss_width // 2).astype(int)
ss_width: int = np.shape(ss)[1]
ss_half_width: int = np.floor(ss_width / 2).astype(int)
y_start_pos: int = Ny // 2 - ss_half_width - 1
y_end_pos: int = y_start_pos + ss_width
z_start_pos: int = Nz // 2 - ss_half_width - 1
z_end_pos: int = z_start_pos + ss_width

source.s_mask = np.zeros((Nx, Ny, Nz), dtype=int)

source.s_mask[offset:s_height + offset, y_start_pos:y_end_pos, z_start_pos:z_end_pos] = ss.astype(int)
source.s_mask[:, :, Nz // 2 - 1:] = 0
source.s_mask[:, :, Nz // 2 - 1:] = int(0)

# print("diff_mat_pml_smask:", np.shape(mpml_smask), np.shape(source.s_mask),
# np.max(np.abs(mpml_smask - source.s_mask)), np.argmax(np.abs(mpml_smask - source.s_mask)),
# np.nonzero(np.abs(mpml_smask - source.s_mask)),
# np.max(np.abs(source.s_mask)), np.max(np.abs(mpml_smask)))

# define the source signal
source.sxx = tone_burst(1.0 / kgrid.dt, 1e6, 3)
source.syy = source.sxx
source.szz = source.sxx
fs = 1.0 / kgrid.dt
source.sxx = tone_burst(sample_freq=fs, signal_freq=1e6, num_cycles=3)

# print(source.sxx)

source.syy = deepcopy(source.sxx)
source.szz = deepcopy(source.sxx)

# define sensor
sensor = kSensor()
Expand All @@ -88,7 +101,7 @@ def test_pstd_elastic_3d_check_mpml_stability():
multi_axial_PML_ratio=0.0)

# run the simulations
sensor_data_pml = pstd_elastic_3d(deepcopy(kgrid),
sensor_data_pml = pstd_elastic_3d(kgrid=deepcopy(kgrid),
medium=deepcopy(medium),
source=deepcopy(source),
sensor=deepcopy(sensor),
Expand All @@ -103,70 +116,74 @@ def test_pstd_elastic_3d_check_mpml_stability():
binary_sensor_mask=True,
multi_axial_PML_ratio=0.1)

sensor_data_mpml = pstd_elastic_3d(deepcopy(kgrid),
sensor_data_mpml = pstd_elastic_3d(kgrid=deepcopy(kgrid),
medium=deepcopy(medium),
source=deepcopy(source),
sensor=deepcopy(sensor),
simulation_options=deepcopy(simulation_options_mpml))

# check magnitudes
pml_max = np.max([sensor_data_pml.ux_final, sensor_data_pml.uy_final, sensor_data_pml.uz_final])
mpml_max = np.max([sensor_data_mpml.ux_final,sensor_data_mpml.uy_final, sensor_data_mpml.uz_final])
mpml_max = np.max([sensor_data_mpml.ux_final, sensor_data_mpml.uy_final, sensor_data_mpml.uz_final])

# set reference magnitude (initial source)
ref_max = 1.0 / np.max(medium.sound_speed_shear * medium.density)


mat_contents = sio.loadmat('mpml_stability.mat')

pml_ux_final = mat_contents['pml_ux_final']
pml_uy_final = mat_contents['pml_uy_final']
pml_uz_final = mat_contents['pml_uz_final']

mpml_ux_final = mat_contents['mpml_ux_final']
mpml_uy_final = mat_contents['mpml_uy_final']
mpml_uz_final = mat_contents['mpml_uz_final']

diff_mat_pml_ux_final = np.max(np.abs(sensor_data_pml.ux_final - pml_ux_final))
diff_mat_pml_uy_final = np.max(np.abs(sensor_data_pml.uy_final - pml_uy_final))
diff_mat_pml_uz_final = np.max(np.abs(sensor_data_pml.uz_final - pml_uz_final))

diff_mat_mpml_ux_final = np.max(np.abs(sensor_data_mpml.ux_final - mpml_ux_final))
diff_mat_mpml_uz_final = np.max(np.abs(sensor_data_mpml.uy_final - mpml_uy_final))
diff_mat_mpml_uy_final = np.max(np.abs(sensor_data_mpml.uz_final - mpml_uz_final))

print("diff_mat_pml_ux_final:", diff_mat_pml_ux_final,
np.max(np.abs(sensor_data_pml.ux_final)), np.argmax(np.abs(sensor_data_pml.ux_final)),
np.max(np.abs(pml_ux_final)), np.argmax(np.abs(pml_ux_final)))
print("diff_mat_pml_uy_final:", diff_mat_pml_uy_final,
np.max(np.abs(sensor_data_pml.uy_final)), np.argmax(np.abs(sensor_data_pml.uy_final)),
np.max(np.abs(pml_uy_final)), np.argmax(np.abs(pml_uy_final)))
print("diff_mat_pml_uz_final:", diff_mat_pml_uz_final,
np.max(np.abs(sensor_data_pml.uz_final)), np.argmax(np.abs(sensor_data_pml.uz_final)),
np.max(np.abs(pml_uz_final)), np.argmax(np.abs(pml_uz_final)))

print("diff_mat_mpml_ux_final:", diff_mat_mpml_ux_final,
np.max(np.abs(sensor_data_mpml.ux_final)), np.argmax(np.abs(sensor_data_mpml.ux_final)),
np.max(np.abs(mpml_ux_final)), np.argmax(np.abs(mpml_ux_final)))
print("diff_mat_mpml_uy_final:", diff_mat_mpml_uy_final,
np.max(np.abs(sensor_data_mpml.uy_final)), np.argmax(np.abs(sensor_data_mpml.uy_final)),
np.max(np.abs(mpml_uy_final)), np.argmax(np.abs(mpml_uy_final)))
print("diff_mat_mpml_uz_final:", diff_mat_mpml_uz_final,
np.max(np.abs(sensor_data_mpml.uz_final)), np.argmax(np.abs(sensor_data_mpml.uz_final)),
np.max(np.abs(mpml_uz_final)), np.argmax(np.abs(mpml_uz_final)))


ref_max = 1.0 / (np.max(medium.sound_speed_shear) * np.max(medium.density))

# pml_ux_final = mat_contents['pml_ux_final']
# pml_uy_final = mat_contents['pml_uy_final']
# pml_uz_final = mat_contents['pml_uz_final']

# mpml_ux_final = mat_contents['mpml_ux_final']
# mpml_uy_final = mat_contents['mpml_uy_final']
# mpml_uz_final = mat_contents['mpml_uz_final']

# diff_mat_pml_ux_final = np.max(np.abs(sensor_data_pml.ux_final - pml_ux_final))
# diff_mat_pml_uy_final = np.max(np.abs(sensor_data_pml.uy_final - pml_uy_final))
# diff_mat_pml_uz_final = np.max(np.abs(sensor_data_pml.uz_final - pml_uz_final))

# diff_mat_mpml_ux_final = np.max(np.abs(sensor_data_mpml.ux_final - mpml_ux_final))
# diff_mat_mpml_uz_final = np.max(np.abs(sensor_data_mpml.uy_final - mpml_uy_final))
# diff_mat_mpml_uy_final = np.max(np.abs(sensor_data_mpml.uz_final - mpml_uz_final))

# print("diff_mat_pml_ux_final:", diff_mat_pml_ux_final,
# np.max(np.abs(sensor_data_pml.ux_final)), np.argmax(np.abs(sensor_data_pml.ux_final)),
# np.max(np.abs(pml_ux_final)), np.argmax(np.abs(pml_ux_final)))
# print("diff_mat_pml_uy_final:", diff_mat_pml_uy_final,
# np.max(np.abs(sensor_data_pml.uy_final)), np.argmax(np.abs(sensor_data_pml.uy_final)),
# np.max(np.abs(pml_uy_final)), np.argmax(np.abs(pml_uy_final)))
# print("diff_mat_pml_uz_final:", diff_mat_pml_uz_final,
# np.max(np.abs(sensor_data_pml.uz_final)), np.argmax(np.abs(sensor_data_pml.uz_final)),
# np.max(np.abs(pml_uz_final)), np.argmax(np.abs(pml_uz_final)))

# print("diff_mat_mpml_ux_final:", diff_mat_mpml_ux_final,
# np.max(np.abs(sensor_data_mpml.ux_final)), np.argmax(np.abs(sensor_data_mpml.ux_final)),
# np.max(np.abs(mpml_ux_final)), np.argmax(np.abs(mpml_ux_final)))
# print("diff_mat_mpml_uy_final:", diff_mat_mpml_uy_final,
# np.max(np.abs(sensor_data_mpml.uy_final)), np.argmax(np.abs(sensor_data_mpml.uy_final)),
# np.max(np.abs(mpml_uy_final)), np.argmax(np.abs(mpml_uy_final)))
# print("diff_mat_mpml_uz_final:", diff_mat_mpml_uz_final,
# np.max(np.abs(sensor_data_mpml.uz_final)), np.argmax(np.abs(sensor_data_mpml.uz_final)),
# np.max(np.abs(mpml_uz_final)), np.argmax(np.abs(mpml_uz_final)))

# mat_pml_max = np.max([pml_ux_final, pml_uy_final, pml_uz_final])
# mat_mpml_max = np.max([mpml_ux_final, mpml_uy_final, mpml_uz_final])

# print("mat_pml_max < ref_max " + str(mat_pml_max < ref_max) + ", mat_pml_max: " + str(mat_pml_max) + ", ref_max: " + str(ref_max))
# print("mat_mpml_max > ref_max " + str(mat_mpml_max > ref_max) + ", mpml_max: " + str(mat_mpml_max) + ", ref_max: " + str(ref_max))

# print("pml_max < ref_max " + str(pml_max < ref_max) + ", pml_max: " + str(pml_max) + ", ref_max: " + str(ref_max))
# print("mpml_max > ref_max " + str(mpml_max > ref_max) + ", mpml_max: " + str(mpml_max) + ", ref_max: " + str(ref_max))

# check results - the test should fail if the pml DOES work (i.e., it
# doesn't become unstable), or if the m-pml DOESN'T work (i.e., it does
# become unstable)
# become unstable). The pml should not work and the mpml should.
if pml_max < ref_max:
test_pass = False
assert test_pass, "pml_max < ref_max " + str(pml_max < ref_max) + ", pml_max: " + str(pml_max) + ", ref_max: " + str(ref_max)

if mpml_max > ref_max:
test_pass = False
assert test_pass, "mpml_max > ref_max " + str(mpml_max > ref_max) + ", pml_max: " + str(pml_max) + ", ref_max: " + str(ref_max)
assert test_pass, "mpml_max > ref_max " + str(mpml_max > ref_max) + ", mpml_max: " + str(mpml_max) + ", ref_max: " + str(ref_max)



0 comments on commit 936dfb3

Please sign in to comment.