Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix oops #1383

Merged
merged 2 commits into from
Jan 16, 2024
Merged

Fix oops #1383

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
865720d
fix how we locate particles with terrain
asalmgren Jan 16, 2024
daf595f
forgot to fix the other particle advection routine
asalmgren Jan 16, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 22 additions & 11 deletions Source/Particles/ERFPCEvolve.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,8 +98,7 @@ void ERFPC::AdvectWithFlow ( MultiFab* a_umac,
bool use_terrain = (a_z_height != nullptr);
auto zheight = use_terrain ? (*a_z_height)[grid].array() : Array4<Real>{};
Copy link
Contributor

@debog debog Jan 20, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@asalmgren @AMLattanzi how does this make zheight available on the GPU?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is due to C++ lambda capture. If a C++ lambda uses [=], it will make a copy of all the variables used in the lambda's function body. These variables are part of the lambda function object. In the case of GPU, the CUDA runtime will copy the lambda function to GPU at kernel launch.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, @WeiqunZhang! Makes sense now.


ParallelFor(n,
[=] AMREX_GPU_DEVICE (int i)
ParallelFor(n, [=] AMREX_GPU_DEVICE (int i)
{
ParticleType& p = p_pbox[i];
if (p.id() <= 0) { return; }
Expand All @@ -121,7 +120,7 @@ void ERFPC::AdvectWithFlow ( MultiFab* a_umac,
p.rdata(dim) = v[dim];
}

// also update z-coordinate here
// Update z-coordinate carried by the particle
IntVect iv(
AMREX_D_DECL(int(amrex::Math::floor((p.pos(0)-plo[0])*dxi[0])),
int(amrex::Math::floor((p.pos(1)-plo[1])*dxi[1])),
Expand All @@ -130,10 +129,16 @@ void ERFPC::AdvectWithFlow ( MultiFab* a_umac,
iv[1] += domain.smallEnd()[1];
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@AMLattanzi @asalmgren Can lines 124-129 be replaced by ERFParticlesAssignor(p, plo, dxi, domain) (defined in ERFPC.H)? Same with lines 226-231 below?

ParticleReal zlo, zhi;
if (use_terrain) {
zlo = Real(0.25) * (zheight(iv[0],iv[1] ,iv[2] ) + zheight(iv[0]+1,iv[1] ,iv[2] ) +
zheight(iv[0],iv[1]+1,iv[2] ) + zheight(iv[0]+1,iv[1]+1,iv[2] ));
zhi = Real(0.25) * (zheight(iv[0],iv[1] ,iv[2]+1) + zheight(iv[0]+1,iv[1] ,iv[2]+1) +
zheight(iv[0],iv[1]+1,iv[2]+1) + zheight(iv[0]+1,iv[1]+1,iv[2]+1));
Real lx = (p.pos(0)-plo[0])*dxi[0] - static_cast<Real>(iv[0]-domain.smallEnd()[0]);
Real ly = (p.pos(1)-plo[1])*dxi[1] - static_cast<Real>(iv[1]-domain.smallEnd()[1]);
zlo = zheight(iv[0] ,iv[1] ,iv[2] ) * (1.0-lx) * (1.0-ly) +
zheight(iv[0]+1,iv[1] ,iv[2] ) * lx * (1.0-ly) +
zheight(iv[0] ,iv[1]+1,iv[2] ) * (1.0-lx) * ly +
zheight(iv[0]+1,iv[1]+1,iv[2] ) * lx * ly;
zhi = zheight(iv[0] ,iv[1] ,iv[2]+1) * (1.0-lx) * (1.0-ly) +
zheight(iv[0]+1,iv[1] ,iv[2]+1) * lx * (1.0-ly) +
zheight(iv[0] ,iv[1]+1,iv[2]+1) * (1.0-lx) * ly +
zheight(iv[0]+1,iv[1]+1,iv[2]+1) * lx * ly;
} else {
zlo = iv[2] * dx[2];
zhi = (iv[2]+1) * dx[2];
Comment on lines 130 to 144
Copy link
Contributor

@debog debog Jan 20, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@AMLattanzi @asalmgren Can lines 130-145 (or 130-151) be written as a separate function/functor? If functor, it could store zheight and use_terrain as member data? And we could use it for lines 232-252 as well.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually I want to move this functionality into a routine in amrex -- it is used by multiple codes. I was just going to make a new function called something like update_particle_idata that takes one particle and the zheight array. What do you think?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That would be nice since this is a particle functionality, not an ERF-specific thing. But terrain and zheight are specific to atmospheric flows. To make it more general, could we have xcoord, ycoord, zcoord (defaulted to nullptr?) instead of just zheight that allows mapping in all dimensions? ERF would then call this function with only zheight.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We're working on the generalized functionality ...

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we could parallel what is done in the BinMapper functor. This would move the functionality to AMReX, clean up the ERF code since the above snippet now lives in a functor, and would be similar to what we do without curvilinear coordinates.

Expand Down Expand Up @@ -226,10 +231,16 @@ void ERFPC::AdvectWithGravity ( int a_lev,
iv[1] += domain.smallEnd()[1];
ParticleReal zlo, zhi;
if (use_terrain) {
zlo = Real(0.25) * (zheight(iv[0],iv[1] ,iv[2] ) + zheight(iv[0]+1,iv[1] ,iv[2] ) +
zheight(iv[0],iv[1]+1,iv[2] ) + zheight(iv[0]+1,iv[1]+1,iv[2] ));
zhi = Real(0.25) * (zheight(iv[0],iv[1] ,iv[2]+1) + zheight(iv[0]+1,iv[1] ,iv[2]+1) +
zheight(iv[0],iv[1]+1,iv[2]+1) + zheight(iv[0]+1,iv[1]+1,iv[2]+1));
Real lx = (p.pos(0)-plo[0])*dxi[0] - static_cast<Real>(iv[0]-domain.smallEnd()[0]);
Real ly = (p.pos(1)-plo[1])*dxi[1] - static_cast<Real>(iv[1]-domain.smallEnd()[1]);
zlo = zheight(iv[0] ,iv[1] ,iv[2] ) * (1.0-lx) * (1.0-ly) +
zheight(iv[0]+1,iv[1] ,iv[2] ) * lx * (1.0-ly) +
zheight(iv[0] ,iv[1]+1,iv[2] ) * (1.0-lx) * ly +
zheight(iv[0]+1,iv[1]+1,iv[2] ) * lx * ly;
zhi = zheight(iv[0] ,iv[1] ,iv[2]+1) * (1.0-lx) * (1.0-ly) +
zheight(iv[0]+1,iv[1] ,iv[2]+1) * lx * (1.0-ly) +
zheight(iv[0] ,iv[1]+1,iv[2]+1) * (1.0-lx) * ly +
zheight(iv[0]+1,iv[1]+1,iv[2]+1) * lx * ly;
} else {
zlo = iv[2] * dx[2];
zhi = (iv[2]+1) * dx[2];
Expand Down
Loading