utestChebyshevDirectPeriodic.cpp 13.6 KB
Newer Older
1
// ===================================================================================
COULAUD Olivier's avatar
COULAUD Olivier committed
2
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas, Matthias Messner
3 4 5 6 7 8 9 10 11 12 13 14
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// 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
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info". 
// "http://www.gnu.org/licenses".
15 16 17 18 19 20
// ===================================================================================

// ==== CMAKE =====
// @FUSE_BLAS
// ================

21
#include "Utils/FGlobal.hpp"
22

23 24
#include "Containers/FOctree.hpp"
#include "Containers/FVector.hpp"
25

26
#include "Files/FFmaGenericLoader.hpp"
27 28
#include "Files/FRandomLoader.hpp"
#include "Files/FTreeIO.hpp"
29

30
#include "Core/FFmmAlgorithmPeriodic.hpp"
31 32 33

#include "FUTester.hpp"

34

35 36 37 38
#include "Kernels/Chebyshev/FChebCell.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Chebyshev/FChebKernel.hpp"
#include "Kernels/Chebyshev/FChebSymKernel.hpp"
39

40 41
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
42

43
/*
COULAUD Olivier's avatar
COULAUD Olivier committed
44 45
  In this test we compare the Chebyshev fmm results and the direct results.
 */
46 47 48 49 50 51


/** the test class
 *
 */
class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
COULAUD Olivier's avatar
COULAUD Olivier committed
52

53 54 55
	///////////////////////////////////////////////////////////
	// The tests!
	///////////////////////////////////////////////////////////
COULAUD Olivier's avatar
COULAUD Olivier committed
56

57
    template <class FReal, class CellClass, class ContainerClass, class KernelClass, class MatrixKernelClass,
COULAUD Olivier's avatar
COULAUD Olivier committed
58
	class LeafClass, class OctreeClass, class FmmClass>
59
	void RunTest(bool readFile=false)	{
60 61 62
		// Warning in make test the exec dir it Build/UTests
		// Load particles

COULAUD Olivier's avatar
COULAUD Olivier committed
63 64
		const int NbLevels        = 4;
		const int SizeSubLevels = 2;
65
		const int PeriodicDeep  = 3;
COULAUD Olivier's avatar
COULAUD Olivier committed
66

67 68 69
		int NbParticles     = 150;
		FReal BoxWidth;
		FPoint<FReal> CenterOfBox ;
70
		FRandomLoader<FReal> loader(NbParticles);
COULAUD Olivier's avatar
COULAUD Olivier committed
71

72 73 74 75
		BoxWidth      = loader.getBoxWidth();
		CenterOfBox = loader.getCenterOfBox();
		NbParticles   = loader.getNumberOfParticles();

COULAUD Olivier's avatar
COULAUD Olivier committed
76 77 78 79
		Print("Number of particles:");
		Print(loader.getNumberOfParticles());

		// Create octree
80
		OctreeClass tree(NbLevels, SizeSubLevels,BoxWidth,CenterOfBox);
COULAUD Olivier's avatar
COULAUD Olivier committed
81

82 83
        // interaction kernel evaluator
        const MatrixKernelClass MatrixKernel;
84

COULAUD Olivier's avatar
COULAUD Olivier committed
85
		struct TestParticle{
86
            FPoint<FReal> position;
COULAUD Olivier's avatar
COULAUD Olivier committed
87 88
			FReal physicalValue;
			FReal potential;
89
			FReal forces[3];
COULAUD Olivier's avatar
COULAUD Olivier committed
90 91 92 93
		};
		FReal coeff = -1.0, value = 0.10, sum = 0.0;
		TestParticle* const particles = new TestParticle[loader.getNumberOfParticles()];
		for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
94
            FPoint<FReal> position;
COULAUD Olivier's avatar
COULAUD Olivier committed
95 96 97 98 99 100 101 102 103 104 105 106 107 108
			loader.fillParticle(&position);
			// put in tree
			value *= coeff ;
			sum += value ;
			// put in tree
			tree.insert(position, idxPart, value);
			// get copy
			particles[idxPart].position         = position;
			particles[idxPart].physicalValue = value;
			particles[idxPart].potential        = 0.0;
			particles[idxPart].forces[0]        = 0.0;
			particles[idxPart].forces[1]        = 0.0;
			particles[idxPart].forces[2]        = 0.0;
		}
