FAdaptChebKernel.hpp

#ifndef FADAPTCHEBKERNEL_HPP
#define FADAPTCHEBKERNEL_HPP

#include "Kernels/Chebyshev/FChebSymKernel.hpp"

#include "FComputeClassDescriptor.hpp"

#ifdef _OPENMP
#include <omp.h>
#endif

template<
    class FReal,
    class CellClass,
    class ContainerClass,
    class MatrixKernelClass,
    int ORDER,
    int NVALS = 1
    >
class FAdaptChebKernel : public FChebSymKernel<FReal, CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>{
    using FBase =  FChebSymKernel<FReal, CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>;

public:
    using FBase::FBase;

    using FBase::M2M;
    using FBase::M2L;
    using FBase::L2L;
    using FBase::P2P;

    void setup() {

        #ifdef _OPENMP
        if(omp_get_thread_num() == 0)
        #endif
        {
            std::cout << "Symetric Chebyshev adaptive kernel" << '\n';
            std::cout << "Compute class: "
                      << typeid(typename ComputeClassDescriptor<FReal>::type).name()
                      << '\n';
        }
    }

    template<class SymbolicData>
    void P2M(typename CellClass::multipole_t* const LeafCell,
             const SymbolicData * const LeafSymbData,
             const ContainerClass* const SourceParticles)
    {
        FReal leafBoxWidth = FBase::BoxWidth / FReal(1 <<  LeafSymbData->getLevel());
        const FPoint<FReal> LeafCellCenter =
            FBase::getCellCenter(LeafSymbData->getCoordinate(), static_cast<int>(LeafSymbData->getLevel()));
        FBase::Interpolator->applyP2M(LeafCellCenter, leafBoxWidth,
                                      LeafCell->get(0), SourceParticles);
    }

    template<class SymbolicData>
    void L2P(const typename CellClass::local_expansion_t* const LeafCell,
             const SymbolicData * const LeafSymbData,
             ContainerClass* const TargetParticles)
    {
        FReal leafBoxWidth = FBase::BoxWidth / FReal(1 <<  LeafSymbData->getLevel());
        const FPoint<FReal> LeafCellCenter =
            FBase::getCellCenter(LeafSymbData->getCoordinate(), static_cast<int>(LeafSymbData->getLevel()));
        // apply Sx and Px (grad Sx)
        FBase::Interpolator->applyL2PTotal(LeafCellCenter, leafBoxWidth,
                                           LeafCell->get(0), TargetParticles);
    }

