From 7c51dadaee3dbe8dbcc5921194750a46e5c69bd3 Mon Sep 17 00:00:00 2001 From: Severin Diederichs Date: Thu, 10 Jun 2021 15:53:39 +0200 Subject: [PATCH 01/10] debug printouts --- src/Hipace.cpp | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/Hipace.cpp b/src/Hipace.cpp index da4bc1f314..e1ef61e8f8 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -246,7 +246,7 @@ void Hipace::MakeNewLevelFromScratch ( int lev, amrex::Real /*time*/, const amrex::BoxArray& ba, const amrex::DistributionMapping&) { - AMREX_ALWAYS_ASSERT(lev == 0); + // AMREX_ALWAYS_ASSERT(lev == 0); // We are going to ignore the DistributionMapping argument and build our own. amrex::DistributionMapping dm; @@ -256,6 +256,16 @@ Hipace::MakeNewLevelFromScratch ( const int nboxes_x = m_numprocs_x; const int nboxes_y = m_numprocs_y; const int nboxes_z = (m_boxes_in_z == 1) ? ncells_global[2] / box_size[2] : m_boxes_in_z; + + amrex::Print() << " level " << lev << "\n"; + amrex::Print() << " box array " << ba << "\n"; + amrex::Print() << " dm " << dm << "\n"; + amrex::Print() << "ncells_global[0] " << ncells_global[0] << " ncells_global[1] " << ncells_global[1] << " ncells_global[2] " << ncells_global[2] << "\n"; + amrex::Print() << " box_size[0] " << box_size[0] << " box_size[1] " << box_size[1] << " box_size[2] " << box_size[2] << "\n"; + amrex::Print() << " nboxes_x " << nboxes_x << " nboxes_y " << nboxes_y << " nboxes_z " << nboxes_z << " ba.size() " << ba.size() << "\n"; + + + AMREX_ALWAYS_ASSERT(static_cast(nboxes_x) * static_cast(nboxes_y) * static_cast(nboxes_z) == ba.size()); From 3f003dbf45cf0c41ec24e410cc8349528b904013 Mon Sep 17 00:00:00 2001 From: Severin Diederichs Date: Fri, 11 Jun 2021 12:39:12 +0200 Subject: [PATCH 02/10] add IO for refined level --- src/Hipace.H | 3 +- src/Hipace.cpp | 30 ++++++++---- src/diagnostics/OpenPMDWriter.H | 32 ++++++++----- src/diagnostics/OpenPMDWriter.cpp | 77 ++++++++++++++++++++----------- src/fields/Fields.cpp | 2 + 5 files changed, 95 insertions(+), 49 deletions(-) diff --git a/src/Hipace.H b/src/Hipace.H index 7b8a73db7a..6207adc0ab 100644 --- a/src/Hipace.H +++ b/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& getDiagComps () { return m_diags.getComps(); }; /** \brief get diagnostics Component names of Fields to output */ diff --git a/src/Hipace.cpp b/src/Hipace.cpp index e1ef61e8f8..3bca77ae2f 100644 --- a/src/Hipace.cpp +++ b/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 loc_array; @@ -370,7 +372,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, maxLevel()+1); #endif if (m_verbose>=1) std::cout<<"Rank "< beamnames = getDiagBeamNames(); #ifdef HIPACE_USE_OPENPMD - constexpr int lev = 0; + const int nlev = maxLevel()+1; 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, 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 diff --git a/src/diagnostics/OpenPMDWriter.H b/src/diagnostics/OpenPMDWriter.H index f033240031..18a340369f 100644 --- a/src/diagnostics/OpenPMDWriter.H +++ b/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& 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, + void WriteFieldData (amrex::FArrayBox & 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 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> 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 m_last_output_dumped; /** vector of length nbeams with the numbers of particles already written to file */ amrex::Vector 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 @@ -134,15 +138,19 @@ public: * \param[in] call_type whether the beams or the fields should be written to file */ void WriteDiagnostics( - amrex::Vector const& a_mf, MultiBeam& a_multi_beam, + amrex::Vector & a_mf, MultiBeam& a_multi_beam, amrex::Vector 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& a_box_sorter_vec, const amrex::Geometry& geom3D, + const amrex::Vector& a_box_sorter_vec, amrex::Vector const& geom3D, const OpenPMDWriterCallType call_type); - void reset () {m_outputSeries.reset();}; + /** \brief Resets the openPMD series of all levels + * + * \param[in] nlev number of mesh refinement levels + */ + void reset (const int nlev); /** Prefix/path for the output files */ std::string m_file_prefix = "diags/hdf5"; diff --git a/src/diagnostics/OpenPMDWriter.cpp b/src/diagnostics/OpenPMDWriter.cpp index 2697894f6c..17c74f091a 100644 --- a/src/diagnostics/OpenPMDWriter.cpp +++ b/src/diagnostics/OpenPMDWriter.cpp @@ -20,7 +20,8 @@ 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()"); @@ -28,8 +29,8 @@ OpenPMDWriter::InitDiagnostics (const int output_step, const int output_period, 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 @@ -37,46 +38,65 @@ OpenPMDWriter::InitDiagnostics (const int output_step, const int output_period, #else m_openpmd_backend = "json"; #endif - } + } + + if (nlev > 0) { + for (int lev=0; lev( + 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 } void OpenPMDWriter::WriteDiagnostics ( - amrex::Vector const& a_mf, MultiBeam& a_multi_beam, + amrex::Vector & a_mf, MultiBeam& a_multi_beam, amrex::Vector 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& a_box_sorter_vec, const amrex::Geometry& geom3D, + const amrex::Vector& a_box_sorter_vec, amrex::Vector const& geom3D, const OpenPMDWriterCallType call_type) { - openPMD::Iteration iteration = m_outputSeries->iterations[output_step]; + for (int lev=0; leviterations[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, + amrex::FArrayBox & 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 +138,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 +149,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 +167,7 @@ OpenPMDWriter::WriteBeamParticleData (MultiBeam& beams, openPMD::Iteration itera const int output_step, const int it, const amrex::Vector& 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 +184,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 +384,10 @@ OpenPMDWriter::SaveRealProperty (BeamParticleContainer& pc, } } +void OpenPMDWriter::reset (const int nlev) +{ + for (int lev = 0; lev0) 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. From d9d718979b18e3dbe6e98464cc355cb12ad8f647 Mon Sep 17 00:00:00 2001 From: Severin Diederichs Date: Fri, 11 Jun 2021 12:44:36 +0200 Subject: [PATCH 03/10] cleaning --- src/Hipace.cpp | 21 ++++----------------- 1 file changed, 4 insertions(+), 17 deletions(-) diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 3bca77ae2f..07dfd36661 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -248,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; @@ -258,16 +257,6 @@ Hipace::MakeNewLevelFromScratch ( const int nboxes_x = m_numprocs_x; const int nboxes_y = m_numprocs_y; const int nboxes_z = (m_boxes_in_z == 1) ? ncells_global[2] / box_size[2] : m_boxes_in_z; - - amrex::Print() << " level " << lev << "\n"; - amrex::Print() << " box array " << ba << "\n"; - amrex::Print() << " dm " << dm << "\n"; - amrex::Print() << "ncells_global[0] " << ncells_global[0] << " ncells_global[1] " << ncells_global[1] << " ncells_global[2] " << ncells_global[2] << "\n"; - amrex::Print() << " box_size[0] " << box_size[0] << " box_size[1] " << box_size[1] << " box_size[2] " << box_size[2] << "\n"; - amrex::Print() << " nboxes_x " << nboxes_x << " nboxes_y " << nboxes_y << " nboxes_z " << nboxes_z << " ba.size() " << ba.size() << "\n"; - - - AMREX_ALWAYS_ASSERT(static_cast(nboxes_x) * static_cast(nboxes_y) * static_cast(nboxes_z) == ba.size()); @@ -372,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, maxLevel()+1); + m_openpmd_writer.InitDiagnostics(step, m_output_period, m_max_step, m_nlev); #endif if (m_verbose>=1) std::cout<<"Rank "< Date: Fri, 11 Jun 2021 12:47:19 +0200 Subject: [PATCH 04/10] clearer explanation for workaround --- src/Hipace.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 07dfd36661..8f3defa534 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -391,7 +391,8 @@ Hipace::Evolve () const amrex::Box& bx = boxArray(lev)[it]; - // FIXME: dirty workaround to not touch the box in general but only for IO. + // 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); From 6b7619a31c8d944badbc7440affaed64a912c1eb Mon Sep 17 00:00:00 2001 From: Severin Diederichs Date: Fri, 11 Jun 2021 12:49:02 +0200 Subject: [PATCH 05/10] add const for fabs in IO, which was removed for testing purposes --- src/diagnostics/OpenPMDWriter.H | 4 ++-- src/diagnostics/OpenPMDWriter.cpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/diagnostics/OpenPMDWriter.H b/src/diagnostics/OpenPMDWriter.H index 18a340369f..d3fbf427ed 100644 --- a/src/diagnostics/OpenPMDWriter.H +++ b/src/diagnostics/OpenPMDWriter.H @@ -85,7 +85,7 @@ private: * \param[in] output_step current time step to dump * \param[in] lev MR level */ - void WriteFieldData (amrex::FArrayBox & fab, amrex::Geometry const& geom, + 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, const int lev); @@ -138,7 +138,7 @@ public: * \param[in] call_type whether the beams or the fields should be written to file */ void WriteDiagnostics( - amrex::Vector & a_mf, MultiBeam& a_multi_beam, + amrex::Vector const& a_mf, MultiBeam& a_multi_beam, amrex::Vector const& geom, const amrex::Real physical_time, const int output_step, const int nlev, const int slice_dir, const amrex::Vector< std::string > varnames, diff --git a/src/diagnostics/OpenPMDWriter.cpp b/src/diagnostics/OpenPMDWriter.cpp index 17c74f091a..8849fa378e 100644 --- a/src/diagnostics/OpenPMDWriter.cpp +++ b/src/diagnostics/OpenPMDWriter.cpp @@ -66,7 +66,7 @@ OpenPMDWriter::InitDiagnostics (const int output_step, const int output_period, void OpenPMDWriter::WriteDiagnostics ( - amrex::Vector & a_mf, MultiBeam& a_multi_beam, + amrex::Vector const& a_mf, MultiBeam& a_multi_beam, amrex::Vector const& geom, const amrex::Real physical_time, const int output_step, const int nlev, const int slice_dir, const amrex::Vector< std::string > varnames, @@ -94,7 +94,7 @@ OpenPMDWriter::WriteDiagnostics ( void OpenPMDWriter::WriteFieldData ( - amrex::FArrayBox & fab, amrex::Geometry const& geom, + 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, const int lev) { From 13f2edb966a3fa8f06cda4412532b95243b8a2bb Mon Sep 17 00:00:00 2001 From: Severin Diederichs Date: Fri, 11 Jun 2021 12:55:31 +0200 Subject: [PATCH 06/10] more cleaning --- src/diagnostics/OpenPMDWriter.cpp | 5 +---- src/fields/Fields.cpp | 2 +- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/src/diagnostics/OpenPMDWriter.cpp b/src/diagnostics/OpenPMDWriter.cpp index 8849fa378e..6106711308 100644 --- a/src/diagnostics/OpenPMDWriter.cpp +++ b/src/diagnostics/OpenPMDWriter.cpp @@ -43,11 +43,10 @@ OpenPMDWriter::InitDiagnostics (const int output_step, const int output_period, if (nlev > 0) { for (int lev=0; lev( filename, openPMD::Access::CREATE) ); - m_last_output_dumped.push_back(-1); } } else { @@ -55,7 +54,6 @@ OpenPMDWriter::InitDiagnostics (const int output_step, const int output_period, m_outputSeries.push_back(std::make_unique< openPMD::Series >( filename, openPMD::Access::CREATE) ); - m_last_output_dumped.push_back(-1); } @@ -389,5 +387,4 @@ void OpenPMDWriter::reset (const int nlev) for (int lev = 0; lev0) return; // The Poisson solver operates on transverse slices only. // The constructor takes the BoxArray and the DistributionMap of a slice, From c2dcab6219a473b9159e7fd4afaed1c421f4b8c2 Mon Sep 17 00:00:00 2001 From: Severin Diederichs Date: Fri, 11 Jun 2021 13:27:47 +0200 Subject: [PATCH 07/10] fix typo to enable untouched IO without MR --- src/diagnostics/OpenPMDWriter.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/diagnostics/OpenPMDWriter.cpp b/src/diagnostics/OpenPMDWriter.cpp index 6106711308..3c5371d849 100644 --- a/src/diagnostics/OpenPMDWriter.cpp +++ b/src/diagnostics/OpenPMDWriter.cpp @@ -40,7 +40,7 @@ OpenPMDWriter::InitDiagnostics (const int output_step, const int output_period, #endif } - if (nlev > 0) { + if (nlev > 1) { for (int lev=0; lev Date: Fri, 11 Jun 2021 21:34:14 +0200 Subject: [PATCH 08/10] call reset based on length of m_openPMD_series array --- src/Hipace.cpp | 2 +- src/diagnostics/OpenPMDWriter.H | 7 ++----- src/diagnostics/OpenPMDWriter.cpp | 6 ++++-- 3 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 8f3defa534..ad47d53e62 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -440,7 +440,7 @@ Hipace::Evolve () } #ifdef HIPACE_USE_OPENPMD - if (m_output_period > 0) m_openpmd_writer.reset(m_nlev); + if (m_output_period > 0) m_openpmd_writer.reset(); #endif } diff --git a/src/diagnostics/OpenPMDWriter.H b/src/diagnostics/OpenPMDWriter.H index d3fbf427ed..eb371ed00d 100644 --- a/src/diagnostics/OpenPMDWriter.H +++ b/src/diagnostics/OpenPMDWriter.H @@ -146,11 +146,8 @@ public: const amrex::Vector& a_box_sorter_vec, amrex::Vector const& geom3D, const OpenPMDWriterCallType call_type); - /** \brief Resets the openPMD series of all levels - * - * \param[in] nlev number of mesh refinement levels - */ - void reset (const int nlev); + /** \brief Resets the openPMD series of all levels */ + void reset (); /** Prefix/path for the output files */ std::string m_file_prefix = "diags/hdf5"; diff --git a/src/diagnostics/OpenPMDWriter.cpp b/src/diagnostics/OpenPMDWriter.cpp index 3c5371d849..7c5c1063df 100644 --- a/src/diagnostics/OpenPMDWriter.cpp +++ b/src/diagnostics/OpenPMDWriter.cpp @@ -382,9 +382,11 @@ OpenPMDWriter::SaveRealProperty (BeamParticleContainer& pc, } } -void OpenPMDWriter::reset (const int nlev) +void OpenPMDWriter::reset () { - for (int lev = 0; lev Date: Sat, 12 Jun 2021 21:58:39 +0200 Subject: [PATCH 09/10] using multiple Poisson solvers --- src/Hipace.cpp | 2 +- src/fields/Fields.H | 4 ++-- src/fields/Fields.cpp | 47 +++++++++++++++++++++---------------------- 3 files changed, 26 insertions(+), 27 deletions(-) diff --git a/src/Hipace.cpp b/src/Hipace.cpp index ad47d53e62..4ef2390ca1 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -281,7 +281,7 @@ Hipace::MakeNewLevelFromScratch ( DefineSliceGDB(ba, dm); // Note: we pass ba[0] as a dummy box, it will be resized properly in the loop over boxes in Evolve m_diags.AllocData(lev, ba[0], Comps[WhichSlice::This]["N"], Geom(lev)); - m_fields.AllocData(lev, ba, dm, Geom(lev), m_slice_ba, m_slice_dm); + m_fields.AllocData(lev, ba, dm, Geom(), m_slice_ba, m_slice_dm); } void diff --git a/src/fields/Fields.H b/src/fields/Fields.H index 0c5b028894..9b081eb34f 100644 --- a/src/fields/Fields.H +++ b/src/fields/Fields.H @@ -82,11 +82,11 @@ public: */ void AllocData ( int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm, - amrex::Geometry const& geom, const amrex::BoxArray& slice_ba, + amrex::Vector const& geom, const amrex::BoxArray& slice_ba, const amrex::DistributionMapping& slice_dm); /** Class to handle transverse FFT Poisson solver on 1 slice */ - std::unique_ptr m_poisson_solver; + amrex::Vector> m_poisson_solver; /** get function for the 2D slices */ amrex::Vector>& getSlices () {return m_slices; } /** get function for the 2D slices diff --git a/src/fields/Fields.cpp b/src/fields/Fields.cpp index 708df89911..81f88be653 100644 --- a/src/fields/Fields.cpp +++ b/src/fields/Fields.cpp @@ -15,7 +15,8 @@ Fields::Fields (Hipace const* a_hipace) void Fields::AllocData ( int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm, - amrex::Geometry const& geom, const amrex::BoxArray& slice_ba, const amrex::DistributionMapping& slice_dm) + amrex::Vector const& geom, const amrex::BoxArray& slice_ba, + const amrex::DistributionMapping& slice_dm) { HIPACE_PROFILE("Fields::AllocData()"); // Need at least 1 guard cell transversally for transverse derivative @@ -32,21 +33,19 @@ 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. if (m_do_dirichlet_poisson){ - m_poisson_solver = std::unique_ptr( + m_poisson_solver.push_back(std::unique_ptr( new FFTPoissonSolverDirichlet(getSlices(lev, WhichSlice::This).boxArray(), getSlices(lev, WhichSlice::This).DistributionMap(), - geom)); + geom[lev])) ); } else { - m_poisson_solver = std::unique_ptr( + m_poisson_solver.push_back(std::unique_ptr( new FFTPoissonSolverPeriodic(getSlices(lev, WhichSlice::This).boxArray(), getSlices(lev, WhichSlice::This).DistributionMap(), - geom)); + geom[lev])) ); } } @@ -255,14 +254,14 @@ Fields::SolvePoissonExmByAndEypBx (amrex::Geometry const& geom, const MPI_Comm& Comps[WhichSlice::This]["Psi"], 1); // calculating the right-hand side 1/episilon0 * -(rho-Jz/c) - amrex::MultiFab::Copy(m_poisson_solver->StagingArea(), getSlices(lev, WhichSlice::This), + amrex::MultiFab::Copy(m_poisson_solver[lev]->StagingArea(), getSlices(lev, WhichSlice::This), Comps[WhichSlice::This]["jz"], 0, 1, 0); - m_poisson_solver->StagingArea().mult(-1./phys_const.c); - amrex::MultiFab::Add(m_poisson_solver->StagingArea(), getSlices(lev, WhichSlice::This), + m_poisson_solver[lev]->StagingArea().mult(-1./phys_const.c); + amrex::MultiFab::Add(m_poisson_solver[lev]->StagingArea(), getSlices(lev, WhichSlice::This), Comps[WhichSlice::This]["rho"], 0, 1, 0); - m_poisson_solver->StagingArea().mult(-1./phys_const.ep0); + m_poisson_solver[lev]->StagingArea().mult(-1./phys_const.ep0); - m_poisson_solver->SolvePoissonEquation(lhs); + m_poisson_solver[lev]->SolvePoissonEquation(lhs); /* ---------- Transverse FillBoundary Psi ---------- */ amrex::ParallelContext::push(m_comm_xy); @@ -306,7 +305,7 @@ Fields::SolvePoissonEz (amrex::Geometry const& geom, const int lev) // from the slice MF, and store in the staging area of poisson_solver TransverseDerivative( getSlices(lev, WhichSlice::This), - m_poisson_solver->StagingArea(), + m_poisson_solver[lev]->StagingArea(), Direction::x, geom.CellSize(Direction::x), 1./(phys_const.ep0*phys_const.c), @@ -315,7 +314,7 @@ Fields::SolvePoissonEz (amrex::Geometry const& geom, const int lev) TransverseDerivative( getSlices(lev, WhichSlice::This), - m_poisson_solver->StagingArea(), + m_poisson_solver[lev]->StagingArea(), Direction::y, geom.CellSize(Direction::y), 1./(phys_const.ep0*phys_const.c), @@ -324,7 +323,7 @@ Fields::SolvePoissonEz (amrex::Geometry const& geom, const int lev) // Solve Poisson equation. // The RHS is in the staging area of poisson_solver. // The LHS will be returned as lhs. - m_poisson_solver->SolvePoissonEquation(lhs); + m_poisson_solver[lev]->SolvePoissonEquation(lhs); } void @@ -338,7 +337,7 @@ Fields::SolvePoissonBx (amrex::MultiFab& Bx_iter, amrex::Geometry const& geom, c // and store in the staging area of poisson_solver TransverseDerivative( getSlices(lev, WhichSlice::This), - m_poisson_solver->StagingArea(), + m_poisson_solver[lev]->StagingArea(), Direction::y, geom.CellSize(Direction::y), -phys_const.mu0, @@ -348,7 +347,7 @@ Fields::SolvePoissonBx (amrex::MultiFab& Bx_iter, amrex::Geometry const& geom, c LongitudinalDerivative( getSlices(lev, WhichSlice::Previous1), getSlices(lev, WhichSlice::Next), - m_poisson_solver->StagingArea(), + m_poisson_solver[lev]->StagingArea(), geom.CellSize(Direction::z), phys_const.mu0, SliceOperatorType::Add, @@ -357,7 +356,7 @@ Fields::SolvePoissonBx (amrex::MultiFab& Bx_iter, amrex::Geometry const& geom, c // Solve Poisson equation. // The RHS is in the staging area of poisson_solver. // The LHS will be returned as lhs. - m_poisson_solver->SolvePoissonEquation(Bx_iter); + m_poisson_solver[lev]->SolvePoissonEquation(Bx_iter); } void @@ -371,7 +370,7 @@ Fields::SolvePoissonBy (amrex::MultiFab& By_iter, amrex::Geometry const& geom, c // and store in the staging area of poisson_solver TransverseDerivative( getSlices(lev, WhichSlice::This), - m_poisson_solver->StagingArea(), + m_poisson_solver[lev]->StagingArea(), Direction::x, geom.CellSize(Direction::x), phys_const.mu0, @@ -381,7 +380,7 @@ Fields::SolvePoissonBy (amrex::MultiFab& By_iter, amrex::Geometry const& geom, c LongitudinalDerivative( getSlices(lev, WhichSlice::Previous1), getSlices(lev, WhichSlice::Next), - m_poisson_solver->StagingArea(), + m_poisson_solver[lev]->StagingArea(), geom.CellSize(Direction::z), -phys_const.mu0, SliceOperatorType::Add, @@ -390,7 +389,7 @@ Fields::SolvePoissonBy (amrex::MultiFab& By_iter, amrex::Geometry const& geom, c // Solve Poisson equation. // The RHS is in the staging area of poisson_solver. // The LHS will be returned as lhs. - m_poisson_solver->SolvePoissonEquation(By_iter); + m_poisson_solver[lev]->SolvePoissonEquation(By_iter); } void @@ -407,7 +406,7 @@ Fields::SolvePoissonBz (amrex::Geometry const& geom, const int lev) // from the slice MF, and store in the staging area of m_poisson_solver TransverseDerivative( getSlices(lev, WhichSlice::This), - m_poisson_solver->StagingArea(), + m_poisson_solver[lev]->StagingArea(), Direction::y, geom.CellSize(Direction::y), phys_const.mu0, @@ -416,7 +415,7 @@ Fields::SolvePoissonBz (amrex::Geometry const& geom, const int lev) TransverseDerivative( getSlices(lev, WhichSlice::This), - m_poisson_solver->StagingArea(), + m_poisson_solver[lev]->StagingArea(), Direction::x, geom.CellSize(Direction::x), -phys_const.mu0, @@ -425,7 +424,7 @@ Fields::SolvePoissonBz (amrex::Geometry const& geom, const int lev) // Solve Poisson equation. // The RHS is in the staging area of m_poisson_solver. // The LHS will be returned as lhs. - m_poisson_solver->SolvePoissonEquation(lhs); + m_poisson_solver[lev]->SolvePoissonEquation(lhs); } void From 822c06cf054e738e8d0a66e7827942c79239ea0d Mon Sep 17 00:00:00 2001 From: Severin Diederichs Date: Mon, 14 Jun 2021 13:22:37 +0200 Subject: [PATCH 10/10] forgot to remove unused parameter in merge conflict --- src/fields/Fields.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/fields/Fields.cpp b/src/fields/Fields.cpp index 43575ad7af..b29126a88e 100644 --- a/src/fields/Fields.cpp +++ b/src/fields/Fields.cpp @@ -22,9 +22,6 @@ Fields::AllocData ( int nguards_xy = std::max(1, Hipace::m_depos_order_xy); m_slices_nguards = {nguards_xy, nguards_xy, 0}; - // Only xy slices need guard cells, there is no deposition to/gather from the output array F. - amrex::IntVect nguards_F = amrex::IntVect(0,0,0); - for (int islice=0; islice