FFTOpenPoissonSolver.hpp

//
// Class FFTOpenPoissonSolver
//   FFT-based Poisson Solver for open boundaries.
//   Solves laplace(phi) = -rho, and E = -grad(phi).
//
//
// Communication helpers (pack / unpack / solver_send / solver_recv) now live in
// Field/FieldBufferOps.hpp and are included via FFTOpenPoissonSolver.h. All call
// sites below use the ippl::detail:: qualified names.
 
namespace ippl {
 
    using detail::pack;
    using detail::solver_recv;
    using detail::solver_send;
    using detail::unpack;
 
    // constructor and destructor
    template <typename FieldLHS, typename FieldRHS>
    FFTOpenPoissonSolver<FieldLHS, FieldRHS>::FFTOpenPoissonSolver()
        : Base()
        , mesh_mp(nullptr)
        , layout_mp(nullptr)
        , mesh2_m(nullptr)
        , layout2_m(nullptr)
        , meshComplex_m(nullptr)
        , layoutComplex_m(nullptr)
        , mesh4_m(nullptr)
        , layout4_m(nullptr)
        , mesh2n1_m(nullptr)
        , layout2n1_m(nullptr)
        , isGradFD_m(false) {
        setDefaultParameters();
    }
 
    template <typename FieldLHS, typename FieldRHS>
    FFTOpenPoissonSolver<FieldLHS, FieldRHS>::FFTOpenPoissonSolver(rhs_type& rhs,
                                                                   ParameterList& params)
        : mesh_mp(nullptr)
        , layout_mp(nullptr)
        , mesh2_m(nullptr)
        , layout2_m(nullptr)
        , meshComplex_m(nullptr)
        , layoutComplex_m(nullptr)
        , mesh4_m(nullptr)
        , layout4_m(nullptr)
        , mesh2n1_m(nullptr)
        , layout2n1_m(nullptr)
        , isGradFD_m(false) {
        using T = typename FieldLHS::value_type::value_type;
        static_assert(std::is_floating_point<T>::value, "Not a floating point type");
 
        setDefaultParameters();
        this->params_m.merge(params);
        this->params_m.update("output_type", Base::SOL);
 
        this->setRhs(rhs);
    }
 
    template <typename FieldLHS, typename FieldRHS>
    FFTOpenPoissonSolver<FieldLHS, FieldRHS>::FFTOpenPoissonSolver(lhs_type& lhs, rhs_type& rhs,
                                                                   ParameterList& params)
        : mesh_mp(nullptr)
        , layout_mp(nullptr)
        , mesh2_m(nullptr)
        , layout2_m(nullptr)
        , meshComplex_m(nullptr)
        , layoutComplex_m(nullptr)
        , mesh4_m(nullptr)
        , layout4_m(nullptr)
        , mesh2n1_m(nullptr)
        , layout2n1_m(nullptr)
        , isGradFD_m(false) {
        using T = typename FieldLHS::value_type::value_type;
        static_assert(std::is_floating_point<T>::value, "Not a floating point type");
 
        setDefaultParameters();
        this->params_m.merge(params);
 
        this->setLhs(lhs);
        this->setRhs(rhs);
    }
 
    // override setRhs to call class-specific initialization
    template <typename FieldLHS, typename FieldRHS>
    void FFTOpenPoissonSolver<FieldLHS, FieldRHS>::setRhs(rhs_type& rhs) {
        Base::setRhs(rhs);
 
        // start a timer
        static IpplTimings::TimerRef initialize = IpplTimings::getTimer("Initialize");
        IpplTimings::startTimer(initialize);
 
        initializeFields();
 
        IpplTimings::stopTimer(initialize);
    }
 
    // allows user to set gradient of phi = Efield instead of spectral
    // calculation of Efield (which uses FFTs)
 
    template <typename FieldLHS, typename FieldRHS>
    void FFTOpenPoissonSolver<FieldLHS, FieldRHS>::setGradFD() {
        // get the output type (sol, grad, or sol & grad)
        const int out = this->params_m.template get<int>("output_type");
 
        if (out != Base::SOL_AND_GRAD) {
            throw IpplException(
                "FFTOpenPoissonSolver::setGradFD()",
                "Cannot use gradient for Efield computation unless output type is SOL_AND_GRAD");
        } else {
            isGradFD_m = true;
        }
    }
 
    // initializeFields method, called in constructor
 
    template <typename FieldLHS, typename FieldRHS>
    void FFTOpenPoissonSolver<FieldLHS, FieldRHS>::initializeFields() {
        // get algorithm and hessian flag from parameter list
        const int alg                = this->params_m.template get<int>("algorithm");
        const int greensFunctionType = this->params_m.template get<int>("greens_function");
        const bool hessian           = this->params_m.template get<bool>("hessian");
 
        // first check if valid algorithm choice
        if ((alg != Algorithm::VICO) && (alg != Algorithm::HOCKNEY)
            && (alg != Algorithm::BIHARMONIC) && (alg != Algorithm::DCT_VICO)) {
            throw IpplException("FFTOpenPoissonSolver::initializeFields()",
                                "Currently only HOCKNEY, VICO, DCT_VICO, and BIHARMONIC are "
                                "supported for open BCs");
        }
 
        if ((greensFunctionType != GreenFunction::STANDARD)
            && (greensFunctionType != GreenFunction::INTEGRATED)) {
            throw IpplException("FFTOpenPoissonSolver::initializeFields()",
                                "Currently only STANDARD and INTEGRATED Green's functions are "
                                "supported for open BCs");
        }
 
        if ((greensFunctionType == GreenFunction::INTEGRATED)
            && ((alg != Algorithm::HOCKNEY) || (Dim != 3))) {
            throw IpplException("FFTOpenPoissonSolver::initializeFields()",
                                "The integrated Green's function is currently only implemented "
                                "for 3D HOCKNEY open BCs");
        }
 
        // check dimension
        if ((Dim == 2) && (alg != Algorithm::HOCKNEY)) {
            throw IpplException("FFTOpenPoissonSolver::initializeFields()",
                                "Only HOCKNEY supports 2D!");
        } else if (Dim == 1) {
            throw IpplException("FFTOpenPoissonSolver::initializeFields()",
                                "Only 2D and 3D solves supported!");
        }
 
        // get layout and mesh
        layout_mp              = &(this->rhs_mp->getLayout());
        mesh_mp                = &(this->rhs_mp->get_mesh());
        mpi::Communicator comm = layout_mp->comm;
 
        // get mesh spacing and origin
        hr_m               = mesh_mp->getMeshSpacing();
        vector_type origin = mesh_mp->getOrigin();
 
        // create domain for the real fields
        domain_m = layout_mp->getDomain();
 
        // get the mesh spacings and sizes for each dimension
        for (unsigned int i = 0; i < Dim; ++i) {
            nr_m[i] = domain_m[i].length();
 
            // create the doubled domain for the FFT procedure
            domain2_m[i] = Index(2 * nr_m[i]);
        }
 
        // define decomposition (parallel / serial)
        std::array<bool, Dim> isParallel = layout_mp->isParallel();
 
        // create double sized mesh and layout objects using the previously defined domain2_m
        using mesh_type = typename lhs_type::Mesh_t;
        mesh2_m         = std::unique_ptr<mesh_type>(new mesh_type(domain2_m, hr_m, origin));
        layout2_m = std::unique_ptr<FieldLayout_t>(new FieldLayout_t(comm, domain2_m, isParallel));
 
        // create the domain for the transformed (complex) fields
        // since we use HeFFTe for the transforms it doesn't require permuting to the right
        // one of the dimensions has only (n/2 +1) as our original fields are fully real
        // the dimension is given by the user via r2c_direction
        unsigned int RCDirection = this->params_m.template get<int>("r2c_direction");
        for (unsigned int i = 0; i < Dim; ++i) {
            if (i == RCDirection) {
                domainComplex_m[RCDirection] = Index(nr_m[RCDirection] + 1);
            } else {
                domainComplex_m[i] = Index(2 * nr_m[i]);
            }
        }
 
        // create mesh and layout for the real to complex FFT transformed fields
        meshComplex_m = std::unique_ptr<mesh_type>(new mesh_type(domainComplex_m, hr_m, origin));
        layoutComplex_m =
            std::unique_ptr<FieldLayout_t>(new FieldLayout_t(comm, domainComplex_m, isParallel));
 
        // initialize fields
        storage_field.initialize(*mesh2_m, *layout2_m);
        rho2tr_m.initialize(*meshComplex_m, *layoutComplex_m);
        grntr_m.initialize(*meshComplex_m, *layoutComplex_m);
 
        int out = this->params_m.template get<int>("output_type");
        if (((out == Base::GRAD || out == Base::SOL_AND_GRAD) && !isGradFD_m) || hessian) {
            temp_m.initialize(*meshComplex_m, *layoutComplex_m);
        }
 
        if (hessian) {
            hess_m.initialize(*mesh_mp, *layout_mp);
        }
 
        // create the FFT object
        fft_m = std::make_unique<FFT_t>(*layout2_m, *layoutComplex_m, this->params_m);
 
        // if Vico, also need to create mesh and layout for 4N Fourier domain
        // on this domain, the truncated Green's function is defined
        // also need to create the 4N complex grid, on which precomputation step done
        if (alg == Algorithm::VICO || alg == Algorithm::BIHARMONIC) {
            // start a timer
            static IpplTimings::TimerRef initialize_vico =
                IpplTimings::getTimer("Initialize: extra Vico");
            IpplTimings::startTimer(initialize_vico);
 
            for (unsigned int i = 0; i < Dim; ++i) {
                domain4_m[i] = Index(4 * nr_m[i]);
            }
 
            // 4N grid
            using mesh_type = typename lhs_type::Mesh_t;
            mesh4_m         = std::unique_ptr<mesh_type>(new mesh_type(domain4_m, hr_m, origin));
            layout4_m =
                std::unique_ptr<FieldLayout_t>(new FieldLayout_t(comm, domain4_m, isParallel));
 
            // initialize fields
            grnL_m.initialize(*mesh4_m, *layout4_m);
 
            // create a Complex-to-Complex FFT object to transform for layout4
            fft4n_m = std::make_unique<FFT<CCTransform, CxField_gt>>(*layout4_m, this->params_m);
 
            IpplTimings::stopTimer(initialize_vico);
        }
 
        // if DCT_VICO, need 2N+1 mesh, layout, domain, and green's function field for
        // precomputation
        if (alg == Algorithm::DCT_VICO) {
            // start a timer
            static IpplTimings::TimerRef initialize_vico =
                IpplTimings::getTimer("Initialize: extra Vico");
            IpplTimings::startTimer(initialize_vico);
 
            // 2N+1 domain for DCT
            for (unsigned int i = 0; i < Dim; ++i) {
                domain2n1_m[i] = Index(2 * nr_m[i] + 1);
            }
 
            // 2N+1 grid
            mesh2n1_m = std::unique_ptr<mesh_type>(new mesh_type(domain2n1_m, hr_m, origin));
            layout2n1_m =
                std::unique_ptr<FieldLayout_t>(new FieldLayout_t(comm, domain2n1_m, isParallel));
 
            // initialize fields
            grn2n1_m.initialize(*mesh2n1_m, *layout2n1_m);  // 2N+1 grnL
 
            // create real to real FFT object for 2N+1 grid
            fft2n1_m = std::make_unique<FFT<Cos1Transform, Field_t>>(*layout2n1_m, this->params_m);
 
            IpplTimings::stopTimer(initialize_vico);
        }
 
        // these are fields that are used for calculating the standard Green's function for Hockney
        if ((alg == Algorithm::HOCKNEY) && (greensFunctionType == GreenFunction::STANDARD)) {
            // start a timer
            static IpplTimings::TimerRef initialize_hockney =
                IpplTimings::getTimer("Initialize: extra Hockney");
            IpplTimings::startTimer(initialize_hockney);
 
            for (unsigned int gd = 0; gd < Dim; ++gd) {
                grnIField_m[gd].initialize(*mesh2_m, *layout2_m);
 
                // get number of ghost points and the Kokkos view to iterate over field
                auto view        = grnIField_m[gd].getView();
                const int nghost = grnIField_m[gd].getNghost();
                const auto& ldom = layout2_m->getLocalNDIndex();
 
                // the length of the physical domain
                const int size = nr_m[gd];
 
                // Kokkos parallel for loop to initialize grnIField[d]
                using index_array_type = typename RangePolicy<Dim>::index_array_type;
                ippl::parallel_for(
                    "Helper index Green field initialization",
                    grnIField_m[gd].getFieldRangePolicy(),
                    KOKKOS_LAMBDA(const index_array_type& args) {
                        scalar_type checkVal = 0.0;
 
                        // go from local indices to global
                        Vector<int, Dim> igVec = args - nghost;
                        for (unsigned d = 0; d < Dim; ++d) {
                            igVec[d] += ldom[d].first();
                            checkVal += igVec[d];
                        }
 
                        // assign (index)^2 if 0 <= index < N, and (2N-index)^2 elsewhere
                        const bool outsideN = (igVec[gd] >= size);
                        apply(view, args) =
                            (2 * size * outsideN - igVec[gd]) * (2 * size * outsideN - igVec[gd]);
 
                        // add 1.0 if at (0,0,0) to avoid singularity
                        const bool isOrig = (checkVal == 0);
                        apply(view, args) += isOrig * 1.0;
                    });
            }
            IpplTimings::stopTimer(initialize_hockney);
        }
 
        static IpplTimings::TimerRef warmup = IpplTimings::getTimer("Warmup");
        IpplTimings::startTimer(warmup);
 
        // "empty" transforms to warmup all the FFTs
        fft_m->warmup(rho2_mr, rho2tr_m);
        if (alg == Algorithm::VICO || alg == Algorithm::BIHARMONIC) {
            fft4n_m->warmup(grnL_m);
        }
        if (alg == Algorithm::DCT_VICO) {
            fft2n1_m->warmup(grn2n1_m);
        }
 
        IpplTimings::stopTimer(warmup);
 
        rho2_mr  = 0.0;
        rho2tr_m = 0.0;
        if (alg == Algorithm::VICO || alg == Algorithm::BIHARMONIC) {
            grnL_m = 0.0;
        }
        if (alg == Algorithm::DCT_VICO) {
            grn2n1_m = 0.0;
        }
 
        // call greensFunction and we will get the transformed G in the class attribute
        // this is done in initialization so that we already have the precomputed fct
        // for all timesteps (green's fct will only change if mesh size changes)
        static IpplTimings::TimerRef ginit = IpplTimings::getTimer("Green Init");
        IpplTimings::startTimer(ginit);
        greensFunction();
        IpplTimings::stopTimer(ginit);
    };
 
