diff options
Diffstat (limited to 'Source/FieldSolver/SpectralSolver')
14 files changed, 458 insertions, 164 deletions
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H index 45a35098b..9d3facf3f 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H @@ -27,6 +27,22 @@ class AvgGalileanAlgorithm : public SpectralBaseAlgorithm const amrex::Real dt); /** + * \brief Virtual function for current correction in Fourier space + * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>). + * This function overrides the virtual function \c CurrentCorrection in the + * base class \c SpectralBaseAlgorithm and cannot be overridden by further + * derived classes. + * + * \param[in,out] field_data All fields in Fourier space + * \param[in,out] current Array of unique pointers to \c MultiFab storing + * the three components of the current density + * \param[in] rho Unique pointer to \c MultiFab storing the charge density + */ + virtual void CurrentCorrection (SpectralFieldData& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current, + const std::unique_ptr<amrex::MultiFab>& rho) override final; + + /** * \brief Virtual function for Vay current deposition in Fourier space * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>). * This function overrides the virtual function \c VayDeposition in the diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp index 0ba603031..e5fbc8261 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp @@ -371,6 +371,14 @@ AvgGalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ }; void +AvgGalileanAlgorithm::CurrentCorrection (SpectralFieldData& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current, + const std::unique_ptr<amrex::MultiFab>& rho) +{ + amrex::Abort("Current correction not implemented for averaged Galilean PSATD"); +} + +void AvgGalileanAlgorithm::VayDeposition (SpectralFieldData& field_data, std::array<std::unique_ptr<amrex::MultiFab>,3>& current) { diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H index 799eca8b1..83ae83d69 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H @@ -10,21 +10,37 @@ class GalileanAlgorithm : public SpectralBaseAlgorithm { public: - GalileanAlgorithm(const SpectralKSpace& spectral_kspace, - const amrex::DistributionMapping& dm, - const int norder_x, const int norder_y, - const int norder_z, const bool nodal, - const amrex::Array<amrex::Real,3>& v_galilean, - const amrex::Real dt); + GalileanAlgorithm (const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const int norder_x, const int norder_y, + const int norder_z, const bool nodal, + const amrex::Array<amrex::Real,3>& v_galilean, + const amrex::Real dt, + const bool update_with_rho); // Redefine update equation from base class - virtual void pushSpectralFields(SpectralFieldData& f) const override final; - virtual int getRequiredNumberOfFields() const override final { + virtual void pushSpectralFields (SpectralFieldData& f) const override final; + virtual int getRequiredNumberOfFields () const override final { return SpectralFieldIndex::n_fields; }; - void InitializeSpectralCoefficients(const SpectralKSpace& spectral_kspace, - const amrex::DistributionMapping& dm, - const amrex::Array<amrex::Real, 3>& v_galilean, - const amrex::Real dt); + void InitializeSpectralCoefficients (const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const amrex::Real dt); + + /** + * \brief Virtual function for current correction in Fourier space + * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>). + * This function overrides the virtual function \c CurrentCorrection in the + * base class \c SpectralBaseAlgorithm and cannot be overridden by further + * derived classes. + * + * \param[in,out] field_data All fields in Fourier space + * \param[in,out] current Array of unique pointers to \c MultiFab storing + * the three components of the current density + * \param[in] rho Unique pointer to \c MultiFab storing the charge density + */ + virtual void CurrentCorrection (SpectralFieldData& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current, + const std::unique_ptr<amrex::MultiFab>& rho) override final; /** * \brief Virtual function for Vay current deposition in Fourier space @@ -43,6 +59,9 @@ class GalileanAlgorithm : public SpectralBaseAlgorithm private: SpectralRealCoefficients C_coef, S_ck_coef; SpectralComplexCoefficients Theta2_coef, X1_coef, X2_coef, X3_coef, X4_coef; + amrex::Array<amrex::Real,3> m_v_galilean; + amrex::Real m_dt; + bool m_update_with_rho; }; #endif // WARPX_USE_PSATD #endif // WARPX_GALILEAN_ALGORITHM_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp index e606b3232..e969fd6cf 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp @@ -14,10 +14,13 @@ GalileanAlgorithm::GalileanAlgorithm(const SpectralKSpace& spectral_kspace, const int norder_x, const int norder_y, const int norder_z, const bool nodal, const Array<Real, 3>& v_galilean, - const Real dt) + const Real dt, + const bool update_with_rho) // Initialize members of base class - : SpectralBaseAlgorithm( spectral_kspace, dm, - norder_x, norder_y, norder_z, nodal ) + : m_v_galilean(v_galilean), + m_dt(dt), + m_update_with_rho(update_with_rho), + SpectralBaseAlgorithm(spectral_kspace, dm, norder_x, norder_y, norder_z, nodal) { const BoxArray& ba = spectral_kspace.spectralspace_ba; @@ -30,14 +33,14 @@ GalileanAlgorithm::GalileanAlgorithm(const SpectralKSpace& spectral_kspace, X4_coef = SpectralComplexCoefficients(ba, dm, 1, 0); Theta2_coef = SpectralComplexCoefficients(ba, dm, 1, 0); - InitializeSpectralCoefficients(spectral_kspace, dm, v_galilean, dt); - + InitializeSpectralCoefficients(spectral_kspace, dm, dt); }; -/* Advance the E and B field in spectral space (stored in `f`) - * over one time step */ +/* Advance the E and B field in spectral space (stored in `f`) over one time step */ void -GalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ +GalileanAlgorithm::pushSpectralFields (SpectralFieldData& f) const +{ + const bool update_with_rho = m_update_with_rho; // Loop over boxes for (MFIter mfi(f.fields); mfi.isValid(); ++mfi){ @@ -46,6 +49,7 @@ GalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ // Extract arrays for the fields to be updated Array4<Complex> fields = f.fields[mfi].array(); + // Extract arrays for the coefficients Array4<const Real> C_arr = C_coef[mfi].array(); Array4<const Real> S_ck_arr = S_ck_coef[mfi].array(); @@ -63,8 +67,7 @@ GalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ const Real* modified_kz_arr = modified_kz_vec[mfi].dataPtr(); // Loop over indices within one box - ParallelFor(bx, - [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Record old values of the fields to be updated using Idx = SpectralFieldIndex; @@ -74,13 +77,15 @@ GalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ const Complex Bx_old = fields(i,j,k,Idx::Bx); const Complex By_old = fields(i,j,k,Idx::By); const Complex Bz_old = fields(i,j,k,Idx::Bz); - // Shortcut for the values of J and rho + + // 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); - // k vector values, and coefficients + + // k vector values const Real kx = modified_kx_arr[i]; #if (AMREX_SPACEDIM==3) const Real ky = modified_ky_arr[j]; @@ -89,8 +94,12 @@ GalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ constexpr Real ky = 0; const Real kz = modified_kz_arr[j]; #endif - constexpr Real c2 = PhysConst::c*PhysConst::c; - constexpr Complex I = Complex{0,1}; + // Physical constant c**2 and imaginary unit + constexpr Real c2 = PhysConst::c*PhysConst::c; + constexpr Complex I = Complex{0._rt,1._rt}; + + // The definition of these coefficients is explained in more detail + // in the function InitializeSpectralCoefficients below const Real C = C_arr(i,j,k); const Real S_ck = S_ck_arr(i,j,k); const Complex X1 = X1_arr(i,j,k); @@ -99,43 +108,62 @@ GalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ const Complex X4 = X4_arr(i,j,k); const Complex T2 = Theta2_arr(i,j,k); - // Update E (see the original Galilean article) - fields(i,j,k,Idx::Ex) = T2*C*Ex_old - + T2*S_ck*c2*I*(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 - + T2*S_ck*c2*I*(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 - + T2*S_ck*c2*I*(kx*By_old - ky*Bx_old) - + X4*Jz - I*(X2*rho_new - T2*X3*rho_old)*kz; - // Update B (see the original Galilean article) - // Note: here X1 is T2*x1/(ep0*c*c*k_norm*k_norm), where - // x1 has the same definition as in the original paper - fields(i,j,k,Idx::Bx) = T2*C*Bx_old - - T2*S_ck*I*(ky*Ez_old - kz*Ey_old) - + X1*I*(ky*Jz - kz*Jy); - fields(i,j,k,Idx::By) = T2*C*By_old - - T2*S_ck*I*(kz*Ex_old - kx*Ez_old) - + X1*I*(kz*Jx - kx*Jz); - fields(i,j,k,Idx::Bz) = T2*C*Bz_old - - T2*S_ck*I*(kx*Ey_old - ky*Ex_old) - + X1*I*(kx*Jy - ky*Jx); + // The equations in the following are the update equations for B and E, + // equations (11a) and (11b) of (Lehe et al, PRE 94, 2016), respectively, + // (or their rho-free formulation) + + // Update E (equation (11b) or its rho-free formulation): + if (update_with_rho) { + + // Ex + fields(i,j,k,Idx::Ex) = T2*C*Ex_old + + T2*S_ck*c2*I*(ky*Bz_old - kz*By_old) + + X4*Jx - I*(X2*rho_new - T2*X3*rho_old)*kx; + // Ey + fields(i,j,k,Idx::Ey) = T2*C*Ey_old + + T2*S_ck*c2*I*(kz*Bx_old - kx*Bz_old) + + X4*Jy - I*(X2*rho_new - T2*X3*rho_old)*ky; + // Ez + fields(i,j,k,Idx::Ez) = T2*C*Ez_old + + T2*S_ck*c2*I*(kx*By_old - ky*Bx_old) + + X4*Jz - I*(X2*rho_new - T2*X3*rho_old)*kz; + } else { + + Complex k_dot_J = kx * Jx + ky * Jy + kz * Jz; + Complex k_dot_E = kx * Ex_old + ky * Ey_old + kz * Ez_old; + + // Ex + fields(i,j,k,Idx::Ex) = T2 * C * Ex_old + I * T2 * S_ck * c2 * (ky * Bz_old - kz * By_old) + + X4 * Jx + X2 * k_dot_E * kx + X3 * k_dot_J * kx; + // Ey + fields(i,j,k,Idx::Ey) = T2 * C * Ey_old + I * T2 * S_ck * c2 * (kz * Bx_old - kx * Bz_old) + + X4 * Jy + X2 * k_dot_E * ky + X3 * k_dot_J * ky; + // Ez + fields(i,j,k,Idx::Ez) = T2 * C * Ez_old + I * T2 * S_ck * c2 * (kx * By_old - ky * Bx_old) + + X4 * Jz + X2 * k_dot_E * kz + X3 * k_dot_J * kz; + } + + // Update B (equation (11a) with X1 rescaled by theta/(epsilon_0*c**2*k**2)): + // Bx + fields(i,j,k,Idx::Bx) = T2*C*Bx_old - T2*S_ck*I*(ky*Ez_old - kz*Ey_old) + X1*I*(ky*Jz - kz*Jy); + // By + fields(i,j,k,Idx::By) = T2*C*By_old - T2*S_ck*I*(kz*Ex_old - kx*Ez_old) + X1*I*(kz*Jx - kx*Jz); + // Bz + fields(i,j,k,Idx::Bz) = T2*C*Bz_old - T2*S_ck*I*(kx*Ey_old - ky*Ex_old) + X1*I*(kx*Jy - ky*Jx); }); } }; - -void GalileanAlgorithm::InitializeSpectralCoefficients(const SpectralKSpace& spectral_kspace, - const amrex::DistributionMapping& dm, - const Array<Real, 3>& v_galilean, - const amrex::Real dt) +void GalileanAlgorithm::InitializeSpectralCoefficients (const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const amrex::Real dt) { + const bool update_with_rho = m_update_with_rho; + const BoxArray& ba = spectral_kspace.spectralspace_ba; - // Fill them with the right values: - // Loop over boxes and allocate the corresponding coefficients - // for each box owned by the local MPI proc - for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi){ + + // Loop over boxes and allocate the corresponding coefficients for each box + for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi) { const Box& bx = ba[mfi]; @@ -145,6 +173,7 @@ void GalileanAlgorithm::InitializeSpectralCoefficients(const SpectralKSpace& spe const Real* modified_ky = modified_ky_vec[mfi].dataPtr(); #endif const Real* modified_kz = modified_kz_vec[mfi].dataPtr(); + // Extract arrays for the coefficients Array4<Real> C = C_coef[mfi].array(); Array4<Real> S_ck = S_ck_coef[mfi].array(); @@ -152,17 +181,17 @@ void GalileanAlgorithm::InitializeSpectralCoefficients(const SpectralKSpace& spe Array4<Complex> X2 = X2_coef[mfi].array(); Array4<Complex> X3 = X3_coef[mfi].array(); Array4<Complex> X4 = X4_coef[mfi].array(); - Array4<Complex> Theta2 = Theta2_coef[mfi].array(); - // Extract reals (for portability on GPU) - Real vx = v_galilean[0]; + Array4<Complex> T2 = Theta2_coef[mfi].array(); + + // Extract Galilean velocity + Real vx = m_v_galilean[0]; #if (AMREX_SPACEDIM==3) - Real vy = v_galilean[1]; + Real vy = m_v_galilean[1]; #endif - Real vz = v_galilean[2]; + Real vz = m_v_galilean[2]; // Loop over indices within one box - ParallelFor(bx, - [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Calculate norm of vector const Real k_norm = std::sqrt( @@ -173,75 +202,277 @@ void GalileanAlgorithm::InitializeSpectralCoefficients(const SpectralKSpace& spe #else std::pow(modified_kz[j], 2)); #endif + // Physical constants c, c**2, and epsilon_0, and imaginary unit + constexpr Real c = PhysConst::c; + constexpr Real c2 = c*c; + constexpr Real ep0 = PhysConst::ep0; + constexpr Complex I = Complex{0._rt,1._rt}; - // Calculate coefficients - constexpr Real c = PhysConst::c; - constexpr Real ep0 = PhysConst::ep0; - const Complex I{0.,1.}; - if (k_norm != 0){ + // Auxiliary coefficients used when update_with_rho=false + const Real dt2 = dt * dt; + const Real dt3 = dt * dt2; + Complex X2_old, X3_old; - C(i,j,k) = std::cos(c*k_norm*dt); - S_ck(i,j,k) = std::sin(c*k_norm*dt)/(c*k_norm); - - // Calculate dot product with galilean velocity - const Real kv = modified_kx[i]*vx + + // Calculate dot product of k vector with Galilean velocity + const Real kv = modified_kx[i]*vx + #if (AMREX_SPACEDIM==3) - modified_ky[j]*vy + - modified_kz[k]*vz; + modified_ky[j]*vy + modified_kz[k]*vz; #else - modified_kz[j]*vz; + modified_kz[j]*vz; #endif + // The coefficients in the following refer to the ones given in equations + // (12a)-(12d) of (Lehe et al, PRE 94, 2016), used to update B and E + // (equations (11a) and (11b) of the same reference, respectively) + + if (k_norm != 0.) { + + // Auxiliary coefficients + const Real k2 = k_norm * k_norm; + const Real ck = c * k_norm; + const Real ckdt = ck * dt; + const Complex tmp1 = amrex::exp( I * ckdt); // limit of T2 for nu = 1 + const Complex tmp2 = amrex::exp(- I * ckdt); // limit of T2 for nu = -1 + + // See equation (12a) + C (i,j,k) = std::cos(ckdt); + S_ck(i,j,k) = std::sin(ckdt) / ck; + + // See equation (12b) + const Real nu = kv / ck; + const Complex theta = amrex::exp( I * 0.5_rt * kv * dt); + const Complex theta_star = amrex::exp(- I * 0.5_rt * kv * dt); - const Real nu = kv/(k_norm*c); - const Complex theta = amrex::exp( 0.5_rt*I*kv*dt ); - const Complex theta_star = amrex::exp( -0.5_rt*I*kv*dt ); - const Complex e_theta = amrex::exp( I*c*k_norm*dt ); - - Theta2(i,j,k) = theta*theta; - - if ( (nu != 1.) && (nu != 0) ) { - - // Note: the coefficients X1, X2, X3 do not correspond - // exactly to the original Galilean paper, but the - // update equation have been modified accordingly so that - // the expressions/ below (with the update equations) - // are mathematically equivalent to those of the paper. - Complex x1 = 1._rt/(1._rt-nu*nu) * - (theta_star - C(i,j,k)*theta + I*kv*S_ck(i,j,k)*theta); - // x1, above, is identical to the original paper - X1(i,j,k) = theta*x1/(ep0*c*c*k_norm*k_norm); - // The difference betwen X2 and X3 below, and those - // from the original paper is the factor ep0*k_norm*k_norm - X2(i,j,k) = (x1 - theta*(1._rt - C(i,j,k))) - /(theta_star-theta)/(ep0*k_norm*k_norm); - X3(i,j,k) = (x1 - theta_star*(1._rt - C(i,j,k))) - /(theta_star-theta)/(ep0*k_norm*k_norm); - X4(i,j,k) = I*kv*X1(i,j,k) - theta*theta*S_ck(i,j,k)/ep0; + // This is exp(i*(k \dot v_gal)*dt) + T2(i,j,k) = theta * theta; + + if ( (nu != 1.) && (nu != 0.) ) { + + // x1 is the coefficient chi_1 in equation (12c) + Complex x1 = 1._rt / (1._rt - nu*nu) + * (theta_star - C(i,j,k) * theta + I * kv * S_ck(i,j,k) * theta); + + // X1 multiplies i*(k \times J) in the update equation for B + X1(i,j,k) = theta * x1 / (ep0 * c2 * k2); + + if (update_with_rho) { + // X2 multiplies rho_new in the update equation for E + // X3 multiplies rho_old in the update equation for E + X2(i,j,k) = (x1 - theta * (1._rt - C(i,j,k))) / (theta_star - theta) / (ep0 * k2); + X3(i,j,k) = (x1 - theta_star * (1._rt - C(i,j,k))) / (theta_star - theta) / (ep0 * k2); + } else { + // X2_old is the coefficient chi_2 in equation (12d) + // X3_old is the coefficient chi_3 in equation (12d) + // X2 multiplies (k \dot E) in the update equation for E + // X3 multiplies (k \dot J) in the update equation for E + X2_old = (x1 - theta * (1._rt - C(i,j,k))) / (theta_star - theta); + X3_old = (x1 - theta_star * (1._rt - C(i,j,k))) / (theta_star - theta); + X2(i,j,k) = T2(i,j,k) * (X2_old - X3_old) / k2; + X3(i,j,k) = I * X2_old * (T2(i,j,k) - 1._rt) / (ep0 * k2 * kv); + } + + // X4 multiplies J in the update equation for E + X4(i,j,k) = I * kv * X1(i,j,k) - T2(i,j,k) * S_ck(i,j,k) / ep0; } - if ( nu == 0) { - X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0*c*c*k_norm*k_norm); - X2(i,j,k) = (1._rt - S_ck(i,j,k)/dt) / (ep0*k_norm*k_norm); - X3(i,j,k) = (C(i,j,k) - S_ck(i,j,k)/dt) / (ep0*k_norm*k_norm); - X4(i,j,k) = -S_ck(i,j,k)/ep0; + + // Limits for nu = 0 + if (nu == 0.) { + + // X1 multiplies i*(k \times J) in the update equation for B + X1(i,j,k) = (1._rt - C(i,j,k)) / (ep0 * c2 * k2); + + if (update_with_rho) { + // X2 multiplies rho_new in the update equation for E + // X3 multiplies rho_old in the update equation for E + X2(i,j,k) = (1._rt - S_ck(i,j,k) / dt) / (ep0 * k2); + X3(i,j,k) = (C(i,j,k) - S_ck(i,j,k) / dt) / (ep0 * k2); + } else { + // X2 multiplies (k \dot E) in the update equation for E + // X3 multiplies (k \dot J) in the update equation for E + X2(i,j,k) = (1._rt - C(i,j,k)) / k2; + X3(i,j,k) = (S_ck(i,j,k) / dt - 1._rt) * dt / (ep0 * k2); + } + + // Coefficient multiplying J in update equation for E + X4(i,j,k) = - S_ck(i,j,k) / ep0; } - if ( nu == 1.) { - X1(i,j,k) = (1._rt - e_theta*e_theta + 2._rt*I*c*k_norm*dt) / (4._rt*c*c*ep0*k_norm*k_norm); - X2(i,j,k) = (3._rt - 4._rt*e_theta + e_theta*e_theta + 2._rt*I*c*k_norm*dt) / (4._rt*ep0*k_norm*k_norm*(1._rt - e_theta)); - X3(i,j,k) = (3._rt - 2._rt/e_theta - 2._rt*e_theta + e_theta*e_theta - 2._rt*I*c*k_norm*dt) / (4._rt*ep0*(e_theta - 1._rt)*k_norm*k_norm); - X4(i,j,k) = I*(-1._rt + e_theta*e_theta + 2._rt*I*c*k_norm*dt) / (4._rt*ep0*c*k_norm); + + // Limits for nu = 1 + if (nu == 1.) { + + // X1 multiplies i*(k \times J) in the update equation for B + X1(i,j,k) = (1._rt - tmp1 * tmp1 + 2._rt * I * ckdt) / (4._rt * ep0 * c2 * k2); + + if (update_with_rho) { + // X2 multiplies rho_new in the update equation for E + // X3 multiplies rho_old in the update equation for E + X2(i,j,k) = (- 3._rt + 4._rt * tmp1 - tmp1 * tmp1 - 2._rt * I * ckdt) + / (4._rt * ep0 * k2 * (tmp1 - 1._rt)); + X3(i,j,k) = (3._rt - 2._rt / tmp1 - 2._rt * tmp1 + tmp1 * tmp1 - 2._rt * I * ckdt) + / (4._rt * ep0 * k2 * (tmp1 - 1._rt)); + } else { + // X2 multiplies (k \dot E) in the update equation for E + // X3 multiplies (k \dot J) in the update equation for E + X2(i,j,k) = (1._rt - C(i,j,k)) * tmp1 / k2; + X3(i,j,k) = (2._rt * ckdt - I * tmp1 * tmp1 + 4._rt * I * tmp1 - 3._rt * I) + / (4._rt * ep0 * ck * k2); + } + + // Coefficient multiplying J in update equation for E + X4(i,j,k) = (- I + I * tmp1 * tmp1 - 2._rt * ckdt) / (4._rt * ep0 * ck); } - } else { // Handle k_norm = 0, by using the analytical limit + // Limits for nu = -1 + if (nu == -1.) { + + // X1 multiplies i*(k \times J) in the update equation for B + X1(i,j,k) = (1._rt - tmp2 * tmp2 - 2._rt * I * ckdt) / (4._rt * ep0 * c2 * k2); + + if (update_with_rho) { + // X2 multiplies rho_new in the update equation for E + // X3 multiplies rho_old in the update equation for E + X2(i,j,k) = (- 4._rt + 3._rt * tmp1 + tmp2 - 2._rt * I * ckdt * tmp1) + / (4._rt * ep0 * k2 * (tmp1 - 1._rt)); + X3(i,j,k) = (2._rt - tmp2 - 3._rt * tmp1 + 2._rt * tmp1 * tmp1 - 2._rt * I * ckdt * tmp1) + / (4._rt * ep0 * k2 * (tmp1 - 1._rt)); + } else { + // X2 multiplies (k \dot E) in the update equation for E + // X3 multiplies (k \dot J) in the update equation for E + X2(i,j,k) = (1._rt - C(i,j,k)) * tmp2 / k2; + X3(i,j,k) = (2._rt * ckdt + I * tmp2 * tmp2 - 4._rt * I * tmp2 + 3._rt * I) + / (4._rt * ep0 * ck * k2); + } + + // Coefficient multiplying J in update equation for E + X4(i,j,k) = (I - I * tmp2 * tmp2 - 2._rt * ckdt) / (4._rt * ep0 * ck); + } + } + + // Limits for k = 0 + else { + + // Limits of cos(c*k*dt) and sin(c*k*dt)/(c*k) C(i,j,k) = 1._rt; S_ck(i,j,k) = dt; - X1(i,j,k) = dt*dt/(2._rt * ep0); - X2(i,j,k) = c*c*dt*dt/(6._rt * ep0); - X3(i,j,k) = - c*c*dt*dt/(3._rt * ep0); - X4(i,j,k) = -dt/ep0; - Theta2(i,j,k) = 1._rt; + + // X1 multiplies i*(k \times J) in the update equation for B + X1(i,j,k) = dt2 / (2._rt * ep0); + + if (update_with_rho) { + // X2 multiplies rho_new in the update equation for E + // X3 multiplies rho_old in the update equation for E + X2(i,j,k) = c2 * dt2 / (6._rt * ep0); + X3(i,j,k) = - c2 * dt2 / (3._rt * ep0); + } else { + // X2 multiplies (k \dot E) in the update equation for E + // X3 multiplies (k \dot J) in the update equation for E + X2(i,j,k) = c2 * dt2 * 0.5_rt; + X3(i,j,k) = - c2 * dt3 / (6._rt * ep0); + } + + // Coefficient multiplying J in update equation for E + X4(i,j,k) = -dt / ep0; + + // Limit of exp(I*(k \dot v_gal)*dt) + T2(i,j,k) = 1._rt; + } + }); + } +} + +void +GalileanAlgorithm::CurrentCorrection (SpectralFieldData& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current, + const std::unique_ptr<amrex::MultiFab>& rho) { + // Profiling + BL_PROFILE("GalileanAlgorithm::CurrentCorrection"); + + using Idx = SpectralFieldIndex; + + // Forward Fourier transform of J and rho + field_data.ForwardTransform(*current[0], Idx::Jx, 0); + field_data.ForwardTransform(*current[1], Idx::Jy, 0); + field_data.ForwardTransform(*current[2], Idx::Jz, 0); + field_data.ForwardTransform(*rho, Idx::rho_old, 0); + field_data.ForwardTransform(*rho, Idx::rho_new, 1); + + // Loop over boxes + for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){ + + const amrex::Box& bx = field_data.fields[mfi].box(); + + // Extract arrays for the fields to be updated + amrex::Array4<Complex> fields = field_data.fields[mfi].array(); + + // Extract pointers for the k vectors + const amrex::Real* const modified_kx_arr = modified_kx_vec[mfi].dataPtr(); +#if (AMREX_SPACEDIM==3) + const amrex::Real* const modified_ky_arr = modified_ky_vec[mfi].dataPtr(); +#endif + const amrex::Real* const modified_kz_arr = modified_kz_vec[mfi].dataPtr(); + + // Local copy of member variables before GPU loop + const amrex::Real dt = m_dt; + + // Galilean velocity + const amrex::Real vgx = m_v_galilean[0]; + const amrex::Real vgy = m_v_galilean[1]; + const amrex::Real vgz = m_v_galilean[2]; + + // Loop over indices within one box + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Record old values of the fields to be updated + using Idx = SpectralFieldIndex; + + // 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); + + // k vector values, and coefficients + const amrex::Real kx = modified_kx_arr[i]; +#if (AMREX_SPACEDIM==3) + const amrex::Real ky = modified_ky_arr[j]; + const amrex::Real kz = modified_kz_arr[k]; +#else + constexpr amrex::Real ky = 0._rt; + const amrex::Real kz = modified_kz_arr[j]; +#endif + constexpr Complex I = Complex{0._rt,1._rt}; + + const amrex::Real k_norm = std::sqrt(kx * kx + ky * ky + kz * kz); + + // Correct J + if (k_norm != 0._rt) + { + const Complex k_dot_J = kx * Jx + ky * Jy + kz * Jz; + const amrex::Real k_dot_vg = kx * vgx + ky * vgy + kz * vgz; + + if ( k_dot_vg != 0._rt ) { + + 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) * kx / (k_norm * k_norm); + 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) * kz / (k_norm * k_norm); + + } else { + + 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) * ky / (k_norm * k_norm); + fields(i,j,k,Idx::Jz) = Jz - (k_dot_J - I * (rho_new - rho_old) / dt) * kz / (k_norm * k_norm); + } } }); } + + // Backward Fourier transform of J + field_data.BackwardTransform(*current[0], Idx::Jx, 0); + field_data.BackwardTransform(*current[1], Idx::Jy, 0); + field_data.BackwardTransform(*current[2], Idx::Jz, 0); } void diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H index 347012d5a..348282dce 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H @@ -35,6 +35,22 @@ class PMLPsatdAlgorithm : public SpectralBaseAlgorithm } /** + * \brief Virtual function for current correction in Fourier space + * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>). + * This function overrides the virtual function \c CurrentCorrection in the + * base class \c SpectralBaseAlgorithm and cannot be overridden by further + * derived classes. + * + * \param[in,out] field_data All fields in Fourier space + * \param[in,out] current Array of unique pointers to \c MultiFab storing + * the three components of the current density + * \param[in] rho Unique pointer to \c MultiFab storing the charge density + */ + virtual void CurrentCorrection (SpectralFieldData& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current, + const std::unique_ptr<amrex::MultiFab>& rho) override final; + + /** * \brief Virtual function for Vay current deposition in Fourier space * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>). * This function overrides the virtual function \c VayDeposition in the diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp index d2f087706..4bc147cd8 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp @@ -156,6 +156,14 @@ void PMLPsatdAlgorithm::InitializeSpectralCoefficients ( }; void +PMLPsatdAlgorithm::CurrentCorrection (SpectralFieldData& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current, + const std::unique_ptr<amrex::MultiFab>& rho) +{ + amrex::Abort("Current correction not implemented for PML PSATD"); +} + +void PMLPsatdAlgorithm::VayDeposition (SpectralFieldData& field_data, std::array<std::unique_ptr<amrex::MultiFab>,3>& current) { diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H index bfa9283e6..4bfcd9fc9 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H @@ -39,19 +39,19 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm /** * \brief Virtual function for current correction in Fourier space - * (equation (19) of https://doi.org/10.1016/j.jcp.2013.03.010). + * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>). * This function overrides the virtual function \c CurrentCorrection in the - * base class \c SpectralBaseAlgorithm (qualifier \c override) and cannot be - * overridden by further derived classes (qualifier \c final). + * base class \c SpectralBaseAlgorithm and cannot be overridden by further + * derived classes. * - * \param[in,out] field_data all fields in Fourier space - * \param[in,out] current three-dimensional array of unique pointers to MultiFab - * storing the three components of the current density - * \param[in] rho unique pointer to MultiFab storing the charge density + * \param[in,out] field_data All fields in Fourier space + * \param[in,out] current Array of unique pointers to \c MultiFab storing + * the three components of the current density + * \param[in] rho Unique pointer to \c MultiFab storing the charge density */ - virtual void CurrentCorrection ( SpectralFieldData& field_data, - std::array<std::unique_ptr<amrex::MultiFab>,3>& current, - const std::unique_ptr<amrex::MultiFab>& rho ) override final; + virtual void CurrentCorrection (SpectralFieldData& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current, + const std::unique_ptr<amrex::MultiFab>& rho) override final; /** * \brief Virtual function for Vay current deposition in Fourier space diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp index 7f9fd3edb..89138641f 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp @@ -220,9 +220,9 @@ void PsatdAlgorithm::InitializeSpectralCoefficients(const SpectralKSpace& spectr } void -PsatdAlgorithm::CurrentCorrection( SpectralFieldData& field_data, +PsatdAlgorithm::CurrentCorrection (SpectralFieldData& field_data, std::array<std::unique_ptr<amrex::MultiFab>,3>& current, - const std::unique_ptr<amrex::MultiFab>& rho ) { + const std::unique_ptr<amrex::MultiFab>& rho) { // Profiling WARPX_PROFILE( "PsatdAlgorithm::CurrentCorrection" ); diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H index fc03f95f2..b4c2597a8 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H @@ -30,19 +30,19 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ /** * \brief Virtual function for current correction in Fourier space - * (equation (19) of https://doi.org/10.1016/j.jcp.2013.03.010). + * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>). * This function overrides the virtual function \c CurrentCorrection in the - * base class \c SpectralBaseAlgorithmRZ (qualifier \c override) and cannot be - * overridden by further derived classes (qualifier \c final). + * base class \c SpectralBaseAlgorithmRZ and cannot be overridden by further + * derived classes. * - * \param[in,out] field_data all fields in Fourier space - * \param[in,out] current two-dimensional array of unique pointers to MultiFab - * storing the three components of the current density - * \param[in] rho unique pointer to MultiFab storing the charge density + * \param[in,out] field_data All fields in Fourier space + * \param[in,out] current Array of unique pointers to \c MultiFab storing + * the three components of the current density + * \param[in] rho Unique pointer to \c MultiFab storing the charge density */ - virtual void CurrentCorrection ( SpectralFieldDataRZ& field_data, - std::array<std::unique_ptr<amrex::MultiFab>,3>& current, - const std::unique_ptr<amrex::MultiFab>& rho ) override final; + virtual void CurrentCorrection (SpectralFieldDataRZ& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current, + const std::unique_ptr<amrex::MultiFab>& rho) override final; /** * \brief Virtual function for Vay current deposition in Fourier space diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp index 35a478f2c..349f266d1 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp @@ -200,9 +200,9 @@ void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const void PsatdAlgorithmRZ::CurrentCorrection (SpectralFieldDataRZ& field_data, std::array<std::unique_ptr<amrex::MultiFab>,3>& current, - const std::unique_ptr<amrex::MultiFab>& rho ) + const std::unique_ptr<amrex::MultiFab>& rho) { - amrex::Abort("Current correction not implemented in RZ"); + amrex::Abort("Current correction not implemented in RZ geometry"); } void diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H index 2a34d21ab..73ad773f3 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H @@ -33,19 +33,17 @@ class SpectralBaseAlgorithm /** * \brief Virtual function for current correction in Fourier space - * (equation (19) of https://doi.org/10.1016/j.jcp.2013.03.010). - * This virtual function is not pure: it can be overridden in derived - * classes (e.g. PsatdAlgorithm, PMLPsatdAlgorithm), but a base-class - * implementation is provided (empty implementation in this case). + * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>). + * This virtual function is pure and must be defined in derived classes. * - * \param[in,out] field_data all fields in Fourier space - * \param[in,out] current three-dimensional array of unique pointers to MultiFab - * storing the three components of the current density - * \param[in] rho unique pointer to MultiFab storing the charge density + * \param[in,out] field_data All fields in Fourier space + * \param[in,out] current Array of unique pointers to \c MultiFab storing + * the three components of the current density + * \param[in] rho Unique pointer to \c MultiFab storing the charge density */ - virtual void CurrentCorrection ( SpectralFieldData& field_data, - std::array<std::unique_ptr<amrex::MultiFab>,3>& current, - const std::unique_ptr<amrex::MultiFab>& rho ) {}; + virtual void CurrentCorrection (SpectralFieldData& field_data, + std::array<std::unique_ptr<amrex::MultiFab>,3>& current, + const std::unique_ptr<amrex::MultiFab>& rho) = 0; /** * \brief Virtual function for Vay current deposition in Fourier space diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H index 833a61aec..b0a30de50 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H @@ -30,19 +30,17 @@ class SpectralBaseAlgorithmRZ /** * \brief Virtual function for current correction in Fourier space - * (equation (19) of https://doi.org/10.1016/j.jcp.2013.03.010). - * This virtual function is not pure: it can be overridden in derived - * classes (e.g. PsatdAlgorithmRZ), but a base-class - * implementation is provided (empty implementation in this case). + * (<a href="https://doi.org/10.1016/j.jcp.2013.03.010"> Vay et al, 2013</a>). + * This virtual function is pure and must be defined in derived classes. * - * \param[in,out] field_data all fields in Fourier space - * \param[in,out] current two-dimensional array of unique pointers to MultiFab - * storing the three components of the current density - * \param[in] rho unique pointer to MultiFab storing the charge density + * \param[in,out] field_data All fields in Fourier space + * \param[in,out] current Array of unique pointers to \c MultiFab storing + * the three components of the current density + * \param[in] rho Unique pointer to \c MultiFab storing the charge density */ virtual void CurrentCorrection ( SpectralFieldDataRZ& field_data, std::array<std::unique_ptr<amrex::MultiFab>,3>& current, - const std::unique_ptr<amrex::MultiFab>& rho ) {}; + const std::unique_ptr<amrex::MultiFab>& rho ) = 0; /** * \brief Compute spectral divergence of E diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index ba805c790..e7c336eb8 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -72,7 +72,7 @@ class SpectralSolver /** * \brief Public interface to call the virtual function \c CurrentCorrection, * defined in the base class SpectralBaseAlgorithm and possibly overridden - * by its derived classes (e.g. PsatdAlgorithm, PMLPsatdAlgorithm), from + * by its derived classes (e.g. PsatdAlgorithm, GalileanAlgorithm), from * objects of class SpectralSolver through the private unique pointer \c algorithm * * \param[in,out] current three-dimensional array of unique pointers to MultiFab diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp index 7146c7d4c..5d2c954ce 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -69,7 +69,7 @@ SpectralSolver::SpectralSolver( } else { algorithm = std::unique_ptr<GalileanAlgorithm>( new GalileanAlgorithm( - k_space, dm, norder_x, norder_y, norder_z, nodal, v_galilean, dt ) ); + k_space, dm, norder_x, norder_y, norder_z, nodal, v_galilean, dt, update_with_rho ) ); } } } |