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 16 17
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/Spherical/FSphericalKernel.hpp"
#include "Kernels/Spherical/FSphericalRotationKernel.hpp"
#include "Kernels/Spherical/FSphericalBlasKernel.hpp"
#include "Kernels/Spherical/FSphericalBlockBlasKernel.hpp"
18

19
#include "Files/FFmaGenericLoader.hpp"
20

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

#include "FUTester.hpp"

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

berenger-bramas's avatar
berenger-bramas committed
29
/** the test class
30 31
 *
 */
32
class TestSphericalDirect : public FUTester<TestSphericalDirect> {
33
	/** The test method to factorize all the test based on different kernels */
34
    template < class FReal, class CellClass, class ContainerClass, class KernelClass, class LeafClass,
35
		   class OctreeClass, class FmmClass, int ORDER>
36 37
	void RunTest( const bool isBlasKernel){
		//
38 39
		const int DevP = ORDER;
		printf("---------------------------------------- \n ------------------ %d ------------------ \n ---------------------------------------- ",ORDER);
40 41 42 43 44 45 46 47 48 49 50 51 52
		//
		// 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);
		//
53
		FFmaGenericLoader<FReal> loader(filename);
54 55 56 57 58 59 60 61 62 63 64 65
		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() ;
66
		FmaRWParticle<FReal, 8,8>* const particles = new FmaRWParticle<FReal, 8,8>[nbParticles];
67 68 69 70 71

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



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

97
		for(FSize idx = 0 ; idx < loader.getNumberOfParticles()  ; ++idx){
98
		    energyD +=  particles[idx].getPotential()*particles[idx].getPhysicalValue() ;
99 100 101 102 103
		}
		/////////////////////////////////////////////////////////////////////////////////////////////////
		// Compare
		/////////////////////////////////////////////////////////////////////////////////////////////////
		Print("Compute Diff...");
104 105
		FMath::FAccurater<FReal> potentialDiff;
		FMath::FAccurater<FReal> fx, fy, fz;
106 107 108 109 110 111 112 113
		{ // 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();
114 115
				const FSize nbParticlesInLeaf           = leaf->getTargets()->getNbParticles();
				const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
116

117 118
				for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
					const FSize indexPartOrig = indexes[idxPart];
119 120 121 122
					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]);
123 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
					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
155
		double epsilon = 1.0/FMath::pow2(DevP);
156 157 158
		const FReal MaximumDiffPotential = FReal(epsilon);
		const FReal MaximumDiffForces     = FReal(10*epsilon);
		printf(" Criteria error - Epsilon  %e\n",epsilon);
159 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

		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(){
198
        typedef double FReal;
199
		typedef FSphericalCell<FReal>            CellClass;
200
		typedef FP2PParticleContainerIndexed<FReal>  ContainerClass;
201

202
		typedef FSphericalKernel< FReal, CellClass, ContainerClass >          KernelClass;
203

204 205
		typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
		typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
206 207 208

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

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


240 241 242
	}


berenger-bramas's avatar
berenger-bramas committed
243

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

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

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

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

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

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

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

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

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

296 297

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

329
	}
330
#endif
berenger-bramas's avatar
berenger-bramas committed
331

332

333 334 335
	///////////////////////////////////////////////////////////
	// Set the tests!
	///////////////////////////////////////////////////////////
berenger-bramas's avatar
berenger-bramas committed
336

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


// You must do this
349
TestClass(TestSphericalDirect)