k-Wave

1. wgrissom
   

   Member
   Joined: May '16
   Posts: 1

   Hi, I am running a relatively small 3D simulation, but it generates enormous (> 100 GB) hd5 files in my /tmp/ folder, and I run out of disk space. I am using k-Wave with the multi-threaded CUDA executable on Ubuntu 16.04, calling it from MATLAB R2015a. I have about 140 GB free on the drive, and I have 96 GB of RAM and dual 6-core Xeons.

   The output is below. My question is, should the files really be that big for this problem size? I am trying to record 'p_max', 'p_min', 'p_rms', 'I_avg' at each location in the simulated volume. Could that be my issue (i.e., I need to reduce my number of sensors, and perhaps repeat with multiple sensor positions)? Thank you for any help.

   Running k-Wave simulation...
   start time: 21-Sep-2016 22:37:19
   reference sound speed: 1484m/s
   dt: 80.8625ns, t_end: 110.0539us, time steps: 1362
   input grid size: 316 by 183 by 183 grid points (126.4 by 73.2 by 73.2mm)
   maximum supported frequency: 1.855MHz by 1.8449MHz by 1.8449MHz
   kspaceFirstOrder_saveToDisk
   precomputation completed in 26.3772s
   saving input files to disk...
   completed in 1min 14.5148s
   ┌───────────────────────────────────────────────────────────────┐
   │ kspaceFirstOrder3D-CUDA v1.1 │
   ├───────────────────────────────────────────────────────────────┤
   │ Reading simulation configuration: Done │
   │ Selected GPU device id: 0 │
   │ GPU device name: Tesla K20c │
   │ Number of CPU threads: 24 │
   ├───────────────────────────────────────────────────────────────┤
   │ Simulation details │
   ├───────────────────────────────────────────────────────────────┤
   │ Domain dimensions: 316 x 183 x 183 │
   │ Simulation time steps: 1362 │
   ├───────────────────────────────────────────────────────────────┤
   │ Initialization │
   ├───────────────────────────────────────────────────────────────┤
   │ Memory allocation: Done │
   │ Data loading: Done │
   │ Elapsed time: 0.44s │
   ├───────────────────────────────────────────────────────────────┤
   │ FFT plans creation: Done │
   │ Pre-processing phase: Done │
   │ Elapsed time: 0.11s │
   ├───────────────────────────────────────────────────────────────┤
   │ Computational resources │
   ├───────────────────────────────────────────────────────────────┤
   │ Current host memory in use: 1647MB │
   │ Current device memory in use: 3649MB │
   ├───────────────────────────────────────────────────────────────┤
   │ Simulation │
   ├──────────┬────────────────┬──────────────┬────────────────────┤
   │ Progress │ Elapsed time │ Time to go │ Est. finish time │
   ├──────────┼────────────────┼──────────────┼────────────────────┤
   │ 0% │ 1.665s │ 1131.926s │ 21/09/16 22:58:00 │
   │ 5% │ 104.034s │ 1920.166s │ 21/09/16 23:12:51 │
   │ 10% │ 197.664s │ 1753.197s │ 21/09/16 23:11:38 │
   │ 15% │ 314.691s │ 1765.936s │ 21/09/16 23:13:47 │
   │ 20% │ 426.971s │ 1695.417s │ 21/09/16 23:14:29 │
   │ 25% │ 540.264s │ 1611.315s │ 21/09/16 23:14:59 │
   │ 30% │ 655.427s │ 1521.870s │ 21/09/16 23:15:24 │
   │ 35% │ 777.675s │ 1438.210s │ 21/09/16 23:16:03 │
   │ 40% │ 908.100s │ 1357.161s │ 21/09/16 23:16:52 │
   │ 45% │ 1019.975s │ 1242.575s │ 21/09/16 23:16:49 │
   │ 50% │ 1135.512s │ 1128.862s │ 21/09/16 23:16:51 │
   └──────────┴────────────────┴──────────────┴────────────────────┘
   ┌───────────────────────────────────────────────────────────────┐
   │ !!! K-Wave experienced a fatal error !!! │
   ├───────────────────────────────────────────────────────────────┤
   │ Error: Could not write into "" dataset. │
   ├───────────────────────────────────────────────────────────────┤
   │ Execution terminated │
   └───────────────────────────────────────────────────────────────┘
   /bin/bash: line 1: 15377 Segmentation fault (core dumped) ./kspaceFirstOrder3D-CUDA -i /tmp/kwave_input_data21-Sep-2016-22-37-19.h5 -o /tmp/kwave_output_data21-Sep-2016-22-37-19.h5 --p_max --p_min --p_rms --u_non_staggered_raw --p_raw
   export LD_LIBRARY_PATH=; cd /buffyexport/home/grissowa/Dropbox/code/zebrography/sim/kWave/k-wave-toolbox-version-1.1.1/k-Wave/binaries/; ./kspaceFirstOrder3D-CUDA -i /tmp/kwave_input_data21-Sep-2016-22-37-19.h5 -o /tmp/kwave_output_data21-Sep-2016-22-37-19.h5 --p_max --p_min --p_rms --u_non_staggered_raw --p_raw: Segmentation fault
   Error using h5readc
   The HDF5 library encountered an error and produced the following stack trace information:

   H5F_super_read file signature not found
   H5F_open unable to read superblock
   H5Fopen unable to open file

   Posted 8 years ago #

