utestSphericalDirect.cpp 14.8 KB
Newer Older
1
// See LICENCE file at project root
2

3

4
#include "Utils/FGlobal.hpp"
5
#include "Utils/FTic.hpp"
6

7
8
#include "Containers/FOctree.hpp"
#include "Containers/FVector.hpp"
9

10
11
#include "Kernels/Spherical/FSphericalCell.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
12

13
14
15
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/Spherical/FSphericalKernel.hpp"
#include "Kernels/Spherical/FSphericalRotationKernel.hpp"
16
#ifdef SCALFMM_USE_BLAS
17
18
#include "Kernels/Spherical/FSphericalBlasKernel.hpp"
#include "Kernels/Spherical/FSphericalBlockBlasKernel.hpp"
19
#endif
20
#include "Files/FFmaGenericLoader.hpp"
21

22
#include "Core/FFmmAlgorithm.hpp"
23
24
25

#include "FUTester.hpp"

berenger-bramas's avatar
berenger-bramas committed
26
/*
berenger-bramas's avatar
berenger-bramas committed
27
  In this test we compare the spherical fmm results and the direct results.
COULAUD Olivier's avatar
COULAUD Olivier committed
28
 */
berenger-bramas's avatar
berenger-bramas committed
29

berenger-bramas's avatar
berenger-bramas committed
30
/** the test class
COULAUD Olivier's avatar
COULAUD Olivier committed
31
32
 *
 */
berenger-bramas's avatar
berenger-bramas committed
33
class TestSphericalDirect : public FUTester<TestSphericalDirect> {
COULAUD Olivier's avatar
COULAUD Olivier committed
34
	/** The test method to factorize all the test based on different kernels */
35
    template < class FReal, class CellClass, class ContainerClass, class KernelClass, class LeafClass,
36
		   class OctreeClass, class FmmClass, int ORDER>
COULAUD Olivier's avatar
COULAUD Olivier committed
37
38
	void RunTest( const bool isBlasKernel){
		//
39
40
		const int DevP = ORDER;
		printf("---------------------------------------- \n ------------------ %d ------------------ \n ---------------------------------------- ",ORDER);
COULAUD Olivier's avatar
COULAUD Olivier committed
41
42
43
44
45
46
47
48
49
50
51
52
53
		//
		// Load particles
		//
		if(sizeof(FReal) == sizeof(float) ) {
			std::cerr << "No input data available for Float "<< std::endl;
			exit(EXIT_FAILURE);
		}
		const std::string parFile( (sizeof(FReal) == sizeof(float))?
				"Test/DirectFloat.bfma":
				"UTest/DirectDouble.bfma");
		//
		std::string filename(SCALFMMDataPath+parFile);
		//
54
		FFmaGenericLoader<FReal> loader(filename);
COULAUD Olivier's avatar
COULAUD Olivier committed
55
56
57
58
59
60
61
62
63
64
65
66
		if(!loader.isOpen()){
			Print("Cannot open particles file.");
			uassert(false);
			return;
		}
		Print("Number of particles:");
		Print(loader.getNumberOfParticles());

		const int NbLevels      = 4;
		const int SizeSubLevels = 2;
		//
		FSize nbParticles = loader.getNumberOfParticles() ;
67
		FmaRWParticle<FReal, 8,8>* const particles = new FmaRWParticle<FReal, 8,8>[nbParticles];
COULAUD Olivier's avatar
COULAUD Olivier committed
68
69
70
71
72

		loader.fillParticle(particles,nbParticles);
		//
		// Create octree
		//
73
		FSphericalCell<FReal>::Init(DevP);
COULAUD Olivier's avatar
COULAUD Olivier committed
74
75
76
		OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
		//   Insert particle in the tree
		//
77
		for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
78
		    tree.insert(particles[idxPart].getPosition() , idxPart, particles[idxPart].getPhysicalValue() );
COULAUD Olivier's avatar
COULAUD Olivier committed
79
80
81
82
83
84
85
86
87
		}



