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"

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

berenger-bramas's avatar
berenger-bramas committed
30
/** the test class
31 32
 *
 */
33
class TestSphericalDirect : public FUTester<TestSphericalDirect> {
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>
37 38
	void RunTest( const bool isBlasKernel){
		//
39 40
		const int DevP = ORDER;
		printf("---------------------------------------- \n ------------------ %d ------------------ \n ---------------------------------------- ",ORDER);
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);
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];
68 69 70 71 72

		loader.fillParticle(particles,nbParticles);
		//
		// Create octree
		//
73
		FSphericalCell<FReal>::Init(DevP);
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() );
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;
89
		algo.execute();
90 91
		timer.tac();
		std::cout << "Computation Time : " << ORDER <<" ; "<< timer.elapsed() << std::endl;
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() ;
100 101 102 103 104
		}
		/////////////////////////////////////////////////////////////////////////////////////////////////
		// Compare
		/////////////////////////////////////////////////////////////////////////////////////////////////
		Print("Compute Diff...");
105 106
		FMath::FAccurater<FReal> potentialDiff;
		FMath::FAccurater<FReal> fx, fy, fz;
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();
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]);
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
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);
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;
202

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

205 206
		typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
		typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
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);


241 242 243
	}


berenger-bramas's avatar
berenger-bramas committed
244

245
#ifdef SCALFMM_USE_BLAS
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

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);
282
	}
berenger-bramas's avatar
berenger-bramas committed
283

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

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

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

333

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

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
345
	}
346 347 348 349
};


// You must do this
350
TestClass(TestSphericalDirect)