Commit c22ccd0a authored by BRAMAS Berenger's avatar BRAMAS Berenger

move to new periodic

parent b1d39835
......@@ -9,7 +9,7 @@ project(Lib_scalfmm)
set(LIBRARY_OUTPUT_PATH ../lib/${CMAKE_BUILD_TYPE})
ADD_DEFINITIONS(${ScaLFMM_CXX_FLAGS})
MESSAGE( ${ScaLFMM_CXX_FLAGS})
# Searching all cpp file
file(
GLOB_RECURSE
......
......@@ -29,111 +29,103 @@
* finally use data pointer as you like
*/
class FMpiBufferReader : public FAbstractBufferReader {
const MPI_Comm comm; //< Communicator needed by MPI_Pack functions
const int arrayCapacity; //< Allocated space
std::unique_ptr<char[]> array; //< Allocated Array
int currentIndex;
/** Test and exit if not enought space */
void assertRemainingSpace(const size_t requestedSpace) const {
if(int(currentIndex + requestedSpace) > arrayCapacity){
printf("Error FMpiBufferReader has not enough space\n");
exit(0);
}
}
const MPI_Comm comm; //< Communicator needed by MPI_Pack functions
const int arrayCapacity; //< Allocated space
std::unique_ptr<char[]> array; //< Allocated Array
int currentIndex;
public :
/*Constructor with a default arrayCapacity of 512 bytes */
FMpiBufferReader(const MPI_Comm inComm, const int inCapacity = 512):
comm(inComm),
arrayCapacity(inCapacity),
array(new char[inCapacity]),
currentIndex(0)
{}
/** Destructor
/*Constructor with a default arrayCapacity of 512 bytes */
FMpiBufferReader(const MPI_Comm inComm, const int inCapacity = 512):
comm(inComm),
arrayCapacity(inCapacity),
array(new char[inCapacity]),
currentIndex(0)
{}
/** Destructor
*/
virtual ~FMpiBufferReader(){
}
/** Get allocated memory pointer */
char* data(){
return array.get();
}
/** Get allocated memory pointer */
const char* data() const {
return array.get();
}
/** get the filled space */
int getSize() const{
return currentIndex;
}
/** Size of the memory initialized */
int getCapacity() const{
return arrayCapacity;
}
/** Move the read index to a position */
void seek(const int inIndex){
if(inIndex > arrayCapacity){
printf("FMpiBufferReader :: Aborting :: Can't move index because buffer isn't long enough");
exit(0);
virtual ~FMpiBufferReader(){
}
/** Get allocated memory pointer */
char* data(){
return array.get();
}
/** Get allocated memory pointer */
const char* data() const {
return array.get();
}
/** get the filled space */
int getSize() const{
return currentIndex;
}
/** Size of the memory initialized */
int getCapacity() const{
return arrayCapacity;
}
/** Move the read index to a position */
void seek(const int inIndex){
if(inIndex > arrayCapacity){
printf("FMpiBufferReader :: Aborting :: Can't move index because buffer isn't long enough");
exit(0);
}
else{
currentIndex = inIndex;
}
}
else{
currentIndex = inIndex;
/** Get the read position */
int tell() const {
return currentIndex;
}
/** Get a value with memory cast */
template <class ClassType>
ClassType getValue(){
ClassType value;
int previousIndex = currentIndex;
seek(int(sizeof(value) + previousIndex));
MPI_Unpack(array.get(),arrayCapacity,&previousIndex,&value,1,FMpi::GetType(value),comm);
return value;
}
/** Get a value with memory cast at a specified index */
template <class ClassType>
ClassType getValue(const int ind){
ClassType value;
int previousIndex = ind;
seek(int(sizeof(value)+ind));
MPI_Unpack(array.get(),arrayCapacity,&previousIndex,&value,1,FMpi::GetType(value),comm);
return value;
}
/** Fill a value with memory cast */
template <class ClassType>
void fillValue(ClassType* const inValue){
int previousIndex = currentIndex;
seek(int(sizeof(ClassType) + previousIndex));
MPI_Unpack(array.get(),arrayCapacity,&previousIndex,inValue,1,FMpi::GetType(*inValue),comm);
}
/** Fill one/many value(s) with memcpy */
template <class ClassType>
void fillArray(ClassType* const inArray, const int inSize){
int previousIndex = currentIndex;
seek(int(sizeof(ClassType) * inSize + previousIndex));
MPI_Unpack(array.get(),arrayCapacity,&previousIndex,inArray,inSize,FMpi::GetType(*inArray),comm);
}
/** Same as fillValue */
template <class ClassType>
FMpiBufferReader& operator>>(ClassType& object){
fillValue(&object);
return *this;
}
}
/** Get the read position */
int tell() const {
return currentIndex;
}
/** Get a value with memory cast */
template <class ClassType>
ClassType getValue(){
ClassType value;
int previousIndex = currentIndex;
seek(sizeof(value) + previousIndex);
MPI_Unpack(array.get(),arrayCapacity,&previousIndex,&value,1,FMpi::GetType(value),comm);
return value;
}
/** Get a value with memory cast at a specified index */
template <class ClassType>
ClassType getValue(const int ind){
ClassType value;
int previousIndex = ind;
seek(sizeof(value)+ind);
MPI_Unpack(array.get(),arrayCapacity,&previousIndex,&value,1,FMpi::GetType(value),comm);
return value;
}
/** Fill a value with memory cast */
template <class ClassType>
void fillValue(ClassType* const inValue){
int previousIndex = currentIndex;
seek(sizeof(ClassType) + previousIndex);
MPI_Unpack(array.get(),arrayCapacity,&previousIndex,inValue,1,FMpi::GetType(*inValue),comm);
}
/** Fill one/many value(s) with memcpy */
template <class ClassType>
void fillArray(ClassType* const inArray, const int inSize){
int previousIndex = currentIndex;
seek(sizeof(ClassType) * inSize + previousIndex);
MPI_Unpack(array.get(),arrayCapacity,&previousIndex,inArray,inSize,FMpi::GetType(*inArray),comm);
}
/** Same as fillValue */
template <class ClassType>
FMpiBufferReader& operator>>(ClassType& object){
fillValue(&object);
return *this;
}
};
#endif
......
......@@ -28,102 +28,99 @@
* finally use data pointer as you like
*/
class FMpiBufferWriter : public FAbstractBufferWriter {
const MPI_Comm mpiComm; //< Communicator needed by MPI_Pack functions
const int arrayCapacity; //< Allocated Space
std::unique_ptr<char[]> array; //< Allocated Array
int currentIndex; //< Currently filled space
/** Test and exit if not enought space */
void assertRemainingSpace(const size_t requestedSpace) const {
if(int(currentIndex + requestedSpace) > arrayCapacity){
printf("Error FMpiBufferWriter has not enough space\n");
exit(0);
const MPI_Comm mpiComm; //< Communicator needed by MPI_Pack functions
const int arrayCapacity; //< Allocated Space
std::unique_ptr<char[]> array; //< Allocated Array
int currentIndex; //< Currently filled space
/** Test and exit if not enought space */
void assertRemainingSpace(const size_t requestedSpace) const {
if(int(currentIndex + requestedSpace) > arrayCapacity){
printf("Error FMpiBufferWriter has not enough space\n");
exit(0);
}
}
}
public:
/** Constructor with a default arrayCapacity of 512 bytes */
FMpiBufferWriter(const MPI_Comm inComm, const int inCapacity = 1024):
mpiComm(inComm),
arrayCapacity(inCapacity),
array(new char[inCapacity]),
currentIndex(0)
{}
/** Destructor */
virtual ~FMpiBufferWriter(){
}
/** Get allocated memory pointer */
char* data(){
return array.get();
}
/** Get allocated memory pointer */
const char* data() const {
return array.get();
}
/** Get the filled space */
int getSize() const {
return currentIndex;
}
/** Get the allocated space */
int getCapacity() const {
return arrayCapacity;
}
/** 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);
}
/**
/** Constructor with a default arrayCapacity of 512 bytes */
FMpiBufferWriter(const MPI_Comm inComm, const int inCapacity = 1024):
mpiComm(inComm),
arrayCapacity(inCapacity),
array(new char[inCapacity]),
currentIndex(0)
{}
/** Destructor */
virtual ~FMpiBufferWriter(){
}
/** Get allocated memory pointer */
char* data(){
return array.get();
}
/** Get allocated memory pointer */
const char* data() const {
return array.get();
}
/** Get the filled space */
int getSize() const {
return currentIndex;
}
/** Get the allocated space */
int getCapacity() const {
return arrayCapacity;
}
/** Write data by packing cpy */
template <class ClassType>
void write(const ClassType& object){
assertRemainingSpace(sizeof(ClassType));
MPI_Pack(const_cast<ClassType*>(&object), 1, FMpi::GetType(object), array.get(), arrayCapacity, &currentIndex, mpiComm);
}
/**
* Allow to pass rvalue to write
*/
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);
}
/** Write back, position + sizeof(object) has to be < size */
template <class ClassType>
void writeAt(const int position, const ClassType& object){
if(position + (int) sizeof(ClassType) > currentIndex){
printf("Not enought space\n");
exit(0);
template <class ClassType>
void write(const ClassType&& object){
assertRemainingSpace(sizeof(ClassType));
MPI_Pack(const_cast<ClassType*>(&object), 1, FMpi::GetType(object), array.get(), arrayCapacity, &currentIndex, mpiComm);
}
int noConstPosition = position;
MPI_Pack(const_cast<ClassType*>(&object), 1, FMpi::GetType(object), array.get(), arrayCapacity, &noConstPosition, mpiComm);
}
/** Write an array
/** Write back, position + sizeof(object) has to be < size */
template <class ClassType>
void writeAt(const int position, const ClassType& object){
if(position + (int) sizeof(ClassType) > currentIndex){
printf("Not enought space\n");
exit(0);
}
int noConstPosition = position;
MPI_Pack(const_cast<ClassType*>(&object), 1, FMpi::GetType(object), array.get(), arrayCapacity, &noConstPosition, mpiComm);
}
/** Write an array
* Warning : inSize is a number of ClassType object to write, not a size in bytes
*/
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);
}
/** Equivalent to write */
template <class ClassType>
FMpiBufferWriter& operator<<(const ClassType& object){
write(object);
return *this;
}
/** Reset the writing index, but do not change the arrayCapacity */
void reset(){
currentIndex = 0;
}
template <class ClassType>
void write(const ClassType* const objects, const int inSize){
assertRemainingSpace(sizeof(ClassType) * inSize);
MPI_Pack( const_cast<ClassType*>(objects), inSize, FMpi::GetType(*objects), array.get(), arrayCapacity, &currentIndex, mpiComm);
}
/** Equivalent to write */
template <class ClassType>
FMpiBufferWriter& operator<<(const ClassType& object){
write(object);
return *this;
}
/** Reset the writing index, but do not change the arrayCapacity */
void reset(){
currentIndex = 0;
}
};
......
......@@ -69,13 +69,13 @@ public:
template<class OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass>
static FAbstractAlgorithm* BuildAlgorithm(OctreeClass*const tree, KernelClass*const kernel,
const MPI_Comm mpiComm = (MPI_Comm)0, const bool isPeriodic = false,
const int periodicUpperlevel = 0, const int inPeriodicDirections = AllDirs){
const int periodicUpperlevel = 0){
#ifdef ScalFMM_USE_MPI
if(isPeriodic == false){
return new FFmmAlgorithmThreadProc<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass>(FMpi::FComm(mpiComm), tree, kernel);
}
else{
auto algo = new FFmmAlgorithmThreadProcPeriodic<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass>(FMpi::FComm(mpiComm), tree, periodicUpperlevel, inPeriodicDirections);
auto algo = new FFmmAlgorithmThreadProcPeriodic<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass>(FMpi::FComm(mpiComm), tree, periodicUpperlevel);
algo->setKernel(kernel);
return algo;
}
......@@ -84,7 +84,7 @@ public:
return new FFmmAlgorithmThread<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass>(tree, kernel);
}
else{
auto algo = new FFmmAlgorithmPeriodic<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass>(tree, periodicUpperlevel, inPeriodicDirections);
auto algo = new FFmmAlgorithmPeriodic<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass>(tree, periodicUpperlevel);
algo->setKernel(kernel);
return algo;
}
......
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
// 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".
// ===================================================================================
#ifndef FFMMALGORITHMPERIODIC_HPP
#define FFMMALGORITHMPERIODIC_HPP
......@@ -50,7 +35,6 @@ class FFmmAlgorithmPeriodic : public FAbstractAlgorithm{
const int OctreeHeight; //< The height of the octree (real height)
const int nbLevelsAboveRoot; //< The nb of level the user ask to go above the tree (>= -1)
const int offsetRealTree; //< nbLevelsAboveRoot GetFackLevel
const int periodicDirections;
public:
......@@ -61,10 +45,9 @@ public:
* @param inUpperLevel this parameter defines the behavior of the periodicity refer to the main doc
*
*/
FFmmAlgorithmPeriodic(OctreeClass* const inTree, const int inUpperLevel = 0, const int inPeriodicDirections = AllDirs)
FFmmAlgorithmPeriodic(OctreeClass* const inTree, const int inUpperLevel = 0)
: tree(inTree) , kernels(0), OctreeHeight(tree->getHeight()),
nbLevelsAboveRoot(inUpperLevel), offsetRealTree(inUpperLevel + 3),
periodicDirections(inPeriodicDirections) {
nbLevelsAboveRoot(inUpperLevel), offsetRealTree(inUpperLevel + 3) {
FAssertLF(tree, "tree cannot be null");
FAssertLF(-1 <= inUpperLevel, "inUpperLevel cannot be < -1");
......@@ -86,7 +69,7 @@ public:
*/
void execute(const unsigned operationsToProceed = FFmmNearAndFarFields){
FAssertLF(kernels, "kernels cannot be null");
FTRACE( FTrace::FFunction functionTrace(__FUNCTION__, "Fmm" , __FILE__ , __LINE__) );
FTRACE( FTrace::FFunction functionTrace(__FUNCTION__, "Fmm" , __FILE__ , __LINE__) );
if(operationsToProceed & FFmmP2M) bottomPass();
......@@ -195,7 +178,7 @@ public:
const int fackLevel = idxLevel + offsetRealTree;
// for each cells
do{
const int counter = tree->getPeriodicInteractionNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(), idxLevel, periodicDirections);
const int counter = tree->getPeriodicInteractionNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(), idxLevel, AllDirs);
FLOG(computationCounter.tic());
if(counter) kernels->M2L( octreeIterator.getCurrentCell() , neighbors, counter, fackLevel);
FLOG(computationCounter.tac());
......@@ -279,7 +262,7 @@ public:
// need the current particles and neighbors particles
const FTreeCoordinate centerOfLeaf = octreeIterator.getCurrentGlobalCoordinate();
const int counter = tree->getPeriodicLeafsNeighbors( neighbors, offsets, &hasPeriodicLeaves, centerOfLeaf, heightMinusOne, periodicDirections);
const int counter = tree->getPeriodicLeafsNeighbors( neighbors, offsets, &hasPeriodicLeaves, centerOfLeaf, heightMinusOne, AllDirs);
int periodicNeighborsCounter = 0;
if(hasPeriodicLeaves){
......@@ -311,7 +294,7 @@ public:
FLOG(computationCounterP2P.tac());
for(int idxNeig = 0 ; idxNeig < 27 ; ++idxNeig){
if( periodicNeighbors[idxNeig] ){
if( periodicNeighbors[idxNeig] ){
FReal*const positionsX = periodicNeighbors[idxNeig]->getWPositions()[0];
FReal*const positionsY = periodicNeighbors[idxNeig]->getWPositions()[1];
FReal*const positionsZ = periodicNeighbors[idxNeig]->getWPositions()[2];
......@@ -455,49 +438,38 @@ public:
return (((x+3)*7) + (y+3))*7 + (z + 3);
}
/** To know how many times the box is repeated in each direction
* -x +x -y +y -z +z
* If the periodicy is not in all direction this number is unchanged
* because it contains the theorical periodic box width for the
* nbLevelsAboveRoot choosen
*/
long long int theoricalRepetition() const {
return nbLevelsAboveRoot == -1 ? 3 : 3LL * (1LL<<(nbLevelsAboveRoot+1)) + 1LL;
if( nbLevelsAboveRoot == -1 ){
// we know it is 3 (-1;+1)
return 3;
}
// Else we find the repetition in one dir and double it
const long long int oneDirectionRepetition = (1<<(nbLevelsAboveRoot+2)); // 2^nbLevelsAboveRoot in each dim
const long long int fullRepetition = 2 * oneDirectionRepetition;
return fullRepetition;
}
/** To know the number of box repeated in each direction
* @param min the number of repetition for -x,-y,-z
* @param max the number of repetition for x,y,z
* The mins value are contains between [-(theoricalRepetition-1 / 2); 0]
* The maxs value are contains between [0;(theoricalRepetition-1 / 2)]
*/
void repetitionsIntervals(FTreeCoordinate*const min, FTreeCoordinate*const max) const {
const int halfRepeated = int((theoricalRepetition()-1)/2);
min->setPosition(-ifDir(DirMinusX,halfRepeated,0),-ifDir(DirMinusY,halfRepeated,0),
-ifDir(DirMinusZ,halfRepeated,0));
max->setPosition(ifDir(DirPlusX,halfRepeated,0),ifDir(DirPlusY,halfRepeated,0),
ifDir(DirPlusZ,halfRepeated,0));
if( nbLevelsAboveRoot == -1 ){
// We know it is (-1;1)
min->setPosition(-1,-1,-1);
max->setPosition(1,1,1);
}
else{
const int halfRepeated = int(theoricalRepetition()/2);
min->setPosition(-halfRepeated,-halfRepeated,-halfRepeated);
// if we repeat the box 8 times, we go from [-4 to 3]
max->setPosition(halfRepeated-1,halfRepeated-1,halfRepeated-1);
}
}
/** To get the number of repetition in all direction (x,y,z)
* @return the number of box in all directions
* Each value is between [1;theoricalRepetition]
*/
FTreeCoordinate repetitions() const {
const int halfRepeated = int((theoricalRepetition()-1)/2);
return FTreeCoordinate(ifDir(DirMinusX,halfRepeated,0) + ifDir(DirPlusX,halfRepeated,0) + 1,
ifDir(DirMinusY,halfRepeated,0) + ifDir(DirPlusY,halfRepeated,0) + 1,
ifDir(DirMinusZ,halfRepeated,0) + ifDir(DirPlusZ,halfRepeated,0) + 1);
}
/** This function has to be used to init the kernel with correct args
* it return the box seen from a kernel point of view from the periodicity the user ask for
* this is computed using the originalBoxWidth given in parameter
* @param originalBoxWidth the real system size
* @return the size the kernel should use
*/
FReal extendedBoxWidth() const {
return tree->getBoxWidth() * FReal(1<<offsetRealTree);
// The simulation box is repeated is repeated 4 times if nbLevelsAboveRoot==-1
// And then it doubles by two
return tree->getBoxWidth() * FReal(1<<(nbLevelsAboveRoot+3));
}
/** This function has to be used to init the kernel with correct args
......@@ -508,11 +480,11 @@ public:
* @return the center the kernel should use
*/
FPoint extendedBoxCenter() const {
const FReal originalBoxWidth = tree->getBoxWidth();
const FReal originalBoxWidth = tree->getBoxWidth();
const FReal originalBoxWidthDiv2 = originalBoxWidth/2.0;
const FPoint originalBoxCenter = tree->getBoxCenter();
const FPoint originalBoxCenter = tree->getBoxCenter();
const FReal offset = extendedBoxWidth()/2;
const FReal offset = extendedBoxWidth()/FReal(2.0);
return FPoint( originalBoxCenter.getX() - originalBoxWidthDiv2 + offset,
originalBoxCenter.getY() - originalBoxWidthDiv2 + offset,
originalBoxCenter.getZ() - originalBoxWidthDiv2 + offset);
......@@ -529,25 +501,6 @@ public:
return OctreeHeight + offsetRealTree;
}
/** To know if a direction is used
* @param testDir a direction to test
* @return true if the direction is used else false
*/
bool usePerDir(const int testDir) const{
return TestPeriodicCondition(periodicDirections , PeriodicCondition(testDir));
}
/** To enable quick test of the direction
* @param testDir the direction to test
* @param correctValue the value to return if direction is use
* @param wrongValue the value to return if direction is not use
* @return correctValue if testDir is used, else wrongValue
*/
template <class T>
const T& ifDir(const PeriodicCondition testDir, const T& correctValue, const T& wrongValue) const {
return (periodicDirections & testDir ? correctValue : wrongValue);
}
/** Periodicity Core
* This function is split in several part:
* 1 - special case managment
......@@ -563,638 +516,84 @@ public:
FTRACE( FTrace::FFunction functionTrace(__FUNCTION__, "Fmm" , __FILE__ , __LINE__) );
FLOG( FLog::Controller.write("\tStart Periodic Pass\n").write(FLog::Flush); );
FLOG(FTic counterTime);
/////////////////////////////////////////////////////
// If nb level == -1 nothing to do
if( nbLevelsAboveRoot == -1 ){
FLOG( FLog::Controller << "\tFinished (@Periodic = " << counterTime.tacAndElapsed() << "s)\n" );
return;
}
/////////////////////////////////////////////////////
// if nb level == 0 only M2L at real root level
if( nbLevelsAboveRoot == 0 ){
CellClass rootUp;
// compute the root