diff --git a/Src/Particle/AMReX_TracerParticle_mod_K.H b/Src/Particle/AMReX_TracerParticle_mod_K.H index 126f980a1a2..5b1262eab08 100644 --- a/Src/Particle/AMReX_TracerParticle_mod_K.H +++ b/Src/Particle/AMReX_TracerParticle_mod_K.H @@ -58,54 +58,59 @@ void cic_interpolate (const P& p, } } + +// AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +// int compute_particle_weights (amrex::Real* s, amrex::Real x) { +// amrex::Real l = (x - plo[0]) * dxi[0] - Real(0.5); + +// } + template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void cic_interpolate_mapped_z (const P& p, - amrex::GpuArray const& plo, - amrex::GpuArray const& dxi, - const amrex::Array4& uccarr, - const amrex::Array4& zheight, - bool use_terrain, amrex::ParticleReal* val, int M = AMREX_SPACEDIM) +void cic_interpolate_mapped (const P& p, + amrex::GpuArray const& plo, + amrex::GpuArray const& dxi, + const amrex::Array4& ucc_arr, + const amrex::Array4& loc_arr, + amrex::ParticleReal* val, int M = AMREX_SPACEDIM) { AMREX_ASSERT(val != nullptr); - // handle x and y normally - amrex::Real lx = (Real(p.pos(0)) - plo[0]) * dxi[0] - Real(0.5); - amrex::Real ly = (Real(p.pos(1)) - plo[1]) * dxi[1] - Real(0.5); + AMREX_D_TERM(int i = p.idata(0);, + int j = p.idata(1);, + int k = p.idata(2)); - int const i = static_cast(amrex::Math::floor(lx)); - int const j = static_cast(amrex::Math::floor(ly)); + AMREX_D_TERM(amrex::Real lx = amrex::Real(p.pos(0)) - 0.5*(loc_arr(AMREX_D_DECL(i-1,j,k),0) + loc_arr(AMREX_D_DECL(i,j,k),0));, + amrex::Real ly = amrex::Real(p.pos(1)) - 0.5*(loc_arr(AMREX_D_DECL(i,j-1,k),1) + loc_arr(AMREX_D_DECL(i,j,k),1));, + amrex::Real lz = amrex::Real(p.pos(2)) - 0.5*(loc_arr(AMREX_D_DECL(i,j,k-1),2) + loc_arr(AMREX_D_DECL(i,j,k),2))); - amrex::Real const xint = lx - static_cast(i); - amrex::Real const yint = ly - static_cast(j); + AMREX_D_TERM(int i0 = i + amrex::Math::floor(lx/(loc_arr(AMREX_D_DECL(i+1,j,k),0) - loc_arr(AMREX_D_DECL(i,j,k),0)));, + int j0 = j + amrex::Math::floor(ly/(loc_arr(AMREX_D_DECL(i,j+1,k),1) - loc_arr(AMREX_D_DECL(i,j,k),1)));, + int k0 = k + amrex::Math::floor(lz/(loc_arr(AMREX_D_DECL(i,j,k+1),2) - loc_arr(AMREX_D_DECL(i,j,k),2)))); - amrex::Real sy[] = {Real(1.0) - yint, yint}; - amrex::Real sx[] = {Real(1.0) - xint, xint}; + AMREX_D_TERM(amrex::Real const xint = 2.0*lx/(loc_arr(AMREX_D_DECL(i0+2,j,k),0) - loc_arr(AMREX_D_DECL(i0,j,k),0));, + amrex::Real const yint = 2.0*ly/(loc_arr(AMREX_D_DECL(i,j0+2,k),1) - loc_arr(AMREX_D_DECL(i,j0,k),1));, + amrex::Real const zint = 2.0*lz/(loc_arr(AMREX_D_DECL(i,j,k0+2),2) - loc_arr(AMREX_D_DECL(i,j,k0),2))); - // handle mapped z if in use - int const k = p.idata(0); - ParticleReal zlo, zhi; - if (use_terrain) { - zlo = zheight(i, j, k ); - zhi = zheight(i, j, k+1); - } else { - zlo = k / dxi[2]; - zhi = (k+1) / dxi[2]; - } - - amrex::Real lz = Real((p.pos(2) - zlo) / (zhi - zlo)) - Real(0.5); - amrex::Real const zint = lz; +#if AMREX_SPACEDIM > 2 amrex::Real sz[] = {Real(1.0) - zint, zint}; +#endif +#if AMREX_SPACEDIM > 1 + amrex::Real sy[] = {Real(1.0) - yint, yint}; +#endif + amrex::Real sx[] = {Real(1.0) - xint, xint}; for (int d=0; d < M; ++d) { val[d] = ParticleReal(0.0); +#if AMREX_SPACEDIM > 2 for (int kk = 0; kk<=1; ++kk) { +#endif +#if AMREX_SPACEDIM > 1 for (int jj = 0; jj <= 1; ++jj) { +#endif for (int ii = 0; ii <= 1; ++ii) { - val[d] += static_cast(sx[ii]*sy[jj]*sz[kk]*uccarr(i+ii,j+jj,k+kk),d); - } - } - } + val[d] += static_cast(AMREX_D_TERM(sx[ii],*sy[jj],*sz[kk])*ucc_arr(IntVect(AMREX_D_DECL(i0+ii,j0+jj,k0+kk)),d)); + AMREX_D_TERM(},},}); } }