Mesh refinement: Adding IO

src/Hipace.H

@@ -290,6 +290,7 @@ public:
      */
     void CheckGhostSlice (int it);
 
+    static int m_nlev; /**< number of MR levels */
     amrex::RealVect patch_lo {0., 0., 0.}; /**< 3D array with lower ends of the refined grid */
     amrex::RealVect patch_hi {0., 0., 0.}; /**< 3D array with upper ends of the refined grid */
 private:
@@ -313,7 +314,7 @@ private:
     Diagnostic m_diags;
 
     void ResizeFDiagFAB (const amrex::Box box, const int lev) { m_diags.ResizeFDiagFAB(box, lev); };
-    void FillDiagnostics (int lev, int i_slice);
+    void FillDiagnostics (int i_slice);
     /** \brief get diagnostics Component names of Fields to output */
     amrex::Vector<std::string>& getDiagComps () { return m_diags.getComps(); };
     /** \brief get diagnostics Component names of Fields to output */

src/Hipace.cpp

@@ -45,6 +45,7 @@ amrex::Real Hipace::m_external_Ez_slope = 0.;
 amrex::Real Hipace::m_external_Ez_uniform = 0.;
 amrex::Real Hipace::m_MG_tolerance_rel = 1.e-4;
 amrex::Real Hipace::m_MG_tolerance_abs = 0.;
+int Hipace::m_nlev = 0;
 
 Hipace&
 Hipace::GetInstance ()
@@ -119,6 +120,7 @@ Hipace::Hipace () :
     pph.query("MG_tolerance_rel", m_MG_tolerance_rel);
     pph.query("MG_tolerance_abs", m_MG_tolerance_abs);
 