    // compute electric potential by solving Poisson's eq given a field rho and mesh spacings hr
    template <typename FieldLHS, typename FieldRHS>
    void FFTOpenPoissonSolver<FieldLHS, FieldRHS>::solve() {
        // start a timer
        static IpplTimings::TimerRef solve = IpplTimings::getTimer("Solve");
        IpplTimings::startTimer(solve);
 
        // get the output type (sol, grad, or sol & grad)
        const int out = this->params_m.template get<int>("output_type");
 
        // get the algorithm (hockney, vico, or biharmonic)
        const int alg = this->params_m.template get<int>("algorithm");
 
        // get hessian flag (if true, we compute the Hessian)
        const bool hessian = this->params_m.template get<bool>("hessian");
 
        // set the mesh & spacing, which may change each timestep
        mesh_mp = &(this->rhs_mp->get_mesh());
 
        // check whether the mesh spacing has changed with respect to the old one
        // if yes, update and set green flag to true
        bool green = false;
        for (unsigned int i = 0; i < Dim; ++i) {
            if (hr_m[i] != mesh_mp->getMeshSpacing(i)) {
                hr_m[i] = mesh_mp->getMeshSpacing(i);
                green   = true;
            }
        }
 
        // set mesh spacing on the other grids again
        mesh2_m->setMeshSpacing(hr_m);
        meshComplex_m->setMeshSpacing(hr_m);
 
        // field object on the doubled grid; zero-padded
        rho2_mr = 0.0;
 
        // start a timer
        static IpplTimings::TimerRef stod = IpplTimings::getTimer("Solve: Physical to double");
        IpplTimings::startTimer(stod);
 
        // store rho (RHS) in the lower left quadrant of the doubled grid
        // with or without communication (if only 1 rank)
 
        const int ranks = Comm->size();
 
        auto view2 = rho2_mr.getView();
        auto view1 = this->rhs_mp->getView();
 
        const int nghost2 = rho2_mr.getNghost();
        const int nghost1 = this->rhs_mp->getNghost();
 
        const auto& ldom2 = layout2_m->getLocalNDIndex();
        const auto& ldom1 = layout_mp->getLocalNDIndex();
 
        using index_array_type = typename RangePolicy<Dim>::index_array_type;
 
        if (ranks > 1) {
            // COMMUNICATION
            const auto& lDomains2 = layout2_m->getHostLocalDomains();
 
            // send
            std::vector<MPI_Request> requests(0);
 
            for (int i = 0; i < ranks; ++i) {
                if (lDomains2[i].touches(ldom1)) {
                    auto intersection = lDomains2[i].intersect(ldom1);
 
                    solver_send(mpi::tag::OPEN_SOLVER, 0, i, intersection, ldom1, nghost1, view1,
                                fd_m, requests);
                }
            }
 
            // receive
            const auto& lDomains1 = layout_mp->getHostLocalDomains();
 
            for (int i = 0; i < ranks; ++i) {
                if (lDomains1[i].touches(ldom2)) {
                    auto intersection = lDomains1[i].intersect(ldom2);
 
                    mpi::Communicator::size_type nrecvs;
                    nrecvs = intersection.size();
 
                    buffer_type buf = Comm->getBuffer<memory_space, Trhs>(nrecvs);
 
                    Comm->recv(i, mpi::tag::OPEN_SOLVER, fd_m, *buf, nrecvs * sizeof(Trhs), nrecvs);
                    buf->resetReadPos();
 
                    unpack(intersection, view2, fd_m, nghost2, ldom2);
                }
            }
 
            // wait for all messages to be received
            if (requests.size() > 0) {
                MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
            }
            ippl::Comm->freeAllBuffers();
        } else {
            ippl::parallel_for(
                "Write rho on the doubled grid", this->rhs_mp->getFieldRangePolicy(),
                KOKKOS_LAMBDA(const index_array_type& args) {
                    scalar_type checkVal = 0;
 
                    Vector<int, Dim> igVec1 = args - nghost1;
                    Vector<int, Dim> igVec2 = args - nghost2;
 
                    for (unsigned d = 0; d < Dim; ++d) {
                        igVec1[d] += ldom1[d].first();
                        igVec2[d] += ldom2[d].first();
 
                        checkVal += Kokkos::abs(igVec1[d] - igVec2[d]);
                    }
 
                    // Write physical rho on [0,N-1] of doubled field.
                    // Check whether we are in the 1st quadrant by checking whether
                    // the global indices for view1 and view2 are equal.
                    // This is done using checkVal, which should be 0 if ig1 = ig2.
                    const bool isQuadrant1 = (checkVal == 0);
                    apply(view2, args)     = apply(view1, args) * isQuadrant1;
                });
        }
 
        IpplTimings::stopTimer(stod);
 
        // start a timer
        static IpplTimings::TimerRef fftrho = IpplTimings::getTimer("FFT: Rho");
        IpplTimings::startTimer(fftrho);
 
        // forward FFT of the charge density field on doubled grid
        fft_m->transform(FORWARD, rho2_mr, rho2tr_m);
 
        IpplTimings::stopTimer(fftrho);
 
        // call greensFunction to recompute if the mesh spacing has changed
        if (green) {
            greensFunction();
        }
 
        // multiply FFT(rho2)*FFT(green)
        // convolution becomes multiplication in FFT
        // minus sign since we are solving laplace(phi) = -rho
        rho2tr_m = -rho2tr_m * grntr_m;
 
        // if output_type is SOL or SOL_AND_GRAD, we caculate solution
        if ((out == Base::SOL) || (out == Base::SOL_AND_GRAD)) {
            // start a timer
            static IpplTimings::TimerRef fftc = IpplTimings::getTimer("FFT: Convolution");
            IpplTimings::startTimer(fftc);
 
            // inverse FFT of the product and store the electrostatic potential in rho2_mr
            fft_m->transform(BACKWARD, rho2_mr, rho2tr_m);
 
            IpplTimings::stopTimer(fftc);
            // Hockney: multiply the rho2_mr field by the total number of points to account for
            // double counting (rho and green) of normalization factor in forward transform
            // also multiply by the mesh spacing^3 (to account for discretization)
            // Vico: need to multiply by normalization factor of 1/4N^3,
            // since only backward transform was performed on the 4N grid
            // DCT_VICO: need to multiply by a factor of (2N)^3 to match the normalization factor in
            // the transform.
            for (unsigned int i = 0; i < Dim; ++i) {
                switch (alg) {
                    case Algorithm::HOCKNEY:
                        rho2_mr = rho2_mr * 2.0 * nr_m[i] * hr_m[i];
                        break;
                    case Algorithm::VICO:
                    case Algorithm::BIHARMONIC:
                        rho2_mr = rho2_mr * 2.0 * (1.0 / 4.0);
                        break;
                    case Algorithm::DCT_VICO:
                        rho2_mr = rho2_mr * (1.0 / 4.0);
                        break;
                    default:
                        throw IpplException(
                            "FFTOpenPoissonSolver::initializeFields()",
                            "Currently only HOCKNEY, VICO, DCT_VICO, and BIHARMONIC are "
                            "supported for open BCs");
                }
            }
 
            // start a timer
            static IpplTimings::TimerRef dtos = IpplTimings::getTimer("Solve: Double to physical");
            IpplTimings::startTimer(dtos);
 
            // get the physical part only --> physical electrostatic potential is now given in RHS
            // need communication if more than one rank
 
            if (ranks > 1) {
                // COMMUNICATION
 
                // send
                const auto& lDomains1 = layout_mp->getHostLocalDomains();
 
                std::vector<MPI_Request> requests(0);
 
                for (int i = 0; i < ranks; ++i) {
                    if (lDomains1[i].touches(ldom2)) {
                        auto intersection = lDomains1[i].intersect(ldom2);
 
                        solver_send(mpi::tag::OPEN_SOLVER, 0, i, intersection, ldom2, nghost2,
                                    view2, fd_m, requests);
                    }
                }
 
                // receive
                const auto& lDomains2 = layout2_m->getHostLocalDomains();
 
                for (int i = 0; i < ranks; ++i) {
                    if (ldom1.touches(lDomains2[i])) {
                        auto intersection = ldom1.intersect(lDomains2[i]);
 
                        mpi::Communicator::size_type nrecvs;
                        nrecvs = intersection.size();
 
                        buffer_type buf = Comm->getBuffer<memory_space, Trhs>(nrecvs);
 
                        Comm->recv(i, mpi::tag::OPEN_SOLVER, fd_m, *buf, nrecvs * sizeof(Trhs),
                                   nrecvs);
                        buf->resetReadPos();
 
                        unpack(intersection, view1, fd_m, nghost1, ldom1);
                    }
                }
 
                // wait for all messages to be received
                if (requests.size() > 0) {
                    MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
                }
                ippl::Comm->freeAllBuffers();
 
            } else {
                ippl::parallel_for(
                    "Write the solution into the LHS on physical grid",
                    this->rhs_mp->getFieldRangePolicy(),
                    KOKKOS_LAMBDA(const index_array_type& args) {
                        scalar_type checkVal = 0;
 
                        Vector<int, Dim> igVec1 = args - nghost1;
                        Vector<int, Dim> igVec2 = args - nghost2;
 
                        for (unsigned d = 0; d < Dim; ++d) {
                            igVec1[d] += ldom1[d].first();
                            igVec2[d] += ldom2[d].first();
 
                            checkVal += Kokkos::abs(igVec1[d] - igVec2[d]);
                        }
 
                        // Take [0,N-1] quadrant as physical solution.
                        // Check whether we are in the 1st quadrant by checking whether
                        // the global indices for view1 and view2 are equal.
                        // This is done using checkVal, which should be 0 if ig1 = ig2.
                        const bool isQuadrant1 = (checkVal == 0);
                        apply(view1, args)     = apply(view2, args) * isQuadrant1;
                    });
            }
            IpplTimings::stopTimer(dtos);
        }
 
        // if we want finite differences Efield = -grad(phi)
        // instead of computing it in Fourier domain
        // this is only possible if SOL_AND_GRAD is the output type
        if (isGradFD_m && (out == Base::SOL_AND_GRAD)) {
            *(this->lhs_mp) = -grad(*this->rhs_mp);
        }
 
        // if output_type is GRAD or SOL_AND_GRAD, we calculate E-field (gradient in Fourier domain)
        if (((out == Base::GRAD) || (out == Base::SOL_AND_GRAD)) && (!isGradFD_m)) {
            // start a timer
            static IpplTimings::TimerRef efield = IpplTimings::getTimer("Solve: Electric field");
            IpplTimings::startTimer(efield);
 
            // get E field view (LHS)
            auto viewL        = this->lhs_mp->getView();
            const int nghostL = this->lhs_mp->getNghost();
 
            // get rho2tr_m view (as we want to multiply by ik then transform)
            auto viewR        = rho2tr_m.getView();
            const int nghostR = rho2tr_m.getNghost();
            const auto& ldomR = layoutComplex_m->getLocalNDIndex();
 
            // use temp_m as a temporary complex field
            auto view_g = temp_m.getView();
 
            // define some constants
            const scalar_type pi          = Kokkos::numbers::pi_v<scalar_type>;
            const Kokkos::complex<Trhs> I = {0.0, 1.0};
 
            // define some member variables in local scope for the parallel_for
            vector_type hsize  = hr_m;
            Vector<int, Dim> N = nr_m;
 
            // loop over each component (E = vector field)
            for (size_t gd = 0; gd < Dim; ++gd) {
                // loop over rho2tr_m to multiply by -ik (gradient in Fourier space)
                ippl::parallel_for(
                    "Gradient - E field", rho2tr_m.getFieldRangePolicy(),
                    KOKKOS_LAMBDA(const index_array_type& args) {
                        // global indices for 2N rhotr_m
                        Vector<int, Dim> igVec = args - nghostR;
                        for (unsigned d = 0; d < Dim; ++d) {
                            igVec[d] += ldomR[d].first();
                        }
 
                        scalar_type k_gd;
                        const scalar_type Len = N[gd] * hsize[gd];
                        const bool shift      = (igVec[gd] > N[gd]);
                        const bool notMid     = (igVec[gd] != N[gd]);
 
                        k_gd = notMid * (pi / Len) * (igVec[gd] - shift * 2 * N[gd]);
 
                        apply(view_g, args) = -(I * k_gd) * apply(viewR, args);
                    });
 
                // start a timer
                static IpplTimings::TimerRef ffte = IpplTimings::getTimer("FFT: Efield");
                IpplTimings::startTimer(ffte);
 
                // transform to get E-field
                fft_m->transform(BACKWARD, rho2_mr, temp_m);
 
                IpplTimings::stopTimer(ffte);
 
                // apply proper normalization
                for (unsigned int i = 0; i < Dim; ++i) {
                    switch (alg) {
                        case Algorithm::HOCKNEY:
                            rho2_mr = rho2_mr * 2.0 * nr_m[i] * hr_m[i];
                            break;
                        case Algorithm::VICO:
                        case Algorithm::BIHARMONIC:
                            rho2_mr = rho2_mr * 2.0 * (1.0 / 4.0);
                            break;
                        case Algorithm::DCT_VICO:
                            rho2_mr = rho2_mr * (1.0 / 4.0);
                            break;
                        default:
                            throw IpplException(
                                "FFTOpenPoissonSolver::initializeFields()",
                                "Currently only HOCKNEY, VICO, DCT_VICO, and BIHARMONIC are "
                                "supported for open BCs");
                    }
                }
 
                // start a timer
                static IpplTimings::TimerRef edtos =
                    IpplTimings::getTimer("Efield: double to phys.");
                IpplTimings::startTimer(edtos);
 
                // restrict to physical grid (N^3) and assign to LHS (E-field)
                // communication needed if more than one rank
                if (ranks > 1) {
                    // COMMUNICATION
 
                    // send
                    const auto& lDomains1 = layout_mp->getHostLocalDomains();
                    std::vector<MPI_Request> requests(0);
 
                    for (int i = 0; i < ranks; ++i) {
                        if (lDomains1[i].touches(ldom2)) {
                            auto intersection = lDomains1[i].intersect(ldom2);
 
                            solver_send(mpi::tag::OPEN_SOLVER, 0, i, intersection, ldom2, nghost2,
                                        view2, fd_m, requests);
                        }
                    }
 
                    // receive
                    const auto& lDomains2 = layout2_m->getHostLocalDomains();
 
                    for (int i = 0; i < ranks; ++i) {
                        if (ldom1.touches(lDomains2[i])) {
                            auto intersection = ldom1.intersect(lDomains2[i]);
 
                            mpi::Communicator::size_type nrecvs;
                            nrecvs = intersection.size();
 
                            buffer_type buf = Comm->getBuffer<memory_space, Trhs>(nrecvs);
 
                            Comm->recv(i, mpi::tag::OPEN_SOLVER, fd_m, *buf, nrecvs * sizeof(Trhs),
                                       nrecvs);
                            buf->resetReadPos();
 
                            unpack(intersection, viewL, gd, fd_m, nghostL, ldom1);
                        }
                    }
 
                    // wait for all messages to be received
                    if (requests.size() > 0) {
                        MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
                    }
                    ippl::Comm->freeAllBuffers();
 
                } else {
                    ippl::parallel_for(
                        "Write the E-field on physical grid", this->lhs_mp->getFieldRangePolicy(),
                        KOKKOS_LAMBDA(const index_array_type& args) {
                            scalar_type checkVal = 0;
 
                            Vector<int, Dim> igVec1 = args - nghostL;
                            Vector<int, Dim> igVec2 = args - nghost2;
 
                            for (unsigned d = 0; d < Dim; ++d) {
                                igVec1[d] += ldom1[d].first();
                                igVec2[d] += ldom2[d].first();
 
                                checkVal += Kokkos::abs(igVec1[d] - igVec2[d]);
                            }
 
                            // Take [0,N-1] quadrant as physical solution.
                            // Check whether we are in the 1st quadrant by checking whether
                            // the global indices for view1 and view2 are equal.
                            // This is done using checkVal, which should be 0 if ig1 = ig2.
                            const bool isQuadrant1 = (checkVal == 0);
                            apply(viewL, args)[gd] = apply(view2, args) * isQuadrant1;
                        });
                }
                IpplTimings::stopTimer(edtos);
            }
            IpplTimings::stopTimer(efield);
        }
 
        // if user asked for Hessian, compute in Fourier domain (-k^2 multiplication)
        if (hessian) {
            // start a timer
            static IpplTimings::TimerRef hess = IpplTimings::getTimer("Solve: Hessian");
            IpplTimings::startTimer(hess);
 
            // get Hessian matrix view (LHS)
            auto viewH        = hess_m.getView();
            const int nghostH = hess_m.getNghost();
 
            // get rho2tr_m view (as we want to multiply by -k^2 then transform)
            auto viewR        = rho2tr_m.getView();
            const int nghostR = rho2tr_m.getNghost();
            const auto& ldomR = layoutComplex_m->getLocalNDIndex();
 
            // use temp_m as a temporary complex field
            auto view_g = temp_m.getView();
 
            // define some constants
            const scalar_type pi = Kokkos::numbers::pi_v<scalar_type>;
 
            // define some member variables in local scope for the parallel_for
            vector_type hsize  = hr_m;
            Vector<int, Dim> N = nr_m;
 
            // loop over each component (Hessian = Matrix field)
            for (size_t row = 0; row < Dim; ++row) {
                for (size_t col = 0; col < Dim; ++col) {
                    // loop over rho2tr_m to multiply by -k^2 (second derivative in Fourier space)
                    // if diagonal element (row = col), do not need N/2 term = 0
                    // else, if mixed derivative, need kVec = 0 at N/2
 
                    ippl::parallel_for(
                        "Hessian", rho2tr_m.getFieldRangePolicy(),
                        KOKKOS_LAMBDA(const index_array_type& args) {
                            // global indices for 2N rhotr_m
                            Vector<int, Dim> igVec = args - nghostR;
                            for (unsigned d = 0; d < Dim; ++d) {
                                igVec[d] += ldomR[d].first();
                            }
 
                            // compute wave-vector
                            Vector_t kVec;
                            for (size_t d = 0; d < Dim; ++d) {
                                const scalar_type Len = N[d] * hsize[d];
                                const bool shift      = (igVec[d] > N[d]);
                                const bool isMid      = (igVec[d] == N[d]);
                                const bool notDiag    = (row != col);
 
                                kVec[d] = (1 - (notDiag * isMid)) * (pi / Len)
                                          * (igVec[d] - shift * 2 * N[d]);
                            }
 
                            apply(view_g, args) = -(kVec[col] * kVec[row]) * apply(viewR, args);
                        });
 
                    // start a timer
                    static IpplTimings::TimerRef ffth = IpplTimings::getTimer("FFT: Hessian");
                    IpplTimings::startTimer(ffth);
 
                    // transform to get Hessian
                    fft_m->transform(BACKWARD, rho2_mr, temp_m);
 
                    IpplTimings::stopTimer(ffth);
 
                    // apply proper normalization
                    for (unsigned int i = 0; i < Dim; ++i) {
                        switch (alg) {
                            case Algorithm::HOCKNEY:
                                rho2_mr = rho2_mr * 2.0 * nr_m[i] * hr_m[i];
                                break;
                            case Algorithm::VICO:
                            case Algorithm::BIHARMONIC:
                                rho2_mr = rho2_mr * 2.0 * (1.0 / 4.0);
                                break;
                            case Algorithm::DCT_VICO:
                                rho2_mr = rho2_mr * (1.0 / 4.0);
                                break;
                            default:
                                throw IpplException(
                                    "FFTOpenPoissonSolver::initializeFields()",
                                    "Currently only HOCKNEY, VICO, DCT_VICO, and BIHARMONIC are "
                                    "supported for open BCs");
                        }
                    }
 
                    // start a timer
                    static IpplTimings::TimerRef hdtos =
                        IpplTimings::getTimer("Hessian: double to phys.");
                    IpplTimings::startTimer(hdtos);
 
                    // restrict to physical grid (N^3) and assign to Matrix field (Hessian)
                    // communication needed if more than one rank
                    if (ranks > 1) {
                        // COMMUNICATION
 
                        // send
                        const auto& lDomains1 = layout_mp->getHostLocalDomains();
                        std::vector<MPI_Request> requests(0);
 
                        for (int i = 0; i < ranks; ++i) {
                            if (lDomains1[i].touches(ldom2)) {
                                auto intersection = lDomains1[i].intersect(ldom2);
 
                                solver_send(mpi::tag::OPEN_SOLVER, 0, i, intersection, ldom2,
                                            nghost2, view2, fd_m, requests);
                            }
                        }
 
                        // receive
                        const auto& lDomains2 = layout2_m->getHostLocalDomains();
 
                        for (int i = 0; i < ranks; ++i) {
                            if (ldom1.touches(lDomains2[i])) {
                                auto intersection = ldom1.intersect(lDomains2[i]);
 
                                mpi::Communicator::size_type nrecvs;
                                nrecvs = intersection.size();
 
                                buffer_type buf = Comm->getBuffer<memory_space, Trhs>(nrecvs);
 
                                Comm->recv(i, mpi::tag::OPEN_SOLVER, fd_m, *buf,
                                           nrecvs * sizeof(Trhs), nrecvs);
                                buf->resetReadPos();
 
                                unpack(intersection, viewH, fd_m, nghostH, ldom1, row, col);
                            }
                        }
 
                        // wait for all messages to be received
                        if (requests.size() > 0) {
                            MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
                        }
                        ippl::Comm->freeAllBuffers();
 
                    } else {
                        ippl::parallel_for(
                            "Write Hessian on physical grid", hess_m.getFieldRangePolicy(),
                            KOKKOS_LAMBDA(const index_array_type& args) {
                                scalar_type checkVal = 0;
 
                                Vector<int, Dim> igVec1 = args - nghostH;
                                Vector<int, Dim> igVec2 = args - nghost2;
 
                                for (unsigned d = 0; d < Dim; ++d) {
                                    igVec1[d] += ldom1[d].first();
                                    igVec2[d] += ldom2[d].first();
 
                                    checkVal += Kokkos::abs(igVec1[d] - igVec2[d]);
                                }
 
                                // Take [0,N-1] quadrant as physical solution.
                                // Check whether we are in the 1st quadrant by checking whether
                                // the global indices for view1 and view2 are equal.
                                // This is done using checkVal, which should be 0 if ig1 = ig2.
                                const bool isQuadrant1       = (checkVal == 0);
                                apply(viewH, args)[row][col] = apply(view2, args) * isQuadrant1;
                            });
                    }
                    IpplTimings::stopTimer(hdtos);
                }
            }
            IpplTimings::stopTimer(hess);
        }
        IpplTimings::stopTimer(solve);
    };
 
