k-Wave User Forum » Topic: Calculation in a loop
http://www.k-wave.org/forum/topic/calculation-in-a-loop
Support for the k-Wave MATLAB toolboxen-USSun, 03 Nov 2024 15:17:45 +0000http://bbpress.org/?v=1.0.2<![CDATA[Search]]>q
http://www.k-wave.org/forum/search.php
lfura on "Calculation in a loop"
http://www.k-wave.org/forum/topic/calculation-in-a-loop#post-9135
Thu, 29 Aug 2024 07:11:40 +0000lfura9135@http://www.k-wave.org/forum/<p>I would like to simulate the oscillations of the microbubbles exposed to ultrasound. So I thought that I will simulate the distribution of acoustic pressure in k-wave in every step and then take this pressure and put it into Rayleigh-Plesset equation as the external acoustic pressure. The result would be a change in the radius of the microbubble, which could be treated as a wave source and introduced as a wave source in the next step of the k-wave simulation. In this way, I wanted to link the k-wave simulations to the Rayleigh-Plesset equation in a loop for the entire simulation space. Of course, there is still the problem of scattering by the microbubble, but I wanted to leave that out for now.
</p>bencox on "Calculation in a loop"
http://www.k-wave.org/forum/topic/calculation-in-a-loop#post-9131
Wed, 28 Aug 2024 17:15:34 +0000bencox9131@http://www.k-wave.org/forum/<p>Hi Ifura,</p>
<p>No, it's not the same. What is it that you would like to simulate?</p>
<p>Best wishes<br />
Ben
</p>lfura on "Calculation in a loop"
http://www.k-wave.org/forum/topic/calculation-in-a-loop#post-9122
Tue, 06 Aug 2024 12:54:14 +0000lfura9122@http://www.k-wave.org/forum/<p>Hi, </p>
<p>Is it possible in k-wave to calculate the sound pressure distribution after N time steps by running simulations in a loop after 1 time step? The loop runs kspaceFirstOrder where at each grid node the wave source is the value from the previous time step in the loop. It seems to me that this should work, but maybe I'm missing something and it's not working. I have attached my code below.</p>
<p>Kind regards</p>
<pre><code>pml_size = 20; % [grid points]
% define the grid parameters
Nx = 512 - 2 * pml_size; % [grid points]
Ny = 512 - 2 * pml_size; % [grid points]
dx = 10e-6; % [m]
dy = 10e-6; % [m]
% create the computational grid
kgrid = kWaveGrid(Nx, dx, Ny, dy);
%%
% define the properties of the propagation medium
medium.sound_speed = 1480; % [m/s]
medium.density = 1000; % [kg/m^3]
medium.alpha_coeff = 0.02; % [dB/(MHz^y cm)]
medium.alpha_power = 1.1;
% define the source parameters
diameter = 4e-3; % [m]
radius = 3e-3; % [m]
freq = 1.05e6; % [Hz]
amp = 0.1e6; % [Pa]
mask= logical(zeros(Nx,Ny));
mask(Nx/2,Ny/2) = 1;
source.p_mask = zeros(Nx,Ny);
source.p_mask(mask) = 1;
% calculate the time step using an integer number of points per period
ppw = medium.sound_speed / (freq * dx); % points per wavelength
cfl = 0.3; % cfl number
ppp = ceil(ppw / cfl); % points per period
T = 1 / freq; % period [s]
dt = T / ppp; % time step [s]
% calculate the number of time steps to reach steady state
t_end = sqrt( kgrid.x_size.^2 + kgrid.y_size.^2 ) / medium.sound_speed;
Nt = round(t_end / dt);
kgrid.setTime(Nt, dt);
% define the input signal
input_signal = createCWSignals(kgrid.t_array, freq, amp, 0);
% create the time array
kgrid.setTime(1, dt);
source.p = zeros(1,1);
source.p = input_signal(:,1);
% set the sensor mask to cover the entire grid
sensor.mask = ones(Nx, Ny);
sensor.record = {'p'};
% record one point
sensor.record_start_index = 1;
% set the input arguements
input_args = {'PMLInside', false, 'PlotPML', false, 'DisplayMask', ...
'off', 'PlotScale', [-1, 1] * amp, ...
'PlotSim', false, 'DataCast', 'gpuArray-single'};
% run the acoustic simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
% % reshape the data
p = reshape(sensor_data.p, Nx, Ny,1);
%%
source.p_mask = ones(Nx,Ny);
f = waitbar(0);
for i=2:length(input_signal)
source.p = p(:,:,i-1);
source.p = source.p(:);
source.p(mask(:)) = input_signal(:,i);
% run the acoustic simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
% % reshape the data
p (:,:,i) = reshape(sensor_data.p, Nx, Ny,1);
waitbar(i/(length(input_signal)),f)
end
close(f)</code></pre>