diff options
Diffstat (limited to 'Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp')
| -rw-r--r-- | Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp | 206 |
1 files changed, 102 insertions, 104 deletions
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp index e7459d0f9..14ffc2d91 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp @@ -31,12 +31,14 @@ using namespace amrex; /* \brief Initialize coefficients for the update equation */ PMLPsatdAlgorithm::PMLPsatdAlgorithm(const SpectralKSpace& spectral_kspace, const DistributionMapping& dm, + const SpectralFieldIndex& spectral_index, const int norder_x, const int norder_y, const int norder_z, const bool nodal, const amrex::IntVect& fill_guards, const Real dt, const bool dive_cleaning, const bool divb_cleaning) // Initialize members of base class - : SpectralBaseAlgorithm(spectral_kspace, dm, norder_x, norder_y, norder_z, nodal, fill_guards), + : SpectralBaseAlgorithm(spectral_kspace, dm, spectral_index, norder_x, norder_y, norder_z, nodal, fill_guards), + m_spectral_index(spectral_index), m_dt(dt), m_dive_cleaning(dive_cleaning), m_divb_cleaning(divb_cleaning) @@ -59,6 +61,8 @@ PMLPsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const { const bool dive_cleaning = m_dive_cleaning; const bool divb_cleaning = m_divb_cleaning; + const SpectralFieldIndex& Idx = m_spectral_index; + // Loop over boxes for (MFIter mfi(f.fields); mfi.isValid(); ++mfi){ @@ -85,21 +89,19 @@ PMLPsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const { 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 Exy = fields(i,j,k,Idx::Exy); - const Complex Exz = fields(i,j,k,Idx::Exz); - const Complex Eyx = fields(i,j,k,Idx::Eyx); - const Complex Eyz = fields(i,j,k,Idx::Eyz); - const Complex Ezx = fields(i,j,k,Idx::Ezx); - const Complex Ezy = fields(i,j,k,Idx::Ezy); - - const Complex Bxy = fields(i,j,k,Idx::Bxy); - const Complex Bxz = fields(i,j,k,Idx::Bxz); - const Complex Byx = fields(i,j,k,Idx::Byx); - const Complex Byz = fields(i,j,k,Idx::Byz); - const Complex Bzx = fields(i,j,k,Idx::Bzx); - const Complex Bzy = fields(i,j,k,Idx::Bzy); + const Complex Exy = fields(i,j,k,Idx.Exy); + const Complex Exz = fields(i,j,k,Idx.Exz); + const Complex Eyx = fields(i,j,k,Idx.Eyx); + const Complex Eyz = fields(i,j,k,Idx.Eyz); + const Complex Ezx = fields(i,j,k,Idx.Ezx); + const Complex Ezy = fields(i,j,k,Idx.Ezy); + + const Complex Bxy = fields(i,j,k,Idx.Bxy); + const Complex Bxz = fields(i,j,k,Idx.Bxz); + const Complex Byx = fields(i,j,k,Idx.Byx); + const Complex Byz = fields(i,j,k,Idx.Byz); + const Complex Bzx = fields(i,j,k,Idx.Bzx); + const Complex Bzy = fields(i,j,k,Idx.Bzy); Complex Ex, Ey, Ez; Complex Bx, By, Bz; @@ -123,21 +125,21 @@ PMLPsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const { } else if (dive_cleaning && divb_cleaning) { - Exx = fields(i,j,k,Idx::Exx); - Eyy = fields(i,j,k,Idx::Eyy); - Ezz = fields(i,j,k,Idx::Ezz); + Exx = fields(i,j,k,Idx.Exx); + Eyy = fields(i,j,k,Idx.Eyy); + Ezz = fields(i,j,k,Idx.Ezz); - Bxx = fields(i,j,k,Idx::Bxx); - Byy = fields(i,j,k,Idx::Byy); - Bzz = fields(i,j,k,Idx::Bzz); + Bxx = fields(i,j,k,Idx.Bxx); + Byy = fields(i,j,k,Idx.Byy); + Bzz = fields(i,j,k,Idx.Bzz); - Fx = fields(i,j,k,Idx::Fx); - Fy = fields(i,j,k,Idx::Fy); - Fz = fields(i,j,k,Idx::Fz); + Fx = fields(i,j,k,Idx.Fx); + Fy = fields(i,j,k,Idx.Fy); + Fz = fields(i,j,k,Idx.Fz); - Gx = fields(i,j,k,Idx::Gx); - Gy = fields(i,j,k,Idx::Gy); - Gz = fields(i,j,k,Idx::Gz); + Gx = fields(i,j,k,Idx.Gx); + Gy = fields(i,j,k,Idx.Gy); + Gz = fields(i,j,k,Idx.Gz); Ex = Exx + Exy + Exz; Ey = Eyx + Eyy + Eyz; @@ -216,42 +218,42 @@ PMLPsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const { const Complex C22_c2 = C22 / c2; // Update E - fields(i,j,k,Idx::Exy) = C2 * Exy + C5 * Exz + C9 * Ey - + C10 * Bx + C11 * By + C19 * Bz; + fields(i,j,k,Idx.Exy) = C2 * Exy + C5 * Exz + C9 * Ey + + C10 * Bx + C11 * By + C19 * Bz; - fields(i,j,k,Idx::Exz) = C6 * Exy + C3 * Exz + C8 * Ez - - C10 * Bx - C22 * By - C12 * Bz; + fields(i,j,k,Idx.Exz) = C6 * Exy + C3 * Exz + C8 * Ez + - C10 * Bx - C22 * By - C12 * Bz; - fields(i,j,k,Idx::Eyz) = C3 * Eyz + C6 * Eyx + C7 * Ez - + C21 * Bx + C10 * By + C13 * Bz; + fields(i,j,k,Idx.Eyz) = C3 * Eyz + C6 * Eyx + C7 * Ez + + C21 * Bx + C10 * By + C13 * Bz; - fields(i,j,k,Idx::Eyx) = C9 * Ex + C4 * Eyz + C1 * Eyx - - C14 * Bx - C10 * By - C18 * Bz; + fields(i,j,k,Idx.Eyx) = C9 * Ex + C4 * Eyz + C1 * Eyx + - C14 * Bx - C10 * By - C18 * Bz; - fields(i,j,k,Idx::Ezx) = C8 * Ex + C1 * Ezx + C4 * Ezy - + C15 * Bx + C17 * By + C10 * Bz; + fields(i,j,k,Idx.Ezx) = C8 * Ex + C1 * Ezx + C4 * Ezy + + C15 * Bx + C17 * By + C10 * Bz; - fields(i,j,k,Idx::Ezy) = C7 * Ey + C5 * Ezx + C2 * Ezy - - C20 * Bx - C16 * By - C10 * Bz; + fields(i,j,k,Idx.Ezy) = C7 * Ey + C5 * Ezx + C2 * Ezy + - C20 * Bx - C16 * By - C10 * Bz; // Update B - fields(i,j,k,Idx::Bxy) = C2 * Bxy + C5 * Bxz + C9 * By - - C10_c2 * Ex - C11_c2 * Ey - C19_c2 * Ez; + fields(i,j,k,Idx.Bxy) = C2 * Bxy + C5 * Bxz + C9 * By + - C10_c2 * Ex - C11_c2 * Ey - C19_c2 * Ez; - fields(i,j,k,Idx::Bxz) = C6 * Bxy + C3 * Bxz + C8 * Bz - + C10_c2 * Ex + C22_c2 * Ey + C12_c2 * Ez; + fields(i,j,k,Idx.Bxz) = C6 * Bxy + C3 * Bxz + C8 * Bz + + C10_c2 * Ex + C22_c2 * Ey + C12_c2 * Ez; - fields(i,j,k,Idx::Byz) = C3 * Byz + C6 * Byx + C7 * Bz - - C21_c2 * Ex - C10_c2 * Ey - C13_c2 * Ez; + fields(i,j,k,Idx.Byz) = C3 * Byz + C6 * Byx + C7 * Bz + - C21_c2 * Ex - C10_c2 * Ey - C13_c2 * Ez; - fields(i,j,k,Idx::Byx) = C9 * Bx + C4 * Byz + C1 * Byx - + C14_c2 * Ex + C10_c2 * Ey + C18_c2 * Ez; + fields(i,j,k,Idx.Byx) = C9 * Bx + C4 * Byz + C1 * Byx + + C14_c2 * Ex + C10_c2 * Ey + C18_c2 * Ez; - fields(i,j,k,Idx::Bzx) = C8 * Bx + C1 * Bzx + C4 * Bzy - - C15_c2 * Ex - C17_c2 * Ey - C10_c2 * Ez; + fields(i,j,k,Idx.Bzx) = C8 * Bx + C1 * Bzx + C4 * Bzy + - C15_c2 * Ex - C17_c2 * Ey - C10_c2 * Ez; - fields(i,j,k,Idx::Bzy) = C7 * By + C5 * Bzx + C2 * Bzy - + C20_c2 * Ex + C16_c2 * Ey + C10_c2 * Ez; + fields(i,j,k,Idx.Bzy) = C7 * By + C5 * Bzx + C2 * Bzy + + C20_c2 * Ex + C16_c2 * Ey + C10_c2 * Ez; } else if (dive_cleaning && divb_cleaning) { @@ -264,80 +266,80 @@ PMLPsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const { const Complex C25_c2 = C25 / c2; // Update E - fields(i,j,k,Idx::Exx) = C1 * Exx + C4 * Exy + C4 * Exz - - C9 * Ey - C8 * Ez + C23 * F; + fields(i,j,k,Idx.Exx) = C1 * Exx + C4 * Exy + C4 * Exz + - C9 * Ey - C8 * Ez + C23 * F; - fields(i,j,k,Idx::Exy) = C5 * Exx + C2 * Exy + C5 * Exz - + C9 * Ey + C24 * Bz - C7 * G; + fields(i,j,k,Idx.Exy) = C5 * Exx + C2 * Exy + C5 * Exz + + C9 * Ey + C24 * Bz - C7 * G; - fields(i,j,k,Idx::Exz) = C6 * Exx + C6 * Exy + C3 * Exz - + C8 * Ez - C25 * By + C7 * G; + fields(i,j,k,Idx.Exz) = C6 * Exx + C6 * Exy + C3 * Exz + + C8 * Ez - C25 * By + C7 * G; - fields(i,j,k,Idx::Eyx) = C9 * Ex + C1 * Eyx + C4 * Eyy - + C4 * Eyz - C23 * Bz + C8 * G; + fields(i,j,k,Idx.Eyx) = C9 * Ex + C1 * Eyx + C4 * Eyy + + C4 * Eyz - C23 * Bz + C8 * G; - fields(i,j,k,Idx::Eyy) = - C9 * Ex + C5 * Eyx + C2 * Eyy - + C5 * Eyz - C7 * Ez + C24 * F; + fields(i,j,k,Idx.Eyy) = - C9 * Ex + C5 * Eyx + C2 * Eyy + + C5 * Eyz - C7 * Ez + C24 * F; - fields(i,j,k,Idx::Eyz) = C6 * Eyx + C6 * Eyy + C3 * Eyz - + C7 * Ez + C25 * Bx - C8 * G; + fields(i,j,k,Idx.Eyz) = C6 * Eyx + C6 * Eyy + C3 * Eyz + + C7 * Ez + C25 * Bx - C8 * G; - fields(i,j,k,Idx::Ezx) = C8 * Ex + C1 * Ezx + C4 * Ezy - + C4 * Ezz + C23 * By - C9 * G; + fields(i,j,k,Idx.Ezx) = C8 * Ex + C1 * Ezx + C4 * Ezy + + C4 * Ezz + C23 * By - C9 * G; - fields(i,j,k,Idx::Ezy) = C7 * Ey + C5 * Ezx + C2 * Ezy - + C5 * Ezz - C24 * Bx + C9 * G; + fields(i,j,k,Idx.Ezy) = C7 * Ey + C5 * Ezx + C2 * Ezy + + C5 * Ezz - C24 * Bx + C9 * G; - fields(i,j,k,Idx::Ezz) = - C8 * Ex - C7 * Ey + C6 * Ezx - + C6 * Ezy + C3 * Ezz + C25 * F; + fields(i,j,k,Idx.Ezz) = - C8 * Ex - C7 * Ey + C6 * Ezx + + C6 * Ezy + C3 * Ezz + C25 * F; // Update B - fields(i,j,k,Idx::Bxx) = C1 * Bxx + C4 * Bxy + C4 * Bxz - - C9 * By - C8 * Bz + C23_c2 * G; + fields(i,j,k,Idx.Bxx) = C1 * Bxx + C4 * Bxy + C4 * Bxz + - C9 * By - C8 * Bz + C23_c2 * G; - fields(i,j,k,Idx::Bxy) = - C24_c2 * Ez + C5 * Bxx + C2 * Bxy - + C5 * Bxz + C9 * By + C7 * F; + fields(i,j,k,Idx.Bxy) = - C24_c2 * Ez + C5 * Bxx + C2 * Bxy + + C5 * Bxz + C9 * By + C7 * F; - fields(i,j,k,Idx::Bxz) = C25_c2 * Ey + C6 * Bxx + C6 * Bxy - + C3 * Bxz + C8 * Bz - C7 * F; + fields(i,j,k,Idx.Bxz) = C25_c2 * Ey + C6 * Bxx + C6 * Bxy + + C3 * Bxz + C8 * Bz - C7 * F; - fields(i,j,k,Idx::Byx) = C23_c2 * Ez + C9 * Bx + C1 * Byx - + C4 * Byy + C4 * Byz - C8 * F; + fields(i,j,k,Idx.Byx) = C23_c2 * Ez + C9 * Bx + C1 * Byx + + C4 * Byy + C4 * Byz - C8 * F; - fields(i,j,k,Idx::Byy) = - C9 * Bx + C5 * Byx + C2 * Byy - + C5 * Byz - C7 * Bz + C24_c2 * G; + fields(i,j,k,Idx.Byy) = - C9 * Bx + C5 * Byx + C2 * Byy + + C5 * Byz - C7 * Bz + C24_c2 * G; - fields(i,j,k,Idx::Byz) = - C25_c2 * Ex + C6 * Byx + C6 * Byy - + C3 * Byz + C7 * Bz + C8 * F; + fields(i,j,k,Idx.Byz) = - C25_c2 * Ex + C6 * Byx + C6 * Byy + + C3 * Byz + C7 * Bz + C8 * F; - fields(i,j,k,Idx::Bzx) = - C23_c2 * Ey + C8 * Bx + C1 * Bzx - + C4 * Bzy + C4 * Bzz + C9 * F; + fields(i,j,k,Idx.Bzx) = - C23_c2 * Ey + C8 * Bx + C1 * Bzx + + C4 * Bzy + C4 * Bzz + C9 * F; - fields(i,j,k,Idx::Bzy) = C24_c2 * Ex + C7 * By + C5 * Bzx - + C2 * Bzy + C5 * Bzz - C9 * F; + fields(i,j,k,Idx.Bzy) = C24_c2 * Ex + C7 * By + C5 * Bzx + + C2 * Bzy + C5 * Bzz - C9 * F; - fields(i,j,k,Idx::Bzz) = - C8 * Bx - C7 * By + C6 * Bzx - + C6 * Bzy + C3 * Bzz + C25_c2 * G; + fields(i,j,k,Idx.Bzz) = - C8 * Bx - C7 * By + C6 * Bzx + + C6 * Bzy + C3 * Bzz + C25_c2 * G; // Update F - fields(i,j,k,Idx::Fx) = C23_c2 * Ex + C8 * By - C9 * Bz - + C1 * Fx + C4 * Fy + C4 * Fz; + fields(i,j,k,Idx.Fx) = C23_c2 * Ex + C8 * By - C9 * Bz + + C1 * Fx + C4 * Fy + C4 * Fz; - fields(i,j,k,Idx::Fy) = C24_c2 * Ey - C7 * Bx + C9 * Bz - + C5 * Fx + C2 * Fy + C5 * Fz; + fields(i,j,k,Idx.Fy) = C24_c2 * Ey - C7 * Bx + C9 * Bz + + C5 * Fx + C2 * Fy + C5 * Fz; - fields(i,j,k,Idx::Fz) = C25_c2 * Ez + C7 * Bx - C8 * By - + C6 * Fx + C6 * Fy + C3 * Fz; + fields(i,j,k,Idx.Fz) = C25_c2 * Ez + C7 * Bx - C8 * By + + C6 * Fx + C6 * Fy + C3 * Fz; // Update G - fields(i,j,k,Idx::Gx) = - C8 * Ey + C9 * Ez + C23 * Bx - + C1 * Gx + C4 * Gy + C4 * Gz; + fields(i,j,k,Idx.Gx) = - C8 * Ey + C9 * Ez + C23 * Bx + + C1 * Gx + C4 * Gy + C4 * Gz; - fields(i,j,k,Idx::Gy) = C7 * Ex - C9 * Ez + C24 * By - + C5 * Gx + C2 * Gy + C5 * Gz; + fields(i,j,k,Idx.Gy) = C7 * Ex - C9 * Ez + C24 * By + + C5 * Gx + C2 * Gy + C5 * Gz; - fields(i,j,k,Idx::Gz) = - C7 * Ex + C8 * Ey + C25 * Bz - + C6 * Gx + C6 * Gy + C3 * Gz; + fields(i,j,k,Idx.Gz) = - C7 * Ex + C8 * Ey + C25 * Bz + + C6 * Gx + C6 * Gy + C3 * Gz; } } }); @@ -416,8 +418,4 @@ PMLPsatdAlgorithm::VayDeposition (const int /*lev*/, amrex::Abort("Vay deposition not implemented for PML PSATD"); } -int PMLPsatdAlgorithm::getRequiredNumberOfFields() const { - return SpectralPMLIndex::n_fields; -} - #endif // WARPX_USE_PSATD |