k-Wave

1. user5137
   

   Member
   Joined: Jul '13
   Posts: 3

   Hi there.

   I'm getting started with trying to simulate the propagation of a simple acoustic wave through water. I am in an exploration stage at the moment and the majority of my code is an adaptation of the given examples. I was wondering if you could help me understand the following.

   I have put a line source (with absorbing boundaries switched off) at an angle of pi/4 wrt the grid axes. The grid (which I copied from your example) is as follows

   % create the computational grid
   Nx = 256;           % number of grid points in the x (row) direction
   Ny = 256;           % number of grid points in the y (column) direction
   dx = 500e-6/Nx;      	% grid point spacing in the x direction [m]

   The 'video' playback looks reasonable - I have a screenshot of it here.

   

   When I increase the number of points in the grid the behaviour of the acoustic field becomes less believable. At Nx = Ny > 470 the propagating wave starts off looking more and more granular. By my estimation increasing Nx,Ny should increase the number of points per wavelength. An example where this granularity is quite clear is

   Nx = 512;           % number of grid points in the x (row) direction
   Ny = 512;           % number of grid points in the y (column) direction
   dx = 500e-6/Nx;      	% grid point spacing in the x direction [m]

   The playback changes to something like

   

   I expect this issue is borne out of my lack of experience using this toolbox. My question is whether this looks like a familiar error, eg spatial / temporal dispersion?

   The source is a single cycle sinusoid of frequency 5MHz.

   Many thanks,

   Posted 11 years ago #

2. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   I certainly agree it looks a bit suspicious. If you are running a linear simulation in a lossless and homogeneous medium, there shouldn't be any numerical dispersion as the numerical scheme is exact. The fact that it starts to happen after a particular grid size is puzzling. What version of MATLAB are you using, and do have the 'DataCast' parameter set?

   Without seeing your code it's difficult to guess exactly what the problem might be. If you post it here I can check whether I see the same error or spot something wrong with your script.

   Brad.

   Posted 11 years ago #

3. user5137
   

   Member
   Joined: Jul '13
   Posts: 3

   Hi Brad,

   Thank you for replying so promptly. It's wonderful to see a toolbox so actively developed.

   In answer to your questions, I am using Matlab 2011b (v7.13). I have not set the DataCast parameter - I see in the kspacefirstorder2D function that it is off by default. Also, I know you didn't ask, but I am using kwave v1.0 (13 Nov 2012 release).

   I have modified my script to make it a little easier to parse, so the simulation will look slightly different to the screenshots posted above, however the effect of changing the grid remains the same.

   clear all;
   close all;
   close all hidden;
   
   % create the computational grid
   Nx = 512;           % number of grid points in the x (row) direction
   Ny = 512;           % number of grid points in the y (column) direction
   dx = 500e-6/Nx;      	% grid point spacing in the x direction [m]
   dy = dx;            % grid point spacing in the y direction [m]
   kgrid = makeGrid(Nx, dx, Ny, dy);
   
   % define the properties of the propagation medium
   medium.sound_speed = 1480;  % [m/s]
   
   %define source.
   source_frequency = 5e6;
   source_mag = 1;
   
   % create the time array
   Nt = 250;  % time steps
   dt = 1e-9; % size of time step
   kgrid.t_array = (1:Nt)*dt;
   
   % remove the detector
   sensor.mask = zeros(Nx, Ny);
   
   % angle of line source wrt grid
   angle_of_source = pi/4;
   
   length_of_line = Nx;
   startpoint = [Nx/2 Ny/2];
   source_line = makeLine(Nx, Ny, startpoint, angle_of_source, length_of_line);
   source_positions = find(source_line >= 1);
   
   % define a time varying sinusoidal source
   source.p = source_mag*sin(2*pi*source_frequency*kgrid.t_array);
   
   % make the source of only 1 cycle
   time_for_one_cycle = 1 / source_frequency;
   x = (find (kgrid.t_array > time_for_one_cycle));
   source.p(x(1):length(kgrid.t_array)) = 0;
   
   % filter the source to remove high frequencies not supported by the grid
    source.p = filterTimeSeries(kgrid, medium, source.p);
   
   % select a source
   source.p_mask = zeros(Nx, Ny);
   source.p_mask(source_positions) = 1;
   
   % create a display mask to display the transducer
   display_mask = source.p_mask + sensor.mask;
   
   % run the simulation with PML turned off
   input_args = {'DisplayMask', display_mask, 'PlotScale', [-0.5 0.5], 'PMLAlpha', [0 0]};
   
   sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

   I do look forward to resolving this so I can continue to explore the toolbox and thank you in advance for your time :)

   Posted 11 years ago #

4. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   The problem is related to the size of the time step. If you set this as dt = CFL * dx / c, where CFL <= 1, this should solve the problem. The CFL (Courant-Friedrichs-Lewy) number can be interpreted as the distance a wave can travel in a single time step relative to the spatial grid spacing. If you use the makeTime function, the CFL number is set to 0.3 by default.

   Brad.

   Posted 11 years ago #

5. user5137
   

   Member
   Joined: Jul '13
   Posts: 3

   Hi Brad,

   Thank you for taking the time to look through this code and establish the source of the error. As I suspected, it was borne out of my ignorance rather than the fault of the toolbox. Instead of setting the duration of the timestep myself, I've now taken to using the makeTime function and the simulation executes as one would expect.

   A couple of suggestions, while I am here.

   1. Have you considered using stackoverflow rather than a custom forum for these types of questions? I've seen quite a few open source projects moving to this platform and there is a strong community presence.

   2. I've seen others mention using a code repository (e.g. github or bitbucket) so I'd just add my vote as it would make contributing to the project much easier :)

   Once again, thank you for your input - I'm looking forward to further exploring this toolbox.

   Posted 11 years ago #

RSS feed for this topic

Reply

You must log in to post.