Commit 25e7da9d authored by BLANCHARD Pierre's avatar BLANCHARD Pierre

Solved merge conflicts in FCheb/UnifCell, fixed and updated FCheb/UnifCell...

Solved merge conflicts in FCheb/UnifCell, fixed and updated FCheb/UnifCell (templates are now <ORDER,NRHS,NLHS,NVALS>), some more cleanup (USE_FFT, ordering of args in new kernels ctor).
parents 12fa0584 33701cb2
......@@ -35,6 +35,11 @@ protected:
void restore(BufferReaderClass&){
static_assert(sizeof(BufferReaderClass) == 0 , "Your class should implement restore");
}
static int GetSize(){
return 0;
}
};
#endif // FABSTRACTSERIALIZABLE_HPP
......@@ -75,6 +75,10 @@ public:
buffer >> dataDown >> dataUp;
}
static int GetSize(){
return sizeof(long long int)*2;
}
/////////////////////////////////////////////////
/** Serialize only up data in a buffer */
......
......@@ -62,7 +62,7 @@
template<class OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass>
class FFmmAlgorithmThreadProc : public FAbstractAlgorithm {
static const int MaxSizePerCell = 1024;
const int MaxSizePerCell = CellClass::GetSize();
OctreeClass* const tree; //< The octree to work on
KernelClass** kernels; //< The kernels
......@@ -191,7 +191,7 @@ public:
currentLimit >>= 3;
}
}
printf("Proc::%d From leaf %lld to leaf %lld\n",idProcess,myLastInterval.min,myLastInterval.max);
// We get the min/max indexes from each procs
FMpi::MpiAssert( MPI_Allgather( myIntervals, int(sizeof(Interval)) * OctreeHeight, MPI_BYTE,
workingIntervalsPerLevel, int(sizeof(Interval)) * OctreeHeight, MPI_BYTE, comm.getComm()), __LINE__ );
......@@ -325,43 +325,44 @@ private:
if( child[idxChild] && getWorkingInterval((idxLevel+1), idProcess).min <= child[idxChild]->getMortonIndex() ){
child[idxChild]->serializeUp(sendBuffer);
state = char(state | (0x1 << idxChild));
}
}
sendBuffer.writeAt(0,state);
while( sendToProc && iterArray[cellsToSend].getCurrentGlobalIndex() == getWorkingInterval(idxLevel , sendToProc - 1).max){
while( sendToProc && iterArray[cellsToSend].getCurrentGlobalIndex() <= getWorkingInterval(idxLevel , sendToProc - 1).max){
--sendToProc;
}
MPI_Isend(sendBuffer.data(), sendBuffer.getSize(), MPI_PACKED, sendToProc,
FMpi::TagFmmM2M, comm.getComm(), &requests[iterRequests++]);
}
// We may need to receive something
bool hasToReceive = false;
int endProcThatSend = firstProcThatSend;
if(idProcess != nbProcess - 1){ // if I'm the last one (idProcess == nbProcess-1), I shall not receive anything in a M2M
while(firstProcThatSend < nbProcess
&& (getWorkingInterval((idxLevel+1), firstProcThatSend).max) < (getWorkingInterval((idxLevel+1), idProcess).max)){
&& (getWorkingInterval((idxLevel+1), firstProcThatSend).max) <= (getWorkingInterval((idxLevel+1), idProcess).max)){
// Second condition :: while firstProcThatSend max morton index is < to myself max interval
++firstProcThatSend;
}
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;
}
if(firstProcThatSend != endProcThatSend){
hasToReceive = true;
for(int idxProc = firstProcThatSend ; idxProc < endProcThatSend ; ++idxProc ){
MPI_Irecv(&recvBuffer.data()[idxProc * recvBufferOffset], recvBufferOffset, MPI_PACKED,
idxProc, FMpi::TagFmmM2M, comm.getComm(), &requests[iterRequests++]);
......@@ -369,6 +370,7 @@ private:
}
}
}
FLOG(prepareCounter.tac());
FTRACE( regionTrace.end() );
......@@ -395,7 +397,7 @@ private:
if( hasToReceive ){
CellClass* currentChild[8];
memcpy(currentChild, iterArray[numberOfCells - 1].getCurrentChild(), 8 * sizeof(CellClass*));
// retreive data and merge my child and the child from others
for(int idxProc = firstProcThatSend ; idxProc < endProcThatSend ; ++idxProc){
recvBuffer.seek(idxProc * recvBufferOffset);
......@@ -407,12 +409,12 @@ private:
state >>= 1;
++position;
}
FAssertLF(!currentChild[position], "Already has a cell here");
FAssertLF(!currentChild[position], "Already has a cell here");
recvBufferCells[position].deserializeUp(recvBuffer);
currentChild[position] = (CellClass*) &recvBufferCells[position];
state >>= 1;
++position;
}
......@@ -423,7 +425,6 @@ private:
(*kernels[0]).M2M( iterArray[numberOfCells - 1].getCurrentCell() , currentChild, idxLevel);
FLOG(computationCounter.tac());
firstProcThatSend = endProcThatSend - 1;
}
}
......@@ -815,7 +816,7 @@ private:
const int heightMinusOne = OctreeHeight - 1;
FMpiBufferWriter sendBuffer(comm.getComm());
FMpiBufferWriter sendBuffer(comm.getComm(),MaxSizePerCell);
FMpiBufferReader recvBuffer(comm.getComm(),MaxSizePerCell);
// for each levels exepted leaf level
......
......@@ -35,10 +35,11 @@
class FEwalLoader : public FAbstractLoader {
protected:
std::ifstream file; //< The file to read
FPoint centerOfBox; //< The center of box read from file
FReal boxWidth; //< the box width read from file
int nbParticles; //< the number of particles read from file
int levcfg ; //< DL_POLY CONFIG file key. 0,1 or 2
FPoint centerOfBox; //< The center of box read from file
FReal boxWidth; //< the box width read from file
int nbParticles; //< the number of particles read from file
int levcfg ; //< DL_POLY CONFIG file key. 0,1 or 2
FReal energy ;
public:
enum Type{
OW,
......@@ -66,7 +67,7 @@ public:
int imcon ;
//int tempi(0);
FReal tempf(0);
file >> levcfg >> imcon >> this->nbParticles;
file >> levcfg >> imcon >> this->nbParticles >> this->energy;
// Periodic case
if( imcon > 0 ) {
FReal widthx, widthy, widthz;
......@@ -83,7 +84,7 @@ public:
this->centerOfBox.setPosition(0.0,0.0,0.0);
}
else {
this->boxWidth = 0;
this->boxWidth = 0;
this->nbParticles = 0;
}
}
......@@ -127,7 +128,7 @@ public:
}
FReal getEnergy() const{
return static_cast<FReal>(0.0);
return this->energy;
}
/**
......@@ -232,7 +233,7 @@ public:
centerOfBox.setPosition(0.0,0.0,0.0);
}
else {
this->boxWidth = 0;
this->boxWidth = 0;
this->nbParticles = 0;
}
}
......
......@@ -45,95 +45,95 @@ class FAbstractChebKernel : public FAbstractKernels< CellClass, ContainerClass>
{
protected:
enum {nnodes = TensorTraits<ORDER>::nnodes};
typedef FChebInterpolator<ORDER> InterpolatorClass;
/// Needed for P2M, M2M, L2L and L2P operators
const FSmartPointer<InterpolatorClass,FSmartPointerMemory> Interpolator;
/// Needed for P2P operator
const FSmartPointer<MatrixKernelClass,FSmartPointerMemory> MatrixKernel;
/// Height of the entire oct-tree
const unsigned int TreeHeight;
/// Corner of oct-tree box
const FPoint BoxCorner;
/// Width of oct-tree box
const FReal BoxWidth;
/// Width of a leaf cell box
const FReal BoxWidthLeaf;
/**
* Compute center of leaf cell from its tree coordinate.
* @param[in] Coordinate tree coordinate
* @return center of leaf cell
*/
const FPoint getLeafCellCenter(const FTreeCoordinate& Coordinate) const
{
return FPoint(BoxCorner.getX() + (FReal(Coordinate.getX()) + FReal(.5)) * BoxWidthLeaf,
BoxCorner.getY() + (FReal(Coordinate.getY()) + FReal(.5)) * BoxWidthLeaf,
BoxCorner.getZ() + (FReal(Coordinate.getZ()) + FReal(.5)) * BoxWidthLeaf);
}
typedef FChebInterpolator<ORDER> InterpolatorClass;
/// Needed for P2M, M2M, L2L and L2P operators
const FSmartPointer<InterpolatorClass,FSmartPointerMemory> Interpolator;
/// Needed for P2P operator
const FSmartPointer<MatrixKernelClass,FSmartPointerMemory> MatrixKernel;
/// Height of the entire oct-tree
const unsigned int TreeHeight;
/// Corner of oct-tree box
const FPoint BoxCorner;
/// Width of oct-tree box
const FReal BoxWidth;
/// Width of a leaf cell box
const FReal BoxWidthLeaf;
/**
* Compute center of leaf cell from its tree coordinate.
* @param[in] Coordinate tree coordinate
* @return center of leaf cell
*/
const FPoint getLeafCellCenter(const FTreeCoordinate& Coordinate) const
{
return FPoint(BoxCorner.getX() + (FReal(Coordinate.getX()) + FReal(.5)) * BoxWidthLeaf,
BoxCorner.getY() + (FReal(Coordinate.getY()) + FReal(.5)) * BoxWidthLeaf,
BoxCorner.getZ() + (FReal(Coordinate.getZ()) + FReal(.5)) * BoxWidthLeaf);
}
public:
/**
* The constructor initializes all constant attributes and it reads the
* precomputed and compressed M2L operators from a binary file (an
* runtime_error is thrown if the required file is not valid).
*/
FAbstractChebKernel(const int inTreeHeight,
const FPoint& inBoxCenter,
const FReal inBoxWidth)
: Interpolator(new InterpolatorClass()),
MatrixKernel(new MatrixKernelClass()),
TreeHeight(inTreeHeight),
BoxCorner(inBoxCenter - inBoxWidth / FReal(2.)),
BoxWidth(inBoxWidth),
BoxWidthLeaf(BoxWidth / FReal(FMath::pow(2, inTreeHeight - 1)))
{
/* empty */
}
/**
* The constructor initializes all constant attributes and it reads the
* precomputed and compressed M2L operators from a binary file (an
* runtime_error is thrown if the required file is not valid).
*/
FAbstractChebKernel(const int inTreeHeight,
const FReal inBoxWidth,
const FPoint& inBoxCenter)
: Interpolator(new InterpolatorClass()),
MatrixKernel(new MatrixKernelClass()),
TreeHeight(inTreeHeight),
BoxCorner(inBoxCenter - inBoxWidth / FReal(2.)),
BoxWidth(inBoxWidth),
BoxWidthLeaf(BoxWidth / FReal(FMath::pow(2, inTreeHeight - 1)))
{
/* empty */
}
virtual ~FAbstractChebKernel(){
// should not be used
}
virtual ~FAbstractChebKernel(){
// should not be used
}
const InterpolatorClass *const getPtrToInterpolator() const
{ return Interpolator.getPtr(); }
const InterpolatorClass *const getPtrToInterpolator() const
{ return Interpolator.getPtr(); }
virtual void P2M(CellClass* const LeafCell,
const ContainerClass* const SourceParticles) = 0;
virtual void P2M(CellClass* const LeafCell,
const ContainerClass* const SourceParticles) = 0;
virtual void M2M(CellClass* const FRestrict ParentCell,
const CellClass*const FRestrict *const FRestrict ChildCells,
const int TreeLevel) = 0;
virtual void M2M(CellClass* const FRestrict ParentCell,
const CellClass*const FRestrict *const FRestrict ChildCells,
const int TreeLevel) = 0;
virtual void M2L(CellClass* const FRestrict TargetCell,
const CellClass* SourceCells[343],
const int NumSourceCells,
const int TreeLevel) = 0;
virtual void M2L(CellClass* const FRestrict TargetCell,
const CellClass* SourceCells[343],
const int NumSourceCells,
const int TreeLevel) = 0;
virtual void L2L(const CellClass* const FRestrict ParentCell,
CellClass* FRestrict *const FRestrict ChildCells,
const int TreeLevel) = 0;
virtual void L2L(const CellClass* const FRestrict ParentCell,
CellClass* FRestrict *const FRestrict ChildCells,
const int TreeLevel) = 0;
virtual void L2P(const CellClass* const LeafCell,
ContainerClass* const TargetParticles) = 0;
virtual void L2P(const CellClass* const LeafCell,
ContainerClass* const TargetParticles) = 0;
virtual void P2P(const FTreeCoordinate& /* LeafCellCoordinate */, // needed for periodic boundary conditions
ContainerClass* const FRestrict TargetParticles,
const ContainerClass* const FRestrict /*SourceParticles*/,
ContainerClass* const NeighborSourceParticles[27],
const int /* size */) = 0;
virtual void P2P(const FTreeCoordinate& /* LeafCellCoordinate */, // needed for periodic boundary conditions
ContainerClass* const FRestrict TargetParticles,
const ContainerClass* const FRestrict /*SourceParticles*/,
ContainerClass* const NeighborSourceParticles[27],
const int /* size */) = 0;
virtual void P2PRemote(const FTreeCoordinate& /*inPosition*/,
ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
ContainerClass* const inNeighbors[27], const int /*inSize*/) = 0;
virtual void P2PRemote(const FTreeCoordinate& /*inPosition*/,
ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
ContainerClass* const inNeighbors[27], const int /*inSize*/) = 0;
};
......
......@@ -16,7 +16,7 @@
#ifndef FCHEBCELL_HPP
#define FCHEBCELL_HPP
#include "../../Components/FBasicCell.hpp"
#include "../../Extensions/FExtendMortonIndex.hpp"
#include "../../Extensions/FExtendCoordinate.hpp"
......@@ -24,76 +24,118 @@
#include "../../Extensions/FExtendCellType.hpp"
/**
* @author Matthias Messner (matthias.messner@inria.fr)
* @class FChebCell
* Please read the license
*
* This class defines a cell used in the Chebyshev based FMM.
* @param NVALS is the number of right hand side.
*/
template <int ORDER, int NMUL = 1, int NLOC = 1>
class FChebCell : public FExtendMortonIndex, public FExtendCoordinate
* @author Matthias Messner (matthias.messner@inria.fr)
* @class FChebCell
* Please read the license
*
* This class defines a cell used in the Chebyshev based FMM.
* @param NVALS is the number of right hand side.
*/
template <int ORDER, int NRHS = 1, int NLHS = 1, int NVALS = 1>
class FChebCell : public FBasicCell
{
static const int VectorSize = TensorTraits<ORDER>::nnodes * 2;
static const int VectorSize = TensorTraits<ORDER>::nnodes * 2;
FReal multipole_exp[NRHS * NVALS * VectorSize]; //< Multipole expansion
FReal local_exp[NLHS * NVALS * VectorSize]; //< Local expansion
FReal multipole_exp[NMUL * VectorSize]; //< Multipole expansion
FReal local_exp[NLOC * VectorSize]; //< Local expansion
public:
FChebCell(){
memset(multipole_exp, 0, sizeof(FReal) * NMUL * VectorSize);
memset(local_exp, 0, sizeof(FReal) * NLOC * VectorSize);
}
FChebCell(){
memset(multipole_exp, 0, sizeof(FReal) * NRHS * NVALS * VectorSize);
memset(local_exp, 0, sizeof(FReal) * NLHS * NVALS * VectorSize);
}
~FChebCell() {}
~FChebCell() {}
/** Get Multipole */
const FReal* getMultipole(const int inRhs) const
{ return this->multipole_exp + inRhs*VectorSize;
}
/** Get Local */
const FReal* getLocal(const int inRhs) const{
return this->local_exp + inRhs*VectorSize;
}
/** Get Multipole */
const FReal* getMultipole(const int inRhs) const
{ return this->multipole_exp + inRhs*VectorSize;
}
/** Get Local */
const FReal* getLocal(const int inRhs) const{
return this->local_exp + inRhs*VectorSize;
}
/** Get Multipole */
FReal* getMultipole(const int inRhs){
return this->multipole_exp + inRhs*VectorSize;
}
/** Get Local */
FReal* getLocal(const int inRhs){
return this->local_exp + inRhs*VectorSize;
}
/** Get Multipole */
FReal* getMultipole(const int inRhs){
return this->multipole_exp + inRhs*VectorSize;
}
/** Get Local */
FReal* getLocal(const int inRhs){
return this->local_exp + inRhs*VectorSize;
}
/** To get the leading dim of a vec */
int getVectorSize() const{
return VectorSize;
}
/** To get the leading dim of a vec */
int getVectorSize() const{
return VectorSize;
}
/** Make it like the begining */
void resetToInitialState(){
memset(multipole_exp, 0, sizeof(FReal) * NRHS * NVALS * VectorSize);
memset(local_exp, 0, sizeof(FReal) * NLHS * NVALS * VectorSize);
}
///////////////////////////////////////////////////////
// to extend FAbstractSendable
///////////////////////////////////////////////////////
template <class BufferWriterClass>
void serializeUp(BufferWriterClass& buffer) const{
buffer.write(multipole_exp, VectorSize*NVALS*NRHS);
}
template <class BufferReaderClass>
void deserializeUp(BufferReaderClass& buffer){
buffer.fillArray(multipole_exp, VectorSize*NVALS*NRHS);
}
template <class BufferWriterClass>
void serializeDown(BufferWriterClass& buffer) const{
buffer.write(local_exp, VectorSize*NVALS*NLHS);
}
template <class BufferReaderClass>
void deserializeDown(BufferReaderClass& buffer){
buffer.fillArray(local_exp, VectorSize*NVALS*NLHS);
}
///////////////////////////////////////////////////////
// to extend Serializable
///////////////////////////////////////////////////////
template <class BufferWriterClass>
void save(BufferWriterClass& buffer) const{
FBasicCell::save(buffer);
buffer.write(multipole_exp, VectorSize*NVALS*NRHS);
buffer.write(local_exp, VectorSize*NVALS*NLHS);
}
template <class BufferReaderClass>
void restore(BufferReaderClass& buffer){
FBasicCell::restore(buffer);
buffer.fillArray(multipole_exp, VectorSize*NVALS*NRHS);
buffer.fillArray(local_exp, VectorSize*NVALS*NLHS);
}
static int GetSize(){
return sizeof(FReal) * VectorSize*(NRHS+NLHS)*NVALS;
}
/** Make it like the begining */
void resetToInitialState(){
memset(multipole_exp, 0, sizeof(FReal) * NMUL * VectorSize);
memset(local_exp, 0, sizeof(FReal) * NLOC * VectorSize);
}
};
template <int ORDER, int NMUL = 1, int NLOC = 1>
class FTypedChebCell : public FChebCell<ORDER,NMUL,NLOC>, public FExtendCellType {
template <int ORDER, int NRHS = 1, int NLHS = 1, int NVALS = 1>
class FTypedChebCell : public FChebCell<ORDER,NRHS,NLHS,NVALS>, public FExtendCellType {
public:
template <class BufferWriterClass>
void save(BufferWriterClass& buffer) const{
FChebCell<ORDER,NMUL,NLOC>::save(buffer);
FExtendCellType::save(buffer);
}
template <class BufferReaderClass>
void restore(BufferReaderClass& buffer){
FChebCell<ORDER,NMUL,NLOC>::restore(buffer);
FExtendCellType::restore(buffer);
}
void resetToInitialState(){
FChebCell<ORDER,NMUL,NLOC>::resetToInitialState();
FExtendCellType::resetToInitialState();
}
template <class BufferWriterClass>
void save(BufferWriterClass& buffer) const{
FChebCell<ORDER,NRHS,NLHS,NVALS>::save(buffer);
FExtendCellType::save(buffer);
}
template <class BufferReaderClass>
void restore(BufferReaderClass& buffer){
FChebCell<ORDER,NRHS,NLHS,NVALS>::restore(buffer);
FExtendCellType::restore(buffer);
}
void resetToInitialState(){
FChebCell<ORDER,NRHS,NLHS,NVALS>::resetToInitialState();
FExtendCellType::resetToInitialState();
}
};
#endif //FCHEBCELL_HPP
......
......@@ -105,9 +105,9 @@ public:
* runtime_error is thrown if the required file is not valid).
*/
FChebFlopsSymKernel(const int _inTreeHeight,
const FPoint& inBoxCenter,
const FReal inBoxWidth,
const FReal Epsilon)
const FReal inBoxWidth,
const FPoint& inBoxCenter,
const FReal Epsilon)
: MatrixKernel(new MatrixKernelClass()),
SymHandler(new SymmetryHandler(MatrixKernel.getPtr(), Epsilon)), inTreeHeight(_inTreeHeight),
flopsP2M(0), flopsM2M(0), flopsM2L(0), flopsL2L(0), flopsL2P(0), flopsP2P(0),
......
......@@ -61,13 +61,13 @@ public:
* runtime_error is thrown if the required file is not valid).
*/
FChebKernel(const int inTreeHeight,
const FPoint& inBoxCenter,
const FReal inBoxWidth,
const FReal Epsilon)
const FReal inBoxWidth,
const FPoint& inBoxCenter,
const FReal Epsilon)
: FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,
inBoxCenter,
inBoxWidth),
M2LHandler(new M2LHandlerClass(Epsilon))
inBoxWidth,
inBoxCenter),
M2LHandler(new M2LHandlerClass(Epsilon))
{
// read precomputed compressed m2l operators from binary file
//M2LHandler->ReadFromBinaryFileAndSet();
......
......@@ -107,10 +107,10 @@ public:
* runtime_error is thrown if the required file is not valid).
*/
FChebSymKernel(const int inTreeHeight,
const FPoint& inBoxCenter,
const FReal inBoxWidth,
const FReal Epsilon)
: AbstractBaseClass(inTreeHeight, inBoxCenter, inBoxWidth),
const FReal inBoxWidth,
const FPoint& inBoxCenter,
const FReal Epsilon)
: AbstractBaseClass(inTreeHeight, inBoxWidth, inBoxCenter),
SymHandler(new SymmetryHandlerClass(AbstractBaseClass::MatrixKernel.getPtr(), Epsilon, inBoxWidth, inTreeHeight)),
Loc(NULL), Mul(NULL), countExp(NULL)
{
......
......@@ -69,12 +69,12 @@ public:
* runtime_error is thrown if the required file is not valid).
*/
FChebTensorialKernel(const int inTreeHeight,
const FPoint& inBoxCenter,
const FReal inBoxWidth,
const FReal Epsilon)
: FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,