Hi,
I've got a simple 1D simulation I'm running to examine ultrasound quarter-wavelength matching layers, and basically it works for some grid spacings (5um and above) and goes unstable for smaller ones (specifically 1um), which is a bit counter-intuitive to me. I am likely doing something wrong, but I can't seem to find how changing grid spacing would cause this.
Note, the variable to change in the below code is dx. For a value of 5e-6 the sim works, and for 1e-6 it is unstable.
%% Create 1D Line Array
dx = 5e-6; %grid point spacing
freq = 5e6; %drive frequency of pressure pulse in Hz
%Define Aluminum Properties
alum.c = 6299;
alum.d = 2698;
alum.lambda = alum.c/freq;
%Define Matching Layer Properties
match_freq = 6.8e6;
pary.c = 2135;
pary.d = 1241.5;
pary.lambda = pary.c/match_freq;
%Define Water Properties
water.c = 1540;
water.d = 1000;
water.lambda = pary.c/freq;
%Define length of domains in wavelengths
alum.num_wave = 10;
water.num_wave = 10;
%Calculate number of elements per domain, given dx and the above domain
%lengths defined
alum.n = floor(alum.num_wave*alum.lambda/dx);
pary.n = floor(pary.lambda/(4*dx));
water.n = floor(water.num_wave*water.lambda/dx);
Nx = alum.n+pary.n+water.n;
%% Define Grid for 1D sim
kgrid = makeGrid(alum.n+pary.n+water.n, dx);
medium.density = [alum.d*ones(1,alum.n), pary.d*ones(1, pary.n), water.d*ones(1, water.n)];
medium.sound_speed = [alum.c*ones(1,alum.n), pary.c*ones(1, pary.n), water.c*ones(1, water.n)];
dt = 1e-9; %time step in seconds
t_end = alum.num_wave+water.num_wave; %Time for wave to travel in both mediums
t_end = t_end/freq;
%Define time array
t_array = 0:dt:t_end;
kgrid.t_array = t_array;
%Definine signal properties
source_freq = freq; %Hz
source_mag = 1; %Pa
source_cycles = 4; %number of cycles
%Define time array for source
t_length = floor((source_cycles*(1/source_freq))/dt);
t_length = min(t_length, length(kgrid.t_array));
%Define source values for pressure source
source.p = source_mag*sin(2*pi*source_freq*kgrid.t_array(1:t_length));
source.p = [source.p zeros(1, floor((2/source_freq)/dt))];
source.p = filterTimeSeries(kgrid, medium, source.p);
%Define mask where pressure occurs
source.p_mask = [zeros(1, floor(alum.n/2)),1];
source.p_mask = [source.p_mask, zeros(1, Nx-length(source.p_mask))]';
sensor.mask = [floor(alum.n/2),alum.n+floor(water.n/2)].';
sensor.record = {'p_final', 'p_max', 'p_rms', 'p'};
%% Run the Simultion
display_mask = source.p_mask;
Arraytype = 'single'; %use gpuArray-single or gpuArray-double for gpu calcs
%Arraytype = 'gpuArray-single'; %use gpuArray-single or gpuArray-double for gpu calcs
% assign the input options
input_args = {'DataCast', Arraytype,'PlotSim', true, 'DisplayMask', display_mask, 'PMLInside', false, 'PlotPML', false};
%input_args = {'DataCast', 'single','PlotSim', true, 'DisplayMask', display_mask, 'PMLInside', false, 'PlotPML', false, 'PlotScale', [-1 1]};
% run the simulation
sensor_data = kspaceFirstOrder1D(kgrid, medium, source, sensor, input_args{:});