Commit 89e17b15 authored by COULAUD Olivier's avatar COULAUD Olivier

Bugs with periodic conditions seem fixed. More tests should be done.

parent d303c1c9
......@@ -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<FReal> position; // Spatial position of the particle.
FReal physicalValue; // Physical value of the particle.
/* Returns the particle position. */
const FPoint<FReal>& 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<FReal> position; // Spatial position of the particle.
FReal physicalValue; // Physical value of the particle.
/* Returns the particle position. */
const FPoint<FReal>& 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<TestParticle> 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<TestParticle> 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<MortonIndex> 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 <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl;
std::cout << std::scientific;
std::cout.precision(15) ;
FReal TotalPhysicalValue=0.0 ;
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];
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<FSize>& 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 "<<posX[idxPart]<<" "<<posY[idxPart]<<" "<<posZ[idxPart]
<< " Forces: " << forcesX[idxPart] << " " << forcesY[idxPart] << " "<< forcesZ[idxPart] <<std::endl;
}
energy += potentials[idxPart]*physicalValues[idxPart] ;
TotalPhysicalValue += physicalValues[idxPart] ;
}
});
FReal gloEnergy = app.global().reduceSum(energy);
if(0 == app.global().processId()){
std::cout <<std::endl<<"aboveRoot: " << aboveTree << " Energy: "<< energy<<" TotalPhysicalValue: " << TotalPhysicalValue<< std::endl; std::cout <<std::endl << "Proc "<< app.global().processId() << " Energy: "<< gloEnergy <<std::endl;
std::cout <<std::endl <<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl<<std::endl;
}
}
// -----------------------------------------------------
if(FParameters::existParameter(argc, argv, FParameterDefinitions::OutputFile.options)){
algorithm->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 <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl;
std::cout << std::scientific;
std::cout.precision(15) ;
FReal locTotalPhysicalValue=0.0 ;
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];
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<FSize>& 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 "<<posX[idxPart]<<" "<<posY[idxPart]<<" "<<posZ[idxPart]
<< " Forces: " << forcesX[idxPart] << " " << forcesY[idxPart] << " "<< forcesZ[idxPart] <<std::endl;
}
energy += potentials[idxPart]*physicalValues[idxPart] ;
locTotalPhysicalValue += physicalValues[idxPart] ;
}
});
FReal gloEnergy = app.global().reduceSum(energy);
FReal TotalPhysicalValue = app.global().reduceSum(locTotalPhysicalValue);
if(0 == app.global().processId()){
std::cout <<std::endl<<"aboveRoot: " << aboveTree << " Energy: "<< energy<<" TotalPhysicalValue: " << TotalPhysicalValue<< std::endl; std::cout <<std::endl << "Proc "<< app.global().processId() << " Energy: "<< gloEnergy <<std::endl;
std::cout <<std::endl <<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl<<std::endl;
}
}
// -----------------------------------------------------
if(FParameters::existParameter(argc, argv, FParameterDefinitions::OutputFile.options)){
std::string name(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options, "output.fma"));
FMpiFmaGenericWriter<FReal> 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<FReal> paraWriter(name,app);
paraWriter.writeDistributionOfParticlesFromOctree(tree,loader.getNumberOfParticles(),localParticlesNumber,mortonLeafDistribution);
}
return 0;
return 0;
}
......@@ -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<MortonIndex> & mortonLeafDistribution) {} ;
#endif
/// [mortonLeafDistribution[2*p], mortonLeafDistribution[2*p+1] is the morton index shared by process p
virtual void getMortonLeafDistribution(std::vector<MortonIndex> & mortonLeafDistribution) {
mortonLeafDistribution.resize(2) ;
mortonLeafDistribution[0] = 0 ;
mortonLeafDistribution[1] = 8 << (nbLevelsInTree-1) ;
};
};
......
......@@ -2,6 +2,8 @@
#ifndef FFMMALGORITHMTHREADPROC_HPP
#define FFMMALGORITHMTHREADPROC_HPP
#define SCALFMM_DISTRIBUTED_ALGORITHM
#include <omp.h>
#include <algorithm>
//
......@@ -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<MortonIndex> & mortonLeafDistribution) final {
for (int p=0 ; p< comm.processCount() ; ++p ){
auto inter = getWorkingInterval(OctreeHeight-1, p );
void getMortonLeafDistribution(std::vector<MortonIndex> & 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;
}
......
......@@ -2,6 +2,8 @@
#ifndef FFMMALGORITHMTHREADPROCPPERIODIC_HPP
#define FFMMALGORITHMTHREADPROCPPERIODIC_HPP
#define SCALFMM_DISTRIBUTED_ALGORITHM
#include <algorithm>
#include <array>
#include <vector>
......@@ -28,6 +30,7 @@
#ifdef _OPENMP
#include <omp.h>
#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<MortonIndex> & 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<MortonIndex> & 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 " <<leaf->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<FTreeCoordinate,maxPeriodicNeighbors> 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] <<std::endl;
// std::cout << " Offset " << xoffset << " " << yoffset << " " << zoffset <<std::endl;
positionsX[idxPart] += xoffset;
positionsY[idxPart] += yoffset;
positionsZ[idxPart] += zoffset;
}
// std::cout << " Xn " << positionsX[idxPart] << " " << positionsY[idxPart] << " " << positionsZ[idxPart] <<std::endl;
}
// periodicOffsets[periodicNeighborsCounter] = offsets[idxNeig];
periodicNeighborPositions[periodicNeighborsCounter] = neighborPositions[idxNeig];
......@@ -1962,7 +1996,7 @@ protected:
}
else{ // set non periodic cells in neighbors and in neighborPositions
// nonPeriodicCounter versus idxNeig-periodicNeighborsCounter
std::cout << nonPeriodicCounter <<" "<< idxNeig-periodicNeighborsCounter <<std::endl;
// std::cout << nonPeriodicCounter <<" "<< idxNeig-periodicNeighborsCounter <<std::endl;
neighbors[idxNeig-periodicNeighborsCounter] = neighbors[idxNeig];
neighborPositions[idxNeig-periodicNeighborsCounter] = neighborPositions[idxNeig];
++nonPeriodicCounter ;
......@@ -2022,7 +2056,7 @@ protected:
const int limitem1 = limite -1 ;
FAssertLF(leafsNeedOtherData.getSize() < std::numeric_limits<int>::max());
const int nbLeafToProceed = int(leafsNeedOtherData.getSize());
// std::cout << "Treat communications " << nbLeafToProceed <<std::endl;
//#pragma omp for schedule(dynamic, userChunkSize)
#pragma omp single
for(int idxLeafs = 0 ; idxLeafs < nbLeafToProceed ; ++idxLeafs){
......@@ -2036,10 +2070,11 @@ protected:
const int nbNeigh = this->getNeighborsIndexesPeriodic(currentIter.coord,limite,
indexesNeighbors,indexArray,AllDirs);
// std::cout << currentIter.coord.getMortonIndex() << " Cells comming from another process " << nbNeigh <<std::endl;
std::array<ContainerClass*,26> periodicNeighbors{};
std::array<ContainerClass*,19> periodicNeighbors{};
int nbPeriodicCells=0 ;
// Pourquoi pas de decallage de boite ??,
// std::cout <<idProcess<< " LEAF " << currentIter.cell->getMortonIndex() << " " << currentIter.coord << " " << nbNeigh <<std::endl;
//
//
for(int idxNeigh = 0 ; idxNeigh < nbNeigh ; ++idxNeigh){
if(indexesNeighbors[idxNeigh] < (intervals[idProcess].leftIndex) || (intervals[idProcess].rightIndex) < indexesNeighbors[idxNeigh]){
ContainerClass*const hypotheticNeighbor = otherP2Ptree.getLeafSrc(indexesNeighbors[idxNeigh]);
......@@ -2069,7 +2104,7 @@ protected:
if(movePos){
/// NEW NEW NEW NEW NEW NEW
// Put periodic cell into another cell (copy data)
periodicNeighbors[nbPeriodicCells] = new ContainerClass(*(hypotheticNeighbor)) ;
periodicNeighbors[nbPeriodicCells] = new ContainerClass(*hypotheticNeighbor) ;
// newPos = pos + ofsset
FReal*const positionsX = periodicNeighbors[nbPeriodicCells]->getPositions()[0];
FReal*const positionsY = periodicNeighbors[nbPeriodicCells]->getPositions()[1];
......@@ -2080,26 +2115,37 @@ protected:
positionsZ[idxPart] += offSet[2];
}
#ifdef toto
std::cout << " MOVE CELLS "<<std::endl;
std::cout << " MOVE CELL "<<idxNeigh<<std::endl;
FReal*const positionsXX = hypotheticNeighbor->getPositions()[0];
FReal*const positionsYX = hypotheticNeighbor->getPositions()[1];
FReal*const positionsZX = hypotheticNeighbor->getPositions()[2];
auto idxPart = 0 ;
std::cout <<neighborPositions[idxNeigh] << " Shift " << positionsX[idxPart] << " "<< positionsY[idxPart] << " "<< positionsZ[idxPart] <<std::endl;
std::cout <<neighborPositions[idxNeigh] << " ORI " << positionsXX[idxPart] << " "<< positionsYX[idxPart] << " "<< positionsZX[idxPart] <<std::endl;
std::cout <<neighborPositions[idxNeigh] << " OFSET " << offSet[0] << " "<< offSet[1] << " "<< offSet[2] <<std::endl;
std::cout <<idProcess << " ORI " << positionsXX[idxPart] << " "<< positionsYX[idxPart] << " "<< positionsZX[idxPart] <<std::endl;
std::cout <<idProcess << " OFSET " << offSet[0] << " "<< offSet[1] << " "<< offSet[2] <<std::endl;
std::cout <<idProcess << " Shift " << positionsX[idxPart] << " "<< positionsY[idxPart] << " "<< positionsZ[idxPart] <<std::endl;
#endif
++nbPeriodicCells;
neighbors[counter] = periodicNeighbors[nbPeriodicCells];
// std::cout << " hypotheticNeighbor " << hypotheticNeighbor << " " <<neighbors[counter] << " " << periodicNeighbors[nbPeriodicCells] <<std::endl;
++nbPeriodicCells;
}
else {
// std::cout << " NO MOVE CELL "<<idxNeigh<<std::endl;
}
++counter;
}
}
}
// std::cout << " counter: "<< counter <<std::endl << " -->";
if(counter){ // neighborPositions doesn't use)
if(counter> 0){ // neighborPositions doesn't use)
// std::cout << " counter: "<< counter <<" nbPeriodicCells " << nbPeriodicCells<<std::endl << " --> P2PRemote" <<std::endl;;
// for (int i = 0 ; i < counter ; ++i){
// std::cout << " " << neighbors[i] ;
// }
// std::cout <<std::endl;
myThreadkernels.P2PRemote( currentIter.coord, currentIter.targets,
currentIter.sources, neighbors, neighborPositions, counter);
nullptr /*currentIter.sources*/, neighbors, neighborPositions, counter);
if(nbPeriodicCells >0){
//to Do
for(int i=0 ; i < nbPeriodicCells ; ++i){
......
......@@ -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
......
// See LICENCE file at project root
// ==== CMAKE =====
// ================
//
// make testPeriodicAlgorithm
// mpirun --oversubscribe -np 4 --tag-output Tests/Release/testPeriodicAlgorithm -h 3 -sh 1
#include <iostream>
#include <vector>
#include <array>
#include <algorithm>
#include <stdexcept>
//#include <cstdio>
#include <cstdlib>
//
// Specific part
//
#include <string>
//
// 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<FReal_t> position; // Spatial position of the particle.
FReal_t physicalValue; // Physical value of the particle.
/* Returns the particle position. */
const FPoint<FReal_t>& getPosition(){
return position;
}
};
void setOneParticleInEachCell(const FReal_t &BoxWidth , const int &Height, std::vector<TestParticle> &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 <<" [" <<std::endl ;
for(MortonIndex idx = startIdx, pos=0 ; idx < endIdx ; ++idx,++pos ) {
FTreeCoordinate coor ;
coor.setPositionFromMorton(idx) ;
particles[pos].position = { (0.5+coor[0])*BoxCellWidth,(0.5+coor[1])*BoxCellWidth,(0.5+coor[2])*BoxCellWidth } ;
particles[pos].index = idx ;
particles[pos].physicalValue =val ;
val *= -1.0 ;
std::cout << idx << " " << particles[pos].index << " " << particles[pos].position << " " << particles[pos].physicalValue << std::endl;
}
}
void setOneParticleInInCellsLinesZ(const std::vector<std::array<MortonIndex,2> >& Lines, const FReal_t &BoxWidth , const int &Height, std::vector<TestParticle> &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<MortonIndex> cells(nb1DLeaves*nbLines) ;
// fill cells