Remove code for hybrid PSATD, as well as remaining Fortran interface file (#927)

* Clean-up hybrid code

* Remove Fortran interface

18 files changed, 87 insertions, 1156 deletions

diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp
index b73447cab..770c7cadd 100644
--- a/Source/Diagnostics/FieldIO.cpp
+++ b/Source/Diagnostics/FieldIO.cpp
@@ -578,8 +578,6 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames,
                 const BoxArray& ba = amrex::convert(boxArray(lev),IntVect::TheUnitVector());
                 MultiFab divE( ba, DistributionMap(lev), 1, 0 );
 #ifdef WARPX_USE_PSATD
-                AMREX_ALWAYS_ASSERT_WITH_MESSAGE( fft_hybrid_mpi_decomposition == false,
-                    "Spectral computation of div(E) not implemented with PICSAR spectral solver" );
                 spectral_solver_fp[lev]->ComputeSpectralDivE( Efield_aux[lev], divE );
 #else
                 m_fdtd_solver_fp[lev]->ComputeDivE( Efield_aux[lev], divE );
diff --git a/Source/Diagnostics/WarpXOpenPMD.cpp b/Source/Diagnostics/WarpXOpenPMD.cpp
index dc5342a78..caaac36ff 100644
--- a/Source/Diagnostics/WarpXOpenPMD.cpp
+++ b/Source/Diagnostics/WarpXOpenPMD.cpp
@@ -586,7 +586,7 @@ WarpXOpenPMDPlot::WriteOpenPMDFields( //const std::string& filename,
   auto meshes = series_iteration.meshes;
   meshes.setAttribute( "fieldSolver", [](){
 #ifdef WARPX_USE_PSATD
-      return "PSATD"; // TODO double-check if WARPX_USE_PSATD_HYBRID is covered
+      return "PSATD";
 #else
       switch( WarpX::particle_pusher_algo ) {
           case MaxwellSolverAlgo::Yee : return "Yee";
diff --git a/Source/FieldSolver/Make.package b/Source/FieldSolver/Make.package
index 03c2f830a..ba5eb1ac3 100644
--- a/Source/FieldSolver/Make.package
+++ b/Source/FieldSolver/Make.package
@@ -3,9 +3,6 @@ CEXE_sources += ElectrostaticSolver.cpp
 CEXE_sources += WarpX_QED_Field_Pushers.cpp
 ifeq ($(USE_PSATD),TRUE)
   include $(WARPX_HOME)/Source/FieldSolver/SpectralSolver/Make.package
-  ifeq ($(USE_PSATD_PICSAR),TRUE)
-    include $(WARPX_HOME)/Source/FieldSolver/PicsarHybridSpectralSolver/Make.package
-  endif
 endif
 include $(WARPX_HOME)/Source/FieldSolver/FiniteDifferenceSolver/Make.package
 
diff --git a/Source/FieldSolver/PicsarHybridSpectralSolver/Make.package b/Source/FieldSolver/PicsarHybridSpectralSolver/Make.package
deleted file mode 100644
index 22d38025c..000000000
--- a/Source/FieldSolver/PicsarHybridSpectralSolver/Make.package
+++ /dev/null
@@ -1,4 +0,0 @@
-CEXE_sources += PicsarHybridSpectralSolver.cpp
-F90EXE_sources += picsar_hybrid_spectral.F90
-
-VPATH_LOCATIONS   += $(WARPX_HOME)/Source/FieldSolver/PicsarHybridSpectralSolver
diff --git a/Source/FieldSolver/PicsarHybridSpectralSolver/PicsarHybridFFTData.H b/Source/FieldSolver/PicsarHybridSpectralSolver/PicsarHybridFFTData.H
deleted file mode 100644
index 7e037cf8d..000000000
--- a/Source/FieldSolver/PicsarHybridSpectralSolver/PicsarHybridFFTData.H
+++ /dev/null
@@ -1,43 +0,0 @@
-/* Copyright 2019 Remi Lehe
- *
- * This file is part of WarpX.
- *
- * License: BSD-3-Clause-LBNL
- */
-#ifndef WARPX_PICSAR_HYBRID_FFTDATA_H_
-#define WARPX_PICSAR_HYBRID_FFTDATA_H_
-
-#ifdef WARPX_USE_PSATD
-#   include <fftw3.h>
-
-
-// FFTData is a stuct containing a 22 pointers to auxiliary arrays
-// 1-11: padded arrays in real space (required by FFTW); 12-22: arrays in spectral space
-struct FFTData
-{
-    static constexpr int N = 22;
-    void* data[N] = { nullptr };
-
-    ~FFTData () {
-        for (int i = 0; i < N; ++i) { // The memory is allocated with fftw_alloc.
-            fftw_free(data[i]);
-            data[i] = nullptr;
-        }
-    }
-
-    FFTData () = default;
-
-    FFTData (FFTData && rhs) noexcept {
-        for (int i = 0; i < N; ++i) {
-            data[i] = rhs.data[i];
-            rhs.data[i] = nullptr;
-        }
-    }
-
-    FFTData (FFTData const&) = delete;
-    void operator= (FFTData const&) = delete;
-    void operator= (FFTData&&) = delete;
-};
-
-#endif // WARPX_USE_PSATD
-#endif // WARPX_PICSAR_HYBRID_FFTDATA_H_
diff --git a/Source/FieldSolver/PicsarHybridSpectralSolver/PicsarHybridSpectralSolver.cpp b/Source/FieldSolver/PicsarHybridSpectralSolver/PicsarHybridSpectralSolver.cpp
deleted file mode 100644
index db6942229..000000000
--- a/Source/FieldSolver/PicsarHybridSpectralSolver/PicsarHybridSpectralSolver.cpp
+++ /dev/null
@@ -1,459 +0,0 @@
-/* Copyright 2019 Axel Huebl, Maxence Thevenet, Remi Lehe
- * Revathi Jambunathan, Weiqun Zhang
- *
- * This file is part of WarpX.
- *
- * License: BSD-3-Clause-LBNL
- */
-#include "PicsarHybridFFTData.H"
-#include "WarpX.H"
-#include "FortranInterface/WarpX_f.H"
-#include <AMReX_iMultiFab.H>
-
-
-#ifdef WARPX_USE_PSATD
-
-using namespace amrex;
-
-constexpr int FFTData::N;
-
-namespace {
-static std::unique_ptr<FFTData> nullfftdata; // This for process with nz_fft=0
-
-/** \brief Returns an "owner mask" which 1 for all cells, except
- *  for the duplicated (physical) cells of a nodal grid.
- *
- *  More precisely, for these cells (which are represented on several grids)
- *  the owner mask is 1 only if these cells are at the lower left end of
- *  the local grid - or if these cells are at the end of the physical domain
- *  Therefore, there for these cells, there will be only one grid for
- *  which the owner mask is non-zero.
- */
-static iMultiFab
-BuildFFTOwnerMask (const MultiFab& mf, const Geometry& geom)
-{
-    const BoxArray& ba = mf.boxArray();
-    const DistributionMapping& dm = mf.DistributionMap();
-    iMultiFab mask(ba, dm, 1, 0);
-    const int owner = 1;
-    const int nonowner = 0;
-    mask.setVal(owner);
-
-    const Box& domain_box = amrex::convert(geom.Domain(), ba.ixType());
-
-    AMREX_ASSERT(ba.complementIn(domain_box).isEmpty());
-
-#ifdef _OPENMP
-#pragma omp parallel
-#endif
-    for (MFIter mfi(mask); mfi.isValid(); ++mfi)
-    {
-        IArrayBox& fab = mask[mfi];
-        const Box& bx = fab.box();
-        Box bx2 = bx;
-        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
-            // Detect nodal dimensions
-            if (bx2.type(idim) == IndexType::NODE) {
-                // Make sure that this grid does not touch the end of
-                // the physical domain.
-                if (bx2.bigEnd(idim) < domain_box.bigEnd(idim)) {
-                    bx2.growHi(idim, -1);
-                }
-            }
-        }
-        const BoxList& bl = amrex::boxDiff(bx, bx2);
-        // Set owner mask in these cells
-        for (const auto& b : bl) {
-            fab.setVal(nonowner, b, 0, 1);
-        }
-
-    }
-
-    return mask;
-}
-
-/** \brief Copy the data from the FFT grid to the regular grid
- *
- * Because, for nodal grid, some cells are duplicated on several boxes,
- * special care has to be taken in order to have consistent values on
- * each boxes when copying this data. Here this is done by setting a
- * mask, where, for these duplicated cells, the mask is non-zero on only
- * one box.
- */
-static void
-CopyDataFromFFTToValid (MultiFab& mf, const MultiFab& mf_fft, const BoxArray& ba_valid_fft, const Geometry& geom)
-{
-    auto idx_type = mf_fft.ixType();
-    MultiFab mftmp(amrex::convert(ba_valid_fft,idx_type), mf_fft.DistributionMap(), 1, 0);
-
-    const iMultiFab& mask = BuildFFTOwnerMask(mftmp, geom);
-
-    // Local copy: whenever an MPI rank owns both the data from the FFT
-    // grid and from the regular grid, for overlapping region, copy it locally
-#ifdef _OPENMP
-#pragma omp parallel
-#endif
-    for (MFIter mfi(mftmp,true); mfi.isValid(); ++mfi)
-    {
-        const Box& bx = mfi.tilebox();
-        FArrayBox& dstfab = mftmp[mfi];
-
-        const FArrayBox& srcfab = mf_fft[mfi];
-        const Box& srcbox = srcfab.box();
-
-        if (srcbox.contains(bx))
-        {
-            // Copy the interior region (without guard cells)
-            dstfab.copy(srcfab, bx, 0, bx, 0, 1);
-            // Set the value to 0 whenever the mask is 0
-            // (i.e. for nodal duplicated cells, there is a single box
-            // for which the mask is different than 0)
-            // if mask == 0, set value to zero
-            dstfab.setValIfNot(0.0, bx, mask[mfi], 0, 1);
-        }
-    }
-
-    // Global copy: Get the remaining the data from other procs
-    // Use ParallelAdd instead of ParallelCopy, so that the value from
-    // the cell that has non-zero mask is the one which is retained.
-    mf.setVal(0.0, 0);
-    mf.ParallelAdd(mftmp);
-
-
-}
-
-}
-
-void
-WarpX::AllocLevelDataFFT (int lev)
-{
-    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(lev == 0, "PSATD doesn't work with mesh refinement yet");
-
-    static_assert(std::is_standard_layout<FFTData>::value, "FFTData must have standard layout");
-    static_assert(sizeof(FFTData) == sizeof(void*)*FFTData::N, "sizeof FFTData is wrong");
-
-
-
-    InitFFTComm(lev);
-
-    BoxArray ba_fp_fft;
-    DistributionMapping dm_fp_fft;
-    FFTDomainDecomposition(lev, ba_fp_fft, dm_fp_fft, ba_valid_fp_fft[lev], domain_fp_fft[lev],
-                           geom[lev].Domain());
-
-    // rho2 has one extra ghost cell, so that it's safe to deposit charge density after
-    // pushing particle.
-
-    Efield_fp_fft[lev][0].reset(new MultiFab(amrex::convert(ba_fp_fft,Ex_nodal_flag),
-                                             dm_fp_fft, 1, 0));
-    Efield_fp_fft[lev][1].reset(new MultiFab(amrex::convert(ba_fp_fft,Ey_nodal_flag),
-                                             dm_fp_fft, 1, 0));
-    Efield_fp_fft[lev][2].reset(new MultiFab(amrex::convert(ba_fp_fft,Ez_nodal_flag),
-                                             dm_fp_fft, 1, 0));
-    Bfield_fp_fft[lev][0].reset(new MultiFab(amrex::convert(ba_fp_fft,Bx_nodal_flag),
-                                             dm_fp_fft, 1, 0));
-    Bfield_fp_fft[lev][1].reset(new MultiFab(amrex::convert(ba_fp_fft,By_nodal_flag),
-                                             dm_fp_fft, 1, 0));
-    Bfield_fp_fft[lev][2].reset(new MultiFab(amrex::convert(ba_fp_fft,Bz_nodal_flag),
-                                             dm_fp_fft, 1, 0));
-    current_fp_fft[lev][0].reset(new MultiFab(amrex::convert(ba_fp_fft,jx_nodal_flag),
-                                              dm_fp_fft, 1, 0));
-    current_fp_fft[lev][1].reset(new MultiFab(amrex::convert(ba_fp_fft,jy_nodal_flag),
-                                              dm_fp_fft, 1, 0));
-    current_fp_fft[lev][2].reset(new MultiFab(amrex::convert(ba_fp_fft,jz_nodal_flag),
-                                              dm_fp_fft, 1, 0));
-    rho_fp_fft[lev].reset(new MultiFab(amrex::convert(ba_fp_fft,IntVect::TheNodeVector()),
-                                       dm_fp_fft, 2, 0));
-
-    dataptr_fp_fft[lev].reset(new LayoutData<FFTData>(ba_fp_fft, dm_fp_fft));
-
-    if (lev > 0)
-    {
-        BoxArray ba_cp_fft;
-        DistributionMapping dm_cp_fft;
-        FFTDomainDecomposition(lev, ba_cp_fft, dm_cp_fft, ba_valid_cp_fft[lev], domain_cp_fft[lev],
-                               amrex::coarsen(geom[lev].Domain(),2));
-
-        Efield_cp_fft[lev][0].reset(new MultiFab(amrex::convert(ba_cp_fft,Ex_nodal_flag),
-                                                 dm_cp_fft, 1, 0));
-        Efield_cp_fft[lev][1].reset(new MultiFab(amrex::convert(ba_cp_fft,Ey_nodal_flag),
-                                                 dm_cp_fft, 1, 0));
-        Efield_cp_fft[lev][2].reset(new MultiFab(amrex::convert(ba_cp_fft,Ez_nodal_flag),
-                                                 dm_cp_fft, 1, 0));
-        Bfield_cp_fft[lev][0].reset(new MultiFab(amrex::convert(ba_cp_fft,Bx_nodal_flag),
-                                                 dm_cp_fft, 1, 0));
-        Bfield_cp_fft[lev][1].reset(new MultiFab(amrex::convert(ba_cp_fft,By_nodal_flag),
-                                                 dm_cp_fft, 1, 0));
-        Bfield_cp_fft[lev][2].reset(new MultiFab(amrex::convert(ba_cp_fft,Bz_nodal_flag),
-                                                 dm_cp_fft, 1, 0));
-        current_cp_fft[lev][0].reset(new MultiFab(amrex::convert(ba_cp_fft,jx_nodal_flag),
-                                                  dm_cp_fft, 1, 0));
-        current_cp_fft[lev][1].reset(new MultiFab(amrex::convert(ba_cp_fft,jy_nodal_flag),
-                                                  dm_cp_fft, 1, 0));
-        current_cp_fft[lev][2].reset(new MultiFab(amrex::convert(ba_cp_fft,jz_nodal_flag),
-                                                  dm_cp_fft, 1, 0));
-        rho_cp_fft[lev].reset(new MultiFab(amrex::convert(ba_cp_fft,IntVect::TheNodeVector()),
-                                           dm_cp_fft, 2, 0));
-
-        dataptr_cp_fft[lev].reset(new LayoutData<FFTData>(ba_cp_fft, dm_cp_fft));
-    }
-
-    InitFFTDataPlan(lev);
-}
-
-/** \brief Create MPI sub-communicators for each FFT group,
- *         and put them in PICSAR module
- *
- * These communicators are passed to the parallel FFTW library, in order
- * to perform a global FFT within each FFT group.
- */
-void
-WarpX::InitFFTComm (int lev)
-{
-    int nprocs = ParallelDescriptor::NProcs();
-    ngroups_fft = std::min(ngroups_fft, nprocs);
-
-    // # of processes in the subcommunicator
-    int np_fft = nprocs / ngroups_fft;
-    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(np_fft*ngroups_fft == nprocs,
-        "Number of processes must be divisible by number of FFT groups");
-
-    int myproc = ParallelDescriptor::MyProc();
-    // my color in ngroups_fft subcommunicators.  0 <= color_fft < ngroups_fft
-    color_fft[lev] = myproc / np_fft;
-    MPI_Comm_split(ParallelDescriptor::Communicator(), color_fft[lev], myproc, &comm_fft[lev]);
-
-    int fcomm = MPI_Comm_c2f(comm_fft[lev]);
-    // Set the communicator of the PICSAR module to the one we just created
-    warpx_fft_mpi_init(fcomm);
-}
-
-/** \brief Perform domain decomposition for the FFTW
- *
- *  Attribute one (unique) box to each proc, in such a way that:
- *    - The global domain is divided among FFT groups,
- *      with additional guard cells around each FFT group
- *    - The domain associated to an FFT group (with its guard cells)
- *      is further divided in sub-subdomains along z, so as to distribute
- *      it among the procs within an FFT group
- *
- *  The attribution is done by setting (within this function):
- *  - ba_fft: the BoxArray representing the final set of sub-domains for the FFT
- *            (includes/covers the guard cells of the FFT groups)
- *  - dm_fft: the mapping between these sub-domains and the corresponding proc
- *              (imposes one unique box for each proc)
- *  - ba_valid: the BoxArray that contains valid part of the sub-domains of ba_fft
- *            (i.e. does not include/cover the guard cells of the FFT groups)
- *  - domain_fft: a Box that represent the domain of the FFT group for the current proc
- */
-void
-WarpX::FFTDomainDecomposition (int lev, BoxArray& ba_fft, DistributionMapping& dm_fft,
-                               BoxArray& ba_valid, Box& domain_fft, const Box& domain)
-{
-
-    IntVect nguards_fft(AMREX_D_DECL(nox_fft/2,noy_fft/2,noz_fft/2));
-
-    int nprocs = ParallelDescriptor::NProcs();
-
-    BoxList bl(domain, ngroups_fft);  // This does a multi-D domain decomposition for groups
-    AMREX_ALWAYS_ASSERT(bl.size() == ngroups_fft);
-    const Vector<Box>& bldata = bl.data();
-
-    // This is the domain for the FFT sub-group (including guard cells)
-    domain_fft = amrex::grow(bldata[color_fft[lev]], nguards_fft);
-    // Ask FFTW to chop the current FFT sub-group domain in the z-direction
-    // and give a chunk to each MPI rank in the current sub-group.
-    int nz_fft, z0_fft;
-
-    warpx_fft_domain_decomp(&nz_fft, &z0_fft, WARPX_TO_FORTRAN_BOX(domain_fft));
-    // Each MPI rank adds a box with its chunk of the FFT grid
-    // (given by the above decomposition) to the list `bx_fft`,
-    // then list is shared among all MPI ranks via AllGather
-    Vector<Box> bx_fft;
-    if (nz_fft > 0) {
-        Box b = domain_fft;
-        b.setRange(AMREX_SPACEDIM-1, z0_fft+domain_fft.smallEnd(AMREX_SPACEDIM-1), nz_fft);
-        bx_fft.push_back(b);
-    } else {
-        // Add empty box for the AllGather call
-        bx_fft.push_back(Box());
-    }
-    amrex::AllGatherBoxes(bx_fft);
-    AMREX_ASSERT(bx_fft.size() == ParallelDescriptor::NProcs());
-    // Build pmap and bx_fft without the empty boxes
-    Vector<int> pmap;
-    for (int i = 0; i < bx_fft.size(); ++i) {
-        if (bx_fft[i].ok()) {
-            pmap.push_back(i);
-        }
-    }
-    bx_fft.erase(std::remove_if(bx_fft.begin(),bx_fft.end(),
-                                [](Box const& b) { return b.isEmpty(); }),
-                 bx_fft.end());
-    AMREX_ASSERT(bx_fft.size() == pmap.size());
-
-    // Define the AMReX objects for the FFT grid: BoxArray and DistributionMapping
-    ba_fft.define(BoxList(std::move(bx_fft)));
-    dm_fft.define(std::move(pmap));
-
-    // For communication between WarpX normal domain and FFT domain, we need to create a
-    // special BoxArray ba_valid
-    const Box foobox(-nguards_fft-2, -nguards_fft-2);
-
-    BoxList bl_valid; // List of boxes: will be filled by the valid part of the subdomains of ba_fft
-    bl_valid.reserve(ba_fft.size());
-    int np_fft = nprocs / ngroups_fft;
-    for (int i = 0; i < ba_fft.size(); ++i)
-    {
-        int igroup = dm_fft[i] / np_fft; // This should be consistent with InitFFTComm
-        const Box& bx = ba_fft[i] & bldata[igroup]; // Intersection with the domain of
-                                                    // the FFT group *without* guard cells
-        if (bx.ok())
-        {
-            bl_valid.push_back(bx);
-        }
-        else
-        {
-            bl_valid.push_back(foobox);
-        }
-    }
-
-    ba_valid.define(std::move(bl_valid));
-}
-
-/** /brief Set all the flags and metadata of the PICSAR FFT module.
- *         Allocate the auxiliary arrays of `fft_data`
- *
- * Note: dataptr_data is a stuct containing 22 pointers to arrays
- * 1-11: padded arrays in real space ; 12-22 arrays for the fields in Fourier space
- */
-void
-WarpX::InitFFTDataPlan (int lev)
-{
-    auto dx_fp = CellSize(lev);
-
-    if (Efield_fp_fft[lev][0]->local_size() == 1)
-      //Only one FFT patch on this MPI
-    {
-        for (MFIter mfi(*Efield_fp_fft[lev][0]); mfi.isValid(); ++mfi)
-        {
-            warpx_fft_dataplan_init(&nox_fft, &noy_fft, &noz_fft,
-                                    (*dataptr_fp_fft[lev])[mfi].data, &FFTData::N,
-                                    dx_fp.data(), &dt[lev], &fftw_plan_measure, &WarpX::do_nodal );
-        }
-    }
-    else if (Efield_fp_fft[lev][0]->local_size() == 0)
-      // No FFT patch on this MPI rank (may happen with FFTW)
-      // Still need to call the MPI-FFT initialization routines
-    {
-        nullfftdata.reset(new FFTData());
-        warpx_fft_dataplan_init(&nox_fft, &noy_fft, &noz_fft,
-                                nullfftdata->data, &FFTData::N,
-                                dx_fp.data(), &dt[lev], &fftw_plan_measure,
-                                &WarpX::do_nodal );
-    }
-    else
-    {
-      // Multiple FFT patches on this MPI rank
-        amrex::Abort("WarpX::InitFFTDataPlan: TODO");
-    }
-
-    if (lev > 0)
-    {
-        amrex::Abort("WarpX::InitFFTDataPlan: TODO");
-    }
-}
-
-void
-WarpX::FreeFFT (int lev)
-{
-    nullfftdata.reset();
-
-    warpx_fft_nullify();
-
-    if (comm_fft[lev] != MPI_COMM_NULL) {
-        MPI_Comm_free(&comm_fft[lev]);
-    }
-    comm_fft[lev] = MPI_COMM_NULL;
-}
-
-void
-WarpX::PushPSATD_hybridFFT (int lev, amrex::Real /* dt */)
-{
-    WARPX_PROFILE_VAR_NS("WarpXFFT::CopyDualGrid", blp_copy);
-    WARPX_PROFILE_VAR_NS("PICSAR::FftPushEB", blp_push_eb);
-
-    auto period_fp = geom[lev].periodicity();
-
-    WARPX_PROFILE_VAR_START(blp_copy);
-    Efield_fp_fft[lev][0]->ParallelCopy(*Efield_fp[lev][0], 0, 0, 1, Efield_fp[lev][0]->nGrow(), 0, period_fp);
-    Efield_fp_fft[lev][1]->ParallelCopy(*Efield_fp[lev][1], 0, 0, 1, Efield_fp[lev][1]->nGrow(), 0, period_fp);
-    Efield_fp_fft[lev][2]->ParallelCopy(*Efield_fp[lev][2], 0, 0, 1, Efield_fp[lev][2]->nGrow(), 0, period_fp);
-    Bfield_fp_fft[lev][0]->ParallelCopy(*Bfield_fp[lev][0], 0, 0, 1, Bfield_fp[lev][0]->nGrow(), 0, period_fp);
-    Bfield_fp_fft[lev][1]->ParallelCopy(*Bfield_fp[lev][1], 0, 0, 1, Bfield_fp[lev][1]->nGrow(), 0, period_fp);
-    Bfield_fp_fft[lev][2]->ParallelCopy(*Bfield_fp[lev][2], 0, 0, 1, Bfield_fp[lev][2]->nGrow(), 0, period_fp);
-    current_fp_fft[lev][0]->ParallelCopy(*current_fp[lev][0], 0, 0, 1, current_fp[lev][0]->nGrow(), 0, period_fp);
-    current_fp_fft[lev][1]->ParallelCopy(*current_fp[lev][1], 0, 0, 1, current_fp[lev][1]->nGrow(), 0, period_fp);
-    current_fp_fft[lev][2]->ParallelCopy(*current_fp[lev][2], 0, 0, 1, current_fp[lev][2]->nGrow(), 0, period_fp);
-    rho_fp_fft[lev]->ParallelCopy(*rho_fp[lev], 0, 0, 2, rho_fp[lev]->nGrow(), 0, period_fp);
-    WARPX_PROFILE_VAR_STOP(blp_copy);
-
-    WARPX_PROFILE_VAR_START(blp_push_eb);
-    if (Efield_fp_fft[lev][0]->local_size() == 1)
-       //Only one FFT patch on this MPI
-    {
-        for (MFIter mfi(*Efield_fp_fft[lev][0]); mfi.isValid(); ++mfi)
-        {
-                warpx_fft_push_eb(WARPX_TO_FORTRAN_ANYD((*Efield_fp_fft[lev][0])[mfi]),
-                                  WARPX_TO_FORTRAN_ANYD((*Efield_fp_fft[lev][1])[mfi]),
-                                  WARPX_TO_FORTRAN_ANYD((*Efield_fp_fft[lev][2])[mfi]),
-                                  WARPX_TO_FORTRAN_ANYD((*Bfield_fp_fft[lev][0])[mfi]),
-                                  WARPX_TO_FORTRAN_ANYD((*Bfield_fp_fft[lev][1])[mfi]),
-                                  WARPX_TO_FORTRAN_ANYD((*Bfield_fp_fft[lev][2])[mfi]),
-                                  WARPX_TO_FORTRAN_ANYD((*current_fp_fft[lev][0])[mfi]),
-                                  WARPX_TO_FORTRAN_ANYD((*current_fp_fft[lev][1])[mfi]),
-                                  WARPX_TO_FORTRAN_ANYD((*current_fp_fft[lev][2])[mfi]),
-                                  WARPX_TO_FORTRAN_N_ANYD((*rho_fp_fft[lev])[mfi],0),
-                                  WARPX_TO_FORTRAN_N_ANYD((*rho_fp_fft[lev])[mfi],1));
-        }
-    }
-    else if (Efield_fp_fft[lev][0]->local_size() == 0)
-      // No FFT patch on this MPI rank
-      // Still need to call the MPI-FFT routine.
-    {
-        FArrayBox fab(Box(IntVect::TheZeroVector(), IntVect::TheUnitVector()));
-        warpx_fft_push_eb(WARPX_TO_FORTRAN_ANYD(fab),
-                          WARPX_TO_FORTRAN_ANYD(fab),
-                          WARPX_TO_FORTRAN_ANYD(fab),
-                          WARPX_TO_FORTRAN_ANYD(fab),
-                          WARPX_TO_FORTRAN_ANYD(fab),
-                          WARPX_TO_FORTRAN_ANYD(fab),
-                          WARPX_TO_FORTRAN_ANYD(fab),
-                          WARPX_TO_FORTRAN_ANYD(fab),
-                          WARPX_TO_FORTRAN_ANYD(fab),
-                          WARPX_TO_FORTRAN_ANYD(fab),
-                          WARPX_TO_FORTRAN_ANYD(fab));
-    }
-    else
-      // Multiple FFT patches on this MPI rank
-    {
-        amrex::Abort("WarpX::PushPSATD: TODO");
-    }
-    WARPX_PROFILE_VAR_STOP(blp_push_eb);
-
-    WARPX_PROFILE_VAR_START(blp_copy);
-    CopyDataFromFFTToValid(*Efield_fp[lev][0], *Efield_fp_fft[lev][0], ba_valid_fp_fft[lev], geom[lev]);
-    CopyDataFromFFTToValid(*Efield_fp[lev][1], *Efield_fp_fft[lev][1], ba_valid_fp_fft[lev], geom[lev]);
-    CopyDataFromFFTToValid(*Efield_fp[lev][2], *Efield_fp_fft[lev][2], ba_valid_fp_fft[lev], geom[lev]);
-    CopyDataFromFFTToValid(*Bfield_fp[lev][0], *Bfield_fp_fft[lev][0], ba_valid_fp_fft[lev], geom[lev]);
-    CopyDataFromFFTToValid(*Bfield_fp[lev][1], *Bfield_fp_fft[lev][1], ba_valid_fp_fft[lev], geom[lev]);
-    CopyDataFromFFTToValid(*Bfield_fp[lev][2], *Bfield_fp_fft[lev][2], ba_valid_fp_fft[lev], geom[lev]);
-    WARPX_PROFILE_VAR_STOP(blp_copy);
-
-    if (lev > 0)
-    {
-        amrex::Abort("WarpX::PushPSATD: TODO");
-    }
-}
-
-#endif // WARPX_USE_PSATD
diff --git a/Source/FieldSolver/PicsarHybridSpectralSolver/picsar_hybrid_spectral.F90 b/Source/FieldSolver/PicsarHybridSpectralSolver/picsar_hybrid_spectral.F90
deleted file mode 100644
index 96bfb6111..000000000
--- a/Source/FieldSolver/PicsarHybridSpectralSolver/picsar_hybrid_spectral.F90
+++ /dev/null
@@ -1,326 +0,0 @@
-! Copyright 2019 Andrew Myers, Maxence Thevenet, Remi Lehe
-! Weiqun Zhang
-!
-! This file is part of WarpX.
-!
-! License: BSD-3-Clause-LBNL
-
-
-module warpx_fft_module
-  use amrex_error_module, only : amrex_error, amrex_abort
-  use amrex_fort_module, only : amrex_real
-  use iso_c_binding
-  implicit none
-
-  include 'fftw3-mpi.f03'
-
-  private
-  public :: warpx_fft_mpi_init, warpx_fft_domain_decomp, warpx_fft_dataplan_init, warpx_fft_nullify, &
-       warpx_fft_push_eb
-
-contains
-
-!> @brief
-!! Set the communicator of the PICSAR module to the one that is passed in argument
-  subroutine warpx_fft_mpi_init (fcomm) bind(c,name='warpx_fft_mpi_init')
-    use shared_data, only : comm, rank, nproc
-    integer, intent(in), value :: fcomm
-
-    integer :: ierr, lnproc, lrank
-
-    comm = fcomm
-
-    call mpi_comm_size(comm, lnproc, ierr)
-    nproc = lnproc
-
-    call mpi_comm_rank(comm, lrank, ierr)
-    rank = lrank
-
-#ifdef _OPENMP
-    ierr = fftw_init_threads()
-    if (ierr.eq.0) call amrex_error("fftw_init_threads failed")
-#endif
-    call fftw_mpi_init()
-#ifdef _OPENMP
-    call dfftw_init_threads(ierr)
-    if (ierr.eq.0) call amrex_error("dfftw_init_threads failed")
-#endif
-  end subroutine warpx_fft_mpi_init
-
-!> @brief
-!! Ask FFTW to do domain decomposition.
-!
-! This is always a 1d domain decomposition along z ; it is typically
-! done on the *FFT sub-groups*, not the all domain
-  subroutine warpx_fft_domain_decomp (warpx_local_nz, warpx_local_z0, global_lo, global_hi) &
-       bind(c,name='warpx_fft_domain_decomp')
-    use picsar_precision, only : idp
-    use shared_data, only : comm, &
-         nx_global, ny_global, nz_global, & ! size of global FFT
-         nx, ny, nz ! size of local subdomains
-    use mpi_fftw3, only : local_nz, local_z0, fftw_mpi_local_size_3d, alloc_local
-
-    integer, intent(out) :: warpx_local_nz, warpx_local_z0
-    integer, dimension(3), intent(in) :: global_lo, global_hi
-
-    nx_global = INT(global_hi(1)-global_lo(1)+1,idp)
-    ny_global = INT(global_hi(2)-global_lo(2)+1,idp)
-    nz_global = INT(global_hi(3)-global_lo(3)+1,idp)
-
-    alloc_local = fftw_mpi_local_size_3d( &
-         INT(nz_global,C_INTPTR_T), &
-         INT(ny_global,C_INTPTR_T), &
-         INT(nx_global,C_INTPTR_T)/2+1, &
-         comm, local_nz, local_z0)
-
-    nx = nx_global
-    ny = ny_global
-    nz = local_nz
-
-    warpx_local_nz = local_nz
-    warpx_local_z0 = local_z0
-  end subroutine warpx_fft_domain_decomp
-
-
-!> @brief
-!! Set all the flags and metadata of the PICSAR FFT module.
-!! Allocate the auxiliary arrays of `fft_data`
-!!
-!! Note: fft_data is a stuct containing 22 pointers to arrays
-!! 1-11: padded arrays in real space ; 12-22 arrays for the fields in Fourier space
-  subroutine warpx_fft_dataplan_init (nox, noy, noz, fft_data, ndata, dx_wrpx, dt_wrpx, fftw_measure, do_nodal) &
-       bind(c,name='warpx_fft_dataplan_init')
-    USE picsar_precision, only: idp
-    use shared_data, only : c_dim,  p3dfft_flag, fftw_plan_measure, &
-         fftw_with_mpi, fftw_threads_ok, fftw_hybrid, fftw_mpi_transpose, &
-         nx, ny, nz, & ! size of local subdomains
-         nkx, nky, nkz, & ! size of local ffts
-         iz_min_r, iz_max_r, iy_min_r, iy_max_r, ix_min_r, ix_max_r, & ! loop bounds
-         dx, dy, dz
-    use fields, only : nxguards, nyguards, nzguards,  & ! size of guard regions
-         ex_r, ey_r, ez_r, bx_r, by_r, bz_r, &
-         jx_r, jy_r, jz_r, rho_r, rhoold_r, &
-         exf, eyf, ezf, bxf, byf, bzf, &
-         jxf, jyf, jzf, rhof, rhooldf, &
-         l_spectral, l_staggered, norderx, nordery, norderz
-    use mpi_fftw3, only : alloc_local
-    use omp_lib, only: omp_get_max_threads
-    USE gpstd_solver, only: init_gpstd
-    USE fourier_psaotd, only: init_plans_fourier_mpi
-    use params, only : dt
-
-    integer, intent(in) :: nox, noy, noz, ndata
-    integer, intent(in) :: fftw_measure
-    integer, intent(in) :: do_nodal
-    type(c_ptr), intent(inout) :: fft_data(ndata)
-    real(c_double), intent(in) :: dx_wrpx(3), dt_wrpx
-
-    integer(idp) :: nopenmp
-    integer :: nx_padded
-    integer, dimension(3) :: shp
-    integer(kind=c_size_t) :: sz
-
-    ! No need to distinguish physical and guard cells for the global FFT;
-    ! only nx+2*nxguards counts. Thus we declare 0 guard cells for simplicity
-    nxguards = 0_idp
-    nyguards = 0_idp
-    nzguards = 0_idp
-
-    ! For the calculation of the modified [k] vectors
-    if (do_nodal == 0) then
-       l_staggered = .TRUE.
-    else
-       l_staggered = .FALSE.
-    endif
-    norderx = int(nox, idp)
-    nordery = int(noy, idp)
-    norderz = int(noz, idp)
-    ! Define parameters of FFT plans
-    c_dim = INT(AMREX_SPACEDIM,idp)   ! Dimensionality of the simulation (2d/3d)
-    fftw_with_mpi = .TRUE. ! Activate MPI FFTW
-    fftw_hybrid = .FALSE.   ! FFT per MPI subgroup (instead of global)
-    fftw_mpi_transpose = .FALSE. ! Do not transpose the data
-    fftw_plan_measure = (fftw_measure .ne. 0)
-    p3dfft_flag = .FALSE.
-    l_spectral  = .TRUE.   ! Activate spectral Solver, using FFT
-#ifdef _OPENMP
-    fftw_threads_ok = .TRUE.
-    nopenmp = OMP_GET_MAX_THREADS()
-#else
-    fftw_threads_ok = .FALSE.
-    nopenmp = 1
-#endif
-
-    ! Allocate padded arrays for MPI FFTW
-    nx_padded = 2*(nx/2 + 1)
-    shp = [nx_padded, int(ny), int(nz)]
-    sz = 2*alloc_local
-    fft_data(1) = fftw_alloc_real(sz)
-    call c_f_pointer(fft_data(1), ex_r, shp)
-    fft_data(2) = fftw_alloc_real(sz)
-    call c_f_pointer(fft_data(2), ey_r, shp)
-    fft_data(3) = fftw_alloc_real(sz)
-    call c_f_pointer(fft_data(3), ez_r, shp)
-    fft_data(4) = fftw_alloc_real(sz)
-    call c_f_pointer(fft_data(4), bx_r, shp)
-    fft_data(5) = fftw_alloc_real(sz)
-    call c_f_pointer(fft_data(5), by_r, shp)
-    fft_data(6) = fftw_alloc_real(sz)
-    call c_f_pointer(fft_data(6), bz_r, shp)
-    fft_data(7) = fftw_alloc_real(sz)
-    call c_f_pointer(fft_data(7), jx_r, shp)
-    fft_data(8) = fftw_alloc_real(sz)
-    call c_f_pointer(fft_data(8), jy_r, shp)
-    fft_data(9) = fftw_alloc_real(sz)
-    call c_f_pointer(fft_data(9), jz_r, shp)
-    fft_data(10) = fftw_alloc_real(sz)
-    call c_f_pointer(fft_data(10), rho_r, shp)
-    fft_data(11) = fftw_alloc_real(sz)
-    call c_f_pointer(fft_data(11), rhoold_r, shp)
-
-    ! Set array bounds when copying ex to ex_r in PICSAR
-    ix_min_r = 1; ix_max_r = nx
-    iy_min_r = 1; iy_max_r = ny
-    iz_min_r = 1; iz_max_r = nz
-    ! Allocate Fourier space fields of the same size
-    nkx = nx/2 + 1
-    nky = ny
-    nkz = nz
-    shp = [int(nkx), int(nky), int(nkz)]
-    sz = alloc_local
-    fft_data(12) = fftw_alloc_complex(sz)
-    call c_f_pointer(fft_data(12), exf, shp)
-    fft_data(13) = fftw_alloc_complex(sz)
-    call c_f_pointer(fft_data(13), eyf, shp)
-    fft_data(14) = fftw_alloc_complex(sz)
-    call c_f_pointer(fft_data(14), ezf, shp)
-    fft_data(15) = fftw_alloc_complex(sz)
-    call c_f_pointer(fft_data(15), bxf, shp)
-    fft_data(16) = fftw_alloc_complex(sz)
-    call c_f_pointer(fft_data(16), byf, shp)
-    fft_data(17) = fftw_alloc_complex(sz)
-    call c_f_pointer(fft_data(17), bzf, shp)
-    fft_data(18) = fftw_alloc_complex(sz)
-    call c_f_pointer(fft_data(18), jxf, shp)
-    fft_data(19) = fftw_alloc_complex(sz)
-    call c_f_pointer(fft_data(19), jyf, shp)
-    fft_data(20) = fftw_alloc_complex(sz)
-    call c_f_pointer(fft_data(20), jzf, shp)
-    fft_data(21) = fftw_alloc_complex(sz)
-    call c_f_pointer(fft_data(21), rhof, shp)
-    fft_data(22) = fftw_alloc_complex(sz)
-    call c_f_pointer(fft_data(22), rhooldf, shp)
-!$acc enter data create (exf,eyf,ezf,bxf,byf,bzf,jxf,jyf,jzf,rhof,rhooldf)
-    if (ndata < 22) then
-       call amrex_abort("size of fft_data is too small")
-    end if
-
-    dx = dx_wrpx(1)
-    dy = dx_wrpx(2)
-    dz = dx_wrpx(3)
-    dt = dt_wrpx
-
-    ! Initialize the matrix blocks for the PSATD solver
-    CALL init_gpstd()
-    ! Initialize the plans for fftw with MPI
-    CALL init_plans_fourier_mpi(nopenmp)
-
-  end subroutine warpx_fft_dataplan_init
-
-
-  subroutine warpx_fft_nullify () bind(c,name='warpx_fft_nullify')
-    use fields, only: ex_r, ey_r, ez_r, bx_r, by_r, bz_r, &
-         jx_r, jy_r, jz_r, rho_r, rhoold_r, &
-         exf, eyf, ezf, bxf, byf, bzf, &
-         jxf, jyf, jzf, rhof, rhooldf
-    use mpi_fftw3, only : plan_r2c_mpi, plan_c2r_mpi
-    nullify(ex_r)
-    nullify(ey_r)
-    nullify(ez_r)
-    nullify(bx_r)
-    nullify(by_r)
-    nullify(bz_r)
-    nullify(jx_r)
-    nullify(jy_r)
-    nullify(jz_r)
-    nullify(rho_r)
-    nullify(rhoold_r)
-    nullify(exf)
-    nullify(eyf)
-    nullify(ezf)
-    nullify(bxf)
-    nullify(byf)
-    nullify(bzf)
-    nullify(jxf)
-    nullify(jyf)
-    nullify(jzf)
-    nullify(rhof)
-    nullify(rhooldf)
-    call fftw_destroy_plan(plan_r2c_mpi)
-    call fftw_destroy_plan(plan_c2r_mpi)
-    call fftw_mpi_cleanup()
-  end subroutine warpx_fft_nullify
-
-
-  subroutine warpx_fft_push_eb ( &
-       ex_wrpx, exlo, exhi, &
-       ey_wrpx, eylo, eyhi, &
-       ez_wrpx, ezlo, ezhi, &
-       bx_wrpx, bxlo, bxhi, &
-       by_wrpx, bylo, byhi, &
-       bz_wrpx, bzlo, bzhi, &
-       jx_wrpx, jxlo, jxhi, &
-       jy_wrpx, jylo, jyhi, &
-       jz_wrpx, jzlo, jzhi, &
-       rhoold_wrpx, r1lo, r1hi, &
-       rho_wrpx, r2lo, r2hi) &
-       bind(c,name='warpx_fft_push_eb')
-
-    use fields, only: ex, ey, ez, bx, by, bz, jx, jy, jz
-    use shared_data, only: rhoold, rho
-    use constants, only: num
-
-    integer, dimension(3), intent(in) :: exlo, exhi, eylo, eyhi, ezlo, ezhi, bxlo, bxhi, &
-         bylo, byhi, bzlo, bzhi, jxlo, jxhi, jylo, jyhi, jzlo, jzhi, r1lo, r1hi, r2lo, r2hi
-    REAL(num), INTENT(INOUT), TARGET :: ex_wrpx(0:exhi(1)-exlo(1),0:exhi(2)-exlo(2),0:exhi(3)-exlo(3))
-    REAL(num), INTENT(INOUT), TARGET :: ey_wrpx(0:eyhi(1)-eylo(1),0:eyhi(2)-eylo(2),0:eyhi(3)-eylo(3))
-    REAL(num), INTENT(INOUT), TARGET :: ez_wrpx(0:ezhi(1)-ezlo(1),0:ezhi(2)-ezlo(2),0:ezhi(3)-ezlo(3))
-    REAL(num), INTENT(INOUT), TARGET :: bx_wrpx(0:bxhi(1)-bxlo(1),0:bxhi(2)-bxlo(2),0:bxhi(3)-bxlo(3))
-    REAL(num), INTENT(INOUT), TARGET :: by_wrpx(0:byhi(1)-bylo(1),0:byhi(2)-bylo(2),0:byhi(3)-bylo(3))
-    REAL(num), INTENT(INOUT), TARGET :: bz_wrpx(0:bzhi(1)-bzlo(1),0:bzhi(2)-bzlo(2),0:bzhi(3)-bzlo(3))
-    REAL(num), INTENT(INOUT), TARGET :: jx_wrpx(0:jxhi(1)-jxlo(1),0:jxhi(2)-jxlo(2),0:jxhi(3)-jxlo(3))
-    REAL(num), INTENT(INOUT), TARGET :: jy_wrpx(0:jyhi(1)-jylo(1),0:jyhi(2)-jylo(2),0:jyhi(3)-jylo(3))
-    REAL(num), INTENT(INOUT), TARGET :: jz_wrpx(0:jzhi(1)-jzlo(1),0:jzhi(2)-jzlo(2),0:jzhi(3)-jzlo(3))
-    REAL(num), INTENT(INOUT), TARGET :: rhoold_wrpx(0:r1hi(1)-r1lo(1),0:r1hi(2)-r1lo(2),0:r1hi(3)-r1lo(3))
-    REAL(num), INTENT(INOUT), TARGET :: rho_wrpx(0:r2hi(1)-r2lo(1),0:r2hi(2)-r2lo(2),0:r2hi(3)-r2lo(3))
-
-    ! Point the fields in the PICSAR modules to the fields provided by WarpX
-    ex => ex_wrpx
-    ey => ey_wrpx
-    ez => ez_wrpx
-    bx => bx_wrpx
-    by => by_wrpx
-    bz => bz_wrpx
-    jx => jx_wrpx
-    jy => jy_wrpx
-    jz => jz_wrpx
-    rho => rho_wrpx
-    rhoold => rhoold_wrpx
-
-    ! Call the corresponding PICSAR function
-    CALL push_psatd_ebfield()
-
-    ex => null()
-    ey => null()
-    ez => null()
-    bx => null()
-    by => null()
-    bz => null()
-    jx => null()
-    jy => null()
-    jz => null()
-    rho => null()
-    rhoold => null()
-  end subroutine warpx_fft_push_eb
-
-end module warpx_fft_module
diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp
index c1e0c59d0..c53070827 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp
@@ -37,7 +37,6 @@ SpectralKSpace::SpectralKSpace( const BoxArray& realspace_ba,
         // For local FFTs, boxes in spectral space start at 0 in
         // each direction and have the same number of points as the
         // (cell-centered) real space box
-        // TODO: this will be different for the hybrid FFT scheme
         Box realspace_bx = realspace_ba[i];
         IntVect fft_size = realspace_bx.length();
         // Because the spectral solver uses real-to-complex FFTs, we only
@@ -115,7 +114,6 @@ SpectralKSpace::getKComponent( const DistributionMapping& dm,
                 k[i] = (i-N)*dk;
             }
         }
-        // TODO: this will be different for the hybrid FFT scheme
     }
     return k_comp;
 }
diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp
index 6bf412b3c..09b1b7dbb 100644
--- a/Source/FieldSolver/WarpXPushFieldsEM.cpp
+++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp
@@ -95,13 +95,7 @@ WarpX::PushPSATD (amrex::Real a_dt)
 {
     for (int lev = 0; lev <= finest_level; ++lev) {
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(dt[lev] == a_dt, "dt must be consistent");
-        if (fft_hybrid_mpi_decomposition){
-#ifdef WARPX_USE_PSATD_HYBRID
-            PushPSATD_hybridFFT(lev, a_dt);
-#endif
-        } else {
-            PushPSATD_localFFT(lev, a_dt);
-        }
+        PushPSATD(lev, a_dt);
 
         // Evolve the fields in the PML boxes
         if (do_pml && pml[lev]->ok()) {
@@ -111,7 +105,7 @@ WarpX::PushPSATD (amrex::Real a_dt)
 }
 
 void
-WarpX::PushPSATD_localFFT (int lev, amrex::Real /* dt */)
+WarpX::PushPSATD (int lev, amrex::Real /* dt */)
 {
     // Update the fields on the fine and coarse patch
     PushPSATDSinglePatch( *spectral_solver_fp[lev],
diff --git a/Source/FortranInterface/Make.package b/Source/FortranInterface/Make.package
deleted file mode 100644
index a1301fb7d..000000000
--- a/Source/FortranInterface/Make.package
+++ /dev/null
@@ -1 +0,0 @@
-VPATH_LOCATIONS   += $(WARPX_HOME)/Source/FortranInterface
diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H
deleted file mode 100644
index a264780a3..000000000
--- a/Source/FortranInterface/WarpX_f.H
+++ /dev/null
@@ -1,62 +0,0 @@
-/* Copyright 2019 Andrew Myers, Aurore Blelly, Axel Huebl
- * David Grote, Jean-Luc Vay, Ligia Diana Amorim
- * Luca Fedeli, Mathieu Lobet, Maxence Thevenet
- * Remi Lehe, Revathi Jambunathan, Weiqun Zhang
- *
- *
- * This file is part of WarpX.
- *
- * License: BSD-3-Clause-LBNL
- */
-#ifndef WARPX_F_H_
-#define WARPX_F_H_
-
-#include <AMReX_BLFort.H>
-
-#ifdef __cplusplus
-
-#if AMREX_SPACEDIM==2
-#define WARPX_ARLIM_ANYD(x) std::array<int,3>{(x)[0], 0, (x)[1]}.data()
-#else
-#define WARPX_ARLIM_ANYD(x) x
-#endif
-
-#define WARPX_TO_FORTRAN_BOX(x) WARPX_ARLIM_ANYD((x).loVect()), WARPX_ARLIM_ANYD((x).hiVect())
-#define WARPX_TO_FORTRAN_ANYD(x) (x).dataPtr(), WARPX_ARLIM_ANYD((x).loVect()), WARPX_ARLIM_ANYD((x).hiVect())
-#define WARPX_TO_FORTRAN_N_ANYD(x,n) (x).dataPtr(n), WARPX_ARLIM_ANYD((x).loVect()), WARPX_ARLIM_ANYD((x).hiVect())
-
-#endif
-
-#ifdef __cplusplus
-extern "C"
-{
-#endif
-
-#ifdef WARPX_USE_PSATD
-    void warpx_fft_mpi_init (int fcomm);
-    void warpx_fft_domain_decomp (int* warpx_local_nz, int* warpx_local_z0,
-                                  const int* global_lo, const int* global_hi);
-    void warpx_fft_dataplan_init (const int* nox, const int* noy, const int* noz,
-                                  void* fft_data, const int* ndata,
-                                  const amrex_real* dx_w, const amrex_real* dt_w,
-                                  const int* fftw_plan_measure, const int* do_nodal );
-    void warpx_fft_nullify ();
-    void warpx_fft_push_eb (amrex_real* ex_w, const int* exlo, const int* exhi,
-                            amrex_real* ey_w, const int* eylo, const int* eyhi,
-                            amrex_real* ez_w, const int* ezlo, const int* ezhi,
-                            amrex_real* bx_w, const int* bxlo, const int* bxhi,
-                            amrex_real* by_w, const int* bylo, const int* byhi,
-                            amrex_real* bz_w, const int* bzlo, const int* bzhi,
-                            amrex_real* jx_w, const int* jxlo, const int* jxhi,
-                            amrex_real* jy_w, const int* jylo, const int* jyhi,
-                            amrex_real* jz_w, const int* jzlo, const int* jzhi,
-                            amrex_real* rhoold_w, const int* r1lo, const int* r1hi,
-                            amrex_real* rho_w, const int* r2lo, const int* r2hi);
-
-#endif
-
-#ifdef __cplusplus
-}
-#endif
-
-#endif //WARPX_F_H_
diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp
index 1f242f761..8f6dd9004 100644
--- a/Source/Initialization/WarpXInitData.cpp
+++ b/Source/Initialization/WarpXInitData.cpp
@@ -422,41 +422,6 @@ WarpX::InitLevelData (int lev, Real /*time*/)
     }
 }
 
-#ifdef WARPX_USE_PSATD_HYBRID
-
-void
-WarpX::InitLevelDataFFT (int lev, Real time)
-{
-
-    Efield_fp_fft[lev][0]->setVal(0.0);
-    Efield_fp_fft[lev][1]->setVal(0.0);
-    Efield_fp_fft[lev][2]->setVal(0.0);
-    Bfield_fp_fft[lev][0]->setVal(0.0);
-    Bfield_fp_fft[lev][1]->setVal(0.0);
-    Bfield_fp_fft[lev][2]->setVal(0.0);
-    current_fp_fft[lev][0]->setVal(0.0);
-    current_fp_fft[lev][1]->setVal(0.0);
-    current_fp_fft[lev][2]->setVal(0.0);
-    rho_fp_fft[lev]->setVal(0.0);
-
-    if (lev > 0)
-    {
-        Efield_cp_fft[lev][0]->setVal(0.0);
-        Efield_cp_fft[lev][1]->setVal(0.0);
-        Efield_cp_fft[lev][2]->setVal(0.0);
-        Bfield_cp_fft[lev][0]->setVal(0.0);
-        Bfield_cp_fft[lev][1]->setVal(0.0);
-        Bfield_cp_fft[lev][2]->setVal(0.0);
-        current_cp_fft[lev][0]->setVal(0.0);
-        current_cp_fft[lev][1]->setVal(0.0);
-        current_cp_fft[lev][2]->setVal(0.0);
-        rho_cp_fft[lev]->setVal(0.0);
-    }
-
-}
-
-#endif
-
 void
 WarpX::InitializeExternalFieldsOnGridUsingParser (
        MultiFab *mfx, MultiFab *mfy, MultiFab *mfz,
diff --git a/Source/Make.WarpX b/Source/Make.WarpX
index 72dce3480..226dc7500 100644
--- a/Source/Make.WarpX
+++ b/Source/Make.WarpX
@@ -63,7 +63,6 @@ include $(WARPX_HOME)/Source/BoundaryConditions/Make.package
 include $(WARPX_HOME)/Source/Diagnostics/Make.package
 include $(WARPX_HOME)/Source/FieldSolver/Make.package
 include $(WARPX_HOME)/Source/Filter/Make.package
-include $(WARPX_HOME)/Source/FortranInterface/Make.package
 include $(WARPX_HOME)/Source/Initialization/Make.package
 include $(WARPX_HOME)/Source/Laser/Make.package
 include $(WARPX_HOME)/Source/Parallelization/Make.package
@@ -166,11 +165,6 @@ ifeq ($(USE_PSATD),TRUE)
        INCLUDE_LOCATIONS += $(FFTW_HOME)/include
        LIBRARY_LOCATIONS += $(FFTW_HOME)/lib
      endif
-     ifeq ($(USE_PSATD_PICSAR),TRUE)
-        # Compile with PICSAR's Hybrid-MPI PSATD
-        DEFINES += -DWARPX_USE_PSATD_HYBRID
-        DEFINES += -DFFTW # PICSAR uses it
-     endif
   else
     # Use cuFFT
     libraries += -lcufft
diff --git a/Source/Parallelization/GuardCellManager.H b/Source/Parallelization/GuardCellManager.H
index 70948baa8..c2b9227d2 100644
--- a/Source/Parallelization/GuardCellManager.H
+++ b/Source/Parallelization/GuardCellManager.H
@@ -24,7 +24,6 @@ public:
      * \param do_fdtd_nci_corr bool, whether to use Godfrey NCI corrector
      * \param do_nodal bool, whether the field solver is nodal
      * \param do_moving_window bool, whether to use moving window
-     * \param do_fft_mpi_dec bool, whether to do parallel FFTs for PSATD
      * \param aux_is_nodal bool, true if the aux grid is nodal
      * \param moving_window_dir direction of moving window
      * \param nox order of current deposition
@@ -40,7 +39,6 @@ public:
         const bool do_fdtd_nci_corr,
         const bool do_nodal,
         const bool do_moving_window,
-        const bool do_fft_mpi_dec,
         const bool aux_is_nodal,
         const int moving_window_dir,
         const int nox,
diff --git a/Source/Parallelization/GuardCellManager.cpp b/Source/Parallelization/GuardCellManager.cpp
index f36c64214..c75b7426f 100644
--- a/Source/Parallelization/GuardCellManager.cpp
+++ b/Source/Parallelization/GuardCellManager.cpp
@@ -17,7 +17,6 @@ guardCellManager::Init(
     const bool do_fdtd_nci_corr,
     const bool do_nodal,
     const bool do_moving_window,
-    const bool do_fft_mpi_dec,
     const bool aux_is_nodal,
     const int moving_window_dir,
     const int nox,
@@ -95,40 +94,38 @@ guardCellManager::Init(
     ng_alloc_F = IntVect(AMREX_D_DECL(ng_alloc_F_int, ng_alloc_F_int, ng_alloc_F_int));
 
 #ifdef WARPX_USE_PSATD
-    if (do_fft_mpi_dec == false){
-        // All boxes should have the same number of guard cells
-        // (to avoid temporary parallel copies)
-        // Thus take the max of the required number of guards for each field
-        // Also: the number of guard cell should be enough to contain
-        // the stencil of the FFT solver. Here, this number (`ngFFT`)
-        // is determined *empirically* to be the order of the solver
-        // for nodal, and half the order of the solver for staggered.
-
-        int ngFFt_x = do_nodal ? nox_fft : nox_fft/2.;
-        int ngFFt_y = do_nodal ? noy_fft : noy_fft/2.;
-        int ngFFt_z = do_nodal ? noz_fft : noz_fft/2.;
-
-        ParmParse pp("psatd");
-        pp.query("nx_guard", ngFFt_x);
-        pp.query("ny_guard", ngFFt_y);
-        pp.query("nz_guard", ngFFt_z);
-
-        IntVect ngFFT = IntVect(AMREX_D_DECL(ngFFt_x, ngFFt_y, ngFFt_z));
-
-        for (int i_dim=0; i_dim<AMREX_SPACEDIM; i_dim++ ){
-            int ng_required = ngFFT[i_dim];
-            // Get the max
-            ng_required = std::max( ng_required, ng_alloc_EB[i_dim] );
-            ng_required = std::max( ng_required, ng_alloc_J[i_dim] );
-            ng_required = std::max( ng_required, ng_alloc_Rho[i_dim] );
-            ng_required = std::max( ng_required, ng_alloc_F[i_dim] );
-            // Set the guard cells to this max
-            ng_alloc_EB[i_dim] = ng_required;
-            ng_alloc_J[i_dim] = ng_required;
-            ng_alloc_F[i_dim] = ng_required;
-            ng_alloc_Rho[i_dim] = ng_required;
-            ng_alloc_F_int = ng_required;
-        }
+    // All boxes should have the same number of guard cells
+    // (to avoid temporary parallel copies)
+    // Thus take the max of the required number of guards for each field
+    // Also: the number of guard cell should be enough to contain
+    // the stencil of the FFT solver. Here, this number (`ngFFT`)
+    // is determined *empirically* to be the order of the solver
+    // for nodal, and half the order of the solver for staggered.
+
+    int ngFFt_x = do_nodal ? nox_fft : nox_fft/2.;
+    int ngFFt_y = do_nodal ? noy_fft : noy_fft/2.;
+    int ngFFt_z = do_nodal ? noz_fft : noz_fft/2.;
+
+    ParmParse pp("psatd");
+    pp.query("nx_guard", ngFFt_x);
+    pp.query("ny_guard", ngFFt_y);
+    pp.query("nz_guard", ngFFt_z);
+
+    IntVect ngFFT = IntVect(AMREX_D_DECL(ngFFt_x, ngFFt_y, ngFFt_z));
+
+    for (int i_dim=0; i_dim<AMREX_SPACEDIM; i_dim++ ){
+        int ng_required = ngFFT[i_dim];
+        // Get the max
+        ng_required = std::max( ng_required, ng_alloc_EB[i_dim] );
+        ng_required = std::max( ng_required, ng_alloc_J[i_dim] );
+        ng_required = std::max( ng_required, ng_alloc_Rho[i_dim] );
+        ng_required = std::max( ng_required, ng_alloc_F[i_dim] );
+        // Set the guard cells to this max
+        ng_alloc_EB[i_dim] = ng_required;
+        ng_alloc_J[i_dim] = ng_required;
+        ng_alloc_F[i_dim] = ng_required;
+        ng_alloc_Rho[i_dim] = ng_required;
+        ng_alloc_F_int = ng_required;
     }
     ng_alloc_F = IntVect(AMREX_D_DECL(ng_alloc_F_int, ng_alloc_F_int, ng_alloc_F_int));
 #endif
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index c63b94a27..a085c440d 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -11,7 +11,6 @@
 #include "PhysicalParticleContainer.H"
 
 #include "MultiParticleContainer.H"
-#include "FortranInterface/WarpX_f.H"
 #include "WarpX.H"
 #include "Utils/WarpXConst.H"
 #include "Utils/WarpXUtil.H"
diff --git a/Source/WarpX.H b/Source/WarpX.H
index d343528d7..1a344a785 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -31,9 +31,6 @@
 #       include "FieldSolver/SpectralSolver/SpectralSolver.H"
 #   endif
 #endif
-#ifdef WARPX_USE_PSATD_HYBRID
-#   include "FieldSolver/PicsarHybridSpectralSolver/PicsarHybridFFTData.H"
-#endif
 
 #include "Parallelization/GuardCellManager.H"
 
@@ -809,21 +806,6 @@ private:
     amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_slice;
     amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_slice;
 
-#ifdef WARPX_USE_PSATD_HYBRID
-    // Store fields in real space on the dual grid (i.e. the grid for the FFT push of the fields)
-    // This includes data for the FFT guard cells (between FFT groups)
-    amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_fp_fft;
-    amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_fp_fft;
-    amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_fp_fft;
-    amrex::Vector<            std::unique_ptr<amrex::MultiFab>      > rho_fp_fft;
-
-    amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_cp_fft;
-    amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_cp_fft;
-    amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_cp_fft;
-    amrex::Vector<            std::unique_ptr<amrex::MultiFab>      > rho_cp_fft;
-#endif
-
-    bool fft_hybrid_mpi_decomposition = false;
     bool fft_periodic_single_box = false;
     int nox_fft = 16;
     int noy_fft = 16;
@@ -833,9 +815,8 @@ private:
 private:
     void EvolvePSATD (int numsteps);
     void PushPSATD (amrex::Real dt);
-    void PushPSATD_localFFT (int lev, amrex::Real dt);
+    void PushPSATD (int lev, amrex::Real dt);
 
-    int ngroups_fft = 4;
     int fftw_plan_measure = 1;
 
 #   ifdef WARPX_DIM_RZ
@@ -849,33 +830,6 @@ private:
     amrex::Vector<std::unique_ptr<FiniteDifferenceSolver>> m_fdtd_solver_fp;
     amrex::Vector<std::unique_ptr<FiniteDifferenceSolver>> m_fdtd_solver_cp;
 
-#ifdef WARPX_USE_PSATD_HYBRID
-private:
-    amrex::Vector<std::unique_ptr<amrex::LayoutData<FFTData> > > dataptr_fp_fft;
-    amrex::Vector<std::unique_ptr<amrex::LayoutData<FFTData> > > dataptr_cp_fft;
-
-    // Domain decomposition containing the valid part of the dual grids (i.e. without FFT guard cells)
-    amrex::Vector<amrex::BoxArray> ba_valid_fp_fft;
-    amrex::Vector<amrex::BoxArray> ba_valid_cp_fft;
-
-    amrex::Vector<amrex::Box> domain_fp_fft;  //  "global" domain for the group this process belongs to
-    amrex::Vector<amrex::Box> domain_cp_fft;
-
-    amrex::Vector<MPI_Comm> comm_fft;
-    amrex::Vector<int>      color_fft;
-
-    void AllocLevelDataFFT (int lev);
-    void InitLevelDataFFT (int lev, amrex::Real time);
-
-    void InitFFTComm (int lev);
-    void FFTDomainDecomposition (int lev, amrex::BoxArray& ba_fft, amrex::DistributionMapping& dm_fft,
-                                 amrex::BoxArray& ba_valid, amrex::Box& domain_fft,
-                                 const amrex::Box& domain);
-    void InitFFTDataPlan (int lev);
-    void FreeFFT (int lev);
-    void PushPSATD_hybridFFT (int lev, amrex::Real dt);
-#endif
-
 #if defined(BL_USE_SENSEI_INSITU)
     amrex::AmrMeshInSituBridge *insitu_bridge;
 #endif
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 54deb59d4..1130a060f 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -303,29 +303,6 @@ WarpX::WarpX ()
 #endif
     m_fdtd_solver_fp.resize(nlevs_max);
     m_fdtd_solver_cp.resize(nlevs_max);
-#ifdef WARPX_USE_PSATD_HYBRID
-    Efield_fp_fft.resize(nlevs_max);
-    Bfield_fp_fft.resize(nlevs_max);
-    current_fp_fft.resize(nlevs_max);
-    rho_fp_fft.resize(nlevs_max);
-
-    Efield_cp_fft.resize(nlevs_max);
-    Bfield_cp_fft.resize(nlevs_max);
-    current_cp_fft.resize(nlevs_max);
-    rho_cp_fft.resize(nlevs_max);
-
-    dataptr_fp_fft.resize(nlevs_max);
-    dataptr_cp_fft.resize(nlevs_max);
-
-    ba_valid_fp_fft.resize(nlevs_max);
-    ba_valid_cp_fft.resize(nlevs_max);
-
-    domain_fp_fft.resize(nlevs_max);
-    domain_cp_fft.resize(nlevs_max);
-
-    comm_fft.resize(nlevs_max,MPI_COMM_NULL);
-    color_fft.resize(nlevs_max,-1);
-#endif
 
 #ifdef BL_USE_SENSEI_INSITU
     insitu_bridge = nullptr;
@@ -738,9 +715,7 @@ WarpX::ReadParameters ()
 #ifdef WARPX_USE_PSATD
     {
         ParmParse pp("psatd");
-        pp.query("hybrid_mpi_decomposition", fft_hybrid_mpi_decomposition);
         pp.query("periodic_single_box_fft", fft_periodic_single_box);
-        pp.query("ngroups_fft", ngroups_fft);
         pp.query("fftw_plan_measure", fftw_plan_measure);
         pp.query("nox", nox_fft);
         pp.query("noy", noy_fft);
@@ -809,17 +784,6 @@ WarpX::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& new_grids,
 {
     AllocLevelData(lev, new_grids, new_dmap);
     InitLevelData(lev, time);
-
-#ifdef WARPX_USE_PSATD
-    if (fft_hybrid_mpi_decomposition){
-#ifdef WARPX_USE_PSATD_HYBRID
-        AllocLevelDataFFT(lev);
-        InitLevelDataFFT(lev, time);
-#else
-    amrex::Abort("The option `psatd.fft_hybrid_mpi_decomposition` does not work on GPU.");
-#endif
-    }
-#endif
 }
 
 void
@@ -855,30 +819,6 @@ WarpX::ClearLevel (int lev)
     rho_cp[lev].reset();
 
     costs[lev].reset();
-
-
-#ifdef WARPX_USE_PSATD_HYBRID
-    for (int i = 0; i < 3; ++i) {
-        Efield_fp_fft[lev][i].reset();
-        Bfield_fp_fft[lev][i].reset();
-        current_fp_fft[lev][i].reset();
-
-        Efield_cp_fft[lev][i].reset();
-        Bfield_cp_fft[lev][i].reset();
-        current_cp_fft[lev][i].reset();
-    }
-
-    rho_fp_fft[lev].reset();
-    rho_cp_fft[lev].reset();
-
-    dataptr_fp_fft[lev].reset();
-    dataptr_cp_fft[lev].reset();
-
-    ba_valid_fp_fft[lev] = BoxArray();
-    ba_valid_cp_fft[lev] = BoxArray();
-
-    FreeFFT(lev);
-#endif
 }
 
 void
@@ -892,7 +832,6 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d
         WarpX::use_fdtd_nci_corr,
         do_nodal,
         do_moving_window,
-        fft_hybrid_mpi_decomposition,
         aux_is_nodal,
         moving_window_dir,
         WarpX::nox,
@@ -942,6 +881,8 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
     //
     // The fine patch
     //
+    std::array<Real,3> dx = CellSize(lev);
+
     Bfield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Bx_nodal_flag),dm,ncomps,ngE+ngextra));
     Bfield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,By_nodal_flag),dm,ncomps,ngE+ngextra));
     Bfield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Bz_nodal_flag),dm,ncomps,ngE+ngextra));
