always use bathymetry when defining z_phys_nd and particles

seahorce-scidac · Dec 15, 2023 · 696fbb3 · 696fbb3
1 parent f27bdf2
commit 696fbb3
Show file tree

Hide file tree

Showing 5 changed files with 16 additions and 83 deletions.
diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp
@@ -198,7 +198,6 @@ ROMSX::WritePlotFile (int which, Vector<std::string> plot_var_names)
         {
             MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 0);
             temp_dat.setVal(0);
-            particleData.tracer_particles->Redistribute();
             particleData.tracer_particles->Increment(temp_dat, lev);
             MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0);
             mf_comp += 1;

diff --git a/Source/Particles/ParticleData.H b/Source/Particles/ParticleData.H
@@ -32,12 +32,7 @@ struct ParticleData {
         // Initialize tracer particles
         if (use_tracer_particles) {
             tracer_particles = std::make_unique<TracerPC>(a_gdb);
-            if (z_phys_nd[0] != nullptr) {
-                tracer_particles->InitParticles(*z_phys_nd[0]);
-            } else {
-                tracer_particles->InitParticles();
-            }
-            tracer_particles->Redistribute();
+            tracer_particles->InitParticles(*z_phys_nd[0]);
             amrex::Print() << "Initialized " << tracer_particles->TotalNumberOfParticles() << " tracer particles." << std::endl;
         }
     }

diff --git a/Source/Particles/TracerPC.H b/Source/Particles/TracerPC.H
@@ -59,7 +59,6 @@ public:
                                    amrex::DefaultAllocator, TracerAssignor>(geom, dmap, ba)
         {}
 
-    void InitParticles ();
     void InitParticles (const amrex::MultiFab& a_z_height);
 
     void AdvectWithUmac (amrex::Array<amrex::MultiFab const*, AMREX_SPACEDIM> umac,

diff --git a/Source/Particles/TracerPC.cpp b/Source/Particles/TracerPC.cpp
@@ -4,58 +4,6 @@
 
 using namespace amrex;
 
-void
-TracerPC::
-InitParticles ()
-{
-    BL_PROFILE("TracerPC::InitParticles");
-
-    const int lev = 0;
-    const Real* dx = Geom(lev).CellSize();
-    const Real* plo = Geom(lev).ProbLo();
-
-    for(MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi)
-    {
-        const Box& tile_box  = mfi.tilebox();
-        Gpu::HostVector<ParticleType> host_particles;
-        for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) {
-            if (iv[0] == 3) {
-                Real r[3] = {0.5, 0.5, 0.5};  // this means place at cell center
-
-                Real x = plo[0] + (iv[0] + r[0])*dx[0];
-                Real y = plo[1] + (iv[1] + r[1])*dx[1];
-                Real z = plo[2] + (iv[2] + r[2])*dx[2];
-
-                ParticleType p;
-                p.id()  = ParticleType::NextID();
-                p.cpu() = ParallelDescriptor::MyProc();
-                p.pos(0) = x;
-                p.pos(1) = y;
-                p.pos(2) = z;
-
-                p.rdata(TracerRealIdx::old_x) = p.pos(0);
-                p.rdata(TracerRealIdx::old_y) = p.pos(1);
-                p.rdata(TracerRealIdx::old_z) = p.pos(2);
-
-                p.idata(TracerIntIdx::k) = iv[2];  // particles carry their z-index
-
-                host_particles.push_back(p);
-           }
-        }
-
-        auto& particles = GetParticles(lev);
-        auto& particle_tile = particles[std::make_pair(mfi.index(), mfi.LocalTileIndex())];
-        auto old_size = particle_tile.GetArrayOfStructs().size();
-        auto new_size = old_size + host_particles.size();
-        particle_tile.resize(new_size);
-
-        Gpu::copy(Gpu::hostToDevice,
-                  host_particles.begin(),
-                  host_particles.end(),
-                  particle_tile.GetArrayOfStructs().begin() + old_size);
-    }
-}
-
 void
 TracerPC::
 InitParticles (const MultiFab& a_z_height)

diff --git a/Source/Utils/DepthStretchTransform.H b/Source/Utils/DepthStretchTransform.H
@@ -212,30 +212,22 @@ ROMSX::stretch_transform (int lev)
       // WARNING: z_w(i,j,k)       refers to the face on the HIGH side of cell (i,j,k)
       // WARNING: z_phys_nd(i,j,k) refers to the node on the LOW  side of cell (i,j,k)
       //
-      if (solverChoice.flat_bathymetry) {
-          Real dz = geom[lev].CellSize(2);
-          ParallelFor(Box(z_phys_nd), [=] AMREX_GPU_DEVICE (int i, int j, int k)
-          {
-             z_phys_nd(i,j,k) = k * dz;
-          });
-      } else {
-          ParallelFor(Box(z_phys_nd), [=] AMREX_GPU_DEVICE (int i, int j, int k)
+      ParallelFor(Box(z_phys_nd), [=] AMREX_GPU_DEVICE (int i, int j, int k)
+      {
+          // For now assume all boundaries are constant height --
+          //     we will enforce periodicity below
+          int kk = (k == 0) ? hi.z : k-1;
+          if ( i >= lo.x && i <= hi.x-1 && j >= lo.y && j <= hi.y-1 )
           {
-              // For now assume all boundaries are constant height --
-              //     we will enforce periodicity below
-              int kk = (k == 0) ? hi.z : k-1;
-              if ( i >= lo.x && i <= hi.x-1 && j >= lo.y && j <= hi.y-1 )
-              {
-                  z_phys_nd(i,j,k)=0.25*( z_w(i,j  ,kk) + z_w(i+1,j  ,kk) +
-                                          z_w(i,j+1,kk) + z_w(i+1,j+1,kk) );
-              } else {
-                 int ii = std::min(std::max(i, lo.x), hi.x);
-                 int jj = std::min(std::max(j, lo.y), hi.y);
-                 z_phys_nd(i,j,k) = z_w(ii,jj,kk);
-              }
-              if (k == 0)  z_phys_nd(i,j,k) *= -1.;
-          });
-      } // if not flat
+              z_phys_nd(i,j,k)=0.25*( z_w(i,j  ,kk) + z_w(i+1,j  ,kk) +
+                                      z_w(i,j+1,kk) + z_w(i+1,j+1,kk) );
+          } else {
+             int ii = std::min(std::max(i, lo.x), hi.x);
+             int jj = std::min(std::max(j, lo.y), hi.y);
+             z_phys_nd(i,j,k) = z_w(ii,jj,kk);
+          }
+          if (k == 0)  z_phys_nd(i,j,k) *= -1.;
+      });
     } // mf
 
     // Note that we do *not* want to do a multilevel fill here -- we have