From 4dcee1f248e18c5a16da5beec287c119258d416e Mon Sep 17 00:00:00 2001 From: Edoardo Zoni <59625522+EZoni@users.noreply.github.com> Date: Mon, 20 Jul 2020 08:33:30 -0700 Subject: Vay current deposition (#1051) * Added stub for current correction in RZ spectral solver * Start implementation of Vay deposition * Continue implementation of Vay deposition * Correct deposition of D * Add phase shift for staggered currents * Small clean-up * Fix units in deposition of D * Implement average of cumulative sum (needs bug fix) * Start fixing bug in average of cumulative sum * Still debugging * Cumulative sums should be correct now * Subtract averages of cumulative sums: - current implementation: cumulative sums, inverse Fourier transform, subtraction of averages - needs to be tested (including units of D after Vay deposition) - needs to be shortened (too many loops over boxes and ParallelFors) * [skip CI] Clean up and fix units * Still fixing units * [skip CI] Remove temporarily averages of cumulative sums * [skip CI] Remove distinction between staggered and nodal * Vay and Esirkepov similar results on periodic single box: TODO: - debug (charge not conserved); - try using compute_shifted_shape_factor as in Esirkepov deposition; - clean up; - try on multiple boxes and with correction of mode at 0 frequency. * [skip CI] Clean up * Fix bug in 3D deposition * [skip CI] Clean up * Fix 2D and 3D implementation: - simulation results agree between direct and Vay deposition in both 2D and 3D - Travis CI tests should pass except for check of charge conservation (debug) * Small clean-up * Fix bug when compiling in RZ geometry * Add benchmark json files (will be reset later) * Do not set zero current at zero frequency * [skip CI] Revert last commit and clean up * Fix small bug after reverting commit * Set nodal test first on Travis * Fix benchmark for nodal test in 3D * Fix particle output for nodal test in 3D * Fix bugs due to staggering * Rename current nodal Travis tests * Add Travis tests staggered in 2D and 3D * Further clean-up after bug fix * Abort when using Vay deposition with domain decomposition * Add optional argument of index type to forward FFT * Fourier shifts can be private members as before * Small clean-up * Clean up and improve Doxygen documentation * Fix small bug in analysis script for 2D tests * Fix tests (remove E and B fields from particle diags) * Add option to fill guard cells and docs * Fix value of last guard cell by enforcing periodicity * Revert changes merged from #1121 * Clean up style * Improve docs * Fix forgotten alignment * Improve docs * Make base class functions VayDeposition pure Co-authored-by: Dave Grote --- .../SpectralAlgorithms/PsatdAlgorithm.cpp | 77 +++++++++++++++++++++- 1 file changed, 76 insertions(+), 1 deletion(-) (limited to 'Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp') diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp index b8b64b80c..7f9fd3edb 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp @@ -256,7 +256,6 @@ PsatdAlgorithm::CurrentCorrection( SpectralFieldData& field_data, // 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 @@ -297,4 +296,80 @@ PsatdAlgorithm::CurrentCorrection( SpectralFieldData& field_data, field_data.BackwardTransform( *current[1], Idx::Jy, 0 ); field_data.BackwardTransform( *current[2], Idx::Jz, 0 ); } + +void +PsatdAlgorithm::VayDeposition (SpectralFieldData& field_data, + std::array,3>& current) { + // Profiling + WARPX_PROFILE("PsatdAlgorithm::VayDeposition"); + + using Idx = SpectralFieldIndex; + + // Forward Fourier transform of D (temporarily stored in current): + // D is nodal and does not match the staggering of J, therefore we pass the + // actual staggering of D (IntVect(1)) to the ForwardTransform function + field_data.ForwardTransform(*current[0], Idx::Jx, 0, IntVect(1)); + field_data.ForwardTransform(*current[1], Idx::Jy, 0, IntVect(1)); + field_data.ForwardTransform(*current[2], Idx::Jz, 0, IntVect(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 fields = field_data.fields[mfi].array(); + + // Extract pointers for the modified 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(); + + // Loop over indices within one box + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + using Idx = SpectralFieldIndex; + + // Shortcuts for the values of D + const Complex Dx = fields(i,j,k,Idx::Jx); + const Complex Dy = fields(i,j,k,Idx::Jy); + const Complex Dz = fields(i,j,k,Idx::Jz); + + // Imaginary unit + constexpr Complex I = Complex{0._rt, 1._rt}; + + // Modified k vector values + const amrex::Real kx_mod = modified_kx_arr[i]; +#if (AMREX_SPACEDIM==3) + const amrex::Real ky_mod = modified_ky_arr[j]; + const amrex::Real kz_mod = modified_kz_arr[k]; +#else + constexpr amrex::Real ky_mod = 0._rt; + const amrex::Real kz_mod = modified_kz_arr[j]; +#endif + + // Compute Jx + if (kx_mod != 0._rt) fields(i,j,k,Idx::Jx) = I * Dx / kx_mod; + else fields(i,j,k,Idx::Jx) = 0._rt; + +#if (AMREX_SPACEDIM==3) + // Compute Jy + if (ky_mod != 0._rt) fields(i,j,k,Idx::Jy) = I * Dy / ky_mod; + else fields(i,j,k,Idx::Jy) = 0._rt; +#endif + + // Compute Jz + if (kz_mod != 0._rt) fields(i,j,k,Idx::Jz) = I * Dz / kz_mod; + else fields(i,j,k,Idx::Jz) = 0._rt; + + }); + } + + // 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); +} #endif // WARPX_USE_PSATD -- cgit v1.2.3