k-Wave

1. zak
   

   Member
   Joined: Oct '12
   Posts: 3

   Hi, Dr. Cox and Dr. Treeby,

   I am trying to model the photoacoustic wave propagation in heterogeneous media, where the heterogeneity is only a small part but very strong (in terms of both sound speed and mass density) compared to rest of the media, e.g. bone or skull in soft tissue.

   In heterogeneous media, this paper says that the k-space method is expected to be unconditionally stable if the reference sound speed c_ref in the k-space operator is chosen to be the maximum sound speed in the media (c_max), and this is also the default choice of the k-Wave.

   However, in your recent paper of nonlinear wave propagation, it says the phase errors will be introduced in the simulation if the c_ref doesn't match the local sound speed. The k-space operator will over compensate for the actual phase error introduced by the finite difference time step if c_ref is much larger than the local sound speed.

   So I am not sure how to choose the proper c_ref in my case. On one hand, it seems I should choose c_max as c_ref given what the first paper says; on the other hand, this choice is problematic according to your paper since the sound speed in most part of the media (c_0) is much smaller than c_max; it indicates that I should set c_ref = c_0.

   Should I choose c_ref be to c_max or c_0? Or is there any other choice?

   Thank you very much!

   Posted 12 years ago #

2. bencox
   

   Developer
   Joined: Mar '10
   Posts: 292

   Hi Zak,

   Good question. The k-space model is certainly unconditionally stable if you choose c_ref > max(c), in other words the solution won't 'blow up'. Unfortunately the phase errors can become very large, so while the solution will be stable it might also be completely wrong.

   The alternative is to choose c_ref<min(c) in which case the phase errors are bounded (they will never be worse than if the k-space correction term sinc(ckdt/2) were just left out) but now the model could become unstable and it is necessary for the timesteps to satisfy the stability criterion:

   dt < (2/(c_ref * max(k))) * asin(c_ref/max(c))

   Hope that helps you choose,

   Ben

   Posted 12 years ago #

3. zak
   

   Member
   Joined: Oct '12
   Posts: 3

   Hi Dr. Cox

   Thanks for your help.

   From your suggestion, I can set c_ref = c_min, then estimate the upper bound of the time step. As such, the phase errors is not worse than the k-space adjustment being left out. But I wonder if we can just neglect the k-space adjustment when doing simulations in heterogeneous media so there is no c_ref to be determined in the first place.

   Zak

   Posted 12 years ago #

4. bencox
   

   Developer
   Joined: Mar '10
   Posts: 292

   Hi Zak,

   You could neglect the k-space adjustment, but there would be no advantage to doing so (other than no having to make a decision about c_ref!) There would be disadvantages though: there would be no phase error correction at all, and you would have to satisfy a stability requirement which is always going to be more stringent:

   dt < 2/(k_max * c_max)

   Ben

   Posted 12 years ago #

5. zak
   

   Member
   Joined: Oct '12
   Posts: 3

   Hi Dr. Cox,

   I tried your suggestion by setting c_ref = c_min and dt < (2/(c_ref * max(k))) * asin(c_ref/max(c)), but the code prompted a warning saying "time step may be too large for a stable simulation" and the simulated results blowed up.

   Here are my codes:

   Nx=512;
   dx=0.5e-3;
   kgrid=makeGrid(Nx,dx,Nx,dx);

   medium.sound_speed=1500*ones(Nx,Nx);
   medium.sound_speed(:,200:300)=3000;
   medium.density=1000*ones(Nx,Nx);
   medium.density(:,200:300)=2000;
   medium.sound_speed_ref=1500;

   (2/(medium.sound_speed_ref*kgrid.k_max))*asin(medium.sound_speed_ref/max(medium.sound_speed(:)))
   ans =
   1.111111111111111e-07

   kgrid.t_array=(0:2000)*1e-7;

   source.p0=makeDisc(Nx,Ny,Nx/2,100,10);

   sensor.mask=zeros(Nx,Nx);
   sensor.mask(Nx/2,400)=1;

   input_args = {'Smooth', true, 'PMLInside', true, 'PMLSize', 10, 'DataCast', 'gpuArray-single', 'DataRecast', true,'PlotPML', false, 'PlotSim', false};

   p=kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

   However, if I used the default dt=5e-8 which is calculated by the criteria dt = CFL * dx_min / c_max, then there's no problem.

   Posted 11 years ago #

6. bencox
   

   Developer
   Joined: Mar '10
   Posts: 292

   Hi Zak,

   Thanks for this; it's raised an interesting point. If you replace kgrid.k_max by max(kgrid.k) it should work.

   The reason that these are not the same is because we have defined kgrid.k_max as the maximum wavenumber supported by the grid in every direction, in other words it is defied as kgrid.k_max = min([kgrid.kx_max, kgrid.ky_max]). However, there will be components whose wavenumbers in both directions are maximum, and so the absolute value sqrt(kgrid.kx_max.^2 + kgrid.ky_max.^2) will be bigger than kgrid.k_max.

   Ben

   Posted 11 years ago #

RSS feed for this topic

Reply

You must log in to post.