k-Wave User Forum » Topic: Instability in shear wave absorption simulation

k-Wave User Forum » Topic: Instability in shear wave absorption simulation http://www.k-wave.org/forum/topic/instability-in-shear-wave-absorption-simulation Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 00:04:28 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php genny_p on "Instability in shear wave absorption simulation" http://www.k-wave.org/forum/topic/instability-in-shear-wave-absorption-simulation#post-6939 Fri, 28 Jun 2019 11:07:46 +0000 genny_p 6939@http://www.k-wave.org/forum/ Hi Brad, Thanks for the suggestion. The simulation seems to actually only be stable with longer time steps (sparser grid spacing and higher CFL number) and when I decrease the time step then the simulation becomes unstable (you can see this by changing the factors in the sample code I provided). I take it then this is not a result of a numerical instability. What other types of instabilities could be possible? Thanks, Genny Bradley Treeby on "Instability in shear wave absorption simulation" http://www.k-wave.org/forum/topic/instability-in-shear-wave-absorption-simulation#post-6894 Sat, 22 Jun 2019 20:18:49 +0000 Bradley Treeby 6894@http://www.k-wave.org/forum/ Hi Genny, For the shear code, we don't have an out of the box function to test for stability with an absorbing medium. For the fluid code, we have an approximate way of checking for the stable time step size (see the <code>checkStability</code> function). It might be possible to find a relationship for the elastic code, but I haven't looked at it. If it is a numerical instability, then keep making the CFL smaller, it should eventually work. Let me know if it doesn't. Brad. genny_p on "Instability in shear wave absorption simulation" http://www.k-wave.org/forum/topic/instability-in-shear-wave-absorption-simulation#post-6851 Wed, 15 May 2019 12:05:54 +0000 genny_p 6851@http://www.k-wave.org/forum/ Hi Brad, Ben, and the k-Wave community, As I mentioned at Photonics West, I am applying k-Wave to model scattering of planar shear waves in heterogeneous tissue, and so far the toolbox has been very easy to use and is working quite well for elastic scattering problems. I recently started trying to model shear wave absorption using the medium.alpha_coeff_shear factor to account for viscoelasticity. The results look quite reasonable for certain cases, but when I increase the absorption factor enough, the result seems to become unstable. I am using a value on the order of 3e6 dB/(MHz^2 cm), since this is around what we measure for shear wave absorption in tissue-mimicking phantoms with wave frequency on the order of 600 Hz, i.e. 1.08 dB/cm. You can see the instability forming in the example code below, where the edges of the plane wave start to increase with an amplitude much larger than the initial wave amplitude. I saw in the Modelling Power Law Absorption Example that you discuss numerical errors (<a href="http://www.k-wave.org/documentation/example_na_modelling_absorption.php#heading3)" rel="nofollow">http://www.k-wave.org/documentation/example_na_modelling_absorption.php#heading3)</a>, but there it is mentioned that the errors can be minimized by reducing the time step. I found the instability I see seems to worsen with a decreased time step (which I implemented by reducing the CFL number). The instability is also worse when I reduce the grid spacing (which in turn also reduced the time step). See the comments in my code below for the factors that influence the instability that I investigated. All the parameters used in the example code are typical for our experiments, except maybe the inital velocity u_0, but there doesn't seem to be any dependence of the instability on this value. This leads me to the following questions: 1. Is this instability related to the numerical errors discussed in the Modelling Power Law Absorption Example? 2. Is there some relationship that tells me for which conditions (up to which absorption) the code works for before becoming unstable? 3. If I have a result that appears stable (i.e. with a low enough medium.alpha_coeff_shear value), how can I check whether the results are reliable, since when I decrease the grid size to check convergence the simulation then blows up? I would be grateful if someone had any insight to any of these questions. Thanks, Genny % ------------------------------------------------------------------------------------- % Modified from Shear Waves And Critical Angle Reflection Example % Simplified code that shows instability during absorption of planar shear wave % Heterogeneities removed to focus on instability % Grid size cropped for quicker simulation % The parameters coded below result in an amplitude that appears to blow up % May 2019 % % If one of the factors is changed to as decribed below (with all other parameters constant): % * decreasing alpha0_s1 to 2e6 the solution appears stable % * increasing alpha0_s1 to 4e6 the solution blows up more % * increasing t_end to 5000 the solution blows up more % * decreasing t_end to 2000 the solution appears stable % * decreasing cfl to 0.1 (decreases dt by a factor of 2) the solution blows up more % * increasing scale to 1.25 (decreases dx & dt by a factor of 0.8) the solution blows up much more % * decreasing scale to 0.75 (increases dx & dt by a factor of 1.3) the solution appears stable % * increasing cp1 to 3000 (decreases dt by a factor of 2) the solution blows up more % * decreasing cp1 to 750 (increases dt by a factor of 2) the solution appears stable % * increasing cs1 to 4.5 the solution blows up more % * decreasing cs1 to 3.5 the solution appears stable % * instability appears insensitive to frequency (300-1800 Hz tested) % * instability appears insensitive to u_0 (0.5 to 1000 tested) close all clearvars % define the medium properties for the outise material cp1 = 1500; % compressional wave speed [m/s] cs1 = 4.35; % shear wave speed [m/s] rho1 = 1000; % density [kg/m^3] alpha0_p1 = 0; % compressional absorption [dB/(MHz^2 cm)] alpha0_s1 = 3e6; % shear absorption [dB/(MHz^2 cm)] (measured approx 1.08 dB/cm at 600 Hz) % create the computational grid scale = 1; % grid spacing factor PML_size = 10; % [grid points] Nx = 40*scale ; % [grid points] Ny = 20*scale ; % [grid points] dx = 0.4e-3/scale; % [m] dy = 0.4e-3/scale; % [m] kgrid = kWaveGrid(Nx, dx, Ny, dy); % create the time array to run until wave passes through entire grid cfl = 0.2; t_end = 4500*1e-6; kgrid.makeTime(max([cp1 cs1]), cfl, t_end); % define the source source_freq = 600; % [Hz] u_0 = 500; % initial velocity source_mask = zeros(Nx, Ny); source_mask(:,1) = 1; source.u_mask = source_mask; source.ux = u_0 * sin(2 * pi * source_freq * kgrid.t_array); source.uy = 0 * sin(2 * pi * source_freq * kgrid.t_array); % define sensor sensor.record = {'u_max_all','u_min_all','u_final'}; % define the medium properties clear medium medium.sound_speed_compression = cp1*ones(Nx, Ny); medium.sound_speed_shear = cs1*ones(Nx, Ny); medium.density = rho1*ones(Nx, Ny); medium.alpha_coeff_compression = alpha0_p1*ones(Nx, Ny); medium.alpha_coeff_shear = alpha0_s1*ones(Nx, Ny); % set the input arguments input_args = {'PMLSize', PML_size, 'PMLAlpha', 2, 'PlotPML', false, ... 'PMLInside', false, 'PlotScale', 'auto', ... 'DataCast', 'single'}; % run the elastic simulation sensor_data_elastic = pstdElastic2D(kgrid, medium, source, sensor, input_args{:}); % visualization ui = u_0*cp1/cs1; % initial velocity figure(1) subplot(2,2,1); imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, sensor_data_elastic.ux_max_all/ui); caxis([0 1.5]) xlabel('y mm');ylabel('x mm');axis image;h = colorbar;ylabel(h, 'Normalized u_x');title('ux max all') subplot(2,2,2); imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, sensor_data_elastic.ux_final/ui); caxis([-1.5 1.5]) xlabel('y mm');ylabel('x mm');axis image;h = colorbar;ylabel(h, 'Normalized u_x');title('ux final') subplot(2,2,3); imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, sensor_data_elastic.uy_max_all/ui); caxis([0 1.5]) xlabel('y mm');ylabel('x mm');axis image;h = colorbar;ylabel(h, 'Normalized u_y');title('uy max all') subplot(2,2,4); imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, sensor_data_elastic.uy_final/ui); caxis([-1.5 1.5]) xlabel('y mm');ylabel('x mm');axis image;h = colorbar;ylabel(h, 'Normalized u_y');title('uy final')