Galilean PSATD with shift (#704)

* 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

* 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

* Preparation for merging dev into galilean.

* 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

* Converted tab to spaces

* Removed EOL white space

* Add documentation and automated tests

* Fix compilation error

* Add automated test

* Update automated test

* Removed temporary used galilean shift

* Removed temporary used particle's push for Galilean PSATD

* Removed unused statement

* Remove EOL white space.

* Added zero shift for LowerCorner in RZ geometry

* Minor changes to Galilean implementation

* Modifications for GPU

* Fix typo

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

10 files changed, 305 insertions, 14 deletions

diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H
new file mode 100644
index 000000000..e59b1902c
--- /dev/null
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.H
@@ -0,0 +1,33 @@
+#ifndef WARPX_GALILEAN_ALGORITHM_H_
+#define WARPX_GALILEAN_ALGORITHM_H_
+
+#include <SpectralBaseAlgorithm.H>
+
+/* \brief Class that updates the field in spectral space
+ * and stores the coefficients of the corresponding update equation.
+ */
+class GalileanAlgorithm : public SpectralBaseAlgorithm
+{
+    public:
+        GalileanAlgorithm(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);
+        // Redefine update equation 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::Array<amrex::Real, 3>& v_galilean,
+                                    const amrex::Real dt);
+
+    private:
+        SpectralRealCoefficients C_coef, S_ck_coef;
+        SpectralComplexCoefficients Theta2_coef, X1_coef, X2_coef, X3_coef, X4_coef;
+};
+
+#endif // WARPX_GALILEAN_ALGORITHM_H_
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp
new file mode 100644
index 000000000..f869da90c
--- /dev/null
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanAlgorithm.cpp
@@ -0,0 +1,239 @@
+#include <GalileanAlgorithm.H>
+#include <WarpXConst.H>
+#include <cmath>
+
+using namespace amrex;
+
+/* \brief Initialize coefficients for the update equation */
+GalileanAlgorithm::GalileanAlgorithm(const SpectralKSpace& spectral_kspace,
+                         const DistributionMapping& dm,
+                         const int norder_x, const int norder_y,
+                         const int norder_z, const bool nodal,
+                         const Array<Real, 3>& v_galilean,
+                         const Real dt)
+     // Initialize members of base class
+     : SpectralBaseAlgorithm( spectral_kspace, dm,
+                              norder_x, norder_y, norder_z, nodal )
+{
+    const BoxArray& ba = spectral_kspace.spectralspace_ba;
+
+    // Allocate the arrays of coefficients
+    C_coef = SpectralRealCoefficients(ba, dm, 1, 0);
+    S_ck_coef = SpectralRealCoefficients(ba, dm, 1, 0);
+    X1_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+    X2_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+    X3_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+    X4_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+    Theta2_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
+
+    InitializeSpectralCoefficients(spectral_kspace, dm, v_galilean, dt);
+
+};
+
+/* Advance the E and B field in spectral space (stored in `f`)
+ * over one time step */
+void
+GalileanAlgorithm::pushSpectralFields(SpectralFieldData& f) const{
+
+    // Loop over boxes
+    for (MFIter mfi(f.fields); mfi.isValid(); ++mfi){
+
+        const Box& bx = f.fields[mfi].box();
+
+        // Extract arrays for the fields to be updated
+        Array4<Complex> fields = f.fields[mfi].array();
+        // Extract arrays for the coefficients
+        Array4<const Real> C_arr = C_coef[mfi].array();
+        Array4<const Real> S_ck_arr = S_ck_coef[mfi].array();
+        Array4<const Complex> X1_arr = X1_coef[mfi].array();
+        Array4<const Complex> X2_arr = X2_coef[mfi].array();
+        Array4<const Complex> X3_arr = X3_coef[mfi].array();
+        Array4<const Complex> X4_arr = X4_coef[mfi].array();
+        Array4<const Complex> Theta2_arr = Theta2_coef[mfi].array();
+
+        // Extract pointers for the k vectors
+        const Real* modified_kx_arr = modified_kx_vec[mfi].dataPtr();
+#if (AMREX_SPACEDIM==3)
+        const Real* modified_ky_arr = modified_ky_vec[mfi].dataPtr();
+#endif
+        const Real* 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
+        {
+            // Record old values of the fields to be updated
+            using Idx = SpectralFieldIndex;
+            const Complex Ex_old = fields(i,j,k,Idx::Ex);
+            const Complex Ey_old = fields(i,j,k,Idx::Ey);
+            const Complex Ez_old = fields(i,j,k,Idx::Ez);
+            const Complex Bx_old = fields(i,j,k,Idx::Bx);
+            const Complex By_old = fields(i,j,k,Idx::By);
+            const Complex Bz_old = fields(i,j,k,Idx::Bz);
+            // Shortcut for the values of J and rho
+            const Complex Jx = fields(i,j,k,Idx::Jx);
+            const Complex Jy = fields(i,j,k,Idx::Jy);
+            const Complex Jz = fields(i,j,k,Idx::Jz);
+            const Complex rho_old = fields(i,j,k,Idx::rho_old);
+            const Complex rho_new = fields(i,j,k,Idx::rho_new);
+            // k vector values, and coefficients
+            const Real kx = modified_kx_arr[i];
+#if (AMREX_SPACEDIM==3)
+            const Real ky = modified_ky_arr[j];
+            const Real kz = modified_kz_arr[k];
+#else
+            constexpr Real ky = 0;
+            const Real kz = modified_kz_arr[j];
+#endif
+            constexpr Real c2 = PhysConst::c*PhysConst::c;
+            constexpr Complex I = Complex{0,1};
+            const Real C = C_arr(i,j,k);
+            const Real S_ck = S_ck_arr(i,j,k);
+            const Complex X1 = X1_arr(i,j,k);
+            const Complex X2 = X2_arr(i,j,k);
+            const Complex X3 = X3_arr(i,j,k);
+            const Complex X4 = X4_arr(i,j,k);
+            const Complex T2 = Theta2_arr(i,j,k);
+
+            // Update E (see the original Galilean article)
+            fields(i,j,k,Idx::Ex) = T2*C*Ex_old
+                        + T2*S_ck*c2*I*(ky*Bz_old - kz*By_old)
+                        + X4*Jx - I*(X2*rho_new - T2*X3*rho_old)*kx;
+            fields(i,j,k,Idx::Ey) = T2*C*Ey_old
+                        + T2*S_ck*c2*I*(kz*Bx_old - kx*Bz_old)
+                        + X4*Jy - I*(X2*rho_new - T2*X3*rho_old)*ky;
+            fields(i,j,k,Idx::Ez) = T2*C*Ez_old
+                        + T2*S_ck*c2*I*(kx*By_old - ky*Bx_old)
+                        + X4*Jz - I*(X2*rho_new - T2*X3*rho_old)*kz;
+            // Update B (see the original Galilean article)
+            // Note: here X1 is T2*x1/(ep0*c*c*k_norm*k_norm), where
+            // x1 has the same definition as in the original paper
+            fields(i,j,k,Idx::Bx) = T2*C*Bx_old
+                        - T2*S_ck*I*(ky*Ez_old - kz*Ey_old)
+                        +      X1*I*(ky*Jz     - kz*Jy);
+            fields(i,j,k,Idx::By) = T2*C*By_old
+                        - T2*S_ck*I*(kz*Ex_old - kx*Ez_old)
+                        +      X1*I*(kz*Jx     - kx*Jz);
+            fields(i,j,k,Idx::Bz) = T2*C*Bz_old
+                        - T2*S_ck*I*(kx*Ey_old - ky*Ex_old)
+                        +      X1*I*(kx*Jy     - ky*Jx);
+        });
+    }
+};
+
+
+void GalileanAlgorithm::InitializeSpectralCoefficients(const SpectralKSpace& spectral_kspace,
+                                    const amrex::DistributionMapping& dm,
+                                    const Array<Real, 3>& v_galilean,
+                                    const amrex::Real dt)
+{
+    const BoxArray& ba = spectral_kspace.spectralspace_ba;
+    // Fill them with the right values:
+    // Loop over boxes and allocate the corresponding coefficients
+    // for each box owned by the local MPI proc
+    for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi){
+
+        const Box& bx = ba[mfi];
+
+        // Extract pointers for the k vectors
+        const Real* modified_kx = modified_kx_vec[mfi].dataPtr();
+    #if (AMREX_SPACEDIM==3)
+        const Real* modified_ky = modified_ky_vec[mfi].dataPtr();
+    #endif
+        const Real* modified_kz = modified_kz_vec[mfi].dataPtr();
+        // Extract arrays for the coefficients
+        Array4<Real> C = C_coef[mfi].array();
+        Array4<Real> S_ck = S_ck_coef[mfi].array();
+        Array4<Complex> X1 = X1_coef[mfi].array();
+        Array4<Complex> X2 = X2_coef[mfi].array();
+        Array4<Complex> X3 = X3_coef[mfi].array();
+        Array4<Complex> X4 = X4_coef[mfi].array();
+        Array4<Complex> Theta2 = Theta2_coef[mfi].array();
+        // Extract reals (for portability on GPU)
+        Real vx = v_galilean[0];
+        Real vy = v_galilean[1];
+        Real vz = v_galilean[2];
+
+        // Loop over indices within one box
+        ParallelFor(bx,
+        [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+        {
+            // Calculate norm of vector
+            const Real k_norm = std::sqrt(
+                std::pow(modified_kx[i], 2) +
+    #if (AMREX_SPACEDIM==3)
+                std::pow(modified_ky[j], 2) +
+                std::pow(modified_kz[k], 2));
+    #else
+                std::pow(modified_kz[j], 2));
+    #endif
+
+            // Calculate coefficients
+            constexpr Real c = PhysConst::c;
+            constexpr Real ep0 = PhysConst::ep0;
+            const Complex I{0.,1.};
+            if (k_norm != 0){
+
+                C(i,j,k) = std::cos(c*k_norm*dt);
+                S_ck(i,j,k) = std::sin(c*k_norm*dt)/(c*k_norm);
+
+                // Calculate dot product with galilean velocity
+                const Real kv = modified_kx[i]*vx +
+    #if (AMREX_SPACEDIM==3)
+                                 modified_ky[j]*vy +
+                                 modified_kz[k]*vz;
+    #else
+                                 modified_kz[j]*vz;
+    #endif
+
+                const Real nu = kv/(k_norm*c);
+                const Complex theta = MathFunc::exp( 0.5*I*kv*dt );
+                const Complex theta_star = MathFunc::exp( -0.5*I*kv*dt );
+                const Complex e_theta = MathFunc::exp( I*c*k_norm*dt );
+
+                Theta2(i,j,k) = theta*theta;
+
+                if ( (nu != 1.) && (nu != 0) ) {
+
+                    // Note: the coefficients X1, X2, X3 do not correspond
+                    // exactly to the original Galilean paper, but the
+                    // update equation have been modified accordingly so that
+                    // the expressions/ below (with the update equations)
+                    // are mathematically equivalent to those of the paper.
+                    Complex x1 = 1./(1.-nu*nu) *
+                        (theta_star - C(i,j,k)*theta + I*kv*S_ck(i,j,k)*theta);
+                    // x1, above, is identical to the original paper
+                    X1(i,j,k) = theta*x1/(ep0*c*c*k_norm*k_norm);
+                    // The difference betwen X2 and X3 below, and those
+                    // from the original paper is the factor ep0*k_norm*k_norm
+                    X2(i,j,k) = (x1 - theta*(1 - C(i,j,k)))
+                                /(theta_star-theta)/(ep0*k_norm*k_norm);
+                    X3(i,j,k) = (x1 - theta_star*(1 - C(i,j,k)))
+                                /(theta_star-theta)/(ep0*k_norm*k_norm);
+                    X4(i,j,k) = I*kv*X1(i,j,k) - theta*theta*S_ck(i,j,k)/ep0;
+                }
+                if ( nu == 0) {
+                    X1(i,j,k) = (1. - C(i,j,k)) / (ep0*c*c*k_norm*k_norm);
+                    X2(i,j,k) = (1. - S_ck(i,j,k)/dt) / (ep0*k_norm*k_norm);
+                    X3(i,j,k) = (C(i,j,k) - S_ck(i,j,k)/dt) / (ep0*k_norm*k_norm);
+                    X4(i,j,k) = -S_ck(i,j,k)/ep0;
+                }
+                if ( nu == 1.) {
+                    X1(i,j,k) = (1. - e_theta*e_theta + 2.*I*c*k_norm*dt) / (4.*c*c*ep0*k_norm*k_norm);
+                    X2(i,j,k) = (3. - 4.*e_theta + e_theta*e_theta + 2.*I*c*k_norm*dt) / (4.*ep0*k_norm*k_norm*(1.- e_theta));
+                    X3(i,j,k) = (3. - 2./e_theta - 2.*e_theta + e_theta*e_theta - 2.*I*c*k_norm*dt) / (4.*ep0*(e_theta - 1.)*k_norm*k_norm);
+                    X4(i,j,k) = I*(-1. + e_theta*e_theta + 2.*I*c*k_norm*dt) / (4.*ep0*c*k_norm);
+                }
+
+            } else { // Handle k_norm = 0, by using the analytical limit
+                C(i,j,k) = 1.;
+                S_ck(i,j,k) = dt;
+                X1(i,j,k) = dt*dt/(2. * ep0);
+                X2(i,j,k) = c*c*dt*dt/(6. * ep0);
+                X3(i,j,k) = - c*c*dt*dt/(3. * ep0);
+                X4(i,j,k) = -dt/ep0;
+                Theta2(i,j,k) = 1.;
+            }
+        });
+    }
+}
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
index ee8376865..82763f501 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
@@ -1,8 +1,13 @@
 CEXE_headers += SpectralBaseAlgorithm.H
 CEXE_headers += PsatdAlgorithm.H
 CEXE_sources += PsatdAlgorithm.cpp