2. Jiri Jaros
   

   Developer
   Joined: Feb '12
   Posts: 115

   Hi wgrissom,

   The problem is in the intensity 'I_avg'. The problem is that pressure and velocity are spatially and temporally staggered. So to calculate avarage intensity, we need to shift velocity and pressure in space (easy), but unfortunately in time as well. Since the GPU doesn't have enough memory to buffer more timesteps, we need to flush data on disk.

   And this is the core of the problem. You sample the whole simulation domain, over 1362 timestep. A single timestep costs you 316 x 183 x 183 * 4B * 4 matrices ~ 160MB. Since you're sampling over 1362 timesteps, it sums up to 214 GB.

   There are two questions:
   (1) Do you really have to sample intensity over the whole domain? If not, make the sensor mask smaller (just to cover the focus area).
   (2) Do you really have to sample intensity over the whole time period? For most of the time there's nothing in the domain (until the wave comes...). I'd recommend to sample 1 period of the signal after it reaches a steady state.

   Jiri

   Posted 8 years ago #

3. bastg
   

   Member
   Joined: Jun '19
   Posts: 7

   A follow up on this thread: I understand that the exact computation of the acoustic intensity requires dumping the entire P and U solutions (in time & space) on disk. For realistic simulations in 3D, this results in files well over 1TB (as mentioned many times on this forum), which is not feasible.

   I wonder if there is an approximate solution to this problem. I just tried, on a small example, saving p_rms and u_rms, and then computed the averaged intensity as I_rms = abs(p_rms ) .* sqrt( ux_rms.^2 + uy_rms.^2 + uz_rms.^2 ). This is strictly speaking incorrect as this assumes that one can permute the temporal averaging with the multiplication, i.e. I am assuming that <P*U> = <P>*<U>. That being said, the error is not that bad. For F=100kHz, dt=10us and a focused transducer, I am finding a mean error of 2% and a peak error of 18%. Honestly, given that the acoustic properties are probably within ~30% of reality (I am simulating the human skull), this kind of error seems acceptable.

   The other source of error is that P and U are staggered in space. I tried correcting this shift using the following code:
   >> p2 = fourierShift( p , 1/2 , 1 );
   >> p3 = fourierShift( p2 , 1/2 , 2 );
   >> p4 = fourierShift( p3 , 1/2 , 3 );
   and using p4 instead of p in the calculation of I_rms, but the error goes up instead of going down (it also goes up when using a -1/2 shift). I may not be correcting the spatial staggering properly...

   My overall question is this: What is the best approximate estimate of the acoustic intensity that one can obtain using only the averaged variables (in time) that are outputted by K-WAVE?

   Posted 4 years ago #

4. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi bastg,

   If your simulations are linear, the most memory efficient way would be using a plane wave approximation with the pressure amplitude extracted using p_max_all:

   I = sensor_data.p_max_all.^2 ./ (2 * medium.density .* medium.sound_speed);

   Provided the time step is small enough, taking the maximum pressure will be a reasonable approximation to the actual spectral amplitude.

   If you want to record the particle velocity on the regular grid (instead of spatially staggered), you can record u_non_staggered. The particle velocity components are all staggered in different directions (i.e., ux is staggered in x, uy in y, etc), which is why the error goes up using your code.

   One of Jiri's students has been working on a memory efficient way to extract the intensity directly using the C++ code, but I think it's a little way from the production version at this point.

   Hope that helps,

   Brad

   Posted 4 years ago #

5. bastg
   

   Member
   Joined: Jun '19
   Posts: 7

   A follow-up: The "plane wave method", which computes I_avg ~ <P>^2/(2*rho*c) gives an error of 0.0120 whereas the "separable method", which computes I_avg ~ <P>*<U>, has an error of 0.0014 (after correctly de-staggering ux, uy and uz as mentioned in previous post). This is for a 100kHz simulation, curved transducer with R=80mm and aperture=61mm in a uniform water bath.

   This means that the separable method is ~8X more accurate than the plane wave method in this particular situation. In addition, the error of the separable method is mostly in the very near field, close to the focus the error is only ~2%.

   I have also ran a comparison @ 500kHz, and I found an even more pronounced difference: The plane wave method error at the focus is ~48%. The error of the separable approximation at the focus is ~1%.

   This is all anecdotal, and such a comparison should be done in a more rigorous way, but I thought it 1) could be interesting to other folks and 2) spark a discussion about why the separable method seem to work so well. I would be interested in hearing what the developers of k-wave think of this (and maybe see if they can reproduce this).

   Thanks -- Bastien

   Posted 4 years ago #

6. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi Bastien,

   Thanks for your patience - it took me a while to find time to take a look. Here's my explanation.

   To calculate the components of the average acoustic intensity, you want integrate p*u over an acoustic period. For a monochromatic source of frequency w the period is T = 2*pi/w. The pressure and the components of the particle velocity can be written as:

   p = a1*sin(w*t + p1)
   u = a2*sin(w*t + p2)

   where a1 and a2 are the amplitudes, and p1 and p2 are the phases (similarly for all the vector components of u). Multiplying these together, integrating over time from 0 to T, and dividing by T to take the average then gives:

   I_avg = integral(p*u) / T
   I_avg = integral(a1*sin(w*t + p1)*a2*sin(w*t + p2)) / T
   I_avg = a1*a2*cos(p1 - p2)/2

   Now if the wave is fairly planar, the phase components p1 and p2 will be similar, and the cos(p1 - p2) term will be approximately one. This leaves:

   I_avg ≈ a1*a2 / 2 = p_rms * u_rms

   which is what you've calculated.

   I coded this up for a simple focused field example in 2D, and the errors are on the order of 5% for this approach, compared to 10% for the plane wave approximation. I've also plotted the cosine terms, which for the p*ux components (which is the dominant one) is close to 1 in central part of the field. You can download my code here. Actually, it surprised me that it worked quite a bit better than the conventional plane wave approximation.

   Let me know any queries or comments.

   Brad

   Posted 4 years ago #

RSS feed for this topic

Reply

You must log in to post.