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!