109 110 111 112
        if (FMath::Abs(sum)> 0.00001){
        		std::cerr << "Sum of charges is not equal zero!!! (sum=<<"<<sum<<" )"<<std::endl;
        		exit(-1);
        }
COULAUD Olivier's avatar
COULAUD Olivier committed
113 114 115
		/////////////////////////////////////////////////////////////////////////////////////////////////
		// Run FMM computation
		/////////////////////////////////////////////////////////////////////////////////////////////////
COULAUD Olivier's avatar
COULAUD Olivier committed
116 117
		Print("Fmm...");
		FmmClass algo(&tree,PeriodicDeep );
118 119
		KernelClass *kernels = new KernelClass(algo.extendedTreeHeight(), algo.extendedBoxWidth(), algo.extendedBoxCenter(),&MatrixKernel);
		algo.setKernel(kernels);
COULAUD Olivier's avatar
COULAUD Olivier committed
120
		algo.execute();
121
		delete kernels ;
COULAUD Olivier's avatar
COULAUD Olivier committed
122
		/////////////////////////////////////////////////////////////////////////////////////////////////
COULAUD Olivier's avatar
COULAUD Olivier committed
123
		// Run direct computation
COULAUD Olivier's avatar
COULAUD Olivier committed
124
		/////////////////////////////////////////////////////////////////////////////////////////////////
COULAUD Olivier's avatar
COULAUD Olivier committed
125 126 127 128
		Print("Direct...");

		FTreeCoordinate min, max;
		algo.repetitionsIntervals(&min, &max);
COULAUD Olivier's avatar
COULAUD Olivier committed
129
		FReal energy= 0.0 , energyD = 0.0 ;
COULAUD Olivier's avatar
COULAUD Olivier committed
130 131 132 133 134 135 136 137 138 139

		for(int idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
			for(int idxOther = idxTarget + 1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
				FP2P::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
						particles[idxTarget].position.getZ(),particles[idxTarget].physicalValue,
						&particles[idxTarget].forces[0],&particles[idxTarget].forces[1],
						&particles[idxTarget].forces[2],&particles[idxTarget].potential,
						particles[idxOther].position.getX(), particles[idxOther].position.getY(),
						particles[idxOther].position.getZ(),particles[idxOther].physicalValue,
						&particles[idxOther].forces[0],&particles[idxOther].forces[1],
140
                              &particles[idxOther].forces[2],&particles[idxOther].potential,&MatrixKernel);
COULAUD Olivier's avatar
COULAUD Olivier committed
141 142 143 144 145 146 147 148

			}
			for(int idxX = min.getX() ; idxX <= max.getX() ; ++idxX){
				for(int idxY = min.getY() ; idxY <= max.getY() ; ++idxY){
					for(int idxZ = min.getZ() ; idxZ <= max.getZ() ; ++idxZ){
						if(idxX ==0 && idxY == 0 && idxZ == 0) continue;
						// next lines for test

149
                        const FPoint<FReal> offset(loader.getBoxWidth() * FReal(idxX),
COULAUD Olivier's avatar
COULAUD Olivier committed
150 151 152 153 154 155 156 157 158 159 160 161 162
								loader.getBoxWidth() * FReal(idxY),
								loader.getBoxWidth() * FReal(idxZ));

						for(int idxSource = 0 ; idxSource < NbParticles ; ++idxSource){
							TestParticle source = particles[idxSource];
							source.position += offset;

							FP2P::NonMutualParticles(
									source.position.getX(), source.position.getY(),
									source.position.getZ(),source.physicalValue,
									particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
									particles[idxTarget].position.getZ(),particles[idxTarget].physicalValue,
									&particles[idxTarget].forces[0],&particles[idxTarget].forces[1],
163
									&particles[idxTarget].forces[2],&particles[idxTarget].potential,&MatrixKernel);
COULAUD Olivier's avatar
COULAUD Olivier committed
164 165 166 167 168
						}
					}
				}
			}
		}
