utestSpherical.cpp 10.8 KB
Newer Older
COULAUD Olivier's avatar
COULAUD Olivier committed
1
// ===================================================================================
2
3
// Copyright ScalFmm 2016 INRIA
//
COULAUD Olivier's avatar
COULAUD Olivier committed
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
COULAUD Olivier's avatar
COULAUD Olivier committed
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/
COULAUD Olivier's avatar
COULAUD Olivier committed
14
15
// ===================================================================================

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

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

BRAMAS Berenger's avatar
BRAMAS Berenger committed
21
22
#include "Kernels/Spherical/FSphericalCell.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
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"
COULAUD Olivier's avatar
COULAUD Olivier committed
29

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

BRAMAS Berenger's avatar
BRAMAS Berenger committed
32
#include "Core/FFmmAlgorithm.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
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,
COULAUD Olivier's avatar
COULAUD Olivier committed
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);
COULAUD Olivier's avatar
COULAUD Olivier committed
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];
COULAUD Olivier's avatar
COULAUD Olivier committed
77
78
79
80
81

		loader.fillParticle(particles,nbParticles);
		//
		// Create octree
		//
82
		FSphericalCell<FReal>::Init(DevP);
COULAUD Olivier's avatar
COULAUD Olivier committed
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
		    tree.insert(particles[idxPart].getPosition() , idxPart, particles[idxPart].getPhysicalValue() );
COULAUD Olivier's avatar
COULAUD Olivier committed
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
		}



		// 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
		    energyD +=  particles[idx].getPotential()*particles[idx].getPhysicalValue() ;
COULAUD Olivier's avatar
COULAUD Olivier committed
106
107
108
109
110
		}
		/////////////////////////////////////////////////////////////////////////////////////////////////
		// Compare
		/////////////////////////////////////////////////////////////////////////////////////////////////
		Print("Compute Diff...");
111
112
		FMath::FAccurater<FReal> potentialDiff;
		FMath::FAccurater<FReal> fx, fy, fz;
COULAUD Olivier's avatar
COULAUD Olivier committed
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();
COULAUD Olivier's avatar
COULAUD Olivier committed
123

124
125
				for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
					const FSize indexPartOrig = indexes[idxPart];
126
127
128
129
					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
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
					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;
COULAUD Olivier's avatar
COULAUD Olivier committed
206

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

209
210
		typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
		typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
COULAUD Olivier's avatar
COULAUD Olivier committed
211
212
213

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

214
        RunTest<FReal, CellClass, ContainerClass, KernelClass, LeafClass,
COULAUD Olivier's avatar
COULAUD Olivier committed
215
216
217
218
219
		OctreeClass, FmmClass>(false);
	}



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

227
		typedef FSphericalBlasKernel<FReal, CellClass, ContainerClass >          KernelClass;
COULAUD Olivier's avatar
COULAUD Olivier committed
228

229
230
		typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
		typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
COULAUD Olivier's avatar
COULAUD Olivier committed
231
232
233

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

234
        RunTest<FReal, CellClass, ContainerClass, KernelClass, LeafClass,
COULAUD Olivier's avatar
COULAUD Olivier committed
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;
COULAUD Olivier's avatar
COULAUD Olivier committed
243

244
		typedef FSphericalBlockBlasKernel< FReal, CellClass, ContainerClass >          KernelClass;
COULAUD Olivier's avatar
COULAUD Olivier committed
245

246
247
		typedef FSimpleLeaf<FReal, ContainerClass >                     LeafClass;
		typedef FOctree<FReal, CellClass, ContainerClass , LeafClass >  OctreeClass;
COULAUD Olivier's avatar
COULAUD Olivier committed
248
249
250

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

251
        RunTest<FReal, CellClass, ContainerClass, KernelClass, LeafClass,
COULAUD Olivier's avatar
COULAUD Olivier committed
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
COULAUD Olivier's avatar
COULAUD Olivier committed
264
265
266
267
268
269
270
271
272
273
274
275
		AddTest(&TestSphericalDirect::TestSphericalBlas,"Test Spherical Blas Kernel");
		AddTest(&TestSphericalDirect::TestSphericalBlockBlas,"Test Spherical Block Blas Kernel");
#endif
	}
};


// You must do this
TestClass(TestSphericalDirect)