FTestKernels.hpp 7.86 KB
Newer Older
1
// ===================================================================================
2 3 4 5 6 7 8 9 10 11 12 13 14
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.  
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info". 
// "http://www.gnu.org/licenses".
15
// ===================================================================================
16 17
#ifndef FTESTKERNELS_HPP
#define FTESTKERNELS_HPP
18

19 20 21 22 23

#include <iostream>

#include "FAbstractKernels.hpp"
#include "../Containers/FOctree.hpp"
24
#include "../Utils/FGlobal.hpp"
25

26

BRAMAS Berenger's avatar
BRAMAS Berenger committed
27 28 29 30 31 32 33

/**
 * @example testFmmAlgorithm
 * It illustrates how the FTestKernels to validate a FMM core algorithm.
 */


34 35 36
/**
* @author Berenger Bramas (berenger.bramas@inria.fr)
* @class AbstractKernels
BRAMAS Berenger's avatar
BRAMAS Berenger committed
37
* @brief This kernels is a virtual kernels to validate that the fmm algorithm is correctly done on particles.
38
*
39
* It should use FTestCell and FTestParticle.
BRAMAS Berenger's avatar
BRAMAS Berenger committed
40
* A the end of a computation, the particles should host then number of particles in the simulation (-1).
41
*/
42 43
template< class CellClass, class ContainerClass>
class FTestKernels  : public FAbstractKernels<CellClass,ContainerClass> {
44 45 46 47 48
public:
    /** Default destructor */
    virtual ~FTestKernels(){
    }

berenger-bramas's avatar
berenger-bramas committed
49
    /** Before upward */
berenger-bramas's avatar
berenger-bramas committed
50
    void P2M(CellClass* const pole, const ContainerClass* const particles) {
51
        // the pole represents all particles under
52
        pole->setDataUp(pole->getDataUp() + particles->getNbParticles());
53
    }
berenger-bramas's avatar
berenger-bramas committed
54

berenger-bramas's avatar
berenger-bramas committed
55
    /** During upward */
56
    void M2M(CellClass* const FRestrict pole, const CellClass *const FRestrict *const FRestrict child, const int /*level*/) {
57 58 59 60 61 62 63
        // A parent represents the sum of the child
        for(int idx = 0 ; idx < 8 ; ++idx){
            if(child[idx]){
                pole->setDataUp(pole->getDataUp() + child[idx]->getDataUp());
            }
        }
    }
berenger-bramas's avatar
berenger-bramas committed
64

berenger-bramas's avatar
berenger-bramas committed
65
    /** Before Downward */
66
    void M2L(CellClass* const FRestrict pole, const CellClass* distantNeighbors[343], const int /*size*/, const int /*level*/) {
67
        // The pole is impacted by what represent other poles
68
        for(int idx = 0 ; idx < 343 ; ++idx){
69 70 71
            if(distantNeighbors[idx]){
                pole->setDataDown(pole->getDataDown() + distantNeighbors[idx]->getDataUp());
            }
72 73
        }
    }
berenger-bramas's avatar
berenger-bramas committed
74

berenger-bramas's avatar
berenger-bramas committed
75
    /** During Downward */
76
    void L2L(const CellClass*const FRestrict local, CellClass* FRestrict *const FRestrict child, const int /*level*/) {
77 78 79 80 81 82
        // Each child is impacted by the father
        for(int idx = 0 ; idx < 8 ; ++idx){
            if(child[idx]){
                child[idx]->setDataDown(local->getDataDown() + child[idx]->getDataDown());
            }
        }
83

84
    }
berenger-bramas's avatar
berenger-bramas committed
85

berenger-bramas's avatar
berenger-bramas committed
86
    /** After Downward */
berenger-bramas's avatar
berenger-bramas committed
87
    void L2P(const CellClass* const  local, ContainerClass*const particles){
88
        // The particles is impacted by the parent cell      
89 90
        long long int*const particlesAttributes = particles->getDataDown();
        for(int idxPart = 0 ; idxPart < particles->getNbParticles() ; ++idxPart){
91
            particlesAttributes[idxPart] += local->getDataDown();
92 93
        }
    }
berenger-bramas's avatar
berenger-bramas committed
94 95


berenger-bramas's avatar
berenger-bramas committed
96
    /** After Downward */
berenger-bramas's avatar
berenger-bramas committed
97 98 99
    void P2P(const FTreeCoordinate& ,
                 ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
                 ContainerClass* const directNeighborsParticles[27], const int ){
berenger-bramas's avatar
berenger-bramas committed
100
        // Each particles targeted is impacted by the particles sources
101
        long long int inc = sources->getNbParticles();
berenger-bramas's avatar
berenger-bramas committed
102 103 104
        if(targets == sources){
            inc -= 1;
        }
berenger-bramas's avatar
berenger-bramas committed
105 106
        for(int idx = 0 ; idx < 27 ; ++idx){
            if( directNeighborsParticles[idx] ){
107
                inc += directNeighborsParticles[idx]->getNbParticles();
berenger-bramas's avatar
berenger-bramas committed
108
            }
berenger-bramas's avatar
berenger-bramas committed
109 110
        }

111 112
        long long int*const particlesAttributes = targets->getDataDown();
        for(int idxPart = 0 ; idxPart < targets->getNbParticles() ; ++idxPart){
113
            particlesAttributes[idxPart] += inc;
114 115 116 117 118 119 120 121 122 123 124
        }
    }

