FAdaptChebKernel.hpp 7.26 KB
Newer Older
1 2 3 4
#ifndef FADAPTCHEBKERNEL_HPP
#define FADAPTCHEBKERNEL_HPP

#include "Kernels/Chebyshev/FChebSymKernel.hpp"
5
#include "InastempCompileConfig.h"
6

7 8 9 10
#ifdef _OPENMP
#include <omp.h>
#endif

11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
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() {
31 32 33 34 35

        #ifdef _OPENMP
        if(omp_get_thread_num() == 0)
        #endif
        {
36
            std::cout << "Symetric Chebyshev adaptive kernel" << '\n';
37
        }
38 39
    }

40 41 42
    template<class SymbolicData>
    void P2M(typename CellClass::multipole_t* const LeafCell,
             const SymbolicData * const LeafSymbData,
43 44
             const ContainerClass* const SourceParticles)
    {
45
        FReal leafBoxWidth = FBase::BoxWidth / FReal(1 <<  LeafSymbData->getLevel());
46
        const FPoint<FReal> LeafCellCenter =
47
            FBase::getCellCenter(LeafSymbData->getCoordinate(), static_cast<int>(LeafSymbData->getLevel()));
48
        FBase::Interpolator->applyP2M(LeafCellCenter, leafBoxWidth,
49
                                      LeafCell->get(0), SourceParticles);
50 51
    }

52 53 54
    template<class SymbolicData>
    void L2P(const typename CellClass::local_expansion_t* const LeafCell,
             const SymbolicData * const LeafSymbData,
55 56
             ContainerClass* const TargetParticles)
    {
57
        FReal leafBoxWidth = FBase::BoxWidth / FReal(1 <<  LeafSymbData->getLevel());
58
        const FPoint<FReal> LeafCellCenter =
59
            FBase::getCellCenter(LeafSymbData->getCoordinate(), static_cast<int>(LeafSymbData->getLevel()));
60 61
        // apply Sx and Px (grad Sx)
        FBase::Interpolator->applyL2PTotal(LeafCellCenter, leafBoxWidth,
62
                                           LeafCell->get(0), TargetParticles);
63 64
    }

65 66 67 68 69 70
    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*/)
    {
71 72
        using ComputeClass = InaVecBestType<FReal>;
        constexpr int FRealCount = ComputeClass::VecLength;
73 74

        // Target cell: local
75
        const FReal localCellWidth(FBase::BoxWidth / FReal(1 << symb->getLevel()));
76
        const FPoint<FReal> localCellCenter(
77
            FBase::getCellCenter(symb->getCoordinate(),static_cast<int>(symb->getLevel())));
78 79 80 81 82 83

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

        // Particles attributes
84 85 86 87
        const FReal * const posX = particles->getPositions()[0];
        const FReal * const posY = particles->getPositions()[1];
        const FReal * const posZ = particles->getPositions()[2];
        const FReal * const physicalValues = particles->getPhysicalValues();
88

89 90 91 92
//        const FReal* pX = particles->getPositions()[0];
//        const FReal* pY = particles->getPositions()[1];
//        const FReal* pZ = particles->getPositions()[2];
//        const FReal* pV = particles->getPhysicalValues();
93

94 95
        // apply P2L
        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
96
            for (unsigned int m = 0; m<FBase::nnodes; ++m) {
Berenger Bramas's avatar
Berenger Bramas committed
97 98 99
                ComputeClass XX = ComputeClass(X[m].getX());
                ComputeClass XY = ComputeClass(X[m].getY());
                ComputeClass XZ = ComputeClass(X[m].getZ());
100

101
                // Compute using vectorization for all but the last array elements
102 103
                ComputeClass tmpLocalExp = ComputeClass::GetZero();
                for (std::size_t idxPart = 0 ; idxPart < particles->getNbParticles() ; idxPart += FRealCount)
104 105 106 107
                {
                    tmpLocalExp +=
                        FBase::MatrixKernel->evaluate(
                            XX, XY, XZ,
108
                                ComputeClass(&posX[idxPart]), ComputeClass(&posY[idxPart]), ComputeClass(&posZ[idxPart]))
109 110 111
                        * physicalValues[idxPart];
                }

112
                local->get(idxRhs)[m] += (tmpLocalExp.horizontalSum());
113 114 115 116
            }
        }// NVALS
    }

117 118 119 120 121 122
    template<class SymbolicData>
    void M2P(const typename CellClass::multipole_t* const pole,
             const SymbolicData * const symb,
             ContainerClass* const particles,
             const SymbolicData * const /*source_symb*/)
    {
123
        using ComputeClass = InaVecBestType<FReal>;
124
        constexpr int FRealCount = ComputeClass::VecLength;
125 126

        // Source cell: pole
127
        const FReal poleCellWidth(FBase::BoxWidth / FReal(1 << symb->getLevel()));
128
        const FPoint<FReal> poleCellCenter(
129
            FBase::getCellCenter(symb->getCoordinate(),static_cast<int>(symb->getLevel())));
130 131 132 133 134 135

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

        // read positions
136 137 138
        const FReal* const posX = (particles->getPositions()[0]);
        const FReal* const posY = (particles->getPositions()[1]);
        const FReal* const posZ = (particles->getPositions()[2]);
139 140

        // get potential
141 142 143 144 145
        FReal* const physVal = (particles->getPhysicalValues());
        FReal* const potentials = (particles->getPotentials());
        FReal* const fx = (particles->getForcesX());
        FReal* const fy = (particles->getForcesY());
        FReal* const fz = (particles->getForcesZ());
146 147 148 149 150 151

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

                ComputeClass  MultipoleExpansion =
Berenger Bramas's avatar
Berenger Bramas committed
152
                    ComputeClass(pole->get(idxRhs)[n]);
153

Berenger Bramas's avatar
Berenger Bramas committed
154 155 156
                ComputeClass YX = ComputeClass(Y[n].getX());
                ComputeClass YY = ComputeClass(Y[n].getY());
                ComputeClass YZ = ComputeClass(Y[n].getZ());
157 158

                for(std::size_t idxPart = 0;
159 160
                    idxPart < particles->getNbParticles();
                    idxPart += FRealCount)
161 162 163 164
                {
                    ComputeClass Kxy[1];
                    ComputeClass dKxy[3];
                    FBase::MatrixKernel->evaluateBlockAndDerivative(
165 166 167
                                ComputeClass(&posX[idxPart]),
                                ComputeClass(&posY[idxPart]),
                                ComputeClass(&posZ[idxPart]),
168 169 170
                        YX, YY, YZ,
                        Kxy,dKxy);

171 172 173 174
                    (ComputeClass(&potentials[idxPart]) + Kxy[0] * MultipoleExpansion).storeInArray(&potentials[idxPart]);
                    (ComputeClass(&fx[idxPart]) + dKxy[0] * physVal[idxPart] * MultipoleExpansion).storeInArray(&fx[idxPart]);
                    (ComputeClass(&fy[idxPart]) + dKxy[1] * physVal[idxPart] * MultipoleExpansion).storeInArray(&fy[idxPart]);
                    (ComputeClass(&fz[idxPart]) + dKxy[2] * physVal[idxPart] * MultipoleExpansion).storeInArray(&fz[idxPart]);
175 176 177 178 179 180 181 182 183 184
                }
            }// Particles
        }// NVALS
    }



};

#endif /* FADAPTCHEBKERNEL_HPP */