drop end of domain « k-Wave User Forum

k-Wave User Forum » General Acoustics » Acoustic Simulations

drop end of domain

(6 posts) (2 voices)

Started 8 years ago by mathiasb
Latest reply from Bradley Treeby

mathiasb

Member
Joined: Apr '16
Posts: 3

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

https://dl.dropboxusercontent.com/u/581753/1wave.png

Posted 8 years ago #
Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,643

Hi Mathias,

It certainly sounds a bit like a PML effect. What happens if you turn the PML off by setting 'PMLAlpha', 0? 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.

Posted 8 years ago #

Member
Joined: Apr '16
Posts: 3

Hi Brad

thanks for the tips.
Here I plot with PMLAlpha set to 0,
https://dl.dropboxusercontent.com/u/581753/onewave_PMLAlphastudy.png
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

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)

Posted 8 years ago #

Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,643
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 makeTime line to:
```
t_end = 15e-3;
[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed, [], t_end);
```
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:
```
CFL = 0.14;
t_end = 15e-3;
[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed, CFL, t_end);
```
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 here 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.
Posted 8 years ago #
mathiasb

Member
Joined: Apr '16
Posts: 3
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:
```
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);
```
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
Posted 8 years ago #
Bradley Treeby

Developer
Joined: Mar '10
Posts: 1,643
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 spect function will not be exact, and will gradually increase as you increase the total time (eventually converging for very large t_end).

To avoid this, you could set the time step dt using something like the following:
```
% 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);
```
Then only record the wave field for the last few periods in steady state
```
% only record the last few periods in steady state
sensor.record_start_index = Nt - NUM_PERIODS_TO_RECORD*PPP + 1
```
This way the recorded time series will be exactly periodic.

Brad.
Posted 8 years ago #

RSS feed for this topic

Reply

You must log in to post.

k-Wave

A MATLAB toolbox for the time-domain simulation of acoustic wave fields

drop end of domain

Reply

A MATLAB toolbox for the time-domain
simulation of acoustic wave fields