testSphericalDebug.cpp 18 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger 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".
// ===================================================================================

17
18
// Keep in private GIT
// @SCALFMM_PRIVATE
19
// @FUSE_BLAS
20
#define DEBUG_SPHERICAL_M2L
COULAUD Olivier's avatar
COULAUD Olivier committed
21
22
#define  BLAS_SPHERICAL_COMPRESS
#define  BLAS_M2L_P
23
24
25
26
27
28
29
30
#include <iostream>

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

#include "../Src/Containers/FVector.hpp"

#include "../Src/Kernels/Spherical/FSphericalCell.hpp"
#include "../Src/Kernels/P2P/FP2PParticleContainerIndexed.hpp"
31
#include "../Src/Kernels/Interpolation/FInterpMatrixKernel.hpp"
32
33
34
35
36
37
38
39
40
41

#include "../Src/Components/FSimpleLeaf.hpp"
#include "../Src/Kernels/Spherical/FSphericalKernel.hpp"
#include "../Src/Kernels/Spherical/FSphericalRotationKernel.hpp"
#include "../Src/Kernels/Spherical/FSphericalBlasKernel.hpp"
#include "../Src/Kernels/Spherical/FSphericalBlockBlasKernel.hpp"

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

42
43
#include "../../Src/Utils/FParameterNames.hpp"

44
45
46
47
48
49
50
51
52
53
54
#include "../UTests/FUTester.hpp"

/*
  In this test we compare the spherical fmm results and the direct results.
  */

/** the test class
  *
  */