+
+CEXE_headers += GalileanAlgorithm.H
+CEXE_sources += GalileanAlgorithm.cpp
+
 CEXE_headers += PMLPsatdAlgorithm.H
 CEXE_sources += PMLPsatdAlgorithm.cpp
 
+
 INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms
 VPATH_LOCATIONS   += $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/SpectralAlgorithms
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H
index 4452002d1..a2b85262d 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H
@@ -33,7 +33,7 @@ class PMLPsatdAlgorithm : public SpectralBaseAlgorithm
         }
 
     private:
-        SpectralCoefficients C_coef, S_ck_coef;
+        SpectralRealCoefficients C_coef, S_ck_coef;
 
 };
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp
index bae74daf6..0c4c4d41a 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp
@@ -23,8 +23,8 @@ PMLPsatdAlgorithm::PMLPsatdAlgorithm(
     const BoxArray& ba = spectral_kspace.spectralspace_ba;
 
     // Allocate the arrays of coefficients
-    C_coef = SpectralCoefficients(ba, dm, 1, 0);
-    S_ck_coef = SpectralCoefficients(ba, dm, 1, 0);
+    C_coef = SpectralRealCoefficients(ba, dm, 1, 0);
+    S_ck_coef = SpectralRealCoefficients(ba, dm, 1, 0);
 
     InitializeSpectralCoefficients(spectral_kspace, dm, dt);
 }
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
index 4abe89d9d..5e9b3e7bf 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
@@ -34,7 +34,7 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm
                                     const amrex::Real dt);
 
     private:
