I was simulating a structure suspended in air with varying density (o.5, 1, 1.5, 2 gm/cm3) and speed (1.2 power of density) . I was changing value of CFL to make simulations stable when value density value was at 1.5 and 2 as it started to break, without changing it. However, I noticed that I could also change the reference sound speed to do the same thing, to make the simulation more stable. Interestingly, it decreases computational time drastically.
1. Changing Ref sound speed: 1500 m/s, computing time: 9 mins.
2. Changing CFL to 0.04, computing time: 1hr 1 minutes (when CFL is changed, reference speed changes automatically to 577 m/s.)
I calculated transmission loss in both cases. In case 1, I am getting 80 dB with formula used in the code. In case 2, I am getting the same value 0f 80 dB.
Does this mean changing ref speed is more efficient than CFL?
Any insight on this would be much appreciated. Thank you.
Here is my code.
clearvars;
close all;
% create the computational grid
Nx = 3032 ;
Ny = 1496;
dx = 6.6e-06;
dy = dx;
kgrid = kWaveGrid(Nx, dx, Ny, dy);
%load image
img1 = loadImage('1.bmp'); %https://i.imgur.com/hmBLu0q.png
%resize
img1 = imbinarize(resize(img1, [Nx, Ny]));
img_indices1 = find(img1 ==1);
% define the medium properties
%Medium: Air
medium.density = 10* ones(Nx, Ny); % [kg/m^3]
medium.sound_speed = 330*ones(Nx, Ny); % [m/s]
medium.alpha_coeff = 1.64 * ones(Nx,Ny);
medium.alpha_power = 2;
medium.BonA = 20* ones(Nx,Ny);
% create the time array
cfl = 0.05;
t_end = 0.5e-4; % [s]
kgrid.makeTime(max(medium.sound_speed(:)), cfl);
% define source mask for a linear transducer with an odd number of elements
x_offset = 50; % [grid points]
sensor_start = round(Ny *0.1);
sensor_end = round(Ny*0.9);
source.p_mask = makeLine(Nx, Ny, [1 sensor_start], [1 sensor_end]);
% define the properties of the tone burst used to drive the transducer
sampling_freq = 1/kgrid.dt; % [Hz]
tone_burst_freq = 0.5e6; % [Hz]
tone_burst_cycles = 3;
source_strength = 1; %[Pa]
%Appending the source
source.p = source_strength * toneBurst(sampling_freq, tone_burst_freq, tone_burst_cycles);
%Sensor
x_offset1 = round(Nx*0.4); % [grid points]
x_offset2 = round(Nx*0.6);
width = 10; % [grid points]
sensor.mask = zeros(Nx, Ny);
sensor.mask(x_offset1, Ny/2 - width/2 + 1:Ny/2 + width/2) = 1;
sensor.mask(x_offset2, Ny/2 - width/2 + 1:Ny/2 + width/2) = 1;
% define the frequency response of the sensor elements
center_freq = 0.5e6; % [Hz]
bandwidth = 80; % [%]
sensor.frequency_response = [center_freq, bandwidth];
% assign the input options
input_args = { 'PlotPML', false,'PlotLayout', false...
'DataCast', 'gpuArray-single', 'PMLInside',false,...
'RecordMovie', false, 'MovieName', 'example_movie',...
'MovieProfile', 'MPEG-4', 'MovieArgs', {'FrameRate', 10} };
counter = 0;
for i = 0.5:0.5:2
clear sensor_data;
%Defining properties of structure as scaling function
medium.density(img_indices1) = i*100;
a = medium.density(img_indices1);
medium.sound_speed(img_indices1)= a .^ 1.2;
b = mean(medium.sound_speed(img_indices1));
medium.alpha_coeff(img_indices1) = 10;
medium.BonA(img_indices1) = 400;
%reference.sound_speed = 1500;
%run simulations
[sensor_data] = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
% compute the amplitude spectra of the recorded time series
[~, sensor_data1] = spect(sensor_data(1,:), 1/kgrid.dt);
[f, sensor_data2] = spect(sensor_data(2,:), 1/kgrid.dt);
% applying gaussian filter to the data
sigma = 3; % pick sigma value for the gaussian
gaussFilter = gausswin(6*sigma + 1)';
gaussFilter = gaussFilter / sum(gaussFilter); % normalize
filteredY1 = conv(sensor_data1, gaussFilter, 'same');
filteredY2 = conv(sensor_data2, gaussFilter, 'same');
max_Y1 = max(filteredY1);
max_Y2 = max(filteredY2);
loss1 = -10* log((max_Y2)/ max_Y1) ;
loss = round(loss1);
disp(loss);
counter = counter +1;
%graphs
subplot(2,2,counter);
plot(f*1e-6, filteredY1, '-k',...
f*1e-6, filteredY2, '-r', 'linewidth', 1);
xlim([0 4]);
xlabel('Frequency (MHz)', 'FontWeight', 'bold');
ylabel('Amplitude (a.u)','FontWeight', 'bold');
legend('Sensor-1', 'Sensor-2', 'Fontsize', 7,'FontWeight', 'bold','FontSize', 9);
title(['Loss = ', num2str(loss), ' dB'], 'FontSize', 9, 'FontWeight', 'bold');
end