ChebyshevInterpolationMPIFMM.cpp 10.2 KB
Newer Older
COULAUD Olivier's avatar
COULAUD Olivier committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
// ===================================================================================
// 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 =====
// @FUSE_MPI
// @FUSE_BLAS
// ================

#include <iostream>
#include <stdexcept>
#include <cstdio>
#include <cstdlib>


#include "ScalFmmConfig.h"
29 30 31
#include "Containers/FOctree.hpp"
#include "Utils/FMpi.hpp"
#include "Core/FFmmAlgorithmThreadProc.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
32

33 34 35
#include "Files/FFmaGenericLoader.hpp"
#include "Files/FMpiFmaGenericLoader.hpp"
#include "Files/FMpiTreeBuilder.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
36

37
#include "BalanceTree/FLeafBalance.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
38

39 40 41
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Chebyshev/FChebSymKernel.hpp"
#include "Kernels/Chebyshev/FChebCell.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
42

43 44
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
45

46 47
#include "Utils/FParameters.hpp"
#include "Utils/FParameterNames.hpp"
COULAUD Olivier's avatar
COULAUD Olivier committed
48 49 50 51 52 53 54 55 56 57 58 59 60


/// \file ChebyshevInterpolationMPIFMM
//!
//! \brief This program runs the MPI FMM with Chebyshev interpolation of 1/r kernel
//!  \authors B. Bramas, O. Coulaud
//!
//!  This code is a short example to use the FMM Algorithm Proc with Chebyshev Interpolation for the 1/r kernel


// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
61
    FHelpDescribeAndExit(argc, argv,
62 63 64 65 66
			 "Driver for Chebyshev Interpolation kernel using MPI  (1/r kernel). "
			 "Usully run using : mpirun -np nb_proc_needed ./ChebyshevInterpolationAlgorithm [params].",
			 FParameterDefinitions::InputFile, FParameterDefinitions::OctreeHeight,
			 FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::InputFile,
			 FParameterDefinitions::NbThreads);
67

68
    const std::string defaultFile(/*SCALFMMDataPath+*/"../Data/test20k.fma");
69 70 71 72 73 74
    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);


COULAUD Olivier's avatar
COULAUD Olivier committed
75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93
#ifdef _OPENMP
	omp_set_num_threads(NbThreads);
	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 values for MPI
	FMpi app(argc,argv);
	//
	// init timer
	FTic time;

94 95
    typedef double FReal;
    FMpiFmaGenericLoader<FReal> loader(filename,app.global());
COULAUD Olivier's avatar
COULAUD Olivier committed
96 97 98 99 100 101

	if(!loader.isOpen()) throw std::runtime_error("Particle file couldn't be opened!") ;
	////////////////////////////////////////////////////////////////////



102
    // begin spherical kernel
COULAUD Olivier's avatar
COULAUD Olivier committed
103 104 105
	// accuracy
	const unsigned int ORDER = 7;
	// typedefs
106 107 108 109 110
    typedef FP2PParticleContainerIndexed<FReal>                      ContainerClass;
    typedef FSimpleLeaf<FReal, ContainerClass >                       LeafClass;
    typedef FChebCell<FReal,ORDER>                                    CellClass;
    typedef FOctree<FReal,CellClass,ContainerClass,LeafClass>         OctreeClass;
    typedef FInterpMatrixKernelR<FReal>                                MatrixKernelClass;
111
	const MatrixKernelClass MatrixKernel;
112
    typedef FChebSymKernel<FReal, CellClass,ContainerClass,MatrixKernelClass,ORDER>  KernelClass;
COULAUD Olivier's avatar
COULAUD Olivier committed
113 114 115 116 117 118

	//
	typedef FFmmAlgorithmThreadProc<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClassProc;

	// init oct-tree
	OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
119

