[New Miniapp] Time-varying pressure source for acoustic2D

[alexpattyn]

Not sure if you want to avoid domain specific miniapps, but I am working towards modifying the acoustic 2D miniapp for ultrasound and photoacoustic modeling. Currently I just played around with getting a time-varying source working.

My current implementation just sets the pressure at a given pixel based on some source signal vector.

freq = 1/dt/5 # divide by 5 to stay below the nyquist limit of the grid's sampling frequency
trange = range(0,2pi,length=nt)
amp = 1
signal = amp*sin.(freq*trange)

# Time loop
for it = 1:nt
        P[nx÷2,ny÷2]=signal[it] # Set source signal value
        @parallel compute_V!(Vx, Vy, P, dt, ρ, dx, dy)
        @parallel compute_P!(P, Vx, Vy, dt, k, dx, dy)
        
        # Visualisation
 end

I can clean my code up a bit and submit it as an acoustic2D with time varying source miniapp, but I wanted to see if this would be a recommended way of implementing a time-varying source? Or if this poses obvious performance penalties that users shouldn't be exposed to.

[luraess]

Thanks @alexpattyn for your suggestion. I'd say the goal of the miniapps is exactly to feature domain specific examples on how to implement and leverage ParallelStencil (as this may often be good documentation forgetting other folks started using it).

Regarding your addition suggestion: form my point of view it would be valuable to have a version of the acoustic miniapp featuring time varying source. I was actually thinking if it would make sense to add a new version of the acoustic apps that would include both the time-varying source and absorbing boundaries (since you suggested working on an implementation in Discourse)

The question I'd like @omlins thoughts on are:

• shall the addition also be propagated to 3D and 3D multi-xpu implementations
• shall the operation setting the source be implemented in a @parallel_indices kernel like (not tested):

@parallel_indices (ix,iy) function add_source!(P::Data.Array, Signal::Data.Array, it::Int)
    if (ix<=size(P, 1) && iy<=size(P, 2)) P[nx÷2, ny÷2] = Signal[it]  end # Set source signal value
    return
end
# [...]
freq   = 1.0 / dt / 5.0 # divide by 5 to stay below the nyquist limit of the grid's sampling frequency
trange = range(0,2pi,length=nt)
amp    = 1.0
Signal = Data.Array(amp*sin.(freq*trange))
# [...]
# Time loop
for it = 1:nt
        @parallel add_source!(P, Signal, it)
        @parallel compute_V!(Vx, Vy, P, dt, ρ, dx, dy)
        @parallel compute_P!(P, Vx, Vy, dt, k, dx, dy)
        # Visualisation
 end

• define (xsource,ysource) as variable and pass it as argument

Upon @omlins feedback, you could draft a PR adding the miniapp version(s) we decided upon for review (including the corresponding entries in the README).

[omlins]

Thanks @alexpattyn for your intentions to create a mini-app. I share @luraess view that mini-apps should solve relevant problems from different domains and mini-apps will often be specific to some domains. @luraess noted:

I was actually thinking if it would make sense to add a new version of the acoustic apps that would include both the time-varying source and absorbing boundaries

I totally believe that it would make sense to have these two things in the same code as we certainly don't want to have a large amount of mini-apps that differ only slightly; that would not be very useful. That said, it would be great if your mini-app would include different commonly used "elements" typically found in your domain!

shall the addition also be propagated to 3D and 3D multi-xpu implementations

Yes, absolutely, that will make them much more relevant. I think this is fundamental.

shall the operation setting the source be implemented in a @parallel_indices kernel

That should probably not matter for performance as long as you set the source only in grid point.

To conclude, @alexpattyn , we would very much appreciate a PR with a mini-app (in 2-D, 3-D and multi-XPU) as discussed above. :)

[alexpattyn]

Sounds good! I would think two good examples from my field would be:

1. Photoacoustic simulation with a PML and a sensor array (e.g. a linear or ring geometry). This would basically mean we have some initial pressure map at t=0 and then sample the grid at particular locations as the wave moves forward in time.
2. Ultrasound sim with a PML, sensor array, and some time-varying source (e.g. a gaussian pulse with a given center frequency and bandwidth).

The time-varying source and sensors are easy enough to implement, but I don't have much experience with integrating a PML. It looks like I can modify the existing acoustic2D mini-app and base it on k-Wave: MATLAB toolbox for the simulation
and reconstruction of photoacoustic wave fields, since they described the first-order acoustic equations and added on an attenuation term. I'll take a look at getting a 2D example working!

[luraess]