    // calculate FFT of the Green's function
 
    template <typename FieldLHS, typename FieldRHS>
    KOKKOS_INLINE_FUNCTION typename FFTOpenPoissonSolver<FieldLHS, FieldRHS>::scalar_type
    FFTOpenPoissonSolver<FieldLHS, FieldRHS>::integratedGreenAtanTerm(
        const scalar_type& numerator, const scalar_type& denominator) {
        const scalar_type pi = Kokkos::numbers::pi_v<scalar_type>;
        if (denominator == scalar_type(0)) {
            if (numerator > scalar_type(0)) {
                return scalar_type(0.5) * pi;
            }
            if (numerator < scalar_type(0)) {
                return -scalar_type(0.5) * pi;
            }
            return scalar_type(0);
        }
        return Kokkos::atan(numerator / denominator);
    }
 
    template <typename FieldLHS, typename FieldRHS>
    KOKKOS_INLINE_FUNCTION typename FFTOpenPoissonSolver<FieldLHS, FieldRHS>::scalar_type
    FFTOpenPoissonSolver<FieldLHS, FieldRHS>::integratedGreenLogTerm(const scalar_type& a,
                                                                     const scalar_type& b,
                                                                     const scalar_type& c,
                                                                     const scalar_type& r) {
        const scalar_type coeff = b * c;
        return (coeff == scalar_type(0)) ? scalar_type(0) : coeff * Kokkos::log(a + r);
    }
 
