utestSphericalDirect.cpp 10.4 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

#include "../Src/Utils/FGlobal.hpp"
18

19 20 21
#include "../Src/Containers/FOctree.hpp"
#include "../Src/Containers/FVector.hpp"

22 23
#include "../Src/Kernels/Spherical/FSphericalCell.hpp"
#include "../Src/Kernels/Spherical/FSphericalParticle.hpp"
24

25
#include "../Src/Components/FSimpleLeaf.hpp"
26 27 28 29
#include "../Src/Kernels/Spherical/FSphericalKernel.hpp"
#include "../Src/Kernels/Spherical/FSphericalRotationKernel.hpp"
#include "../Src/Kernels/Spherical/FSphericalBlasKernel.hpp"
#include "../Src/Kernels/Spherical/FSphericalBlockBlasKernel.hpp"
30

31 32 33 34 35 36 37 38
#include "../Src/Files/FFmaBinLoader.hpp"
#include "../Src/Files/FTreeIO.hpp"

#include "../Src/Core/FFmmAlgorithmThread.hpp"
#include "../Src/Core/FFmmAlgorithm.hpp"

#include "FUTester.hpp"

berenger-bramas's avatar
berenger-bramas committed
39
/*
berenger-bramas's avatar
berenger-bramas committed
40
  In this test we compare the spherical fmm results and the direct results.
berenger-bramas's avatar
berenger-bramas committed
41 42
  */

43
/** We need to know the position of the particle in the array */
44
class IndexedParticle : public FSphericalParticle {
45 46 47 48 49 50 51 52 53 54 55 56
    int index;
public:
    IndexedParticle(): index(-1){}

    int getIndex() const{
        return index;
    }
    void setIndex( const int inIndex ){
        index = inIndex;
    }
};

berenger-bramas's avatar
berenger-bramas committed
57 58 59
/** the test class
  *
  */
berenger-bramas's avatar
berenger-bramas committed
60
class TestSphericalDirect : public FUTester<TestSphericalDirect> {
berenger-bramas's avatar
berenger-bramas committed
61 62 63 64
    /** The test method to factorize all the test based on different kernels */
    template <class ParticleClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass,
              class OctreeClass, class FmmClass>
    void RunTest(const bool isBlasKernel){
65
        // Warning in make test the exec dir it Build/UTests
66
        // Load particles
67 68 69 70
        const char* const filename = (sizeof(FReal) == sizeof(float))?
                                        "../../Data/utestDirect.bin.fma.single":
                                        "../../Data/utestDirect.bin.fma.double";
        FFmaBinLoader<ParticleClass> loader(filename);
71 72
        if(!loader.isOpen()){
            Print("Cannot open particles file.");
73
            uassert(false);
74 75 76 77 78
            return;
        }
        Print("Number of particles:");
        Print(loader.getNumberOfParticles());

79 80
        const int NbLevels      = 4;
        const int SizeSubLevels = 2;
berenger-bramas's avatar
berenger-bramas committed
81 82
        const int DevP = 9;
        FSphericalCell::Init(DevP, isBlasKernel);
83

84 85
        // Create octree
        OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
86 87 88 89 90 91 92 93
        ParticleClass* const particles = new ParticleClass[loader.getNumberOfParticles()];
        for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
            loader.fillParticle(particles[idxPart]);
            particles[idxPart].setIndex( idxPart );
            tree.insert(particles[idxPart]);
        }


94
        // Run FMM
95
        Print("Fmm...");
96
        //KernelClass kernels(NbLevels,loader.getBoxWidth());
97
        KernelClass kernels(DevP,NbLevels,loader.getBoxWidth(), loader.getCenterOfBox());
98 99 100
        FmmClass algo(&tree,&kernels);
        algo.execute();

101
        // Run direct computation
102 103 104
        Print("Direct...");
        for(int idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
            for(int idxOther = idxTarget + 1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
105
                kernels.directInteractionMutual(&particles[idxTarget], &particles[idxOther]);
106 107 108
            }
        }

109
        // Compare
110
        Print("Compute Diff...");
111 112
        FMath::FAccurater potentialDiff;
        FMath::FAccurater fx, fy, fz;
113
        { // Check that each particle has been summed with all other
berenger-bramas's avatar
berenger-bramas committed
114
            typename OctreeClass::Iterator octreeIterator(&tree);
115 116 117 118 119 120
            octreeIterator.gotoBottomLeft();

            do{
                typename ContainerClass::BasicIterator leafIter(*octreeIterator.getCurrentListTargets());

                while( leafIter.hasNotFinished() ){
berenger-bramas's avatar
berenger-bramas committed
121
                    const ParticleClass& other = particles[leafIter.data().getIndex()];
122

123
                    potentialDiff.add(other.getPotential(),leafIter.data().getPotential());
124

125
                    fx.add(other.getForces().getX(),leafIter.data().getForces().getX());
126

127
                    fy.add(other.getForces().getY(),leafIter.data().getForces().getY());
128

129
                    fz.add(other.getForces().getZ(),leafIter.data().getForces().getZ());
130 131 132 133 134 135 136 137

                    leafIter.gotoNext();
                }
            } while(octreeIterator.moveRight());
        }

