Commit 3b1d2118 authored by BRAMAS Berenger's avatar BRAMAS Berenger

Change the kernel interface to have more efficient iteration in sparse fmm

parent f1c82ed0
......@@ -320,7 +320,7 @@ public:
}
void P2P(ContainerClass* target, const ContainerClass* sources) override {
ContainerClass* sourcesArray[27] = { const_cast<ContainerClass*> (sources) };
ContainerClass* sourcesArray[1] = { const_cast<ContainerClass*> (sources) };
DirectInteractionComputer<FReal, MatrixKernelClass::NCMP, NVALS>::template P2PRemote(target,sourcesArray,1,MatrixKernel);
}
......
This diff is collapsed.
......@@ -34,12 +34,10 @@ public:
}
}
void M2L(CellClass* const FRestrict local, const CellClass* interactions[], const int counter, const int level) override {
void M2L(CellClass* const FRestrict local, const CellClass* interactions[], const int /*positions*/[], const int counter, const int level) override {
std::cout << "Usual] M2L Idx = " << local->getMortonIndex() << " at level " << level << " with\n";
for(int idxInter = 0 ; idxInter < 8 ; ++idxInter){
if(interactions[idxInter]){
for(int idxInter = 0 ; idxInter < counter ; ++idxInter){
std::cout << "\t Idx " << interactions[idxInter]->getMortonIndex() << "\n";
}
}
}
......@@ -58,23 +56,19 @@ public:
void P2P(const FTreeCoordinate& ,
ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict ,
ContainerClass* const neighs[27], const int ) override {
ContainerClass* const neighs[], const int /*positions*/[], const int size) override {
std::cout << "Usual] P2P @" << targets << " has " << targets->getNbParticles() << " with\n";
for(int idxNeigh = 0 ; idxNeigh < 27 ; ++idxNeigh){
if(neighs[idxNeigh]){
for(int idxNeigh = 0 ; idxNeigh < size ; ++idxNeigh){
std::cout << "\t @" << neighs[idxNeigh]<< " has " << neighs[idxNeigh]->getNbParticles() << "\n";
}
}
}
void P2PRemote(const FTreeCoordinate& ,
ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict ,
ContainerClass* const neighs[27], const int ) override {
ContainerClass* const neighs[], const int /*positions*/[], const int size) override {
std::cout << "Usual] P2P remote @" << targets << " has " << targets->getNbParticles() << " with\n";
for(int idxNeigh = 0 ; idxNeigh < 27 ; ++idxNeigh){
if(neighs[idxNeigh]){
for(int idxNeigh = 0 ; idxNeigh < size ; ++idxNeigh){
std::cout << "\t @" << neighs[idxNeigh]<< " has " << neighs[idxNeigh]->getNbParticles() << "\n";
}
}
}
......
......@@ -197,7 +197,7 @@ public:
}
void P2M(CellClass* const cell, const ContainerClass* const SourceParticles) {
void P2M(CellClass* const cell, const ContainerClass* const SourceParticles) override {
FSize tmpCost = countFlopsP2M(SourceParticles->getNbParticles());
flopsP2M += tmpCost;
cell->addCost(tmpCost);
......@@ -208,7 +208,7 @@ public:
void M2M(CellClass* const FRestrict cell,
const CellClass*const FRestrict *const FRestrict ChildCells,
const int /*TreeLevel*/) {
const int /*TreeLevel*/) override {
FSize flops = 0;
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex)
if (ChildCells[ChildIndex]) flops += countFlopsM2MorL2L();
......@@ -221,17 +221,18 @@ public:
void M2L(CellClass* const FRestrict cell,
const CellClass* SourceCells[343],
const int /* not needed */,
const int /* TreeLevel */)
const CellClass* /*SourceCells*/[],
const int positions[],
const int size,
const int /* TreeLevel */) override
{
FSize flops = 0;
// count how ofter each of the 16 interactions is used
memset(countExp, 0, sizeof(int) * 343);
for (unsigned int idx=0; idx<343; ++idx)
if (SourceCells[idx]) countExp[SymHandler->pindices[idx]]++;
for (int idx=0; idx<size; ++idx)
countExp[SymHandler->pindices[positions[idx]]]++;
// multiply (mat-mat-mul)
for (unsigned int pidx=0; pidx<343; ++pidx)
for (int pidx=0; pidx<343; ++pidx)
if (countExp[pidx])
flops += countFlopsM2L(countExp[pidx], SymHandler->LowRank[pidx])
+ countExp[pidx]*nnodes;
......@@ -243,7 +244,7 @@ public:
void L2L(const CellClass* const FRestrict /* not needed */,
CellClass* FRestrict *const FRestrict ChildCells,
const int /* TreeLevel*/) {
const int /* TreeLevel*/) override {
FSize flops = 0;
FSize tmpCost = 0;
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex)
......@@ -260,7 +261,7 @@ public:
void L2P(const CellClass* const cell,
ContainerClass* const TargetParticles) {
ContainerClass* const TargetParticles) override {
//// 1.a) apply Sx
//flopsL2P += countFlopsP2MorL2P(TargetParticlesParticles->getNbParticles()) + TargetParticles->getNbParticles();
//// 1.b) apply Px (grad Sx)
......@@ -281,8 +282,9 @@ public:
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 */) {
ContainerClass* const NeighborSourceParticles[],
const int positions[],
const int size) override {
FSize tmpCost = 0;
FSize srcPartCount = SourceParticles->getNbParticles();
FSize tgtPartCount = TargetParticles->getNbParticles();
......@@ -290,13 +292,11 @@ public:
if ( TargetParticles != SourceParticles ) {
tmpCost += countFlopsP2P() * tgtPartCount * srcPartCount;
for ( unsigned int idx = 0; idx < 27; ++idx ) {
if (NeighborSourceParticles[idx]) {
for ( int idx = 0; idx < size; ++idx ) {
tmpCost +=
countFlopsP2P()
* tgtPartCount
* NeighborSourceParticles[idx]->getNbParticles();
}
* NeighborSourceParticles[idx]->getNbParticles();
}
} else {
tmpCost +=
......@@ -304,15 +304,13 @@ public:
* ((tgtPartCount * tgtPartCount
+ tgtPartCount)
/ 2);
for (unsigned int idx=0; idx<=13; ++idx)
for (int idx=0; idx < size && positions[idx]<=13; ++idx)
{
if (NeighborSourceParticles[idx]) {
tmpCost +=
countFlopsP2Pmutual()
* tgtPartCount
* NeighborSourceParticles[idx]->getNbParticles();
}
}
}
flopsP2P += tmpCost;
......
......@@ -113,7 +113,7 @@ public:
// delete [] _countExp;
}
void P2M(CellClass* const cell, const ContainerClass* const SourceParticles) {
void P2M(CellClass* const cell, const ContainerClass* const SourceParticles) override {
FSize tmpCost = SourceParticles->getNbParticles();
flopsP2M += tmpCost;
cell->addCost(tmpCost);
......@@ -123,7 +123,7 @@ public:
void M2M(CellClass* const FRestrict cell,
const CellClass* const FRestrict* const FRestrict ChildrenCells,
const int /*TreeLevel*/) {
const int /*TreeLevel*/) override {
FSize flops = 0;
for ( unsigned int childIdx = 0; childIdx < 8; ++childIdx )
if ( ChildrenCells[childIdx] )
......@@ -135,16 +135,15 @@ public:
void M2L(CellClass* const FRestrict cell,
const CellClass* SourceCells[343],
const int /* not needed */,
const int /* TreeLevel */) {
const CellClass* /*SourceCells*/[],
const int /*positions*/[],
const int size,
const int /* TreeLevel */) override {
FSize flops = 0;
// count how often each of the 16 interactions is used
memset(_countExp, 0, sizeof(int) * 343);
for (unsigned int idx = 0; idx < 343; ++idx) {
if ( SourceCells[idx] ) {
flops += 1;
}
memset(_countExp, 0, sizeof(int) * 343);
for (int idx = 0; idx < size; ++idx) {
flops += 1;
}
flopsM2L += flops;
......@@ -155,7 +154,7 @@ public:
void L2L(const CellClass* const FRestrict /* not needed */,
CellClass* FRestrict *const FRestrict ChildCells,
const int /* TreeLevel*/) {
const int /* TreeLevel*/) override {
FSize flops = 0;
FSize tmpCost = 0;
for (unsigned int childIdx = 0; childIdx < 8; ++childIdx ) {
......@@ -173,7 +172,7 @@ public:
void L2P(const CellClass* const cell,
ContainerClass* const TargetParticles) {
ContainerClass* const TargetParticles) override {
FSize tmpCost = TargetParticles->getNbParticles();
flopsL2P += tmpCost;
cell->addCost(tmpCost);
......@@ -185,8 +184,9 @@ public:
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 */) {
ContainerClass* const /*NeighborSourceParticles*/[],
const int /*positions*/[],
const int size) override {
FSize srcPartCount = SourceParticles->getNbParticles();
FSize tgtPartCount = TargetParticles->getNbParticles();
......
......@@ -320,19 +320,21 @@ protected:
// Cell count in zone (for this level)
int zoneCellCount = costzones.at(threadNum).at(level).second;
// Array for cell neighbours
const CellClass* neighbours[343];
const CellClass* neighbours[342];
int neighborPositions[342];
// Call M2L kernel on cells
while(zoneCellCount-- > 0) {
const int counter =
tree->getInteractionNeighbors(
neighbours,
neighbours, neighborPositions,
zoneIterator.getCurrentGlobalCoordinate(),
level);
if(counter)
myThreadkernels->M2L(
zoneIterator.getCurrentCell(),
neighbours,
neighborPositions,
counter,
level);
zoneIterator.moveRight();
......@@ -439,7 +441,8 @@ protected:
// Cell count in zone (for this level)
int zoneCellCount;
// Cell neighbours
ContainerClass* neighbours[27];
ContainerClass* neighbours[26];
int neighborPositions[26];
// Cell data
LeafData leafdata;
......@@ -474,7 +477,7 @@ protected:
if( p2pEnabled ) {
// needs the current particles and neighbours particles
const int counter =
tree->getLeafsNeighbors(neighbours,
tree->getLeafsNeighbors(neighbours,neighborPositions,
leafdata.cell->getCoordinate(),
OctreeHeight - 1); // Leaf level
......@@ -482,6 +485,7 @@ protected:
leafdata.targets,
leafdata.sources,
neighbours,
neighborPositions,
counter);
}
......
......@@ -63,11 +63,25 @@ public:
* M2L
* Multipole to local
* @param local the element to fill using distant neighbors
* @param distantNeighbors is an array containing fathers's direct neighbors's child - direct neigbors
* @param distantNeighbors is an array containing fathers's direct neighbors's child - direct neigbors (max 189)
* @param neighborPositions the relative position of the neighbor (between O and 342)
* @param size the number of neighbors
* @param inLevel the current level of the computation
*
* @code // the code was:
* @code for(int idxNeigh = 0 ; idxNeigh < 343 ; ++idxNeigh){
* @code if(distantNeighbors[idxNeigh]){
* @code ...
* @code }
* @code }
* @code // Now it must be :
* @code for(int idxExistingNeigh = 0 ; idxExistingNeigh < size ; ++idxExistingNeigh){
* @code const int idxNeigh = neighborPositions[idxExistingNeigh]
* @code distantNeighbors[idxExistingNeigh]...
* @code }
*/
virtual void M2L(CellClass* const FRestrict local, const CellClass* distantNeighbors[343],
virtual void M2L(CellClass* const FRestrict local, const CellClass* distantNeighbors[],
const int neighborPositions[],
const int size, const int inLevel) = 0;
......@@ -113,12 +127,14 @@ public:
* @param inLeafPosition tree coordinate of the leaf
* @param targets current boxe targets particles
* @param sources current boxe sources particles (can be == to targets)
* @param directNeighborsParticles the particles from direct neighbors (this is an array of list)
* @param directNeighborsParticles the particles from direct neighbors (this is an array of list) max length = 26
* @param neighborPositions the relative position of the neighbors (between O and 25) max length = 26
* @param size the number of direct neighbors
*/
virtual void P2P(const FTreeCoordinate& inLeafPosition,
ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
ContainerClass* const directNeighborsParticles[27], const int size) = 0;
ContainerClass* const directNeighborsParticles[], const int neighborPositions[],
const int size) = 0;
/**
* P2P
......@@ -126,7 +142,8 @@ public:
* @param inLeafPosition tree coordinate of the leaf
* @param targets current boxe targets particles
* @param sources current boxe sources particles (can be == to targets)
* @param directNeighborsParticles the particles from direct neighbors (this is an array of list)
* @param directNeighborsParticles the particles from direct neighbors (this is an array of list) max length = 26
* @param neighborPositions the relative position of the neighbors (between 0 and 25) max length = 26
* @param size the number of direct neighbors
*
* This method is called by the MPI algorithm with leaves from other hosts.
......@@ -135,7 +152,8 @@ public:
*/
virtual void P2PRemote(const FTreeCoordinate& /*inLeafPosition*/,
ContainerClass* const FRestrict /*targets*/, const ContainerClass* const FRestrict /*sources*/,
ContainerClass* const /*directNeighborsParticles*/[27], const int /*size*/) {
ContainerClass* const /*directNeighborsParticles*/[],
const int /*neighborPositions*/[], const int /*size*/) {
FLOG( FLog::Controller.write("Warning, P2P remote is used but not implemented!").write(FLog::Flush) );
}
......
......@@ -35,42 +35,43 @@ public:
}
/** Do nothing */
virtual void P2M(CellClass* const , const ContainerClass* const ) {
virtual void P2M(CellClass* const /*targetCell*/, const ContainerClass* const /*sourceParticles*/) override {
}
/** Do nothing */
virtual void M2M(CellClass* const FRestrict , const CellClass*const FRestrict *const FRestrict , const int ) {
virtual void M2M(CellClass* const FRestrict /*parentCell*/, const CellClass*const FRestrict *const FRestrict /*children*/, const int /*level*/) override {
}
/** Do nothing */
virtual void M2L(CellClass* const FRestrict , const CellClass* [], const int , const int ) {
virtual void M2L(CellClass* const FRestrict /*targetLocal*/, const CellClass* /*sourceMultipoles*/[],
const int /*relativePostions*/[], const int /*nbInteractions*/, const int /*level*/) override {
}
/** Do nothing */
virtual void L2L(const CellClass* const FRestrict , CellClass* FRestrict *const FRestrict , const int ) {
virtual void L2L(const CellClass* const FRestrict /*parentCell*/, CellClass* FRestrict *const FRestrict /*children*/, const int /*level*/) override {
}
/** Do nothing */
virtual void L2P(const CellClass* const , ContainerClass* const ){
virtual void L2P(const CellClass* const /*sourceCell*/, ContainerClass* const /*targetPaticles*/) override {
}
/** Do nothing */
virtual void P2P(const FTreeCoordinate& ,
ContainerClass* const FRestrict , const ContainerClass* const FRestrict ,
ContainerClass* const [27], const int ){
virtual void P2P(const FTreeCoordinate& /*treeCoord*/,
ContainerClass* const FRestrict /*targetParticles*/, const ContainerClass* const FRestrict /*sourceParticles*/,
ContainerClass* const /*neigbhorsParticles*/[], const int /*neighborPositions*/[], const int /*nbNeighbors*/) override {
}
/** Do nothing */
virtual void P2PRemote(const FTreeCoordinate& ,
ContainerClass* const FRestrict , const ContainerClass* const FRestrict ,
ContainerClass* const [27], const int ){
virtual void P2PRemote(const FTreeCoordinate& /*treeCoord*/,
ContainerClass* const FRestrict /*targetParticles*/, const ContainerClass* const FRestrict /*sourceParticles*/,
ContainerClass* const /*neigbhorsParticles*/[], const int /*neighborPositions*/[], const int /*nbNeighbors*/) override {
}
......
......@@ -47,13 +47,13 @@ public:
}
/** Before upward */
void P2M(CellClass* const pole, const ContainerClass* const particles) {
void P2M(CellClass* const pole, const ContainerClass* const particles) override {
// the pole represents all particles under
pole->setDataUp(pole->getDataUp() + particles->getNbParticles());
}
/** During upward */
void M2M(CellClass* const FRestrict pole, const CellClass *const FRestrict *const FRestrict child, const int /*level*/) {
void M2M(CellClass* const FRestrict pole, const CellClass *const FRestrict *const FRestrict child, const int /*level*/) override {
// A parent represents the sum of the child
for(int idx = 0 ; idx < 8 ; ++idx){
if(child[idx]){
......@@ -63,17 +63,16 @@ public:
}
/** Before Downward */
void M2L(CellClass* const FRestrict pole, const CellClass* distantNeighbors[343], const int /*size*/, const int /*level*/) {
void M2L(CellClass* const FRestrict pole, const CellClass* distantNeighbors[], const int /*position*/[],
const int size, const int /*level*/) override {
// The pole is impacted by what represent other poles
for(int idx = 0 ; idx < 343 ; ++idx){
if(distantNeighbors[idx]){
for(int idx = 0 ; idx < size ; ++idx){
pole->setDataDown(pole->getDataDown() + distantNeighbors[idx]->getDataUp());
}
}
}
/** During Downward */
void L2L(const CellClass*const FRestrict local, CellClass* FRestrict *const FRestrict child, const int /*level*/) {
void L2L(const CellClass*const FRestrict local, CellClass* FRestrict *const FRestrict child, const int /*level*/) override {
// Each child is impacted by the father
for(int idx = 0 ; idx < 8 ; ++idx){
if(child[idx]){
......@@ -84,7 +83,7 @@ public:
}
/** After Downward */
void L2P(const CellClass* const local, ContainerClass*const particles){
void L2P(const CellClass* const local, ContainerClass*const particles) override {
// The particles is impacted by the parent cell
long long int*const particlesAttributes = particles->getDataDown();
for(FSize idxPart = 0 ; idxPart < particles->getNbParticles() ; ++idxPart){
......@@ -96,16 +95,14 @@ public:
/** After Downward */
void P2P(const FTreeCoordinate& ,
ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict sources,
ContainerClass* const directNeighborsParticles[27], const int ){
ContainerClass* const directNeighborsParticles[], const int /*positions*/[], const int inSize) override {
// Each particles targeted is impacted by the particles sources
long long int inc = sources->getNbParticles();
if(targets == sources){
inc -= 1;
}
for(int idx = 0 ; idx < 27 ; ++idx){
if( directNeighborsParticles[idx] ){
for(int idx = 0 ; idx < inSize ; ++idx){
inc += directNeighborsParticles[idx]->getNbParticles();
}
}
long long int*const particlesAttributes = targets->getDataDown();
......@@ -117,13 +114,11 @@ public:
/** After Downward */
void P2PRemote(const FTreeCoordinate& ,
ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict /*sources*/,
ContainerClass* const directNeighborsParticles[27], const int ){
ContainerClass* const directNeighborsParticles[], const int /*positions*/[], const int inSize) override {
// Each particles targeted is impacted by the particles sources
long long int inc = 0;
for(int idx = 0 ; idx < 27 ; ++idx){
if( directNeighborsParticles[idx] ){
for(int idx = 0 ; idx < inSize ; ++idx){
inc += directNeighborsParticles[idx]->getNbParticles();
}
}
long long int*const particlesAttributes = targets->getDataDown();
......
This diff is collapsed.
......@@ -190,7 +190,8 @@ protected:
typename OctreeClass::Iterator avoidGotoLeftIterator(octreeIterator);
const CellClass* neighbors[343];
const CellClass* neighbors[342];
int neighborPositions[342];
// for each levels
for(int idxLevel = FAbstractAlgorithm::upperWorkingLevel ; idxLevel < FAbstractAlgorithm::lowerWorkingLevel ; ++idxLevel ){
......@@ -200,9 +201,9 @@ protected:
// for each cells
do{
const int counter = tree->getInteractionNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(), idxLevel, separationCriteria);
const int counter = tree->getInteractionNeighbors(neighbors, neighborPositions, octreeIterator.getCurrentGlobalCoordinate(), idxLevel, separationCriteria);
FLOG(computationCounter.tic());
if(counter) kernels->M2L( octreeIterator.getCurrentCell() , neighbors, counter, idxLevel);
if(counter) kernels->M2L( octreeIterator.getCurrentCell() , neighbors, neighborPositions, counter, idxLevel);
FLOG(computationCounter.tac());
} while(octreeIterator.moveRight());
......@@ -282,7 +283,8 @@ protected:
typename OctreeClass::Iterator octreeIterator(tree);
octreeIterator.gotoBottomLeft();
// There is a maximum of 26 neighbors
ContainerClass* neighbors[27];
ContainerClass* neighbors[26];
int neighborPositions[26];
// for each leafs
do{
if(l2pEnabled){
......@@ -292,10 +294,10 @@ protected:
}
if(p2pEnabled){
// need the current particles and neighbors particles
const int counter = tree->getLeafsNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(),heightMinusOne);
const int counter = tree->getLeafsNeighbors(neighbors, neighborPositions, octreeIterator.getCurrentGlobalCoordinate(),heightMinusOne);
FLOG(computationCounterP2P.tic());
kernels->P2P(octreeIterator.getCurrentGlobalCoordinate(),octreeIterator.getCurrentListTargets(),
octreeIterator.getCurrentListSrc(), neighbors, counter);
octreeIterator.getCurrentListSrc(), neighbors, neighborPositions, counter);
FLOG(computationCounterP2P.tac());
}
} while(octreeIterator.moveRight());
......
This diff is collapsed.
......@@ -247,7 +247,8 @@ protected:
typename OctreeClass::Iterator octreeIterator(tree);
typename OctreeClass::Iterator avoidGotoLeftIterator(octreeIterator);
const CellClass* neighbors[343];
const CellClass* neighbors[342];
int neighborPositions[342];
// for each levels
for(int idxLevel = 1 ; idxLevel < OctreeHeight ; ++idxLevel ){
......@@ -256,9 +257,9 @@ protected:
const int separationCriteria = (idxLevel != OctreeHeight-1 ? 1 : leafLevelSeperationCriteria);
// for each cells
do{
const int counter = tree->getPeriodicInteractionNeighbors(neighbors, octreeIterator.getCurrentGlobalCoordinate(), idxLevel, AllDirs, separationCriteria);
const int counter = tree->getPeriodicInteractionNeighbors(neighbors, neighborPositions, octreeIterator.getCurrentGlobalCoordinate(), idxLevel, AllDirs, separationCriteria);
FLOG(computationCounter.tic());
if(counter) kernels->M2L( octreeIterator.getCurrentCell() , neighbors, counter, fackLevel);
if(counter) kernels->M2L( octreeIterator.getCurrentCell() , neighbors, neighborPositions, counter, fackLevel);
FLOG(computationCounter.tac());
} while(octreeIterator.moveRight());
avoidGotoLeftIterator.moveDown();
......@@ -327,8 +328,9 @@ protected:
typename OctreeClass::Iterator octreeIterator(tree);
octreeIterator.gotoBottomLeft();
// There is a maximum of 26 neighbors
ContainerClass* neighbors[27];
FTreeCoordinate offsets[27];
ContainerClass* neighbors[26];
FTreeCoordinate offsets[26];
int neighborPositions[26];
bool hasPeriodicLeaves;
// for each leafs
do{
......@@ -340,18 +342,18 @@ protected:
if(p2pEnabled){
// need the current particles and neighbors particles
const FTreeCoordinate centerOfLeaf = octreeIterator.getCurrentGlobalCoordinate();
const int counter = tree->getPeriodicLeafsNeighbors( neighbors, offsets, &hasPeriodicLeaves, centerOfLeaf, heightMinusOne, AllDirs);
const int counter = tree->getPeriodicLeafsNeighbors( neighbors, neighborPositions, offsets, &hasPeriodicLeaves, centerOfLeaf, heightMinusOne, AllDirs);
int periodicNeighborsCounter = 0;
if(hasPeriodicLeaves){
ContainerClass* periodicNeighbors[27];
memset(periodicNeighbors, 0, 27 * sizeof(ContainerClass*));
ContainerClass* periodicNeighbors[26];
int periodicNeighborPositions[26];
for(int idxNeig = 0 ; idxNeig < 27 ; ++idxNeig){
if( neighbors[idxNeig] && !offsets[idxNeig].equals(0,0,0) ){
for(int idxNeig = 0 ; idxNeig < counter ; ++idxNeig){
if( !offsets[idxNeig].equals(0,0,0) ){
// Put periodic neighbors into other array
periodicNeighbors[idxNeig] = neighbors[idxNeig];
neighbors[idxNeig] = nullptr;
periodicNeighbors[periodicNeighborsCounter] = neighbors[idxNeig];
periodicNeighborPositions[periodicNeighborsCounter] = neighborPositions[idxNeig];
++periodicNeighborsCounter;
FReal*const positionsX = periodicNeighbors[idxNeig]->getWPositions()[0];
......@@ -364,31 +366,33 @@ protected:
positionsZ[idxPart] += boxWidth * FReal(offsets[idxNeig].getZ());
}
}
else{
neighbors[idxNeig-periodicNeighborsCounter] = neighbors[idxNeig];
neighborPositions[idxNeig-periodicNeighborsCounter] = neighborPositions[idxNeig];
}
}
FLOG(computationCounterP2P.tic());
kernels->P2PRemote(octreeIterator.getCurrentGlobalCoordinate(),octreeIterator.getCurrentListTargets(),
octreeIterator.getCurrentListSrc(), periodicNeighbors, periodicNeighborsCounter);
octreeIterator.getCurrentListSrc(), periodicNeighbors, periodicNeighborPositions, periodicNeighborsCounter);
FLOG(computationCounterP2P.tac());
for(int idxNeig = 0 ; idxNeig < 27 ; ++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];
for(int idxNeig = 0 ; idxNeig < periodicNeighborsCounter ; ++idxNeig){
FReal*const positionsX = periodicNeighbors[idxNeig]->getWPositions()[0];
FReal*const positionsY = periodicNeighbors[idxNeig]->getWPositions()[1];
FReal*const positionsZ = periodicNeighbors[idxNeig]->getWPositions()[2];
for(FSize idxPart = 0; idxPart < periodicNeighbors[idxNeig]->getNbParticles() ; ++idxPart){
positionsX[idxPart] -= boxWidth * FReal(offsets[idxNeig].getX());
positionsY[idxPart] -= boxWidth * FReal(offsets[idxNeig].getY());
positionsZ[idxPart] -= boxWidth * FReal(offsets[idxNeig].getZ());
}
for(FSize idxPart = 0; idxPart < periodicNeighbors[idxNeig]->getNbParticles() ; ++idxPart){
positionsX[idxPart] -= boxWidth * FReal(offsets[idxNeig].getX());
positionsY[idxPart] -= boxWidth * FReal(offsets[idxNeig].getY());
positionsZ[idxPart] -= boxWidth * FReal(offsets[idxNeig].getZ());
}
}
}
FLOG(computationCounterP2P.tic());
kernels->P2P(octreeIterator.getCurrentGlobalCoordinate(),octreeIterator.getCurrentListTargets(),
octreeIterator.getCurrentListSrc(), neighbors, counter - periodicNeighborsCounter);
octreeIterator.getCurrentListSrc(), neighbors, neighborPositions, counter - periodicNeighborsCounter);
FLOG(computationCounterP2P.tac());
}
} while(octreeIterator.moveRight());
......@@ -450,32 +454,34 @@ protected:
{
const int idxUpperLevel = 2;
const CellClass* neighbors[343];
memset(neighbors, 0, sizeof(CellClass*) * 343);
const CellClass* neighbors[342];
int neighborPositions[342];
int counter = 0;
for(int idxX = -3 ; idxX <= 2 ; ++idxX){
for(int idxY = -3 ; idxY <= 2 ; ++idxY){
for(int idxZ = -3 ; idxZ <= 2 ; ++idxZ){
if( FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
neighbors[neighIndex(idxX,idxY,idxZ)] = &upperCells[idxUpperLevel-1];
neighbors[counter] = &upperCells[idxUpperLevel-1];
neighborPositions[counter] = neighIndex(idxX,idxY,idxZ);
++counter;
}
}
}
}
// compute M2L
kernels->M2L( &downerCells[idxUpperLevel-1] , neighbors, counter, idxUpperLevel);
kernels->M2L( &downerCells[idxUpperLevel-1] , neighbors, neighborPositions, counter, idxUpperLevel);
}
for(int idxUpperLevel = 3 ; idxUpperLevel <= offsetRealTree ; ++idxUpperLevel){
const CellClass* neighbors[343];
memset(neighbors, 0, sizeof(CellClass*) * 343);
const CellClass* neighbors[342];
int neighborPositions[342];
int counter = 0;
for(int idxX = -2 ; idxX <= 3 ; ++idxX){
for(int idxY = -2 ; idxY <= 3 ; ++idxY){
for(int idxZ = -2 ; idxZ <= 3 ; ++idxZ){
if( FMath::Abs(idxX) > 1 || FMath::Abs(idxY) > 1 || FMath::Abs(idxZ) > 1){
neighbors[neighIndex(idxX,idxY,idxZ)] = &upperCells[idxUpperLevel-1];
neighbors[counter] = &upperCells[idxUpperLevel-1];
neighborPositions[counter] = neighIndex(idxX,idxY,idxZ);
++counter;
}