diff options
Diffstat (limited to 'Source')
26 files changed, 728 insertions, 28 deletions
diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 4e59ee1a1..c210c6f0a 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -113,12 +113,16 @@ WarpX::Evolve (int numsteps) FillBoundaryB(guard_cells.ng_FieldGather, guard_cells.ng_Extra); // E and B: enough guard cells to update Aux or call Field Gather in fp and cp // Need to update Aux on lower levels, to interpolate to higher levels. + if (fft_do_time_averaging) + { + FillBoundaryE_avg(guard_cells.ng_FieldGather, guard_cells.ng_Extra); + FillBoundaryB_avg(guard_cells.ng_FieldGather, guard_cells.ng_Extra); + } #ifndef WARPX_USE_PSATD FillBoundaryAux(guard_cells.ng_UpdateAux); #endif UpdateAuxilaryData(); } - if (do_subcycling == 0 || finest_level == 0) { OneStep_nosub(cur_time); // E : guard cells are up-to-date @@ -593,6 +597,8 @@ WarpX::PushParticlesandDepose (int lev, amrex::Real cur_time, DtType a_dt_type) mypc->Evolve(lev, *Efield_aux[lev][0],*Efield_aux[lev][1],*Efield_aux[lev][2], *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2], + *Efield_avg_aux[lev][0],*Efield_avg_aux[lev][1],*Efield_avg_aux[lev][2], + *Bfield_avg_aux[lev][0],*Bfield_avg_aux[lev][1],*Bfield_avg_aux[lev][2], *current_fp[lev][0],*current_fp[lev][1],*current_fp[lev][2], current_buf[lev][0].get(), current_buf[lev][1].get(), current_buf[lev][2].get(), rho_fp[lev].get(), charge_buf[lev].get(), 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 diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 1925bd883..a6487cf13 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -36,10 +36,12 @@ namespace { #endif std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield, std::array<std::unique_ptr<amrex::MultiFab>,3>& Bfield, + std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield_avg, + std::array<std::unique_ptr<amrex::MultiFab>,3>& Bfield_avg, std::array<std::unique_ptr<amrex::MultiFab>,3>& current, std::unique_ptr<amrex::MultiFab>& rho ) { - using Idx = SpectralFieldIndex; + using Idx = SpectralAvgFieldIndex; // Perform forward Fourier transform #ifdef WARPX_DIM_RZ @@ -87,6 +89,18 @@ namespace { solver.BackwardTransform(*Bfield[1], Idx::By); #endif solver.BackwardTransform(*Bfield[2], Idx::Bz); + +#ifndef WARPX_DIM_RZ + if (solver.fft_do_time_averaging){ + solver.BackwardTransform(*Efield_avg[0], Idx::Ex_avg); + solver.BackwardTransform(*Efield_avg[1], Idx::Ey_avg); + solver.BackwardTransform(*Efield_avg[2], Idx::Ez_avg); + + solver.BackwardTransform(*Bfield_avg[0], Idx::Bx_avg); + solver.BackwardTransform(*Bfield_avg[1], Idx::By_avg); + solver.BackwardTransform(*Bfield_avg[2], Idx::Bz_avg); + } +#endif } } @@ -109,10 +123,10 @@ WarpX::PushPSATD (int lev, amrex::Real /* dt */) { // Update the fields on the fine and coarse patch PushPSATDSinglePatch( *spectral_solver_fp[lev], - Efield_fp[lev], Bfield_fp[lev], current_fp[lev], rho_fp[lev] ); + Efield_fp[lev], Bfield_fp[lev], Efield_avg_fp[lev], Bfield_avg_fp[lev], current_fp[lev], rho_fp[lev] ); if (spectral_solver_cp[lev]) { PushPSATDSinglePatch( *spectral_solver_cp[lev], - Efield_cp[lev], Bfield_cp[lev], current_cp[lev], rho_cp[lev] ); + Efield_cp[lev], Bfield_cp[lev], Efield_avg_cp[lev], Bfield_avg_cp[lev], current_cp[lev], rho_cp[lev] ); } } #endif diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index f9f95e15e..60121b3f4 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -247,6 +247,11 @@ WarpX::InitLevelData (int lev, Real /*time*/) if (E_ext_grid_s == "constant") pp.getarr("E_external_grid", E_external_grid); + // initialize the averaged fields only if the averaged algorithm + // is activated ('psatd.do_time_averaging=1') + ParmParse ppsatd("psatd"); + ppsatd.query("do_time_averaging", fft_do_time_averaging ); + for (int i = 0; i < 3; ++i) { current_fp[lev][i]->setVal(0.0); if (lev > 0) @@ -254,16 +259,33 @@ WarpX::InitLevelData (int lev, Real /*time*/) if (B_ext_grid_s == "constant" || B_ext_grid_s == "default") { Bfield_fp[lev][i]->setVal(B_external_grid[i]); + if (fft_do_time_averaging) { + Bfield_avg_fp[lev][i]->setVal(B_external_grid[i]); + } + if (lev > 0) { Bfield_aux[lev][i]->setVal(B_external_grid[i]); Bfield_cp[lev][i]->setVal(B_external_grid[i]); + if (fft_do_time_averaging) { + Bfield_avg_aux[lev][i]->setVal(B_external_grid[i]); + Bfield_avg_cp[lev][i]->setVal(B_external_grid[i]); + } } } if (E_ext_grid_s == "constant" || E_ext_grid_s == "default") { Efield_fp[lev][i]->setVal(E_external_grid[i]); + if (fft_do_time_averaging) { + Efield_avg_fp[lev][i]->setVal(E_external_grid[i]); + } + if (lev > 0) { Efield_aux[lev][i]->setVal(E_external_grid[i]); Efield_cp[lev][i]->setVal(E_external_grid[i]); + if (fft_do_time_averaging) { + Efield_avg_aux[lev][i]->setVal(E_external_grid[i]); + Efield_avg_cp[lev][i]->setVal(E_external_grid[i]); + } + } } } diff --git a/Source/Laser/LaserParticleContainer.H b/Source/Laser/LaserParticleContainer.H index 0ac494cc6..95dc05094 100644 --- a/Source/Laser/LaserParticleContainer.H +++ b/Source/Laser/LaserParticleContainer.H @@ -41,6 +41,8 @@ public: virtual void Evolve (int lev, const amrex::MultiFab&, const amrex::MultiFab&, const amrex::MultiFab&, const amrex::MultiFab&, const amrex::MultiFab&, const amrex::MultiFab&, + const amrex::MultiFab&, const amrex::MultiFab&, const amrex::MultiFab&, + const amrex::MultiFab&, const amrex::MultiFab&, const amrex::MultiFab&, amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, amrex::MultiFab*, amrex::MultiFab*, amrex::MultiFab*, amrex::MultiFab* rho, amrex::MultiFab* crho, diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index 6100a4408..c80f64d2e 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -381,6 +381,8 @@ void LaserParticleContainer::Evolve (int lev, const MultiFab&, const MultiFab&, const MultiFab&, const MultiFab&, const MultiFab&, const MultiFab&, + const MultiFab&, const MultiFab&, const MultiFab&, + const MultiFab&, const MultiFab&, const MultiFab&, MultiFab& jx, MultiFab& jy, MultiFab& jz, MultiFab* cjx, MultiFab* cjy, MultiFab* cjz, MultiFab* rho, MultiFab* crho, diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 95ed848f1..79bcce2f4 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -353,6 +353,25 @@ WarpX::FillBoundaryF (IntVect ng) } void +WarpX::FillBoundaryB_avg (IntVect ng, IntVect ng_extra_fine) +{ + for (int lev = 0; lev <= finest_level; ++lev) + { + FillBoundaryB_avg(lev, ng, ng_extra_fine); + } +} + +void +WarpX::FillBoundaryE_avg (IntVect ng, IntVect ng_extra_fine) +{ + for (int lev = 0; lev <= finest_level; ++lev) + { + FillBoundaryE_avg(lev, ng, ng_extra_fine); + } +} + + +void WarpX::FillBoundaryE(int lev, IntVect ng, IntVect ng_extra_fine) { FillBoundaryE(lev, PatchType::fine, ng+ng_extra_fine); @@ -475,6 +494,112 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type, IntVect ng) } void +WarpX::FillBoundaryE_avg(int lev, IntVect ng, IntVect ng_extra_fine) +{ + FillBoundaryE_avg(lev, PatchType::fine, ng+ng_extra_fine); + if (lev > 0) FillBoundaryE_avg(lev, PatchType::coarse, ng); +} + +void +WarpX::FillBoundaryE_avg (int lev, PatchType patch_type, IntVect ng) +{ + if (patch_type == PatchType::fine) + { + if (do_pml && pml[lev]->ok()) + { + amrex::Abort("Averaged Galilean PSATD with PML is not yet implemented"); + } + + const auto& period = Geom(lev).periodicity(); + if ( safe_guard_cells ){ + Vector<MultiFab*> mf{Efield_avg_fp[lev][0].get(),Efield_avg_fp[lev][1].get(),Efield_avg_fp[lev][2].get()}; + amrex::FillBoundary(mf, period); + } else { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + ng <= Efield_avg_fp[lev][0]->nGrowVect(), + "Error: in FillBoundaryE_avg, requested more guard cells than allocated"); + Efield_avg_fp[lev][0]->FillBoundary(ng, period); + Efield_avg_fp[lev][1]->FillBoundary(ng, period); + Efield_avg_fp[lev][2]->FillBoundary(ng, period); + } + } + else if (patch_type == PatchType::coarse) + { + if (do_pml && pml[lev]->ok()) + { + amrex::Abort("Averaged Galilean PSATD with PML is not yet implemented"); + } + + const auto& cperiod = Geom(lev-1).periodicity(); + if ( safe_guard_cells ) { + Vector<MultiFab*> mf{Efield_avg_cp[lev][0].get(),Efield_avg_cp[lev][1].get(),Efield_avg_cp[lev][2].get()}; + amrex::FillBoundary(mf, cperiod); + + } else { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + ng <= Efield_avg_cp[lev][0]->nGrowVect(), + "Error: in FillBoundaryE, requested more guard cells than allocated"); + Efield_avg_cp[lev][0]->FillBoundary(ng, cperiod); + Efield_avg_cp[lev][1]->FillBoundary(ng, cperiod); + Efield_avg_cp[lev][2]->FillBoundary(ng, cperiod); + } + } +} + + +void +WarpX::FillBoundaryB_avg (int lev, IntVect ng, IntVect ng_extra_fine) +{ + FillBoundaryB_avg(lev, PatchType::fine, ng + ng_extra_fine); + if (lev > 0) FillBoundaryB_avg(lev, PatchType::coarse, ng); +} + +void +WarpX::FillBoundaryB_avg (int lev, PatchType patch_type, IntVect ng) +{ + if (patch_type == PatchType::fine) + { + if (do_pml && pml[lev]->ok()) + { + amrex::Abort("Averaged Galilean PSATD with PML is not yet implemented"); + } + const auto& period = Geom(lev).periodicity(); + if ( safe_guard_cells ) { + Vector<MultiFab*> mf{Bfield_avg_fp[lev][0].get(),Bfield_avg_fp[lev][1].get(),Bfield_avg_fp[lev][2].get()}; + amrex::FillBoundary(mf, period); + } else { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + ng <= Bfield_fp[lev][0]->nGrowVect(), + "Error: in FillBoundaryB, requested more guard cells than allocated"); + Bfield_avg_fp[lev][0]->FillBoundary(ng, period); + Bfield_avg_fp[lev][1]->FillBoundary(ng, period); + Bfield_avg_fp[lev][2]->FillBoundary(ng, period); + } + } + else if (patch_type == PatchType::coarse) + { + if (do_pml && pml[lev]->ok()) + { + amrex::Abort("Averaged Galilean PSATD with PML is not yet implemented"); + } + + const auto& cperiod = Geom(lev-1).periodicity(); + if ( safe_guard_cells ){ + Vector<MultiFab*> mf{Bfield_avg_cp[lev][0].get(),Bfield_avg_cp[lev][1].get(),Bfield_avg_cp[lev][2].get()}; + amrex::FillBoundary(mf, cperiod); + } else { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + ng <= Bfield_avg_cp[lev][0]->nGrowVect(), + "Error: in FillBoundaryB_avg, requested more guard cells than allocated"); + Bfield_avg_cp[lev][0]->FillBoundary(ng, cperiod); + Bfield_avg_cp[lev][1]->FillBoundary(ng, cperiod); + Bfield_avg_cp[lev][2]->FillBoundary(ng, cperiod); + } + } +} + + +void WarpX::FillBoundaryF (int lev, IntVect ng) { FillBoundaryF(lev, PatchType::fine, ng); diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index ecd6717c0..20af8126d 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -86,6 +86,8 @@ public: void Evolve (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, + const amrex::MultiFab& Ex_avg, const amrex::MultiFab& Ey_avg, const amrex::MultiFab& Ez_avg, + const amrex::MultiFab& Bx_avg, const amrex::MultiFab& By_avg, const amrex::MultiFab& Bz_avg, amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, amrex::MultiFab* cjx, amrex::MultiFab* cjy, amrex::MultiFab* cjz, amrex::MultiFab* rho, amrex::MultiFab* crho, diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index fd0980d45..9f015124b 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -316,6 +316,8 @@ void MultiParticleContainer::Evolve (int lev, const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz, + const MultiFab& Ex_avg, const MultiFab& Ey_avg, const MultiFab& Ez_avg, + const MultiFab& Bx_avg, const MultiFab& By_avg, const MultiFab& Bz_avg, MultiFab& jx, MultiFab& jy, MultiFab& jz, MultiFab* cjx, MultiFab* cjy, MultiFab* cjz, MultiFab* rho, MultiFab* crho, @@ -332,7 +334,7 @@ MultiParticleContainer::Evolve (int lev, if (rho) rho->setVal(0.0); if (crho) crho->setVal(0.0); for (auto& pc : allcontainers) { - pc->Evolve(lev, Ex, Ey, Ez, Bx, By, Bz, jx, jy, jz, cjx, cjy, cjz, + pc->Evolve(lev, Ex, Ey, Ez, Bx, By, Bz, Ex_avg, Ey_avg, Ez_avg, Bx_avg, By_avg, Bz_avg, jx, jy, jz, cjx, cjy, cjz, rho, crho, cEx, cEy, cEz, cBx, cBy, cBz, t, dt, a_dt_type); } } diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H index ad89ca0a5..0bfe38f42 100644 --- a/Source/Particles/PhotonParticleContainer.H +++ b/Source/Particles/PhotonParticleContainer.H @@ -36,6 +36,12 @@ public: const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, + const amrex::MultiFab& Ex_avg, + const amrex::MultiFab& Ey_avg, + const amrex::MultiFab& Ez_avg, + const amrex::MultiFab& Bx_avg, + const amrex::MultiFab& By_avg, + const amrex::MultiFab& Bz_avg, amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index abf56dd7c..ff1107eea 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -196,6 +196,8 @@ void PhotonParticleContainer::Evolve (int lev, const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz, + const MultiFab& Ex_avg, const MultiFab& Ey_avg, const MultiFab& Ez_avg, + const MultiFab& Bx_avg, const MultiFab& By_avg, const MultiFab& Bz_avg, MultiFab& jx, MultiFab& jy, MultiFab& jz, MultiFab* cjx, MultiFab* cjy, MultiFab* cjz, MultiFab* rho, MultiFab* crho, @@ -209,6 +211,8 @@ PhotonParticleContainer::Evolve (int lev, PhysicalParticleContainer::Evolve (lev, Ex, Ey, Ez, Bx, By, Bz, + Ex_avg, Ey_avg, Ez_avg, + Bx_avg, By_avg, Bz_avg, jx, jy, jz, cjx, cjy, cjz, rho, crho, diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 0427807b8..c259e1e12 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -104,6 +104,12 @@ public: const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, + const amrex::MultiFab& Ex_avg, + const amrex::MultiFab& Ey_avg, + const amrex::MultiFab& Ez_avg, + const amrex::MultiFab& Bx_avg, + const amrex::MultiFab& By_avg, + const amrex::MultiFab& Bz_avg, amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 94766d48a..4d8dc67cf 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -919,6 +919,8 @@ void PhysicalParticleContainer::Evolve (int lev, const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz, + const MultiFab& Ex_avg, const MultiFab& Ey_avg, const MultiFab& Ez_avg, + const MultiFab& Bx_avg, const MultiFab& By_avg, const MultiFab& Bz_avg, MultiFab& jx, MultiFab& jy, MultiFab& jz, MultiFab* cjx, MultiFab* cjy, MultiFab* cjz, MultiFab* rho, MultiFab* crho, @@ -926,6 +928,11 @@ PhysicalParticleContainer::Evolve (int lev, const MultiFab* cBx, const MultiFab* cBy, const MultiFab* cBz, Real /*t*/, Real dt, DtType a_dt_type) { + + bool fft_do_time_averaging = false; + ParmParse pp("psatd"); + pp.query("do_time_averaging", fft_do_time_averaging); + WARPX_PROFILE("PPC::Evolve()"); WARPX_PROFILE_VAR_NS("PPC::GatherAndPush", blp_fg); @@ -984,12 +991,12 @@ PhysicalParticleContainer::Evolve (int lev, const long np = pti.numParticles(); // Data on the grid - FArrayBox const* exfab = &(Ex[pti]); - FArrayBox const* eyfab = &(Ey[pti]); - FArrayBox const* ezfab = &(Ez[pti]); - FArrayBox const* bxfab = &(Bx[pti]); - FArrayBox const* byfab = &(By[pti]); - FArrayBox const* bzfab = &(Bz[pti]); + FArrayBox const* exfab = fft_do_time_averaging ? &(Ex_avg[pti]) : &(Ex[pti]); + FArrayBox const* eyfab = fft_do_time_averaging ? &(Ey_avg[pti]) : &(Ey[pti]); + FArrayBox const* ezfab = fft_do_time_averaging ? &(Ez_avg[pti]) : &(Ez[pti]); + FArrayBox const* bxfab = fft_do_time_averaging ? &(Bx_avg[pti]) : &(Bx[pti]); + FArrayBox const* byfab = fft_do_time_averaging ? &(By_avg[pti]) : &(By[pti]); + FArrayBox const* bzfab = fft_do_time_averaging ? &(Bz_avg[pti]) : &(Bz[pti]); Elixir exeli, eyeli, ezeli, bxeli, byeli, bzeli; diff --git a/Source/Particles/RigidInjectedParticleContainer.H b/Source/Particles/RigidInjectedParticleContainer.H index 9d4257702..068b67bed 100644 --- a/Source/Particles/RigidInjectedParticleContainer.H +++ b/Source/Particles/RigidInjectedParticleContainer.H @@ -50,6 +50,12 @@ public: const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, + const amrex::MultiFab& Ex_avg, + const amrex::MultiFab& Ey_avg, + const amrex::MultiFab& Ez_avg, + const amrex::MultiFab& Bx_avg, + const amrex::MultiFab& By_avg, + const amrex::MultiFab& Bz_avg, amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 6d5565f2e..110d304bb 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -340,6 +340,8 @@ void RigidInjectedParticleContainer::Evolve (int lev, const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz, + const MultiFab& Ex_avg, const MultiFab& Ey_avg, const MultiFab& Ez_avg, + const MultiFab& Bx_avg, const MultiFab& By_avg, const MultiFab& Bz_avg, MultiFab& jx, MultiFab& jy, MultiFab& jz, MultiFab* cjx, MultiFab* cjy, MultiFab* cjz, MultiFab* rho, MultiFab* crho, @@ -367,6 +369,8 @@ RigidInjectedParticleContainer::Evolve (int lev, PhysicalParticleContainer::Evolve (lev, Ex, Ey, Ez, Bx, By, Bz, + Ex_avg, Ey_avg, Ez_avg, + Bx_avg, By_avg, Bz_avg, jx, jy, jz, cjx, cjy, cjz, rho, crho, diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 33769c61c..4c493b21f 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -157,6 +157,8 @@ public: virtual void Evolve (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, + const amrex::MultiFab& Ex_avg, const amrex::MultiFab& Ey_avg, const amrex::MultiFab& Ez_avg, + const amrex::MultiFab& Bx_avg, const amrex::MultiFab& By_avg, const amrex::MultiFab& Bz_avg, amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, amrex::MultiFab* cjx, amrex::MultiFab* cjy, amrex::MultiFab* cjz, amrex::MultiFab* rho, amrex::MultiFab* crho, diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp index 9f96512d9..8435864b8 100644 --- a/Source/Utils/WarpXMovingWindow.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp @@ -130,6 +130,10 @@ WarpX::MoveWindow (bool move_j) } shiftMF(*Bfield_fp[lev][dim], geom[lev], num_shift, dir, ng_extra, B_external_grid[dim], use_Bparser, Bfield_parser); shiftMF(*Efield_fp[lev][dim], geom[lev], num_shift, dir, ng_extra, E_external_grid[dim], use_Eparser, Efield_parser); + if (fft_do_time_averaging) { + shiftMF(*Bfield_avg_fp[lev][dim], geom[lev], num_shift, dir, ng_extra, B_external_grid[dim], use_Bparser, Bfield_parser); + shiftMF(*Efield_avg_fp[lev][dim], geom[lev], num_shift, dir, ng_extra, E_external_grid[dim], use_Eparser, Efield_parser); + } if (move_j) { shiftMF(*current_fp[lev][dim], geom[lev], num_shift, dir, ng_zero); } @@ -146,6 +150,12 @@ WarpX::MoveWindow (bool move_j) shiftMF(*Efield_cp[lev][dim], geom[lev-1], num_shift_crse, dir, ng_zero, E_external_grid[dim], use_Eparser, Efield_parser); shiftMF(*Bfield_aux[lev][dim], geom[lev], num_shift, dir, ng_zero); shiftMF(*Efield_aux[lev][dim], geom[lev], num_shift, dir, ng_zero); + if (fft_do_time_averaging) { + shiftMF(*Bfield_avg_cp[lev][dim], geom[lev-1], num_shift_crse, dir, ng_zero, B_external_grid[dim], use_Bparser, Bfield_parser); + shiftMF(*Efield_avg_cp[lev][dim], geom[lev-1], num_shift_crse, dir, ng_zero, E_external_grid[dim], use_Eparser, Efield_parser); + shiftMF(*Bfield_avg_aux[lev][dim], geom[lev], num_shift, dir, ng_zero); + shiftMF(*Efield_avg_aux[lev][dim], geom[lev], num_shift, dir, ng_zero); + } if (move_j) { shiftMF(*current_cp[lev][dim], geom[lev-1], num_shift_crse, dir, ng_zero); } diff --git a/Source/WarpX.H b/Source/WarpX.H index 401ed5b40..be02066e7 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -84,6 +84,7 @@ public: int Verbose () const { return verbose; } + void InitData (); void Evolve (int numsteps = -1); @@ -267,6 +268,7 @@ public: amrex::Real time_of_last_gal_shift = 0; amrex::Array<amrex::Real,3> v_galilean = {{0}}; + static int num_mirrors; amrex::Vector<amrex::Real> mirror_z; amrex::Vector<amrex::Real> mirror_z_width; @@ -385,10 +387,16 @@ public: // Fill boundary cells including coarse/fine boundaries void FillBoundaryB (amrex::IntVect ng, amrex::IntVect ng_extra_fine=amrex::IntVect::TheZeroVector()); void FillBoundaryE (amrex::IntVect ng, amrex::IntVect ng_extra_fine=amrex::IntVect::TheZeroVector()); + void FillBoundaryB_avg (amrex::IntVect ng, amrex::IntVect ng_extra_fine=amrex::IntVect::TheZeroVector()); + void FillBoundaryE_avg (amrex::IntVect ng, amrex::IntVect ng_extra_fine=amrex::IntVect::TheZeroVector()); + void FillBoundaryF (amrex::IntVect ng); void FillBoundaryAux (amrex::IntVect ng); void FillBoundaryE (int lev, amrex::IntVect ng, amrex::IntVect ng_extra_fine=amrex::IntVect::TheZeroVector()); void FillBoundaryB (int lev, amrex::IntVect ng, amrex::IntVect ng_extra_fine=amrex::IntVect::TheZeroVector()); + void FillBoundaryE_avg (int lev, amrex::IntVect ng, amrex::IntVect ng_extra_fine=amrex::IntVect::TheZeroVector()); + void FillBoundaryB_avg (int lev, amrex::IntVect ng, amrex::IntVect ng_extra_fine=amrex::IntVect::TheZeroVector()); + void FillBoundaryF (int lev, amrex::IntVect ng); void FillBoundaryAux (int lev, amrex::IntVect ng); @@ -575,6 +583,9 @@ private: void FillBoundaryE (int lev, PatchType patch_type, amrex::IntVect ng); void FillBoundaryF (int lev, PatchType patch_type, amrex::IntVect ng); + void FillBoundaryB_avg (int lev, PatchType patch_type, amrex::IntVect ng); + void FillBoundaryE_avg (int lev, PatchType patch_type, amrex::IntVect ng); + void OneStep_nosub (amrex::Real t); void OneStep_sub1 (amrex::Real t); @@ -664,6 +675,8 @@ private: // Full solution amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_aux; amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_aux; + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_avg_aux; + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_avg_aux; // Fine patch amrex::Vector< std::unique_ptr<amrex::MultiFab> > F_fp; @@ -671,7 +684,8 @@ private: amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_fp; amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_fp; amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_fp; - + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_avg_fp; + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_avg_fp; // store fine patch amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_store; @@ -681,6 +695,8 @@ private: amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_cp; amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_cp; amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_cp; + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_avg_cp; + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_avg_cp; // Copy of the coarse aux amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3 > > Efield_cax; @@ -797,6 +813,7 @@ private: amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_slice; amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_slice; + bool fft_do_time_averaging = false; bool fft_periodic_single_box = false; int nox_fft = 16; int noy_fft = 16; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 22debe82e..8c1dded61 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -142,7 +142,6 @@ WarpX::ResetInstance () WarpX::WarpX () { m_instance = this; - ReadParameters(); BackwardCompatibility(); @@ -190,11 +189,16 @@ WarpX::WarpX () Efield_aux.resize(nlevs_max); Bfield_aux.resize(nlevs_max); + Efield_avg_aux.resize(nlevs_max); + Bfield_avg_aux.resize(nlevs_max); + F_fp.resize(nlevs_max); rho_fp.resize(nlevs_max); current_fp.resize(nlevs_max); Efield_fp.resize(nlevs_max); Bfield_fp.resize(nlevs_max); + Efield_avg_fp.resize(nlevs_max); + Bfield_avg_fp.resize(nlevs_max); current_store.resize(nlevs_max); @@ -203,6 +207,8 @@ WarpX::WarpX () current_cp.resize(nlevs_max); Efield_cp.resize(nlevs_max); Bfield_cp.resize(nlevs_max); + Efield_avg_cp.resize(nlevs_max); + Bfield_avg_cp.resize(nlevs_max); Efield_cax.resize(nlevs_max); Bfield_cax.resize(nlevs_max); @@ -608,6 +614,7 @@ WarpX::ReadParameters () pp.query("current_correction", current_correction); pp.query("update_with_rho", update_with_rho); pp.query("v_galilean", v_galilean); + pp.query("do_time_averaging", fft_do_time_averaging); // Scale the velocity by the speed of light for (int i=0; i<3; i++) v_galilean[i] *= PhysConst::c; } @@ -860,6 +867,15 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm current_fp[lev][1].reset( new MultiFab(amrex::convert(ba,jy_nodal_flag),dm,ncomps,ngJ)); current_fp[lev][2].reset( new MultiFab(amrex::convert(ba,jz_nodal_flag),dm,ncomps,ngJ)); + + Bfield_avg_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Bx_nodal_flag),dm,ncomps,ngE)); + Bfield_avg_fp[lev][1].reset( new MultiFab(amrex::convert(ba,By_nodal_flag),dm,ncomps,ngE)); + Bfield_avg_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Bz_nodal_flag),dm,ncomps,ngE)); + + Efield_avg_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Ex_nodal_flag),dm,ncomps,ngE)); + Efield_avg_fp[lev][1].reset( new MultiFab(amrex::convert(ba,Ey_nodal_flag),dm,ncomps,ngE)); + Efield_avg_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Ez_nodal_flag),dm,ncomps,ngE)); + if (do_dive_cleaning || (plot_rho && do_back_transformed_diagnostics)) { rho_fp[lev].reset(new MultiFab(amrex::convert(ba,rho_nodal_flag),dm,2*ncomps,ngRho)); @@ -932,6 +948,9 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm for (int idir = 0; idir < 3; ++idir) { Efield_aux[lev][idir].reset(new MultiFab(*Efield_fp[lev][idir], amrex::make_alias, 0, ncomps)); Bfield_aux[lev][idir].reset(new MultiFab(*Bfield_fp[lev][idir], amrex::make_alias, 0, ncomps)); + + Efield_avg_aux[lev][idir].reset(new MultiFab(*Efield_avg_fp[lev][idir], amrex::make_alias, 0, ncomps)); + Bfield_avg_aux[lev][idir].reset(new MultiFab(*Bfield_avg_fp[lev][idir], amrex::make_alias, 0, ncomps)); } } else @@ -943,6 +962,16 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm Efield_aux[lev][0].reset( new MultiFab(amrex::convert(ba,Ex_nodal_flag),dm,ncomps,ngE)); Efield_aux[lev][1].reset( new MultiFab(amrex::convert(ba,Ey_nodal_flag),dm,ncomps,ngE)); Efield_aux[lev][2].reset( new MultiFab(amrex::convert(ba,Ez_nodal_flag),dm,ncomps,ngE)); + + + Bfield_avg_aux[lev][0].reset( new MultiFab(amrex::convert(ba,Bx_nodal_flag),dm,ncomps,ngE)); + Bfield_avg_aux[lev][1].reset( new MultiFab(amrex::convert(ba,By_nodal_flag),dm,ncomps,ngE)); + Bfield_avg_aux[lev][2].reset( new MultiFab(amrex::convert(ba,Bz_nodal_flag),dm,ncomps,ngE)); + + Efield_avg_aux[lev][0].reset( new MultiFab(amrex::convert(ba,Ex_nodal_flag),dm,ncomps,ngE)); + Efield_avg_aux[lev][1].reset( new MultiFab(amrex::convert(ba,Ey_nodal_flag),dm,ncomps,ngE)); + Efield_avg_aux[lev][2].reset( new MultiFab(amrex::convert(ba,Ez_nodal_flag),dm,ncomps,ngE)); + } // @@ -964,6 +993,16 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm Efield_cp[lev][1].reset( new MultiFab(amrex::convert(cba,Ey_nodal_flag),dm,ncomps,ngE)); Efield_cp[lev][2].reset( new MultiFab(amrex::convert(cba,Ez_nodal_flag),dm,ncomps,ngE)); + // Create the MultiFabs for B_avg + Bfield_avg_cp[lev][0].reset( new MultiFab(amrex::convert(cba,Bx_nodal_flag),dm,ncomps,ngE)); + Bfield_avg_cp[lev][1].reset( new MultiFab(amrex::convert(cba,By_nodal_flag),dm,ncomps,ngE)); + Bfield_avg_cp[lev][2].reset( new MultiFab(amrex::convert(cba,Bz_nodal_flag),dm,ncomps,ngE)); + + // Create the MultiFabs for E_avg + Efield_avg_cp[lev][0].reset( new MultiFab(amrex::convert(cba,Ex_nodal_flag),dm,ncomps,ngE)); + Efield_avg_cp[lev][1].reset( new MultiFab(amrex::convert(cba,Ey_nodal_flag),dm,ncomps,ngE)); + Efield_avg_cp[lev][2].reset( new MultiFab(amrex::convert(cba,Ez_nodal_flag),dm,ncomps,ngE)); + // Create the MultiFabs for the current current_cp[lev][0].reset( new MultiFab(amrex::convert(cba,jx_nodal_flag),dm,ncomps,ngJ)); current_cp[lev][1].reset( new MultiFab(amrex::convert(cba,jy_nodal_flag),dm,ncomps,ngJ)); |