		// Run FMM
		Print("Fmm...");
		//KernelClass kernels(NbLevels,loader.getBoxWidth());
		KernelClass kernels(DevP,NbLevels,loader.getBoxWidth(), loader.getCenterOfBox());
		FmmClass algo(&tree,&kernels);
88
		FTic timer;
COULAUD Olivier's avatar
COULAUD Olivier committed
89
		algo.execute();
90
91
		timer.tac();
		std::cout << "Computation Time : " << ORDER <<" ; "<< timer.elapsed() << std::endl;
COULAUD Olivier's avatar
COULAUD Olivier committed
92
93
94
95
96
97
		//
		FReal energy= 0.0 , energyD = 0.0 ;
		/////////////////////////////////////////////////////////////////////////////////////////////////
		// Compute direct energy
		/////////////////////////////////////////////////////////////////////////////////////////////////

98
		for(FSize idx = 0 ; idx < loader.getNumberOfParticles()  ; ++idx){
99
		    energyD +=  particles[idx].getPotential()*particles[idx].getPhysicalValue() ;
COULAUD Olivier's avatar
COULAUD Olivier committed
100
101
102
103
104
		}
		/////////////////////////////////////////////////////////////////////////////////////////////////
		// Compare
		/////////////////////////////////////////////////////////////////////////////////////////////////
		Print("Compute Diff...");
105
106
		FMath::FAccurater<FReal> potentialDiff;
		FMath::FAccurater<FReal> fx, fy, fz;
COULAUD Olivier's avatar
COULAUD Olivier committed
107
108
109
110
111
112
113
114
		{ // Check that each particle has been summed with all other

			tree.forEachLeaf([&](LeafClass* leaf){
				const FReal*const potentials        = leaf->getTargets()->getPotentials();
				const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
				const FReal*const forcesX            = leaf->getTargets()->getForcesX();
				const FReal*const forcesY            = leaf->getTargets()->getForcesY();
				const FReal*const forcesZ            = leaf->getTargets()->getForcesZ();
115
116
				const FSize nbParticlesInLeaf           = leaf->getTargets()->getNbParticles();
				const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
COULAUD Olivier's avatar
COULAUD Olivier committed
117

118
119
				for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
					const FSize indexPartOrig = indexes[idxPart];
120
121
122
123
					potentialDiff.add(particles[indexPartOrig].getPotential(),potentials[idxPart]);
					fx.add(particles[indexPartOrig].getForces()[0],forcesX[idxPart]);
					fy.add(particles[indexPartOrig].getForces()[1],forcesY[idxPart]);
					fz.add(particles[indexPartOrig].getForces()[2],forcesZ[idxPart]);
COULAUD Olivier's avatar
COULAUD Olivier committed
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
					energy   += potentials[idxPart]*physicalValues[idxPart];
				}
			});
		}

		delete[] particles;

		// Print for information

		Print("Potential diff is = ");
		printf("         Pot L2Norm     %e\n",potentialDiff.getL2Norm());
		printf("         Pot RL2Norm   %e\n",potentialDiff.getRelativeL2Norm());
		printf("         Pot RMSError   %e\n",potentialDiff.getRMSError());
		Print("Fx diff is = ");
		printf("         Fx L2Norm     %e\n",fx.getL2Norm());
		printf("         Fx RL2Norm   %e\n",fx.getRelativeL2Norm());
		printf("         Fx RMSError   %e\n",fx.getRMSError());
		Print("Fy diff is = ");
		printf("        Fy L2Norm     %e\n",fy.getL2Norm());
		printf("        Fy RL2Norm   %e\n",fy.getRelativeL2Norm());
		printf("        Fy RMSError   %e\n",fy.getRMSError());
		Print("Fz diff is = ");
		printf("        Fz L2Norm     %e\n",fz.getL2Norm());
		printf("        Fz RL2Norm   %e\n",fz.getRelativeL2Norm());
		printf("        Fz RMSError   %e\n",fz.getRMSError());
		FReal L2error = (fx.getRelativeL2Norm()*fx.getRelativeL2Norm() + fy.getRelativeL2Norm()*fy.getRelativeL2Norm()  + fz.getRelativeL2Norm() *fz.getRelativeL2Norm()  );
		printf(" Total L2 Force Error= %e\n",FMath::Sqrt(L2error)) ;
		printf("  Energy Error  =   %.12e\n",FMath::Abs(energy-energyD));
		printf("  Energy FMM    =   %.12e\n",FMath::Abs(energy));
		printf("  Energy DIRECT =   %.12e\n",FMath::Abs(energyD));

		// Assert
COULAUD Olivier's avatar
COULAUD Olivier committed
156
		double epsilon = 1.0/FMath::pow2(DevP);
