From d528e233a7710dc497174c825620f1de7323faa5 Mon Sep 17 00:00:00 2001 From: Dave Grote Date: Fri, 26 Apr 2019 21:30:35 -0700 Subject: Implemented dive cleaning for RZ geometry --- Source/Particles/WarpXParticleContainer.cpp | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 6a882be99..46f6fc1fc 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -645,6 +645,11 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, &ngRho, &ngRho, &ngRho, &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect, &WarpX::charge_deposition_algo); +#ifdef WARPX_RZ + warpx_charge_deposition_rz_volume_scaling( + data_ptr, &ngRho, rholen.getVect(), + &xyzmin[0], &dx[0]); +#endif BL_PROFILE_VAR_STOP(blp_pxr_chd); BL_PROFILE_VAR_START(blp_accumulate); @@ -703,6 +708,11 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, &ngRho, &ngRho, &ngRho, &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect, &WarpX::charge_deposition_algo); +#ifdef WARPX_RZ + warpx_charge_deposition_rz_volume_scaling( + data_ptr, &ngRho, rholen.getVect(), + &cxyzmin_tile[0], &cdx[0]); +#endif BL_PROFILE_VAR_STOP(blp_pxr_chd); BL_PROFILE_VAR_START(blp_accumulate); @@ -863,6 +873,12 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) &dx[0], &dx[1], &dx[2], &nx, &ny, &nz, &nxg, &nyg, &nzg, &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect, &WarpX::charge_deposition_algo); +#ifdef WARPX_RZ + long ngRho = WarpX::nox; + warpx_charge_deposition_rz_volume_scaling( + data_ptr, &ngRho, rholen.getVect(), + &xyzmin[0], &dx[0]); +#endif #ifdef _OPENMP rhofab.atomicAdd(local_rho); -- cgit v1.2.3 From 9f99d4f819c7760d174c8888134496bdee5b03f0 Mon Sep 17 00:00:00 2001 From: Revathi Jambunathan Date: Mon, 6 May 2019 20:54:02 -0400 Subject: printoutstatement to debug the swap in rho_old and rho_new in the gpu vs cpu --- Source/Evolve/WarpXEvolveEM.cpp | 53 +++++++++++++++++++++++++++++ Source/FieldSolver/WarpXFFT.cpp | 14 ++++++++ Source/Initialization/WarpXInitData.cpp | 9 +++++ Source/Particles/WarpXParticleContainer.cpp | 2 +- 4 files changed, 77 insertions(+), 1 deletion(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 60b0b5fa3..025ff993a 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -80,12 +80,32 @@ WarpX::EvolveEM (int numsteps) *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2]); } is_synchronized = false; + + for (MFIter mfi(*rho_fp[0]); mfi.isValid(); ++mfi) + { + MultiFab &mf = *rho_fp[0]; + Box realspace_bx = mf[mfi].box(); // Copy the box + Array4 mf_arr = mf[mfi].array(); + amrex::Print() << " at if synchronized rho " << mf_arr(0,0,0,0) ; + amrex::Print() << " new rho " << mf_arr(0,0,0,1) << "\n"; + } + + } else { // Beyond one step, we have E^{n} and B^{n}. // Particles have p^{n-1/2} and x^{n}. FillBoundaryE(); FillBoundaryB(); UpdateAuxilaryData(); + + for (MFIter mfi(*rho_fp[0]); mfi.isValid(); ++mfi) + { + MultiFab &mf = *rho_fp[0]; + Box realspace_bx = mf[mfi].box(); // Copy the box + Array4 mf_arr = mf[mfi].array(); + amrex::Print() << " at if synchronized rho " << mf_arr(0,0,0,0) ; + amrex::Print() << " new rho " << mf_arr(0,0,0,1) << "\n"; + } } if (do_subcycling == 0 || finest_level == 0) { @@ -258,15 +278,48 @@ WarpX::OneStep_nosub (Real cur_time) if (warpx_py_particlescraper) warpx_py_particlescraper(); if (warpx_py_beforedeposition) warpx_py_beforedeposition(); #endif + + for (MFIter mfi(*rho_fp[0]); mfi.isValid(); ++mfi) + { + MultiFab &mf = *rho_fp[0]; + Box realspace_bx = mf[mfi].box(); // Copy the box + Array4 mf_arr = mf[mfi].array(); + amrex::Print() << " at before pushf particles and depose rho " << mf_arr(0,0,0,0) ; + amrex::Print() << " new rho " << mf_arr(0,0,0,1) << "\n"; + } PushParticlesandDepose(cur_time); + for (MFIter mfi(*rho_fp[0]); mfi.isValid(); ++mfi) + { + MultiFab &mf = *rho_fp[0]; + Box realspace_bx = mf[mfi].box(); // Copy the box + Array4 mf_arr = mf[mfi].array(); + amrex::Print() << " at after pushf particles and depose rho " << mf_arr(0,0,0,0) ; + amrex::Print() << " new rho " << mf_arr(0,0,0,1) << "\n"; + } #ifdef WARPX_USE_PY if (warpx_py_afterdeposition) warpx_py_afterdeposition(); #endif SyncCurrent(); + for (MFIter mfi(*rho_fp[0]); mfi.isValid(); ++mfi) + { + MultiFab &mf = *rho_fp[0]; + Box realspace_bx = mf[mfi].box(); // Copy the box + Array4 mf_arr = mf[mfi].array(); + amrex::Print() << " before sync rho " << mf_arr(0,0,0,0) ; + amrex::Print() << " new rho " << mf_arr(0,0,0,1) << "\n"; + } SyncRho(rho_fp, rho_cp); + for (MFIter mfi(*rho_fp[0]); mfi.isValid(); ++mfi) + { + MultiFab &mf = *rho_fp[0]; + Box realspace_bx = mf[mfi].box(); // Copy the box + Array4 mf_arr = mf[mfi].array(); + amrex::Print() << " after sync rho " << mf_arr(0,0,0,0) ; + amrex::Print() << " new rho " << mf_arr(0,0,0,1) << "\n"; + } // Push E and B from {n} to {n+1} // (And update guard cells immediately afterwards) #ifdef WARPX_USE_PSATD diff --git a/Source/FieldSolver/WarpXFFT.cpp b/Source/FieldSolver/WarpXFFT.cpp index e8b1442bb..ffeb237c0 100644 --- a/Source/FieldSolver/WarpXFFT.cpp +++ b/Source/FieldSolver/WarpXFFT.cpp @@ -405,6 +405,20 @@ WarpX::PushPSATD (int lev, amrex::Real /* dt */) rho_fp_fft[lev]->ParallelCopy(*rho_fp[lev], 0, 0, 2, 0, 0, period_fp); BL_PROFILE_VAR_STOP(blp_copy); + + for (MFIter mfi(*rho_fp_fft[lev]); mfi.isValid(); ++mfi) + { + MultiFab &mf = *rho_fp_fft[lev]; + Box realspace_bx = mf[mfi].box(); // Copy the box + Array4 mf_arr = mf[mfi].array(); + amrex::Print() << " at begin Push PSATD rho " << mf_arr(0,0,0,0) ; + amrex::Print() << " new rho " << mf_arr(0,0,0,1) << "\n"; + MultiFab &mf_orig = *rho_fp[lev]; + Array4 mforig_arr = mf_orig[mfi].array(); + amrex::Print() << " at begin Push PSATD rho " << mforig_arr(0,0,0,0) ; + amrex::Print() << " new rho " << mforig_arr(0,0,0,1) << "\n"; + } + BL_PROFILE_VAR_START(blp_push_eb); if (fft_hybrid_mpi_decomposition){ if (Efield_fp_fft[lev][0]->local_size() == 1) diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 23637ec97..589ee168c 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -342,6 +342,15 @@ WarpX::InitLevelDataFFT (int lev, Real time) current_cp_fft[lev][2]->setVal(0.0); rho_cp_fft[lev]->setVal(0.0); } + + for (MFIter mfi(*rho_fp_fft[lev]); mfi.isValid(); ++mfi) + { + MultiFab &mf = *rho_fp_fft[lev]; + Box realspace_bx = mf[mfi].box(); // Copy the box + Array4 mf_arr = mf[mfi].array(); + amrex::Print() << " at initialization rho " << mf_arr(0,0,0,0) ; + amrex::Print() << " new rho " << mf_arr(0,0,0,1) << "\n"; + } } #endif diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 2edd3c636..695faaa62 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -608,7 +608,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, const std::array& xyzmin = xyzmin_tile; #ifdef AMREX_USE_GPU - data_ptr = (*rhomf)[pti].dataPtr(); + data_ptr = (*rhomf)[pti].dataPtr(icomp); auto rholen = (*rhomf)[pti].length(); #else tile_box.grow(ngRho); -- cgit v1.2.3 From baed8d31df7b9dd0013aaa1c1382801eba918bbd Mon Sep 17 00:00:00 2001 From: Revathi Jambunathan Date: Mon, 6 May 2019 20:58:57 -0400 Subject: Passing icomp as a variable to the dataPtr() function in the charge deposit subroutine -- ensuring that the rho_old values are not over-written by rho_new --- Source/Particles/WarpXParticleContainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 2edd3c636..695faaa62 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -608,7 +608,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, const std::array& xyzmin = xyzmin_tile; #ifdef AMREX_USE_GPU - data_ptr = (*rhomf)[pti].dataPtr(); + data_ptr = (*rhomf)[pti].dataPtr(icomp); auto rholen = (*rhomf)[pti].length(); #else tile_box.grow(ngRho); -- cgit v1.2.3 From eb76cadbc5e76111679f6189e6728e462af41f66 Mon Sep 17 00:00:00 2001 From: Revathi Jambunathan Date: Wed, 8 May 2019 16:21:18 -0400 Subject: temporary print files for debugging --- GNUmakefile | 2 +- Source/Evolve/WarpXEvolveEM.cpp | 58 +++--------------- .../SpectralSolver/SpectralFieldData.cpp | 52 ++++++---------- Source/FieldSolver/WarpXFFT.cpp | 70 +++++++--------------- Source/Initialization/WarpXInitData.cpp | 20 ++++--- Source/Parallelization/WarpXComm.cpp | 7 ++- Source/Particles/WarpXParticleContainer.cpp | 2 +- Source/WarpX.cpp | 2 + Source/main.cpp | 3 + 9 files changed, 72 insertions(+), 144 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/GNUmakefile b/GNUmakefile index c662281b3..e3ed55d80 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -29,7 +29,7 @@ USE_ASCENT_INSITU = FALSE WarpxBinDir = Bin -USE_PSATD = TRUE +USE_PSATD = FALSE USE_RZ = FALSE DO_ELECTROSTATIC = FALSE diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 025ff993a..1c5dc0104 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -69,27 +69,20 @@ WarpX::EvolveEM (int numsteps) // At the beginning, we have B^{n} and E^{n}. // Particles have p^{n} and x^{n}. // is_synchronized is true. + amrex::Print() << " in evolve before fill boundary \n"; if (is_synchronized) { FillBoundaryE(); FillBoundaryB(); UpdateAuxilaryData(); // on first step, push p by -0.5*dt + amrex::Print() << " in evolve before pushP \n"; for (int lev = 0; lev <= finest_level; ++lev) { mypc->PushP(lev, -0.5*dt[lev], *Efield_aux[lev][0],*Efield_aux[lev][1],*Efield_aux[lev][2], *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2]); } is_synchronized = false; - - for (MFIter mfi(*rho_fp[0]); mfi.isValid(); ++mfi) - { - MultiFab &mf = *rho_fp[0]; - Box realspace_bx = mf[mfi].box(); // Copy the box - Array4 mf_arr = mf[mfi].array(); - amrex::Print() << " at if synchronized rho " << mf_arr(0,0,0,0) ; - amrex::Print() << " new rho " << mf_arr(0,0,0,1) << "\n"; - } - + amrex::Print() << " in evolve after pushP \n"; } else { // Beyond one step, we have E^{n} and B^{n}. @@ -98,15 +91,8 @@ WarpX::EvolveEM (int numsteps) FillBoundaryB(); UpdateAuxilaryData(); - for (MFIter mfi(*rho_fp[0]); mfi.isValid(); ++mfi) - { - MultiFab &mf = *rho_fp[0]; - Box realspace_bx = mf[mfi].box(); // Copy the box - Array4 mf_arr = mf[mfi].array(); - amrex::Print() << " at if synchronized rho " << mf_arr(0,0,0,0) ; - amrex::Print() << " new rho " << mf_arr(0,0,0,1) << "\n"; - } } + amrex::Print() << " in evolve after fill boundary \n"; if (do_subcycling == 0 || finest_level == 0) { OneStep_nosub(cur_time); @@ -116,6 +102,7 @@ WarpX::EvolveEM (int numsteps) amrex::Print() << "Error: do_subcycling = " << do_subcycling << std::endl; amrex::Abort("Unsupported do_subcycling type"); } + amrex::Print() << " in evolve after onestep no sub \n"; #ifdef WARPX_USE_PY if (warpx_py_beforeEsolve) warpx_py_beforeEsolve(); @@ -279,52 +266,23 @@ WarpX::OneStep_nosub (Real cur_time) if (warpx_py_beforedeposition) warpx_py_beforedeposition(); #endif - for (MFIter mfi(*rho_fp[0]); mfi.isValid(); ++mfi) - { - MultiFab &mf = *rho_fp[0]; - Box realspace_bx = mf[mfi].box(); // Copy the box - Array4 mf_arr = mf[mfi].array(); - amrex::Print() << " at before pushf particles and depose rho " << mf_arr(0,0,0,0) ; - amrex::Print() << " new rho " << mf_arr(0,0,0,1) << "\n"; - } PushParticlesandDepose(cur_time); - for (MFIter mfi(*rho_fp[0]); mfi.isValid(); ++mfi) - { - MultiFab &mf = *rho_fp[0]; - Box realspace_bx = mf[mfi].box(); // Copy the box - Array4 mf_arr = mf[mfi].array(); - amrex::Print() << " at after pushf particles and depose rho " << mf_arr(0,0,0,0) ; - amrex::Print() << " new rho " << mf_arr(0,0,0,1) << "\n"; - } + #ifdef WARPX_USE_PY if (warpx_py_afterdeposition) warpx_py_afterdeposition(); #endif SyncCurrent(); - for (MFIter mfi(*rho_fp[0]); mfi.isValid(); ++mfi) - { - MultiFab &mf = *rho_fp[0]; - Box realspace_bx = mf[mfi].box(); // Copy the box - Array4 mf_arr = mf[mfi].array(); - amrex::Print() << " before sync rho " << mf_arr(0,0,0,0) ; - amrex::Print() << " new rho " << mf_arr(0,0,0,1) << "\n"; - } SyncRho(rho_fp, rho_cp); - for (MFIter mfi(*rho_fp[0]); mfi.isValid(); ++mfi) - { - MultiFab &mf = *rho_fp[0]; - Box realspace_bx = mf[mfi].box(); // Copy the box - Array4 mf_arr = mf[mfi].array(); - amrex::Print() << " after sync rho " << mf_arr(0,0,0,0) ; - amrex::Print() << " new rho " << mf_arr(0,0,0,1) << "\n"; - } // Push E and B from {n} to {n+1} // (And update guard cells immediately afterwards) #ifdef WARPX_USE_PSATD PushPSATD(dt[0]); + amrex::Print() << " before fill bndry E \n"; FillBoundaryE(); + amrex::Print() << " before fill bndry B \n"; FillBoundaryB(); #else EvolveF(0.5*dt[0], DtType::FirstHalf); diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 7c2061f8d..5998bdd2b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -55,15 +55,13 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, #ifdef AMREX_USE_GPU // Add cuFFT-specific code // Creating 3D plan for real to complex -- double precision + cudaDeviceSynchronize(); cufftResult result; result = cufftPlan3d( &forward_plan[mfi], fft_size[2], fft_size[1],fft_size[0], CUFFT_D2Z); if ( result != CUFFT_SUCCESS ) { amrex::Print() << " cufftplan3d forward failed! \n"; } - if ( result == CUFFT_SUCCESS ) { - amrex::Print() << " created cufft forward plan\n"; - } // Add 2D cuFFT-spacific code for D2Z // Note that D2Z is inherently forward plan @@ -73,9 +71,7 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, if ( result != CUFFT_SUCCESS ) { amrex::Print() << " cufftplan3d backward failed! \n"; } - if ( result == CUFFT_SUCCESS ) { - amrex::Print() << " created cufft backward plan\n"; - } + cudaDeviceSynchronize(); #else // Create FFTW plans @@ -156,15 +152,9 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { tmp_arr(i,j,k) = mf_arr(i,j,k,i_comp); }); -#ifdef AMREX_USE_GPU - cudaDeviceSynchronize(); -#endif - amrex::Print() << " in forward trans icomp " << i_comp << " " << tmp_arr(0,0,0) << " mf arr " ; - amrex::Print() << " " << mf_arr(0,0,0,0); - amrex::Print() << " " << mf_arr(15,15,15,0); - amrex::Print() << " " << mf_arr(0,0,0,1); - amrex::Print() << " " << mf_arr(15,15,15,1); - amrex::Print() << "\n"; +//#ifdef AMREX_USE_GPU +// cudaDeviceSynchronize(); +//#endif } // Perform Fourier transform from `tmpRealField` to `tmpSpectralField` @@ -173,8 +163,9 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, // GPU stream as the above copy cudaDeviceSynchronize(); cufftResult result; - //cudaStream_t stream = amrex::Cuda::Device::cudaStream(); - //cufftSetStream ( forward_plan[mfi], stream); + cudaStream_t stream = amrex::Cuda::Device::cudaStream(); + amrex::Print() << " stream is " << stream << "\n"; + cufftSetStream ( forward_plan[mfi], stream); result = cufftExecD2Z( forward_plan[mfi], tmpRealField[mfi].dataPtr(), reinterpret_cast( @@ -182,13 +173,9 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, if ( result != CUFFT_SUCCESS ) { amrex::Print() << " cufftplan3d execute failed ! \n"; } - if ( result == CUFFT_SUCCESS ) { - amrex::Print() << " created cufft forward transform\n"; - } cudaDeviceSynchronize(); #else fftw_execute( forward_plan[mfi] ); - amrex::Print() << " forward fft on cpu\n"; #endif // Copy the spectral-space field `tmpSpectralField` to the appropriate @@ -219,8 +206,6 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, // Copy field into the right index fields_arr(i,j,k,field_index) = spectral_field_value; }); -// amrex::Print() << " in forward trans after D2Z" << fields_arr(0,0,0,0) ; - amrex::Print() << "\n"; } } } @@ -279,10 +264,11 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, #ifdef AMREX_USE_GPU // Add cuFFT-specific code ; make sure that this is done on the same // GPU stream as the above copy - cudaDeviceSynchronize(); + cudaDeviceSynchronize(); cufftResult result; - //cudaStream_t stream = amrex::Cuda::Device::cudaStream(); - //cufftSetStream ( backward_plan[mfi], stream); + cudaStream_t stream = amrex::Cuda::Device::cudaStream(); + amrex::Print() << " stream is " << stream << "\n"; + cufftSetStream ( backward_plan[mfi], stream); result = cufftExecZ2D( backward_plan[mfi], reinterpret_cast( tmpSpectralField[mfi].dataPtr()), @@ -293,10 +279,9 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, if ( result == CUFFT_SUCCESS ) { amrex::Print() << " created cufft inverse transform\n"; } - cudaDeviceSynchronize(); + cudaDeviceSynchronize(); #else fftw_execute( backward_plan[mfi] ); - amrex::Print() << " cpu inverse done\n"; #endif // Copy the temporary field `tmpRealField` to the real-space field `mf` @@ -313,15 +298,14 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, // Copy and normalize field mf_arr(i,j,k,i_comp) = inv_N*tmp_arr(i,j,k); }); +//#ifdef AMREX_USE_GPU +// cudaDeviceSynchronize(); +//#endif + #ifdef AMREX_USE_GPU cudaDeviceSynchronize(); #endif - amrex::Print() << " after backward plan in real space 0,0,0 " << mf_arr(0,0,0,0) << " tmp " << tmp_arr(0,0,0) << "\n"; - amrex::Print() << " after backward plan in real space 15, 15, 15 " << mf_arr(15,15,15,0) << " tmp " << tmp_arr(0,0,0) << "\n"; - amrex::Print() << "\n"; -#ifdef AMREX_USE_GPU - cudaDeviceSynchronize(); -#endif + amrex::Print() << " divided by 1/N \n"; } } } diff --git a/Source/FieldSolver/WarpXFFT.cpp b/Source/FieldSolver/WarpXFFT.cpp index ffeb237c0..c4e0461f9 100644 --- a/Source/FieldSolver/WarpXFFT.cpp +++ b/Source/FieldSolver/WarpXFFT.cpp @@ -405,20 +405,6 @@ WarpX::PushPSATD (int lev, amrex::Real /* dt */) rho_fp_fft[lev]->ParallelCopy(*rho_fp[lev], 0, 0, 2, 0, 0, period_fp); BL_PROFILE_VAR_STOP(blp_copy); - - for (MFIter mfi(*rho_fp_fft[lev]); mfi.isValid(); ++mfi) - { - MultiFab &mf = *rho_fp_fft[lev]; - Box realspace_bx = mf[mfi].box(); // Copy the box - Array4 mf_arr = mf[mfi].array(); - amrex::Print() << " at begin Push PSATD rho " << mf_arr(0,0,0,0) ; - amrex::Print() << " new rho " << mf_arr(0,0,0,1) << "\n"; - MultiFab &mf_orig = *rho_fp[lev]; - Array4 mforig_arr = mf_orig[mfi].array(); - amrex::Print() << " at begin Push PSATD rho " << mforig_arr(0,0,0,0) ; - amrex::Print() << " new rho " << mforig_arr(0,0,0,1) << "\n"; - } - BL_PROFILE_VAR_START(blp_push_eb); if (fft_hybrid_mpi_decomposition){ if (Efield_fp_fft[lev][0]->local_size() == 1) @@ -465,44 +451,29 @@ WarpX::PushPSATD (int lev, amrex::Real /* dt */) // Not using the hybrid decomposition auto& solver = *spectral_solver_fp[lev]; - //// Perform forward Fourier transform - //amrex::Print() << " FTT of Ex \n"; - //solver.ForwardTransform(*Efield_fp_fft[lev][0], SpectralFieldIndex::Ex); - //amrex::Print() << " FTT of Ey \n"; - //solver.ForwardTransform(*Efield_fp_fft[lev][1], SpectralFieldIndex::Ey); - //amrex::Print() << " FTT of Ez \n"; - //solver.ForwardTransform(*Efield_fp_fft[lev][2], SpectralFieldIndex::Ez); - //amrex::Print() << " FTT of Bx \n"; - //solver.ForwardTransform(*Bfield_fp_fft[lev][0], SpectralFieldIndex::Bx); - //amrex::Print() << " FTT of By \n"; - //solver.ForwardTransform(*Bfield_fp_fft[lev][1], SpectralFieldIndex::By); - //amrex::Print() << " FTT of Bz \n"; - //solver.ForwardTransform(*Bfield_fp_fft[lev][2], SpectralFieldIndex::Bz); - //amrex::Print() << " FTT of Jx \n"; - //solver.ForwardTransform(*current_fp_fft[lev][0], SpectralFieldIndex::Jx); - //amrex::Print() << " FTT of Jy \n"; - //solver.ForwardTransform(*current_fp_fft[lev][1], SpectralFieldIndex::Jy); - //amrex::Print() << " FTT of Jz \n"; - //solver.ForwardTransform(*current_fp_fft[lev][2], SpectralFieldIndex::Jz); + // Perform forward Fourier transform + solver.ForwardTransform(*Efield_fp_fft[lev][0], SpectralFieldIndex::Ex); + solver.ForwardTransform(*Efield_fp_fft[lev][1], SpectralFieldIndex::Ey); + solver.ForwardTransform(*Efield_fp_fft[lev][2], SpectralFieldIndex::Ez); + solver.ForwardTransform(*Bfield_fp_fft[lev][0], SpectralFieldIndex::Bx); + solver.ForwardTransform(*Bfield_fp_fft[lev][1], SpectralFieldIndex::By); + solver.ForwardTransform(*Bfield_fp_fft[lev][2], SpectralFieldIndex::Bz); + solver.ForwardTransform(*current_fp_fft[lev][0], SpectralFieldIndex::Jx); + solver.ForwardTransform(*current_fp_fft[lev][1], SpectralFieldIndex::Jy); + solver.ForwardTransform(*current_fp_fft[lev][2], SpectralFieldIndex::Jz); solver.ForwardTransform(*rho_fp_fft[lev], SpectralFieldIndex::rho_old, 0); solver.ForwardTransform(*rho_fp_fft[lev], SpectralFieldIndex::rho_new, 1); - ////// Advance fields in spectral space - //solver.pushSpectralFields(); - - ////// Perform backward Fourier Transform - //amrex::Print() << " BT of Ex \n"; - //solver.BackwardTransform(*Efield_fp_fft[lev][0], SpectralFieldIndex::Ex); - //amrex::Print() << " BT of Ey \n"; - //solver.BackwardTransform(*Efield_fp_fft[lev][1], SpectralFieldIndex::Ey); - //amrex::Print() << " BT of Ez \n"; - //solver.BackwardTransform(*Efield_fp_fft[lev][2], SpectralFieldIndex::Ez); - //amrex::Print() << " BT of Bx \n"; - //solver.BackwardTransform(*Bfield_fp_fft[lev][0], SpectralFieldIndex::Bx); - //amrex::Print() << " BT of By \n"; - //solver.BackwardTransform(*Bfield_fp_fft[lev][1], SpectralFieldIndex::By); - //amrex::Print() << " BT of Bz \n"; - //solver.BackwardTransform(*Bfield_fp_fft[lev][2], SpectralFieldIndex::Bz); + // Advance fields in spectral space + solver.pushSpectralFields(); + + // Perform backward Fourier Transform + solver.BackwardTransform(*Efield_fp_fft[lev][0], SpectralFieldIndex::Ex); + solver.BackwardTransform(*Efield_fp_fft[lev][1], SpectralFieldIndex::Ey); + solver.BackwardTransform(*Efield_fp_fft[lev][2], SpectralFieldIndex::Ez); + solver.BackwardTransform(*Bfield_fp_fft[lev][0], SpectralFieldIndex::Bx); + solver.BackwardTransform(*Bfield_fp_fft[lev][1], SpectralFieldIndex::By); + solver.BackwardTransform(*Bfield_fp_fft[lev][2], SpectralFieldIndex::Bz); } BL_PROFILE_VAR_STOP(blp_push_eb); @@ -519,4 +490,5 @@ WarpX::PushPSATD (int lev, amrex::Real /* dt */) { amrex::Abort("WarpX::PushPSATD: TODO"); } + amrex::Print() << " coped data from fft to valid \n"; } diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 589ee168c..b9f27f07e 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -318,16 +318,28 @@ WarpX::InitLevelData (int lev, Real time) void WarpX::InitLevelDataFFT (int lev, Real time) { + + amrex::Print() << " print out pointer " << &Efield_fp_fft[lev][0] << "\n"; Efield_fp_fft[lev][0]->setVal(0.0); + amrex::Print() << " ex set \n"; Efield_fp_fft[lev][1]->setVal(0.0); + amrex::Print() << " ey set \n"; Efield_fp_fft[lev][2]->setVal(0.0); + amrex::Print() << " ez set \n"; Bfield_fp_fft[lev][0]->setVal(0.0); + amrex::Print() << " bx set \n"; Bfield_fp_fft[lev][1]->setVal(0.0); + amrex::Print() << " by set \n"; Bfield_fp_fft[lev][2]->setVal(0.0); + amrex::Print() << " bz set \n"; current_fp_fft[lev][0]->setVal(0.0); + amrex::Print() << " jx set \n"; current_fp_fft[lev][1]->setVal(0.0); + amrex::Print() << " jy set \n"; current_fp_fft[lev][2]->setVal(0.0); + amrex::Print() << " jz set \n"; rho_fp_fft[lev]->setVal(0.0); + amrex::Print() << " rhofp set \n"; if (lev > 0) { @@ -343,14 +355,6 @@ WarpX::InitLevelDataFFT (int lev, Real time) rho_cp_fft[lev]->setVal(0.0); } - for (MFIter mfi(*rho_fp_fft[lev]); mfi.isValid(); ++mfi) - { - MultiFab &mf = *rho_fp_fft[lev]; - Box realspace_bx = mf[mfi].box(); // Copy the box - Array4 mf_arr = mf[mfi].array(); - amrex::Print() << " at initialization rho " << mf_arr(0,0,0,0) ; - amrex::Print() << " new rho " << mf_arr(0,0,0,1) << "\n"; - } } #endif diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 28971eb0c..23d97133b 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -221,7 +221,9 @@ WarpX::FillBoundaryE () { for (int lev = 0; lev <= finest_level; ++lev) { + amrex::Print() << " lev " << lev << " in fill bndry E \n"; FillBoundaryE(lev); + amrex::Print() << " lev " << lev << " after fill bndry E \n"; } } @@ -245,7 +247,7 @@ void WarpX::FillBoundaryE (int lev, PatchType patch_type) { if (patch_type == PatchType::fine) - { + { if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeE(patch_type, @@ -270,9 +272,12 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type) pml[lev]->FillBoundaryE(patch_type); } + amrex::Print() << " before defining multifab \n"; const auto& cperiod = Geom(lev-1).periodicity(); Vector mf{Efield_cp[lev][0].get(),Efield_cp[lev][1].get(),Efield_cp[lev][2].get()}; + amrex::Print() << " before amrex fill bndry \n"; amrex::FillBoundary(mf, cperiod); + amrex::Print() << " after amrex fill bndry \n"; } } diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 695faaa62..94b58ca07 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -660,7 +660,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, const std::array& cxyzmin_tile = WarpX::LowerCorner(ctilebox, lev-1); #ifdef AMREX_USE_GPU - data_ptr = (*crhomf)[pti].dataPtr(); + data_ptr = (*crhomf)[pti].dataPtr(icomp); auto rholen = (*crhomf)[pti].length(); #else tile_box = amrex::convert(ctilebox, IntVect::TheUnitVector()); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 532858556..5814d6829 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -558,7 +558,9 @@ WarpX::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& new_grids, InitLevelData(lev, time); #ifdef WARPX_USE_PSATD + amrex::Print() << " alloc level data fft \n"; AllocLevelDataFFT(lev); + amrex::Print() << " init level data fft \n"; InitLevelDataFFT(lev, time); #endif } diff --git a/Source/main.cpp b/Source/main.cpp index d89e89e47..2e1a5e9cf 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -30,10 +30,13 @@ int main(int argc, char* argv[]) const Real strt_total = amrex::second(); { + amrex::Print() << " about to construct warpx \n"; WarpX warpx; + amrex::Print() << " call warpx init data \n"; warpx.InitData(); + amrex::Print() << " call warpx evolve \n"; warpx.Evolve(); Real end_total = amrex::second() - strt_total; -- cgit v1.2.3 From 5cefcae9e6ae84553211ecece151ac8108430cf4 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 8 May 2019 14:35:30 -0700 Subject: Convert the particle position pusher to C++ --- Source/FortranInterface/WarpX_f.H | 6 --- Source/FortranInterface/WarpX_picsar.F90 | 32 ------------ Source/Particles/WarpXParticleContainer.cpp | 80 ++++++++++++++++------------- 3 files changed, 44 insertions(+), 74 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index 3d92b7651..4457e34d8 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -176,12 +176,6 @@ extern "C" const amrex::Real* charge, const amrex::Real* mass, const amrex::Real* dt, const long* particle_pusher_algo); - // Particle pusher (position) - - void warpx_particle_pusher_positions(const long* np, - amrex::Real* xp, amrex::Real* yp, amrex::Real* zp, - amrex::Real* uxp, amrex::Real* uyp, amrex::Real* uzp, amrex::Real* gaminv, - const amrex::Real* dt); // Laser pusher diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90 index 6d6a2e411..a4cc99926 100644 --- a/Source/FortranInterface/WarpX_picsar.F90 +++ b/Source/FortranInterface/WarpX_picsar.F90 @@ -519,38 +519,6 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n end subroutine warpx_particle_pusher_momenta - ! _________________________________________________________________ - !> - !> @brief - !> Main subroutine for the particle pusher of positions - !> - !> @param[in] np number of super-particles - !> @param[in] xp,yp,zp particle position arrays - !> @param[in] uxp,uyp,uzp normalized momentum in each direction - !> @param[in] gaminv particle Lorentz factors - !> @param[in] dt time step - !> @param[in] particle_pusher_algo Particle pusher algorithm - subroutine warpx_particle_pusher_positions(np,xp,yp,zp,uxp,uyp,uzp, & - gaminv,dt) & - bind(C, name="warpx_particle_pusher_positions") - - INTEGER(c_long), INTENT(IN) :: np - REAL(amrex_real),INTENT(INOUT) :: gaminv(np) - REAL(amrex_real),INTENT(INOUT) :: xp(np),yp(np),zp(np) - REAL(amrex_real),INTENT(INOUT) :: uxp(np),uyp(np),uzp(np) - REAL(amrex_real),INTENT(IN) :: dt - - CALL pxr_set_gamma(np,uxp,uyp,uzp,gaminv) - - !!!! --- push particle species positions a time step -#if (AMREX_SPACEDIM == 3) || (defined WARPX_RZ) - CALL pxr_pushxyz(np,xp,yp,zp,uxp,uyp,uzp,gaminv,dt) -#elif (AMREX_SPACEDIM == 2) - CALL pxr_pushxz(np,xp,zp,uxp,uzp,gaminv,dt) -#endif - - end subroutine warpx_particle_pusher_positions - ! _________________________________________________________________ !> !> @brief diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 80f3882a0..ab6e71d40 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -86,9 +86,9 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) 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 @@ -174,7 +174,7 @@ WarpXParticleContainer::AddNParticles (int lev, int n, const Real* x, const Real* y, const Real* z, const Real* vx, const Real* vy, const Real* vz, int nattr, const Real* attr, int uniqueparticles, int id) -{ +{ BL_ASSERT(nattr == 1); const Real* weight = attr; @@ -230,15 +230,15 @@ WarpXParticleContainer::AddNParticles (int lev, #endif p.pos(1) = z[i]; #endif - + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["xold"], x[i]); particle_tile.push_back_real(particle_comps["yold"], y[i]); - particle_tile.push_back_real(particle_comps["zold"], z[i]); + particle_tile.push_back_real(particle_comps["zold"], z[i]); } - + particle_tile.push_back(p); } @@ -254,9 +254,9 @@ WarpXParticleContainer::AddNParticles (int lev, auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); particle_tile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend); - particle_tile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); + particle_tile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); } - + for (int comp = PIdx::uz+1; comp < PIdx::nattribs; ++comp) { #ifdef WARPX_RZ @@ -327,7 +327,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, jx_ptr = local_jx[thread_num].dataPtr(); jy_ptr = local_jy[thread_num].dataPtr(); jz_ptr = local_jz[thread_num].dataPtr(); - + local_jx[thread_num].setVal(0.0); local_jy[thread_num].setVal(0.0); local_jz[thread_num].setVal(0.0); @@ -426,7 +426,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, BL_PROFILE_VAR_STOP(blp_pxr_cd); -#ifndef AMREX_USE_GPU +#ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); jx[pti].atomicAdd(local_jx[thread_num], tbx, tbx, 0, 0, 1); @@ -477,7 +477,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, auto jyntot = local_jy[thread_num].length(); auto jzntot = local_jz[thread_num].length(); #endif - + long ncrse = np - np_current; BL_PROFILE_VAR_START(blp_pxr_cd); if (j_is_nodal) { @@ -648,7 +648,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, #endif BL_PROFILE_VAR_STOP(blp_pxr_chd); -#ifndef AMREX_USE_GPU +#ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); (*rhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp, 1); @@ -708,7 +708,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, #endif BL_PROFILE_VAR_STOP(blp_pxr_chd); -#ifndef AMREX_USE_GPU +#ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); (*crhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp, 1); @@ -1024,7 +1024,7 @@ WarpXParticleContainer::PushX (int lev, Real dt) { BL_PROFILE("WPC::PushX()"); BL_PROFILE_VAR_NS("WPC::PushX::Copy", blp_copy); - BL_PROFILE_VAR_NS("WPC:PushX::Push", blp_pxr_pp); + BL_PROFILE_VAR_NS("WPC:PushX::Push", blp_pp); if (do_not_push) return; @@ -1040,16 +1040,6 @@ WarpXParticleContainer::PushX (int lev, Real dt) { Real wt = amrex::second(); - auto& attribs = pti.GetAttribs(); - - auto& uxp = attribs[PIdx::ux]; - auto& uyp = attribs[PIdx::uy]; - auto& uzp = attribs[PIdx::uz]; - - const long np = pti.numParticles(); - - giv.resize(np); - // // copy data from particle container to temp arrays // @@ -1060,14 +1050,32 @@ WarpXParticleContainer::PushX (int lev, Real dt) // // Particle Push // - BL_PROFILE_VAR_START(blp_pxr_pp); - warpx_particle_pusher_positions(&np, - xp.dataPtr(), - yp.dataPtr(), - zp.dataPtr(), - uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), - giv.dataPtr(), &dt); - BL_PROFILE_VAR_STOP(blp_pxr_pp); + BL_PROFILE_VAR_START(blp_pp); + // Extract pointers to particle position and momenta, for this particle tile + Real* AMREX_RESTRICT x = xp.dataPtr(); + Real* AMREX_RESTRICT y = yp.dataPtr(); + Real* AMREX_RESTRICT z = zp.dataPtr(); + auto& attribs = pti.GetAttribs(); + Real* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + Real* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + Real* AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); + // Loop over the particles and update the position + const long np = pti.numParticles(); + const Real inv_c2 = 1./(PhysConst::c*PhysConst::c); + amrex::ParallelFor( np, + [=] AMREX_GPU_DEVICE (long i) { + // Compute inverse Lorentz factor + const Real inv_gamma = 1./std::sqrt( + 1. + (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])*inv_c2); + // Update positions over one time step + x[i] += ux[i] * inv_gamma * dt; +#if (AMREX_SPACEDIM == 3) + y[i] += uy[i] * inv_gamma * dt; +#endif + z[i] += uz[i] * inv_gamma * dt; + } + ); + BL_PROFILE_VAR_STOP(blp_pp); // // copy particle data back @@ -1090,15 +1098,15 @@ WarpXParticleContainer::PushX (int lev, Real dt) } // This function is called in Redistribute, just after locate -void -WarpXParticleContainer::particlePostLocate(ParticleType& p, +void +WarpXParticleContainer::particlePostLocate(ParticleType& p, const ParticleLocData& pld, const int lev) { // Tag particle if goes to higher level. // It will be split later in the loop - if (pld.m_lev == lev+1 - and p.m_idata.id != NoSplitParticleID + if (pld.m_lev == lev+1 + and p.m_idata.id != NoSplitParticleID and p.m_idata.id >= 0) { p.m_idata.id = DoSplitParticleID; -- cgit v1.2.3 From f7784a2641b9ef8f5c045cf282dd0069f903997d Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 8 May 2019 15:01:26 -0700 Subject: Fix RZ bug --- Source/Particles/WarpXParticleContainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index ab6e71d40..4bcd43472 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -1069,7 +1069,7 @@ WarpXParticleContainer::PushX (int lev, Real dt) 1. + (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])*inv_c2); // Update positions over one time step x[i] += ux[i] * inv_gamma * dt; -#if (AMREX_SPACEDIM == 3) +#if (AMREX_SPACEDIM == 3) || (defined WARPX_RZ) // RZ pushes particles in 3D y[i] += uy[i] * inv_gamma * dt; #endif z[i] += uz[i] * inv_gamma * dt; -- cgit v1.2.3 From 5b43b8220d4d817615c2d360cd8e647277e8d273 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 8 May 2019 15:39:29 -0700 Subject: Use array of struct in pusher for positions --- Source/Particles/WarpXParticleContainer.cpp | 31 ++++++----------------------- 1 file changed, 6 insertions(+), 25 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 4bcd43472..1abbd747d 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -1023,8 +1023,6 @@ void WarpXParticleContainer::PushX (int lev, Real dt) { BL_PROFILE("WPC::PushX()"); - BL_PROFILE_VAR_NS("WPC::PushX::Copy", blp_copy); - BL_PROFILE_VAR_NS("WPC:PushX::Push", blp_pp); if (do_not_push) return; @@ -1034,27 +1032,18 @@ WarpXParticleContainer::PushX (int lev, Real dt) #pragma omp parallel #endif { - Cuda::ManagedDeviceVector xp, yp, zp, giv; for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { Real wt = amrex::second(); - // - // copy data from particle container to temp arrays - // - BL_PROFILE_VAR_START(blp_copy); - pti.GetPosition(xp, yp, zp); - BL_PROFILE_VAR_STOP(blp_copy); - // // Particle Push // - BL_PROFILE_VAR_START(blp_pp); // Extract pointers to particle position and momenta, for this particle tile - Real* AMREX_RESTRICT x = xp.dataPtr(); - Real* AMREX_RESTRICT y = yp.dataPtr(); - Real* AMREX_RESTRICT z = zp.dataPtr(); + // - positions are stored as an array of struct, in `ParticleType` + ParticleType * AMREX_RESTRICT pstructs = &(pti.GetArrayOfStructs()[0]); + // - momenta are stored as a struct of array, in `attribs` auto& attribs = pti.GetAttribs(); Real* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); Real* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); @@ -1068,21 +1057,13 @@ WarpXParticleContainer::PushX (int lev, Real dt) const Real inv_gamma = 1./std::sqrt( 1. + (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])*inv_c2); // Update positions over one time step - x[i] += ux[i] * inv_gamma * dt; + pstructs[i].pos(0) += ux[i] * inv_gamma * dt; #if (AMREX_SPACEDIM == 3) || (defined WARPX_RZ) // RZ pushes particles in 3D - y[i] += uy[i] * inv_gamma * dt; + pstructs[i].pos(1) += uy[i] * inv_gamma * dt; #endif - z[i] += uz[i] * inv_gamma * dt; + pstructs[i].pos(2) += uz[i] * inv_gamma * dt; } ); - BL_PROFILE_VAR_STOP(blp_pp); - - // - // copy particle data back - // - BL_PROFILE_VAR_START(blp_copy); - pti.SetPosition(xp, yp, zp); - BL_PROFILE_VAR_STOP(blp_copy); if (cost) { const Box& tbx = pti.tilebox(); -- cgit v1.2.3 From 271eb52eafffdefb8dfc3bb43ee521094fbd5a74 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Fri, 10 May 2019 17:49:29 -0700 Subject: lab diags is a species argument. add comments --- Source/Diagnostics/BoostedFrameDiagnostic.H | 4 +++ Source/Diagnostics/BoostedFrameDiagnostic.cpp | 25 ++++++++------- Source/Particles/MultiParticleContainer.H | 5 +++ Source/Particles/MultiParticleContainer.cpp | 36 ++++++++++++++++++++-- Source/Particles/PhysicalParticleContainer.cpp | 15 +++++---- .../Particles/RigidInjectedParticleContainer.cpp | 2 +- Source/Particles/WarpXParticleContainer.H | 8 +++++ Source/Particles/WarpXParticleContainer.cpp | 6 ++-- Source/WarpX.cpp | 8 +++-- 9 files changed, 84 insertions(+), 25 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.H b/Source/Diagnostics/BoostedFrameDiagnostic.H index e35d307a6..ef4bd2ec1 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.H +++ b/Source/Diagnostics/BoostedFrameDiagnostic.H @@ -75,6 +75,10 @@ class BoostedFrameDiagnostic { int boost_direction_; amrex::Vector > data_buffer_; + // particles_buffer_ is current blind to refinement level. + // particles_buffer_[i][j] is a WarpXParticleContainer::DiagnosticParticleData where + // - i is the back-transformed snapshot number + // - j is the species number amrex::Vector > particles_buffer_; int num_buffer_ = 256; int max_box_size_ = 256; diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.cpp b/Source/Diagnostics/BoostedFrameDiagnostic.cpp index 13972075d..5d85fc8f8 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.cpp +++ b/Source/Diagnostics/BoostedFrameDiagnostic.cpp @@ -540,13 +540,14 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) if (WarpX::do_boosted_frame_particles) { for (int j = 0; j < mypc.nSpecies(); ++j) { + std::string species_name = species_names[mypc.map_species_lab_diags[j]]; #ifdef WARPX_USE_HDF5 writeParticleDataHDF5(particles_buffer_[i][j], snapshots_[i].file_name, - species_names[j]); + species_name); #else std::stringstream part_ss; - part_ss << snapshots_[i].file_name + "/" + species_names[j] + "/"; + part_ss << snapshots_[i].file_name + "/" + species_name + "/"; writeParticleData(particles_buffer_[i][j], part_ss.str(), i_lab); #endif } @@ -600,7 +601,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, DistributionMapping buff_dm(buff_ba); data_buffer_[i].reset( new MultiFab(buff_ba, buff_dm, ncomp, 0) ); } - if (WarpX::do_boosted_frame_particles) particles_buffer_[i].resize(mypc.nSpecies()); + if (WarpX::do_boosted_frame_particles) particles_buffer_[i].resize(mypc.nspecies_lab_frame_diags); } if (WarpX::do_boosted_frame_fields) { @@ -666,14 +667,15 @@ writeLabFrameData(const MultiFab* cell_centered_data, } if (WarpX::do_boosted_frame_particles) { - for (int j = 0; j < mypc.nSpecies(); ++j) { + for (int j = 0; j < mypc.nspecies_lab_frame_diags; ++j) { + const std::string species_name = species_names[mypc.map_species_lab_diags[j]]; #ifdef WARPX_USE_HDF5 writeParticleDataHDF5(particles_buffer_[i][j], snapshots_[i].file_name, - species_names[j]); + species_name); #else std::stringstream part_ss; - part_ss << snapshots_[i].file_name + "/" + species_names[j] + "/"; + part_ss << snapshots_[i].file_name + "/" + species_name + "/"; writeParticleData(particles_buffer_[i][j], part_ss.str(), i_lab); #endif } @@ -855,16 +857,16 @@ LabSnapShot(Real t_lab_in, Real t_boost, Real zmin_lab_in, ParallelDescriptor::Barrier(); - if (WarpX::do_boosted_frame_particles) - { + if (WarpX::do_boosted_frame_particles){ auto & mypc = WarpX::GetInstance().GetPartContainer(); const std::vector species_names = mypc.GetSpeciesNames(); for (int j = 0; j < mypc.nSpecies(); ++j) { - output_create_species_group(file_name, species_names[j]); + std::string species_name = species_names[mypc.map_species_lab_diags[j]]; + output_create_species_group(file_name, species_name); for (int k = 0; k < static_cast(particle_field_names.size()); ++k) { - std::string field_path = species_names[j] + "/" + particle_field_names[k]; + std::string field_path = species_name + "/" + particle_field_names[k]; output_create_particle_field(file_name, field_path); } } @@ -888,7 +890,8 @@ LabSnapShot(Real t_lab_in, Real t_boost, Real zmin_lab_in, const std::string particles_prefix = "particle"; for(int i = 0; i < nspecies; ++i) { - const std::string fullpath = file_name + "/" + species_names[i]; + std::string species_name = species_names[mypc.map_species_lab_diags[i]]; + const std::string fullpath = file_name + "/" + species_name; if (!UtilCreateDirectory(fullpath, 0755)) CreateDirectoryFailed(fullpath); } diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 9291e0358..5a79443d0 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -153,6 +153,8 @@ public: void SetParticleBoxArray (int lev, amrex::BoxArray& new_ba); void SetParticleDistributionMap (int lev, amrex::DistributionMapping& new_dm); + void setSpeciesLabFrameDiags() const; + int nSpecies() const {return nspecies;} int nSpeciesDepositOnMainGrid () const { @@ -184,6 +186,9 @@ public: // Number of coefficients for the stencil of the NCI corrector. // The stencil is applied in the z direction only. static constexpr int nstencilz_fdtd_nci_corr=5; + int nspecies_lab_frame_diags = 0; + std::vector map_species_lab_diags; + int do_boosted_frame_diags = 0; amrex::Vector > fdtd_nci_stencilz_ex; amrex::Vector > fdtd_nci_stencilz_by; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 440906348..bde8d244e 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -31,7 +31,21 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) pc_tmp.reset(new PhysicalParticleContainer(amr_core)); - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + // Compute the number of species for which mab frame data is dumped + // nspecies_lab_frame_diags, and map their ID to MultiParticleContainer + // particle IDs in map_species_lab_diags. + map_species_lab_diags.resize(nspecies); + nspecies_lab_frame_diags = 0; + for (int i=0; ido_boosted_frame_diags){ + map_species_lab_diags[nspecies_lab_frame_diags] = i; + nspecies_lab_frame_diags += 1; + do_boosted_frame_diags = 1; + } + } + + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { for (int i = 0; i < nspecies + nlasers; ++i) { @@ -376,13 +390,24 @@ MultiParticleContainer BL_PROFILE("MultiParticleContainer::GetLabFrameData"); + // Loop over particle species for (int i = 0; i < nspecies; ++i){ WarpXParticleContainer* pc = allcontainers[i].get(); WarpXParticleContainer::DiagnosticParticles diagnostic_particles; pc->GetParticleSlice(direction, z_old, z_new, t_boost, t_lab, dt, diagnostic_particles); - + // Here, diagnostic_particles[lev][index] is a WarpXParticleContainer::DiagnosticParticleData + // where "lev" is the AMR level and "index" is a [grid index][tile index] pair. + + // Loop over AMR levels for (int lev = 0; lev <= pc->finestLevel(); ++lev){ + // Loop over [grid index][tile index] pairs + // and Fills parts[species number i] with particle data from all grids and + // tiles in diagnostic_particles. parts contains particles from all + // AMR levels indistinctly. for (auto it = diagnostic_particles[lev].begin(); it != diagnostic_particles[lev].end(); ++it){ + // it->first is the [grid index][tile index] key + // it->second is the corresponding + // WarpXParticleContainer::DiagnosticParticleData value parts[i].GetRealData(DiagIdx::w).insert( parts[i].GetRealData(DiagIdx::w ).end(), it->second.GetRealData(DiagIdx::w ).begin(), it->second.GetRealData(DiagIdx::w ).end()); @@ -459,3 +484,10 @@ MultiParticleContainer::doContinuousInjection() const } return warpx_do_continuous_injection; } + +// Set number of species for which lab frame data is dumped +// and maps their ID to MultiParticleContainer IDs. +//void +//MultiParticleContainer::setSpeciesLabFrameDiags() const +//{ +//} diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index b3c598c22..d99cc9c66 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -82,6 +82,9 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_splitting", do_splitting); pp.query("split_type", split_type); pp.query("do_continuous_injection", do_continuous_injection); + // Whether to plot back-transformed (lab-frame) diagnostics + // for this species. + pp.query("do_boosted_frame_diags", do_boosted_frame_diags); pp.query("plot_species", plot_species); int do_user_plot_vars; @@ -90,14 +93,14 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp // By default, all particle variables are dumped to plotfiles, // including {x,y,z,ux,uy,uz}old variables when running in a // boosted frame - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles){ + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){ plot_flags.resize(PIdx::nattribs + 6, 1); } else { plot_flags.resize(PIdx::nattribs, 1); } } else { // Set plot_flag to 0 for all attribs - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles){ + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){ plot_flags.resize(PIdx::nattribs + 6, 0); } else { plot_flags.resize(PIdx::nattribs, 0); @@ -216,7 +219,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, attribs[PIdx::uz] = u[2]; attribs[PIdx::w ] = weight; - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["xold"], x); @@ -500,7 +503,7 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) attribs[PIdx::uy] = u[1]; attribs[PIdx::uz] = u[2]; - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); particle_tile.push_back_real(particle_comps["xold"], x); @@ -742,7 +745,7 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) attribs[PIdx::uz] = u[2]; // note - this will be slow on the GPU, need to revisit - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); particle_tile.push_back_real(particle_comps["xold"], x); @@ -1715,7 +1718,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, auto& Bzp = attribs[PIdx::Bz]; const long np = pti.numParticles(); - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& xpold = pti.GetAttribs(particle_comps["xold"]); auto& ypold = pti.GetAttribs(particle_comps["yold"]); diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index a5acca281..fd1b2dfb5 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -225,7 +225,7 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, auto& Bzp = attribs[PIdx::Bz]; const long np = pti.numParticles(); - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& xpold = pti.GetAttribs(particle_comps["xold"]); auto& ypold = pti.GetAttribs(particle_comps["yold"]); diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 51dc5ec05..600061e8d 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -85,7 +85,13 @@ class WarpXParticleContainer public: friend MultiParticleContainer; + // amrex::StructOfArrays with DiagIdx::nattribs amrex::Real components + // for the particle data and 0 int components. using DiagnosticParticleData = amrex::StructOfArrays; + // DiagnosticParticles is a vector, with one element per MR level. + // DiagnosticParticles[lev] is typically a key-value pair where the key is + // a pair [grid_index, tile_index], and the value is the corresponding + // DiagnosticParticleData (see above) on this tile. using DiagnosticParticles = amrex::Vector, DiagnosticParticleData> >; WarpXParticleContainer (amrex::AmrCore* amr_core, int ispecies); @@ -265,6 +271,8 @@ protected: // support all features allowed by direct injection. int do_continuous_injection = 0; + int do_boosted_frame_diags = 0; + amrex::Vector local_rho; amrex::Vector local_jx; amrex::Vector local_jy; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 1abbd747d..0cf5c10b4 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -78,7 +78,7 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) particle_comps["theta"] = PIdx::theta; #endif - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { particle_comps["xold"] = PIdx::nattribs; particle_comps["yold"] = PIdx::nattribs+1; @@ -231,7 +231,7 @@ WarpXParticleContainer::AddNParticles (int lev, p.pos(1) = z[i]; #endif - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["xold"], x[i]); @@ -249,7 +249,7 @@ WarpXParticleContainer::AddNParticles (int lev, particle_tile.push_back_real(PIdx::uy, vy + ibegin, vy + iend); particle_tile.push_back_real(PIdx::uz, vz + ibegin, vz + iend); - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index eb84af2c7..c863e9a4c 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -56,9 +56,12 @@ int WarpX::num_mirrors = 0; int WarpX::sort_int = -1; -bool WarpX::do_boosted_frame_diagnostic = false; int WarpX::num_snapshots_lab = std::numeric_limits::lowest(); Real WarpX::dt_snapshots_lab = std::numeric_limits::lowest(); +// bool WarpX::do_boosted_frame_diagnostic = false; +// bool WarpX::do_boosted_frame_fields = true; +// bool WarpX::do_boosted_frame_particles = true; +bool WarpX::do_boosted_frame_diagnostic = false; bool WarpX::do_boosted_frame_fields = true; bool WarpX::do_boosted_frame_particles = true; @@ -117,7 +120,7 @@ WarpX::ResetInstance () { delete m_instance; m_instance = nullptr; -} +} WarpX::WarpX () { @@ -157,6 +160,7 @@ WarpX::WarpX () current_injection_position = geom[0].ProbLo(moving_window_dir); } } + do_boosted_frame_particles = mypc->do_boosted_frame_diags; Efield_aux.resize(nlevs_max); Bfield_aux.resize(nlevs_max); -- cgit v1.2.3 From c952dbe2e3d9a2c7bab2774e8036e0e9ff72e0ed Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Sun, 12 May 2019 07:58:38 -0700 Subject: only selected species BFD-dumped, but all old attribs initialized --- Source/Diagnostics/BoostedFrameDiagnostic.cpp | 30 ++++---------------------- Source/Particles/MultiParticleContainer.cpp | 19 ++++++++++------ Source/Particles/PhysicalParticleContainer.cpp | 9 +++++--- Source/Particles/WarpXParticleContainer.cpp | 9 +++++--- Source/WarpX.cpp | 2 +- 5 files changed, 29 insertions(+), 40 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.cpp b/Source/Diagnostics/BoostedFrameDiagnostic.cpp index 934610b9e..7a44bc66a 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.cpp +++ b/Source/Diagnostics/BoostedFrameDiagnostic.cpp @@ -582,27 +582,19 @@ writeLabFrameData(const MultiFab* cell_centered_data, const std::vector species_names = mypc.GetSpeciesNames(); - Print()<<"in BoostedFrameDiagnostic::writeLabFrameData 1\n"; - for (int i = 0; i < N_snapshots_; ++i) { - std::cout<<"in BoostedFrameDiagnostic::writeLabFrameData 2\n"; - std::cout<<"i "< zhi_boost) or (snapshots_[i].current_z_lab < snapshots_[i].zmin_lab) or (snapshots_[i].current_z_lab > snapshots_[i].zmax_lab) ) continue; - std::cout<<"toto 3\n"; int i_lab = (snapshots_[i].current_z_lab - snapshots_[i].zmin_lab) / dz_lab_; - std::cout<<"in BoostedFrameDiagnostic::writeLabFrameData 3\n"; if (buff_counter_[i] == 0) { if (WarpX::do_boosted_frame_fields) { @@ -617,7 +609,6 @@ writeLabFrameData(const MultiFab* cell_centered_data, } if (WarpX::do_boosted_frame_particles) particles_buffer_[i].resize(mypc.nspecies_lab_frame_diags); } - std::cout<<"in BoostedFrameDiagnostic::writeLabFrameData 4\n"; if (WarpX::do_boosted_frame_fields) { const int ncomp = cell_centered_data->nComp(); @@ -657,14 +648,11 @@ writeLabFrameData(const MultiFab* cell_centered_data, &ncomp, &i_boost, &i_lab); } } - std::cout<<"in BoostedFrameDiagnostic::writeLabFrameData 5\n"; if (WarpX::do_boosted_frame_particles) { - std::cout<<"in BoostedFrameDiagnostic::writeLabFrameData 6\n"; - //mypc.GetLabFrameData(snapshots_[i].file_name, i_lab, boost_direction_, - // old_z_boost, snapshots_[i].current_z_boost, - // t_boost, snapshots_[i].t_lab, dt, particles_buffer_[i]); - std::cout<<"in BoostedFrameDiagnostic::writeLabFrameData 7\n"; + mypc.GetLabFrameData(snapshots_[i].file_name, i_lab, boost_direction_, + old_z_boost, snapshots_[i].current_z_boost, + t_boost, snapshots_[i].t_lab, dt, particles_buffer_[i]); } @@ -686,28 +674,21 @@ writeLabFrameData(const MultiFab* cell_centered_data, } if (WarpX::do_boosted_frame_particles) { - std::cout<<"in BoostedFrameDiagnostic::writeLabFrameData 8\n"; for (int j = 0; j < mypc.nspecies_lab_frame_diags; ++j) { - std::cout<<"in BoostedFrameDiagnostic::writeLabFrameData 9\n"; const std::string species_name = species_names[mypc.map_species_lab_diags[j]]; #ifdef WARPX_USE_HDF5 - std::cout<<"particles_buffer_ j "<(particle_field_names.size()); ++k) @@ -919,10 +899,8 @@ LabSnapShot(Real t_lab_in, Real t_boost, Real zmin_lab_in, auto & mypc = WarpX::GetInstance().GetPartContainer(); const std::vector species_names = mypc.GetSpeciesNames(); -// int nspecies = mypc.nSpecies(); const std::string particles_prefix = "particle"; -// for(int i = 0; i < nspecies; ++i) { for(int i = 0; i < mypc.nspecies_lab_frame_diags; ++i) { int is = mypc.map_species_lab_diags[i]; std::string species_name = species_names[is]; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 841c41835..337005dc3 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -55,7 +55,19 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) allcontainers[i]->AddRealComp("uxold"); allcontainers[i]->AddRealComp("uyold"); allcontainers[i]->AddRealComp("uzold"); + } + /* + for (int i = 0; i < nspecies_lab_frame_diags; ++i) + { + int is = map_species_lab_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"); @@ -389,17 +401,13 @@ MultiParticleContainer { BL_PROFILE("MultiParticleContainer::GetLabFrameData"); - std::cout<<"GetLabFrameData 1\n"; // Loop over particle species for (int i = 0; i < nspecies_lab_frame_diags; ++i){ int isp = map_species_lab_diags[i]; -std::cout<<"GetLabFrameData 2\n"; WarpXParticleContainer* pc = allcontainers[isp].get(); - std::cout<<"getparticleslice "<< isp<GetParticleSlice(direction, z_old, z_new, t_boost, t_lab, dt, diagnostic_particles); -std::cout<<"GetLabFrameData 3\n"; // Here, diagnostic_particles[lev][index] is a WarpXParticleContainer::DiagnosticParticleData // where "lev" is the AMR level and "index" is a [grid index][tile index] pair. @@ -409,16 +417,13 @@ std::cout<<"GetLabFrameData 3\n"; // and Fills parts[species number i] with particle data from all grids and // tiles in diagnostic_particles. parts contains particles from all // AMR levels indistinctly. -std::cout<<"GetLabFrameData 4\n"; for (auto it = diagnostic_particles[lev].begin(); it != diagnostic_particles[lev].end(); ++it){ // it->first is the [grid index][tile index] key // it->second is the corresponding // WarpXParticleContainer::DiagnosticParticleData value -std::cout<<"GetLabFrameData 5\n"; parts[i].GetRealData(DiagIdx::w).insert( parts[i].GetRealData(DiagIdx::w ).end(), it->second.GetRealData(DiagIdx::w ).begin(), it->second.GetRealData(DiagIdx::w ).end()); -std::cout<<"GetLabFrameData 6\n"; parts[i].GetRealData(DiagIdx::x).insert( parts[i].GetRealData(DiagIdx::x ).end(), it->second.GetRealData(DiagIdx::x ).begin(), diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index c4d4abfd6..a3f555bd0 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -219,7 +219,8 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, attribs[PIdx::uz] = u[2]; attribs[PIdx::w ] = weight; - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + if (WarpX::do_boosted_frame_diagnostic) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["xold"], x); @@ -503,7 +504,8 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) attribs[PIdx::uy] = u[1]; attribs[PIdx::uz] = u[2]; - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + if (WarpX::do_boosted_frame_diagnostic) { auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); particle_tile.push_back_real(particle_comps["xold"], x); @@ -745,7 +747,8 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) attribs[PIdx::uz] = u[2]; // note - this will be slow on the GPU, need to revisit - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + if (WarpX::do_boosted_frame_diagnostic) { auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); particle_tile.push_back_real(particle_comps["xold"], x); diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 0cf5c10b4..9bee0031a 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -78,7 +78,8 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) particle_comps["theta"] = PIdx::theta; #endif - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + if (WarpX::do_boosted_frame_diagnostic) { particle_comps["xold"] = PIdx::nattribs; particle_comps["yold"] = PIdx::nattribs+1; @@ -231,7 +232,8 @@ WarpXParticleContainer::AddNParticles (int lev, p.pos(1) = z[i]; #endif - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + if (WarpX::do_boosted_frame_diagnostic) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["xold"], x[i]); @@ -249,7 +251,8 @@ WarpXParticleContainer::AddNParticles (int lev, particle_tile.push_back_real(PIdx::uy, vy + ibegin, vy + iend); particle_tile.push_back_real(PIdx::uz, vz + ibegin, vz + iend); - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + if (WarpX::do_boosted_frame_diagnostic) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index c863e9a4c..7c026e105 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -329,7 +329,7 @@ WarpX::ReadParameters () pp.get("gamma_boost", gamma_boost); pp.query("do_boosted_frame_fields", do_boosted_frame_fields); - pp.query("do_boosted_frame_particles", do_boosted_frame_particles); + // pp.query("do_boosted_frame_particles", do_boosted_frame_particles); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_moving_window, -- cgit v1.2.3 From 9a26a71845fde091c7840772bef1b23dbc46d6ac Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Sun, 12 May 2019 10:17:28 -0700 Subject: old attribs not allocated if species not BFD --- Source/Diagnostics/BoostedFrameDiagnostic.cpp | 22 ++++++++++------------ Source/Evolve/WarpXEvolveEM.cpp | 2 -- Source/Laser/LaserParticleContainer.cpp | 5 +++-- Source/Particles/MultiParticleContainer.H | 11 ++++++++--- Source/Particles/MultiParticleContainer.cpp | 7 +++++-- Source/Particles/PhysicalParticleContainer.cpp | 12 +++++------- Source/Particles/WarpXParticleContainer.H | 1 - Source/Particles/WarpXParticleContainer.cpp | 19 +++++++++++++------ Source/WarpX.cpp | 2 +- 9 files changed, 45 insertions(+), 36 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.cpp b/Source/Diagnostics/BoostedFrameDiagnostic.cpp index 7a44bc66a..413730f91 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.cpp +++ b/Source/Diagnostics/BoostedFrameDiagnostic.cpp @@ -497,8 +497,6 @@ BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, void BoostedFrameDiagnostic::Flush(const Geometry& geom) { BL_PROFILE("BoostedFrameDiagnostic::Flush"); - - std::cout<<"in Flush\n"; VisMF::Header::Version current_version = VisMF::GetHeaderVersion(); VisMF::SetHeaderVersion(amrex::VisMF::Header::NoFabHeader_v1); @@ -542,11 +540,10 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) if (WarpX::do_boosted_frame_particles) { // for (int j = 0; j < mypc.nSpecies(); ++j) { - for (int j = 0; j < mypc.nspecies_lab_frame_diags; ++j) { - int js = mypc.map_species_lab_diags[j]; + for (int j = 0; j < mypc.nSpeciesLabFrameDiags(); ++j) { + int js = mypc.mapSpeciesLabDiags(j); std::string species_name = species_names[js]; #ifdef WARPX_USE_HDF5 - std::cout<<"particles_buffer_ j "<(particle_field_names.size()); ++k) { std::string field_path = species_name + "/" + particle_field_names[k]; @@ -872,9 +870,9 @@ LabSnapShot(Real t_lab_in, Real t_boost, Real zmin_lab_in, auto & mypc = WarpX::GetInstance().GetPartContainer(); const std::vector species_names = mypc.GetSpeciesNames(); // for (int j = 0; j < mypc.nSpecies(); ++j) - for (int j = 0; j < mypc.nspecies_lab_frame_diags; ++j) + for (int j = 0; j < mypc.nSpeciesLabFrameDiags(); ++j) { - int js = mypc.map_species_lab_diags[j]; + int js = mypc.mapSpeciesLabDiags(j); std::string species_name = species_names[js]; output_create_species_group(file_name, species_name); for (int k = 0; k < static_cast(particle_field_names.size()); ++k) @@ -901,8 +899,8 @@ LabSnapShot(Real t_lab_in, Real t_boost, Real zmin_lab_in, const std::vector species_names = mypc.GetSpeciesNames(); const std::string particles_prefix = "particle"; - for(int i = 0; i < mypc.nspecies_lab_frame_diags; ++i) { - int is = mypc.map_species_lab_diags[i]; + for(int i = 0; i < mypc.nSpeciesLabFrameDiags(); ++i) { + int is = mypc.mapSpeciesLabDiags(i); std::string species_name = species_names[is]; const std::string fullpath = file_name + "/" + species_name; if (!UtilCreateDirectory(fullpath, 0755)) diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 5f6760f22..d88f178df 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -173,9 +173,7 @@ WarpX::EvolveEM (int numsteps) if (WarpX::do_boosted_frame_fields) { cell_centered_data = GetCellCenteredData(); } - std::cout<<"before myBFD->writeLabFrameData"<writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); - std::cout<<"after myBFD->writeLabFrameData"<::max(); - - ParmParse pp(laser_name); + do_boosted_frame_diags = 0; + + ParmParse pp(laser_name); // Parse the type of laser profile and set the corresponding flag `profile` std::string laser_type_s; diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 5a79443d0..217f727b0 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -157,6 +157,10 @@ public: int nSpecies() const {return nspecies;} + int nSpeciesLabFrameDiags() const {return nspecies_lab_frame_diags;} + int mapSpeciesLabDiags(int i) const {return map_species_lab_diags[i];} + int doBoostedFrameDiags() const {return do_boosted_frame_diags;} + int nSpeciesDepositOnMainGrid () const { int r = 0; for (int i : deposit_on_main_grid) { @@ -186,9 +190,6 @@ public: // Number of coefficients for the stencil of the NCI corrector. // The stencil is applied in the z direction only. static constexpr int nstencilz_fdtd_nci_corr=5; - int nspecies_lab_frame_diags = 0; - std::vector map_species_lab_diags; - int do_boosted_frame_diags = 0; amrex::Vector > fdtd_nci_stencilz_ex; amrex::Vector > fdtd_nci_stencilz_by; @@ -219,6 +220,10 @@ private: void ReadParameters (); + int nspecies_lab_frame_diags = 0; + std::vector map_species_lab_diags; + int do_boosted_frame_diags = 0; + // runtime parameters int nlasers = 0; int nspecies = 1; // physical particles only. nspecies+nlasers == allcontainers.size(). diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 337005dc3..508f3f606 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -47,6 +47,8 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { + //maxoldattribs + /* for (int i = 0; i < nspecies + nlasers; ++i) { allcontainers[i]->AddRealComp("xold"); @@ -56,7 +58,7 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) allcontainers[i]->AddRealComp("uyold"); allcontainers[i]->AddRealComp("uzold"); } - /* + */ for (int i = 0; i < nspecies_lab_frame_diags; ++i) { int is = map_species_lab_diags[i]; @@ -67,13 +69,14 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) 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"); + */ } } diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index a3f555bd0..dd167ba41 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -219,8 +219,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, attribs[PIdx::uz] = u[2]; attribs[PIdx::w ] = weight; - // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - if (WarpX::do_boosted_frame_diagnostic) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["xold"], x); @@ -504,8 +503,7 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) attribs[PIdx::uy] = u[1]; attribs[PIdx::uz] = u[2]; - // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - if (WarpX::do_boosted_frame_diagnostic) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); particle_tile.push_back_real(particle_comps["xold"], x); @@ -747,8 +745,7 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) attribs[PIdx::uz] = u[2]; // note - this will be slow on the GPU, need to revisit - // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - if (WarpX::do_boosted_frame_diagnostic) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); particle_tile.push_back_real(particle_comps["xold"], x); @@ -838,8 +835,9 @@ FieldGatherES (const amrex::Vector, const auto& particles = pti.GetArrayOfStructs(); int nstride = particles.dataShape().first; const long np = pti.numParticles(); - + std::cout<<"start 1 GetAttribs\n"; auto& attribs = pti.GetAttribs(); + std::cout<<"end 2 GetAttribs\n"; auto& Exp = attribs[PIdx::Ex]; auto& Eyp = attribs[PIdx::Ey]; #if AMREX_SPACEDIM == 3 diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index d1e25f3ad..6968f2a61 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -59,7 +59,6 @@ public: const amrex::Cuda::ManagedDeviceVector& y, const amrex::Cuda::ManagedDeviceVector& z); #endif - const std::array& GetAttribs () const { return GetStructOfArrays().GetRealData(); } diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 9bee0031a..66f0dfb5c 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -6,6 +6,8 @@ #include #include #include +#include +#include using namespace amrex; @@ -22,7 +24,9 @@ WarpXParIter::GetPosition (Cuda::ManagedDeviceVector& x, Cuda::ManagedDevi { amrex::ParIter<0,0,PIdx::nattribs>::GetPosition(x, z); #ifdef WARPX_RZ + std::cout<<"start 2 GetAttribs()\n"; const auto& attribs = GetAttribs(); + std::cout<<"stop 2 GetAttribs()\n"; const auto& theta = attribs[PIdx::theta]; y.resize(x.size()); for (unsigned int i=0 ; i < x.size() ; i++) { @@ -39,7 +43,9 @@ void WarpXParIter::SetPosition (const Cuda::ManagedDeviceVector& x, const Cuda::ManagedDeviceVector& y, const Cuda::ManagedDeviceVector& z) { #ifdef WARPX_RZ + std::cout<<"start 3 GetAttribs()\n"; auto& attribs = GetAttribs(); + std::cout<<"stop 3 GetAttribs()\n"; auto& theta = attribs[PIdx::theta]; Cuda::DeviceVector r(x.size()); for (unsigned int i=0 ; i < x.size() ; i++) { @@ -78,8 +84,7 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) particle_comps["theta"] = PIdx::theta; #endif - // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - if (WarpX::do_boosted_frame_diagnostic) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { particle_comps["xold"] = PIdx::nattribs; particle_comps["yold"] = PIdx::nattribs+1; @@ -232,8 +237,7 @@ WarpXParticleContainer::AddNParticles (int lev, p.pos(1) = z[i]; #endif - // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - if (WarpX::do_boosted_frame_diagnostic) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["xold"], x[i]); @@ -251,8 +255,7 @@ WarpXParticleContainer::AddNParticles (int lev, particle_tile.push_back_real(PIdx::uy, vy + ibegin, vy + iend); particle_tile.push_back_real(PIdx::uz, vz + ibegin, vz + iend); - // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - if (WarpX::do_boosted_frame_diagnostic) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); @@ -998,7 +1001,9 @@ WarpXParticleContainer::PushXES (Real dt) int nstride = particles.dataShape().first; const long np = pti.numParticles(); + std::cout<<"start 4 GetAttribs()\n"; auto& attribs = pti.GetAttribs(); + std::cout<<"stop 4 GetAttribs()\n"; auto& uxp = attribs[PIdx::ux]; auto& uyp = attribs[PIdx::uy]; auto& uzp = attribs[PIdx::uz]; @@ -1047,7 +1052,9 @@ WarpXParticleContainer::PushX (int lev, Real dt) // - positions are stored as an array of struct, in `ParticleType` ParticleType * AMREX_RESTRICT pstructs = &(pti.GetArrayOfStructs()[0]); // - momenta are stored as a struct of array, in `attribs` + std::cout<<"start 5 GetAttribs()\n"; auto& attribs = pti.GetAttribs(); + std::cout<<"stop 5 GetAttribs()\n"; Real* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); Real* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); Real* AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 7c026e105..e4c54cc05 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -160,7 +160,7 @@ WarpX::WarpX () current_injection_position = geom[0].ProbLo(moving_window_dir); } } - do_boosted_frame_particles = mypc->do_boosted_frame_diags; + do_boosted_frame_particles = mypc->doBoostedFrameDiags(); Efield_aux.resize(nlevs_max); Bfield_aux.resize(nlevs_max); -- cgit v1.2.3 From 449c1254dfd126b7b2b6291b05177cb22ed38f5c Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Sun, 12 May 2019 10:35:02 -0700 Subject: cleaning (remove print statements etc.) --- Source/Diagnostics/BoostedFrameDiagnostic.cpp | 5 +---- Source/Evolve/WarpXEvolveEM.cpp | 2 -- Source/Particles/MultiParticleContainer.H | 2 -- Source/Particles/MultiParticleContainer.cpp | 21 --------------------- Source/Particles/PhysicalParticleContainer.cpp | 2 -- Source/Particles/WarpXParticleContainer.cpp | 10 ---------- Source/WarpX.cpp | 7 +------ 7 files changed, 2 insertions(+), 47 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.cpp b/Source/Diagnostics/BoostedFrameDiagnostic.cpp index 413730f91..49038b9e2 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.cpp +++ b/Source/Diagnostics/BoostedFrameDiagnostic.cpp @@ -539,7 +539,6 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) } if (WarpX::do_boosted_frame_particles) { - // for (int j = 0; j < mypc.nSpecies(); ++j) { for (int j = 0; j < mypc.nSpeciesLabFrameDiags(); ++j) { int js = mypc.mapSpeciesLabDiags(j); std::string species_name = species_names[js]; @@ -645,8 +644,8 @@ writeLabFrameData(const MultiFab* cell_centered_data, &ncomp, &i_boost, &i_lab); } } - if (WarpX::do_boosted_frame_particles) { + if (WarpX::do_boosted_frame_particles) { mypc.GetLabFrameData(snapshots_[i].file_name, i_lab, boost_direction_, old_z_boost, snapshots_[i].current_z_boost, t_boost, snapshots_[i].t_lab, dt, particles_buffer_[i]); @@ -732,7 +731,6 @@ writeParticleDataHDF5(const WarpXParticleContainer::DiagnosticParticleData& pdat ParallelDescriptor::ReduceLongMax(old_np); // Write data here - Print()<<"particle_field_names.size()"<(particle_field_names.size()); ++k) { std::string field_path = species_name + "/" + particle_field_names[k]; @@ -869,7 +867,6 @@ LabSnapShot(Real t_lab_in, Real t_boost, Real zmin_lab_in, if (WarpX::do_boosted_frame_particles){ auto & mypc = WarpX::GetInstance().GetPartContainer(); const std::vector species_names = mypc.GetSpeciesNames(); - // for (int j = 0; j < mypc.nSpecies(); ++j) for (int j = 0; j < mypc.nSpeciesLabFrameDiags(); ++j) { int js = mypc.mapSpeciesLabDiags(j); diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index d88f178df..dab58f95b 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -247,9 +247,7 @@ WarpX::EvolveEM (int numsteps) } if (do_boosted_frame_diagnostic) { - std::cout<<"before myBFD->Flush"<Flush(geom[0]); - std::cout<<"after myBFD->Flush"<AddRealComp("xold"); - allcontainers[i]->AddRealComp("yold"); - allcontainers[i]->AddRealComp("zold"); - allcontainers[i]->AddRealComp("uxold"); - allcontainers[i]->AddRealComp("uyold"); - allcontainers[i]->AddRealComp("uzold"); - } - */ for (int i = 0; i < nspecies_lab_frame_diags; ++i) { int is = map_species_lab_diags[i]; @@ -69,14 +57,12 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) 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"); - */ } } @@ -500,10 +486,3 @@ MultiParticleContainer::doContinuousInjection() const } return warpx_do_continuous_injection; } - -// Set number of species for which lab frame data is dumped -// and maps their ID to MultiParticleContainer IDs. -//void -//MultiParticleContainer::setSpeciesLabFrameDiags() const -//{ -//} diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index dd167ba41..37c136a3d 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -835,9 +835,7 @@ FieldGatherES (const amrex::Vector, const auto& particles = pti.GetArrayOfStructs(); int nstride = particles.dataShape().first; const long np = pti.numParticles(); - std::cout<<"start 1 GetAttribs\n"; auto& attribs = pti.GetAttribs(); - std::cout<<"end 2 GetAttribs\n"; auto& Exp = attribs[PIdx::Ex]; auto& Eyp = attribs[PIdx::Ey]; #if AMREX_SPACEDIM == 3 diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 66f0dfb5c..0cf5c10b4 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -6,8 +6,6 @@ #include #include #include -#include -#include using namespace amrex; @@ -24,9 +22,7 @@ WarpXParIter::GetPosition (Cuda::ManagedDeviceVector& x, Cuda::ManagedDevi { amrex::ParIter<0,0,PIdx::nattribs>::GetPosition(x, z); #ifdef WARPX_RZ - std::cout<<"start 2 GetAttribs()\n"; const auto& attribs = GetAttribs(); - std::cout<<"stop 2 GetAttribs()\n"; const auto& theta = attribs[PIdx::theta]; y.resize(x.size()); for (unsigned int i=0 ; i < x.size() ; i++) { @@ -43,9 +39,7 @@ void WarpXParIter::SetPosition (const Cuda::ManagedDeviceVector& x, const Cuda::ManagedDeviceVector& y, const Cuda::ManagedDeviceVector& z) { #ifdef WARPX_RZ - std::cout<<"start 3 GetAttribs()\n"; auto& attribs = GetAttribs(); - std::cout<<"stop 3 GetAttribs()\n"; auto& theta = attribs[PIdx::theta]; Cuda::DeviceVector r(x.size()); for (unsigned int i=0 ; i < x.size() ; i++) { @@ -1001,9 +995,7 @@ WarpXParticleContainer::PushXES (Real dt) int nstride = particles.dataShape().first; const long np = pti.numParticles(); - std::cout<<"start 4 GetAttribs()\n"; auto& attribs = pti.GetAttribs(); - std::cout<<"stop 4 GetAttribs()\n"; auto& uxp = attribs[PIdx::ux]; auto& uyp = attribs[PIdx::uy]; auto& uzp = attribs[PIdx::uz]; @@ -1052,9 +1044,7 @@ WarpXParticleContainer::PushX (int lev, Real dt) // - positions are stored as an array of struct, in `ParticleType` ParticleType * AMREX_RESTRICT pstructs = &(pti.GetArrayOfStructs()[0]); // - momenta are stored as a struct of array, in `attribs` - std::cout<<"start 5 GetAttribs()\n"; auto& attribs = pti.GetAttribs(); - std::cout<<"stop 5 GetAttribs()\n"; Real* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); Real* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); Real* AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index e4c54cc05..6b1edaccd 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -56,12 +56,9 @@ int WarpX::num_mirrors = 0; int WarpX::sort_int = -1; +bool WarpX::do_boosted_frame_diagnostic = false; int WarpX::num_snapshots_lab = std::numeric_limits::lowest(); Real WarpX::dt_snapshots_lab = std::numeric_limits::lowest(); -// bool WarpX::do_boosted_frame_diagnostic = false; -// bool WarpX::do_boosted_frame_fields = true; -// bool WarpX::do_boosted_frame_particles = true; -bool WarpX::do_boosted_frame_diagnostic = false; bool WarpX::do_boosted_frame_fields = true; bool WarpX::do_boosted_frame_particles = true; @@ -329,8 +326,6 @@ WarpX::ReadParameters () pp.get("gamma_boost", gamma_boost); pp.query("do_boosted_frame_fields", do_boosted_frame_fields); - // pp.query("do_boosted_frame_particles", do_boosted_frame_particles); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_moving_window, "The moving window should be on if using the boosted frame diagnostic."); -- cgit v1.2.3 From 002cdf8fd7fb445055956f8de7e6b7bf678ab3ce Mon Sep 17 00:00:00 2001 From: Revathi Jambunathan Date: Tue, 14 May 2019 18:55:53 -0400 Subject: changes made after merging with lastest dev version --- Docs/source/building/building.rst | 2 +- Docs/source/running_cpp/parameters.rst | 56 ++++++++++-- Examples/Tests/Langmuir/inputs.multi.rt | 3 +- Python/pywarpx/picmi.py | 6 +- Source/Diagnostics/BoostedFrameDiagnostic.H | 9 ++ Source/Diagnostics/BoostedFrameDiagnostic.cpp | 72 ++++++++++++---- Source/Diagnostics/FieldIO.cpp | 29 ++++--- Source/Diagnostics/ParticleIO.cpp | 12 ++- Source/Diagnostics/WarpXIO.cpp | 4 +- Source/Evolve/WarpXEvolveEM.cpp | 48 +++++++++++ Source/FortranInterface/WarpX_f.H | 6 -- Source/FortranInterface/WarpX_picsar.F90 | 32 ------- Source/Laser/LaserParticleContainer.cpp | 5 +- Source/Particles/MultiParticleContainer.H | 12 ++- Source/Particles/MultiParticleContainer.cpp | 49 ++++++++--- Source/Particles/PhysicalParticleContainer.cpp | 45 ++++++++-- .../Particles/RigidInjectedParticleContainer.cpp | 2 +- Source/Particles/WarpXParticleContainer.H | 18 +++- Source/Particles/WarpXParticleContainer.cpp | 99 ++++++++++------------ Source/Utils/WarpXUtil.cpp | 3 +- Source/WarpX.H | 10 ++- Source/WarpX.cpp | 49 +++-------- 22 files changed, 373 insertions(+), 198 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Docs/source/building/building.rst b/Docs/source/building/building.rst index 01787545b..a9ea4e755 100644 --- a/Docs/source/building/building.rst +++ b/Docs/source/building/building.rst @@ -19,7 +19,7 @@ single directory (e.g. ``warpx_directory``): cd warpx_directory git clone --branch dev https://github.com/ECP-WarpX/WarpX.git git clone --branch master https://bitbucket.org/berkeleylab/picsar.git - git clone --branch developement https://github.com/AMReX-Codes/amrex.git + git clone --branch development https://github.com/AMReX-Codes/amrex.git Basic compilation ----------------- diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 076e39ce6..b22cf745c 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -28,6 +28,15 @@ Overall simulation parameters The direction of the Lorentz-transform for boosted-frame simulations (The direction ``y`` cannot be used in 2D simulations.) +* ``warpx.zmax_plasma_to_compute_max_step`` (`float`) optional + Can be useful when running in a boosted frame. If specified, automatically + calculates the number of iterations required in the boosted frame for the + lower `z` end of the simulation domain to reach + ``warpx.zmax_plasma_to_compute_max_step`` (typically the plasma end, + given in the lab frame). The value of ``max_step`` is overwritten, and + printed to standard output. Currently only works if the Lorentz boost and + the moving window are along the z direction. + * ``warpx.verbose`` (`0` or `1`) Controls how much information is printed to the terminal, when running WarpX. @@ -283,6 +292,25 @@ Particle initialization axes (4 particles in 2D, 6 particles in 3D). When `1`, particles are split along the diagonals (4 particles in 2D, 8 particles in 3D). +* ``.plot_species`` (`0` or `1` optional; default `1`) + Whether to plot particle quantities for this species. + +* ``.plot_vars`` (list of `strings` separated by spaces, optional) + List of particle quantities to write to `plotfiles`. By defaults, all + quantities are written to file. Choices are + * ``w`` for the particle weight, + * ``ux`` ``uy`` ``uz`` for the particle momentum, + * ``Ex`` ``Ey`` ``Ez`` for the electric field on particles, + * ``Bx`` ``By`` ``Bz`` for the magnetic field on particles. + The particle positions are always included. Use + ``.plot_vars = none`` to plot no particle data, except + particle position. + +* ``.do_boosted_frame_diags`` (`0` or `1` optional, default `1`) + Only used when ``warpx.do_boosted_frame_diagnostic=1``. When running in a + boosted frame, whether or not to plot back-transformed diagnostics for + this species. + * ``warpx.serialize_ics`` (`0 or 1`) Whether or not to use OpenMP threading for particle initialization. @@ -547,6 +575,20 @@ Numerics and algorithms same points in space) or a staggered grid (i.e. Yee grid ; different fields are defined at different points in space) +* ``warpx.do_subcycling`` (`0` or `1`; default: 0) + Whether or not to use sub-cycling. Different refinement levels have a + different cell size, which results in different Courant–Friedrichs–Lewy + (CFL) limits for the time step. By default, when using mesh refinement, + the same time step is used for all levels. This time step is + taken as the CFL limit of the finest level. Hence, for coarser + levels, the timestep is only a fraction of the CFL limit for this + level, which may lead to numerical artifacts. With sub-cycling, each level + evolves with its own time step, set to its own CFL limit. In practice, it + means that when level 0 performs one iteration, level 1 performs two + iterations. Currently, this option is only supported when + ``amr.max_level = 1``. More information can be found at + https://ieeexplore.ieee.org/document/8659392. + * ``psatd.nox``, ``psatd.noy``, ``pstad.noz`` (`integer`) optional (default `16` for all) The order of accuracy of the spatial derivatives, when using the code compiled with a PSATD solver. @@ -638,14 +680,18 @@ Diagnostics and output (This is done by averaging the field.) ``plot_coarsening_ratio`` should be an integer divisor of ``blocking_factor``. -* ``warpx.particle_plot_vars`` (`strings`, separated by spaces ; default: all) - Control which particle variables get written to the plot file. Choices are: - `w`, `ux`, `uy`, `uz`, `Ex`, `Ey`, `Ez`, `Bx`, `By`, and `Bz`. - The particle positions and ids are always included. - * ``amr.plot_file`` (`string`) Root for output file names. Supports sub-directories. Default `diags/plotfiles/plt` +* ``warpx.plot_J_field`` (`0` or `1` optional; default `1`) + Whether to plot the current density. + +* ``warpx.plot_E_field`` (`0` or `1` optional; default `1`) + Whether to plot the electric field. + +* ``warpx.plot_B_field`` (`0` or `1` optional; default `1`) + Whether to plot the magnetic field. + Checkpoints and restart ----------------------- WarpX supports checkpoints/restart via AMReX. diff --git a/Examples/Tests/Langmuir/inputs.multi.rt b/Examples/Tests/Langmuir/inputs.multi.rt index 334d4ae94..c9969d39a 100644 --- a/Examples/Tests/Langmuir/inputs.multi.rt +++ b/Examples/Tests/Langmuir/inputs.multi.rt @@ -13,7 +13,6 @@ amr.max_level = 0 amr.plot_int = 20 # How often to write plotfiles. "<= 0" means no plotfiles. warpx.plot_rho = 1 -warpx.particle_plot_vars = w ux Bz Ey # Geometry geometry.coord_sys = 0 # 0: Cartesian @@ -61,6 +60,7 @@ electrons.ymin = -20.e-6 electrons.ymax = 20.e-6 electrons.zmin = -20.e-6 electrons.zmax = 20.e-6 +electrons.plot_vars = w ux Bz Ey electrons.profile = constant electrons.density = 2.e24 # number of electrons per m^3 @@ -79,6 +79,7 @@ positrons.ymin = -20.e-6 positrons.ymax = 20.e-6 positrons.zmin = -20.e-6 positrons.zmax = 20.e-6 +positrons.plot_vars = w ux Bz Ey positrons.profile = constant positrons.density = 2.e24 # number of positrons per m^3 diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index a439205a7..f242f8589 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -113,9 +113,6 @@ class GaussianBunchDistribution(picmistandard.PICMI_GaussianBunchDistribution): species.uy = self.centroid_velocity[1] species.uz = self.centroid_velocity[2] - if self.fill_in: - species.do_continuous_injection = 1 - class UniformDistribution(picmistandard.PICMI_UniformDistribution): def initialize_inputs(self, species_number, layout, species): @@ -156,6 +153,9 @@ class UniformDistribution(picmistandard.PICMI_UniformDistribution): species.uy = self.directed_velocity[1] species.uz = self.directed_velocity[2] + if self.fill_in: + species.do_continuous_injection = 1 + class AnalyticDistribution(picmistandard.PICMI_AnalyticDistribution): def initialize_inputs(self, species_number, layout, species): diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.H b/Source/Diagnostics/BoostedFrameDiagnostic.H index e35d307a6..5b447689a 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.H +++ b/Source/Diagnostics/BoostedFrameDiagnostic.H @@ -74,10 +74,19 @@ class BoostedFrameDiagnostic { unsigned Nz_lab_; int boost_direction_; + // For back-transformed diagnostics of grid fields, data_buffer_[i] + // stores a buffer of the fields in the lab frame (in a MultiFab, i.e. + // with all box data etc.). When the buffer if full, dump to file. amrex::Vector > data_buffer_; + // particles_buffer_ is currently blind to refinement level. + // particles_buffer_[i][j] is a WarpXParticleContainer::DiagnosticParticleData where + // - i is the back-transformed snapshot number + // - j is the species number amrex::Vector > particles_buffer_; int num_buffer_ = 256; int max_box_size_ = 256; + // buff_counter_[i] is the number of z slices in data_buffer_[i] + // for snapshot number i. amrex::Vector buff_counter_; amrex::Vector snapshots_; diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.cpp b/Source/Diagnostics/BoostedFrameDiagnostic.cpp index 13972075d..11789579c 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.cpp +++ b/Source/Diagnostics/BoostedFrameDiagnostic.cpp @@ -504,6 +504,7 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) auto & mypc = WarpX::GetInstance().GetPartContainer(); const std::vector species_names = mypc.GetSpeciesNames(); + // Loop over BFD snapshots for (int i = 0; i < N_snapshots_; ++i) { int i_lab = (snapshots_[i].current_z_lab - snapshots_[i].zmin_lab) / dz_lab_; @@ -539,14 +540,20 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) } if (WarpX::do_boosted_frame_particles) { - for (int j = 0; j < mypc.nSpecies(); ++j) { + // Loop over species to be dumped to BFD + for (int j = 0; j < mypc.nSpeciesBoostedFrameDiags(); ++j) { + // Get species name + std::string species_name = + species_names[mypc.mapSpeciesBoostedFrameDiags(j)]; #ifdef WARPX_USE_HDF5 + // Dump species data writeParticleDataHDF5(particles_buffer_[i][j], snapshots_[i].file_name, - species_names[j]); + species_name); #else std::stringstream part_ss; - part_ss << snapshots_[i].file_name + "/" + species_names[j] + "/"; + part_ss << snapshots_[i].file_name + "/" + species_name + "/"; + // Dump species data writeParticleData(particles_buffer_[i][j], part_ss.str(), i_lab); #endif } @@ -575,21 +582,30 @@ writeLabFrameData(const MultiFab* cell_centered_data, const Real zhi_boost = domain_z_boost.hi(boost_direction_); const std::vector species_names = mypc.GetSpeciesNames(); - + + // Loop over snapshots for (int i = 0; i < N_snapshots_; ++i) { + + // Get updated z position of snapshot const Real old_z_boost = snapshots_[i].current_z_boost; snapshots_[i].updateCurrentZPositions(t_boost, inv_gamma_boost_, inv_beta_boost_); + // If snapshot out of the domain, nothing to do if ( (snapshots_[i].current_z_boost < zlo_boost) or (snapshots_[i].current_z_boost > zhi_boost) or (snapshots_[i].current_z_lab < snapshots_[i].zmin_lab) or (snapshots_[i].current_z_lab > snapshots_[i].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. int i_lab = (snapshots_[i].current_z_lab - snapshots_[i].zmin_lab) / dz_lab_; + // If buffer of snapshot i is empty... if (buff_counter_[i] == 0) { + // ... reset fields buffer data_buffer_[i] if (WarpX::do_boosted_frame_fields) { const int ncomp = cell_centered_data->nComp(); Box buff_box = geom.Domain(); @@ -600,13 +616,16 @@ writeLabFrameData(const MultiFab* cell_centered_data, DistributionMapping buff_dm(buff_ba); data_buffer_[i].reset( new MultiFab(buff_ba, buff_dm, ncomp, 0) ); } - if (WarpX::do_boosted_frame_particles) particles_buffer_[i].resize(mypc.nSpecies()); + // ... reset particle buffer particles_buffer_[i] + if (WarpX::do_boosted_frame_particles) + particles_buffer_[i].resize(mypc.nSpeciesBoostedFrameDiags()); } if (WarpX::do_boosted_frame_fields) { const int ncomp = cell_centered_data->nComp(); const int start_comp = 0; const bool interpolate = true; + // Get slice in the boosted frame std::unique_ptr slice = amrex::get_slice_data(boost_direction_, snapshots_[i].current_z_boost, *cell_centered_data, geom, @@ -620,19 +639,23 @@ writeLabFrameData(const MultiFab* cell_centered_data, &gamma_boost_, &beta_boost_); } + // Create a 2D box for the slice in the boosted frame Real dx = geom.CellSize(boost_direction_); int i_boost = (snapshots_[i].current_z_boost - geom.ProbLo(boost_direction_))/dx; Box slice_box = geom.Domain(); slice_box.setSmall(boost_direction_, i_boost); slice_box.setBig(boost_direction_, i_boost); + // Make it a BoxArray slice_ba BoxArray slice_ba(slice_box); slice_ba.maxSize(max_box_size_); + // Create MultiFab tmp on slice_ba with data from slice MultiFab tmp(slice_ba, data_buffer_[i]->DistributionMap(), ncomp, 0); tmp.copy(*slice, 0, 0, ncomp); #ifdef _OPENMP #pragma omp parallel #endif + // Copy data from MultiFab tmp to MultiDab data_buffer[i] for (MFIter mfi(tmp, true); mfi.isValid(); ++mfi) { const Box& tile_box = mfi.tilebox(); WRPX_COPY_SLICE(BL_TO_FORTRAN_BOX(tile_box), @@ -641,16 +664,18 @@ writeLabFrameData(const MultiFab* cell_centered_data, &ncomp, &i_boost, &i_lab); } } - + if (WarpX::do_boosted_frame_particles) { mypc.GetLabFrameData(snapshots_[i].file_name, i_lab, boost_direction_, old_z_boost, snapshots_[i].current_z_boost, t_boost, snapshots_[i].t_lab, dt, particles_buffer_[i]); + } ++buff_counter_[i]; + // If buffer full, write to disk. if (buff_counter_[i] == num_buffer_) { if (WarpX::do_boosted_frame_fields) { @@ -666,14 +691,21 @@ writeLabFrameData(const MultiFab* cell_centered_data, } if (WarpX::do_boosted_frame_particles) { - for (int j = 0; j < mypc.nSpecies(); ++j) { + // Loop over species to be dumped to BFD + for (int j = 0; j < mypc.nSpeciesBoostedFrameDiags(); ++j) { + // Get species name + const std::string species_name = species_names[mypc.mapSpeciesBoostedFrameDiags(j)]; #ifdef WARPX_USE_HDF5 + // Write data to disk (HDF5) writeParticleDataHDF5(particles_buffer_[i][j], snapshots_[i].file_name, - species_names[j]); + species_name); #else std::stringstream part_ss; - part_ss << snapshots_[i].file_name + "/" + species_names[j] + "/"; + + part_ss << snapshots_[i].file_name + "/" + species_name + "/"; + + // Write data to disk (custom) writeParticleData(particles_buffer_[i][j], part_ss.str(), i_lab); #endif } @@ -855,16 +887,19 @@ LabSnapShot(Real t_lab_in, Real t_boost, Real zmin_lab_in, ParallelDescriptor::Barrier(); - if (WarpX::do_boosted_frame_particles) - { + if (WarpX::do_boosted_frame_particles){ auto & mypc = WarpX::GetInstance().GetPartContainer(); const std::vector species_names = mypc.GetSpeciesNames(); - for (int j = 0; j < mypc.nSpecies(); ++j) + // Loop over species to be dumped to BFD + for (int j = 0; j < mypc.nSpeciesBoostedFrameDiags(); ++j) { - output_create_species_group(file_name, species_names[j]); + // Loop over species to be dumped to BFD + std::string species_name = + species_names[mypc.mapSpeciesBoostedFrameDiags(j)]; + output_create_species_group(file_name, species_name); for (int k = 0; k < static_cast(particle_field_names.size()); ++k) { - std::string field_path = species_names[j] + "/" + particle_field_names[k]; + std::string field_path = species_name + "/" + particle_field_names[k]; output_create_particle_field(file_name, field_path); } } @@ -884,11 +919,14 @@ LabSnapShot(Real t_lab_in, Real t_boost, Real zmin_lab_in, auto & mypc = WarpX::GetInstance().GetPartContainer(); const std::vector species_names = mypc.GetSpeciesNames(); - int nspecies = mypc.nSpecies(); const std::string particles_prefix = "particle"; - for(int i = 0; i < nspecies; ++i) { - const std::string fullpath = file_name + "/" + species_names[i]; + // Loop over species to be dumped to BFD + for(int i = 0; i < mypc.nSpeciesBoostedFrameDiags(); ++i) { + // Get species name + std::string species_name = + species_names[mypc.mapSpeciesBoostedFrameDiags(i)]; + const std::string fullpath = file_name + "/" + species_name; if (!UtilCreateDirectory(fullpath, 0755)) CreateDirectoryFailed(fullpath); } diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp index 209d8e9b4..ced1f8ea3 100644 --- a/Source/Diagnostics/FieldIO.cpp +++ b/Source/Diagnostics/FieldIO.cpp @@ -299,7 +299,10 @@ WarpX::AverageAndPackFields ( Vector& varnames, amrex::Vector& mf_avg, const int ngrow) const { // Count how many different fields should be written (ncomp) - const int ncomp = 3*3 + const int ncomp = 0 + + static_cast(plot_E_field)*3 + + static_cast(plot_B_field)*3 + + static_cast(plot_J_field)*3 + static_cast(plot_part_per_cell) + static_cast(plot_part_per_grid) + static_cast(plot_part_per_proc) @@ -321,15 +324,21 @@ WarpX::AverageAndPackFields ( Vector& varnames, // Go through the different fields, pack them into mf_avg[lev], // add the corresponding names to `varnames` and increment dcomp int dcomp = 0; - AverageAndPackVectorField(mf_avg[lev], current_fp[lev], dcomp, ngrow); - if(lev==0) for(auto name:{"jx","jy","jz"}) varnames.push_back(name); - dcomp += 3; - AverageAndPackVectorField(mf_avg[lev], Efield_aux[lev], dcomp, ngrow); - if(lev==0) for(auto name:{"Ex","Ey","Ez"}) varnames.push_back(name); - dcomp += 3; - AverageAndPackVectorField(mf_avg[lev], Bfield_aux[lev], dcomp, ngrow); - if(lev==0) for(auto name:{"Bx","By","Bz"}) varnames.push_back(name); - dcomp += 3; + if (plot_J_field) { + AverageAndPackVectorField(mf_avg[lev], current_fp[lev], dcomp, ngrow); + if(lev==0) for(auto name:{"jx","jy","jz"}) varnames.push_back(name); + dcomp += 3; + } + if (plot_E_field) { + AverageAndPackVectorField(mf_avg[lev], Efield_aux[lev], dcomp, ngrow); + if(lev==0) for(auto name:{"Ex","Ey","Ez"}) varnames.push_back(name); + dcomp += 3; + } + if (plot_B_field) { + AverageAndPackVectorField(mf_avg[lev], Bfield_aux[lev], dcomp, ngrow); + if(lev==0) for(auto name:{"Bx","By","Bz"}) varnames.push_back(name); + dcomp += 3; + } if (plot_part_per_cell) { diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index f15c084a0..258069a37 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -28,16 +28,20 @@ MultiParticleContainer::Checkpoint (const std::string& dir) const void MultiParticleContainer::WritePlotFile (const std::string& dir, - const Vector& real_flags, const Vector& real_names) const { Vector int_names; Vector int_flags; for (unsigned i = 0, n = species_names.size(); i < n; ++i) { - allcontainers[i]->WritePlotFile(dir, species_names[i], - real_flags, int_flags, - real_names, int_names); + auto& pc = allcontainers[i]; + if (pc->plot_species) { + // real_names contains a list of all particle attributes. + // pc->plot_flags is 1 or 0, whether quantity is dumped or not. + pc->WritePlotFile(dir, species_names[i], + pc->plot_flags, int_flags, + real_names, int_names); + } } } diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index 186a8d45e..1762a9b79 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -422,7 +422,7 @@ WarpX::GetCellCenteredData() { AverageAndPackVectorField( *cc[lev], Efield_aux[lev], dcomp, ng ); dcomp += 3; // then the magnetic field - AverageAndPackVectorField( *cc[lev], Efield_aux[lev], dcomp, ng ); + AverageAndPackVectorField( *cc[lev], Bfield_aux[lev], dcomp, ng ); dcomp += 3; // then the current density AverageAndPackVectorField( *cc[lev], current_fp[lev], dcomp, ng ); @@ -642,7 +642,7 @@ WarpX::WritePlotFile () const particle_varnames.push_back("uzold"); } - mypc->WritePlotFile(plotfilename, particle_plot_flags, particle_varnames); + mypc->WritePlotFile(plotfilename, particle_varnames); WriteJobInfo(plotfilename); diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 6541dd4be..1c7729a83 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -25,6 +25,10 @@ WarpX::EvolveEM (int numsteps) static int last_check_file_step = 0; static int last_insitu_step = 0; + if (do_compute_max_step_from_zmax){ + computeMaxStepBoostAccelerator(geom[0]); + } + int numsteps_max; if (numsteps < 0) { // Note that the default argument is numsteps = -1 numsteps_max = max_step; @@ -506,6 +510,50 @@ WarpX::ComputeDt () } } +/* \brief computes max_step for wakefield simulation in boosted frame. + * \param geom: Geometry object that contains simulation domain. + * + * max_step is set so that the simulation stop when the lower corner of the + * simulation box passes input parameter zmax_plasma_to_compute_max_step. + */ +void +WarpX::computeMaxStepBoostAccelerator(amrex::Geometry geom){ + // Sanity checks: can use zmax_plasma_to_compute_max_step only if + // the moving window and the boost are all in z direction. + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + WarpX::moving_window_dir == AMREX_SPACEDIM-1, + "Can use zmax_plasma_to_compute_max_step only if " + + "moving window along z. TODO: all directions."); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + maxLevel() == 0, + "Can use zmax_plasma_to_compute_max_step only if " + + "max level = 0."); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + (WarpX::boost_direction[0]-0)*(WarpX::boost_direction[0]-0) + + (WarpX::boost_direction[1]-0)*(WarpX::boost_direction[1]-0) + + (WarpX::boost_direction[2]-1)*(WarpX::boost_direction[2]-1) < 1.e-12, + "Can use zmax_plasma_to_compute_max_step only if " + + "warpx.boost_direction = z. TODO: all directions."); + + // Lower end of the simulation domain. All quantities are given in boosted + // frame except zmax_plasma_to_compute_max_step. + const Real zmin_domain_boost = geom.ProbLo(AMREX_SPACEDIM-1); + // End of the plasma: Transform input argument + // zmax_plasma_to_compute_max_step to boosted frame. + const Real len_plasma_boost = zmax_plasma_to_compute_max_step/gamma_boost; + // Plasma velocity + const Real v_plasma_boost = -beta_boost * PhysConst::c; + // Get time at which the lower end of the simulation domain passes the + // upper end of the plasma (in the z direction). + const Real interaction_time_boost = (len_plasma_boost-zmin_domain_boost)/ + (moving_window_v-v_plasma_boost); + // Divide by dt, and update value of max_step. + const int computed_max_step = interaction_time_boost/dt[0]; + max_step = computed_max_step; + Print()<<"max_step computed in computeMaxStepBoostAccelerator: " + < - !> @brief - !> Main subroutine for the particle pusher of positions - !> - !> @param[in] np number of super-particles - !> @param[in] xp,yp,zp particle position arrays - !> @param[in] uxp,uyp,uzp normalized momentum in each direction - !> @param[in] gaminv particle Lorentz factors - !> @param[in] dt time step - !> @param[in] particle_pusher_algo Particle pusher algorithm - subroutine warpx_particle_pusher_positions(np,xp,yp,zp,uxp,uyp,uzp, & - gaminv,dt) & - bind(C, name="warpx_particle_pusher_positions") - - INTEGER(c_long), INTENT(IN) :: np - REAL(amrex_real),INTENT(INOUT) :: gaminv(np) - REAL(amrex_real),INTENT(INOUT) :: xp(np),yp(np),zp(np) - REAL(amrex_real),INTENT(INOUT) :: uxp(np),uyp(np),uzp(np) - REAL(amrex_real),INTENT(IN) :: dt - - CALL pxr_set_gamma(np,uxp,uyp,uzp,gaminv) - - !!!! --- push particle species positions a time step -#if (AMREX_SPACEDIM == 3) || (defined WARPX_RZ) - CALL pxr_pushxyz(np,xp,yp,zp,uxp,uyp,uzp,gaminv,dt) -#elif (AMREX_SPACEDIM == 2) - CALL pxr_pushxz(np,xp,zp,uxp,uzp,gaminv,dt) -#endif - - end subroutine warpx_particle_pusher_positions - ! _________________________________________________________________ !> !> @brief diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index 8a2590e5e..bda5dcd86 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -25,8 +25,9 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, { charge = 1.0; mass = std::numeric_limits::max(); - - ParmParse pp(laser_name); + do_boosted_frame_diags = 0; + + ParmParse pp(laser_name); // Parse the type of laser profile and set the corresponding flag `profile` std::string laser_type_s; diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 0c5e49c04..513f30b99 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -130,7 +130,6 @@ public: void Checkpoint (const std::string& dir) const; void WritePlotFile( const std::string& dir, - const amrex::Vector& real_flags, const amrex::Vector& real_names) const; void Restart (const std::string& dir); @@ -156,6 +155,10 @@ public: int nSpecies() const {return nspecies;} + int nSpeciesBoostedFrameDiags() const {return nspecies_boosted_frame_diags;} + int mapSpeciesBoostedFrameDiags(int i) const {return map_species_boosted_frame_diags[i];} + int doBoostedFrameDiags() const {return do_boosted_frame_diags;} + int nSpeciesDepositOnMainGrid () const { int r = 0; for (int i : deposit_on_main_grid) { @@ -215,6 +218,13 @@ private: void ReadParameters (); + // Number of species dumped in BoostedFrameDiagnostics + int nspecies_boosted_frame_diags = 0; + // map_species_boosted_frame_diags[i] is the species ID in + // MultiParticleContainer for 0 map_species_boosted_frame_diags; + int do_boosted_frame_diags = 0; + // runtime parameters int nlasers = 0; int nspecies = 1; // physical particles only. nspecies+nlasers == allcontainers.size(). diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 440906348..6d618c096 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -31,16 +31,31 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) pc_tmp.reset(new PhysicalParticleContainer(amr_core)); - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + // Compute the number of species for which lab-frame data is dumped + // nspecies_lab_frame_diags, and map their ID to MultiParticleContainer + // particle IDs in map_species_lab_diags. + map_species_boosted_frame_diags.resize(nspecies); + nspecies_boosted_frame_diags = 0; + for (int i=0; ido_boosted_frame_diags){ + map_species_boosted_frame_diags[nspecies_boosted_frame_diags] = i; + do_boosted_frame_diags = 1; + nspecies_boosted_frame_diags += 1; + } + } + + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - for (int i = 0; i < nspecies + nlasers; ++i) + for (int i = 0; i < nspecies_boosted_frame_diags; ++i) { - allcontainers[i]->AddRealComp("xold"); - allcontainers[i]->AddRealComp("yold"); - allcontainers[i]->AddRealComp("zold"); - allcontainers[i]->AddRealComp("uxold"); - allcontainers[i]->AddRealComp("uyold"); - allcontainers[i]->AddRealComp("uzold"); + 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"); @@ -376,13 +391,25 @@ MultiParticleContainer BL_PROFILE("MultiParticleContainer::GetLabFrameData"); - for (int i = 0; i < nspecies; ++i){ - WarpXParticleContainer* pc = allcontainers[i].get(); + // Loop over particle species + for (int i = 0; i < nspecies_boosted_frame_diags; ++i){ + int isp = map_species_boosted_frame_diags[i]; + WarpXParticleContainer* pc = allcontainers[isp].get(); WarpXParticleContainer::DiagnosticParticles diagnostic_particles; pc->GetParticleSlice(direction, z_old, z_new, t_boost, t_lab, dt, diagnostic_particles); - + // Here, diagnostic_particles[lev][index] is a WarpXParticleContainer::DiagnosticParticleData + // where "lev" is the AMR level and "index" is a [grid index][tile index] pair. + + // Loop over AMR levels for (int lev = 0; lev <= pc->finestLevel(); ++lev){ + // Loop over [grid index][tile index] pairs + // and Fills parts[species number i] with particle data from all grids and + // tiles in diagnostic_particles. parts contains particles from all + // AMR levels indistinctly. for (auto it = diagnostic_particles[lev].begin(); it != diagnostic_particles[lev].end(); ++it){ + // it->first is the [grid index][tile index] key + // it->second is the corresponding + // WarpXParticleContainer::DiagnosticParticleData value parts[i].GetRealData(DiagIdx::w).insert( parts[i].GetRealData(DiagIdx::w ).end(), it->second.GetRealData(DiagIdx::w ).begin(), it->second.GetRealData(DiagIdx::w ).end()); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index a1d4e1319..8ef4c13a9 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -82,6 +82,40 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_splitting", do_splitting); pp.query("split_type", split_type); pp.query("do_continuous_injection", do_continuous_injection); + // Whether to plot back-transformed (lab-frame) diagnostics + // for this species. + pp.query("do_boosted_frame_diags", do_boosted_frame_diags); + + pp.query("plot_species", plot_species); + int do_user_plot_vars; + do_user_plot_vars = pp.queryarr("plot_vars", plot_vars); + if (not do_user_plot_vars){ + // By default, all particle variables are dumped to plotfiles, + // including {x,y,z,ux,uy,uz}old variables when running in a + // boosted frame + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){ + plot_flags.resize(PIdx::nattribs + 6, 1); + } else { + plot_flags.resize(PIdx::nattribs, 1); + } + } else { + // Set plot_flag to 0 for all attribs + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){ + plot_flags.resize(PIdx::nattribs + 6, 0); + } else { + plot_flags.resize(PIdx::nattribs, 0); + } + // If not none, set plot_flags values to 1 for elements in plot_vars. + if (plot_vars[0] != "none"){ + for (const auto& var : plot_vars){ + // Return error if var not in PIdx. + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + ParticleStringNames::to_index.count(var), + "plot_vars argument not in ParticleStringNames"); + plot_flags[ParticleStringNames::to_index.at(var)] = 1; + } + } + } } PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) @@ -185,7 +219,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, attribs[PIdx::uz] = u[2]; attribs[PIdx::w ] = weight; - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["xold"], x); @@ -469,7 +503,7 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) attribs[PIdx::uy] = u[1]; attribs[PIdx::uz] = u[2]; - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); particle_tile.push_back_real(particle_comps["xold"], x); @@ -711,7 +745,7 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) attribs[PIdx::uz] = u[2]; // note - this will be slow on the GPU, need to revisit - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); particle_tile.push_back_real(particle_comps["xold"], x); @@ -801,7 +835,6 @@ FieldGatherES (const amrex::Vector, const auto& particles = pti.GetArrayOfStructs(); int nstride = particles.dataShape().first; const long np = pti.numParticles(); - auto& attribs = pti.GetAttribs(); auto& Exp = attribs[PIdx::Ex]; auto& Eyp = attribs[PIdx::Ey]; @@ -1686,7 +1719,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, auto& Bzp = attribs[PIdx::Bz]; const long np = pti.numParticles(); - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& xpold = pti.GetAttribs(particle_comps["xold"]); auto& ypold = pti.GetAttribs(particle_comps["yold"]); @@ -1831,6 +1864,8 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real // Note the the slice should always move in the negative boost direction. AMREX_ALWAYS_ASSERT(z_new < z_old); + AMREX_ALWAYS_ASSERT(do_boosted_frame_diags == 1); + const int nlevs = std::max(0, finestLevel()+1); // we figure out a box for coarse-grained rejection. If the RealBox corresponding to a diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index a5acca281..fd1b2dfb5 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -225,7 +225,7 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, auto& Bzp = attribs[PIdx::Bz]; const long np = pti.numParticles(); - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& xpold = pti.GetAttribs(particle_comps["xold"]); auto& ypold = pti.GetAttribs(particle_comps["yold"]); diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 6fa02b824..89244237b 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -59,7 +59,6 @@ public: const amrex::Cuda::ManagedDeviceVector& y, const amrex::Cuda::ManagedDeviceVector& z); #endif - const std::array& GetAttribs () const { return GetStructOfArrays().GetRealData(); } @@ -85,7 +84,13 @@ class WarpXParticleContainer public: friend MultiParticleContainer; + // amrex::StructOfArrays with DiagIdx::nattribs amrex::Real components + // and 0 int components for the particle data. using DiagnosticParticleData = amrex::StructOfArrays; + // DiagnosticParticles is a vector, with one element per MR level. + // DiagnosticParticles[lev] is typically a key-value pair where the key is + // a pair [grid_index, tile_index], and the value is the corresponding + // DiagnosticParticleData (see above) on this tile. using DiagnosticParticles = amrex::Vector, DiagnosticParticleData> >; WarpXParticleContainer (amrex::AmrCore* amr_core, int ispecies); @@ -265,6 +270,8 @@ protected: // support all features allowed by direct injection. int do_continuous_injection = 0; + int do_boosted_frame_diags = 1; + amrex::Vector local_rho; amrex::Vector local_jx; amrex::Vector local_jy; @@ -272,6 +279,15 @@ protected: amrex::Vector > m_xp, m_yp, m_zp, m_giv; + // Whether to dump particle quantities. + // If true, particle position is always dumped. + int plot_species = 1; + // For all particle attribs (execept position), whether or not + // to dump to file. + amrex::Vector plot_flags; + // list of names of attributes to dump. + amrex::Vector plot_vars; + 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 89e21df1c..32816ba5f 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -78,7 +78,7 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) particle_comps["theta"] = PIdx::theta; #endif - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { particle_comps["xold"] = PIdx::nattribs; particle_comps["yold"] = PIdx::nattribs+1; @@ -86,9 +86,9 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) 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 @@ -174,7 +174,7 @@ WarpXParticleContainer::AddNParticles (int lev, int n, const Real* x, const Real* y, const Real* z, const Real* vx, const Real* vy, const Real* vz, int nattr, const Real* attr, int uniqueparticles, int id) -{ +{ BL_ASSERT(nattr == 1); const Real* weight = attr; @@ -230,15 +230,15 @@ WarpXParticleContainer::AddNParticles (int lev, #endif p.pos(1) = z[i]; #endif - - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["xold"], x[i]); particle_tile.push_back_real(particle_comps["yold"], y[i]); - particle_tile.push_back_real(particle_comps["zold"], z[i]); + particle_tile.push_back_real(particle_comps["zold"], z[i]); } - + particle_tile.push_back(p); } @@ -249,14 +249,14 @@ WarpXParticleContainer::AddNParticles (int lev, particle_tile.push_back_real(PIdx::uy, vy + ibegin, vy + iend); particle_tile.push_back_real(PIdx::uz, vz + ibegin, vz + iend); - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); particle_tile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend); - particle_tile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); + particle_tile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); } - + for (int comp = PIdx::uz+1; comp < PIdx::nattribs; ++comp) { #ifdef WARPX_RZ @@ -327,7 +327,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, jx_ptr = local_jx[thread_num].dataPtr(); jy_ptr = local_jy[thread_num].dataPtr(); jz_ptr = local_jz[thread_num].dataPtr(); - + local_jx[thread_num].setVal(0.0); local_jy[thread_num].setVal(0.0); local_jz[thread_num].setVal(0.0); @@ -426,7 +426,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, BL_PROFILE_VAR_STOP(blp_pxr_cd); -#ifndef AMREX_USE_GPU +#ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); jx[pti].atomicAdd(local_jx[thread_num], tbx, tbx, 0, 0, 1); @@ -477,7 +477,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, auto jyntot = local_jy[thread_num].length(); auto jzntot = local_jz[thread_num].length(); #endif - + long ncrse = np - np_current; BL_PROFILE_VAR_START(blp_pxr_cd); if (j_is_nodal) { @@ -650,7 +650,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, #endif BL_PROFILE_VAR_STOP(blp_pxr_chd); -#ifndef AMREX_USE_GPU +#ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); (*rhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp, 1); @@ -710,7 +710,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, #endif BL_PROFILE_VAR_STOP(blp_pxr_chd); -#ifndef AMREX_USE_GPU +#ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); (*crhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp, 1); @@ -1025,8 +1025,6 @@ void WarpXParticleContainer::PushX (int lev, Real dt) { BL_PROFILE("WPC::PushX()"); - BL_PROFILE_VAR_NS("WPC::PushX::Copy", blp_copy); - BL_PROFILE_VAR_NS("WPC:PushX::Push", blp_pxr_pp); if (do_not_push) return; @@ -1036,47 +1034,38 @@ WarpXParticleContainer::PushX (int lev, Real dt) #pragma omp parallel #endif { - Cuda::ManagedDeviceVector xp, yp, zp, giv; for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { Real wt = amrex::second(); - auto& attribs = pti.GetAttribs(); - - auto& uxp = attribs[PIdx::ux]; - auto& uyp = attribs[PIdx::uy]; - auto& uzp = attribs[PIdx::uz]; - - const long np = pti.numParticles(); - - giv.resize(np); - - // - // copy data from particle container to temp arrays - // - BL_PROFILE_VAR_START(blp_copy); - pti.GetPosition(xp, yp, zp); - BL_PROFILE_VAR_STOP(blp_copy); - // // Particle Push // - BL_PROFILE_VAR_START(blp_pxr_pp); - warpx_particle_pusher_positions(&np, - xp.dataPtr(), - yp.dataPtr(), - zp.dataPtr(), - uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), - giv.dataPtr(), &dt); - BL_PROFILE_VAR_STOP(blp_pxr_pp); - - // - // copy particle data back - // - BL_PROFILE_VAR_START(blp_copy); - pti.SetPosition(xp, yp, zp); - BL_PROFILE_VAR_STOP(blp_copy); + // Extract pointers to particle position and momenta, for this particle tile + // - positions are stored as an array of struct, in `ParticleType` + ParticleType * AMREX_RESTRICT pstructs = &(pti.GetArrayOfStructs()[0]); + // - momenta are stored as a struct of array, in `attribs` + auto& attribs = pti.GetAttribs(); + Real* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + Real* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + Real* AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); + // Loop over the particles and update the position + const long np = pti.numParticles(); + const Real inv_c2 = 1./(PhysConst::c*PhysConst::c); + amrex::ParallelFor( np, + [=] AMREX_GPU_DEVICE (long i) { + // Compute inverse Lorentz factor + const Real inv_gamma = 1./std::sqrt( + 1. + (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])*inv_c2); + // Update positions over one time step + pstructs[i].pos(0) += ux[i] * inv_gamma * dt; +#if (AMREX_SPACEDIM == 3) || (defined WARPX_RZ) // RZ pushes particles in 3D + pstructs[i].pos(1) += uy[i] * inv_gamma * dt; +#endif + pstructs[i].pos(2) += uz[i] * inv_gamma * dt; + } + ); if (cost) { const Box& tbx = pti.tilebox(); @@ -1092,15 +1081,15 @@ WarpXParticleContainer::PushX (int lev, Real dt) } // This function is called in Redistribute, just after locate -void -WarpXParticleContainer::particlePostLocate(ParticleType& p, +void +WarpXParticleContainer::particlePostLocate(ParticleType& p, const ParticleLocData& pld, const int lev) { // Tag particle if goes to higher level. // It will be split later in the loop - if (pld.m_lev == lev+1 - and p.m_idata.id != NoSplitParticleID + if (pld.m_lev == lev+1 + and p.m_idata.id != NoSplitParticleID and p.m_idata.id >= 0) { p.m_idata.id = DoSplitParticleID; diff --git a/Source/Utils/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp index 4ec7ebb51..a5ea6d75a 100644 --- a/Source/Utils/WarpXUtil.cpp +++ b/Source/Utils/WarpXUtil.cpp @@ -118,8 +118,7 @@ void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, amrex::Real zmax) const int lo_ind = bx.loVect()[1]; #endif // Check if box intersect with [zmin, zmax] - if ( (zmin>zmin_box && zmin<=zmax_box) || - (zmax>zmin_box && zmax<=zmax_box) ){ + if ( (zmax>zmin_box && zmin<=zmax_box) ){ Array4 arr = mf[mfi].array(); // Set field to 0 between zmin and zmax ParallelFor(bx, diff --git a/Source/WarpX.H b/Source/WarpX.H index 44e84033f..71fc54af2 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -108,6 +108,8 @@ public: static amrex::Real gamma_boost; static amrex::Real beta_boost; static amrex::Vector boost_direction; + static amrex::Real zmax_plasma_to_compute_max_step; + static int do_compute_max_step_from_zmax; static bool do_dynamic_scheduling; static bool refine_plasma; @@ -152,6 +154,8 @@ public: void applyMirrors(amrex::Real time); void ComputeDt (); + // Compute max_step automatically for simulations in a boosted frame. + void computeMaxStepBoostAccelerator(amrex::Geometry geom); int MoveWindow (bool move_j); void UpdatePlasmaInjectionPosition (amrex::Real dt); @@ -515,6 +519,9 @@ private: bool dump_plotfiles = true; bool dump_openpmd = false; #endif + bool plot_E_field = true; + bool plot_B_field = true; + bool plot_J_field = true; bool plot_part_per_cell = true; bool plot_part_per_grid = false; bool plot_part_per_proc = false; @@ -543,9 +550,6 @@ private: bool is_synchronized = true; - amrex::Vector particle_plot_vars; - amrex::Vector particle_plot_flags; - #ifdef WARPX_USE_PSATD // Store fields in real space on the dual grid (i.e. the grid for the FFT push of the fields) // This includes data for the FFT guard cells (between FFT groups) diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index a3a24897a..6b1edaccd 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -32,6 +32,8 @@ int WarpX::moving_window_dir = -1; Real WarpX::gamma_boost = 1.; Real WarpX::beta_boost = 0.; Vector WarpX::boost_direction = {0,0,0}; +int WarpX::do_compute_max_step_from_zmax = 0; +Real WarpX::zmax_plasma_to_compute_max_step = 0.; long WarpX::current_deposition_algo = 3; long WarpX::charge_deposition_algo = 0; @@ -115,7 +117,7 @@ WarpX::ResetInstance () { delete m_instance; m_instance = nullptr; -} +} WarpX::WarpX () { @@ -155,6 +157,7 @@ WarpX::WarpX () current_injection_position = geom[0].ProbLo(moving_window_dir); } } + do_boosted_frame_particles = mypc->doBoostedFrameDiags(); Efield_aux.resize(nlevs_max); Bfield_aux.resize(nlevs_max); @@ -272,6 +275,12 @@ WarpX::ReadParameters () ReadBoostedFrameParameters(gamma_boost, beta_boost, boost_direction); + // pp.query returns 1 if argument zmax_plasma_to_compute_max_step is + // specified by the user, 0 otherwise. + do_compute_max_step_from_zmax = + pp.query("zmax_plasma_to_compute_max_step", + zmax_plasma_to_compute_max_step); + pp.queryarr("B_external", B_external); pp.query("do_moving_window", do_moving_window); @@ -317,8 +326,6 @@ WarpX::ReadParameters () pp.get("gamma_boost", gamma_boost); pp.query("do_boosted_frame_fields", do_boosted_frame_fields); - pp.query("do_boosted_frame_particles", do_boosted_frame_particles); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_moving_window, "The moving window should be on if using the boosted frame diagnostic."); @@ -371,6 +378,9 @@ WarpX::ReadParameters () if (ParallelDescriptor::NProcs() == 1) { plot_proc_number = false; } + pp.query("plot_E_field" , plot_E_field); + pp.query("plot_B_field" , plot_B_field); + pp.query("plot_J_field" , plot_J_field); pp.query("plot_part_per_cell", plot_part_per_cell); pp.query("plot_part_per_grid", plot_part_per_grid); pp.query("plot_part_per_proc", plot_part_per_proc); @@ -428,39 +438,6 @@ WarpX::ReadParameters () fine_tag_hi = RealVect{hi}; } - // select which particle comps to write - { - pp.queryarr("particle_plot_vars", particle_plot_vars); - - if (particle_plot_vars.size() == 0) - { - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) - { - particle_plot_flags.resize(PIdx::nattribs + 6, 1); - } - else - { - particle_plot_flags.resize(PIdx::nattribs, 1); - } - } - else - { - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) - { - particle_plot_flags.resize(PIdx::nattribs + 6, 0); - } - else - { - particle_plot_flags.resize(PIdx::nattribs, 0); - } - - for (const auto& var : particle_plot_vars) - { - particle_plot_flags[ParticleStringNames::to_index.at(var)] = 1; - } - } - } - pp.query("load_balance_int", load_balance_int); pp.query("load_balance_with_sfc", load_balance_with_sfc); pp.query("load_balance_knapsack_factor", load_balance_knapsack_factor); -- cgit v1.2.3 From 33a2926da6fa6f32ff9ec0f40108d99e7a1492b9 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Tue, 21 May 2019 16:15:18 -0700 Subject: Added new file for pusher kernel --- Source/Particles/Make.package | 2 ++ Source/Particles/Pusher/GetAndSetPosition.H | 40 +++++++++++++++++++++++++++++ Source/Particles/Pusher/Make.package | 3 +++ Source/Particles/WarpXParticleContainer.cpp | 2 ++ 4 files changed, 47 insertions(+) create mode 100644 Source/Particles/Pusher/GetAndSetPosition.H create mode 100644 Source/Particles/Pusher/Make.package (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package index acbfe3390..f33095738 100644 --- a/Source/Particles/Make.package +++ b/Source/Particles/Make.package @@ -10,5 +10,7 @@ CEXE_headers += WarpXParticleContainer.H CEXE_headers += RigidInjectedParticleContainer.H CEXE_headers += PhysicalParticleContainer.H +include $(WARPX_HOME)/Source/Particles/Pusher/Make.package + INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles diff --git a/Source/Particles/Pusher/GetAndSetPosition.H b/Source/Particles/Pusher/GetAndSetPosition.H new file mode 100644 index 000000000..8eaa99dd8 --- /dev/null +++ b/Source/Particles/Pusher/GetAndSetPosition.H @@ -0,0 +1,40 @@ +#ifndef WARPX_PARTICLES_PUSHER_GETANDSETPOSITION_H_ +#define WARPX_PARTICLES_PUSHER_GETANDSETPOSITION_H_ + +#include +#include +#include + +#ifndef WARPX_RZ + +void +GetPosition( + amrex::Real& x, amrex::Real& y, amrex::Real& z, + const WarpXParticleContainer::ParticleType& p) +{ +#if (AMREX_SPACEDIM==3) + x = p.pos(0); + y = p.pos(1); + z = p.pos(2); +#else + x = p.pos(0); + z = p.pos(1); +#endif +} + +# else // if WARPX_RZ is True + +void +GetCartesianPositionFromCylindrical( + amrex::Real& x, amrex::Real& y, amrex::Real& z, + const WarpXParticleContainer::ParticleType& p, const amrex::Real theta) +{ + const amrex::Real r = p.pos(0); + x = r*std::cos(theta); + y = r*std::sin(theta); + z = p.pos(1); +} + +#endif // WARPX_RZ + +#endif // WARPX_PARTICLES_PUSHER_GETANDSETPOSITION_H_ diff --git a/Source/Particles/Pusher/Make.package b/Source/Particles/Pusher/Make.package new file mode 100644 index 000000000..deb0f5260 --- /dev/null +++ b/Source/Particles/Pusher/Make.package @@ -0,0 +1,3 @@ +CEXE_headers += GetAndSetPosition.H +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/Pusher +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Pusher diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 0cf5c10b4..332760ebf 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -7,6 +7,8 @@ #include #include +#include + using namespace amrex; int WarpXParticleContainer::do_not_push = 0; -- cgit v1.2.3 From b78301705e3802e164b32bd09ededc9139baabda Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 22 May 2019 10:28:22 -0700 Subject: Refactor particle push --- Source/Particles/Pusher/GetAndSetPosition.H | 44 ++++++++++++++++++++++++++--- Source/Particles/Pusher/Make.package | 1 + Source/Particles/Pusher/UpdatePosition.H | 28 ++++++++++++++++++ Source/Particles/WarpXParticleContainer.cpp | 32 +++++++++++++-------- 4 files changed, 89 insertions(+), 16 deletions(-) create mode 100644 Source/Particles/Pusher/UpdatePosition.H (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/Pusher/GetAndSetPosition.H b/Source/Particles/Pusher/GetAndSetPosition.H index 8eaa99dd8..a08000860 100644 --- a/Source/Particles/Pusher/GetAndSetPosition.H +++ b/Source/Particles/Pusher/GetAndSetPosition.H @@ -7,8 +7,10 @@ #ifndef WARPX_RZ -void -GetPosition( +/* \brief Extract the particle's coordinates from the ParticleType struct `p`, + * and stores them in the variables `x`, `y`, `z`. */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void GetPosition( amrex::Real& x, amrex::Real& y, amrex::Real& z, const WarpXParticleContainer::ParticleType& p) { @@ -18,14 +20,35 @@ GetPosition( z = p.pos(2); #else x = p.pos(0); + y = std::numeric_limits::quiet_NaN(); z = p.pos(1); #endif } +/* \brief Set the particle's coordinates in the ParticleType struct `p`, + * from their values in the variables `x`, `y`, `z`. */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void SetPosition( + WarpXParticleContainer::ParticleType& p, + const amrex::Real x, const amrex::Real y, const amrex::Real z) +{ +#if (AMREX_SPACEDIM==3) + p.pos(0) = x; + p.pos(1) = y; + p.pos(2) = z; +#else + p.pos(0) = x; + p.pos(1) = z; +#endif +} + # else // if WARPX_RZ is True -void -GetCartesianPositionFromCylindrical( +/* \brief Extract the particle's coordinates from `theta` and the attributes + * of the ParticleType struct `p` (which contains the radius), + * and store them in the variables `x`, `y`, `z` */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void GetCartesianPositionFromCylindrical( amrex::Real& x, amrex::Real& y, amrex::Real& z, const WarpXParticleContainer::ParticleType& p, const amrex::Real theta) { @@ -35,6 +58,19 @@ GetCartesianPositionFromCylindrical( z = p.pos(1); } +/* \brief Set the particle's cylindrical coordinates by setting `theta` + * and the attributes of the ParticleType struct `p` (which stores the radius), + * from the values of `x`, `y`, `z` */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void SetCylindricalPositionFromCartesian( + WarpXParticleContainer::ParticleType& p, amrex::Real& theta, + const amrex::Real x, const amrex::Real y, const amrex::Real z ) +{ + theta = std::atan2(y, x); + p.pos(0) = std::sqrt(x*x + y*y); + p.pos(1) = z; +} + #endif // WARPX_RZ #endif // WARPX_PARTICLES_PUSHER_GETANDSETPOSITION_H_ diff --git a/Source/Particles/Pusher/Make.package b/Source/Particles/Pusher/Make.package index deb0f5260..8c8e77905 100644 --- a/Source/Particles/Pusher/Make.package +++ b/Source/Particles/Pusher/Make.package @@ -1,3 +1,4 @@ CEXE_headers += GetAndSetPosition.H +CEXE_headers += UpdatePosition.H INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/Pusher VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Pusher diff --git a/Source/Particles/Pusher/UpdatePosition.H b/Source/Particles/Pusher/UpdatePosition.H new file mode 100644 index 000000000..1442ac28e --- /dev/null +++ b/Source/Particles/Pusher/UpdatePosition.H @@ -0,0 +1,28 @@ +#ifndef WARPX_PARTICLES_PUSHER_UPDATEPOSITION_H_ +#define WARPX_PARTICLES_PUSHER_UPDATEPOSITION_H_ + +#include + +/* \brief Push the particle's positions over one timestep, + * given the value of its momenta `ux`, `uy`, `uz` */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void UpdatePosition( + amrex::Real& x, amrex::Real& y, amrex::Real& z, + const amrex::Real ux, const amrex::Real uy, const amrex::Real uz, + const amrex::Real dt ) +{ + + constexpr amrex::Real inv_c2 = 1./(PhysConst::c*PhysConst::c); + + // Compute inverse Lorentz factor + const amrex::Real inv_gamma = 1./std::sqrt(1. + (ux*ux + uy*uy + uz*uz)*inv_c2); + // Update positions over one time step + x += ux * inv_gamma * dt; +#if (AMREX_SPACEDIM == 3) || (defined WARPX_RZ) // RZ pushes particles in 3D + y += uy * inv_gamma * dt; +#endif + z += uz * inv_gamma * dt; + +} + +#endif // WARPX_PARTICLES_PUSHER_UPDATEPOSITION_H_ diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 332760ebf..e04d6aa69 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -7,7 +7,9 @@ #include #include +// Import low-level single-particle kernels #include +#include using namespace amrex; @@ -1050,20 +1052,26 @@ WarpXParticleContainer::PushX (int lev, Real dt) Real* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); Real* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); Real* AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); - // Loop over the particles and update the position - const long np = pti.numParticles(); - const Real inv_c2 = 1./(PhysConst::c*PhysConst::c); - amrex::ParallelFor( np, +#ifdef WARPX_RZ + Real* AMREX_RESTRICT theta = attribs[PIdx::theta].dataPtr(); +#endif + // Loop over the particles and update their position + amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - // Compute inverse Lorentz factor - const Real inv_gamma = 1./std::sqrt( - 1. + (ux[i]*ux[i] + uy[i]*uy[i] + uz[i]*uz[i])*inv_c2); - // Update positions over one time step - pstructs[i].pos(0) += ux[i] * inv_gamma * dt; -#if (AMREX_SPACEDIM == 3) || (defined WARPX_RZ) // RZ pushes particles in 3D - pstructs[i].pos(1) += uy[i] * inv_gamma * dt; + ParticleType& p = pstructs[i]; // Particle object that gets updated + Real x, y, z; // Temporary variables +#ifndef WARPX_RZ + GetPosition( x, y, z, p ); // Initialize x, y, z +#else + GetCartesianPositionFromCylindrical( x, y, z, p, theta ); +#endif + // Even for RZ, the particles are pushed in 3D Cartesian + UpdatePosition( x, y, z, ux[i], uy[i], uz[i], dt); +#ifndef WARPX_RZ + SetPosition( p, x, y, z ); // Update the object p +#else + SetCylindricalPositionFromCartesian( p, theta, x, y, z ); #endif - pstructs[i].pos(2) += uz[i] * inv_gamma * dt; } ); -- cgit v1.2.3 From 8ccd9bb67fa5fe100e59d999f5698d738e75d93c Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 22 May 2019 15:23:06 -0700 Subject: Implement style comments --- Source/Particles/Pusher/UpdatePosition.H | 2 ++ Source/Particles/WarpXParticleContainer.cpp | 8 +++----- 2 files changed, 5 insertions(+), 5 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/Pusher/UpdatePosition.H b/Source/Particles/Pusher/UpdatePosition.H index 1442ac28e..0a4f579f4 100644 --- a/Source/Particles/Pusher/UpdatePosition.H +++ b/Source/Particles/Pusher/UpdatePosition.H @@ -1,6 +1,8 @@ #ifndef WARPX_PARTICLES_PUSHER_UPDATEPOSITION_H_ #define WARPX_PARTICLES_PUSHER_UPDATEPOSITION_H_ +#include +#include #include /* \brief Push the particle's positions over one timestep, diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index e04d6aa69..fa7fda233 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -1062,14 +1062,12 @@ WarpXParticleContainer::PushX (int lev, Real dt) Real x, y, z; // Temporary variables #ifndef WARPX_RZ GetPosition( x, y, z, p ); // Initialize x, y, z -#else - GetCartesianPositionFromCylindrical( x, y, z, p, theta ); -#endif - // Even for RZ, the particles are pushed in 3D Cartesian UpdatePosition( x, y, z, ux[i], uy[i], uz[i], dt); -#ifndef WARPX_RZ SetPosition( p, x, y, z ); // Update the object p #else + // For WARPX_RZ, the particles are still pushed in 3D Cartesian + GetCartesianPositionFromCylindrical( x, y, z, p, theta ); + UpdatePosition( x, y, z, ux[i], uy[i], uz[i], dt); SetCylindricalPositionFromCartesian( p, theta, x, y, z ); #endif } -- cgit v1.2.3 From c394ff331f484fb04d0266f82a0ca41124af6f98 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 22 May 2019 16:11:26 -0700 Subject: Fix compilation error for cylindrical --- Source/Particles/WarpXParticleContainer.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index fa7fda233..c7ffe956b 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -1066,9 +1066,9 @@ WarpXParticleContainer::PushX (int lev, Real dt) SetPosition( p, x, y, z ); // Update the object p #else // For WARPX_RZ, the particles are still pushed in 3D Cartesian - GetCartesianPositionFromCylindrical( x, y, z, p, theta ); + GetCartesianPositionFromCylindrical( x, y, z, p, theta[i] ); UpdatePosition( x, y, z, ux[i], uy[i], uz[i], dt); - SetCylindricalPositionFromCartesian( p, theta, x, y, z ); + SetCylindricalPositionFromCartesian( p, theta[i], x, y, z ); #endif } ); -- cgit v1.2.3 From 5308cbefa8fbde30bed8a88943d58ec646731d9e Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Thu, 23 May 2019 15:27:17 -0700 Subject: fix a bunch of unused variable / parameter shadowing warnings --- Source/Diagnostics/BoostedFrameDiagnostic.H | 2 +- Source/Diagnostics/BoostedFrameDiagnostic.cpp | 5 --- Source/Diagnostics/FieldIO.cpp | 1 - Source/Evolve/WarpXEvolveEM.cpp | 28 +++++++------- Source/FieldSolver/WarpXPushFieldsEM.cpp | 48 ++++++++++++------------ Source/FortranInterface/WarpX_picsar.F90 | 7 +++- Source/Laser/LaserParticleContainer.cpp | 6 +-- Source/Parallelization/WarpXComm.cpp | 52 +++++++++++++------------- Source/Particles/PhysicalParticleContainer.cpp | 3 -- Source/Particles/WarpXParticleContainer.cpp | 44 +++++++++++----------- Source/Utils/WarpXMovingWindow.cpp | 6 +-- Source/WarpX.cpp | 1 - 12 files changed, 96 insertions(+), 107 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.H b/Source/Diagnostics/BoostedFrameDiagnostic.H index 55a124c51..3a09259b0 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.H +++ b/Source/Diagnostics/BoostedFrameDiagnostic.H @@ -40,8 +40,8 @@ class BoostedFrameDiagnostic { amrex::Real current_z_lab; amrex::Real current_z_boost; - std::vector mesh_field_names_; int ncomp_to_dump_; + std::vector mesh_field_names_; int file_num; int initial_i; diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.cpp b/Source/Diagnostics/BoostedFrameDiagnostic.cpp index f62369091..57872790f 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.cpp +++ b/Source/Diagnostics/BoostedFrameDiagnostic.cpp @@ -497,8 +497,6 @@ LorentzTransformZ(MultiFab& data, Real gamma_boost, Real beta_boost, int ncomp) // Transform the transverse E and B fields. Note that ez and bz are not // changed by the tranform. Real e_lab, b_lab, j_lab, r_lab; - int i0 = 0; - int i4 = 4; e_lab = gamma_boost * (arr(i, j, k, 0) + beta_boost*clight*arr(i, j, k, 4)); b_lab = gamma_boost * (arr(i, j, k, 4) + beta_boost*arr(i, j, k, 0)/clight); @@ -718,7 +716,6 @@ writeLabFrameData(const MultiFab* cell_centered_data, if (buff_counter_[i] == 0) { // ... reset fields buffer data_buffer_[i] if (WarpX::do_boosted_frame_fields) { - const int ncomp = cell_centered_data->nComp(); Box buff_box = geom.Domain(); buff_box.setSmall(boost_direction_, i_lab - num_buffer_ + 1); buff_box.setBig(boost_direction_, i_lab); @@ -726,7 +723,6 @@ writeLabFrameData(const MultiFab* cell_centered_data, buff_ba.maxSize(max_box_size_); DistributionMapping buff_dm(buff_ba); data_buffer_[i].reset( new MultiFab(buff_ba, buff_dm, ncomp_to_dump, 0) ); - // data_buffer_[i].reset( new MultiFab(buff_ba, buff_dm, ncomp, 0) ); } // ... reset particle buffer particles_buffer_[i] if (WarpX::do_boosted_frame_particles) @@ -961,7 +957,6 @@ LabSnapShot(Real t_lab_in, Real t_boost, RealBox prob_domain_lab, my_bfd(bfd) { Real zmin_lab = prob_domain_lab_.lo(AMREX_SPACEDIM-1); - Real zmax_lab = prob_domain_lab_.hi(AMREX_SPACEDIM-1); current_z_lab = 0.0; current_z_boost = 0.0; updateCurrentZPositions(t_boost, my_bfd.inv_gamma_boost_, my_bfd.inv_beta_boost_); diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp index ced1f8ea3..e3d44d1fc 100644 --- a/Source/Diagnostics/FieldIO.cpp +++ b/Source/Diagnostics/FieldIO.cpp @@ -660,7 +660,6 @@ getInterpolatedScalar( interpolated_F->setVal(0.); // Loop through the boxes and interpolate the values from the _cp data - const int use_limiter = 0; #ifdef _OPEMP #pragma omp parallel #endif diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index dab58f95b..ad7c7d840 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -514,7 +514,7 @@ WarpX::ComputeDt () * simulation box passes input parameter zmax_plasma_to_compute_max_step. */ void -WarpX::computeMaxStepBoostAccelerator(amrex::Geometry geom){ +WarpX::computeMaxStepBoostAccelerator(amrex::Geometry a_geom){ // Sanity checks: can use zmax_plasma_to_compute_max_step only if // the moving window and the boost are all in z direction. AMREX_ALWAYS_ASSERT_WITH_MESSAGE( @@ -534,7 +534,7 @@ WarpX::computeMaxStepBoostAccelerator(amrex::Geometry geom){ // Lower end of the simulation domain. All quantities are given in boosted // frame except zmax_plasma_to_compute_max_step. - const Real zmin_domain_boost = geom.ProbLo(AMREX_SPACEDIM-1); + const Real zmin_domain_boost = a_geom.ProbLo(AMREX_SPACEDIM-1); // End of the plasma: Transform input argument // zmax_plasma_to_compute_max_step to boosted frame. const Real len_plasma_boost = zmax_plasma_to_compute_max_step/gamma_boost; @@ -591,19 +591,19 @@ WarpX::applyMirrors(Real time){ NullifyMF(Bz, lev, z_min, z_max); if (lev>0){ // Get coarse patch field MultiFabs - MultiFab& Ex = *Efield_cp[lev][0].get(); - MultiFab& Ey = *Efield_cp[lev][1].get(); - MultiFab& Ez = *Efield_cp[lev][2].get(); - MultiFab& Bx = *Bfield_cp[lev][0].get(); - MultiFab& By = *Bfield_cp[lev][1].get(); - MultiFab& Bz = *Bfield_cp[lev][2].get(); + MultiFab& cEx = *Efield_cp[lev][0].get(); + MultiFab& cEy = *Efield_cp[lev][1].get(); + MultiFab& cEz = *Efield_cp[lev][2].get(); + MultiFab& cBx = *Bfield_cp[lev][0].get(); + MultiFab& cBy = *Bfield_cp[lev][1].get(); + MultiFab& cBz = *Bfield_cp[lev][2].get(); // Set each field to zero between z_min and z_max - NullifyMF(Ex, lev, z_min, z_max); - NullifyMF(Ey, lev, z_min, z_max); - NullifyMF(Ez, lev, z_min, z_max); - NullifyMF(Bx, lev, z_min, z_max); - NullifyMF(By, lev, z_min, z_max); - NullifyMF(Bz, lev, z_min, z_max); + NullifyMF(cEx, lev, z_min, z_max); + NullifyMF(cEy, lev, z_min, z_max); + NullifyMF(cEz, lev, z_min, z_max); + NullifyMF(cBx, lev, z_min, z_max); + NullifyMF(cBy, lev, z_min, z_max); + NullifyMF(cBz, lev, z_min, z_max); } } } diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index a66d14980..c53e13f8f 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -17,30 +17,30 @@ using namespace amrex; void -WarpX::EvolveB (Real dt) +WarpX::EvolveB (Real a_dt) { for (int lev = 0; lev <= finest_level; ++lev) { - EvolveB(lev, dt); + EvolveB(lev, a_dt); } } void -WarpX::EvolveB (int lev, Real dt) +WarpX::EvolveB (int lev, Real a_dt) { BL_PROFILE("WarpX::EvolveB()"); - EvolveB(lev, PatchType::fine, dt); + EvolveB(lev, PatchType::fine, a_dt); if (lev > 0) { - EvolveB(lev, PatchType::coarse, dt); + EvolveB(lev, PatchType::coarse, a_dt); } } void -WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) +WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) { const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const std::array& dx = WarpX::CellSize(patch_level); - Real dtsdx = dt/dx[0], dtsdy = dt/dx[1], dtsdz = dt/dx[2]; + Real dtsdx = a_dt/dx[0], dtsdy = a_dt/dx[1], dtsdz = a_dt/dx[2]; MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz; if (patch_type == PatchType::fine) @@ -164,30 +164,30 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) } void -WarpX::EvolveE (Real dt) +WarpX::EvolveE (Real a_dt) { for (int lev = 0; lev <= finest_level; ++lev) { - EvolveE(lev, dt); + EvolveE(lev, a_dt); } } void -WarpX::EvolveE (int lev, Real dt) +WarpX::EvolveE (int lev, Real a_dt) { BL_PROFILE("WarpX::EvolveE()"); - EvolveE(lev, PatchType::fine, dt); + EvolveE(lev, PatchType::fine, a_dt); if (lev > 0) { - EvolveE(lev, PatchType::coarse, dt); + EvolveE(lev, PatchType::coarse, a_dt); } } void -WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) +WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) { - const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * dt; - const Real c2dt = (PhysConst::c*PhysConst::c) * dt; + const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * a_dt; + const Real c2dt = (PhysConst::c*PhysConst::c) * a_dt; int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const std::array& dx = WarpX::CellSize(patch_level); @@ -358,27 +358,27 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) } void -WarpX::EvolveF (Real dt, DtType dt_type) +WarpX::EvolveF (Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; for (int lev = 0; lev <= finest_level; ++lev) { - EvolveF(lev, dt, dt_type); + EvolveF(lev, a_dt, a_dt_type); } } void -WarpX::EvolveF (int lev, Real dt, DtType dt_type) +WarpX::EvolveF (int lev, Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; - EvolveF(lev, PatchType::fine, dt, dt_type); - if (lev > 0) EvolveF(lev, PatchType::coarse, dt, dt_type); + EvolveF(lev, PatchType::fine, a_dt, a_dt_type); + if (lev > 0) EvolveF(lev, PatchType::coarse, a_dt, a_dt_type); } void -WarpX::EvolveF (int lev, PatchType patch_type, Real dt, DtType dt_type) +WarpX::EvolveF (int lev, PatchType patch_type, Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; @@ -388,7 +388,7 @@ WarpX::EvolveF (int lev, PatchType patch_type, Real dt, DtType dt_type) int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const auto& dx = WarpX::CellSize(patch_level); - const std::array dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]}; + const std::array dtsdx {a_dt/dx[0], a_dt/dx[1], a_dt/dx[2]}; MultiFab *Ex, *Ey, *Ez, *rho, *F; if (patch_type == PatchType::fine) @@ -408,12 +408,12 @@ WarpX::EvolveF (int lev, PatchType patch_type, Real dt, DtType dt_type) F = F_cp[lev].get(); } - const int rhocomp = (dt_type == DtType::FirstHalf) ? 0 : 1; + const int rhocomp = (a_dt_type == DtType::FirstHalf) ? 0 : 1; MultiFab src(rho->boxArray(), rho->DistributionMap(), 1, 0); ComputeDivE(src, 0, {Ex,Ey,Ez}, dx); MultiFab::Saxpy(src, -mu_c2, *rho, rhocomp, 0, 1, 0); - MultiFab::Saxpy(*F, dt, src, 0, 0, 1, 0); + MultiFab::Saxpy(*F, a_dt, src, 0, 0, 1, 0); if (do_pml && pml[lev]->ok()) { diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90 index a4cc99926..c17e8861b 100644 --- a/Source/FortranInterface/WarpX_picsar.F90 +++ b/Source/FortranInterface/WarpX_picsar.F90 @@ -265,8 +265,10 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n real(amrex_real), intent(IN OUT):: rho(*) real(amrex_real), intent(IN) :: rmin, dr +#ifdef WARPX_RZ integer(c_long) :: type_rz_depose = 1 - +#endif + ! Compute the number of valid cells and guard cells integer(c_long) :: rho_nvalid(AMREX_SPACEDIM), rho_nguards(AMREX_SPACEDIM) rho_nvalid = rho_ntot - 2*rho_ng @@ -385,8 +387,9 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n real(amrex_real), intent(IN OUT):: jx(*), jy(*), jz(*) real(amrex_real), intent(IN) :: rmin, dr +#ifdef WARPX_RZ integer(c_long) :: type_rz_depose = 1 - +#endif ! Compute the number of valid cells and guard cells integer(c_long) :: jx_nvalid(AMREX_SPACEDIM), jy_nvalid(AMREX_SPACEDIM), jz_nvalid(AMREX_SPACEDIM), & jx_nguards(AMREX_SPACEDIM), jy_nguards(AMREX_SPACEDIM), jz_nguards(AMREX_SPACEDIM) diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index 52d506bb8..2f964b6f3 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -78,14 +78,14 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, parser.define(field_function); parser.registerVariables({"X","Y","t"}); - ParmParse pp("my_constants"); + ParmParse ppc("my_constants"); std::set symbols = parser.symbols(); symbols.erase("X"); symbols.erase("Y"); symbols.erase("t"); // after removing variables, we are left with constants for (auto it = symbols.begin(); it != symbols.end(); ) { Real v; - if (pp.query(it->c_str(), v)) { + if (ppc.query(it->c_str(), v)) { parser.setConstant(*it, v); it = symbols.erase(it); } else { @@ -429,8 +429,6 @@ LaserParticleContainer::Evolve (int lev, { Real wt = amrex::second(); - const Box& box = pti.validbox(); - auto& attribs = pti.GetAttribs(); auto& wp = attribs[PIdx::w ]; diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 348b31c8b..1cf347420 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -79,7 +79,7 @@ WarpX::UpdateAuxilaryData () MultiFab::Subtract(dBz, *Bfield_cp[lev][2], 0, 0, 1, ng); const Real* dx = Geom(lev-1).CellSize(); - const int ref_ratio = refRatio(lev-1)[0]; + const int rr = refRatio(lev-1)[0]; #ifdef _OPENMP #pragma omp parallel #endif @@ -89,7 +89,7 @@ WarpX::UpdateAuxilaryData () { Box ccbx = mfi.fabbox(); ccbx.enclosedCells(); - ccbx.coarsen(ref_ratio).refine(ref_ratio); // so that ccbx is coarsenable + ccbx.coarsen(rr).refine(rr); // so that ccbx is coarsenable const FArrayBox& cxfab = dBx[mfi]; const FArrayBox& cyfab = dBy[mfi]; @@ -106,18 +106,18 @@ WarpX::UpdateAuxilaryData () BL_TO_FORTRAN_ANYD(cxfab), BL_TO_FORTRAN_ANYD(cyfab), BL_TO_FORTRAN_ANYD(czfab), - dx, &ref_ratio,&use_limiter); + dx, &rr,&use_limiter); #else amrex_interp_div_free_bfield(ccbx.loVect(), ccbx.hiVect(), BL_TO_FORTRAN_ANYD(bfab[0]), BL_TO_FORTRAN_ANYD(bfab[2]), BL_TO_FORTRAN_ANYD(cxfab), BL_TO_FORTRAN_ANYD(czfab), - dx, &ref_ratio,&use_limiter); + dx, &rr,&use_limiter); amrex_interp_cc_bfield(ccbx.loVect(), ccbx.hiVect(), BL_TO_FORTRAN_ANYD(bfab[1]), BL_TO_FORTRAN_ANYD(cyfab), - &ref_ratio,&use_limiter); + &rr,&use_limiter); #endif for (int idim = 0; idim < 3; ++idim) @@ -153,7 +153,7 @@ WarpX::UpdateAuxilaryData () MultiFab::Subtract(dEy, *Efield_cp[lev][1], 0, 0, 1, ng); MultiFab::Subtract(dEz, *Efield_cp[lev][2], 0, 0, 1, ng); - const int ref_ratio = refRatio(lev-1)[0]; + const int rr = refRatio(lev-1)[0]; #ifdef _OPEMP #pragma omp parallel #endif @@ -163,7 +163,7 @@ WarpX::UpdateAuxilaryData () { Box ccbx = mfi.fabbox(); ccbx.enclosedCells(); - ccbx.coarsen(ref_ratio).refine(ref_ratio); // so that ccbx is coarsenable + ccbx.coarsen(rr).refine(rr); // so that ccbx is coarsenable const FArrayBox& cxfab = dEx[mfi]; const FArrayBox& cyfab = dEy[mfi]; @@ -180,18 +180,18 @@ WarpX::UpdateAuxilaryData () BL_TO_FORTRAN_ANYD(cxfab), BL_TO_FORTRAN_ANYD(cyfab), BL_TO_FORTRAN_ANYD(czfab), - &ref_ratio,&use_limiter); + &rr,&use_limiter); #else amrex_interp_efield(ccbx.loVect(), ccbx.hiVect(), BL_TO_FORTRAN_ANYD(efab[0]), BL_TO_FORTRAN_ANYD(efab[2]), BL_TO_FORTRAN_ANYD(cxfab), BL_TO_FORTRAN_ANYD(czfab), - &ref_ratio,&use_limiter); + &rr,&use_limiter); amrex_interp_nd_efield(ccbx.loVect(), ccbx.hiVect(), BL_TO_FORTRAN_ANYD(efab[1]), BL_TO_FORTRAN_ANYD(cyfab), - &ref_ratio); + &rr); #endif for (int idim = 0; idim < 3; ++idim) @@ -364,7 +364,7 @@ WarpX::SyncCurrent () current_cp[lev][1]->setVal(0.0); current_cp[lev][2]->setVal(0.0); - const IntVect& ref_ratio = refRatio(lev-1); + const IntVect& rr = refRatio(lev-1); std::array fine { current_fp[lev][0].get(), current_fp[lev][1].get(), @@ -372,7 +372,7 @@ WarpX::SyncCurrent () std::array< MultiFab*,3> crse { current_cp[lev][0].get(), current_cp[lev][1].get(), current_cp[lev][2].get() }; - SyncCurrent(fine, crse, ref_ratio[0]); + SyncCurrent(fine, crse, rr[0]); } Vector,3> > j_fp(finest_level+1); @@ -524,10 +524,10 @@ WarpX::SyncCurrent () void WarpX::SyncCurrent (const std::array& fine, const std::array< amrex::MultiFab*,3>& crse, - int ref_ratio) + int rr) { - BL_ASSERT(ref_ratio == 2); - const IntVect& ng = (fine[0]->nGrowVect() + 1) /ref_ratio; + BL_ASSERT(rr == 2); + const IntVect& ng = (fine[0]->nGrowVect() + 1) /rr; #ifdef _OPEMP #pragma omp parallel @@ -539,7 +539,7 @@ WarpX::SyncCurrent (const std::array& fine, for (MFIter mfi(*crse[idim],true); mfi.isValid(); ++mfi) { const Box& bx = mfi.growntilebox(ng); - Box fbx = amrex::grow(amrex::refine(bx,ref_ratio),1); + Box fbx = amrex::grow(amrex::refine(bx,rr),1); ffab.resize(fbx); fbx &= (*fine[idim])[mfi].box(); ffab.setVal(0.0); @@ -564,8 +564,8 @@ WarpX::SyncRho (amrex::Vector >& rhof, for (int lev = 1; lev <= finest_level; ++lev) { rhoc[lev]->setVal(0.0); - const IntVect& ref_ratio = refRatio(lev-1); - SyncRho(*rhof[lev], *rhoc[lev], ref_ratio[0]); + const IntVect& rr = refRatio(lev-1); + SyncRho(*rhof[lev], *rhoc[lev], rr[0]); } Vector > rho_f_g(finest_level+1); @@ -673,10 +673,10 @@ WarpX::SyncRho (amrex::Vector >& rhof, * averaging the values of the charge density of the fine patch (on the same level). */ void -WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int ref_ratio) +WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int rr) { - BL_ASSERT(ref_ratio == 2); - const IntVect& ng = (fine.nGrowVect()+1)/ref_ratio; + BL_ASSERT(rr == 2); + const IntVect& ng = (fine.nGrowVect()+1)/rr; const int nc = fine.nComp(); #ifdef _OPEMP @@ -687,7 +687,7 @@ WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int ref_ratio) for (MFIter mfi(crse,true); mfi.isValid(); ++mfi) { const Box& bx = mfi.growntilebox(ng); - Box fbx = amrex::grow(amrex::refine(bx,ref_ratio),1); + Box fbx = amrex::grow(amrex::refine(bx,rr),1); ffab.resize(fbx, nc); fbx &= fine[mfi].box(); ffab.setVal(0.0); @@ -710,7 +710,7 @@ WarpX::RestrictCurrentFromFineToCoarsePatch (int lev) current_cp[lev][1]->setVal(0.0); current_cp[lev][2]->setVal(0.0); - const IntVect& ref_ratio = refRatio(lev-1); + const IntVect& rr = refRatio(lev-1); std::array fine { current_fp[lev][0].get(), current_fp[lev][1].get(), @@ -718,7 +718,7 @@ WarpX::RestrictCurrentFromFineToCoarsePatch (int lev) std::array< MultiFab*,3> crse { current_cp[lev][0].get(), current_cp[lev][1].get(), current_cp[lev][2].get() }; - SyncCurrent(fine, crse, ref_ratio[0]); + SyncCurrent(fine, crse, rr[0]); } void @@ -824,8 +824,8 @@ WarpX::RestrictRhoFromFineToCoarsePatch (int lev) { if (rho_fp[lev]) { rho_cp[lev]->setVal(0.0); - const IntVect& ref_ratio = refRatio(lev-1); - SyncRho(*rho_fp[lev], *rho_cp[lev], ref_ratio[0]); + const IntVect& rr = refRatio(lev-1); + SyncRho(*rho_fp[lev], *rho_cp[lev], rr[0]); } } diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 37c136a3d..212084e64 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1024,9 +1024,6 @@ PhysicalParticleContainer::FieldGather (int lev, { const std::array& dx = WarpX::CellSize(lev); - // WarpX assumes the same number of guard cells for Ex, Ey, Ez, Bx, By, Bz - long ng = Ex.nGrow(); - BL_ASSERT(OnSameGrids(lev,Ex)); MultiFab* cost = WarpX::getCosts(lev); diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index c7ffe956b..9791eee80 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -237,10 +237,10 @@ WarpXParticleContainer::AddNParticles (int lev, if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - particle_tile.push_back_real(particle_comps["xold"], x[i]); - particle_tile.push_back_real(particle_comps["yold"], y[i]); - particle_tile.push_back_real(particle_comps["zold"], z[i]); + auto& ptile = DefineAndReturnParticleTile(0, 0, 0); + 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); @@ -255,10 +255,10 @@ WarpXParticleContainer::AddNParticles (int lev, if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - particle_tile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); - particle_tile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend); - particle_tile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); + auto& ptile = DefineAndReturnParticleTile(0, 0, 0); + 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) @@ -737,7 +737,6 @@ WarpXParticleContainer::DepositCharge (Vector >& rho, const auto& gm = m_gdb->Geom(lev); const auto& ba = m_gdb->ParticleBoxArray(lev); - const auto& dm = m_gdb->DistributionMap(lev); const Real* dx = gm.CellSize(); const Real* plo = gm.ProbLo(); @@ -807,36 +806,36 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) #ifdef _OPENMP #pragma omp parallel -#endif { +#endif Cuda::ManagedDeviceVector xp, yp, zp; - FArrayBox local_rho; +#ifdef _OPENMP + FArrayBox rho_loc; +#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { - const Box& box = pti.validbox(); - auto& wp = pti.GetAttribs(PIdx::w); const long np = pti.numParticles(); pti.GetPosition(xp, yp, zp); - const std::array& xyzmin_tile = WarpX::LowerCorner(pti.tilebox(), lev); - const std::array& xyzmin_grid = WarpX::LowerCorner(box, lev); - // Data on the grid Real* data_ptr; FArrayBox& rhofab = (*rho)[pti]; #ifdef _OPENMP + const std::array& xyzmin_tile = WarpX::LowerCorner(pti.tilebox(), lev); Box tile_box = convert(pti.tilebox(), IntVect::TheUnitVector()); const std::array& xyzmin = xyzmin_tile; tile_box.grow(ng); - local_rho.resize(tile_box); - local_rho = 0.0; - data_ptr = local_rho.dataPtr(); - auto rholen = local_rho.length(); + rho_loc.resize(tile_box); + rho_loc = 0.0; + data_ptr = rho_loc.dataPtr(); + auto rholen = rho_loc.length(); #else + const Box& box = pti.validbox(); + const std::array& xyzmin_grid = WarpX::LowerCorner(box, lev); const std::array& xyzmin = xyzmin_grid; data_ptr = rhofab.dataPtr(); auto rholen = rhofab.length(); @@ -874,10 +873,9 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) #endif #ifdef _OPENMP - rhofab.atomicAdd(local_rho); -#endif + rhofab.atomicAdd(rho_loc); } - +#endif } if (!local) rho->SumBoundary(gm.periodicity()); diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp index 5f9e36bf6..a0ab1f26f 100644 --- a/Source/Utils/WarpXMovingWindow.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp @@ -5,7 +5,7 @@ using namespace amrex; void -WarpX::UpdatePlasmaInjectionPosition (Real dt) +WarpX::UpdatePlasmaInjectionPosition (Real a_dt) { int dir = moving_window_dir; // Continuously inject plasma in new cells (by default only on level 0) @@ -14,12 +14,12 @@ WarpX::UpdatePlasmaInjectionPosition (Real dt) // call to this function, and injection position needs to be updated current_injection_position -= WarpX::beta_boost * #if ( AMREX_SPACEDIM == 3 ) - WarpX::boost_direction[dir] * PhysConst::c * dt; + WarpX::boost_direction[dir] * PhysConst::c * a_dt; #elif ( AMREX_SPACEDIM == 2 ) // In 2D, dir=0 corresponds to x and dir=1 corresponds to z // This needs to be converted in order to index `boost_direction` // which has 3 components, for both 2D and 3D simulations. - WarpX::boost_direction[2*dir] * PhysConst::c * dt; + WarpX::boost_direction[2*dir] * PhysConst::c * a_dt; #endif } } diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 6b6752bf1..3d7f7dcc5 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -1001,7 +1001,6 @@ WarpX::ComputeDivE (MultiFab& divE, int dcomp, void WarpX::BuildBufferMasks () { - int ngbuffer = std::max(n_field_gather_buffer, n_current_deposition_buffer); for (int lev = 1; lev <= maxLevel(); ++lev) { for (int ipass = 0; ipass < 2; ++ipass) -- cgit v1.2.3 From e7bc3c34fc838c4ab6ee845de2e3f6354c3cbdec Mon Sep 17 00:00:00 2001 From: Revathi Jambunathan Date: Thu, 6 Jun 2019 17:23:22 -0400 Subject: adding 2D cufft plans --- .../SpectralSolver/SpectralFieldData.cpp | 21 +++++++++++++-------- Source/Parallelization/WarpXComm.cpp | 5 ----- Source/Particles/PhysicalParticleContainer.cpp | 2 -- Source/Particles/WarpXParticleContainer.cpp | 2 -- 4 files changed, 13 insertions(+), 17 deletions(-) (limited to 'Source/Particles/WarpXParticleContainer.cpp') diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index cb26b19b7..6eeb266e7 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -55,7 +55,7 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, #ifdef AMREX_USE_GPU // Add cuFFT-specific code // Creating 3D plan for real to complex -- double precision - + // Assuming CUDA is used for programming GPU cufftResult result; #if (AMREX_SPACEDIM == 3) result = cufftPlan3d( &forward_plan[mfi], fft_size[2], @@ -63,21 +63,26 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, if ( result != CUFFT_SUCCESS ) { amrex::Print() << " cufftplan3d forward failed! \n"; } -#else - // Add 2D cuFFT-spacific code for D2Z - // Note that D2Z is inherently forward plan -#endif -#if (AMREX_SPACEDIM == 3) result = cufftPlan3d( &backward_plan[mfi], fft_size[2], fft_size[1], fft_size[0], CUFFT_Z2D); if ( result != CUFFT_SUCCESS ) { amrex::Print() << " cufftplan3d backward failed! \n"; } #else - // Add 2D cuFFT-specific code for Z2D - // Note that Z2D is inherently backward plan + // Add 2D cuFFT-spacific code for D2Z + // Note that D2Z is inherently forward plan + result = cufftPlan2d( &forward_plan[mfi], fft_size[1], + fft_size[0], CUFFT_D2Z ); + if ( result != CUFFT_SUCCESS ) { + amrex::Print() << " cufftplan2d forward failed! \n"; + } + result = cufftPlan2d( &backward_plan[mfi], fft_size[1], + fft_size[0], CUFFT_Z2D ); + if ( result != CUFFT_SUCCESS ) { + amrex::Print() << " cufftplan2d backward failed! \n"; + } #endif #else diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 998fc4de1..00dcb85d0 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -222,9 +222,7 @@ WarpX::FillBoundaryE () { for (int lev = 0; lev <= finest_level; ++lev) { - amrex::Print() << " lev " << lev << " in fill bndry E \n"; FillBoundaryE(lev); - amrex::Print() << " lev " << lev << " after fill bndry E \n"; } } @@ -273,12 +271,9 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type) pml[lev]->FillBoundaryE(patch_type); } - amrex::Print() << " before defining multifab \n"; const auto& cperiod = Geom(lev-1).periodicity(); Vector mf{Efield_cp[lev][0].get(),Efield_cp[lev][1].get(),Efield_cp[lev][2].get()}; - amrex::Print() << " before amrex fill bndry \n"; amrex::FillBoundary(mf, cperiod); - amrex::Print() << " after amrex fill bndry \n"; } } diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 0826605ec..1f517fccb 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1388,9 +1388,7 @@ PhysicalParticleContainer::Evolve (int lev, pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); BL_PROFILE_VAR_STOP(blp_copy); - amrex::Print() << " before deposit chage \n"; if (rho) DepositCharge(pti, wp, rho, crho, 0, np_current, np, thread_num, lev); - amrex::Print() << " after deposit chage \n"; if (! do_not_push) { diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 576c53319..47d57294d 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -612,7 +612,6 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, const std::array& xyzmin = xyzmin_tile; #ifdef AMREX_USE_GPU - amrex::Print() << " before icomp data ptr " << icomp << "\n"; data_ptr = (*rhomf)[pti].dataPtr(icomp); auto rholen = (*rhomf)[pti].length(); #else @@ -635,7 +634,6 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, const long nz = rholen[1]-1-2*ngRho; #endif BL_PROFILE_VAR_START(blp_pxr_chd); - amrex::Print() << " before warpxcharge deposition \n"; warpx_charge_deposition(data_ptr, &np_current, m_xp[thread_num].dataPtr(), m_yp[thread_num].dataPtr(), -- cgit v1.2.3