Current correction in Fourier space (#675)

* Start implementing PSATD push without rho.

TODO: 1) fix unit test pml_x_psatd;
      2) try new PSATD push in PML;
      3) avoid intro of new derived class?

* Correct PSATD push to fix PML test.

* Few improvements on new PSATD push:

- new class name is 'PsatdAlgorithmMixed' (both rho and J are used);
- new algorithm parameter to choose between available implementations:
  parameter name is 'psatd_push', possible values are 'standard' for
  old implementation using rho via Gauss law and continuity equation
  and 'mixed' for new implementation using rho via Gauss law and J for
  all remaining terms.

* Fix style error (tabs vs four white spaces).

* Improve comments for available PSATD algorithms.

* Correct few typos in latest comments.

* Implement first current correction:

- new member function 'CurrentCorrection' in class SpectralSolver;
- correction applied only without subcycling (in function 'OneStep_nosub');
- TODO: add correction when subcycling is used (in function 'OneStep_sub1');
- back to old implementation of PSATD push of E and B (class 'PsatdAlgorithmMixed' removed);
- TODO: PML unit test 'pml_x_psatd' does not pass.

* Small cleanup:

- remove residual option for choice of PSATD push algorithm (only one choice now);
- improve comments.

* Implement div(E) diagnostics for spectral case.

* split travis tests in bigger matrix

* split more TravisCI tests, add electrostatic, use defaults values

* typo

* Move computation of div(E) to base class SpectralBaseAlgorithm.

* need to split psatd too

* consistent variable names and use function to avoid duplication

* fix typo for qed tests

* typo

* also need to update run_tests.sg

* Update copyright tags.

* change matrix

* Add test of div(E) vs rho/epsilon_0 in PML test.

* Small clean-up.

* Small clean-up

* Remove option for current correction from input files of two new tests (not used)

* Small clean-up: remove unnecessary references

* More clean-up (minimize style changes to keep PR simple).

* Add specific 2D/3D tests for current correction:

- 'Langmuir_multi_2d_psatd_cc': same input file as 'Langmuir_multi_2d_psatd',
  except for current correction and output of divE;
- 'Langmuir_multi_psatd_cc': same input file (3D) as 'Langmuir_multi_psatd',
  except for current correction and output of divE;
- add corresponding Python scripts for analysis: same as previous ones,
  except for check on L-infinity spatial norm of rho/epsilon_0 vs div(E);
- revert changes on old tests: do not use current correction in old tests
  (benchmarks on Battra do not need to be updated).

* Improve comments.

* Start implementation of new averaging with staggering:

- face-to-cell-center and edge-to-cell-center replaced so far;
- TODO: node-to-cell-center and 1D behavior (AMREX_SPACEDIM=1).

* Remove unrelated style changes (cleaner PR).

* Avoid duplication of input files and analysis scripts for new tests

* Improve comments for Doxygen documentation.

* Improve comments for Doxygen documentation.

* Small clean-up

* first implementation of Diags base classes

* Small clean-up

* Small clean-up

* Fix erroneous non-ASCII character

* Small clean-up

* Auxiliary function for current correction in class WarpX to keep OneStep_nosub clean

* Remove unrelated style changes (cleaner PR)

* Improve comments

* Instrument virtual function 'CurrentCorrection'

* Trying to fix build error detected by LGTM analysis only

* add example, temporarily

* Continue implementation of new averaging with staggering:

- new function takes reference to single MultiFab (no vector);
- TODO: node-to-cell-center still in progress.

* Fix small bug and clean up

* Fix bug in loop over n=0,...,ncomp-1 and clean up

* add more functions

* Add Doxygen documentation and clean up

* Small clean-up in Doxygen documentation

* Compile in single precision: add _rt suffix to avoid unnecessary conversions

* Avoid accessing staggering index directly from IntVect in innermost loops

* Implement periodic-single box option for spectral

* Fix out-of-bound in the periodic, single-box case

* Replace do-while loop with for loop (default ncomp=1)

* Remove temporary pointer and pass reference to MultiFab (instead of MultiFab*)

* Replace AMREX_LAUNCH_HOST_DEVICE_LAMBDA with ParallelFor

* cleaning and initialize output mf

* use general average routine

* move flush in new class, and implemented the Plotfile derived class

* add comments

* eol

* free memory in destructor

* typo

* typo

* no need to clear MF pointers there

* though shalt not break existing tests

* FlushRaw doesnt have to be virtual for now

* The importance of being constant

* Capability to select fields in output files

* EOL

* revert to old inputs

* const in right place

* avoid brace initializer there

* oops, fix logic error in is_in

* Use old name for output image of new 2D test

* user can choose flush interval, same behavior as plot_int

* Small clean-up in Doxygen documentation

* Add option to plot raw fields

* eol

* replace ter flush with dump to avoid confusion

* add options

* Diagnostics stores a vector of functors to compute diags on the fly

* eol