157
158
159
		const FReal MaximumDiffPotential = FReal(epsilon);
		const FReal MaximumDiffForces     = FReal(10*epsilon);
		printf(" Criteria error - Epsilon  %e\n",epsilon);
COULAUD Olivier's avatar
COULAUD Olivier committed
160
161
162
163
164
165
166
167
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
198

		Print("Test1 - Error Relative L2 norm Potential ");
		uassert(potentialDiff.getRelativeL2Norm() < MaximumDiffPotential);    //1
		Print("Test2 - Error RMS L2 norm Potential ");
		uassert(potentialDiff.getRMSError() < MaximumDiffPotential);  //2
		Print("Test3 - Error Relative L2 norm FX ");
		uassert(fx.getRelativeL2Norm()  < MaximumDiffForces);                       //3
		Print("Test4 - Error RMS L2 norm FX ");
		uassert(fx.getRMSError() < MaximumDiffForces);                      //4
		Print("Test5 - Error Relative L2 norm FY ");
		uassert(fy.getRelativeL2Norm()  < MaximumDiffForces);                       //5
		Print("Test6 - Error RMS L2 norm FY ");
		uassert(fy.getRMSError() < MaximumDiffForces);                      //6
		Print("Test7 - Error Relative L2 norm FZ ");
		uassert(fz.getRelativeL2Norm()  < MaximumDiffForces);                      //8
		Print("Test8 - Error RMS L2 norm FZ ");
		uassert(fz.getRMSError() < MaximumDiffForces);                                           //8
		Print("Test9 - Error Relative L2 norm F ");
		uassert(L2error              < MaximumDiffForces);                                            //9   Total Force
		Print("Test10 - Relative error Energy ");
		uassert(FMath::Abs(energy-energyD) /energyD< MaximumDiffPotential);                     //10  Total Energy

	}

	/** 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(){
199
        typedef double FReal;
200
		typedef FSphericalCell<FReal>            CellClass;
201
		typedef FP2PParticleContainerIndexed<FReal>  ContainerClass;
COULAUD Olivier's avatar
COULAUD Olivier committed
202

203
		typedef FSphericalKernel< FReal, CellClass, ContainerClass >          KernelClass;
COULAUD Olivier's avatar
COULAUD Olivier committed
204

205
206
		typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
		typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
COULAUD Olivier's avatar
COULAUD Olivier committed
207
208
209

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

210
		printf("Spherical\n \n");
211
        // RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
212
		//	 OctreeClass, FmmClass, 4>(false);
213
        // RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
214
		//	 OctreeClass, FmmClass, 6>(false);
215
        // RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
216
		//	 OctreeClass, FmmClass, 8>(false);
217
        // RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
218
		//	 OctreeClass, FmmClass, 10>(false);
219
        // RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
220
		//	 OctreeClass, FmmClass, 12>(false);
221
        // RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
222
		//	 OctreeClass, FmmClass, 14>(false);
223
        // RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
224
		//	 OctreeClass, FmmClass, 16>(false);
225
        // RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
226
		//	 OctreeClass, FmmClass, 18>(false);
227
        // RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
228
		//	 OctreeClass, FmmClass, 20>(true);
229
        // RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
230
		//	 OctreeClass, FmmClass, 22>(true);
231
        // RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
232
		//	 OctreeClass, FmmClass, 24>(true);
233
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
234
			 OctreeClass, FmmClass, 50>(true);
235
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
236
			 OctreeClass, FmmClass, 55>(true);
237
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
238
239
240
			 OctreeClass, FmmClass, 60>(true);


COULAUD Olivier's avatar
COULAUD Olivier committed
241
242
243
	}


berenger-bramas's avatar
berenger-bramas committed
244

245
#ifdef SCALFMM_USE_BLAS
COULAUD Olivier's avatar
COULAUD Olivier committed
246
247
	/** Blas */
	void TestSphericalBlas(){
248
        typedef double FReal;
249
		typedef FSphericalCell<FReal>            CellClass;
250
		typedef FP2PParticleContainerIndexed<FReal>  ContainerClass;
berenger-bramas's avatar
berenger-bramas committed
251

252
		typedef FSphericalBlasKernel<FReal, CellClass, ContainerClass >          KernelClass;
berenger-bramas's avatar
berenger-bramas committed
253

254
255
		typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
		typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
berenger-bramas's avatar
berenger-bramas committed
256

COULAUD Olivier's avatar
COULAUD Olivier committed
257
		typedef FFmmAlgorithm<OctreeClass,  CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
berenger-bramas's avatar
berenger-bramas committed
258

259
		printf("Spherical BLAS\n \n");
260
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
261
			 OctreeClass, FmmClass, 4>(true);
262
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
263
			 OctreeClass, FmmClass, 6>(true);
264
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
265
			 OctreeClass, FmmClass, 8>(true);
266
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
267
			 OctreeClass, FmmClass, 10>(true);
268
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
269
			 OctreeClass, FmmClass, 12>(true);
270
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
271
			 OctreeClass, FmmClass, 14>(true);
272
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
273
			 OctreeClass, FmmClass, 16>(true);
274
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
275
			 OctreeClass, FmmClass, 18>(true);
276
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
277
			 OctreeClass, FmmClass, 20>(true);
278
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
279
			 OctreeClass, FmmClass, 22>(true);
280
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
281
			 OctreeClass, FmmClass, 24>(true);
COULAUD Olivier's avatar
COULAUD Olivier committed
282
	}