+    m_nlev = maxLevel()+1;
     if (maxLevel() > 0) {
         AMREX_ALWAYS_ASSERT(maxLevel() < 2);
         amrex::Array<amrex::Real, AMREX_SPACEDIM> loc_array;
@@ -246,7 +248,6 @@ void
 Hipace::MakeNewLevelFromScratch (
     int lev, amrex::Real /*time*/, const amrex::BoxArray& ba, const amrex::DistributionMapping&)
 {
-    AMREX_ALWAYS_ASSERT(lev == 0);
 
     // We are going to ignore the DistributionMapping argument and build our own.
     amrex::DistributionMapping dm;
@@ -360,7 +361,7 @@ Hipace::Evolve ()
     for (int step = m_numprocs_z - 1 - m_rank_z; step <= m_max_step; step += m_numprocs_z)
     {
 #ifdef HIPACE_USE_OPENPMD
-        m_openpmd_writer.InitDiagnostics(step, m_output_period, m_max_step);
+        m_openpmd_writer.InitDiagnostics(step, m_output_period, m_max_step, m_nlev);
 #endif
 
         if (m_verbose>=1) std::cout<<"Rank "<<rank<<" started  step "<<step<<" with dt = "<<m_dt<<'\n';
@@ -389,7 +390,13 @@ Hipace::Evolve ()
             if (it>0) m_multi_beam.PackLocalGhostParticles(it-1, m_box_sorters);
 
             const amrex::Box& bx = boxArray(lev)[it];
-            ResizeFDiagFAB(bx, lev);
+
+            // FIXME: dirty workaround to not touch the box in general
+            // but resize the boxes on all levels for IO.
+            for (int lev_loc = 0; lev_loc < m_nlev; ++lev_loc) {
+                const amrex::Box& bx_loc = boxArray(lev_loc)[it];
+                ResizeFDiagFAB(bx_loc, lev_loc);
+            }
 
             amrex::Vector<BeamBins> bins;
             bins = m_multi_beam.findParticlesInEachSlice(lev, it, bx, geom[lev], m_box_sorters);
@@ -524,7 +531,7 @@ Hipace::SolveOneSlice (int islice, int lev, const int ibox,
     // Push beam particles
     m_multi_beam.AdvanceBeamParticlesSlice(m_fields, geom[lev], lev, islice, bx, bins, m_box_sorters, ibox);
 
-    FillDiagnostics(lev, islice);
+    FillDiagnostics(islice);
 
     m_fields.ShiftSlices(lev);
 
@@ -1227,10 +1234,12 @@ Hipace::NotifyFinish (const int it, bool only_ghost)
 }
 
 void
-Hipace::FillDiagnostics (int lev, int i_slice)
+Hipace::FillDiagnostics (int i_slice)
 {
-    m_fields.Copy(lev, i_slice, FieldCopyType::StoF, 0, 0, Comps[WhichSlice::This]["N"],
-                  m_diags.getF(lev), m_diags.sliceDir(), Geom(lev));
+    for (int lev = 0; lev < m_nlev; ++lev) {
+        m_fields.Copy(lev, i_slice, FieldCopyType::StoF, 0, 0, Comps[WhichSlice::This]["N"],
+                      m_diags.getF(lev), m_diags.sliceDir(), Geom(lev));
+    }
 }
 
 void
@@ -1247,10 +1256,9 @@ Hipace::WriteDiagnostics (int output_step, const int it, const OpenPMDWriterCall
     const amrex::Vector< std::string > beamnames = getDiagBeamNames();
 
 #ifdef HIPACE_USE_OPENPMD
-    constexpr int lev = 0;
     m_openpmd_writer.WriteDiagnostics(getDiagF(), m_multi_beam, getDiagGeom(),
-                        m_physical_time, output_step, lev, getDiagSliceDir(), varnames, beamnames,
-                        it, m_box_sorters, geom[lev], call_type);
+                        m_physical_time, output_step, m_nlev, getDiagSliceDir(), varnames, beamnames,
+                        it, m_box_sorters, geom, call_type);
 #else
     amrex::Print()<<"WARNING: HiPACE++ compiled without openPMD support, the simulation has no I/O.\n";
 #endif

src/diagnostics/OpenPMDWriter.H

@@ -67,12 +67,13 @@ private:
      * \param[in] a_box_sorter_vec Vector (over species) of particles sorted by box
      * \param[in] geom Geometry of the simulation, to get the cell size etc.
      * \param[in] beamnames list of the names of the beam to be written to file
+     * \param[in] lev MR level
      */
     void WriteBeamParticleData (MultiBeam& beams, openPMD::Iteration iteration,
                                 const int output_step, const int it,
                                 const amrex::Vector<BoxSorter>& a_box_sorter_vec,
                                 const amrex::Geometry& geom,
-                                const amrex::Vector< std::string > beamnames);
+                                const amrex::Vector< std::string > beamnames, const int lev);
 
     /** \brief writing openPMD field data
      *
@@ -82,24 +83,25 @@ private:
      * \param[in] varnames list of variable names for the fields (ExmBy, EypBx, Ey, ...)
      * \param[in,out] iteration openPMD iteration to which the data is written
      * \param[in] output_step current time step to dump
+     * \param[in] lev MR level
      */
     void WriteFieldData (amrex::FArrayBox const& fab, amrex::Geometry const& geom,
                          const int slice_dir, const amrex::Vector< std::string > varnames,
-                         openPMD::Iteration iteration, const int output_step);
+                         openPMD::Iteration iteration, const int output_step, const int lev);
 
     /** Named Beam SoA attributes per particle as defined in BeamIdx
      */
     amrex::Vector<std::string> m_real_names {"weighting", "momentum_x", "momentum_y", "momentum_z"};
 
-    /** Parallel openPMD-api Series object for output */
-    std::unique_ptr< openPMD::Series > m_outputSeries;
+    /** vector over levels of openPMD-api Series object for output */
+    amrex::Vector<std::unique_ptr< openPMD::Series >> m_outputSeries;
 
     /** openPMD backend: h5, bp, or json. Default depends on what is available */
     std::string m_openpmd_backend = "default";
 
     /** Last iteration that was written to file.
      * This is stored to make sure we don't write the last iteration multiple times. */
-    int m_last_output_dumped = -1;
+    amrex::Vector<int> m_last_output_dumped;
 
     /** vector of length nbeams with the numbers of particles already written to file */
     amrex::Vector<uint64_t> m_offset;
@@ -114,8 +116,10 @@ public:
      * \param[in] output_step current iteration
      * \param[in] output_period output period
      * \param[in] max_step maximum time step of the simulation
+     * \param[in] nlev number of mesh refinement levels
      */
-    void InitDiagnostics (const int output_step, const int output_period, const int max_step);
+    void InitDiagnostics (const int output_step, const int output_period, const int max_step,
+                          const int nlev);
 
     /** \brief writing openPMD data
      *
@@ -124,7 +128,7 @@ public:
      * \param[in] geom Geometry of the simulation, to get the cell size etc.
      * \param[in] physical_time Physical time of the currenerationt it.
      * \param[in] output_step current iteration to be written to file
-     * \param[in] lev MR level
+     * \param[in] nlev number of MR levels
      * \param[in] slice_dir direction of slicing. 0=x, 1=y, -1=no slicing (3D array)
      * \param[in] varnames list of variable names for the fields (ExmBy, EypBx, Ey, ...)
      * \param[in] beamnames list of the names of the beam to be written to file
@@ -136,13 +140,14 @@ public:
     void WriteDiagnostics(
         amrex::Vector<amrex::FArrayBox> const& a_mf, MultiBeam& a_multi_beam,
         amrex::Vector<amrex::Geometry> const& geom,
-        const amrex::Real physical_time, const int output_step, const int lev,
+        const amrex::Real physical_time, const int output_step, const int nlev,
         const int slice_dir, const amrex::Vector< std::string > varnames,
         const amrex::Vector< std::string > beamnames, const int it,
-        const amrex::Vector<BoxSorter>& a_box_sorter_vec, const amrex::Geometry& geom3D,
+        const amrex::Vector<BoxSorter>& a_box_sorter_vec, amrex::Vector<amrex::Geometry> const& geom3D,
         const OpenPMDWriterCallType call_type);
 
-    void reset () {m_outputSeries.reset();};
+    /** \brief Resets the openPMD series of all levels */
+    void reset ();
 
     /** Prefix/path for the output files */
     std::string m_file_prefix = "diags/hdf5";

src/diagnostics/OpenPMDWriter.cpp

@@ -20,29 +20,44 @@ OpenPMDWriter::OpenPMDWriter ()
 }
 
 void