        delete[] particles;

berenger-bramas's avatar
berenger-bramas committed
138
        // Print for information
139
        Print("Potential diff is = ");
140 141
        Print(potentialDiff.getL2Norm());
        Print(potentialDiff.getInfNorm());
142
        Print("Fx diff is = ");
143 144
        Print(fx.getL2Norm());
        Print(fx.getInfNorm());
145
        Print("Fy diff is = ");
146 147
        Print(fy.getL2Norm());
        Print(fy.getInfNorm());
148
        Print("Fz diff is = ");
149 150 151
        Print(fz.getL2Norm());
        Print(fz.getInfNorm());

berenger-bramas's avatar
berenger-bramas committed
152
        // Assert
153
        const FReal MaximumDiff = FReal(0.0001);
154 155 156 157 158 159 160 161
        uassert(potentialDiff.getL2Norm() < MaximumDiff);
        uassert(potentialDiff.getInfNorm() < MaximumDiff);
        uassert(fx.getL2Norm()  < MaximumDiff);
        uassert(fx.getInfNorm() < MaximumDiff);
        uassert(fy.getL2Norm()  < MaximumDiff);
        uassert(fy.getInfNorm() < MaximumDiff);
        uassert(fz.getL2Norm()  < MaximumDiff);
        uassert(fz.getInfNorm() < MaximumDiff);
162 163
    }

berenger-bramas's avatar
berenger-bramas committed
164 165
    /** If memstas is running print the memory used */
    void PostTest() {
166 167 168 169 170 171 172
        if( FMemStats::controler.isUsed() ){
            std::cout << "Memory used at the end " << FMemStats::controler.getCurrentAllocated() << " Bytes (" << FMemStats::controler.getCurrentAllocatedMB() << "MB)\n";
            std::cout << "Max memory used " << FMemStats::controler.getMaxAllocated() << " Bytes (" << FMemStats::controler.getMaxAllocatedMB() << "MB)\n";
            std::cout << "Total memory used " << FMemStats::controler.getTotalAllocated() << " Bytes (" << FMemStats::controler.getTotalAllocatedMB() << "MB)\n";
        }
    }

berenger-bramas's avatar
berenger-bramas committed
173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210
    ///////////////////////////////////////////////////////////
    // The tests!
    ///////////////////////////////////////////////////////////

    /** Classic */
    void TestSpherical(){
        typedef IndexedParticle         ParticleClass;
        typedef FSphericalCell            CellClass;
        typedef FVector<ParticleClass>  ContainerClass;

        typedef FSphericalKernel<ParticleClass, CellClass, ContainerClass >          KernelClass;

        typedef FSimpleLeaf<ParticleClass, ContainerClass >                     LeafClass;
        typedef FOctree<ParticleClass, CellClass, ContainerClass , LeafClass >  OctreeClass;

        typedef FFmmAlgorithm<OctreeClass, ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;

        RunTest<ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass,
                OctreeClass, FmmClass>(false);
    }

    /** Rotation */
    void TestRotation(){
        typedef IndexedParticle         ParticleClass;
        typedef FSphericalCell            CellClass;
        typedef FVector<ParticleClass>  ContainerClass;

        typedef FSphericalRotationKernel<ParticleClass, CellClass, ContainerClass >          KernelClass;

        typedef FSimpleLeaf<ParticleClass, ContainerClass >                     LeafClass;
        typedef FOctree<ParticleClass, CellClass, ContainerClass , LeafClass >  OctreeClass;

        typedef FFmmAlgorithm<OctreeClass, ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;

        RunTest<ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass,
                OctreeClass, FmmClass>(false);
    }

211
#ifdef SCALFMM_USE_BLAS
berenger-bramas's avatar
berenger-bramas committed
212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244
    /** Blas */
    void TestSphericalBlas(){
        typedef IndexedParticle         ParticleClass;
        typedef FSphericalCell            CellClass;
        typedef FVector<ParticleClass>  ContainerClass;

        typedef FSphericalBlasKernel<ParticleClass, CellClass, ContainerClass >          KernelClass;

        typedef FSimpleLeaf<ParticleClass, ContainerClass >                     LeafClass;
        typedef FOctree<ParticleClass, CellClass, ContainerClass , LeafClass >  OctreeClass;

        typedef FFmmAlgorithm<OctreeClass, ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;

        RunTest<ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass,
                OctreeClass, FmmClass>(true);
    }

    /** Block blas */
    void TestSphericalBlockBlas(){
        typedef IndexedParticle         ParticleClass;
        typedef FSphericalCell            CellClass;
        typedef FVector<ParticleClass>  ContainerClass;

        typedef FSphericalBlockBlasKernel<ParticleClass, CellClass, ContainerClass >          KernelClass;

        typedef FSimpleLeaf<ParticleClass, ContainerClass >                     LeafClass;
        typedef FOctree<ParticleClass, CellClass, ContainerClass , LeafClass >  OctreeClass;

        typedef FFmmAlgorithm<OctreeClass, ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;

        RunTest<ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass,
                OctreeClass, FmmClass>(true);
    }
245
#endif
berenger-bramas's avatar
berenger-bramas committed
246 247 248 249 250 251

    ///////////////////////////////////////////////////////////
    // Set the tests!
    ///////////////////////////////////////////////////////////

    /** set test */
252
    void SetTests(){
berenger-bramas's avatar
berenger-bramas committed
253 254
        AddTest(&TestSphericalDirect::TestSpherical,"Test Spherical Kernel");
        AddTest(&TestSphericalDirect::TestRotation,"Test Rotation Spherical Kernel");
255
#ifdef SCALFMM_USE_BLAS
berenger-bramas's avatar
berenger-bramas committed
256 257
        AddTest(&TestSphericalDirect::TestSphericalBlas,"Test Spherical Blas Kernel");
        AddTest(&TestSphericalDirect::TestSphericalBlockBlas,"Test Spherical Block Blas Kernel");
258
#endif
259 260 261 262 263
    }
};


// You must do this
berenger-bramas's avatar
berenger-bramas committed
264
TestClass(TestSphericalDirect)
265 266


267