RotationFMM.cpp 7.76 KB
Newer Older
1
// ===================================================================================
2 3
// Copyright ScalFmm 2016 INRIA
//
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
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/
14 15 16 17 18 19 20
// ===================================================================================

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

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

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

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

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


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

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

45
#include "Utils/FParameterNames.hpp"
46

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


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

75
    const std::string defaultFile(/*SCALFMMDataPath+*/"../Data/test20k.fma");
76 77 78 79 80
    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);

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


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

			// 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
146
				<< time.elapsed() << "  s) ." << std::endl;
147 148 149 150 151 152 153 154 155
	} // -----------------------------------------------------

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

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

190 191
            for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
                const FSize indexPartOrig = indexes[idxPart];
192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207
				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;
}