    template<class SymbolicData>
    void P2L(typename CellClass::local_expansion_t* const local,
             const SymbolicData * const symb,
             const ContainerClass* const particles,
             const SymbolicData * const /*source_symb*/)
    {
        using ComputeClass = typename ComputeClassDescriptor<FReal>::type;
        constexpr std::size_t FRealCount = ComputeClassDescriptor<FReal>::count;

        // Target cell: local
        const FReal localCellWidth(FBase::BoxWidth / FReal(1 << symb->getLevel()));
        const FPoint<FReal> localCellCenter(
            FBase::getCellCenter(symb->getCoordinate(),static_cast<int>(symb->getLevel())));

        // interpolation points of target (X) cell
        FPoint<FReal> X[FBase::nnodes];
        FChebTensor<FReal,ORDER>::setRoots(localCellCenter, localCellWidth, X);

        // Particles attributes
        const ComputeClass * const posX = (const ComputeClass * const) particles->getPositions()[0];
        const ComputeClass * const posY = (const ComputeClass * const) particles->getPositions()[1];
        const ComputeClass * const posZ = (const ComputeClass * const) particles->getPositions()[2];
        const ComputeClass * const physicalValues = (const ComputeClass * const) particles->getPhysicalValues();

        const FReal* pX = particles->getPositions()[0];
        const FReal* pY = particles->getPositions()[1];
        const FReal* pZ = particles->getPositions()[2];
        const FReal* pV = particles->getPhysicalValues();

        // apply P2L
        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
            for (unsigned int m = 0; m<FBase::nnodes; ++m) {
                ComputeClass XX = FMath::ConvertTo<ComputeClass>(X[m].getX());
                ComputeClass XY = FMath::ConvertTo<ComputeClass>(X[m].getY());
                ComputeClass XZ = FMath::ConvertTo<ComputeClass>(X[m].getZ());

                std::size_t idxPart = 0;
                // Compute using vectorization for all but the last array elements
                ComputeClass tmpLocalExp = FMath::Zero<ComputeClass>();
                for (;
                     idxPart < ((particles->getNbParticles())
                                / FRealCount);
                     ++idxPart)
                {
                    tmpLocalExp +=
                        FBase::MatrixKernel->evaluate(
                            XX, XY, XZ,
                            posX[idxPart], posY[idxPart], posZ[idxPart])
                        * physicalValues[idxPart];
                }

                local->get(idxRhs)[m] += FMath::ConvertTo<FReal>(tmpLocalExp);

                // Compute the last array elements one by one if they exist
                if(idxPart < ((particles->getNbParticles() + FRealCount - 1) / FRealCount)) {
                    auto Xx = X[m].getX();
                    auto Xy = X[m].getY();
                    auto Xz = X[m].getZ();

                    for(idxPart = FRealCount * (particles->getNbParticles() / FRealCount);
                        idxPart < static_cast<std::size_t>(particles->getNbParticles());
                        ++idxPart)
                    {
                        local->get(idxRhs)[m] +=
                            FBase::MatrixKernel->evaluate(
                                Xx, Xy, Xz,
                                pX[idxPart], pY[idxPart], pZ[idxPart])
                            * pV[idxPart];
                    }
                }
            }
        }// NVALS
    }

    template<class SymbolicData>
    void M2P(const typename CellClass::multipole_t* const pole,
             const SymbolicData * const symb,
             ContainerClass* const particles,
             const SymbolicData * const /*source_symb*/)
    {
        using ComputeClass = typename ComputeClassDescriptor<FReal>::type;
        constexpr std::size_t FRealCount = ComputeClassDescriptor<FReal>::count;

        // Source cell: pole
        const FReal poleCellWidth(FBase::BoxWidth / FReal(1 << symb->getLevel()));
        const FPoint<FReal> poleCellCenter(
            FBase::getCellCenter(symb->getCoordinate(),static_cast<int>(symb->getLevel())));

        // interpolation points of source (Y) cell
        FPoint<FReal> Y[FBase::nnodes];
        FChebTensor<FReal,ORDER>::setRoots(poleCellCenter, poleCellWidth, Y);

        // read positions
        const ComputeClass* const posX = (const ComputeClass* const)(particles->getPositions()[0]);
        const ComputeClass* const posY = (const ComputeClass* const)(particles->getPositions()[1]);
        const ComputeClass* const posZ = (const ComputeClass* const)(particles->getPositions()[2]);

        // get potential
        ComputeClass* const physVal = (ComputeClass* const)(particles->getPhysicalValues());
        ComputeClass* const potentials = (ComputeClass* const)(particles->getPotentials());
        ComputeClass* const fx = (ComputeClass* const)(particles->getForcesX());
        ComputeClass* const fy = (ComputeClass* const)(particles->getForcesY());
        ComputeClass* const fz = (ComputeClass* const)(particles->getForcesZ());

        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
            // apply M2P
            for (unsigned int n=0; n<FBase::nnodes; ++n){

                ComputeClass  MultipoleExpansion =
                    FMath::ConvertTo<ComputeClass, FReal>(pole->get(idxRhs)[n]);

                ComputeClass YX = FMath::ConvertTo<ComputeClass, FReal>(Y[n].getX());
                ComputeClass YY = FMath::ConvertTo<ComputeClass, FReal>(Y[n].getY());
                ComputeClass YZ = FMath::ConvertTo<ComputeClass, FReal>(Y[n].getZ());

                for(std::size_t idxPart = 0;
                    idxPart < ( (particles->getNbParticles() + FRealCount - 1)
                                / FRealCount);
                    ++idxPart)
                {
                    ComputeClass Kxy[1];
                    ComputeClass dKxy[3];
                    FBase::MatrixKernel->evaluateBlockAndDerivative(
                        posX[idxPart], posY[idxPart], posZ[idxPart],
                        YX, YY, YZ,
                        Kxy,dKxy);

                    potentials[idxPart] += Kxy[0] * MultipoleExpansion;
                    fx[idxPart] += dKxy[0] * physVal[idxPart] * MultipoleExpansion;
                    fy[idxPart] += dKxy[1] * physVal[idxPart] * MultipoleExpansion;
                    fz[idxPart] += dKxy[2] * physVal[idxPart] * MultipoleExpansion;
                }
            }// Particles
        }// NVALS
    }



};

#endif /* FADAPTCHEBKERNEL_HPP */