-OpenPMDWriter::InitDiagnostics (const int output_step, const int output_period, const int max_step)
+OpenPMDWriter::InitDiagnostics (const int output_step, const int output_period, const int max_step,
+                                const int nlev)
 {
     HIPACE_PROFILE("OpenPMDWriter::InitDiagnostics()");
 
     // Dump every m_output_period steps and after last step
     if (output_period < 0 ||
        (!(output_step == max_step) && output_step % output_period != 0)) return;
 
-     // pick first available backend if default is chosen
-     if( m_openpmd_backend == "default" ) {
+    // pick first available backend if default is chosen
+    if( m_openpmd_backend == "default" ) {
 #if openPMD_HAVE_HDF5==1
         m_openpmd_backend = "h5";
 #elif openPMD_HAVE_ADIOS2==1
         m_openpmd_backend = "bp";
 #else
         m_openpmd_backend = "json";
 #endif
-     }
+    }
+
+    if (nlev > 1) {
+        for (int lev=0; lev<nlev; ++lev) {
+            std::string filename = m_file_prefix + "/lev_" + std::to_string(lev) +  "/openpmd_%06T."
+                                   + m_openpmd_backend;
+
+            m_outputSeries.push_back(std::make_unique< openPMD::Series >(
+                filename, openPMD::Access::CREATE) );
+            m_last_output_dumped.push_back(-1);
+        }
+    } else {
+        std::string filename = m_file_prefix + "/openpmd_%06T." + m_openpmd_backend;
+
+        m_outputSeries.push_back(std::make_unique< openPMD::Series >(
+            filename, openPMD::Access::CREATE) );
+        m_last_output_dumped.push_back(-1);
+    }
 
-    std::string filename = m_file_prefix + "/openpmd_%06T." + m_openpmd_backend;
 
-    m_outputSeries = std::make_unique< openPMD::Series >(
-        filename, openPMD::Access::CREATE);
 
     // TODO: meta-data: author, mesh path, extensions, software
 }
