Bugfix for particle advection without terrain

LLNL · Jan 31, 2024 · c62d827 · c62d827
1 parent 4e88708
commit c62d827
Show file tree

Hide file tree

Showing 2 changed files with 19 additions and 16 deletions.
diff --git a/Source/Particles/ERFPC.H b/Source/Particles/ERFPC.H
@@ -74,21 +74,23 @@ void update_location_idata (P& p,
     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;
+    if (height_arr) {
+        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;
+        }
     }
 }
 

diff --git a/Source/Particles/ERFPCEvolve.cpp b/Source/Particles/ERFPCEvolve.cpp
@@ -103,7 +103,8 @@ void ERFPC::AdvectWithFlow ( MultiFab*                           a_umac,
                 ParticleType& p = p_pbox[i];
                 if (p.id() <= 0) { return; }
                 ParticleReal v[AMREX_SPACEDIM];
-                mac_interpolate_mapped_z(p, plo, dxi, umacarr, zheight, v);
+                if (use_terrain) mac_interpolate_mapped_z(p, plo, dxi, umacarr, zheight, v);
+                else             mac_interpolate(p, plo, dxi, umacarr, v);
                 if (ipass == 0)
                 {
                     for (int dim=0; dim < AMREX_SPACEDIM; dim++)