COULAUD Olivier's avatar
COULAUD Olivier committed
120 121

	{ // -----------------------------------------------------
122 123 124 125 126
	  if(app.global().processId() == 0){
	    std::cout << "Creating & Inserting " << loader.getNumberOfParticles()
		      << " particles ..." << std::endl;
	    std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
	  }
COULAUD Olivier's avatar
COULAUD Olivier committed
127 128
	  time.tic();
	  //
129

COULAUD Olivier's avatar
COULAUD Olivier committed
130
	  struct TestParticle{
131
	    FSize index;
132
        FPoint<FReal> position;
COULAUD Olivier's avatar
COULAUD Olivier committed
133
	    FReal physicalValue;
134
        const FPoint<FReal>& getPosition(){
COULAUD Olivier's avatar
COULAUD Olivier committed
135 136 137
	      return position;
	    }
	  };
138

139 140
	  TestParticle* particles = new TestParticle[loader.getMyNumberOfParticles()];
	  memset(particles, 0, (unsigned int) (sizeof(TestParticle) * loader.getMyNumberOfParticles()));
141 142

	  //idx (in file) of the first part that will be used by this proc.
143 144
	  FSize idxStart = loader.getStart();
	  printf("Proc %d idxStart %lld \n",app.global().processId(),idxStart);
145

146
	  for(FSize idxPart = 0 ; idxPart < loader.getMyNumberOfParticles() ; ++idxPart){
COULAUD Olivier's avatar
COULAUD Olivier committed
147 148 149 150 151
	    //Storage of the index (in the original file) of each part.
	    particles[idxPart].index = idxPart + idxStart;
	    // Read particles from file
	    loader.fillParticle(&particles[idxPart].position,&particles[idxPart].physicalValue);
	  }
152

COULAUD Olivier's avatar
COULAUD Olivier committed
153
	  FVector<TestParticle> finalParticles;
154 155
      FLeafBalance balancer;
	  // FMpiTreeBuilder< FReal,TestParticle >::ArrayToTree(app.global(), particles, loader.getMyNumberOfParticles(),
156 157 158
	  //					 tree.getBoxCenter(),
	  //					 tree.getBoxWidth(),
	  //					 tree.getHeight(), &finalParticles,&balancer);
159
      FMpiTreeBuilder< FReal, TestParticle >::DistributeArrayToContainer(app.global(),particles,
PIACIBELLO Cyrille's avatar
PIACIBELLO Cyrille committed
160 161 162 163
								      loader.getMyNumberOfParticles(),
								      tree.getBoxCenter(),
								      tree.getBoxWidth(),tree.getHeight(),
								      &finalParticles, &balancer);
164

COULAUD Olivier's avatar
COULAUD Olivier committed
165 166 167 168 169
	  for(int idx = 0 ; idx < finalParticles.getSize(); ++idx){
	    tree.insert(finalParticles[idx].position,finalParticles[idx].index,finalParticles[idx].physicalValue);
	  }
	  printf("%d parts have been inserted in Tree \n",finalParticles.getSize());
	  delete[] particles;
170

COULAUD Olivier's avatar
COULAUD Olivier committed
171
	  time.tac();
172 173 174 175 176 177 178 179
	  double timeUsed = time.elapsed();
	  double minTime,maxTime;
	  std::cout << "Done  " << "(@Reading and Inserting Particles = "  << time.elapsed() << "s)." << std::endl;
	  MPI_Reduce(&timeUsed,&minTime,1,MPI_DOUBLE,MPI_MIN,0,app.global().getComm());
	  MPI_Reduce(&timeUsed,&maxTime,1,MPI_DOUBLE,MPI_MAX,0,app.global().getComm());
	  if(app.global().processId() == 0){
	      printf("Reading and Inserting Particles Time used : \t MIN : %f \t MAX %f in s\n",minTime,maxTime);
	  }
COULAUD Olivier's avatar
COULAUD Olivier committed
180 181 182 183
	} // -----------------------------------------------------

	{ // -----------------------------------------------------
	  std::cout << "\nChebyshev Interpolation  FMM Proc (P="<< ORDER << ") ... " << std::endl;
184

COULAUD Olivier's avatar
COULAUD Olivier committed
185 186 187 188
	  time.tic();
	  //
	  // Here we use a pointer due to the limited size of the stack
	  //
189
	  KernelClass *kernels = new KernelClass(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel);
COULAUD Olivier's avatar
COULAUD Olivier committed
190 191 192 193 194 195
	  //
	  FmmClassProc algorithm(app.global(),&tree, kernels);
	  //
	  algorithm.execute();   // Here the call of the FMM algorithm
	  //
	  time.tac();
196 197
	  double timeUsed = time.elapsed();
	  double minTime,maxTime;
COULAUD Olivier's avatar
COULAUD Olivier committed
198
	  std::cout << "Done  " << "(@Algorithm = " << time.elapsed() << " s)." << std::endl;
199 200 201 202 203
	  MPI_Reduce(&timeUsed,&minTime,1,MPI_DOUBLE,MPI_MIN,0,app.global().getComm());
	  MPI_Reduce(&timeUsed,&maxTime,1,MPI_DOUBLE,MPI_MAX,0,app.global().getComm());
	  if(app.global().processId() == 0){
	      printf("Time used : \t MIN : %f \t MAX %f \n",minTime,maxTime);
	  }
COULAUD Olivier's avatar
COULAUD Olivier committed
204 205 206 207 208 209 210
	}
	// -----------------------------------------------------
	//
	// Some output
	//
	//
	{ // -----------------------------------------------------
211
	  long int N1=0, N2= loader.getNumberOfParticles()/2, N3= (loader.getNumberOfParticles()-1); ;
COULAUD Olivier's avatar
COULAUD Olivier committed
212 213 214 215 216 217 218 219 220 221 222 223
	  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 posX = leaf->getTargets()->getPositions()[0];
	      const FReal*const posY = leaf->getTargets()->getPositions()[1];
	      const FReal*const posZ = leaf->getTargets()->getPositions()[2];
224

COULAUD Olivier's avatar
COULAUD Olivier committed
225 226 227 228 229 230
	      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 int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
	      const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
231

COULAUD Olivier's avatar
COULAUD Olivier committed
232
	      const FVector<int>& indexes = leaf->getTargets()->getIndexes();
233

COULAUD Olivier's avatar
COULAUD Olivier committed
234 235 236 237
	      for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
		const int indexPartOrig = indexes[idxPart];
		if ((indexPartOrig == N1) || (indexPartOrig == N2) || (indexPartOrig == N3)  ) {
		  std::cout << "Proc "<< app.global().processId() << " Index "<< indexPartOrig <<"  potential  " << potentials[idxPart]
238
			    << " Pos "<<posX[idxPart]<<" "<<posY[idxPart]<<" "<<posZ[idxPart]
COULAUD Olivier's avatar
COULAUD Olivier committed
239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254
			    << "   Forces: " << forcesX[idxPart] << " " << forcesY[idxPart] << " "<< forcesZ[idxPart] <<std::endl;
		}
		energy += potentials[idxPart]*physicalValues[idxPart] ;
	      }
	    });
	  FReal gloEnergy = app.global().reduceSum(energy);
	  if(0 == app.global().processId()){
	    std::cout <<std::endl << "Proc "<< app.global().processId() << " Energy: "<< gloEnergy <<std::endl;
	    std::cout <<std::endl <<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl<<std::endl;
	  }
	}
	// -----------------------------------------------------


	return 0;
}