Skip to content

Commit

Permalink
added update_location_idata() function to compute particle cell index…
Browse files Browse the repository at this point in the history 
… from position
  • Loading branch information
@debog
debog committed Jan 25, 2024
1 parent 7687f0e commit c109aca
Show file tree
Hide file tree
Showing 2 changed files with 37 additions and 56 deletions.
32 changes: 32 additions & 0 deletions Source/Particles/ERFPC.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,38 @@ struct ERFParticlesAssignor
}
};

template <typename P>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void update_location_idata (P& p,
amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
const amrex::Array4<amrex::Real const>& height_arr)
{
amrex::IntVect iv( int(amrex::Math::floor((p.pos(0)-plo[0])*dxi[0])),
int(amrex::Math::floor((p.pos(1)-plo[1])*dxi[1])),
p.idata(ERFParticlesIntIdx::k) );

p.idata(ERFParticlesIntIdx::i) = iv[0];
p.idata(ERFParticlesIntIdx::j) = iv[1];

amrex::Real lx = (p.pos(0)-plo[0])*dxi[0] - static_cast<amrex::Real>(iv[0]);
amrex::Real ly = (p.pos(1)-plo[1])*dxi[1] - static_cast<amrex::Real>(iv[1]);
auto zlo = height_arr(iv[0] ,iv[1] ,iv[2] ) * (1.0-lx) * (1.0-ly) +
height_arr(iv[0]+1,iv[1] ,iv[2] ) * lx * (1.0-ly) +
height_arr(iv[0] ,iv[1]+1,iv[2] ) * (1.0-lx) * ly +
height_arr(iv[0]+1,iv[1]+1,iv[2] ) * lx * ly;
auto zhi = height_arr(iv[0] ,iv[1] ,iv[2]+1) * (1.0-lx) * (1.0-ly) +
height_arr(iv[0]+1,iv[1] ,iv[2]+1) * lx * (1.0-ly) +
height_arr(iv[0] ,iv[1]+1,iv[2]+1) * (1.0-lx) * ly +
height_arr(iv[0]+1,iv[1]+1,iv[2]+1) * lx * ly;

if (p.pos(2) > zhi) {
p.idata(ERFParticlesIntIdx::k) += 1;
} else if (p.pos(2) <= zlo) {
p.idata(ERFParticlesIntIdx::k) -= 1;
}
}

class ERFPC : public amrex::ParticleContainer< ERFParticlesRealIdx::ncomps,
ERFParticlesIntIdx::ncomps,
0,
Expand Down
61 changes: 5 additions & 56 deletions Source/Particles/ERFPCEvolve.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,14 +103,16 @@ void ERFPC::AdvectWithFlow ( MultiFab* a_umac,
ParticleType& p = p_pbox[i];
if (p.id() <= 0) { return; }
ParticleReal v[AMREX_SPACEDIM];
mac_interpolate(p, plo, dxi, umacarr, v);
mac_interpolate_mapped_z(p, plo, dxi, umacarr, zheight, v);
if (ipass == 0)
{
for (int dim=0; dim < AMREX_SPACEDIM; dim++)
{
p.rdata(dim) = p.pos(dim);
p.pos(dim) += static_cast<ParticleReal>(ParticleReal(0.5)*a_dt*v[dim]);
}
// Update z-coordinate carried by the particle
update_location_idata(p,plo,dxi,zheight);
}
else
{
Expand All @@ -119,35 +121,8 @@ void ERFPC::AdvectWithFlow ( MultiFab* a_umac,
p.pos(dim) = p.rdata(dim) + static_cast<ParticleReal>(a_dt*v[dim]);
p.rdata(dim) = v[dim];
}

// 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])),
p.idata(ERFParticlesIntIdx::k)));
iv[0] += domain.smallEnd()[0];
iv[1] += domain.smallEnd()[1];
ParticleReal zlo, zhi;
if (use_terrain) {
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];
}
if (p.pos(2) > zhi) { // need to be careful here
p.idata(ERFParticlesIntIdx::k) += 1;
} else if (p.pos(2) <= zlo) {
p.idata(ERFParticlesIntIdx::k) -= 1;
}
update_location_idata(p,plo,dxi,zheight);
}
});
}
Expand Down Expand Up @@ -223,33 +198,7 @@ void ERFPC::AdvectWithGravity ( int a_lev,
p.rdata(ERFParticlesRealIdx::vz) -= (grav - drag) * half_dt;

// also update z-coordinate here
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])),
p.idata(ERFParticlesIntIdx::k)));
iv[0] += domain.smallEnd()[0];
iv[1] += domain.smallEnd()[1];
ParticleReal zlo, zhi;
if (use_terrain) {
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];
}
if (p.pos(2) > zhi) { // need to be careful here
p.idata(ERFParticlesIntIdx::k) += 1;
} else if (p.pos(2) <= zlo) {
p.idata(ERFParticlesIntIdx::k) -= 1;
}
update_location_idata(p,plo,dxi,zheight);
});
}

Expand Down

0 comments on commit c109aca

Please sign in to comment.