From 89e17b154b47777e00e06c209288a7e314c7b221 Mon Sep 17 00:00:00 2001 From: Olivier Coulaud Date: Tue, 29 May 2018 10:38:43 +0200 Subject: [PATCH] Bugs with periodic conditions seem fixed. More tests should be done. --- Examples/MPIInterpolationFMM.hpp | 304 ++++++++-------- Src/Core/FCoreCommon.hpp | 11 +- Src/Core/FFmmAlgorithmThreadProc.hpp | 10 +- Src/Core/FFmmAlgorithmThreadProcPeriodic.hpp | 96 +++-- Src/Kernels/P2P/FP2PTensorialKij.hpp | 2 +- Tests/Utils/generateDistrib.cpp | 356 +++++++++++++++++++ Tests/Utils/testFmmAlgorithmProcPeriodic.cpp | 2 +- 7 files changed, 600 insertions(+), 181 deletions(-) create mode 100644 Tests/Utils/generateDistrib.cpp diff --git a/Examples/MPIInterpolationFMM.hpp b/Examples/MPIInterpolationFMM.hpp index 353f9760..8455f593 100644 --- a/Examples/MPIInterpolationFMM.hpp +++ b/Examples/MPIInterpolationFMM.hpp @@ -80,9 +80,9 @@ int main(int argc, char* argv[]) FParameterDefinitions::OutputFile, FParameterDefinitions::NbThreads, FParameterDefinitions::PeriodicityNbLevels, - localIncreaseBox - ) -; + localIncreaseBox + ) + ; const std::string defaultFile("../Data/test20k.fma"); const std::string filename = FParameters::getStr(argc,argv, FParameterDefinitions::InputFile.options, defaultFile.c_str()); @@ -127,98 +127,98 @@ int main(int argc, char* argv[]) auto boxWidth = loader.getBoxWidth() ; // if(FParameters::existParameter(argc, argv, localIncreaseBox.options)){ - FReal ratio= FParameters::getValue(argc, argv, localIncreaseBox.options, 1.0); - boxWidth *= ratio; - } - + FReal ratio= FParameters::getValue(argc, argv, localIncreaseBox.options, 1.0); + boxWidth *= ratio; + } + // Initialize empty oct-tree OctreeClass tree(TreeHeight, SubTreeHeight, boxWidth, loader.getCenterOfBox()); FSize localParticlesNumber = 0 ; // ----------------------------------------------------- - if(app.global().processId() == 0){ - std::cout << "Loading & Inserting " << loader.getNumberOfParticles() - << " particles ..." << std::endl - <<" Box: "<< std::endl - << " width " << boxWidth << std::endl - << " Centre " << loader.getCenterOfBox()<< std::endl; - std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl; - } - time.tic(); + if(app.global().processId() == 0){ + std::cout << "Loading & Inserting " << loader.getNumberOfParticles() + << " particles ..." << std::endl + <<" Box: "<< std::endl + << " width " << boxWidth << std::endl + << " Centre " << loader.getCenterOfBox()<< std::endl; + std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl; + } + time.tic(); + + /* Mock particle structure to balance the tree over the processes. */ + struct TestParticle{ + FSize index; // Index of the particle in the original file. + FPoint position; // Spatial position of the particle. + FReal physicalValue; // Physical value of the particle. + /* Returns the particle position. */ + const FPoint& getPosition(){ + return position; + } + }; + + // Temporary array of particles read by this process. + TestParticle* particles = new TestParticle[loader.getMyNumberOfParticles()]; + memset(particles, 0, (sizeof(TestParticle) * loader.getMyNumberOfParticles())); + + // Index (in file) of the first particle that will be read by this process. + FSize idxStart = loader.getStart(); + std::cout << "Proc:" << app.global().processId() << " start-index: " << idxStart << std::endl; + + // Read particles from parts. + for(FSize idxPart = 0 ; idxPart < loader.getMyNumberOfParticles() ; ++idxPart){ + // Store the index (in the original file) the particle. + particles[idxPart].index = idxPart + idxStart; + // Read particle from file + loader.fillParticle(&particles[idxPart].position, + &particles[idxPart].physicalValue); + } - /* Mock particle structure to balance the tree over the processes. */ - struct TestParticle{ - FSize index; // Index of the particle in the original file. - FPoint position; // Spatial position of the particle. - FReal physicalValue; // Physical value of the particle. - /* Returns the particle position. */ - const FPoint& getPosition(){ - return position; - } - }; - - // Temporary array of particles read by this process. - TestParticle* particles = new TestParticle[loader.getMyNumberOfParticles()]; - memset(particles, 0, (sizeof(TestParticle) * loader.getMyNumberOfParticles())); - - // Index (in file) of the first particle that will be read by this process. - FSize idxStart = loader.getStart(); - std::cout << "Proc:" << app.global().processId() << " start-index: " << idxStart << std::endl; - - // Read particles from parts. - for(FSize idxPart = 0 ; idxPart < loader.getMyNumberOfParticles() ; ++idxPart){ - // Store the index (in the original file) the particle. - particles[idxPart].index = idxPart + idxStart; - // Read particle from file - loader.fillParticle(&particles[idxPart].position, - &particles[idxPart].physicalValue); - } + // Final vector of particles + FVector finalParticles; + FLeafBalance balancer; + // Redistribute particules between processes + FMpiTreeBuilder< FReal, TestParticle >:: + DistributeArrayToContainer(app.global(), + particles, + loader.getMyNumberOfParticles(), + tree.getBoxCenter(), + tree.getBoxWidth(), + tree.getHeight(), + &finalParticles, + &balancer); + + // Free temporary array memory. + delete[] particles; + + // Insert final particles into tree. + + for(FSize idx = 0 ; idx < finalParticles.getSize(); ++idx){ + tree.insert(finalParticles[idx].position, + finalParticles[idx].index, + finalParticles[idx].physicalValue); + } - // Final vector of particles - FVector finalParticles; - FLeafBalance balancer; - // Redistribute particules between processes - FMpiTreeBuilder< FReal, TestParticle >:: - DistributeArrayToContainer(app.global(), - particles, - loader.getMyNumberOfParticles(), - tree.getBoxCenter(), - tree.getBoxWidth(), - tree.getHeight(), - &finalParticles, - &balancer); - - // Free temporary array memory. - delete[] particles; - - // Insert final particles into tree. - - for(FSize idx = 0 ; idx < finalParticles.getSize(); ++idx){ - tree.insert(finalParticles[idx].position, - finalParticles[idx].index, - finalParticles[idx].physicalValue); - } + time.tac(); - time.tac(); + localParticlesNumber = finalParticles.getSize() ; - localParticlesNumber = finalParticles.getSize() ; + double timeUsed = time.elapsed(); + double minTime,maxTime; + std::cout << "Proc:" << app.global().processId() + << " " << finalParticles.getSize() + << " particles have been inserted in the tree. (@Reading and Inserting Particles = " + << time.elapsed() << " s)." + << std::endl; - double timeUsed = time.elapsed(); - double minTime,maxTime; - std::cout << "Proc:" << app.global().processId() - << " " << finalParticles.getSize() - << " particles have been inserted in the tree. (@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){ - std::cout << "readinsert-time-min:" << minTime - << " readinsert-time-max:" << maxTime - << 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){ + std::cout << "readinsert-time-min:" << minTime + << " readinsert-time-max:" << maxTime + << std::endl; + } // ----------------------------------------------------- std::vector mortonLeafDistribution(2*app.global().processCount()); @@ -255,7 +255,25 @@ int main(int argc, char* argv[]) // FMM exectution FFmmFarField FFmmNearField algorithm->execute(FFmmNearField); // + algorithm->getMortonLeafDistribution(mortonLeafDistribution); +// auto nbProcess = app.global().processCount() ; +// auto level = tree.getHeight()- 1; +// for (int i =0; i < nbProcess ; ++i){ +// auto I = algoPer.getWorkingInterval(level, i); +// mortonLeafDistribution[2*i] = I.leftIndex; +// mortonLeafDistribution[2*i+1] = I.rightIndex; +// } + + // if(app.global().processId() == 0) + { + std::cout << app.global().processId() <<" Morton distribution "<< std::endl ; + for(auto v : mortonLeafDistribution) + std::cout << v << " "; + + std::cout << std::endl; + } + app.global().barrier(); time.tac(); timeUsed = time.elapsed(); std::cout << "Done " << "(@Algorithm = " << time.elapsed() << " s)." << std::endl; @@ -268,70 +286,62 @@ int main(int argc, char* argv[]) } } -// ----------------------------------------------------- -// -// Some output -// -// -{ // ----------------------------------------------------- -FSize N1=0, N2= loader.getNumberOfParticles()/2, N3= (loader.getNumberOfParticles()-1); ; -FReal energy =0.0 ; -// -// Loop over all leaves -// -std::cout <getTargets()->getPositions()[0]; - const FReal*const posY = leaf->getTargets()->getPositions()[1]; - const FReal*const posZ = leaf->getTargets()->getPositions()[2]; - - 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(); - const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues(); - - const FVector& indexes = leaf->getTargets()->getIndexes(); - - for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){ - const FSize indexPartOrig = indexes[idxPart]; - if ((indexPartOrig == N1) || (indexPartOrig == N2) || (indexPartOrig == N3) ) { - std::cout << "Proc "<< app.global().processId() << " Index "<< indexPartOrig <<" potential " << potentials[idxPart] - << " Pos "<getAndBuildMortonIndexAtLeaf(mortonLeafDistribution); + // ----------------------------------------------------- // - if(app.global().processId() == 0) - { - std::cout << " Morton distribution "<< std::endl ; - for(auto v : mortonLeafDistribution) - std::cout << v << " "; + // Some output + // + // + { // ----------------------------------------------------- + FSize N1=0, N2= loader.getNumberOfParticles()/2, N3= (loader.getNumberOfParticles()-1); ; + FReal energy =0.0 ; + // + // Loop over all leaves + // + std::cout <getTargets()->getPositions()[0]; + const FReal*const posY = leaf->getTargets()->getPositions()[1]; + const FReal*const posZ = leaf->getTargets()->getPositions()[2]; + + 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(); + const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues(); + + const FVector& indexes = leaf->getTargets()->getIndexes(); + + for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){ + const FSize indexPartOrig = indexes[idxPart]; + if ((indexPartOrig == N1) || (indexPartOrig == N2) || (indexPartOrig == N3) ) + { + std::cout << "Proc "<< app.global().processId() << " Index "<< indexPartOrig <<" potential " << potentials[idxPart] + << " Pos "< paraWriter(name,app); + paraWriter.writeDistributionOfParticlesFromOctree(tree,loader.getNumberOfParticles(),localParticlesNumber,mortonLeafDistribution); - std::cout << std::endl; } - std::string name(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options, "output.fma")); - FMpiFmaGenericWriter paraWriter(name,app); - paraWriter.writeDistributionOfParticlesFromOctree(tree,loader.getNumberOfParticles(),localParticlesNumber,mortonLeafDistribution); - -} -return 0; + return 0; } diff --git a/Src/Core/FCoreCommon.hpp b/Src/Core/FCoreCommon.hpp index 706b1b71..171b74b9 100644 --- a/Src/Core/FCoreCommon.hpp +++ b/Src/Core/FCoreCommon.hpp @@ -151,12 +151,15 @@ public: executeCore(operationsToProceed); } - #ifdef SCALFMM_USE_MPI /// Build and dill vector of the MortonIndex Distribution at Leaf level /// p = mpi process id then - /// [mortonLeafDistribution[2*p], mortonLeafDistribution[2*p+1] is the morton index shared by process p - virtual void getAndBuildMortonIndexAtLeaf(std::vector & mortonLeafDistribution) {} ; -#endif + /// [mortonLeafDistribution[2*p], mortonLeafDistribution[2*p+1] is the morton index shared by process p + virtual void getMortonLeafDistribution(std::vector & mortonLeafDistribution) { + mortonLeafDistribution.resize(2) ; + mortonLeafDistribution[0] = 0 ; + mortonLeafDistribution[1] = 8 << (nbLevelsInTree-1) ; + + }; }; diff --git a/Src/Core/FFmmAlgorithmThreadProc.hpp b/Src/Core/FFmmAlgorithmThreadProc.hpp index 6f646803..891e4284 100644 --- a/Src/Core/FFmmAlgorithmThreadProc.hpp +++ b/Src/Core/FFmmAlgorithmThreadProc.hpp @@ -2,6 +2,8 @@ #ifndef FFMMALGORITHMTHREADPROC_HPP #define FFMMALGORITHMTHREADPROC_HPP +#define SCALFMM_DISTRIBUTED_ALGORITHM + #include #include // @@ -138,9 +140,11 @@ public: /// Build and dill vector of the MortonIndex Distribution at Leaf level /// p = mpi process id then /// [mortonLeafDistribution[2*p], mortonLeafDistribution[2*p+1] is the morton index shared by process p - void getAndBuildMortonIndexAtLeaf(std::vector & mortonLeafDistribution) final { - for (int p=0 ; p< comm.processCount() ; ++p ){ - auto inter = getWorkingInterval(OctreeHeight-1, p ); + void getMortonLeafDistribution(std::vector & mortonLeafDistribution) final { + mortonLeafDistribution.resize(2*nbProcess) ; + auto level = OctreeHeight - 1; + for (int p=0 ; p< nbProcess ; ++p ){ + auto inter = this->getWorkingInterval(level, p ); mortonLeafDistribution[2*p] = inter.leftIndex; mortonLeafDistribution[2*p+1] = inter.rightIndex; } diff --git a/Src/Core/FFmmAlgorithmThreadProcPeriodic.hpp b/Src/Core/FFmmAlgorithmThreadProcPeriodic.hpp index 26340663..70d863b7 100644 --- a/Src/Core/FFmmAlgorithmThreadProcPeriodic.hpp +++ b/Src/Core/FFmmAlgorithmThreadProcPeriodic.hpp @@ -2,6 +2,8 @@ #ifndef FFMMALGORITHMTHREADPROCPPERIODIC_HPP #define FFMMALGORITHMTHREADPROCPPERIODIC_HPP +#define SCALFMM_DISTRIBUTED_ALGORITHM + #include #include #include @@ -28,6 +30,7 @@ #ifdef _OPENMP #include #endif + #include "FCoreCommon.hpp" #include "FP2PExclusion.hpp" @@ -95,6 +98,25 @@ public: using local_expansion_t = typename CellClass::local_expansion_t; using symbolic_data_t = CellClass; + + /// Get the Morton index Distribution at the leaf level + /// + /// Fill the vector mortonDistribution + /// + /// p = mpi process id then + /// Processor p owns indexes between [mortonLeafDistribution[2*p], mortonLeafDistribution[2*p]+1] + /// + /// parameter[out] mortonLeafDistribution + /// + void getMortonLeafDistribution(std::vector & mortonLeafDistribution) final { + mortonLeafDistribution.resize(2*nbProcess) ; + auto level = OctreeHeight - 1; + for (int p=0 ; p< nbProcess ; ++p ){ + auto inter = this->getWorkingInterval(level, p ); + mortonLeafDistribution[2*p] = inter.leftIndex; + mortonLeafDistribution[2*p+1] = inter.rightIndex; + } + } private: /// Get an interval from a process id and tree level Interval& getWorkingInterval( int level, int proc){ @@ -105,7 +127,6 @@ private: const Interval& getWorkingInterval( int level, int proc) const { return workingIntervalsPerLevel[OctreeHeight * proc + level]; } - /// Does \a procIdx have work at given \a idxLevel /** i.e. does it hold cells and is responsible of them ? */ bool procHasWorkAtLevel(const int idxLevel , const int idxProc) const { @@ -134,7 +155,14 @@ public: } } } +// /// Build and dill vector of the MortonIndex Distribution at Leaf level +// /// p = mpi process id then +// /// [mortonLeafDistribution[2*p], mortonLeafDistribution[2*p+1] is the morton index shared by process p +// virtual void getAndBuildMortonIndexAtLeaf(std::vector & mortonLeafDistribution) { +// mortonLeafDistribution.resize(2*nbProcess) ; + +// } ; Interval& getWorkingInterval(const int level){ return getWorkingInterval(level, idProcess); @@ -1697,13 +1725,15 @@ protected: if(needOther){ //means that something need to be sent (or received) leafsNeedOther.set(idxLeaf,true); ++countNeedOther; - // std::cout << " Leaf "<< iterArray[idxLeaf].getCurrentGlobalCoordinate().getMortonIndex() << std::endl ; - // for(int i = 0 ; i < nbProcess ; ++i ){ - // std::cout << " "<< i << " alreadySent " << alreadySent[i]<< " " << toSend[i].size() <<" : "; - // for(int j = 0 ; j < toSend[i].size() ; ++j ){ - // std::cout << " " << toSend[i][j].getCurrentGlobalCoordinate().getMortonIndex() ; - // } - // std::cout << std::endl; +// std::cout << " Leaf "<< iterArray[idxLeaf].getCurrentGlobalCoordinate().getMortonIndex() << std::endl ; +// for(int i = 0 ; i < nbProcess ; ++i ){ +// std::cout << " "<< i << " alreadySent " << alreadySent[i]<< " " << toSend[i].size() <<" : "; +// for(int j = 0 ; j < toSend[i].size() ; ++j ){ +// std::cout << " " << toSend[i][j].getCurrentGlobalCoordinate().getMortonIndex() ; +// } +// std::cout << std::endl; +// } + } } @@ -1818,7 +1848,7 @@ protected: // Check the new tree // std::cout <<"otherP2Ptree " << std::endl; -// otherP2Ptree.forEachCellLeaf([&](CellClass* LeafCell, LeafClass* leaf){ +// otherP2Ptree.forEachCellLeaf([&](CellClass* LeafCell, LeafClass* leaf){ // std::cout << LeafCell->getMortonIndex() <<" leaf " << leaf << " source " <getSrc() << std::endl; // } ); @@ -1932,7 +1962,7 @@ protected: // if(hasPeriodicLeaves){ // Current cell has periodic cells - constexpr int maxPeriodicNeighbors = 26; + constexpr int maxPeriodicNeighbors = 19; ContainerClass* periodicNeighbors[maxPeriodicNeighbors]{}; int periodicNeighborPositions[maxPeriodicNeighbors]; // std::array periodicOffsets{} ; @@ -1951,10 +1981,14 @@ protected: FReal yoffset= boxWidth * FReal(offsets[idxNeig].getY()) ; FReal zoffset= boxWidth * FReal(offsets[idxNeig].getZ()) ; for(FSize idxPart = 0; idxPart < periodicNeighbors[periodicNeighborsCounter]->getNbParticles() ; ++idxPart){ + // std::cout << " Xo " << positionsX[idxPart] << " " << positionsY[idxPart] << " " << positionsZ[idxPart] <::max()); const int nbLeafToProceed = int(leafsNeedOtherData.getSize()); - + // std::cout << "Treat communications " << nbLeafToProceed <getNeighborsIndexesPeriodic(currentIter.coord,limite, indexesNeighbors,indexArray,AllDirs); // std::cout << currentIter.coord.getMortonIndex() << " Cells comming from another process " << nbNeigh < periodicNeighbors{}; + std::array periodicNeighbors{}; int nbPeriodicCells=0 ; - - // Pourquoi pas de decallage de boite ??, + // std::cout <getMortonIndex() << " " << currentIter.coord << " " << nbNeigh <getPositions()[0]; FReal*const positionsY = periodicNeighbors[nbPeriodicCells]->getPositions()[1]; @@ -2080,26 +2115,37 @@ protected: positionsZ[idxPart] += offSet[2]; } #ifdef toto - std::cout << " MOVE CELLS "<getPositions()[0]; FReal*const positionsYX = hypotheticNeighbor->getPositions()[1]; FReal*const positionsZX = hypotheticNeighbor->getPositions()[2]; auto idxPart = 0 ; - std::cout <"; - if(counter){ // neighborPositions doesn't use) + if(counter> 0){ // neighborPositions doesn't use) + // std::cout << " counter: "<< counter <<" nbPeriodicCells " << nbPeriodicCells< P2PRemote" <0){ //to Do for(int i=0 ; i < nbPeriodicCells ; ++i){ diff --git a/Src/Kernels/P2P/FP2PTensorialKij.hpp b/Src/Kernels/P2P/FP2PTensorialKij.hpp index 73d4b2b7..05c576f3 100644 --- a/Src/Kernels/P2P/FP2PTensorialKij.hpp +++ b/Src/Kernels/P2P/FP2PTensorialKij.hpp @@ -227,7 +227,7 @@ inline void FullRemoteKIJ(ContainerClass* const FRestrict inTargets, const Conta // get information on tensorial aspect of matrix kernel const int ncmp = MatrixKernelClass::NCMP; - const int npv = MatrixKernelClass::NPV; + const int npv = MatrixKernelClass::NPV; const int npot = MatrixKernelClass::NPOT; // get info on targets diff --git a/Tests/Utils/generateDistrib.cpp b/Tests/Utils/generateDistrib.cpp new file mode 100644 index 00000000..7fc6b1c1 --- /dev/null +++ b/Tests/Utils/generateDistrib.cpp @@ -0,0 +1,356 @@ +// See LICENCE file at project root + +// ==== CMAKE ===== + +// ================ +// +// make testPeriodicAlgorithm +// mpirun --oversubscribe -np 4 --tag-output Tests/Release/testPeriodicAlgorithm -h 3 -sh 1 + + +#include +#include +#include +#include +#include +//#include +#include + +// +// Specific part +// +#include +// +// Generic part +// +#include "ScalFmmConfig.h" + +#include "Files/FFmaGenericLoader.hpp" + + +#include "Utils/FParameters.hpp" +#include "Utils/FParameterNames.hpp" + +// +// Order of the Interpolation approximation +using FReal_t = double; + +/** This program show an example of use of + * the fmm basic algo + * it also check that each particles is impacted each other particles + */ + +/* Mock particle structure to balance the tree over the processes. */ +struct TestParticle{ + FSize index; // Index of the particle in the original file. + FPoint position; // Spatial position of the particle. + FReal_t physicalValue; // Physical value of the particle. + /* Returns the particle position. */ + const FPoint& getPosition(){ + return position; + } +}; + + +void setOneParticleInEachCell(const FReal_t &BoxWidth , const int &Height, std::vector &particles, + MortonIndex & startIdx, MortonIndex & endIdx, int nbProc = 1, int myId = 0){ + // + auto leavesLevel = Height- 1 ; + auto BoxCellWidth = BoxWidth/ (2 << (leavesLevel-1)) ; + auto nbLeaves = 2 << (leavesLevel-1) ; + nbLeaves = nbLeaves*nbLeaves*nbLeaves ; + std::cout << " H: " << leavesLevel << " nbCells: " << nbLeaves << " cellwidth " << BoxCellWidth << std::endl ; + auto bloc = nbLeaves / nbProc ; + startIdx = bloc*myId ; + endIdx = startIdx+bloc; + if(myId == nbProc-1){ + endIdx = nbLeaves ;} + particles.resize(endIdx-startIdx) ; + FReal_t val = 0.1 *( myId %2 ==0 ? 1 : -1 ); + std::cout << "procId " << myId << " MortonIdx: [ " << startIdx << " , " << endIdx <<" [" < >& Lines, const FReal_t &BoxWidth , const int &Height, std::vector &particles, + MortonIndex & startIdx, MortonIndex & endIdx, int nbProc = 1, int myId = 0){ + // + auto leavesLevel = Height - 1 ; + auto BoxCellWidth = BoxWidth/ (2 << (leavesLevel-1)) ; + auto nb1DLeaves = 2 << (leavesLevel-1) ; + auto nbLines = Lines.size(); + std::vector cells(nb1DLeaves*nbLines) ; + // fill cells + std::cout << "Cells " << std::endl ; + + for (auto v : Lines){ + std::cout << " ("<< v[0]<< " "<< v[1] << ") "; + } + std::cout << std::endl; + std::size_t pos = 0 ; + for(std::size_t i = 0 ; i < nbLines ; ++i ) { + auto & index = Lines[i] ; + for(MortonIndex idx = 0; idx < nb1DLeaves ; ++idx ) { + FTreeCoordinate coor(index[0],index[1],idx) ; + cells[pos] = coor.getMortonIndex() ; + ++pos ; + } + } + std::cout << "Cells " << std::endl ; + + for (auto v : cells){ + std::cout << " "<< v; + } + std::cout << std::endl ; + std::sort(cells.begin(), cells.end()) ; + std::cout << "Sorted Cells " << std::endl ; + for (auto v : cells){ + std::cout << " "<< v; + } + std::cout << std::endl ; + // nbLeaves = nbLeaves*nbLeaves*nbLeaves ; + std::cout << " H: " << leavesLevel << " nb1DLeaves "<< nb1DLeaves + <<" nbCells: " << cells.size() << " cellwidth " << BoxCellWidth << std::endl ; + auto nbLeaves = cells.size() ; + auto bloc = nbLeaves / nbProc ; + startIdx = bloc*myId ; + endIdx = startIdx+bloc; + if(myId == nbProc-1){ + endIdx = nbLeaves ;} + particles.resize(endIdx-startIdx) ; + FReal_t val ; + if(particles.size() ==1) { + val = 0.1*( myId %2 ==0 ? 1 : -1 ); + } + else { + val = 0.1; + + } + std::cout << "procId " << myId << " MortonIdx: [ " << startIdx << " , " << endIdx <<" [" < &particles, + MortonIndex & startIdx, MortonIndex & endIdx, int nbProc = 1, int myId = 0){ + // + auto leavesLevel = Height - 1 ; + auto BoxCellWidth = BoxWidth/ (2 << (leavesLevel-1)) ; + auto nbLeaves = 2 << (leavesLevel-1) ; + // nbLeaves = nbLeaves*nbLeaves*nbLeaves ; + std::cout << " H: " << leavesLevel << " nbCells: " << nbLeaves << " cellwidth " << BoxCellWidth << std::endl ; + auto bloc = nbLeaves / nbProc ; + startIdx = bloc*myId ; + endIdx = startIdx+bloc; + if(myId == nbProc-1){ + endIdx = nbLeaves ;} + particles.resize(endIdx-startIdx) ; + FReal_t val ; + if(particles.size() ==1) { + val = 0.1*( myId %2 ==0 ? 1 : -1 ); + } + else { + val = 0.1; + + } + std::cout << "procId " << myId << " MortonIdx: [ " << startIdx << " , " << endIdx <<" [" < &particles, + MortonIndex & startIdx, MortonIndex & endIdx, int nbProc = 1, int myId = 0){ + // + auto leavesLevel = tree.getHeight() - 1 ; + auto BoxCellWidth = tree .getBoxWidth()/ (2 << (leavesLevel-1)) ; + auto nbLeaves = 2 << (leavesLevel-1) ; + // nbLeaves = nbLeaves*nbLeaves*nbLeaves ; + std::cout << " H: " << leavesLevel << " nbCells: " << nbLeaves << " cellwidth " << BoxCellWidth << std::endl ; + auto bloc = nbLeaves / nbProc ; + startIdx = bloc*myId ; + endIdx = startIdx+bloc; + if(myId == nbProc-1){ + endIdx = nbLeaves ;} + particles.resize(endIdx-startIdx) ; + FReal_t val ; + if(particles.size() ==1) { + val = 0.1*( myId %2 ==0 ? 1 : -1 ); + } + else { + val = 0.1; + + } + std::cout << "procId " << myId << " MortonIdx: [ " << startIdx << " , " << endIdx <<" [" < &particles, + MortonIndex & startIdx, MortonIndex & endIdx, int nbProc = 1, int myId = 0){ + // + auto leavesLevel = tree.getHeight() - 1 ; + auto BoxCellWidth = tree .getBoxWidth()/ (2 << (leavesLevel-1)) ; + auto nbLeaves = 2 << (leavesLevel-1) ; + // nbLeaves = nbLeaves*nbLeaves*nbLeaves ; + std::cout << " H: " << leavesLevel << " nbCells: " << nbLeaves << " cellwidth " << BoxCellWidth << std::endl ; + auto bloc = nbLeaves / nbProc ; + startIdx = bloc*myId ; + endIdx = startIdx+bloc; + if(myId == nbProc-1){ + endIdx = nbLeaves ;} + particles.resize(endIdx-startIdx) ; + FReal_t val ; + if(particles.size() ==1) { + val = 0.1*( myId %2 ==0 ? 1 : -1 ); + } + else { + val = 0.1; + + } + std::cout << "procId " << myId << " MortonIdx: [ " << startIdx << " , " << endIdx <<" [" < centerOfBox = {0.5,0.5,0.5} ; + + std::vector particles; + + ////////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////// + + // + MortonIndex startIdx=0,endIdx=0; +//#ifdef SLICES + // Check a slice (Y,Z) for a given X = 0 + auto nb1DLeaves = 2 << (TreeHeight-2) ; + MortonIndex mPos = 0 , mPos1 = nb1DLeaves-1 ; +#ifdef ss + std::vector > Lines ; + + for( MortonIndex i = 0 ; i < nb1DLeaves; ++i){ + // first slice + auto tmpPos = mPos; + std::array tmp = {mPos,i}; + // if( i==0 || i== nb1DLeaves-1) + Lines.emplace_back(tmp); + // second slice + tmp[0] = ++tmpPos; +// if( i==0|| i== nb1DLeaves-1 + Lines.emplace_back(tmp); +// tmp[0] = ++tmpPos ; +// // if( i==0 || i== nb1DLeaves-1) +// Lines.emplace_back(tmp); + // tmp[0] = mPos1; + // Lines.emplace_back(tmp); + + } +#else + std::vector > Lines ={{1,1}} ; + +#endif + for ( auto v : Lines){ + std::cout << " {" << v[0] << ", " << v[1]<< "}" ; + } + std::cout << std::endl; + setOneParticleInInCellsLinesZ(Lines, boxWidth,TreeHeight, particles,startIdx, endIdx); +//#endif +// setOneParticleInInCellsLineZ(boxWidth,TreeHeight, particles,startIdx, endIdx); + // setOneParticleInInCellsLineY(tree, particles,startIdx); + // setOneParticleInInCellsLineY(tree, particles,startIdx); + //setOneParticleInInCellsLineX(tree, particles,startIdx, endIdx); +// setOneParticleInEachCell(boxWidth,TreeHeight, particles,startIdx, endIdx); + FSize NbPoints = particles.size(); + + if(FParameters::existParameter(argc, argv, "-zeromean")){ + // set the total charge to zero in order to have a solution + FReal_t charge = 0 ; + for ( auto v : particles){ + charge += v.physicalValue ; + } + charge /= static_cast(NbPoints) ; + for ( auto v : particles){ + v.physicalValue -= charge ; + } + charge = 0.0 ; + for ( auto v : particles){ + charge += v.physicalValue ; + } + + std::cout << " GlobalCharge: "<< charge << std::endl; + } + + ////////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////// + // ----------------------------------------------------- + std::cout << "Write "<< NbPoints <<" particles to '" << filename << "'" << std::endl; + FFmaGenericWriter writer(filename); + + constexpr int nbData = 4 ; + writer.writeHeader(centerOfBox, boxWidth, NbPoints, sizeof(FReal_t),nbData); + + std::vector > finalPart(NbPoints); + for (int i=0 ; i < NbPoints ; ++i){ + finalPart[i].setPosition(particles[i].position) ; + finalPart[i].setPhysicalValue(particles[i].physicalValue) ; + } + writer.writeArrayOfReal(finalPart[0].getPtrFirstData(), nbData, NbPoints); + std::cout << "End of writing" <