add fully mapped version of cic interpolation

AMReX-Codes · Jan 4, 2024 · ba7c8d6 · ba7c8d6
1 parent a032bd3
commit ba7c8d6
Showing 1 changed file with 37 additions and 32 deletions.
diff --git 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 <typename P>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-void cic_interpolate_mapped_z (const P& p,
-                               amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
-                               amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
-                               const amrex::Array4<amrex::Real const>&  uccarr,
-                               const amrex::Array4<amrex::Real const>&  zheight,
-                               bool use_terrain, amrex::ParticleReal* val, int M = AMREX_SPACEDIM)
+void cic_interpolate_mapped (const P& p,
+                             amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
+                             amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
+                             const amrex::Array4<amrex::Real const>&  ucc_arr,
+                             const amrex::Array4<amrex::Real const>&  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<int>(amrex::Math::floor(lx));
-    int const j = static_cast<int>(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<Real>(i);
-    amrex::Real const yint = ly - static_cast<Real>(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<ParticleReal>(sx[ii]*sy[jj]*sz[kk]*uccarr(i+ii,j+jj,k+kk),d);
-                }
-            }
-        }
+                    val[d] += static_cast<ParticleReal>(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(},},});
     }
 }