Skip to content

GitLab

S
ScalFMM
Commit b0830fbe authored by BRAMAS Berenger's avatar BRAMAS Berenger
Browse files
Options

Merge branch 'master' of git+ssh://scm.gforge.inria.fr//gitroot//scalfmm/scalfmm

parents 0b18cdb9 b21d1352
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 966 additions and 522 deletions
Doc/ParallelDetails.pdf
View file @ b0830fbe
No preview for this file type
Src/Containers/FMpiBufferWriter.hpp
View file @ b0830fbe
......@@ -43,7 +43,7 @@ class FMpiBufferWriter : public FAbstractBufferWriter {
public:
/** Constructor with a default arrayCapacity of 512 bytes */
FMpiBufferWriter(const MPI_Comm inComm, const int inCapacity = 512):
FMpiBufferWriter(const MPI_Comm inComm, const int inCapacity = 1024):
mpiComm(inComm),
arrayCapacity(inCapacity),
array(new char[inCapacity]),
......@@ -77,6 +77,7 @@ public:
/** Write data by packing cpy */
template <class ClassType>
void write(const ClassType& object){
//printf("Space need in the write : %d, index set on %d \n",sizeof(ClassType),currentIndex);
assertRemainingSpace(sizeof(ClassType));
MPI_Pack(const_cast<ClassType*>(&object), 1, FMpi::GetType(object), array.get(), arrayCapacity, &currentIndex, mpiComm);
}
......@@ -86,6 +87,7 @@ public:
*/
template <class ClassType>
void write(const ClassType&& object){
// printf("Space need in the write : %d \n",sizeof(ClassType));
assertRemainingSpace(sizeof(ClassType));
MPI_Pack(const_cast<ClassType*>(&object), 1, FMpi::GetType(object), array.get(), arrayCapacity, &currentIndex, mpiComm);
}
......@@ -93,7 +95,7 @@ public:
/** Write back, position + sizeof(object) has to be < size */
template <class ClassType>
void writeAt(const int position, const ClassType& object){
if(position + sizeof(ClassType) > currentIndex){
if(position + (int) sizeof(ClassType) > currentIndex){
printf("Not enought space\n");
exit(0);
}
......@@ -106,6 +108,7 @@ public:
*/
template <class ClassType>
void write(const ClassType* const objects, const int inSize){
// printf("Space need in the write : %d, index set on %d, and capacity is %d \n",sizeof(ClassType)*inSize,currentIndex,arrayCapacity);
assertRemainingSpace(sizeof(ClassType) * inSize);
MPI_Pack( const_cast<ClassType*>(objects), inSize, FMpi::GetType(*objects), array.get(), arrayCapacity, &currentIndex, mpiComm);
}
......
Src/Core/FFmmAlgorithmThreadProc.hpp
View file @ b0830fbe
......@@ -246,10 +246,8 @@ private:
}
FLOG(computationCounter.tac());
FLOG( FLog::Controller << "\tFinished (@Bottom Pass (P2M) = " << counterTime.tacAndElapsed() << "s)\n" );
FLOG( FLog::Controller << "\t\t Computation : " << computationCounter.elapsed() << " s\n" );
}
......@@ -281,11 +279,11 @@ private:
MPI_Status status[14];
// Maximum data per message is:
FBufferWriter sendBuffer/*(comm.getComm())*/;
FMpiBufferWriter sendBuffer(comm.getComm(),7*MaxSizePerCell);
const int recvBufferOffset = (8 * MaxSizePerCell + 1);
FBufferReader recvBuffer(/*comm.getComm(),*/ nbProcess*recvBufferOffset);
FMpiBufferReader recvBuffer(comm.getComm(), nbProcess*recvBufferOffset);
CellClass recvBufferCells[8];
int firstProcThatSend = idProcess + 1;
// for each levels
......@@ -318,28 +316,25 @@ private:
FLOG(prepareCounter.tic());
if(idProcess != 0
&& (getWorkingInterval((idxLevel+1), idProcess).min >>3) <= (getWorkingInterval((idxLevel+1), idProcess - 1).max >>3)){
char state = 0;
sendBuffer.write(state);
printf("Index position after writing first state(%d) :: %d\n",state,sendBuffer.getSize());
const CellClass* const* const child = iterArray[cellsToSend].getCurrentChild();
for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
if( child[idxChild] && getWorkingInterval((idxLevel+1), idProcess).min <= child[idxChild]->getMortonIndex() ){
child[idxChild]->serializeUp(sendBuffer);
state = char(state | (0x1 << idxChild));
}
}
printf("state just before Send :: %d \n",state);
sendBuffer.writeAt(0,state);
printf("what is it in sendBuffer :: %d\n", sendBuffer.data()[0]);
while( sendToProc && iterArray[cellsToSend].getCurrentGlobalIndex() == getWorkingInterval(idxLevel , sendToProc - 1).max){
--sendToProc;
}
MPI_Isend(sendBuffer.data(), sendBuffer.getSize(), MPI_BYTE, sendToProc, FMpi::TagFmmM2M, comm.getComm(), &requests[iterRequests++]);
MPI_Isend(sendBuffer.data(), sendBuffer.getSize(), MPI_PACKED, sendToProc,
FMpi::TagFmmM2M, comm.getComm(), &requests[iterRequests++]);
}
// We may need to receive something
......@@ -351,16 +346,15 @@ private:
&& (getWorkingInterval((idxLevel+1), firstProcThatSend).max) < (getWorkingInterval((idxLevel+1), idProcess).max)){
// Second condition :: while firstProcThatSend max morton index is < to myself max interval
++firstProcThatSend;
printf("\n \n PLOP \n \n");
}
if(firstProcThatSend < nbProcess &&
(getWorkingInterval((idxLevel+1), firstProcThatSend).min >>3) <= (getWorkingInterval((idxLevel+1) , idProcess).max>>3) ){
(getWorkingInterval((idxLevel+1), firstProcThatSend).min >>3) == (getWorkingInterval((idxLevel+1) , idProcess).max>>3) ){
endProcThatSend = firstProcThatSend;
while( endProcThatSend < nbProcess &&
(getWorkingInterval((idxLevel+1) ,endProcThatSend).min >>3) <= (getWorkingInterval((idxLevel+1) , idProcess).max>>3)){
(getWorkingInterval((idxLevel+1) ,endProcThatSend).min >>3) == (getWorkingInterval((idxLevel+1) , idProcess).max>>3)){
++endProcThatSend;
}
......@@ -369,7 +363,7 @@ private:
hasToReceive = true;
for(int idxProc = firstProcThatSend ; idxProc < endProcThatSend ; ++idxProc ){
MPI_Irecv(&recvBuffer.data()[idxProc * recvBufferOffset], recvBufferOffset, MPI_BYTE,
MPI_Irecv(&recvBuffer.data()[idxProc * recvBufferOffset], recvBufferOffset, MPI_PACKED,
idxProc, FMpi::TagFmmM2M, comm.getComm(), &requests[iterRequests++]);
}
}
......@@ -406,7 +400,7 @@ private:
for(int idxProc = firstProcThatSend ; idxProc < endProcThatSend ; ++idxProc){
recvBuffer.seek(idxProc * recvBufferOffset);
int state = int(recvBuffer.getValue<char>());
printf("state read in recv :: %d \n",state);
int position = 0;
while( state && position < 8){
while(!(state & 0x1)){
......@@ -583,35 +577,37 @@ private:
const int SizeOfCellToSend = sizeof(MortonIndex) + sizeof(int) + MaxSizePerCell;
FBufferWriter**const sendBuffer = new FBufferWriter*[nbProcess * OctreeHeight];
memset(sendBuffer, 0, sizeof(FBufferWriter*) * nbProcess * OctreeHeight);
FMpiBufferWriter**const sendBuffer = new FMpiBufferWriter*[nbProcess * OctreeHeight];
memset(sendBuffer, 0, sizeof(FMpiBufferWriter*) * nbProcess * OctreeHeight);
FBufferReader**const recvBuffer = new FBufferReader*[nbProcess * OctreeHeight];
memset(recvBuffer, 0, sizeof(FBufferReader*) * nbProcess * OctreeHeight);
FMpiBufferReader**const recvBuffer = new FMpiBufferReader*[nbProcess * OctreeHeight];
memset(recvBuffer, 0, sizeof(FMpiBufferReader*) * nbProcess * OctreeHeight);
for(int idxLevel = 2 ; idxLevel < OctreeHeight ; ++idxLevel ){
for(int idxProc = 0 ; idxProc < nbProcess ; ++idxProc){
const int toSendAtProcAtLevel = indexToSend[idxLevel * nbProcess + idxProc];
if(toSendAtProcAtLevel != 0){
sendBuffer[idxLevel * nbProcess + idxProc] = new FBufferWriter(toSendAtProcAtLevel * SizeOfCellToSend);
sendBuffer[idxLevel * nbProcess + idxProc] = new FMpiBufferWriter(comm.getComm(),toSendAtProcAtLevel * SizeOfCellToSend);
for(int idxLeaf = 0 ; idxLeaf < toSendAtProcAtLevel; ++idxLeaf){
const MortonIndex cellIndex = toSend[idxLevel * nbProcess + idxProc][idxLeaf].getCurrentGlobalIndex();
sendBuffer[idxLevel * nbProcess + idxProc]->write(cellIndex);
toSend[idxLevel * nbProcess + idxProc][idxLeaf].getCurrentCell()->serializeUp(*sendBuffer[idxLevel * nbProcess + idxProc]);
}
FMpi::MpiAssert( MPI_Isend( sendBuffer[idxLevel * nbProcess + idxProc]->data(), sendBuffer[idxLevel * nbProcess + idxProc]->getSize()
, MPI_BYTE , idxProc, FMpi::TagLast + idxLevel, comm.getComm(), &requests[iterRequest++]) , __LINE__ );
FMpi::MpiAssert( MPI_Isend( sendBuffer[idxLevel * nbProcess + idxProc]->data(),
sendBuffer[idxLevel * nbProcess + idxProc]->getSize(),MPI_PACKED, idxProc,
FMpi::TagLast + idxLevel, comm.getComm(), &requests[iterRequest++]) , __LINE__ );
}
const int toReceiveFromProcAtLevel = globalReceiveMap[(idxProc * nbProcess * OctreeHeight) + idxLevel * nbProcess + idProcess];
if(toReceiveFromProcAtLevel){
recvBuffer[idxLevel * nbProcess + idxProc] = new FBufferReader(toReceiveFromProcAtLevel * SizeOfCellToSend);
recvBuffer[idxLevel * nbProcess + idxProc] = new FMpiBufferReader(comm.getComm(),toReceiveFromProcAtLevel * SizeOfCellToSend);
FMpi::MpiAssert( MPI_Irecv(recvBuffer[idxLevel * nbProcess + idxProc]->data(), recvBuffer[idxLevel * nbProcess + idxProc]->getSize(), MPI_BYTE,
idxProc, FMpi::TagLast + idxLevel, comm.getComm(), &requests[iterRequest++]) , __LINE__ );
FMpi::MpiAssert( MPI_Irecv(recvBuffer[idxLevel * nbProcess + idxProc]->data(),
recvBuffer[idxLevel * nbProcess + idxProc]->getCapacity(), MPI_PACKED,idxProc,
FMpi::TagLast + idxLevel, comm.getComm(), &requests[iterRequest++]) , __LINE__ );
}
}
}
......@@ -699,12 +695,12 @@ private:
const int toReceiveFromProcAtLevel = globalReceiveMap[(idxProc * nbProcess * OctreeHeight) + idxLevel * nbProcess + idProcess];
for(int idxCell = 0 ; idxCell < toReceiveFromProcAtLevel ; ++idxCell){
const MortonIndex cellIndex = recvBuffer[idxLevel * nbProcess + idxProc]->FBufferReader::getValue<MortonIndex>();
const MortonIndex cellIndex = recvBuffer[idxLevel * nbProcess + idxProc]->FMpiBufferReader::getValue<MortonIndex>();
CellClass* const newCell = new CellClass;
newCell->setMortonIndex(cellIndex);
newCell->deserializeUp(*recvBuffer[idxLevel * nbProcess + idxProc]);
tempTree.insertCell(cellIndex, idxLevel, newCell);
}
}
......@@ -819,8 +815,8 @@ private:
const int heightMinusOne = OctreeHeight - 1;
FBufferWriter sendBuffer;
FBufferReader recvBuffer(MaxSizePerCell);
FMpiBufferWriter sendBuffer(comm.getComm());
FMpiBufferReader recvBuffer(comm.getComm(),MaxSizePerCell);
// for each levels exepted leaf level
for(int idxLevel = 2 ; idxLevel < heightMinusOne ; ++idxLevel ){
......@@ -858,7 +854,7 @@ private:
needToRecv = true;
MPI_Irecv( recvBuffer.data(), recvBuffer.getSize(), MPI_BYTE, MPI_ANY_SOURCE,
MPI_Irecv( recvBuffer.data(), recvBuffer.getCapacity(), MPI_PACKED, MPI_ANY_SOURCE,
FMpi::TagFmmL2L, comm.getComm(), &requests[iterRequests++]);
}
......@@ -881,7 +877,7 @@ private:
for(int idxProc = firstProcThatRecv ; idxProc < endProcThatRecv ; ++idxProc ){
MPI_Isend(sendBuffer.data(), sendBuffer.getSize(), MPI_BYTE, idxProc,
MPI_Isend(sendBuffer.data(), sendBuffer.getSize(), MPI_PACKED, idxProc,
FMpi::TagFmmL2L, comm.getComm(), &requests[iterRequests++]);
}
......@@ -961,11 +957,11 @@ private:
int iterRequest = 0;
int nbMessagesToRecv = 0;
FBufferWriter**const sendBuffer = new FBufferWriter*[nbProcess];
memset(sendBuffer, 0, sizeof(FBufferWriter*) * nbProcess);
FMpiBufferWriter**const sendBuffer = new FMpiBufferWriter*[nbProcess];
memset(sendBuffer, 0, sizeof(FMpiBufferWriter*) * nbProcess);
FBufferReader**const recvBuffer = new FBufferReader*[nbProcess];
memset(recvBuffer, 0, sizeof(FBufferReader*) * nbProcess);
FMpiBufferReader**const recvBuffer = new FMpiBufferReader*[nbProcess];
memset(recvBuffer, 0, sizeof(FMpiBufferReader*) * nbProcess);
/* This a nbProcess x nbProcess matrix of integer
* let U and V be id of processes :
......@@ -1074,8 +1070,8 @@ private:
for(int idxProc = 0 ; idxProc < nbProcess ; ++idxProc){
if(globalReceiveMap[idxProc * nbProcess + idProcess]){ //if idxProc has sth for me.
//allocate buffer of right size
recvBuffer[idxProc] = new FBufferReader(globalReceiveMap[idxProc * nbProcess + idProcess]);
FMpi::MpiAssert( MPI_Irecv(recvBuffer[idxProc]->data(), recvBuffer[idxProc]->getSize(), MPI_BYTE,
recvBuffer[idxProc] = new FMpiBufferReader(comm.getComm(),globalReceiveMap[idxProc * nbProcess + idProcess]);
FMpi::MpiAssert( MPI_Irecv(recvBuffer[idxProc]->data(), recvBuffer[idxProc]->getCapacity(), MPI_PACKED,
idxProc, FMpi::TagFmmP2P, comm.getComm(), &requests[iterRequest++]) , __LINE__ );
}
}
......@@ -1084,8 +1080,7 @@ private:
// Prepare send
for(int idxProc = 0 ; idxProc < nbProcess ; ++idxProc){
if(toSend[idxProc].getSize() != 0){
//sendBuffer[idxProc] = new FBufferWriter(partsToSend[idxProc]); //Could be read out of globalReceiveMap
sendBuffer[idxProc] = new FBufferWriter(globalReceiveMap[idProcess*nbProcess+idxProc]);
sendBuffer[idxProc] = new FMpiBufferWriter(comm.getComm(),globalReceiveMap[idProcess*nbProcess+idxProc]);
// << is equivalent to write().
(*sendBuffer[idxProc]) << toSend[idxProc].getSize();
for(int idxLeaf = 0 ; idxLeaf < toSend[idxProc].getSize() ; ++idxLeaf){
......@@ -1095,7 +1090,7 @@ private:
//TEST BERENGER
//if(sendBuffer[idxProc]->getSize() != partsToSend[idxProc]){
FMpi::MpiAssert( MPI_Isend( sendBuffer[idxProc]->data(), sendBuffer[idxProc]->getSize() , MPI_BYTE ,
FMpi::MpiAssert( MPI_Isend( sendBuffer[idxProc]->data(), sendBuffer[idxProc]->getSize() , MPI_PACKED ,
idxProc, FMpi::TagFmmP2P, comm.getComm(), &requests[iterRequest++]) , __LINE__ );
}
......
Src/Utils/FMpi.hpp
View file @ b0830fbe
......@@ -23,6 +23,8 @@
#include "FNoCopyable.hpp"
#include "FMath.hpp"
//Need that for converting datas
#include "FComplexe.hpp"
/////////////////////////////////////////////////////////////////////////////////////////
......@@ -280,6 +282,12 @@ public:
return MPI_CHAR;
}
static const MPI_Datatype GetType(const FComplexe& a){
MPI_Datatype FMpiComplexe;
MPI_Type_contiguous(2, GetType(a.getReal()) , &FMpiComplexe);
return FMpiComplexe;
}
////////////////////////////////////////////////////////////
// Mpi interface functions
////////////////////////////////////////////////////////////
......
Tests/Kernels/testSphericalProcAlgorithm.cpp
View file @ b0830fbe
......@@ -29,6 +29,9 @@
#include "../../Src/Kernels/Spherical/FSphericalKernel.hpp"
#include "../../Src/Kernels/Spherical/FSphericalCell.hpp"
#include "../../Src/Kernels/Rotation/FRotationKernel.hpp"
#include "../../Src/Kernels/Rotation/FRotationCell.hpp"
#include "../../Src/Core/FFmmAlgorithmThreadProc.hpp"
#include "../../Src/Core/FFmmAlgorithmThread.hpp"
......@@ -73,7 +76,7 @@ bool isEqualPole(const CellClass& me, const CellClass& other, FReal*const cumul)
}
/** To compare data */
bool isEqualLocal(const FSphericalCell& me, const FSphericalCell& other, FReal*const cumul){
bool isEqualLocal(const FSphericalCell& me, const FSphericalCell& other,FReal*const cumul){
FMath::FAccurater accurate;
for(int idx = 0; idx < FSphericalCell::GetLocalSize(); ++idx){
accurate.add(me.getLocal()[idx].getImag(),other.getLocal()[idx].getImag());
......@@ -210,202 +213,206 @@ void ValidateFMMAlgoProc(OctreeClass* const badTree,
// Simply create particles and try the kernels
int main(int argc, char ** argv){
typedef FSphericalCell CellClass;
typedef FP2PParticleContainer ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass;
typedef FOctree< CellClass, ContainerClass , LeafClass > OctreeClass;
typedef FSphericalKernel< CellClass, ContainerClass > KernelClass;
typedef FFmmAlgorithmThreadProc<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
typedef FFmmAlgorithmThread<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClassNoProc;
///////////////////////What we do/////////////////////////////
std::cout << ">> This executable has to be used to test Spherical algorithm.\n";
//////////////////////////////////////////////////////////////
FMpi app( argc, argv);
const int DevP = FParameters::getValue(argc,argv,"-p", 8);
const int NbLevels = FParameters::getValue(argc,argv,"-h", 5);
const int SizeSubLevels = FParameters::getValue(argc,argv,"-sh", 3);
FTic counter;
const char* const defaultFilename = (sizeof(FReal) == sizeof(float))?
"../../Data/test20k.bin.fma.single":
"../../Data/test20k.bin.fma.double";
const char* const filename = FParameters::getStr(argc,argv,"-f", defaultFilename);
std::cout << "Opening : " << filename << "\n";
FMpiFmaLoader loader(filename, app.global());
if(!loader.isOpen()){
std::cout << "Loader Error, " << filename << " is missing\n";
return 1;
}
typedef FSphericalCell CellClass;
typedef FP2PParticleContainer ContainerClass;
// -----------------------------------------------------
CellClass::Init(DevP);
OctreeClass tree(NbLevels, SizeSubLevels,loader.getBoxWidth(),loader.getCenterOfBox());
typedef FSimpleLeaf< ContainerClass > LeafClass;
typedef FOctree< CellClass, ContainerClass , LeafClass > OctreeClass;
typedef FSphericalKernel< CellClass, ContainerClass > KernelClass;
// -----------------------------------------------------
typedef FFmmAlgorithmThreadProc<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
typedef FFmmAlgorithmThread<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClassNoProc;
std::cout << "Creating & Inserting " << loader.getNumberOfParticles() << " particles ..." << std::endl;
std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
counter.tic();
if( app.global().processCount() != 1){
//////////////////////////////////////////////////////////////////////////////////
// Build tree from mpi loader
//////////////////////////////////////////////////////////////////////////////////
std::cout << "Build Tree ..." << std::endl;
counter.tic();
struct TestParticle{
FPoint position;
FReal physicalValue;
const FPoint& getPosition(){
return position;
}
};
///////////////////////What we do/////////////////////////////
std::cout << ">> This executable has to be used to test Spherical algorithm.\n";
//////////////////////////////////////////////////////////////
TestParticle* particles = new TestParticle[loader.getNumberOfParticles()];
memset(particles, 0, sizeof(TestParticle) * loader.getNumberOfParticles());
FMpi app( argc, argv);
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
loader.fillParticle(&particles[idxPart].position,&particles[idxPart].physicalValue);
}
const int DevP = FParameters::getValue(argc,argv,"-p", 8);
const int NbLevels = FParameters::getValue(argc,argv,"-h", 5);
const int SizeSubLevels = FParameters::getValue(argc,argv,"-sh", 3);
FTic counter;
const char* const defaultFilename = (sizeof(FReal) == sizeof(float))?
"../Data/test20k.bin.fma.single":
"../Data/test20k.bin.fma.double";
const char* const filename = FParameters::getStr(argc,argv,"-f", defaultFilename);
FVector<TestParticle> finalParticles;
FMpiTreeBuilder< TestParticle >::ArrayToTree(app.global(), particles, loader.getNumberOfParticles(),
tree.getBoxCenter(),
tree.getBoxWidth(),
tree.getHeight(), &finalParticles);
std::cout << "Opening : " << filename << "\n";
for(int idx = 0 ; idx < finalParticles.getSize(); ++idx){
tree.insert(finalParticles[idx].position);
}
FMpiFmaLoader loader(filename, app.global());
if(!loader.isOpen()){
std::cout << "Loader Error, " << filename << " is missing\n";
return 1;
}
CellClass::Init(DevP);
OctreeClass tree(NbLevels, SizeSubLevels,loader.getBoxWidth(),loader.getCenterOfBox());
delete[] particles;
// -----------------------------------------------------
counter.tac();
std::cout << "Done " << "(" << counter.elapsed() << "s)." << std::endl;
std::cout << "Creating & Inserting " << loader.getNumberOfParticles() << " particles ..." << std::endl;
std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
counter.tic();
if( app.global().processCount() != 1){
//////////////////////////////////////////////////////////////////////////////////
// Build tree from mpi loader
//////////////////////////////////////////////////////////////////////////////////
std::cout << "Build Tree ..." << std::endl;
counter.tic();
//////////////////////////////////////////////////////////////////////////////////
struct TestParticle{
FPoint position;
FReal physicalValue;
const FPoint& getPosition(){
return position;
}
};
TestParticle* particles = new TestParticle[loader.getNumberOfParticles()];
memset(particles, 0, sizeof(TestParticle) * loader.getNumberOfParticles());
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
loader.fillParticle(&particles[idxPart].position,&particles[idxPart].physicalValue);
}
else{
FPoint position;
FReal physicalValue;
for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
loader.fillParticle(&position,&physicalValue);
tree.insert(position, physicalValue);
}
FVector<TestParticle> finalParticles;
FMpiTreeBuilder< TestParticle >::ArrayToTree(app.global(), particles, loader.getNumberOfParticles(),
tree.getBoxCenter(),
tree.getBoxWidth(),
tree.getHeight(), &finalParticles);
for(int idx = 0 ; idx < finalParticles.getSize(); ++idx){
tree.insert(finalParticles[idx].position,finalParticles[idx].physicalValue);
}
delete[] particles;
counter.tac();
std::cout << "Done " << "(@Creating and Inserting Particles = " << counter.elapsed() << "s)." << std::endl;
std::cout << "Done " << "(" << counter.elapsed() << "s)." << std::endl;
//////////////////////////////////////////////////////////////////////////////////
}
else{
FPoint position;
FReal physicalValue;
for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
loader.fillParticle(&position,&physicalValue);
tree.insert(position, physicalValue);
}
}
// -----------------------------------------------------
std::cout << "Create kernel..." << std::endl;
counter.tac();
std::cout << "Done " << "(@Creating and Inserting Particles = " << counter.elapsed() << "s)." << std::endl;
KernelClass kernels(DevP, NbLevels,loader.getBoxWidth(), loader.getCenterOfBox());
// -----------------------------------------------------
std::cout << "Create kernel..." << std::endl;
std::cout << "Done " << " in " << counter.elapsed() << "s)." << std::endl;
KernelClass kernels(DevP, NbLevels,loader.getBoxWidth(), loader.getCenterOfBox());
// -----------------------------------------------------
std::cout << "Done " << " in " << counter.elapsed() << "s)." << std::endl;
std::cout << "Working on particles ..." << std::endl;
// -----------------------------------------------------
FmmClass algo(app.global(),&tree,&kernels);
std::cout << "Working on particles ..." << std::endl;
counter.tic();
algo.execute();
counter.tac();
FmmClass algo(app.global(),&tree,&kernels);
std::cout << "Done " << "(@Algorithm = " << counter.elapsed() << "s)." << std::endl;
counter.tic();
algo.execute();
counter.tac();
{ // get sum forces&potential
FTRACE( FTrace::FFunction functionTrace(__FUNCTION__, "Sum Result" , __FILE__ , __LINE__) );
std::cout << "Done " << "(@Algorithm = " << counter.elapsed() << "s)." << std::endl;
FReal potential = 0;
FReal fx = 0.0, fy = 0.0, fz = 0.0;
{ // get sum forces&potential
FTRACE( FTrace::FFunction functionTrace(__FUNCTION__, "Sum Result" , __FILE__ , __LINE__) );
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();
const int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
FReal potential = 0;
FReal fx = 0.0, fy = 0.0, fz = 0.0;
for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
potential += potentials[idxPart];
fx += forcesX[idxPart];
fy += forcesY[idxPart];
fz += forcesZ[idxPart];
}
});
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();
const int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
std::cout << "My potential is " << potential << std::endl;
for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
potential += potentials[idxPart];
fx += forcesX[idxPart];
fy += forcesY[idxPart];
fz += forcesZ[idxPart];
}
});
potential = app.global().reduceSum(potential);
fx = app.global().reduceSum(fx);
fy = app.global().reduceSum(fy);
fz = app.global().reduceSum(fz);
std::cout << "My potential is " << potential << std::endl;
potential = app.global().reduceSum(potential);
fx = app.global().reduceSum(fx);
fy = app.global().reduceSum(fy);
fz = app.global().reduceSum(fz);
if(app.global().processId() == 0){
std::cout << "Foces Sum x = " << fx << " y = " << fy << " z = " << fz << std::endl;
std::cout << "Potential Sum = " << potential << std::endl;
}
if(app.global().processId() == 0){
std::cout << "Foces Sum x = " << fx << " y = " << fy << " z = " << fz << std::endl;
std::cout << "Potential Sum = " << potential << std::endl;
}
}
#ifdef VALIDATE_FMM
{
OctreeClass treeValide(NbLevels, SizeSubLevels,loader.getBoxWidth(),loader.getCenterOfBox());
{
OctreeClass treeValide(NbLevels, SizeSubLevels,loader.getBoxWidth(),loader.getCenterOfBox());
{
FFmaBinLoader loaderSeq(filename);
FPoint position;
FReal physicalValue;
for(FSize idxPart = 0 ; idxPart < loaderSeq.getNumberOfParticles() ; ++idxPart){
loaderSeq.fillParticle(&position,&physicalValue);
treeValide.insert(position,physicalValue);