RotationFMM.cpp 7.85 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas
// 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".
// ===================================================================================

// ==== CMAKE =====
//
// ================

#include <iostream>
22
#include <stdexcept>
COULAUD Olivier's avatar
COULAUD Olivier committed
23
#include <string>
24
#include <cstdlib>
25
#include <cstdio>
26

COULAUD Olivier's avatar
COULAUD Olivier committed
27
#include "ScalFmmConfig.h"
COULAUD Olivier's avatar
COULAUD Olivier committed
28 29
#include "Utils/FParameters.hpp"
#include "Files/FFmaGenericLoader.hpp"
30

COULAUD Olivier's avatar
COULAUD Olivier committed
31 32
#include "Kernels/Rotation/FRotationKernel.hpp"
#include "Kernels/Rotation/FRotationCell.hpp"
33

COULAUD Olivier's avatar
COULAUD Olivier committed
34 35
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
36 37


COULAUD Olivier's avatar
COULAUD Olivier committed
38
#include "Containers/FOctree.hpp"
39

COULAUD Olivier's avatar
COULAUD Olivier committed
40
#ifdef _OPENMP
COULAUD Olivier's avatar
COULAUD Olivier committed
41
#include "Core/FFmmAlgorithmThread.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
42 43 44
#else
#include "Core/FFmmAlgorithm.hpp"
#endif
45

46
#include "Utils/FParameterNames.hpp"
47

48 49 50 51 52 53 54 55 56 57 58 59 60 61
/// \file  RotationFMM.cpp
//!
//! \brief This program runs the FMM Algorithm with harmonic spherical approximation of 1/r kernel
//!  \authors B. Bramas, O. Coulaud
//!
//!  This code is a short example to use the  rotation harmonic spherical approximation  for the 1/r kernel
//!
//!
//!  <b> General arguments:</b>
//!     \param   -help(-h)      to see the parameters available in this driver
//!     \param   -depth          The depth of the octree
//!     \param   -subdepth     Specifies the size of the sub octree
//!     \param   -t                   The number of threads
//!
62
//!     \param   -f name          Name of the particles file with extension (.fma or .bfma). The data in  file have to be in our FMA format
63 64 65 66 67 68 69
//!
//


// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
70 71 72 73 74 75
    FHelpDescribeAndExit(argc, argv,
                         "Driver for HArmonic Spherical + Rotation  --  kernel  (1/r kernel).",
                         FParameterDefinitions::InputFile, FParameterDefinitions::OctreeHeight,
                         FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::InputFile,
                         FParameterDefinitions::NbThreads);

76
    const std::string defaultFile(/*SCALFMMDataPath+*/"../Data/test20k.fma");
77 78 79 80 81
    const std::string filename(FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, defaultFile.c_str()));
    const unsigned int TreeHeight       = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 5);
    const unsigned int SubTreeHeight  = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeSubHeight.options, 2);
    const unsigned int NbThreads        = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, 1);

82
#ifdef _OPENMP
COULAUD Olivier's avatar
COULAUD Olivier committed
83
	omp_set_num_threads(NbThreads) ;
84 85 86 87 88 89 90 91 92 93 94 95 96 97
	std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl;
#else
	std::cout << "\n>> Sequential version.\n" << std::endl;
#endif
	//
	std::cout <<	 "Parameters  "<< std::endl
			<<     "      Octree Depth      "<< TreeHeight <<std::endl
			<<	  "      SubOctree depth "<< SubTreeHeight <<std::endl
			<<     "      Input file  name: " <<filename <<std::endl
			<<     "      Thread number:  " << NbThreads <<std::endl
			<<std::endl;
	//
	// init timer
	FTic time;
98
	//
99
	// open particle file
100
	//
101
	////////////////////////////////////////////////////////////////////
COULAUD Olivier's avatar
COULAUD Olivier committed
102
	//
103 104
    typedef double FReal;
    FFmaGenericLoader<FReal> loader(filename);
COULAUD Olivier's avatar
COULAUD Olivier committed
105 106 107
	//
	////////////////////////////////////////////////////////////////////
	//
108
    // begin spherical kernel