    template <typename FieldLHS, typename FieldRHS>
    KOKKOS_INLINE_FUNCTION typename FFTOpenPoissonSolver<FieldLHS, FieldRHS>::scalar_type
    FFTOpenPoissonSolver<FieldLHS, FieldRHS>::integratedGreenAntiderivative(const scalar_type& x,
                                                                            const scalar_type& y,
                                                                            const scalar_type& z) {
        const scalar_type r2 = x * x + y * y + z * z;
        if (r2 == scalar_type(0)) {
            return scalar_type(0);
        }
 
        const scalar_type r = Kokkos::sqrt(r2);
        scalar_type value   = scalar_type(0);
        value -= scalar_type(0.5) * z * z * integratedGreenAtanTerm(x * y, z * r);
        value -= scalar_type(0.5) * y * y * integratedGreenAtanTerm(x * z, y * r);
        value -= scalar_type(0.5) * x * x * integratedGreenAtanTerm(y * z, x * r);
        value += integratedGreenLogTerm(x, y, z, r);
        value += integratedGreenLogTerm(y, x, z, r);
        value += integratedGreenLogTerm(z, x, y, r);
        return value;
    }
 
    template <typename FieldLHS, typename FieldRHS>
    KOKKOS_INLINE_FUNCTION typename FFTOpenPoissonSolver<FieldLHS, FieldRHS>::scalar_type
    FFTOpenPoissonSolver<FieldLHS, FieldRHS>::integratedGreenAverage(const scalar_type& x,
                                                                     const scalar_type& y,
                                                                     const scalar_type& z,
                                                                     const vector_type& h) {
        const scalar_type half = scalar_type(0.5);
        const scalar_type x0   = x - half * h[0];
        const scalar_type x1   = x + half * h[0];
        const scalar_type y0   = y - half * h[1];
        const scalar_type y1   = y + half * h[1];
        const scalar_type z0   = z - half * h[2];
        const scalar_type z1   = z + half * h[2];
 
        scalar_type integral = integratedGreenAntiderivative(x1, y1, z1);
        integral -= integratedGreenAntiderivative(x0, y1, z1);
        integral -= integratedGreenAntiderivative(x1, y0, z1);
        integral += integratedGreenAntiderivative(x0, y0, z1);
        integral -= integratedGreenAntiderivative(x1, y1, z0);
        integral += integratedGreenAntiderivative(x0, y1, z0);
        integral += integratedGreenAntiderivative(x1, y0, z0);
        integral -= integratedGreenAntiderivative(x0, y0, z0);
 
        return integral / (h[0] * h[1] * h[2]);
    }
 
