k-Wave User Forum » Topic: Changing CFL or reference speed to make simulation stable.

k-Wave User Forum » Topic: Changing CFL or reference speed to make simulation stable. http://www.k-wave.org/forum/topic/changing-cfl-or-reference-speed-to-make-simulation-stable Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 22:28:58 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php bencox on "Changing CFL or reference speed to make simulation stable." http://www.k-wave.org/forum/topic/changing-cfl-or-reference-speed-to-make-simulation-stable#post-8568 Wed, 13 Jul 2022 16:52:26 +0000 bencox 8568@http://www.k-wave.org/forum/ Hi sgrghimire, I'm not sure what you mean by reference sound speed. <code>medium.sound_speed_ref</code> is an optional parameter in <code>kspaceFirstOrder2D</code> but you don't seem to be using that. If you have a large impedance mismatch, then the simulation may need many points per wavelength for the simulation to converge. In some cases, it may be possible to make the smaller values of density or sound speed unrealistically large without affecting the reflection coefficient very much at all. This can be a useful trick to reduce the number of points per wavelength required. Best wishes Ben sgrghimire on "Changing CFL or reference speed to make simulation stable." http://www.k-wave.org/forum/topic/changing-cfl-or-reference-speed-to-make-simulation-stable#post-8526 Thu, 05 May 2022 19:05:30 +0000 sgrghimire 8526@http://www.k-wave.org/forum/ 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