k-Wave

1. jkiv
   

   Member
   Joined: Sep '10
   Posts: 7

   Good day,

   I'm a student using k-Wave to simulate propagation for high intensity focused ultrasound (HIFU). I'm currently having some issues dealing with the PML.

   I've been defining my heterogeneous material without thought of the PML. That is, if the top half of my system is material A and the bottom is material B, I make all of the top material A and all of the bottom material B including the area where the PML lies. However, if a material with a high alpha_coeff (eg. 2.43 dB/MHz-cm) lies within the PML we get unusual artifacts. Here's an example, where PMLSize is 20 pixels, PMLAlpha is 4 Np/(rad/s-m), and alpha_power is 1.15.

   

   
   (Note: the image above is a plot of the log10 of the pressure)

   If I set the medium's alpha_coeff in the PML to a small value (eg. 0.2 dB/MHz-cm) , the effect is diminished. Trying an alpha_coeff = 0 in the PML leaves only a minimal effect:

   

   
   (Note: the image above is a plot of the log10 of the pressure)

   This leads me to my first question: why might this be happening? Is this a bug, a consequence of my setup, or a consequence of the methods being used?

   My second question is more general that the documentation does not address directly. How does one set up their PML optimally? Given my material, my source frequency, etc, how might it be possible to automate the setup of the PML properties such that a sufficient amount of the propagation is absorbed?

   Thanks in advance,
   Jon

   Posted 13 years ago #

2. bencox
   

   Developer
   Joined: Mar '10
   Posts: 285

   Hi Jon,

   Thanks for using k-Wave. Do you get the same thing if you set the PML to be outside the domain of interest, ie. setting 'PMLInside' to false when calling kspaceFirstOrder2D?

   Just so I can replicate your example, what sort of source/initial condition are you using?

   With regard to the optimal PML conditions, we've not done a systematic study of the ideal PML settings but have set defaults that we found worked for most cases. k-Wave ramps the absorption within the PML with the curve of the form ((x-x0)/(xmax-x0))^4, where x0 is the inside edge of the PML and xmax the outer edge. This is suggested by Tabei, M., T. D. Mast, et al., J. Acoust. Soc. Am., 111(1), 53-63, (2002), the paper which the earliest versions of k-Wave were derived from. Designing PMLs that guarantee a certain level of reduction for a given incident wave sounds like an interesting project, but not one we've done yet, I'm afraid.

   Ben

   PS. For newcomers reading this, see http://www.k-wave.org/documentation/example_na_controlling_the_pml.php for an intro to the PML.

   Posted 13 years ago #

3. jkiv
   

   Member
   Joined: Sep '10
   Posts: 7

   Hey Ben,

   I tried PMLInside to find it fails:

   Running k-space simulation...
   dt: 7.5ns, t_end: 41.97us, time steps: 5597
   input grid size: 216 by 984 pixels (13.5 by 61.5mm)
   maximum supported frequency: 12MHz
   expanding computational grid...
   expanding sound speed grid...
   expanding density grid...
   computational grid size: 256 by 1024 pixels
   smoothing sound speed distribution...
   smoothing density distribution...
   ??? Error using ==> times
   Matrix dimensions must agree.

   Error in ==> kspaceFirstOrder2D at 827
   absorb_param = absorb_tau.*absorb_nabla1.*rho0;

   Error in ==> test_PML at 161
   [sensor_data, p_final] = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

   This happens because absorb_nabla1 and rho0 are 1024x256 whereas absorb_tau is 984x216.

   The same error occurs when I alter example_na_controlling_the_pml.m setting example_number = 3 and adding a medium.alpha_coeff and medium.alpha_power (specifically medium.alpha_coeff = 2*ones(Nz,Nx) and medium.alpha_power = 1.15.)

   Either way, I will eventually be doing 3D. I've only been doing 2D because I get NaNs when I run a 3D version of what I described in the first post. I'll clean up the .m I'm using and post it later today to see if anyone can reproduce what I'm seeing.

   Posted 13 years ago #

4. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,640

   Hi Jon,

   The current release of the toolbox (B.0.2) was coded assuming scalar values of the absorption coefficients (so these are not expanded when PMLInside is set to false). We are currently working on a new release so will take a look into what problems using heterogeneous absorption values creates. In the meantime, perhaps using a single value of absorption will be sufficient for you to move forward in your simulations?

   Thanks for letting us know the issues you are encountering,

   Brad.

   Posted 13 years ago #

5. jkiv
   

   Member
   Joined: Sep '10
   Posts: 7

   Here's the example: https://gist.github.com/dc9aa9351ec1997c8321

   Line 66 (% medium.alpha_coeff(PML_mask == 1) = 0; % [dB/(MHz^y cm)]), if uncommented, sets the medium.alpha_coeff = 0 in the PML.

   In order to see the effects I plot the log of the pressure in kspaceFirstOrder2D.m:

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   p_plot = 20*log10(abs(p_plot));
   plot_scale = [-30 0];
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   % update plot
   imagesc(kgrid.x(1,:)*scale, kgrid.z(:,1)*scale, p_plot, plot_scale);
   colormap(COLOR_MAP);
   ylabel(['z-position [' prefix 'm]']);
   xlabel(['x-position [' prefix 'm]']);
   axis image;
   drawnow;

   Also let me know what you think about the maximum pressure over both. If the high alpha_coeff` is in the PML, there is a lot of attenuation in the entire system:
   

   If it is 0 in the PML there is very little attenuation:
   

   Posted 13 years ago #

6. jkiv
   

   Member
   Joined: Sep '10
   Posts: 7

   Thanks Brad,

   This may be the source of my confusion. I'll see if I can find an appropriate scalar value to use.

   Posted 13 years ago #

7. bencox
   

   Developer
   Joined: Mar '10
   Posts: 285

   Hi Jon,

   Did setting a scalar absorption value sort out your problem?

   Ben

   Posted 13 years ago #

8. jkiv
   

   Member
   Joined: Sep '10
   Posts: 7

   Using a single value behaves as one would expect: there is no ghosting; attenuation is high with a high alpha_coeff, low with a low one; etc.

   I'm not too cozy with the math, so I wouldn't be able to figure out what is needed to have heterogeneous alpha_coeff. So I'm eager to see this addressed in the a future version. Hopefully it won't affect the algorithm too much as I've reimplemented the time-loop in CUDA and it's highly orchestrated (ie, a lot of my kernels are written for partial steps in calculations, storing results in then-unused memory spaces to avoid allocating any more memory than I need, etc.)

   Anyway, I'll keep you updated on anything else I run in to.

   Thanks,
   Jon

   Posted 13 years ago #

RSS feed for this topic

Reply

You must log in to post.