    template <typename FieldLHS, typename FieldRHS>
    void FFTOpenPoissonSolver<FieldLHS, FieldRHS>::greensFunction() {
        const scalar_type pi = Kokkos::numbers::pi_v<scalar_type>;
        grn_mr               = 0.0;
 
        const int alg                = this->params_m.template get<int>("algorithm");
        const int greensFunctionType = this->params_m.template get<int>("greens_function");
 
        using index_array_type = typename RangePolicy<Dim>::index_array_type;
 
        if (alg == Algorithm::VICO || alg == Algorithm::BIHARMONIC) {
            if constexpr (Dim == 3) {
                Vector_t l(hr_m * nr_m);
                Vector_t hs_m;
                double L_sum(0.0);
 
                // compute length of the physical domain
                // compute Fourier domain spacing
                for (unsigned int i = 0; i < Dim; ++i) {
                    hs_m[i] = pi * 0.5 / l[i];
                    L_sum   = L_sum + l[i] * l[i];
                }
 
                // define the origin of the 4N grid
                vector_type origin;
 
                for (unsigned int i = 0; i < Dim; ++i) {
                    origin[i] = -2 * nr_m[i] * pi / l[i];
                }
 
                // set mesh for the 4N mesh
                mesh4_m->setMeshSpacing(hs_m);
 
                // size of truncation window
                L_sum = Kokkos::sqrt(L_sum);
                // we choose a window 10% larger than domain (arbitrary choice)
                L_sum = 1.1 * L_sum;
 
                // initialize grnL_m
                typename CxField_gt::view_type view_g = grnL_m.getView();
                const int nghost_g                    = grnL_m.getNghost();
                const auto& ldom_g                    = layout4_m->getLocalNDIndex();
 
                Vector<int, Dim> size = nr_m;
 
                // Kokkos parallel for loop to assign analytic grnL_m
                if (alg == Algorithm::VICO) {
                    ippl::parallel_for(
                        "Initialize Green's function ", grnL_m.getFieldRangePolicy(),
                        KOKKOS_LAMBDA(const index_array_type& args) {
                            scalar_type checkVal = 0;
 
                            // go from local indices to global
                            Vector<int, Dim> igVec = args - nghost_g;
                            for (unsigned d = 0; d < Dim; ++d) {
                                igVec[d] += ldom_g[d].first();
                                checkVal += igVec[d];
                            }
 
                            // compute s
                            Tg s = 0;
                            for (unsigned d = 0; d < Dim; ++d) {
                                bool isOutside = (igVec[d] > 2 * size[d] - 1);
                                const Tg t     = igVec[d] * hs_m[d] + isOutside * origin[d];
                                s += (t * t);
                            }
                            s = Kokkos::sqrt(s);
 
                            // assign the green's function value
                            // if (0,0,0), assign L^2/2 (analytical limit of sinc)
                            const bool isOrig    = (checkVal == 0);
                            const Tg analyticLim = -L_sum * L_sum * 0.5;
                            const Tg value       = -2.0
                                             * (Kokkos::sin(0.5 * L_sum * s) / (s + isOrig * 1.0))
                                             * (Kokkos::sin(0.5 * L_sum * s) / (s + isOrig * 1.0));
 
                            apply(view_g, args) = (!isOrig) * value + isOrig * analyticLim;
                        });
                } else if (alg == Algorithm::BIHARMONIC) {
                    ippl::parallel_for(
                        "Initialize Green's function ", grnL_m.getFieldRangePolicy(),
                        KOKKOS_LAMBDA(const index_array_type& args) {
                            scalar_type checkVal = 0;
 
                            // go from local indices to global
                            Vector<int, Dim> igVec = args - nghost_g;
                            for (unsigned d = 0; d < Dim; ++d) {
                                igVec[d] += ldom_g[d].first();
                                checkVal += igVec[d];
                            }
 
                            // compute s
                            Tg s = 0;
                            for (unsigned d = 0; d < Dim; ++d) {
                                bool isOutside = (igVec[d] > 2 * size[d] - 1);
                                const Tg t     = igVec[d] * hs_m[d] + isOutside * origin[d];
                                s += (t * t);
                            }
                            s = Kokkos::sqrt(s);
 
                            // assign the green's function value
                            // if (0,0,0), assign L^2/2 (analytical limit of sinc)
                            const bool isOrig    = (checkVal == 0);
                            const Tg analyticLim = -L_sum * L_sum * L_sum * L_sum / 8.0;
                            const Tg value =
                                -((2 - (L_sum * L_sum * s * s)) * Kokkos::cos(L_sum * s)
                                  + 2 * L_sum * s * Kokkos::sin(L_sum * s) - 2)
                                / (2 * s * s * s * s + isOrig * 1.0);
 
                            apply(view_g, args) = (!isOrig) * value + isOrig * analyticLim;
                        });
                }
 
                // start a timer
                static IpplTimings::TimerRef fft4 = IpplTimings::getTimer("FFT: Precomputation");
                IpplTimings::startTimer(fft4);
 
                // inverse Fourier transform of the green's function for precomputation
                fft4n_m->transform(BACKWARD, grnL_m);
 
                IpplTimings::stopTimer(fft4);
 
                // Restrict transformed grnL_m to 2N domain after precomputation step
 
                // get the field data first
                typename Field_t::view_type view = grn_mr.getView();
                const int nghost                 = grn_mr.getNghost();
                const auto& ldom                 = layout2_m->getLocalNDIndex();
 
                // start a timer
                static IpplTimings::TimerRef ifftshift = IpplTimings::getTimer("Vico shift loop");
                IpplTimings::startTimer(ifftshift);
 
                // get number of ranks to see if need communication
                const int ranks = Comm->size();
 
                if (ranks > 1) {
                    communicateVico(size, view_g, ldom_g, nghost_g, view, ldom, nghost);
                } else {
                    // restrict the green's function to a (2N)^3 grid from the (4N)^3 grid
                    using mdrange_type = Kokkos::MDRangePolicy<Kokkos::Rank<3>>;
                    Kokkos::parallel_for(
                        "Restrict domain of Green's function from 4N to 2N",
                        mdrange_type({nghost, nghost, nghost}, {view.extent(0) - nghost - size[0],
                                                                view.extent(1) - nghost - size[1],
                                                                view.extent(2) - nghost - size[2]}),
                        KOKKOS_LAMBDA(const int i, const int j, const int k) {
                            // go from local indices to global
                            const int ig = i + ldom[0].first() - nghost;
                            const int jg = j + ldom[1].first() - nghost;
                            const int kg = k + ldom[2].first() - nghost;
 
                            const int ig2 = i + ldom_g[0].first() - nghost_g;
                            const int jg2 = j + ldom_g[1].first() - nghost_g;
                            const int kg2 = k + ldom_g[2].first() - nghost_g;
 
                            const bool isOrig = (ig == ig2) && (jg == jg2) && (kg == kg2);
                            view(i, j, k)     = isOrig * real(view_g(i, j, k));
 
                            // Now fill the rest of the field
                            const int s = 2 * size[0] - ig - 1 - ldom_g[0].first() + nghost_g;
                            const int p = 2 * size[1] - jg - 1 - ldom_g[1].first() + nghost_g;
                            const int q = 2 * size[2] - kg - 1 - ldom_g[2].first() + nghost_g;
 
                            view(s, j, k) = real(view_g(i + 1, j, k));
                            view(i, p, k) = real(view_g(i, j + 1, k));
                            view(i, j, q) = real(view_g(i, j, k + 1));
                            view(s, j, q) = real(view_g(i + 1, j, k + 1));
                            view(s, p, k) = real(view_g(i + 1, j + 1, k));
                            view(i, p, q) = real(view_g(i, j + 1, k + 1));
                            view(s, p, q) = real(view_g(i + 1, j + 1, k + 1));
                        });
                }
                IpplTimings::stopTimer(ifftshift);
            }
        } else if (alg == Algorithm::DCT_VICO) {
            if constexpr (Dim == 3) {
                Vector_t l(hr_m * nr_m);
                Vector_t hs_m;
                double L_sum(0.0);
 
                // compute length of the physical domain
                // compute Fourier domain spacing
                for (unsigned int i = 0; i < Dim; ++i) {
                    hs_m[i] = Kokkos::numbers::pi_v<Trhs> * 0.5 / l[i];
                    L_sum   = L_sum + l[i] * l[i];
                }
 
                // define the origin of the 2N+1 grid
                Vector_t origin = 0.0;
 
                // set mesh for the 2N+1 mesh
                mesh2n1_m->setMeshSpacing(hs_m);
 
                // size of truncation window
                L_sum = Kokkos::sqrt(L_sum);
                L_sum = 1.1 * L_sum;
 
                // initialize grn2n1_m
                Vector<int, Dim> size = nr_m;
 
                // initialize grn2n1_m
                typename Field_t::view_type view_g2n1 = grn2n1_m.getView();
                const int nghost_g2n1                 = grn2n1_m.getNghost();
                const auto& ldom_g2n1                 = layout2n1_m->getLocalNDIndex();
 
                ippl::parallel_for(
                    "Initialize 2N+1 Green's function ", grn2n1_m.getFieldRangePolicy(),
                    KOKKOS_LAMBDA(const index_array_type& args) {
                        scalar_type checkVal = 0;
 
                        // go from local indices to global
                        Vector<int, Dim> igVec = args - nghost_g2n1;
                        for (unsigned d = 0; d < Dim; ++d) {
                            igVec[d] += ldom_g2n1[d].first();
                            checkVal += igVec[d];
                        }
 
                        // compute s
                        Tg s = 0;
                        for (unsigned d = 0; d < Dim; ++d) {
                            double t = igVec[d] * hs_m[d];
                            s += (t * t);
                        }
                        s = Kokkos::sqrt(s);
 
                        const bool isOrig = (checkVal == 0);
                        const double val  = -2.0
                                           * (Kokkos::sin(0.5 * L_sum * s) / (s + isOrig * 1.0))
                                           * (Kokkos::sin(0.5 * L_sum * s) / (s + isOrig * 1.0));
                        const double analyticLim = -L_sum * L_sum * 0.5;
 
                        apply(view_g2n1, args) = ((!isOrig) * val) + (isOrig * analyticLim);
                    });
 
                // start a timer
                static IpplTimings::TimerRef fft4 = IpplTimings::getTimer("FFT: Precomputation");
                IpplTimings::startTimer(fft4);
 
                // inverse DCT transform of 2N+1 green's function for the precomputation
                fft2n1_m->transform(BACKWARD, grn2n1_m);
 
                IpplTimings::stopTimer(fft4);
 
                // Restrict transformed grn2n1_m to 2N domain after precomputation step
 
                // get the field data first
                typename Field_t::view_type view = grn_mr.getView();
                const int nghost                 = grn_mr.getNghost();
                const auto& ldom                 = layout2_m->getLocalNDIndex();
 
                // start a timer
                static IpplTimings::TimerRef ifftshift = IpplTimings::getTimer("Vico shift loop");
                IpplTimings::startTimer(ifftshift);
 
                // get number of ranks to see if need communication
                int ranks = ippl::Comm->size();
 
                // range type for Kokkos loop
                using mdrange_type = Kokkos::MDRangePolicy<Kokkos::Rank<3>>;
 
                if (ranks > 1) {
                    communicateVico(size, view_g2n1, ldom_g2n1, nghost_g2n1, view, ldom, nghost);
                } else {
                    // restrict the green's function to a (2N)^3 grid from the (2N+1)^3 grid
                    Kokkos::parallel_for(
                        "Restrict domain of Green's function from 2N+1 to 2N",
                        mdrange_type({nghost, nghost, nghost}, {view.extent(0) - nghost - size[0],
                                                                view.extent(1) - nghost - size[1],
                                                                view.extent(2) - nghost - size[2]}),
                        KOKKOS_LAMBDA(const int i, const int j, const int k) {
                            // go from local indices to global
                            const int ig = i + ldom[0].first() - nghost;
                            const int jg = j + ldom[1].first() - nghost;
                            const int kg = k + ldom[2].first() - nghost;
 
                            const int ig2 = i + ldom_g2n1[0].first() - nghost_g2n1;
                            const int jg2 = j + ldom_g2n1[1].first() - nghost_g2n1;
                            const int kg2 = k + ldom_g2n1[2].first() - nghost_g2n1;
 
                            const bool isOrig = ((ig == ig2) && (jg == jg2) && (kg == kg2));
                            view(i, j, k)     = isOrig * view_g2n1(i, j, k);
 
                            // Now fill the rest of the field
                            const int s = 2 * size[0] - ig - 1 - ldom_g2n1[0].first() + nghost_g2n1;
                            const int p = 2 * size[1] - jg - 1 - ldom_g2n1[1].first() + nghost_g2n1;
                            const int q = 2 * size[2] - kg - 1 - ldom_g2n1[2].first() + nghost_g2n1;
 
                            view(s, j, k) = view_g2n1(i + 1, j, k);
                            view(i, p, k) = view_g2n1(i, j + 1, k);
                            view(i, j, q) = view_g2n1(i, j, k + 1);
                            view(s, j, q) = view_g2n1(i + 1, j, k + 1);
                            view(s, p, k) = view_g2n1(i + 1, j + 1, k);
                            view(i, p, q) = view_g2n1(i, j + 1, k + 1);
                            view(s, p, q) = view_g2n1(i + 1, j + 1, k + 1);
                        });
                }
            }
        } else {
            // Hockney case
 
            typename Field_t::view_type view = grn_mr.getView();
            const int nghost                 = grn_mr.getNghost();
            const auto& ldom                 = layout2_m->getLocalNDIndex();
 
            if (greensFunctionType == GreenFunction::INTEGRATED) {
                Vector<int, Dim> nr = nr_m;
                vector_type hs      = hr_m;
 
                ippl::parallel_for(
                    "Integrated Green's function ", grn_mr.getFieldRangePolicy(),
                    KOKKOS_LAMBDA(const index_array_type& args) {
                        Vector<int, Dim> igVec = args - nghost;
                        for (unsigned d = 0; d < Dim; ++d) {
                            igVec[d] += ldom[d].first();
                        }
 
                        scalar_type offset[Dim];
                        for (unsigned int d = 0; d < Dim; ++d) {
                            const int igSigned =
                                (igVec[d] < nr[d]) ? igVec[d] : igVec[d] - 2 * nr[d];
                            offset[d] = static_cast<scalar_type>(igSigned) * hs[d];
                        }
 
                        const scalar_type avg =
                            integratedGreenAverage(offset[0], offset[1], offset[2], hs);
                        apply(view, args) = -avg / (scalar_type(4) * pi);
                    });
            } else {
                // calculate square of the mesh spacing for each dimension
                Vector_t hrsq(hr_m * hr_m);
 
                // use the grnIField_m helper field to compute Green's function
                for (unsigned int i = 0; i < Dim; ++i) {
                    grn_mr = grn_mr + grnIField_m[i] * hrsq[i];
                }
 
                // Formula of Green's function (2D and 3D supported for Hockney)
                if (Dim == 2) {
                    grn_mr = log(sqrt(grn_mr)) / (2 * pi);
                } else if (Dim == 3) {
                    grn_mr = -1.0 / (4.0 * pi * sqrt(grn_mr));
                }
 
                // Kokkos parallel for loop to find (0,0,0) point and regularize
                ippl::parallel_for(
                    "Regularize Green's function ", grn_mr.getFieldRangePolicy(),
                    KOKKOS_LAMBDA(const index_array_type& args) {
                        scalar_type checkVal = 0;
 
                        // go from local indices to global
                        Vector<int, Dim> igVec = args - nghost;
                        for (unsigned d = 0; d < Dim; ++d) {
                            igVec[d] += ldom[d].first();
                            checkVal += igVec[d];
                        }
 
                        // if (0,0,0), assign to it 1/(4*pi)
                        const bool isOrig = (checkVal == 0);
                        apply(view, args) =
                            isOrig * (-1.0 / (4.0 * pi)) + (!isOrig) * apply(view, args);
                    });
            }
        }
 
        // start a timer
        static IpplTimings::TimerRef fftg = IpplTimings::getTimer("FFT: Green");
        IpplTimings::startTimer(fftg);
 
        // perform the FFT of the Green's function for the convolution
        fft_m->transform(FORWARD, grn_mr, grntr_m);
 
        IpplTimings::stopTimer(fftg);
    };
 
