From f60a8212d09abb9a0e487e66c88c12170fe8a251 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 26 Apr 2019 20:30:17 -0700 Subject: Reorganize directories --- .../SpectralSolver/SpectralAlgorithms/Make.package | 5 + .../SpectralAlgorithms/PsatdAlgorithm.H | 32 ++++ .../SpectralAlgorithms/PsatdAlgorithm.cpp | 167 +++++++++++++++++++++ 3 files changed, 204 insertions(+) create mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package create mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H create mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp (limited to 'Source/FieldSolver/SpectralSolver/SpectralAlgorithms') diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package new file mode 100644 index 000000000..1caf69397 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package @@ -0,0 +1,5 @@ +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/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H new file mode 100644 index 000000000..acefcc466 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H @@ -0,0 +1,32 @@ +#ifndef WARPX_PSATD_ALGORITHM_H_ +#define WARPX_PSATD_ALGORITHM_H_ + +#include +#include + +/* \brief Class that updates the field in spectral space + * and stores the coefficients of the corresponding update equation. + */ +class PsatdAlgorithm +{ + using SpectralCoefficients = amrex::FabArray< amrex::BaseFab >; + + 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; + + 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; +}; + +#endif // WARPX_PSATD_ALGORITHM_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp new file mode 100644 index 000000000..ada7506c3 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp @@ -0,0 +1,167 @@ +#include +#include +#include + +using namespace amrex; + +/* \brief Initialize coefficients for the update equation */ +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 +{ + const BoxArray& ba = spectral_kspace.spectralspace_ba; + + // Allocate the arrays of coefficients + C_coef = SpectralCoefficients(ba, dm, 1, 0); + S_ck_coef = SpectralCoefficients(ba, dm, 1, 0); + X1_coef = SpectralCoefficients(ba, dm, 1, 0); + 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 C = C_coef[mfi].array(); + Array4 S_ck = S_ck_coef[mfi].array(); + Array4 X1 = X1_coef[mfi].array(); + Array4 X2 = X2_coef[mfi].array(); + Array4 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); + } + }); + } +}; + +/* Advance the E and B field in spectral space (stored in `f`) + * over one time step */ +void +PsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ + + // Loop over boxes + for (MFIter mfi(f.fields); mfi.isValid(); ++mfi){ + + const Box& bx = f.fields[mfi].box(); + + // Extract arrays for the fields to be updated + Array4 fields = f.fields[mfi].array(); + // Extract arrays for the coefficients + Array4 C_arr = C_coef[mfi].array(); + Array4 S_ck_arr = S_ck_coef[mfi].array(); + Array4 X1_arr = X1_coef[mfi].array(); + Array4 X2_arr = X2_coef[mfi].array(); + Array4 X3_arr = X3_coef[mfi].array(); + // Extract pointers for the k vectors + const Real* modified_kx_arr = modified_kx_vec[mfi].dataPtr(); +#if (AMREX_SPACEDIM==3) + const Real* modified_ky_arr = modified_ky_vec[mfi].dataPtr(); +#endif + const Real* modified_kz_arr = modified_kz_vec[mfi].dataPtr(); + + // Loop over indices within one box + ParallelFor(bx, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Record old values of the fields to be updated + 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 = 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 + constexpr Real c2 = PhysConst::c*PhysConst::c; + constexpr Real inv_ep0 = 1./PhysConst::ep0; + constexpr 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) + 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; + 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; + 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) + fields(i,j,k,Idx::Bx) = C*Bx_old + - S_ck*I*(ky*Ez_old - kz*Ey_old) + + X1*I*(ky*Jz - kz*Jy); + fields(i,j,k,Idx::By) = C*By_old + - S_ck*I*(kz*Ex_old - kx*Ez_old) + + X1*I*(kz*Jx - kx*Jz); + fields(i,j,k,Idx::Bz) = C*Bz_old + - S_ck*I*(kx*Ey_old - ky*Ex_old) + + X1*I*(kx*Jy - ky*Jx); + }); + } +}; -- cgit v1.2.3 From 9e8aca63877e32463a912cd6b3fafb660e898f15 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 26 Apr 2019 20:59:35 -0700 Subject: Use factory method for the solver --- .../SpectralSolver/SpectralAlgorithms/Make.package | 1 + .../SpectralAlgorithms/PsatdAlgorithm.H | 20 +++------- .../SpectralAlgorithms/PsatdAlgorithm.cpp | 11 ++---- .../SpectralAlgorithms/SpectralBaseAlgorithm.H | 44 ++++++++++++++++++++++ Source/FieldSolver/SpectralSolver/SpectralSolver.H | 13 ++++--- 5 files changed, 62 insertions(+), 27 deletions(-) create mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H (limited to 'Source/FieldSolver/SpectralSolver/SpectralAlgorithms') diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package index 1caf69397..c62c21f44 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package @@ -1,3 +1,4 @@ +CEXE_headers += SpectralBaseAlgorithm.H CEXE_headers += PsatdAlgorithm.H CEXE_sources += PsatdAlgorithm.cpp diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H index acefcc466..0487e5226 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H @@ -1,31 +1,23 @@ #ifndef WARPX_PSATD_ALGORITHM_H_ #define WARPX_PSATD_ALGORITHM_H_ -#include -#include +#include /* \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 >; - 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; + const int norder_z, const bool nodal, + const amrex::Real dt); + // Redefine update equation from base class + virtual 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/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp index ada7506c3..37892d35a 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/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; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H new file mode 100644 index 000000000..18d26e0c8 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H @@ -0,0 +1,44 @@ +#ifndef WARPX_SPECTRAL_BASE_ALGORITHM_H_ +#define WARPX_SPECTRAL_BASE_ALGORITHM_H_ + +#include +#include + +/* \brief Class that updates the field in spectral space + * and stores the coefficients of the corresponding update equation. + * TODO: Mention base class + */ +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; + + protected: // Meant to be used in the subclasses + + using SpectralCoefficients = amrex::FabArray< amrex::BaseFab >; + + // 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/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index 7444452af..5266025b6 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -31,8 +31,8 @@ class SpectralSolver // 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 ); + algorithm = std::unique_ptr( new 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 ); }; @@ -59,14 +59,17 @@ class SpectralSolver /* \brief Update the fields in spectral space, over one timestep */ void pushSpectralFields(){ BL_PROFILE("SpectralSolver::pushSpectralFields"); - algorithm.pushSpectralFields( field_data ); + 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 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_ -- cgit v1.2.3 From e4985aa4188ce6f9605171db2a1c621ecc730e26 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 1 May 2019 08:20:32 -0700 Subject: Add comments --- .../SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H | 5 ++++- Source/FieldSolver/SpectralSolver/SpectralSolver.cpp | 6 +++++- 2 files changed, 9 insertions(+), 2 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralAlgorithms') diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H index 18d26e0c8..5c662e533 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H @@ -6,7 +6,10 @@ /* \brief Class that updates the field in spectral space * and stores the coefficients of the corresponding update equation. - * TODO: Mention base class + * + * `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 { diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp index 5da39a7c7..c21c3cfb1 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -2,7 +2,11 @@ #include #include -/* \brief TODO Description +/* \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, -- cgit v1.2.3 From f07d411bc3b6110a48e70a6b69907d5f7100759d Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 1 May 2019 12:07:11 -0700 Subject: Add virtual destructor --- .../SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H | 4 ++++ 1 file changed, 4 insertions(+) (limited to 'Source/FieldSolver/SpectralSolver/SpectralAlgorithms') diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H index 5c662e533..602eb2473 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H @@ -17,6 +17,10 @@ class SpectralBaseAlgorithm // 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 -- cgit v1.2.3 From 4766b39209bf3ed2849e936c4f2dca7e437f991e Mon Sep 17 00:00:00 2001 From: Revathi Jambunathan Date: Tue, 14 May 2019 18:54:24 -0400 Subject: changed made after merging with lastest dev version --- GNUmakefile | 10 +++++----- Source/Evolve/WarpXEvolveEM.cpp | 9 --------- .../SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H | 18 +++++++++--------- .../SpectralAlgorithms/PsatdAlgorithm.cpp | 5 ++--- .../FieldSolver/SpectralSolver/SpectralFieldData.cpp | 18 ++++++++++-------- Source/FieldSolver/WarpXFFT.cpp | 11 +++++++++++ Source/Initialization/WarpXInitData.cpp | 11 ----------- Source/WarpX.cpp | 2 -- Source/main.cpp | 3 --- 9 files changed, 37 insertions(+), 50 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralAlgorithms') diff --git a/GNUmakefile b/GNUmakefile index e3ed55d80..dd80ee161 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -8,17 +8,17 @@ DEBUG = FALSE #DIM = 2 DIM = 3 -#COMP = gcc +COMP = gcc #COMP = intel -COMP = pgi +#COMP = pgi TINY_PROFILE = TRUE #PROFILE = TRUE #COMM_PROFILE = TRUE #TRACE_PROFILE = TRUE -USE_OMP = FALSE -USE_GPU = TRUE +USE_OMP = TRUE +USE_GPU = FALSE EBASE = main @@ -29,7 +29,7 @@ USE_ASCENT_INSITU = FALSE WarpxBinDir = Bin -USE_PSATD = FALSE +USE_PSATD = TRUE USE_RZ = FALSE DO_ELECTROSTATIC = FALSE diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 02c21529b..6541dd4be 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -69,20 +69,17 @@ 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; - amrex::Print() << " in evolve after pushP \n"; } else { // Beyond one step, we have E^{n} and B^{n}. @@ -92,7 +89,6 @@ WarpX::EvolveEM (int numsteps) UpdateAuxilaryData(); } - amrex::Print() << " in evolve after fill boundary \n"; if (do_subcycling == 0 || finest_level == 0) { OneStep_nosub(cur_time); @@ -102,7 +98,6 @@ 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"; if (num_mirrors>0){ applyMirrors(cur_time); @@ -274,9 +269,7 @@ WarpX::OneStep_nosub (Real cur_time) if (warpx_py_particlescraper) warpx_py_particlescraper(); if (warpx_py_beforedeposition) warpx_py_beforedeposition(); #endif - amrex::Print() << " before push \n"; PushParticlesandDepose(cur_time); - amrex::Print() << " after push \n"; #ifdef WARPX_USE_PY if (warpx_py_afterdeposition) warpx_py_afterdeposition(); @@ -290,9 +283,7 @@ WarpX::OneStep_nosub (Real cur_time) // (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/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H index 34743525e..36d5782e8 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H @@ -1,8 +1,8 @@ #ifndef WARPX_PSATD_ALGORITHM_H_ #define WARPX_PSATD_ALGORITHM_H_ -#include -#include +//#include +//#include #include /* \brief Class that updates the field in spectral space @@ -17,19 +17,19 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm 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; + //PsatdAlgorithm() = default; // Default constructor + //PsatdAlgorithm& operator=(PsatdAlgorithm&& algorithm) = default; void InitializeCoefficience(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 +// 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/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp index 8dd2a830f..3da0ef453 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp @@ -23,6 +23,7 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, X3_coef = SpectralCoefficients(ba, dm, 1, 0); InitializeCoefficience(spectral_kspace, dm, dt); +} // // Fill them with the right values: // // Loop over boxes and allocate the corresponding coefficients // // for each box owned by the local MPI proc @@ -76,7 +77,6 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, // } // }); // } -}; /* Advance the E and B field in spectral space (stored in `f`) * over one time step */ @@ -173,8 +173,7 @@ void PsatdAlgorithm::InitializeCoefficience(const SpectralKSpace& spectral_kspac // for each box owned by the local MPI proc for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi){ - //const Box& bx = ba[mfi]; - const Box bx = ba[mfi]; + const Box& bx = ba[mfi]; // Extract pointers for the k vectors const Real* modified_kx = modified_kx_vec[mfi].dataPtr(); diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 5998bdd2b..ca9e87f0f 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -152,9 +152,10 @@ 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 +#ifdef AMREX_USE_GPU + cudaDeviceSynchronize(); +#endif + amrex::Print() << " before forward transform " << mf_arr(15,15,15,0) << " " << mf_arr(15,15,15,1); } // Perform Fourier transform from `tmpRealField` to `tmpSpectralField` @@ -165,7 +166,7 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, cufftResult result; cudaStream_t stream = amrex::Cuda::Device::cudaStream(); amrex::Print() << " stream is " << stream << "\n"; - cufftSetStream ( forward_plan[mfi], stream); +// cufftSetStream ( forward_plan[mfi], stream); result = cufftExecD2Z( forward_plan[mfi], tmpRealField[mfi].dataPtr(), reinterpret_cast( @@ -268,7 +269,7 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, cufftResult result; cudaStream_t stream = amrex::Cuda::Device::cudaStream(); amrex::Print() << " stream is " << stream << "\n"; - cufftSetStream ( backward_plan[mfi], stream); +// cufftSetStream ( backward_plan[mfi], stream); result = cufftExecZ2D( backward_plan[mfi], reinterpret_cast( tmpSpectralField[mfi].dataPtr()), @@ -298,9 +299,10 @@ 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() << " mf_arr after BT " << mf_arr(15,15,15,0) << " " << mf_arr(15,15,15,1) << "\n"; #ifdef AMREX_USE_GPU cudaDeviceSynchronize(); diff --git a/Source/FieldSolver/WarpXFFT.cpp b/Source/FieldSolver/WarpXFFT.cpp index c4e0461f9..3c7127d24 100644 --- a/Source/FieldSolver/WarpXFFT.cpp +++ b/Source/FieldSolver/WarpXFFT.cpp @@ -452,16 +452,27 @@ WarpX::PushPSATD (int lev, amrex::Real /* dt */) auto& solver = *spectral_solver_fp[lev]; // Perform forward Fourier transform + amrex::Print() << " FT Ex \n"; solver.ForwardTransform(*Efield_fp_fft[lev][0], SpectralFieldIndex::Ex); + amrex::Print() << " FT Ey \n"; solver.ForwardTransform(*Efield_fp_fft[lev][1], SpectralFieldIndex::Ey); + amrex::Print() << " FT Ez \n"; solver.ForwardTransform(*Efield_fp_fft[lev][2], SpectralFieldIndex::Ez); + amrex::Print() << " FT Bx \n"; solver.ForwardTransform(*Bfield_fp_fft[lev][0], SpectralFieldIndex::Bx); + amrex::Print() << " FT By \n"; solver.ForwardTransform(*Bfield_fp_fft[lev][1], SpectralFieldIndex::By); + amrex::Print() << " FT Bz \n"; solver.ForwardTransform(*Bfield_fp_fft[lev][2], SpectralFieldIndex::Bz); + amrex::Print() << " FT jx \n"; solver.ForwardTransform(*current_fp_fft[lev][0], SpectralFieldIndex::Jx); + amrex::Print() << " FT jy \n"; solver.ForwardTransform(*current_fp_fft[lev][1], SpectralFieldIndex::Jy); + amrex::Print() << " FT jz \n"; solver.ForwardTransform(*current_fp_fft[lev][2], SpectralFieldIndex::Jz); + amrex::Print() << " FT rho old \n"; solver.ForwardTransform(*rho_fp_fft[lev], SpectralFieldIndex::rho_old, 0); + amrex::Print() << " FT rho new \n"; solver.ForwardTransform(*rho_fp_fft[lev], SpectralFieldIndex::rho_new, 1); // Advance fields in spectral space diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index b9f27f07e..6a7e79924 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -319,27 +319,16 @@ 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) { diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 83788d5b9..a3a24897a 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -556,9 +556,7 @@ 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 2e1a5e9cf..d89e89e47 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -30,13 +30,10 @@ 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 7cad9dae5a2acf76e2352d428cd2a03c88453a3f Mon Sep 17 00:00:00 2001 From: Revathi Jambunathan Date: Thu, 6 Jun 2019 14:38:11 -0400 Subject: cleaning code for spectral cufft for 3D --- .../SpectralAlgorithms/PsatdAlgorithm.H | 10 ++-- .../SpectralAlgorithms/PsatdAlgorithm.cpp | 58 ++-------------------- .../SpectralSolver/SpectralFieldData.cpp | 50 +++++++------------ .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 1 + Source/FieldSolver/WarpXFFT.cpp | 23 ++++----- 5 files changed, 35 insertions(+), 107 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralAlgorithms') diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H index 36d5782e8..0a8614e54 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H @@ -17,19 +17,15 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm 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 InitializeCoefficience(const SpectralKSpace& spectral_kspace, + + 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/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp index 3da0ef453..d45b01bda 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp @@ -22,61 +22,8 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, X2_coef = SpectralCoefficients(ba, dm, 1, 0); X3_coef = SpectralCoefficients(ba, dm, 1, 0); - InitializeCoefficience(spectral_kspace, dm, dt); + InitializeSpectralCoefficients(spectral_kspace, dm, dt); } -// // 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]; -// 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 C = C_coef[mfi].array(); -// Array4 S_ck = S_ck_coef[mfi].array(); -// Array4 X1 = X1_coef[mfi].array(); -// Array4 X2 = X2_coef[mfi].array(); -// Array4 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); -// } -// }); -// } /* Advance the E and B field in spectral space (stored in `f`) * over one time step */ @@ -139,6 +86,7 @@ PsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ const Real X2 = X2_arr(i,j,k); const Real X3 = X3_arr(i,j,k); + // Update E (see WarpX online documentation: theory section) fields(i,j,k,Idx::Ex) = C*Ex_old + S_ck*(c2*I*(ky*Bz_old - kz*By_old) - inv_ep0*Jx) @@ -163,7 +111,7 @@ PsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ } }; -void PsatdAlgorithm::InitializeCoefficience(const SpectralKSpace& spectral_kspace, +void PsatdAlgorithm::InitializeSpectralCoefficients(const SpectralKSpace& spectral_kspace, const amrex::DistributionMapping& dm, const amrex::Real dt) { diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index ca9e87f0f..cb26b19b7 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -55,23 +55,30 @@ 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; +#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"; } +#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); - // Add 2D cuFFT-specific code for Z2D if ( result != CUFFT_SUCCESS ) { amrex::Print() << " cufftplan3d backward failed! \n"; } - cudaDeviceSynchronize(); +#else + // Add 2D cuFFT-specific code for Z2D + // Note that Z2D is inherently backward plan + +#endif #else // Create FFTW plans @@ -152,29 +159,22 @@ 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() << " before forward transform " << mf_arr(15,15,15,0) << " " << mf_arr(15,15,15,1); } // 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 - cudaDeviceSynchronize(); cufftResult result; - cudaStream_t stream = amrex::Cuda::Device::cudaStream(); - amrex::Print() << " stream is " << stream << "\n"; -// cufftSetStream ( forward_plan[mfi], stream); + cudaStream_t stream = amrex::Gpu::Device::cudaStream(); + cufftSetStream ( forward_plan[mfi], stream); result = cufftExecD2Z( forward_plan[mfi], tmpRealField[mfi].dataPtr(), reinterpret_cast( tmpSpectralField[mfi].dataPtr()) ); if ( result != CUFFT_SUCCESS ) { - amrex::Print() << " cufftplan3d execute failed ! \n"; + amrex::Print() << " forward transform using cufftExecD2Z failed ! \n"; } - cudaDeviceSynchronize(); #else fftw_execute( forward_plan[mfi] ); #endif @@ -193,6 +193,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); @@ -245,6 +246,7 @@ 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(); + ParallelFor( spectralspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { Complex spectral_field_value = field_arr(i,j,k,field_index); @@ -265,22 +267,16 @@ 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(); cufftResult result; - cudaStream_t stream = amrex::Cuda::Device::cudaStream(); - amrex::Print() << " stream is " << stream << "\n"; -// cufftSetStream ( backward_plan[mfi], stream); + cudaStream_t stream = amrex::Gpu::Device::cudaStream(); + cufftSetStream ( backward_plan[mfi], stream); result = cufftExecZ2D( backward_plan[mfi], reinterpret_cast( tmpSpectralField[mfi].dataPtr()), tmpRealField[mfi].dataPtr() ); if ( result != CUFFT_SUCCESS ) { - amrex::Print() << " cufftplan3d execute inverse failed ! \n"; + amrex::Print() << " Backward transform using cufftexecZ2D failed! \n"; } - if ( result == CUFFT_SUCCESS ) { - amrex::Print() << " created cufft inverse transform\n"; - } - cudaDeviceSynchronize(); #else fftw_execute( backward_plan[mfi] ); #endif @@ -294,20 +290,12 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, Array4 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 { // Copy and normalize field mf_arr(i,j,k,i_comp) = inv_N*tmp_arr(i,j,k); }); -#ifdef AMREX_USE_GPU - cudaDeviceSynchronize(); -#endif - amrex::Print() << " mf_arr after BT " << mf_arr(15,15,15,0) << " " << mf_arr(15,15,15,1) << "\n"; - -#ifdef AMREX_USE_GPU - cudaDeviceSynchronize(); -#endif - amrex::Print() << " divided by 1/N \n"; } } } diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 6a88a52a3..6fe5e3939 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -149,6 +149,7 @@ SpectralKSpace::getSpectralShiftFactor( const DistributionMapping& dm, #else shift[i] = std::exp( I*sign*k[i]*0.5*dx[i_dim] ); #endif + } } return shift_factor; diff --git a/Source/FieldSolver/WarpXFFT.cpp b/Source/FieldSolver/WarpXFFT.cpp index d2cb83e60..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); + + } } @@ -449,34 +452,23 @@ 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.") + 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]; // Perform forward Fourier transform - amrex::Print() << " FT Ex \n"; solver.ForwardTransform(*Efield_fp_fft[lev][0], SpectralFieldIndex::Ex); - amrex::Print() << " FT Ey \n"; solver.ForwardTransform(*Efield_fp_fft[lev][1], SpectralFieldIndex::Ey); - amrex::Print() << " FT Ez \n"; solver.ForwardTransform(*Efield_fp_fft[lev][2], SpectralFieldIndex::Ez); - amrex::Print() << " FT Bx \n"; solver.ForwardTransform(*Bfield_fp_fft[lev][0], SpectralFieldIndex::Bx); - amrex::Print() << " FT By \n"; solver.ForwardTransform(*Bfield_fp_fft[lev][1], SpectralFieldIndex::By); - amrex::Print() << " FT Bz \n"; solver.ForwardTransform(*Bfield_fp_fft[lev][2], SpectralFieldIndex::Bz); - amrex::Print() << " FT jx \n"; solver.ForwardTransform(*current_fp_fft[lev][0], SpectralFieldIndex::Jx); - amrex::Print() << " FT jy \n"; solver.ForwardTransform(*current_fp_fft[lev][1], SpectralFieldIndex::Jy); - amrex::Print() << " FT jz \n"; solver.ForwardTransform(*current_fp_fft[lev][2], SpectralFieldIndex::Jz); - amrex::Print() << " FT rho old \n"; solver.ForwardTransform(*rho_fp_fft[lev], SpectralFieldIndex::rho_old, 0); - amrex::Print() << " FT rho new \n"; solver.ForwardTransform(*rho_fp_fft[lev], SpectralFieldIndex::rho_new, 1); // Advance fields in spectral space @@ -489,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); @@ -505,5 +498,7 @@ WarpX::PushPSATD (int lev, amrex::Real /* dt */) { amrex::Abort("WarpX::PushPSATD: TODO"); } - amrex::Print() << " coped data from fft to valid \n"; + } + + -- cgit v1.2.3 From 2569bcd08921227bedc7ccdb1018a5614ab31610 Mon Sep 17 00:00:00 2001 From: Revathi Jambunathan Date: Thu, 6 Jun 2019 19:19:57 -0400 Subject: Included revisions as suggested in the PR --- .../SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H | 4 ---- Source/FieldSolver/SpectralSolver/SpectralFieldData.H | 1 - Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp | 16 +++++++++------- 3 files changed, 9 insertions(+), 12 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralAlgorithms') diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H index 0a8614e54..12718e38b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H @@ -1,8 +1,6 @@ #ifndef WARPX_PSATD_ALGORITHM_H_ #define WARPX_PSATD_ALGORITHM_H_ -//#include -//#include #include /* \brief Class that updates the field in spectral space @@ -10,7 +8,6 @@ */ class PsatdAlgorithm : public SpectralBaseAlgorithm { - using SpectralCoefficients = amrex::FabArray< amrex::BaseFab >; public: PsatdAlgorithm(const SpectralKSpace& spectral_kspace, @@ -25,7 +22,6 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm void pushSpectralFields(SpectralFieldData& f) const override final; private: - // Modified finite-order vectors SpectralCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef; }; diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index 1d64817ef..7954414b8 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -25,7 +25,6 @@ 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; #else using FFTplans = amrex::LayoutData; diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 6eeb266e7..a2b695568 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -53,9 +53,11 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, // 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], @@ -70,8 +72,6 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, amrex::Print() << " cufftplan3d backward failed! \n"; } #else - // 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 ) { @@ -117,7 +117,7 @@ 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 @@ -168,8 +168,9 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, // 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); @@ -270,7 +271,8 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, // 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(); -- cgit v1.2.3