COULAUD Olivier's avatar
COULAUD Olivier committed
169 170 171 172
		for(int idx = 0 ; idx <  loader.getNumberOfParticles()  ; ++idx){
			energyD +=  particles[idx].potential*particles[idx].physicalValue ;
		}
		/////////////////////////////////////////////////////////////////////////////////////////////////
COULAUD Olivier's avatar
COULAUD Olivier committed
173
		// Compare
COULAUD Olivier's avatar
COULAUD Olivier committed
174 175
		/////////////////////////////////////////////////////////////////////////////////////////////////

COULAUD Olivier's avatar
COULAUD Olivier committed
176
		Print("Compute Diff...");
177 178
		FMath::FAccurater<FReal> potentialDiff;
		FMath::FAccurater<FReal> fx, fy, fz;
COULAUD Olivier's avatar
COULAUD Olivier committed
179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196

		{ // 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();
				const int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
				const FVector<int>& indexes = leaf->getTargets()->getIndexes();

				for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
					const int indexPartOrig = indexes[idxPart];
					potentialDiff.add(particles[indexPartOrig].potential,potentials[idxPart]);
					fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]);
					fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
					fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);
COULAUD Olivier's avatar
COULAUD Olivier committed
197
					energy+=potentials[idxPart]*physicalValues[idxPart];
COULAUD Olivier's avatar
COULAUD Olivier committed
198 199 200 201 202
				}
			});
		}
		// Print for information

COULAUD Olivier's avatar
COULAUD Olivier committed
203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219
        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()  );
COULAUD Olivier's avatar
COULAUD Olivier committed
220
		printf(" Total L2 Force Error= %e\n",FMath::Sqrt(L2error)) ;
COULAUD Olivier's avatar
COULAUD Olivier committed
221 222
		printf("  Energy Error  =   %.12e\n",FMath::Abs(energy-energyD));
		printf("  Energy FMM    =   %.12e\n",FMath::Abs(energy));
COULAUD Olivier's avatar
COULAUD Olivier committed
223
		printf("  Energy DIRECT =   %.12e\n",FMath::Abs(energyD));
COULAUD Olivier's avatar
COULAUD Olivier committed
224

COULAUD Olivier's avatar
COULAUD Olivier committed
225
		// Assert
COULAUD Olivier's avatar
COULAUD Olivier committed
226 227 228 229 230
        const FReal MaximumDiffPotential = FReal(9e-3);
        const FReal MaximumDiffForces     = FReal(9e-2);


        uassert(potentialDiff.getL2Norm() < MaximumDiffPotential);    //1
231
   //     uassert(potentialDiff.getRMSError() < MaximumDiffPotential);  //2
COULAUD Olivier's avatar
COULAUD Olivier committed
232
        uassert(fx.getL2Norm()  < MaximumDiffForces);                       //3
233
  //      uassert(fx.getRMSError() < MaximumDiffForces);                      //4
COULAUD Olivier's avatar
COULAUD Olivier committed
234
        uassert(fy.getL2Norm()  < MaximumDiffForces);                       //5
235
 //       uassert(fy.getRMSError() < MaximumDiffForces);                      //6
COULAUD Olivier's avatar
COULAUD Olivier committed
236
        uassert(fz.getL2Norm()  < MaximumDiffForces);                      //8
237
  //      uassert(fz.getRMSError() < MaximumDiffForces);                                           //8
COULAUD Olivier's avatar
COULAUD Olivier committed
238 239
        uassert(L2error              < MaximumDiffForces);                                            //9   Total Force
        uassert(FMath::Abs(energy-energyD) < 10*MaximumDiffPotential);                     //10  Total Energy
240

241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261


	    //
	    // ----------------------------------------------------------------
	    //  Save  computation in binary format
	    //
	    // 2 is for PeriodicDeep
