diff options
Diffstat (limited to 'Source/FieldSolver')
| -rw-r--r-- | Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp | 52 | ||||
| -rw-r--r-- | Source/FieldSolver/WarpXFFT.cpp | 70 |
2 files changed, 39 insertions, 83 deletions
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<cuDoubleComplex*>( @@ -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<cuDoubleComplex*>( 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<const Real> 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<const Real> 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"; } |