-        SpectralCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef;
+        SpectralRealCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef;
 };
 
 #endif // WARPX_PSATD_ALGORITHM_H_
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
index 4f4963e95..b2675ff91 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
@@ -22,11 +22,11 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace,
     const BoxArray& ba = spectral_kspace.spectralspace_ba;
 
     // Allocate the arrays of coefficients
-    C_coef = SpectralCoefficients(ba, dm, 1, 0);
-    S_ck_coef = SpectralCoefficients(ba, dm, 1, 0);
-    X1_coef = SpectralCoefficients(ba, dm, 1, 0);
-    X2_coef = SpectralCoefficients(ba, dm, 1, 0);
-    X3_coef = SpectralCoefficients(ba, dm, 1, 0);
+    C_coef = SpectralRealCoefficients(ba, dm, 1, 0);
+    S_ck_coef = SpectralRealCoefficients(ba, dm, 1, 0);
+    X1_coef = SpectralRealCoefficients(ba, dm, 1, 0);
+    X2_coef = SpectralRealCoefficients(ba, dm, 1, 0);
+    X3_coef = SpectralRealCoefficients(ba, dm, 1, 0);
 
     InitializeSpectralCoefficients(spectral_kspace, dm, dt);
 }
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
index 2c4946190..ed5b98183 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
@@ -30,7 +30,10 @@ class SpectralBaseAlgorithm
 
     protected: // Meant to be used in the subclasses
 