berenger-bramas's avatar
berenger-bramas committed
283

COULAUD Olivier's avatar
COULAUD Olivier committed
284
285
	/** Block blas */
	void TestSphericalBlockBlas(){
286
        typedef double FReal;
287
		typedef FSphericalCell<FReal>            CellClass;
288
		typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
berenger-bramas's avatar
berenger-bramas committed
289

290
		typedef FSphericalBlockBlasKernel< FReal, CellClass, ContainerClass >          KernelClass;
berenger-bramas's avatar
berenger-bramas committed
291

292
293
		typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
		typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
berenger-bramas's avatar
berenger-bramas committed
294

COULAUD Olivier's avatar
COULAUD Olivier committed
295
		typedef FFmmAlgorithm<OctreeClass,  CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
berenger-bramas's avatar
berenger-bramas committed
296

297
298

		printf("Spherical BLOCK BLAS\n \n");
299
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
300
			 OctreeClass, FmmClass, 4>(true);
301
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
302
			 OctreeClass, FmmClass, 6>(true);
303
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
304
			 OctreeClass, FmmClass, 8>(true);
305
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
306
			 OctreeClass, FmmClass, 10>(true);
307
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
308
			 OctreeClass, FmmClass, 12>(true);
309
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
310
			 OctreeClass, FmmClass, 14>(true);
311
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
312
			 OctreeClass, FmmClass, 16>(true);
313
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
314
			 OctreeClass, FmmClass, 18>(true);
315
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
316
			 OctreeClass, FmmClass, 20>(true);
317
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
318
			 OctreeClass, FmmClass, 22>(true);
319
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
320
			 OctreeClass, FmmClass, 24>(true);
321
        RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
322
			 OctreeClass, FmmClass, 26>(true);
323
//		RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
324
//			 OctreeClass, FmmClass, 28>(true);
325
//		RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
326
//			 OctreeClass, FmmClass, 30>(true);
327
//		RunTest< FReal, CellClass, ContainerClass, KernelClass, LeafClass,
328
//			 OctreeClass, FmmClass, 32>(true);
329

COULAUD Olivier's avatar
COULAUD Olivier committed
330
	}
331
#endif
berenger-bramas's avatar
berenger-bramas committed
332

333

COULAUD Olivier's avatar
COULAUD Olivier committed
334
335
336
	///////////////////////////////////////////////////////////
	// Set the tests!
	///////////////////////////////////////////////////////////
berenger-bramas's avatar
berenger-bramas committed
337

COULAUD Olivier's avatar
COULAUD Olivier committed
338
339
	/** set test */
	void SetTests(){
340
	    //AddTest(&TestSphericalDirect::TestSpherical,"Test Spherical Kernel");
341
#ifdef SCALFMM_USE_BLAS
342
343
	    //AddTest(&TestSphericalDirect::TestSphericalBlas,"Test Spherical Blas Kernel");
	    AddTest(&TestSphericalDirect::TestSphericalBlockBlas,"Test Spherical Block Blas Kernel");
344
#endif
COULAUD Olivier's avatar
COULAUD Olivier committed
345
	}
346
347
348
349
};


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