RotationFMM.cpp 8.05 KB
Newer Older
1
// ===================================================================================
2 3 4 5
// Copyright ScalFmm 2016 INRIA, Olivier Coulaud, Bérenger Bramas,
// Matthias Messner olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the
// FMM.
6
//
7
// This software is governed by the CeCILL-C and LGPL licenses and
8
// abiding by the rules of distribution of free software.
9 10 11
// An extension to the license is given to allow static linking of scalfmm
// inside a proprietary application (no matter its license).
// See the main license file for more details.
12 13 14 15
//
// 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
16 17 18
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
19 20 21 22 23 24 25
// ===================================================================================

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

#include <iostream>
26
#include <stdexcept>
COULAUD Olivier's avatar
COULAUD Olivier committed
27
#include <string>
28
#include <cstdlib>
29
#include <cstdio>
30

COULAUD Olivier's avatar
COULAUD Olivier committed
31
#include "ScalFmmConfig.h"
COULAUD Olivier's avatar
COULAUD Olivier committed
32 33
#include "Utils/FParameters.hpp"
#include "Files/FFmaGenericLoader.hpp"
34

COULAUD Olivier's avatar
COULAUD Olivier committed
35 36
#include "Kernels/Rotation/FRotationKernel.hpp"
#include "Kernels/Rotation/FRotationCell.hpp"
37

COULAUD Olivier's avatar
COULAUD Olivier committed
38 39
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
40 41


COULAUD Olivier's avatar
COULAUD Olivier committed
42
#include "Containers/FOctree.hpp"
43

COULAUD Olivier's avatar
COULAUD Olivier committed
44
#ifdef _OPENMP
COULAUD Olivier's avatar
COULAUD Olivier committed
45
#include "Core/FFmmAlgorithmThread.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
46 47 48
#else
#include "Core/FFmmAlgorithm.hpp"
#endif
49

50
#include "Utils/FParameterNames.hpp"
51

52 53 54 55 56 57 58 59 60 61 62 63 64 65
/// \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
//!
66
//!     \param   -f name          Name of the particles file with extension (.fma or .bfma). The data in  file have to be in our FMA format
67 68 69 70 71 72 73
//!
//


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

80
    const std::string defaultFile(/*SCALFMMDataPath+*/"../Data/test20k.fma");
81 82 83 84 85
    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);

86
#ifdef _OPENMP
COULAUD Olivier's avatar
COULAUD Olivier committed
87
	omp_set_num_threads(NbThreads) ;
88 89 90 91 92 93 94 95 96 97 98 99 100 101
	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;
102
	//
103
	// open particle file
104
	//
105
	////////////////////////////////////////////////////////////////////
COULAUD Olivier's avatar
COULAUD Olivier committed
106
	//
107 108
    typedef double FReal;
    FFmaGenericLoader<FReal> loader(filename);
COULAUD Olivier's avatar
COULAUD Olivier committed
109 110 111
	//
	////////////////////////////////////////////////////////////////////
	//
112
    // begin spherical kernel
113
	// accuracy
COULAUD Olivier's avatar
COULAUD Olivier committed
114
	const unsigned int P = 22;
115
	// typedefs
116 117 118 119
    typedef FP2PParticleContainerIndexed<FReal>                     ContainerClass;
    typedef FSimpleLeaf< FReal, ContainerClass >                       LeafClass;
    typedef FRotationCell<FReal, P>                                             CellClass;
    typedef FOctree<FReal, CellClass,ContainerClass,LeafClass>  OctreeClass;
120
	//
121
    typedef FRotationKernel< FReal, CellClass, ContainerClass , P>   KernelClass;
122
	//
COULAUD Olivier's avatar
COULAUD Olivier committed
123
#ifdef _OPENMP
124
	typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
COULAUD Olivier's avatar
COULAUD Olivier committed
125 126 127
#else
	typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
#endif
128
	// init oct-tree
129
	OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
130 131 132


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

			// 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
151
				<< time.elapsed() << "  s) ." << std::endl;
152 153 154 155 156 157 158 159 160
	} // -----------------------------------------------------

	{ // -----------------------------------------------------
		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
		//
161
		KernelClass *kernels = new KernelClass(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
162
		//
COULAUD Olivier's avatar
COULAUD Olivier committed
163
		FmmClass algorithm(&tree, kernels) ;
164 165 166 167
		//
		algorithm.execute();   // Here the call of the FMM algorithm
		//
		time.tac();
COULAUD Olivier's avatar
COULAUD Olivier committed
168
		std::cout << "Done  " << "(@Algorithm = " << time.elapsed() << " s) ." << std::endl;
169 170 171 172 173 174 175
	}
	// -----------------------------------------------------
	//
	// Some output
	//
	//
	{ // -----------------------------------------------------
176
        FSize N1=0, N2= loader.getNumberOfParticles()/2, N3= loader.getNumberOfParticles() -1; ;
177 178 179 180 181 182 183 184 185 186 187 188 189
		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();
190
            const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
191 192
			const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();

193
            const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
194

195 196
            for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                const FSize indexPartOrig = indexes[idxPart];
197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212
				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;
}