-        using SpectralCoefficients = amrex::FabArray< amrex::BaseFab <amrex::Real> >;
+        using SpectralRealCoefficients = \
+            amrex::FabArray< amrex::BaseFab <amrex::Real> >;
+        using SpectralComplexCoefficients = \
+            amrex::FabArray< amrex::BaseFab <Complex> >;
 
         // Constructor
         SpectralBaseAlgorithm(const SpectralKSpace& spectral_kspace,
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
index a6375e483..65f975682 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
@@ -30,6 +30,7 @@ class SpectralSolver
                         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::RealVect dx, const amrex::Real dt,
                         const bool pml=false );
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
index ca7bd06a0..c24a7af69 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
@@ -7,8 +7,10 @@
 #include <SpectralKSpace.H>
 #include <SpectralSolver.H>
 #include <PsatdAlgorithm.H>
+#include <GalileanAlgorithm.H>
 #include <PMLPsatdAlgorithm.H>
 
+
 /* \brief Initialize the spectral Maxwell solver
  *
  * This function selects the spectral algorithm to be used, allocates the
@@ -28,6 +30,7 @@ SpectralSolver::SpectralSolver(
                 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::RealVect dx, const amrex::Real dt,
                 const bool pml ) {
 
@@ -40,13 +43,20 @@ SpectralSolver::SpectralSolver(
 
     // - Select the algorithm depending on the input parameters
     //   Initialize the corresponding coefficients over k space
+
     if (pml) {
         algorithm = std::unique_ptr<PMLPsatdAlgorithm>( new PMLPsatdAlgorithm(
             k_space, dm, norder_x, norder_y, norder_z, nodal, dt ) );
-    } else {
-        algorithm = std::unique_ptr<PsatdAlgorithm>( new PsatdAlgorithm(
-            k_space, dm, norder_x, norder_y, norder_z, nodal, dt ) );
-    }
+    } else if ((v_galilean[0]==0) && (v_galilean[1]==0) && (v_galilean[2]==0)){
+         // v_galilean is 0: use standard PSATD algorithm
+         algorithm = std::unique_ptr<PsatdAlgorithm>( new PsatdAlgorithm(
+             k_space, dm, norder_x, norder_y, norder_z, nodal, dt ) );
+      } else {
+          // Otherwise: use the Galilean algorithm
+          algorithm = std::unique_ptr<GalileanAlgorithm>( new GalileanAlgorithm(
+              k_space, dm, norder_x, norder_y, norder_z, nodal, v_galilean, dt ));
+       }
+
 
     // - Initialize arrays for fields in spectral space + FFT plans
     field_data = SpectralFieldData( realspace_ba, k_space, dm,