    /** After Downward */
    void P2PRemote(const FTreeCoordinate& ,
                 ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
                 ContainerClass* const directNeighborsParticles[27], const int ){
        // Each particles targeted is impacted by the particles sources
        long long int inc = 0;
        for(int idx = 0 ; idx < 27 ; ++idx){
            if( directNeighborsParticles[idx] ){
125
                inc += directNeighborsParticles[idx]->getNbParticles();
126 127 128
            }
        }

129 130
        long long int*const particlesAttributes = targets->getDataDown();
        for(int idxPart = 0 ; idxPart < targets->getNbParticles() ; ++idxPart){
131
            particlesAttributes[idxPart] += inc;
berenger-bramas's avatar
berenger-bramas committed
132 133
        }
    }
134 135 136 137 138 139
};


/** This function test the octree to be sure that the fmm algorithm
  * has worked completly.
  */
140
template< class OctreeClass, class CellClass, class ContainerClass, class LeafClass>
berenger-bramas's avatar
berenger-bramas committed
141
void ValidateFMMAlgo(OctreeClass* const tree){
142
    std::cout << "Check Result\n";
143
    const int TreeHeight = tree->getHeight();
144
    long long int NbPart = 0;
145
    { // Check that each particle has been summed with all other
146 147 148 149
        tree->forEachCellLeaf([&](CellClass* cell, LeafClass* leaf){
            if(cell->getDataUp() != leaf->getSrc()->getNbParticles() ){
                    std::cout << "Problem P2M : " << cell->getDataUp() <<
                                 " (should be " << leaf->getSrc()->getNbParticles() << ")\n";
150
            }
151 152
            NbPart += leaf->getSrc()->getNbParticles();
        });
153 154
    }
    { // Ceck if there is number of NbPart summed at level 1
berenger-bramas's avatar
berenger-bramas committed
155
        typename OctreeClass::Iterator octreeIterator(tree);
156
        octreeIterator.moveDown();
157
        long long int res = 0;
158 159 160 161 162 163 164 165
        do{
            res += octreeIterator.getCurrentCell()->getDataUp();
        } while(octreeIterator.moveRight());
        if(res != NbPart){
            std::cout << "Problem M2M at level 1 : " << res << "\n";
        }
    }
    { // Ceck if there is number of NbPart summed at level 1
berenger-bramas's avatar
berenger-bramas committed
166
        typename OctreeClass::Iterator octreeIterator(tree);
167 168
        octreeIterator.gotoBottomLeft();
        for(int idxLevel = TreeHeight - 1 ; idxLevel > 1 ; --idxLevel ){
169
            long long int res = 0;
170 171 172 173 174 175 176 177 178 179 180
            do{
                res += octreeIterator.getCurrentCell()->getDataUp();
            } while(octreeIterator.moveRight());
            if(res != NbPart){
                std::cout << "Problem M2M at level " << idxLevel << " : " << res << "\n";
            }
            octreeIterator.moveUp();
            octreeIterator.gotoLeft();
        }
    }
    { // Check that each particle has been summed with all other
berenger-bramas's avatar
berenger-bramas committed
181
        typename OctreeClass::Iterator octreeIterator(tree);
182 183 184
        octreeIterator.gotoBottomLeft();
        do{
            const bool isUsingTsm = (octreeIterator.getCurrentListTargets() != octreeIterator.getCurrentListSrc());
185 186 187 188
            const long long int* dataDown = octreeIterator.getCurrentListTargets()->getDataDown();
            for(int idxPart = 0 ; idxPart < octreeIterator.getCurrentListTargets()->getNbParticles() ; ++idxPart){
                if( (!isUsingTsm && dataDown[idxPart] != NbPart - 1) ||
                    (isUsingTsm && dataDown[idxPart] != NbPart) ){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
189 190 191
                    std::cout << "Problem L2P + P2P : " << dataDown[idxPart] << ", " <<
                                 " NbPart : " << NbPart << ", " <<
                                 " ( Index " << octreeIterator.getCurrentGlobalIndex() << ")\n";
192 193 194 195 196 197 198 199 200 201 202 203
                }
            }
        } while(octreeIterator.moveRight());
    }

    std::cout << "Done\n";
}



#endif //FTESTKERNELS_HPP

204