http://www.k-wave.org/forum/topic/spurious-wave-generation-in-absorbing-media Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 23:09:20 +0000 http://bbpress.org/?v=1.0.2 q http://www.k-wave.org/forum/search.php http://www.k-wave.org/forum/topic/spurious-wave-generation-in-absorbing-media#post-4984 Mon, 09 Feb 2015 15:04:54 +0000 DNF 4984@http://www.k-wave.org/forum/ <p>Thank you, Bradley. At this point, it isn't terribly urgent to get the files quickly, but I am very interested in the content. </p> http://www.k-wave.org/forum/topic/spurious-wave-generation-in-absorbing-media#post-4975 Wed, 04 Feb 2015 23:46:05 +0000 Bradley Treeby 4975@http://www.k-wave.org/forum/ <p>Hi DNF,</p> <p>If you capture the time-varying pressure over a 2D plane, and then enforce this as a Dirichlet boundary condition in the next simulation, this is sufficient (you don't also need particle velocity). However, the problem is that if you are enforcing a boundary condition, any waves that are reflected from heterogeneities will again be reflected from the plane over which you are enforcing the boundary condition (with an amplitude dependent on the exact value of pressure you happen to be enforcing at that time). This means your simulations results will not match the non-split case.</p> <p>I'm currently away for a meeting, but when I'm back at my desk, I'll tidy up the codes and email them to you.</p> <p>Brad. </p> http://www.k-wave.org/forum/topic/spurious-wave-generation-in-absorbing-media#post-4963 Wed, 28 Jan 2015 09:38:27 +0000 DNF 4963@http://www.k-wave.org/forum/ <p>Thank you, Bradley.</p> <p>I am interested both in the elastic code, and in the SimSonic wrapper you mentioned previously.</p> <p>As for splitting the simulation, I would have thought that recording the pressure field a slightly before the first interface, then re-transmitting it into the layered structure (with much higher sampling), and then recording the pressure again, I would be able to capture multiple reflections. Anyway, is capturing the pressure sufficient? Are not the particle velocities required? </p> http://www.k-wave.org/forum/topic/spurious-wave-generation-in-absorbing-media#post-4960 Fri, 23 Jan 2015 17:46:44 +0000 Bradley Treeby 4960@http://www.k-wave.org/forum/ <p>Hi DNF,</p> <p>The function <code>kspaceElastic2D</code> is a k-space code for elastic media with power law absorption, as discussed in <a href="http://www.homepages.ucl.ac.uk/~rmapbtr/papers/JOURN_23_2014_Treeby_JASA_PowerLawAbsorptionElastic.pdf">this paper</a>. This model is still under development, so is not in the current release code. At the moment, there is only a 2D version and no PML is implemented. However, I'd be happy to send you a copy to try out if that's helpful.</p> <p>If you are only interested in the first transmitted wave, then dividing your simulation into chunks in this way will work. If you record the pressure field over a plane in the first simulation, you can enforce this as a Dirichlet boundary condition in the next simulation. However, keep in mind that you will not be able to model multiple reflections, etc.</p> <p>If you are using the C++ code, there is a checkpoint-restart facility that can be easily used (see the <a href="http://www.k-wave.org/doxygen/index.html">doxygen documentation</a> for more details). Unfortunately, the same facility doesn't exist in the MATLAB code at the moment</p> <p>Hope that helps,</p> <p>Brad. </p> http://www.k-wave.org/forum/topic/spurious-wave-generation-in-absorbing-media#post-4956 Fri, 23 Jan 2015 09:01:06 +0000 DNF 4956@http://www.k-wave.org/forum/ <p>Thanks for looking into this.</p> <p>Now, <code>pstdElastic2D</code> is not not a k-space method, as you say, but I recently came over a reference to a function <code>kspaceElastic2D</code>, inside the script <code>kspaceFirstOrder_inputChecking</code> (line 45). Is that code which is under development?</p> <p>I am considering splitting up the calculation into several chunks. Since I have a layered structure, I thought it might be possible to propagate in one homogeneous layer, save the state, propagate across a boundary with higher sampling, save the state, etc. Might this be a viable approach? </p> <p>Is it easy or even possible to save the state of the simulation? This would also enable me to pause and resume the simulation, or interrupt/fail while keeping the results up to that point, which is something I miss a bit. </p> http://www.k-wave.org/forum/topic/spurious-wave-generation-in-absorbing-media#post-4945 Thu, 22 Jan 2015 16:28:49 +0000 Bradley Treeby 4945@http://www.k-wave.org/forum/ <p>Hi DNF,</p> <p>Unfortunately, I haven't managed to get to the bottom of the behaviour you have observed. It seems to be related the fractional absorption operator in the case <code>medium.alpha_power</code> is not equal to 2, and <code>medium.alpha_coeff</code> is high. Using <code>medium.alpha_mode = &#39;no_absorption&#39;</code> causes the effect to disappear, so it doesn't appear to be related to the dispersion term. It's possible it's related the particle velocity used in the absorption term being temporally offset by dt/2 from the acoustic density to which it's added. I will keep digging.</p> <p>In the elastic case, my guess is you're just seeing classic instability. Unfortunately, we have not derived a closed form expression for stability in the absorbing elastic case using the PSTD model. However, remember the elastic code <code>pstdElastic2D</code> does not use a k-space corrected time-stepping scheme, so is not unconditionally stable. You will need to use a smaller time step than in the associated fluid case, both to achieve stability, and to avoid numerical dispersion from the time-stepping scheme.</p> <p>One point to note, if you are interested in modelling very strong impedance discontinuities, spectral methods might not be your best choice. In <a href="http://www.k-wave.org/papers/2014-Robertson-IEEEIUS.pdf">this paper</a> there is a plot of the accuracy in the reflection coefficient versus the strength of the discontinuity. Roughly speaking, as the code stands, for contrasts stronger than 4:1, finite difference methods will give you more accurate reflection coefficients. If your domain size isn't very large (so the accumulation of dispersion error isn't a problem), you might be better served using a tool like <a href="http://www.simsonic.fr">SimSonic</a>.</p> <p>The only other thing to mention is that when viewing the simulations with a log scale, small-scale behaviour becomes magnified. Given the perfectly matched layer will only be accurate to a few decimal places at best, in some cases the spurious waves you have observed might not be a problem.</p> <p>Hope that helps,</p> <p>Brad. </p> http://www.k-wave.org/forum/topic/spurious-wave-generation-in-absorbing-media#post-4938 Fri, 16 Jan 2015 09:32:59 +0000 DNF 4938@http://www.k-wave.org/forum/ <p>Yes, I am aware of the basic reason why this happens. But absorption and the choice of alpha_power can increase the effect dramatically. The effect is very noticeable. Even using 30-40 points per wavelength, it is still clearly present. I am in a situation where a 1D simulation can take up to 20 hours on a new and powerful pc.</p> <p>And in the elastic simulations, as I said, the effect grows without bounds, eventually completely drowning out the real solution. </p> http://www.k-wave.org/forum/topic/spurious-wave-generation-in-absorbing-media#post-4937 Thu, 15 Jan 2015 16:09:15 +0000 ouelletn 4937@http://www.k-wave.org/forum/ <p>Run the 1D code for a very short time and stop it. If you zoom into the actual wave you should notice that the wave is undergoing a ringing effect that extends into the second media. This ringing effect is inherent of the high discontinuities in your simulation and a limitation of expressing the wave as a Fourier series. </p> <p>The only way to minimize this is to increase the number of grid points within the simulations.Changing the frequency distribution of the impulse may also help reduce this effect. You should consider mostly if this is necessary for your simulation, as the amplitude of this effect is fairly low and may not influence any statements made about the effect of the layers.<br /> -Andrew </p> http://www.k-wave.org/forum/topic/spurious-wave-generation-in-absorbing-media#post-4934 Wed, 14 Jan 2015 11:06:52 +0000 DNF 4934@http://www.k-wave.org/forum/ <p>This effect is also visible in the 2D version of the code.</p> <p>Worse, when using <code>pstdElastic2D.m</code> with strong impedance contrasts, spurious waves are generated that in some cases grow exponentially and without bound, to the point that they reach a value of <code>inf</code>, causing the visualization to give an error. Both spurious shear and pressure waves are generated, and grow wildly.</p> <p>This happens even when I use no absorption at all. </p> http://www.k-wave.org/forum/topic/spurious-wave-generation-in-absorbing-media#post-4925 Mon, 05 Jan 2015 10:42:52 +0000 DNF 4925@http://www.k-wave.org/forum/ <p>I am using k-Wave to simulate wave propagation in strongly contrasting media. A strange effect shows up during propagation when there is strong absorption. The attached figure illustrates the effect: Shortly after starting the excitation, and before the pulse arrives, a wave is generated inside the neighbouring layer.</p> <p>The magnitude of this effect is apparently influenced by several factors:<br /> * The level of absorption — more absorption _increases_ the effect, contrary to my intuition.<br /> * The distance from the source to the medium interface.<br /> * The number of layers in the total medium (four layers is worse than three).<br /> * When the wave reaches the PML it gets even worse.<br /> * Using a power law exponent of 2 seems to remove the problem (this was pointed out to me by Dr. Treeby.)</p> <p>I can reduce the effect by increasing the temporal and spatial sampling rates, but some times even using 40 points per wavelength and a timestep less than a 1ns is not good enough for a transmit frequency of 500kHz.</p> <p>I have enclosed an example script that illustrates the effect. It seems that the effect gets worse at some ’sweet spot’ of combined impedances and layer thicknesses, which I cannot quite re-create.</p> <p>My questions:<br /> * Do you see why this is happening? Could there be a bug in the absorption implementation, or is this more or less expected behaviour?<br /> * What can I do to reduce this effect? For certain geometries, cranking up the sampling rates leads to 20h+ long 1D simulation.</p> <p>Here's an example script:</p> <pre><code>%% % This example shows a strange effect where a spurious wave is generated inside a % medium layer before the excitation wave arrives in that layer. % % This happens when there is strong absorption inside one of the layers (e.g. 10dB/cm/MHz), and the % effect is smaller when the absorption is reduced (e.g. to 1dB/cm/MHz). Tweaking the % geometry in various ways also seems to change the effect, having four layers seems % to be worse than having three layers, and the thickness of the various layers also % make a difference. % % The effect can be reduced by reducing the timestep, but sometimes to extreme levels. %% Geometry padding = 5e-3; Lx = [6,2,10,20]*1e-3; %% Material properties c = [600, 3000, 6000, 3000]; rho = [200, 2000, 8000, 2000]; attn = [0.2, 10, 0.05, 0.1]; % attn = [0.2, 1.0, 0.05, 0.1]; % reducing the absorption reduces the effect attn_power = 1.0; % changing it to 2.0 removes the problem useattn = true; %% Transmit pulse Npulse = 2; fc = 0.5e6; %% Computational grid ppw = 20; % points per wavelength cfl = 0.1; % Courant-Friedrich-Lewy % grid size dx = min(c)/(ppw*fc); Nx = ceil(sum([Lx,padding])/dx); kgrid = makeGrid(Nx,dx); % index positions of interfaces ind = round(cumsum([padding,Lx(1:end-1)])/dx); medium = struct(); % sound velocity medium.sound_speed = c(1) + zeros(Nx,1); for jj = 2:length(ind), medium.sound_speed(ind(jj):end) = c(jj); end % density medium.density = rho(1) + zeros(Nx,1); for jj = 2:length(ind), medium.density(ind(jj):end) = rho(2); end % attenuation if useattn, medium.alpha_coeff = attn(1) + zeros(Nx,1); for jj = 2:length(ind), medium.alpha_coeff(ind(jj):end) = attn(jj); end medium.alpha_power = attn_power; medium.alpha_mode = &#39;no_dispersion&#39;; end [kgrid.t_array,dt] = makeTime(kgrid,medium.sound_speed,cfl); %% Transducer Amp = 1e-6; source = struct(); source.u_mask = zeros(Nx,1); source.u_mask(ind(1)) = 1; source.ux = -Amp*toneBurst(1/kgrid.dt,fc,Npulse,&#39;Envelope&#39;,&#39;Gaussian&#39;); sensor = struct(); sensorind = ind(1):10:Nx; sensor.mask = zeros(Nx,1); sensor.mask(sensorind) = 1; sensor.record = {&#39;p&#39;}; % display display_mask = zeros(Nx,1); display_mask(ind) = 1; %% 1D sim args = {&#39;DataCast&#39;,&#39;double&#39;, ... &#39;LogScale&#39;,true, ... &#39;DisplayMask&#39;,display_mask, ... &#39;PMLSize&#39;,ppw*2, ... % 2 wavelengths &#39;PMLInside&#39;,false, ... &#39;PlotPML&#39;,false, ... &#39;PlotSim&#39;,true}; data = kspaceFirstOrder1D(kgrid,medium,source,sensor,args{:});</code></pre>