That sounds great @alexpattyn. Most valuable would certainly be to have time-varying source and PML within a single code. Output could be a double panel figure with forward simulation evolution on one panel and sensor arrays output on the other panel.

The time-varying source and sensors are easy enough to implement, but I don't have much experience with integrating a PML.

Indeed, PML may be the tricky bit but would also be a very nice addition. Feel free to ping us when you have a 2D prototype.

[alexpattyn]

A brief update:

Updated the two equations

@parallel function compute_V!(Vx::Data.Array, Vy::Data.Array, P::Data.Array, αx::Data.Array, αy::Data.Array, dt::Data.Number, ρ::Data.Number, dx::Data.Number, dy::Data.Number)
    @inn(Vx) = @inn(Vx) - dt*(1/ρ*@d_xi(P)/dx + @inn(αx)*@inn(Vx))
    @inn(Vy) = @inn(Vy) - dt*(1/ρ*@d_yi(P)/dy + @inn(αy)*@inn(Vy))
    return
end

@parallel function compute_P!(P::Data.Array, Vx::Data.Array, Vy::Data.Array, dt::Data.Number, ρ::Data.Number, dx::Data.Number, dy::Data.Number)
    @all(P) = @all(P) - dt*(ρ*(@d_xa(Vx)/dx + @d_ya(Vy)/dy)+(@all(αx) + @all(αy))*@all(P)) 
    return
end

I got a little confused on how to break the pressure field into x and y components as described in the paper I linked above. I tried taking the Px = P*cos(theta)^2 and Py = P*sin(theta)^2 so P = Px + Py which was a condition in the paper above, however this keep leading to numerical dispersion. So I found this paper by some of the same authors and there equation (31) shows that the attenuation terms are just added together for the first-order pressure equation. By following this second paper I got the below:

It seems that the PML causes some back-reflection itself, so I will have to look into this a little further. I know sometimes people define the PML as a gradient of attenuation to limit that back-reflection but I have to test that out. As far as I am aware k-wave doesn't do this, and I would prefer to avoid that additional complexity.

Once that is finished I have to clean my 2D code up and can submit it for review!

@omlins @luraess One thing I wanted to ask was is there a reason to use @all(...) vs @inn(...)? Since I don't care about what happens at the boundaries - due to the PML - are they effectively the same? I ask since in the default acoustic2D mini app we use @inn() for the particle velocity and @all() for the pressure field.

[luraess]

Thanks @alexpattyn for sharing those details. Your absorbing boundaries implementation looks promising. As far as I remember from early experiments, the PML approach should ideally not generate any reflections (while the less formal "sponge" approach would). In the Cox et al. paper you refer to, it seems they use the single parameter α (as you in the code snippet) to absorb at the boundaries. If you didn't yet, it may be worth trying to define α as specially increasing within the absorbing layer (in e.g. the log space). Also, there may be quite some literature on PML coming from geophysics.

Also, I wouldn't split the pressure update as it's not a vector field; there is only one pressure.

is there a reason to use @all(...) vs @inn(...)

The reason is that we are using a staggered grid; pressure is defined in the cell centre, while velocity is defined, like a flux, at cell faces. For this reason, one needs different array sizes (independent from the boundary condition) to get sizes to match upon taking the derivatives.

Once that is finished I have to clean my 2D code up and can submit it for review!

Sure! And then we should do the 3D multi-XPU version based on that 🙂

[alexpattyn]

Another Update

I made the PML a gradient between 0 and 1 and this seems to yield better results. I still get some back-reflections, but they are orders of magnitude lower in amplitude. I will continue to take a closer look at this and then compare it against k-wave to see if I get similar results - I can also benchmark it against their OpenMP C++ code.

Acoustic 2D with PML

Acoustic 2D without PML

Sure! And then we should do the 3D multi-XPU version based on that 🙂

Definitely! I have a couple of NVIDIA cards I can test it against to see how the GPU/XPU version runs.

[luraess]

Thanks for the update. I guess the making the absorbing layer gradually absorbing is indeed one way to go. I am still wondering wether this approach is then not rather spine boundaries as from what I recall PML should be, as the name suggests, perfectly matched (i.e. provide the exact correction).

I can also benchmark it [...]

I'd guess doing some kind of benchmarking activity would be good just to ensure one has a solid approach.

I have a couple of NVIDIA cards I can test it against to see how the GPU/XPU version runs

That's be great 🙂

[alexpattyn]

I am still wondering wether this approach is then not rather spine boundaries ...

From the k-wave documentation:

Consider the case of a wave propagating in the x direction. If αx is constant, between the edge of the PML and one grid point inside,the wave will be forced to decrease by a factor of exp(−αx∆x/c0). If αx is large then the PML will impose a large
gradient across the PML boundary, which will cause a reflection of the incoming wave.

