Comments (19)
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. :)
from parallelstencil.jl.
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.
from parallelstencil.jl.
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).
from parallelstencil.jl.
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 🙂
from parallelstencil.jl.
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.
from parallelstencil.jl.
Sounds good! I would think two good examples from my field would be:
- 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.
- 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!
from parallelstencil.jl.
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.
from parallelstencil.jl.
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.
from parallelstencil.jl.
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 🙂
from parallelstencil.jl.
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.
from parallelstencil.jl.
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.
from parallelstencil.jl.
@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.
from parallelstencil.jl.
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.
from parallelstencil.jl.
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.
from parallelstencil.jl.
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!
from parallelstencil.jl.
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? 😀
from parallelstencil.jl.
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 😄
from parallelstencil.jl.
... 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. 😄
from parallelstencil.jl.
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!
from parallelstencil.jl.
Related Issues (20)
- AMDGPU v0.5.0 compat HOT 1
- Add device_sync
- sync issues on AMDGPU backend
- Make CellArrays mutable HOT 4
- finite volume method HOT 3
- [JuliaCon/proceedings-review] @parallel keyword argument `loopopt` deprecated? HOT 1
- ParallelStencil on 1.10 HOT 6
- [JuliaCon/proceedings-review] DOI of paper by Besard et al. HOT 2
- [JuliaCon/proceedings-review] Community guidelines HOT 1
- [JuliaCon/proceedings-review] Performance metrics HOT 4
- Type unstable Data.Number HOT 2
- GPU memory management issue when running multi-GPU code HOT 10
- Add support for Polyester's `@batch` HOT 20
- Generalize loopopt
- Create and update GPU unit tests
- Thread (CPU) Float32/Float64 performance comparison on miniapp acoustic2D HOT 12
- Example for init_global_grid_usage HOT 3
- How to implement custom finite differencing operators HOT 8
- CUDA Crash with julia 1.9.0 HOT 8
- Non cartesian gather! HOT 2
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from parallelstencil.jl.