diff options
Diffstat (limited to 'Source')
| -rw-r--r-- | Source/FieldSolver/Make.package | 5 | ||||
| -rw-r--r-- | Source/FieldSolver/PicsarHybridSpectralSolver/Make.package | 6 | ||||
| -rw-r--r-- | Source/FieldSolver/PicsarHybridSpectralSolver/PicsarHybridFFTData.H | 32 | ||||
| -rw-r--r-- | Source/FieldSolver/PicsarHybridSpectralSolver/PicsarHybridSpectralSolver.cpp (renamed from Source/FieldSolver/WarpXFFT.cpp) | 54 | ||||
| -rw-r--r-- | Source/FieldSolver/PicsarHybridSpectralSolver/picsar_hybrid_spectral.F90 (renamed from Source/FieldSolver/WarpX_fft.F90) | 0 | ||||
| -rw-r--r-- | Source/FieldSolver/WarpXPushFieldsEM.cpp | 48 | ||||
| -rw-r--r-- | Source/Initialization/WarpXInitData.cpp | 2 | ||||
| -rw-r--r-- | Source/Make.WarpX | 25 | ||||
| -rw-r--r-- | Source/WarpX.H | 64 | ||||
| -rw-r--r-- | Source/WarpX.cpp | 15 |
10 files changed, 137 insertions, 114 deletions
diff --git a/Source/FieldSolver/Make.package b/Source/FieldSolver/Make.package index 7eea52ee7..25e21c146 100644 --- a/Source/FieldSolver/Make.package +++ b/Source/FieldSolver/Make.package @@ -2,9 +2,10 @@ CEXE_headers += WarpX_K.H CEXE_headers += WarpX_FDTD.H CEXE_sources += WarpXPushFieldsEM.cpp ifeq ($(USE_PSATD),TRUE) - CEXE_sources += WarpXFFT.cpp - F90EXE_sources += WarpX_fft.F90 include $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/Make.package + ifeq ($(USE_CUDA),FALSE) + include $(WARPX_HOME)/Source/FieldSolver/PicsarHybridSpectralSolver/Make.package + endif endif ifeq ($(USE_OPENBC_POISSON),TRUE) F90EXE_sources += openbc_poisson_solver.F90 diff --git a/Source/FieldSolver/PicsarHybridSpectralSolver/Make.package b/Source/FieldSolver/PicsarHybridSpectralSolver/Make.package new file mode 100644 index 000000000..e185525e6 --- /dev/null +++ b/Source/FieldSolver/PicsarHybridSpectralSolver/Make.package @@ -0,0 +1,6 @@ +CEXE_headers += PicsarHybridFFTData.H +CEXE_sources += PicsarHybridSpectralSolver.cpp +F90EXE_sources += picsar_hybrid_spectral.F90 + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/PicsarHybridSpectralSolver +VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/PicsarHybridSpectralSolver
\ No newline at end of file diff --git a/Source/FieldSolver/PicsarHybridSpectralSolver/PicsarHybridFFTData.H b/Source/FieldSolver/PicsarHybridSpectralSolver/PicsarHybridFFTData.H new file mode 100644 index 000000000..4e97ab675 --- /dev/null +++ b/Source/FieldSolver/PicsarHybridSpectralSolver/PicsarHybridFFTData.H @@ -0,0 +1,32 @@ +#ifndef WARPX_PICSAR_HYBRID_FFTDATA_H_ +#define WARPX_PICSAR_HYBRID_FFTDATA_H_ + +// FFTData is a stuct containing a 22 pointers to auxiliary arrays +// 1-11: padded arrays in real space (required by FFTW); 12-22: arrays in spectral space +struct FFTData +{ + static constexpr int N = 22; + void* data[N] = { nullptr }; + + ~FFTData () { + for (int i = 0; i < N; ++i) { // The memory is allocated with fftw_alloc. + fftw_free(data[i]); + data[i] = nullptr; + } + } + + FFTData () = default; + + FFTData (FFTData && rhs) noexcept { + for (int i = 0; i < N; ++i) { + data[i] = rhs.data[i]; + rhs.data[i] = nullptr; + } + } + + FFTData (FFTData const&) = delete; + void operator= (FFTData const&) = delete; + void operator= (FFTData&&) = delete; +}; + +#endif // WARPX_PICSAR_HYBRID_FFTDATA_H_
\ No newline at end of file diff --git a/Source/FieldSolver/WarpXFFT.cpp b/Source/FieldSolver/PicsarHybridSpectralSolver/PicsarHybridSpectralSolver.cpp index f9cd7570a..26c93086a 100644 --- a/Source/FieldSolver/WarpXFFT.cpp +++ b/Source/FieldSolver/PicsarHybridSpectralSolver/PicsarHybridSpectralSolver.cpp @@ -5,10 +5,10 @@ using namespace amrex; -constexpr int WarpX::FFTData::N; +constexpr int FFTData::N; namespace { -static std::unique_ptr<WarpX::FFTData> nullfftdata; // This for process with nz_fft=0 +static std::unique_ptr<FFTData> nullfftdata; // This for process with nz_fft=0 /** \brief Returns an "owner mask" which 1 for all cells, except * for the duplicated (physical) cells of a nodal grid. @@ -122,6 +122,8 @@ WarpX::AllocLevelDataFFT (int lev) static_assert(std::is_standard_layout<FFTData>::value, "FFTData must have standard layout"); static_assert(sizeof(FFTData) == sizeof(void*)*FFTData::N, "sizeof FFTData is wrong"); + + InitFFTComm(lev); BoxArray ba_fp_fft; @@ -366,52 +368,8 @@ WarpX::FreeFFT (int lev) } void -WarpX::PushPSATD (amrex::Real a_dt) -{ - for (int lev = 0; lev <= finest_level; ++lev) { - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(dt[lev] == a_dt, "dt must be consistent"); - if (fft_hybrid_mpi_decomposition){ - PushPSATD_hybridFFT(lev, a_dt); - } else { - PushPSATD_localFFT(lev, a_dt); - } - } -} - -void WarpX::PushPSATD_localFFT (int lev, amrex::Real /* dt */) -{ - auto& solver = *spectral_solver_fp[lev]; - - // Perform forward Fourier transform - solver.ForwardTransform(*Efield_fp[lev][0], SpectralFieldIndex::Ex); - solver.ForwardTransform(*Efield_fp[lev][1], SpectralFieldIndex::Ey); - solver.ForwardTransform(*Efield_fp[lev][2], SpectralFieldIndex::Ez); - solver.ForwardTransform(*Bfield_fp[lev][0], SpectralFieldIndex::Bx); - solver.ForwardTransform(*Bfield_fp[lev][1], SpectralFieldIndex::By); - solver.ForwardTransform(*Bfield_fp[lev][2], SpectralFieldIndex::Bz); - solver.ForwardTransform(*current_fp[lev][0], SpectralFieldIndex::Jx); - solver.ForwardTransform(*current_fp[lev][1], SpectralFieldIndex::Jy); - solver.ForwardTransform(*current_fp[lev][2], SpectralFieldIndex::Jz); - solver.ForwardTransform(*rho_fp[lev], SpectralFieldIndex::rho_old, 0); - solver.ForwardTransform(*rho_fp[lev], SpectralFieldIndex::rho_new, 1); - - // Advance fields in spectral space - solver.pushSpectralFields(); - - // Perform backward Fourier Transform - solver.BackwardTransform(*Efield_fp[lev][0], SpectralFieldIndex::Ex); - solver.BackwardTransform(*Efield_fp[lev][1], SpectralFieldIndex::Ey); - solver.BackwardTransform(*Efield_fp[lev][2], SpectralFieldIndex::Ez); - solver.BackwardTransform(*Bfield_fp[lev][0], SpectralFieldIndex::Bx); - solver.BackwardTransform(*Bfield_fp[lev][1], SpectralFieldIndex::By); - solver.BackwardTransform(*Bfield_fp[lev][2], SpectralFieldIndex::Bz); -} - -void WarpX::PushPSATD_hybridFFT (int lev, amrex::Real /* dt */) { -#ifndef AMREX_USE_CUDA // Running on CPU ; use PICSAR code for the hybrid FFT - BL_PROFILE_VAR_NS("WarpXFFT::CopyDualGrid", blp_copy); BL_PROFILE_VAR_NS("PICSAR::FftPushEB", blp_push_eb); @@ -486,8 +444,4 @@ WarpX::PushPSATD_hybridFFT (int lev, amrex::Real /* dt */) { amrex::Abort("WarpX::PushPSATD: TODO"); } -#else // AMREX_USE_CUDA is defined ; running on GPU - amrex::Abort("The option `psatd.fft_hybrid_mpi_decomposition` does not work on GPU."); -#endif - } diff --git a/Source/FieldSolver/WarpX_fft.F90 b/Source/FieldSolver/PicsarHybridSpectralSolver/picsar_hybrid_spectral.F90 index 35357a63f..35357a63f 100644 --- a/Source/FieldSolver/WarpX_fft.F90 +++ b/Source/FieldSolver/PicsarHybridSpectralSolver/picsar_hybrid_spectral.F90 diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 0887fa6a2..4fce4717b 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -17,6 +17,54 @@ using namespace amrex; +#ifdef WARPX_USE_PSATD + +void +WarpX::PushPSATD (amrex::Real a_dt) +{ + for (int lev = 0; lev <= finest_level; ++lev) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(dt[lev] == a_dt, "dt must be consistent"); + if (fft_hybrid_mpi_decomposition){ +#ifndef AMREX_USE_CUDA // Only available on CPU + PushPSATD_hybridFFT(lev, a_dt); +#endif + } else { + PushPSATD_localFFT(lev, a_dt); + } + } +} + +void WarpX::PushPSATD_localFFT (int lev, amrex::Real /* dt */) +{ + auto& solver = *spectral_solver_fp[lev]; + + // Perform forward Fourier transform + solver.ForwardTransform(*Efield_fp[lev][0], SpectralFieldIndex::Ex); + solver.ForwardTransform(*Efield_fp[lev][1], SpectralFieldIndex::Ey); + solver.ForwardTransform(*Efield_fp[lev][2], SpectralFieldIndex::Ez); + solver.ForwardTransform(*Bfield_fp[lev][0], SpectralFieldIndex::Bx); + solver.ForwardTransform(*Bfield_fp[lev][1], SpectralFieldIndex::By); + solver.ForwardTransform(*Bfield_fp[lev][2], SpectralFieldIndex::Bz); + solver.ForwardTransform(*current_fp[lev][0], SpectralFieldIndex::Jx); + solver.ForwardTransform(*current_fp[lev][1], SpectralFieldIndex::Jy); + solver.ForwardTransform(*current_fp[lev][2], SpectralFieldIndex::Jz); + solver.ForwardTransform(*rho_fp[lev], SpectralFieldIndex::rho_old, 0); + solver.ForwardTransform(*rho_fp[lev], SpectralFieldIndex::rho_new, 1); + + // Advance fields in spectral space + solver.pushSpectralFields(); + + // Perform backward Fourier Transform + solver.BackwardTransform(*Efield_fp[lev][0], SpectralFieldIndex::Ex); + solver.BackwardTransform(*Efield_fp[lev][1], SpectralFieldIndex::Ey); + solver.BackwardTransform(*Efield_fp[lev][2], SpectralFieldIndex::Ez); + solver.BackwardTransform(*Bfield_fp[lev][0], SpectralFieldIndex::Bx); + solver.BackwardTransform(*Bfield_fp[lev][1], SpectralFieldIndex::By); + solver.BackwardTransform(*Bfield_fp[lev][2], SpectralFieldIndex::Bz); +} + +#endif + void WarpX::EvolveB (Real a_dt) { diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index d583b2b0f..2442e0205 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -317,7 +317,7 @@ WarpX::InitLevelData (int lev, Real time) } } -#ifdef WARPX_USE_PSATD +#ifdef WARPX_USE_PSATD_HYBRID void WarpX::InitLevelDataFFT (int lev, Real time) diff --git a/Source/Make.WarpX b/Source/Make.WarpX index bc837c797..3060ae8f0 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -111,20 +111,23 @@ endif ifeq ($(USE_PSATD),TRUE) - FFTW_HOME ?= NOT_SET - ifneq ($(FFTW_HOME),NOT_SET) - VPATH_LOCATIONS += $(FFTW_HOME)/include - INCLUDE_LOCATIONS += $(FFTW_HOME)/include - LIBRARY_LOCATIONS += $(FFTW_HOME)/lib - endif USERSuffix := $(USERSuffix).PSATD DEFINES += -DWARPX_USE_PSATD - ifeq ($(USE_CUDA),TRUE) - DEFINES += -DFFTW -DCUDA_FFT=1 # PICSAR uses it - LIBRARIES += -lfftw3_mpi -lfftw3 -lfftw3_threads -lcufft + ifeq ($(USE_CUDA),FALSE) # Running on CPU + # Use FFTW + LIBRARIES += -lfftw3_mpi -lfftw3 -lfftw3_threads + FFTW_HOME ?= NOT_SET + ifneq ($(FFTW_HOME),NOT_SET) + VPATH_LOCATIONS += $(FFTW_HOME)/include + INCLUDE_LOCATIONS += $(FFTW_HOME)/include + LIBRARY_LOCATIONS += $(FFTW_HOME)/lib + endif + # Compile with PICSAR's Hybrid-MPI PSATD + DEFINES += -DWARPX_USE_PSATD_HYBRID + DEFINES += -DFFTW # PICSAR uses it else - DEFINES += -DFFTW # PICSAR uses it - LIBRARIES += -lfftw3_mpi -lfftw3 -lfftw3_threads + # Use cuFFT + LIBRARIES += -lcufft endif endif diff --git a/Source/WarpX.H b/Source/WarpX.H index ca7596589..35de5c88d 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -26,9 +26,11 @@ #include <NCIGodfreyFilter.H> #ifdef WARPX_USE_PSATD -#include <fftw3.h> #include <SpectralSolver.H> #endif +#ifdef WARPX_USE_PSATD_HYBRID +#include <PicsarHybridFFTData.H> +#endif #if defined(BL_USE_SENSEI_INSITU) namespace amrex { @@ -578,7 +580,7 @@ private: amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_slice; amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_slice; -#ifdef WARPX_USE_PSATD +#ifdef WARPX_USE_PSATD_HYBRID // 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) amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_fp_fft; @@ -590,44 +592,30 @@ private: amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_cp_fft; amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_cp_fft; amrex::Vector< std::unique_ptr<amrex::MultiFab> > rho_cp_fft; +#endif -public: - // FFTData is a stuct containing a 22 pointers to auxiliary arrays - // 1-11: padded arrays in real space (required by FFTW); 12-22: arrays in spectral space - struct FFTData - { - static constexpr int N = 22; - void* data[N] = { nullptr }; - - ~FFTData () { - for (int i = 0; i < N; ++i) { // The memory is allocated with fftw_alloc. - fftw_free(data[i]); - data[i] = nullptr; - } - } - - FFTData () = default; +#ifdef WARPX_USE_PSATD +private: + void EvolvePSATD (int numsteps); + void PushPSATD (amrex::Real dt); + void PushPSATD_localFFT (int lev, amrex::Real dt); - FFTData (FFTData && rhs) noexcept { - for (int i = 0; i < N; ++i) { - data[i] = rhs.data[i]; - rhs.data[i] = nullptr; - } - } + bool fft_hybrid_mpi_decomposition = false; + int ngroups_fft = 4; + int fftw_plan_measure = 1; + int nox_fft = 16; + int noy_fft = 16; + int noz_fft = 16; - FFTData (FFTData const&) = delete; - void operator= (FFTData const&) = delete; - void operator= (FFTData&&) = delete; - }; + amrex::Vector<std::unique_ptr<SpectralSolver>> spectral_solver_fp; + amrex::Vector<std::unique_ptr<SpectralSolver>> spectral_solver_cp; +#endif +#ifdef WARPX_USE_PSATD_HYBRID private: - amrex::Vector<std::unique_ptr<amrex::LayoutData<FFTData> > > dataptr_fp_fft; amrex::Vector<std::unique_ptr<amrex::LayoutData<FFTData> > > dataptr_cp_fft; - amrex::Vector<std::unique_ptr<SpectralSolver>> spectral_solver_fp; - amrex::Vector<std::unique_ptr<SpectralSolver>> spectral_solver_cp; - // Domain decomposition containing the valid part of the dual grids (i.e. without FFT guard cells) amrex::Vector<amrex::BoxArray> ba_valid_fp_fft; amrex::Vector<amrex::BoxArray> ba_valid_cp_fft; @@ -638,13 +626,6 @@ private: amrex::Vector<MPI_Comm> comm_fft; amrex::Vector<int> color_fft; - bool fft_hybrid_mpi_decomposition = false; - int ngroups_fft = 4; - int fftw_plan_measure = 1; - int nox_fft = 16; - int noy_fft = 16; - int noz_fft = 16; - void AllocLevelDataFFT (int lev); void InitLevelDataFFT (int lev, amrex::Real time); @@ -654,12 +635,7 @@ private: const amrex::Box& domain); void InitFFTDataPlan (int lev); void FreeFFT (int lev); - - void EvolvePSATD (int numsteps); - void PushPSATD (amrex::Real dt); - void PushPSATD_localFFT (int lev, amrex::Real dt); void PushPSATD_hybridFFT (int lev, amrex::Real dt); - #endif #if defined(BL_USE_SENSEI_INSITU) diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 08948227c..f379d5f4c 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -202,6 +202,10 @@ WarpX::WarpX () costs.resize(nlevs_max); #ifdef WARPX_USE_PSATD + spectral_solver_fp.resize(nlevs_max); + spectral_solver_cp.resize(nlevs_max); +#endif +#ifdef WARPX_USE_PSATD_HYBRID Efield_fp_fft.resize(nlevs_max); Bfield_fp_fft.resize(nlevs_max); current_fp_fft.resize(nlevs_max); @@ -215,9 +219,6 @@ WarpX::WarpX () dataptr_fp_fft.resize(nlevs_max); dataptr_cp_fft.resize(nlevs_max); - spectral_solver_fp.resize(nlevs_max); - spectral_solver_cp.resize(nlevs_max); - ba_valid_fp_fft.resize(nlevs_max); ba_valid_cp_fft.resize(nlevs_max); @@ -511,8 +512,6 @@ WarpX::ReadParameters () pp.query("nox", nox_fft); pp.query("noy", noy_fft); pp.query("noz", noz_fft); - // Override value - if (fft_hybrid_mpi_decomposition==false) ngroups_fft=ParallelDescriptor::NProcs(); } #endif @@ -568,8 +567,12 @@ WarpX::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& new_grids, #ifdef WARPX_USE_PSATD if (fft_hybrid_mpi_decomposition){ +#ifdef WARPX_USE_PSATD_HYBRID AllocLevelDataFFT(lev); InitLevelDataFFT(lev, time); +#else + amrex::Abort("The option `psatd.fft_hybrid_mpi_decomposition` does not work on GPU."); +#endif } #endif } @@ -614,7 +617,7 @@ WarpX::ClearLevel (int lev) costs[lev].reset(); -#ifdef WARPX_USE_PSATD +#ifdef WARPX_USE_PSATD_HYBRID for (int i = 0; i < 3; ++i) { Efield_fp_fft[lev][i].reset(); Bfield_fp_fft[lev][i].reset(); |