diff options
Diffstat (limited to 'Source')
| -rw-r--r-- | Source/Diagnostics/FieldIO.cpp | 35 | ||||
| -rw-r--r-- | Source/Diagnostics/ParticleIO.cpp | 12 | ||||
| -rw-r--r-- | Source/Diagnostics/SliceDiagnostic.cpp | 4 | ||||
| -rw-r--r-- | Source/Diagnostics/WarpXIO.cpp | 6 | ||||
| -rw-r--r-- | Source/Evolve/WarpXEvolveEM.cpp | 16 | ||||
| -rw-r--r-- | Source/FieldSolver/WarpX_FDTD.H | 32 | ||||
| -rw-r--r-- | Source/Make.WarpX | 2 | ||||
| -rwxr-xr-x | Source/Particles/Deposition/ChargeDeposition.H | 2 | ||||
| -rw-r--r-- | Source/Particles/Deposition/CurrentDeposition.H | 8 | ||||
| -rw-r--r-- | Source/Particles/MultiParticleContainer.cpp | 20 | ||||
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 110 | ||||
| -rw-r--r-- | Source/Particles/RigidInjectedParticleContainer.cpp | 43 | ||||
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.H | 18 | ||||
| -rw-r--r-- | Source/Particles/WarpXParticleContainer.cpp | 23 | ||||
| -rw-r--r-- | Source/WarpX.H | 1 | ||||
| -rw-r--r-- | Source/WarpX.cpp | 1 |
16 files changed, 138 insertions, 195 deletions
diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp index 5c523804e..540a968a2 100644 --- a/Source/Diagnostics/FieldIO.cpp +++ b/Source/Diagnostics/FieldIO.cpp @@ -13,7 +13,7 @@ using namespace amrex; namespace { // Return true if any element in elems is in vect, false otherwise - static bool is_in_vector(std::vector<std::string> vect, + static bool is_in_vector(std::vector<std::string> vect, std::vector<std::string> elems) { bool value = false; @@ -143,9 +143,16 @@ WriteOpenPMDFields( const std::string& filename, auto dataset = openPMD::Dataset(datatype, global_size); // Create new file and store the time/iteration info - auto series = openPMD::Series( filename, - openPMD::AccessType::CREATE, - MPI_COMM_WORLD ); + auto series = [filename](){ + if( ParallelDescriptor::NProcs() > 1 ) + return openPMD::Series( filename, + openPMD::AccessType::CREATE, + ParallelDescriptor::Communicator() ); + else + return openPMD::Series( filename, + openPMD::AccessType::CREATE ); + }(); + auto series_iteration = series.iterations[iteration]; series_iteration.setTime( time ); @@ -167,7 +174,7 @@ WriteOpenPMDFields( const std::string& filename, } } - // Setup the mesh accordingly + // Setup the mesh record accordingly auto mesh = series_iteration.meshes[field_name]; mesh.setDataOrder(openPMD::Mesh::DataOrder::F); // MultiFab: Fortran order mesh.setAxisLabels( axis_labels ); @@ -175,11 +182,11 @@ WriteOpenPMDFields( const std::string& filename, mesh.setGridGlobalOffset( global_offset ); setOpenPMDUnit( mesh, field_name ); - // Create a new mesh record, and store the associated metadata - auto mesh_record = mesh[comp_name]; - mesh_record.resetDataset( dataset ); + // Create a new mesh record component, and store the associated metadata + auto mesh_comp = mesh[comp_name]; + mesh_comp.resetDataset( dataset ); // Cell-centered data: position is at 0.5 of a cell size. - mesh_record.setPosition(std::vector<double>{AMREX_D_DECL(0.5, 0.5, 0.5)}); + mesh_comp.setPosition(std::vector<double>{AMREX_D_DECL(0.5, 0.5, 0.5)}); // Loop through the multifab, and store each box as a chunk, // in the openPMD file. @@ -195,8 +202,8 @@ WriteOpenPMDFields( const std::string& filename, // Write local data const double* local_data = fab.dataPtr(icomp); - mesh_record.storeChunk(openPMD::shareRaw(local_data), - chunk_offset, chunk_size); + mesh_comp.storeChunk(openPMD::shareRaw(local_data), + chunk_offset, chunk_size); } } // Flush data to disk after looping over all components @@ -318,14 +325,14 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames, + static_cast<int>(plot_finepatch)*6 + static_cast<int>(plot_crsepatch)*6 + static_cast<int>(costs[0] != nullptr); - + // Loop over levels of refinement for (int lev = 0; lev <= finest_level; ++lev) { // Allocate pointers to the `ncomp` fields that will be added mf_avg.push_back( MultiFab(grids[lev], dmap[lev], ncomp, ngrow)); - // For E, B and J, if at least one component is requested, + // For E, B and J, if at least one component is requested, // build cell-centered temporary MultiFab with 3 comps MultiFab mf_tmp_E, mf_tmp_B, mf_tmp_J; // Build mf_tmp_E is at least one component of E is requested @@ -420,7 +427,7 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames, {Bfield_aux[lev][0].get(), Bfield_aux[lev][1].get(), Bfield_aux[lev][2].get()}, - WarpX::CellSize(lev) ); + WarpX::CellSize(lev) ); } else if (fieldname == "divE"){ if (do_nodal) amrex::Abort("TODO: do_nodal && plot dive"); const BoxArray& ba = amrex::convert(boxArray(lev),IntVect::TheUnitVector()); diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index 743df4ce6..8bfa45a59 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -103,17 +103,6 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const real_names.push_back("theta"); #endif - if (WarpX::do_boosted_frame_diagnostic && pc->DoBoostedFrameDiags()) - { - real_names.push_back("xold"); - real_names.push_back("yold"); - real_names.push_back("zold"); - - real_names.push_back("uxold"); - real_names.push_back("uyold"); - real_names.push_back("uzold"); - } - if(pc->do_field_ionization){ int_names.push_back("ionization_level"); // int_flags specifies, for each integer attribs, whether it is @@ -122,6 +111,7 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const // when ionization is on. int_flags.resize(1, 1); } + // Convert momentum to SI pc->ConvertUnits(ConvertDirection::WarpX_to_SI); // real_names contains a list of all particle attributes. diff --git a/Source/Diagnostics/SliceDiagnostic.cpp b/Source/Diagnostics/SliceDiagnostic.cpp index 9f365b39d..0ec528e32 100644 --- a/Source/Diagnostics/SliceDiagnostic.cpp +++ b/Source/Diagnostics/SliceDiagnostic.cpp @@ -140,10 +140,10 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, MFIter mfi_dst(mfDst); for (MFIter mfi(mfSrc); mfi.isValid(); ++mfi) { - FArrayBox& Src_fabox = mfSrc[mfi]; + Array4<Real const> const& Src_fabox = mfSrc.const_array(mfi); const Box& Dst_bx = mfi_dst.validbox(); - FArrayBox& Dst_fabox = mfDst[mfi_dst]; + Array4<Real> const& Dst_fabox = mfDst.array(mfi_dst); int scomp = 0; int dcomp = 0; diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index 38399bf9e..9da0348d6 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -525,8 +525,10 @@ WarpX::WritePlotFile () const #ifdef WARPX_USE_OPENPMD if (dump_openpmd){ // Write openPMD format: only for level 0 - std::string filename = amrex::Concatenate("diags/hdf5/data", istep[0]); - filename += ".h5"; + std::string filename = std::string("diags/"); + filename.append(openpmd_backend); + filename += amrex::Concatenate("/data", istep[0]); + filename += std::string(".") + openpmd_backend; WriteOpenPMDFields( filename, varnames, *output_mf[0], output_geom[0], istep[0], t_new[0] ); } diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 49d855d0e..75643d748 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -140,6 +140,14 @@ WarpX::EvolveEM (int numsteps) bool do_insitu = ((step+1) >= insitu_start) && (insitu_int > 0) && ((step+1) % insitu_int == 0); + if (do_boosted_frame_diagnostic) { + std::unique_ptr<MultiFab> cell_centered_data = nullptr; + if (WarpX::do_boosted_frame_fields) { + cell_centered_data = GetCellCenteredData(); + } + myBFD->writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); + } + bool move_j = is_synchronized || to_make_plot || do_insitu; // 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. @@ -173,14 +181,6 @@ WarpX::EvolveEM (int numsteps) t_new[i] = cur_time; } - if (do_boosted_frame_diagnostic) { - std::unique_ptr<MultiFab> cell_centered_data = nullptr; - if (WarpX::do_boosted_frame_fields) { - cell_centered_data = GetCellCenteredData(); - } - myBFD->writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); - } - // slice gen // if (to_make_plot || do_insitu || to_make_slice_plot) { diff --git a/Source/FieldSolver/WarpX_FDTD.H b/Source/FieldSolver/WarpX_FDTD.H index e2d40f4c2..8f91036cc 100644 --- a/Source/FieldSolver/WarpX_FDTD.H +++ b/Source/FieldSolver/WarpX_FDTD.H @@ -13,7 +13,7 @@ void warpx_push_bx_yee(int i, int j, int k, Array4<Real> const& Bx, #if defined WARPX_DIM_3D Bx(i,j,k) += - dtsdy * (Ez(i,j+1,k ) - Ez(i,j,k)) + dtsdz * (Ey(i,j ,k+1) - Ey(i,j,k)); -#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ) +#elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) // Note that the 2D Cartesian and RZ versions are the same Bx(i,j,0) += + dtsdz * (Ey(i,j+1,0) - Ey(i,j,0)); #endif @@ -27,7 +27,7 @@ void warpx_push_by_yee(int i, int j, int k, Array4<Real> const& By, #if defined WARPX_DIM_3D By(i,j,k) += + dtsdx * (Ez(i+1,j,k ) - Ez(i,j,k)) - dtsdz * (Ex(i ,j,k+1) - Ex(i,j,k)); -#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ) +#elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) // Note that the 2D Cartesian and RZ versions are the same By(i,j,0) += + dtsdx * (Ez(i+1,j ,0) - Ez(i,j,0)) - dtsdz * (Ex(i ,j+1,0) - Ex(i,j,0)); @@ -42,7 +42,7 @@ void warpx_push_bz_yee(int i, int j, int k, Array4<Real> const& Bz, #if defined WARPX_DIM_3D Bz(i,j,k) += - dtsdx * (Ey(i+1,j ,k) - Ey(i,j,k)) + dtsdy * (Ex(i ,j+1,k) - Ex(i,j,k)); -#elif defined WARPX_DIM_2D +#elif defined WARPX_DIM_XZ Bz(i,j,0) += - dtsdx * (Ey(i+1,j,0) - Ey(i,j,0)); #elif defined WARPX_DIM_RZ const Real ru = 1. + 0.5/(rmin*dxinv + i + 0.5); @@ -60,7 +60,7 @@ void warpx_push_ex_yee(int i, int j, int k, Array4<Real> const& Ex, Ex(i,j,k) += + dtsdy_c2 * (Bz(i,j,k) - Bz(i,j-1,k )) - dtsdz_c2 * (By(i,j,k) - By(i,j ,k-1)) - mu_c2_dt * Jx(i,j,k); -#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ) +#elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) // Note that the 2D Cartesian and RZ versions are the same Ex(i,j,0) += - dtsdz_c2 * (By(i,j,0) - By(i,j-1,0)) - mu_c2_dt * Jx(i,j,0); @@ -76,7 +76,7 @@ void warpx_push_ey_yee(int i, int j, int k, Array4<Real> const& Ey, Ey(i,j,k) += - dtsdx_c2 * (Bz(i,j,k) - Bz(i-1,j,k)) + dtsdz_c2 * (Bx(i,j,k) - Bx(i,j,k-1)) - mu_c2_dt * Jy(i,j,k); -#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ) +#elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) // 2D Cartesian and RZ are the same, except that the axis is skipped with RZ #ifdef WARPX_DIM_RZ if (i != 0 || rmin != 0.) @@ -98,7 +98,7 @@ void warpx_push_ez_yee(int i, int j, int k, Array4<Real> const& Ez, Ez(i,j,k) += + dtsdx_c2 * (By(i,j,k) - By(i-1,j ,k)) - dtsdy_c2 * (Bx(i,j,k) - Bx(i ,j-1,k)) - mu_c2_dt * Jz(i,j,k); -#elif defined WARPX_DIM_2D +#elif defined WARPX_DIM_XZ Ez(i,j,0) += + dtsdx_c2 * (By(i,j,0) - By(i-1,j,0)) - mu_c2_dt * Jz(i,j,0); #elif defined WARPX_DIM_RZ @@ -120,7 +120,7 @@ void warpx_push_ex_f_yee(int j, int k, int l, Array4<Real> const& Ex, { #if defined WARPX_DIM_3D Ex(j,k,l) += + dtsdx_c2 * (F(j+1,k,l) - F(j,k,l)); -#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ) +#elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) Ex(j,k,0) += + dtsdx_c2 * (F(j+1,k,0) - F(j ,k,0)); #endif } @@ -140,7 +140,7 @@ void warpx_push_ez_f_yee(int j, int k, int l, Array4<Real> const& Ez, { #if defined WARPX_DIM_3D Ez(j,k,l) += + dtsdz_c2 * (F(j,k,l+1) - F(j,k,l)); -#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ) +#elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) Ez(j,k,0) += + dtsdz_c2 * (F(j,k+1,0) - F(j,k,0)); #endif } @@ -183,7 +183,7 @@ static void warpx_calculate_ckc_coefficients(Real dtsdx, Real dtsdy, Real dtsdz, gammax *= dtsdx; gammay *= dtsdy; gammaz *= dtsdz; -#elif defined WARPX_DIM_2D +#elif defined WARPX_DIM_XZ const Real delta = std::max(dtsdx,dtsdz); const Real rx = (dtsdx/delta)*(dtsdx/delta); const Real rz = (dtsdz/delta)*(dtsdz/delta); @@ -225,7 +225,7 @@ void warpx_push_bx_ckc(int j, int k, int l, Array4<Real> const& Bx, + Ey(j-1,k+1,l+1) - Ey(j-1,k+1,l ) + Ey(j+1,k-1,l+1) - Ey(j+1,k-1,l ) + Ey(j-1,k-1,l+1) - Ey(j-1,k-1,l )); -#elif defined WARPX_DIM_2D +#elif defined WARPX_DIM_XZ Bx(j,k,0) += + alphaz * (Ey(j ,k+1,0) - Ey(j, k,0)) + betazx * (Ey(j+1,k+1,0) - Ey(j+1,k,0) + Ey(j-1,k+1,0) - Ey(j-1,k,0)); @@ -258,7 +258,7 @@ void warpx_push_by_ckc(int j, int k, int l, Array4<Real> const& By, + Ex(j-1,k+1,l+1) - Ex(j-1,k+1,l ) + Ex(j+1,k-1,l+1) - Ex(j+1,k-1,l ) + Ex(j-1,k-1,l+1) - Ex(j-1,k-1,l )); -#elif defined WARPX_DIM_2D +#elif defined WARPX_DIM_XZ By(j,k,0) += + alphax * (Ez(j+1,k ,0) - Ez(j,k ,0)) + betaxz * (Ez(j+1,k+1,0) - Ez(j,k+1,0) + Ez(j+1,k-1,0) - Ez(j,k-1,0)) @@ -294,7 +294,7 @@ void warpx_push_bz_ckc(int j, int k, int l, Array4<Real> const& Bz, + Ex(j-1,k+1,l+1) - Ex(j-1,k ,l+1) + Ex(j+1,k+1,l-1) - Ex(j+1,k ,l-1) + Ex(j-1,k+1,l-1) - Ex(j-1,k ,l-1)); -#elif defined WARPX_DIM_2D +#elif defined WARPX_DIM_XZ Bz(j,k,0) += - alphax * (Ey(j+1,k ,0) - Ey(j,k ,0)) - betaxz * (Ey(j+1,k+1,0) - Ey(j,k+1,0) + Ey(j+1,k-1,0) - Ey(j,k-1,0)); @@ -318,7 +318,7 @@ void warpx_push_ex_f_ckc(int j, int k, int l, Array4<Real> const& Ex, + F(j+1,k-1,l+1) - F(j,k-1,l+1) + F(j+1,k+1,l-1) - F(j,k+1,l-1) + F(j+1,k-1,l-1) - F(j,k-1,l-1)); -#elif defined WARPX_DIM_2D +#elif defined WARPX_DIM_XZ Ex(j,k,0) += + alphax * (F(j+1,k ,0) - F(j,k ,0)) + betaxz * (F(j+1,k+1,0) - F(j,k+1,0) + F(j+1,k-1,0) - F(j,k-1,0)); @@ -362,7 +362,7 @@ void warpx_push_ez_f_ckc(int j, int k, int l, Array4<Real> const& Ez, + F(j-1,k+1,l+1) - F(j-1,k+1,l) + F(j+1,k-1,l+1) - F(j+1,k-1,l) + F(j-1,k-1,l+1) - F(j-1,k-1,l)); -#elif defined WARPX_DIM_2D +#elif defined WARPX_DIM_XZ Ez(j,k,0) += + alphaz * (F(j ,k+1,0) - F(j ,k,0)) + betazx * (F(j+1,k+1,0) - F(j+1,k,0) + F(j-1,k+1,0) - F(j-1,k,0)); @@ -382,7 +382,7 @@ void warpx_computedivb(int i, int j, int k, int dcomp, Array4<Real> const& divB, divB(i,j,k,dcomp) = (Bx(i+1,j ,k ) - Bx(i,j,k))*dxinv + (By(i ,j+1,k ) - By(i,j,k))*dyinv + (Bz(i ,j ,k+1) - Bz(i,j,k))*dzinv; -#elif defined WARPX_DIM_2D +#elif defined WARPX_DIM_XZ divB(i,j,0,dcomp) = (Bx(i+1,j ,0) - Bx(i,j,0))*dxinv + (Bz(i ,j+1,0) - Bz(i,j,0))*dzinv; #elif defined WARPX_DIM_RZ @@ -406,7 +406,7 @@ void warpx_computedive(int i, int j, int k, int dcomp, Array4<Real> const& divE, divE(i,j,k,dcomp) = (Ex(i,j,k) - Ex(i-1,j,k))*dxinv + (Ey(i,j,k) - Ey(i,j-1,k))*dyinv + (Ez(i,j,k) - Ez(i,j,k-1))*dzinv; -#elif defined WARPX_DIM_2D +#elif defined WARPX_DIM_XZ divE(i,j,0,dcomp) = (Ex(i,j,0) - Ex(i-1,j,0))*dxinv + (Ez(i,j,0) - Ez(i,j-1,0))*dzinv; #elif defined WARPX_DIM_RZ diff --git a/Source/Make.WarpX b/Source/Make.WarpX index 00d5d7bf5..ea06f8e06 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -51,7 +51,7 @@ else ifeq ($(USE_RZ),TRUE) DEFINES += -DWARPX_DIM_RZ else - DEFINES += -DWARPX_DIM_2D + DEFINES += -DWARPX_DIM_XZ endif endif diff --git a/Source/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H index 4253f5e76..f02eb1033 100755 --- a/Source/Particles/Deposition/ChargeDeposition.H +++ b/Source/Particles/Deposition/ChargeDeposition.H @@ -82,7 +82,7 @@ void doChargeDepositionShapeN(const amrex::Real * const xp, const int k = compute_shape_factor<depos_order>(sz, z); // Deposit charge into rho_arr -#if (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ) +#if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order; ix++){ amrex::Gpu::Atomic::Add( diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index e8e7fab0a..b879970ef 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -133,7 +133,7 @@ void doDepositionShapeN(const amrex::Real * const xp, const int l0 = compute_shape_factor<depos_order>(sz0, zmid-stagger_shift); // Deposit current into jx_arr, jy_arr and jz_arr -#if (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ) +#if (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order; ix++){ amrex::Gpu::Atomic::Add( @@ -224,7 +224,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, const amrex::Real invdtdx = 1.0/(dt*dx[1]*dx[2]); const amrex::Real invdtdy = 1.0/(dt*dx[0]*dx[2]); const amrex::Real invdtdz = 1.0/(dt*dx[0]*dx[1]); -#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ) +#elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) const amrex::Real invdtdx = 1.0/(dt*dx[2]); const amrex::Real invdtdz = 1.0/(dt*dx[0]); const amrex::Real invvol = 1.0/(dx[0]*dx[2]); @@ -282,7 +282,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, sintheta = 0.; } const amrex::Real vy = (-uxp[ip]*sintheta + uyp[ip]*costheta)*gaminv; -#elif (defined WARPX_DIM_2D) +#elif (defined WARPX_DIM_XZ) const amrex::Real vy = uyp[ip]*gaminv; #endif @@ -357,7 +357,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, } } -#elif (defined WARPX_DIM_2D) || (defined WARPX_DIM_RZ) +#elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) for (int k=dkl; k<=depos_order+2-dku; k++) { amrex::Real sdxi = 0.; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 9e2cd097f..143f4b7f9 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -42,26 +42,6 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) nspecies_boosted_frame_diags += 1; } } - - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - for (int i = 0; i < nspecies_boosted_frame_diags; ++i) - { - int is = map_species_boosted_frame_diags[i]; - allcontainers[is]->AddRealComp("xold"); - allcontainers[is]->AddRealComp("yold"); - allcontainers[is]->AddRealComp("zold"); - allcontainers[is]->AddRealComp("uxold"); - allcontainers[is]->AddRealComp("uyold"); - allcontainers[is]->AddRealComp("uzold"); - } - pc_tmp->AddRealComp("xold"); - pc_tmp->AddRealComp("yold"); - pc_tmp->AddRealComp("zold"); - pc_tmp->AddRealComp("uxold"); - pc_tmp->AddRealComp("uyold"); - pc_tmp->AddRealComp("uzold"); - } } void diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 9750e7e59..015482e3f 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -212,17 +212,8 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, { // need to create old values auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - particle_tile.push_back_real(particle_comps["xold"], x); - particle_tile.push_back_real(particle_comps["yold"], y); - particle_tile.push_back_real(particle_comps["zold"], z); - - particle_tile.push_back_real(particle_comps["uxold"], u[0]); - particle_tile.push_back_real(particle_comps["uyold"], u[1]); - particle_tile.push_back_real(particle_comps["uzold"], u[2]); - } } + // add particle AddOneParticle(0, 0, 0, x, y, z, attribs); } @@ -301,11 +292,13 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) #ifdef _OPENMP // First touch all tiles in the map in serial for (MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) { - const int grid_id = mfi.index(); - const int tile_id = mfi.LocalTileIndex(); - GetParticles(lev)[std::make_pair(grid_id, tile_id)]; + auto index = std::make_pair(mfi.index(), mfi.LocalTileIndex()); + GetParticles(lev)[index]; + tmp_particle_data.resize(finestLevel()+1); + // Create map entry if not there + tmp_particle_data[lev][index]; if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ) { - DefineAndReturnParticleTile(lev, grid_id, tile_id); + DefineAndReturnParticleTile(lev, mfi.index(), mfi.LocalTileIndex()); } } #endif @@ -405,14 +398,16 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) // then how many new particles will be injected is not that simple // We have to shift fine_injection_box because overlap_box has been shifted. Box fine_overlap_box = overlap_box & amrex::shift(fine_injection_box,shifted); - max_new_particles += fine_overlap_box.numPts() * num_ppc - * (AMREX_D_TERM(rrfac,*rrfac,*rrfac)-1); - for (int icell = 0, ncells = overlap_box.numPts(); icell < ncells; ++icell) { - IntVect iv = overlap_box.atOffset(icell); - int r = (fine_overlap_box.contains(iv)) ? AMREX_D_TERM(rrfac,*rrfac,*rrfac) : 1; - for (int ipart = 0; ipart < r; ++ipart) { - cellid_v.push_back(icell); - cellid_v.push_back(ipart); + if (fine_overlap_box.ok()) { + max_new_particles += fine_overlap_box.numPts() * num_ppc + * (AMREX_D_TERM(rrfac,*rrfac,*rrfac)-1); + for (int icell = 0, ncells = overlap_box.numPts(); icell < ncells; ++icell) { + IntVect iv = overlap_box.atOffset(icell); + int r = (fine_overlap_box.contains(iv)) ? AMREX_D_TERM(rrfac,*rrfac,*rrfac) : 1; + for (int ipart = 0; ipart < r; ++ipart) { + cellid_v.push_back(icell); + cellid_v.push_back(ipart); + } } } } @@ -430,12 +425,10 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) const int cpuid = ParallelDescriptor::MyProc(); auto& particle_tile = GetParticles(lev)[std::make_pair(grid_id,tile_id)]; - bool do_boosted = false; + if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ) { DefineAndReturnParticleTile(lev, grid_id, tile_id); } - do_boosted = WarpX::do_boosted_frame_diagnostic - && do_boosted_frame_diags; auto old_size = particle_tile.GetArrayOfStructs().size(); auto new_size = old_size + max_new_particles; @@ -447,15 +440,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) for (int ia = 0; ia < PIdx::nattribs; ++ia) { pa[ia] = soa.GetRealData(ia).data() + old_size; } - GpuArray<Real*,6> pb; - if (do_boosted) { - pb[0] = soa.GetRealData(particle_comps[ "xold"]).data() + old_size; - pb[1] = soa.GetRealData(particle_comps[ "yold"]).data() + old_size; - pb[2] = soa.GetRealData(particle_comps[ "zold"]).data() + old_size; - pb[3] = soa.GetRealData(particle_comps["uxold"]).data() + old_size; - pb[4] = soa.GetRealData(particle_comps["uyold"]).data() + old_size; - pb[5] = soa.GetRealData(particle_comps["uzold"]).data() + old_size; - } + int* pi; if (do_field_ionization) { pi = soa.GetIntData(particle_icomps["ionization_level"]).data() + old_size; @@ -617,15 +602,6 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) pa[PIdx::uy][ip] = u.y; pa[PIdx::uz][ip] = u.z; - if (do_boosted) { - pb[0][ip] = x; - pb[1][ip] = y; - pb[2][ip] = z; - pb[3][ip] = u.x; - pb[4][ip] = u.y; - pb[5][ip] = u.z; - } - #if (AMREX_SPACEDIM == 3) p.pos(0) = x; p.pos(1) = y; @@ -976,6 +952,19 @@ PhysicalParticleContainer::Evolve (int lev, bool has_buffer = cEx || cjx; + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + { + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) + { + const auto np = pti.numParticles(); + const auto lev = pti.GetLevel(); + const auto index = pti.GetPairIndex(); + tmp_particle_data.resize(finestLevel()+1); + for (int i = 0; i < TmpIdx::nattribs; ++i) + tmp_particle_data[lev][index][i].resize(np); + } + } + #ifdef _OPENMP #pragma omp parallel #endif @@ -1688,22 +1677,21 @@ PhysicalParticleContainer::PushP (int lev, Real dt, void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, const Real* yp, const Real* zp) { - auto& attribs = pti.GetAttribs(); - Real* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); Real* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr(); Real* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr(); - Real* AMREX_RESTRICT xpold = pti.GetAttribs(particle_comps["xold"]).dataPtr(); - Real* AMREX_RESTRICT ypold = pti.GetAttribs(particle_comps["yold"]).dataPtr(); - Real* AMREX_RESTRICT zpold = pti.GetAttribs(particle_comps["zold"]).dataPtr(); - Real* AMREX_RESTRICT uxpold = pti.GetAttribs(particle_comps["uxold"]).dataPtr(); - Real* AMREX_RESTRICT uypold = pti.GetAttribs(particle_comps["uyold"]).dataPtr(); - Real* AMREX_RESTRICT uzpold = pti.GetAttribs(particle_comps["uzold"]).dataPtr(); - - const long np = pti.numParticles(); - + const auto np = pti.numParticles(); + const auto lev = pti.GetLevel(); + const auto index = pti.GetPairIndex(); + Real* AMREX_RESTRICT xpold = tmp_particle_data[lev][index][TmpIdx::xold ].dataPtr(); + Real* AMREX_RESTRICT ypold = tmp_particle_data[lev][index][TmpIdx::yold ].dataPtr(); + Real* AMREX_RESTRICT zpold = tmp_particle_data[lev][index][TmpIdx::zold ].dataPtr(); + Real* AMREX_RESTRICT uxpold = tmp_particle_data[lev][index][TmpIdx::uxold].dataPtr(); + Real* AMREX_RESTRICT uypold = tmp_particle_data[lev][index][TmpIdx::uyold].dataPtr(); + Real* AMREX_RESTRICT uzpold = tmp_particle_data[lev][index][TmpIdx::uzold].dataPtr(); + ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { xpold[i]=xp[i]; @@ -1788,15 +1776,15 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real auto& uyp_new = attribs[PIdx::uy ]; auto& uzp_new = attribs[PIdx::uz ]; - auto& xp_old = pti.GetAttribs(particle_comps["xold"]); - auto& yp_old = pti.GetAttribs(particle_comps["yold"]); - auto& zp_old = pti.GetAttribs(particle_comps["zold"]); - auto& uxp_old = pti.GetAttribs(particle_comps["uxold"]); - auto& uyp_old = pti.GetAttribs(particle_comps["uyold"]); - auto& uzp_old = pti.GetAttribs(particle_comps["uzold"]); + auto& xp_old = tmp_particle_data[lev][index][TmpIdx::xold]; + auto& yp_old = tmp_particle_data[lev][index][TmpIdx::yold]; + auto& zp_old = tmp_particle_data[lev][index][TmpIdx::zold]; + auto& uxp_old = tmp_particle_data[lev][index][TmpIdx::uxold]; + auto& uyp_old = tmp_particle_data[lev][index][TmpIdx::uyold]; + auto& uzp_old = tmp_particle_data[lev][index][TmpIdx::uzold]; const long np = pti.numParticles(); - + Real uzfrm = -WarpX::gamma_boost*WarpX::beta_boost*PhysConst::c; Real inv_c2 = 1.0/PhysConst::c/PhysConst::c; diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 36cb9d224..0c94d1e33 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -239,15 +239,13 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, Real* const AMREX_RESTRICT Bzp = attribs[PIdx::Bz].dataPtr(); if (!done_injecting_lev) { - if (!(WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)) { - // If the old values are not already saved, create copies here. - xp_save = xp; - yp_save = yp; - zp_save = zp; - uxp_save = uxp; - uyp_save = uyp; - uzp_save = uzp; - } + // If the old values are not already saved, create copies here. + xp_save = xp; + yp_save = yp; + zp_save = zp; + uxp_save = uxp; + uyp_save = uyp; + uzp_save = uzp; // Scale the fields of particles about to cross the injection plane. // This only approximates what should be happening. The particles @@ -275,27 +273,12 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, if (!done_injecting_lev) { - Real* AMREX_RESTRICT x_save; - Real* AMREX_RESTRICT y_save; - Real* AMREX_RESTRICT z_save; - Real* AMREX_RESTRICT ux_save; - Real* AMREX_RESTRICT uy_save; - Real* AMREX_RESTRICT uz_save; - if (!(WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)) { - x_save = xp_save.dataPtr(); - y_save = yp_save.dataPtr(); - z_save = zp_save.dataPtr(); - ux_save = uxp_save.dataPtr(); - uy_save = uyp_save.dataPtr(); - uz_save = uzp_save.dataPtr(); - } else { - x_save = pti.GetAttribs(particle_comps["xold"]).dataPtr(); - y_save = pti.GetAttribs(particle_comps["yold"]).dataPtr(); - z_save = pti.GetAttribs(particle_comps["zold"]).dataPtr(); - ux_save = pti.GetAttribs(particle_comps["uxold"]).dataPtr(); - uy_save = pti.GetAttribs(particle_comps["uyold"]).dataPtr(); - uz_save = pti.GetAttribs(particle_comps["uzold"]).dataPtr(); - } + Real* AMREX_RESTRICT x_save = xp_save.dataPtr(); + Real* AMREX_RESTRICT y_save = yp_save.dataPtr(); + Real* AMREX_RESTRICT z_save = zp_save.dataPtr(); + Real* AMREX_RESTRICT ux_save = uxp_save.dataPtr(); + Real* AMREX_RESTRICT uy_save = uyp_save.dataPtr(); + Real* AMREX_RESTRICT uz_save = uzp_save.dataPtr(); // Undo the push for particles not injected yet. // The zp are advanced a fixed amount. diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index cc553f9a7..7393f7301 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -29,6 +29,15 @@ struct DiagIdx }; }; +struct TmpIdx +{ + enum { + xold = 0, + yold, zold, uxold, uyold, uzold, + nattribs + }; +}; + namespace ParticleStringNames { const std::map<std::string, int> to_index = { @@ -308,7 +317,10 @@ protected: amrex::Vector<amrex::FArrayBox> local_jy; amrex::Vector<amrex::FArrayBox> local_jz; - amrex::Vector<amrex::Cuda::ManagedDeviceVector<amrex::Real> > m_xp, m_yp, m_zp, m_giv; + using DataContainer = amrex::Gpu::ManagedDeviceVector<amrex::Real>; + using PairIndex = std::pair<int, int>; + + amrex::Vector<DataContainer> m_xp, m_yp, m_zp, m_giv; // Whether to dump particle quantities. // If true, particle position is always dumped. @@ -318,7 +330,9 @@ protected: amrex::Vector<int> plot_flags; // list of names of attributes to dump. amrex::Vector<std::string> plot_vars; - + + amrex::Vector<std::map<PairIndex, std::array<DataContainer, TmpIdx::nattribs> > > tmp_particle_data; + private: virtual void particlePostLocate(ParticleType& p, const amrex::ParticleLocData& pld, const int lev) override; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 646644f94..80d1caa0f 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -85,17 +85,6 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) particle_comps["theta"] = PIdx::theta; #endif - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - particle_comps["xold"] = PIdx::nattribs; - particle_comps["yold"] = PIdx::nattribs+1; - particle_comps["zold"] = PIdx::nattribs+2; - particle_comps["uxold"] = PIdx::nattribs+3; - particle_comps["uyold"] = PIdx::nattribs+4; - particle_comps["uzold"] = PIdx::nattribs+5; - - } - // Initialize temporary local arrays for charge/current deposition int num_threads = 1; #ifdef _OPENMP @@ -240,12 +229,6 @@ WarpXParticleContainer::AddNParticles (int lev, if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ){ auto& ptile = DefineAndReturnParticleTile(0, 0, 0); - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - ptile.push_back_real(particle_comps["xold"], x[i]); - ptile.push_back_real(particle_comps["yold"], y[i]); - ptile.push_back_real(particle_comps["zold"], z[i]); - } } particle_tile.push_back(p); @@ -260,12 +243,6 @@ WarpXParticleContainer::AddNParticles (int lev, if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ){ auto& ptile = DefineAndReturnParticleTile(0, 0, 0); - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - ptile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); - ptile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend); - ptile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); - } } for (int comp = PIdx::uz+1; comp < PIdx::nattribs; ++comp) diff --git a/Source/WarpX.H b/Source/WarpX.H index 0f78f22f3..f1a3b1479 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -540,6 +540,7 @@ private: int check_int = -1; int plot_int = -1; + std::string openpmd_backend {"h5"}; #ifdef WARPX_USE_OPENPMD bool dump_plotfiles = false; bool dump_openpmd = true; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index b3cb2fbc9..4723626ab 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -412,6 +412,7 @@ WarpX::ReadParameters () pp.query("dump_openpmd", dump_openpmd); + pp.query("openpmd_backend", openpmd_backend); pp.query("dump_plotfiles", dump_plotfiles); pp.query("plot_raw_fields", plot_raw_fields); pp.query("plot_raw_fields_guards", plot_raw_fields_guards); |