k-Wave User Forum » Topic: Stability of elastic simulations

k-Wave User Forum » Topic: Stability of elastic simulations http://www.k-wave.org/forum/topic/stability-of-elastic-simulations Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 12:11:49 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "Stability of elastic simulations" http://www.k-wave.org/forum/topic/stability-of-elastic-simulations#post-5681 Thu, 15 Sep 2016 16:43:32 +0000 Bradley Treeby 5681@http://www.k-wave.org/forum/ Hi Bellamy, The instability seems to be related to the very high values of absorption, and the PML. If you turn either the absorption or the PML off, your simulation seems to run ok. It's possibly because the M-PML that the elastic code uses is derived for lossless media. Unfortunately I don't have time right now to dig any deeper, but I will definitely add it to the to-do list. As an aside, a couple of other tips for your code. First, if you set your grid size to have small prime factors, it will be faster (e.g., use <code>Nx = 128*scale+448- 2*PML_size;</code>). Second, you will use less memory if you just define one source input (this is used for all the source elements), without the trailing zeros. Brad. Bellamy on "Stability of elastic simulations" http://www.k-wave.org/forum/topic/stability-of-elastic-simulations#post-5649 Mon, 29 Aug 2016 14:47:42 +0000 Bellamy 5649@http://www.k-wave.org/forum/ This is an example of one of the simulations I've tested: When plotting the pressure record at one of the points, there is a growing ringing recorded as time goes on, after a certain amount of time. clear all; % ========================================================================= % SIMULATION PARAMETERS % ========================================================================= % change scale to 2 to reproduce the higher resolution figures used in the % help file scale = 1; lambda=1540/.8E6; % create the computational grid PML_size = 10; % size of the PML in grid points Nx = 128*scale+450- 2*PML_size; % number of grid points in the x direction Ny = 192*scale - 2*PML_size; % number of grid points in the y direction dx = lambda/8; % grid point spacing in the x direction [m] dy = dx; % grid point spacing in the y direction [m] kgrid = makeGrid(Nx, dx, Ny, dy); % define the medium properties cp1 = 1540; % compressional wave speed [m/s] cs1 = 0; % shear wave speed [m/s] rho1 = 1000; % density [kg/m^3] alpha0_p1 = 0.1; % compressional wave absorption [dB/(MHz^2 cm)] alpha0_s1 = 0.1; % shear wave absorption [dB/(MHz^2 cm)] cp2 = 3250; % compressional wave speed [m/s] cs2 = 1685; % shear wave speed [m/s] rho2 = 2533; % density [kg/m^3] alpha0_p2 = 38; % compressional wave absorption [dB/(MHz^2 cm)] alpha0_s2 = 38; % shear wave absorption [dB/(MHz^2 cm)] cp3 = 2139; % compressional wave speed [m/s] cs3 = 1100; % shear wave speed [m/s] rho3 = 1500; % density [kg/m^3] alpha0_p3 = 38; % compressional wave absorption [dB/(MHz^2 cm)] alpha0_s3 = 38; % shear wave absorption [dB/(MHz^2 cm)] % create the time array cfl = 0.08; t_end = 1.4e-4; kgrid.t_array= makeTime(kgrid, cp2, cfl, t_end); % define position of heterogeneous slab slab = zeros(Nx, Ny); slab(20:Nx, :) = 1; slab3 = zeros(Nx, Ny); slab3(40:60, :) = 1; % define the source properties source_freq = .8E6; % [Hz] source_cycles = 3; % number of tone burst cycles source_focus_dist = 5*scale; % position of focus inside slab [grid points] source_slab_dist = 5*scale; % distance between source and slab [grid points] %source_mask = makeCircle(Nx, Ny, Nx/2 + source_focus_dist, Ny/2, 50*scale, pi/3); source_mask(1:Nx, 1:Ny) = 0; source_mask(10, 20:30) = 1; sensor.mask=ones(Nx,Ny); % define the sensor to record the maximum particle velocity everywhere sensor.record = {'p','p_max_all'}; % set the input arguments input_args = {'PMLSize', PML_size, 'PMLAlpha', 2, 'PlotPML', false, ... 'PMLInside', false, 'PlotScale', 'auto', ... 'DisplayMask', 'off', 'DataCast', 'single'}; % assign the source source.p_mask = source_mask; source.p = toneBurst(1/kgrid.dt, source_freq, source_cycles); % ========================================================================= % ELASTIC SIMULATION % ========================================================================= % define the medium properties clear medium medium.sound_speed_compression = cp1*ones(Nx, Ny); medium.sound_speed_compression(slab == 1) = cp2; medium.sound_speed_compression(slab3 == 1) = cp3; medium.sound_speed_shear = cs1*ones(Nx, Ny); medium.sound_speed_shear(slab == 1) = cs2; medium.sound_speed_shear(slab3 == 1) = cs3; medium.density = rho1*ones(Nx, Ny); medium.density(slab == 1) = rho2; medium.density(slab3 == 1) = rho3; medium.alpha_coeff_compression = alpha0_p1*ones(Nx, Ny); medium.alpha_coeff_compression(slab == 1) = alpha0_p2; medium.alpha_coeff_compression(slab3 == 1) = alpha0_p3; medium.alpha_coeff_shear = alpha0_s1*ones(Nx, Ny); medium.alpha_coeff_shear(slab == 1) = alpha0_s2; medium.alpha_coeff_shear(slab3 == 1) = alpha0_s3; % assign the source clear source source.u_mask = source_mask; signal = toneBurst(1/kgrid.dt, source_freq, 3,'SignalLength',5000); signal= filterTimeSeries(kgrid, medium, signal); for i=1:sum(source.u_mask(:)) source.ux(i,:)=signal; end % run the elastic simulation sensor_data_elastic = pstdElastic2D(kgrid, medium, source, sensor, input_args{:}); p(:,:,:) = reshape(sensor_data_elastic.p, Nx,Ny,[]); plot(squeeze(p(10,25,:))) I've tried reducing the cfl until I run into out of memory errors, and I still see this instability forming (albeit with a smaller Nx). Does anyone know the source of this instability, and what I can do to get rid of it? Thanks! Bellamy on "Stability of elastic simulations" http://www.k-wave.org/forum/topic/stability-of-elastic-simulations#post-5647 Tue, 23 Aug 2016 20:31:34 +0000 Bellamy 5647@http://www.k-wave.org/forum/ Hi, I'm having some trouble with instabilities in my grid when I simulate a strongly heterogeneous medium, with some high density (max 2700 kg/m^3) and high attenuation values (~38 dB/cm MHz^2). I am using the PSTD elastic code. I've tried reducing the grid size to lambda/12 (actually in testing using the shear example in 2D I found that for the same cfl that it's more stable for lambda/8 than lambda/12 for some reason..?) and I've tried it with lamba/8 grid spacing and cfl=0.03 using the max compressional speed of sound (3350m/s) as the reference. I've been testing out some simple situations, modifying the shear wave example and I find that it looks ok if I keep the duration of the simulation to be the time it takes for the ultrasound to travel from one side of the grid to the other, but if I make it much longer than that, the simulation becomes unstable. I've tried additive and dirichlet modes, and it was still unstable. Is this an issue with the PML, or am I missing something? Thanks!