diff options
Diffstat (limited to 'Source')
12 files changed, 332 insertions, 41 deletions
diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 0cf367284..469805843 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -6,6 +6,10 @@ #include <AMReX_MultiFab.H> #include <AMReX_Geometry.H> +#ifdef WARPX_USE_PSATD +#include <SpectralSolver.H> +#endif + #if (AMREX_SPACEDIM == 3) #define WRPX_PML_TO_FORTRAN(x) \ @@ -93,7 +97,9 @@ class PML public: PML (const amrex::BoxArray& ba, const amrex::DistributionMapping& dm, const amrex::Geometry* geom, const amrex::Geometry* cgeom, - int ncell, int delta, int ref_ratio, int do_dive_cleaning, int do_moving_window); + int ncell, int delta, int ref_ratio, amrex::Real dt, + int nox_fft, int noy_fft, int noz_fft, bool do_nodal, + int do_dive_cleaning, int do_moving_window); void ComputePMLFactors (amrex::Real dt); @@ -111,6 +117,10 @@ public: const MultiSigmaBox& GetMultiSigmaBox_cp () const { return *sigba_cp; } +#ifdef WARPX_USE_PSATD + void PushPSATD (); +#endif + void ExchangeB (const std::array<amrex::MultiFab*,3>& B_fp, const std::array<amrex::MultiFab*,3>& B_cp); void ExchangeE (const std::array<amrex::MultiFab*,3>& E_fp, @@ -154,6 +164,11 @@ private: std::unique_ptr<MultiSigmaBox> sigba_fp; std::unique_ptr<MultiSigmaBox> sigba_cp; +#ifdef WARPX_USE_PSATD + std::unique_ptr<SpectralSolver> spectral_solver_fp; + std::unique_ptr<SpectralSolver> spectral_solver_cp; +#endif + static amrex::BoxArray MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, int ncell); diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index f780f335c..5bcec5065 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -329,7 +329,9 @@ MultiSigmaBox::ComputePMLFactorsE (const Real* dx, Real dt) PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, const Geometry* geom, const Geometry* cgeom, - int ncell, int delta, int ref_ratio, int do_dive_cleaning, int do_moving_window) + int ncell, int delta, int ref_ratio, Real dt, + int nox_fft, int noy_fft, int noz_fft, bool do_nodal, + int do_dive_cleaning, int do_moving_window) : m_geom(geom), m_cgeom(cgeom) { @@ -343,10 +345,30 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, DistributionMapping dm{ba}; - int nge = 2; - int ngb = 2; - int ngf = (do_moving_window) ? 2 : 0; - if (WarpX::maxwell_fdtd_solver_id == 1) ngf = std::max( ngf, 1 ); + // Define the number of guard cells in each direction, for E, B, and F + IntVect nge = IntVect(AMREX_D_DECL(2, 2, 2)); + IntVect ngb = IntVect(AMREX_D_DECL(2, 2, 2)); + int ngf_int = (do_moving_window) ? 2 : 0; + if (WarpX::maxwell_fdtd_solver_id == 1) ngf_int = std::max( ngf_int, 1 ); + IntVect ngf = IntVect(AMREX_D_DECL(ngf_int, ngf_int, ngf_int)); +#ifdef WARPX_USE_PSATD + // Increase the number of guard cells, in order to fit the extent + // of the stencil for the spectral solver + IntVect ngFFT; + if (do_nodal) { + ngFFT = IntVect(AMREX_D_DECL(nox_fft, noy_fft, noz_fft)); + } else { + ngFFT = IntVect(AMREX_D_DECL(nox_fft/2, noy_fft/2, noz_fft/2)); + } + // Set the number of guard cells to the maximum of each field + // (all fields should have the same number of guard cells) + ngFFT = ngFFT.max(nge); + ngFFT = ngFFT.max(ngb); + ngFFT = ngFFT.max(ngf); + nge = ngFFT; + ngb = ngFFT; + ngf = ngFFT; + #endif pml_E_fp[0].reset(new MultiFab(amrex::convert(ba,WarpX::Ex_nodal_flag), dm, 3, nge)); pml_E_fp[1].reset(new MultiFab(amrex::convert(ba,WarpX::Ey_nodal_flag), dm, 3, nge)); @@ -370,11 +392,22 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, sigba_fp.reset(new MultiSigmaBox(ba, dm, grid_ba, geom->CellSize(), ncell, delta)); +#ifdef WARPX_USE_PSATD + const bool in_pml = true; // Tells spectral solver to use split-PML equations + const RealVect dx{AMREX_D_DECL(geom->CellSize(0), geom->CellSize(1), geom->CellSize(2))}; + // Get the cell-centered box, with guard cells + BoxArray realspace_ba = ba; // Copy box + realspace_ba.enclosedCells().grow(nge); // cell-centered + guard cells + spectral_solver_fp.reset( new SpectralSolver( realspace_ba, dm, + nox_fft, noy_fft, noz_fft, do_nodal, dx, dt, in_pml ) ); +#endif + if (cgeom) { - - nge = 1; - ngb = 1; +#ifndef WARPX_USE_PSATD + nge = IntVect(AMREX_D_DECL(1, 1, 1)); + ngb = IntVect(AMREX_D_DECL(1, 1, 1)); +#endif BoxArray grid_cba = grid_ba; grid_cba.coarsen(ref_ratio); @@ -403,8 +436,17 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, } sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta)); - } +#ifdef WARPX_USE_PSATD + const bool in_pml = true; // Tells spectral solver to use split-PML equations + const RealVect cdx{AMREX_D_DECL(cgeom->CellSize(0), cgeom->CellSize(1), cgeom->CellSize(2))}; + // Get the cell-centered box, with guard cells + BoxArray realspace_ba = cba; // Copy box + realspace_ba.enclosedCells().grow(nge); // cell-centered + guard cells + spectral_solver_cp.reset( new SpectralSolver( realspace_ba, cdm, + nox_fft, noy_fft, noz_fft, do_nodal, cdx, dt, in_pml ) ); +#endif + } } BoxArray @@ -753,3 +795,41 @@ PML::Restart (const std::string& dir) VisMF::Read(*pml_B_cp[2], dir+"_Bz_cp"); } } + +#ifdef WARPX_USE_PSATD +void +PML::PushPSATD() { + SpectralSolver& solver = *(spectral_solver_fp); + + using Idx = SpectralPMLIndex; + + // Perform forward Fourier transform + solver.ForwardTransform(*pml_E_fp[0], Idx::Exy, 0); + solver.ForwardTransform(*pml_E_fp[0], Idx::Exz, 1); + solver.ForwardTransform(*pml_E_fp[1], Idx::Eyx, 0); + solver.ForwardTransform(*pml_E_fp[1], Idx::Eyz, 1); + solver.ForwardTransform(*pml_E_fp[2], Idx::Ezx, 0); + solver.ForwardTransform(*pml_E_fp[2], Idx::Ezy, 1); + solver.ForwardTransform(*pml_B_fp[0], Idx::Bxy, 0); + solver.ForwardTransform(*pml_B_fp[0], Idx::Bxz, 1); + solver.ForwardTransform(*pml_B_fp[1], Idx::Byx, 0); + solver.ForwardTransform(*pml_B_fp[1], Idx::Byz, 1); + solver.ForwardTransform(*pml_B_fp[2], Idx::Bzx, 0); + solver.ForwardTransform(*pml_B_fp[2], Idx::Bzy, 1); + // Advance fields in spectral space + solver.pushSpectralFields(); + // Perform backward Fourier Transform + solver.BackwardTransform(*pml_E_fp[0], Idx::Exy, 0); + solver.BackwardTransform(*pml_E_fp[0], Idx::Exz, 1); + solver.BackwardTransform(*pml_E_fp[1], Idx::Eyx, 0); + solver.BackwardTransform(*pml_E_fp[1], Idx::Eyz, 1); + solver.BackwardTransform(*pml_E_fp[2], Idx::Ezx, 0); + solver.BackwardTransform(*pml_E_fp[2], Idx::Ezy, 1); + solver.BackwardTransform(*pml_B_fp[0], Idx::Bxy, 0); + solver.BackwardTransform(*pml_B_fp[0], Idx::Bxz, 1); + solver.BackwardTransform(*pml_B_fp[1], Idx::Byx, 0); + solver.BackwardTransform(*pml_B_fp[1], Idx::Byz, 1); + solver.BackwardTransform(*pml_B_fp[2], Idx::Bzx, 0); + solver.BackwardTransform(*pml_B_fp[2], Idx::Bzy, 1); +} +#endif diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package index c62c21f44..ee8376865 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package @@ -1,6 +1,8 @@ CEXE_headers += SpectralBaseAlgorithm.H CEXE_headers += PsatdAlgorithm.H CEXE_sources += PsatdAlgorithm.cpp +CEXE_headers += PMLPsatdAlgorithm.H +CEXE_sources += PMLPsatdAlgorithm.cpp INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms VPATH_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H new file mode 100644 index 000000000..d77597d53 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H @@ -0,0 +1,27 @@ +#ifndef WARPX_PML_PSATD_ALGORITHM_H_ +#define WARPX_PML_PSATD_ALGORITHM_H_ + +#include <SpectralBaseAlgorithm.H> + +/* \brief Class that updates the field in spectral space + * and stores the coefficients of the corresponding update equation. + */ +class PMLPsatdAlgorithm : public SpectralBaseAlgorithm +{ + public: + PMLPsatdAlgorithm(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); + // Redefine functions from base class + virtual void pushSpectralFields(SpectralFieldData& f) const override final; + virtual int getRequiredNumberOfFields() const override final { + return SpectralPMLIndex::n_fields; + } + + private: + SpectralCoefficients C_coef, S_ck_coef; +}; + +#endif // WARPX_PML_PSATD_ALGORITHM_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp new file mode 100644 index 000000000..46ac851d2 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp @@ -0,0 +1,138 @@ +#include <PMLPsatdAlgorithm.H> +#include <WarpXConst.H> +#include <cmath> + +using namespace amrex; + +/* \brief Initialize coefficients for the update equation */ +PMLPsatdAlgorithm::PMLPsatdAlgorithm( + 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) + // Initialize members of base class + : SpectralBaseAlgorithm( spectral_kspace, dm, + norder_x, norder_y, norder_z, nodal ) +{ + 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); + + // Fill them with the right values: + // Loop over boxes and allocate the corresponding coefficients + // for each box owned by the local MPI proc + for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi){ + + const Box& bx = ba[mfi]; + + // Extract pointers for the k vectors + const Real* modified_kx = modified_kx_vec[mfi].dataPtr(); +#if (AMREX_SPACEDIM==3) + const Real* modified_ky = modified_ky_vec[mfi].dataPtr(); +#endif + const Real* modified_kz = modified_kz_vec[mfi].dataPtr(); + // Extract arrays for the coefficients + Array4<Real> C = C_coef[mfi].array(); + Array4<Real> S_ck = S_ck_coef[mfi].array(); + + // 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); + } else { // Handle k_norm = 0, by using the analytical limit + C(i,j,k) = 1.; + S_ck(i,j,k) = dt; + } + }); + } +}; + +/* Advance the E and B field in spectral space (stored in `f`) + * over one time step */ +void +PMLPsatdAlgorithm::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<Complex> fields = f.fields[mfi].array(); + // Extract arrays for the coefficients + Array4<const Real> C_arr = C_coef[mfi].array(); + Array4<const Real> S_ck_arr = S_ck_coef[mfi].array(); + // 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 = SpectralPMLIndex; + const Complex Ex_old = fields(i,j,k,Idx::Exy) \ + + fields(i,j,k,Idx::Exz); + const Complex Ey_old = fields(i,j,k,Idx::Eyx) \ + + fields(i,j,k,Idx::Eyz); + const Complex Ez_old = fields(i,j,k,Idx::Ezx) \ + + fields(i,j,k,Idx::Ezy); + const Complex Bx_old = fields(i,j,k,Idx::Bxy) \ + + fields(i,j,k,Idx::Bxz); + const Complex By_old = fields(i,j,k,Idx::Byx) \ + + fields(i,j,k,Idx::Byz); + const Complex Bz_old = fields(i,j,k,Idx::Bzx) \ + + fields(i,j,k,Idx::Bzy); + // 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 Complex I = Complex{0,1}; + const Real C = C_arr(i,j,k); + const Real S_ck = S_ck_arr(i,j,k); + + // Update E + fields(i,j,k,Idx::Exy) = C*fields(i,j,k,Idx::Exy) + S_ck*c2*I*ky*Bz_old; + fields(i,j,k,Idx::Exz) = C*fields(i,j,k,Idx::Exz) - S_ck*c2*I*kz*By_old; + fields(i,j,k,Idx::Eyz) = C*fields(i,j,k,Idx::Eyz) + S_ck*c2*I*kz*Bx_old; + fields(i,j,k,Idx::Eyx) = C*fields(i,j,k,Idx::Eyx) - S_ck*c2*I*kx*Bz_old; + fields(i,j,k,Idx::Ezx) = C*fields(i,j,k,Idx::Ezx) + S_ck*c2*I*kx*By_old; + fields(i,j,k,Idx::Ezy) = C*fields(i,j,k,Idx::Ezy) - S_ck*c2*I*ky*Bx_old; + // Update B + fields(i,j,k,Idx::Bxy) = C*fields(i,j,k,Idx::Bxy) - S_ck*I*ky*Ez_old; + fields(i,j,k,Idx::Bxz) = C*fields(i,j,k,Idx::Bxz) + S_ck*I*kz*Ey_old; + fields(i,j,k,Idx::Byz) = C*fields(i,j,k,Idx::Byz) - S_ck*I*kz*Ex_old; + fields(i,j,k,Idx::Byx) = C*fields(i,j,k,Idx::Byx) + S_ck*I*kx*Ez_old; + fields(i,j,k,Idx::Bzx) = C*fields(i,j,k,Idx::Bzx) - S_ck*I*kx*Ey_old; + fields(i,j,k,Idx::Bzy) = C*fields(i,j,k,Idx::Bzy) + S_ck*I*ky*Ex_old; + }); + } +}; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H index 12718e38b..825d04dc2 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H @@ -13,14 +13,18 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm 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); - + const int norder_z, const bool nodal, + const amrex::Real dt); + // Redefine functions from base class + virtual void pushSpectralFields(SpectralFieldData& f) const override final; + virtual int getRequiredNumberOfFields() const override final { + return SpectralFieldIndex::n_fields; + } + void InitializeSpectralCoefficients(const SpectralKSpace& spectral_kspace, - const amrex::DistributionMapping& dm, + const amrex::DistributionMapping& dm, const amrex::Real dt); - void pushSpectralFields(SpectralFieldData& f) const override final; - private: SpectralCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef; }; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H index 602eb2473..5d5e376c1 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H @@ -14,9 +14,9 @@ class SpectralBaseAlgorithm { public: - // Member function that updates the fields in spectral space ; - // meant to be overridden in subclasses + // Virtual member function ; meant to be overridden in subclasses virtual void pushSpectralFields(SpectralFieldData& f) const = 0; + virtual int getRequiredNumberOfFields() const = 0; // The destructor should also be a virtual function, so that // a pointer to subclass of `SpectraBaseAlgorithm` actually // calls the subclass's destructor. diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index 7954414b8..6a2446981 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -8,18 +8,24 @@ // Declare type for spectral fields using SpectralField = amrex::FabArray< amrex::BaseFab <Complex> >; -/* Index for the fields that will be stored in spectral space */ +/* Index for the regular fields, when stored in spectral space */ struct SpectralFieldIndex { enum { Ex=0, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz, rho_old, rho_new, n_fields }; // n_fields is automatically the total number of fields }; +/* Index for the PML fields, when stored in spectral space */ +struct SpectralPMLIndex { + enum { Exy=0, Exz, Eyx, Eyz, Ezx, Ezy, + Bxy, Bxz, Byx, Byz, Bzx, Bzy, n_fields }; + // n_fields is automatically the total number of fields +}; + /* \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 @@ -32,8 +38,9 @@ class SpectralFieldData public: SpectralFieldData( const amrex::BoxArray& realspace_ba, - const SpectralKSpace& k_space, - const amrex::DistributionMapping& dm ); + const SpectralKSpace& k_space, + const amrex::DistributionMapping& dm, + const int n_field_required ); SpectralFieldData() = default; // Default constructor SpectralFieldData& operator=(SpectralFieldData&& field_data) = default; ~SpectralFieldData(); @@ -41,10 +48,10 @@ class SpectralFieldData const int field_index, const int i_comp); void BackwardTransform( amrex::MultiFab& mf, const int field_index, const int i_comp); - - private: // `fields` stores fields in spectral space, as multicomponent FabArray SpectralField fields; + + private: // tmpRealField and tmpSpectralField store fields // right before/after the Fourier transform SpectralField tmpSpectralField; // contains Complexs diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 948baf0a6..8f0853484 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -5,14 +5,14 @@ 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 ) + const DistributionMapping& dm, + const int n_field_required ) { const BoxArray& spectralspace_ba = k_space.spectralspace_ba; // Allocate the arrays that contain the fields in spectral space // (one component per field) - fields = SpectralField(spectralspace_ba, dm, - SpectralFieldIndex::n_fields, 0); + fields = SpectralField(spectralspace_ba, dm, n_field_required, 0); // Allocate temporary arrays - in real space and spectral space // These arrays will store the data just before/after the FFT diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index d4019a9a3..c570b017b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -23,7 +23,8 @@ class SpectralSolver const amrex::DistributionMapping& dm, const int norder_x, const int norder_y, const int norder_z, const bool nodal, - const amrex::RealVect dx, const amrex::Real dt ); + const amrex::RealVect dx, const amrex::Real dt, + const bool pml=false ); /* \brief Transform the component `i_comp` of MultiFab `mf` * to spectral space, and store the corresponding result internally diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp index c21c3cfb1..4b9def013 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -1,19 +1,29 @@ #include <SpectralKSpace.H> #include <SpectralSolver.H> #include <PsatdAlgorithm.H> +#include <PMLPsatdAlgorithm.H> /* \brief Initialize the spectral Maxwell solver * * This function selects the spectral algorithm to be used, allocates the * corresponding coefficients for the discretized field update equation, * and prepares the structures that store the fields in spectral space. + * + * \param norder_x Order of accuracy of the spatial derivatives along x + * \param norder_y Order of accuracy of the spatial derivatives along y + * \param norder_z Order of accuracy of the spatial derivatives along z + * \param nodal Whether the solver is applied to a nodal or staggered grid + * \param dx Cell size along each dimension + * \param dt Time step + * \param pml Whether the boxes in which the solver is applied are PML boxes */ SpectralSolver::SpectralSolver( const amrex::BoxArray& realspace_ba, const amrex::DistributionMapping& dm, const int norder_x, const int norder_y, const int norder_z, const bool nodal, - const amrex::RealVect dx, const amrex::Real dt ) { + const amrex::RealVect dx, const amrex::Real dt, + const bool pml ) { // Initialize all structures using the same distribution mapping dm @@ -24,12 +34,16 @@ SpectralSolver::SpectralSolver( // - Select the algorithm depending on the input parameters // Initialize the corresponding coefficients over k space - // TODO: Add more algorithms + selection depending on input parameters - // For the moment, this only uses the standard PsatdAlgorithm - algorithm = std::unique_ptr<PsatdAlgorithm>( new PsatdAlgorithm( + if (pml) { + algorithm = std::unique_ptr<PMLPsatdAlgorithm>( new PMLPsatdAlgorithm( + k_space, dm, norder_x, norder_y, norder_z, nodal, dt ) ); + } else { + algorithm = std::unique_ptr<PsatdAlgorithm>( new PsatdAlgorithm( k_space, dm, norder_x, norder_y, norder_z, nodal, dt ) ); + } // - Initialize arrays for fields in spectral space + FFT plans - field_data = SpectralFieldData( realspace_ba, k_space, dm ); + field_data = SpectralFieldData( realspace_ba, k_space, dm, + algorithm->getRequiredNumberOfFields() ); }; diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 2442e0205..6f27787e8 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -88,7 +88,7 @@ WarpX::InitDiagnostics () { const Real* current_lo = geom[0].ProbLo(); const Real* current_hi = geom[0].ProbHi(); Real dt_boost = dt[0]; - + // Find the positions of the lab-frame box that corresponds to the boosted-frame box at t=0 Real zmin_lab = current_lo[moving_window_dir]/( (1.+beta_boost)*gamma_boost ); Real zmax_lab = current_hi[moving_window_dir]/( (1.+beta_boost)*gamma_boost ); @@ -97,7 +97,7 @@ WarpX::InitDiagnostics () { zmax_lab, moving_window_v, dt_snapshots_lab, num_snapshots_lab, gamma_boost, - t_new[0], dt_boost, + t_new[0], dt_boost, moving_window_dir, geom[0])); } } @@ -118,10 +118,10 @@ WarpX::InitFromScratch () InitPML(); -#ifdef WARPX_DO_ELECTROSTATIC +#ifdef WARPX_DO_ELECTROSTATIC if (do_electrostatic) { getLevelMasks(masks); - + // the plus one is to convert from num_cells to num_nodes getLevelMasks(gather_masks, n_buffer + 1); } @@ -134,13 +134,16 @@ WarpX::InitPML () if (do_pml) { pml[0].reset(new PML(boxArray(0), DistributionMap(0), &Geom(0), nullptr, - pml_ncell, pml_delta, 0, do_dive_cleaning, do_moving_window)); + pml_ncell, pml_delta, 0, dt[0], + nox_fft, noy_fft, noz_fft, do_nodal, + do_dive_cleaning, do_moving_window)); for (int lev = 1; lev <= finest_level; ++lev) { pml[lev].reset(new PML(boxArray(lev), DistributionMap(lev), &Geom(lev), &Geom(lev-1), - pml_ncell, pml_delta, refRatio(lev-1)[0], do_dive_cleaning, - do_moving_window)); + pml_ncell, pml_delta, refRatio(lev-1)[0], dt[lev], + nox_fft, noy_fft, noz_fft, do_nodal, + do_dive_cleaning, do_moving_window)); } } } @@ -226,7 +229,7 @@ WarpX::InitOpenbc () Vector<int> alllohi(6*nprocs,100000); MPI_Allgather(lohi, 6, MPI_INT, alllohi.data(), 6, MPI_INT, ParallelDescriptor::Communicator()); - + BoxList bl{IndexType::TheNodeType()}; for (int i = 0; i < nprocs; ++i) { @@ -252,7 +255,7 @@ WarpX::InitOpenbc () rho_openbc.copy(*rho, 0, 0, 1, rho->nGrow(), 0, gm.periodicity(), FabArrayBase::ADD); const Real* dx = gm.CellSize(); - + warpx_openbc_potential(rho_openbc[myproc].dataPtr(), phi_openbc[myproc].dataPtr(), dx); BoxArray nba = boxArray(lev); |