@@ -975,40 +916,35 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
     {
         rho_fp[lev].reset(new MultiFab(amrex::convert(ba,rho_nodal_flag),dm,2*ncomps,ngRho));
     }
-    if (fft_hybrid_mpi_decomposition == false){
-        // Allocate and initialize the spectral solver
-        std::array<Real,3> dx = CellSize(lev);
-#if (AMREX_SPACEDIM == 3)
-        RealVect dx_vect(dx[0], dx[1], dx[2]);
-#elif (AMREX_SPACEDIM == 2)
-        RealVect dx_vect(dx[0], dx[2]);
-#endif
-
-        // Check whether the option periodic, single box is valid here
-        if (fft_periodic_single_box) {
-            AMREX_ALWAYS_ASSERT_WITH_MESSAGE( geom[0].isAllPeriodic() && ba.size()==1 && lev==0,
-            "The option `psatd.periodic_single_box_fft` can only be used for a periodic domain, decomposed in a single box.");
-        }
-        // Get the cell-centered box
-        BoxArray realspace_ba = ba;  // Copy box
-        realspace_ba.enclosedCells(); // Make it cell-centered
-        // Define spectral solver
-#ifdef WARPX_DIM_RZ
-        realspace_ba.grow(1, ngE[1]); // add guard cells only in z
-        spectral_solver_fp[lev].reset( new SpectralSolverRZ( realspace_ba, dm,
-            n_rz_azimuthal_modes, noz_fft, do_nodal, dx_vect, dt[lev], lev ) );
-#else
-        if ( fft_periodic_single_box == false ) {
-            realspace_ba.grow(ngE); // add guard cells
-        }
-        bool const pml=false;
-        spectral_solver_fp[lev].reset( new SpectralSolver( realspace_ba, dm,
-            nox_fft, noy_fft, noz_fft, do_nodal, v_galilean, dx_vect, dt[lev],
-            pml, fft_periodic_single_box ) );
-#endif
+    // Allocate and initialize the spectral solver
+#   if (AMREX_SPACEDIM == 3)
+    RealVect dx_vect(dx[0], dx[1], dx[2]);
+#   elif (AMREX_SPACEDIM == 2)
+    RealVect dx_vect(dx[0], dx[2]);
+#   endif
+    // Check whether the option periodic, single box is valid here
+    if (fft_periodic_single_box) {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE( geom[0].isAllPeriodic() && ba.size()==1 && lev==0,
+        "The option `psatd.periodic_single_box_fft` can only be used for a periodic domain, decomposed in a single box.");
+    }
+    // Get the cell-centered box
+    BoxArray realspace_ba = ba;  // Copy box
+    realspace_ba.enclosedCells(); // Make it cell-centered
+    // Define spectral solver
+#   ifdef WARPX_DIM_RZ
+    realspace_ba.grow(1, ngE[1]); // add guard cells only in z
+    spectral_solver_fp[lev].reset( new SpectralSolverRZ( realspace_ba, dm,
+        n_rz_azimuthal_modes, noz_fft, do_nodal, dx_vect, dt[lev], lev ) );
+#   else
+    if ( fft_periodic_single_box == false ) {
+        realspace_ba.grow(ngE); // add guard cells
     }
+    bool const pml=false;
+    spectral_solver_fp[lev].reset( new SpectralSolver( realspace_ba, dm,
+        nox_fft, noy_fft, noz_fft, do_nodal, v_galilean, dx_vect, dt[lev],
+        pml, fft_periodic_single_box ) );
+#   endif
 #endif
-    std::array<Real,3> const dx = CellSize(lev);
     m_fdtd_solver_fp[lev].reset(
         new FiniteDifferenceSolver(maxwell_fdtd_solver_id, dx, do_nodal) );
     //
@@ -1051,6 +987,7 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
     {
         BoxArray cba = ba;
         cba.coarsen(refRatio(lev-1));
+        std::array<Real,3> cdx = CellSize(lev-1);
 
         // Create the MultiFabs for B
         Bfield_cp[lev][0].reset( new MultiFab(amrex::convert(cba,Bx_nodal_flag),dm,ncomps,ngE));
@@ -1079,33 +1016,28 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
         {
             rho_cp[lev].reset(new MultiFab(amrex::convert(cba,rho_nodal_flag),dm,2*ncomps,ngRho));
         }
-        if (fft_hybrid_mpi_decomposition == false){
-            // Allocate and initialize the spectral solver
-            std::array<Real,3> cdx = CellSize(lev-1);
-#if (AMREX_SPACEDIM == 3)
-            RealVect cdx_vect(cdx[0], cdx[1], cdx[2]);
-#elif (AMREX_SPACEDIM == 2)
-            RealVect cdx_vect(cdx[0], cdx[2]);
-#endif
-            // Get the cell-centered box, with guard cells
-            BoxArray realspace_ba = cba;// Copy box
-            realspace_ba.enclosedCells(); // Make it cell-centered
-            // Define spectral solver
-#ifdef WARPX_DIM_RZ
-            realspace_ba.grow(1, ngE[1]); // add guard cells only in z
-            spectral_solver_cp[lev].reset( new SpectralSolverRZ( realspace_ba, dm,
-                n_rz_azimuthal_modes, noz_fft, do_nodal, cdx_vect, dt[lev], lev ) );
-
-#else
-            realspace_ba.grow(ngE); // add guard cells
-            spectral_solver_cp[lev].reset( new SpectralSolver( realspace_ba, dm,
-                nox_fft, noy_fft, noz_fft, do_nodal, v_galilean, cdx_vect, dt[lev] ) );
-#endif
-        }
+        // Allocate and initialize the spectral solver
+#   if (AMREX_SPACEDIM == 3)
+        RealVect cdx_vect(cdx[0], cdx[1], cdx[2]);
+#   elif (AMREX_SPACEDIM == 2)
+        RealVect cdx_vect(cdx[0], cdx[2]);
+#   endif
+        // Get the cell-centered box, with guard cells
+        BoxArray realspace_ba = cba;// Copy box
+        realspace_ba.enclosedCells(); // Make it cell-centered
+        // Define spectral solver
+#   ifdef WARPX_DIM_RZ
+        realspace_ba.grow(1, ngE[1]); // add guard cells only in z
+        spectral_solver_cp[lev].reset( new SpectralSolverRZ( realspace_ba, dm,
+            n_rz_azimuthal_modes, noz_fft, do_nodal, cdx_vect, dt[lev], lev ) );
+#   else
+        realspace_ba.grow(ngE); // add guard cells
+        spectral_solver_cp[lev].reset( new SpectralSolver( realspace_ba, dm,
+            nox_fft, noy_fft, noz_fft, do_nodal, v_galilean, cdx_vect, dt[lev] ) );
+#   endif
 #endif
-    std::array<Real,3> cdx = CellSize(lev-1);
-    m_fdtd_solver_cp[lev].reset(
-        new FiniteDifferenceSolver( maxwell_fdtd_solver_id, cdx, do_nodal ) );
+        m_fdtd_solver_cp[lev].reset(
+            new FiniteDifferenceSolver( maxwell_fdtd_solver_id, cdx, do_nodal ) );
     }
 
     //