From 936dfb33f51d953aaf2ea085bcdaecbe0ac7a478 Mon Sep 17 00:00:00 2001 From: "Sinden, David" Date: Tue, 10 Sep 2024 22:46:12 +0200 Subject: [PATCH] passing tests and fixes --- .../create_storage_variables.py | 41 +++--- .../scale_source_terms_func.py | 49 ++++-- ...st_pstd_elastic_3d_check_mpml_stability.py | 139 ++++++++++-------- 3 files changed, 136 insertions(+), 93 deletions(-) diff --git a/kwave/kWaveSimulation_helper/create_storage_variables.py b/kwave/kWaveSimulation_helper/create_storage_variables.py index 81f2e6b6..a8c1a3f2 100644 --- a/kwave/kWaveSimulation_helper/create_storage_variables.py +++ b/kwave/kWaveSimulation_helper/create_storage_variables.py @@ -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 @@ -446,22 +455,8 @@ 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) @@ -469,6 +464,18 @@ def compute_triangulation_points(flags, kgrid, record, mask): 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...') diff --git a/kwave/kWaveSimulation_helper/scale_source_terms_func.py b/kwave/kWaveSimulation_helper/scale_source_terms_func.py index 4f980fd5..ecc23333 100644 --- a/kwave/kWaveSimulation_helper/scale_source_terms_func.py +++ b/kwave/kWaveSimulation_helper/scale_source_terms_func.py @@ -1,5 +1,7 @@ import logging +from copy import deepcopy + import numpy as np from kwave.kgrid import kWaveGrid @@ -7,8 +9,10 @@ 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. @@ -211,20 +215,31 @@ 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)) @@ -232,10 +247,14 @@ def scale_stress_source(source, c0, is_source_exists, is_p0_exists, source_s, dt 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 @@ -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): diff --git a/tests/test_pstd_elastic_3d_check_mpml_stability.py b/tests/test_pstd_elastic_3d_check_mpml_stability.py index db830b0d..30eb813f 100644 --- a/tests/test_pstd_elastic_3d_check_mpml_stability.py +++ b/tests/test_pstd_elastic_3d_check_mpml_stability.py @@ -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 @@ -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() @@ -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), @@ -103,7 +116,7 @@ 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), @@ -111,62 +124,66 @@ def test_pstd_elastic_3d_check_mpml_stability(): # 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)