Implement averaged algo on staggered grids & merge spectral classes (#1544)

* Refactor and clean up some spectral classes

* Abort when current correction or Vay deposition are not implemented

* Implement general equations for averaged Galilean

* Allocate averaged MultiFabs also when aux_is_nodal=1 and do_nodal=0

* Allocate +ngextra guard cells also for averaged MultiFabs

* Make alias MultiFabs for averaged aux data

* With averaging, interpolate from avg_fp (not fp) to aux

* Fix some limits of the coefficients

* Fix bug causing NaNs in spectral coefficients

* Add 2D CI test with same analysis as nodal test

* Add 3D CI test with same analysis as nodal test

* Add limit that was not covered (knorm=0 && knorm_c!=0 && nu=0)

* Allocate T2_coef only if Galilean algorithm is used

* Allocate X4_coef only if Galilean algorithm is used

* Remove extra ghost cell from 'avg_fp' MultiFabs

1 files changed, 61 insertions, 28 deletions

diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
index 4bfcd9fc9..1ca4aaf31 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
@@ -1,5 +1,4 @@
-/* Copyright 2019 Maxence Thevenet, Remi Lehe, Revathi Jambunathan, Edoardo Zoni
- *
+/* Copyright 2019
  *
  * This file is part of WarpX.
  *
@@ -10,32 +9,45 @@
 
 #include "SpectralBaseAlgorithm.H"
 
-
 #if WARPX_USE_PSATD
-
-/**
- * \brief Class that updates the field in spectral space
+/* \brief Class that updates the field in spectral space
  * and stores the coefficients of the corresponding update equation.
  */
 class PsatdAlgorithm : public SpectralBaseAlgorithm
 {
-
     public:
-        PsatdAlgorithm(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::Real dt,
-                         const bool update_with_rho);
-        // Redefine functions from base class
-        virtual void pushSpectralFields(SpectralFieldData& f) const override final;
-        virtual int getRequiredNumberOfFields() const override final {
-            return SpectralFieldIndex::n_fields;
-        }
 
-        void InitializeSpectralCoefficients(const SpectralKSpace& spectral_kspace,
-                                    const amrex::DistributionMapping& dm,
-                                    const amrex::Real dt);
+        // TODO Add Doxygen docs
+        PsatdAlgorithm (
+            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,
+            const bool update_with_rho,
+            const bool time_averaging);
+
+        // TODO Add Doxygen docs
+        virtual void pushSpectralFields (SpectralFieldData& f) const override final;
+
+        // TODO Add Doxygen docs
+        virtual int getRequiredNumberOfFields () const override final
+        {
+            if (m_time_averaging) {
+                return SpectralAvgFieldIndex::n_fields;
+            } else {
+                return SpectralFieldIndex::n_fields;
+            }
+        };
+
+        // TODO Add Doxygen docs
+        void InitializeSpectralCoefficients (
+            const SpectralKSpace& spectral_kspace,
+            const amrex::DistributionMapping& dm,
+            const amrex::Real dt);
 
         /**
          * \brief Virtual function for current correction in Fourier space
@@ -49,9 +61,10 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm
          *                           the three components of the current density
          * \param[in]     rho        Unique pointer to \c MultiFab storing the charge density
          */
-        virtual void CurrentCorrection (SpectralFieldData& field_data,
-                                        std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
-                                        const std::unique_ptr<amrex::MultiFab>& rho) override final;
+        virtual void CurrentCorrection (
+            SpectralFieldData& field_data,
+            std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
+            const std::unique_ptr<amrex::MultiFab>& rho) override final;
 
         /**
          * \brief Virtual function for Vay current deposition in Fourier space
@@ -64,14 +77,34 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm
          * \param[in,out] current    Array of unique pointers to \c MultiFab storing
          *                           the three components of the current density
          */
-        virtual void VayDeposition (SpectralFieldData& field_data,
-                                    std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final;
+        virtual void VayDeposition (
+            SpectralFieldData& field_data,
+            std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final;
 
     private:
-        SpectralRealCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef;
+
+        // These real and complex coefficients are always allocated
+        SpectralRealCoefficients C_coef, S_ck_coef;
+        SpectralComplexCoefficients T2_coef, X1_coef, X2_coef, X3_coef, X4_coef;
+
+        // These real and complex coefficients are allocated only with averaged Galilean PSATD
+        SpectralRealCoefficients C1_coef, C3_coef, S1_coef, S3_coef;
+        SpectralComplexCoefficients Psi1_coef, Psi2_coef, Psi3_coef, Psi4_coef,
+                                    A1_coef, A2_coef, Rhoold_coef, Rhonew_coef, Jcoef_coef;
+
+        // Centered modified finite-order k vectors
+        KVectorComponent modified_kx_vec_centered;
+#if (AMREX_SPACEDIM==3)
+        KVectorComponent modified_ky_vec_centered;
+#endif
+        KVectorComponent modified_kz_vec_centered;
+
+        // Other member variables
+        amrex::Array<amrex::Real,3> m_v_galilean;
         amrex::Real m_dt;
         bool m_update_with_rho;
+        bool m_time_averaging;
+        bool m_is_galilean;
 };
-
 #endif // WARPX_USE_PSATD
 #endif // WARPX_PSATD_ALGORITHM_H_