utestSphericalBlasAlgorithm.cpp 10.8 KB
Newer Older
1
// ===================================================================================
2
3
// Copyright ScalFmm 2016 INRIA
//
4
5
// This software is a computer program whose purpose is to compute the FMM.
//
6
// This software is governed by Mozilla Public License Version 2.0 (MPL 2.0) and
7
8
9
10
11
// 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
12
13
// Mozilla Public License Version 2.0 (MPL 2.0) for more details.
// https://www.mozilla.org/en-US/MPL/2.0/
14
15
// ===================================================================================

BRAMAS Berenger's avatar
BRAMAS Berenger committed
16
#include "Utils/FGlobal.hpp"
17

BRAMAS Berenger's avatar
BRAMAS Berenger committed
18
19
#include "Containers/FOctree.hpp"
#include "Containers/FVector.hpp"
20

BRAMAS Berenger's avatar
BRAMAS Berenger committed
21
22
#include "Kernels/Spherical/FSphericalCell.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
23

BRAMAS Berenger's avatar
BRAMAS Berenger committed
24
25
26
27
28
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/Spherical/FSphericalKernel.hpp"
#include "Kernels/Spherical/FSphericalRotationKernel.hpp"
#include "Kernels/Spherical/FSphericalBlasKernel.hpp"
#include "Kernels/Spherical/FSphericalBlockBlasKernel.hpp"
29

BRAMAS Berenger's avatar
BRAMAS Berenger committed
30
#include "Files/FFmaGenericLoader.hpp"
31

BRAMAS Berenger's avatar
BRAMAS Berenger committed
32
#include "Core/FFmmAlgorithm.hpp"
33
34
35
36
37
38
39
40
41
42
43
44

#include "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 */
45
    template < class FReal, class CellClass, class ContainerClass, class KernelClass, class LeafClass,
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
	class OctreeClass, class FmmClass>
	void RunTest( const bool isBlasKernel){
		//
		const int DevP = 9;
		//
		// 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);
		//
63
		FFmaGenericLoader<FReal> loader(filename);
64
65
66
67
68
69
70
71
72
73
74
75
		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() ;
76
		FmaRWParticle<FReal, 8,8>* const particles = new FmaRWParticle<FReal, 8,8>[nbParticles];
77
78
79
80
81

		loader.fillParticle(particles,nbParticles);
		//
		// Create octree
		//
82
		FSphericalCell<FReal>::Init(DevP);
83
84
85
		OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
		//   Insert particle in the tree
		//
86
		for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
		    tree.insert(particles[idxPart].getPosition() , idxPart, particles[idxPart].getPhysicalValue() );
		}



		// Run FMM
		Print("Fmm...");
		//KernelClass kernels(NbLevels,loader.getBoxWidth());
		KernelClass kernels(DevP,NbLevels,loader.getBoxWidth(), loader.getCenterOfBox());
		FmmClass algo(&tree,&kernels);
		algo.execute();
		//
		FReal energy= 0.0 , energyD = 0.0 ;
		/////////////////////////////////////////////////////////////////////////////////////////////////
		// Compute direct energy
		/////////////////////////////////////////////////////////////////////////////////////////////////

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

124
125
				for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
					const FSize indexPartOrig = indexes[idxPart];
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
156
157
158
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
198
199
200
201
202
					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]);
					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
		const FReal MaximumDiffPotential = FReal(9e-3);
		const FReal MaximumDiffForces     = FReal(9e-2);

		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(){
203
        typedef double FReal;
204
		typedef FSphericalCell<FReal>            CellClass;
205
		typedef FP2PParticleContainerIndexed<FReal>  ContainerClass;
206

207
		typedef FSphericalKernel< FReal, CellClass, ContainerClass >          KernelClass;
208

209
210
		typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
		typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
211
212
213

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

214
        RunTest<FReal, CellClass, ContainerClass, KernelClass, LeafClass,
215
216
217
218
219
		OctreeClass, FmmClass>(false);
	}



220
#ifdef SCALFMM_USE_BLAS
221
222
	/** Blas */
	void TestSphericalBlas(){
223
        typedef double FReal;
224
		typedef FSphericalCell<FReal>            CellClass;
225
		typedef FP2PParticleContainerIndexed<FReal>  ContainerClass;
226

227
		typedef FSphericalBlasKernel<FReal, CellClass, ContainerClass >          KernelClass;
228

229
230
		typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
		typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
231
232
233

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

234
        RunTest<FReal, CellClass, ContainerClass, KernelClass, LeafClass,
235
236
237
238
239
		OctreeClass, FmmClass>(true);
	}

	/** Block blas */
	void TestSphericalBlockBlas(){
240
        typedef double FReal;
241
		typedef FSphericalCell<FReal>            CellClass;
242
		typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
243

244
		typedef FSphericalBlockBlasKernel< FReal, CellClass, ContainerClass >          KernelClass;
245

246
247
		typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
		typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
248
249
250

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

251
        RunTest<FReal, CellClass, ContainerClass, KernelClass, LeafClass,
252
253
254
255
256
257
258
259
260
261
262
		OctreeClass, FmmClass>(true);
	}
#endif

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

	/** set test */
	void SetTests(){
		AddTest(&TestSphericalDirect::TestSpherical,"Test Spherical Kernel");
263
#ifdef SCALFMM_USE_BLAS
264
265
266
267
268
269
270
271
272
		AddTest(&TestSphericalDirect::TestSphericalBlas,"Test Spherical Blas Kernel");
		AddTest(&TestSphericalDirect::TestSphericalBlockBlas,"Test Spherical Block Blas Kernel");
#endif
	}
};


// You must do this
TestClass(TestSphericalDirect)