diff options
Diffstat (limited to 'Source/FieldSolver')
24 files changed, 630 insertions, 490 deletions
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H index 577ded61f..40dfeef38 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H @@ -28,6 +28,7 @@ class ComovingPsatdAlgorithm : public SpectralBaseAlgorithm */ ComovingPsatdAlgorithm (const SpectralKSpace& spectral_kspace, const amrex::DistributionMapping& dm, + const SpectralFieldIndex& spectral_index, const int norder_x, const int norder_y, const int norder_z, @@ -42,11 +43,6 @@ class ComovingPsatdAlgorithm : public SpectralBaseAlgorithm */ virtual void pushSpectralFields (SpectralFieldData& f) const override final; - virtual int getRequiredNumberOfFields () const override final - { - return SpectralFieldIndex::n_fields; - } - /* \brief Initialize the coefficients needed in the update equations */ void InitializeSpectralCoefficients (const SpectralKSpace& spectral_kspace, @@ -89,6 +85,8 @@ class ComovingPsatdAlgorithm : public SpectralBaseAlgorithm SpectralRealCoefficients C_coef, S_ck_coef; SpectralComplexCoefficients Theta2_coef, X1_coef, X2_coef, X3_coef, X4_coef; + SpectralFieldIndex m_spectral_index; + // k vectors KVectorComponent kx_vec; #if (AMREX_SPACEDIM==3) diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp index d78ece8f5..f80ee7749 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp @@ -22,6 +22,7 @@ using namespace amrex; ComovingPsatdAlgorithm::ComovingPsatdAlgorithm (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, @@ -29,7 +30,8 @@ ComovingPsatdAlgorithm::ComovingPsatdAlgorithm (const SpectralKSpace& spectral_k const amrex::Real dt, const bool update_with_rho) // Members initialization - : 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), // Initialize the infinite-order k vectors (the argument n_order = -1 selects // the infinite order option, the argument nodal = false is then irrelevant) kx_vec(spectral_kspace.getModifiedKComponent(dm, 0, -1, false)), @@ -64,6 +66,8 @@ ComovingPsatdAlgorithm::ComovingPsatdAlgorithm (const SpectralKSpace& spectral_k void ComovingPsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const { + const SpectralFieldIndex& Idx = m_spectral_index; + // Loop over boxes for (amrex::MFIter mfi(f.fields); mfi.isValid(); ++mfi){ @@ -91,20 +95,19 @@ ComovingPsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const amrex::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); + 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); // Shortcuts 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); + 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 const amrex::Real kx_mod = modified_kx_arr[i]; @@ -130,23 +133,23 @@ ComovingPsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const const Complex X4 = X4_arr(i,j,k); // Update E - fields(i,j,k,Idx::Ex) = C*Ex_old + S_ck*c2*I*(ky_mod*Bz_old - kz_mod*By_old) + fields(i,j,k,Idx.Ex) = C*Ex_old + S_ck*c2*I*(ky_mod*Bz_old - kz_mod*By_old) + X4*Jx - I*(X2*rho_new - X3*rho_old)*kx_mod; - fields(i,j,k,Idx::Ey) = C*Ey_old + S_ck*c2*I*(kz_mod*Bx_old - kx_mod*Bz_old) + fields(i,j,k,Idx.Ey) = C*Ey_old + S_ck*c2*I*(kz_mod*Bx_old - kx_mod*Bz_old) + X4*Jy - I*(X2*rho_new - X3*rho_old)*ky_mod; - fields(i,j,k,Idx::Ez) = C*Ez_old + S_ck*c2*I*(kx_mod*By_old - ky_mod*Bx_old) + fields(i,j,k,Idx.Ez) = C*Ez_old + S_ck*c2*I*(kx_mod*By_old - ky_mod*Bx_old) + X4*Jz - I*(X2*rho_new - X3*rho_old)*kz_mod; // Update B - fields(i,j,k,Idx::Bx) = C*Bx_old - S_ck*I*(ky_mod*Ez_old - kz_mod*Ey_old) + fields(i,j,k,Idx.Bx) = C*Bx_old - S_ck*I*(ky_mod*Ez_old - kz_mod*Ey_old) + X1*I*(ky_mod*Jz - kz_mod*Jy); - fields(i,j,k,Idx::By) = C*By_old - S_ck*I*(kz_mod*Ex_old - kx_mod*Ez_old) + fields(i,j,k,Idx.By) = C*By_old - S_ck*I*(kz_mod*Ex_old - kx_mod*Ez_old) + X1*I*(kz_mod*Jx - kx_mod*Jz); - fields(i,j,k,Idx::Bz) = C*Bz_old - S_ck*I*(kx_mod*Ey_old - ky_mod*Ex_old) + fields(i,j,k,Idx.Bz) = C*Bz_old - S_ck*I*(kx_mod*Ey_old - ky_mod*Ex_old) + X1*I*(kx_mod*Jy - ky_mod*Jx); }); } @@ -416,14 +419,14 @@ ComovingPsatdAlgorithm::CurrentCorrection (const int lev, // Profiling BL_PROFILE("ComovingAlgorithm::CurrentCorrection"); - using Idx = SpectralFieldIndex; + const SpectralFieldIndex& Idx = m_spectral_index; // Forward Fourier transform of J and rho - field_data.ForwardTransform(lev, *current[0], Idx::Jx, 0); - field_data.ForwardTransform(lev, *current[1], Idx::Jy, 0); - field_data.ForwardTransform(lev, *current[2], Idx::Jz, 0); - field_data.ForwardTransform(lev, *rho, Idx::rho_old, 0); - field_data.ForwardTransform(lev, *rho, Idx::rho_new, 1); + field_data.ForwardTransform(lev, *current[0], Idx.Jx, 0); + field_data.ForwardTransform(lev, *current[1], Idx.Jy, 0); + field_data.ForwardTransform(lev, *current[2], Idx.Jz, 0); + field_data.ForwardTransform(lev, *rho, Idx.rho_old, 0); + field_data.ForwardTransform(lev, *rho, Idx.rho_new, 1); const amrex::IntVect& fill_guards = m_fill_guards; @@ -457,11 +460,11 @@ ComovingPsatdAlgorithm::CurrentCorrection (const int lev, amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Shortcuts 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); + 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 amrex::Real kx_mod = modified_kx_arr[i]; @@ -492,24 +495,24 @@ ComovingPsatdAlgorithm::CurrentCorrection (const int lev, const Complex theta = amrex::exp(- I * k_dot_v * dt * 0.5_rt); const Complex den = 1._rt - theta * theta; - fields(i,j,k,Idx::Jx) = Jx - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * kx_mod / (knorm_mod * knorm_mod); - fields(i,j,k,Idx::Jy) = Jy - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * ky_mod / (knorm_mod * knorm_mod); - fields(i,j,k,Idx::Jz) = Jz - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * kz_mod / (knorm_mod * knorm_mod); + fields(i,j,k,Idx.Jx) = Jx - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * kx_mod / (knorm_mod * knorm_mod); + fields(i,j,k,Idx.Jy) = Jy - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * ky_mod / (knorm_mod * knorm_mod); + fields(i,j,k,Idx.Jz) = Jz - (kmod_dot_J + k_dot_v * theta * (rho_new - rho_old) / den) * kz_mod / (knorm_mod * knorm_mod); } else { - fields(i,j,k,Idx::Jx) = Jx - (kmod_dot_J - I * (rho_new - rho_old) / dt) * kx_mod / (knorm_mod * knorm_mod); - fields(i,j,k,Idx::Jy) = Jy - (kmod_dot_J - I * (rho_new - rho_old) / dt) * ky_mod / (knorm_mod * knorm_mod); - fields(i,j,k,Idx::Jz) = Jz - (kmod_dot_J - I * (rho_new - rho_old) / dt) * kz_mod / (knorm_mod * knorm_mod); + fields(i,j,k,Idx.Jx) = Jx - (kmod_dot_J - I * (rho_new - rho_old) / dt) * kx_mod / (knorm_mod * knorm_mod); + fields(i,j,k,Idx.Jy) = Jy - (kmod_dot_J - I * (rho_new - rho_old) / dt) * ky_mod / (knorm_mod * knorm_mod); + fields(i,j,k,Idx.Jz) = Jz - (kmod_dot_J - I * (rho_new - rho_old) / dt) * kz_mod / (knorm_mod * knorm_mod); } } }); } // Backward Fourier transform of J - field_data.BackwardTransform(lev, *current[0], Idx::Jx, 0, fill_guards); - field_data.BackwardTransform(lev, *current[1], Idx::Jy, 0, fill_guards); - field_data.BackwardTransform(lev, *current[2], Idx::Jz, 0, fill_guards); + field_data.BackwardTransform(lev, *current[0], Idx.Jx, 0, fill_guards); + field_data.BackwardTransform(lev, *current[1], Idx.Jy, 0, fill_guards); + field_data.BackwardTransform(lev, *current[2], Idx.Jz, 0, fill_guards); } void diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.H index c5680a4ce..8ac6a1d58 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.H @@ -18,6 +18,7 @@ class GalileanPsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ public: GalileanPsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace, amrex::DistributionMapping const & dm, + const SpectralFieldIndex& spectral_index, int const n_rz_azimuthal_modes, int const norder_z, bool const nodal, const amrex::Array<amrex::Real,3>& v_galilean, @@ -25,9 +26,6 @@ class GalileanPsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ bool const update_with_rho); // Redefine functions from base class virtual void pushSpectralFields (SpectralFieldDataRZ & f) override final; - virtual int getRequiredNumberOfFields () const override final { - return SpectralFieldIndex::n_fields; - } void InitializeSpectralCoefficients (SpectralFieldDataRZ const & f); @@ -61,6 +59,9 @@ class GalileanPsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final; private: + + SpectralFieldIndex m_spectral_index; + bool coefficients_initialized; // Note that dt and v_galilean are saved to use in InitializeSpectralCoefficients amrex::Real const m_dt; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.cpp index ce6d5b933..79fee3c2e 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.cpp @@ -17,13 +17,15 @@ using namespace amrex::literals; /* \brief Initialize coefficients for the update equation */ GalileanPsatdAlgorithmRZ::GalileanPsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace, amrex::DistributionMapping const & dm, + const SpectralFieldIndex& spectral_index, int const n_rz_azimuthal_modes, int const norder_z, bool const nodal, const amrex::Array<amrex::Real,3>& v_galilean, amrex::Real const dt, bool const update_with_rho) // Initialize members of base class - : SpectralBaseAlgorithmRZ(spectral_kspace, dm, norder_z, nodal), + : SpectralBaseAlgorithmRZ(spectral_kspace, dm, spectral_index, norder_z, nodal), + m_spectral_index(spectral_index), m_dt(dt), m_v_galilean(v_galilean), m_update_with_rho(update_with_rho) @@ -60,6 +62,8 @@ GalileanPsatdAlgorithmRZ::pushSpectralFields (SpectralFieldDataRZ & f) coefficients_initialized = true; } + const SpectralFieldIndex& Idx = m_spectral_index; + // Loop over boxes for (amrex::MFIter mfi(f.fields); mfi.isValid(); ++mfi){ @@ -91,18 +95,17 @@ GalileanPsatdAlgorithmRZ::pushSpectralFields (SpectralFieldDataRZ & f) { // All of the fields of each mode are grouped together - using Idx = SpectralFieldIndex; - auto const Ep_m = Idx::Ex + Idx::n_fields*mode; - auto const Em_m = Idx::Ey + Idx::n_fields*mode; - auto const Ez_m = Idx::Ez + Idx::n_fields*mode; - auto const Bp_m = Idx::Bx + Idx::n_fields*mode; - auto const Bm_m = Idx::By + Idx::n_fields*mode; - auto const Bz_m = Idx::Bz + Idx::n_fields*mode; - auto const Jp_m = Idx::Jx + Idx::n_fields*mode; - auto const Jm_m = Idx::Jy + Idx::n_fields*mode; - auto const Jz_m = Idx::Jz + Idx::n_fields*mode; - auto const rho_old_m = Idx::rho_old + Idx::n_fields*mode; - auto const rho_new_m = Idx::rho_new + Idx::n_fields*mode; + auto const Ep_m = Idx.Ex + Idx.n_fields*mode; + auto const Em_m = Idx.Ey + Idx.n_fields*mode; + auto const Ez_m = Idx.Ez + Idx.n_fields*mode; + auto const Bp_m = Idx.Bx + Idx.n_fields*mode; + auto const Bm_m = Idx.By + Idx.n_fields*mode; + auto const Bz_m = Idx.Bz + Idx.n_fields*mode; + auto const Jp_m = Idx.Jx + Idx.n_fields*mode; + auto const Jm_m = Idx.Jy + Idx.n_fields*mode; + auto const Jz_m = Idx.Jz + Idx.n_fields*mode; + auto const rho_old_m = Idx.rho_old + Idx.n_fields*mode; + auto const rho_new_m = Idx.rho_new + Idx.n_fields*mode; // Record old values of the fields to be updated Complex const Ep_old = fields(i,j,k,Ep_m); @@ -295,14 +298,14 @@ GalileanPsatdAlgorithmRZ::CurrentCorrection (const int lev, // Profiling WARPX_PROFILE( "GalileanPsatdAlgorithmRZ::CurrentCorrection" ); - using Idx = SpectralFieldIndex; + const SpectralFieldIndex& Idx = m_spectral_index; // Forward Fourier transform of J and rho - field_data.ForwardTransform( lev, *current[0], Idx::Jx, - *current[1], Idx::Jy); - field_data.ForwardTransform( lev, *current[2], Idx::Jz, 0); - field_data.ForwardTransform( lev, *rho, Idx::rho_old, 0 ); - field_data.ForwardTransform( lev, *rho, Idx::rho_new, 1 ); + field_data.ForwardTransform( lev, *current[0], Idx.Jx, + *current[1], Idx.Jy); + field_data.ForwardTransform( lev, *current[2], Idx.Jz, 0); + field_data.ForwardTransform( lev, *rho, Idx.rho_old, 0 ); + field_data.ForwardTransform( lev, *rho, Idx.rho_new, 1 ); // Loop over boxes for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){ @@ -328,11 +331,11 @@ GalileanPsatdAlgorithmRZ::CurrentCorrection (const int lev, [=] AMREX_GPU_DEVICE(int i, int j, int k, int mode) noexcept { // All of the fields of each mode are grouped together - auto const Jp_m = Idx::Jx + Idx::n_fields*mode; - auto const Jm_m = Idx::Jy + Idx::n_fields*mode; - auto const Jz_m = Idx::Jz + Idx::n_fields*mode; - auto const rho_old_m = Idx::rho_old + Idx::n_fields*mode; - auto const rho_new_m = Idx::rho_new + Idx::n_fields*mode; + auto const Jp_m = Idx.Jx + Idx.n_fields*mode; + auto const Jm_m = Idx.Jy + Idx.n_fields*mode; + auto const Jz_m = Idx.Jz + Idx.n_fields*mode; + auto const rho_old_m = Idx.rho_old + Idx.n_fields*mode; + auto const rho_new_m = Idx.rho_new + Idx.n_fields*mode; // Shortcuts for the values of J and rho Complex const Jp = fields(i,j,k,Jp_m); @@ -367,9 +370,9 @@ GalileanPsatdAlgorithmRZ::CurrentCorrection (const int lev, // Backward Fourier transform of J field_data.BackwardTransform( lev, - *current[0], Idx::Jx, - *current[1], Idx::Jy); - field_data.BackwardTransform( lev, *current[2], Idx::Jz, 0 ); + *current[0], Idx.Jx, + *current[1], Idx.Jy); + field_data.BackwardTransform( lev, *current[2], Idx.Jz, 0 ); } diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H index 83c93b25d..ce71d5a37 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H @@ -29,6 +29,7 @@ class PMLPsatdAlgorithm : public SpectralBaseAlgorithm public: PMLPsatdAlgorithm(const SpectralKSpace& spectral_kspace, const amrex::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, @@ -43,7 +44,6 @@ class PMLPsatdAlgorithm : public SpectralBaseAlgorithm // Redefine functions from base class virtual void pushSpectralFields(SpectralFieldData& f) const override final; - virtual int getRequiredNumberOfFields() const override final; /** * \brief Virtual function for current correction in Fourier space @@ -78,6 +78,7 @@ class PMLPsatdAlgorithm : public SpectralBaseAlgorithm std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final; private: + SpectralFieldIndex m_spectral_index; SpectralRealCoefficients C_coef, S_ck_coef, inv_k2_coef; amrex::Real m_dt; bool m_dive_cleaning; 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 diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H index a53319327..f7cd06edf 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H @@ -33,6 +33,7 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm * * \param[in] spectral_kspace spectral space * \param[in] dm distribution mapping + * \param[in] spectral_index object containing indices to access data in spectral space * \param[in] norder_x order of the spectral solver along x * \param[in] norder_y order of the spectral solver along y * \param[in] norder_z order of the spectral solver along z @@ -47,6 +48,7 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm PsatdAlgorithm ( const SpectralKSpace& spectral_kspace, const amrex::DistributionMapping& dm, + const SpectralFieldIndex& spectral_index, const int norder_x, const int norder_y, const int norder_z, @@ -66,28 +68,6 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm virtual void pushSpectralFields (SpectralFieldData& f) const override final; /** - * \brief Returns the number of fields stored in spectral space - */ - virtual int getRequiredNumberOfFields () const override final - { - if (m_J_linear_in_time) - { - return SpectralFieldIndexJLinearInTime::n_fields; - } - else - { - if (m_time_averaging) - { - return SpectralFieldIndexTimeAveraging::n_fields; - } - else - { - return SpectralFieldIndex::n_fields; - } - } - } - - /** * \brief Initializes the coefficients used in \c pushSpectralFields to update the E and B fields * * \param[in] spectral_kspace spectral space @@ -170,6 +150,8 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm // These real and complex coefficients are allocated only with averaged Galilean PSATD SpectralComplexCoefficients Psi1_coef, Psi2_coef, Y1_coef, Y2_coef, Y3_coef, Y4_coef; + SpectralFieldIndex m_spectral_index; + // Centered modified finite-order k vectors KVectorComponent modified_kx_vec_centered; #if (AMREX_SPACEDIM==3) diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp index bc74a36a2..36e113555 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp @@ -29,6 +29,7 @@ using namespace amrex; PsatdAlgorithm::PsatdAlgorithm( const SpectralKSpace& spectral_kspace, const DistributionMapping& dm, + const SpectralFieldIndex& spectral_index, const int norder_x, const int norder_y, const int norder_z, @@ -40,7 +41,8 @@ PsatdAlgorithm::PsatdAlgorithm( const bool time_averaging, const bool J_linear_in_time) // Initializer list - : 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), // Initialize the centered finite-order modified k vectors: // these are computed always with the assumption of centered grids // (argument nodal = true), for both nodal and staggered simulations @@ -108,6 +110,8 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const const amrex::Real dt = m_dt; + const SpectralFieldIndex& Idx = m_spectral_index; + // Loop over boxes for (amrex::MFIter mfi(f.fields); mfi.isValid(); ++mfi) { @@ -167,24 +171,20 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const // Loop over indices within one box ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - using Idx = SpectralFieldIndex; - using IdxAvg = SpectralFieldIndexTimeAveraging; - using IdxLin = SpectralFieldIndexJLinearInTime; - // Record old values of the fields to be updated - 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); + 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); // Shortcuts 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); + 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 const amrex::Real kx = modified_kx_arr[i]; @@ -215,17 +215,17 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const if (update_with_rho) { - fields(i,j,k,Idx::Ex) = T2 * C * Ex_old - + I * c2 * T2 * S_ck * (ky * Bz_old - kz * By_old) - + X4 * Jx - I * (X2 * rho_new - T2 * X3 * rho_old) * kx; + fields(i,j,k,Idx.Ex) = T2 * C * Ex_old + + I * c2 * T2 * S_ck * (ky * Bz_old - kz * By_old) + + X4 * Jx - I * (X2 * rho_new - T2 * X3 * rho_old) * kx; - fields(i,j,k,Idx::Ey) = T2 * C * Ey_old - + I * c2 * T2 * S_ck * (kz * Bx_old - kx * Bz_old) - + X4 * Jy - I * (X2 * rho_new - T2 * X3 * rho_old) * ky; + fields(i,j,k,Idx.Ey) = T2 * C * Ey_old + + I * c2 * T2 * S_ck * (kz * Bx_old - kx * Bz_old) + + X4 * Jy - I * (X2 * rho_new - T2 * X3 * rho_old) * ky; - fields(i,j,k,Idx::Ez) = T2 * C * Ez_old - + I * c2 * T2 * S_ck * (kx * By_old - ky * Bx_old) - + X4 * Jz - I * (X2 * rho_new - T2 * X3 * rho_old) * kz; + fields(i,j,k,Idx.Ez) = T2 * C * Ez_old + + I * c2 * T2 * S_ck * (kx * By_old - ky * Bx_old) + + X4 * Jz - I * (X2 * rho_new - T2 * X3 * rho_old) * kz; } // Update equations for E in the formulation without rho @@ -236,56 +236,56 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const Complex k_dot_J = kx * Jx + ky * Jy + kz * Jz; Complex k_dot_E = kx * Ex_old + ky * Ey_old + kz * Ez_old; - fields(i,j,k,Idx::Ex) = T2 * C * Ex_old - + I * c2 * T2 * S_ck * (ky * Bz_old - kz * By_old) - + X4 * Jx + X2 * k_dot_E * kx + X3 * k_dot_J * kx; + fields(i,j,k,Idx.Ex) = T2 * C * Ex_old + + I * c2 * T2 * S_ck * (ky * Bz_old - kz * By_old) + + X4 * Jx + X2 * k_dot_E * kx + X3 * k_dot_J * kx; - fields(i,j,k,Idx::Ey) = T2 * C * Ey_old - + I * c2 * T2 * S_ck * (kz * Bx_old - kx * Bz_old) - + X4 * Jy + X2 * k_dot_E * ky + X3 * k_dot_J * ky; + fields(i,j,k,Idx.Ey) = T2 * C * Ey_old + + I * c2 * T2 * S_ck * (kz * Bx_old - kx * Bz_old) + + X4 * Jy + X2 * k_dot_E * ky + X3 * k_dot_J * ky; - fields(i,j,k,Idx::Ez) = T2 * C * Ez_old - + I * c2 * T2 * S_ck * (kx * By_old - ky * Bx_old) - + X4 * Jz + X2 * k_dot_E * kz + X3 * k_dot_J * kz; + fields(i,j,k,Idx.Ez) = T2 * C * Ez_old + + I * c2 * T2 * S_ck * (kx * By_old - ky * Bx_old) + + X4 * Jz + X2 * k_dot_E * kz + X3 * k_dot_J * kz; } // Update equations for B // T2 = 1 always with standard PSATD (zero Galilean velocity) - fields(i,j,k,Idx::Bx) = T2 * C * Bx_old - - I * T2 * S_ck * (ky * Ez_old - kz * Ey_old) - + I * X1 * (ky * Jz - kz * Jy); + fields(i,j,k,Idx.Bx) = T2 * C * Bx_old + - I * T2 * S_ck * (ky * Ez_old - kz * Ey_old) + + I * X1 * (ky * Jz - kz * Jy); - fields(i,j,k,Idx::By) = T2 * C * By_old - - I * T2 * S_ck * (kz * Ex_old - kx * Ez_old) - + I * X1 * (kz * Jx - kx * Jz); + fields(i,j,k,Idx.By) = T2 * C * By_old + - I * T2 * S_ck * (kz * Ex_old - kx * Ez_old) + + I * X1 * (kz * Jx - kx * Jz); - fields(i,j,k,Idx::Bz) = T2 * C * Bz_old - - I * T2 * S_ck * (kx * Ey_old - ky * Ex_old) - + I * X1 * (kx * Jy - ky * Jx); + fields(i,j,k,Idx.Bz) = T2 * C * Bz_old + - I * T2 * S_ck * (kx * Ey_old - ky * Ex_old) + + I * X1 * (kx * Jy - ky * Jx); if (J_linear_in_time) { - const Complex Jx_new = fields(i,j,k,IdxLin::Jx_new); - const Complex Jy_new = fields(i,j,k,IdxLin::Jy_new); - const Complex Jz_new = fields(i,j,k,IdxLin::Jz_new); + const Complex Jx_new = fields(i,j,k,Idx.Jx_new); + const Complex Jy_new = fields(i,j,k,Idx.Jy_new); + const Complex Jz_new = fields(i,j,k,Idx.Jz_new); - const Complex F_old = fields(i,j,k,IdxLin::F); - const Complex G_old = fields(i,j,k,IdxLin::G); + const Complex F_old = fields(i,j,k,Idx.F); + const Complex G_old = fields(i,j,k,Idx.G); - fields(i,j,k,Idx::Ex) += -X1 * (Jx_new - Jx) / dt + I * c2 * S_ck * F_old * kx; + fields(i,j,k,Idx.Ex) += -X1 * (Jx_new - Jx) / dt + I * c2 * S_ck * F_old * kx; - fields(i,j,k,Idx::Ey) += -X1 * (Jy_new - Jy) / dt + I * c2 * S_ck * F_old * ky; + fields(i,j,k,Idx.Ey) += -X1 * (Jy_new - Jy) / dt + I * c2 * S_ck * F_old * ky; - fields(i,j,k,Idx::Ez) += -X1 * (Jz_new - Jz) / dt + I * c2 * S_ck * F_old * kz; + fields(i,j,k,Idx.Ez) += -X1 * (Jz_new - Jz) / dt + I * c2 * S_ck * F_old * kz; - fields(i,j,k,Idx::Bx) += I * X2/c2 * (ky * (Jz_new - Jz) - kz * (Jy_new - Jy)); + fields(i,j,k,Idx.Bx) += I * X2/c2 * (ky * (Jz_new - Jz) - kz * (Jy_new - Jy)); + I * c2 * S_ck * G_old * kx; - fields(i,j,k,Idx::By) += I * X2/c2 * (kz * (Jx_new - Jx) - kx * (Jz_new - Jz)); + fields(i,j,k,Idx.By) += I * X2/c2 * (kz * (Jx_new - Jx) - kx * (Jz_new - Jz)); + I * c2 * S_ck * G_old * ky; - fields(i,j,k,Idx::Bz) += I * X2/c2 * (kx * (Jy_new - Jy) - ky * (Jx_new - Jx)); + fields(i,j,k,Idx.Bz) += I * X2/c2 * (kx * (Jy_new - Jy) - ky * (Jx_new - Jx)); + I * c2 * S_ck * G_old * kz; const Complex k_dot_J = kx * Jx + ky * Jy + kz * Jz; @@ -293,10 +293,10 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const const Complex k_dot_E = kx * Ex_old + ky * Ey_old + kz * Ez_old; const Complex k_dot_B = kx * Bx_old + ky * By_old + kz * Bz_old; - fields(i,j,k,IdxLin::F) = C * F_old + S_ck * (I * k_dot_E - rho_old * inv_ep0) + fields(i,j,k,Idx.F) = C * F_old + S_ck * (I * k_dot_E - rho_old * inv_ep0) - X1 * ((rho_new - rho_old) / dt + I * k_dot_J) - I * X2/c2 * k_dot_dJ; - fields(i,j,k,IdxLin::G) = C * G_old + I * S_ck * k_dot_B; + fields(i,j,k,Idx.G) = C * G_old + I * S_ck * k_dot_B; if (time_averaging) { @@ -307,32 +307,32 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const // because it is meant to be used with sub-cycling // maybe this should be made more generic - fields(i,j,k,IdxLin::Ex_avg) += S_ck * Ex_old + fields(i,j,k,Idx.Ex_avg) += S_ck * Ex_old + I * c2 * ep0 * X1 * (ky * Bz_old - kz * By_old) + I * X5 * rho_old * kx + I * X6 * rho_new * kx + X3/c2 * Jx - X2/c2 * Jx_new + I * c2 * ep0 * X1 * F_old * kx; - fields(i,j,k,IdxLin::Ey_avg) += S_ck * Ey_old + fields(i,j,k,Idx.Ey_avg) += S_ck * Ey_old + I * c2 * ep0 * X1 * (kz * Bx_old - kx * Bz_old) + I * X5 * rho_old * ky + I * X6 * rho_new * ky + X3/c2 * Jy - X2/c2 * Jy_new + I * c2 * ep0 * X1 * F_old * ky; - fields(i,j,k,IdxLin::Ez_avg) += S_ck * Ez_old + fields(i,j,k,Idx.Ez_avg) += S_ck * Ez_old + I * c2 * ep0 * X1 * (kx * By_old - ky * Bx_old) + I * X5 * rho_old * kz + I * X6 * rho_new * kz + X3/c2 * Jz - X2/c2 * Jz_new + I * c2 * ep0 * X1 * F_old * kz; - fields(i,j,k,IdxLin::Bx_avg) += S_ck * Bx_old + fields(i,j,k,Idx.Bx_avg) += S_ck * Bx_old - I * ep0 * X1 * (ky * Ez_old - kz * Ey_old) - I * X5/c2 * (ky * Jz - kz * Jy) - I * X6/c2 * (ky * Jz_new - kz * Jy_new); + I * c2 * ep0 * X1 * G_old * kx; - fields(i,j,k,IdxLin::By_avg) += S_ck * By_old + fields(i,j,k,Idx.By_avg) += S_ck * By_old - I * ep0 * X1 * (kz * Ex_old - kx * Ez_old) - I * X5/c2 * (kz * Jx - kx * Jz) - I * X6/c2 * (kz * Jx_new - kx * Jz_new); + I * c2 * ep0 * X1 * G_old * ky; - fields(i,j,k,IdxLin::Bz_avg) += S_ck * Bz_old + fields(i,j,k,Idx.Bz_avg) += S_ck * Bz_old - I * ep0 * X1 * (kx * Ey_old - ky * Ex_old) - I * X5/c2 * (kx * Jy - ky * Jx) - I * X6/c2 * (kx * Jy_new - ky * Jx_new); + I * c2 * ep0 * X1 * G_old * kz; @@ -350,29 +350,29 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const const Complex Y2 = Y2_arr(i,j,k); const Complex Y4 = Y4_arr(i,j,k); - fields(i,j,k,IdxAvg::Ex_avg) = Psi1 * Ex_old - - I * c2 * Psi2 * (ky * Bz_old - kz * By_old) - + Y4 * Jx + (Y2 * rho_new + Y3 * rho_old) * kx; + fields(i,j,k,Idx.Ex_avg) = Psi1 * Ex_old + - I * c2 * Psi2 * (ky * Bz_old - kz * By_old) + + Y4 * Jx + (Y2 * rho_new + Y3 * rho_old) * kx; - fields(i,j,k,IdxAvg::Ey_avg) = Psi1 * Ey_old - - I * c2 * Psi2 * (kz * Bx_old - kx * Bz_old) - + Y4 * Jy + (Y2 * rho_new + Y3 * rho_old) * ky; + fields(i,j,k,Idx.Ey_avg) = Psi1 * Ey_old + - I * c2 * Psi2 * (kz * Bx_old - kx * Bz_old) + + Y4 * Jy + (Y2 * rho_new + Y3 * rho_old) * ky; - fields(i,j,k,IdxAvg::Ez_avg) = Psi1 * Ez_old - - I * c2 * Psi2 * (kx * By_old - ky * Bx_old) - + Y4 * Jz + (Y2 * rho_new + Y3 * rho_old) * kz; + fields(i,j,k,Idx.Ez_avg) = Psi1 * Ez_old + - I * c2 * Psi2 * (kx * By_old - ky * Bx_old) + + Y4 * Jz + (Y2 * rho_new + Y3 * rho_old) * kz; - fields(i,j,k,IdxAvg::Bx_avg) = Psi1 * Bx_old - + I * Psi2 * (ky * Ez_old - kz * Ey_old) - + I * Y1 * (ky * Jz - kz * Jy); + fields(i,j,k,Idx.Bx_avg) = Psi1 * Bx_old + + I * Psi2 * (ky * Ez_old - kz * Ey_old) + + I * Y1 * (ky * Jz - kz * Jy); - fields(i,j,k,IdxAvg::By_avg) = Psi1 * By_old - + I * Psi2 * (kz * Ex_old - kx * Ez_old) - + I * Y1 * (kz * Jx - kx * Jz); + fields(i,j,k,Idx.By_avg) = Psi1 * By_old + + I * Psi2 * (kz * Ex_old - kx * Ez_old) + + I * Y1 * (kz * Jx - kx * Jz); - fields(i,j,k,IdxAvg::Bz_avg) = Psi1 * Bz_old - + I * Psi2 * (kx * Ey_old - ky * Ex_old) - + I * Y1 * (kx * Jy - ky * Jx); + fields(i,j,k,Idx.Bz_avg) = Psi1 * Bz_old + + I * Psi2 * (kx * Ey_old - ky * Ex_old) + + I * Y1 * (kx * Jy - ky * Jx); } }); } @@ -851,14 +851,14 @@ PsatdAlgorithm::CurrentCorrection ( // Profiling BL_PROFILE("PsatdAlgorithm::CurrentCorrection"); - using Idx = SpectralFieldIndex; + const SpectralFieldIndex& Idx = m_spectral_index; // Forward Fourier transform of J and rho - field_data.ForwardTransform(lev, *current[0], Idx::Jx, 0); - field_data.ForwardTransform(lev, *current[1], Idx::Jy, 0); - field_data.ForwardTransform(lev, *current[2], Idx::Jz, 0); - field_data.ForwardTransform(lev, *rho, Idx::rho_old, 0); - field_data.ForwardTransform(lev, *rho, Idx::rho_new, 1); + field_data.ForwardTransform(lev, *current[0], Idx.Jx, 0); + field_data.ForwardTransform(lev, *current[1], Idx.Jy, 0); + field_data.ForwardTransform(lev, *current[2], Idx.Jz, 0); + field_data.ForwardTransform(lev, *rho, Idx.rho_old, 0); + field_data.ForwardTransform(lev, *rho, Idx.rho_new, 1); const amrex::IntVect& fill_guards = m_fill_guards; @@ -892,11 +892,11 @@ PsatdAlgorithm::CurrentCorrection ( ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Shortcuts 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); + 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 amrex::Real kx = modified_kx_arr[i]; @@ -928,25 +928,25 @@ PsatdAlgorithm::CurrentCorrection ( const Complex rho_old_mod = rho_old * amrex::exp(I * k_dot_vg * dt); const Complex den = 1._rt - amrex::exp(I * k_dot_vg * dt); - fields(i,j,k,Idx::Jx) = Jx - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den) + fields(i,j,k,Idx.Jx) = Jx - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den) * kx / (k_norm * k_norm); - fields(i,j,k,Idx::Jy) = Jy - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den) + fields(i,j,k,Idx.Jy) = Jy - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den) * ky / (k_norm * k_norm); - fields(i,j,k,Idx::Jz) = Jz - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den) + fields(i,j,k,Idx.Jz) = Jz - (k_dot_J - k_dot_vg * (rho_new - rho_old_mod) / den) * kz / (k_norm * k_norm); } else { - fields(i,j,k,Idx::Jx) = Jx - (k_dot_J - I * (rho_new - rho_old) / dt) + fields(i,j,k,Idx.Jx) = Jx - (k_dot_J - I * (rho_new - rho_old) / dt) * kx / (k_norm * k_norm); - fields(i,j,k,Idx::Jy) = Jy - (k_dot_J - I * (rho_new - rho_old) / dt) + fields(i,j,k,Idx.Jy) = Jy - (k_dot_J - I * (rho_new - rho_old) / dt) * ky / (k_norm * k_norm); - fields(i,j,k,Idx::Jz) = Jz - (k_dot_J - I * (rho_new - rho_old) / dt) + fields(i,j,k,Idx.Jz) = Jz - (k_dot_J - I * (rho_new - rho_old) / dt) * kz / (k_norm * k_norm); } } @@ -954,9 +954,9 @@ PsatdAlgorithm::CurrentCorrection ( } // Backward Fourier transform of J - field_data.BackwardTransform(lev, *current[0], Idx::Jx, 0, fill_guards); - field_data.BackwardTransform(lev, *current[1], Idx::Jy, 0, fill_guards); - field_data.BackwardTransform(lev, *current[2], Idx::Jz, 0, fill_guards); + field_data.BackwardTransform(lev, *current[0], Idx.Jx, 0, fill_guards); + field_data.BackwardTransform(lev, *current[1], Idx.Jy, 0, fill_guards); + field_data.BackwardTransform(lev, *current[2], Idx.Jz, 0, fill_guards); } void @@ -968,14 +968,14 @@ PsatdAlgorithm::VayDeposition ( // Profiling BL_PROFILE("PsatdAlgorithm::VayDeposition()"); - using Idx = SpectralFieldIndex; + const SpectralFieldIndex& Idx = m_spectral_index; // Forward Fourier transform of D (temporarily stored in current): // D is nodal and does not match the staggering of J, therefore we pass the // actual staggering of D (IntVect(1)) to the ForwardTransform function - field_data.ForwardTransform(lev, *current[0], Idx::Jx, 0, IntVect(1)); - field_data.ForwardTransform(lev, *current[1], Idx::Jy, 0, IntVect(1)); - field_data.ForwardTransform(lev, *current[2], Idx::Jz, 0, IntVect(1)); + field_data.ForwardTransform(lev, *current[0], Idx.Jx, 0, IntVect(1)); + field_data.ForwardTransform(lev, *current[1], Idx.Jy, 0, IntVect(1)); + field_data.ForwardTransform(lev, *current[2], Idx.Jz, 0, IntVect(1)); const amrex::IntVect& fill_guards = m_fill_guards; @@ -998,11 +998,11 @@ PsatdAlgorithm::VayDeposition ( ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Shortcuts for the values of D - const Complex Dx = fields(i,j,k,Idx::Jx); + const Complex Dx = fields(i,j,k,Idx.Jx); #if (AMREX_SPACEDIM == 3) - const Complex Dy = fields(i,j,k,Idx::Jy); + const Complex Dy = fields(i,j,k,Idx.Jy); #endif - const Complex Dz = fields(i,j,k,Idx::Jz); + const Complex Dz = fields(i,j,k,Idx.Jz); // Imaginary unit constexpr Complex I = Complex{0._rt, 1._rt}; @@ -1017,25 +1017,25 @@ PsatdAlgorithm::VayDeposition ( #endif // Compute Jx - if (kx_mod != 0._rt) fields(i,j,k,Idx::Jx) = I * Dx / kx_mod; - else fields(i,j,k,Idx::Jx) = 0._rt; + if (kx_mod != 0._rt) fields(i,j,k,Idx.Jx) = I * Dx / kx_mod; + else fields(i,j,k,Idx.Jx) = 0._rt; #if (AMREX_SPACEDIM == 3) // Compute Jy - if (ky_mod != 0._rt) fields(i,j,k,Idx::Jy) = I * Dy / ky_mod; - else fields(i,j,k,Idx::Jy) = 0._rt; + if (ky_mod != 0._rt) fields(i,j,k,Idx.Jy) = I * Dy / ky_mod; + else fields(i,j,k,Idx.Jy) = 0._rt; #endif // Compute Jz - if (kz_mod != 0._rt) fields(i,j,k,Idx::Jz) = I * Dz / kz_mod; - else fields(i,j,k,Idx::Jz) = 0._rt; + if (kz_mod != 0._rt) fields(i,j,k,Idx.Jz) = I * Dz / kz_mod; + else fields(i,j,k,Idx.Jz) = 0._rt; }); } // Backward Fourier transform of J - field_data.BackwardTransform(lev, *current[0], Idx::Jx, 0, fill_guards); - field_data.BackwardTransform(lev, *current[1], Idx::Jy, 0, fill_guards); - field_data.BackwardTransform(lev, *current[2], Idx::Jz, 0, fill_guards); + field_data.BackwardTransform(lev, *current[0], Idx.Jx, 0, fill_guards); + field_data.BackwardTransform(lev, *current[1], Idx.Jy, 0, fill_guards); + field_data.BackwardTransform(lev, *current[2], Idx.Jz, 0, fill_guards); } #endif // WARPX_USE_PSATD diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H index 9d6ed9cfa..4987e88d7 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H @@ -18,14 +18,12 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ public: PsatdAlgorithmRZ(SpectralKSpaceRZ const & spectral_kspace, amrex::DistributionMapping const & dm, + const SpectralFieldIndex& spectral_index, int const n_rz_azimuthal_modes, int const norder_z, bool const nodal, amrex::Real const dt_step, bool const update_with_rho); // Redefine functions from base class virtual void pushSpectralFields(SpectralFieldDataRZ & f) override final; - virtual int getRequiredNumberOfFields() const override final { - return SpectralFieldIndex::n_fields; - } void InitializeSpectralCoefficients(SpectralFieldDataRZ const & f); @@ -62,6 +60,9 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final; private: + + SpectralFieldIndex m_spectral_index; + bool coefficients_initialized; // Note that dt is saved to use in InitializeSpectralCoefficients amrex::Real m_dt; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp index cb055a27b..6a3c86bf6 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp @@ -17,12 +17,13 @@ using amrex::operator""_rt; /* \brief Initialize coefficients for the update equation */ PsatdAlgorithmRZ::PsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace, amrex::DistributionMapping const & dm, + const SpectralFieldIndex& spectral_index, int const n_rz_azimuthal_modes, int const norder_z, bool const nodal, amrex::Real const dt, bool const update_with_rho) // Initialize members of base class - : SpectralBaseAlgorithmRZ(spectral_kspace, dm, - norder_z, nodal), + : SpectralBaseAlgorithmRZ(spectral_kspace, dm, spectral_index, norder_z, nodal), + m_spectral_index(spectral_index), m_dt(dt), m_update_with_rho(update_with_rho) { @@ -53,6 +54,8 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) coefficients_initialized = true; } + const SpectralFieldIndex& Idx = m_spectral_index; + // Loop over boxes for (amrex::MFIter mfi(f.fields); mfi.isValid(); ++mfi){ @@ -82,18 +85,17 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f) { // All of the fields of each mode are grouped together - using Idx = SpectralFieldIndex; - auto const Ep_m = Idx::Ex + Idx::n_fields*mode; - auto const Em_m = Idx::Ey + Idx::n_fields*mode; - auto const Ez_m = Idx::Ez + Idx::n_fields*mode; - auto const Bp_m = Idx::Bx + Idx::n_fields*mode; - auto const Bm_m = Idx::By + Idx::n_fields*mode; - auto const Bz_m = Idx::Bz + Idx::n_fields*mode; - auto const Jp_m = Idx::Jx + Idx::n_fields*mode; - auto const Jm_m = Idx::Jy + Idx::n_fields*mode; - auto const Jz_m = Idx::Jz + Idx::n_fields*mode; - auto const rho_old_m = Idx::rho_old + Idx::n_fields*mode; - auto const rho_new_m = Idx::rho_new + Idx::n_fields*mode; + auto const Ep_m = Idx.Ex + Idx.n_fields*mode; + auto const Em_m = Idx.Ey + Idx.n_fields*mode; + auto const Ez_m = Idx.Ez + Idx.n_fields*mode; + auto const Bp_m = Idx.Bx + Idx.n_fields*mode; + auto const Bm_m = Idx.By + Idx.n_fields*mode; + auto const Bz_m = Idx.Bz + Idx.n_fields*mode; + auto const Jp_m = Idx.Jx + Idx.n_fields*mode; + auto const Jm_m = Idx.Jy + Idx.n_fields*mode; + auto const Jz_m = Idx.Jz + Idx.n_fields*mode; + auto const rho_old_m = Idx.rho_old + Idx.n_fields*mode; + auto const rho_new_m = Idx.rho_new + Idx.n_fields*mode; // Record old values of the fields to be updated Complex const Ep_old = fields(i,j,k,Ep_m); @@ -223,15 +225,15 @@ PsatdAlgorithmRZ::CurrentCorrection (const int lev, // Profiling WARPX_PROFILE( "PsatdAlgorithmRZ::CurrentCorrection" ); - using Idx = SpectralFieldIndex; + const SpectralFieldIndex& Idx = m_spectral_index; // Forward Fourier transform of J and rho field_data.ForwardTransform( lev, - *current[0], Idx::Jx, - *current[1], Idx::Jy); - field_data.ForwardTransform( lev, *current[2], Idx::Jz, 0); - field_data.ForwardTransform( lev, *rho, Idx::rho_old, 0 ); - field_data.ForwardTransform( lev, *rho, Idx::rho_new, 1 ); + *current[0], Idx.Jx, + *current[1], Idx.Jy); + field_data.ForwardTransform( lev, *current[2], Idx.Jz, 0); + field_data.ForwardTransform( lev, *rho, Idx.rho_old, 0 ); + field_data.ForwardTransform( lev, *rho, Idx.rho_new, 1 ); // Loop over boxes for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){ @@ -256,11 +258,11 @@ PsatdAlgorithmRZ::CurrentCorrection (const int lev, [=] AMREX_GPU_DEVICE(int i, int j, int k, int mode) noexcept { // All of the fields of each mode are grouped together - auto const Jp_m = Idx::Jx + Idx::n_fields*mode; - auto const Jm_m = Idx::Jy + Idx::n_fields*mode; - auto const Jz_m = Idx::Jz + Idx::n_fields*mode; - auto const rho_old_m = Idx::rho_old + Idx::n_fields*mode; - auto const rho_new_m = Idx::rho_new + Idx::n_fields*mode; + auto const Jp_m = Idx.Jx + Idx.n_fields*mode; + auto const Jm_m = Idx.Jy + Idx.n_fields*mode; + auto const Jz_m = Idx.Jz + Idx.n_fields*mode; + auto const rho_old_m = Idx.rho_old + Idx.n_fields*mode; + auto const rho_new_m = Idx.rho_new + Idx.n_fields*mode; // Shortcuts for the values of J and rho Complex const Jp = fields(i,j,k,Jp_m); @@ -292,10 +294,10 @@ PsatdAlgorithmRZ::CurrentCorrection (const int lev, // Backward Fourier transform of J field_data.BackwardTransform( lev, - *current[0], Idx::Jx, - *current[1], Idx::Jy); + *current[0], Idx.Jx, + *current[1], Idx.Jy); field_data.BackwardTransform( lev, - *current[2], Idx::Jz, 0 ); + *current[2], Idx.Jz, 0 ); } void diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H index f412231b7..4a35871cd 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H @@ -11,6 +11,7 @@ #include "Utils/WarpX_Complex.H" #include "FieldSolver/SpectralSolver/SpectralFieldData_fwd.H" +#include "FieldSolver/SpectralSolver/SpectralFieldData.H" #include <AMReX_BaseFab.H> #include <AMReX_Config.H> @@ -36,7 +37,7 @@ class SpectralBaseAlgorithm public: // 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. @@ -92,10 +93,13 @@ class SpectralBaseAlgorithm */ SpectralBaseAlgorithm(const SpectralKSpace& spectral_kspace, const amrex::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); + SpectralFieldIndex m_spectral_index; + // Modified finite-order vectors KVectorComponent modified_kx_vec; #if (AMREX_SPACEDIM==3) diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.cpp index 4445705cb..92e3fcf94 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.cpp @@ -29,11 +29,13 @@ using namespace amrex; */ SpectralBaseAlgorithm::SpectralBaseAlgorithm(const SpectralKSpace& spectral_kspace, const amrex::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): - // Compute and assign the modified k vectors m_fill_guards(fill_guards), + m_spectral_index(spectral_index), + // 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)), @@ -57,12 +59,12 @@ SpectralBaseAlgorithm::ComputeSpectralDivE ( const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield, amrex::MultiFab& divE ) { - using Idx = SpectralFieldIndex; + const SpectralFieldIndex& Idx = m_spectral_index; // Forward Fourier transform of E - field_data.ForwardTransform(lev, *Efield[0], Idx::Ex, 0 ); - field_data.ForwardTransform(lev, *Efield[1], Idx::Ey, 0 ); - field_data.ForwardTransform(lev, *Efield[2], Idx::Ez, 0 ); + field_data.ForwardTransform(lev, *Efield[0], Idx.Ex, 0 ); + field_data.ForwardTransform(lev, *Efield[1], Idx.Ey, 0 ); + field_data.ForwardTransform(lev, *Efield[2], Idx.Ez, 0 ); const amrex::IntVect& fill_guards = m_fill_guards; @@ -85,9 +87,9 @@ SpectralBaseAlgorithm::ComputeSpectralDivE ( [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Shortcuts for the components of E - const Complex Ex = fields(i,j,k,Idx::Ex); - const Complex Ey = fields(i,j,k,Idx::Ey); - const Complex Ez = fields(i,j,k,Idx::Ez); + const Complex Ex = fields(i,j,k,Idx.Ex); + const Complex Ey = fields(i,j,k,Idx.Ey); + const Complex Ez = fields(i,j,k,Idx.Ez); // k vector values const Real kx = modified_kx_arr[i]; #if (AMREX_SPACEDIM==3) @@ -100,10 +102,10 @@ SpectralBaseAlgorithm::ComputeSpectralDivE ( const Complex I = Complex{0,1}; // div(E) in Fourier space - fields(i,j,k,Idx::divE) = I*(kx*Ex+ky*Ey+kz*Ez); + fields(i,j,k,Idx.divE) = I*(kx*Ex+ky*Ey+kz*Ez); }); } // Backward Fourier transform - field_data.BackwardTransform(lev, divE, Idx::divE, 0, fill_guards); + field_data.BackwardTransform(lev, divE, Idx.divE, 0, fill_guards); } diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H index 51bba5b87..a494a3291 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H @@ -22,7 +22,7 @@ class SpectralBaseAlgorithmRZ public: // Virtual member function ; meant to be overridden in subclasses virtual void pushSpectralFields(SpectralFieldDataRZ & f) = 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. @@ -71,11 +71,15 @@ class SpectralBaseAlgorithmRZ // Constructor SpectralBaseAlgorithmRZ(SpectralKSpaceRZ const & spectral_kspace, amrex::DistributionMapping const & dm, + const SpectralFieldIndex& spectral_index, int const norder_z, bool const nodal) // Compute and assign the modified k vectors - : modified_kz_vec(spectral_kspace.getModifiedKComponent(dm, 1, norder_z, nodal)) + : m_spectral_index(spectral_index), + modified_kz_vec(spectral_kspace.getModifiedKComponent(dm, 1, norder_z, nodal)) {} + SpectralFieldIndex m_spectral_index; + // Modified finite-order vectors KVectorComponent modified_kz_vec; }; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.cpp index f80f54208..7e9794d31 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.cpp @@ -21,14 +21,15 @@ SpectralBaseAlgorithmRZ::ComputeSpectralDivE ( amrex::MultiFab& divE ) { using amrex::operator""_rt; - using Idx = SpectralFieldIndex; + + const SpectralFieldIndex& Idx = m_spectral_index; // Forward Fourier transform of E field_data.ForwardTransform( lev, - *Efield[0], Idx::Ex, - *Efield[1], Idx::Ey ); + *Efield[0], Idx.Ex, + *Efield[1], Idx.Ey ); field_data.ForwardTransform( lev, - *Efield[2], Idx::Ez, 0 ); + *Efield[2], Idx.Ez, 0 ); // Loop over boxes for (MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){ @@ -44,15 +45,15 @@ SpectralBaseAlgorithmRZ::ComputeSpectralDivE ( int const nr = bx.length(0); int const modes = field_data.n_rz_azimuthal_modes; - constexpr int n_fields = SpectralFieldIndex::n_fields; + const int n_fields = m_spectral_index.n_fields; // Loop over indices within one box ParallelFor(bx, modes, [=] AMREX_GPU_DEVICE(int i, int j, int /*k*/, int mode) noexcept { - int const ic1 = Idx::Ex + mode*n_fields; - int const ic2 = Idx::Ey + mode*n_fields; - int const ic3 = Idx::Ez + mode*n_fields; + int const ic1 = Idx.Ex + mode*n_fields; + int const ic2 = Idx.Ey + mode*n_fields; + int const ic3 = Idx.Ez + mode*n_fields; // Shortcuts for the components of E Complex const Ep = fields(i,j,0,ic1); @@ -66,11 +67,11 @@ SpectralBaseAlgorithmRZ::ComputeSpectralDivE ( Complex const I = Complex{0._rt,1._rt}; // div(E) in Fourier space - int const ic = Idx::divE + mode*n_fields; + int const ic = Idx.divE + mode*n_fields; fields(i,j,0,ic) = kr*(Ep - Em) + I*kz*Ez; }); } // Backward Fourier transform - field_data.BackwardTransform( lev, divE, Idx::divE, 0 ); + field_data.BackwardTransform( lev, divE, Idx.divE, 0 ); } diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index e7764627b..aac4035a1 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -29,33 +29,74 @@ // Declare type for spectral fields using SpectralField = amrex::FabArray< amrex::BaseFab <Complex> >; -/** Index for the regular fields, when stored in spectral space: - * - n_fields is automatically the total number of fields - * - divE reuses the memory slot for Bx, since Bx is not used when computing divE - */ -struct SpectralFieldIndex { - enum { Ex=0, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz, rho_old, rho_new, n_fields, divE=3 }; -}; - -struct SpectralFieldIndexJLinearInTime { - enum { Ex=0, Ey, Ez, Bx, By, Bz, Jx_old, Jy_old, Jz_old, rho_old, rho_new, - Ex_avg, Ey_avg, Ez_avg, Bx_avg, By_avg, Bz_avg, - Jx_new, Jy_new, Jz_new, F, G, n_fields, divE=3 }; -}; - -/* Index for the regular fields + averaged fields, when stored in spectral space */ -struct SpectralFieldIndexTimeAveraging { - enum { Ex=0, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz, rho_old, rho_new, Ex_avg, Ey_avg, Ez_avg, Bx_avg, By_avg, Bz_avg,n_fields }; - // n_fields is automatically the total number of fields -}; +class SpectralFieldIndex +{ + public: -/** Index for the PML fields, when stored in spectral space, - * (n_fields is automatically set to the total number of fields) - * TODO How to include the diagonal components only when needed? - */ -struct SpectralPMLIndex { - enum {Exx=0, Exy, Exz, Eyx, Eyy, Eyz, Ezx, Ezy, Ezz, - Bxx , Bxy, Bxz, Byx, Byy, Byz, Bzx, Bzy, Bzz, Fx, Fy, Fz, Gx, Gy, Gz, n_fields}; + /** + * \brief Constructor of the class SpectralFieldIndex + * + * Set integer indices to access data in spectral space + * and total number of fields to be stored. + * + * \param[in] update_with_rho whether rho is used in the field update equations + * \param[in] time_averaging whether the time averaging algorithm is used + * \param[in] J_linear_in_time whether to use two currents computed at the beginning and + * the end of the time interval (instead of using one current + * computed at half time) + * \param[in] dive_cleaning whether to use div(E) cleaning to account for errors in + * Gauss law (new field F in the update equations) + * \param[in] divb_cleaning whether to use div(B) cleaning to account for errors in + * div(B) = 0 law (new field G in the update equations) + * \param[in] pml whether the indices are used to access spectral data + * for the PML spectral solver + */ + SpectralFieldIndex (const bool update_with_rho, + const bool time_averaging, + const bool J_linear_in_time, + const bool dive_cleaning, + const bool divb_cleaning, + const bool pml); + + /** + * \brief Default constructor + */ + SpectralFieldIndex () = default; + + /** + * \brief Default destructor + */ + ~SpectralFieldIndex () = default; + + // Total number of fields that are actually allocated + int n_fields; + + // Indices overwritten in the constructor, for the fields that are actually allocated + // (index -1 will never be used, unless there is some bug in the code implementation, + // which would result in a runtime crash due to out-of-bound accesses that can be detected + // by running the code in DEBUG mode) + + // Always + int Ex = -1, Ey = -1, Ez = -1; + int Bx = -1, By = -1, Bz = -1; + int Jx = -1, Jy = -1, Jz = -1; + int rho_old = -1, rho_new = -1, divE = -1; + + // Time averaging + int Ex_avg = -1, Ey_avg = -1, Ez_avg = -1; + int Bx_avg = -1, By_avg = -1, Bz_avg = -1; + + // J linear in time + int Jx_new = -1, Jy_new = -1, Jz_new = -1; + int F = -1, G = -1; + + // PML + int Exy = -1, Exz = -1, Eyx = -1, Eyz = -1, Ezx = -1, Ezy = -1; + int Bxy = -1, Bxz = -1, Byx = -1, Byz = -1, Bzx = -1, Bzy = -1; + + // PML with div(E) and/or div(B) cleaning + int Exx = -1, Eyy = -1, Ezz = -1, Bxx = -1, Byy = -1, Bzz = -1; + int Fx = -1, Fy = -1, Fz = -1, Gx = -1, Gy = -1, Gz = -1; }; /** \brief Class that stores the fields in spectral space, and performs the diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index b26ac4943..e50eae92e 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -32,6 +32,66 @@ using namespace amrex; +SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho, + const bool time_averaging, + const bool J_linear_in_time, + const bool dive_cleaning, + const bool divb_cleaning, + const bool pml) +{ + // TODO Use these to allocate rho_old, rho_new, F, and G only when needed + amrex::ignore_unused(update_with_rho); + + int c = 0; + + if (pml == false) + { + Ex = c++; Ey = c++; Ez = c++; + Bx = c++; By = c++; Bz = c++; + Jx = c++; Jy = c++; Jz = c++; + + // TODO Allocate rho_old and rho_new only when needed + rho_old = c++; rho_new = c++; + + // Reuse data corresponding to index Bx = 3 to avoid storing extra memory + divE = 3; + + if (time_averaging) + { + Ex_avg = c++; Ey_avg = c++; Ez_avg = c++; + Bx_avg = c++; By_avg = c++; Bz_avg = c++; + } + + if (J_linear_in_time) + { + Jx_new = c++; Jy_new = c++; Jz_new = c++; + + // TODO Allocate F and G only when needed + F = c++; G = c++; + } + } + else // PML + { + Exy = c++; Exz = c++; Eyx = c++; Eyz = c++; Ezx = c++; Ezy = c++; + Bxy = c++; Bxz = c++; Byx = c++; Byz = c++; Bzx = c++; Bzy = c++; + + if (dive_cleaning) + { + Exx = c++; Eyy = c++; Ezz = c++; + Fx = c++; Fy = c++; Fz = c++; + } + + if (divb_cleaning) + { + Bxx = c++; Byy = c++; Bzz = c++; + Gx = c++; Gy = c++; Gz = c++; + } + } + + // This is the number of arrays that will be actually allocated in spectral space + n_fields = c; +} + /* \brief Initialize fields in spectral space, and FFT plans */ SpectralFieldData::SpectralFieldData( const int lev, const amrex::BoxArray& realspace_ba, diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.H b/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.H index 87158314d..bf0d6b780 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.H @@ -83,6 +83,9 @@ class SpectralFieldDataRZ private: + SpectralFieldIndex m_spectral_index; + int m_n_fields; + // tempHTransformed and tmpSpectralField store fields // right before/after the z Fourier transform SpectralField tempHTransformed; // contains Complexes diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp index 80760afb3..9604ccb44 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp @@ -25,7 +25,8 @@ SpectralFieldDataRZ::SpectralFieldDataRZ (const int lev, amrex::DistributionMapping const & dm, int const n_field_required, int const n_modes) - : n_rz_azimuthal_modes(n_modes) + : n_rz_azimuthal_modes(n_modes), + m_n_fields(n_field_required) { amrex::BoxArray const & spectralspace_ba = k_space.spectralspace_ba; @@ -281,7 +282,7 @@ SpectralFieldDataRZ::FABZForwardTransform (amrex::MFIter const & mfi, amrex::Box amrex::Box const& spectralspace_bx = tmpSpectralField[mfi].box(); int const nz = spectralspace_bx.length(1); amrex::Real inv_nz = 1._rt/nz; - constexpr int n_fields = SpectralFieldIndex::n_fields; + const int n_fields = m_n_fields; ParallelFor(spectralspace_bx, modes, [=] AMREX_GPU_DEVICE(int i, int j, int k, int mode) noexcept { @@ -319,7 +320,7 @@ SpectralFieldDataRZ::FABZBackwardTransform (amrex::MFIter const & mfi, amrex::Bo amrex::Box const& spectralspace_bx = tmpSpectralField[mfi].box(); int const modes = n_rz_azimuthal_modes; - constexpr int n_fields = SpectralFieldIndex::n_fields; + const int n_fields = m_n_fields; ParallelFor(spectralspace_bx, modes, [=] AMREX_GPU_DEVICE(int i, int j, int k, int mode) noexcept { int const ic = field_index + mode*n_fields; @@ -735,7 +736,7 @@ SpectralFieldDataRZ::ApplyFilter (const int lev, int const field_index) amrex::Array4<Complex> const& fields_arr = fields[mfi].array(); int const modes = n_rz_azimuthal_modes; - constexpr int n_fields = SpectralFieldIndex::n_fields; + const int n_fields = m_n_fields; amrex::Box const& spectralspace_bx = fields[mfi].box(); int const nr = spectralspace_bx.length(0); @@ -779,7 +780,7 @@ SpectralFieldDataRZ::ApplyFilter (const int lev, int const field_index1, amrex::Array4<Complex> const& fields_arr = fields[mfi].array(); int const modes = n_rz_azimuthal_modes; - constexpr int n_fields = SpectralFieldIndex::n_fields; + const int n_fields = m_n_fields; amrex::Box const& spectralspace_bx = fields[mfi].box(); int const nr = spectralspace_bx.length(0); diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData_fwd.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData_fwd.H index 631dd3b06..d01cc0a83 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData_fwd.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData_fwd.H @@ -5,8 +5,5 @@ * License: BSD-3-Clause-LBNL */ -struct SpectralFieldIndex; -struct SpectralAvgFieldIndex; -struct SpectralPMLIndex; - +class SpectralFieldIndex; class SpectralFieldData; diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index 45ba1a193..50b9978f4 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -60,11 +60,11 @@ class SpectralSolver * \param[in] fft_do_time_averaging whether the time averaging algorithm is used * \param[in] J_linear_in_time whether to use two currents computed at the beginning and * the end of the time interval (instead of using one current - * compute at half time) - * \param[in] dive_cleaning whether to use div(E) cleaning to account for errors in Gauss law - * (new field F in the field update equations) - * \param[in] divb_cleaning whether to use div(B) cleaning to account for errors in magnetic Gauss law - * (new field G in the field update equations) + * computed at half time) + * \param[in] dive_cleaning whether to use div(E) cleaning to account for errors in + * Gauss law (new field F in the update equations) + * \param[in] divb_cleaning whether to use div(B) cleaning to account for errors in + * div(B) = 0 law (new field G in the update equations) */ SpectralSolver (const int lev, const amrex::BoxArray& realspace_ba, @@ -76,13 +76,13 @@ class SpectralSolver const amrex::Array<amrex::Real,3>& v_comoving, const amrex::RealVect dx, const amrex::Real dt, - const bool pml = false, - const bool periodic_single_box = false, - const bool update_with_rho = false, - const bool fft_do_time_averaging = false, - const bool J_linear_in_time = false, - const bool dive_cleaning = false, - const bool divb_cleaning = false); + const bool pml, + const bool periodic_single_box, + const bool update_with_rho, + const bool fft_do_time_averaging, + const bool J_linear_in_time, + const bool dive_cleaning, + const bool divb_cleaning); /** * \brief Transform the component `i_comp` of MultiFab `mf` @@ -185,6 +185,8 @@ class SpectralSolver field_data.fields.mult(scale_factor, icomp, 1); } + SpectralFieldIndex m_spectral_index; + protected: amrex::IntVect m_fill_guards; diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp index 113ea97c3..29f2b7256 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -41,32 +41,35 @@ SpectralSolver::SpectralSolver( // as well as the value of the corresponding k coordinates) const SpectralKSpace k_space= SpectralKSpace(realspace_ba, dm, dx); + m_spectral_index = SpectralFieldIndex(update_with_rho, fft_do_time_averaging, + J_linear_in_time, dive_cleaning, divb_cleaning, pml); + // - Select the algorithm depending on the input parameters // Initialize the corresponding coefficients over k space if (pml) { algorithm = std::make_unique<PMLPsatdAlgorithm>( - k_space, dm, norder_x, norder_y, norder_z, nodal, + k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal, fill_guards, dt, dive_cleaning, divb_cleaning); } else { // Comoving PSATD algorithm if (v_comoving[0] != 0. || v_comoving[1] != 0. || v_comoving[2] != 0.) { algorithm = std::make_unique<ComovingPsatdAlgorithm>( - k_space, dm, norder_x, norder_y, norder_z, nodal, + k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal, fill_guards, v_comoving, dt, update_with_rho); } // PSATD algorithms: standard, Galilean, or averaged Galilean else { algorithm = std::make_unique<PsatdAlgorithm>( - k_space, dm, norder_x, norder_y, norder_z, nodal, fill_guards, + k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal, fill_guards, v_galilean, dt, update_with_rho, fft_do_time_averaging, J_linear_in_time); } } // - Initialize arrays for fields in spectral space + FFT plans - field_data = SpectralFieldData( lev, realspace_ba, k_space, dm, - algorithm->getRequiredNumberOfFields(), periodic_single_box); + field_data = SpectralFieldData(lev, realspace_ba, k_space, dm, + m_spectral_index.n_fields, periodic_single_box); m_fill_guards = fill_guards; } diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H index 1c6b3d613..8e84f7c00 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H @@ -112,6 +112,8 @@ class SpectralSolverRZ */ void VayDeposition (const int lev, std::array<std::unique_ptr<amrex::MultiFab>,3>& current); + SpectralFieldIndex m_spectral_index; + private: SpectralKSpaceRZ k_space; // Save the instance to initialize filtering diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp index ee4e5fd6d..8f46d21cd 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp @@ -41,22 +41,30 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev, // the spectral space corresponding to each box in `realspace_ba`, // as well as the value of the corresponding k coordinates. + const bool time_averaging = false; + const bool J_linear_in_time = false; + const bool dive_cleaning = false; + const bool divb_cleaning = false; + const bool pml = false; + m_spectral_index = SpectralFieldIndex(update_with_rho, time_averaging, J_linear_in_time, + dive_cleaning, divb_cleaning, pml); + // - Select the algorithm depending on the input parameters // Initialize the corresponding coefficients over k space // PML is not supported. if (v_galilean[2] == 0) { // v_galilean is 0: use standard PSATD algorithm algorithm = std::make_unique<PsatdAlgorithmRZ>( - k_space, dm, n_rz_azimuthal_modes, norder_z, nodal, dt, update_with_rho); + k_space, dm, m_spectral_index, n_rz_azimuthal_modes, norder_z, nodal, dt, update_with_rho); } else { // Otherwise: use the Galilean algorithm algorithm = std::make_unique<GalileanPsatdAlgorithmRZ>( - k_space, dm, n_rz_azimuthal_modes, norder_z, nodal, v_galilean, dt, update_with_rho); + k_space, dm, m_spectral_index, n_rz_azimuthal_modes, norder_z, nodal, v_galilean, dt, update_with_rho); } // - Initialize arrays for fields in spectral space + FFT plans field_data = SpectralFieldDataRZ(lev, realspace_ba, k_space, dm, - algorithm->getRequiredNumberOfFields(), + m_spectral_index.n_fields, n_rz_azimuthal_modes); } diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 29d244d12..0715e90c8 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -95,21 +95,20 @@ namespace { } } -using IdxAvg = SpectralFieldIndexTimeAveraging; -using IdxLin = SpectralFieldIndexJLinearInTime; - void WarpX::PSATDForwardTransformEB () { + const SpectralFieldIndex& Idx = spectral_solver_fp[0]->m_spectral_index; + for (int lev = 0; lev <= finest_level; ++lev) { - ForwardTransformVect(lev, *spectral_solver_fp[lev], Efield_fp[lev], IdxAvg::Ex, IdxAvg::Ey, IdxAvg::Ez); - ForwardTransformVect(lev, *spectral_solver_fp[lev], Bfield_fp[lev], IdxAvg::Bx, IdxAvg::By, IdxAvg::Bz); + ForwardTransformVect(lev, *spectral_solver_fp[lev], Efield_fp[lev], Idx.Ex, Idx.Ey, Idx.Ez); + ForwardTransformVect(lev, *spectral_solver_fp[lev], Bfield_fp[lev], Idx.Bx, Idx.By, Idx.Bz); if (spectral_solver_cp[lev]) { - ForwardTransformVect(lev, *spectral_solver_cp[lev], Efield_cp[lev], IdxAvg::Ex, IdxAvg::Ey, IdxAvg::Ez); - ForwardTransformVect(lev, *spectral_solver_cp[lev], Bfield_cp[lev], IdxAvg::Bx, IdxAvg::By, IdxAvg::Bz); + ForwardTransformVect(lev, *spectral_solver_cp[lev], Efield_cp[lev], Idx.Ex, Idx.Ey, Idx.Ez); + ForwardTransformVect(lev, *spectral_solver_cp[lev], Bfield_cp[lev], Idx.Bx, Idx.By, Idx.Bz); } } } @@ -117,15 +116,17 @@ WarpX::PSATDForwardTransformEB () void WarpX::PSATDBackwardTransformEB () { + const SpectralFieldIndex& Idx = spectral_solver_fp[0]->m_spectral_index; + for (int lev = 0; lev <= finest_level; ++lev) { - BackwardTransformVect(lev, *spectral_solver_fp[lev], Efield_fp[lev], IdxAvg::Ex, IdxAvg::Ey, IdxAvg::Ez); - BackwardTransformVect(lev, *spectral_solver_fp[lev], Bfield_fp[lev], IdxAvg::Bx, IdxAvg::By, IdxAvg::Bz); + BackwardTransformVect(lev, *spectral_solver_fp[lev], Efield_fp[lev], Idx.Ex, Idx.Ey, Idx.Ez); + BackwardTransformVect(lev, *spectral_solver_fp[lev], Bfield_fp[lev], Idx.Bx, Idx.By, Idx.Bz); if (spectral_solver_cp[lev]) { - BackwardTransformVect(lev, *spectral_solver_cp[lev], Efield_cp[lev], IdxAvg::Ex, IdxAvg::Ey, IdxAvg::Ez); - BackwardTransformVect(lev, *spectral_solver_cp[lev], Bfield_cp[lev], IdxAvg::Bx, IdxAvg::By, IdxAvg::Bz); + BackwardTransformVect(lev, *spectral_solver_cp[lev], Efield_cp[lev], Idx.Ex, Idx.Ey, Idx.Ez); + BackwardTransformVect(lev, *spectral_solver_cp[lev], Bfield_cp[lev], Idx.Bx, Idx.By, Idx.Bz); } } @@ -144,15 +145,17 @@ WarpX::PSATDBackwardTransformEB () void WarpX::PSATDBackwardTransformEBavg () { + const SpectralFieldIndex& Idx = spectral_solver_fp[0]->m_spectral_index; + for (int lev = 0; lev <= finest_level; ++lev) { - BackwardTransformVect(lev, *spectral_solver_fp[lev], Efield_avg_fp[lev], IdxAvg::Ex_avg, IdxAvg::Ey_avg, IdxAvg::Ez_avg); - BackwardTransformVect(lev, *spectral_solver_fp[lev], Bfield_avg_fp[lev], IdxAvg::Bx_avg, IdxAvg::By_avg, IdxAvg::Bz_avg); + BackwardTransformVect(lev, *spectral_solver_fp[lev], Efield_avg_fp[lev], Idx.Ex_avg, Idx.Ey_avg, Idx.Ez_avg); + BackwardTransformVect(lev, *spectral_solver_fp[lev], Bfield_avg_fp[lev], Idx.Bx_avg, Idx.By_avg, Idx.Bz_avg); if (spectral_solver_cp[lev]) { - BackwardTransformVect(lev, *spectral_solver_cp[lev], Efield_avg_cp[lev], IdxAvg::Ex_avg, IdxAvg::Ey_avg, IdxAvg::Ez_avg); - BackwardTransformVect(lev, *spectral_solver_cp[lev], Bfield_avg_cp[lev], IdxAvg::Bx_avg, IdxAvg::By_avg, IdxAvg::Bz_avg); + BackwardTransformVect(lev, *spectral_solver_cp[lev], Efield_avg_cp[lev], Idx.Ex_avg, Idx.Ey_avg, Idx.Ez_avg); + BackwardTransformVect(lev, *spectral_solver_cp[lev], Bfield_avg_cp[lev], Idx.Bx_avg, Idx.By_avg, Idx.Bz_avg); } } } @@ -160,13 +163,15 @@ WarpX::PSATDBackwardTransformEBavg () void WarpX::PSATDForwardTransformF () { + const SpectralFieldIndex& Idx = spectral_solver_fp[0]->m_spectral_index; + for (int lev = 0; lev <= finest_level; ++lev) { - if (F_fp[lev]) spectral_solver_fp[lev]->ForwardTransform(lev, *F_fp[lev], IdxLin::F); + if (F_fp[lev]) spectral_solver_fp[lev]->ForwardTransform(lev, *F_fp[lev], Idx.F); if (spectral_solver_cp[lev]) { - if (F_cp[lev]) spectral_solver_cp[lev]->ForwardTransform(lev, *F_cp[lev], IdxLin::F); + if (F_cp[lev]) spectral_solver_cp[lev]->ForwardTransform(lev, *F_cp[lev], Idx.F); } } } @@ -174,13 +179,15 @@ WarpX::PSATDForwardTransformF () void WarpX::PSATDBackwardTransformF () { + const SpectralFieldIndex& Idx = spectral_solver_fp[0]->m_spectral_index; + for (int lev = 0; lev <= finest_level; ++lev) { - if (F_fp[lev]) spectral_solver_fp[lev]->BackwardTransform(lev, *F_fp[lev], IdxLin::F); + if (F_fp[lev]) spectral_solver_fp[lev]->BackwardTransform(lev, *F_fp[lev], Idx.F); if (spectral_solver_cp[lev]) { - if (F_cp[lev]) spectral_solver_cp[lev]->BackwardTransform(lev, *F_cp[lev], IdxLin::F); + if (F_cp[lev]) spectral_solver_cp[lev]->BackwardTransform(lev, *F_cp[lev], Idx.F); } } } @@ -188,13 +195,15 @@ WarpX::PSATDBackwardTransformF () void WarpX::PSATDForwardTransformG () { + const SpectralFieldIndex& Idx = spectral_solver_fp[0]->m_spectral_index; + for (int lev = 0; lev <= finest_level; ++lev) { - if (G_fp[lev]) spectral_solver_fp[lev]->ForwardTransform(lev, *G_fp[lev], IdxLin::G); + if (G_fp[lev]) spectral_solver_fp[lev]->ForwardTransform(lev, *G_fp[lev], Idx.G); if (spectral_solver_cp[lev]) { - if (G_cp[lev]) spectral_solver_cp[lev]->ForwardTransform(lev, *G_cp[lev], IdxLin::G); + if (G_cp[lev]) spectral_solver_cp[lev]->ForwardTransform(lev, *G_cp[lev], Idx.G); } } } @@ -202,13 +211,15 @@ WarpX::PSATDForwardTransformG () void WarpX::PSATDBackwardTransformG () { + const SpectralFieldIndex& Idx = spectral_solver_fp[0]->m_spectral_index; + for (int lev = 0; lev <= finest_level; ++lev) { - if (G_fp[lev]) spectral_solver_fp[lev]->BackwardTransform(lev, *G_fp[lev], IdxLin::G); + if (G_fp[lev]) spectral_solver_fp[lev]->BackwardTransform(lev, *G_fp[lev], Idx.G); if (spectral_solver_cp[lev]) { - if (G_cp[lev]) spectral_solver_cp[lev]->BackwardTransform(lev, *G_cp[lev], IdxLin::G); + if (G_cp[lev]) spectral_solver_cp[lev]->BackwardTransform(lev, *G_cp[lev], Idx.G); } } } @@ -216,12 +227,14 @@ WarpX::PSATDBackwardTransformG () void WarpX::PSATDForwardTransformJ () { - const int idx_jx = (WarpX::J_linear_in_time) ? static_cast<int>(IdxLin::Jx_new) - : static_cast<int>(IdxAvg::Jx); - const int idx_jy = (WarpX::J_linear_in_time) ? static_cast<int>(IdxLin::Jy_new) - : static_cast<int>(IdxAvg::Jy); - const int idx_jz = (WarpX::J_linear_in_time) ? static_cast<int>(IdxLin::Jz_new) - : static_cast<int>(IdxAvg::Jz); + const SpectralFieldIndex& Idx = spectral_solver_fp[0]->m_spectral_index; + + const int idx_jx = (WarpX::J_linear_in_time) ? static_cast<int>(Idx.Jx_new) + : static_cast<int>(Idx.Jx); + const int idx_jy = (WarpX::J_linear_in_time) ? static_cast<int>(Idx.Jy_new) + : static_cast<int>(Idx.Jy); + const int idx_jz = (WarpX::J_linear_in_time) ? static_cast<int>(Idx.Jz_new) + : static_cast<int>(Idx.Jz); for (int lev = 0; lev <= finest_level; ++lev) { @@ -239,11 +252,11 @@ WarpX::PSATDForwardTransformJ () { for (int lev = 0; lev <= finest_level; ++lev) { - spectral_solver_fp[lev]->ApplyFilter(lev, IdxAvg::Jx, IdxAvg::Jy, IdxAvg::Jz); + spectral_solver_fp[lev]->ApplyFilter(lev, Idx.Jx, Idx.Jy, Idx.Jz); if (spectral_solver_cp[lev]) { - spectral_solver_cp[lev]->ApplyFilter(lev, IdxAvg::Jx, IdxAvg::Jy, IdxAvg::Jz); + spectral_solver_cp[lev]->ApplyFilter(lev, Idx.Jx, Idx.Jy, Idx.Jz); } } } @@ -253,8 +266,10 @@ WarpX::PSATDForwardTransformJ () void WarpX::PSATDForwardTransformRho (const int icomp) { + const SpectralFieldIndex& Idx = spectral_solver_fp[0]->m_spectral_index; + // Select index in k space - const int dst_comp = (icomp == 0) ? IdxAvg::rho_old : IdxAvg::rho_new; + const int dst_comp = (icomp == 0) ? Idx.rho_old : Idx.rho_new; for (int lev = 0; lev <= finest_level; ++lev) { @@ -301,13 +316,15 @@ WarpX::PSATDPushSpectralFields () void WarpX::PSATDMoveRhoNewToRhoOld () { + const SpectralFieldIndex& Idx = spectral_solver_fp[0]->m_spectral_index; + for (int lev = 0; lev <= finest_level; ++lev) { - spectral_solver_fp[lev]->CopySpectralDataComp(IdxAvg::rho_new, IdxAvg::rho_old); + spectral_solver_fp[lev]->CopySpectralDataComp(Idx.rho_new, Idx.rho_old); if (spectral_solver_cp[lev]) { - spectral_solver_cp[lev]->CopySpectralDataComp(IdxAvg::rho_new, IdxAvg::rho_old); + spectral_solver_cp[lev]->CopySpectralDataComp(Idx.rho_new, Idx.rho_old); } } } @@ -315,17 +332,19 @@ WarpX::PSATDMoveRhoNewToRhoOld () void WarpX::PSATDMoveJNewToJOld () { + const SpectralFieldIndex& Idx = spectral_solver_fp[0]->m_spectral_index; + for (int lev = 0; lev <= finest_level; ++lev) { - spectral_solver_fp[lev]->CopySpectralDataComp(IdxLin::Jx_new, IdxLin::Jx_old); - spectral_solver_fp[lev]->CopySpectralDataComp(IdxLin::Jy_new, IdxLin::Jy_old); - spectral_solver_fp[lev]->CopySpectralDataComp(IdxLin::Jz_new, IdxLin::Jz_old); + spectral_solver_fp[lev]->CopySpectralDataComp(Idx.Jx_new, Idx.Jx); + spectral_solver_fp[lev]->CopySpectralDataComp(Idx.Jy_new, Idx.Jy); + spectral_solver_fp[lev]->CopySpectralDataComp(Idx.Jz_new, Idx.Jz); if (spectral_solver_cp[lev]) { - spectral_solver_cp[lev]->CopySpectralDataComp(IdxLin::Jx_new, IdxLin::Jx_old); - spectral_solver_cp[lev]->CopySpectralDataComp(IdxLin::Jy_new, IdxLin::Jy_old); - spectral_solver_cp[lev]->CopySpectralDataComp(IdxLin::Jz_new, IdxLin::Jz_old); + spectral_solver_cp[lev]->CopySpectralDataComp(Idx.Jx_new, Idx.Jx); + spectral_solver_cp[lev]->CopySpectralDataComp(Idx.Jy_new, Idx.Jy); + spectral_solver_cp[lev]->CopySpectralDataComp(Idx.Jz_new, Idx.Jz); } } } @@ -333,23 +352,25 @@ WarpX::PSATDMoveJNewToJOld () void WarpX::PSATDEraseAverageFields () { + const SpectralFieldIndex& Idx = spectral_solver_fp[0]->m_spectral_index; + for (int lev = 0; lev <= finest_level; ++lev) { - spectral_solver_fp[lev]->ZeroOutDataComp(IdxAvg::Ex_avg); - spectral_solver_fp[lev]->ZeroOutDataComp(IdxAvg::Ey_avg); - spectral_solver_fp[lev]->ZeroOutDataComp(IdxAvg::Ez_avg); - spectral_solver_fp[lev]->ZeroOutDataComp(IdxAvg::Bx_avg); - spectral_solver_fp[lev]->ZeroOutDataComp(IdxAvg::By_avg); - spectral_solver_fp[lev]->ZeroOutDataComp(IdxAvg::Bz_avg); + spectral_solver_fp[lev]->ZeroOutDataComp(Idx.Ex_avg); + spectral_solver_fp[lev]->ZeroOutDataComp(Idx.Ey_avg); + spectral_solver_fp[lev]->ZeroOutDataComp(Idx.Ez_avg); + spectral_solver_fp[lev]->ZeroOutDataComp(Idx.Bx_avg); + spectral_solver_fp[lev]->ZeroOutDataComp(Idx.By_avg); + spectral_solver_fp[lev]->ZeroOutDataComp(Idx.Bz_avg); if (spectral_solver_cp[lev]) { - spectral_solver_cp[lev]->ZeroOutDataComp(IdxAvg::Ex_avg); - spectral_solver_cp[lev]->ZeroOutDataComp(IdxAvg::Ey_avg); - spectral_solver_cp[lev]->ZeroOutDataComp(IdxAvg::Ez_avg); - spectral_solver_cp[lev]->ZeroOutDataComp(IdxAvg::Bx_avg); - spectral_solver_cp[lev]->ZeroOutDataComp(IdxAvg::By_avg); - spectral_solver_cp[lev]->ZeroOutDataComp(IdxAvg::Bz_avg); + spectral_solver_cp[lev]->ZeroOutDataComp(Idx.Ex_avg); + spectral_solver_cp[lev]->ZeroOutDataComp(Idx.Ey_avg); + spectral_solver_cp[lev]->ZeroOutDataComp(Idx.Ez_avg); + spectral_solver_cp[lev]->ZeroOutDataComp(Idx.Bx_avg); + spectral_solver_cp[lev]->ZeroOutDataComp(Idx.By_avg); + spectral_solver_cp[lev]->ZeroOutDataComp(Idx.Bz_avg); } } } @@ -357,23 +378,25 @@ WarpX::PSATDEraseAverageFields () void WarpX::PSATDScaleAverageFields (const amrex::Real scale_factor) { + const SpectralFieldIndex& Idx = spectral_solver_fp[0]->m_spectral_index; + for (int lev = 0; lev <= finest_level; ++lev) { - spectral_solver_fp[lev]->ScaleDataComp(IdxAvg::Ex_avg, scale_factor); - spectral_solver_fp[lev]->ScaleDataComp(IdxAvg::Ey_avg, scale_factor); - spectral_solver_fp[lev]->ScaleDataComp(IdxAvg::Ez_avg, scale_factor); - spectral_solver_fp[lev]->ScaleDataComp(IdxAvg::Bx_avg, scale_factor); - spectral_solver_fp[lev]->ScaleDataComp(IdxAvg::By_avg, scale_factor); - spectral_solver_fp[lev]->ScaleDataComp(IdxAvg::Bz_avg, scale_factor); + spectral_solver_fp[lev]->ScaleDataComp(Idx.Ex_avg, scale_factor); + spectral_solver_fp[lev]->ScaleDataComp(Idx.Ey_avg, scale_factor); + spectral_solver_fp[lev]->ScaleDataComp(Idx.Ez_avg, scale_factor); + spectral_solver_fp[lev]->ScaleDataComp(Idx.Bx_avg, scale_factor); + spectral_solver_fp[lev]->ScaleDataComp(Idx.By_avg, scale_factor); + spectral_solver_fp[lev]->ScaleDataComp(Idx.Bz_avg, scale_factor); if (spectral_solver_cp[lev]) { - spectral_solver_cp[lev]->ScaleDataComp(IdxAvg::Ex_avg, scale_factor); - spectral_solver_cp[lev]->ScaleDataComp(IdxAvg::Ey_avg, scale_factor); - spectral_solver_cp[lev]->ScaleDataComp(IdxAvg::Ez_avg, scale_factor); - spectral_solver_cp[lev]->ScaleDataComp(IdxAvg::Bx_avg, scale_factor); - spectral_solver_cp[lev]->ScaleDataComp(IdxAvg::By_avg, scale_factor); - spectral_solver_cp[lev]->ScaleDataComp(IdxAvg::Bz_avg, scale_factor); + spectral_solver_cp[lev]->ScaleDataComp(Idx.Ex_avg, scale_factor); + spectral_solver_cp[lev]->ScaleDataComp(Idx.Ey_avg, scale_factor); + spectral_solver_cp[lev]->ScaleDataComp(Idx.Ez_avg, scale_factor); + spectral_solver_cp[lev]->ScaleDataComp(Idx.Bx_avg, scale_factor); + spectral_solver_cp[lev]->ScaleDataComp(Idx.By_avg, scale_factor); + spectral_solver_cp[lev]->ScaleDataComp(Idx.Bz_avg, scale_factor); } } } |