From 94bb0b90d8ed4bf8571f472e33310159ffb7b6a3 Mon Sep 17 00:00:00 2001 From: Tobias Hienzsch Date: Thu, 22 Aug 2024 21:37:34 +0200 Subject: [PATCH] [python] Add tests for sim3d engine --- data/materials/sabine_03.h5 | 2 +- data/materials/sabine_04.h5 | 2 +- src/python/pffdtd/analysis/room_modes.py | 126 +++++++++++---------- src/python/pffdtd/materials/build.py | 4 +- src/python/pffdtd/sim3d/process_outputs.py | 91 ++++++++++----- src/python/pffdtd/sim3d/sim_fdtd.py | 11 +- src/python/test/test_sim3d.py | 84 ++++++++++++++ 7 files changed, 222 insertions(+), 98 deletions(-) create mode 100644 src/python/test/test_sim3d.py diff --git a/data/materials/sabine_03.h5 b/data/materials/sabine_03.h5 index 68dcc95..9314609 100644 --- a/data/materials/sabine_03.h5 +++ b/data/materials/sabine_03.h5 @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:2453cdc3bd6455d93da182bb5b5e5de501da90e8cb207ac2c11a94d3c95127f2 +oid sha256:4022ba75e7cb9fc0d514365f59c95428e74efdd31e0978be044d42838c9b6618 size 2072 diff --git a/data/materials/sabine_04.h5 b/data/materials/sabine_04.h5 index 68dcc95..f852e9e 100644 --- a/data/materials/sabine_04.h5 +++ b/data/materials/sabine_04.h5 @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:2453cdc3bd6455d93da182bb5b5e5de501da90e8cb207ac2c11a94d3c95127f2 +oid sha256:d2f10df62936b59865d1ab9414c3926f916c4c9f1498394974aef4eb3026cf0d size 2072 diff --git a/src/python/pffdtd/analysis/room_modes.py b/src/python/pffdtd/analysis/room_modes.py index 09a7159..8160189 100644 --- a/src/python/pffdtd/analysis/room_modes.py +++ b/src/python/pffdtd/analysis/room_modes.py @@ -12,7 +12,7 @@ from pffdtd.common.plot import plot_styles -def _find_nearest(array, value): +def find_nearest(array, value): array = np.asarray(array) idx = (np.abs(array - value)).argmin() return array[idx] @@ -61,6 +61,22 @@ def room_mode_kind(m, n, p): return "oblique" +def room_modes(L, W, H, max_order=6): + modes = [] + for m in range(max_order+1): + for n in range(max_order+1): + for p in range(max_order+1): + if m+n+p > 0: + modes.append({ + "m": m, + "n": n, + "p": p, + "frequency": room_mode(L, W, H, m, n, p) + }) + + return sorted(modes, key=lambda x: x['frequency']) + + def frequency_spacing_index(modes): psi = 0.0 num = len(modes) @@ -76,59 +92,33 @@ def frequency_spacing_index(modes): return psi / (num-1) -def main(): - parser = argparse.ArgumentParser() - parser.add_argument('filename', nargs='*') - parser.add_argument('--data_dir', type=str, help='run directory') - parser.add_argument('--fmax', type=float, default=200.0) - parser.add_argument('--fmin', type=float, default=1.0) - parser.add_argument('--modes', type=int, default=10) - args = parser.parse_args() - - directory = args.data_dir - paths = args.filename +def main( + *, + data_dir=None, + filename=None, + fmin=None, + fmax=None, + width=None, + length=None, + height=None, + num_modes=None, + plot=True, +): + directory = data_dir + paths = filename if not paths: paths = collect_wav_paths(directory, "*_out_normalised.wav") - scale = 0.4 - - # Tobi - L = 6.0 * scale - W = 3.65 * scale - H = 3.12 * scale - - # Optimal ratio A - # L = 6 * scale - # W = 4.96 * scale - # H = 4.14 * scale - - # Optimal ratio B - L = 7 * scale - W = 5.19 * scale - H = 3.70 * scale - - # Worst ratio (Cube) - # W = L - # H = L + L = length + W = width + H = height A = L*W V = A*H S = 2*(L*W+L*H+W*H) - modes = [] - max_order = 6 - for m in range(max_order+1): - for n in range(max_order+1): - for p in range(max_order+1): - if m+n+p > 0: - modes.append({ - "m": m, - "n": n, - "p": p, - "frequency": room_mode(L, W, H, m, n, p) - }) - - modes = sorted(modes, key=lambda x: x['frequency'])[:25] + max_order = 8 + modes = room_modes(L, W, H, max_order=max_order)[:25] print(f"{L=:.3f}m {W=:.3f}m {H=:.3f}m") print(f"{A=:.2f}m^2 {S=:.2f}m^2 {V=:.2f}m^3") @@ -150,8 +140,8 @@ def main(): for file in paths: file = pathlib.Path(file) fs, buf = wavfile.read(file) - fmin = args.fmin if args.fmin > 0 else 1 - fmax = args.fmax if args.fmax > 0 else fs/2 + fmin = fmin if fmin > 0 else 1 + fmax = fmax if fmax > 0 else fs/2 nfft = (2**iceil(np.log2(buf.shape[0])))*2 spectrum = np.fft.rfft(buf, nfft) @@ -163,22 +153,23 @@ def main(): dB_max = np.max(dB) peaks, _ = find_peaks(dB, width=2) - mode_peaks = freqs[peaks] + measured_mode_freqs = freqs[peaks] - print(mode_peaks[:10]) + print(measured_mode_freqs[:10]) plt.plot(freqs, dB, linestyle='-', label=f'{file.stem[:4]}') - if args.modes > 0: - mode_freqs = [mode['frequency'] for mode in modes][:args.modes] - plt.vlines(mode_freqs, dB_max-80, dB_max+10, + if num_modes > 0: + calculated_mode_freqs = [mode['frequency'] + for mode in modes][:num_modes] + plt.vlines(calculated_mode_freqs, dB_max-80, dB_max+10, colors='#AAAAAA', linestyles='--', label="Modes") if len(paths) == 1: - plt.plot(mode_peaks, dB[peaks], 'r.', + plt.plot(measured_mode_freqs, dB[peaks], 'r.', markersize=10, label='Peaks') - for mode in mode_freqs: - nearest = _find_nearest(mode_peaks, mode) + for mode in calculated_mode_freqs: + nearest = find_nearest(measured_mode_freqs, mode) print(f"{mode=:03.3f} Hz - {nearest=:03.3f} Hz = {mode-nearest:03.3f} Hz") plt.title("") @@ -190,8 +181,27 @@ def main(): plt.grid(which='minor', color='#DDDDDD', linestyle=':', linewidth=0.5) plt.minorticks_on() plt.legend() - plt.show() + + if plot: + plt.show() + + return calculated_mode_freqs, measured_mode_freqs if __name__ == "__main__": - main() + parser = argparse.ArgumentParser() + parser.add_argument('filename', nargs='*') + parser.add_argument('--data_dir', type=str, help='run directory') + parser.add_argument('--fmax', type=float, default=200.0) + parser.add_argument('--fmin', type=float, default=1.0) + parser.add_argument('--modes', type=int, default=10) + args = parser.parse_args() + + main( + data_dir=args.data_dir, + filename=args.filename, + fmin=args.fmin, + fmax=args.fmax, + num_modes=args.modes, + plot=True + ) diff --git a/src/python/pffdtd/materials/build.py b/src/python/pffdtd/materials/build.py index 85c3ad0..5fac740 100644 --- a/src/python/pffdtd/materials/build.py +++ b/src/python/pffdtd/materials/build.py @@ -28,8 +28,8 @@ def main(): #freq-independent impedance from Sabine abs coefficient write_freq_ind_mat_from_Yn(convert_Sabs_to_Yn(0.01),filename=Path(write_folder / 'sabine_01.h5')) write_freq_ind_mat_from_Yn(convert_Sabs_to_Yn(0.02),filename=Path(write_folder / 'sabine_02.h5')) - write_freq_ind_mat_from_Yn(convert_Sabs_to_Yn(0.02),filename=Path(write_folder / 'sabine_03.h5')) - write_freq_ind_mat_from_Yn(convert_Sabs_to_Yn(0.02),filename=Path(write_folder / 'sabine_04.h5')) + write_freq_ind_mat_from_Yn(convert_Sabs_to_Yn(0.03),filename=Path(write_folder / 'sabine_03.h5')) + write_freq_ind_mat_from_Yn(convert_Sabs_to_Yn(0.04),filename=Path(write_folder / 'sabine_04.h5')) write_freq_ind_mat_from_Yn(convert_Sabs_to_Yn(0.05),filename=Path(write_folder / 'sabine_05.h5')) write_freq_ind_mat_from_Yn(convert_Sabs_to_Yn(0.1),filename=Path(write_folder / 'sabine_1.h5')) write_freq_ind_mat_from_Yn(convert_Sabs_to_Yn(0.2),filename=Path(write_folder / 'sabine_2.h5')) diff --git a/src/python/pffdtd/sim3d/process_outputs.py b/src/python/pffdtd/sim3d/process_outputs.py index 144c907..bd5dc5c 100644 --- a/src/python/pffdtd/sim3d/process_outputs.py +++ b/src/python/pffdtd/sim3d/process_outputs.py @@ -303,6 +303,53 @@ def save_h5(self): h5f.create_dataset('Fs_f', data=self.Fs_f) h5f.close() + +def process_outputs( + *, + data_dir=None, + resample_Fs=None, + fcut_lowcut=None, + N_order_lowcut=None, + fcut_lowpass=None, + N_order_lowpass=None, + symmetric_lowpass=None, + air_abs_filter=None, + save_wav=None, + plot_raw=None, + plot=None, +): + po = ProcessOutputs(data_dir) + + po.initial_process(fcut=fcut_lowcut,N_order=N_order_lowcut) + + if resample_Fs: + po.resample(resample_Fs) + + if fcut_lowpass>0: + po.apply_lowpass(fcut=fcut_lowpass,N_order=N_order_lowpass,symmetric=symmetric_lowpass) + + #these are only needed if you're simulating with fmax >1kHz, but generally fine to use + if air_abs_filter.lower() == 'modal': #best, but slowest + po.apply_modal_filter() + elif air_abs_filter.lower() == 'stokes': #generally fine for indoor air + po.apply_stokes_filter() + elif air_abs_filter.lower() == 'ola': #fastest, but not as recommended + po.apply_ola_filter() + + po.save_h5() + + if save_wav: + po.save_wav() + + plt.rcParams.update(plot_styles) + + if plot_raw: + po.plot_raw_outputs() + + if plot or plot_raw: + po.plot_filtered_outputs() + po.show_plots() + def main(): import argparse parser = argparse.ArgumentParser() @@ -331,37 +378,19 @@ def main(): args = parser.parse_args() - po = ProcessOutputs(args.data_dir) - - po.initial_process(fcut=args.fcut_lowcut,N_order=args.N_order_lowcut) - - if args.resample_Fs: - po.resample(args.resample_Fs) - - if args.fcut_lowpass>0: - po.apply_lowpass(fcut=args.fcut_lowpass,N_order=args.N_order_lowpass,symmetric=args.symmetric_lowpass) - - #these are only needed if you're simulating with fmax >1kHz, but generally fine to use - if args.air_abs_filter.lower() == 'modal': #best, but slowest - po.apply_modal_filter() - elif args.air_abs_filter.lower() == 'stokes': #generally fine for indoor air - po.apply_stokes_filter() - elif args.air_abs_filter.lower() == 'ola': #fastest, but not as recommended - po.apply_ola_filter() - - po.save_h5() - - if args.save_wav: - po.save_wav() - - plt.rcParams.update(plot_styles) - - if args.plot_raw: - po.plot_raw_outputs() - - if args.plot or args.plot_raw: - po.plot_filtered_outputs() - po.show_plots() + process_outputs( + data_dir=args.data_dir, + resample_Fs=args.resample_Fs, + fcut_lowcut=args.fcut_lowcut, + N_order_lowcut=args.N_order_lowcut, + fcut_lowpass=args.fcut_lowpass, + N_order_lowpass=args.N_order_lowpass, + symmetric_lowpass=args.symmetric_lowpass, + air_abs_filter=args.air_abs_filter, + save_wav=args.save_wav, + plot_raw=args.plot_raw, + plot=args.plot, + ) if __name__ == '__main__': main() diff --git a/src/python/pffdtd/sim3d/sim_fdtd.py b/src/python/pffdtd/sim3d/sim_fdtd.py index bbe7644..a3d36b8 100644 --- a/src/python/pffdtd/sim3d/sim_fdtd.py +++ b/src/python/pffdtd/sim3d/sim_fdtd.py @@ -44,6 +44,12 @@ def __init__(self,data_dir,energy_on=False,nthreads=None): self.print(f'numba set for {nthreads=}') nb.set_num_threads(nthreads) + self.load_h5_data() + self.setup_mask() + self.allocate_mem() + self.set_coeffs() + self.checks() + def print(self,fstring): print(f'--ENGINE: {fstring}') @@ -921,11 +927,6 @@ def main(): assert args.draw_backend=='mayavi' eng = SimEngine(args.data_dir,energy_on=args.energy,nthreads=args.nthreads) - eng.load_h5_data() - eng.setup_mask() - eng.allocate_mem() - eng.set_coeffs() - eng.checks() if args.plot: eng.run_plot(draw_backend=args.draw_backend,json_model=args.json_model) else: diff --git a/src/python/test/test_sim3d.py b/src/python/test/test_sim3d.py new file mode 100644 index 0000000..3df3128 --- /dev/null +++ b/src/python/test/test_sim3d.py @@ -0,0 +1,84 @@ +import pytest + +from pffdtd.analysis.room_modes import find_nearest +from pffdtd.analysis.room_modes import main as room_modes +from pffdtd.materials.adm_funcs import write_freq_ind_mat_from_Yn, convert_Sabs_to_Yn +from pffdtd.sim3d.room_builder import RoomBuilder +from pffdtd.sim3d.sim_setup import sim_setup +from pffdtd.sim3d.sim_fdtd import SimEngine +from pffdtd.sim3d.process_outputs import process_outputs + + +@pytest.mark.parametrize("fcc,dx_scale,tolerance", [(True, 3, 7), (False, 2, 3.9)]) +def test_sim3d(tmp_path, fcc, dx_scale, tolerance): + fmin = 20 + fmax = 300 + ppw = 10.5 + dx = 343/(fmax*ppw) + sim_dir = tmp_path + cpu_dir = sim_dir/'cpu' + model_file = sim_dir/'model.json' + material = 'sabine_02.h5' + num_modes = 25 + + L = 2.8 + W = 2.076 + H = 1.48 + + offset = dx*dx_scale + room = RoomBuilder(W, L, H, wall_color=[255, 255, 255]) + room.add_source("S1", [offset, offset, offset]) + room.add_receiver("R1", [W-offset, L-offset, H-offset]) + room.build(model_file) + + write_freq_ind_mat_from_Yn(convert_Sabs_to_Yn(0.02), sim_dir / material) + + sim_setup( + model_json_file=model_file, + mat_folder=sim_dir, + mat_files_dict={ + 'Ceiling': material, + 'Floor': material, + 'Walls': material, + }, + diff_source=True, + duration=3.75, + fcc_flag=fcc, + fmax=fmax, + PPW=ppw, + insig_type='impulse', + save_folder=cpu_dir, + ) + + eng = SimEngine(cpu_dir) + eng.run_all(1) + eng.save_outputs() + + process_outputs( + data_dir=cpu_dir, + resample_Fs=48_000, + fcut_lowcut=fmin, + N_order_lowcut=4, + fcut_lowpass=fmax, + N_order_lowpass=8, + symmetric_lowpass=True, + air_abs_filter="none", + save_wav=True, + plot_raw=False, + plot=False, + ) + + actual, measured = room_modes( + data_dir=cpu_dir, + fmin=fmin, + fmax=fmax, + num_modes=num_modes, + plot=False, + width=W, + length=L, + height=H, + ) + + for mode in actual[:num_modes]: + nearest = find_nearest(measured, mode) + assert abs(mode-nearest) < tolerance