http://www.k-wave.org/forum/topic/analysis-of-the-power-loss-along-propagation-in-lossless-medium Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 23:31:55 +0000 http://bbpress.org/?v=1.0.2 q http://www.k-wave.org/forum/search.php http://www.k-wave.org/forum/topic/analysis-of-the-power-loss-along-propagation-in-lossless-medium#post-5132 Wed, 08 Jul 2015 16:35:03 +0000 Bradley Treeby 5132@http://www.k-wave.org/forum/ <p>Great! Nice to see the non-staggered velocity calculations working in practice. </p> http://www.k-wave.org/forum/topic/analysis-of-the-power-loss-along-propagation-in-lossless-medium#post-5130 Wed, 08 Jul 2015 15:13:59 +0000 Anthony 5130@http://www.k-wave.org/forum/ <p>Sorry, I made a stupid mistake in my 2D example: everything is perfect. </p> <pre><code>clear all close all clc %% PARAMETERS TO CHANGE medium.alpha_coeff = 0; % liver: 0.37 medium.alpha_power = 0; % liver: 1.266 % create the computational grid Nx = 128-40; Ny = 128-40; dx = 0.1e-3; dy = 0.1e-3; kgrid = makeGrid(Nx, dx, Ny, dy); % define the properties of the propagation medium medium.sound_speed = 1500; medium.density = 1000; medium.BonA = 0; % create time array kgrid.t_array = makeTime(kgrid, medium.sound_speed); % create initial pressure distribution using makeDisc R = Ny*0.45; Rap = 28/38*R; source.p_mask = zeros(Nx, Ny); xc = round(Nx/2); yc = round(Ny/2); xgrid = (1:Nx)-xc; ygrid = ((1:Ny)-yc)&#39;; rgrid = sqrt(bsxfun(@plus,xgrid.^2,ygrid.^2)); source.p_mask(logical(bsxfun(@times,(abs(rgrid-R)&lt;1/2), bsxfun(@times,ygrid&lt;0,abs(xgrid)&lt;Rap)))) = 1; source.p_mode = &#39;dirichlet&#39;; figure(1) clf imagesc(source.p_mask) axis equal % define source source_mag = 0.1e6; source_f0 = 3e6; source.p = source_mag*cos(2*pi*kgrid.t_array*source_f0); % zero pad and filter source.p = filterTimeSeries(kgrid, medium, [source.p, zeros(1, length(source.p))]); % define sensor sensor.mask = ones(Nx, Ny); sensor.record = {&#39;p&#39;, &#39;u_non_staggered&#39;}; % run the simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor,... &#39;DisplayMask&#39;, &#39;off&#39;, &#39;PlotScale&#39;, [-0.5, 0.5]*source_mag,&#39;PMLInside&#39;,false,&#39;PlotPML&#39;,false); % time shift ux = timeShift(sensor_data.ux_non_staggered, kgrid.dt); uy = timeShift(sensor_data.uy_non_staggered, kgrid.dt); % Intensity NpointsToKeep = 50; % in order to have an integer number of points) Ix = reshape(sensor_data.p.*ux,[Nx,Ny,size(ux,2)]); Iy = reshape(sensor_data.p.*uy,[Nx,Ny,size(uy,2)]); Ix = mean(Ix(:,:,end-NpointsToKeep+1:end),3); Iy = mean(Iy(:,:,end-NpointsToKeep+1:end),3); % Display figure(2) clf sidePowerLoss = ( cumsum(-Iy(:,1),1)+cumsum(Iy(:,end),1) ) * dx; compensatedPower = squeeze(sum(Ix,2)) * dy + sidePowerLoss; plot(compensatedPower, &#39;-&#39;, &#39;LineWidth&#39;,2) line([0 numel(compensatedPower)],[max(compensatedPower) max(compensatedPower)],&#39;Color&#39;,[0 0 0],&#39;LineStyle&#39;,&#39;--&#39;) ylim([0 1.1*max(compensatedPower)]) ylabel(&#39;Integrated intensity [W]&#39;) xlabel(&#39;Main propagation axis [pixels]&#39;) set(gca,&#39;FontSize&#39;,18) title([&#39;Total power loss: &#39; num2str(max(compensatedPower)-compensatedPower(end),2) &#39;W&#39;])</code></pre> http://www.k-wave.org/forum/topic/analysis-of-the-power-loss-along-propagation-in-lossless-medium#post-5129 Wed, 08 Jul 2015 14:42:41 +0000 Anthony 5129@http://www.k-wave.org/forum/ <p>Hi Brad,</p> <p>My goal is to calculate the heat deposition in the context of HIFU treatments. </p> <p>I don't know if you remember an old topic I created about heat source term calculation. I would like to avoid the plane wave approximation, so I tried to compute -div(&lt;I&gt;), without success. Then I moved to a compromise:<br /> - Let z denote the main propagation axis (in 3D).<br /> - I want to compute how much power is lost by the acoustic wave within each xy plane. To do that, I compute \int Iz dx dy and differentiate it with respect to z to get the deposited power per meter (of course, I also compensate for the power leaving the domain through the sides).<br /> - Then, I use the plane wave approximation (or anything else) to distribute this deposited power within each xy plane. </p> <p>At least, I am sure that the energy deposited within each plane exactly corresponds to the energy that the acoustic wave has lost. Then, the exact repartition within each xy plane is approximated but it should not change much my simulated lesions. </p> <p>My problem is:<br /> when I test this method with lossless propagation, there should be no power power loss and therefore no heating. However, in the 2D example I reported, there seems to be some significant power loss along the propagation (Note: the main propagation axis is noted x instead of z in my 2D example). </p> <p>I am not 100% sure it comes from k-wave: maybe the way I reason or calculate is wrong and yet it seems simple in the 2D example...</p> <p>Thanks for taking the time to answer!<br /> Anthony </p> http://www.k-wave.org/forum/topic/analysis-of-the-power-loss-along-propagation-in-lossless-medium#post-5128 Wed, 08 Jul 2015 13:30:39 +0000 Bradley Treeby 5128@http://www.k-wave.org/forum/ <p>Ok, no problems.</p> <p>Could you walk me through exactly what you are trying to calculate? </p> http://www.k-wave.org/forum/topic/analysis-of-the-power-loss-along-propagation-in-lossless-medium#post-5123 Wed, 08 Jul 2015 10:17:01 +0000 Anthony 5123@http://www.k-wave.org/forum/ <p>Hi Brad,</p> <p>I usually do that also in my real codes. Here I simply adjusted the number of points to keep: as there were 16.66666667 points per wavelength with the default CFL, I simply kept 50 points to have 3 entire periods. I admit this is very unelegant ;-)</p> <p>Anthony </p> http://www.k-wave.org/forum/topic/analysis-of-the-power-loss-along-propagation-in-lossless-medium#post-5120 Wed, 08 Jul 2015 09:11:45 +0000 Bradley Treeby 5120@http://www.k-wave.org/forum/ <p>Hi Anthony,</p> <p>Sorry for not coming back to you sooner.</p> <p>I haven't run your code, however, if you want to calculate time average intensity, then it's important to make sure that you are using an integer number of time points per period (remember the formula for time average intensity is defined as the integral over a period). Otherwise, your time average intensity will be incorrect.</p> <p>I normally use some code like this:</p> <pre><code>% compute points per temporal period PPP = round(points_per_wavelength / CFL); % compute corresponding time spacing dt = 1/(PPP*source_freq); % create the time array using an integer number of points per period t_end = 30e-6; kgrid.t_array = 0:dt:t_end; % only record the final few periods when the field is in steady state RECORD_PERIODS = 3; sensor.record_start_index = kgrid.Nt - RECORD_PERIODS*PPP + 1;</code></pre> <p>Let me know if that helps.</p> <p>Brad. </p> http://www.k-wave.org/forum/topic/analysis-of-the-power-loss-along-propagation-in-lossless-medium#post-5101 Tue, 09 Jun 2015 10:13:33 +0000 Anthony 5101@http://www.k-wave.org/forum/ <p>Hi Brad,</p> <p>Thanks for your very convincing answer. However, I still have a problem of energy conservation. I slightly modified your code in order to show you what I mean. There is probably something wrong about my calculation but, even with alpha = 0, I calculate a power loss of 2.6 W along the beam path, which represents 30 percents of the power loss in the case of a classical absorbing material (liver) in the same situation... As you can see, I take into account the power leaving the domain through the sides. By the way, you are right, nonlinearity does not change anything.</p> <p>Best regards,<br /> Anthony</p> <pre><code>clear all close all clc %% PARAMETERS TO CHANGE medium.alpha_coeff = 0; % liver: 0.37 medium.alpha_power = 0; % liver: 1.266 % create the computational grid Nx = 128-40; Ny = 128-40; dx = 0.1e-3; dy = 0.1e-3; kgrid = makeGrid(Nx, dx, Ny, dy); % define the properties of the propagation medium medium.sound_speed = 1500; medium.density = 1000; medium.BonA = 0; % create time array kgrid.t_array = makeTime(kgrid, medium.sound_speed); % create initial pressure distribution using makeDisc R = Ny*0.45; Rap = 28/38*R; source.p_mask = zeros(Nx, Ny); xc = round(Nx/2); yc = round(Ny/2); xgrid = (1:Nx)-xc; ygrid = ((1:Ny)-yc)&#39;; rgrid = sqrt(bsxfun(@plus,xgrid.^2,ygrid.^2)); source.p_mask(logical(bsxfun(@times,(abs(rgrid-R)&lt;1/2), bsxfun(@times,ygrid&lt;0,abs(xgrid)&lt;Rap)))) = 1; source.p_mode = &#39;dirichlet&#39;; figure(1) clf imagesc(source.p_mask) axis equal % define source source_mag = 0.1e6; source_f0 = 3e6; source.p = source_mag*cos(2*pi*kgrid.t_array*source_f0); % zero pad and filter source.p = filterTimeSeries(kgrid, medium, [source.p, zeros(1, length(source.p))]); % define sensor sensor.mask = ones(Nx, Ny); sensor.record = {&#39;p&#39;, &#39;u_non_staggered&#39;}; % run the simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor,... &#39;DisplayMask&#39;, &#39;off&#39;, &#39;PlotScale&#39;, [-0.5, 0.5]*source_mag,&#39;PMLInside&#39;,false,&#39;PlotPML&#39;,false); % time shift ux = timeShift(sensor_data.ux_non_staggered, kgrid.dt); uy = timeShift(sensor_data.uy_non_staggered, kgrid.dt); % compute velocity magnitude u = sqrt(ux.^2 + uy.^2); % Intensity NpointsToKeep = 50; % in order to have an integer number of points) Ix = reshape(sensor_data.p.*ux,[Nx,Ny,size(ux,2)]); Iy = reshape(sensor_data.p.*uy,[Nx,Ny,size(uy,2)]); Ix = mean(Ix(:,:,end-NpointsToKeep+1:end),3); Iy = mean(Iy(:,:,end-NpointsToKeep+1:end),3); % Display figure(2) clf sidePowerLoss = ( cumsum(-Iy(1,:),2)&#39;+cumsum(Iy(end,:),2)&#39; ) * dx; compensatedPower = squeeze(sum(Ix,2)) * dy + sidePowerLoss; plot(compensatedPower, &#39;-&#39;, &#39;LineWidth&#39;,2) line([0 numel(compensatedPower)],[max(compensatedPower) max(compensatedPower)],&#39;Color&#39;,[0 0 0],&#39;LineStyle&#39;,&#39;--&#39;) ylim([0 1.1*max(compensatedPower)]) ylabel(&#39;Integrated intensity [W]&#39;) xlabel(&#39;Main propagation axis [pixels]&#39;) set(gca,&#39;FontSize&#39;,18) title([&#39;Total power loss: &#39; num2str(max(compensatedPower)-compensatedPower(end),2) &#39;W&#39;])</code></pre> http://www.k-wave.org/forum/topic/analysis-of-the-power-loss-along-propagation-in-lossless-medium#post-5096 Tue, 02 Jun 2015 18:25:56 +0000 Bradley Treeby 5096@http://www.k-wave.org/forum/ <p>Hi Anthony,</p> <p>I'm not sure if I follow exactly what you are trying to do? If you take a simple measure of the kinetic and potential energy, does that work out? I've written a simple script to demonstrate what I mean. This calculates the kinetic and potential energy, and the Lagrangian density. Once the source has been added, the energy is constant (there is no dissipation by the numerical scheme until the waves reach the PML), and the Lagrangian density is zero as expected. I would expect the same in nonlinear case. </p> <p>Hope that helps,</p> <p>Brad.</p> <pre><code>clear all; % create the computational grid Nx = 128; Ny = 128; dx = 0.1e-3; dy = 0.1e-3; kgrid = makeGrid(Nx, dx, Ny, dy); % define the properties of the propagation medium medium.sound_speed = 1500; medium.density = 1000; % create time array kgrid.t_array = makeTime(kgrid, medium.sound_speed); % create initial pressure distribution using makeDisc source.p_mask = zeros(Nx, Ny); source.p_mask(Nx/2, 50:60) = 1; % define source source_mag = 10e6; source.p = source_mag*toneBurst(1/kgrid.dt, 3e6, 3); % zero pad and filter source.p = filterTimeSeries(kgrid, medium, [source.p, zeros(1, length(source.p))]); % define sensor sensor.mask = ones(Nx, Ny); sensor.record = {&#39;p&#39;, &#39;u_non_staggered&#39;}; % run the simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor,... &#39;DisplayMask&#39;, &#39;off&#39;, &#39;PlotScale&#39;, [-0.5, 0.5]*source_mag); % time shift ux = timeShift(sensor_data.ux_non_staggered, kgrid.dt); uy = timeShift(sensor_data.uy_non_staggered, kgrid.dt); % compute velocity magnitude u = sqrt(ux.^2 + uy.^2); % compute Lagrangian density kinetic = sum(u.^2, 1) * 0.5 * medium.density; potential = sum(sensor_data.p.^2, 1) / (2*medium.density*medium.sound_speed.^2); Lagrangian_density = kinetic - potential; % plot change in energy over time figure; subplot(2, 1, 1); plot(kgrid.t_array, kinetic, &#39;b-&#39;); hold on; plot(kgrid.t_array, potential, &#39;r-&#39;); legend(&#39;kinetic&#39;, &#39;potential&#39;) title(&#39;Energy Components&#39;); subplot(2, 1, 2); plot(kgrid.t_array, Lagrangian_density, &#39;k-&#39;); title(&#39;Lagrangian Density&#39;);</code></pre> http://www.k-wave.org/forum/topic/analysis-of-the-power-loss-along-propagation-in-lossless-medium#post-5090 Wed, 27 May 2015 10:11:21 +0000 Anthony 5090@http://www.k-wave.org/forum/ <p>Hi,</p> <p>I am working on the simulation of HIFU treatments and I wanted to check the power which is dissipated by the numerical scheme itself (in the case of nonlinear propagation) in order to correct the deposited energy in the absorbing case.</p> <p>Let z denote the main propagation of my field (represented here for information: <a href="https://www.dropbox.com/s/089etiubtkhkkwt/XZ-field.png?dl=0)" rel="nofollow">https://www.dropbox.com/s/089etiubtkhkkwt/XZ-field.png?dl=0)</a>.<br /> I computed the integral of I_z on each xy-plane and plotted the derivative of the power with respect to z (see <a href="https://www.dropbox.com/s/lgynn4yaa3esg6o/deposited_power.png?dl=0)" rel="nofollow">https://www.dropbox.com/s/lgynn4yaa3esg6o/deposited_power.png?dl=0)</a>.</p> <p>My question is: before the focus, the energy conservation is approximately OK but downstream to it, there is some significant power loss: could you tell me why?</p> <p>A few precisions:<br /> - of course, I removed the power which leaves the domain by the sides (using I_x and I_y).<br /> - my source frequency is 3 MHz and my grid step is 40 microns, which is OK to properly propagate several harmonics.<br /> --&gt; At first, I thought that the energy loss was due to the generation of harmonics which are not supported by the grid but, even at the end of my domain, the high harmonics are not carrying that much energy (a few percents). More info about this in PS.<br /> - I set my source as a Dirichlet boundary condition. I first simulated my acoustic field with the linear equations up to 6 mm upstream to the focus and then used this as my source term for the nonlinear simulation because I couldn't simulate the whole field in the nonlinear case due to memory limitations. By the way, why do I get these power oscillations? (my medium is homogeneous).<br /> - I didn't forget to unstagger the velocity fields in time using timeShift before computing the intensity.</p> <p>Don't hesitate to ask me for more details or figures if I am not clear. I could also send you the acoustic field itself (18 GB) or only the intensity values I computed (~100MB along each direction).</p> <p>Best regards and sorry to nitpick ;-)<br /> Anthony</p> <p>PS: in order to assess the harmonic content of the field, I computed the module of the time-fft of my pressure field. Then for each voxel, I divided the module of each harmonics by the sum of the modules (taking the square gives approximately the same results). Finally, I took the mean of this proportion over each xy-plane and plotted this curve: <a href="https://www.dropbox.com/s/9ee9vfy4vd6w8lo/harmonic_content.png?dl=0" rel="nofollow">https://www.dropbox.com/s/9ee9vfy4vd6w8lo/harmonic_content.png?dl=0</a><br /> Little detail: in each xy-plane, I only considered the pixels where the rms pressure was greater than -6dB of the max in this plane i.e. only the pixels "inside the beam".<br /> --&gt; As you can see, the high harmonics are not very important... (well, I admit that my metric is not wonderful, if you want me to measure this using another formula, feel free to ask). </p>