diff options
Diffstat (limited to 'Source/FieldSolver')
13 files changed, 425 insertions, 256 deletions
diff --git a/Source/FieldSolver/SpectralSolver/Make.package b/Source/FieldSolver/SpectralSolver/Make.package index 50127914d..b0ee54095 100644 --- a/Source/FieldSolver/SpectralSolver/Make.package +++ b/Source/FieldSolver/SpectralSolver/Make.package @@ -1,11 +1,12 @@ CEXE_headers += WarpX_Complex.H CEXE_headers += SpectralSolver.H +CEXE_sources += SpectralSolver.cpp CEXE_headers += SpectralFieldData.H CEXE_sources += SpectralFieldData.cpp -CEXE_headers += PsatdAlgorithm.H -CEXE_sources += PsatdAlgorithm.cpp CEXE_headers += SpectralKSpace.H CEXE_sources += SpectralKSpace.cpp +include $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package + INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package new file mode 100644 index 000000000..c62c21f44 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package @@ -0,0 +1,6 @@ +CEXE_headers += SpectralBaseAlgorithm.H +CEXE_headers += PsatdAlgorithm.H +CEXE_sources += PsatdAlgorithm.cpp + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms +VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H index acefcc466..12718e38b 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H @@ -1,31 +1,27 @@ #ifndef WARPX_PSATD_ALGORITHM_H_ #define WARPX_PSATD_ALGORITHM_H_ -#include <SpectralKSpace.H> -#include <SpectralFieldData.H> +#include <SpectralBaseAlgorithm.H> /* \brief Class that updates the field in spectral space * and stores the coefficients of the corresponding update equation. */ -class PsatdAlgorithm +class PsatdAlgorithm : public SpectralBaseAlgorithm { - using SpectralCoefficients = amrex::FabArray< amrex::BaseFab <amrex::Real> >; public: PsatdAlgorithm(const SpectralKSpace& spectral_kspace, const amrex::DistributionMapping& dm, const int norder_x, const int norder_y, const int norder_z, const bool nodal, const amrex::Real dt); - PsatdAlgorithm() = default; // Default constructor - PsatdAlgorithm& operator=(PsatdAlgorithm&& algorithm) = default; - void pushSpectralFields(SpectralFieldData& f) const; + + void InitializeSpectralCoefficients(const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const amrex::Real dt); + + void pushSpectralFields(SpectralFieldData& f) const override final; private: - // Modified finite-order vectors - KVectorComponent modified_kx_vec, modified_kz_vec; -#if (AMREX_SPACEDIM==3) - KVectorComponent modified_ky_vec; -#endif SpectralCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef; }; diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp index 56e58bcc4..d45b01bda 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp @@ -9,14 +9,9 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, const DistributionMapping& dm, const int norder_x, const int norder_y, const int norder_z, const bool nodal, const Real dt) -// Compute and assign the modified k vectors -: modified_kx_vec(spectral_kspace.getModifiedKComponent(dm,0,norder_x,nodal)), -#if (AMREX_SPACEDIM==3) - modified_ky_vec(spectral_kspace.getModifiedKComponent(dm,1,norder_y,nodal)), - modified_kz_vec(spectral_kspace.getModifiedKComponent(dm,2,norder_z,nodal)) -#else - modified_kz_vec(spectral_kspace.getModifiedKComponent(dm,1,norder_z,nodal)) -#endif + // Initialize members of base class + : SpectralBaseAlgorithm( spectral_kspace, dm, + norder_x, norder_y, norder_z, nodal ) { const BoxArray& ba = spectral_kspace.spectralspace_ba; @@ -27,57 +22,8 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, X2_coef = SpectralCoefficients(ba, dm, 1, 0); X3_coef = SpectralCoefficients(ba, dm, 1, 0); - // Fill them with the right values: - // Loop over boxes and allocate the corresponding coefficients - // for each box owned by the local MPI proc - for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi){ - - const Box& bx = ba[mfi]; - - // Extract pointers for the k vectors - const Real* modified_kx = modified_kx_vec[mfi].dataPtr(); -#if (AMREX_SPACEDIM==3) - const Real* modified_ky = modified_ky_vec[mfi].dataPtr(); -#endif - const Real* modified_kz = modified_kz_vec[mfi].dataPtr(); - // Extract arrays for the coefficients - Array4<Real> C = C_coef[mfi].array(); - Array4<Real> S_ck = S_ck_coef[mfi].array(); - Array4<Real> X1 = X1_coef[mfi].array(); - Array4<Real> X2 = X2_coef[mfi].array(); - Array4<Real> X3 = X3_coef[mfi].array(); - - // Loop over indices within one box - ParallelFor(bx, - [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - // Calculate norm of vector - const Real k_norm = std::sqrt( - std::pow(modified_kx[i], 2) + -#if (AMREX_SPACEDIM==3) - std::pow(modified_ky[j], 2) + -#endif - std::pow(modified_kz[k], 2)); - - // Calculate coefficients - constexpr Real c = PhysConst::c; - constexpr Real ep0 = PhysConst::ep0; - if (k_norm != 0){ - C(i,j,k) = std::cos(c*k_norm*dt); - S_ck(i,j,k) = std::sin(c*k_norm*dt)/(c*k_norm); - X1(i,j,k) = (1. - C(i,j,k))/(ep0 * c*c * k_norm*k_norm); - X2(i,j,k) = (1. - S_ck(i,j,k)/dt)/(ep0 * k_norm*k_norm); - X3(i,j,k) = (C(i,j,k) - S_ck(i,j,k)/dt)/(ep0 * k_norm*k_norm); - } else { // Handle k_norm = 0, by using the analytical limit - C(i,j,k) = 1.; - S_ck(i,j,k) = dt; - X1(i,j,k) = 0.5 * dt*dt / ep0; - X2(i,j,k) = c*c * dt*dt / (6.*ep0); - X3(i,j,k) = - c*c * dt*dt / (3.*ep0); - } - }); - } -}; + InitializeSpectralCoefficients(spectral_kspace, dm, dt); +} /* Advance the E and B field in spectral space (stored in `f`) * over one time step */ @@ -85,23 +31,12 @@ void PsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ // Loop over boxes - for (MFIter mfi(f.Ex); mfi.isValid(); ++mfi){ + for (MFIter mfi(f.fields); mfi.isValid(); ++mfi){ - const Box& bx = f.Ex[mfi].box(); + const Box& bx = f.fields[mfi].box(); // Extract arrays for the fields to be updated - Array4<Complex> Ex_arr = f.Ex[mfi].array(); - Array4<Complex> Ey_arr = f.Ey[mfi].array(); - Array4<Complex> Ez_arr = f.Ez[mfi].array(); - Array4<Complex> Bx_arr = f.Bx[mfi].array(); - Array4<Complex> By_arr = f.By[mfi].array(); - Array4<Complex> Bz_arr = f.Bz[mfi].array(); - // Extract arrays for J and rho - Array4<const Complex> Jx_arr = f.Jx[mfi].array(); - Array4<const Complex> Jy_arr = f.Jy[mfi].array(); - Array4<const Complex> Jz_arr = f.Jz[mfi].array(); - Array4<const Complex> rho_old_arr = f.rho_old[mfi].array(); - Array4<const Complex> rho_new_arr = f.rho_new[mfi].array(); + Array4<Complex> fields = f.fields[mfi].array(); // Extract arrays for the coefficients Array4<const Real> C_arr = C_coef[mfi].array(); Array4<const Real> S_ck_arr = S_ck_coef[mfi].array(); @@ -120,55 +55,118 @@ PsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Record old values of the fields to be updated - const Complex Ex_old = Ex_arr(i,j,k); - const Complex Ey_old = Ey_arr(i,j,k); - const Complex Ez_old = Ez_arr(i,j,k); - const Complex Bx_old = Bx_arr(i,j,k); - const Complex By_old = By_arr(i,j,k); - const Complex Bz_old = Bz_arr(i,j,k); + using Idx = SpectralFieldIndex; + const Complex Ex_old = fields(i,j,k,Idx::Ex); + const Complex Ey_old = fields(i,j,k,Idx::Ey); + const Complex Ez_old = fields(i,j,k,Idx::Ez); + const Complex Bx_old = fields(i,j,k,Idx::Bx); + const Complex By_old = fields(i,j,k,Idx::By); + const Complex Bz_old = fields(i,j,k,Idx::Bz); // Shortcut for the values of J and rho - const Complex Jx = Jx_arr(i,j,k); - const Complex Jy = Jy_arr(i,j,k); - const Complex Jz = Jz_arr(i,j,k); - const Complex rho_old = rho_old_arr(i,j,k); - const Complex rho_new = rho_new_arr(i,j,k); + const Complex Jx = fields(i,j,k,Idx::Jx); + const Complex Jy = fields(i,j,k,Idx::Jy); + const Complex Jz = fields(i,j,k,Idx::Jz); + const Complex rho_old = fields(i,j,k,Idx::rho_old); + const Complex rho_new = fields(i,j,k,Idx::rho_new); // k vector values, and coefficients const Real kx = modified_kx_arr[i]; #if (AMREX_SPACEDIM==3) const Real ky = modified_ky_arr[j]; + const Real kz = modified_kz_arr[k]; #else constexpr Real ky = 0; + const Real kz = modified_kz_arr[j]; #endif - const Real kz = modified_kz_arr[k]; constexpr Real c2 = PhysConst::c*PhysConst::c; constexpr Real inv_ep0 = 1./PhysConst::ep0; - constexpr Complex I = Complex{0,1}; + const Complex I = Complex{0,1}; const Real C = C_arr(i,j,k); const Real S_ck = S_ck_arr(i,j,k); const Real X1 = X1_arr(i,j,k); const Real X2 = X2_arr(i,j,k); const Real X3 = X3_arr(i,j,k); + // Update E (see WarpX online documentation: theory section) - Ex_arr(i,j,k) = C*Ex_old + fields(i,j,k,Idx::Ex) = C*Ex_old + S_ck*(c2*I*(ky*Bz_old - kz*By_old) - inv_ep0*Jx) - I*(X2*rho_new - X3*rho_old)*kx; - Ey_arr(i,j,k) = C*Ey_old + fields(i,j,k,Idx::Ey) = C*Ey_old + S_ck*(c2*I*(kz*Bx_old - kx*Bz_old) - inv_ep0*Jy) - I*(X2*rho_new - X3*rho_old)*ky; - Ez_arr(i,j,k) = C*Ez_old + fields(i,j,k,Idx::Ez) = C*Ez_old + S_ck*(c2*I*(kx*By_old - ky*Bx_old) - inv_ep0*Jz) - I*(X2*rho_new - X3*rho_old)*kz; // Update B (see WarpX online documentation: theory section) - Bx_arr(i,j,k) = C*Bx_old + fields(i,j,k,Idx::Bx) = C*Bx_old - S_ck*I*(ky*Ez_old - kz*Ey_old) + X1*I*(ky*Jz - kz*Jy); - By_arr(i,j,k) = C*By_old + fields(i,j,k,Idx::By) = C*By_old - S_ck*I*(kz*Ex_old - kx*Ez_old) + X1*I*(kz*Jx - kx*Jz); - Bz_arr(i,j,k) = C*Bz_old + fields(i,j,k,Idx::Bz) = C*Bz_old - S_ck*I*(kx*Ey_old - ky*Ex_old) + X1*I*(kx*Jy - ky*Jx); }); } }; + +void PsatdAlgorithm::InitializeSpectralCoefficients(const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const amrex::Real dt) +{ + const BoxArray& ba = spectral_kspace.spectralspace_ba; + // Fill them with the right values: + // Loop over boxes and allocate the corresponding coefficients + // for each box owned by the local MPI proc + for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi){ + + const Box& bx = ba[mfi]; + + // Extract pointers for the k vectors + const Real* modified_kx = modified_kx_vec[mfi].dataPtr(); +#if (AMREX_SPACEDIM==3) + const Real* modified_ky = modified_ky_vec[mfi].dataPtr(); +#endif + const Real* modified_kz = modified_kz_vec[mfi].dataPtr(); + // Extract arrays for the coefficients + Array4<Real> C = C_coef[mfi].array(); + Array4<Real> S_ck = S_ck_coef[mfi].array(); + Array4<Real> X1 = X1_coef[mfi].array(); + Array4<Real> X2 = X2_coef[mfi].array(); + Array4<Real> X3 = X3_coef[mfi].array(); + + // Loop over indices within one box + ParallelFor(bx, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Calculate norm of vector + const Real k_norm = std::sqrt( + std::pow(modified_kx[i], 2) + +#if (AMREX_SPACEDIM==3) + std::pow(modified_ky[j], 2) + + std::pow(modified_kz[k], 2)); +#else + std::pow(modified_kz[j], 2)); +#endif + + + // Calculate coefficients + constexpr Real c = PhysConst::c; + constexpr Real ep0 = PhysConst::ep0; + if (k_norm != 0){ + C(i,j,k) = std::cos(c*k_norm*dt); + S_ck(i,j,k) = std::sin(c*k_norm*dt)/(c*k_norm); + X1(i,j,k) = (1. - C(i,j,k))/(ep0 * c*c * k_norm*k_norm); + X2(i,j,k) = (1. - S_ck(i,j,k)/dt)/(ep0 * k_norm*k_norm); + X3(i,j,k) = (C(i,j,k) - S_ck(i,j,k)/dt)/(ep0 * k_norm*k_norm); + } else { // Handle k_norm = 0, by using the analytical limit + C(i,j,k) = 1.; + S_ck(i,j,k) = dt; + X1(i,j,k) = 0.5 * dt*dt / ep0; + X2(i,j,k) = c*c * dt*dt / (6.*ep0); + X3(i,j,k) = - c*c * dt*dt / (3.*ep0); + } + }); + } +} diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H new file mode 100644 index 000000000..602eb2473 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H @@ -0,0 +1,51 @@ +#ifndef WARPX_SPECTRAL_BASE_ALGORITHM_H_ +#define WARPX_SPECTRAL_BASE_ALGORITHM_H_ + +#include <SpectralKSpace.H> +#include <SpectralFieldData.H> + +/* \brief Class that updates the field in spectral space + * and stores the coefficients of the corresponding update equation. + * + * `SpectralBaseAlgorithm` is only a base class and cannot be used directly. + * Instead use its subclasses, which implement the specific field update + * equations for a given spectral algorithm. + */ +class SpectralBaseAlgorithm +{ + public: + // Member function that updates the fields in spectral space ; + // meant to be overridden in subclasses + virtual void pushSpectralFields(SpectralFieldData& f) const = 0; + // The destructor should also be a virtual function, so that + // a pointer to subclass of `SpectraBaseAlgorithm` actually + // calls the subclass's destructor. + virtual ~SpectralBaseAlgorithm() {}; + + protected: // Meant to be used in the subclasses + + using SpectralCoefficients = amrex::FabArray< amrex::BaseFab <amrex::Real> >; + + // Constructor + SpectralBaseAlgorithm(const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const int norder_x, const int norder_y, + const int norder_z, const bool nodal) + // Compute and assign the modified k vectors + : modified_kx_vec(spectral_kspace.getModifiedKComponent(dm,0,norder_x,nodal)), +#if (AMREX_SPACEDIM==3) + modified_ky_vec(spectral_kspace.getModifiedKComponent(dm,1,norder_y,nodal)), + modified_kz_vec(spectral_kspace.getModifiedKComponent(dm,2,norder_z,nodal)) +#else + modified_kz_vec(spectral_kspace.getModifiedKComponent(dm,1,norder_z,nodal)) +#endif + {}; + + // Modified finite-order vectors + KVectorComponent modified_kx_vec, modified_kz_vec; +#if (AMREX_SPACEDIM==3) + KVectorComponent modified_ky_vec; +#endif +}; + +#endif // WARPX_SPECTRAL_BASE_ALGORITHM_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index c62513de6..7954414b8 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -9,8 +9,9 @@ using SpectralField = amrex::FabArray< amrex::BaseFab <Complex> >; /* Index for the fields that will be stored in spectral space */ -struct SpectralFieldIndex{ - enum { Ex=0, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz, rho_old, rho_new }; +struct SpectralFieldIndex { + enum { Ex=0, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz, rho_old, rho_new, n_fields }; + // n_fields is automatically the total number of fields }; /* \brief Class that stores the fields in spectral space, and performs the @@ -24,7 +25,7 @@ class SpectralFieldData // (plans are only initialized for the boxes that are owned by // the local MPI rank) #ifdef AMREX_USE_GPU - // Add cuFFT-specific code + using FFTplans = amrex::LayoutData<cufftHandle>; #else using FFTplans = amrex::LayoutData<fftw_plan>; #endif @@ -37,15 +38,17 @@ class SpectralFieldData SpectralFieldData& operator=(SpectralFieldData&& field_data) = default; ~SpectralFieldData(); void ForwardTransform( const amrex::MultiFab& mf, - const int field_index, const int i_comp ); + const int field_index, const int i_comp); void BackwardTransform( amrex::MultiFab& mf, - const int field_index, const int i_comp ); + const int field_index, const int i_comp); private: - SpectralField Ex, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz, rho_old, rho_new; + // `fields` stores fields in spectral space, as multicomponent FabArray + SpectralField fields; // tmpRealField and tmpSpectralField store fields // right before/after the Fourier transform - SpectralField tmpRealField, tmpSpectralField; + SpectralField tmpSpectralField; // contains Complexs + amrex::MultiFab tmpRealField; // contains Reals FFTplans forward_plan, backward_plan; // Correcting "shift" factors when performing FFT from/to // a cell-centered grid in real space, instead of a nodal grid @@ -54,7 +57,6 @@ class SpectralFieldData #if (AMREX_SPACEDIM==3) SpectralShiftFactor yshift_FFTfromCell, yshift_FFTtoCell; #endif - SpectralField& getSpectralField( const int field_index ); }; #endif // WARPX_SPECTRAL_FIELD_DATA_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 291fe945e..a2b695568 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -10,21 +10,13 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, const BoxArray& spectralspace_ba = k_space.spectralspace_ba; // Allocate the arrays that contain the fields in spectral space - Ex = SpectralField(spectralspace_ba, dm, 1, 0); - Ey = SpectralField(spectralspace_ba, dm, 1, 0); - Ez = SpectralField(spectralspace_ba, dm, 1, 0); - Bx = SpectralField(spectralspace_ba, dm, 1, 0); - By = SpectralField(spectralspace_ba, dm, 1, 0); - Bz = SpectralField(spectralspace_ba, dm, 1, 0); - Jx = SpectralField(spectralspace_ba, dm, 1, 0); - Jy = SpectralField(spectralspace_ba, dm, 1, 0); - Jz = SpectralField(spectralspace_ba, dm, 1, 0); - rho_old = SpectralField(spectralspace_ba, dm, 1, 0); - rho_new = SpectralField(spectralspace_ba, dm, 1, 0); + // (one component per field) + fields = SpectralField(spectralspace_ba, dm, + SpectralFieldIndex::n_fields, 0); // Allocate temporary arrays - in real space and spectral space // These arrays will store the data just before/after the FFT - tmpRealField = SpectralField(realspace_ba, dm, 1, 0); + tmpRealField = MultiFab(realspace_ba, dm, 1, 0); tmpSpectralField = SpectralField(spectralspace_ba, dm, 1, 0); // By default, we assume the FFT is done from/to a nodal grid in real space @@ -56,31 +48,65 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, // Loop over boxes and allocate the corresponding plan // for each box owned by the local MPI proc for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ - Box bx = spectralspace_ba[mfi]; + // Note: the size of the real-space box and spectral-space box + // differ when using real-to-complex FFT. When initializing + // the FFT plan, the valid dimensions are those of the real-space box. + IntVect fft_size = realspace_ba[mfi].length(); #ifdef AMREX_USE_GPU - // Add cuFFT-specific code + // Create cuFFT plans + // Creating 3D plan for real to complex -- double precision + // Assuming CUDA is used for programming GPU + // Note that D2Z is inherently forward plan + // and Z2D is inherently backward plan + cufftResult result; +#if (AMREX_SPACEDIM == 3) + 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"; + } + + 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 + 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 // Create FFTW plans forward_plan[mfi] = // Swap dimensions: AMReX FAB are Fortran-order but FFTW is C-order #if (AMREX_SPACEDIM == 3) - fftw_plan_dft_3d( bx.length(2), bx.length(1), bx.length(0), + fftw_plan_dft_r2c_3d( fft_size[2], fft_size[1], fft_size[0], #else - fftw_plan_dft_2d( bx.length(1), bx.length(0), + fftw_plan_dft_r2c_2d( fft_size[1], fft_size[0], #endif - reinterpret_cast<fftw_complex*>( tmpRealField[mfi].dataPtr() ), + tmpRealField[mfi].dataPtr(), reinterpret_cast<fftw_complex*>( tmpSpectralField[mfi].dataPtr() ), - FFTW_FORWARD, FFTW_ESTIMATE ); + FFTW_ESTIMATE ); backward_plan[mfi] = // Swap dimensions: AMReX FAB are Fortran-order but FFTW is C-order #if (AMREX_SPACEDIM == 3) - fftw_plan_dft_3d( bx.length(2), bx.length(1), bx.length(0), + fftw_plan_dft_c2r_3d( fft_size[2], fft_size[1], fft_size[0], #else - fftw_plan_dft_2d( bx.length(1), bx.length(0), + fftw_plan_dft_c2r_2d( fft_size[1], fft_size[0], #endif reinterpret_cast<fftw_complex*>( tmpSpectralField[mfi].dataPtr() ), - reinterpret_cast<fftw_complex*>( tmpRealField[mfi].dataPtr() ), - FFTW_BACKWARD, FFTW_ESTIMATE ); + tmpRealField[mfi].dataPtr(), + FFTW_ESTIMATE ); #endif } } @@ -91,7 +117,9 @@ SpectralFieldData::~SpectralFieldData() if (tmpRealField.size() > 0){ for ( MFIter mfi(tmpRealField); mfi.isValid(); ++mfi ){ #ifdef AMREX_USE_GPU - // Add cuFFT-specific code + // Destroy cuFFT plans + cufftDestroy( forward_plan[mfi] ); + cufftDestroy( backward_plan[mfi] ); #else // Destroy FFTW plans fftw_destroy_plan( forward_plan[mfi] ); @@ -131,28 +159,38 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, realspace_bx.enclosedCells(); // Discard last point in nodal direction AMREX_ALWAYS_ASSERT( realspace_bx == tmpRealField[mfi].box() ); Array4<const Real> mf_arr = mf[mfi].array(); - Array4<Complex> tmp_arr = tmpRealField[mfi].array(); + Array4<Real> tmp_arr = tmpRealField[mfi].array(); ParallelFor( realspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - tmp_arr(i,j,k) = mf_arr(i,j,k,i_comp); + tmp_arr(i,j,k) = mf_arr(i,j,k,i_comp); }); } // Perform Fourier transform from `tmpRealField` to `tmpSpectralField` #ifdef AMREX_USE_GPU - // Add cuFFT-specific code ; make sure that this is done on the same - // GPU stream as the above copy + // Perform Fast Fourier Transform on GPU using cuFFT + // make sure that this is done on the same + // GPU stream as the above copy + cufftResult result; + cudaStream_t stream = amrex::Gpu::Device::cudaStream(); + cufftSetStream ( forward_plan[mfi], stream); + result = cufftExecD2Z( forward_plan[mfi], + tmpRealField[mfi].dataPtr(), + reinterpret_cast<cuDoubleComplex*>( + tmpSpectralField[mfi].dataPtr()) ); + if ( result != CUFFT_SUCCESS ) { + amrex::Print() << " forward transform using cufftExecD2Z failed ! \n"; + } #else fftw_execute( forward_plan[mfi] ); #endif // Copy the spectral-space field `tmpSpectralField` to the appropriate - // field (specified by the input argument field_index ) + // index of the FabArray `fields` (specified by `field_index`) // and apply correcting shift factor if the real space data comes // from a cell-centered grid in real space instead of a nodal grid. { - SpectralField& field = getSpectralField( field_index ); - Array4<Complex> field_arr = field[mfi].array(); + Array4<Complex> fields_arr = SpectralFieldData::fields[mfi].array(); Array4<const Complex> tmp_arr = tmpSpectralField[mfi].array(); const Complex* xshift_arr = xshift_FFTfromCell[mfi].dataPtr(); #if (AMREX_SPACEDIM == 3) @@ -161,6 +199,7 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, const Complex* zshift_arr = zshift_FFTfromCell[mfi].dataPtr(); // Loop over indices within one box const Box spectralspace_bx = tmpSpectralField[mfi].box(); + ParallelFor( spectralspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { Complex spectral_field_value = tmp_arr(i,j,k); @@ -168,10 +207,12 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, if (is_nodal_x==false) spectral_field_value *= xshift_arr[i]; #if (AMREX_SPACEDIM == 3) if (is_nodal_y==false) spectral_field_value *= yshift_arr[j]; -#endif if (is_nodal_z==false) spectral_field_value *= zshift_arr[k]; - // Copy field into temporary array - field_arr(i,j,k) = spectral_field_value; +#elif (AMREX_SPACEDIM == 2) + if (is_nodal_z==false) spectral_field_value *= zshift_arr[j]; +#endif + // Copy field into the right index + fields_arr(i,j,k,field_index) = spectral_field_value; }); } } @@ -182,7 +223,8 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, * real space, and store it in the component `i_comp` of `mf` */ void SpectralFieldData::BackwardTransform( MultiFab& mf, - const int field_index, const int i_comp ) + const int field_index, + const int i_comp ) { // Check field index type, in order to apply proper shift in spectral space const bool is_nodal_x = mf.is_nodal(0); @@ -200,10 +242,8 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, // field (specified by the input argument field_index) // and apply correcting shift factor if the field is to be transformed // to a cell-centered grid in real space instead of a nodal grid. - // Normalize (divide by 1/N) since the FFT+IFFT results in a factor N { - SpectralField& field = getSpectralField( field_index ); - Array4<const Complex> field_arr = field[mfi].array(); + Array4<const Complex> field_arr = SpectralFieldData::fields[mfi].array(); Array4<Complex> tmp_arr = tmpSpectralField[mfi].array(); const Complex* xshift_arr = xshift_FFTtoCell[mfi].dataPtr(); #if (AMREX_SPACEDIM == 3) @@ -212,60 +252,57 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, const Complex* zshift_arr = zshift_FFTtoCell[mfi].dataPtr(); // Loop over indices within one box const Box spectralspace_bx = tmpSpectralField[mfi].box(); - // For normalization: divide by the number of points in the box - const Real inv_N = 1./spectralspace_bx.numPts(); + ParallelFor( spectralspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - Complex spectral_field_value = field_arr(i,j,k); + Complex spectral_field_value = field_arr(i,j,k,field_index); // Apply proper shift in each dimension if (is_nodal_x==false) spectral_field_value *= xshift_arr[i]; #if (AMREX_SPACEDIM == 3) if (is_nodal_y==false) spectral_field_value *= yshift_arr[j]; -#endif if (is_nodal_z==false) spectral_field_value *= zshift_arr[k]; - // Copy field into temporary array (after normalization) - tmp_arr(i,j,k) = inv_N*spectral_field_value; +#elif (AMREX_SPACEDIM == 2) + if (is_nodal_z==false) spectral_field_value *= zshift_arr[j]; +#endif + // Copy field into temporary array + tmp_arr(i,j,k) = spectral_field_value; }); } // Perform Fourier transform from `tmpSpectralField` to `tmpRealField` #ifdef AMREX_USE_GPU - // Add cuFFT-specific code ; make sure that this is done on the same + // Perform Fast Fourier Transform on GPU using cuFFT. + // make sure that this is done on the same // GPU stream as the above copy + cufftResult result; + cudaStream_t stream = amrex::Gpu::Device::cudaStream(); + cufftSetStream ( backward_plan[mfi], stream); + result = cufftExecZ2D( backward_plan[mfi], + reinterpret_cast<cuDoubleComplex*>( + tmpSpectralField[mfi].dataPtr()), + tmpRealField[mfi].dataPtr() ); + if ( result != CUFFT_SUCCESS ) { + amrex::Print() << " Backward transform using cufftexecZ2D failed! \n"; + } #else fftw_execute( backward_plan[mfi] ); #endif // Copy the temporary field `tmpRealField` to the real-space field `mf` + + // Normalize (divide by 1/N) since the FFT+IFFT results in a factor N { const Box realspace_bx = tmpRealField[mfi].box(); Array4<Real> mf_arr = mf[mfi].array(); - Array4<const Complex> tmp_arr = tmpRealField[mfi].array(); + Array4<const Real> tmp_arr = tmpRealField[mfi].array(); + // Normalization: divide by the number of points in realspace + const Real inv_N = 1./realspace_bx.numPts(); + ParallelFor( realspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - mf_arr(i,j,k,i_comp) = tmp_arr(i,j,k).real(); + // Copy and normalize field + mf_arr(i,j,k,i_comp) = inv_N*tmp_arr(i,j,k); }); } } } - - -SpectralField& -SpectralFieldData::getSpectralField( const int field_index ) -{ - switch(field_index) - { - case SpectralFieldIndex::Ex : return Ex; break; - case SpectralFieldIndex::Ey : return Ey; break; - case SpectralFieldIndex::Ez : return Ez; break; - case SpectralFieldIndex::Bx : return Bx; break; - case SpectralFieldIndex::By : return By; break; - case SpectralFieldIndex::Bz : return Bz; break; - case SpectralFieldIndex::Jx : return Jx; break; - case SpectralFieldIndex::Jy : return Jy; break; - case SpectralFieldIndex::Jz : return Jz; break; - case SpectralFieldIndex::rho_old : return rho_old; break; - case SpectralFieldIndex::rho_new : return rho_new; break; - default : return tmpSpectralField; // For synthax; should not occur in practice - } -} diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H index ad17e6423..ae7124b2e 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H @@ -33,7 +33,9 @@ class SpectralKSpace const amrex::DistributionMapping& dm, const amrex::RealVect realspace_dx ); KVectorComponent getKComponent( - const amrex::DistributionMapping& dm, const int i_dim ) const; + const amrex::DistributionMapping& dm, + const amrex::BoxArray& realspace_ba, + const int i_dim, const bool only_positive_k ) const; KVectorComponent getModifiedKComponent( const amrex::DistributionMapping& dm, const int i_dim, const int n_order, const bool nodal ) const; diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index d91891a30..6fe5e3939 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -28,19 +28,33 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, // For local FFTs, boxes in spectral space start at 0 in // each direction and have the same number of points as the // (cell-centered) real space box - // TODO: this will be different for the real-to-complex FFT // TODO: this will be different for the hybrid FFT scheme Box realspace_bx = realspace_ba[i]; - Print() << realspace_bx.smallEnd() << " " << realspace_bx.bigEnd() << std::endl; - Box bx = Box( IntVect::TheZeroVector(), - realspace_bx.bigEnd() - realspace_bx.smallEnd() ); - spectral_bl.push_back( bx ); + IntVect fft_size = realspace_bx.length(); + // Because the spectral solver uses real-to-complex FFTs, we only + // need the positive k values along the fastest axis + // (first axis for AMReX Fortran-order arrays) in spectral space. + // This effectively reduces the size of the spectral space by half + // see e.g. the FFTW documentation for real-to-complex FFTs + IntVect spectral_bx_size = fft_size; + spectral_bx_size[0] = fft_size[0]/2 + 1; + // Define the corresponding box + Box spectral_bx = Box( IntVect::TheZeroVector(), + spectral_bx_size - IntVect::TheUnitVector() ); + spectral_bl.push_back( spectral_bx ); } spectralspace_ba.define( spectral_bl ); // Allocate the components of the k vector: kx, ky (only in 3D), kz + bool only_positive_k; for (int i_dim=0; i_dim<AMREX_SPACEDIM; i_dim++) { - k_vec[i_dim] = getKComponent( dm, i_dim ); + if (i_dim==0) { + // Real-to-complex FFTs: first axis contains only the positive k + only_positive_k = true; + } else { + only_positive_k = false; + } + k_vec[i_dim] = getKComponent(dm, realspace_ba, i_dim, only_positive_k); } } @@ -50,7 +64,9 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, */ KVectorComponent SpectralKSpace::getKComponent( const DistributionMapping& dm, - const int i_dim ) const + const BoxArray& realspace_ba, + const int i_dim, + const bool only_positive_k ) const { // Initialize an empty ManagedVector in each box KVectorComponent k_comp(spectralspace_ba, dm); @@ -65,21 +81,31 @@ SpectralKSpace::getKComponent( const DistributionMapping& dm, k.resize( N ); // Fill the k vector - const Real dk = 2*MathConst::pi/(N*dx[i_dim]); + IntVect fft_size = realspace_ba[mfi].length(); + const Real dk = 2*MathConst::pi/(fft_size[i_dim]*dx[i_dim]); AMREX_ALWAYS_ASSERT_WITH_MESSAGE( bx.smallEnd(i_dim) == 0, "Expected box to start at 0, in spectral space."); AMREX_ALWAYS_ASSERT_WITH_MESSAGE( bx.bigEnd(i_dim) == N-1, "Expected different box end index in spectral space."); - const int mid_point = (N+1)/2; - // Fill positive values of k (FFT conventions: first half is positive) - for (int i=0; i<mid_point; i++ ){ - k[i] = i*dk; - } - // Fill negative values of k (FFT conventions: second half is negative) - for (int i=mid_point; i<N; i++){ - k[i] = (i-N)*dk; + if (only_positive_k){ + // Fill the full axis with positive k values + // (typically: first axis, in a real-to-complex FFT) + for (int i=0; i<N; i++ ){ + k[i] = i*dk; + } + } else { + const int mid_point = (N+1)/2; + // Fill positive values of k + // (FFT conventions: first half is positive) + for (int i=0; i<mid_point; i++ ){ + k[i] = i*dk; + } + // Fill negative values of k + // (FFT conventions: second half is negative) + for (int i=mid_point; i<N; i++){ + k[i] = (i-N)*dk; + } } - // TODO: this will be different for the real-to-complex transform // TODO: this will be different for the hybrid FFT scheme } return k_comp; @@ -116,9 +142,14 @@ SpectralKSpace::getSpectralShiftFactor( const DistributionMapping& dm, case ShiftType::TransformFromCellCentered: sign = -1.; break; case ShiftType::TransformToCellCentered: sign = 1.; } - constexpr Complex I{0,1}; + const Complex I{0,1}; for (int i=0; i<k.size(); i++ ){ +#ifdef AMREX_USE_GPU + shift[i] = thrust::exp( I*sign*k[i]*0.5*dx[i_dim] ); +#else shift[i] = std::exp( I*sign*k[i]*0.5*dx[i_dim] ); +#endif + } } return shift_factor; @@ -161,12 +192,13 @@ SpectralKSpace::getModifiedKComponent( const DistributionMapping& dm, // Fill the modified k vector for (int i=0; i<k.size(); i++ ){ + modified_k[i] = 0; for (int n=1; n<stencil_coef.size(); n++){ if (nodal){ - modified_k[i] = stencil_coef[n]* \ + modified_k[i] += stencil_coef[n]* \ std::sin( k[i]*n*delta_x )/( n*delta_x ); } else { - modified_k[i] = stencil_coef[n]* \ + modified_k[i] += stencil_coef[n]* \ std::sin( k[i]*(n-0.5)*delta_x )/( (n-0.5)*delta_x ); } } diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index 7444452af..d4019a9a3 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -1,8 +1,7 @@ #ifndef WARPX_SPECTRAL_SOLVER_H_ #define WARPX_SPECTRAL_SOLVER_H_ -#include <SpectralKSpace.H> -#include <PsatdAlgorithm.H> +#include <SpectralBaseAlgorithm.H> #include <SpectralFieldData.H> /* \brief Top-level class for the electromagnetic spectral solver @@ -14,28 +13,17 @@ class SpectralSolver { public: - // Inline definition of the member functions of `SpectralSolver` + // Inline definition of the member functions of `SpectralSolver`, + // except the constructor (see `SpectralSolver.cpp`) // The body of these functions is short, since the work is done in the // underlying classes `SpectralFieldData` and `PsatdAlgorithm` - /* \brief Initialize the spectral solver */ + // Constructor SpectralSolver( const amrex::BoxArray& realspace_ba, const amrex::DistributionMapping& dm, const int norder_x, const int norder_y, const int norder_z, const bool nodal, - const amrex::RealVect dx, const amrex::Real dt ) { - // Initialize all structures using the same distribution mapping dm - - // - Initialize k space object (Contains info about the size of - // the spectral space corresponding to each box in `realspace_ba`, - // as well as the value of the corresponding k coordinates) - const SpectralKSpace k_space= SpectralKSpace(realspace_ba, dm, dx); - // - Initialize the algorithm (coefficients) over this space - algorithm = PsatdAlgorithm( k_space, dm, norder_x, norder_y, - norder_z, nodal, dt ); - // - Initialize arrays for fields in Fourier space + FFT plans - field_data = SpectralFieldData( realspace_ba, k_space, dm ); - }; + const amrex::RealVect dx, const amrex::Real dt ); /* \brief Transform the component `i_comp` of MultiFab `mf` * to spectral space, and store the corresponding result internally @@ -59,14 +47,20 @@ class SpectralSolver /* \brief Update the fields in spectral space, over one timestep */ void pushSpectralFields(){ BL_PROFILE("SpectralSolver::pushSpectralFields"); - algorithm.pushSpectralFields( field_data ); + // Virtual function: the actual function used here depends + // on the sub-class of `SpectralBaseAlgorithm` that was + // initialized in the constructor of `SpectralSolver` + algorithm->pushSpectralFields( field_data ); }; private: SpectralFieldData field_data; // Store field in spectral space // and perform the Fourier transforms - PsatdAlgorithm algorithm; // Contains Psatd coefficients - // and field update equation + std::unique_ptr<SpectralBaseAlgorithm> algorithm; + // Defines field update equation in spectral space, + // and the associated coefficients. + // SpectralBaseAlgorithm is a base class ; this pointer is meant + // to point an instance of a *sub-class* defining a specific algorithm }; #endif // WARPX_SPECTRAL_SOLVER_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp new file mode 100644 index 000000000..c21c3cfb1 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -0,0 +1,35 @@ +#include <SpectralKSpace.H> +#include <SpectralSolver.H> +#include <PsatdAlgorithm.H> + +/* \brief Initialize the spectral Maxwell solver + * + * This function selects the spectral algorithm to be used, allocates the + * corresponding coefficients for the discretized field update equation, + * and prepares the structures that store the fields in spectral space. + */ +SpectralSolver::SpectralSolver( + const amrex::BoxArray& realspace_ba, + const amrex::DistributionMapping& dm, + const int norder_x, const int norder_y, + const int norder_z, const bool nodal, + const amrex::RealVect dx, const amrex::Real dt ) { + + // Initialize all structures using the same distribution mapping dm + + // - Initialize k space object (Contains info about the size of + // the spectral space corresponding to each box in `realspace_ba`, + // as well as the value of the corresponding k coordinates) + const SpectralKSpace k_space= SpectralKSpace(realspace_ba, dm, dx); + + // - Select the algorithm depending on the input parameters + // Initialize the corresponding coefficients over k space + // TODO: Add more algorithms + selection depending on input parameters + // For the moment, this only uses the standard PsatdAlgorithm + algorithm = std::unique_ptr<PsatdAlgorithm>( new PsatdAlgorithm( + k_space, dm, norder_x, norder_y, norder_z, nodal, dt ) ); + + // - Initialize arrays for fields in spectral space + FFT plans + field_data = SpectralFieldData( realspace_ba, k_space, dm ); + +}; diff --git a/Source/FieldSolver/WarpXFFT.cpp b/Source/FieldSolver/WarpXFFT.cpp index 4bccc956b..13d92f6f3 100644 --- a/Source/FieldSolver/WarpXFFT.cpp +++ b/Source/FieldSolver/WarpXFFT.cpp @@ -56,6 +56,7 @@ BuildFFTOwnerMask (const MultiFab& mf, const Geometry& geom) for (const auto& b : bl) { fab.setVal(nonowner, b, 0, 1); } + } return mask; @@ -89,7 +90,7 @@ CopyDataFromFFTToValid (MultiFab& mf, const MultiFab& mf_fft, const BoxArray& ba const FArrayBox& srcfab = mf_fft[mfi]; const Box& srcbox = srcfab.box(); - + if (srcbox.contains(bx)) { // Copy the interior region (without guard cells) @@ -107,6 +108,8 @@ CopyDataFromFFTToValid (MultiFab& mf, const MultiFab& mf_fft, const BoxArray& ba // the cell that has non-zero mask is the one which is retained. mf.setVal(0.0, 0); mf.ParallelAdd(mftmp); + + } } @@ -130,7 +133,11 @@ WarpX::AllocLevelDataFFT (int lev) // Allocate and initialize objects for the spectral solver // (all use the same distribution mapping) std::array<Real,3> dx = CellSize(lev); - RealVect dx_vect = RealVect( AMREX_D_DECL(dx[0], dx[1], dx[2]) ); +#if (AMREX_SPACEDIM == 3) + RealVect dx_vect(dx[0], dx[1], dx[2]); +#elif (AMREX_SPACEDIM == 2) + RealVect dx_vect(dx[0], dx[2]); +#endif spectral_solver_fp[lev].reset( new SpectralSolver( ba_fp_fft, dm_fp_fft, nox_fft, noy_fft, noz_fft, do_nodal, dx_vect, dt[lev] ) ); } @@ -403,6 +410,7 @@ WarpX::PushPSATD (int lev, amrex::Real /* dt */) BL_PROFILE_VAR_START(blp_push_eb); if (fft_hybrid_mpi_decomposition){ +#ifndef AMREX_USE_CUDA // When running on CPU: use PICSAR code if (Efield_fp_fft[lev][0]->local_size() == 1) //Only one FFT patch on this MPI { @@ -443,6 +451,9 @@ WarpX::PushPSATD (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 } else { // Not using the hybrid decomposition auto& solver = *spectral_solver_fp[lev]; @@ -470,6 +481,7 @@ WarpX::PushPSATD (int lev, amrex::Real /* dt */) 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); @@ -486,4 +498,7 @@ WarpX::PushPSATD (int lev, amrex::Real /* dt */) { amrex::Abort("WarpX::PushPSATD: TODO"); } + } + + diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 6a9ceae5a..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<Real,3>& 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) @@ -65,6 +65,10 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) MultiFab* cost = costs[lev].get(); const IntVect& rr = (lev > 0) ? refRatio(lev-1) : IntVect::TheUnitVector(); + // xmin is only used by the picsar kernel with cylindrical geometry, + // in which case it is actually rmin. + const Real xmin = Geom(0).ProbLo(0); + // Loop through the grids, and over the tiles within each grid #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) @@ -77,10 +81,6 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) const Box& tby = mfi.tilebox(By_nodal_flag); const Box& tbz = mfi.tilebox(Bz_nodal_flag); - // xmin is only used by the picsar kernel with cylindrical geometry, - // in which case it is actually rmin. - const Real xmin = mfi.tilebox().smallEnd(0)*dx[0]; - if (do_nodal) { auto const& Bxfab = Bx->array(mfi); auto const& Byfab = By->array(mfi); @@ -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<Real,3>& dx = WarpX::CellSize(patch_level); @@ -224,6 +224,10 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) MultiFab* cost = costs[lev].get(); const IntVect& rr = (lev > 0) ? refRatio(lev-1) : IntVect::TheUnitVector(); + // xmin is only used by the picsar kernel with cylindrical geometry, + // in which case it is actually rmin. + const Real xmin = Geom(0).ProbLo(0); + // Loop through the grids, and over the tiles within each grid #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) @@ -236,10 +240,6 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) const Box& tey = mfi.tilebox(Ey_nodal_flag); const Box& tez = mfi.tilebox(Ez_nodal_flag); - // xmin is only used by the picsar kernel with cylindrical geometry, - // in which case it is actually rmin. - const Real xmin = mfi.tilebox().smallEnd(0)*dx[0]; - if (do_nodal) { auto const& Exfab = Ex->array(mfi); auto const& Eyfab = Ey->array(mfi); @@ -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<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]}; + const std::array<Real,3> 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()) { |