From 438590e9a37556ec5f0a029f0f88469e7ba3565a Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 19 Apr 2019 05:39:49 -0700 Subject: Change name of PsatdSolver to PsatdAlgorithm --- .../FieldSolver/SpectralSolver/PsatdAlgorithm.cpp | 170 +++++++++++++++++++++ 1 file changed, 170 insertions(+) create mode 100644 Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp (limited to 'Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp') diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp new file mode 100644 index 000000000..78a4ba191 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp @@ -0,0 +1,170 @@ +#include +#include +#include + +using namespace amrex; + +PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, + const DistributionMapping& dm, + const int norder_x, const int norder_y, + const int norder_z, const Real dt) +{ + const BoxArray& ba = spectral_kspace.spectralspace_ba; + + // Allocate the 1D vectors + modified_kx_vec = SpectralKVector( ba, dm ); + modified_ky_vec = SpectralKVector( ba, dm ); + modified_kz_vec = SpectralKVector( ba, dm ); + // Allocate and fill them by computing the modified vector + for ( MFIter mfi(ba, dm); mfi.isValid(); ++mfi ){ + Box bx = ba[mfi]; + ComputeModifiedKVector( + modified_kx_vec[mfi], spectral_kspace.kx_vec[mfi], + bx, spectral_kspace.dx[0], norder_x ); + ComputeModifiedKVector( + modified_ky_vec[mfi], spectral_kspace.ky_vec[mfi], + bx, spectral_kspace.dx[1], norder_y ); + ComputeModifiedKVector( + modified_kz_vec[mfi], spectral_kspace.kz_vec[mfi], + bx, spectral_kspace.dx[2], norder_z ); + } + + // 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 + 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(); + const Real* modified_ky = modified_ky_vec[mfi].dataPtr(); + 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 ) + + std::pow( modified_ky[j], 2 ) + + 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); + } + }); + } +}; + +void +PsatdAlgorithm::pushSpectralFields( SpectralData& f ) const{ + + // Loop over boxes + for ( MFIter mfi(f.Ex); mfi.isValid(); ++mfi ){ + + const Box& bx = f.Ex[mfi].box(); + + // Extract arrays for the fields to be updated + Array4 Ex_arr = f.Ex[mfi].array(); + Array4 Ey_arr = f.Ey[mfi].array(); + Array4 Ez_arr = f.Ez[mfi].array(); + Array4 Bx_arr = f.Bx[mfi].array(); + Array4 By_arr = f.By[mfi].array(); + Array4 Bz_arr = f.Bz[mfi].array(); + // Extract arrays for J and rho + Array4 Jx_arr = f.Jx[mfi].array(); + Array4 Jy_arr = f.Jy[mfi].array(); + Array4 Jz_arr = f.Jz[mfi].array(); + Array4 rho_old_arr = f.rho_old[mfi].array(); + Array4 rho_new_arr = f.rho_new[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(); + const Real* modified_ky_arr = modified_ky_vec[mfi].dataPtr(); + 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 + 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); + // 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); + // k vector values, and coefficients + const Real kx = modified_kx_arr[i]; + const Real ky = modified_ky_arr[j]; + 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 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 + + 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 + + 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 + + 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 + - S_ck*I*(ky*Ez_old - kz*Ey_old) + + X1*I*(ky*Jz - kz*Jy ); + By_arr(i,j,k) = 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 + - S_ck*I*(kx*Ey_old - ky*Ex_old) + + X1*I*(kx*Jy - ky*Jx ); + }); + } +}; -- cgit v1.2.3 From 0e9c0209744cf032da4713a1041dd0ccbaa0b1c4 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 19 Apr 2019 06:04:40 -0700 Subject: Define SpectralSolver class --- Source/FieldSolver/SpectralSolver/Make.package | 5 +- Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H | 4 +- .../FieldSolver/SpectralSolver/PsatdAlgorithm.cpp | 2 +- Source/FieldSolver/SpectralSolver/SpectralData.H | 45 ----- Source/FieldSolver/SpectralSolver/SpectralData.cpp | 184 -------------------- .../FieldSolver/SpectralSolver/SpectralFieldData.H | 46 +++++ .../SpectralSolver/SpectralFieldData.cpp | 186 +++++++++++++++++++++ 7 files changed, 238 insertions(+), 234 deletions(-) delete mode 100644 Source/FieldSolver/SpectralSolver/SpectralData.H delete mode 100644 Source/FieldSolver/SpectralSolver/SpectralData.cpp create mode 100644 Source/FieldSolver/SpectralSolver/SpectralFieldData.H create mode 100644 Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp (limited to 'Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp') diff --git a/Source/FieldSolver/SpectralSolver/Make.package b/Source/FieldSolver/SpectralSolver/Make.package index 52328b259..50127914d 100644 --- a/Source/FieldSolver/SpectralSolver/Make.package +++ b/Source/FieldSolver/SpectralSolver/Make.package @@ -1,6 +1,7 @@ CEXE_headers += WarpX_Complex.H -CEXE_headers += SpectralData.H -CEXE_sources += SpectralData.cpp +CEXE_headers += SpectralSolver.H +CEXE_headers += SpectralFieldData.H +CEXE_sources += SpectralFieldData.cpp CEXE_headers += PsatdAlgorithm.H CEXE_sources += PsatdAlgorithm.cpp CEXE_headers += SpectralKSpace.H diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H index a235315ae..7e94bdec6 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H @@ -2,7 +2,7 @@ #define WARPX_PSATD_ALGORITHM_H_ #include -#include +#include using namespace amrex; using namespace Gpu; @@ -18,7 +18,7 @@ class PsatdAlgorithm const DistributionMapping& dm, const int norder_x, const int norder_y, const int norder_z, const Real dt); - void pushSpectralFields( SpectralData& f ) const; + void pushSpectralFields( SpectralFieldData& f ) const; private: // Modified finite-order vectors diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp index 78a4ba191..acc74148e 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp @@ -84,7 +84,7 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, }; void -PsatdAlgorithm::pushSpectralFields( SpectralData& f ) const{ +PsatdAlgorithm::pushSpectralFields( SpectralFieldData& f ) const{ // Loop over boxes for ( MFIter mfi(f.Ex); mfi.isValid(); ++mfi ){ diff --git a/Source/FieldSolver/SpectralSolver/SpectralData.H b/Source/FieldSolver/SpectralSolver/SpectralData.H deleted file mode 100644 index 49764b96b..000000000 --- a/Source/FieldSolver/SpectralSolver/SpectralData.H +++ /dev/null @@ -1,45 +0,0 @@ -#ifndef WARPX_PSATD_DATA_H_ -#define WARPX_PSATD_DATA_H_ - -#include -#include - -using namespace amrex; - -// Declare type for spectral fields -using SpectralField = FabArray>; - -/* 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 }; -}; - -/* \brief Class that stores the fields in spectral space - * and performs the spectral transforms to/from real space - */ -class SpectralData -{ - friend class PsatdAlgorithm; - -#ifdef AMREX_USE_GPU -// Add cuFFT-specific code -#else - using FFTplans = LayoutData; -#endif - - public: - SpectralData( const BoxArray& realspace_ba, - const BoxArray& spectralspace_ba, - const DistributionMapping& dm ); - ~SpectralData(); - void ForwardTransform( const MultiFab& mf, const int field_index ); - void BackwardTransform( MultiFab& mf, const int field_index ); - - private: - SpectralField Ex, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz, rho_old, rho_new; - SpectralField tmpRealField, tmpSpectralField; // Store fields before/after transform - FFTplans forward_plan, backward_plan; - SpectralField& getSpectralField( const int field_index ); -}; - -#endif // WARPX_PSATD_DATA_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralData.cpp b/Source/FieldSolver/SpectralSolver/SpectralData.cpp deleted file mode 100644 index b3c902b20..000000000 --- a/Source/FieldSolver/SpectralSolver/SpectralData.cpp +++ /dev/null @@ -1,184 +0,0 @@ -#include - -using namespace amrex; - -SpectralData::SpectralData( const BoxArray& realspace_ba, - const BoxArray& spectralspace_ba, - const DistributionMapping& dm ) -{ - // 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); - - // Allocate temporary arrays - over different boxarrays - tmpRealField = SpectralField(realspace_ba, dm, 1, 0); - tmpSpectralField = SpectralField(spectralspace_ba, dm, 1, 0); - - // Allocate and initialize the FFT plans - forward_plan = FFTplans(spectralspace_ba, dm); - backward_plan = FFTplans(spectralspace_ba, dm); - for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ - Box bx = spectralspace_ba[mfi]; -#ifdef AMREX_USE_GPU - // Add cuFFT-specific code -#else - // Create FFTW plans - forward_plan[mfi] = fftw_plan_dft_3d( - // Swap dimensions: AMReX data is Fortran-order, but FFTW is C-order - bx.length(2), bx.length(1), bx.length(0), - reinterpret_cast( tmpRealField[mfi].dataPtr() ), - reinterpret_cast( tmpSpectralField[mfi].dataPtr() ), - FFTW_FORWARD, FFTW_ESTIMATE ); - backward_plan[mfi] = fftw_plan_dft_3d( - // Swap dimensions: AMReX data is Fortran-order, but FFTW is C-order - bx.length(2), bx.length(1), bx.length(0), - reinterpret_cast( tmpSpectralField[mfi].dataPtr() ), - reinterpret_cast( tmpRealField[mfi].dataPtr() ), - FFTW_BACKWARD, FFTW_ESTIMATE ); - // TODO: Add 2D code - // TODO: Do real-to-complex transform instead of complex-to-complex - // TODO: Let the user decide whether to use FFTW_ESTIMATE or FFTW_MEASURE -#endif - } -} - - -SpectralData::~SpectralData() -{ - for ( MFIter mfi(tmpRealField); mfi.isValid(); ++mfi ){ -#ifdef AMREX_USE_GPU - // Add cuFFT-specific code -#else - // Destroy FFTW plans - fftw_destroy_plan( forward_plan[mfi] ); - fftw_destroy_plan( backward_plan[mfi] ); -#endif - } -} - -/* TODO: Documentation - * Example: ForwardTransform( Efield_cp[0], SpectralFieldIndex::Ex ) - */ -void -SpectralData::ForwardTransform( const MultiFab& mf, const int field_index ) -{ - // Loop over boxes - for ( MFIter mfi(mf); mfi.isValid(); ++mfi ){ - - // Copy the real-space field `mf` to the temporary field `tmpRealField` - // This ensures that all fields have the same number of points - // before the Fourier transform. - // As a consequence, the copy discards the *last* point of `mf` - // in any direction that has *nodal* index type. - { - Box bx = mf[mfi].box(); - const Box realspace_bx = bx.enclosedCells(); // discards last point in each nodal direction - AMREX_ALWAYS_ASSERT( realspace_bx == tmpRealField[mfi].box() ); - Array4 mf_arr = mf[mfi].array(); - Array4 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); - }); - } - - // 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 -#else - fftw_execute( forward_plan[mfi] ); -#endif - - // Copy the spectral-space field `tmpSpectralField` to the appropriate field - // (specified by the input argument field_index ) - { - SpectralField& field = getSpectralField( field_index ); - Array4 field_arr = field[mfi].array(); - Array4 tmp_arr = tmpSpectralField[mfi].array(); - const Box spectralspace_bx = tmpSpectralField[mfi].box(); - ParallelFor( spectralspace_bx, - [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - field_arr(i,j,k) = tmp_arr(i,j,k); - }); - } - } -} - - -/* TODO: Documentation - */ -void -SpectralData::BackwardTransform( MultiFab& mf, const int field_index ) -{ - // Loop over boxes - for ( MFIter mfi(mf); mfi.isValid(); ++mfi ){ - - // Copy the appropriate field (specified by the input argument field_index) - // to the spectral-space field `tmpSpectralField` - { - SpectralField& field = getSpectralField( field_index ); - Array4 field_arr = field[mfi].array(); - Array4 tmp_arr = tmpSpectralField[mfi].array(); - const Box spectralspace_bx = tmpSpectralField[mfi].box(); - ParallelFor( spectralspace_bx, - [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - tmp_arr(i,j,k) = field_arr(i,j,k); - }); - } - - // Perform Fourier transform from `tmpSpectralField` to `tmpRealField` -#ifdef AMREX_USE_GPU - // Add cuFFT-specific code ; make sure that this is done on the same - // GPU stream as the above copy -#else - fftw_execute( backward_plan[mfi] ); -#endif - - // Copy the temporary field `tmpRealField` to the real-space field `mf` - // The copy does *not* fill the *last* point of `mf` - // in any direction that has *nodal* index type (but this point is - // in the guard cells and will be filled by guard cell exchange) - { - Box bx = mf[mfi].box(); - const Box realspace_bx = bx.enclosedCells(); // discards last point in each nodal direction - AMREX_ALWAYS_ASSERT( realspace_bx == tmpRealField[mfi].box() ); - Array4 mf_arr = mf[mfi].array(); - Array4 tmp_arr = tmpRealField[mfi].array(); - ParallelFor( realspace_bx, - [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - mf_arr(i,j,k) = tmp_arr(i,j,k).real(); - }); - } - } -} - - -SpectralField& -SpectralData::getSpectralField( const int field_index ) -{ - switch(field_index) - { - case SpectralFieldIndex::Ex : return Ex; - case SpectralFieldIndex::Ey : return Ey; - case SpectralFieldIndex::Ez : return Ez; - case SpectralFieldIndex::Bx : return Bx; - case SpectralFieldIndex::By : return By; - case SpectralFieldIndex::Bz : return Bz; - case SpectralFieldIndex::Jx : return Jx; - case SpectralFieldIndex::Jy : return Jy; - case SpectralFieldIndex::Jz : return Jz; - case SpectralFieldIndex::rho_old : return rho_old; - case SpectralFieldIndex::rho_new : return rho_new; - default : return tmpSpectralField; // For synthax; should not occur in practice - } -} diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H new file mode 100644 index 000000000..2f6274e40 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -0,0 +1,46 @@ +#ifndef WARPX_SPECTRAL_FIELD_DATA_H_ +#define WARPX_SPECTRAL_FIELD_DATA_H_ + +#include +#include +#include + +using namespace amrex; + +// Declare type for spectral fields +using SpectralField = FabArray>; + +/* 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 }; +}; + +/* \brief Class that stores the fields in spectral space + * and performs the spectral transforms to/from real space + */ +class SpectralFieldData +{ + friend class PsatdAlgorithm; + +#ifdef AMREX_USE_GPU +// Add cuFFT-specific code +#else + using FFTplans = LayoutData; +#endif + + public: + SpectralFieldData( const BoxArray& realspace_ba, + const SpectralKSpace& k_space, + const DistributionMapping& dm ); + ~SpectralFieldData(); + void ForwardTransform( const MultiFab& mf, const int field_index ); + void BackwardTransform( MultiFab& mf, const int field_index ); + + private: + SpectralField Ex, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz, rho_old, rho_new; + SpectralField tmpRealField, tmpSpectralField; // Store fields before/after transform + FFTplans forward_plan, backward_plan; + 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 new file mode 100644 index 000000000..25c8f438b --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -0,0 +1,186 @@ +#include + +using namespace amrex; + +SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, + const SpectralKSpace& k_space, + const DistributionMapping& dm ) +{ + 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); + + // Allocate temporary arrays - over different boxarrays + tmpRealField = SpectralField(realspace_ba, dm, 1, 0); + tmpSpectralField = SpectralField(spectralspace_ba, dm, 1, 0); + + // Allocate and initialize the FFT plans + forward_plan = FFTplans(spectralspace_ba, dm); + backward_plan = FFTplans(spectralspace_ba, dm); + for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ + Box bx = spectralspace_ba[mfi]; +#ifdef AMREX_USE_GPU + // Add cuFFT-specific code +#else + // Create FFTW plans + forward_plan[mfi] = fftw_plan_dft_3d( + // Swap dimensions: AMReX data is Fortran-order, but FFTW is C-order + bx.length(2), bx.length(1), bx.length(0), + reinterpret_cast( tmpRealField[mfi].dataPtr() ), + reinterpret_cast( tmpSpectralField[mfi].dataPtr() ), + FFTW_FORWARD, FFTW_ESTIMATE ); + backward_plan[mfi] = fftw_plan_dft_3d( + // Swap dimensions: AMReX data is Fortran-order, but FFTW is C-order + bx.length(2), bx.length(1), bx.length(0), + reinterpret_cast( tmpSpectralField[mfi].dataPtr() ), + reinterpret_cast( tmpRealField[mfi].dataPtr() ), + FFTW_BACKWARD, FFTW_ESTIMATE ); + // TODO: Add 2D code + // TODO: Do real-to-complex transform instead of complex-to-complex + // TODO: Let the user decide whether to use FFTW_ESTIMATE or FFTW_MEASURE +#endif + } +} + + +SpectralFieldData::~SpectralFieldData() +{ + for ( MFIter mfi(tmpRealField); mfi.isValid(); ++mfi ){ +#ifdef AMREX_USE_GPU + // Add cuFFT-specific code +#else + // Destroy FFTW plans + fftw_destroy_plan( forward_plan[mfi] ); + fftw_destroy_plan( backward_plan[mfi] ); +#endif + } +} + +/* TODO: Documentation + * Example: ForwardTransform( Efield_cp[0], SpectralFieldIndex::Ex ) + */ +void +SpectralFieldData::ForwardTransform( const MultiFab& mf, const int field_index ) +{ + // Loop over boxes + for ( MFIter mfi(mf); mfi.isValid(); ++mfi ){ + + // Copy the real-space field `mf` to the temporary field `tmpRealField` + // This ensures that all fields have the same number of points + // before the Fourier transform. + // As a consequence, the copy discards the *last* point of `mf` + // in any direction that has *nodal* index type. + { + Box bx = mf[mfi].box(); + const Box realspace_bx = bx.enclosedCells(); // discards last point in each nodal direction + AMREX_ALWAYS_ASSERT( realspace_bx == tmpRealField[mfi].box() ); + Array4 mf_arr = mf[mfi].array(); + Array4 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); + }); + } + + // 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 +#else + fftw_execute( forward_plan[mfi] ); +#endif + + // Copy the spectral-space field `tmpSpectralField` to the appropriate field + // (specified by the input argument field_index ) + { + SpectralField& field = getSpectralField( field_index ); + Array4 field_arr = field[mfi].array(); + Array4 tmp_arr = tmpSpectralField[mfi].array(); + const Box spectralspace_bx = tmpSpectralField[mfi].box(); + ParallelFor( spectralspace_bx, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + field_arr(i,j,k) = tmp_arr(i,j,k); + }); + } + } +} + + +/* TODO: Documentation + */ +void +SpectralFieldData::BackwardTransform( MultiFab& mf, const int field_index ) +{ + // Loop over boxes + for ( MFIter mfi(mf); mfi.isValid(); ++mfi ){ + + // Copy the appropriate field (specified by the input argument field_index) + // to the spectral-space field `tmpSpectralField` + { + SpectralField& field = getSpectralField( field_index ); + Array4 field_arr = field[mfi].array(); + Array4 tmp_arr = tmpSpectralField[mfi].array(); + const Box spectralspace_bx = tmpSpectralField[mfi].box(); + ParallelFor( spectralspace_bx, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + tmp_arr(i,j,k) = field_arr(i,j,k); + }); + } + + // Perform Fourier transform from `tmpSpectralField` to `tmpRealField` +#ifdef AMREX_USE_GPU + // Add cuFFT-specific code ; make sure that this is done on the same + // GPU stream as the above copy +#else + fftw_execute( backward_plan[mfi] ); +#endif + + // Copy the temporary field `tmpRealField` to the real-space field `mf` + // The copy does *not* fill the *last* point of `mf` + // in any direction that has *nodal* index type (but this point is + // in the guard cells and will be filled by guard cell exchange) + { + Box bx = mf[mfi].box(); + const Box realspace_bx = bx.enclosedCells(); // discards last point in each nodal direction + AMREX_ALWAYS_ASSERT( realspace_bx == tmpRealField[mfi].box() ); + Array4 mf_arr = mf[mfi].array(); + Array4 tmp_arr = tmpRealField[mfi].array(); + ParallelFor( realspace_bx, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + mf_arr(i,j,k) = tmp_arr(i,j,k).real(); + }); + } + } +} + + +SpectralField& +SpectralFieldData::getSpectralField( const int field_index ) +{ + switch(field_index) + { + case SpectralFieldIndex::Ex : return Ex; + case SpectralFieldIndex::Ey : return Ey; + case SpectralFieldIndex::Ez : return Ez; + case SpectralFieldIndex::Bx : return Bx; + case SpectralFieldIndex::By : return By; + case SpectralFieldIndex::Bz : return Bz; + case SpectralFieldIndex::Jx : return Jx; + case SpectralFieldIndex::Jy : return Jy; + case SpectralFieldIndex::Jz : return Jz; + case SpectralFieldIndex::rho_old : return rho_old; + case SpectralFieldIndex::rho_new : return rho_new; + default : return tmpSpectralField; // For synthax; should not occur in practice + } +} -- cgit v1.2.3 From d96e23d67f6b181657169453daef7c37cb62e635 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 19 Apr 2019 21:50:52 -0700 Subject: Implement spectral solver in 2D --- Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H | 5 ++++- Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp | 18 ++++++++++++++++++ .../FieldSolver/SpectralSolver/SpectralFieldData.cpp | 19 ++++++++++++------- Source/FieldSolver/SpectralSolver/SpectralKSpace.H | 11 +++++++---- Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp | 10 ++++++++-- Source/FieldSolver/SpectralSolver/SpectralSolver.H | 2 +- Source/FieldSolver/WarpXFFT.cpp | 4 +++- 7 files changed, 53 insertions(+), 16 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp') diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H index 375fa48bb..b54f5c143 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H @@ -21,7 +21,10 @@ class PsatdAlgorithm private: // Modified finite-order vectors - SpectralKVector modified_kx_vec, modified_ky_vec, modified_kz_vec; + SpectralKVector modified_kx_vec, modified_kz_vec; +#if (AMREX_SPACEDIM==3) + SpectralKVector 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/PsatdAlgorithm.cpp index acc74148e..75df01343 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp @@ -13,7 +13,9 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, // Allocate the 1D vectors modified_kx_vec = SpectralKVector( ba, dm ); +#if (AMREX_SPACEDIM==3) modified_ky_vec = SpectralKVector( ba, dm ); +#endif modified_kz_vec = SpectralKVector( ba, dm ); // Allocate and fill them by computing the modified vector for ( MFIter mfi(ba, dm); mfi.isValid(); ++mfi ){ @@ -21,12 +23,18 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, ComputeModifiedKVector( modified_kx_vec[mfi], spectral_kspace.kx_vec[mfi], bx, spectral_kspace.dx[0], norder_x ); +#if (AMREX_SPACEDIM==3) ComputeModifiedKVector( modified_ky_vec[mfi], spectral_kspace.ky_vec[mfi], bx, spectral_kspace.dx[1], norder_y ); ComputeModifiedKVector( modified_kz_vec[mfi], spectral_kspace.kz_vec[mfi], bx, spectral_kspace.dx[2], norder_z ); +#else + ComputeModifiedKVector( + modified_kz_vec[mfi], spectral_kspace.kz_vec[mfi], + bx, spectral_kspace.dx[1], norder_z ); +#endif } // Allocate the arrays of coefficients @@ -44,7 +52,9 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, // 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(); @@ -60,7 +70,9 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, // 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 @@ -112,7 +124,9 @@ PsatdAlgorithm::pushSpectralFields( SpectralFieldData& f ) const{ 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 @@ -134,7 +148,11 @@ PsatdAlgorithm::pushSpectralFields( SpectralFieldData& f ) const{ const Complex rho_new = rho_new_arr(i,j,k); // k vector values, and coefficients const Real kx = modified_kx_arr[i]; +#if (AMREX_SPACEDIM==3) const Real ky = modified_ky_arr[j]; +#else + constexpr Real ky = 0; +#endif const Real kz = modified_kz_arr[k]; constexpr Real c2 = PhysConst::c*PhysConst::c; constexpr Real inv_ep0 = 1./PhysConst::ep0; diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 595cbac16..25c08a367 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -34,19 +34,24 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, // Add cuFFT-specific code #else // Create FFTW plans - forward_plan[mfi] = fftw_plan_dft_3d( - // Swap dimensions: AMReX data is Fortran-order, but FFTW is C-order - bx.length(2), bx.length(1), bx.length(0), + forward_plan[mfi] = +#if (AMREX_SPACEDIM == 3) // Swap dimensions: AMReX data is Fortran-order, but FFTW is C-order + fftw_plan_dft_3d( bx.length(2), bx.length(1), bx.length(0), +#else + fftw_plan_dft_2d( bx.length(1), bx.length(0), +#endif reinterpret_cast( tmpRealField[mfi].dataPtr() ), reinterpret_cast( tmpSpectralField[mfi].dataPtr() ), FFTW_FORWARD, FFTW_ESTIMATE ); - backward_plan[mfi] = fftw_plan_dft_3d( - // Swap dimensions: AMReX data is Fortran-order, but FFTW is C-order - bx.length(2), bx.length(1), bx.length(0), + backward_plan[mfi] = +#if (AMREX_SPACEDIM == 3) // Swap dimensions: AMReX data is Fortran-order, but FFTW is C-order + fftw_plan_dft_3d( bx.length(2), bx.length(1), bx.length(0), +#else + fftw_plan_dft_2d( bx.length(1), bx.length(0), +#endif reinterpret_cast( tmpSpectralField[mfi].dataPtr() ), reinterpret_cast( tmpRealField[mfi].dataPtr() ), FFTW_BACKWARD, FFTW_ESTIMATE ); - // TODO: Add 2D code // TODO: Do real-to-complex transform instead of complex-to-complex // TODO: Let the user decide whether to use FFTW_ESTIMATE or FFTW_MEASURE #endif diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H index a790be94f..16ec95cf5 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H @@ -14,15 +14,18 @@ class SpectralKSpace public: SpectralKSpace( const amrex::BoxArray& realspace_ba, const amrex::DistributionMapping& dm, - const amrex::Array dx ); + const amrex::RealVect dx ); amrex::BoxArray spectralspace_ba; - SpectralKVector kx_vec, ky_vec, kz_vec; - amrex::Array dx; + SpectralKVector kx_vec, kz_vec; +#if (AMREX_SPACEDIM==3) + SpectralKVector ky_vec; +#endif + amrex::RealVect dx; }; void AllocateAndFillKvector( amrex::Gpu::ManagedVector& k, - const amrex::Box& bx, const amrex::Array dx, const int i_dim ); + const amrex::Box& bx, const amrex::RealVect dx, const int i_dim ); void ComputeModifiedKVector( amrex::Gpu::ManagedVector& modified_k, diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index d05748192..76e8aef15 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -6,7 +6,7 @@ using namespace Gpu; SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, const DistributionMapping& dm, - const Array realspace_dx ) + const RealVect realspace_dx ) { // Store the cell size dx = realspace_dx; @@ -25,20 +25,26 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, // Allocate the 1D vectors kx_vec = SpectralKVector( spectralspace_ba, dm ); +#if (AMREX_SPACEDIM == 3) ky_vec = SpectralKVector( spectralspace_ba, dm ); +#endif kz_vec = SpectralKVector( spectralspace_ba, dm ); // Initialize the values on each box for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ Box bx = spectralspace_ba[mfi]; AllocateAndFillKvector( kx_vec[mfi], bx, dx, 0 ); +#if (AMREX_SPACEDIM == 3) AllocateAndFillKvector( ky_vec[mfi], bx, dx, 1 ); AllocateAndFillKvector( kz_vec[mfi], bx, dx, 2 ); +#else + AllocateAndFillKvector( kz_vec[mfi], bx, dx, 1 ); +#endif } } void AllocateAndFillKvector( ManagedVector& k, const Box& bx, - const Array dx, const int i_dim ) + const RealVect dx, const int i_dim ) { // Alllocate k to the right size int N = bx.length( i_dim ); diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index a471666b9..06d53a2a9 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -19,7 +19,7 @@ class SpectralSolver SpectralSolver( const amrex::BoxArray& realspace_ba, const amrex::DistributionMapping& dm, const int norder_x, const int norder_y, const int norder_z, - const amrex::Array dx, const amrex::Real dt ) + const amrex::RealVect dx, const amrex::Real dt ) { // Initialize all structures using the same distribution mapping dm // - Initialize k space (Contains size of each box in spectral space, diff --git a/Source/FieldSolver/WarpXFFT.cpp b/Source/FieldSolver/WarpXFFT.cpp index 1f7fe6cd0..e655bd11f 100644 --- a/Source/FieldSolver/WarpXFFT.cpp +++ b/Source/FieldSolver/WarpXFFT.cpp @@ -129,8 +129,10 @@ WarpX::AllocLevelDataFFT (int lev) if (fft_hybrid_mpi_decomposition == false){ // Allocate and initialize objects for the spectral solver // (all use the same distribution mapping) + std::array dx = CellSize(lev); + RealVect dx_vect = RealVect( AMREX_D_DECL(dx[0], dx[1], dx[2]) ); spectral_solver_fp[lev].reset( new SpectralSolver( ba_fp_fft, dm_fp_fft, - nox_fft, noy_fft, noz_fft, CellSize(lev), dt[lev] ) ); + nox_fft, noy_fft, noz_fft, dx_vect, dt[lev] ) ); } // rho2 has one extra ghost cell, so that it's safe to deposit charge density after -- cgit v1.2.3 From 7104eee32f6f909d8a0b3abd11528b9e059a36d4 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sat, 20 Apr 2019 21:29:07 -0700 Subject: Change function interface --- .../FieldSolver/SpectralSolver/PsatdAlgorithm.cpp | 25 +---- Source/FieldSolver/SpectralSolver/SpectralKSpace.H | 23 ++--- .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 106 +++++++++++---------- 3 files changed, 67 insertions(+), 87 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp') diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp index 75df01343..0b1dbc4da 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp @@ -12,30 +12,13 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, const BoxArray& ba = spectral_kspace.spectralspace_ba; // Allocate the 1D vectors - modified_kx_vec = SpectralKVector( ba, dm ); + modified_kx_vec = spectral_kspace.AllocateAndFillModifiedKVector( dm, 0, norder_x ); #if (AMREX_SPACEDIM==3) - modified_ky_vec = SpectralKVector( ba, dm ); -#endif - modified_kz_vec = SpectralKVector( ba, dm ); - // Allocate and fill them by computing the modified vector - for ( MFIter mfi(ba, dm); mfi.isValid(); ++mfi ){ - Box bx = ba[mfi]; - ComputeModifiedKVector( - modified_kx_vec[mfi], spectral_kspace.kx_vec[mfi], - bx, spectral_kspace.dx[0], norder_x ); -#if (AMREX_SPACEDIM==3) - ComputeModifiedKVector( - modified_ky_vec[mfi], spectral_kspace.ky_vec[mfi], - bx, spectral_kspace.dx[1], norder_y ); - ComputeModifiedKVector( - modified_kz_vec[mfi], spectral_kspace.kz_vec[mfi], - bx, spectral_kspace.dx[2], norder_z ); + modified_ky_vec = spectral_kspace.AllocateAndFillModifiedKVector( dm, 1, norder_y ); + modified_kz_vec = spectral_kspace.AllocateAndFillModifiedKVector( dm, 2, norder_z ); #else - ComputeModifiedKVector( - modified_kz_vec[mfi], spectral_kspace.kz_vec[mfi], - bx, spectral_kspace.dx[1], norder_z ); + modified_kz_vec = spectral_kspace.AllocateAndFillModifiedKVector( dm, 1, norder_z ); #endif - } // Allocate the arrays of coefficients C_coef = SpectralCoefficients( ba, dm, 1, 0 ); diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H index 16ec95cf5..6589c6346 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H @@ -15,22 +15,17 @@ class SpectralKSpace SpectralKSpace( const amrex::BoxArray& realspace_ba, const amrex::DistributionMapping& dm, const amrex::RealVect dx ); + SpectralKVector& AllocateAndFillKVector( + const amrex::DistributionMapping& dm, const int i_dim ) const; + SpectralKVector& AllocateAndFillModifiedKVector( + const amrex::DistributionMapping& dm, const int i_dim, const int order ) const; amrex::BoxArray spectralspace_ba; - SpectralKVector kx_vec, kz_vec; -#if (AMREX_SPACEDIM==3) - SpectralKVector ky_vec; -#endif + + private: + amrex::Array k_vec; + // 3D: array of 3 components, corresponding to kx, ky, kz + // 2D: array of 2 components, corresponding to kx, kz amrex::RealVect dx; }; -void -AllocateAndFillKvector( amrex::Gpu::ManagedVector& k, - const amrex::Box& bx, const amrex::RealVect dx, const int i_dim ); - -void -ComputeModifiedKVector( amrex::Gpu::ManagedVector& modified_k, - const amrex::Gpu::ManagedVector& k, - const amrex::Box& bx, const amrex::Real dx, - const int norder ); - #endif diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 76e8aef15..3c600221e 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -11,7 +11,7 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, // Store the cell size dx = realspace_dx; - // Create the box array that corresponds to spectral space + // Create the box array that corresponds to spectral space BoxList spectral_bl; // Create empty box list // Loop over boxes and fill the box list for (int i=0; i < realspace_ba.size(); i++ ) { @@ -23,66 +23,68 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, } spectralspace_ba.define( spectral_bl ); - // Allocate the 1D vectors - kx_vec = SpectralKVector( spectralspace_ba, dm ); -#if (AMREX_SPACEDIM == 3) - ky_vec = SpectralKVector( spectralspace_ba, dm ); -#endif - kz_vec = SpectralKVector( spectralspace_ba, dm ); - // Initialize the values on each box - for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ - Box bx = spectralspace_ba[mfi]; - AllocateAndFillKvector( kx_vec[mfi], bx, dx, 0 ); -#if (AMREX_SPACEDIM == 3) - AllocateAndFillKvector( ky_vec[mfi], bx, dx, 1 ); - AllocateAndFillKvector( kz_vec[mfi], bx, dx, 2 ); -#else - AllocateAndFillKvector( kz_vec[mfi], bx, dx, 1 ); -#endif + // Allocate the components of the k vector: kx, ky (only in 3D), kz + for (int i_dim=0; i_dim& k, const Box& bx, - const RealVect dx, const int i_dim ) +SpectralKVector& +SpectralKSpace::AllocateAndFillKVector( const DistributionMapping& dm, const int i_dim ) const { - // Alllocate k to the right size - int N = bx.length( i_dim ); - k.resize( N ); + // Initialize an empty vector in each box + SpectralKVector k_vec = SpectralKVector(spectralspace_ba, dm); + // Loop over boxes + for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ + Box bx = spectralspace_ba[mfi]; + ManagedVector& k = k_vec[mfi]; - // Fill the k vector - const Real dk = 2*MathConst::pi/(N*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."); - // Fill positive values of k (FFT conventions: first half is positive) - for (int i=0; i<(N+1)/2; i++ ){ - k[i] = i*dk; - } - // Fill negative values of k (FFT conventions: second half is negative) - for (int i=(N+1)/2; i& modified_k, - const ManagedVector& k, - const Box& bx, const Real dx, const int norder ) +SpectralKVector& +SpectralKSpace::AllocateAndFillModifiedKVector( + const DistributionMapping& dm, const int i_dim, const int order ) const { - // Allocate modified_k to the right size - int N = k.size(); - modified_k.resize( N ); + // Initialize an empty vector in each box + SpectralKVector modified_k_vec = SpectralKVector( spectralspace_ba, dm ); + // Loop over boxes + for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ + const ManagedVector& k = k_vec[i_dim][mfi]; + ManagedVector& modified_k = modified_k_vec[mfi]; - // For now, this simply copies the infinite order k - for (int i=0; i Date: Sat, 20 Apr 2019 22:53:00 -0700 Subject: Added shift function --- .../FieldSolver/SpectralSolver/PsatdAlgorithm.cpp | 8 ++-- .../FieldSolver/SpectralSolver/SpectralFieldData.H | 3 ++ .../SpectralSolver/SpectralFieldData.cpp | 8 ++++ Source/FieldSolver/SpectralSolver/SpectralKSpace.H | 15 ++++-- .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 55 +++++++++++++++++----- 5 files changed, 69 insertions(+), 20 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp') diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp index 0b1dbc4da..95c637907 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp @@ -12,12 +12,12 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, const BoxArray& ba = spectral_kspace.spectralspace_ba; // Allocate the 1D vectors - modified_kx_vec = spectral_kspace.AllocateAndFillModifiedKVector( dm, 0, norder_x ); + modified_kx_vec = spectral_kspace.AllocateAndFillModifiedKComponent( dm, 0, norder_x ); #if (AMREX_SPACEDIM==3) - modified_ky_vec = spectral_kspace.AllocateAndFillModifiedKVector( dm, 1, norder_y ); - modified_kz_vec = spectral_kspace.AllocateAndFillModifiedKVector( dm, 2, norder_z ); + modified_ky_vec = spectral_kspace.AllocateAndFillModifiedKComponent( dm, 1, norder_y ); + modified_kz_vec = spectral_kspace.AllocateAndFillModifiedKComponent( dm, 2, norder_z ); #else - modified_kz_vec = spectral_kspace.AllocateAndFillModifiedKVector( dm, 1, norder_z ); + modified_kz_vec = spectral_kspace.AllocateAndFillModifiedKComponent( dm, 1, norder_z ); #endif // Allocate the arrays of coefficients diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index 2abe81889..052fc045e 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -42,6 +42,9 @@ class SpectralFieldData SpectralField Ex, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz, rho_old, rho_new; SpectralField tmpRealField, tmpSpectralField; // Store fields before/after transform FFTplans forward_plan, backward_plan; + // Factors that shift the fields between nodal and cell-centered, in spectral space + // (for each dimension of space) + amrex::Array shift_N2C, shift_C2N; SpectralField& getSpectralField( const int field_index ); }; diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 6b39e9fa4..78b7d7c67 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -25,6 +25,14 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, tmpRealField = SpectralField(realspace_ba, dm, 1, 0); tmpSpectralField = SpectralField(spectralspace_ba, dm, 1, 0); + // Allocate the vectors that allow to shift between nodal and cell-centered + for (int i_dim=0; i_dim #include #include using KVectorComponent = amrex::LayoutData< amrex::Gpu::ManagedVector >; +using SpectralShiftFactor = amrex::LayoutData< amrex::Gpu::ManagedVector >; + +struct ShiftType { + enum{ NodalToCentered=0, CenteredToNodal=1 }; +}; /* TODO Documentation: class represent k space, and how it is distribution * for local FFT on each MPI: k spaces are not connected. @@ -12,14 +18,17 @@ using KVectorComponent = amrex::LayoutData< amrex::Gpu::ManagedVector k_vec; diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 6e6115cec..94b1486af 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -1,5 +1,6 @@ #include #include +#include using namespace amrex; using namespace Gpu; @@ -25,19 +26,19 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, // Allocate the components of the k vector: kx, ky (only in 3D), kz for (int i_dim=0; i_dim& k = k_vec[mfi]; + ManagedVector& k = k_comp[mfi]; // Allocate k to the right size int N = bx.length( i_dim ); @@ -62,19 +63,19 @@ SpectralKSpace::AllocateAndFillKVector( const DistributionMapping& dm, const int // and distinguish the size of the local box and that of the global FFT // This will also be different for the real-to-complex transform } - return k_vec; + return k_comp; } -KVectorComponent& -SpectralKSpace::AllocateAndFillModifiedKVector( +KVectorComponent +SpectralKSpace::AllocateAndFillModifiedKComponent( const DistributionMapping& dm, const int i_dim, const int order ) const { - // Initialize an empty vector in each box - KVectorComponent modified_k_vec = KVectorComponent( spectralspace_ba, dm ); + // Initialize an empty ManagedVector in each box + KVectorComponent modified_k_comp = KVectorComponent( spectralspace_ba, dm ); // Loop over boxes for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ const ManagedVector& k = k_vec[i_dim][mfi]; - ManagedVector& modified_k = modified_k_vec[mfi]; + ManagedVector& modified_k = modified_k_comp[mfi]; // Allocate modified_k to the same size as k modified_k.resize( k.size() ); @@ -86,5 +87,33 @@ SpectralKSpace::AllocateAndFillModifiedKVector( modified_k[i] = k[i]; } } - return modified_k_vec; + return modified_k_comp; +} + +SpectralShiftFactor +SpectralKSpace::AllocateAndFillSpectralShiftFactor( + const DistributionMapping& dm, const int i_dim, const int shift_type ) const +{ + // Initialize an empty ManagedVector in each box + SpectralShiftFactor shift_factor = SpectralShiftFactor( spectralspace_ba, dm ); + // Loop over boxes + for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ + const ManagedVector& k = k_vec[i_dim][mfi]; + ManagedVector& shift = shift_factor[mfi]; + + // Allocate shift coefficients + shift.resize( k.size() ); + + // Fill the shift coefficients + Real sign = 0; + switch (shift_type){ + case ShiftType::CenteredToNodal: sign = -1.; + case ShiftType::NodalToCentered: sign = 1.; + } + constexpr Complex I{0,1}; + for (int i=0; i Date: Sun, 21 Apr 2019 07:00:59 -0700 Subject: Apply shift factors --- .../FieldSolver/SpectralSolver/PsatdAlgorithm.cpp | 8 ++-- .../FieldSolver/SpectralSolver/SpectralFieldData.H | 6 ++- .../SpectralSolver/SpectralFieldData.cpp | 45 +++++++++++++++++----- Source/FieldSolver/SpectralSolver/SpectralKSpace.H | 9 ++--- .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 8 ++-- 5 files changed, 52 insertions(+), 24 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp') diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp index 95c637907..37d3b33ee 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp @@ -12,12 +12,12 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, const BoxArray& ba = spectral_kspace.spectralspace_ba; // Allocate the 1D vectors - modified_kx_vec = spectral_kspace.AllocateAndFillModifiedKComponent( dm, 0, norder_x ); + modified_kx_vec = spectral_kspace.getModifiedKComponent( dm, 0, norder_x ); #if (AMREX_SPACEDIM==3) - modified_ky_vec = spectral_kspace.AllocateAndFillModifiedKComponent( dm, 1, norder_y ); - modified_kz_vec = spectral_kspace.AllocateAndFillModifiedKComponent( dm, 2, norder_z ); + modified_ky_vec = spectral_kspace.getModifiedKComponent( dm, 1, norder_y ); + modified_kz_vec = spectral_kspace.getModifiedKComponent( dm, 2, norder_z ); #else - modified_kz_vec = spectral_kspace.AllocateAndFillModifiedKComponent( dm, 1, norder_z ); + modified_kz_vec = spectral_kspace.getModifiedKComponent( dm, 1, norder_z ); #endif // Allocate the arrays of coefficients diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index 052fc045e..63a7c7520 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -43,8 +43,10 @@ class SpectralFieldData SpectralField tmpRealField, tmpSpectralField; // Store fields before/after transform FFTplans forward_plan, backward_plan; // Factors that shift the fields between nodal and cell-centered, in spectral space - // (for each dimension of space) - amrex::Array shift_N2C, shift_C2N; + SpectralShiftFactor xshift_N2C, xshift_C2N, zshift_N2C, zshift_C2N; +#if (AMREX_SPACEDIM==3) + SpectralShiftFactor yshift_N2C, yshift_C2N; +#endif SpectralField& getSpectralField( const int field_index ); }; diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 3fd7177e5..7d712ba03 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -25,13 +25,18 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, tmpRealField = SpectralField(realspace_ba, dm, 1, 0); tmpSpectralField = SpectralField(spectralspace_ba, dm, 1, 0); - // Allocate the vectors that allow to shift between nodal and cell-centered - for (int i_dim=0; i_dim field_arr = field[mfi].array(); Array4 tmp_arr = tmpSpectralField[mfi].array(); + const Complex* xshift_C2N_arr = xshift_C2N[mfi].dataPtr(); + const Complex* yshift_C2N_arr = yshift_C2N[mfi].dataPtr(); + const Complex* zshift_C2N_arr = zshift_C2N[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 { - tmp_arr(i,j,k) = field_arr(i,j,k); + Complex spectral_field_value = field_arr(i,j,k); + // Apply proper shift in each dimension + if (is_nodal_x==false) spectral_field_value *= xshift_C2N_arr[i]; +#if (AMREX_SPACEDIM == 3) + if (is_nodal_y==false) spectral_field_value *= yshift_C2N_arr[j]; +#endif + if (is_nodal_z==false) spectral_field_value *= zshift_C2N_arr[k]; + // Copy field into temporary array + tmp_arr(i,j,k) = spectral_field_value; }); } @@ -170,7 +196,8 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, // Normalize (divide by 1/N) since the FFT result in a factor N { Box bx = mf[mfi].box(); - const Box realspace_bx = bx.enclosedCells(); // discards last point in each nodal direction + const Box realspace_bx = bx.enclosedCells(); + // `enclosedells` discards last point in each nodal direction AMREX_ALWAYS_ASSERT( realspace_bx == tmpRealField[mfi].box() ); Array4 mf_arr = mf[mfi].array(); Array4 tmp_arr = tmpRealField[mfi].array(); diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H index 90755bc0d..2f0681c3b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H @@ -22,13 +22,12 @@ class SpectralKSpace SpectralKSpace( const amrex::BoxArray& realspace_ba, const amrex::DistributionMapping& dm, const amrex::RealVect dx ); - KVectorComponent AllocateAndFillKComponent( + KVectorComponent getKComponent( const amrex::DistributionMapping& dm, const int i_dim ) const; - KVectorComponent AllocateAndFillModifiedKComponent( + KVectorComponent getModifiedKComponent( const amrex::DistributionMapping& dm, const int i_dim, const int order ) const; - SpectralShiftFactor AllocateAndFillSpectralShiftFactor( - const amrex::DistributionMapping& dm, const int i_dim, const int order ) const; - + SpectralShiftFactor getSpectralShiftFactor( + const amrex::DistributionMapping& dm, const int i_dim, const int shift_type ) const; private: amrex::Array k_vec; diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index a7f726c8a..90c20b01d 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -26,12 +26,12 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, // Allocate the components of the k vector: kx, ky (only in 3D), kz for (int i_dim=0; i_dim Date: Mon, 22 Apr 2019 11:28:08 -0700 Subject: Use Fonberg coefficients to calculate the modified k --- Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H | 6 +- .../FieldSolver/SpectralSolver/PsatdAlgorithm.cpp | 64 +++++++++++----------- Source/FieldSolver/SpectralSolver/SpectralKSpace.H | 6 +- .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 21 +++++-- Source/FieldSolver/SpectralSolver/SpectralSolver.H | 4 +- Source/FieldSolver/WarpXFFT.cpp | 2 +- 6 files changed, 59 insertions(+), 44 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp') diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H index 2c3d8abfd..16d27cab2 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H @@ -11,13 +11,13 @@ class PsatdAlgorithm using SpectralCoefficients = amrex::FabArray< amrex::BaseFab >; public: - PsatdAlgorithm( const SpectralKSpace& spectral_kspace, + PsatdAlgorithm(const SpectralKSpace& spectral_kspace, const amrex::DistributionMapping& dm, const int norder_x, const int norder_y, - const int norder_z, const amrex::Real dt); + const int norder_z, const bool nodal, const amrex::Real dt); PsatdAlgorithm() = default; // Default constructor PsatdAlgorithm& operator=(PsatdAlgorithm&& algorithm) = default; // Default move assignment - void pushSpectralFields( SpectralFieldData& f ) const; + void pushSpectralFields(SpectralFieldData& f) const; private: // Modified finite-order vectors diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp index 37d3b33ee..5ebb7144d 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp @@ -7,29 +7,29 @@ using namespace amrex; PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, const DistributionMapping& dm, const int norder_x, const int norder_y, - const int norder_z, const Real dt) + const int norder_z, const bool nodal, const Real dt) { const BoxArray& ba = spectral_kspace.spectralspace_ba; // Allocate the 1D vectors - modified_kx_vec = spectral_kspace.getModifiedKComponent( dm, 0, norder_x ); + 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 ); - modified_kz_vec = spectral_kspace.getModifiedKComponent( dm, 2, norder_z ); + 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 ); + modified_kz_vec = spectral_kspace.getModifiedKComponent(dm, 1, norder_z, nodal); #endif // 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 ); + 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 - for ( MFIter mfi(ba, dm); mfi.isValid(); ++mfi ){ + for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi){ const Box& bx = ba[mfi]; @@ -47,26 +47,26 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, Array4 X3 = X3_coef[mfi].array(); // Loop over indices within one box - ParallelFor( bx, + 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 ) + + std::pow(modified_kx[i], 2) + #if (AMREX_SPACEDIM==3) - std::pow( modified_ky[j], 2 ) + + std::pow(modified_ky[j], 2) + #endif - std::pow( modified_kz[k], 2 ) ); + 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 ); + 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); + 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; @@ -79,10 +79,10 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, }; void -PsatdAlgorithm::pushSpectralFields( SpectralFieldData& f ) const{ +PsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ // Loop over boxes - for ( MFIter mfi(f.Ex); mfi.isValid(); ++mfi ){ + for (MFIter mfi(f.Ex); mfi.isValid(); ++mfi){ const Box& bx = f.Ex[mfi].box(); @@ -113,7 +113,7 @@ PsatdAlgorithm::pushSpectralFields( SpectralFieldData& f ) const{ const Real* modified_kz_arr = modified_kz_vec[mfi].dataPtr(); // Loop over indices within one box - ParallelFor( bx, + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Record old values of the fields to be updated @@ -148,24 +148,24 @@ PsatdAlgorithm::pushSpectralFields( SpectralFieldData& f ) const{ // Update E (see WarpX online documentation: theory section) Ex_arr(i,j,k) = C*Ex_old - + S_ck*( c2*I*(ky*Bz_old - kz*By_old) - inv_ep0*Jx ) - - I*( X2*rho_new - X3*rho_old )*kx; + + 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 - + S_ck*( c2*I*(kz*Bx_old - kx*Bz_old) - inv_ep0*Jy ) - - I*( X2*rho_new - X3*rho_old )*ky; + + 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 - + S_ck*( c2*I*(kx*By_old - ky*Bx_old) - inv_ep0*Jz ) - - I*( X2*rho_new - X3*rho_old )*kz; + + 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 - S_ck*I*(ky*Ez_old - kz*Ey_old) - + X1*I*(ky*Jz - kz*Jy ); + + X1*I*(ky*Jz - kz*Jy); By_arr(i,j,k) = C*By_old - S_ck*I*(kz*Ex_old - kx*Ez_old) - + X1*I*(kz*Jx - kx*Jz ); + + X1*I*(kz*Jx - kx*Jz); Bz_arr(i,j,k) = C*Bz_old - S_ck*I*(kx*Ey_old - ky*Ex_old) - + X1*I*(kx*Jy - ky*Jx ); + + X1*I*(kx*Jy - ky*Jx); }); } }; diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H index 2f0681c3b..fd77bb7b8 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H @@ -25,7 +25,8 @@ class SpectralKSpace KVectorComponent getKComponent( const amrex::DistributionMapping& dm, const int i_dim ) const; KVectorComponent getModifiedKComponent( - const amrex::DistributionMapping& dm, const int i_dim, const int order ) const; + const amrex::DistributionMapping& dm, const int i_dim, + const int n_order, const bool nodal ) const; SpectralShiftFactor getSpectralShiftFactor( const amrex::DistributionMapping& dm, const int i_dim, const int shift_type ) const; @@ -36,4 +37,7 @@ class SpectralKSpace amrex::RealVect dx; }; +amrex::Vector +getFonbergStencilCoefficients( const int n_order, const bool nodal ); + #endif diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 77aa7342f..111643197 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -101,8 +101,13 @@ SpectralKSpace::getModifiedKComponent( { // Initialize an empty ManagedVector in each box KVectorComponent modified_k_comp = KVectorComponent( spectralspace_ba, dm ); + + // Compute real-space stencil coefficients + Vector stencil_coef = getFonbergStencilCoefficients(n_order, nodal); + // Loop over boxes for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ + Real delta_x = dx[i_dim]; const ManagedVector& k = k_vec[i_dim][mfi]; ManagedVector& modified_k = modified_k_comp[mfi]; @@ -111,9 +116,15 @@ SpectralKSpace::getModifiedKComponent( // Fill the modified k vector for (int i=0; i +Vector getFonbergStencilCoefficients( const int n_order, const bool nodal ) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE( n_order%2 == 0, "n_order should be even."); const int m = n_order/2; - Array stencil_coef; + Vector stencil_coef; stencil_coef.resize( m+1 ); // Coefficients for nodal (a.k.a. centered) finite-difference diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index 06d53a2a9..363ac7bd8 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -19,14 +19,14 @@ class SpectralSolver SpectralSolver( const amrex::BoxArray& realspace_ba, const amrex::DistributionMapping& dm, const int norder_x, const int norder_y, const int norder_z, - const amrex::RealVect dx, const amrex::Real dt ) + const bool nodal, const amrex::RealVect dx, const amrex::Real dt ) { // Initialize all structures using the same distribution mapping dm // - Initialize k space (Contains size of each box in spectral space, // and corresponding values of the k vectors) const SpectralKSpace k_space = SpectralKSpace( realspace_ba, dm, dx ); // - Initialize algorithm (coefficients) over this space - algorithm = PsatdAlgorithm( k_space, dm, norder_x, norder_y, norder_z, dt ); + 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 ); }; diff --git a/Source/FieldSolver/WarpXFFT.cpp b/Source/FieldSolver/WarpXFFT.cpp index e655bd11f..0997643ab 100644 --- a/Source/FieldSolver/WarpXFFT.cpp +++ b/Source/FieldSolver/WarpXFFT.cpp @@ -132,7 +132,7 @@ WarpX::AllocLevelDataFFT (int lev) std::array dx = CellSize(lev); RealVect dx_vect = RealVect( AMREX_D_DECL(dx[0], dx[1], dx[2]) ); spectral_solver_fp[lev].reset( new SpectralSolver( ba_fp_fft, dm_fp_fft, - nox_fft, noy_fft, noz_fft, dx_vect, dt[lev] ) ); + nox_fft, noy_fft, noz_fft, do_nodal, dx_vect, dt[lev] ) ); } // rho2 has one extra ghost cell, so that it's safe to deposit charge density after -- cgit v1.2.3 From 4a3265daf691f0a8fa7461322c7299016fd706bf Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Tue, 23 Apr 2019 16:12:11 -0700 Subject: Add comments ; copy only to valid box --- Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H | 7 +-- .../FieldSolver/SpectralSolver/PsatdAlgorithm.cpp | 6 ++- .../FieldSolver/SpectralSolver/SpectralFieldData.H | 15 ++++-- .../SpectralSolver/SpectralFieldData.cpp | 54 +++++++++++----------- .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 51 ++++++++++++++++---- 5 files changed, 88 insertions(+), 45 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp') diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H index 16d27cab2..acefcc466 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.H @@ -4,8 +4,9 @@ #include #include -/* TODO: Write documentation -*/ +/* \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 >; @@ -16,7 +17,7 @@ class PsatdAlgorithm 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; // Default move assignment + PsatdAlgorithm& operator=(PsatdAlgorithm&& algorithm) = default; void pushSpectralFields(SpectralFieldData& f) const; private: diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp index 5ebb7144d..60e9d58c0 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp @@ -4,6 +4,7 @@ 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, @@ -28,7 +29,8 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, X3_coef = SpectralCoefficients(ba, dm, 1, 0); // Fill them with the right values: - // Loop over boxes + // 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]; @@ -78,6 +80,8 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, } }; +/* Advance the E and B field in spectral space (stored in `f`) + * over one time step */ void PsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index be1765ca0..c62513de6 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -8,20 +8,23 @@ // Declare type for spectral fields using SpectralField = amrex::FabArray< amrex::BaseFab >; -/* Fields that will be stored in spectral space */ +/* 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 }; }; -/* \brief Class that stores the fields in spectral space - * and performs the spectral transforms to/from real space +/* \brief Class that stores the fields in spectral space, and performs the + * Fourier transforms between real space and spectral space */ class SpectralFieldData { friend class PsatdAlgorithm; + // Define the FFTplans type, which holds one fft plan per box + // (plans are only initialized for the boxes that are owned by + // the local MPI rank) #ifdef AMREX_USE_GPU -// Add cuFFT-specific code + // Add cuFFT-specific code #else using FFTplans = amrex::LayoutData; #endif @@ -40,7 +43,9 @@ class SpectralFieldData private: SpectralField Ex, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz, rho_old, rho_new; - SpectralField tmpRealField, tmpSpectralField; // Store fields before/after transform + // tmpRealField and tmpSpectralField store fields + // right before/after the Fourier transform + SpectralField tmpRealField, tmpSpectralField; 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 diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index ca5e338e5..15652addc 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -2,6 +2,7 @@ using namespace amrex; +/* \brief Initialize fields in spectral space, and FFT plans */ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, const SpectralKSpace& k_space, const DistributionMapping& dm ) @@ -52,6 +53,8 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, // Allocate and initialize the FFT plans forward_plan = FFTplans(spectralspace_ba, dm); backward_plan = FFTplans(spectralspace_ba, dm); + // 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]; #ifdef AMREX_USE_GPU @@ -59,7 +62,8 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, #else // Create FFTW plans forward_plan[mfi] = -#if (AMREX_SPACEDIM == 3) // Swap dimensions: AMReX data is Fortran-order, but FFTW is C-order + // 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), #else fftw_plan_dft_2d( bx.length(1), bx.length(0), @@ -68,7 +72,8 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, reinterpret_cast( tmpSpectralField[mfi].dataPtr() ), FFTW_FORWARD, FFTW_ESTIMATE ); backward_plan[mfi] = -#if (AMREX_SPACEDIM == 3) // Swap dimensions: AMReX data is Fortran-order, but FFTW is C-order + // 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), #else fftw_plan_dft_2d( bx.length(1), bx.length(0), @@ -76,8 +81,6 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, reinterpret_cast( tmpSpectralField[mfi].dataPtr() ), reinterpret_cast( tmpRealField[mfi].dataPtr() ), FFTW_BACKWARD, FFTW_ESTIMATE ); - // TODO: Do real-to-complex transform instead of complex-to-complex - // TODO: Let the user decide whether to use FFTW_ESTIMATE or FFTW_MEASURE #endif } } @@ -98,12 +101,13 @@ SpectralFieldData::~SpectralFieldData() } } -/* TODO: Documentation - * Example: ForwardTransform( Efield_cp[0], SpectralFieldIndex::Ex ) - */ +/* \brief Transform the component `i_comp` of MultiFab `mf` + * to spectral space, and store the corresponding result internally + * (in the spectral field specified by `field_index`) */ void SpectralFieldData::ForwardTransform( const 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); @@ -123,8 +127,8 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, // As a consequence, the copy discards the *last* point of `mf` // in any direction that has *nodal* index type. { - Box bx = mf[mfi].box(); - const Box realspace_bx = bx.enclosedCells(); // discards last point in each nodal direction + Box realspace_bx = mf[mfi].box(); // Copy the box + realspace_bx.enclosedCells(); // Discard last point in nodal direction AMREX_ALWAYS_ASSERT( realspace_bx == tmpRealField[mfi].box() ); Array4 mf_arr = mf[mfi].array(); Array4 tmp_arr = tmpRealField[mfi].array(); @@ -174,8 +178,8 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, } -/* TODO: Documentation - */ +/* \brief Transform spectral field specified by `field_index` back to + * 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 ) @@ -193,9 +197,10 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, for ( MFIter mfi(mf); mfi.isValid(); ++mfi ){ // Copy the spectral-space field `tmpSpectralField` to the appropriate - // field (specified by the input argument field_index ) + // 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 field_arr = field[mfi].array(); @@ -207,6 +212,8 @@ 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); @@ -216,8 +223,8 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, 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 - tmp_arr(i,j,k) = spectral_field_value; + // Copy field into temporary array (after normalization) + tmp_arr(i,j,k) = inv_N*spectral_field_value; }); } @@ -230,22 +237,15 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, #endif // Copy the temporary field `tmpRealField` to the real-space field `mf` - // The copy does *not* fill the *last* point of `mf` - // in any direction that has *nodal* index type (but this point is - // in the guard cells and will be filled by guard cell exchange) - // Normalize (divide by 1/N) since the FFT result in a factor N + // in the *valid* part of the domain (guard cells are filled later, + // through guard cell exchange) { - Box bx = mf[mfi].box(); - const Box realspace_bx = bx.enclosedCells(); - // `enclosedells` discards last point in each nodal direction - AMREX_ALWAYS_ASSERT( realspace_bx == tmpRealField[mfi].box() ); + const Box realspace_valid_bx = mfi.validbox(); Array4 mf_arr = mf[mfi].array(); Array4 tmp_arr = tmpRealField[mfi].array(); - // For normalization: divide by the number of points in the box - const Real inv_N = 1./bx.numPts(); - ParallelFor( realspace_bx, + ParallelFor( realspace_valid_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - mf_arr(i,j,k,i_comp) = inv_N*tmp_arr(i,j,k).real(); + mf_arr(i,j,k,i_comp) = tmp_arr(i,j,k).real(); }); } } diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 90866cd3b..c82e45c5e 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -45,10 +45,9 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba, } } -/* For each box, in `spectralspace_ba`, which is owned - * by the local MPI proc (as indicated by the argument `dm`), - * compute the values of the corresponding k coordinate - * along the dimension specified by `i_dim` +/* For each box, in `spectralspace_ba`, which is owned by the local MPI rank + * (as indicated by the argument `dm`), compute the values of the + * corresponding k coordinate along the dimension specified by `i_dim` */ KVectorComponent SpectralKSpace::getKComponent( const DistributionMapping& dm, @@ -57,7 +56,7 @@ SpectralKSpace::getKComponent( const DistributionMapping& dm, // Initialize an empty ManagedVector in each box KVectorComponent k_comp = KVectorComponent(spectralspace_ba, dm); // Loop over boxes and allocate the corresponding ManagedVector - // for each box owned by the local MPI proc ("mfi.isValid") + // for each box owned by the local MPI proc for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ Box bx = spectralspace_ba[mfi]; ManagedVector& k = k_comp[mfi]; @@ -87,6 +86,15 @@ SpectralKSpace::getKComponent( const DistributionMapping& dm, return k_comp; } +/* For each box, in `spectralspace_ba`, which is owned by the local MPI rank + * (as indicated by the argument `dm`), compute the values of the + * corresponding correcting "shift" factor, along the dimension + * specified by `i_dim`. + * + * (By default, we assume the FFT is done from/to a nodal grid in real space + * It the FFT is performed from/to a cell-centered grid in real space, + * a correcting "shift" factor must be applied in spectral space.) + */ SpectralShiftFactor SpectralKSpace::getSpectralShiftFactor( const DistributionMapping& dm, const int i_dim, @@ -94,7 +102,8 @@ SpectralKSpace::getSpectralShiftFactor( const DistributionMapping& dm, { // Initialize an empty ManagedVector in each box SpectralShiftFactor shift_factor = SpectralShiftFactor( spectralspace_ba, dm ); - // Loop over boxes + // Loop over boxes and allocate the corresponding ManagedVector + // for each box owned by the local MPI proc for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ const ManagedVector& k = k_vec[i_dim][mfi]; ManagedVector& shift = shift_factor[mfi]; @@ -116,6 +125,19 @@ SpectralKSpace::getSpectralShiftFactor( const DistributionMapping& dm, return shift_factor; } +/* \brief For each box, in `spectralspace_ba`, which is owned by the local MPI + * rank (as indicated by the argument `dm`), compute the values of the + * corresponding finite-order modified k vector, along the + * dimension specified by `i_dim` + * + * The finite-order modified k vector is the spectral-space representation + * of a finite-order stencil in real space. + * + * \param n_order Order of accuracy of the stencil, in discretizing + * a spatial derivative + * \param nodal Whether the stencil is to be applied to a nodal or + staggered set of fields + */ KVectorComponent SpectralKSpace::getModifiedKComponent( const DistributionMapping& dm, const int i_dim, @@ -123,12 +145,13 @@ SpectralKSpace::getModifiedKComponent( const DistributionMapping& dm, const bool nodal ) const { // Initialize an empty ManagedVector in each box - KVectorComponent modified_k_comp = KVectorComponent( spectralspace_ba, dm ); + KVectorComponent modified_k_comp = KVectorComponent(spectralspace_ba, dm); // Compute real-space stencil coefficients Vector stencil_coef = getFonbergStencilCoefficients(n_order, nodal); - // Loop over boxes + // Loop over boxes and allocate the corresponding ManagedVector + // for each box owned by the local MPI proc for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ Real delta_x = dx[i_dim]; const ManagedVector& k = k_vec[i_dim][mfi]; @@ -153,7 +176,12 @@ SpectralKSpace::getModifiedKComponent( const DistributionMapping& dm, return modified_k_comp; } -/* TODO: Documentation: point to Fonberg paper ; explain recurrence relation +/* Returns an array of coefficients, corresponding to the weight + * of each point in a finite-difference approximation (to order `n_order`) + * of a derivative. + * + * `nodal` indicates whether this finite-difference approximation is + * taken on a nodal grid or a staggered grid. */ Vector getFonbergStencilCoefficients( const int n_order, const bool nodal ) @@ -164,6 +192,11 @@ getFonbergStencilCoefficients( const int n_order, const bool nodal ) Vector coefs; coefs.resize( m+1 ); + // Note: there are closed-form formula for these coefficients, + // but they result in an overflow when evaluated numerically. + // One way to avoid the overflow is to calculate the coefficients + // by recurrence. + // Coefficients for nodal (a.k.a. centered) finite-difference if (nodal == true) { coefs[0] = -2.; // First coefficient -- cgit v1.2.3 From 759f7e1bae3db1c2f436800189a3780a8447841e Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 26 Apr 2019 19:48:42 -0700 Subject: Cosmetic changes --- .../FieldSolver/SpectralSolver/PsatdAlgorithm.cpp | 15 ++++++------ .../FieldSolver/SpectralSolver/SpectralKSpace.cpp | 6 ++--- Source/FieldSolver/SpectralSolver/SpectralSolver.H | 11 +++++---- Source/FieldSolver/SpectralSolver/WarpX_Complex.H | 27 ---------------------- 4 files changed, 16 insertions(+), 43 deletions(-) delete mode 100644 Source/FieldSolver/SpectralSolver/WarpX_Complex.H (limited to 'Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp') diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp index 60e9d58c0..900b542be 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp @@ -9,17 +9,16 @@ 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) -{ - const BoxArray& ba = spectral_kspace.spectralspace_ba; - - // Allocate the 1D vectors - modified_kx_vec = spectral_kspace.getModifiedKComponent(dm, 0, norder_x, 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); + 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); + 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); diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 6df3ea8fa..d91891a30 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -53,7 +53,7 @@ SpectralKSpace::getKComponent( const DistributionMapping& dm, const int i_dim ) const { // Initialize an empty ManagedVector in each box - KVectorComponent k_comp = KVectorComponent(spectralspace_ba, dm); + KVectorComponent k_comp(spectralspace_ba, dm); // Loop over boxes and allocate the corresponding ManagedVector // for each box owned by the local MPI proc for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ @@ -100,7 +100,7 @@ SpectralKSpace::getSpectralShiftFactor( const DistributionMapping& dm, const int shift_type ) const { // Initialize an empty ManagedVector in each box - SpectralShiftFactor shift_factor = SpectralShiftFactor( spectralspace_ba, dm ); + SpectralShiftFactor shift_factor( spectralspace_ba, dm ); // Loop over boxes and allocate the corresponding ManagedVector // for each box owned by the local MPI proc for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){ @@ -144,7 +144,7 @@ SpectralKSpace::getModifiedKComponent( const DistributionMapping& dm, const bool nodal ) const { // Initialize an empty ManagedVector in each box - KVectorComponent modified_k_comp = KVectorComponent(spectralspace_ba, dm); + KVectorComponent modified_k_comp(spectralspace_ba, dm); // Compute real-space stencil coefficients Vector stencil_coef = getFonbergStencilCoefficients(n_order, nodal); diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index 4ed44cb9f..7444452af 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -13,11 +13,6 @@ */ class SpectralSolver { - private: - SpectralFieldData field_data; // Store field in spectral space - // and perform the Fourier transforms - PsatdAlgorithm algorithm; // Contains Psatd coefficients - // and field update equation public: // Inline definition of the member functions of `SpectralSolver` // The body of these functions is short, since the work is done in the @@ -66,6 +61,12 @@ class SpectralSolver BL_PROFILE("SpectralSolver::pushSpectralFields"); 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 }; #endif // WARPX_SPECTRAL_SOLVER_H_ diff --git a/Source/FieldSolver/SpectralSolver/WarpX_Complex.H b/Source/FieldSolver/SpectralSolver/WarpX_Complex.H deleted file mode 100644 index c898c5baa..000000000 --- a/Source/FieldSolver/SpectralSolver/WarpX_Complex.H +++ /dev/null @@ -1,27 +0,0 @@ -#ifndef WARPX_COMPLEX_H_ -#define WARPX_COMPLEX_H_ - -#include - -// Define complex type on GPU/CPU -#ifdef AMREX_USE_GPU - -#include -#include -using Complex = thrust::complex; -static_assert( sizeof(Complex) == sizeof(cuDoubleComplex), - "The complex types in WarpX and cuFFT do not match."); - -#else - -#include -#include -using Complex = std::complex; -static_assert( sizeof(Complex) == sizeof(fftw_complex), - "The complex types in WarpX and FFTW do not match."); - -#endif -static_assert(sizeof(Complex) == sizeof(amrex::Real[2]), - "Unexpected complex type."); - -#endif //WARPX_COMPLEX_H_ -- cgit v1.2.3 From 07c720803c5dbc671164dd4ea6acf7175ff6f2ce Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 26 Apr 2019 20:25:34 -0700 Subject: Fix typos --- Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) (limited to 'Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp') diff --git a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp index 900b542be..56e58bcc4 100644 --- a/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/PsatdAlgorithm.cpp @@ -10,12 +10,12 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, 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)) +: 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)) + 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)) + modified_kz_vec(spectral_kspace.getModifiedKComponent(dm,1,norder_z,nodal)) #endif { const BoxArray& ba = spectral_kspace.spectralspace_ba; -- cgit v1.2.3