class TestSphericalDirect : public FUTester<TestSphericalDirect> {
    /** The test method to factorize all the test based on different kernels */
55
56
    template < class FReal, class CellClass, class ContainerClass, class KernelClass, class LeafClass,
               class OctreeClass, class FmmClass>
57
    void RunTest(const bool isBlasKernel){
COULAUD Olivier's avatar
COULAUD Olivier committed
58
        const int DevP = 2;
59
60
        // Warning in make test the exec dir it Build/UTests
        // Load particles
61
        const FSize nbParticles = 2;
62
        const FReal boxWidth = 1.0;
63
        const FPoint<FReal> boxCenter(boxWidth/2.0,boxWidth/2.0,boxWidth/2.0);
64
65

        struct TestParticle{
66
            FPoint<FReal> position;
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
            FReal forces[3];
            FReal physicalValue;
            FReal potential;
        };

        const int NbLevels      = 3;
        const int SizeSubLevels = 2;

        const int dimGrid = (1 << (NbLevels-1));
        const FReal dimLeaf = (boxWidth/FReal(dimGrid));
        const FReal quarterDimLeaf = (dimLeaf/4.0);

        Print("dimGrid:");
        Print(dimGrid);
        Print("dimLeaf:");
        Print(dimLeaf);
        Print("quarterDimLeaf:");
        Print(quarterDimLeaf);

86
        FSphericalCell<FReal>::Init(DevP, isBlasKernel);
87
88

        TestParticle* const particles = new TestParticle[nbParticles];
89
        particles[0].position = FPoint<FReal>(quarterDimLeaf, quarterDimLeaf, quarterDimLeaf);
90
        particles[0].physicalValue = 1;
91
        //        particles[1].position = FPoint<FReal>(2*quarterDimLeaf, quarterDimLeaf, quarterDimLeaf);
92
93
94
95
        particles[1].physicalValue = -1;

        Print("Number of particles:");
        Print(nbParticles);
BRAMAS Berenger's avatar
BRAMAS Berenger committed
96
        //int idxLeafY = 0,idxLeafZ = 0 ;
97
        std::cout << "\n ------  Loop starts ---"<< std::endl ;
98
99
        for(int idxLeafX = 2 ; idxLeafX < dimGrid ; ++idxLeafX){
            /*for(int idxLeafY = 0 ; idxLeafY < dimGrid ; ++idxLeafY)*/{
100
                /*  for(int idxLeafZ = 0 ; idxLeafZ < dimGrid ; ++idxLeafZ)*/{
101
102
                    //std::cout << "Shift : " << idxLeafX << " " << idxLeafY << " " << idxLeafZ << std::endl;

103
                    /*  particles[1].position = FPoint<FReal>(FReal(idxLeafX)*dimLeaf + 3*quarterDimLeaf,
104
                                                   FReal(idxLeafY)*dimLeaf + quarterDimLeaf,
105
                                                   FReal(idxLeafZ)*dimLeaf + quarterDimLeaf);*/
106
107
108
109
110
111
112
                    //                   particles[1].position = FPoint<FReal>(FReal(idxLeafX)*dimLeaf + quarterDimLeaf,
                    //                               quarterDimLeaf,
                    //                               quarterDimLeaf);
                    particles[1].position = FPoint<FReal>(FReal(idxLeafX)*dimLeaf + 2*quarterDimLeaf,
                                                   quarterDimLeaf,
                                                   quarterDimLeaf);
                    /*     particles[1].position = FPoint<FReal>(FReal(idxLeafX)*dimLeaf + 2*quarterDimLeaf,
113
114
115
116
117
                               quarterDimLeaf,
                               quarterDimLeaf);*/

                    // Create octree
                    OctreeClass tree(NbLevels, SizeSubLevels, boxWidth, boxCenter);
118
                    for(FSize idxPart = 0 ; idxPart < nbParticles ; ++idxPart){
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
                        // put in tree
                        tree.insert(particles[idxPart].position, idxPart, particles[idxPart].physicalValue);
                        // get copy
                        particles[idxPart].potential = 0.0;
                        particles[idxPart].forces[0] = 0.0;
                        particles[idxPart].forces[1] = 0.0;
                        particles[idxPart].forces[2] = 0.0;

                        std::cout << idxPart << " " << particles[idxPart].position << std::endl;
                    }


                    // Run FMM
                    //Print("Fmm...");
                    KernelClass kernels(DevP,NbLevels,boxWidth, boxCenter);
                    FmmClass algo(&tree,&kernels);
                    algo.execute();

                    // Run direct computation
                    //Print("Direct...");
                    for(int idxTarget = 0 ; idxTarget < nbParticles ; ++idxTarget){
                        for(int idxOther = idxTarget + 1 ; idxOther < nbParticles ; ++idxOther){
141
                            FP2PR::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
142
143
144
145
146
147
148
                                                   particles[idxTarget].position.getZ(),particles[idxTarget].physicalValue,
                                                   &particles[idxTarget].forces[0],&particles[idxTarget].forces[1],
                                    &particles[idxTarget].forces[2],&particles[idxTarget].potential,
                                    particles[idxOther].position.getX(), particles[idxOther].position.getY(),
                                    particles[idxOther].position.getZ(),particles[idxOther].physicalValue,
                                    &particles[idxOther].forces[0],&particles[idxOther].forces[1],
                                    &particles[idxOther].forces[2],&particles[idxOther].potential);
149
150
151
152
153
                        }
                    }

                    // Compare
                    //Print("Compute Diff...");
154
155
                    FMath::FAccurater<FReal> potentialDiff;
                    FMath::FAccurater<FReal> fx, fy, fz;
156
157
158
159
160
161
162
                    { // Check that each particle has been summed with all other

                        tree.forEachLeaf([&](LeafClass* leaf){
                            const FReal*const potentials = leaf->getTargets()->getPotentials();
                            const FReal*const forcesX = leaf->getTargets()->getForcesX();
                            const FReal*const forcesY = leaf->getTargets()->getForcesY();
                            const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
163
164
                            const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
                            const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
165

166
167
                            for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                                const FSize indexPartOrig = indexes[idxPart];
168
169
170
171
172
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
                                potentialDiff.add(particles[indexPartOrig].potential,potentials[idxPart]);
                                fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]);
                                fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
                                fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);

                                std::cout <<"indx: " << indexPartOrig << " Potential  Direct " << particles[indexPartOrig].potential<<" FMM "<<potentials[idxPart] << std::endl;
                                std::cout <<"indx: " << indexPartOrig << " Fx Direct " << particles[indexPartOrig].forces[0]<<" FMM "<<forcesX[idxPart] << std::endl;
                                std::cout <<"indx: " << indexPartOrig << " Fy Direct " << particles[indexPartOrig].forces[1]<<" FMM "<<forcesY[idxPart] << std::endl;
                                std::cout <<"indx: " << indexPartOrig << " Fz Direct " << particles[indexPartOrig].forces[2]<<" FMM "<<forcesZ[idxPart] << std::endl;
                                printf("Printf : forces x %e y %e z %e\n", forcesX[idxPart],forcesY[idxPart],forcesZ[idxPart]);//TODO delete
                            }
                        });

