diff options
Diffstat (limited to 'Source/FieldSolver/SpectralSolver')
8 files changed, 439 insertions, 15 deletions
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H new file mode 100644 index 000000000..79536cf27 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H @@ -0,0 +1,30 @@ +#ifndef WARPX_AVG_GALILEAN_ALGORITHM_H_ +#define WARPX_AVG_GALILEAN_ALGORITHM_H_ + +#include "FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H" + +/* \brief Class that updates the field in spectral space + * and stores the coefficients of the corresponding update equation. + */ +class AvgGalileanAlgorithm : public SpectralBaseAlgorithm +{ + public: + AvgGalileanAlgorithm (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); + // Redefine update equation from base class + virtual void pushSpectralFields (SpectralFieldData& f) const override final; + virtual int getRequiredNumberOfFields () const override final { + return SpectralAvgFieldIndex::n_fields; + } + + private: + SpectralRealCoefficients C_coef, S_ck_coef, C1_coef, C3_coef, S1_coef,S3_coef; + SpectralComplexCoefficients Theta2_coef, X1_coef, X2_coef, X3_coef, X4_coef, Psi1_coef, Psi2_coef, Psi3_coef, Psi4_coef, A1_coef, A2_coef, Rhoold_coef, Rhonew_coef, Jcoef_coef; + +}; + +#endif // WARPX_GALILEAN_ALGORITHM_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp new file mode 100644 index 000000000..1369bb2e9 --- /dev/null +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.cpp @@ -0,0 +1,371 @@ +#include "FieldSolver/SpectralSolver/SpectralAlgorithms/AvgGalileanAlgorithm.H" +#include "Utils/WarpXConst.H" +#include <cmath> + +using namespace amrex; + +/* \brief Initialize coefficients for the update equation */ +AvgGalileanAlgorithm::AvgGalileanAlgorithm(const SpectralKSpace& spectral_kspace, + const 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 Real dt) + // Initialize members of base classinde + : SpectralBaseAlgorithm( spectral_kspace, dm, + norder_x, norder_y, norder_z, nodal ) +{ + const BoxArray& ba = spectral_kspace.spectralspace_ba; + + // Allocate the arrays of coefficients + C_coef = SpectralRealCoefficients(ba, dm, 1, 0); + S_ck_coef = SpectralRealCoefficients(ba, dm, 1, 0); + + C1_coef = SpectralRealCoefficients(ba, dm, 1, 0); + S1_coef = SpectralRealCoefficients(ba, dm, 1, 0); + C3_coef = SpectralRealCoefficients(ba, dm, 1, 0); + S3_coef = SpectralRealCoefficients(ba, dm, 1, 0); + + Psi1_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + Psi2_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + Psi3_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + + + X1_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + X2_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + X3_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + X4_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + Theta2_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + + A1_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + A2_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + + Rhoold_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + Rhonew_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + Jcoef_coef = SpectralComplexCoefficients(ba, dm, 1, 0); + + // Fill them with the right values: + // Loop over boxes and allocate the corresponding coefficients + // for each box owned by the local MPI proc + for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi){ + + const Box& bx = ba[mfi]; + + // Extract pointers for the k vectors + const Real* modified_kx = modified_kx_vec[mfi].dataPtr(); +#if (AMREX_SPACEDIM==3) + const Real* modified_ky = modified_ky_vec[mfi].dataPtr(); +#endif + const Real* modified_kz = modified_kz_vec[mfi].dataPtr(); + // Extract arrays for the coefficients + Array4<Real> C = C_coef[mfi].array(); + Array4<Real> S_ck = S_ck_coef[mfi].array(); + Array4<Real> C1 = C1_coef[mfi].array(); + Array4<Real> S1 = S1_coef[mfi].array(); + Array4<Real> C3 = C3_coef[mfi].array(); + Array4<Real> S3 = S3_coef[mfi].array(); + + Array4<Complex> Psi1 = Psi1_coef[mfi].array(); + Array4<Complex> Psi2 = Psi2_coef[mfi].array(); + Array4<Complex> Psi3 = Psi3_coef[mfi].array(); + Array4<Complex> X1 = X1_coef[mfi].array(); + 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(); + Array4<Complex> A1 = A1_coef[mfi].array(); + Array4<Complex> A2 = A2_coef[mfi].array(); + + Array4<Complex> CRhoold = Rhoold_coef[mfi].array(); + Array4<Complex> CRhonew = Rhonew_coef[mfi].array(); + Array4<Complex> Jcoef = Jcoef_coef[mfi].array(); + // Extract reals (for portability on GPU) + + Real vx = v_galilean[0]; + Real vy = v_galilean[1]; + Real vz = v_galilean[2]; + + // Loop over indices within one box + ParallelFor(bx, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Calculate norm of vector + const Real k_norm = std::sqrt( + std::pow(modified_kx[i], 2) + +#if (AMREX_SPACEDIM==3) + std::pow(modified_ky[j], 2) + + std::pow(modified_kz[k], 2)); +#else + std::pow(modified_kz[j], 2)); +#endif + + // Calculate coefficients + constexpr Real c = PhysConst::c; + constexpr Real c2 = PhysConst::c*PhysConst::c; + constexpr Real ep0 = PhysConst::ep0; + const Complex I{0.,1.}; + if (k_norm != 0){ + + C(i,j,k) = std::cos(c*k_norm*dt); + S_ck(i,j,k) = std::sin(c*k_norm*dt)/(c*k_norm); + + C1(i,j,k) = std::cos(0.5_rt*c*k_norm*dt); + S1(i,j,k) = std::sin(0.5_rt*c*k_norm*dt); + C3(i,j,k) = std::cos(1.5_rt*c*k_norm*dt); + S3(i,j,k) = std::sin(1.5_rt*c*k_norm*dt); + + // Calculate dot product with galilean velocity + const Real kv = modified_kx[i]*vx + +#if (AMREX_SPACEDIM==3) + modified_ky[j]*vy + + modified_kz[k]*vz; +#else + modified_kz[j]*vz; +#endif + + 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); + + Complex C_rho = I* c2 /( (1._rt-theta*theta) * ep0); + + Psi1(i,j,k) = theta * ((S1(i,j,k) + I*nu*C1(i,j,k)) + - Theta2(i,j,k) * (S3(i,j,k) + I*nu*C3(i,j,k))) /(c*k_norm*dt * (nu*nu - 1._rt)); + Psi2(i,j,k) = theta * ((C1(i,j,k) - I*nu*S1(i,j,k)) + - Theta2(i,j,k) * (C3(i,j,k) - I*nu*S3(i,j,k))) /(c2*k_norm*k_norm*dt * (nu*nu - 1._rt)); + Psi3(i,j,k) = I * theta * (1._rt - theta*theta) /(c*k_norm*dt*nu); + + A1(i,j,k) = (Psi1(i,j,k) - 1._rt + I * kv*Psi2(i,j,k) )/ (c2* k_norm*k_norm * (nu*nu - 1._rt)); + A2(i,j,k) = (Psi3(i,j,k) - Psi1(i,j,k)) / (c2*k_norm*k_norm); + + CRhoold(i,j,k) = C_rho * (theta*theta * A1(i,j,k) - A2(i,j,k)); + CRhonew(i,j,k) = C_rho * (A2(i,j,k) - A1(i,j,k)); + Jcoef(i,j,k) = (I*kv*A1(i,j,k) + Psi2(i,j,k))/ep0; + // 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; + } + 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; + + Psi1(i,j,k) = (-S1(i,j,k) + S3(i,j,k)) / (c*k_norm*dt); + Psi2(i,j,k) = (-C1(i,j,k) + C3(i,j,k)) / (c2*k_norm*k_norm*dt); + Psi3(i,j,k) = 1._rt; + A1(i,j,k) = (c*k_norm*dt + S1(i,j,k) - S3(i,j,k)) / (c*c2 * k_norm*k_norm*k_norm * dt); + A2(i,j,k) = (c*k_norm*dt + S1(i,j,k) - S3(i,j,k)) / (c*c2 * k_norm*k_norm*k_norm * dt); + CRhoold(i,j,k) = 2._rt * I * S1(i,j,k) * ( dt*C(i,j,k) - S_ck(i,j,k)) + / (c*k_norm*k_norm*k_norm*dt*dt*ep0); + CRhonew(i,j,k) = - I * (c2* k_norm*k_norm * dt*dt - C1(i,j,k) + C3(i,j,k)) + / (c2 * k_norm*k_norm*k_norm*k_norm * ep0 * dt*dt); + Jcoef(i,j,k) = (-C1(i,j,k) + C3(i,j,k)) / (c2*ep0*k_norm*k_norm*dt); + } + 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); + } + + } else { // Handle k_norm = 0, by using the analytical limit + C(i,j,k) = 1._rt; + S_ck(i,j,k) = dt; + C1(i,j,k) = 1._rt; + S1(i,j,k) = 0._rt; + C3(i,j,k) = 1._rt; + S3(i,j,k) = 0._rt; + + X1(i,j,k) = dt*dt/(2._rt * ep0); + X2(i,j,k) = c2*dt*dt/(6._rt * ep0); + X3(i,j,k) = - c2*dt*dt/(3._rt * ep0); + X4(i,j,k) = -dt/ep0; + Theta2(i,j,k) = 1._rt; + + Psi1(i,j,k) = 1._rt; + Psi2(i,j,k) = -dt; + Psi3(i,j,k) = 1._rt; + A1(i,j,k) = 13._rt * dt*dt /24._rt; + A2(i,j,k) = 13._rt * dt*dt /24._rt; + CRhoold(i,j,k) = -I*c2 * dt*dt / (3._rt * ep0); + CRhonew(i,j,k) = -5._rt*I*c2 * dt*dt / (24._rt * ep0); + Jcoef(i,j,k) = -dt/ep0; + } + + }); + } +}; + +/* Advance the E and B field in spectral space (stored in `f`) + * over one time step */ +void +AvgGalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ + + // Loop over boxes + for (MFIter mfi(f.fields); mfi.isValid(); ++mfi){ + + const Box& bx = f.fields[mfi].box(); + + // Extract arrays for the fields to be updated + Array4<Complex> fields = f.fields[mfi].array(); + // Extract arrays for the coefficients + Array4<const Real> C_arr = C_coef[mfi].array(); + Array4<const Real> S_ck_arr = S_ck_coef[mfi].array(); + Array4<const Complex> X1_arr = X1_coef[mfi].array(); + Array4<const Complex> X2_arr = X2_coef[mfi].array(); + Array4<const Complex> X3_arr = X3_coef[mfi].array(); + Array4<const Complex> X4_arr = X4_coef[mfi].array(); + Array4<const Complex> Theta2_arr = Theta2_coef[mfi].array(); + Array4<const Complex> Psi1_arr = Psi1_coef[mfi].array(); + Array4<const Complex> Psi2_arr = Psi2_coef[mfi].array(); + Array4<const Complex> Psi3_arr = Psi3_coef[mfi].array(); + Array4<const Real> C1_arr = C1_coef[mfi].array(); + Array4<const Real> S1_arr = S1_coef[mfi].array(); + Array4<const Real> C3_arr = C3_coef[mfi].array(); + Array4<const Real> S3_arr = S3_coef[mfi].array(); + + + Array4<const Complex> A1_arr = A1_coef[mfi].array(); + Array4<const Complex> A2_arr = A2_coef[mfi].array(); + Array4<const Complex> Rhonew_arr = Rhonew_coef[mfi].array(); + Array4<const Complex> Rhoold_arr = Rhoold_coef[mfi].array(); + Array4<const Complex> Jcoef_arr =Jcoef_coef[mfi].array(); + // Extract pointers for the k vectors + const Real* modified_kx_arr = modified_kx_vec[mfi].dataPtr(); +#if (AMREX_SPACEDIM==3) + const Real* modified_ky_arr = modified_ky_vec[mfi].dataPtr(); +#endif + const Real* modified_kz_arr = modified_kz_vec[mfi].dataPtr(); + + // Loop over indices within one box + ParallelFor(bx, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Record old values of the fields to be updated + using Idx = SpectralAvgFieldIndex; + + const Complex Ex_old = fields(i,j,k,Idx::Ex); + const Complex Ey_old = fields(i,j,k,Idx::Ey); + const Complex Ez_old = fields(i,j,k,Idx::Ez); + const Complex Bx_old = fields(i,j,k,Idx::Bx); + const Complex By_old = fields(i,j,k,Idx::By); + const Complex Bz_old = fields(i,j,k,Idx::Bz); + + // Shortcut for the values of J and rho + const Complex Jx = fields(i,j,k,Idx::Jx); + const Complex Jy = fields(i,j,k,Idx::Jy); + const Complex Jz = fields(i,j,k,Idx::Jz); + const Complex rho_old = fields(i,j,k,Idx::rho_old); + const Complex rho_new = fields(i,j,k,Idx::rho_new); + + const Complex Ex_avg = fields(i,j,k,Idx::Ex_avg); + const Complex Ey_avg= fields(i,j,k,Idx::Ey_avg); + const Complex Ez_avg = fields(i,j,k,Idx::Ez_avg); + const Complex Bx_avg = fields(i,j,k,Idx::Bx_avg); + const Complex By_avg = fields(i,j,k,Idx::By_avg); + const Complex Bz_avg = fields(i,j,k,Idx::Bz_avg); + // k vector values, and coefficients + const Real kx = modified_kx_arr[i]; + +#if (AMREX_SPACEDIM==3) + const Real ky = modified_ky_arr[j]; + const Real kz = modified_kz_arr[k]; +#else + constexpr Real ky = 0; + const Real kz = modified_kz_arr[j]; +#endif + constexpr Real c = PhysConst::c; + constexpr Real c2 = PhysConst::c*PhysConst::c; + constexpr Real inv_ep0 = 1._rt/PhysConst::ep0; + constexpr Complex I = Complex{0,1}; + + const Real C = C_arr(i,j,k); + const Real S_ck = S_ck_arr(i,j,k); + + const Real C1 = C1_arr(i,j,k); + const Real C3 = C3_arr(i,j,k); + const Real S1 = S1_arr(i,j,k); + const Real S3 = S3_arr(i,j,k); + + const Complex X1 = X1_arr(i,j,k); + const Complex X2 = X2_arr(i,j,k); + const Complex X3 = X3_arr(i,j,k); + const Complex X4 = X4_arr(i,j,k); + const Complex T2 = Theta2_arr(i,j,k); + + const Complex Psi1 = Psi1_arr(i,j,k); + const Complex Psi2 = Psi2_arr(i,j,k); + const Complex Psi3 = Psi3_arr(i,j,k); + const Complex A1 = A1_arr(i,j,k); + const Complex A2 = A2_arr(i,j,k); + const Complex CRhoold= Rhoold_arr(i,j,k); + const Complex CRhonew= Rhonew_arr(i,j,k); + const Complex Jcoef = Jcoef_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); + +//Update the averaged E,B fields in time on the interval [(n-1/2)dx, (n+1/2)dx] + fields(i,j,k,Idx::Ex_avg) = Psi1*Ex_old + - Psi2*c2*I*(ky*Bz_old - kz*By_old) + + Jcoef*Jx + ( CRhonew * rho_new + CRhoold*rho_old )*kx; + fields(i,j,k,Idx::Ey_avg) = Psi1*Ey_old + - Psi2*c2*I*(kz*Bx_old - kx*Bz_old) + + Jcoef*Jy +( CRhonew * rho_new + CRhoold*rho_old )*ky; + fields(i,j,k,Idx::Ez_avg) = Psi1*Ez_old + - Psi2*c2*I*(kx*By_old - ky*Bx_old) + + Jcoef*Jz + ( CRhonew * rho_new + CRhoold*rho_old )*kz; + + fields(i,j,k,Idx::Bx_avg) = Psi1*Bx_old + + I*Psi2*(ky*Ez_old - kz*Ey_old) + + A1*I*(ky*Jz - kz*Jy)*inv_ep0; + fields(i,j,k,Idx::By_avg) = Psi1*By_old + + I*Psi2*(kz*Ex_old - kx*Ez_old) + + A1*I*(kz*Jx - kx*Jz)*inv_ep0; + fields(i,j,k,Idx::Bz_avg) = Psi1*Bz_old + + I*Psi2*(kx*Ey_old - ky*Ex_old) + + A1*I*(kx*Jy - ky*Jx)*inv_ep0; + }); + } +}; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt index b80091aaf..c5ebf1c3e 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt @@ -4,6 +4,7 @@ target_sources(WarpX PMLPsatdAlgorithm.cpp PsatdAlgorithm.cpp SpectralBaseAlgorithm.cpp + AvgGalileanAlgorithm.cpp ) if(WarpX_DIMS STREQUAL RZ) diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package index 436d99adf..f316abefa 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package @@ -1,6 +1,7 @@ CEXE_sources += SpectralBaseAlgorithm.cpp CEXE_sources += PsatdAlgorithm.cpp CEXE_sources += GalileanAlgorithm.cpp +CEXE_sources += AvgGalileanAlgorithm.cpp CEXE_sources += PMLPsatdAlgorithm.cpp ifeq ($(USE_RZ),TRUE) diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index f618fda35..b445054cc 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -27,7 +27,13 @@ struct SpectralFieldIndex { enum { Ex=0, Ey, Ez, Bx, By, Bz, Jx, Jy, Jz, rho_old, rho_new, n_fields, divE=3 }; }; -/** Index for the PML fields, when stored in spectral space */ +/* Index for the regular fields + averaged fields, when stored in spectral space */ +struct SpectralAvgFieldIndex { + 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 +}; + +/* Index for the PML fields, when stored in spectral space */ struct SpectralPMLIndex { enum { Exy=0, Exz, Eyx, Eyz, Ezx, Ezy, Bxy, Bxz, Byx, Byz, Bzx, Bzy, n_fields }; diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 89771a3f7..8de6de185 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -10,7 +10,6 @@ #include <cmath> - using namespace amrex; using namespace Gpu; @@ -217,11 +216,11 @@ SpectralKSpace::getModifiedKComponent( const DistributionMapping& dm, // contains only the positive k, and the Nyquist frequency is // the last element of the array. modified_k[k.size()-1] = 0.0_rt; - } else { + } else { // The other axes contains both positive and negative k ; // the Nyquist frequency is in the middle of the array. modified_k[k.size()/2] = 0.0_rt; - } + } } } return modified_k_comp; diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index 41523d1b4..bdec4499d 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -83,7 +83,10 @@ class SpectralSolver algorithm->CurrentCorrection( field_data, current, rho ); }; + bool fft_do_time_averaging = false; + private: + void ReadParameters (); // Store field in spectral space and perform the Fourier transforms SpectralFieldData field_data; diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp index 0be623eed..16894da78 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -8,10 +8,11 @@ #include "SpectralSolver.H" #include "SpectralAlgorithms/PsatdAlgorithm.H" #include "SpectralAlgorithms/GalileanAlgorithm.H" +#include "SpectralAlgorithms/AvgGalileanAlgorithm.H" #include "SpectralAlgorithms/PMLPsatdAlgorithm.H" #include "WarpX.H" #include "Utils/WarpXProfilerWrapper.H" - +#include "Utils/WarpXUtil.H" #if WARPX_USE_PSATD @@ -50,19 +51,30 @@ SpectralSolver::SpectralSolver( // - Select the algorithm depending on the input parameters // Initialize the corresponding coefficients over k space + amrex::ParmParse pp("psatd"); + pp.query("do_time_averaging", fft_do_time_averaging); + if (pml) { algorithm = std::unique_ptr<PMLPsatdAlgorithm>( new PMLPsatdAlgorithm( k_space, dm, norder_x, norder_y, norder_z, nodal, dt ) ); - } else if ((v_galilean[0]==0) && (v_galilean[1]==0) && (v_galilean[2]==0)){ - // v_galilean is 0: use standard PSATD algorithm - algorithm = std::unique_ptr<PsatdAlgorithm>( new PsatdAlgorithm( - k_space, dm, norder_x, norder_y, norder_z, nodal, dt, update_with_rho ) ); - } else { - // Otherwise: use the Galilean algorithm - algorithm = std::unique_ptr<GalileanAlgorithm>( new GalileanAlgorithm( - k_space, dm, norder_x, norder_y, norder_z, nodal, v_galilean, dt )); - } - + } + else { + if (fft_do_time_averaging){ + algorithm = std::unique_ptr<AvgGalileanAlgorithm>( new AvgGalileanAlgorithm( + k_space, dm, norder_x, norder_y, norder_z, nodal, v_galilean, dt ) ); + } + else { + if ((v_galilean[0]==0) && (v_galilean[1]==0) && (v_galilean[2]==0)){ + // v_galilean is 0: use standard PSATD algorithm + algorithm = std::unique_ptr<PsatdAlgorithm>( new PsatdAlgorithm( + k_space, dm, norder_x, norder_y, norder_z, nodal, dt, update_with_rho ) ); + } + else { + algorithm = std::unique_ptr<GalileanAlgorithm>( new GalileanAlgorithm( + k_space, dm, norder_x, norder_y, norder_z, nodal, v_galilean, dt ) ); + } + } + } // - Initialize arrays for fields in spectral space + FFT plans field_data = SpectralFieldData( realspace_ba, k_space, dm, @@ -96,4 +108,5 @@ SpectralSolver::pushSpectralFields(){ // initialized in the constructor of `SpectralSolver` algorithm->pushSpectralFields( field_data ); } + #endif // WARPX_USE_PSATD |