k-Wave User Forum » Topic: drop end of domain

k-Wave User Forum » Topic: drop end of domain http://www.k-wave.org/forum/topic/drop-end-of-domain Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 01:16:54 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php Bradley Treeby on "drop end of domain" http://www.k-wave.org/forum/topic/drop-end-of-domain#post-5636 Thu, 11 Aug 2016 14:51:33 +0000 Bradley Treeby 5636@http://www.k-wave.org/forum/ Hi Mathias, That doesn't sound like expected behaviour. Is it related to your processing script? For example, if you don't record an integer number of points per period, the amplitude you extract using the <code>spect</code> function will not be exact, and will gradually increase as you increase the total time (eventually converging for very large <code>t_end</code>). To avoid this, you could set the time step <code>dt</code> using something like the following: <pre><code>% define time step using an integer number of PPP PPW = SOUND_SPEED/(SOURCE_FREQ*dx); PPP = floor(PPW/CFL); dt = 1/(SOURCE_FREQ*PPP); % set the number of time steps Nt = round(t_end / dt);</code></pre> Then only record the wave field for the last few periods in steady state <pre><code>% only record the last few periods in steady state sensor.record_start_index = Nt - NUM_PERIODS_TO_RECORD*PPP + 1</code></pre> This way the recorded time series will be exactly periodic. Brad. mathiasb on "drop end of domain" http://www.k-wave.org/forum/topic/drop-end-of-domain#post-5624 Mon, 08 Aug 2016 14:43:01 +0000 mathiasb 5624@http://www.k-wave.org/forum/ Hi Bradley, thanks for your reply, which helped me to better understand and mitigate the issue. So now I am using something like this for the time vector: <pre><code>t_end = (Nx + 2*PML_X_SIZE)*dx/medium.sound_speed*1.3; %increased by 30%. [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed, CLF, t_end);</code></pre> That is, I increase the time for the wave to go from start to end of domain of 30% (roughly what I think yourself used in your reply) so that I should cover the whole simulation. However, I noticed that longer time array will keep giving me higher pressure levels: the more I increase t_array, the higher the calculated pressure is. I do not see convergence to an optimal length (tried until I increased t_array of 500%). Any suggestion on how to interpret that? Thanks Mathias Bradley Treeby on "drop end of domain" http://www.k-wave.org/forum/topic/drop-end-of-domain#post-5558 Mon, 04 Jul 2016 14:51:32 +0000 Bradley Treeby 5558@http://www.k-wave.org/forum/ Hi Mathias, The drop off at the end of the domain is because the simulation time isn't long enough for the waves to reach steady state at the end of the domain. If you change the <code>makeTime</code> line to: <pre><code>t_end = 15e-3; [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed, [], t_end);</code></pre> This should sort everything out. Note, in your code, the CFL number that you define is not passed to make time. If you want to set the CFL, you should change this to: <pre><code>CFL = 0.14; t_end = 15e-3; [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed, CFL, t_end);</code></pre> In addition, your source frequency is only at ~2.5 points per wavelength. I would recommend using at least 3 PPW for the PML to work correctly (see <a href="http://www.k-wave.org/papers/2014-Robertson-IEEEIUS.pdf">here</a> for example). The ripples you see are just the on-axis response for a piston transducer in 2D (due to the interaction of the planar wave front and the edge waves). They will disappear if you make your source a single point, although I expect you may need to expand your grid laterally to make room for a larger PML and avoid shallow angles of incidence. Hope that helps, Brad. mathiasb on "drop end of domain" http://www.k-wave.org/forum/topic/drop-end-of-domain#post-5543 Thu, 23 Jun 2016 09:29:51 +0000 mathiasb 5543@http://www.k-wave.org/forum/ Hi Brad thanks for the tips. Here I plot with PMLAlpha set to 0, <a href="https://dl.dropboxusercontent.com/u/581753/onewave_PMLAlphastudy.png" rel="nofollow">https://dl.dropboxusercontent.com/u/581753/onewave_PMLAlphastudy.png</a> It looks PML effect, doesn't it? How could I avoid it or at least mitigate it? I am struggling to get my simulation run as further away possible in terms of x distance from the source and with larger domains that drop effect already starts at ca half of the x length (the code below is an example simulation, my actual simulations are a bit more sophisticated generating nonlinear effects). My dt is in the order of microseconds, as seen in the text in the figure; number of time points is 6161. Thanks Mathias <pre><code>clear all; % ========================================================================= % grid % ========================================================================= % long_grid flag to select between a computational with a long x-axis (1) % and a short x-axis (0). PML_X_SIZE = 100; % [grid points] PML_Y_SIZE = 10; % [grid points] Nx = 2048 - 2*PML_X_SIZE; % number of grid points in the x (row) %10368 Ny = 32 - 2*PML_Y_SIZE; % number of grid points in y (column) dx = 2e-3; % grid point spacing in x [m] dy = 2e-3; % grid point spacing in y [m] kgrid = makeGrid(Nx, dx, Ny, dy); % ========================================================================= % define the properties of the propagation medium % ========================================================================= medium.sound_speed = 343; % [m/s] medium.alpha_coeff = 0.01*2.4e-7; % [dB/(MHz^y cm)] medium.alpha_power = 1.44; medium.density = 1.25; % [kg/m3] medium.BonA = 0.4; % ========================================================================= % create the time array % ========================================================================= CLF = 0.14; [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed); % ========================================================================= % select geometry of the source: line % ========================================================================= source.p_mask = zeros(Nx,Ny); source.p_mask(1,Ny/2-5:Ny/2+5) = 1; num_elements = sum(sum(source.p_mask)); % ========================================================================= % define input signal to transducers: sin() % ========================================================================= source_mag1 = 0.01; % [Pa] source_freq1 = 70e3; % [Hz] signal1 = source_mag1*sin(2*pi*source_freq1*kgrid.t_array); input_signal = signal1; source.p = input_signal; % ========================================================================= % options for solver % ========================================================================= display_mask = source.p_mask; % sensor mask sensor.mask = ones(Nx,Ny); % assign the input options input_args = {'DataCast','double','DisplayMask', source.p_mask, ... 'PMLSize', [PML_X_SIZE, PML_Y_SIZE],... 'PMLInside', false, 'PlotPML', true,... 'PlotScale', [-2*source_mag1, 2*source_mag1]}; % ========================================================================= % solver % ========================================================================= sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); %% ======================================================================== % AXIAL PRESSURE, FFT of discrete frequencies along central array axis % ========================================================================= sensor_data = reshape(sensor_data, [Nx, Ny, kgrid.Nt]); sensor_data_axial=sensor_data(:,end/2,:); nfft = 2^nextpow2(1/dt/1e3); [freq, amp_spect] = spect(sensor_data_axial, 1/dt, 'Dim', 3,'FFTLength',nfft); [f1_value, f1_index] = findClosest(freq, source_freq1); % extract the amplitude at the source frequency and store beam_pattern_f1 = amp_spect(:, :, f1_index); figure semilogx((kgrid.x_vec - min(kgrid.x_vec(:)))*1e3,20*log10(beam_pattern_f1./20e-6),'b') ylabel('Pressure magnitude [dB re 20e-6 Pa]'); xlabel('x-position on axis [mm]'); titles=... sprintf... (['f1=%d kHz, mag1=%g dB (%g Pa),'... 'Nx=%d, Ny=%d, \n CLF=%g,'... ' dx=%g mm, dy=%g mm, dt=%g s, PMLx=%d, PMLy=%d.']... ,source_freq1/1e3,round(20*log10(source_mag1/2e-5)),... source_mag1,Nx,Ny,CLF,dx*1e3,dy*1e3,... dt,PML_X_SIZE,PML_Y_SIZE); title(titles)</code></pre> Bradley Treeby on "drop end of domain" http://www.k-wave.org/forum/topic/drop-end-of-domain#post-5541 Wed, 22 Jun 2016 21:14:42 +0000 Bradley Treeby 5541@http://www.k-wave.org/forum/ Hi Mathias, It certainly sounds a bit like a PML effect. What happens if you turn the PML off by setting <code>'PMLAlpha', 0</code>? What does the display look like when the simulation is running? One other thought, is the number of time steps in your simulation large enough for the wave to reach that part of the domain? It's hard to say about the oscillations without seeing your code or knowing more about the source and medium. Brad. mathiasb on "drop end of domain" http://www.k-wave.org/forum/topic/drop-end-of-domain#post-5535 Tue, 21 Jun 2016 16:03:47 +0000 mathiasb 5535@http://www.k-wave.org/forum/ Hello I am simulating various acoustic linear and nonlinear phenomena in 2D. In various simulations, the length of my domain may vary between 128 pixels to 10368, the width between 12 to 512 pixels, with grid point spacing of typically 1 or 2 mm; my signals consist of typically 1 or 2 frequencies between 40 kHz and 100 kHz with magnitude from 50 dB to 150 dB. In all tested cases, I encounter an effect which I do not understand. As seen in the linked plot below, at the end of the domain the amplitude of my wave decays. The effect happens regardless of domain length, PML size, magnitude. I do not think it is related to physics (e.g. actual wave decay) because if the domain is, e.g., increased while keeping everything else the same, then the amplitude drop is simply shifted further away from the source. Therefore I think it is due to computational/numerical effect, but I do not understand why; I thought about reflections at the end of the domain resulting in destructive interference but I tried with PML sizes between 10 and 1000 pixels and not much changes besides resulting in shorter domains (I keep PML inside false). Also, I do not understand the oscillations in amplitude of the wave, e.g. see the dip just before 1000 mm. Any hint? Thank you, Mathias <a href="https://dl.dropboxusercontent.com/u/581753/1wave.png" rel="nofollow">https://dl.dropboxusercontent.com/u/581753/1wave.png</a>