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

* Clean-up hybrid code

* Remove Fortran interface

7 files changed, 2 insertions, 845 deletions

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],