I am interested in such a multi-step approach in 2-D for elastic waves. Is there an update about the amplitude mismatch issue mentioned in this topic?
k-Wave
A MATLAB toolbox for the time-domain
simulation of acoustic wave fields
Simulation "layer by layer"
(32 posts) (7 voices)-
Posted 5 years ago #
-
Hi Brad,
I just read the article you mentioned and I must admit that I am not sure to understand why directly replacing the pressure values in the domain "leads to errors in the imposed spatial gradients" (although I don't doubt it, since the example you give is very convincing). I am especially curious about the example of the simulated field, where measurement errors cannot be responsible for these errors on gradients, since you set the "true" pressure values... By the way, are you talking about inconsistencies in the 3D spatial gradients (including voxels out of the measurement plane) or only inside this 2D plane?
Then, about the interpretation of the method, I understand that it is "restoring" the spatial consistency by allowing some (little) changes in the plane where the measurements were performed to make them physically consistent with the spatial gradients as they should be according to the equations, but how far do you need to put your equivalent source from the measurement plane? Is there a minimal distance?
Finally, in the example of the section III B, how does the "additive" mode in k-Wave behave? Does it give the same errors as the "Dirichlet" mode, since it seems a more "permissive" way to set the source term.
Warmest regards to everyone in this troubled times,
Anthony
Posted 4 years ago # -
Hi Anthony,
I'm facing the same problem you had before ,just like ,I let z denote my main propagation axis. When I compute the integral of Iz (intensity along z) over each xy-plane , and plot this quantity versus z, I get very strong oscillations close to the source plane.How did you solve this problem?Can I add CFL to solve it?Posted 4 years ago # -
Hi jackYANG,
Unfortunately, I did not find a general solution... But there is a simple (yet ugly) trick: you can make your "layers" overlap enough so the oscillations at the beginning of the layer L+1 vanishes before the z value corresponding to the end of the layer L :-) Then, when you merge all the layers in the zone where L and L+1 overlap, take the pressure values coming from the layer L where it overlaps witj L+1... See figure 2 of https://jtultrasound.biomedcentral.com/articles/10.1186/s40349-016-0056-9.
More elegantly, maybe that using the equivalent source computation method described in Brad's paper (see his last post in this very conversation) solves the problem. Please tell my if you try it and it works ;-) I would be very interested!
Hope it helps,
AnthonyPosted 4 years ago # -
Hi Anthony,
Thank you for your reply. I have downloaded the article you mentioned. I will try my best to solve this problem,If there is any progress, I will contact you again.Posted 4 years ago # -
Hi Anthony,
The errors are because the Dirichlet replacement in k-Wave is not exactly setting a boundary condition in a formal sense, just replacing the pressure values at discrete time points. If you run a simulation in 1D and record the pressure everywhere, then re-emit the recorded field from a single point using a Dirichlet source and record the pressure everywhere, and the compare the spatial gradients of the two fields at different time snapshots, you'll see they're not quite right.
We have an idea how to fix this by formally deriving a boundary value propagator (similar to how we have derived solutions for initial value problems and single frequency problems). Practically, your work around of applying the source term over a plane several grid points thick helps to alleviate this.
The key advantage of an additive source is that it doesn't act like a boundary, so if you have a heterogeneous medium, the scattered waves are not reflected from the source plane (they are if you use a Dirichlet source).
For the equivalent source, I've not really played around with the minimum distance between the "measurement" plane and the source plane. In principle it should work even if they're separated by a single grid point, although I'd advise trying it out as there might be strange behaviour if you're this close.
Brad.
Posted 4 years ago # -
Hi Brad,
thank you for your answer, it is getting clearer in my mind.
Warm regards,
AnthonyPosted 4 years ago # -
Hi Brad, Anthony,
I've been doing multi-stage simulations as well.
I have compared my multi-stage simulation (MULTI) results to the results obtained from a 'control' single-stage simulation (SINGLE), but I’m observing something odd in my results. For me, I am struggling with the transition between stages. The result is that the peak pressure is about 12% higher in MULTI than in SINGLE, as well as being offset in space (for MULTI, the focus is closer to the transducer surface than for SINGLE).
The simulation was performed in a nonlinear homogeneous acoustic medium. The single-stage simulation was performed at a spatial resolution of 8.2 PPW, with CFL around 0.3.
For the multi-stage simulation, all stages had the spatial resolution of 8.2 points per wavelength, but each stage was decreased in size as the acoustic beam narrowed. CFL was also around 0.3. Similar to Anthony, I recorded the source planes from the last several rows of a stage, then used them as the source in the following stage (the white bars). Source planes were 9 gridpoints thick, and were placed 10 gridpoints away from the non-PML domain boundary (so if the PML is 20 gridpoints thick, the source plane is 30 gridpoints away from the computational boundary edge). I record the last 4 cycles, and use interpftn to interpolate the recorded sensor pressures into becoming the new source plane.
I looked at reducing the CFL number from 0.3 to 0.1 in MULTI. Peak intensity is closer to results obtained from SINGLE, but still too high (by approx 10% and the MULTI focus is still closer to the transducer surface than SINGLE).
I am wondering if the interpolation from sensor to source is wrong somehow? Should I be using a windowing function, or time-shifting the interpolation somehow?
Daniel
Posted 4 years ago # -
Hi,
Unfortunately, I did not find a general solution on my side, and did not have the time to dig further.
From what I remember, I tweaked the parameters until it gave me appropriate results for my source. In particular, if I remember well, the errors were not decreasing with decreasing CFL (in an example cite above in the discussion, results were better with CFL = 0.5 than 0.3), and it was not monotonous either with increasing thickness of source "plane".
I don't remember, but there might be something related to the staggered grid though, (see user manual for details and maybe on the forum, I think we discussed something on this matter somewhere)... It needs to be accounted for properly.
Otherwise, you may try the equivalent source described by Brad in the paper cited above...
Sorry not to be of more help :-( Maybe Brad's team can help.
Posted 4 years ago # -
Thanks for your experience and advice, Anthony.
Yes, from my own work I noticed that changing the source plane thickness doesn't seem to change the peak intensity in a predictable way.
I will try CFL = 0.5 next, and using the acoustic holography method (calculateMassSource set of functions?) to generate the source planes.Posted 4 years ago # -
In case it's helpful, I dug out some of my old test codes from when Anthony was working on this and tidied one up. You can download it here. This particular example is for a linear CW simulation with a focused transducer, and gives errors on the order of a few percent for a single grid point source plane. There is also a slightly different example in the Equivalent Source Holography Example.
Posted 4 years ago #
Reply
You must log in to post.