[WIP] Allow new levels on ParticleContainer restart

Src/Base/AMReX_PlotFileUtil.cpp

@@ -1,4 +1,5 @@
 
+#include "AMReX_ParallelDescriptor.H"
 #include <AMReX_VisMF.H>
 #include <AMReX_AsyncOut.H>
 #include <AMReX_PlotFileUtil.H>
@@ -236,6 +237,7 @@ WriteMultiLevelPlotfile (const std::string& plotfilename, int nlevels,
             }
             VisMF::Write(*data, MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mfPrefix));
         }
+        ParallelDescriptor::Barrier();
     }
 }
 

Src/Base/AMReX_VisMF.cpp

@@ -1111,7 +1111,7 @@ VisMF::Write (const FabArray<FArrayBox>&    mf,
             nfi.Stream().flush();
             delete [] allFabData;
 
-        } else {    // ---- write fabs individually
+          } else {    // ---- write fabs individually
             for(MFIter mfi(mf); mfi.isValid(); ++mfi) {
                 std::streamoff hLength = 0;
                 const FArrayBox &fab = mf[mfi];
@@ -1123,10 +1123,9 @@ VisMF::Write (const FabArray<FArrayBox>&    mf,
                     hLength = static_cast<std::streamoff>(hss.tellp());
                     auto tstr = hss.str();
                     nfi.Stream().write(tstr.c_str(), hLength);    // ---- the fab header
-                    nfi.Stream().flush();
                 }
                 Real const* fabdata = fab.dataPtr();
-#ifdef AMREX_USE_GPU
+        #ifdef AMREX_USE_GPU
                 std::unique_ptr<FArrayBox> hostfab;
                 if (fab.arena()->isManaged() || fab.arena()->isDevice()) {
                     hostfab = std::make_unique<FArrayBox>(fab.box(), fab.nComp(),
@@ -1136,20 +1135,19 @@ VisMF::Write (const FabArray<FArrayBox>&    mf,
                     Gpu::streamSynchronize();
                     fabdata = hostfab->dataPtr();
                 }
-#endif
+        #endif
                 if(doConvert) {
                     char *cDataPtr = new char[writeDataSize];
                     RealDescriptor::convertFromNativeFormat(static_cast<void *> (cDataPtr),
                                                             writeDataItems,
                                                             fabdata, *whichRD);
                     nfi.Stream().write(cDataPtr, writeDataSize);
-                    nfi.Stream().flush();
                     delete [] cDataPtr;
                 } else {    // ---- copy from the fab
                     nfi.Stream().write((char *) fabdata, writeDataSize);
-                    nfi.Stream().flush();
                 }
             }
+            nfi.Stream().flush();
         }
     }
 

Src/Particle/AMReX_ParticleIO.H

@@ -613,10 +613,10 @@ ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssig
                                     particle_io_flags, tile_map[grid], count[grid], is_checkpoint);
 
         writeIntData(istuff.dataPtr(), istuff.size(), ofs);
-        ofs.flush();  // Some systems require this flush() (probably due to a bug)
+        //ofs.flush();  // Some systems require this flush() (probably due to a bug)
 
         WriteParticleRealData(rstuff.dataPtr(), rstuff.size(), ofs);
-        ofs.flush();  // Some systems require this flush() (probably due to a bug)
+        //ofs.flush();  // Some systems require this flush() (probably due to a bug)
     }
 }
 
@@ -746,6 +746,9 @@ ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssig
     HdrFile >> finest_level_in_file;
     AMREX_ASSERT(finest_level_in_file >= 0);
 
+    // handle the case finestLevel() > finest_level_in_file
+    //int finest_level_overall = std::max(finest_level_in_file, finestLevel());
+
     // Determine whether this is a dual-grid restart or not.
     Vector<DistributionMapping> old_dms(finest_level_in_file + 1);
     Vector<BoxArray> old_bas(finest_level_in_file + 1);
@@ -768,9 +771,9 @@ ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssig
 
     if (have_pheaders)
     {
-        for (int lev = 0; lev <= finestLevel(); lev++)
+        for (int lev = 0; lev <= finest_level_in_file; lev++)
         {
-            old_dms[lev] = ParticleDistributionMap(lev);
+            old_dms[lev] = ParticleDistributionMap(lev); // FIXME: segfault here!!!
             old_bas[lev] = ParticleBoxArray(lev);
             std::string phdr_name = fullname;
             phdr_name += "/Level_";
@@ -793,7 +796,7 @@ ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssig
     }
 
     if (dual_grid) {
-        for (int lev = 0; lev <= finestLevel(); lev++) {
+        for (int lev = 0; lev <= finest_level_in_file; lev++) {
             // this can happen if there are no particles at a given level in the checkpoint
             if (particle_box_arrays[lev].empty()) {
                 particle_box_arrays[lev] = BoxArray(Geom(lev).Domain());
@@ -817,6 +820,7 @@ ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssig
     }
 
     for (int lev = 0; lev <= finest_level_in_file; lev++) {
+        amrex::Print() << "ParticleContainer::Restart(): Reading level " << lev << std::endl;
         Vector<int>  which(ngrids[lev]);
         Vector<int>  count(ngrids[lev]);
         Vector<Long> where(ngrids[lev]);
@@ -914,12 +918,19 @@ ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssig
     }
 
     if (dual_grid) {
-        for (int lev = 0; lev <= finestLevel(); lev++) {
+        for (int lev = 0; lev <= finest_level_in_file; lev++) {
             SetParticleBoxArray(lev, old_bas[lev]);
             SetParticleDistributionMap(lev, old_dms[lev]);
         }
     }
 
+    // Clear particles for levels higher than finest_level_in_file but present in current container
+    for (int lev = finest_level_in_file + 1; lev <= finestLevel(); ++lev) {
+        if (lev < m_particles.size()) {
+            m_particles[lev].clear();
+        }
+    }
+
     Redistribute();
 
     AMREX_ASSERT(OK());

Src/Particle/AMReX_WriteBinaryParticleData.H

@@ -673,6 +673,7 @@ void WriteBinaryParticleDataSync (PC const& pc,
                 }
             }
         }
+        ParallelDescriptor::Barrier();
     }
 
     if (ParallelDescriptor::IOProcessor())