k-Wave User Forum » Topic: Simulation of knife-edge diffraction?

k-Wave User Forum » Topic: Simulation of knife-edge diffraction? http://www.k-wave.org/forum/topic/simulation-of-knife-edge-diffraction Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 17:20:48 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "Simulation of knife-edge diffraction?" http://www.k-wave.org/forum/topic/simulation-of-knife-edge-diffraction#post-6496 Sun, 24 Jun 2018 19:04:58 +0000 Bradley Treeby 6496@http://www.k-wave.org/forum/ HI jkusumaFB, Thanks for the interesting puzzle - unfortunately, it's got us stumped! It's a stability issue related to the heterogeneous density, but we have no idea why it depends on the spatial distribution. If you have a single grid point on the left or top, it is also unstable. Fortunately, it can be solved by reducing the size of the time step, or only having a heterogeneous sound speed - the formula for stability used in the example doesn't account for the heterogeneous density. We'll keep digging and let you know if we figure it out. Brad. jkusumaFB on "Simulation of knife-edge diffraction?" http://www.k-wave.org/forum/topic/simulation-of-knife-edge-diffraction#post-6468 Thu, 17 May 2018 16:43:06 +0000 jkusumaFB 6468@http://www.k-wave.org/forum/ I think I've figured this out -- sort of -- enough to do a small numerical experiment anyway. If I place the single barrier on the right, it works. If I place the single barrier on the left, it crashes, unless I place a tiny barrier on the right, with a wide open slit. It seems to me that this has to do with the solver, which is unfortunately beyond my ability to figure out at this point. :-( jkusumaFB on "Simulation of knife-edge diffraction?" http://www.k-wave.org/forum/topic/simulation-of-knife-edge-diffraction#post-6463 Tue, 15 May 2018 17:19:25 +0000 jkusumaFB 6463@http://www.k-wave.org/forum/ Hi everybody, I'm new to k-wave and I am trying to learn by example -- in this case I'm interested in simulating a single-edge (knife-edge) diffraction simulation. I started with the built-in example "example_tvsp_slit_diffraction.m" that Bradley Treeby wrote, and simply removed one of the two sides of the barrier. When I do so, the simulation seems to become divergent. I tried removing the "Datacast = single" option and it still doesn't work. If I remove all the barriers, the simulation runs fine ... basically a homogeneous media. Can somebody help me understand why the one-sided simulation doesn't run correctly, please? Thanks in advance! <blockquote> % Diffraction Through A Slit Example Example % % This example illustrates the diffraction of a plane acoustic wave through % a slit. It builds on the Monopole Point Source In A Homogeneous % Propagation Medium and Simulating Transducer Field Patterns examples. % % author: Bradley Treeby % date: 1st November 2010 % last update: 20th June 2017 % % This function is part of the k-Wave Toolbox (<a href="http://www.k-wave.org" rel="nofollow">http://www.k-wave.org</a>) % Copyright (C) 2010-2017 Bradley Treeby % This file is part of k-Wave. k-Wave is free software: you can % redistribute it and/or modify it under the terms of the GNU Lesser % General Public License as published by the Free Software Foundation, % either version 3 of the License, or (at your option) any later version. % % k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY % WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS % FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for % more details. % % You should have received a copy of the GNU Lesser General Public License % along with k-Wave. If not, see <http://www.gnu.org/licenses/>. clearvars; % ========================================================================= % SIMULATION % ========================================================================= % create the computational grid and define the bulk medium properties scale = 1; % change to 2 to produce higher resolution images PML_size = 10; % size of the perfectly matched layer [grid points] Nx = scale * 128 - 2 * PML_size; % number of grid points in the x (row) direction Ny = Nx; % number of grid points in the y (column) direction dx = 50e-3 / Nx; % grid point spacing in the x direction [m] dy = dx; % grid point spacing in the y direction [m] c0 = 1500; % [m/s] speed of sound in water rho0 = 1000; % [kg/m^3] kgrid = kWaveGrid(Nx, dx, Ny, dy); % define the ratio between the barrier and background sound speed and % density %barrier_scale = 20; barrier_scale = 20; % lower contrast values = better diffraction contrast %% define the barrier and the source wavelength depending on the example % create a mask of a barrier with a slit slit_thickness = scale * 2; % [grid points] slit_width = scale * 50; % [grid points] slit_x_pos = Nx - Nx/4; % [grid points] slit_y_end = Ny/2; % [grid points] slit_offset = Ny/2 - slit_width/2 - 1; % [grid points] slit_mask = zeros(Nx, Ny); %%% one sided -- crashes? slit_mask(slit_x_pos:slit_x_pos + slit_thickness, 1:slit_y_end) = 1; % left side only %%% two slits -- works % slit_mask(slit_x_pos:slit_x_pos + slit_thickness, 1:slit_offset) = 1; % left side % slit_mask(slit_x_pos:slit_x_pos + slit_thickness, end - 1:end) = 1; % right side % define the source wavelength to be a quarter of the slit size source_wavelength = 0.25 * slit_width * dx; % [m] source_freq = 500e3; % [Hz] source_wavelength = c0/source_freq; % [m] % source_wavelength = 0.0058; % assign the slit to the properties of the propagation medium medium.sound_speed = c0 * ones(Nx, Ny); medium.density = rho0 * ones(Nx, Ny); medium.sound_speed(slit_mask == 1) = barrier_scale * c0; medium.density(slit_mask == 1) = barrier_scale * rho0; % assign the reference sound speed to the background medium medium.sound_speed_ref = c0; %% prepare for simulation % find the time step at the stability limit c_ref = medium.sound_speed_ref; c_max = barrier_scale * c0; k_max = max(kgrid.k(:)); dt_limit = 2 / (c_ref * k_max) * asin(c_ref / c_max); % create the time array, with the time step just below the stability limit dt = 0.95 * dt_limit; % [s] t_end = 40e-6; % [s] kgrid.setTime(round(t_end / dt) + 1, dt); % create a source mask of a single line source.p_mask = zeros(Nx, Ny); source.p_mask(end, :) = 1; % create and filter the time varying sinusoidal source source_mag = 2; source_freq = c0 / source_wavelength; source.p = source_mag * sin(2 * pi * source_freq * kgrid.t_array); source.p = filterTimeSeries(kgrid, medium, source.p); % define the field parameters to record sensor.mask = ones(Nx, Ny); sensor.record = {'u_final', 'p_final'}; %% run simulation and display results % set the input options % input_args = {'PMLInside', false, 'PMLSize', PML_size, 'PlotPML', false, ... % 'DisplayMask', slit_mask, 'DataCast', 'single'}; input_args = {'PMLInside', false, 'PMLSize', PML_size, 'PlotPML', false, ... 'DisplayMask', slit_mask}; % run the simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); % ========================================================================= % VISUALISATION % ========================================================================= % plot the final wave-field figure; mx = max(abs(sensor_data.p_final(:))); sensor_data.p_final(slit_mask == 1) = mx; imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, sensor_data.p_final, [-mx, mx]); colormap(getColorMap); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('p'); % plot the final wave-field figure; mx = max(abs(sensor_data.ux_final(:))); sensor_data.ux_final(slit_mask == 1) = mx; imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, sensor_data.ux_final, [-mx, mx]); colormap(getColorMap); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('ux'); % plot the final wave-field figure; mx = max(abs(sensor_data.uy_final(:))); sensor_data.uy_final(slit_mask == 1) = mx; imagesc(kgrid.y_vec * 1e3, kgrid.x_vec * 1e3, sensor_data.uy_final, [-mx, mx]); colormap(getColorMap); ylabel('x-position [mm]'); xlabel('y-position [mm]'); axis image; title('uy'); </blockquote>