* Field gather from diags to sync particle quantities

* New diagnostics handle RZ with same behavior as old ones

* cleaning and doc

* const ref for string

* smarter for loop from Axel and typo fix from Reva

* Functors to compute some fields

* simplify code following Dave's comments

* Create subfolders and add more output options (divE etc.)

* eol

* Add documentation for periodic_single_box_fft

* For periodic, single-box, apply current correction after guard cell exchange

* rename mode_avg to convertRZmodes2cartesian

* Update CellCenterFunctor.H

* fill varnames and vector of functors at the same time

* output rho_new, not rho_old

* WarpX instance not needed here

* add const

* little bit more of reorganization

* Travis CI: force 2 MPI processes only for numprocs > 2

* Use special FFT (PR #834) and new diagnostics (PR #844) in new tests

* Improve Doxygen documentation

* Move option do_current_correction from warpx to psatd

* Fix path to output files for tests using new diagnostics

* Fix additional paths to output files for new diagnostics

* Add input paramter do_current_correction to documentation

* Fix test Langmuir_multi_psatd_hybrid: do not plot divE

* Remove input parameter amr.plot_int in tests using new diagnostics

* Trigger failing source/style checks on Travis CI

* Fix build error due to public include

* Add missing const keywords

* Change test names and corresponding analysis scripts

* Improve Python script for analysis

* Do not rename output files in old CI tests (without current correction)

* Fix output file name prefix for some tests

* Trigger Travis CI build after AMReX bug fix

* Void commit: trigger Travis CI build

* Fix some tests failing due to recent changes in master

* Use new diagnostics for particle output correctly

Co-authored-by: MaxThevenet <mthevenet@lbl.gov>

* Print tolerance and error in Python analysis

Co-authored-by: MaxThevenet <mthevenet@lbl.gov>

* Print tolerance and error in Python analysis

Co-authored-by: MaxThevenet <mthevenet@lbl.gov>

* Improve documentation

Co-authored-by: MaxThevenet <mthevenet@lbl.gov>

* Fix EOL white spaces

* Fix name of particle output variables

Co-authored-by: MaxThevenet <mthevenet@lbl.gov>
Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

1 files changed, 42 insertions, 1 deletions

diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp
index 8ca059184..7cd1abd8a 100644
--- a/Source/Evolve/WarpXEvolve.cpp
+++ b/Source/Evolve/WarpXEvolve.cpp
@@ -16,6 +16,9 @@
 #ifdef WARPX_USE_PY
 #   include "Python/WarpX_py.H"
 #endif
+#ifdef WARPX_USE_PSATD
+#include "FieldSolver/SpectralSolver/SpectralSolver.H"
+#endif
 
 #ifdef BL_USE_SENSEI_INSITU
 #   include <AMReX_AmrMeshInSituBridge.H>
@@ -341,15 +344,36 @@ WarpX::OneStep_nosub (Real cur_time)
     if (warpx_py_afterdeposition) warpx_py_afterdeposition();
 #endif
 
+#ifdef WARPX_USE_PSATD
+    // Apply current correction in Fourier space
+    // (equation (19) of https://doi.org/10.1016/j.jcp.2013.03.010)
+    if ( fft_periodic_single_box == false ) {
+        // For domain decomposition with local FFT over guard cells,
+        // apply this before `SyncCurrent`, i.e. before exchanging guard cells for J
+        if ( do_current_correction ) CurrentCorrection();
+    }
+#endif
+
 #ifdef WARPX_QED
     //Do QED processes
     mypc->doQedEvents();
 #endif
 
+    // Synchronize J and rho
     SyncCurrent();
-
     SyncRho();
 
+#ifdef WARPX_USE_PSATD
+    // Apply current correction in Fourier space
+    // (equation (19) of https://doi.org/10.1016/j.jcp.2013.03.010)
+    if ( fft_periodic_single_box == true ) {
+        // For periodic, single-box FFT (FFT without guard cells)
+        // apply this after `SyncCurrent`, i.e. after exchanging guard cells for J
+        if ( do_current_correction ) CurrentCorrection();
+    }
+#endif
+
+
     // At this point, J is up-to-date inside the domain, and E and B are
     // up-to-date including enough guard cells for first step of the field
     // solve.
@@ -785,3 +809,20 @@ WarpX::applyMirrors(Real time){
         }
     }
 }
+
+#ifdef WARPX_USE_PSATD
+void
+WarpX::CurrentCorrection ()
+{
+    for ( int lev = 0; lev <= finest_level; ++lev )
+    {
+        // Apply correction on fine patch
+        spectral_solver_fp[lev]->CurrentCorrection( current_fp[lev], rho_fp[lev] );
+        if ( spectral_solver_cp[lev] )
+        {
+            // Apply correction on coarse patch
+            spectral_solver_cp[lev]->CurrentCorrection( current_cp[lev], rho_cp[lev] );
+        }
+    }
+}
+#endif