From 865720d6e745b92379259c751d60847234581509 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Tue, 16 Jan 2024 10:45:26 -0800 Subject: [PATCH 1/2] fix how we locate particles with terrain --- Source/Particles/ERFPCEvolve.cpp | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/Source/Particles/ERFPCEvolve.cpp b/Source/Particles/ERFPCEvolve.cpp index dc7e7a941..0aed01002 100644 --- a/Source/Particles/ERFPCEvolve.cpp +++ b/Source/Particles/ERFPCEvolve.cpp @@ -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{}; - ParallelFor(n, - [=] AMREX_GPU_DEVICE (int i) + ParallelFor(n, [=] AMREX_GPU_DEVICE (int i) { ParticleType& p = p_pbox[i]; if (p.id() <= 0) { return; } @@ -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])), @@ -130,10 +129,16 @@ void ERFPC::AdvectWithFlow ( MultiFab* a_umac, 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(iv[0]-domain.smallEnd()[0]); + Real ly = (p.pos(1)-plo[1])*dxi[1] - static_cast(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]; From daf595feb69698aa09003dd37a98c163bc256276 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Tue, 16 Jan 2024 11:10:37 -0800 Subject: [PATCH 2/2] forgot to fix the other particle advection routine --- Source/Particles/ERFPCEvolve.cpp | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/Source/Particles/ERFPCEvolve.cpp b/Source/Particles/ERFPCEvolve.cpp index 0aed01002..47a0efd27 100644 --- a/Source/Particles/ERFPCEvolve.cpp +++ b/Source/Particles/ERFPCEvolve.cpp @@ -231,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(iv[0]-domain.smallEnd()[0]); + Real ly = (p.pos(1)-plo[1])*dxi[1] - static_cast(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];