From f60a8212d09abb9a0e487e66c88c12170fe8a251 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 26 Apr 2019 20:30:17 -0700 Subject: Reorganize directories --- .../SpectralAlgorithms/PsatdAlgorithm.cpp | 167 +++++++++++++++++++++ 1 file changed, 167 insertions(+) create mode 100644 Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp (limited to 'Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp') 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/PsatdAlgorithm.cpp') 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