Commit fa137a72 authored by BRAMAS Berenger's avatar BRAMAS Berenger
Browse files
parents 51a1acac ea61f86d
......@@ -380,11 +380,20 @@ typedef void (*Callback_M2M)(int level, void* parentCell, int childPosition, voi
* @param targetCell pointer to cell to be filled
* @param sourceCellPosition number of source cell (source cell
* position in the tree can inferred from its number (refer to doc))
* @param sourceCell array of cell to be read
* @param sourceCell cell to be read
* @param userData datas specific to the user's kernel
*/
typedef void (*Callback_M2L)(int level, void* targetCell, int sourceCellPosition, void* sourceCell, void* userData);
/**
* @brief Function to be filled by user's M2L
* @param level current level in the tree
* @param targetCell pointer to cell to be filled
* @param sourceCell array of cell to be read
* @param userData datas specific to the user's kernel
*/
typedef void (*Callback_M2LFull)(int level, void* targetCell, void* sourceCell[343], void* userData);
/**
* @brief Function to be filled by user's L2L
* @param level current level in the tree
......@@ -414,6 +423,19 @@ typedef void (*Callback_L2P)(void* leafCell, int nbParticles,const int* particle
*/
typedef void (*Callback_P2P)(int nbParticles, const int* particleIndexes, int nbSourceParticles, const int* sourceParticleIndexes, void* userData);
/**
* @brief Function to be filled by user's P2P
* @param nbParticles number of particle in current leaf
* @param particleIndexes indexes of particles currently computed
* @param sourceCell array of the neighbors source cells
* @param sourceParticleIndexes array of indices of source particles currently computed
* @param sourceNbPart array containing the number of part in each neighbors
* @param userData datas specific to the user's kernel
*/
typedef void (*Callback_P2PFull)(int nbParticles, const int* particleIndexes,
const int * sourceParticleIndexes[27],int sourceNbPart[27], void* userData);
/**
* @brief Function to be filled by user's P2P inside the leaf
* @param nbParticles number of particle in current leaf
......@@ -433,9 +455,11 @@ typedef struct User_Scalfmm_Kernel_Descriptor {
Callback_P2M p2m;
Callback_M2M m2m;
Callback_M2L m2l;
Callback_M2LFull m2l_full;
Callback_L2L l2l;
Callback_L2P l2p;
Callback_P2P p2p;
Callback_P2PFull p2p_full;
Callback_P2PInner p2pinner;
}Scalfmm_Kernel_Descriptor;
......
......@@ -187,8 +187,6 @@ public:
forcesToFill[indexes[idxPart]*3+0] = sources->getForcesX()[idxPart];
forcesToFill[indexes[idxPart]*3+1] = sources->getForcesY()[idxPart];
forcesToFill[indexes[idxPart]*3+2] = sources->getForcesZ()[idxPart];
printf("forces found : %e,%e,%e\n",
sources->getForcesX()[idxPart],sources->getForcesY()[idxPart],sources->getForcesZ()[idxPart]);
}
});
}
......
......@@ -31,10 +31,12 @@
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
//For lagrange interpolation
#include "Kernels/Uniform/FUnifCell.hpp"
//For interpolation
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Uniform/FUnifKernel.hpp"
//For lagrange interpolation
// #include "Kernels/Uniform/FUnifCell.hpp"
// #include "Kernels/Uniform/FUnifKernel.hpp"
//For chebyshev Interpolation
#include "Kernels/Chebyshev/FChebCell.hpp"
......
......@@ -24,16 +24,16 @@ extern "C" scalfmm_handle scalfmm_init(/*int TreeHeight,double BoxWidth,double*
handle->engine = new FInterEngine<ChebCell,ChebKernel>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
break;
case 2:
//TODO typedefs
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FUnifCell<7> UnifCell;
// case 2:
// //TODO typedefs
// typedef FP2PParticleContainerIndexed<> ContainerClass;
// typedef FUnifCell<7> UnifCell;
typedef FInterpMatrixKernelR MatrixKernelClass;
typedef FUnifKernel<UnifCell,ContainerClass,MatrixKernelClass,7> UnifKernel;
// typedef FInterpMatrixKernelR MatrixKernelClass;
// typedef FUnifKernel<UnifCell,ContainerClass,MatrixKernelClass,7> UnifKernel;
handle->engine = new FInterEngine<UnifCell,UnifKernel>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
break;
// handle->engine = new FInterEngine<UnifCell,UnifKernel>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
// break;
default:
std::cout<< "Kernel type unsupported" << std::endl;
......@@ -48,3 +48,263 @@ extern "C" void scalfmm_dealloc_handle(scalfmm_handle handle, Callback_free_cell
delete ((ScalFmmCoreHandle *) handle)->engine ;
delete (ScalFmmCoreHandle *) handle;
}
/**
* This parts implements all the function defined in FChebInterface.h
* using the Chebyshev classes
*/
#ifndef CHEBINTERFACE_HPP
#define CHEBINTERFACE_HPP
extern "C" {
#include "Kernels/Chebyshev/FChebInterface.h"
}
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
#include "Components/FSimpleLeaf.hpp"
typedef struct FChebCell_struct{
//Store what's needed
FChebCell<7> * cell;
long long int index;
}ChebCellStruct;
//How to create/destroy cells
extern "C" ChebCellStruct * ChebCellStruct_create(long long int inIndex,int * position){
ChebCellStruct * newCell = new ChebCellStruct();
newCell->index = inIndex;
newCell->cell = new FChebCell<7>();
newCell->cell->setMortonIndex(inIndex);
newCell->cell->setCoordinate(position[0],position[1],position[2]);
return newCell;
}
extern "C" void ChebCellStruct_free(ChebCellStruct * inCell){
if(inCell->cell) {
delete inCell->cell;
}
delete inCell;
}
typedef struct FChebKernel_struct{
//Be ready full duplication go there !!!
FChebSymKernel<FChebCell<7>,FP2PParticleContainerIndexed<>,FInterpMatrixKernelR,7> * kernel;
FInterpMatrixKernelR * matrix;
} ChebKernelStruct;
//Kernel functions
extern "C" ChebKernelStruct * ChebKernelStruct_create(int inTreeHeight,
double inBoxWidth,
double* inBoxCenter){
ChebKernelStruct * newKernel = new ChebKernelStruct();
newKernel->matrix= new FInterpMatrixKernelR();
newKernel->kernel =
new FChebSymKernel<FChebCell<7>,FP2PParticleContainerIndexed<>,
FInterpMatrixKernelR,7>(inTreeHeight,
inBoxWidth,
FPoint(inBoxCenter[0],inBoxCenter[1],inBoxCenter[2]),
newKernel->matrix);
return newKernel;
}
extern "C" void ChebKernelStruct_free(void *inKernel){
delete reinterpret_cast<ChebKernelStruct *>(inKernel)->matrix;
delete reinterpret_cast<ChebKernelStruct *>(inKernel)->kernel;
delete reinterpret_cast<ChebKernelStruct *>(inKernel);
}
extern "C" void ChebKernel_P2M(void * leafCell, int nbParticles, const int *particleIndexes, void * inKernel){
//make temporary array of parts
FP2PParticleContainerIndexed<>* tempContainer = new FP2PParticleContainerIndexed<>();
tempContainer->reserve(nbParticles);
FPoint pos;
for(int i=0 ; i<nbParticles ; ++i){
pos = FPoint(reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3 ],
reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3+1],
reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3+2]);
double Phi = reinterpret_cast<UserData *>(inKernel)->myPhyValues[particleIndexes[i]];
tempContainer->push(pos,particleIndexes[i],Phi);
}
//get the real cell struct
ChebCellStruct * realCellStruct = reinterpret_cast<ChebCellStruct *>(leafCell);
FChebCell<7> * realCell = realCellStruct->cell;
//get the real chebyshev struct
UserData * userDataKernel = reinterpret_cast<UserData *>(inKernel);
ChebKernelStruct * realKernel = userDataKernel->kernelStruct;
realKernel->kernel->P2M(realCell, tempContainer);
delete tempContainer;
}
extern "C" void ChebKernel_M2M(int level, void* parentCell, int childPosition, void *childCell, void *inKernel){
//Get our structures
ChebCellStruct * parentCellStruct = reinterpret_cast<ChebCellStruct *>(parentCell);
ChebCellStruct * childCellStruct = reinterpret_cast<ChebCellStruct *>(childCell);
//get real cheb cell
FChebCell<7>* parentChebCell = parentCellStruct->cell;
FChebCell<7>* childChebCell = childCellStruct->cell;
//Get the kernel
ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;
inKernelStruct->kernel->getPtrToInterpolator()->applyM2M(childPosition,
childChebCell->getMultipole(0),
parentChebCell->getMultipole(0));
}
extern "C" void ChebKernel_M2L(int level, void* targetCell, void* sourceCell[343], void* inKernel){
//Get our structures
ChebCellStruct * targetCellStruct = reinterpret_cast<ChebCellStruct *>(targetCell);
//get real cheb cell
FChebCell<7>* const targetChebCell = targetCellStruct->cell;
//copy to an array of FChebCell
const FChebCell<7>* arrayOfChebCell[343];
for(int i=0; i<343 ; ++i){
if(sourceCell[i] != nullptr){
arrayOfChebCell[i] = reinterpret_cast<ChebCellStruct*>(sourceCell[i])->cell;
}
else{
arrayOfChebCell[i] = nullptr;
}
}
//Get the kernel
ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;
inKernelStruct->kernel->M2L(targetChebCell,arrayOfChebCell,0,level);
}
extern "C" void ChebKernel_L2L(int level, void* parentCell, int childPosition, void* childCell, void* inKernel){
//Get our structures
ChebCellStruct * parentCellStruct = reinterpret_cast<ChebCellStruct *>(parentCell);
ChebCellStruct * childCellStruct = reinterpret_cast<ChebCellStruct *>(childCell);
//get real cheb cell
FChebCell<7>* parentChebCell = parentCellStruct->cell;
FChebCell<7>* childChebCell = childCellStruct->cell;
//Get the kernel
ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;
inKernelStruct->kernel->getPtrToInterpolator()->applyL2L(childPosition,
parentChebCell->getLocal(0),
childChebCell->getLocal(0));
}
extern "C" void ChebKernel_L2P(void* leafCell, int nbParticles, const int* particleIndexes, void* inKernel){
//Create temporary FSimpleLeaf
FP2PParticleContainerIndexed<>* tempContainer = new FP2PParticleContainerIndexed<>();
tempContainer->reserve(nbParticles);
FPoint pos;
for(int i=0 ; i<nbParticles ; ++i){
pos = FPoint(reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3 ],
reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3+1],
reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3+2]);
double Phi = reinterpret_cast<UserData *>(inKernel)->myPhyValues[particleIndexes[i]];
tempContainer->push(pos,particleIndexes[i],Phi);
}
//Get our structures
ChebCellStruct * leafCellStruct = reinterpret_cast<ChebCellStruct *>(leafCell);
//get real cheb cell
FChebCell<7>* leafChebCell = leafCellStruct->cell;
//Get the kernel
ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;
inKernelStruct->kernel->L2P(leafChebCell,tempContainer);
//Then retrieve the results
double * forcesToFill = reinterpret_cast<UserData *>(inKernel)->forcesComputed;
const FVector<int>& indexes = tempContainer->getIndexes();
for(int idxPart = 0 ; idxPart<nbParticles ; ++idxPart){
forcesToFill[indexes[idxPart]*3+0] += tempContainer->getForcesX()[idxPart];
forcesToFill[indexes[idxPart]*3+1] += tempContainer->getForcesY()[idxPart];
forcesToFill[indexes[idxPart]*3+2] += tempContainer->getForcesZ()[idxPart];
}
delete tempContainer;
tempContainer=nullptr;
}
void ChebKernel_P2P(int nbParticles, const int* particleIndexes,
const int * sourceParticleIndexes[27],int sourceNbPart[27],void* inKernel){
//Create temporary FSimpleLeaf for target
FP2PParticleContainerIndexed<>* tempContTarget = new FP2PParticleContainerIndexed<>();
tempContTarget->reserve(nbParticles);
for(int i=0 ; i<nbParticles ; ++i){
FPoint pos = FPoint(reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3 ],
reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3+1],
reinterpret_cast<UserData *>(inKernel)->insertedPositions[particleIndexes[i]*3+2]);
double Phi = reinterpret_cast<UserData *>(inKernel)->myPhyValues[particleIndexes[i]];
tempContTarget->push(pos,particleIndexes[i],Phi);
}
//Create 27 FSimpleLeaf for 27 sources
FP2PParticleContainerIndexed<>* tempContSources[27];
for(int idSource=0; idSource<27 ; ++idSource){
if(sourceNbPart[idSource] != 0){
//Create container
tempContSources[idSource] = new FP2PParticleContainerIndexed<>();
//Allocate memory
tempContSources[idSource]->reserve(sourceNbPart[idSource]);
//Store a ptr to the indices of that source leaf
const int * indSource = sourceParticleIndexes[idSource];
//Then, for each part in this source
for(int i=0 ; i<sourceNbPart[idSource] ; ++i){
FPoint pos = FPoint(reinterpret_cast<UserData *>(inKernel)->insertedPositions[indSource[i]*3 ],
reinterpret_cast<UserData *>(inKernel)->insertedPositions[indSource[i]*3+1],
reinterpret_cast<UserData *>(inKernel)->insertedPositions[indSource[i]*3+2]);
double Phi = reinterpret_cast<UserData *>(inKernel)->myPhyValues[indSource[i]];
tempContSources[idSource]->push(pos,indSource[i],Phi);
}
}
else{
tempContSources[idSource] = nullptr;
}
}
//Everything is fine, now, call Chebyshev P2P
//Get the kernel
ChebKernelStruct * inKernelStruct = reinterpret_cast<UserData*>(inKernel)->kernelStruct;
//Empty tree coordinate
int coord[3] = {0,0,0};
inKernelStruct->kernel->P2P(FTreeCoordinate(coord),tempContTarget,nullptr,tempContSources,0);
//get back forces
double * forcesToFill = reinterpret_cast<UserData *>(inKernel)->forcesComputed;
const FVector<int>& indexes = tempContTarget->getIndexes();
for(int idxPart = 0 ; idxPart<nbParticles ; ++idxPart){
forcesToFill[indexes[idxPart]*3+0] += tempContTarget->getForcesX()[idxPart];
forcesToFill[indexes[idxPart]*3+1] += tempContTarget->getForcesY()[idxPart];
forcesToFill[indexes[idxPart]*3+2] += tempContTarget->getForcesZ()[idxPart];
}
//Note that sources are also modified.
//get back sources forces
for(int idSource = 0 ; idSource < 27 ; ++idSource){
const FVector<int>& indexes = tempContSources[idSource]->getIndexes();
const int nbPartInSource = sourceNbPart[idSource];
for(int idxSourcePart = 0; idxSourcePart < nbPartInSource ; ++idxSourcePart){
forcesToFill[indexes[idxSourcePart]*3+0] += tempContSources[idSource]->getForcesX()[idxSourcePart];
forcesToFill[indexes[idxSourcePart]*3+1] += tempContSources[idSource]->getForcesY()[idxSourcePart];
forcesToFill[indexes[idxSourcePart]*3+2] += tempContSources[idSource]->getForcesZ()[idxSourcePart];
}
}
//Release memory
for(int idSource=0; idSource<27 ; ++idSource){
if(tempContSources[idSource]) delete tempContSources[idSource];
}
delete tempContTarget;
}
#endif
......@@ -72,6 +72,7 @@ public:
void* getContainer() const {
return userData;
}
};
/**
......@@ -115,10 +116,25 @@ public:
/** Do nothing */
virtual void M2L(CellClass* const FRestrict cell, const CellClass* interactions[], const int , const int level) {
if(kernel.m2l){
for(int idx = 0 ; idx < 343 ; ++idx){
if( interactions[idx] ){
kernel.m2l(level, cell->getContainer(), idx, interactions[idx]->getContainer(), userData);
if(kernel.m2l_full){//all 343 interactions will be computed directly
//First, copy the fmm cell inside an array of user cells
void * userCellArray[343];
for(int i=0 ; i<343 ; ++i){
if(interactions[i] != nullptr){
userCellArray[i] = interactions[i]->getContainer();
}
else{
userCellArray[i] = nullptr;
}
}
kernel.m2l_full(level,cell->getContainer(),userCellArray,userData);
}
else{
if(kernel.m2l){
for(int idx = 0 ; idx < 343 ; ++idx){
if( interactions[idx] ){
kernel.m2l(level, cell->getContainer(), idx, interactions[idx]->getContainer(), userData);
}
}
}
}
......@@ -147,11 +163,27 @@ public:
ContainerClass* const neighbors[27], const int ){
if(kernel.p2pinner) kernel.p2pinner(targets->getNbParticles(), targets->getIndexes().data(), userData);
if(kernel.p2p_full){
//Create the arrays of size and indexes
int nbPartPerNeighbors[27];
const int * indicesPerNeighbors[27];
for(int idx=0 ; idx<27 ; ++idx){
if(neighbors[idx]){
nbPartPerNeighbors[idx] = neighbors[idx]->getNbParticles();
indicesPerNeighbors[idx] = neighbors[idx]->getIndexes().data();
}
else{
nbPartPerNeighbors[idx] = 0;
indicesPerNeighbors[idx] = nullptr;
}
}
kernel.p2p_full(targets->getNbParticles(),targets->getIndexes().data(),indicesPerNeighbors,nbPartPerNeighbors,userData);
}
if(kernel.p2p){
for(int idx = 0 ; idx < 27 ; ++idx){
if( neighbors[idx] ){
kernel.p2p(targets->getNbParticles(), targets->getIndexes().data(),
neighbors[idx]->getNbParticles(), neighbors[idx]->getIndexes().data(), userData);
neighbors[idx]->getNbParticles(), neighbors[idx]->getIndexes().data(), userData);
}
}
}
......@@ -200,12 +232,14 @@ public:
~FUserKernelEngine(){
delete octree;
octree=nullptr;
if(arranger){
delete arranger;
arranger=nullptr;
}
if(kernel){
delete kernel;
kernel=nullptr;
}
}
......@@ -233,6 +267,7 @@ public:
octree->insert(FPoint(&XYZ[3*idPart]),idPart);
}
nbPart += NbPositions;
this->init_cell();
}
/**
......@@ -389,8 +424,6 @@ public:
}
}
/*
* Call the user allocator on userDatas member field of each cell
*/
......@@ -422,6 +455,7 @@ public:
octree->forEachCell([&](CoreCell * currCell){
if(currCell->getContainer()){
user_cell_deallocator(currCell->getContainer());
currCell->setContainer(nullptr);
}
});
}
......
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
//For timing monitoring
#include <time.h>
#include <sys/time.h>
#include "../Src/CScalfmmApi.h"
#include "../../Src/Kernels/Chebyshev/FChebInterface.h"
/**
* This file is an example of the user defined kernel API, with link
* to our Chebyshev Kernel
**/
/**
* @brief Wrapper to init internal ChebCell
*/
void* cheb_init_cell(int level, long long morton_index, int* tree_position, double* spatial_position){
return ChebCellStruct_create(morton_index,tree_position);
}
/**
* @brief Wrapper to free internal ChebCell
*/
void cheb_free_cell(void * inCell){
ChebCellStruct_free(inCell);
}
/**
* @brief Wrapper to FMM operators (refer to CScalfmmApi.h to get the
* detailed descriptions)
*/
void cheb_p2m(void* cellData, int nbParticlesInLeaf, const int* particleIndexes, void* userData){
ChebKernel_P2M(cellData,nbParticlesInLeaf,particleIndexes,userData);
}
void cheb_m2m(int level, void* parentCell, int childPosition, void* childCell, void* userData){
ChebKernel_M2M(level,parentCell,childPosition,childCell,userData);
}
void cheb_m2l_full(int level, void* targetCell, void* sourceCell[343], void* userData){
ChebKernel_M2L(level, targetCell, sourceCell, userData);
}
void cheb_l2l(int level, void* parentCell, int childPosition, void* childCell, void* userData){
ChebKernel_L2L( level, parentCell, childPosition, childCell, userData);
}
void cheb_l2p(void* leafCell, int nbParticles, const int* particleIndexes, void* userData){
ChebKernel_L2P( leafCell, nbParticles, particleIndexes, userData);
}
void cheb_p2pFull(int nbParticles, const int* particleIndexes,
const int * sourceParticleIndexes[27], int sourceNbPart[27],void* userData) {
ChebKernel_P2P(nbParticles, particleIndexes, sourceParticleIndexes, sourceNbPart, userData);
}
/**
* @brief Wrapper on timeval struct
*/
typedef struct timer{
struct timeval start;
struct timeval end;
}Timer;
/**
* @brief monitoring function (init start)
*/
void tic(Timer* inTimer){
gettimeofday(&inTimer->start, NULL);
}
/**
* @brief monitoring function (init end)
*/
void tac(Timer* inTimer){
gettimeofday(&inTimer->end, NULL);
}
/**
* @brief monitoring function (return elapsed time in micro second)
*/
long int get_elapsed(Timer* inTimer){
return ((inTimer->end.tv_sec * 1000000 + inTimer->end.tv_usec)
- (inTimer->start.tv_sec * 1000000 + inTimer->start.tv_usec));
}
/**
* @brief monitoring function (print elapsed)
*/
void print_elapsed(Timer* inTimer){
long int elapsed = get_elapsed(inTimer);
printf("Elapsed : %ld us (or %f seconds)\n",elapsed,elapsed/1000000.0);
}
/**
* @brief monitoring function (print difference :: First-Second)
*/
void print_difference_elapsed(Timer* inTimer1,Timer*inTimer2){
long int diff = get_elapsed(inTimer1)-get_elapsed(inTimer2);
printf("Timer Difference : %ld us (or %f seconds)\n",diff,(double)diff/1000000.0);
}
/**
* @brief Do everything
* @param number of particle (no default value)
*/
int main(int argc, char ** av){
omp_set_num_threads(1);
printf("Start\n");
if(argc<2){
printf("Use : %s nb_part\n exiting\n",av[0]);
exit(0);
}
int nbPart= atoi(av[1]);
double* particleXYZ =