109
	// accuracy
COULAUD Olivier's avatar
COULAUD Olivier committed
110
	const unsigned int P = 22;
111
	// typedefs
112 113 114 115
    typedef FP2PParticleContainerIndexed<FReal>                     ContainerClass;
    typedef FSimpleLeaf< FReal, ContainerClass >                       LeafClass;
    typedef FRotationCell<FReal, P>                                             CellClass;
    typedef FOctree<FReal, CellClass,ContainerClass,LeafClass>  OctreeClass;
116
	//
117
    typedef FRotationKernel< FReal, CellClass, ContainerClass , P>   KernelClass;
118
	//
COULAUD Olivier's avatar
COULAUD Olivier committed
119
#ifdef _OPENMP
120
	typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
COULAUD Olivier's avatar
COULAUD Olivier committed
121 122 123
#else
	typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
#endif
124
	// init oct-tree
125
	OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
126 127 128


	{ // -----------------------------------------------------
129
		std::cout << "Creating & Inserting " << loader.getNumberOfParticles()
COULAUD Olivier's avatar
COULAUD Olivier committed
130
                																<< " particles ..." << std::endl;
131 132 133
		std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
		time.tic();
		//
134
        FPoint<FReal> position;
135 136
		FReal physicalValue = 0.0;
		//
137
        for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
138
			// Read particles from file
139
			loader.fillParticle(&position,&physicalValue);
140 141 142 143 144 145 146

			// put particles in octree
			tree.insert(position, idxPart, physicalValue);
		}

		time.tac();
		std::cout << "Done  " << "(@Creating and Inserting Particles = "
COULAUD Olivier's avatar
COULAUD Olivier committed
147
				<< time.elapsed() << "  s) ." << std::endl;
148 149 150 151 152 153 154 155 156
	} // -----------------------------------------------------

	{ // -----------------------------------------------------
		std::cout << "\nRotation harmonic Spherical FMM (P="<< P << ") ... " << std::endl;

		time.tic();
		//
		// Here we use a pointer due to the limited size of the stack
		//
157
		KernelClass *kernels = new KernelClass(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
158
		//
COULAUD Olivier's avatar
COULAUD Olivier committed
159
		FmmClass algorithm(&tree, kernels) ;
160 161 162 163
		//
		algorithm.execute();   // Here the call of the FMM algorithm
		//
		time.tac();
COULAUD Olivier's avatar
COULAUD Olivier committed
164
		std::cout << "Done  " << "(@Algorithm = " << time.elapsed() << " s) ." << std::endl;
165 166 167 168 169 170 171
	}
	// -----------------------------------------------------
	//
	// Some output
	//
	//
	{ // -----------------------------------------------------
172
        FSize N1=0, N2= loader.getNumberOfParticles()/2, N3= loader.getNumberOfParticles() -1; ;
173 174 175 176 177 178 179 180 181 182 183 184 185
		FReal energy =0.0 ;
		//
		//   Loop over all leaves
		//
		std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl;
		std::cout << std::scientific;
		std::cout.precision(10) ;

		tree.forEachLeaf([&](LeafClass* leaf){
			const FReal*const potentials = leaf->getTargets()->getPotentials();
			const FReal*const forcesX = leaf->getTargets()->getForcesX();
			const FReal*const forcesY = leaf->getTargets()->getForcesY();
			const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
186
            const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
187 188
			const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();

189
            const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
190

191 192
            for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                const FSize indexPartOrig = indexes[idxPart];
193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208
				if ((indexPartOrig == N1) || (indexPartOrig == N2) || (indexPartOrig == N3)  ) {
					std::cout << "Index "<< indexPartOrig <<"  potential  " << potentials[idxPart]
					                                                                      << "   Forces: " << forcesX[idxPart] << " " << forcesY[idxPart] << " "<< forcesZ[idxPart] <<std::endl;
				}
				energy += potentials[idxPart]*physicalValues[idxPart] ;
			}
		});
		std::cout <<std::endl<<"Energy: "<< energy<<std::endl;
		std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl<<std::endl;

	}
	// -----------------------------------------------------


	return 0;
}