RotationFMM.cpp 9.2 KB
Newer Older
1
// See LICENCE file at project root
2 3 4 5 6 7

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

#include <iostream>
8
#include <stdexcept>
COULAUD Olivier's avatar
COULAUD Olivier committed
9
#include <string>
10
#include <cstdlib>
11
#include <cstdio>
12

COULAUD Olivier's avatar
COULAUD Olivier committed
13
#include "ScalFmmConfig.h"
14 15
#include "Utils/FParameters.hpp"
#include "Files/FFmaGenericLoader.hpp"
16

17 18
#include "Kernels/Rotation/FRotationKernel.hpp"
#include "Kernels/Rotation/FRotationCell.hpp"
19

20 21
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
22 23


24
#include "Containers/FOctree.hpp"
25

COULAUD Olivier's avatar
COULAUD Olivier committed
26
#ifdef _OPENMP
27
#include "Core/FFmmAlgorithmThread.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
28 29 30
#else
#include "Core/FFmmAlgorithm.hpp"
#endif
31

32
#include "Utils/FParameterNames.hpp"
33

34 35 36 37 38 39 40 41 42 43 44 45 46 47
/// \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
//!
48
//!     \param   -f name          Name of the particles file with extension (.fma or .bfma). The data in  file have to be in our FMA format
49 50 51 52 53 54 55
//!
//


// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
56 57 58
    FHelpDescribeAndExit(argc, argv,
                         "Driver for HArmonic Spherical + Rotation  --  kernel  (1/r kernel).",
                         FParameterDefinitions::InputFile, FParameterDefinitions::OctreeHeight,
59
                         FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::OutputFile,
60 61
                         FParameterDefinitions::NbThreads);

62
    const std::string defaultFile(/*SCALFMMDataPath+*/"../Data/test20k.fma");
63 64 65 66 67
    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);

68
#ifdef _OPENMP
69
	omp_set_num_threads(NbThreads) ;
70 71 72 73 74 75 76 77 78 79 80 81 82 83
	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;
84
	//
85
	// open particle file
86
	//
87
	////////////////////////////////////////////////////////////////////
COULAUD Olivier's avatar
COULAUD Olivier committed
88
	//
89 90
    typedef double FReal;
    FFmaGenericLoader<FReal> loader(filename);
COULAUD Olivier's avatar
COULAUD Olivier committed
91 92 93
	//
	////////////////////////////////////////////////////////////////////
	//
94
    // begin spherical kernel
95
	// accuracy
COULAUD Olivier's avatar
COULAUD Olivier committed
96
	const unsigned int P = 22;
97
	// typedefs
98 99 100 101
    typedef FP2PParticleContainerIndexed<FReal>                     ContainerClass;
    typedef FSimpleLeaf< FReal, ContainerClass >                       LeafClass;
    typedef FRotationCell<FReal, P>                                             CellClass;
    typedef FOctree<FReal, CellClass,ContainerClass,LeafClass>  OctreeClass;
102
	//
103
    typedef FRotationKernel< FReal, CellClass, ContainerClass , P>   KernelClass;
104
	//
COULAUD Olivier's avatar
COULAUD Olivier committed
105
#ifdef _OPENMP
106
	typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
COULAUD Olivier's avatar
COULAUD Olivier committed
107 108 109
#else
	typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
#endif
110
	// init oct-tree
111
	OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
112 113 114


	{ // -----------------------------------------------------
115
		std::cout << "Creating & Inserting " << loader.getNumberOfParticles()
COULAUD Olivier's avatar
COULAUD Olivier committed
116
                																<< " particles ..." << std::endl;
117 118 119
		std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
		time.tic();
		//
120
        FPoint<FReal> position;
121 122
		FReal physicalValue = 0.0;
		//
123
        for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
124
			// Read particles from file
125
			loader.fillParticle(&position,&physicalValue);
126 127 128 129 130 131 132

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

		time.tac();
		std::cout << "Done  " << "(@Creating and Inserting Particles = "
133
				<< time.elapsed() << "  s) ." << std::endl;
134 135 136 137 138 139 140 141 142
	} // -----------------------------------------------------

	{ // -----------------------------------------------------
		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
		//
143
		KernelClass *kernels = new KernelClass(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
144
		//
145
		FmmClass algorithm(&tree, kernels) ;
146 147 148 149
		//
		algorithm.execute();   // Here the call of the FMM algorithm
		//
		time.tac();
150
		std::cout << "Done  " << "(@Algorithm = " << time.elapsed() << " s) ." << std::endl;
151 152 153 154 155 156 157
	}
	// -----------------------------------------------------
	//
	// Some output
	//
	//
	{ // -----------------------------------------------------
158
        FSize N1=0, N2= loader.getNumberOfParticles()/2, N3= loader.getNumberOfParticles() -1; ;
159 160 161 162 163 164 165 166 167 168 169 170 171
		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();
172
            const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
173 174
			const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();

175
            const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
176

177 178
            for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                const FSize indexPartOrig = indexes[idxPart];
179 180 181 182 183 184 185 186 187 188 189 190
				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;

	}
	// -----------------------------------------------------
191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
	// -----------------------------------------------------
	if(FParameters::existParameter(argc, argv, FParameterDefinitions::OutputFile.options)){
	  std::string name(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options,   "output.fma"));
	  FFmaGenericWriter<FReal> writer(name) ;
	  //
	  FSize NbPoints = loader.getNumberOfParticles();
	  FReal * particles ;
	  particles = new FReal[8*NbPoints] ;
	  memset(particles,0,8*NbPoints*sizeof(FReal));
	  FSize j = 0 ;
	  tree.forEachLeaf([&](LeafClass* leaf){
	      //
	      // Input
	      const FReal*const posX = leaf->getTargets()->getPositions()[0];
	      const FReal*const posY = leaf->getTargets()->getPositions()[1];
	      const FReal*const posZ = leaf->getTargets()->getPositions()[2];
	      const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
	      const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
	      //
	      // Computed data
	      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();
	      //
	      const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
	      for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
		j = 8*indexes[idxPart];
		particles[j]      = posX[idxPart] ;
		particles[j+1]  = posY[idxPart] ;
		particles[j+2]  = posZ[idxPart] ;
		particles[j+3]  = physicalValues[idxPart] ;
		particles[j+4]  = potentials[idxPart] ;
		particles[j+5]  =  forcesX[idxPart] ;
		particles[j+6]  =  forcesY[idxPart] ;
		particles[j+7]  =  forcesZ[idxPart] ;
	      }
	    });

          writer.writeHeader( loader.getCenterOfBox(), loader.getBoxWidth() ,  NbPoints, sizeof(FReal), 8) ;
          writer.writeArrayOfReal(particles,  8 , NbPoints);

          delete[] particles;

          //
          std::string name1( "output.fma");
          //
          FFmaGenericWriter<FReal> writer1(name1) ;
          writer1.writeDistributionOfParticlesFromOctree(&tree,NbPoints) ;
        }
241 242 243

	return 0;
}