                        tree.forEachCell([&](CellClass* cell){

                            std::cout << "cell " << cell->getMortonIndex() << "\n   Multipole: \n";
                            int index_j_k = 0;
                            for (int j = 0 ; j <= DevP ; ++j ){
                                std::cout <<"[" << j << "] ----- level\n";
                                for (int k=0; k<=j ;++k, ++index_j_k){
                                    std::cout << "[" << k << "] ( " << cell->getMultipole()[index_j_k].getReal() << " , " << cell->getMultipole()[index_j_k].getImag() << " )   ";
                                }
                                std::cout << "\n";
                            }
                            std::cout << "\n";
                            std::cout << "   Local:\n";
                            index_j_k = 0;
                            for (int j = 0 ; j <= DevP ; ++j ){
                                std::cout <<"[" << j << "] ----- level \n";
                                for (int k=0; k<=j ;++k, ++index_j_k){
198
                                    std::cout << "[" << k << "] ( " << cell->getLocal()[index_j_k].getReal() << " , " << cell->getLocal()[index_j_k].getImag() << " )   ";
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
                                }
                                std::cout << "\n";
                            }
                            std::cout << "\n\n";
                        });
                    }

                    // Print for information
                    Print("Potential diff is = ");
                    Print(potentialDiff.getL2Norm());
                    Print(potentialDiff.getInfNorm());
                    Print("Fx diff is = ");
                    Print(fx.getL2Norm());
                    Print(fx.getInfNorm());
                    Print("Fy diff is = ");
                    Print(fy.getL2Norm());
                    Print(fy.getInfNorm());
                    Print("Fz diff is = ");
                    Print(fz.getL2Norm());
                    Print(fz.getInfNorm());

                    // Assert
                    // Assert
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
                    const FReal MaximumDiff = FReal(0.0001);
                    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);
                    //                   const FReal MaximumDiff = FReal(0.0001);
                    //                    if(fx.getL2Norm() > MaximumDiff || fx.getInfNorm() > MaximumDiff){
                    //                        std::cout << "Error in X " << fx.getL2Norm() << " " << fx.getInfNorm() << std::endl;
                    //                    }
                    //                    if(fy.getL2Norm() > MaximumDiff || fy.getInfNorm() > MaximumDiff){
                    //                        std::cout << "Error in Y " << fy.getL2Norm() << " " << fy.getInfNorm() << std::endl;
                    //                    }
                    //                    if(fz.getL2Norm() > MaximumDiff || fz.getInfNorm() > MaximumDiff){
                    //                        std::cout << "Error in Z " << fz.getL2Norm() << " " << fz.getInfNorm() << std::endl;
                    //                    }
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
                }
            }
        }

        delete[] particles;
    }

    /** If memstas is running print the memory used */
    void PostTest() {
        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";
        }
    }

    ///////////////////////////////////////////////////////////
    // The tests!
    ///////////////////////////////////////////////////////////

    /** Classic */
    void TestSpherical(){
263
264
        typedef double FReal;

265
        typedef FSphericalCell<FReal>            CellClass;
266
        typedef FP2PParticleContainerIndexed<FReal>  ContainerClass;
267

268
        typedef FSphericalKernel< FReal, CellClass, ContainerClass >          KernelClass;
269

270
271
        typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
        typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
272
273

        typedef FFmmAlgorithm<OctreeClass,  CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
274
        std::cout << std::endl << std::endl << std::endl
275
276
277
                  << " $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<<std::endl
                  << " $$$$$$$$$                                                    TestSpherical                                             $$$$$$$$$$$$$$$$"<<std::endl
                  << " $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<<std::endl;
278

279
        RunTest< FReal,  CellClass, ContainerClass, KernelClass, LeafClass,
280
281
282
283
284
                OctreeClass, FmmClass>(false);
    }

    /** Rotation */
    void TestRotation(){
285
286
        typedef double FReal;

287
        typedef FSphericalCell<FReal>            CellClass;
288
        typedef FP2PParticleContainerIndexed<FReal>  ContainerClass;
289

290
        typedef FSphericalRotationKernel<  FReal, CellClass, ContainerClass >          KernelClass;
291

292
293
        typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
        typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
294
295
296

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

297
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
298
299
300
                OctreeClass, FmmClass>(false);
    }

301
#ifdef SCALFMM_USE_BLAS
302
303
    /** Blas */
    void TestSphericalBlas(){
304
305
        typedef double FReal;

306
        typedef FSphericalCell<FReal>            CellClass;
307
        typedef FP2PParticleContainerIndexed<FReal>  ContainerClass;
308

309
        typedef FSphericalBlasKernel<FReal, CellClass, ContainerClass >          KernelClass;
310

311
312
        typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
        typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
313
314
315

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

316
        RunTest< FReal,  CellClass, ContainerClass, KernelClass, LeafClass,
317
318
319
320
321
                OctreeClass, FmmClass>(true);
    }

    /** Block blas */
    void TestSphericalBlockBlas(){
322
323
        typedef double FReal;

324
        typedef FSphericalCell<FReal>            CellClass;
325
        typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
326

327
        typedef FSphericalBlockBlasKernel< FReal, CellClass, ContainerClass >          KernelClass;
328

329
330
        typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
        typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
331
332

        typedef FFmmAlgorithm<OctreeClass,  CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
333
334
        //
        std::cout << std::endl << std::endl << std::endl
335
336
337
                  << " $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<<std::endl
                  << " $$$$$$$$$                                                    TestSphericalBlockBlas                                $$$$$$$$$$$$$$$$"<<std::endl
                  << " $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<<std::endl;
338
        RunTest<  FReal, CellClass, ContainerClass, KernelClass, LeafClass,
339
340
341
342
343
344
345
346
347
348
                OctreeClass, FmmClass>(true);
    }
#endif

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

    /** set test */
    void SetTests(){
349
        //    AddTest(&TestSphericalDirect::TestSpherical,"Test Spherical Kernel");
350
        //AddTest(&TestSphericalDirect::TestRotation,"Test Rotation Spherical Kernel");
351
#ifdef SCALFMM_USE_BLAS
352
353
        AddTest(&TestSphericalDirect::TestSphericalBlas,"Test Spherical Blas Kernel");
        //        AddTest(&TestSphericalDirect::TestSphericalBlockBlas,"Test Spherical Block Blas Kernel");
354
355
356
357
358
359
360
361
362
363
#endif
    }
};


// You must do this
TestClass(TestSphericalDirect)