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 <dpgrote@lbl.gov>

1 files changed, 76 insertions, 1 deletions

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<std::unique_ptr<amrex::MultiFab>,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<Complex> 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