diff options
Diffstat (limited to 'Source')
38 files changed, 1562 insertions, 212 deletions
diff --git a/Source/Diagnostics/BTD_Plotfile_Header_Impl.H b/Source/Diagnostics/BTD_Plotfile_Header_Impl.H index 37540a5f1..a7fc9fd1e 100644 --- a/Source/Diagnostics/BTD_Plotfile_Header_Impl.H +++ b/Source/Diagnostics/BTD_Plotfile_Header_Impl.H @@ -1,3 +1,9 @@ +/* Copyright 2021 Revathi Jambunathan + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ #ifndef BTD_PLOTFILE_HEADER_IMPL_H #define BTD_PLOTFILE_HEADER_IMPL_H @@ -229,4 +235,113 @@ class BTDMultiFabHeaderImpl amrex::Vector<amrex::Real> src); }; +/** + * \brief Class to read, modify, and write species header when + * back-transformed diag format is selected as plotfile. + * This class enables multiple BTD particle buffers to be interweaved and stitched into + * a single plotfile with a single Header. + */ +class BTDSpeciesHeaderImpl +{ +public: + /** Constructor. + * \param[in] string containing path of Headerfile + */ + BTDSpeciesHeaderImpl (std::string const& Headerfile_path, std::string const& species_name); + ~BTDSpeciesHeaderImpl () = default; + /** Reads the Header file for BTD species*/ + void ReadHeader (); + /** Writes the meta-data of species Header file, with path, m_Header_path*/ + void WriteHeader (); + /** Set data Index of the data-file, DATAXXXXX, that the particles belong to*/ + void set_DataIndex (const int lev, const int box_id, const int data_index); + /** Add new_particles to existing to obtain the total number of particles of the species. + \param[in] new_particles, total particles in the new buffer + */ + void AddTotalParticles (const int new_particles) { m_total_particles += new_particles;} + /** Increment number of boxes in a box array by 1, with every flush. */ + void IncrementParticleBoxArraySize () { m_particleBoxArray_size[m_finestLevel]++; } + /** Append particle info of the newly added box, namely, + m_which_data, m_particles_per_box, m_offset_per_box. + */ + void AppendParticleInfoForNewBox (const int data_index, const int particles_per_box, const int offset); + /** Vector of data indices of all particle data for each level*/ + amrex::Vector< amrex::Vector<int> > m_which_data; + /** Vector of particles per box for each level*/ + amrex::Vector< amrex::Vector<int> > m_particles_per_box; + /** Vector of particle offset per box of each particle data file for each level*/ + amrex::Vector< amrex::Vector<int> > m_offset_per_box; + /** Vector of box array size per for each level*/ + amrex::Vector<int> m_particleBoxArray_size; + /** Total number of particles of species, m_species_name, in the particle output*/ + int m_total_particles; + /** Id of the next particle*/ + int m_nextid; + /** finest level of the particle output*/ + int m_finestLevel; + private: + /** Path of the headerfile */ + std::string m_Header_path; + /** Species name of the particles flushed out */ + std::string m_species_name; + /** String containing file version of the species output */ + std::string m_file_version; + /** Number of dimensions stored in the plotfile, should be same as AMREX_SPACEDIM */ + int m_spacedim; + /** Number of real attributes*/ + int m_num_output_real; + /** Number of integer attributes*/ + int m_num_output_int; + /** The real attribute names*/ + amrex::Vector<std::string> m_real_comp_names; + /** The real integer names*/ + amrex::Vector<std::string> m_int_comp_names; + /** if the output is a checkpoint */ + bool m_is_checkpoint; +}; + +/** + * \brief Class to read, modify, and write particle header file, Particle_H, when + * back-transformed diag format is selected as plotfile. + * This class enables multiple particle buffers to be interweaved and stitched into + * a single plotfile with a single Particle_H file. + */ +class BTDParticleDataHeaderImpl +{ +public: + /** Constructor. + * \param[in] string containing path of Headerfile + */ + BTDParticleDataHeaderImpl (std::string const& Headerfile_path); + /** Destructor */ + ~BTDParticleDataHeaderImpl () = default; + /** Reads the particle header file at m_Header_path and stores its data*/ + void ReadHeader (); + /** Writes the meta-data of particle box array in header file, with path, m_Header_path*/ + void WriteHeader (); + /** Returns the size of the box array, m_ba_size */ + int ba_size () {return m_ba_size; } + /** Increases Box array size, m_ba_size, by add_size + * \param[in] add_size + */ + void IncreaseBoxArraySize ( const int add_size) { m_ba_size += add_size;} + /** Returns box corresponding to the ith box in the BoxArray, m_ba. + * \param[in] int ibox, index of the box in the BoxArray. + */ + amrex::Box ba_box (int ibox) {return m_ba[ibox]; } + /** Resize boxArray, m_ba, to size, m_ba_size. */ + void ResizeBoxArray () { m_ba.resize(m_ba_size); } + /** Set Box indices of the ith-box in Box Array, m_ba, to the new Box, ba_box. + * \param[in] int ibox, index of the ith box in BoxArray, m_ba. + * \param[in] amrex::Box box dimensions corresponding to the ith Fab. + */ + void SetBox (int ibox, amrex::Box ba_box) { m_ba.set(ibox, ba_box); } + /** Size of BoxArray, m_ba*/ + int m_ba_size; + /** BoxArray for particle output*/ + amrex::BoxArray m_ba; + /** string containing path of the particle output of species, species_name.*/ + std::string m_Header_path; +}; + #endif diff --git a/Source/Diagnostics/BTD_Plotfile_Header_Impl.cpp b/Source/Diagnostics/BTD_Plotfile_Header_Impl.cpp index 6b02760c4..bf9f8e877 100644 --- a/Source/Diagnostics/BTD_Plotfile_Header_Impl.cpp +++ b/Source/Diagnostics/BTD_Plotfile_Header_Impl.cpp @@ -1,3 +1,9 @@ +/* Copyright 2021 Revathi Jambunathan + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ #include "BTD_Plotfile_Header_Impl.H" #include "WarpX.H" @@ -341,3 +347,182 @@ BTDMultiFabHeaderImpl::CopyVec(amrex::Vector<amrex::Real>& dst, dst[i] = src[i]; } } + + +BTDSpeciesHeaderImpl::BTDSpeciesHeaderImpl (std::string const & Headerfile_path, std::string const& species_name) + : m_Header_path(Headerfile_path), m_species_name(species_name) +{ + +} + +void +BTDSpeciesHeaderImpl::ReadHeader () +{ + amrex::Vector<char> HeaderCharPtr; + amrex::Long fileLength(0), fileLengthPadded(0); + std::ifstream iss; + + iss.open(m_Header_path.c_str(), std::ios::in); + iss.seekg(0, std::ios::end); + fileLength = static_cast<std::streamoff>(iss.tellg()); + iss.seekg(0, std::ios::beg); + + fileLengthPadded = fileLength + 1; + HeaderCharPtr.resize(fileLengthPadded); + iss.read(HeaderCharPtr.dataPtr(), fileLength); + iss.close(); + HeaderCharPtr[fileLength] = '\0'; + + std::istringstream is(HeaderCharPtr.dataPtr(), std::istringstream::in); + is.exceptions(std::ios_base::failbit | std::ios_base::badbit); + + is >> m_file_version; + is >> m_spacedim; + is >> m_num_output_real; + m_real_comp_names.resize(m_num_output_real); + for (int i = 0; i < m_real_comp_names.size(); ++i) { + is >> m_real_comp_names[i]; + } + is >> m_num_output_int; + for (int i = 0; i < m_int_comp_names.size(); ++i) { + is >> m_int_comp_names[i]; + } + is >> m_is_checkpoint; + is >> m_total_particles; + is >> m_nextid; + is >> m_finestLevel; + m_particleBoxArray_size.resize(m_finestLevel+1); + m_which_data.resize(m_finestLevel+1); + m_particles_per_box.resize(m_finestLevel+1); + m_offset_per_box.resize(m_finestLevel+1); + for (int lev = 0; lev <= m_finestLevel; ++lev) { + is >> m_particleBoxArray_size[lev]; + m_which_data[lev].resize(m_particleBoxArray_size[lev]); + m_particles_per_box[lev].resize(m_particleBoxArray_size[lev]); + m_offset_per_box[lev].resize(m_particleBoxArray_size[lev]); + for (int i = 0; i < m_particleBoxArray_size[lev]; ++i) { + is >> m_which_data[lev][i] >> m_particles_per_box[lev][i] >> m_offset_per_box[lev][i]; + } + } + +} + +void +BTDSpeciesHeaderImpl::WriteHeader () +{ + if (amrex::FileExists(m_Header_path)) { + amrex::FileSystem::Remove(m_Header_path); + } + std::ofstream HeaderFile; + HeaderFile.open(m_Header_path.c_str(), std::ofstream::out | + std::ofstream::trunc | + std::ofstream::binary); + if ( !HeaderFile.good()) amrex::FileOpenFailed(m_Header_path); + + HeaderFile.precision(17); + + // File Version + HeaderFile << m_file_version << '\n'; + // spacedim + HeaderFile << m_spacedim << '\n'; + // number of real components output + HeaderFile << m_num_output_real << '\n'; + for (int i = 0; i < m_num_output_real; ++i) { + HeaderFile << m_real_comp_names[i] << '\n'; + } + HeaderFile << m_num_output_int << '\n'; + for (int i = 0; i < m_num_output_int; ++i) { + HeaderFile << m_int_comp_names[i] << '\n'; + } + HeaderFile << m_is_checkpoint << '\n'; + HeaderFile << m_total_particles << '\n'; + HeaderFile << m_nextid << '\n'; + HeaderFile << m_finestLevel << '\n'; + for (int lev = 0; lev <= m_finestLevel; ++lev) { + HeaderFile << m_particleBoxArray_size[lev] << '\n'; + for (int i = 0; i < m_particleBoxArray_size[lev]; ++i) { + HeaderFile << m_which_data[lev][i] << ' ' << m_particles_per_box[lev][i] << ' ' << m_offset_per_box[lev][i] << '\n'; + } + } +} + +void +BTDSpeciesHeaderImpl::set_DataIndex(const int lev, const int box_id, const int data_index) +{ + m_which_data[lev][box_id] = data_index; +} + +void +BTDSpeciesHeaderImpl::AppendParticleInfoForNewBox ( const int data_index, + const int particles_per_box, + const int offset) +{ + m_which_data[m_finestLevel].resize(m_particleBoxArray_size[m_finestLevel]); + m_particles_per_box[m_finestLevel].resize(m_particleBoxArray_size[m_finestLevel]); + m_offset_per_box[m_finestLevel].resize(m_particleBoxArray_size[m_finestLevel]); + + const int last_boxId = m_particleBoxArray_size[m_finestLevel] - 1; + m_which_data[m_finestLevel][last_boxId] = data_index; + m_particles_per_box[m_finestLevel][last_boxId] = particles_per_box; + m_offset_per_box[m_finestLevel][last_boxId] = offset; + +} + +BTDParticleDataHeaderImpl::BTDParticleDataHeaderImpl (std::string const & Headerfile_path) + : m_Header_path(Headerfile_path) +{ + +} + +void +BTDParticleDataHeaderImpl::ReadHeader () +{ + // Read existing fab Header first + amrex::Vector<char> HeaderCharPtr; + amrex::Long fileLength(0), fileLengthPadded(0); + std::ifstream iss; + + iss.open(m_Header_path.c_str(), std::ios::in); + iss.seekg(0, std::ios::end); + fileLength = static_cast<std::streamoff>(iss.tellg()); + iss.seekg(0, std::ios::beg); + + fileLengthPadded = fileLength + 1; + HeaderCharPtr.resize(fileLengthPadded); + iss.read(HeaderCharPtr.dataPtr(), fileLength); + iss.close(); + HeaderCharPtr[fileLength] = '\0'; + + std::istringstream is(HeaderCharPtr.dataPtr(), std::istringstream::in); + is.exceptions(std::ios_base::failbit | std::ios_base::badbit); + + + int in_hash; + int bl_ignore_max = 100000; + + is.ignore(bl_ignore_max,'(') >> m_ba_size >> in_hash; + m_ba.resize(m_ba_size); + for (int ibox = 0; ibox < m_ba.size(); ++ibox) { + amrex::Box bx; + is >> bx; + m_ba.set(ibox, bx); + } + is.ignore(bl_ignore_max, ')'); + +} + +void +BTDParticleDataHeaderImpl::WriteHeader () +{ + if (amrex::FileExists(m_Header_path)) { + amrex::FileSystem::Remove(m_Header_path); + } + std::ofstream HeaderFile; + HeaderFile.open(m_Header_path.c_str(), std::ofstream::out | + std::ofstream::trunc | + std::ofstream::binary); + if ( !HeaderFile.good()) amrex::FileOpenFailed(m_Header_path); + + HeaderFile.precision(17); + m_ba.writeOn(HeaderFile); +} diff --git a/Source/Diagnostics/BTDiagnostics.H b/Source/Diagnostics/BTDiagnostics.H index 629fe0192..a016ec777 100644 --- a/Source/Diagnostics/BTDiagnostics.H +++ b/Source/Diagnostics/BTDiagnostics.H @@ -1,3 +1,9 @@ +/* Copyright 2021 Revathi Jambunathan + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ #ifndef WARPX_BTDIAGNOSTICS_H_ #define WARPX_BTDIAGNOSTICS_H_ @@ -69,7 +75,13 @@ private: * This is currently an empty function: * The particle containers required for this must be added to populate this function. */ - void InitializeParticleBuffer () override {} + void InitializeParticleBuffer () override; + /** This function prepares the current z coordinate in the boosted-frame and lab-frame as required by + the particle and fields. + */ + void PrepareBufferData () override; + /** This function increments the buffer counter and identifies if the snapshot is fully populated */ + void UpdateBufferData () override; /** The cell-centered data for all fields, namely, * Ex, Ey, Ez, Bx, By, Bz, jx, jy, jz, and rho is computed and stored in * the multi-level cell-centered multifab, m_mf_cc. This MultiFab extends @@ -78,6 +90,8 @@ private: * is sliced, back-transformed, and stored in the output multifab, m_mf_output. */ void PrepareFieldDataForOutput () override; + /** The Particle Geometry, BoxArray, and RealBox are set for the lab-frame output */ + void PrepareParticleDataForOutput() override; /** whether z-slice that corresponds to the buffer, i_buffer, is within the * boosted-domain and lab-frame domains at level, lev. @@ -92,9 +106,9 @@ private: /** Initialize buffer domain, buffer box and lab-frame parameters such as * m_t_lab, and z-positions for the i^th snapshot, i_buffer, and level, lev. * \param[in] i_buffer i^th snapshot or buffer - * \param[in] lev mesh-refinement level for which the field buffer data is initialized. + * \param[in] lev mesh-refinement level for which the field and/or particle buffer data is initialized. */ - void InitializeFieldBufferData ( int i_buffer , int lev) override; + void InitializeBufferData ( int i_buffer , int lev) override; /** Whether to compute back-tranformed values for field-data. * default value is true. */ @@ -158,6 +172,9 @@ private: /** Vector of boosted-frame z co-ordinate corresponding to each back-transformed snapshot at the current timestep */ amrex::Vector<amrex::Real> m_current_z_boost; + /** Vector of previous boosted-frame z co-ordinate corresponding to each + back-transformed snapshot*/ + amrex::Vector<amrex::Real> m_old_z_boost; /** Geometry that defines the domain attributes corresponding to the * full snapshot in the back-transformed lab-frame. * Specifically, the user-defined physical co-ordinates for the diagnostics @@ -299,9 +316,6 @@ private: "Bx", "By", "Bz", "jx", "jy", "jz", "rho"}; - /** Temporarily clear species output for BTD until particle buffer is added */ - void TMP_ClearSpeciesDataForBTD() override; - /** Merge the lab-frame buffer multifabs so it can be visualized as * a single plotfile */ @@ -317,6 +331,28 @@ private: void InterleaveFabArrayHeader( std::string Buffer_FabHeaderFilename, std::string snapshot_FabHeaderFilename, std::string newsnapshot_FabFilename); -}; + /** Interleave lab-frame metadata of the species header file in the buffers to + * be consistent with the merged plotfile lab-frame data + */ + void InterleaveSpeciesHeader(std::string buffer_species_Header_path, + std::string snapshot_species_Header_path, + std::string species_name, const int new_data_index); + /** Interleave lab-frame metadata of the particle header file in the buffers to + * be consistent with the merged plotfile lab-frame data + */ + void InterleaveParticleDataHeader( std::string buffer_ParticleHdrFilename, + std::string snapshot_ParticleHdrFilename); + /** Initialize particle functors for each species to compute the back-transformed + lab-frame data. */ + void InitializeParticleFunctors () override; + + /** Update total number of particles flushed for all species for ith snapshot */ + void UpdateTotalParticlesFlushed(int i_buffer); + /** Reset total number of particles in the particle buffer to 0 for ith snapshot */ + void ResetTotalParticlesInBuffer(int i_buffer); + /** Clear particle data stored in the particle buffer */ + void ClearParticleBuffer(int i_buffer); + +}; #endif // WARPX_BTDIAGNOSTICS_H_ diff --git a/Source/Diagnostics/BTDiagnostics.cpp b/Source/Diagnostics/BTDiagnostics.cpp index b679cb524..7b61a6adb 100644 --- a/Source/Diagnostics/BTDiagnostics.cpp +++ b/Source/Diagnostics/BTDiagnostics.cpp @@ -1,3 +1,10 @@ +/* Copyright 2021 Revathi Jambunathan + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + #include "BTDiagnostics.H" #include "BTD_Plotfile_Header_Impl.H" #include "ComputeDiagFunctors/BackTransformFunctor.H" @@ -7,6 +14,7 @@ #include "Diagnostics/Diagnostics.H" #include "Diagnostics/FlushFormats/FlushFormat.H" #include "Parallelization/WarpXCommUtil.H" +#include "ComputeDiagFunctors/BackTransformParticleFunctor.H" #include "Utils/CoarsenIO.H" #include "Utils/WarpXConst.H" #include "Utils/WarpXUtil.H" @@ -51,53 +59,60 @@ void BTDiagnostics::DerivedInitData () // The number of levels to be output is nlev_output. nlev_output = 1; - // allocate vector of m_t_lab with m_num_buffers; m_t_lab.resize(m_num_buffers); - // allocate vector of RealBost of the simulation domain in lab-frame m_prob_domain_lab.resize(m_num_buffers); - // allocate vector of RealBox of the diag domain m_snapshot_domain_lab.resize(m_num_buffers); - // allocate vector of RealBox of the buffers that fill the snapshot m_buffer_domain_lab.resize(m_num_buffers); - // define box correctly (one for all snapshots) m_snapshot_box.resize(m_num_buffers); - // define box for each buffer that fills the snapshot m_buffer_box.resize(m_num_buffers); - // allocate vector of m_current_z_lab m_current_z_lab.resize(m_num_buffers); - // allocate vector of m_num_buffers m_current_z_boost.resize(m_num_buffers); - // allocate vector of m_buff_counter to counter number of slices filled in the buffer + m_old_z_boost.resize(m_num_buffers); m_buffer_counter.resize(m_num_buffers); - // allocate vector of num_Cells in the lab-frame m_snapshot_ncells_lab.resize(m_num_buffers); - // allocate vector of cell centered multifabs for nlevels m_cell_centered_data.resize(nmax_lev); - // allocate vector of cell-center functors for nlevels m_cell_center_functors.resize(nmax_lev); - // allocate vector to estimate maximum number of buffer multifabs needed to - // obtain the lab-frame snapshot. m_max_buffer_multifabs.resize(m_num_buffers); - // allocate vector to count number of times the buffer multifab - // has been flushed and refilled m_buffer_flush_counter.resize(m_num_buffers); - // allocate vector of geometry objects corresponding to each snapshot m_geom_snapshot.resize( m_num_buffers ); m_snapshot_full.resize( m_num_buffers ); m_lastValidZSlice.resize( m_num_buffers ); for (int i = 0; i < m_num_buffers; ++i) { m_geom_snapshot[i].resize(nmax_lev); - // initialize snapshot full boolean to false m_snapshot_full[i] = 0; m_lastValidZSlice[i] = 0; } - for (int lev = 0; lev < nmax_lev; ++lev) { // Define cell-centered multifab over the whole domain with // user-defined crse_ratio for nlevels DefineCellCenteredMultiFab(lev); } + /* Allocate vector of particle buffer vectors for each snapshot */ + MultiParticleContainer& mpc = warpx.GetPartContainer(); + // If not specified, and write species is not 0, dump all species + amrex::ParmParse pp_diag_name(m_diag_name); + int write_species = 1; + pp_diag_name.query("write_species", write_species); + if (m_output_species_names.size() == 0 and write_species == 1) + m_output_species_names = mpc.GetSpeciesNames(); + + if (m_output_species_names.size() > 0) { + m_do_back_transformed_particles = true; + } else { + m_do_back_transformed_particles = false; + } + // Turn on do_back_transformed_particles in the particle containers so that + // the tmp_particle_data is allocated and the data of the corresponding species is + // copied and stored in tmp_particle_data before particles are pushed. + for (auto const& species : m_output_species_names){ + mpc.SetDoBackTransformedParticles(m_do_back_transformed_particles); + mpc.SetDoBackTransformedParticles(species, m_do_back_transformed_particles); + } + m_particles_buffer.resize(m_num_buffers); + m_output_species.resize(m_num_buffers); + m_totalParticles_flushed_already.resize(m_num_buffers); + m_totalParticles_in_buffer.resize(m_num_buffers); } void @@ -135,6 +150,8 @@ BTDiagnostics::ReadParameters () pp_diag_name.query("do_back_transformed_fields", m_do_back_transformed_fields); pp_diag_name.query("do_back_transformed_particles", m_do_back_transformed_particles); AMREX_ALWAYS_ASSERT(m_do_back_transformed_fields or m_do_back_transformed_particles); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_do_back_transformed_fields, " fields must be turned on for the new back-transformed diagnostics"); + if (m_do_back_transformed_fields == false) m_varnames.clear(); getWithParser(pp_diag_name, "num_snapshots_lab", m_num_snapshots_lab); m_num_buffers = m_num_snapshots_lab; @@ -181,15 +198,24 @@ BTDiagnostics::DoDump (int step, int i_buffer, bool force_flush) bool -BTDiagnostics::DoComputeAndPack (int step, bool /*force_flush*/) +BTDiagnostics::DoComputeAndPack (int step, bool force_flush) { // always set to true for BTDiagnostics since back-transform buffers are potentially - // computed and packed every timstep, except at initialization when step == -1. - return (step>=0); + // computed and packed every timstep, except at initialization when step == -1, or when + // force_flush is set to true, because we dont need to redundantly re-compute + // buffers when force_flush = true. We only need to dump the buffers when + // force_flush=true. Note that the BTD computation is performed every timestep (step>=0) + if (step < 0 ) { + return false; + } else if (force_flush) { + return false; + } else { + return true; + } } void -BTDiagnostics::InitializeFieldBufferData ( int i_buffer , int lev) +BTDiagnostics::InitializeBufferData ( int i_buffer , int lev) { auto & warpx = WarpX::GetInstance(); // Lab-frame time for the i^th snapshot @@ -281,6 +307,8 @@ BTDiagnostics::InitializeFieldBufferData ( int i_buffer , int lev) m_buffer_counter[i_buffer] = 0; m_current_z_lab[i_buffer] = 0._rt; m_current_z_boost[i_buffer] = 0._rt; + // store old z boost before updated zboost position + m_old_z_boost[i_buffer] = m_current_z_boost[i_buffer]; // Now Update Current Z Positions m_current_z_boost[i_buffer] = UpdateCurrentZBoostCoordinate(m_t_lab[i_buffer], warpx.gett_new(lev) ); @@ -329,6 +357,7 @@ BTDiagnostics::InitializeFieldBufferData ( int i_buffer , int lev) void BTDiagnostics::DefineCellCenteredMultiFab(int lev) { + if (m_do_back_transformed_fields == false) return; // Creating MultiFab to store cell-centered data in boosted-frame for the entire-domain // This MultiFab will store all the user-requested fields in the boosted-frame auto & warpx = WarpX::GetInstance(); @@ -345,6 +374,9 @@ BTDiagnostics::DefineCellCenteredMultiFab(int lev) void BTDiagnostics::InitializeFieldFunctors (int lev) { + // Initialize fields functors only if do_back_transformed_fields is selected + if (m_do_back_transformed_fields == false) return; + auto & warpx = WarpX::GetInstance(); // Clear any pre-existing vector to release stored data // This ensures that when domain is load-balanced, the functors point @@ -398,8 +430,54 @@ BTDiagnostics::InitializeFieldFunctors (int lev) } void +BTDiagnostics::PrepareBufferData () +{ + auto & warpx = WarpX::GetInstance(); + int num_BT_functors = 1; + + for (int lev = 0; lev < nlev_output; ++lev) + { + for (int i = 0; i < num_BT_functors; ++i) + { + for (int i_buffer = 0; i_buffer < m_num_buffers; ++i_buffer ) + { + m_old_z_boost[i_buffer] = m_current_z_boost[i_buffer]; + // Update z-boost and z-lab positions + m_current_z_boost[i_buffer] = UpdateCurrentZBoostCoordinate(m_t_lab[i_buffer], + warpx.gett_new(lev) ); + m_current_z_lab[i_buffer] = UpdateCurrentZLabCoordinate(m_t_lab[i_buffer], + warpx.gett_new(lev) ); + } + } + } +} + +void +BTDiagnostics::UpdateBufferData () +{ + int num_BT_functors = 1; + + for (int lev = 0; lev < nlev_output; ++lev) + { + for (int i = 0; i < num_BT_functors; ++i) + { + for (int i_buffer = 0; i_buffer < m_num_buffers; ++i_buffer ) + { + bool ZSliceInDomain = GetZSliceInDomainFlag (i_buffer, lev); + if (ZSliceInDomain) ++m_buffer_counter[i_buffer]; + // when the 0th z-index is filled, then set lastValidZSlice to 1 + if (k_index_zlab(i_buffer, lev) == 0) m_lastValidZSlice[i_buffer] = 1; + } + } + } +} + +void BTDiagnostics::PrepareFieldDataForOutput () { + // Initialize fields functors only if do_back_transformed_fields is selected + if (m_do_back_transformed_fields == false) return; + auto & warpx = WarpX::GetInstance(); // In this function, we will get cell-centered data for every level, lev, // using the cell-center functors and their respective opeators() @@ -434,11 +512,6 @@ BTDiagnostics::PrepareFieldDataForOutput () { for (int i_buffer = 0; i_buffer < m_num_buffers; ++i_buffer ) { - // Update z-boost and z-lab positions - m_current_z_boost[i_buffer] = UpdateCurrentZBoostCoordinate(m_t_lab[i_buffer], - warpx.gett_new(lev) ); - m_current_z_lab[i_buffer] = UpdateCurrentZLabCoordinate(m_t_lab[i_buffer], - warpx.gett_new(lev) ); // Check if the zslice is in domain bool ZSliceInDomain = GetZSliceInDomainFlag (i_buffer, lev); // Initialize and define field buffer multifab if buffer is empty @@ -459,9 +532,6 @@ BTDiagnostics::PrepareFieldDataForOutput () k_index_zlab(i_buffer, lev), m_max_box_size, m_snapshot_full[i_buffer] ); - if (ZSliceInDomain) ++m_buffer_counter[i_buffer]; - // when the 0th z-index is filled, then set lastValidZSlice to 1 - if (k_index_zlab(i_buffer, lev) == 0) m_lastValidZSlice[i_buffer] = 1; } } } @@ -616,7 +686,7 @@ BTDiagnostics::GetZSliceInDomainFlag (const int i_buffer, const int lev) amrex::Real buffer_zmax_lab = m_snapshot_domain_lab[i_buffer].hi( m_moving_window_dir ); if ( ( m_current_z_boost[i_buffer] < boost_domain.lo(m_moving_window_dir) ) or ( m_current_z_boost[i_buffer] > boost_domain.hi(m_moving_window_dir) ) or - ( m_current_z_lab[i_buffer] < buffer_zmin_lab ) or + ( m_current_z_lab[i_buffer] < buffer_zmin_lab ) or ( m_current_z_lab[i_buffer] > buffer_zmax_lab ) ) { // the slice is not in the boosted domain or lab-frame domain @@ -641,9 +711,10 @@ BTDiagnostics::Flush (int i_buffer) double const labtime = m_t_lab[i_buffer]; m_flush_format->WriteToFile( m_varnames, m_mf_output[i_buffer], m_geom_output[i_buffer], warpx.getistep(), - labtime, m_output_species, nlev_output, file_name, m_file_min_digits, + labtime, m_output_species[i_buffer], nlev_output, file_name, m_file_min_digits, m_plot_raw_fields, m_plot_raw_fields_guards, - isBTD, i_buffer, m_geom_snapshot[i_buffer][0], isLastBTDFlush); + isBTD, i_buffer, m_geom_snapshot[i_buffer][0], isLastBTDFlush, + m_totalParticles_flushed_already[i_buffer]); if (m_format == "plotfile") { MergeBuffersForPlotfile(i_buffer); @@ -652,12 +723,12 @@ BTDiagnostics::Flush (int i_buffer) // Reset the buffer counter to zero after flushing out data stored in the buffer. ResetBufferCounter(i_buffer); IncrementBufferFlushCounter(i_buffer); -} - -void BTDiagnostics::TMP_ClearSpeciesDataForBTD () -{ - m_output_species.clear(); - m_output_species_names.clear(); + // if particles are selected for output then update and reset counters + if (m_output_species_names.size() > 0) { + UpdateTotalParticlesFlushed(i_buffer); + ResetTotalParticlesInBuffer(i_buffer); + ClearParticleBuffer(i_buffer); + } } void BTDiagnostics::MergeBuffersForPlotfile (int i_snapshot) @@ -676,13 +747,6 @@ void BTDiagnostics::MergeBuffersForPlotfile (int i_snapshot) // BTD plotfile have only one level, Level0. std::string snapshot_Level0_path = snapshot_path + "/Level_0"; std::string snapshot_Header_filename = snapshot_path + "/Header"; - // Create directory only when the first buffer is flushed out. - if (m_buffer_flush_counter[i_snapshot] == 0 ) { - // Create Level_0 directory to store all Cell_D and Cell_H files - if (!amrex::UtilCreateDirectory(snapshot_Level0_path, 0755) ) - amrex::CreateDirectoryFailed(snapshot_Level0_path); - } - // Path of the buffer recently flushed std::string BufferPath_prefix = snapshot_path + "/buffer"; const std::string recent_Buffer_filepath = amrex::Concatenate(BufferPath_prefix,iteration[0]); @@ -690,38 +754,104 @@ void BTDiagnostics::MergeBuffersForPlotfile (int i_snapshot) std::string recent_Header_filename = recent_Buffer_filepath+"/Header"; std::string recent_Buffer_Level0_path = recent_Buffer_filepath + "/Level_0"; std::string recent_Buffer_FabHeaderFilename = recent_Buffer_Level0_path + "/Cell_H"; - // Read the header file to get the fab on disk string - BTDMultiFabHeaderImpl Buffer_FabHeader(recent_Buffer_FabHeaderFilename); - Buffer_FabHeader.ReadMultiFabHeader(); - if (Buffer_FabHeader.ba_size() > 1) amrex::Abort("BTD Buffer has more than one fabs."); - // Every buffer that is flushed only has a single fab. - std::string recent_Buffer_FabFilename = recent_Buffer_Level0_path + "/" - + Buffer_FabHeader.FabName(0); - // Existing snapshot Fab Header Filename - std::string snapshot_FabHeaderFilename = snapshot_Level0_path + "/Cell_H"; - std::string snapshot_FabFilename = amrex::Concatenate(snapshot_Level0_path+"/Cell_D_",m_buffer_flush_counter[i_snapshot], 5); - // Name of the newly appended fab in the snapshot - std::string new_snapshotFabFilename = amrex::Concatenate("Cell_D_",m_buffer_flush_counter[i_snapshot],5); - - if ( m_buffer_flush_counter[i_snapshot] == 0) { - std::rename(recent_Header_filename.c_str(), snapshot_Header_filename.c_str()); - Buffer_FabHeader.SetFabName(0, Buffer_FabHeader.fodPrefix(0), - new_snapshotFabFilename, - Buffer_FabHeader.FabHead(0)); - Buffer_FabHeader.WriteMultiFabHeader(); - std::rename(recent_Buffer_FabHeaderFilename.c_str(), - snapshot_FabHeaderFilename.c_str()); - std::rename(recent_Buffer_FabFilename.c_str(), - snapshot_FabFilename.c_str()); - } else { - // Interleave Header file - InterleaveBufferAndSnapshotHeader(recent_Header_filename, - snapshot_Header_filename); - InterleaveFabArrayHeader(recent_Buffer_FabHeaderFilename, - snapshot_FabHeaderFilename, - new_snapshotFabFilename); - std::rename(recent_Buffer_FabFilename.c_str(), - snapshot_FabFilename.c_str()); + // Create directory only when the first buffer is flushed out. + if (m_buffer_flush_counter[i_snapshot] == 0 ) { + // Create Level_0 directory to store all Cell_D and Cell_H files + if (!amrex::UtilCreateDirectory(snapshot_Level0_path, 0755) ) + amrex::CreateDirectoryFailed(snapshot_Level0_path); + // Create directory for each species selected for diagnostic + for (int i = 0; i < m_particles_buffer[i_snapshot].size(); ++i) { + std::string snapshot_species_path = snapshot_path + "/" + m_output_species_names[i]; + if ( !amrex::UtilCreateDirectory(snapshot_species_path, 0755)) + amrex::CreateDirectoryFailed(snapshot_species_path); + // Create Level_0 directory for particles to store Particle_H and DATA files + std::string species_Level0_path = snapshot_species_path + "/Level_0"; + if ( !amrex::UtilCreateDirectory(species_Level0_path, 0755)) + amrex::CreateDirectoryFailed(species_Level0_path); + } + std::string buffer_WarpXHeader_path = recent_Buffer_filepath + "/WarpXHeader"; + std::string snapshot_WarpXHeader_path = snapshot_path + "/WarpXHeader"; + std::string buffer_job_info_path = recent_Buffer_filepath + "/warpx_job_info"; + std::string snapshot_job_info_path = snapshot_path + "/warpx_job_info"; + std::rename(buffer_WarpXHeader_path.c_str(), snapshot_WarpXHeader_path.c_str()); + std::rename(buffer_job_info_path.c_str(), snapshot_job_info_path.c_str()); + } + + if (m_do_back_transformed_fields == true) { + // Read the header file to get the fab on disk string + BTDMultiFabHeaderImpl Buffer_FabHeader(recent_Buffer_FabHeaderFilename); + Buffer_FabHeader.ReadMultiFabHeader(); + if (Buffer_FabHeader.ba_size() > 1) amrex::Abort("BTD Buffer has more than one fabs."); + // Every buffer that is flushed only has a single fab. + std::string recent_Buffer_FabFilename = recent_Buffer_Level0_path + "/" + + Buffer_FabHeader.FabName(0); + // Existing snapshot Fab Header Filename + std::string snapshot_FabHeaderFilename = snapshot_Level0_path + "/Cell_H"; + std::string snapshot_FabFilename = amrex::Concatenate(snapshot_Level0_path+"/Cell_D_",m_buffer_flush_counter[i_snapshot], 5); + // Name of the newly appended fab in the snapshot + std::string new_snapshotFabFilename = amrex::Concatenate("Cell_D_",m_buffer_flush_counter[i_snapshot],5); + + if ( m_buffer_flush_counter[i_snapshot] == 0) { + std::rename(recent_Header_filename.c_str(), snapshot_Header_filename.c_str()); + Buffer_FabHeader.SetFabName(0, Buffer_FabHeader.fodPrefix(0), + new_snapshotFabFilename, + Buffer_FabHeader.FabHead(0)); + Buffer_FabHeader.WriteMultiFabHeader(); + std::rename(recent_Buffer_FabHeaderFilename.c_str(), + snapshot_FabHeaderFilename.c_str()); + std::rename(recent_Buffer_FabFilename.c_str(), + snapshot_FabFilename.c_str()); + } else { + // Interleave Header file + InterleaveBufferAndSnapshotHeader(recent_Header_filename, + snapshot_Header_filename); + InterleaveFabArrayHeader(recent_Buffer_FabHeaderFilename, + snapshot_FabHeaderFilename, + new_snapshotFabFilename); + std::rename(recent_Buffer_FabFilename.c_str(), + snapshot_FabFilename.c_str()); + } + } + for (int i = 0; i < m_particles_buffer[i_snapshot].size(); ++i) { + // species filename of recently flushed buffer + std::string recent_species_prefix = recent_Buffer_filepath+"/"+m_output_species_names[i]; + std::string recent_species_Header = recent_species_prefix + "/Header"; + std::string recent_ParticleHdrFilename = recent_species_prefix + "/Level_0/Particle_H"; + BTDSpeciesHeaderImpl BufferSpeciesHeader(recent_species_Header, + m_output_species_names[i]); + BufferSpeciesHeader.ReadHeader(); + // only one box is flushed out at a time + std::string recent_ParticleDataFilename = amrex::Concatenate(recent_species_prefix + "/Level_0/DATA_",BufferSpeciesHeader.m_which_data[0][0]); + // Path to snapshot particle files + std::string snapshot_species_path = snapshot_path + "/" + m_output_species_names[i]; + std::string snapshot_species_Level0path = snapshot_species_path + "/Level_0"; + std::string snapshot_species_Header = snapshot_species_path + "/Header"; + std::string snapshot_ParticleHdrFilename = snapshot_species_Level0path + "/Particle_H"; + std::string snapshot_ParticleDataFilename = amrex::Concatenate(snapshot_species_Level0path + "/DATA_",m_buffer_flush_counter[i_snapshot],5); + + + if (m_buffer_flush_counter[i_snapshot] == 0) { + BufferSpeciesHeader.set_DataIndex(0,0,m_buffer_flush_counter[i_snapshot]); + BufferSpeciesHeader.WriteHeader(); + + // copy Header file for the species + std::rename(recent_species_Header.c_str(), snapshot_species_Header.c_str()); + if (BufferSpeciesHeader.m_total_particles == 0) continue; + // if finite number of particles in the output, copy ParticleHdr and Data file + std::rename(recent_ParticleHdrFilename.c_str(), snapshot_ParticleHdrFilename.c_str()); + std::rename(recent_ParticleDataFilename.c_str(), snapshot_ParticleDataFilename.c_str()); + } else { + InterleaveSpeciesHeader(recent_species_Header,snapshot_species_Header, + m_output_species_names[i], m_buffer_flush_counter[i_snapshot]); + if (BufferSpeciesHeader.m_total_particles == 0) continue; + if (m_totalParticles_flushed_already[i_snapshot][i]==0) { + std::rename(recent_ParticleHdrFilename.c_str(), snapshot_ParticleHdrFilename.c_str()); + } else { + InterleaveParticleDataHeader(recent_ParticleHdrFilename, + snapshot_ParticleHdrFilename); + } + std::rename(recent_ParticleDataFilename.c_str(), snapshot_ParticleDataFilename.c_str()); + } } // Destroying the recently flushed buffer directory since it is already merged. amrex::FileSystem::RemoveAll(recent_Buffer_filepath); @@ -804,3 +934,164 @@ BTDiagnostics::InterleaveFabArrayHeader(std::string Buffer_FabHeader_path, snapshot_FabHeader.WriteMultiFabHeader(); } + +void +BTDiagnostics::InterleaveSpeciesHeader(std::string buffer_species_Header_path, + std::string snapshot_species_Header_path, + std::string species_name, const int new_data_index) +{ + BTDSpeciesHeaderImpl BufferSpeciesHeader(buffer_species_Header_path, + species_name); + BufferSpeciesHeader.ReadHeader(); + + BTDSpeciesHeaderImpl SnapshotSpeciesHeader(snapshot_species_Header_path, + species_name); + SnapshotSpeciesHeader.ReadHeader(); + SnapshotSpeciesHeader.AddTotalParticles( BufferSpeciesHeader.m_total_particles); + + SnapshotSpeciesHeader.IncrementParticleBoxArraySize(); + const int buffer_finestLevel = BufferSpeciesHeader.m_finestLevel; + const int buffer_boxId = BufferSpeciesHeader.m_particleBoxArray_size[buffer_finestLevel]-1; + SnapshotSpeciesHeader.AppendParticleInfoForNewBox( + new_data_index, + BufferSpeciesHeader.m_particles_per_box[buffer_finestLevel][buffer_boxId], + BufferSpeciesHeader.m_offset_per_box[buffer_finestLevel][buffer_boxId]); + SnapshotSpeciesHeader.WriteHeader(); +} + +void +BTDiagnostics::InterleaveParticleDataHeader(std::string buffer_ParticleHdrFilename, + std::string snapshot_ParticleHdrFilename) +{ + BTDParticleDataHeaderImpl BufferParticleHeader(buffer_ParticleHdrFilename); + BufferParticleHeader.ReadHeader(); + + BTDParticleDataHeaderImpl SnapshotParticleHeader(snapshot_ParticleHdrFilename); + SnapshotParticleHeader.ReadHeader(); + + // Increment BoxArraySize + SnapshotParticleHeader.IncreaseBoxArraySize( BufferParticleHeader.ba_size() ); + // Append New box in snapshot + for (int ibox = 0; ibox < BufferParticleHeader.ba_size(); ++ibox) { + int new_ibox = SnapshotParticleHeader.ba_size() - 1 + ibox; + SnapshotParticleHeader.ResizeBoxArray(); + SnapshotParticleHeader.SetBox(new_ibox, BufferParticleHeader.ba_box(ibox) ); + } + SnapshotParticleHeader.WriteHeader(); +} + +void +BTDiagnostics::InitializeParticleFunctors () +{ + auto & warpx = WarpX::GetInstance(); + const MultiParticleContainer& mpc = warpx.GetPartContainer(); + // allocate with total number of species diagnostics + m_all_particle_functors.resize(m_output_species_names.size()); + // Create an object of class BackTransformParticleFunctor + for (int i = 0; i < m_all_particle_functors.size(); ++i) + { + // species id corresponding to ith diag species + const int idx = mpc.getSpeciesID(m_output_species_names[i]); + m_all_particle_functors[i] = std::make_unique<BackTransformParticleFunctor>(mpc.GetParticleContainerPtr(idx), m_output_species_names[i], m_num_buffers); + } + +} + +void +BTDiagnostics::InitializeParticleBuffer () +{ + auto& warpx = WarpX::GetInstance(); + const MultiParticleContainer& mpc = warpx.GetPartContainer(); + for (int i = 0; i < m_num_buffers; ++i) { + m_particles_buffer[i].resize(m_output_species_names.size()); + m_totalParticles_flushed_already[i].resize(m_output_species_names.size()); + m_totalParticles_in_buffer[i].resize(m_output_species_names.size()); + for (int isp = 0; isp < m_particles_buffer[i].size(); ++isp) { + m_totalParticles_flushed_already[i][isp] = 0; + m_totalParticles_in_buffer[i][isp] = 0; + m_particles_buffer[i][isp] = std::make_unique<PinnedMemoryParticleContainer>(&WarpX::GetInstance()); + const int idx = mpc.getSpeciesID(m_output_species_names[isp]); + m_output_species[i].push_back(ParticleDiag(m_diag_name, + m_output_species_names[isp], + mpc.GetParticleContainerPtr(idx), + m_particles_buffer[i][isp].get())); + } + } + +} + +void +BTDiagnostics::PrepareParticleDataForOutput() +{ + + auto& warpx = WarpX::GetInstance(); + for (int lev = 0; lev < nlev_output; ++lev) { + for (int i = 0; i < m_all_particle_functors.size(); ++i) + { + for (int i_buffer = 0; i_buffer < m_num_buffers; ++i_buffer ) + { + // Check if the zslice is in domain + bool ZSliceInDomain = GetZSliceInDomainFlag (i_buffer, lev); + if (ZSliceInDomain) { + if ( m_totalParticles_in_buffer[i_buffer][i] == 0) { + amrex::Box particle_buffer_box = m_buffer_box[i_buffer]; + particle_buffer_box.setSmall(m_moving_window_dir, + m_buffer_box[i_buffer].smallEnd(m_moving_window_dir)-1); + particle_buffer_box.setBig(m_moving_window_dir, + m_buffer_box[i_buffer].bigEnd(m_moving_window_dir)+1); + amrex::BoxArray buffer_ba( particle_buffer_box ); + //amrex::BoxArray buffer_ba( particle_buffer_box ); + buffer_ba.maxSize(m_max_box_size); + amrex::DistributionMapping buffer_dmap(buffer_ba); + m_particles_buffer[i_buffer][i]->SetParticleBoxArray(lev, buffer_ba); + m_particles_buffer[i_buffer][i]->SetParticleDistributionMap(lev, buffer_dmap); + amrex::IntVect particle_DomBox_lo = m_snapshot_box[i_buffer].smallEnd(); + amrex::IntVect particle_DomBox_hi = m_snapshot_box[i_buffer].bigEnd(); + int zmin = std::max(0, particle_DomBox_lo[m_moving_window_dir] ); + particle_DomBox_lo[m_moving_window_dir] = zmin; + amrex::Box ParticleBox(particle_DomBox_lo, particle_DomBox_hi); + int num_cells = particle_DomBox_hi[m_moving_window_dir] - zmin + 1; + amrex::IntVect ref_ratio = amrex::IntVect(1); + amrex::Real new_lo = m_snapshot_domain_lab[i_buffer].hi(m_moving_window_dir) - + num_cells * dz_lab(warpx.getdt(lev), ref_ratio[m_moving_window_dir]); + amrex::RealBox ParticleRealBox = m_snapshot_domain_lab[i_buffer]; + ParticleRealBox.setLo(m_moving_window_dir, new_lo); + amrex::Vector<int> BTdiag_periodicity(AMREX_SPACEDIM, 0); + amrex::Geometry geom; + geom.define(ParticleBox, &ParticleRealBox, amrex::CoordSys::cartesian, + BTdiag_periodicity.data() ); + m_particles_buffer[i_buffer][i]->SetParticleGeometry(lev, geom); + } + } + m_all_particle_functors[i]->PrepareFunctorData ( + i_buffer, ZSliceInDomain, m_old_z_boost[i_buffer], + m_current_z_boost[i_buffer], m_t_lab[i_buffer], + m_snapshot_full[i_buffer]); + } + } + } +} + +void +BTDiagnostics::UpdateTotalParticlesFlushed(int i_buffer) +{ + for (int isp = 0; isp < m_totalParticles_flushed_already[i_buffer].size(); ++isp) { + m_totalParticles_flushed_already[i_buffer][isp] += m_totalParticles_in_buffer[i_buffer][isp]; + } +} + +void +BTDiagnostics::ResetTotalParticlesInBuffer(int i_buffer) +{ + for (int isp = 0; isp < m_totalParticles_in_buffer[i_buffer].size(); ++isp) { + m_totalParticles_in_buffer[i_buffer][isp] = 0; + } +} + +void +BTDiagnostics::ClearParticleBuffer(int i_buffer) +{ + for (int isp = 0; isp < m_particles_buffer[i_buffer].size(); ++isp) { + m_particles_buffer[i_buffer][isp]->clearParticles(); + } +} diff --git a/Source/Diagnostics/ComputeDiagFunctors/BackTransformFunctor.H b/Source/Diagnostics/ComputeDiagFunctors/BackTransformFunctor.H index a09861225..35c22908c 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/BackTransformFunctor.H +++ b/Source/Diagnostics/ComputeDiagFunctors/BackTransformFunctor.H @@ -1,3 +1,9 @@ +/* Copyright 2021 Revathi Jambunathan + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ #ifndef WARPX_BACKTRANSFORMFUNCTOR_H_ #define WARPX_BACKTRANSFORMFUNCTOR_H_ diff --git a/Source/Diagnostics/ComputeDiagFunctors/BackTransformFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/BackTransformFunctor.cpp index 7eb50aba1..6718c5833 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/BackTransformFunctor.cpp +++ b/Source/Diagnostics/ComputeDiagFunctors/BackTransformFunctor.cpp @@ -1,3 +1,9 @@ +/* Copyright 2021 Revathi Jambunathan + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ #include "BackTransformFunctor.H" #include "Diagnostics/ComputeDiagFunctors/ComputeDiagFunctor.H" @@ -44,13 +50,18 @@ BackTransformFunctor::operator ()(amrex::MultiFab& mf_dst, int /*dcomp*/, const amrex::Real gamma_boost = warpx.gamma_boost; int moving_window_dir = warpx.moving_window_dir; amrex::Real beta_boost = std::sqrt( 1._rt - 1._rt/( gamma_boost * gamma_boost) ); - bool interpolate = true; + const bool interpolate = true; std::unique_ptr< amrex::MultiFab > slice = nullptr; int scomp = 0; // Generate slice of the cell-centered multifab containing boosted-frame field-data // at current z-boost location for the ith buffer - slice = amrex::get_slice_data (moving_window_dir, m_current_z_boost[i_buffer], - *m_mf_src, geom, scomp, m_mf_src->nComp(), interpolate); + slice = amrex::get_slice_data(moving_window_dir, + m_current_z_boost[i_buffer], + *m_mf_src, + geom, + scomp, + m_mf_src->nComp(), + interpolate); // Perform in-place Lorentz-transform of all the fields stored in the slice. LorentzTransformZ( *slice, gamma_boost, beta_boost); diff --git a/Source/Diagnostics/ComputeDiagFunctors/BackTransformParticleFunctor.H b/Source/Diagnostics/ComputeDiagFunctors/BackTransformParticleFunctor.H new file mode 100644 index 000000000..59ceb12cc --- /dev/null +++ b/Source/Diagnostics/ComputeDiagFunctors/BackTransformParticleFunctor.H @@ -0,0 +1,237 @@ +/* Copyright 2021 Revathi Jambunathan + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef WARPX_BACKTRANSFORMPARTICLEFUNCTOR_H_ +#define WARPX_BACKTRANSFORMPARTICLEFUNCTOR_H_ + +#include "ComputeParticleDiagFunctor.H" +#include "Particles/Pusher/GetAndSetPosition.H" + +#include <AMReX.H> +#include <AMReX_AmrParticles.H> + + + + +/** + * \brief Filter to select particles that correspond to a z-slice of the corresponding lab-frame. + */ +struct SelectParticles +{ + using TmpParticles = WarpXParticleContainer::TmpParticles; + + /** + * \brief Constructor of SelectParticles functor. + * + * @param[in] a_pti WarpX particle iterator + * @param[in] tmp_particle_data temporary particle data + * @param[in] current_z_boost current z-position of the slice in boosted frame + * @param[in] old_z_boost previous z-position of the slice in boosted frame + */ + SelectParticles( const WarpXParIter& a_pti, TmpParticles& tmp_particle_data, + amrex::Real current_z_boost, amrex::Real old_z_boost, + int a_offset = 0); + + /** + * \brief Functor call. This method determines if a given particle should be selected + * for Lorentz transformation in obtaining the lab-frame data. The particles that + * with positions that correspond to the specific z-slice in boosted frame are selected. + * + * @param[in] SrcData particle tile data + * @param[in] i particle index + * @return 1 if particles is selected for transformation, else 0 + */ + template <typename SrcData> + AMREX_GPU_HOST_DEVICE + int operator() (const SrcData& src, int i) const noexcept + { + amrex::ignore_unused(src); + amrex::ParticleReal xp, yp, zp; + m_get_position(i, xp, yp, zp); + int Flag = 0; + if ( ( (zp >= m_current_z_boost) && (zpold[i] <= m_old_z_boost) ) || + ( (zp <= m_current_z_boost) && (zpold[i] >= m_old_z_boost) )) + { Flag = 1; + } + return Flag; + } + + /** Object to extract the positions of the macroparticles inside a ParallelFor kernel */ + GetParticlePosition m_get_position; + /** Z coordinate in boosted frame that corresponds to a give snapshot*/ + amrex::Real m_current_z_boost; + /** Previous Z coordinate in boosted frame that corresponds to a give snapshot*/ + amrex::Real m_old_z_boost; + /** Particle z coordinate in boosted frame*/ + amrex::ParticleReal* AMREX_RESTRICT zpold = nullptr; +}; + +/** + * \brief Transform functor to Lorentz-transform particles and obtain lab-frame data. + */ +struct LorentzTransformParticles +{ + + using TmpParticles = WarpXParticleContainer::TmpParticles; + + /** + * \brief Constructor of the LorentzTransformParticles functor. + * + * @param[in] a_pti WarpX particle iterator + * @param[in] tmp_particle_data temporary particle data + * @param[in] t_boost time in boosted frame + * @param[in] dt timestep in boosted-frame + * @param[in] t_lab time in lab-frame + */ + LorentzTransformParticles ( const WarpXParIter& a_pti, TmpParticles& tmp_particle_data, + amrex::Real t_boost, amrex::Real dt, + amrex::Real t_lab, int a_offset = 0); + + /** + * \brief Functor call. This method computes the Lorentz-transform for particle + * attributes to obtain the lab-frame snapshot data. + * + * @param[out] DstData particle tile data that stores the transformed particle data + * @param[in] SrcData particle tile data that is selected for transformation + * @param[in] i_src particle index of the source particles + * @param[in] i_dst particle index of the target particles (transformed data). + */ + template <typename DstData, typename SrcData> + AMREX_GPU_HOST_DEVICE + void operator () (const DstData& dst, const SrcData& src, int i_src, int i_dst) const noexcept + { + amrex::ignore_unused(src); + using namespace amrex::literals; + // get current src position + amrex::ParticleReal xpnew, ypnew, zpnew; + m_get_position(i_src, xpnew, ypnew, zpnew); + const amrex::Real gamma_new_p = std::sqrt(1.0_rt + m_inv_c2* + ( m_uxpnew[i_src] * m_uxpnew[i_src] + + m_uypnew[i_src] * m_uypnew[i_src] + + m_uzpnew[i_src] * m_uzpnew[i_src])); + const amrex::Real gamma_old_p = std::sqrt(1.0_rt + m_inv_c2* + ( m_uxpold[i_src] * m_uxpold[i_src] + + m_uypold[i_src] * m_uypold[i_src] + + m_uzpold[i_src] * m_uzpold[i_src])); + const amrex::Real t_new_p = m_gammaboost * m_t_boost - m_uzfrm * zpnew * m_inv_c2; + const amrex::Real z_new_p = m_gammaboost* ( zpnew + m_betaboost * m_Phys_c * m_t_boost); + const amrex::Real uz_new_p = m_gammaboost * m_uzpnew[i_src] - gamma_new_p * m_uzfrm; + const amrex::Real t_old_p = m_gammaboost * (m_t_boost - m_dt) + - m_uzfrm * m_zpold[i_src] * m_inv_c2; + const amrex::Real z_old_p = m_gammaboost * ( m_zpold[i_src] + m_betaboost + * m_Phys_c * (m_t_boost - m_dt ) ); + const amrex::Real uz_old_p = m_gammaboost * m_uzpold[i_src] - gamma_old_p * m_uzfrm; + // interpolate in time to t_lab + const amrex::Real weight_old = (t_new_p - m_t_lab) + / (t_new_p - t_old_p); + const amrex::Real weight_new = (m_t_lab - t_old_p) + / (t_new_p - t_old_p); + // weighted sum of old and new values + const amrex::ParticleReal xp = m_xpold[i_src] * weight_old + xpnew * weight_new; + const amrex::ParticleReal yp = m_ypold[i_src] * weight_old + ypnew * weight_new; + const amrex::ParticleReal zp = z_old_p * weight_old + z_new_p * weight_new; + const amrex::ParticleReal uxp = m_uxpold[i_src] * weight_old + + m_uxpnew[i_src] * weight_new; + const amrex::ParticleReal uyp = m_uypold[i_src] * weight_old + + m_uypnew[i_src] * weight_new; + const amrex::ParticleReal uzp = uz_old_p * weight_old + + uz_new_p * weight_new; +#if defined (WARPX_DIM_3D) + dst.m_aos[i_dst].pos(0) = xp; + dst.m_aos[i_dst].pos(1) = yp; + dst.m_aos[i_dst].pos(2) = zp; +#elif defined (WARPX_DIM_XZ) + dst.m_aos[i_dst].pos(0) = xp; + dst.m_aos[i_dst].pos(1) = zp; + amrex::ignore_unused(yp); +#elif defined (WARPX_DIM_1D) + dst.m_aos[i_dst].pos(0) = zp; + amrex::ignore_unused(xp, yp); +#else + amrex::ignore_unused(xp, yp, zp); +#endif + dst.m_rdata[PIdx::w][i_dst] = m_wpnew[i_src]; + dst.m_rdata[PIdx::ux][i_dst] = uxp; + dst.m_rdata[PIdx::uy][i_dst] = uyp; + dst.m_rdata[PIdx::uz][i_dst] = uzp; + } + + GetParticlePosition m_get_position; + + amrex::ParticleReal* AMREX_RESTRICT m_xpold = nullptr; + amrex::ParticleReal* AMREX_RESTRICT m_ypold = nullptr; + amrex::ParticleReal* AMREX_RESTRICT m_zpold = nullptr; + + amrex::ParticleReal* AMREX_RESTRICT m_uxpold = nullptr; + amrex::ParticleReal* AMREX_RESTRICT m_uypold = nullptr; + amrex::ParticleReal* AMREX_RESTRICT m_uzpold = nullptr; + + const amrex::ParticleReal* AMREX_RESTRICT m_uxpnew = nullptr; + const amrex::ParticleReal* AMREX_RESTRICT m_uypnew = nullptr; + const amrex::ParticleReal* AMREX_RESTRICT m_uzpnew = nullptr; + const amrex::ParticleReal* AMREX_RESTRICT m_wpnew = nullptr; + + amrex::Real m_gammaboost; + amrex::Real m_betaboost; + amrex::Real m_Phys_c; + amrex::Real m_uzfrm; + amrex::Real m_inv_c2; + amrex::Real m_t_boost; + amrex::Real m_dt; + amrex::Real m_t_lab; +}; + +/** + * \brief BackTransform functor to select particles and Lorentz Transform them + * and store in particle buffers + */ +class +BackTransformParticleFunctor final : public ComputeParticleDiagFunctor +{ +public: + BackTransformParticleFunctor(WarpXParticleContainer *pc_src, std::string species_name, int num_buffers); + /** Computes the Lorentz transform of source particles to obtain lab-frame data in pc_dst*/ + void operator () (ParticleContainer& pc_dst, int &TotalParticleCounter, int i_buffer) const override; + void InitData() override; + /** \brief Prepare data required to back-transform particle attribtutes for lab-frame snapshot, i_buffer + * + * \param[in] i_buffer index of the snapshot + * \param[in] z_slice_in_domain if the z-slice at current_z_boost is within the bounds of + * the boosted-frame and lab-frame domain. The particles are transformed + * only if this value is true. + * \param[in] current_z_boost z co-ordinate of the slice selected in boosted-frame. + * \param[in] t_lab current time in lab-frame for snapshot, i_buffer. + * \param[in] snapshot_full if the current snapshot, with index, i_buffer, is + * full (1) or not (0). If it is full, then Lorentz-transform is not performed + * by setting m_perform_backtransform to 0 for the corresponding ith snapshot. + */ + void PrepareFunctorData ( int i_buffer, bool z_slice_in_domain, amrex::Real old_z_boost, + amrex::Real current_z_boost, amrex::Real t_lab, + int snapshot_full) override; +private: + /** Source particle container */ + WarpXParticleContainer* m_pc_src = nullptr; + /** String containing species name of particles being transformed*/ + std::string m_species_name; + /** Number of buffers or snapshots*/ + int m_num_buffers; + /** Vector of current z co-ordinate in the boosted frame for each snapshot*/ + amrex::Vector<amrex::Real> m_current_z_boost; + /** Vector of previous z co-ordinate in the boosted frame for each snapshot*/ + amrex::Vector<amrex::Real> m_old_z_boost; + /** Vector of lab-frame time for each snapshot*/ + amrex::Vector<amrex::Real> m_t_lab; + /** Vector of 0s and 1s stored to check if back-transformation is to be performed for + * the ith snapshot. The value is set to 0 (false) or 1 (true) using the + * boolean ZSliceInDomain in PrepareFunctorData() + */ + amrex::Vector<int> m_perform_backtransform; +}; + + + + +#endif // diff --git a/Source/Diagnostics/ComputeDiagFunctors/BackTransformParticleFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/BackTransformParticleFunctor.cpp new file mode 100644 index 000000000..b8dac11ab --- /dev/null +++ b/Source/Diagnostics/ComputeDiagFunctors/BackTransformParticleFunctor.cpp @@ -0,0 +1,173 @@ +/* Copyright 2021 Revathi Jambunathan + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#include "BackTransformParticleFunctor.H" +#include "Particles/Pusher/GetAndSetPosition.H" +#include "Particles/WarpXParticleContainer.H" +#include "WarpX.H" +#include <AMReX.H> +#include <AMReX_Print.H> +#include <AMReX_BaseFwd.H> + +SelectParticles::SelectParticles (const WarpXParIter& a_pti, TmpParticles& tmp_particle_data, + amrex::Real current_z_boost, amrex::Real old_z_boost, + int a_offset) + : m_current_z_boost(current_z_boost), m_old_z_boost(old_z_boost) +{ + m_get_position = GetParticlePosition(a_pti, a_offset); + + const auto lev = a_pti.GetLevel(); + const auto index = a_pti.GetPairIndex(); + + zpold = tmp_particle_data[lev][index][TmpIdx::zold].dataPtr(); +} + + +LorentzTransformParticles::LorentzTransformParticles ( const WarpXParIter& a_pti, + TmpParticles& tmp_particle_data, + amrex::Real t_boost, amrex::Real dt, + amrex::Real t_lab, int a_offset) + : m_t_boost(t_boost), m_dt(dt), m_t_lab(t_lab) +{ + using namespace amrex::literals; + + if (tmp_particle_data.size() == 0) return; + m_get_position = GetParticlePosition(a_pti, a_offset); + + auto& attribs = a_pti.GetAttribs(); + m_wpnew = attribs[PIdx::w].dataPtr(); + m_uxpnew = attribs[PIdx::ux].dataPtr(); + m_uypnew = attribs[PIdx::uy].dataPtr(); + m_uzpnew = attribs[PIdx::uz].dataPtr(); + + const auto lev = a_pti.GetLevel(); + const auto index = a_pti.GetPairIndex(); + + m_xpold = tmp_particle_data[lev][index][TmpIdx::xold].dataPtr(); + m_ypold = tmp_particle_data[lev][index][TmpIdx::yold].dataPtr(); + m_zpold = tmp_particle_data[lev][index][TmpIdx::zold].dataPtr(); + m_uxpold = tmp_particle_data[lev][index][TmpIdx::uxold].dataPtr(); + m_uypold = tmp_particle_data[lev][index][TmpIdx::uyold].dataPtr(); + m_uzpold = tmp_particle_data[lev][index][TmpIdx::uzold].dataPtr(); + + m_betaboost = WarpX::beta_boost; + m_gammaboost = WarpX::gamma_boost; + m_Phys_c = PhysConst::c; + m_inv_c2 = 1._rt/(m_Phys_c * m_Phys_c); + m_uzfrm = -m_gammaboost*m_betaboost*m_Phys_c; +} + +/** + * \brief Functor to compute Lorentz Transform and store the selected particles in existing + * particle buffers + */ +BackTransformParticleFunctor::BackTransformParticleFunctor ( + WarpXParticleContainer *pc_src, + std::string species_name, + int num_buffers) + : m_pc_src(pc_src), m_species_name(species_name), m_num_buffers(num_buffers) +{ + InitData(); +} + + +void +BackTransformParticleFunctor::operator () (ParticleContainer& pc_dst, int &totalParticleCounter, int i_buffer) const +{ + if (m_perform_backtransform[i_buffer] == 0) return; + auto &warpx = WarpX::GetInstance(); + // get particle slice + const int nlevs = std::max(0, m_pc_src->finestLevel()+1); + auto tmp_particle_data = m_pc_src->getTmpParticleData(); + int total_particles_added = 0; + for (int lev = 0; lev < nlevs; ++lev) { + amrex::Real t_boost = warpx.gett_new(0); + amrex::Real dt = warpx.getdt(0); + + for (WarpXParIter pti(*m_pc_src, lev); pti.isValid(); ++pti) { + auto ptile_dst = pc_dst.DefineAndReturnParticleTile(lev, pti.index(), pti.LocalTileIndex() ); + } + + auto& particles = m_pc_src->GetParticles(lev); +#ifdef AMREX_USE_OMP +#pragma omp parallel +#endif + { + // Temporary arrays to store copy_flag and copy_index for particles + // that cross the z-slice + amrex::Gpu::DeviceVector<int> FlagForPartCopy; + amrex::Gpu::DeviceVector<int> IndexForPartCopy; + + for (WarpXParIter pti(*m_pc_src, lev); pti.isValid(); ++pti) { + + auto index = std::make_pair(pti.index(), pti.LocalTileIndex()); + + const auto GetParticleFilter = SelectParticles(pti, tmp_particle_data, + m_current_z_boost[i_buffer], + m_old_z_boost[i_buffer]); + const auto GetParticleLorentzTransform = LorentzTransformParticles( + pti, tmp_particle_data, + t_boost, dt, + m_t_lab[i_buffer]); + + long const np = pti.numParticles(); + + FlagForPartCopy.resize(np); + IndexForPartCopy.resize(np); + + int* const AMREX_RESTRICT Flag = FlagForPartCopy.dataPtr(); + int* const AMREX_RESTRICT IndexLocation = IndexForPartCopy.dataPtr(); + + const auto& ptile_src = particles.at(index); + auto src_data = ptile_src.getConstParticleTileData(); + // Flag particles that need to be copied if they cross the z-slice + // setting this to 1 for testing (temporarily) + amrex::ParallelFor(np, + [=] AMREX_GPU_DEVICE(int i) + { + Flag[i] = GetParticleFilter(src_data, i); + }); + + const int total_partdiag_size = amrex::Scan::ExclusiveSum(np,Flag,IndexLocation); + auto& ptile_dst = pc_dst.DefineAndReturnParticleTile(lev, pti.index(), pti.LocalTileIndex() ); + auto old_size = ptile_dst.numParticles(); + ptile_dst.resize(old_size + total_partdiag_size); + auto count = amrex::filterParticles(ptile_dst, ptile_src, GetParticleFilter, 0, old_size, np); + auto dst_data = ptile_dst.getParticleTileData(); + amrex::ParallelFor(np, + [=] AMREX_GPU_DEVICE(int i) + { + if (Flag[i] == 1) GetParticleLorentzTransform(dst_data, src_data, i, + old_size + IndexLocation[i]); + }); + total_particles_added += count; + } + } + } + totalParticleCounter = pc_dst.TotalNumberOfParticles(); +} + + +void +BackTransformParticleFunctor::InitData() +{ + m_current_z_boost.resize(m_num_buffers); + m_old_z_boost.resize(m_num_buffers); + m_t_lab.resize(m_num_buffers); + m_perform_backtransform.resize(m_num_buffers); +} + +void +BackTransformParticleFunctor::PrepareFunctorData ( int i_buffer, bool z_slice_in_domain, + amrex::Real old_z_boost, amrex::Real current_z_boost, + amrex::Real t_lab, int snapshot_full) +{ + m_old_z_boost.at(i_buffer) = old_z_boost; + m_current_z_boost.at(i_buffer) = current_z_boost; + m_t_lab.at(i_buffer) = t_lab; + m_perform_backtransform.at(i_buffer) = 0; + if (z_slice_in_domain == true and snapshot_full == 0) m_perform_backtransform.at(i_buffer) = 1; +} diff --git a/Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt b/Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt index 0439c9d67..5e085db4d 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt +++ b/Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt @@ -7,4 +7,5 @@ target_sources(WarpX PartPerCellFunctor.cpp PartPerGridFunctor.cpp BackTransformFunctor.cpp + BackTransformParticleFunctor.cpp ) diff --git a/Source/Diagnostics/ComputeDiagFunctors/ComputeParticleDiagFunctor.H b/Source/Diagnostics/ComputeDiagFunctors/ComputeParticleDiagFunctor.H new file mode 100644 index 000000000..d7875ed24 --- /dev/null +++ b/Source/Diagnostics/ComputeDiagFunctors/ComputeParticleDiagFunctor.H @@ -0,0 +1,59 @@ +/* Copyright 2021 Revathi Jambunathan + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef WARPX_COMPUTEPARTICLEDIAGFUNCTOR_H_ +#define WARPX_COMPUTEPARTICLEDIAGFUNCTOR_H_ + +#include "Particles/WarpXParticleContainer.H" +#include <AMReX.H> +#include <AMReX_AmrParticles.H> + +/** + * \brief Functor to compute a diagnostic and store the result in existing ParticleContainer. + */ +class +ComputeParticleDiagFunctor +{ +public: + using ParticleContainer = typename amrex::AmrParticleContainer<0, 0, PIdx::nattribs, 0, amrex::PinnedArenaAllocator>; + + ComputeParticleDiagFunctor( ) {} + /** Virtual Destructor to handle clean destruction of derived classes */ + virtual ~ComputeParticleDiagFunctor() = default; + + /** \brief Prepare data required to back-transform particle attribtutes for + * lab-frame snapshot, with index i_buffer. + * Note that this function has parameters that are specific to + * back-transformed diagnostics, that are unused for regular diagnostics. + * + * \param[in] i_buffer index of the snapshot + * \param[in] z_slice_in_domain if the z-slice at current_z_boost is within the bounds of + * the boosted-frame and lab-frame domain. The particles are transformed + * only if this value is true. + * \param[in] current_z_boost z co-ordinate of the slice selected in boosted-frame. + * \param[in] t_lab current time in lab-frame for snapshot, i_buffer. + * \param[in] snapshot_full if the current snapshot, with index, i_buffer, is + * full (1) or not (0). If it is full, then Lorentz-transform is not performed + * by setting m_perform_backtransform to 0 for the corresponding ith snapshot. + */ + virtual void PrepareFunctorData ( int i_buffer, bool ZSliceInDomain, + amrex::Real old_z_boost, + amrex::Real current_z_boost, amrex::Real t_lab, + int snapshot_full) + { + amrex::ignore_unused(i_buffer, + ZSliceInDomain, old_z_boost, + current_z_boost, t_lab, snapshot_full); + } + /** Compute particle attributes and store the result in pc_dst particle container. + * \param[out] pc_dst output particle container where the result is stored. + * \param[in] i_buffer snapshot index for which the particle buffer is processed + */ + virtual void operator () (ParticleContainer& pc_dst, int &totalParticlesInBuffer, int i_buffer = 0) const = 0; + virtual void InitData () {} +}; + +#endif // WARPX_COMPUTEPARTICLEDIAGFUNCTOR_H_ diff --git a/Source/Diagnostics/ComputeDiagFunctors/Make.package b/Source/Diagnostics/ComputeDiagFunctors/Make.package index 297eee576..7a07273e2 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/Make.package +++ b/Source/Diagnostics/ComputeDiagFunctors/Make.package @@ -5,5 +5,6 @@ CEXE_sources += DivBFunctor.cpp CEXE_sources += DivEFunctor.cpp CEXE_sources += RhoFunctor.cpp CEXE_sources += BackTransformFunctor.cpp +CEXE_sources += BackTransformParticleFunctor.cpp VPATH_LOCATIONS += $(WARPX_HOME)/Source/Diagnostics/ComputeDiagFunctors diff --git a/Source/Diagnostics/Diagnostics.H b/Source/Diagnostics/Diagnostics.H index c4435e28c..1301ea6f5 100644 --- a/Source/Diagnostics/Diagnostics.H +++ b/Source/Diagnostics/Diagnostics.H @@ -4,11 +4,15 @@ #include "ParticleDiag/ParticleDiag.H" #include "ComputeDiagFunctors/ComputeDiagFunctor_fwd.H" +#include "ComputeDiagFunctors/ComputeParticleDiagFunctor.H" #include "FlushFormats/FlushFormat_fwd.H" +#include "Particles/WarpXParticleContainer.H" +#include "Particles/PinnedMemoryParticleContainer.H" #include <AMReX_IntVect.H> #include <AMReX_REAL.H> #include <AMReX_Vector.H> +#include <AMReX_AmrParticles.H> #include <AMReX_BaseFwd.H> @@ -58,6 +62,8 @@ public: * \param[in] lev level on which the vector of unique_ptrs to field functors is initialized. */ virtual void InitializeFieldFunctors (int lev) = 0; + /** Initialize functors that store pointers to the species data requested by the user. */ + virtual void InitializeParticleFunctors () {} /** whether to compute and pack data in output buffers at this time step * \param[in] step current time step * \param[in] force_flush if true, return true for any step @@ -92,17 +98,27 @@ protected: * \param[in] i_buffer index of buffer for which the output MultiFab is defined. * \param[in] lev level on which the output MultiFab is defined */ - virtual void InitializeFieldBufferData (int i_buffer, int lev) = 0; + virtual void InitializeBufferData (int i_buffer, int lev) = 0; /** Initialize member variables and arrays specific to the diagnostics in the * derived classes.(FullDiagnostics, BTDiagnostics) */ virtual void DerivedInitData () {} /** This function initialized particle buffers (not implemented in diagnostics, yet) */ virtual void InitializeParticleBuffer () = 0; + /** This function prepares buffer data as required for fields and particles. For back-transformed diagnostics, + * this function prepares the z coordinate in the boosted-frame and lab-frame. + */ + virtual void PrepareBufferData () {} + /** This function updates buffer data and computes the number of buffers filled in the output multifab as well + * as identifies if the last buffer has been filled as needed to close the output files. + */ + virtual void UpdateBufferData () {} /** Prepare data (either fill-boundary or cell-centered data for back-transform diagnostics) to be processed for diagnostics. */ virtual void PrepareFieldDataForOutput () {} + /** The Particle Geometry, BoxArray, and RealBox are set for the lab-frame output */ + virtual void PrepareParticleDataForOutput () {} /** Update the physical extent of the diagnostic domain for moving window and * galilean shift simulations * @@ -154,8 +170,10 @@ protected: std::vector< std::string > m_output_species_names; /** Names of all species in the simulation */ std::vector< std::string > m_all_species_names; - /** Each element of this vector handles output for 1 species */ - amrex::Vector< ParticleDiag > m_output_species; + /** The first vector is for total number of snapshots. (=1 for FullDiagnostics) + The second vector handles output for 1 species + */ + amrex::Vector< amrex::Vector< ParticleDiag> > m_output_species; /** Vector of (pointers to) functors to compute output fields, per level, * per component. This allows for simple operations (averaging to * cell-center for standard EB fields) as well as more involved operations @@ -173,8 +191,22 @@ protected: int m_num_buffers; /** Array of species indices that dump rho per species */ amrex::Vector<int> m_rho_per_species_index; - /** Temporarily clear species output for BTD until particle buffer is added */ - virtual void TMP_ClearSpeciesDataForBTD() {} + /** Vector of particle buffer vectors for each snapshot */ + amrex::Vector< amrex::Vector<std::unique_ptr<PinnedMemoryParticleContainer> > > m_particles_buffer; + /** Vector of pointers to functors to compute particle output per species*/ + amrex::Vector< std::unique_ptr<ComputeParticleDiagFunctor> > m_all_particle_functors; + + /** Vector of total number of particles previously flushed, per species, per snapshot. + * The first vector is for total number of snapshots and second vector loops + * over the total number of species selected for diagnostics. + */ + amrex::Vector< amrex::Vector <int> > m_totalParticles_flushed_already; + /** Vector of total number of particles in the buffer, per species, per snapshot. + * The first vector is for total number of snapshots and second vector loops + * over the total number of species selected for diagnostics. + */ + amrex::Vector< amrex::Vector <int> > m_totalParticles_in_buffer; + }; #endif // WARPX_DIAGNOSTICS_H_ diff --git a/Source/Diagnostics/Diagnostics.cpp b/Source/Diagnostics/Diagnostics.cpp index 903219c31..7792a7c3d 100644 --- a/Source/Diagnostics/Diagnostics.cpp +++ b/Source/Diagnostics/Diagnostics.cpp @@ -1,6 +1,7 @@ #include "Diagnostics.H" #include "Diagnostics/ComputeDiagFunctors/ComputeDiagFunctor.H" +#include "ComputeDiagFunctors/BackTransformParticleFunctor.H" #include "Diagnostics/FlushFormats/FlushFormat.H" #include "Diagnostics/ParticleDiag/ParticleDiag.H" #include "FlushFormats/FlushFormatAscent.H" @@ -197,33 +198,42 @@ Diagnostics::InitData () for (int lev = 0; lev < nmax_lev; ++lev) { // allocate and initialize m_all_field_functors depending on diag type InitializeFieldFunctors(lev); - // Initialize field buffer data, m_mf_output - InitializeFieldBufferData(i_buffer, lev); + // Initialize buffer data required for particle and/or fields + InitializeBufferData(i_buffer, lev); } } - // When particle buffers, m_particle_boundary_buffer are included, - // they will be initialized here - InitializeParticleBuffer(); amrex::ParmParse pp_diag_name(m_diag_name); + // default for writing species output is 1 + int write_species = 1; + pp_diag_name.query("write_species", write_species); + if (write_species == 1) { + // When particle buffers, m_particle_boundary_buffer are included, + // they will be initialized here + InitializeParticleBuffer(); + InitializeParticleFunctors(); + } + amrex::Vector <amrex::Real> dummy_val(AMREX_SPACEDIM); if ( queryArrWithParser(pp_diag_name, "diag_lo", dummy_val, 0, AMREX_SPACEDIM) || queryArrWithParser(pp_diag_name, "diag_hi", dummy_val, 0, AMREX_SPACEDIM) ) { // set geometry filter for particle-diags to true when the diagnostic domain-extent - // is specified by the user - for (int i = 0; i < m_output_species.size(); ++i) { - m_output_species[i].m_do_geom_filter = true; + // is specified by the user. + // Note that the filter is set for every ith snapshot, and the number of snapshots + // for full diagnostics is 1, while for BTD it is user-defined. + for (int i_buffer = 0; i_buffer < m_num_buffers; ++i_buffer ) { + for (int i = 0; i < m_output_species.size(); ++i) { + m_output_species[i_buffer][i].m_do_geom_filter = true; + } + // Disabling particle-io for reduced domain diagnostics by reducing + // the particle-diag vector to zero. + // This is a temporary fix until particle_buffer is supported in diagnostics. + m_output_species[i_buffer].clear(); } - // Disabling particle-io for reduced domain diagnostics by reducing - // the particle-diag vector to zero. - // This is a temporary fix until particle_buffer is supported in diagnostics. m_output_species.clear(); amrex::Print() << " WARNING: For full diagnostics on a reduced domain, particle io is not supported, yet! Therefore, particle-io is disabled for this diag " << m_diag_name << "\n"; } - // default for writing species output is 1 - int write_species = 1; - pp_diag_name.query("write_species", write_species); if (write_species == 0) { if (m_format == "checkpoint"){ amrex::Abort("For checkpoint format, write_species flag must be 1."); @@ -232,8 +242,6 @@ Diagnostics::InitData () m_output_species.clear(); m_output_species_names.clear(); } - // temporarily clear out species output sincce particle buffers are not supported. - TMP_ClearSpeciesDataForBTD(); } @@ -294,13 +302,19 @@ Diagnostics::InitBaseData () for (int i = 0; i < m_num_buffers; ++i) { m_geom_output[i].resize( nmax_lev ); } + } void Diagnostics::ComputeAndPack () { + PrepareBufferData(); // prepare the field-data necessary to compute output data PrepareFieldDataForOutput(); + // Prepare the particle data necessary to compute output data + // Field-data is called first for BTD, since the z-slice location is used to prepare particle data + // to determine if the transform is to be done this step. + PrepareParticleDataForOutput(); auto & warpx = WarpX::GetInstance(); @@ -308,7 +322,7 @@ Diagnostics::ComputeAndPack () for (int i_buffer = 0; i_buffer < m_num_buffers; ++i_buffer) { for(int lev=0; lev<nlev_output; lev++){ int icomp_dst = 0; - for (int icomp=0, n=m_all_field_functors[0].size(); icomp<n; icomp++){ + for (int icomp=0, n=m_all_field_functors[lev].size(); icomp<n; icomp++){ // Call all functors in m_all_field_functors[lev]. Each of them computes // a diagnostics and writes in one or more components of the output // multifab m_mf_output[lev]. @@ -324,7 +338,13 @@ Diagnostics::ComputeAndPack () WarpXCommUtil::FillBoundary(m_mf_output[i_buffer][lev], warpx.Geom(lev).periodicity()); } } + // Call Particle functor + for (int isp = 0; isp < m_all_particle_functors.size(); ++isp) { + m_all_particle_functors[isp]->operator()(*m_particles_buffer[i_buffer][isp], m_totalParticles_in_buffer[i_buffer][isp], i_buffer); + } } + + UpdateBufferData(); } @@ -332,17 +352,16 @@ void Diagnostics::FilterComputePackFlush (int step, bool force_flush) { WARPX_PROFILE("Diagnostics::FilterComputePackFlush()"); - MovingWindowAndGalileanDomainShift (step); if ( DoComputeAndPack (step, force_flush) ) { ComputeAndPack(); + } - for (int i_buffer = 0; i_buffer < m_num_buffers; ++i_buffer) { - if ( !DoDump (step, i_buffer, force_flush) ) continue; - Flush(i_buffer); - } - + for (int i_buffer = 0; i_buffer < m_num_buffers; ++i_buffer) { + if ( !DoDump (step, i_buffer, force_flush) ) continue; + Flush(i_buffer); } + } diff --git a/Source/Diagnostics/FlushFormats/FlushFormat.H b/Source/Diagnostics/FlushFormats/FlushFormat.H index f7b28cf58..4b23796bf 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormat.H +++ b/Source/Diagnostics/FlushFormats/FlushFormat.H @@ -22,7 +22,8 @@ public: bool plot_raw_fields_guards, bool isBTD = false, int snapshotID = -1, const amrex::Geometry& full_BTD_snapshot = amrex::Geometry(), - bool isLastBTDFlush = false) const = 0; + bool isLastBTDFlush = false, + const amrex::Vector<int>& totalParticlesFlushedAlready = amrex::Vector<int>() ) const = 0; virtual ~FlushFormat() {} }; diff --git a/Source/Diagnostics/FlushFormats/FlushFormatAscent.H b/Source/Diagnostics/FlushFormats/FlushFormatAscent.H index 18f0a5e10..ca2339469 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatAscent.H +++ b/Source/Diagnostics/FlushFormats/FlushFormatAscent.H @@ -39,7 +39,8 @@ public: bool plot_raw_fields_guards, bool isBTD = false, int snapshotID = -1, const amrex::Geometry& full_BTD_snapshot = amrex::Geometry(), - bool isLastBTDFlush = false) const override; + bool isLastBTDFlush = false, + const amrex::Vector<int>& totalParticlesFlushedAlready = amrex::Vector<int>() ) const override; /** \brief Do in-situ visualization for particle data. * \param[in] particle_diags Each element of this vector handles output of 1 species. diff --git a/Source/Diagnostics/FlushFormats/FlushFormatAscent.cpp b/Source/Diagnostics/FlushFormats/FlushFormatAscent.cpp index a3f8d593e..a6811f5f7 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatAscent.cpp +++ b/Source/Diagnostics/FlushFormats/FlushFormatAscent.cpp @@ -17,7 +17,8 @@ FlushFormatAscent::WriteToFile ( const amrex::Vector<ParticleDiag>& particle_diags, int nlev, const std::string prefix, int /*file_min_digits*/, bool plot_raw_fields, bool plot_raw_fields_guards, - bool /*isBTD*/, int /*snapshotID*/, const amrex::Geometry& /*full_BTD_snapshot*/, bool /*isLastBTDFlush*/) const + bool /*isBTD*/, int /*snapshotID*/, const amrex::Geometry& /*full_BTD_snapshot*/, + bool /*isLastBTDFlush*/, const amrex::Vector<int>& /* totalParticlesFlushedAlready*/) const { #ifdef AMREX_USE_ASCENT WARPX_PROFILE("FlushFormatAscent::WriteToFile()"); diff --git a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.H b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.H index e057e8854..b91c51862 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.H +++ b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.H @@ -26,7 +26,8 @@ class FlushFormatCheckpoint final : public FlushFormatPlotfile bool plot_raw_fields_guards, bool isBTD = false, int snapshotID = -1, const amrex::Geometry& full_BTD_snapshot = amrex::Geometry(), - bool isLastBTDFlush = false) const override final; + bool isLastBTDFlush = false, + const amrex::Vector<int>& totalParticlesFlushedAlready = amrex::Vector<int>() ) const override final; void CheckpointParticles (const std::string& dir, const amrex::Vector<ParticleDiag>& particle_diags) const; diff --git a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp index 915839683..890d4ceec 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp +++ b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp @@ -33,7 +33,7 @@ FlushFormatCheckpoint::WriteToFile ( bool /*plot_raw_fields_guards*/, bool /*isBTD*/, int /*snapshotID*/, const amrex::Geometry& /*full_BTD_snapshot*/, - bool /*isLastBTDFlush*/) const + bool /*isLastBTDFlush*/, const amrex::Vector<int>& /* totalParticlesFlushedAlready*/) const { WARPX_PROFILE("FlushFormatCheckpoint::WriteToFile()"); diff --git a/Source/Diagnostics/FlushFormats/FlushFormatOpenPMD.H b/Source/Diagnostics/FlushFormats/FlushFormatOpenPMD.H index bdac2147b..3979e7b80 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatOpenPMD.H +++ b/Source/Diagnostics/FlushFormats/FlushFormatOpenPMD.H @@ -38,7 +38,8 @@ public: bool plot_raw_fields_guards, bool isBTD = false, int snapshotID = -1, const amrex::Geometry& full_BTD_snapshot = amrex::Geometry(), - bool isLastBTDFlush = false) const override; + bool isLastBTDFlush = false, + const amrex::Vector<int>& totalParticlesFlushedAlready = amrex::Vector<int>() ) const override; ~FlushFormatOpenPMD () override = default; diff --git a/Source/Diagnostics/FlushFormats/FlushFormatOpenPMD.cpp b/Source/Diagnostics/FlushFormats/FlushFormatOpenPMD.cpp index 1a611680a..401d0e184 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatOpenPMD.cpp +++ b/Source/Diagnostics/FlushFormats/FlushFormatOpenPMD.cpp @@ -82,6 +82,26 @@ FlushFormatOpenPMD::FlushFormatOpenPMD (const std::string& diag_name) operator_type, operator_parameters, warpx.getPMLdirections() ); + + // Temporarily adding Abort for adios filetype if species is selected for BTD output + bool species_output = true; + int write_species = 1; + std::vector< std::string > output_species_names; + bool species_specified = pp_diag_name.queryarr("species", output_species_names); + if (species_specified == true and output_species_names.size() > 0) { + species_output = true; + } else { + // By default species output is computed for all diagnostics, if write_species is not set to 0 + species_output = true; + } + // Check user-defined option to turn off species output + pp_diag_name.query("write_species", write_species); + if (write_species == 0) species_output = false; + if (diag_type_str == "BackTransformed" and species_output == true) { + if (m_OpenPMDPlotWriter->OpenPMDFileType() == "bp") { + amrex::Abort(" Currently BackTransformed diagnostics type does not support species output for ADIOS backend. Please select h5 as openpmd backend"); + } + } } void @@ -94,7 +114,7 @@ FlushFormatOpenPMD::WriteToFile ( const std::string prefix, int file_min_digits, bool plot_raw_fields, bool plot_raw_fields_guards, bool isBTD, int snapshotID, const amrex::Geometry& full_BTD_snapshot, - bool isLastBTDFlush) const + bool isLastBTDFlush, const amrex::Vector<int>& totalParticlesFlushedAlready) const { WARPX_PROFILE("FlushFormatOpenPMD::WriteToFile()"); @@ -117,7 +137,7 @@ FlushFormatOpenPMD::WriteToFile ( varnames, mf, geom, output_iteration, time, isBTD, full_BTD_snapshot); // particles: all (reside only on locally finest level) - m_OpenPMDPlotWriter->WriteOpenPMDParticles(particle_diags, isBTD); + m_OpenPMDPlotWriter->WriteOpenPMDParticles(particle_diags, isBTD, totalParticlesFlushedAlready); // signal that no further updates will be written to this iteration m_OpenPMDPlotWriter->CloseStep(isBTD, isLastBTDFlush); diff --git a/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.H b/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.H index 37bd3c2e0..0d673607b 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.H +++ b/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.H @@ -33,7 +33,8 @@ public: bool plot_raw_fields_guards, bool isBTD = false, int snapshotID = -1, const amrex::Geometry& full_BTD_snapshot = amrex::Geometry(), - bool isLastBTDFlush = false) const override; + bool isLastBTDFlush = false, + const amrex::Vector<int>& totalParticlesFlushedAlready = amrex::Vector<int>() ) const override; /** Write general info of the run into the plotfile */ void WriteJobInfo(const std::string& dir) const; @@ -47,7 +48,8 @@ public: * \param[in] particle_diags Each element of this vector handles output of 1 species. */ void WriteParticles(const std::string& filename, - const amrex::Vector<ParticleDiag>& particle_diags) const; + const amrex::Vector<ParticleDiag>& particle_diags, + bool isBTD = false) const; ~FlushFormatPlotfile() {} }; diff --git a/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp b/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp index b1d2a8723..4e9fe41cb 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp +++ b/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp @@ -4,6 +4,7 @@ #include "Particles/Filter/FilterFunctors.H" #include "Particles/WarpXParticleContainer.H" #include "Particles/ParticleBuffer.H" +#include "Particles/PinnedMemoryParticleContainer.H" #include "Utils/Interpolate.H" #include "Utils/WarpXProfilerWrapper.H" #include "WarpX.H" @@ -58,8 +59,8 @@ FlushFormatPlotfile::WriteToFile ( const amrex::Vector<ParticleDiag>& particle_diags, int nlev, const std::string prefix, int file_min_digits, bool plot_raw_fields, bool plot_raw_fields_guards, - bool /*isBTD*/, int /*snapshotID*/, const amrex::Geometry& /*full_BTD_snapshot*/, - bool /*isLastBTDFlush*/) const + bool isBTD, int /*snapshotID*/, const amrex::Geometry& /*full_BTD_snapshot*/, + bool /*isLastBTDFlush*/, const amrex::Vector<int>& /* totalParticlesFlushedAlready*/) const { WARPX_PROFILE("FlushFormatPlotfile::WriteToFile()"); auto & warpx = WarpX::GetInstance(); @@ -82,7 +83,7 @@ FlushFormatPlotfile::WriteToFile ( WriteAllRawFields(plot_raw_fields, nlev, filename, plot_raw_fields_guards); - WriteParticles(filename, particle_diags); + WriteParticles(filename, particle_diags, isBTD); WriteJobInfo(filename); @@ -294,13 +295,20 @@ FlushFormatPlotfile::WriteWarpXHeader( } void -FlushFormatPlotfile::WriteParticles (const std::string& dir, - const amrex::Vector<ParticleDiag>& particle_diags) const +FlushFormatPlotfile::WriteParticles(const std::string& dir, + const amrex::Vector<ParticleDiag>& particle_diags, + bool isBTD) const { for (unsigned i = 0, n = particle_diags.size(); i < n; ++i) { WarpXParticleContainer* pc = particle_diags[i].getParticleContainer(); auto tmp = ParticleBuffer::getTmpPC<amrex::PinnedArenaAllocator>(pc); + if (isBTD) { + PinnedMemoryParticleContainer* pinned_pc = particle_diags[i].getPinnedParticleContainer(); + tmp.SetParticleGeometry(0,pinned_pc->Geom(0)); + tmp.SetParticleBoxArray(0,pinned_pc->ParticleBoxArray(0)); + tmp.SetParticleDistributionMap(0, pinned_pc->ParticleDistributionMap(0)); + } Vector<std::string> real_names; Vector<std::string> int_names; @@ -348,15 +356,19 @@ FlushFormatPlotfile::WriteParticles (const std::string& dir, GeometryFilter const geometry_filter(particle_diags[i].m_do_geom_filter, particle_diags[i].m_diag_domain); - using SrcData = WarpXParticleContainer::ParticleTileType::ConstParticleTileDataType; - tmp.copyParticles(*pc, - [=] AMREX_GPU_HOST_DEVICE (const SrcData& src, int ip, const amrex::RandomEngine& engine) - { - const SuperParticleType& p = src.getSuperParticle(ip); - return random_filter(p, engine) * uniform_filter(p, engine) - * parser_filter(p, engine) * geometry_filter(p, engine); - }, true); - + if (!isBTD) { + using SrcData = WarpXParticleContainer::ParticleTileType::ConstParticleTileDataType; + tmp.copyParticles(*pc, + [=] AMREX_GPU_HOST_DEVICE (const SrcData& src, int ip, const amrex::RandomEngine& engine) + { + const SuperParticleType& p = src.getSuperParticle(ip); + return random_filter(p, engine) * uniform_filter(p, engine) + * parser_filter(p, engine) * geometry_filter(p, engine); + }, true); + } else { + PinnedMemoryParticleContainer* pinned_pc = particle_diags[i].getPinnedParticleContainer(); + tmp.copyParticles(*pinned_pc); + } // real_names contains a list of all particle attributes. // real_flags & int_flags are 1 or 0, whether quantity is dumped or not. tmp.WritePlotFile( diff --git a/Source/Diagnostics/FlushFormats/FlushFormatSensei.H b/Source/Diagnostics/FlushFormats/FlushFormatSensei.H index 513b16b61..7f4371d3b 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatSensei.H +++ b/Source/Diagnostics/FlushFormats/FlushFormatSensei.H @@ -57,7 +57,8 @@ public: bool plot_raw_fields_guards, bool isBTD = false, int snapshotID = -1, const amrex::Geometry& full_BTD_snapshot = amrex::Geometry(), - bool isLastBTDFlush = false) const override; + bool isLastBTDFlush = false, + const amrex::Vector<int>& totalParticlesFlushedAlready = amrex::Vector<int>() ) const override; /** \brief Do in-situ visualization for particle data. * \param[in] particle_diags Each element of this vector handles output of 1 species. diff --git a/Source/Diagnostics/FlushFormats/FlushFormatSensei.cpp b/Source/Diagnostics/FlushFormats/FlushFormatSensei.cpp index 9ccb1bcaa..443d39e97 100644 --- a/Source/Diagnostics/FlushFormats/FlushFormatSensei.cpp +++ b/Source/Diagnostics/FlushFormats/FlushFormatSensei.cpp @@ -57,13 +57,14 @@ FlushFormatSensei::WriteToFile ( int nlev, const std::string prefix, int file_min_digits, bool plot_raw_fields, bool plot_raw_fields_guards, bool isBTD, int snapshotID, - const amrex::Geometry& full_BTD_snapshot, bool isLastBTDFlush) const + const amrex::Geometry& full_BTD_snapshot, bool isLastBTDFlush, + const amrex::Vector<int>& totalParticlesFlushedAlready) const { amrex::ignore_unused( geom, nlev, prefix, file_min_digits, plot_raw_fields, plot_raw_fields_guards, isBTD, snapshotID, full_BTD_snapshot, - isLastBTDFlush); + isLastBTDFlush, totalParticlesFlushedAlready); #ifndef AMREX_USE_SENSEI_INSITU amrex::ignore_unused(varnames, mf, iteration, time, particle_diags); diff --git a/Source/Diagnostics/FullDiagnostics.H b/Source/Diagnostics/FullDiagnostics.H index 85904a3cf..0bc41051f 100644 --- a/Source/Diagnostics/FullDiagnostics.H +++ b/Source/Diagnostics/FullDiagnostics.H @@ -57,7 +57,7 @@ private: * \param[in] i_buffer index of a back-transformed snapshot * \param[in] lev level on which source multifabs are defined */ - void InitializeFieldBufferData ( int i_buffer, int lev ) override; + void InitializeBufferData ( int i_buffer, int lev ) override; /** Initialize functors that store pointers to the fields requested by the user. * \param[in] lev level on which the vector of unique_ptrs to field functors is initialized. */ @@ -65,6 +65,8 @@ private: void InitializeParticleBuffer () override; /** Prepare field data to be used for diagnostics */ void PrepareFieldDataForOutput () override; + /** Prepare particle data to be used for diagnostics. */ + void PrepareParticleDataForOutput() override {} /** Update the physical extent of the diagnostic domain for moving window and * galilean shift simulations * diff --git a/Source/Diagnostics/FullDiagnostics.cpp b/Source/Diagnostics/FullDiagnostics.cpp index 6b057f574..81f04a330 100644 --- a/Source/Diagnostics/FullDiagnostics.cpp +++ b/Source/Diagnostics/FullDiagnostics.cpp @@ -62,9 +62,12 @@ FullDiagnostics::InitializeParticleBuffer () } } // Initialize one ParticleDiag per species requested - for (auto const& species : m_output_species_names){ - const int idx = mpc.getSpeciesID(species); - m_output_species.push_back(ParticleDiag(m_diag_name, species, mpc.GetParticleContainerPtr(idx))); + m_output_species.resize(m_num_buffers); + for (int i_buffer = 0; i_buffer < m_num_buffers; ++i_buffer) { + for (auto const& species : m_output_species_names){ + const int idx = mpc.getSpeciesID(species); + m_output_species[i_buffer].push_back(ParticleDiag(m_diag_name, species, mpc.GetParticleContainerPtr(idx))); + } } } @@ -123,8 +126,8 @@ FullDiagnostics::Flush ( int i_buffer ) m_flush_format->WriteToFile( m_varnames, m_mf_output[i_buffer], m_geom_output[i_buffer], warpx.getistep(), - warpx.gett_new(0), m_output_species, nlev_output, m_file_prefix, m_file_min_digits, - m_plot_raw_fields, m_plot_raw_fields_guards); + warpx.gett_new(0), m_output_species[i_buffer], nlev_output, m_file_prefix, + m_file_min_digits, m_plot_raw_fields, m_plot_raw_fields_guards); FlushRaw(); } @@ -274,7 +277,7 @@ FullDiagnostics::AddRZModesToOutputNames (const std::string& field, int ncomp){ void -FullDiagnostics::InitializeFieldBufferData (int i_buffer, int lev ) { +FullDiagnostics::InitializeBufferData (int i_buffer, int lev ) { auto & warpx = WarpX::GetInstance(); amrex::RealBox diag_dom; bool use_warpxba = true; @@ -468,8 +471,15 @@ FullDiagnostics::PrepareFieldDataForOutput () warpx.FillBoundaryAux(warpx.getngUpdateAux()); // Update the RealBox used for the geometry filter in particle diags - for (int i = 0; i < m_output_species.size(); ++i) { - m_output_species[i].m_diag_domain = m_geom_output[0][0].ProbDomain(); + // Note that full diagnostics every diag has only one buffer. (m_num_buffers = 1). + // For m_geom_output[i_buffer][lev], the first element is the buffer index, and + // second is level=0 + // The level is set to 0, because the whole physical domain of the simulation is used + // to set the domain dimensions for the output particle container. + for (int i_buffer = 0; i_buffer < m_num_buffers; ++i_buffer) { + for (int i = 0; i < m_output_species[i_buffer].size(); ++i) { + m_output_species[i_buffer][i].m_diag_domain = m_geom_output[i_buffer][0].ProbDomain(); + } } } @@ -481,6 +491,10 @@ FullDiagnostics::MovingWindowAndGalileanDomainShift (int step) // Account for galilean shift amrex::Real new_lo[AMREX_SPACEDIM]; amrex::Real new_hi[AMREX_SPACEDIM]; + // Note that Full diagnostics has only one snapshot, m_num_buffers = 1 + // m_geom_output[i_buffer][lev] below have values 0 and 0, respectively, because + // we need the physical extent from mesh-refinement level = 0, + // and only for the 0th snapshot, since full diagnostics has only one snapshot. const amrex::Real* current_lo = m_geom_output[0][0].ProbLo(); const amrex::Real* current_hi = m_geom_output[0][0].ProbHi(); @@ -502,11 +516,12 @@ FullDiagnostics::MovingWindowAndGalileanDomainShift (int step) new_hi[0] = current_hi[0] + warpx.m_galilean_shift[2]; } #endif - // Update RealBox of geometry with galilean-shifted boundary + // Update RealBox of geometry with galilean-shifted boundary. for (int lev = 0; lev < nmax_lev; ++lev) { + // Note that Full diagnostics has only one snapshot, m_num_buffers = 1 + // Thus here we set the prob domain for the 0th snapshot only. m_geom_output[0][lev].ProbDomain( amrex::RealBox(new_lo, new_hi) ); } - // For Moving Window Shift if (warpx.moving_window_active(step+1)) { int moving_dir = warpx.moving_window_dir; @@ -527,6 +542,8 @@ FullDiagnostics::MovingWindowAndGalileanDomainShift (int step) new_hi[moving_dir] = cur_hi[moving_dir] + num_shift_base*geom_dx[moving_dir]; // Update RealBox of geometry with shifted domain geometry for moving-window for (int lev = 0; lev < nmax_lev; ++lev) { + // Note that Full diagnostics has only one snapshot, m_num_buffers = 1 + // Thus here we set the prob domain for the 0th snapshot only. m_geom_output[0][lev].ProbDomain( amrex::RealBox(new_lo, new_hi) ); } } diff --git a/Source/Diagnostics/MultiDiagnostics.H b/Source/Diagnostics/MultiDiagnostics.H index e4f6dcf0d..d5ccf2113 100644 --- a/Source/Diagnostics/MultiDiagnostics.H +++ b/Source/Diagnostics/MultiDiagnostics.H @@ -26,7 +26,7 @@ public: /** \brief Loop over diags in alldiags and call their InitDiags */ void InitData (); /** \brief Called at each iteration. Compute diags and flush. */ - void FilterComputePackFlush (int step, bool force_flush=false); + void FilterComputePackFlush (int step, bool force_flush=false, bool BackTransform=false); /** \brief Called only at the last iteration. Loop over each diag and if m_dump_last_timestep * is true, compute diags and flush with force_flush=true. */ void FilterComputePackFlushLastTimestep (int step); diff --git a/Source/Diagnostics/MultiDiagnostics.cpp b/Source/Diagnostics/MultiDiagnostics.cpp index 680bd462a..ad6bd500a 100644 --- a/Source/Diagnostics/MultiDiagnostics.cpp +++ b/Source/Diagnostics/MultiDiagnostics.cpp @@ -69,10 +69,18 @@ MultiDiagnostics::ReadParameters () } void -MultiDiagnostics::FilterComputePackFlush (int step, bool force_flush) +MultiDiagnostics::FilterComputePackFlush (int step, bool force_flush, bool BackTransform) { + int i = 0; for (auto& diag : alldiags){ - diag->FilterComputePackFlush (step, force_flush); + if (BackTransform == true) { + if (diags_types[i] == DiagTypes::BackTransformed) + diag->FilterComputePackFlush (step, force_flush); + } else { + if (diags_types[i] != DiagTypes::BackTransformed) + diag->FilterComputePackFlush (step, force_flush); + } + ++i; } } diff --git a/Source/Diagnostics/ParticleDiag/ParticleDiag.H b/Source/Diagnostics/ParticleDiag/ParticleDiag.H index a54377ce5..c625e900d 100644 --- a/Source/Diagnostics/ParticleDiag/ParticleDiag.H +++ b/Source/Diagnostics/ParticleDiag/ParticleDiag.H @@ -4,6 +4,7 @@ #include "ParticleDiag_fwd.H" #include "Particles/WarpXParticleContainer_fwd.H" +#include "Particles/PinnedMemoryParticleContainer.H" #include <AMReX_Parser.H> #include <AMReX_REAL.H> @@ -16,8 +17,9 @@ class ParticleDiag { public: - ParticleDiag(std::string diag_name, std::string name, WarpXParticleContainer* pc); + ParticleDiag(std::string diag_name, std::string name, WarpXParticleContainer* pc, PinnedMemoryParticleContainer *pinned_pc = nullptr); WarpXParticleContainer* getParticleContainer() const { return m_pc; } + PinnedMemoryParticleContainer* getPinnedParticleContainer() const { return m_pinned_pc; } std::string getSpeciesName() const { return m_name; } amrex::Vector<int> plot_flags; @@ -36,6 +38,7 @@ private: std::string m_name; amrex::Vector< std::string > variables; WarpXParticleContainer* m_pc; + PinnedMemoryParticleContainer* m_pinned_pc; }; #endif // WARPX_PARTICLEDIAG_H_ diff --git a/Source/Diagnostics/ParticleDiag/ParticleDiag.cpp b/Source/Diagnostics/ParticleDiag/ParticleDiag.cpp index eb00b6c09..80fc3aabd 100644 --- a/Source/Diagnostics/ParticleDiag/ParticleDiag.cpp +++ b/Source/Diagnostics/ParticleDiag/ParticleDiag.cpp @@ -12,8 +12,8 @@ using namespace amrex; -ParticleDiag::ParticleDiag(std::string diag_name, std::string name, WarpXParticleContainer* pc) - : m_diag_name(diag_name), m_name(name), m_pc(pc) +ParticleDiag::ParticleDiag(std::string diag_name, std::string name, WarpXParticleContainer* pc, PinnedMemoryParticleContainer* pinned_pc) + : m_diag_name(diag_name), m_name(name), m_pc(pc), m_pinned_pc(pinned_pc) { ParmParse pp_diag_name_species_name(diag_name + "." + name); if (!pp_diag_name_species_name.queryarr("variables", variables)){ diff --git a/Source/Diagnostics/WarpXOpenPMD.H b/Source/Diagnostics/WarpXOpenPMD.H index 925e2a6f3..21c52ef10 100644 --- a/Source/Diagnostics/WarpXOpenPMD.H +++ b/Source/Diagnostics/WarpXOpenPMD.H @@ -129,8 +129,10 @@ public: */ void CloseStep (bool isBTD = false, bool isLastBTDFlush = false); - void WriteOpenPMDParticles (const amrex::Vector<ParticleDiag>& particle_diags, - const bool isBTD = false); + void WriteOpenPMDParticles ( + const amrex::Vector<ParticleDiag>& particle_diags, + const bool isBTD = false, + const amrex::Vector<int>& totalParticlesFlushedAlready = amrex::Vector<int>()); void WriteOpenPMDFieldsAll ( const std::vector<std::string>& varnames, @@ -140,6 +142,8 @@ public: bool isBTD = false, const amrex::Geometry& full_BTD_snapshot=amrex::Geometry() ) const; + /** Return OpenPMD File type ("bp" or "h5" or "json")*/ + std::string OpenPMDFileType () { return m_OpenPMDFileType; } private: void Init (openPMD::Access access, bool isBTD); @@ -199,7 +203,7 @@ private: openPMD::ParticleSpecies& currSpecies, const unsigned long long& np, amrex::ParticleReal const charge, - amrex::ParticleReal const mass) const ; + amrex::ParticleReal const mass, bool const isBTD = false); /** This function sets up the entries for particle properties * @@ -215,7 +219,7 @@ private: const amrex::Vector<std::string>& real_comp_names, const amrex::Vector<int>& write_int_comp, const amrex::Vector<std::string>& int_comp_names, - unsigned long long np) const; + const unsigned long long np, bool const isBTD = false) const; /** This function saves the values of the entries for particle properties * @@ -256,8 +260,8 @@ private: const amrex::Vector<std::string>& real_comp_names, const amrex::Vector<std::string>& int_comp_names, amrex::ParticleReal const charge, - amrex::ParticleReal const mass, - const bool isBTD) const; + amrex::ParticleReal const mass, const bool isBTD = false, + int ParticleFlushOffset = 0); /** Get the openPMD-api filename for openPMD::Series * @@ -271,6 +275,9 @@ private: std::unique_ptr<openPMD::Series> m_Series; + /** To set particle meta-data for the openpmd dump. */ + bool m_doParticleSetUp = true; + /** This is the output directory * * This usually does not yet end in a `/`. diff --git a/Source/Diagnostics/WarpXOpenPMD.cpp b/Source/Diagnostics/WarpXOpenPMD.cpp index 6f5c6a36b..daf5fbefc 100644 --- a/Source/Diagnostics/WarpXOpenPMD.cpp +++ b/Source/Diagnostics/WarpXOpenPMD.cpp @@ -489,7 +489,7 @@ WarpXOpenPMDPlot::Init (openPMD::Access access, bool isBTD) void WarpXOpenPMDPlot::WriteOpenPMDParticles (const amrex::Vector<ParticleDiag>& particle_diags, - const bool isBTD) + const bool isBTD, const amrex::Vector<int>& totalParticlesFlushedAlready) { WARPX_PROFILE("WarpXOpenPMDPlot::WriteOpenPMDParticles()"); @@ -552,28 +552,48 @@ WarpXOpenPMDPlot::WriteOpenPMDParticles (const amrex::Vector<ParticleDiag>& part GeometryFilter const geometry_filter(particle_diags[i].m_do_geom_filter, particle_diags[i].m_diag_domain); - using SrcData = WarpXParticleContainer::ParticleTileType::ConstParticleTileDataType; - tmp.copyParticles(*pc, - [=] AMREX_GPU_HOST_DEVICE (const SrcData& src, int ip, const amrex::RandomEngine& engine) - { - const SuperParticleType& p = src.getSuperParticle(ip); - return random_filter(p, engine) * uniform_filter(p, engine) - * parser_filter(p, engine) * geometry_filter(p, engine); - }, true); + if (! isBTD) { + using SrcData = WarpXParticleContainer::ParticleTileType::ConstParticleTileDataType; + tmp.copyParticles(*pc, + [=] AMREX_GPU_HOST_DEVICE (const SrcData& src, int ip, const amrex::RandomEngine& engine) + { + const SuperParticleType& p = src.getSuperParticle(ip); + return random_filter(p, engine) * uniform_filter(p, engine) + * parser_filter(p, engine) * geometry_filter(p, engine); + }, true); + } else if (isBTD) { + PinnedMemoryParticleContainer* pinned_pc = particle_diags[i].getPinnedParticleContainer(); + tmp.SetParticleGeometry(0,pinned_pc->Geom(0)); + tmp.SetParticleBoxArray(0,pinned_pc->ParticleBoxArray(0)); + tmp.SetParticleDistributionMap(0, pinned_pc->ParticleDistributionMap(0)); + tmp.copyParticles(*pinned_pc); + } // real_names contains a list of all real particle attributes. // real_flags is 1 or 0, whether quantity is dumped or not. { - DumpToFile(&tmp, - particle_diags[i].getSpeciesName(), - m_CurrentStep, - real_flags, - int_flags, - real_names, int_names, - pc->getCharge(), pc->getMass(), - isBTD - ); + if (isBTD) { + DumpToFile(&tmp, + particle_diags[i].getSpeciesName(), + m_CurrentStep, + real_flags, + int_flags, + real_names, int_names, + pc->getCharge(), pc->getMass(), + isBTD, totalParticlesFlushedAlready[i] + ); + } else { + DumpToFile(&tmp, + particle_diags[i].getSpeciesName(), + m_CurrentStep, + real_flags, + int_flags, + real_names, int_names, + pc->getCharge(), pc->getMass(), + isBTD, 0 + ); + } } // Convert momentum back to WarpX units @@ -590,12 +610,13 @@ WarpXOpenPMDPlot::DumpToFile (ParticleContainer* pc, const amrex::Vector<std::string>& real_comp_names, const amrex::Vector<std::string>& int_comp_names, amrex::ParticleReal const charge, - amrex::ParticleReal const mass, - const bool isBTD) const + amrex::ParticleReal const mass, const bool isBTD, + int ParticleFlushOffset) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_Series != nullptr, "openPMD: series must be initialized"); WarpXParticleCounter counter(pc); + if (counter.GetTotalNumParticles() == 0) return; openPMD::Iteration currIteration = GetIteration(iteration, isBTD); openPMD::ParticleSpecies currSpecies = currIteration.particles[name]; @@ -640,20 +661,29 @@ WarpXOpenPMDPlot::DumpToFile (ParticleContainer* pc, // // define positions & offsets // - SetupPos(currSpecies, counter.GetTotalNumParticles(), charge, mass); - SetupRealProperties(currSpecies, write_real_comp, real_comp_names, write_int_comp, int_comp_names, counter.GetTotalNumParticles()); - + const unsigned long long NewParticleVectorSize = counter.GetTotalNumParticles() + ParticleFlushOffset; + m_doParticleSetUp = false; + if (counter.GetTotalNumParticles() > 0 and ParticleFlushOffset == 0) { + // This will trigger meta-data flush for particles the first-time non-zero number of particles are flushed. + m_doParticleSetUp = true; + } + SetupPos(currSpecies, NewParticleVectorSize, charge, mass, isBTD); + SetupRealProperties(currSpecies, write_real_comp, real_comp_names, write_int_comp, int_comp_names, NewParticleVectorSize, isBTD); // open files from all processors, in case some will not contribute below m_Series->flush(); - for (auto currentLevel = 0; currentLevel <= pc->finestLevel(); currentLevel++) { uint64_t offset = static_cast<uint64_t>( counter.m_ParticleOffsetAtRank[currentLevel] ); - + // For BTD, the offset include the number of particles already flushed + if (isBTD) offset += ParticleFlushOffset; for (ParticleIter pti(*pc, currentLevel); pti.isValid(); ++pti) { auto const numParticleOnTile = pti.numParticles(); uint64_t const numParticleOnTile64 = static_cast<uint64_t>( numParticleOnTile ); + // Do not call storeChunk() with zero-sized particle tiles: + // https://github.com/openPMD/openPMD-api/issues/1147 + if (numParticleOnTile == 0) continue; + // get position and particle ID from aos // note: this implementation iterates the AoS 4x... // if we flush late as we do now, we can also copy out the data in one go @@ -740,11 +770,12 @@ WarpXOpenPMDPlot::SetupRealProperties (openPMD::ParticleSpecies& currSpecies, const amrex::Vector<std::string>& real_comp_names, const amrex::Vector<int>& write_int_comp, const amrex::Vector<std::string>& int_comp_names, - unsigned long long np) const + const unsigned long long np, bool const isBTD) const { - auto dtype_real = openPMD::Dataset(openPMD::determineDatatype<amrex::ParticleReal>(), {np}); - auto dtype_int = openPMD::Dataset(openPMD::determineDatatype<int>(), {np}); - + std::string options = "{}"; + if (isBTD) options = "{ \"resizable\": true }"; + auto dtype_real = openPMD::Dataset(openPMD::determineDatatype<amrex::ParticleReal>(), {np}, options); + auto dtype_int = openPMD::Dataset(openPMD::determineDatatype<int>(), {np}, options); // // the beam/input3d showed write_real_comp.size() = 16 while only 10 real comp names // so using the min to be safe. @@ -767,6 +798,8 @@ WarpXOpenPMDPlot::SetupRealProperties (openPMD::ParticleSpecies& currSpecies, } } + // attributes need to be set only the first time BTD flush is called for a snapshot + if (isBTD and m_doParticleSetUp == false) return; std::set< std::string > addedRecords; // add meta-data per record only once for (auto idx=0; idx<m_NumSoARealAttributes; idx++) { auto ii = m_NumAoSRealAttributes + idx; // jump over AoS names @@ -832,7 +865,6 @@ WarpXOpenPMDPlot::SaveRealProperty (ParticleIter& pti, uint64_t const numParticleOnTile64 = static_cast<uint64_t>( numParticleOnTile ); auto const& aos = pti.GetArrayOfStructs(); // size = numParticlesOnTile auto const& soa = pti.GetStructOfArrays(); - // first we concatinate the AoS into contiguous arrays { for( auto idx=0; idx<m_NumAoSRealAttributes; idx++ ) { @@ -892,23 +924,31 @@ WarpXOpenPMDPlot::SetupPos ( openPMD::ParticleSpecies& currSpecies, const unsigned long long& np, amrex::ParticleReal const charge, - amrex::ParticleReal const mass) const + amrex::ParticleReal const mass, + bool const isBTD) { - auto const realType = openPMD::Dataset(openPMD::determineDatatype<amrex::ParticleReal>(), {np}); - auto const idType = openPMD::Dataset(openPMD::determineDatatype< uint64_t >(), {np}); + std::string options = "{}"; + if (isBTD) options = "{ \"resizable\": true }"; + auto realType= openPMD::Dataset(openPMD::determineDatatype<amrex::ParticleReal>(), {np}, options); + auto idType = openPMD::Dataset(openPMD::determineDatatype< uint64_t >(), {np}, options); auto const positionComponents = detail::getParticlePositionComponentLabels(); for( auto const& comp : positionComponents ) { currSpecies["positionOffset"][comp].resetDataset( realType ); - currSpecies["positionOffset"][comp].makeConstant( 0. ); currSpecies["position"][comp].resetDataset( realType ); } auto const scalar = openPMD::RecordComponent::SCALAR; currSpecies["id"][scalar].resetDataset( idType ); currSpecies["charge"][scalar].resetDataset( realType ); - currSpecies["charge"][scalar].makeConstant( charge ); currSpecies["mass"][scalar].resetDataset( realType ); + + if (isBTD and m_doParticleSetUp == false) return; + // make constant + for( auto const& comp : positionComponents ) { + currSpecies["positionOffset"][comp].makeConstant( 0. ); + } + currSpecies["charge"][scalar].makeConstant( charge ); currSpecies["mass"][scalar].makeConstant( mass ); // meta data diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 1a6932696..ad3ec8483 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -251,6 +251,13 @@ WarpX::Evolve (int numsteps) myBFD->writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); } + + // sync up time + for (int i = 0; i <= max_level; ++i) { + t_new[i] = cur_time; + } + multi_diags->FilterComputePackFlush( step, false, true ); + bool move_j = is_synchronized; // If is_synchronized we need to shift j too so that next step we can evolve E by dt/2. // We might need to move j because we are going to make a plotfile. @@ -318,11 +325,6 @@ WarpX::Evolve (int numsteps) } } - // sync up time - for (int i = 0; i <= max_level; ++i) { - t_new[i] = cur_time; - } - // warpx_py_afterstep runs with the updated global time. It is included // in the evolve timing. if (warpx_py_afterstep) { @@ -364,7 +366,6 @@ WarpX::Evolve (int numsteps) // End loop on time steps } - multi_diags->FilterComputePackFlushLastTimestep( istep[0] ); if (do_back_transformed_diagnostics) { diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index d45749617..00fb5e8f5 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -246,6 +246,17 @@ public: int nSpeciesBackTransformedDiagnostics() const {return nspecies_back_transformed_diagnostics;} int mapSpeciesBackTransformedDiagnostics(int i) const {return map_species_back_transformed_diagnostics[i];} int doBackTransformedDiagnostics() const {return do_back_transformed_diagnostics;} + /** Whether back-transformed diagnostics need to be performed for any plasma species. + * + * \param[in] do_back_transformed_particles The parameter to set if back-transformed particles are set to true/false + */ + void SetDoBackTransformedParticles (const bool do_back_transformed_particles); + /** Whether back-transformed diagnostics is set for species with species_name. + * + * \param[in] species_name The species for which back-transformed particles is set. + * \param[in] do_back_transformed_particles The parameter to set if back-transformed particles are set to true/false + */ + void SetDoBackTransformedParticles (std::string species_name, const bool do_back_transformed_particles); int nSpeciesDepositOnMainGrid () const { bool const onMainGrid = true; @@ -495,6 +506,7 @@ private: // MultiParticleContainer for 0<i<nspecies_back_transformed_diagnostics std::vector<int> map_species_back_transformed_diagnostics; int do_back_transformed_diagnostics = 0; + bool m_do_back_transformed_particles = false; void MFItInfoCheckTiling(const WarpXParticleContainer& /*pc_src*/) const noexcept { diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 86ff22257..4a5b2ebd8 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -847,6 +847,31 @@ MultiParticleContainer::getSpeciesID (std::string product_str) const } void +MultiParticleContainer::SetDoBackTransformedParticles (const bool do_back_transformed_particles) { + m_do_back_transformed_particles = do_back_transformed_particles; +} + +void +MultiParticleContainer::SetDoBackTransformedParticles (std::string species_name, const bool do_back_transformed_particles) { + auto species_names_list = GetSpeciesNames(); + bool found = 0; + // Loop over species + for (int i = 0; i < static_cast<int>(species_names.size()); ++i) { + // If species name matches, set back-transformed particles parameters + if (species_names_list[i] == species_name) { + found = 1; + auto& pc = allcontainers[i]; + pc->SetDoBackTransformedParticles(do_back_transformed_particles); + } + } + WarpXUtilMsg::AlwaysAssert( + found != 0, + "ERROR: could not find the ID of product species '" + + species_name + "'" + ". Wrong name?" + ); +} + +void MultiParticleContainer::doFieldIonization (int lev, const MultiFab& Ex, const MultiFab& Ey, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index c669d976a..6a001717a 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1605,7 +1605,8 @@ PhysicalParticleContainer::Evolve (int lev, bool has_buffer = cEx || cjx; - if (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics) + if ( (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics) || + (m_do_back_transformed_particles) ) { for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { @@ -2287,6 +2288,8 @@ PhysicalParticleContainer::GetParticleSlice ( uzpold = tmp_particle_data[lev][index][TmpIdx::uzold].dataPtr(); const long np = pti.numParticles(); + amrex::Print() << " np old BTD " << np << "\n"; + amrex::Print() << " tmp particle size : " << tmp_particle_data.size() << "\n"; Real uzfrm = -WarpX::gamma_boost*WarpX::beta_boost*PhysConst::c; Real inv_c2 = 1.0/PhysConst::c/PhysConst::c; @@ -2505,9 +2508,10 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr() + offset; auto copyAttribs = CopyParticleAttribs(pti, tmp_particle_data, offset); - int do_copy = (WarpX::do_back_transformed_diagnostics && - do_back_transformed_diagnostics && - (a_dt_type!=DtType::SecondHalf)); + int do_copy = ( (WarpX::do_back_transformed_diagnostics + && do_back_transformed_diagnostics + && a_dt_type!=DtType::SecondHalf) + || (m_do_back_transformed_particles && (a_dt_type!=DtType::SecondHalf)) ); int* AMREX_RESTRICT ion_lev = nullptr; if (do_field_ionization) { diff --git a/Source/Particles/PinnedMemoryParticleContainer.H b/Source/Particles/PinnedMemoryParticleContainer.H new file mode 100644 index 000000000..f4f3a0b28 --- /dev/null +++ b/Source/Particles/PinnedMemoryParticleContainer.H @@ -0,0 +1,14 @@ +#ifndef PinnedMemoryParticleContainer_H_ +#define PinnedMemoryParticleContainer_H_ + +#include "WarpXParticleContainer.H" +#include <AMReX.H> +#include <AMReX_Vector.H> +#include <AMReX_AmrParticles.H> +#include <AMReX_Particles.H> + +#include <AMReX_ParIter.H> + +using PinnedMemoryParticleContainer = typename amrex::AmrParticleContainer<0,0, PIdx::nattribs, 0, amrex::PinnedArenaAllocator>; + +#endif diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 1561fcfcf..40f352df9 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -317,6 +317,13 @@ public: } int doBackTransformedDiagnostics () const { return do_back_transformed_diagnostics; } + /** Whether back-transformed diagnostics need to be performed for a particular species. + * + * \param[in] do_back_transformed_particles The parameter to set if back-transformed particles are set to true/false + */ + void SetDoBackTransformedParticles(const bool do_back_transformed_particles) { + m_do_back_transformed_particles = do_back_transformed_particles; + } std::map<std::string, int> getParticleComps () const noexcept { return particle_comps;} std::map<std::string, int> getParticleiComps () const noexcept { return particle_icomps;} @@ -417,6 +424,8 @@ protected: int do_resampling = 0; int do_back_transformed_diagnostics = 1; + /** Whether back-transformed diagnostics is turned on for the corresponding species.*/ + bool m_do_back_transformed_particles = false; #ifdef WARPX_QED //Species can receive a shared pointer to a QED engine (species for @@ -445,6 +454,7 @@ public: TmpIdx::nattribs>; using TmpParticles = amrex::Vector<std::map<PairIndex, TmpParticleTile> >; + TmpParticles getTmpParticleData () const noexcept {return tmp_particle_data;} protected: TmpParticles tmp_particle_data; |