//		std::string filenameOut("../Data/UTest/DirectPer2Double.bfma");
//	    std::cout << "Generate " << filenameOut <<"  for output file" << std::endl;
//	    //
//	    std::cout << " nbParticles: " << NbParticles <<"  " << sizeof(NbParticles) <<std::endl;
//	    std::cout << " Box size: " << loader.getBoxWidth() << "  " << sizeof(loader.getBoxWidth())<<std::endl;
//	    //
//	    FFmaGenericWriter<FReal> writer(filenameOut) ;
//
//	    writer.writeHeader( loader.getCenterOfBox(),loader.getBoxWidth(),  NbParticles,8,8) ;
//	//    writer.writeArrayOfParticlesp(articles, NbParticles);
//	     writer.writeArrayOfReal(particles[0].position.getDataValue(), 8, NbParticles);

        delete[] particles;

COULAUD Olivier's avatar
COULAUD Olivier committed
262
	}
263 264 265 266 267

	/** If memstas is running print the memory used */
	void PostTest() {
		if( FMemStats::controler.isUsed() ){
			std::cout << "Memory used at the end " << FMemStats::controler.getCurrentAllocated()
COULAUD Olivier's avatar
COULAUD Olivier committed
268
										<< " Bytes (" << FMemStats::controler.getCurrentAllocatedMB() << "MB)\n";
269
			std::cout << "Max memory used " << FMemStats::controler.getMaxAllocated()
COULAUD Olivier's avatar
COULAUD Olivier committed
270
										<< " Bytes (" << FMemStats::controler.getMaxAllocatedMB() << "MB)\n";
271
			std::cout << "Total memory used " << FMemStats::controler.getTotalAllocated()
COULAUD Olivier's avatar
COULAUD Olivier committed
272
										<< " Bytes (" << FMemStats::controler.getTotalAllocatedMB() << "MB)\n";
273 274 275 276 277 278 279 280 281 282 283
		}
	}


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


	/** TestChebKernel */
	void TestChebKernel(){
COULAUD Olivier's avatar
COULAUD Olivier committed
284
		const unsigned int ORDER = 6;
285
        typedef double FReal;
286 287 288 289 290 291
		typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
		typedef FSimpleLeaf<FReal, ContainerClass> LeafClass;
		typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
		typedef FChebCell<FReal,ORDER> CellClass;
		typedef FOctree<FReal, CellClass,ContainerClass,LeafClass> OctreeClass;
		typedef FChebKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
292
		typedef FFmmAlgorithmPeriodic<FReal,OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
293
		// run test
COULAUD Olivier's avatar
COULAUD Olivier committed
294
		std::cout <<" TEST 1  "<<std::endl;
295
        RunTest<FReal,CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass>();
296 297 298 299
	}

	/** TestChebSymKernel */
	void TestChebSymKernel(){
COULAUD Olivier's avatar
COULAUD Olivier committed
300
		const unsigned int ORDER = 7;
301
        typedef double FReal;
302 303 304 305 306 307
		typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
		typedef FSimpleLeaf<FReal, ContainerClass> LeafClass;
		typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
		typedef FChebCell<FReal,ORDER> CellClass;
		typedef FOctree<FReal, CellClass,ContainerClass,LeafClass> OctreeClass;
		typedef FChebSymKernel<FReal,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
308
		typedef FFmmAlgorithmPeriodic<FReal,OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
309
		// run test
COULAUD Olivier's avatar
COULAUD Olivier committed
310
		std::cout <<std::endl<<" TEST 2 "<<std::endl;
311
        RunTest<FReal,CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass>();
312 313 314 315 316 317 318 319 320 321
	}



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

	/** set test */
	void SetTests(){
COULAUD Olivier's avatar
COULAUD Olivier committed
322
		AddTest(&TestChebyshevDirect::TestChebKernel,"Test Chebyshev Kernel with one big SVD") ;
COULAUD Olivier's avatar
COULAUD Olivier committed
323
		AddTest(&TestChebyshevDirect::TestChebSymKernel,"Test Chebyshev Kernel with 16 small SVDs and symmetries");
324 325 326 327 328 329 330 331 332 333
	}
};


// You must do this
TestClass(TestChebyshevDirect)