k-Wave

1. KHO
   

   Member
   Joined: Nov '12
   Posts: 11

   Hello,
   I have problems understanding the the effect of the temporal staggered grid. restricting this to the pressure table 2.1 in the user manual states that the input pressure is offset by half a time step and the output pressure is offset by a whole time step. As I understand this, it would mean that defining source.p0 will effect that the actual wave starts to propagate half a time ahead and reading out the first entry (or column) of sensor_data would return a wavefield that has actually already travelled for one timestep minus a half time step (input delay), so after half a time step. But the first entry in sensor_data returns the exact distribution of p0. so how do I have to understand the offset between input and output?

   Posted 10 years ago #

2. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,640

   The table you mention actually corresponds to the use of time varying pressure sources (source.p) and time varying velocity sources (source.ux). These are added to the acoustic fields at different points in the code as the injection of mass and force, respectively. In the case of an initial pressure distribution (source.p0), the pressure and the velocity are set to the user values after the end of the first time-step. They are not temporally offset, so the behaviour you have observed is correct. Apologies that this in not very clear from the table!

   Brad.

   Posted 10 years ago #

3. KHO
   

   Member
   Joined: Nov '12
   Posts: 11

   Thanks alot for your detailed explaination! I thought adding a source.p of only one time step would equal making use of source.p0 but that's obviosly not the case. Now everything totally makes sense!
   Best regards

   Posted 10 years ago #

4. Anthony
   

   Member
   Joined: Apr '13
   Posts: 96

   Hi Bradley,

   Your answer is very clear, but I am confused about something. In the mesocentre I am using for my simulations the maximum duration of a job is 24 hours but some of my simulations require more time. In this context, I need to stop my simulation and I would like to "continue" it by setting the appropriate initial conditions.

   I planned to simply run a new simulation with my source.p term and add a source.p0 term corresponding to the end of the previous simulation. However, if I understand well, my source.p0 will be "set to the user values after the end of the first time-step." i.e. a time-step too late for me... Am I right? If yes, is there any workaround?

   Best regards,
   Anthony

   Posted 10 years ago #

5. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,640

   Hi Anthony,

   To restart the simulation, you need to reset the previous values of both the particle velocity and the acoustic density used within the code. The values of the particle velocity can be reset exactly using the velocity source terms, while the values of the acoustic density can be reset approximately using the pressure source terms.

   To explain, when using a pressure source, the user values are first scaled and then added to acoustic density fields in the simulation as a mass source. However, the acoustic density field is artificially split to allow the PML to be applied. The scaled pressure source values are thus divided evenly between the Cartesian components of the density rho_x, rho_y, rho_z.

   To give you an example of how you might restart a simulation, in the first simulation, set the value of sensor.record to include the final pressure and velocity fields (in addition to the fields you are interested in recording), e.g.,

   sensor.record = {'p', 'p_final', 'u_final'};

   In the second simulation, use the final pressure and velocity fields from the previous simulation to set the initial conditions. For example, for a 3D simulation the code would look something like
   

   % set the source using the checkpointing data
   clear source;
   source.p_mask = ones(Nx, Ny, Nz);
   source.p = reshape(sensor_data_1.p_final, [], 1);
   source.p_mode = 'dirichlet';
   source.u_mask = ones(Nx, Ny, Nz);
   source.ux = reshape(sensor_data_1.ux_final, [], 1);
   source.uy = reshape(sensor_data_1.uy_final, [], 1);
   source.uz = reshape(sensor_data_1.uz_final, [], 1);
   source.u_mode = 'dirichlet';

   This will force the entire acoustic density and particle velocity fields to be replaced by the final values from the previous simulation. Note, the second simulation will need to run for one extra time step, as the first time step will be used to reset the fields. You will also need to remove this time step from any output data. For example

   % join the sensor data together
   output_p = [sensor_data_1.p, sensor_data_2.p(:, 2:end)];

   In the linear, lossless case with the PML absorption set to zero, this approach for checkpoint / restart will work exactly. However, if the PML absorption is not zero (which will normally be the case), there will be an error introduced because the reset pressure values are divided evenly between the Cartesian components of the acoustic density. While the total acoustic density will be correct, where rho = rho_x + rho_y + rho_z, the individual components in the restarted simulation will be different to the components in the previous simulation.

   If your waves haven't yet reached the PML, this error will be very small. However, if the waves are propagating through the PML at the point where the simulations are split, you may see a small reflection from the PML. This error might be on the order of a few percent.

   In the nonlinear and absorbing case, there will also be a small error introduced as the calculation of the final pressure values from the acoustic density will include these effects, which are not "undone" when reseting the acoustic density using a pressure source.

   I hope that makes some sense! If a restart capability is critical for your work, we could certainly look at adding it for the next release (by resetting the acoustic density values directly and thus avoid the errors mentioned above). Are you using the C++ codes, or just the MATLAB code?

   Kind regards,

   Brad.

   Posted 10 years ago #

6. Anthony
   

   Member
   Joined: Apr '13
   Posts: 96

   Hi Bradley,

   Thanks for your response: crystal clear, as usual :-)
   I am directly using the C++ code (for memory reasons). I must admit that it would be great for me to have this feature, if it does not represent too much work for you. In general, I always try to split my long simulations, so if something happens, I don't lose everything and I can restart from last "time segment".

   Best regards,
   Anthony

   Posted 10 years ago #

RSS feed for this topic

Reply

You must log in to post.