    template <typename FieldLHS, typename FieldRHS>
    void FFTOpenPoissonSolver<FieldLHS, FieldRHS>::communicateVico(
        Vector<int, Dim> size, typename CxField_gt::view_type view_g,
        const ippl::NDIndex<Dim> ldom_g, const int nghost_g, typename Field_t::view_type view,
        const ippl::NDIndex<Dim> ldom, const int nghost) {
        const auto& lDomains2 = layout2_m->getHostLocalDomains();
        const auto& lDomains4 = layout4_m->getHostLocalDomains();
 
        std::vector<MPI_Request> requests(0);
        const int ranks = Comm->size();
 
        // 1st step: Define 8 domains corresponding to the different quadrants
        ippl::NDIndex<Dim> none;
        for (unsigned i = 0; i < Dim; i++) {
            none[i] = ippl::Index(size[i]);
        }
 
        ippl::NDIndex<Dim> x;
        x[0] = ippl::Index(size[0], 2 * size[0] - 1);
        x[1] = ippl::Index(size[1]);
        x[2] = ippl::Index(size[2]);
 
        ippl::NDIndex<Dim> y;
        y[0] = ippl::Index(size[0]);
        y[1] = ippl::Index(size[1], 2 * size[1] - 1);
        y[2] = ippl::Index(size[2]);
 
        ippl::NDIndex<Dim> z;
        z[0] = ippl::Index(size[0]);
        z[1] = ippl::Index(size[1]);
        z[2] = ippl::Index(size[2], 2 * size[2] - 1);
 
        ippl::NDIndex<Dim> xy;
        xy[0] = ippl::Index(size[0], 2 * size[0] - 1);
        xy[1] = ippl::Index(size[1], 2 * size[1] - 1);
        xy[2] = ippl::Index(size[2]);
 
        ippl::NDIndex<Dim> xz;
        xz[0] = ippl::Index(size[0], 2 * size[0] - 1);
        xz[1] = ippl::Index(size[1]);
        xz[2] = ippl::Index(size[2], 2 * size[2] - 1);
 
        ippl::NDIndex<Dim> yz;
        yz[0] = ippl::Index(size[0]);
        yz[1] = ippl::Index(size[1], 2 * size[1] - 1);
        yz[2] = ippl::Index(size[2], 2 * size[2] - 1);
 
        ippl::NDIndex<Dim> xyz;
        for (unsigned i = 0; i < Dim; i++) {
            xyz[i] = ippl::Index(size[i], 2 * size[i] - 1);
        }
 
        // 2nd step: send
        for (int i = 0; i < ranks; ++i) {
            auto domain2 = lDomains2[i];
 
            if (domain2.touches(none)) {
                auto intersection = domain2.intersect(none);
 
                if (ldom_g.touches(intersection)) {
                    intersection = intersection.intersect(ldom_g);
 
                    solver_send(mpi::tag::VICO_SOLVER, 0, i, intersection, ldom_g, nghost_g, view_g,
                                fd_m, requests);
                }
            }
 
            if (domain2.touches(x)) {
                auto intersection = domain2.intersect(x);
                auto xdom         = ippl::Index((2 * size[0] - intersection[0].first()),
                                                (2 * size[0] - intersection[0].last()), -1);
 
                ippl::NDIndex<Dim> domain4;
                domain4[0] = xdom;
                domain4[1] = intersection[1];
                domain4[2] = intersection[2];
 
                if (ldom_g.touches(domain4)) {
                    intersection = ldom_g.intersect(domain4);
 
                    solver_send(mpi::tag::VICO_SOLVER, 1, i, intersection, ldom_g, nghost_g, view_g,
                                fd_m, requests);
                }
            }
 
            if (domain2.touches(y)) {
                auto intersection = domain2.intersect(y);
                auto ydom         = ippl::Index((2 * size[1] - intersection[1].first()),
                                                (2 * size[1] - intersection[1].last()), -1);
 
                ippl::NDIndex<Dim> domain4;
                domain4[0] = intersection[0];
                domain4[1] = ydom;
                domain4[2] = intersection[2];
 
                if (ldom_g.touches(domain4)) {
                    intersection = ldom_g.intersect(domain4);
 
                    solver_send(mpi::tag::VICO_SOLVER, 2, i, intersection, ldom_g, nghost_g, view_g,
                                fd_m, requests);
                }
            }
 
            if (domain2.touches(z)) {
                auto intersection = domain2.intersect(z);
                auto zdom         = ippl::Index((2 * size[2] - intersection[2].first()),
                                                (2 * size[2] - intersection[2].last()), -1);
 
                ippl::NDIndex<Dim> domain4;
                domain4[0] = intersection[0];
                domain4[1] = intersection[1];
                domain4[2] = zdom;
 
                if (ldom_g.touches(domain4)) {
                    intersection = ldom_g.intersect(domain4);
 
                    solver_send(mpi::tag::VICO_SOLVER, 3, i, intersection, ldom_g, nghost_g, view_g,
                                fd_m, requests);
                }
            }
 
            if (domain2.touches(xy)) {
                auto intersection = domain2.intersect(xy);
                auto xdom         = ippl::Index((2 * size[0] - intersection[0].first()),
                                                (2 * size[0] - intersection[0].last()), -1);
                auto ydom         = ippl::Index((2 * size[1] - intersection[1].first()),
                                                (2 * size[1] - intersection[1].last()), -1);
 
                ippl::NDIndex<Dim> domain4;
                domain4[0] = xdom;
                domain4[1] = ydom;
                domain4[2] = intersection[2];
 
                if (ldom_g.touches(domain4)) {
                    intersection = ldom_g.intersect(domain4);
 
                    solver_send(mpi::tag::VICO_SOLVER, 4, i, intersection, ldom_g, nghost_g, view_g,
                                fd_m, requests);
                }
            }
 
            if (domain2.touches(yz)) {
                auto intersection = domain2.intersect(yz);
                auto ydom         = ippl::Index((2 * size[1] - intersection[1].first()),
                                                (2 * size[1] - intersection[1].last()), -1);
                auto zdom         = ippl::Index((2 * size[2] - intersection[2].first()),
                                                (2 * size[2] - intersection[2].last()), -1);
 
                ippl::NDIndex<Dim> domain4;
                domain4[0] = intersection[0];
                domain4[1] = ydom;
                domain4[2] = zdom;
 
                if (ldom_g.touches(domain4)) {
                    intersection = ldom_g.intersect(domain4);
                    solver_send(mpi::tag::VICO_SOLVER, 5, i, intersection, ldom_g, nghost_g, view_g,
                                fd_m, requests);
                }
            }
 
            if (domain2.touches(xz)) {
                auto intersection = domain2.intersect(xz);
                auto xdom         = ippl::Index((2 * size[0] - intersection[0].first()),
                                                (2 * size[0] - intersection[0].last()), -1);
                auto zdom         = ippl::Index((2 * size[2] - intersection[2].first()),
                                                (2 * size[2] - intersection[2].last()), -1);
 
                ippl::NDIndex<Dim> domain4;
                domain4[0] = xdom;
                domain4[1] = intersection[1];
                domain4[2] = zdom;
 
                if (ldom_g.touches(domain4)) {
                    intersection = ldom_g.intersect(domain4);
                    solver_send(mpi::tag::VICO_SOLVER, 6, i, intersection, ldom_g, nghost_g, view_g,
                                fd_m, requests);
                }
            }
 
            if (domain2.touches(xyz)) {
                auto intersection = domain2.intersect(xyz);
                auto xdom         = ippl::Index((2 * size[0] - intersection[0].first()),
                                                (2 * size[0] - intersection[0].last()), -1);
                auto ydom         = ippl::Index((2 * size[1] - intersection[1].first()),
                                                (2 * size[1] - intersection[1].last()), -1);
                auto zdom         = ippl::Index((2 * size[2] - intersection[2].first()),
                                                (2 * size[2] - intersection[2].last()), -1);
 
                ippl::NDIndex<Dim> domain4;
                domain4[0] = xdom;
                domain4[1] = ydom;
                domain4[2] = zdom;
 
                if (ldom_g.touches(domain4)) {
                    intersection = ldom_g.intersect(domain4);
                    solver_send(mpi::tag::VICO_SOLVER, 7, i, intersection, ldom_g, nghost_g, view_g,
                                fd_m, requests);
                }
            }
        }
 
        // 3rd step: receive
        for (int i = 0; i < ranks; ++i) {
            if (ldom.touches(none)) {
                auto intersection = ldom.intersect(none);
 
                if (lDomains4[i].touches(intersection)) {
                    intersection = intersection.intersect(lDomains4[i]);
 
                    solver_recv(mpi::tag::VICO_SOLVER, 0, i, intersection, ldom, nghost, view,
                                fd_m);
                }
            }
 
            if (ldom.touches(x)) {
                auto intersection = ldom.intersect(x);
 
                auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
                                        (2 * size[0] - intersection[0].last()), -1);
 
                ippl::NDIndex<Dim> domain4;
                domain4[0] = xdom;
                domain4[1] = intersection[1];
                domain4[2] = intersection[2];
 
                if (lDomains4[i].touches(domain4)) {
                    domain4    = lDomains4[i].intersect(domain4);
                    domain4[0] = ippl::Index(2 * size[0] - domain4[0].first(),
                                             2 * size[0] - domain4[0].last(), -1);
 
                    intersection = intersection.intersect(domain4);
 
                    solver_recv(mpi::tag::VICO_SOLVER, 1, i, intersection, ldom, nghost, view, fd_m,
                                {true, false, false});
                }
            }
 
            if (ldom.touches(y)) {
                auto intersection = ldom.intersect(y);
 
                auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
                                        (2 * size[1] - intersection[1].last()), -1);
 
                ippl::NDIndex<Dim> domain4;
                domain4[0] = intersection[0];
                domain4[1] = ydom;
                domain4[2] = intersection[2];
 
                if (lDomains4[i].touches(domain4)) {
                    domain4    = lDomains4[i].intersect(domain4);
                    domain4[1] = ippl::Index(2 * size[1] - domain4[1].first(),
                                             2 * size[1] - domain4[1].last(), -1);
 
                    intersection = intersection.intersect(domain4);
 
                    solver_recv(mpi::tag::VICO_SOLVER, 2, i, intersection, ldom, nghost, view, fd_m,
                                {false, true, false});
                }
            }
 
            if (ldom.touches(z)) {
                auto intersection = ldom.intersect(z);
 
                auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
                                        (2 * size[2] - intersection[2].last()), -1);
 
                ippl::NDIndex<Dim> domain4;
                domain4[0] = intersection[0];
                domain4[1] = intersection[1];
                domain4[2] = zdom;
 
