diff options
Diffstat (limited to 'Source/Diagnostics/BackTransformedDiagnostic.cpp')
| -rw-r--r-- | Source/Diagnostics/BackTransformedDiagnostic.cpp | 633 |
1 files changed, 351 insertions, 282 deletions
diff --git a/Source/Diagnostics/BackTransformedDiagnostic.cpp b/Source/Diagnostics/BackTransformedDiagnostic.cpp index 2880b37b1..452828f02 100644 --- a/Source/Diagnostics/BackTransformedDiagnostic.cpp +++ b/Source/Diagnostics/BackTransformedDiagnostic.cpp @@ -459,7 +459,7 @@ namespace bool compare_tlab_uptr(const std::unique_ptr<LabFrameDiag>&a, const std::unique_ptr<LabFrameDiag>&b) { - return a->t_lab < b->t_lab; + return a->m_t_lab < b->m_t_lab; } namespace @@ -468,7 +468,7 @@ void LorentzTransformZ(MultiFab& data, Real gamma_boost, Real beta_boost, int ncomp) { // Loop over tiles/boxes and in-place convert each slice from boosted - // frame to lab frame. + // frame to back-transformed lab frame. #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -520,14 +520,16 @@ BackTransformedDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, Real dt_slice_snapshots_lab, int N_slice_snapshots, Real gamma_boost, Real t_boost, Real dt_boost, int boost_direction, const Geometry& geom, - amrex::RealBox& slice_realbox) - : gamma_boost_(gamma_boost), - dt_snapshots_lab_(dt_snapshots_lab), - dt_boost_(dt_boost), - N_snapshots_(N_snapshots), - boost_direction_(boost_direction), - N_slice_snapshots_(N_slice_snapshots), - dt_slice_snapshots_lab_(dt_slice_snapshots_lab) + amrex::RealBox& slice_realbox, + amrex::Real particle_slice_width_lab) + : m_gamma_boost_(gamma_boost), + m_dt_snapshots_lab_(dt_snapshots_lab), + m_dt_boost_(dt_boost), + m_N_snapshots_(N_snapshots), + m_boost_direction_(boost_direction), + m_N_slice_snapshots_(N_slice_snapshots), + m_dt_slice_snapshots_lab_(dt_slice_snapshots_lab), + m_particle_slice_width_lab_(particle_slice_width_lab) { @@ -536,13 +538,13 @@ BackTransformedDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, AMREX_ALWAYS_ASSERT(WarpX::do_back_transformed_fields or WarpX::do_back_transformed_particles); - inv_gamma_boost_ = 1.0 / gamma_boost_; - beta_boost_ = std::sqrt(1.0 - inv_gamma_boost_*inv_gamma_boost_); - inv_beta_boost_ = 1.0 / beta_boost_; + m_inv_gamma_boost_ = 1.0 / m_gamma_boost_; + m_beta_boost_ = std::sqrt(1.0 - m_inv_gamma_boost_*m_inv_gamma_boost_); + m_inv_beta_boost_ = 1.0 / m_beta_boost_; - dz_lab_ = PhysConst::c * dt_boost_ * inv_beta_boost_ * inv_gamma_boost_; - inv_dz_lab_ = 1.0 / dz_lab_; - int Nz_lab = static_cast<unsigned>((zmax_lab - zmin_lab) * inv_dz_lab_); + m_dz_lab_ = PhysConst::c * m_dt_boost_ * m_inv_beta_boost_ * m_inv_gamma_boost_; + m_inv_dz_lab_ = 1.0 / m_dz_lab_; + int Nz_lab = static_cast<unsigned>((zmax_lab - zmin_lab) * m_inv_dz_lab_); int Nx_lab = geom.Domain().length(0); #if (AMREX_SPACEDIM == 3) int Ny_lab = geom.Domain().length(1); @@ -565,40 +567,38 @@ BackTransformedDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, map_actual_fields_to_dump.push_back(i); if (do_user_fields){ - ncomp_to_dump = user_fields_to_dump.size(); - map_actual_fields_to_dump.resize(ncomp_to_dump); - mesh_field_names.resize(ncomp_to_dump); - for (int i=0; i<ncomp_to_dump; i++){ + m_ncomp_to_dump = user_fields_to_dump.size(); + map_actual_fields_to_dump.resize(m_ncomp_to_dump); + m_mesh_field_names.resize(m_ncomp_to_dump); + for (int i=0; i<m_ncomp_to_dump; i++){ std::string fieldstr = user_fields_to_dump[i]; - mesh_field_names[i] = fieldstr; - map_actual_fields_to_dump[i] = possible_fields_to_dump[fieldstr]; + m_mesh_field_names[i] = fieldstr; + map_actual_fields_to_dump[i] = m_possible_fields_to_dump[fieldstr]; } } // allocating array with total number of lab frame diags (snapshots+slices) - LabFrameDiags_.resize(N_snapshots+N_slice_snapshots); + m_LabFrameDiags_.resize(N_snapshots+N_slice_snapshots); for (int i = 0; i < N_snapshots; ++i) { - Real t_lab = i * dt_snapshots_lab_; + Real t_lab = i * m_dt_snapshots_lab_; // Get simulation domain physical coordinates (in boosted frame). RealBox prob_domain_lab = geom.ProbDomain(); // Replace z bounds by lab-frame coordinates - // x and y bounds are the same for lab frame and boosted frame + // x and y bounds are the same for back-transformed lab frame and boosted frame prob_domain_lab.setLo(AMREX_SPACEDIM-1, zmin_lab + v_window_lab * t_lab); prob_domain_lab.setHi(AMREX_SPACEDIM-1, zmax_lab + v_window_lab * t_lab); Box diag_box = geom.Domain(); - LabFrameDiags_[i].reset(new LabFrameSnapShot(t_lab, t_boost, - inv_gamma_boost_, inv_beta_boost_, dz_lab_, + m_LabFrameDiags_[i].reset(new LabFrameSnapShot(t_lab, t_boost, + m_inv_gamma_boost_, m_inv_beta_boost_, m_dz_lab_, prob_domain_lab, prob_ncells_lab, - ncomp_to_dump, mesh_field_names, prob_domain_lab, + m_ncomp_to_dump, m_mesh_field_names, prob_domain_lab, diag_box, i)); } for (int i = 0; i < N_slice_snapshots; ++i) { - amrex::Real cell_dx = 0; - amrex::Real cell_dy = 0; IntVect slice_ncells_lab ; // To construct LabFrameSlice(), the location of lo() and hi() of the @@ -611,18 +611,17 @@ BackTransformedDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, const amrex::Real* current_slice_hi = slice_realbox.hi(); const amrex::Real zmin_slice_lab = current_slice_lo[AMREX_SPACEDIM-1] / - ( (1.+beta_boost_)*gamma_boost_); + ( (1.+m_beta_boost_)*m_gamma_boost_); const amrex::Real zmax_slice_lab = current_slice_hi[AMREX_SPACEDIM-1] / - ( (1.+beta_boost_)*gamma_boost_); + ( (1.+m_beta_boost_)*m_gamma_boost_); int Nz_slice_lab = static_cast<unsigned>(( - zmax_slice_lab - zmin_slice_lab) * inv_dz_lab_); + zmax_slice_lab - zmin_slice_lab) * m_inv_dz_lab_); int Nx_slice_lab = ( current_slice_hi[0] - current_slice_lo[0] ) / geom.CellSize(0); if (Nx_slice_lab == 0 ) Nx_slice_lab = 1; // if the x-dimension is reduced, increase total_cells by 1 // to be consistent with the number of cells created for the output. if (Nx_lab != Nx_slice_lab) Nx_slice_lab++; - cell_dx = geom.CellSize(0); #if (AMREX_SPACEDIM == 3) int Ny_slice_lab = ( current_slice_hi[1] - current_slice_lo[1]) / geom.CellSize(1); @@ -631,7 +630,6 @@ BackTransformedDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, // to be consistent with the number of cells created for the output. if (Ny_lab != Ny_slice_lab) Ny_slice_lab++; slice_ncells_lab = {Nx_slice_lab, Ny_slice_lab, Nz_slice_lab}; - cell_dy = geom.CellSize(1); #else slice_ncells_lab = {Nx_slice_lab, Nz_slice_lab}; #endif @@ -653,7 +651,7 @@ BackTransformedDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, Box stmp(slice_lo,slice_hi); Box slicediag_box = stmp; - Real t_slice_lab = i * dt_slice_snapshots_lab_ ; + Real t_slice_lab = i * m_dt_slice_snapshots_lab_ ; RealBox prob_domain_lab = geom.ProbDomain(); // replace z bounds by lab-frame coordinates prob_domain_lab.setLo(AMREX_SPACEDIM-1, zmin_lab + v_window_lab * t_slice_lab); @@ -667,16 +665,16 @@ BackTransformedDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, v_window_lab * t_slice_lab ); // construct labframeslice - LabFrameDiags_[i+N_snapshots].reset(new LabFrameSlice(t_slice_lab, t_boost, - inv_gamma_boost_, inv_beta_boost_, dz_lab_, + m_LabFrameDiags_[i+N_snapshots].reset(new LabFrameSlice(t_slice_lab, t_boost, + m_inv_gamma_boost_, m_inv_beta_boost_, m_dz_lab_, prob_domain_lab, slice_ncells_lab, - ncomp_to_dump, mesh_field_names, slice_dom_lab, - slicediag_box, i, cell_dx, cell_dy)); + m_ncomp_to_dump, m_mesh_field_names, slice_dom_lab, + slicediag_box, i, m_particle_slice_width_lab_)); } // sort diags based on their respective t_lab - std::stable_sort(LabFrameDiags_.begin(), LabFrameDiags_.end(), compare_tlab_uptr); + std::stable_sort(m_LabFrameDiags_.begin(), m_LabFrameDiags_.end(), compare_tlab_uptr); - AMREX_ALWAYS_ASSERT(max_box_size_ >= num_buffer_); + AMREX_ALWAYS_ASSERT(m_max_box_size_ >= m_num_buffer_); } void BackTransformedDiagnostic::Flush(const Geometry& geom) @@ -690,42 +688,42 @@ void BackTransformedDiagnostic::Flush(const Geometry& geom) const std::vector<std::string> species_names = mypc.GetSpeciesNames(); // Loop over BFD snapshots - for (int i = 0; i < LabFrameDiags_.size(); ++i) { + for (int i = 0; i < m_LabFrameDiags_.size(); ++i) { - Real zmin_lab = LabFrameDiags_[i]->prob_domain_lab_.lo(AMREX_SPACEDIM-1); - int i_lab = (LabFrameDiags_[i]->current_z_lab - zmin_lab) / dz_lab_; + Real zmin_lab = m_LabFrameDiags_[i]->m_prob_domain_lab_.lo(AMREX_SPACEDIM-1); + int i_lab = (m_LabFrameDiags_[i]->m_current_z_lab - zmin_lab) / m_dz_lab_; - if (LabFrameDiags_[i]->buff_counter_ != 0) { + if (m_LabFrameDiags_[i]->m_buff_counter_ != 0) { if (WarpX::do_back_transformed_fields) { - const BoxArray& ba = LabFrameDiags_[i]->data_buffer_->boxArray(); - const int hi = ba[0].bigEnd(boost_direction_); - const int lo = hi - LabFrameDiags_[i]->buff_counter_ + 1; + const BoxArray& ba = m_LabFrameDiags_[i]->m_data_buffer_->boxArray(); + const int hi = ba[0].bigEnd(m_boost_direction_); + const int lo = hi - m_LabFrameDiags_[i]->m_buff_counter_ + 1; //Box buff_box = geom.Domain(); - Box buff_box = LabFrameDiags_[i]->buff_box_; - buff_box.setSmall(boost_direction_, lo); - buff_box.setBig(boost_direction_, hi); + Box buff_box = m_LabFrameDiags_[i]->m_buff_box_; + buff_box.setSmall(m_boost_direction_, lo); + buff_box.setBig(m_boost_direction_, hi); BoxArray buff_ba(buff_box); - buff_ba.maxSize(max_box_size_); + buff_ba.maxSize(m_max_box_size_); DistributionMapping buff_dm(buff_ba); - const int ncomp = LabFrameDiags_[i]->data_buffer_->nComp(); + const int ncomp = m_LabFrameDiags_[i]->m_data_buffer_->nComp(); MultiFab tmp(buff_ba, buff_dm, ncomp, 0); - tmp.copy(*LabFrameDiags_[i]->data_buffer_, 0, 0, ncomp); + tmp.copy(*m_LabFrameDiags_[i]->m_data_buffer_, 0, 0, ncomp); #ifdef WARPX_USE_HDF5 for (int comp = 0; comp < ncomp; ++comp) { - output_write_field(LabFrameDiags_[i]->file_name, - mesh_field_names[comp], tmp, comp, + output_write_field(m_LabFrameDiags_[i]->m_file_name, + m_mesh_field_names[comp], tmp, comp, lbound(buff_box).x, lbound(buff_box).y, lbound(buff_box).z); } #else std::stringstream ss; - ss << LabFrameDiags_[i]->file_name << "/Level_0/" + ss << m_LabFrameDiags_[i]->m_file_name << "/Level_0/" << Concatenate("buffer", i_lab, 5); VisMF::Write(tmp, ss.str()); #endif @@ -739,21 +737,21 @@ void BackTransformedDiagnostic::Flush(const Geometry& geom) species_names[mypc.mapSpeciesBackTransformedDiagnostics(j)]; #ifdef WARPX_USE_HDF5 // Dump species data - writeParticleDataHDF5(LabFrameDiags_[i]->particles_buffer_[j], - LabFrameDiags_[i]->file_name, + writeParticleDataHDF5(m_LabFrameDiags_[i]->m_particles_buffer_[j], + m_LabFrameDiags_[i]->m_file_name, species_name); #else std::stringstream part_ss; - part_ss << LabFrameDiags_[i]->file_name + "/" + + part_ss << m_LabFrameDiags_[i]->m_file_name + "/" + species_name + "/"; // Dump species data - writeParticleData(LabFrameDiags_[i]->particles_buffer_[j], + writeParticleData(m_LabFrameDiags_[i]->m_particles_buffer_[j], part_ss.str(), i_lab); #endif } - LabFrameDiags_[i]->particles_buffer_.clear(); + m_LabFrameDiags_[i]->m_particles_buffer_.clear(); } - LabFrameDiags_[i]->buff_counter_ = 0; + m_LabFrameDiags_[i]->m_buff_counter_ = 0; } } @@ -775,8 +773,8 @@ writeLabFrameData(const MultiFab* cell_centered_data, VisMF::SetHeaderVersion(amrex::VisMF::Header::NoFabHeader_v1); const RealBox& domain_z_boost = geom.ProbDomain(); - const Real zlo_boost = domain_z_boost.lo(boost_direction_); - const Real zhi_boost = domain_z_boost.hi(boost_direction_); + const Real zlo_boost = domain_z_boost.lo(m_boost_direction_); + const Real zhi_boost = domain_z_boost.hi(m_boost_direction_); const std::vector<std::string> species_names = mypc.GetSpeciesNames(); Real prev_t_lab = -dt; @@ -785,43 +783,43 @@ writeLabFrameData(const MultiFab* cell_centered_data, amrex::Vector<WarpXParticleContainer::DiagnosticParticleData> tmp_particle_buffer; // Loop over snapshots - for (int i = 0; i < LabFrameDiags_.size(); ++i) { + for (int i = 0; i < m_LabFrameDiags_.size(); ++i) { // Get updated z position of snapshot - const Real old_z_boost = LabFrameDiags_[i]->current_z_boost; - LabFrameDiags_[i]->updateCurrentZPositions(t_boost, - inv_gamma_boost_, - inv_beta_boost_); + const Real old_z_boost = m_LabFrameDiags_[i]->m_current_z_boost; + m_LabFrameDiags_[i]->updateCurrentZPositions(t_boost, + m_inv_gamma_boost_, + m_inv_beta_boost_); - Real diag_zmin_lab = LabFrameDiags_[i]->diag_domain_lab_.lo(AMREX_SPACEDIM-1); - Real diag_zmax_lab = LabFrameDiags_[i]->diag_domain_lab_.hi(AMREX_SPACEDIM-1); + Real diag_zmin_lab = m_LabFrameDiags_[i]->m_diag_domain_lab_.lo(AMREX_SPACEDIM-1); + Real diag_zmax_lab = m_LabFrameDiags_[i]->m_diag_domain_lab_.hi(AMREX_SPACEDIM-1); - if ( ( LabFrameDiags_[i]->current_z_boost < zlo_boost) or - ( LabFrameDiags_[i]->current_z_boost > zhi_boost) or - ( LabFrameDiags_[i]->current_z_lab < diag_zmin_lab) or - ( LabFrameDiags_[i]->current_z_lab > diag_zmax_lab) ) continue; + if ( ( m_LabFrameDiags_[i]->m_current_z_boost < zlo_boost) or + ( m_LabFrameDiags_[i]->m_current_z_boost > zhi_boost) or + ( m_LabFrameDiags_[i]->m_current_z_lab < diag_zmin_lab) or + ( m_LabFrameDiags_[i]->m_current_z_lab > diag_zmax_lab) ) continue; // Get z index of data_buffer_ (i.e. in the lab frame) where // simulation domain (t', [zmin',zmax']), back-transformed to lab // frame, intersects with snapshot. - Real dom_zmin_lab = LabFrameDiags_[i]->prob_domain_lab_.lo(AMREX_SPACEDIM-1); - int i_lab = ( LabFrameDiags_[i]->current_z_lab - dom_zmin_lab) / dz_lab_; + Real dom_zmin_lab = m_LabFrameDiags_[i]->m_prob_domain_lab_.lo(AMREX_SPACEDIM-1); + int i_lab = ( m_LabFrameDiags_[i]->m_current_z_lab - dom_zmin_lab) / m_dz_lab_; // If buffer of snapshot i is empty... - if ( LabFrameDiags_[i]->buff_counter_ == 0) { + if ( m_LabFrameDiags_[i]->m_buff_counter_ == 0) { // ... reset fields buffer data_buffer_ if (WarpX::do_back_transformed_fields) { - LabFrameDiags_[i]->buff_box_.setSmall(boost_direction_, - i_lab - num_buffer_ + 1); - LabFrameDiags_[i]->buff_box_.setBig(boost_direction_, i_lab); + m_LabFrameDiags_[i]->m_buff_box_.setSmall(m_boost_direction_, + i_lab - m_num_buffer_ + 1); + m_LabFrameDiags_[i]->m_buff_box_.setBig(m_boost_direction_, i_lab); - BoxArray buff_ba(LabFrameDiags_[i]->buff_box_); - buff_ba.maxSize(max_box_size_); + BoxArray buff_ba(m_LabFrameDiags_[i]->m_buff_box_); + buff_ba.maxSize(m_max_box_size_); DistributionMapping buff_dm(buff_ba); - LabFrameDiags_[i]->data_buffer_.reset( new MultiFab(buff_ba, - buff_dm, ncomp_to_dump, 0) ); + m_LabFrameDiags_[i]->m_data_buffer_.reset( new MultiFab(buff_ba, + buff_dm, m_ncomp_to_dump, 0) ); } // ... reset particle buffer particles_buffer_[i] if (WarpX::do_back_transformed_particles) - LabFrameDiags_[i]->particles_buffer_.resize( + m_LabFrameDiags_[i]->m_particles_buffer_.resize( mypc.nSpeciesBackTransformedDiagnostics()); } @@ -830,34 +828,34 @@ writeLabFrameData(const MultiFab* cell_centered_data, const int start_comp = 0; const bool interpolate = true; // slice containing back-transformed data is generated only if t_lab != prev_t_lab and is re-used if multiple diags have the same z_lab,t_lab. - if (LabFrameDiags_[i]->t_lab != prev_t_lab ) { + if (m_LabFrameDiags_[i]->m_t_lab != prev_t_lab ) { if (slice) { slice.reset(new MultiFab); slice.reset(nullptr); } - slice = amrex::get_slice_data(boost_direction_, - LabFrameDiags_[i]->current_z_boost, + slice = amrex::get_slice_data(m_boost_direction_, + m_LabFrameDiags_[i]->m_current_z_boost, *cell_centered_data, geom, start_comp, ncomp, interpolate); // Back-transform data to the lab-frame - LorentzTransformZ(*slice, gamma_boost_, beta_boost_, ncomp); + LorentzTransformZ(*slice, m_gamma_boost_, m_beta_boost_, ncomp); } // Create a 2D box for the slice in the boosted frame - Real dx = geom.CellSize(boost_direction_); - int i_boost = ( LabFrameDiags_[i]->current_z_boost - - geom.ProbLo(boost_direction_))/dx; + Real dx = geom.CellSize(m_boost_direction_); + int i_boost = ( m_LabFrameDiags_[i]->m_current_z_boost - + geom.ProbLo(m_boost_direction_))/dx; //Box slice_box = geom.Domain(); - Box slice_box = LabFrameDiags_[i]->buff_box_; - slice_box.setSmall(boost_direction_, i_boost); - slice_box.setBig(boost_direction_, i_boost); + Box slice_box = m_LabFrameDiags_[i]->m_buff_box_; + slice_box.setSmall(m_boost_direction_, i_boost); + slice_box.setBig(m_boost_direction_, i_boost); // Make it a BoxArray slice_ba BoxArray slice_ba(slice_box); - slice_ba.maxSize(max_box_size_); + slice_ba.maxSize(m_max_box_size_); tmp_slice_ptr = std::unique_ptr<MultiFab>(new MultiFab(slice_ba, - LabFrameDiags_[i]->data_buffer_->DistributionMap(), + m_LabFrameDiags_[i]->m_data_buffer_->DistributionMap(), ncomp, 0)); // slice is re-used if the t_lab of a diag is equal to @@ -867,7 +865,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, // tmp_slice_ptr which has the dmap of the // data_buffer that stores the back-transformed data. tmp_slice_ptr->copy(*slice, 0, 0, ncomp); - LabFrameDiags_[i]->AddDataToBuffer(*tmp_slice_ptr, i_lab, + m_LabFrameDiags_[i]->AddDataToBuffer(*tmp_slice_ptr, i_lab, map_actual_fields_to_dump); tmp_slice_ptr.reset(new MultiFab); tmp_slice_ptr.reset(nullptr); @@ -875,44 +873,43 @@ writeLabFrameData(const MultiFab* cell_centered_data, if (WarpX::do_back_transformed_particles) { - if (LabFrameDiags_[i]->t_lab != prev_t_lab ) { + if (m_LabFrameDiags_[i]->m_t_lab != prev_t_lab ) { if (tmp_particle_buffer.size()>0) { tmp_particle_buffer.clear(); tmp_particle_buffer.shrink_to_fit(); } tmp_particle_buffer.resize(mypc.nSpeciesBackTransformedDiagnostics()); - mypc.GetLabFrameData( LabFrameDiags_[i]->file_name, i_lab, - boost_direction_, old_z_boost, - LabFrameDiags_[i]->current_z_boost, - t_boost, LabFrameDiags_[i]->t_lab, dt, - tmp_particle_buffer); + mypc.GetLabFrameData(m_LabFrameDiags_[i]->m_file_name, i_lab, + m_boost_direction_, old_z_boost, + m_LabFrameDiags_[i]->m_current_z_boost, + t_boost, m_LabFrameDiags_[i]->m_t_lab, dt, + tmp_particle_buffer); } - LabFrameDiags_[i]->AddPartDataToParticleBuffer(tmp_particle_buffer, + m_LabFrameDiags_[i]->AddPartDataToParticleBuffer(tmp_particle_buffer, mypc.nSpeciesBackTransformedDiagnostics()); } - ++LabFrameDiags_[i]->buff_counter_; - prev_t_lab = LabFrameDiags_[i]->t_lab; - + ++m_LabFrameDiags_[i]->m_buff_counter_; + prev_t_lab = m_LabFrameDiags_[i]->m_t_lab; // If buffer full, write to disk. - if ( LabFrameDiags_[i]->buff_counter_ == num_buffer_) { + if (m_LabFrameDiags_[i]->m_buff_counter_ == m_num_buffer_) { if (WarpX::do_back_transformed_fields) { #ifdef WARPX_USE_HDF5 - Box buff_box = LabFrameDiags_[i]->buff_box_; - for (int comp = 0; comp < LabFrameDiags_[i]->data_buffer_->nComp(); ++comp) - output_write_field( LabFrameDiags_[i]->file_name, - mesh_field_names[comp], - *LabFrameDiags_[i]->data_buffer_, comp, - lbound(buff_box).x, lbound(buff_box).y, - lbound(buff_box).z); + Box buff_box = m_LabFrameDiags_[i]->m_buff_box_; + for (int comp = 0; comp < m_LabFrameDiags_[i]->m_data_buffer_->nComp(); ++comp) + output_write_field(m_LabFrameDiags_[i]->m_file_name, + m_mesh_field_names[comp], + *m_LabFrameDiags_[i]->m_data_buffer_, comp, + lbound(buff_box).x, lbound(buff_box).y, + lbound(buff_box).z); #else std::stringstream mesh_ss; - mesh_ss << LabFrameDiags_[i]->file_name << "/Level_0/" << + mesh_ss << m_LabFrameDiags_[i]->m_file_name << "/Level_0/" << Concatenate("buffer", i_lab, 5); - VisMF::Write( (*LabFrameDiags_[i]->data_buffer_), mesh_ss.str()); + VisMF::Write( (*m_LabFrameDiags_[i]->m_data_buffer_), mesh_ss.str()); #endif } @@ -924,23 +921,23 @@ writeLabFrameData(const MultiFab* cell_centered_data, mypc.mapSpeciesBackTransformedDiagnostics(j)]; #ifdef WARPX_USE_HDF5 // Write data to disk (HDF5) - writeParticleDataHDF5(LabFrameDiags_[i]->particles_buffer_[j], - LabFrameDiags_[i]->file_name, + writeParticleDataHDF5(m_LabFrameDiags_[i]->m_particles_buffer_[j], + m_LabFrameDiags_[i]->m_file_name, species_name); #else std::stringstream part_ss; - part_ss << LabFrameDiags_[i]->file_name + "/" + + part_ss << m_LabFrameDiags_[i]->m_file_name + "/" + species_name + "/"; // Write data to disk (custom) - writeParticleData(LabFrameDiags_[i]->particles_buffer_[j], + writeParticleData(m_LabFrameDiags_[i]->m_particles_buffer_[j], part_ss.str(), i_lab); #endif } - LabFrameDiags_[i]->particles_buffer_.clear(); + m_LabFrameDiags_[i]->m_particles_buffer_.clear(); } - LabFrameDiags_[i]->buff_counter_ = 0; + m_LabFrameDiags_[i]->m_buff_counter_ = 0; } } @@ -1069,12 +1066,12 @@ writeMetaData () HeaderFile.precision(17); - HeaderFile << N_snapshots_ << "\n"; - HeaderFile << dt_snapshots_lab_ << "\n"; - HeaderFile << gamma_boost_ << "\n"; - HeaderFile << beta_boost_ << "\n"; + HeaderFile << m_N_snapshots_ << "\n"; + HeaderFile << m_dt_snapshots_lab_ << "\n"; + HeaderFile << m_gamma_boost_ << "\n"; + HeaderFile << m_beta_boost_ << "\n"; - if (N_slice_snapshots_ > 0) { + if (m_N_slice_snapshots_ > 0) { const std::string fullpath_slice = WarpX::lab_data_directory + "/slices"; if (!UtilCreateDirectory(fullpath_slice, 0755)) CreateDirectoryFailed(fullpath_slice); @@ -1095,10 +1092,10 @@ writeMetaData () HeaderFile_slice.precision(17); - HeaderFile_slice << N_slice_snapshots_ << "\n"; - HeaderFile_slice << dt_slice_snapshots_lab_ << "\n"; - HeaderFile_slice << gamma_boost_ << "\n"; - HeaderFile_slice << beta_boost_ << "\n"; + HeaderFile_slice << m_N_slice_snapshots_ << "\n"; + HeaderFile_slice << m_dt_slice_snapshots_lab_ << "\n"; + HeaderFile_slice << m_gamma_boost_ << "\n"; + HeaderFile_slice << m_beta_boost_ << "\n"; } @@ -1114,35 +1111,35 @@ LabFrameSnapShot(Real t_lab_in, Real t_boost, Real inv_gamma_boost_in, std::vector<std::string> mesh_field_names, amrex::RealBox diag_domain_lab, Box diag_box, int file_num_in) { - t_lab = t_lab_in; - dz_lab_ = dz_lab_in; - inv_gamma_boost_ = inv_gamma_boost_in; - inv_beta_boost_ = inv_beta_boost_in; - prob_domain_lab_ = prob_domain_lab; - prob_ncells_lab_ = prob_ncells_lab; - diag_domain_lab_ = diag_domain_lab; - buff_box_ = diag_box; - ncomp_to_dump_ = ncomp_to_dump; - mesh_field_names_ = mesh_field_names; - file_num = file_num_in; - current_z_lab = 0.0; - current_z_boost = 0.0; - updateCurrentZPositions(t_boost, inv_gamma_boost_, inv_beta_boost_); - Real zmin_lab = prob_domain_lab_.lo(AMREX_SPACEDIM-1); - initial_i = (current_z_lab - zmin_lab) / dz_lab_ ; - file_name = Concatenate(WarpX::lab_data_directory + "/snapshots/snapshot", - file_num, 5); + m_t_lab = t_lab_in; + m_dz_lab_ = dz_lab_in; + m_inv_gamma_boost_ = inv_gamma_boost_in; + m_inv_beta_boost_ = inv_beta_boost_in; + m_prob_domain_lab_ = prob_domain_lab; + m_prob_ncells_lab_ = prob_ncells_lab; + m_diag_domain_lab_ = diag_domain_lab; + m_buff_box_ = diag_box; + m_ncomp_to_dump_ = ncomp_to_dump; + m_mesh_field_names_ = mesh_field_names; + m_file_num = file_num_in; + m_current_z_lab = 0.0; + m_current_z_boost = 0.0; + updateCurrentZPositions(t_boost, m_inv_gamma_boost_, m_inv_beta_boost_); + Real zmin_lab = m_prob_domain_lab_.lo(AMREX_SPACEDIM-1); + m_initial_i = (m_current_z_lab - zmin_lab) / m_dz_lab_ ; + m_file_name = Concatenate(WarpX::lab_data_directory + "/snapshots/snapshot", + m_file_num, 5); createLabFrameDirectories(); - buff_counter_ = 0; - if (WarpX::do_back_transformed_fields) data_buffer_.reset(nullptr); + m_buff_counter_ = 0; + if (WarpX::do_back_transformed_fields) m_data_buffer_.reset(nullptr); } void LabFrameDiag:: updateCurrentZPositions(Real t_boost, Real inv_gamma, Real inv_beta) { - current_z_boost = (t_lab*inv_gamma - t_boost)*PhysConst::c*inv_beta; - current_z_lab = (t_lab - t_boost*inv_gamma)*PhysConst::c*inv_beta; + m_current_z_boost = (m_t_lab*inv_gamma - t_boost)*PhysConst::c*inv_beta; + m_current_z_lab = (m_t_lab - t_boost*inv_gamma)*PhysConst::c*inv_beta; } void @@ -1151,7 +1148,7 @@ createLabFrameDirectories() { #ifdef WARPX_USE_HDF5 if (ParallelDescriptor::IOProcessor()) { - output_create(file_name); + output_create(m_file_name); } ParallelDescriptor::Barrier(); @@ -1160,16 +1157,16 @@ createLabFrameDirectories() { { if (WarpX::do_back_transformed_fields) { - const auto lo = lbound(buff_box_); - for (int comp = 0; comp < ncomp_to_dump_; ++comp) { - output_create_field(file_name, mesh_field_names_[comp], - prob_ncells_lab_[0], + const auto lo = lbound(m_buff_box_); + for (int comp = 0; comp < m_ncomp_to_dump_; ++comp) { + output_create_field(m_file_name, m_mesh_field_names_[comp], + m_prob_ncells_lab_[0], #if ( AMREX_SPACEDIM == 3 ) - prob_ncells_lab_[1], + m_prob_ncells_lab_[1], #else 1, #endif - prob_ncells_lab_[AMREX_SPACEDIM-1]+1); + m_prob_ncells_lab_[AMREX_SPACEDIM-1]+1); } } } @@ -1185,23 +1182,23 @@ createLabFrameDirectories() { // Loop over species to be dumped to BFD std::string species_name = species_names[mypc.mapSpeciesBackTransformedDiagnostics(j)]; - output_create_species_group(file_name, species_name); + output_create_species_group(m_file_name, species_name); for (int k = 0; k < static_cast<int>(particle_field_names.size()); ++k) { std::string field_path = species_name + "/" + particle_field_names[k]; - output_create_particle_field(file_name, field_path); + output_create_particle_field(m_file_name, field_path); } } } #else if (ParallelDescriptor::IOProcessor()) { - if (!UtilCreateDirectory(file_name, 0755)) - CreateDirectoryFailed(file_name); + if (!UtilCreateDirectory(m_file_name, 0755)) + CreateDirectoryFailed(m_file_name); const int nlevels = 1; for(int i = 0; i < nlevels; ++i) { - const std::string &fullpath = LevelFullPath(i, file_name); + const std::string &fullpath = LevelFullPath(i, m_file_name); if (!UtilCreateDirectory(fullpath, 0755)) CreateDirectoryFailed(fullpath); } @@ -1215,7 +1212,7 @@ createLabFrameDirectories() { // Get species name std::string species_name = species_names[mypc.mapSpeciesBackTransformedDiagnostics(i)]; - const std::string fullpath = file_name + "/" + species_name; + const std::string fullpath = m_file_name + "/" + species_name; if (!UtilCreateDirectory(fullpath, 0755)) CreateDirectoryFailed(fullpath); } @@ -1234,7 +1231,7 @@ writeLabFrameHeader() { VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size); std::ofstream HeaderFile; HeaderFile.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size()); - std::string HeaderFileName(file_name + "/Header"); + std::string HeaderFileName(m_file_name + "/Header"); HeaderFile.open(HeaderFileName.c_str(), std::ofstream::out | std::ofstream::trunc | std::ofstream::binary); @@ -1243,30 +1240,30 @@ writeLabFrameHeader() { HeaderFile.precision(17); - HeaderFile << t_lab << "\n"; + HeaderFile << m_t_lab << "\n"; // Write domain number of cells - HeaderFile << prob_ncells_lab_[0] << ' ' + HeaderFile << m_prob_ncells_lab_[0] << ' ' #if ( AMREX_SPACEDIM==3 ) - << prob_ncells_lab_[1] << ' ' + << m_prob_ncells_lab_[1] << ' ' #endif - << prob_ncells_lab_[AMREX_SPACEDIM-1] <<'\n'; + << m_prob_ncells_lab_[AMREX_SPACEDIM-1] <<'\n'; // Write domain physical boundaries // domain lower bound - HeaderFile << diag_domain_lab_.lo(0) << ' ' + HeaderFile << m_diag_domain_lab_.lo(0) << ' ' #if ( AMREX_SPACEDIM==3 ) - << diag_domain_lab_.lo(1) << ' ' + << m_diag_domain_lab_.lo(1) << ' ' #endif - << diag_domain_lab_.lo(AMREX_SPACEDIM-1) <<'\n'; + << m_diag_domain_lab_.lo(AMREX_SPACEDIM-1) <<'\n'; // domain higher bound - HeaderFile << diag_domain_lab_.hi(0) << ' ' + HeaderFile << m_diag_domain_lab_.hi(0) << ' ' #if ( AMREX_SPACEDIM==3 ) - << diag_domain_lab_.hi(1) << ' ' + << m_diag_domain_lab_.hi(1) << ' ' #endif - << diag_domain_lab_.hi(AMREX_SPACEDIM-1) <<'\n'; + << m_diag_domain_lab_.hi(AMREX_SPACEDIM-1) <<'\n'; // List of fields dumped to file - for (int i=0; i<ncomp_to_dump_; i++) + for (int i=0; i<m_ncomp_to_dump_; i++) { - HeaderFile << mesh_field_names_[i] << ' '; + HeaderFile << m_mesh_field_names_[i] << ' '; } HeaderFile << "\n"; } @@ -1281,28 +1278,30 @@ LabFrameSlice(Real t_lab_in, Real t_boost, Real inv_gamma_boost_in, IntVect prob_ncells_lab, int ncomp_to_dump, std::vector<std::string> mesh_field_names, RealBox diag_domain_lab, Box diag_box, int file_num_in, - amrex::Real cell_dx, amrex::Real cell_dy) + amrex::Real particle_slice_dx_lab) { - t_lab = t_lab_in; - prob_domain_lab_ = prob_domain_lab; - prob_ncells_lab_ = prob_ncells_lab; - diag_domain_lab_ = diag_domain_lab; - buff_box_ = diag_box; - ncomp_to_dump_ = ncomp_to_dump; - mesh_field_names_ = mesh_field_names; - file_num = file_num_in; - current_z_lab = 0.0; - current_z_boost = 0.0; - updateCurrentZPositions(t_boost, inv_gamma_boost_, inv_beta_boost_); - Real zmin_lab = prob_domain_lab_.lo(AMREX_SPACEDIM-1); - initial_i = (current_z_lab - zmin_lab)/dz_lab_; - file_name = Concatenate(WarpX::lab_data_directory+"/slices/slice",file_num,5); + m_t_lab = t_lab_in; + m_dz_lab_ = dz_lab_in; + m_inv_gamma_boost_ = inv_gamma_boost_in; + m_inv_beta_boost_ = inv_beta_boost_in; + m_prob_domain_lab_ = prob_domain_lab; + m_prob_ncells_lab_ = prob_ncells_lab; + m_diag_domain_lab_ = diag_domain_lab; + m_buff_box_ = diag_box; + m_ncomp_to_dump_ = ncomp_to_dump; + m_mesh_field_names_ = mesh_field_names; + m_file_num = file_num_in; + m_current_z_lab = 0.0; + m_current_z_boost = 0.0; + updateCurrentZPositions(t_boost, m_inv_gamma_boost_, m_inv_beta_boost_); + Real zmin_lab = m_prob_domain_lab_.lo(AMREX_SPACEDIM-1); + m_initial_i = (m_current_z_lab - zmin_lab)/m_dz_lab_; + m_file_name = Concatenate(WarpX::lab_data_directory+"/slices/slice",m_file_num,5); createLabFrameDirectories(); - buff_counter_ = 0; - dx_ = cell_dx; - dy_ = cell_dy; + m_buff_counter_ = 0; + m_particle_slice_dx_lab_ = particle_slice_dx_lab; - if (WarpX::do_back_transformed_fields) data_buffer_.reset(nullptr); + if (WarpX::do_back_transformed_fields) m_data_buffer_.reset(nullptr); } void @@ -1311,7 +1310,7 @@ AddDataToBuffer( MultiFab& tmp, int k_lab, amrex::Gpu::ManagedDeviceVector<int> map_actual_fields_to_dump) { const int ncomp_to_dump = map_actual_fields_to_dump.size(); - MultiFab& buf = *data_buffer_; + MultiFab& buf = *m_data_buffer_; for (MFIter mfi(tmp, TilingIfNotGPU()); mfi.isValid(); ++mfi) { Array4<Real> tmp_arr = tmp[mfi].array(); Array4<Real> buf_arr = buf[mfi].array(); @@ -1340,10 +1339,10 @@ AddDataToBuffer( MultiFab& tmp, int k_lab, amrex::Gpu::ManagedDeviceVector<int> map_actual_fields_to_dump) { const int ncomp_to_dump = map_actual_fields_to_dump.size(); - MultiFab& buf = *data_buffer_; + MultiFab& buf = *m_data_buffer_; for (MFIter mfi(tmp, TilingIfNotGPU()); mfi.isValid(); ++mfi) { - Box& bx = buff_box_; + Box& bx = m_buff_box_; const Box& bx_bf = mfi.tilebox(); bx.setSmall(AMREX_SPACEDIM-1,bx_bf.smallEnd(AMREX_SPACEDIM-1)); bx.setBig(AMREX_SPACEDIM-1,bx_bf.bigEnd(AMREX_SPACEDIM-1)); @@ -1374,40 +1373,56 @@ AddPartDataToParticleBuffer( auto np = tmp_particle_buffer[isp].GetRealData(DiagIdx::w).size(); if (np == 0) return; - particles_buffer_[isp].GetRealData(DiagIdx::w).insert( - particles_buffer_[isp].GetRealData(DiagIdx::w).end(), - tmp_particle_buffer[isp].GetRealData(DiagIdx::w).begin(), - tmp_particle_buffer[isp].GetRealData(DiagIdx::w).end()); - - particles_buffer_[isp].GetRealData(DiagIdx::x).insert( - particles_buffer_[isp].GetRealData(DiagIdx::x).end(), - tmp_particle_buffer[isp].GetRealData(DiagIdx::x).begin(), - tmp_particle_buffer[isp].GetRealData(DiagIdx::x).end()); - - particles_buffer_[isp].GetRealData(DiagIdx::y).insert( - particles_buffer_[isp].GetRealData(DiagIdx::y).end(), - tmp_particle_buffer[isp].GetRealData(DiagIdx::y).begin(), - tmp_particle_buffer[isp].GetRealData(DiagIdx::y).end()); - - particles_buffer_[isp].GetRealData(DiagIdx::z).insert( - particles_buffer_[isp].GetRealData(DiagIdx::z).end(), - tmp_particle_buffer[isp].GetRealData(DiagIdx::z).begin(), - tmp_particle_buffer[isp].GetRealData(DiagIdx::z).end()); - - particles_buffer_[isp].GetRealData(DiagIdx::ux).insert( - particles_buffer_[isp].GetRealData(DiagIdx::ux).end(), - tmp_particle_buffer[isp].GetRealData(DiagIdx::ux).begin(), - tmp_particle_buffer[isp].GetRealData(DiagIdx::ux).end()); - - particles_buffer_[isp].GetRealData(DiagIdx::uy).insert( - particles_buffer_[isp].GetRealData(DiagIdx::uy).end(), - tmp_particle_buffer[isp].GetRealData(DiagIdx::uy).begin(), - tmp_particle_buffer[isp].GetRealData(DiagIdx::uy).end()); - - particles_buffer_[isp].GetRealData(DiagIdx::uz).insert( - particles_buffer_[isp].GetRealData(DiagIdx::uz).end(), - tmp_particle_buffer[isp].GetRealData(DiagIdx::uz).begin(), - tmp_particle_buffer[isp].GetRealData(DiagIdx::uz).end()); + // allocate size of particle buffer array to np + // This is a growing array. Each time we add np elements + // to the existing array which has size = init_size + const int init_size = m_particles_buffer_[isp].GetRealData(DiagIdx::w).size(); + const int total_size = init_size + np; + m_particles_buffer_[isp].resize(total_size); + + // Data pointers to particle attributes // + Real* const AMREX_RESTRICT wp_buff = + m_particles_buffer_[isp].GetRealData(DiagIdx::w).data(); + Real* const AMREX_RESTRICT x_buff = + m_particles_buffer_[isp].GetRealData(DiagIdx::x).data(); + Real* const AMREX_RESTRICT y_buff = + m_particles_buffer_[isp].GetRealData(DiagIdx::y).data(); + Real* const AMREX_RESTRICT z_buff = + m_particles_buffer_[isp].GetRealData(DiagIdx::z).data(); + Real* const AMREX_RESTRICT ux_buff = + m_particles_buffer_[isp].GetRealData(DiagIdx::ux).data(); + Real* const AMREX_RESTRICT uy_buff = + m_particles_buffer_[isp].GetRealData(DiagIdx::uy).data(); + Real* const AMREX_RESTRICT uz_buff = + m_particles_buffer_[isp].GetRealData(DiagIdx::uz).data(); + + Real const* const AMREX_RESTRICT wp_temp = + tmp_particle_buffer[isp].GetRealData(DiagIdx::w).data(); + Real const* const AMREX_RESTRICT x_temp = + tmp_particle_buffer[isp].GetRealData(DiagIdx::x).data(); + Real const* const AMREX_RESTRICT y_temp = + tmp_particle_buffer[isp].GetRealData(DiagIdx::y).data(); + Real const* const AMREX_RESTRICT z_temp = + tmp_particle_buffer[isp].GetRealData(DiagIdx::z).data(); + Real const* const AMREX_RESTRICT ux_temp = + tmp_particle_buffer[isp].GetRealData(DiagIdx::ux).data(); + Real const* const AMREX_RESTRICT uy_temp = + tmp_particle_buffer[isp].GetRealData(DiagIdx::uy).data(); + Real const* const AMREX_RESTRICT uz_temp = + tmp_particle_buffer[isp].GetRealData(DiagIdx::uz).data(); + + // copy all the particles from tmp to buffer + amrex::ParallelFor(np, + [=] AMREX_GPU_DEVICE(int i) + { + wp_buff[init_size + i] = wp_temp[i]; + x_buff[init_size + i] = x_temp[i]; + y_buff[init_size + i] = y_temp[i]; + z_buff[init_size + i] = z_temp[i]; + ux_buff[init_size + i] = ux_temp[i]; + uy_buff[init_size + i] = uy_temp[i]; + uz_buff[init_size + i] = uz_temp[i]; + }); } } @@ -1415,53 +1430,107 @@ void LabFrameSlice:: AddPartDataToParticleBuffer( Vector<WarpXParticleContainer::DiagnosticParticleData> tmp_particle_buffer, - int nSpeciesBoostedFrame) { - for (int isp = 0; isp < nSpeciesBoostedFrame; ++isp) { + int nSpeciesBackTransformedDiagnostics) { + + + for (int isp = 0; isp < nSpeciesBackTransformedDiagnostics; ++isp) { auto np = tmp_particle_buffer[isp].GetRealData(DiagIdx::w).size(); if (np == 0) return; - auto const& wpc = tmp_particle_buffer[isp].GetRealData(DiagIdx::w); - auto const& xpc = tmp_particle_buffer[isp].GetRealData(DiagIdx::x); - auto const& ypc = tmp_particle_buffer[isp].GetRealData(DiagIdx::y); - auto const& zpc = tmp_particle_buffer[isp].GetRealData(DiagIdx::z); - auto const& uxpc = tmp_particle_buffer[isp].GetRealData(DiagIdx::ux); - auto const& uypc = tmp_particle_buffer[isp].GetRealData(DiagIdx::uy); - auto const& uzpc = tmp_particle_buffer[isp].GetRealData(DiagIdx::uz); - - particles_buffer_[isp].resize(np); - auto wpc_buff = particles_buffer_[isp].GetRealData(DiagIdx::w); - auto xpc_buff = particles_buffer_[isp].GetRealData(DiagIdx::x); - auto ypc_buff = particles_buffer_[isp].GetRealData(DiagIdx::y); - auto zpc_buff = particles_buffer_[isp].GetRealData(DiagIdx::z); - auto uxpc_buff = particles_buffer_[isp].GetRealData(DiagIdx::ux); - auto uypc_buff = particles_buffer_[isp].GetRealData(DiagIdx::uy); - auto uzpc_buff = particles_buffer_[isp].GetRealData(DiagIdx::uz); - - - int partcounter = 0; - for (int i = 0; i < np; ++i) - { - if( xpc[i] >= (diag_domain_lab_.lo(0)-dx_) && - xpc[i] <= (diag_domain_lab_.hi(0)+dx_) ) { - #if (AMREX_SPACEDIM == 3) - if( ypc[i] >= (diag_domain_lab_.lo(1)-dy_) && - ypc[i] <= (diag_domain_lab_.hi(1) + dy_)) - #endif - { - wpc_buff[partcounter] = wpc[i]; - xpc_buff[partcounter] = xpc[i]; - ypc_buff[partcounter] = ypc[i]; - zpc_buff[partcounter] = zpc[i]; - uxpc_buff[partcounter] = uxpc[i]; - uypc_buff[partcounter] = uypc[i]; - uzpc_buff[partcounter] = uzpc[i]; - ++partcounter; - } - } - } + Real const* const AMREX_RESTRICT wp_temp = + tmp_particle_buffer[isp].GetRealData(DiagIdx::w).data(); + Real const* const AMREX_RESTRICT x_temp = + tmp_particle_buffer[isp].GetRealData(DiagIdx::x).data(); + Real const* const AMREX_RESTRICT y_temp = + tmp_particle_buffer[isp].GetRealData(DiagIdx::y).data(); + Real const* const AMREX_RESTRICT z_temp = + tmp_particle_buffer[isp].GetRealData(DiagIdx::z).data(); + Real const* const AMREX_RESTRICT ux_temp = + tmp_particle_buffer[isp].GetRealData(DiagIdx::ux).data(); + Real const* const AMREX_RESTRICT uy_temp = + tmp_particle_buffer[isp].GetRealData(DiagIdx::uy).data(); + Real const* const AMREX_RESTRICT uz_temp = + tmp_particle_buffer[isp].GetRealData(DiagIdx::uz).data(); + + // temporary arrays to store copy_flag and copy_index + // for particles that cross the reduced domain for diagnostics. + amrex::Gpu::ManagedDeviceVector<int> FlagForPartCopy(np); + amrex::Gpu::ManagedDeviceVector<int> IndexForPartCopy(np); + + int* const AMREX_RESTRICT Flag = FlagForPartCopy.dataPtr(); + int* const AMREX_RESTRICT IndexLocation = IndexForPartCopy.dataPtr(); + + // Compute extent of the reduced domain +/- user-defined physical width + Real const xmin = m_diag_domain_lab_.lo(0)-m_particle_slice_dx_lab_; + Real const xmax = m_diag_domain_lab_.hi(0)+m_particle_slice_dx_lab_; +#if (AMREX_SPACEDIM == 3) + Real const ymin = m_diag_domain_lab_.lo(1)-m_particle_slice_dx_lab_; + Real const ymax = m_diag_domain_lab_.hi(1)+m_particle_slice_dx_lab_; +#endif - particles_buffer_[isp].resize(partcounter); + //Flag particles that need to be copied if they are + // within the reduced slice +/- user-defined physical width + amrex::ParallelFor(np, + [=] AMREX_GPU_DEVICE(int i) + { + Flag[i] = 0; + if ( x_temp[i] >= (xmin) && + x_temp[i] <= (xmax) ) { +#if (AMREX_SPACEDIM == 3) + if (y_temp[i] >= (ymin) && + y_temp[i] <= (ymax) ) +#endif + { + Flag[i] = 1; + } + } + }); + + // Call exclusive scan to obtain location indices using + // flag values. These location indices are used to copy data + // from src to dst when the copy-flag is set to 1. + amrex::Gpu::exclusive_scan(Flag,Flag+np,IndexLocation); + const int copy_size = IndexLocation[np-1] + Flag[np-1]; + const int init_size = m_particles_buffer_[isp].GetRealData(DiagIdx::w).size(); + const int total_reducedDiag_size = copy_size + init_size; + + // allocate array size for reduced diagnostic buffer array + m_particles_buffer_[isp].resize(total_reducedDiag_size); + + // Data pointers to particle attributes // + Real* const AMREX_RESTRICT wp_buff = + m_particles_buffer_[isp].GetRealData(DiagIdx::w).data(); + Real* const AMREX_RESTRICT x_buff = + m_particles_buffer_[isp].GetRealData(DiagIdx::x).data(); + Real* const AMREX_RESTRICT y_buff = + m_particles_buffer_[isp].GetRealData(DiagIdx::y).data(); + Real* const AMREX_RESTRICT z_buff = + m_particles_buffer_[isp].GetRealData(DiagIdx::z).data(); + Real* const AMREX_RESTRICT ux_buff = + m_particles_buffer_[isp].GetRealData(DiagIdx::ux).data(); + Real* const AMREX_RESTRICT uy_buff = + m_particles_buffer_[isp].GetRealData(DiagIdx::uy).data(); + Real* const AMREX_RESTRICT uz_buff = + m_particles_buffer_[isp].GetRealData(DiagIdx::uz).data(); + + // Selective copy of particle data from tmp array to reduced buffer + // array on the GPU using the flag value and index location. + amrex::ParallelFor(np, + [=] AMREX_GPU_DEVICE(int i) + { + if (Flag[i] == 1) + { + const int loc = IndexLocation[i] + init_size; + wp_buff[loc] = wp_temp[i]; + x_buff[loc] = x_temp[i]; + y_buff[loc] = y_temp[i]; + z_buff[loc] = z_temp[i]; + ux_buff[loc] = ux_temp[i]; + uy_buff[loc] = uy_temp[i]; + uz_buff[loc] = uz_temp[i]; + } + }); } } |