The following code is numerically unstable and blows up. I'm not sure why.
The same code 'works' if I change to an additive source. That said, the pressures inside the additive source region are very inhomogeneous, which I do not understand.
This instability doesn't make much sense to me and I'm having a hard time developing confidence in the model. Any feedback would be appreciated.
Best,
Bill Walker
clear
cfl = 0.1; % Courant-Friedrichs-Lewy number
%% Set geometry
N = 64;
Nx = N;
Ny = Nx;
Nz = N;
dx = 0.2e-3;
dy = dx;
dz = dx;
% number of samples for the lossy boundary layers
n_bnd = 6;
%%
x = ((-(Nx-1)/2):((Nx-1)/2))*dx;
y = x;
z = ((-(Nz-1)/2):((Nz-1)/2))*dx;
z = z;
[X,Y,Z] = meshgrid(x,y,z);
R = ((X).^2+(Y).^2).^0.5;
kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);
%% Set pulse characteristics and timing
f0 = 0.5e6;
N_cycles = 8;
t_max = 20e-6;
%% Define the object regions
% PZT 5A Source
r0(1) = 3.29e-3 / 2;
r1(1) = 9e-3 / 2;
z0(1) = -1e-3;
z1(1) = 1e-3 + eps;
material(1) = 1;
%% Fill out the geometry matrix
mat = zeros(size(X));
for j = 1:length(material),
ind = find((r0(j) < R) & (r1(j) >= R) & (z0(j) < Z) & (z1(j) >= Z));
mat(ind) = material(j);
end
%% Set the material properties
% define the properties of the surrounding medium
medium.sound_speed_compression = 2*2^2*330*ones(size(mat)); % [m/s] ----- Made up
medium.sound_speed_shear = 0*ones(size(mat)); % [m/s] ----- Made up
medium.density = 1.225*ones(size(mat)); % [kg/m^3] ----- Made up
% define the properties of the PZT 5A
ind = find(mat == 1);
medium.sound_speed_compression(ind) = 4350; % [m/s] Szabo book
% poisson's ratio = 0.35
medium.sound_speed_shear(ind) = sqrt(1/(2*(1+0.35)))*4350; % [m/s]
medium.density(ind) = 7750; % [kg/m^3] Szabo book
figure(1)
imagesc(squeeze(medium.density(:,Nx/2,:)));
axis equal
axis off
%% Define the source
% create the time array
[kgrid.t_array, dt] = makeTime(kgrid, max(medium.sound_speed_compression(:)), cfl, t_max);
clear kgrid.t_array
kgrid.t_array = 0:dt:t_max;
% define a ring source element
source.u_mask = zeros(size(mat));
ind = find(mat == 1);
source.u_mask(ind) = 1;
tone_burst_freq = f0;
tone_burst_cycles = N_cycles;
source_signal = sin(2*pi*tone_burst_freq*(0:dt:(tone_burst_cycles/tone_burst_freq)));
source_signal = source_signal.*hann(length(source_signal))';
fir_coeff = fir1(128,1.5e6*dt*2,'low');
source_signal = conv(fir_coeff,source_signal);
ind = find(mat == 1);
source.uz = zeros(length(ind),length(source_signal));
for j = 1:length(ind),
source.uz(j,:) = source_signal;
end
source.u_mode = 'dirichlet';
%% Define the sensor points
% define a series of Cartesian points to collect the data
ind = find(abs(X - x(Nx/2)) < dx/2);
sensor.mask = [kgrid.x(ind)'; kgrid.y(ind)'; kgrid.z(ind)'];
%% Clean up some more memory
clear ind mat R
%% compute the pressure field
% define the field parameters to record
sensor.record = {'p','u'};
% input arguments
input_args = {'DisplayMask', source.u_mask, 'DataCast', 'single', 'CartInterp', 'nearest'};
% run the simulation
sensor_data = pstdElastic3D(kgrid, medium, source, sensor, input_args{:});
save(['simple_test'])
close all