@@ -51,32 +66,35 @@ void
 OpenPMDWriter::WriteDiagnostics (
     amrex::Vector<amrex::FArrayBox> const& a_mf, MultiBeam& a_multi_beam,
     amrex::Vector<amrex::Geometry> const& geom,
-    const amrex::Real physical_time, const int output_step, const int lev,
+    const amrex::Real physical_time, const int output_step, const int nlev,
     const int slice_dir, const amrex::Vector< std::string > varnames,
     const amrex::Vector< std::string > beamnames, const int it,
-    const amrex::Vector<BoxSorter>& a_box_sorter_vec, const amrex::Geometry& geom3D,
+    const amrex::Vector<BoxSorter>& a_box_sorter_vec, amrex::Vector<amrex::Geometry> const& geom3D,
     const OpenPMDWriterCallType call_type)
 {
-    openPMD::Iteration iteration = m_outputSeries->iterations[output_step];
+    for (int lev=0; lev<nlev; ++lev) {
+        openPMD::Iteration iteration = m_outputSeries[lev]->iterations[output_step];
 
-    if (call_type == OpenPMDWriterCallType::beams ) {
-        iteration.setTime(physical_time);
-        WriteBeamParticleData(a_multi_beam, iteration, output_step, it, a_box_sorter_vec, geom3D, beamnames);
-        m_outputSeries->flush();
+        if (call_type == OpenPMDWriterCallType::beams ) {
+            iteration.setTime(physical_time);
+            if (lev == 0) {
+                WriteBeamParticleData(a_multi_beam, iteration, output_step, it, a_box_sorter_vec, geom3D[lev], beamnames, lev);
+            }
+            m_outputSeries[lev]->flush();
 
-    } else if (call_type == OpenPMDWriterCallType::fields ) {
-        WriteFieldData(a_mf[lev], geom[lev], slice_dir, varnames, iteration, output_step);
-        m_outputSeries->flush();
-        m_last_output_dumped = output_step;
+        } else if (call_type == OpenPMDWriterCallType::fields ) {
+            WriteFieldData(a_mf[lev], geom[lev], slice_dir, varnames, iteration, output_step, lev);
+            m_outputSeries[lev]->flush();
+            m_last_output_dumped[lev] = output_step;
+        }
     }
-
 }
 
 void
 OpenPMDWriter::WriteFieldData (
     amrex::FArrayBox const& fab, amrex::Geometry const& geom,
     const int slice_dir, const amrex::Vector< std::string > varnames,
-    openPMD::Iteration iteration, const int output_step)
+    openPMD::Iteration iteration, const int output_step, const int lev)
 {
     // todo: periodicity/boundary, field solver, particle pusher, etc.
     auto meshes = iteration.meshes;
@@ -118,7 +136,7 @@ OpenPMDWriter::WriteFieldData (
         // If slicing requested, remove number of points for the slicing direction
         if (slice_dir >= 0) global_size.erase(global_size.begin() + 2-slice_dir);
 
-        if (m_last_output_dumped != output_step) {
+        if (m_last_output_dumped[lev] != output_step) {
             openPMD::Dataset dataset(datatype, global_size);
             field_comp.resetDataset(dataset);
         }
@@ -129,7 +147,6 @@ OpenPMDWriter::WriteFieldData (
 
         data = openPMD::shareRaw( fab.dataPtr( icomp ) ); // non-owning view until flush()
 
-
         // Determine the offset and size of this data chunk in the global output
         amrex::IntVect const box_offset = data_box.smallEnd();
         openPMD::Offset chunk_offset = utils::getReversedVec(box_offset);
@@ -148,7 +165,7 @@ OpenPMDWriter::WriteBeamParticleData (MultiBeam& beams, openPMD::Iteration itera
                                       const int output_step, const int it,
                                       const amrex::Vector<BoxSorter>& a_box_sorter_vec,
                                       const amrex::Geometry& geom,
-                                      const amrex::Vector< std::string > beamnames)
+                                      const amrex::Vector< std::string > beamnames, const int lev)
 {
     HIPACE_PROFILE("WriteBeamParticleData()");
 
@@ -165,7 +182,7 @@ OpenPMDWriter::WriteBeamParticleData (MultiBeam& beams, openPMD::Iteration itera
         auto& beam = beams.getBeam(ibeam);
 
         const unsigned long long np = beams.get_total_num_particles(ibeam);
-        if (m_last_output_dumped != output_step) {
+        if (m_last_output_dumped[lev] != output_step) {
             SetupPos(beam_species, beam, np, geom);
             SetupRealProperties(beam_species, m_real_names, np);
         }
@@ -365,4 +382,11 @@ OpenPMDWriter::SaveRealProperty (BeamParticleContainer& pc,
     }
 }
 
+void OpenPMDWriter::reset ()
+{
+    for (int lev = 0; lev<m_outputSeries.size(); ++lev) {
+        m_outputSeries[lev].reset();
+    }
+}
+
 #endif // HIPACE_USE_OPENPMD

src/fields/Fields.cpp

@@ -32,6 +32,8 @@ Fields::AllocData (
         m_slices[lev][islice].setVal(0.0);
     }
 
+    // FIXME the Poisson solver must be constructed per level, here it is only constructed for lev 0
+    if (lev>0) return;
     // The Poisson solver operates on transverse slices only.
     // The constructor takes the BoxArray and the DistributionMap of a slice,
     // so the FFTPlans are built on a slice.