k-Wave

1. xiong
   

   Member
   Joined: Aug '10
   Posts: 8

   I am doing thermoacoustic research, which is quite similar to photoacoustic. I would like to try k-wave for my project. But I have the following two questions.

   1. In thermoacoustic simulation, the initial value for solving the differential equations is the absorbed heat (not pressure). How to set such initial condition?

   2. The duration of the source pulse (microwave pulse) sometimes cannot be ignored. Take 1 us pulse as an example. The acoustic wave is generated continously and starts to propagate in the first 1 us. Thus the initial condition should be properly set in the first 1 us. Is it possible to use k-wave to do this?

   Thank you very much.

   Xiong Wang

   Posted 13 years ago #

2. bencox
   

   Developer
   Joined: Mar '10
   Posts: 285

   Dear Xiong Wang,

   Thank you for raising these points.

   1) The current version of k-Wave models acoustic pressure propagation, so the initial condition must be given as a pressure. In photoacoustics it is normal to relate the initial acoustic pressure distribution to the absorbed energy density (the heat deposited per unit volume) by the Grueneisen parameter. In other words, the increase in temperature due to the absorbed EM pulse and the accompanying increase in pressure are related via a thermodynamic identity (see Section IIA in our paper http://www.k-wave.org/papers/2005-Cox-JASA.pdf ).

   2) The heating pulse (the absorbed EM pulse) can be treated as a delta function in time as long as the pulse duration is much shorter than the charateristic acoustic travel time, ie. it is much less than d/c where d is a charateristic length of the absorbing region and c is the speed of sound. Under this assumption, the problem reduces to an initial value problem (IVP) and so setting an initial acoustic pressure is appropriate.

   There are two ways to deal with longer pulses which don't satisfy this 'stress confinement condition':
   a) You could run the model using an initial condition (ie. assuming an EM pulse that is a delta function in time) and convolve the actual pulse shape with the time series output afterwards.
   b) Alternatively, you could use a time-varying pressure source. How to do this in k-Wave is described (for a sine wave source) in this example http://www.k-wave.org/documentation/example_us_homogeneous_medium.php, which is also in the Toolbox. You would just need to change the shape of the pulse you are using to whatever temporal shape you think your pulse has. When a thermoacoustically generated acoustic pressure is included as a source term, it has the form of the time derivative of the absorbed energy density. So if your EM pulse is Gaussian, then the pressure source will be the time derivative of a Gaussian.

   3) There is one more thing to note. The equations that k-Wave is based on (linear acoustics) assume that thermal conduction can be neglected. This is reasonable when heat diffusion takes place on a much longer time scale than acoustic propagation, ie. when the thermal relaxation time >> d/c where d is a characteristic length of the absorbing region, and c is the speed of sound. If this is not the case then the temperature and pressure fluctuations become coupled and the current version of k-Wave can't model that.

   I hope that helps - good luck with your modelling,

   Ben

   Posted 13 years ago #

3. xiong
   

   Member
   Joined: Aug '10
   Posts: 8

   Dear Ben,

   Thanks for your reply.

   I think using the time varying pressure source can deal with my project. But I still have a problem with the setting of "source.p". Because the source (known) for my TAI is the "heat energy deposited per unit volume and per unit time (let's note it by H1)", I need to transfer it to pressure. From the paper you mentioned I know how to relate the "initial pressure" with the "heat deposited per unit volume (let's note it by H2)" for a delta funtion pulse. But for a time varying long pulse, how to get the time varying pressure which is to be used by "source.p"? I plan to use the following way to do it. Please take a look at it and tell me if it is right.

   Take a 1us electromagnetic source as an example. I can get the time varying electric field E at all spatical points. Then the H1 is obtained by

   H1(x,y,z,t) = 0.5 * conductivity(x,y,z) * (E(x,y,z,t))^2

   Divide the 1us pulse into 100 equal parts or 100 short pulses (10 ns). The 100 10ns pulses heat the tissue one by one. Thus, the stress confinement is satisfied in each 10ns pulse for my accuracy requirement. The initial pressure for each 10ns pulse is got by

   p0(x,y,z,m) = Gruneisen coefficient(x,y,z) * H2(x,y,z,m)
   (m means the mth 10ns pulse, m=1,2,...,100)

   The H2(x,y,z,m) is calculated by time integration of H1(x,y,z,t,m) in the 10ns (from 10*(m-1) ns to 10*m ns). The 3D k-wave code is then run for 100 times for all the 10ns pulses to get the pressure(x,y,z,t,m). Finally the total pressure(x,y,z,t) is obtained by properly combining all of the pressure(x,y,z,t,m).

   Please tell me if this is right and if you have any other suggestions.

   Thank you for your time.

   Xiong

   Posted 13 years ago #