                if (lDomains4[i].touches(domain4)) {
                    domain4    = lDomains4[i].intersect(domain4);
                    domain4[2] = ippl::Index(2 * size[2] - domain4[2].first(),
                                             2 * size[2] - domain4[2].last(), -1);
 
                    intersection = intersection.intersect(domain4);
 
                    solver_recv(mpi::tag::VICO_SOLVER, 3, i, intersection, ldom, nghost, view, fd_m,
                                {false, false, true});
                }
            }
 
            if (ldom.touches(xy)) {
                auto intersection = ldom.intersect(xy);
 
                auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
                                        (2 * size[0] - intersection[0].last()), -1);
                auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
                                        (2 * size[1] - intersection[1].last()), -1);
 
                ippl::NDIndex<Dim> domain4;
                domain4[0] = xdom;
                domain4[1] = ydom;
                domain4[2] = intersection[2];
 
                if (lDomains4[i].touches(domain4)) {
                    domain4    = lDomains4[i].intersect(domain4);
                    domain4[0] = ippl::Index(2 * size[0] - domain4[0].first(),
                                             2 * size[0] - domain4[0].last(), -1);
                    domain4[1] = ippl::Index(2 * size[1] - domain4[1].first(),
                                             2 * size[1] - domain4[1].last(), -1);
 
                    intersection = intersection.intersect(domain4);
 
                    solver_recv(mpi::tag::VICO_SOLVER, 4, i, intersection, ldom, nghost, view, fd_m,
                                {true, true, false});
                }
            }
 
            if (ldom.touches(yz)) {
                auto intersection = ldom.intersect(yz);
 
                auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
                                        (2 * size[1] - intersection[1].last()), -1);
                auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
                                        (2 * size[2] - intersection[2].last()), -1);
 
                ippl::NDIndex<Dim> domain4;
                domain4[0] = intersection[0];
                domain4[1] = ydom;
                domain4[2] = zdom;
 
                if (lDomains4[i].touches(domain4)) {
                    domain4    = lDomains4[i].intersect(domain4);
                    domain4[1] = ippl::Index(2 * size[1] - domain4[1].first(),
                                             2 * size[1] - domain4[1].last(), -1);
                    domain4[2] = ippl::Index(2 * size[2] - domain4[2].first(),
                                             2 * size[2] - domain4[2].last(), -1);
 
                    intersection = intersection.intersect(domain4);
 
                    solver_recv(mpi::tag::VICO_SOLVER, 5, i, intersection, ldom, nghost, view, fd_m,
                                {false, true, true});
                }
            }
 
            if (ldom.touches(xz)) {
                auto intersection = ldom.intersect(xz);
 
                auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
                                        (2 * size[0] - intersection[0].last()), -1);
                auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
                                        (2 * size[2] - intersection[2].last()), -1);
 
                ippl::NDIndex<Dim> domain4;
                domain4[0] = xdom;
                domain4[1] = intersection[1];
                domain4[2] = zdom;
 
                if (lDomains4[i].touches(domain4)) {
                    domain4    = lDomains4[i].intersect(domain4);
                    domain4[0] = ippl::Index(2 * size[0] - domain4[0].first(),
                                             2 * size[0] - domain4[0].last(), -1);
                    domain4[2] = ippl::Index(2 * size[2] - domain4[2].first(),
                                             2 * size[2] - domain4[2].last(), -1);
 
                    intersection = intersection.intersect(domain4);
 
                    solver_recv(mpi::tag::VICO_SOLVER, 6, i, intersection, ldom, nghost, view, fd_m,
                                {true, false, true});
                }
            }
 
            if (ldom.touches(xyz)) {
                auto intersection = ldom.intersect(xyz);
 
                auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
                                        (2 * size[0] - intersection[0].last()), -1);
                auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
                                        (2 * size[1] - intersection[1].last()), -1);
                auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
                                        (2 * size[2] - intersection[2].last()), -1);
 
                ippl::NDIndex<Dim> domain4;
                domain4[0] = xdom;
                domain4[1] = ydom;
                domain4[2] = zdom;
 
                if (lDomains4[i].touches(domain4)) {
                    domain4    = lDomains4[i].intersect(domain4);
                    domain4[0] = ippl::Index(2 * size[0] - domain4[0].first(),
                                             2 * size[0] - domain4[0].last(), -1);
                    domain4[1] = ippl::Index(2 * size[1] - domain4[1].first(),
                                             2 * size[1] - domain4[1].last(), -1);
                    domain4[2] = ippl::Index(2 * size[2] - domain4[2].first(),
                                             2 * size[2] - domain4[2].last(), -1);
 
                    intersection = intersection.intersect(domain4);
 
                    solver_recv(mpi::tag::VICO_SOLVER, 7, i, intersection, ldom, nghost, view, fd_m,
                                {true, true, true});
                }
            }
        }
 
        if (requests.size() > 0) {
            MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
        }
        ippl::Comm->freeAllBuffers();
    };
 
    // CommunicateVico for DCT_VICO (2N+1 to 2N)
    template <typename FieldLHS, typename FieldRHS>
    void FFTOpenPoissonSolver<FieldLHS, FieldRHS>::communicateVico(
        Vector<int, Dim> size, typename Field_t::view_type view_g, const ippl::NDIndex<Dim> ldom_g,
        const int nghost_g, typename Field_t::view_type view, const ippl::NDIndex<Dim> ldom,
        const int nghost) {
        const auto& lDomains2   = layout2_m->getHostLocalDomains();
        const auto& lDomains2n1 = layout2n1_m->getHostLocalDomains();
 
        std::vector<MPI_Request> requests(0);
        const int ranks = ippl::Comm->size();
 
        // 1st step: Define 8 domains corresponding to the different quadrants
        ippl::NDIndex<Dim> none;
        for (unsigned i = 0; i < Dim; i++) {
            none[i] = ippl::Index(size[i]);
        }
 
        ippl::NDIndex<Dim> x;
        x[0] = ippl::Index(size[0], 2 * size[0] - 1);
        x[1] = ippl::Index(size[1]);
        x[2] = ippl::Index(size[2]);
 
        ippl::NDIndex<Dim> y;
        y[0] = ippl::Index(size[0]);
        y[1] = ippl::Index(size[1], 2 * size[1] - 1);
        y[2] = ippl::Index(size[2]);
 
        ippl::NDIndex<Dim> z;
        z[0] = ippl::Index(size[0]);
        z[1] = ippl::Index(size[1]);
        z[2] = ippl::Index(size[2], 2 * size[2] - 1);
 
        ippl::NDIndex<Dim> xy;
        xy[0] = ippl::Index(size[0], 2 * size[0] - 1);
        xy[1] = ippl::Index(size[1], 2 * size[1] - 1);
        xy[2] = ippl::Index(size[2]);
 
        ippl::NDIndex<Dim> xz;
        xz[0] = ippl::Index(size[0], 2 * size[0] - 1);
        xz[1] = ippl::Index(size[1]);
        xz[2] = ippl::Index(size[2], 2 * size[2] - 1);
 
        ippl::NDIndex<Dim> yz;
        yz[0] = ippl::Index(size[0]);
        yz[1] = ippl::Index(size[1], 2 * size[1] - 1);
        yz[2] = ippl::Index(size[2], 2 * size[2] - 1);
 
        ippl::NDIndex<Dim> xyz;
        for (unsigned i = 0; i < Dim; i++) {
            xyz[i] = ippl::Index(size[i], 2 * size[i] - 1);
        }
 
        // 2nd step: send
        for (int i = 0; i < ranks; ++i) {
            auto domain2 = lDomains2[i];
 
            if (domain2.touches(none)) {
                auto intersection = domain2.intersect(none);
 
                if (ldom_g.touches(intersection)) {
                    intersection = intersection.intersect(ldom_g);
 
                    solver_send(mpi::tag::VICO_SOLVER, 0, i, intersection, ldom_g, nghost_g, view_g,
                                fd_m, requests);
                }
            }
 
            if (domain2.touches(x)) {
                auto intersection = domain2.intersect(x);
                auto xdom         = ippl::Index((2 * size[0] - intersection[0].first()),
                                                (2 * size[0] - intersection[0].last()), -1);
 
                ippl::NDIndex<Dim> domain2n1;
                domain2n1[0] = xdom;
                domain2n1[1] = intersection[1];
                domain2n1[2] = intersection[2];
 
                if (ldom_g.touches(domain2n1)) {
                    intersection = ldom_g.intersect(domain2n1);
 
                    solver_send(mpi::tag::VICO_SOLVER, 1, i, intersection, ldom_g, nghost_g, view_g,
                                fd_m, requests);
                }
            }
 
            if (domain2.touches(y)) {
                auto intersection = domain2.intersect(y);
                auto ydom         = ippl::Index((2 * size[1] - intersection[1].first()),
                                                (2 * size[1] - intersection[1].last()), -1);
 
                ippl::NDIndex<Dim> domain2n1;
                domain2n1[0] = intersection[0];
                domain2n1[1] = ydom;
                domain2n1[2] = intersection[2];
 
                if (ldom_g.touches(domain2n1)) {
                    intersection = ldom_g.intersect(domain2n1);
 
                    solver_send(mpi::tag::VICO_SOLVER, 2, i, intersection, ldom_g, nghost_g, view_g,
                                fd_m, requests);
                }
            }
 
            if (domain2.touches(z)) {
                auto intersection = domain2.intersect(z);
                auto zdom         = ippl::Index((2 * size[2] - intersection[2].first()),
                                                (2 * size[2] - intersection[2].last()), -1);
 
                ippl::NDIndex<Dim> domain2n1;
                domain2n1[0] = intersection[0];
                domain2n1[1] = intersection[1];
                domain2n1[2] = zdom;
 
                if (ldom_g.touches(domain2n1)) {
                    intersection = ldom_g.intersect(domain2n1);
 
                    solver_send(mpi::tag::VICO_SOLVER, 3, i, intersection, ldom_g, nghost_g, view_g,
                                fd_m, requests);
                }
            }
 
            if (domain2.touches(xy)) {
                auto intersection = domain2.intersect(xy);
                auto xdom         = ippl::Index((2 * size[0] - intersection[0].first()),
                                                (2 * size[0] - intersection[0].last()), -1);
                auto ydom         = ippl::Index((2 * size[1] - intersection[1].first()),
                                                (2 * size[1] - intersection[1].last()), -1);
 
                ippl::NDIndex<Dim> domain2n1;
                domain2n1[0] = xdom;
                domain2n1[1] = ydom;
                domain2n1[2] = intersection[2];
 
                if (ldom_g.touches(domain2n1)) {
                    intersection = ldom_g.intersect(domain2n1);
 
                    solver_send(mpi::tag::VICO_SOLVER, 4, i, intersection, ldom_g, nghost_g, view_g,
                                fd_m, requests);
                }
            }
 
            if (domain2.touches(yz)) {
                auto intersection = domain2.intersect(yz);
                auto ydom         = ippl::Index((2 * size[1] - intersection[1].first()),
                                                (2 * size[1] - intersection[1].last()), -1);
                auto zdom         = ippl::Index((2 * size[2] - intersection[2].first()),
                                                (2 * size[2] - intersection[2].last()), -1);
 
                ippl::NDIndex<Dim> domain2n1;
                domain2n1[0] = intersection[0];
                domain2n1[1] = ydom;
                domain2n1[2] = zdom;
 
                if (ldom_g.touches(domain2n1)) {
                    intersection = ldom_g.intersect(domain2n1);
 
                    solver_send(mpi::tag::VICO_SOLVER, 5, i, intersection, ldom_g, nghost_g, view_g,
                                fd_m, requests);
                }
            }
 
            if (domain2.touches(xz)) {
                auto intersection = domain2.intersect(xz);
                auto xdom         = ippl::Index((2 * size[0] - intersection[0].first()),
                                                (2 * size[0] - intersection[0].last()), -1);
                auto zdom         = ippl::Index((2 * size[2] - intersection[2].first()),
                                                (2 * size[2] - intersection[2].last()), -1);
 
                ippl::NDIndex<Dim> domain2n1;
                domain2n1[0] = xdom;
                domain2n1[1] = intersection[1];
                domain2n1[2] = zdom;
 
                if (ldom_g.touches(domain2n1)) {
                    intersection = ldom_g.intersect(domain2n1);
 
                    solver_send(mpi::tag::VICO_SOLVER, 6, i, intersection, ldom_g, nghost_g, view_g,
                                fd_m, requests);
                }
            }
 
            if (domain2.touches(xyz)) {
                auto intersection = domain2.intersect(xyz);
                auto xdom         = ippl::Index((2 * size[0] - intersection[0].first()),
                                                (2 * size[0] - intersection[0].last()), -1);
                auto ydom         = ippl::Index((2 * size[1] - intersection[1].first()),
                                                (2 * size[1] - intersection[1].last()), -1);
                auto zdom         = ippl::Index((2 * size[2] - intersection[2].first()),
                                                (2 * size[2] - intersection[2].last()), -1);
 
                ippl::NDIndex<Dim> domain2n1;
                domain2n1[0] = xdom;
                domain2n1[1] = ydom;
                domain2n1[2] = zdom;
 
                if (ldom_g.touches(domain2n1)) {
                    intersection = ldom_g.intersect(domain2n1);
 
                    solver_send(mpi::tag::VICO_SOLVER, 7, i, intersection, ldom_g, nghost_g, view_g,
                                fd_m, requests);
                }
            }
        }
 
        // 3rd step: receive
        for (int i = 0; i < ranks; ++i) {
            if (ldom.touches(none)) {
                auto intersection = ldom.intersect(none);
 
                if (lDomains2n1[i].touches(intersection)) {
                    intersection = intersection.intersect(lDomains2n1[i]);
 
                    solver_recv(mpi::tag::VICO_SOLVER, 0, i, intersection, ldom, nghost, view,
                                fd_m);
                }
            }
 
            if (ldom.touches(x)) {
                auto intersection = ldom.intersect(x);
 
                auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
                                        (2 * size[0] - intersection[0].last()), -1);
 
                ippl::NDIndex<Dim> domain2n1;
                domain2n1[0] = xdom;
                domain2n1[1] = intersection[1];
                domain2n1[2] = intersection[2];
 
                if (lDomains2n1[i].touches(domain2n1)) {
                    domain2n1    = lDomains2n1[i].intersect(domain2n1);
                    domain2n1[0] = ippl::Index(2 * size[0] - domain2n1[0].first(),
                                               2 * size[0] - domain2n1[0].last(), -1);
 
                    intersection = intersection.intersect(domain2n1);
 
                    solver_recv(mpi::tag::VICO_SOLVER, 1, i, intersection, ldom, nghost, view, fd_m,
                                {true, false, false});
                }
            }
 
            if (ldom.touches(y)) {
                auto intersection = ldom.intersect(y);
 
                auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
                                        (2 * size[1] - intersection[1].last()), -1);
 
                ippl::NDIndex<Dim> domain2n1;
                domain2n1[0] = intersection[0];
                domain2n1[1] = ydom;
                domain2n1[2] = intersection[2];
 
                if (lDomains2n1[i].touches(domain2n1)) {
                    domain2n1    = lDomains2n1[i].intersect(domain2n1);
                    domain2n1[1] = ippl::Index(2 * size[1] - domain2n1[1].first(),
                                               2 * size[1] - domain2n1[1].last(), -1);
 
                    intersection = intersection.intersect(domain2n1);
 
                    solver_recv(mpi::tag::VICO_SOLVER, 2, i, intersection, ldom, nghost, view, fd_m,
                                {false, true, false});
                }
            }
 
            if (ldom.touches(z)) {
                auto intersection = ldom.intersect(z);
 
                auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
                                        (2 * size[2] - intersection[2].last()), -1);
 
                ippl::NDIndex<Dim> domain2n1;
                domain2n1[0] = intersection[0];
                domain2n1[1] = intersection[1];
                domain2n1[2] = zdom;
 
                if (lDomains2n1[i].touches(domain2n1)) {
                    domain2n1    = lDomains2n1[i].intersect(domain2n1);
                    domain2n1[2] = ippl::Index(2 * size[2] - domain2n1[2].first(),
                                               2 * size[2] - domain2n1[2].last(), -1);
 
                    intersection = intersection.intersect(domain2n1);
 
                    solver_recv(mpi::tag::VICO_SOLVER, 3, i, intersection, ldom, nghost, view, fd_m,
                                {false, false, true});
                }
            }
 
            if (ldom.touches(xy)) {
                auto intersection = ldom.intersect(xy);
 
                auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
                                        (2 * size[0] - intersection[0].last()), -1);
                auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
                                        (2 * size[1] - intersection[1].last()), -1);
 
                ippl::NDIndex<Dim> domain2n1;
                domain2n1[0] = xdom;
                domain2n1[1] = ydom;
                domain2n1[2] = intersection[2];
 
                if (lDomains2n1[i].touches(domain2n1)) {
                    domain2n1    = lDomains2n1[i].intersect(domain2n1);
                    domain2n1[0] = ippl::Index(2 * size[0] - domain2n1[0].first(),
                                               2 * size[0] - domain2n1[0].last(), -1);
                    domain2n1[1] = ippl::Index(2 * size[1] - domain2n1[1].first(),
                                               2 * size[1] - domain2n1[1].last(), -1);
 
                    intersection = intersection.intersect(domain2n1);
 
                    solver_recv(mpi::tag::VICO_SOLVER, 4, i, intersection, ldom, nghost, view, fd_m,
                                {true, true, false});
                }
            }
 
            if (ldom.touches(yz)) {
                auto intersection = ldom.intersect(yz);
 
                auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
                                        (2 * size[1] - intersection[1].last()), -1);
                auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
                                        (2 * size[2] - intersection[2].last()), -1);
 
                ippl::NDIndex<Dim> domain2n1;
                domain2n1[0] = intersection[0];
                domain2n1[1] = ydom;
                domain2n1[2] = zdom;
 
                if (lDomains2n1[i].touches(domain2n1)) {
                    domain2n1    = lDomains2n1[i].intersect(domain2n1);
                    domain2n1[1] = ippl::Index(2 * size[1] - domain2n1[1].first(),
                                               2 * size[1] - domain2n1[1].last(), -1);
                    domain2n1[2] = ippl::Index(2 * size[2] - domain2n1[2].first(),
                                               2 * size[2] - domain2n1[2].last(), -1);
 
                    intersection = intersection.intersect(domain2n1);
 
                    solver_recv(mpi::tag::VICO_SOLVER, 5, i, intersection, ldom, nghost, view, fd_m,
                                {false, true, true});
                }
            }
 
            if (ldom.touches(xz)) {
                auto intersection = ldom.intersect(xz);
 
                auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
                                        (2 * size[0] - intersection[0].last()), -1);
                auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
                                        (2 * size[2] - intersection[2].last()), -1);
 
                ippl::NDIndex<Dim> domain2n1;
                domain2n1[0] = xdom;
                domain2n1[1] = intersection[1];
                domain2n1[2] = zdom;
 
                if (lDomains2n1[i].touches(domain2n1)) {
                    domain2n1    = lDomains2n1[i].intersect(domain2n1);
                    domain2n1[0] = ippl::Index(2 * size[0] - domain2n1[0].first(),
                                               2 * size[0] - domain2n1[0].last(), -1);
                    domain2n1[2] = ippl::Index(2 * size[2] - domain2n1[2].first(),
                                               2 * size[2] - domain2n1[2].last(), -1);
 
                    intersection = intersection.intersect(domain2n1);
 
                    solver_recv(mpi::tag::VICO_SOLVER, 6, i, intersection, ldom, nghost, view, fd_m,
                                {true, false, true});
                }
            }
 
            if (ldom.touches(xyz)) {
                auto intersection = ldom.intersect(xyz);
 
                auto xdom = ippl::Index((2 * size[0] - intersection[0].first()),
                                        (2 * size[0] - intersection[0].last()), -1);
                auto ydom = ippl::Index((2 * size[1] - intersection[1].first()),
                                        (2 * size[1] - intersection[1].last()), -1);
                auto zdom = ippl::Index((2 * size[2] - intersection[2].first()),
                                        (2 * size[2] - intersection[2].last()), -1);
 
                ippl::NDIndex<Dim> domain2n1;
                domain2n1[0] = xdom;
                domain2n1[1] = ydom;
                domain2n1[2] = zdom;
 
                if (lDomains2n1[i].touches(domain2n1)) {
                    domain2n1    = lDomains2n1[i].intersect(domain2n1);
                    domain2n1[0] = ippl::Index(2 * size[0] - domain2n1[0].first(),
                                               2 * size[0] - domain2n1[0].last(), -1);
                    domain2n1[1] = ippl::Index(2 * size[1] - domain2n1[1].first(),
                                               2 * size[1] - domain2n1[1].last(), -1);
                    domain2n1[2] = ippl::Index(2 * size[2] - domain2n1[2].first(),
                                               2 * size[2] - domain2n1[2].last(), -1);
 
                    intersection = intersection.intersect(domain2n1);
 
                    solver_recv(mpi::tag::VICO_SOLVER, 7, i, intersection, ldom, nghost, view, fd_m,
                                {true, true, true});
                }
            }
        }
 
        if (requests.size() > 0) {
            MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
        }
        ippl::Comm->freeAllBuffers();
    };
 
    // Shifted Green's function: fills grn_mr with G(r - shift) on the doubled
    // grid and caches FFT into grntr_m. Kokkos-parallel, mirrors the HOCKNEY
    // greensFunction() method in structure.
    //
    // After this call, solve() will convolve rho with the shifted kernel.
    // To restore the standard kernel, call greensFunction() explicitly.
 
    template <typename FieldLHS, typename FieldRHS>
    void FFTOpenPoissonSolver<FieldLHS, FieldRHS>::shiftedGreensFunction(
        const Vector<double, Dim>& shift) {
        const int alg                = this->params_m.template get<int>("algorithm");
        const int greensFunctionType = this->params_m.template get<int>("greens_function");
        if (alg != Algorithm::HOCKNEY) {
            throw IpplException("FFTOpenPoissonSolver::shiftedGreensFunction",
                                "Shifted Green's function is only implemented for HOCKNEY.");
        }
        if ((greensFunctionType == GreenFunction::INTEGRATED) && (Dim != 3)) {
            throw IpplException("FFTOpenPoissonSolver::shiftedGreensFunction",
                                "Shifted integrated Green's function is only implemented for 3D.");
        }
 
        // Sync mesh spacing with the current RHS mesh (same logic as solve()'s
        // mesh-change detection). Without this, two failure modes compound:
        //   1. We would compute the shifted kernel at a STALE hr_m.
        //   2. A subsequent solve() would see hr_m != mesh->getMeshSpacing(),
        //      set green=true, and call greensFunction(), overwriting the
        //      shifted kernel with the standard one.
        // By updating hr_m (and the dependent mesh2_m / meshComplex_m) here,
        // solve()'s mesh check finds no change and leaves grntr_m intact.
        mesh_mp = &(this->rhs_mp->get_mesh());
        for (unsigned int i = 0; i < Dim; ++i) {
            hr_m[i] = mesh_mp->getMeshSpacing(i);
        }
        mesh2_m->setMeshSpacing(hr_m);
        meshComplex_m->setMeshSpacing(hr_m);
 
        const scalar_type pi = Kokkos::numbers::pi_v<scalar_type>;
 
        typename Field_t::view_type view = grn_mr.getView();
        const int nghost                 = grn_mr.getNghost();
        const auto& ldom                 = layout2_m->getLocalNDIndex();
 
        // Capture simple value arrays for the lambda.
        Vector<int, Dim> nr      = nr_m;
        Vector_t hs              = hr_m;
        Vector<double, Dim> shft = shift;
 
        using index_array_type = typename ippl::RangePolicy<Dim>::index_array_type;
        if (greensFunctionType == GreenFunction::INTEGRATED) {
            ippl::parallel_for(
                "Shifted integrated Green's function", grn_mr.getFieldRangePolicy(),
                KOKKOS_LAMBDA(const index_array_type& args) {
                    Vector<int, Dim> igVec = args - nghost;
                    for (unsigned d = 0; d < Dim; ++d) {
                        igVec[d] += ldom[d].first();
                    }
 
                    scalar_type offset[Dim];
                    for (unsigned int d = 0; d < Dim; ++d) {
                        const int igSigned = (igVec[d] < nr[d]) ? igVec[d] : igVec[d] - 2 * nr[d];
                        offset[d]          = static_cast<scalar_type>(igSigned) * hs[d]
                                    - static_cast<scalar_type>(shft[d]);
                    }
 
                    const scalar_type avg =
                        integratedGreenAverage(offset[0], offset[1], offset[2], hs);
                    apply(view, args) = -avg / (scalar_type(4) * pi);
                });
        } else {
            // Regularization threshold (axis-min mesh spacing, squared, quartered).
            scalar_type hmin2 = hs[0] * hs[0];
            for (unsigned int d = 1; d < Dim; ++d) {
                hmin2 = (hs[d] * hs[d] < hmin2) ? (hs[d] * hs[d]) : hmin2;
            }
            const scalar_type regThresh = 0.25 * hmin2;
 
            ippl::parallel_for(
                "Shifted Green's function", grn_mr.getFieldRangePolicy(),
                KOKKOS_LAMBDA(const index_array_type& args) {
                    // local -> global indices
                    Vector<int, Dim> igVec = args - nghost;
                    for (unsigned d = 0; d < Dim; ++d) {
                        igVec[d] += ldom[d].first();
                    }
 
                    // Hockney convention: map doubled-grid indices.
                    //   ig in [0, N)   -> offset =  ig        * h
                    //   ig in [N, 2N)  -> offset = (ig - 2N)  * h
                    // Subtract the shift to evaluate the shifted Green's function G(r-s).
                    double rsq = 0.0;
                    for (unsigned int d = 0; d < Dim; ++d) {
                        const double ig_signed = (igVec[d] < nr[d])
                                                     ? static_cast<double>(igVec[d])
                                                     : static_cast<double>(igVec[d] - 2 * nr[d]);
                        const double xoff      = ig_signed * hs[d] - shft[d];
                        rsq += xoff * xoff;
                    }
 
                    // Sign convention matches greensFunction() (HOCKNEY): grn_mr stores -G0
                    // so that solve()'s `rho2tr_m = -rho2tr_m * grntr_m` yields +conv(rho,G0)
                    // where G0(r) = 1/(4 pi |r|).
                    const bool nearSing = (rsq < regThresh);
                    const scalar_type r = Kokkos::sqrt(rsq + nearSing * regThresh);
                    apply(view, args)   = -1.0 / (4.0 * pi * r);
                });
        }
 
        // Store Fourier transform of shifted Green's function for convolution in solve()
        static IpplTimings::TimerRef fftsg = IpplTimings::getTimer("FFT: Shifted Green");
        IpplTimings::startTimer(fftsg);
 
        fft_m->transform(FORWARD, grn_mr, grntr_m);
 
        IpplTimings::stopTimer(fftsg);
    }
 
}  // namespace ippl