They then go on to say that the attenuation coefficient should vary spatially, and they use the equation below:

So it looks like we are on the right track still!

Edit: They also say after discussing the specifics of the PML:

This corresponds to around 4 or 5 decimal places of accuracy, which should be sufficient for most simulations

So I suppose it should be up to the user to determine their threshold for accuracy - i.e. we will never have a perfect PML ironically. Once I get their implementation of the PML working I would consider that part done.

[luraess]

That's interesting, thanks for reporting. I'd say that indeed, if the errors are sufficiently small we could go with it. Would be good, as you proposed, to quantify that briefly in a kind of benchmark.

[alexpattyn]

@luraess Sorry for disappearing on this. Had a busy end of the fall semester and start of this semester. Interested in completing this issue though. How would you recommend quantifying the error? Taking the average pressure within the non-PML region after the wave enters the PML? Ideally the pressure in the non-PML region would be zero, but it seems this is not possible based on the docs from k-wave.

[luraess]

Thanks @alexpattyn for resuming this. Good question regarding the error quantification. I guess the real accurate approach may be to quantify the energy loss through the PML and ensure that most of the energy is indeed absorbed in the PML-region. I don't know however if this level of complexity is needed here. It may indeed be sufficient to proceed as you suggest

Taking the average pressure within the non-PML region after the wave enters the PML

You could maybe try the above one for various cases (including optimal and less optimal absorbing parameters) to get an idea about how relevant this approach could be. I guess it would make sense the pressure error to be smallest for optimal PML absorbing params.

[alexpattyn]

Reviewing k-waves docs, they quantified the performance of the PML by looking at the change in amplitude of the transmitted and reflected wave. From my earlier comments \alpha_max = A c/dx; where A = maximum attenuation, c = sound speed, and dx = spatial discretization. When A = 4 I am seeing ~37dB drop in the reflected signal, where dB drop = 20*log10(P/P0) | P = max reflected amplitude and P0 = max amplitude.

I'll continue to look into this.

[alexpattyn]

I saw a paper from Yuan et al where to compare the effects of a PML verses Mur's second order boundary condition. They'd defined some error function that looks at the relative difference between the acoustic field with a boundary condition and one without (where the forward model is stopped before reflection occur).

Which I implemented as:

Error += @. 10*log10((Pₜᵣᵤₑ - Pₚₘₗ)^2/Pₜᵣᵤₑ^2)

Given this I got the results below:

So I think we are good with the current PML implementation. I can add a comment in the code that the user should fine-tune the PML parameters for their particular applications.

Unless there are other comments I can start to prepare a PR.

[luraess]

Hi @alexpattyn, thanks for your new addition. I'd say your latest error analysis makes sense and serves the purpose. One thing I was wondering is why the setup is not fully symmetric? Do I miss something?

Unless there are other comments I can start to prepare a PR.

That sounds reasonable to me, unless @omlins has something else to add prior to that. Note however that ParallelStencil is undergoing some refactoring, so I there are two options:

• prepare the PR now and maybe modify some bits soonish
• wait until the refactored version is ready for release and add your contribution in the flow

In all cases, thanks for your contribution!

[alexpattyn]

Hey @luraess,

I'd say your latest error analysis makes sense and serves the purpose. One thing I was wondering is why the setup is not fully symmetric? Do I miss something?

You are referring to the "purple" diagonal in the error heatmap? That was something I was wondering as well. I'll try to investigate that a bit further.

Considering there is a refactoring going on I'll hold off on the PR for now. Additionally, I would also like to get an ultrasound simulation going with realistic physical properties (e.g. domain of 220 mm x 220 mm), however I am getting stability issues that I need to look at.

But no problem! It is helpful for my own research project. I just wish I had more time to spend on this.

But will the refactor include AMD GPU support? 😀

[luraess]

You are referring to the "purple" diagonal in the error heatmap

Yes, it would be good to see that if the entire configuration is symmetric, then an asymmetry may be a hint for a minor onset somewhere.

But will the refactor include AMD GPU support

We are on it, yes 😄

[alexpattyn]

... may be a hint for a minor onset somewhere.

What do you mean here? Like stability issues or a model mismatch?

We are on it, yes smile

Fantastic! I have a Radeon VII ready to be benchmarked when the time comes. 😄

[luraess]

Like stability issues

Not really, rather something like the initial perturbation not being fully centred, or some indexing bug - something maybe super minor that would make the configuration non-symmetric.

I have a Radeon VII ready

Super!