4. xiong
   

   Member
   Joined: Aug '10
   Posts: 8

   Ben,

   Please ignore my last post. I made some mistakes in it. The following is the correct one.

   I think using the time varying pressure source can deal with my project. But I still have a problem with the setting of "source.p". Because the source (known) for my TAI is the "heat energy deposited per unit volume and per unit time (let's note it by H1)", I need to transfer it to pressure. From the paper you mentioned I know how to relate the "initial pressure" with the "heat deposited per unit volume (let's note it by H2)" for a delta funtion pulse. But for a time varying long pulse, how to get the time varying pressure which is to be used by "source.p"? I plan to use the following way to do it. Please take a look at it and tell me if it is right.

   Take a 1us electromagnetic source as an example. I can get the time varying electric field E at all spatical points. Then the H1 is obtained by

   H1(x,y,z,t) = 0.5 * conductivity(x,y,z) * (E(x,y,z,t))^2

   Assume the "dt" set by function "makeTime" is 10ns and the stress confinement is satisfied for a 10ns pulse for my accuracy requirement, divide the 1us pulse into 100 equal parts or 100 shorter pulses (10 ns). The 100 10ns pulses heat the tissue one by one. The initial pressure for each 10ns pulse is got by

   p0(x,y,z,m) = Gruneisen coefficient(x,y,z) * H2(x,y,z,m)
   (m means the mth 10ns pulse, m=1,2,...,100)

   The H2(x,y,z,m) is calculated by time integration of H1(x,y,z,t,m) in the 10ns (from 10*(m-1) ns to 10*m ns). The first 100 terms of source.p will be p0(x,y,z,m) (m from 1 to 100). The rest terms of source.p are 0.

   Please tell me if this is right and you are very welcome to give any other suggestions.

   Thank you for your time.

   Xiong

   Posted 13 years ago #

5. bencox
   

   Developer
   Joined: Mar '10
   Posts: 285

   Dear Xiong,

   Your question has made me think of a few points. I hope they’re helpful.

   1) k-Wave will apply the time integration that you mention automatically. It will multiply each sample of the time-varying source.p by the timestep dt. I think you therefore just need to write:

   source.p = Gruneisen coefficient * H1;

   where this H1 is the sampled version of your known time varying H1(t). See comment (3) below about checking the amplitude is correct.

   2) Just to clear up a possible confusion. I said in my last post that "if your EM pulse is Gaussian, then the pressure source will be the time derivative of a Gaussian". That is not true for k-Wave. The source term for the underlying wave equation is the time derivative, but that is taken care of in the code. The input source.p should be the actual shape of the EM pulse you are simulating. For instance, take the example "Point Source In A Homogeneous Propagation Medium Example" as a starting point, comment out the definition of the sinusoidal source and replace it with a Gaussian. Note the time offset.

   % define a time varying photoacoustic Gaussian pulse
   % with the Gaussian normalised so its integral sum(gaussian)*dt = 1
   source_mag = 1;
   gruneisen = 0.11;
   pulsewidth = 1e-6;
   t_offset = 5e-6;
   gaussian = exp(-((kgrid.t_array-t_offset)/pulsewidth).^2)/(pulsewidth*sqrt(pi));
   source.p = source_mag*gruneisen*gaussian;

   Also adjust the following line to make sure the plot scaling is sensible:

   [sensor_data p_final] = kspaceFirstOrder2D(kgrid, medium, source, sensor,... 'PlotScale',[-1e4 1e4]);

   Just to be tidy, you might also want to delete the instructions that "plot the final wave-field".

   3) It might be a good idea to check that you are getting the amplitude you expect by putting in a plane source (so there’s no geometrical spreading). Adjusting the above example a little:

   % define a line of source elements to generate a plane wave
   source.p_mask = zeros(Nz, Nx);
   source.p_mask(Nz/2, :) = 1;

   Also set PMLAlpha to turn off the PML at the sides so the plane wave remains truly plane:

   [sensor_data p_final] = kspaceFirstOrder2D(kgrid, medium, source, sensor,... 'PlotScale',[-7e4 7e4],'PMLAlpha',[0 4]);

   4) Final thing. There’s a bug, which we hope to fix for the next release, which means that if you want quantitatively accurate results when you are using a time-varying source you should set dx = dz to be equal.

   Kind regards,

   Ben

   Posted 13 years ago #

