Averaged Galilean PSATD (#869)

* Read Galilean velocity

* Prepare structures for Galilean solver

* Started implementing Galilean equations

* Analytical limits for X1, X2, X3, X4 coefficients added

* Slight changes added

* Added Galilean position pusher

* Scale galilean velocity

* Remove unneeded Abort

* Fix Galilean pusher

* Allocate Theta2 array

* Fix definition of coefficients

* Increase guard cells for Galilean

* Add guard cell in particle exchange

* Added modifications for PICSAR galilean branch

* Averaged coefficients added

* Type corrected

* v_gal added to warpx_current_deposition

* v_gal added to WarpXParticleContainer.H

* Bug fixed  - update particle x-position over one time step

* Fix issues with merge from dev

* Averaged fileds allocated on fine patch (Ex)

* Preparation for merging dev into galilean.

* remove TABs by hand

* Removed a tab.

* Redeclared v_galilean as a vector & related changed

* Added an automated test

* Moved v_galilean inside WarpX constructor

* Added analysis script for the automated test

* Changed name of the automated test to galilean_psatd

* Added InitializeSpectralCoefficients method

* Removed temporary comments

* Averaged fields added to FiledGather

* Added infinite order k[i]

* Setting kx_mod( kz_mod) = 0 for index = nx/2(-nz/2)

* Adding galilean shift

* Implemented galilean shift

* Changed method's name from GalileanShift to ShiftGalileanBoundary

* Added doxygen string for ShiftGalileanBoundary

* Removed never used method LowerCornerWithCentering

* Removed temporary comments

* Removed dt as a variable from DepositCharge method and its dependencies

* Changing type of v_galilean from amrex::Vector to amrex::Array

* Changed back ng_fft to be equal nox_fft

* Changed v_galilean's type from amrex::vector to amrex::array

* Removed type

* Removed temporary comments

* Added flag 'do_time_averaging' to swith from unaveraged to averaged Galilean PSATD

* Updated doxygen

* Small fix on using 'do_time_averaging' flag

* Remove some commented Print statements

* [skip ci] Further cleanup

* Fix compilation

* Guard cells update of the averaged E,B fields

* Corrected included header file accordingly

* Removed EOL

* Removed EOL

* Corrected path of the included header file

* Updated choice of the spectral solver.

* Cleanup.

* use amrex::exp instead of std::exp

* no backward FFT for avg fields if avg is off

* Need to shift avg fields in MoveWindow

* Further cleaning

* Added 2D automated test for averaged PSATD

* Added automated 2D and 3D tests for averaged PSATD

* Removed comments

* Added specifications for averaged Galilean PSATD tests.

* Bug fixed - do FillBoundary of the averaged fields only if averaged is activated

* Do shiftMF of the averaged fields only if fft_do_time_averaging=true

* Add checksum json benchmarks for averaged Galilean PSATD

* Add missing indentation

* Add missing indentation

* Updated automated analysis script

* Updated diags period

* Fixed bug: no backtransform of the averaged fields when WARPX_DIM_RZ is set

* Fixed an uninitialized variable

* Clean-up

* Changed permissions for analysis_avg_*d.py

* Compactified variables initialization via conditional assignment

* Initialized averaged E_avg, B_avg fields

* Updated automated test for 2D averaged galilean PSATD

* Updated benchmark for 2D averaged galilean PSATD

* Updated automated 2D averaged PSATD test & benchmark

* Updated automated 3D averaged PSATD test & benchmark

* Fixed typo

* Updated benchmark averaged_galilean_3d_psatd

* Initialize the averaged fields only if the averaged Galilean PSATD algorithm is enabled

* Encreased time step for automated 3D averaged Galilean PSATD (to be equal c*dt=dz while dz/dx=3)

* Updated 3D automated test: specified standard deviations `.ux_th, .uy_th, .uz_th` along each direction

* Typo from the previous commit: enabled the  averaged algorithm

* Trying out another input parameters for 3D averaged automated test with v_gal = v_plasma

* Update benchmark for previously added 3D averaged automated test

* Update 2D averaged automated test and corresponding benchmark

* Revert "Updated diags period"

This reverts commit 7334729b05e96589e020c981efdb430ca095991d.

* Cleanup: removed unwanted comment.

* Added descriptions of the 2D and 3D automated tests for the averaged Galilean PSATD.

* Updated  value calculated via standard Galilean PSATD.

* Encreased relative error tolerance for 3D automated test for the averaged Galilean PSATD.

* Removed to avoid duplication since it already specified in the corresponding 2D/3D input scripts.

* Removed unwanted empty lines

* Added spaces after function's names

* Removed unwanted empty line

* Removed

* Fixed indentation

* Cleanup: removed #include <math.h>

* Cleanup: removed empty lines and fixed indentation

* Added 'AvgGalileanAlgorithm.cpp' to 'CMakeLists.txt'

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

1 files changed, 40 insertions, 1 deletions

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));