6. xiong
   

   Member
   Joined: Aug '10
   Posts: 8

   Dear Ben,

   Thanks for your quick reply.

   This time I just have one quick quiestion. I want to make some change to the length of "t_array" and value of "dt" defined by "makeTime". Since my model is 3D, the t_array set by your code is too long and takes very long time in both Photoacoustic Simulation and Photoacoustic Image Reconstruction. I tried to make t_array=t_array(1:1000), but it does not work. The parameters printed to the command line still show "time steps: 2020". How can I do this?

   Thank you.

   Best regards.

   Xiong

   Posted 13 years ago #

7. bencox
   

   Developer
   Joined: Mar '10
   Posts: 285

   Dear Xiong,

   I tend to do it this way:

   [kgrid.t_array dt] = makeTime(kgrid, medium.sound_speed);
   kgrid.t_array = kgrid.t_array(1:1000);

   or

   [kgrid.t_array dt] = makeTime(kgrid, medium.sound_speed);
   kgrid.t_array = (0:999)*dt;

   Of course, you don't have to use makeTime at all, you could just write:

   dt = 1e-8; %or whatever you want here
   kgrid.t_array = (0:999)*dt;

   The advantage of using makeTime is that dt is chosen to satisfy a certain CFL stability condition. CFL = 0.3 by default, but for greater accuracy you could adjust it:

   cfl = 0.1;
   [kgrid.t_array dt] = makeTime(kgrid, medium.sound_speed,cfl);

   This will return a smaller time step dt. Hope this helps,

   Ben

   Posted 13 years ago #

8. xiong
   

   Member
   Joined: Aug '10
   Posts: 8

   Dear Ben,

   Now I know how to change t_array and dt. Thank you.

   I planned to use the time-varying source. But the array of source.p has a very big size which cannot be handled by my computer memory. I have to change to assuming a delta function first and then do convolution. You siad k-wave will multiply each sample of the time-varying source.p by the timestep dt. So

   source.p = Gruneisen coefficient * H1;

   H1 is the sampled version of my known time varying H1(t). How to set source.p0? Just times the Gruneisen coefficient with H1? If this is right, how can an arbitrarily sampled H1 give the correct results?

   Thank you.

   Best regards,

   Xiong

   Posted 13 years ago #

9. bencox
   

   Developer
   Joined: Mar '10
   Posts: 285

   Dear Xiong,

   If all of your source elements are varying in time in the same way, ie. if they all have the same input, then source.p need only contain one time series, not one for every sensor element in sensor.mask. Will that help reduce the size of it?

   If you're going to convolve the output with your time-varying input, then do you just need a single 1 in source.p?

   We've been running some examples this week and think there may sometimes be a scaling issue with the time-varying pressure input, so if it looks wrong then please let us know. We're planning to fix it and write an example for a time varying thermoacoustic source for the next beta release, which we hope will be the end of next week.

   I hope you can get your example to work.

   Kind regards,

   Ben

   Posted 13 years ago #

RSS feed for this topic

Reply

You must log in to post.