Commit 5c67f64e authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille

Api is almost finished, new test added in order to test the user kernel interface

parent 3d3be250
......@@ -362,7 +362,7 @@ typedef void (*Callback_L2L)(int level, void* parentCell, int childPosition, voi
* @param particleIndexes indexes of particles currently computed
* @param userData datas specific to the user's kernel
*/
typedef void (*Callback_L2P)(void* leafCell, int nbParticles, int* particleIndexes, void* userData);
typedef void (*Callback_L2P)(void* leafCell, int nbParticles,const int* particleIndexes, void* userData);
/**
* @brief Function to be filled by user's P2P
......@@ -380,7 +380,7 @@ typedef void (*Callback_P2P)(int nbParticles, const int* particleIndexes, int nb
* @param particleIndexes indexes of particles currently computed
* @param userData datas specific to the user's kernel
*/
typedef void (*Callback_P2PInner)(int nbParticles, int* particleIndexes, void* userData);
typedef void (*Callback_P2PInner)(int nbParticles, const int* particleIndexes, void* userData);
......@@ -434,6 +434,10 @@ typedef void* (*Callback_init_cell)(int level, long long morton_index, int* tree
typedef void (*Callback_free_cell)(void*);
void scalfmm_init_cell(scalfmm_handle Handle, Callback_init_cell user_cell_initializer);
void scalfmm_free_cell(scalfmm_handle Handle, Callback_free_cell user_cell_deallocator);
///////////////// Common kernel part : /////////////////
......
......@@ -67,11 +67,12 @@ public:
* @param BoxCenter double[3] coordinate of the center of the
* simulation box
*/
FInterEngine(int TreeHeight, double BoxWidth , double * BoxCenter) : kernel(nullptr), matrix(nullptr), octree(nullptr),arranger(nullptr){
FInterEngine(int TreeHeight, double BoxWidth , double * BoxCenter,scalfmm_kernel_type KernelType) :
kernel(nullptr), matrix(nullptr), octree(nullptr),arranger(nullptr){
octree = new OctreeClass(TreeHeight,FMath::Min(3,TreeHeight-1),BoxWidth,FPoint(BoxCenter));
this->matrix = new MatrixKernelClass();
this->kernel = new InterKernel(TreeHeight,BoxWidth,FPoint(BoxCenter),matrix);
kernelType = chebyshev;
kernelType = KernelType;
arranger = nullptr;
}
......
......@@ -56,7 +56,6 @@ protected:
public:
FScalFMMEngine() : Algorithm(multi_thread), progress(), nbPart(0){
//Default behavior in case of out of the box particles is exiting
}
//First function displayed there are common function for every
......@@ -197,6 +196,14 @@ public:
FAssertLF(0,"No user kernel defined, exiting ...\n");
}
virtual void init_cell(Callback_init_cell user_cell_initializer){
FAssertLF(0,"No user kernel defined, exiting ...\n");
}
virtual void free_cell(Callback_free_cell user_cell_initializer){
FAssertLF(0,"No user kernel defined, exiting ...\n");
}
virtual void execute_fmm(){
FAssertLF(0,"No kernel set, cannot execute anything, exiting ...\n");
}
......@@ -371,6 +378,13 @@ extern "C" void scalfmm_execute_fmm(scalfmm_handle Handle){
((ScalFmmCoreHandle * ) Handle)->engine->execute_fmm();
}
extern "C" void scalfmm_init_cell(scalfmm_handle Handle, Callback_init_cell user_cell_initializer){
((ScalFmmCoreHandle * ) Handle)->engine->init_cell(user_cell_initializer);
}
extern "C" void scalfmm_free_cell(scalfmm_handle Handle, Callback_free_cell user_cell_deallocator){
((ScalFmmCoreHandle * ) Handle)->engine->free_cell(user_cell_deallocator);
}
#endif
......@@ -8,6 +8,7 @@ extern "C" {
}
#include "FInterEngine.hpp"
#include "FUserKernelEngine.hpp"
extern "C" scalfmm_handle scalfmm_init(int TreeHeight,double BoxWidth,double* BoxCenter, scalfmm_kernel_type KernelType){
ScalFmmCoreHandle * handle = new ScalFmmCoreHandle();
......@@ -26,7 +27,7 @@ extern "C" scalfmm_handle scalfmm_init(int TreeHeight,double BoxWidth,double* Bo
typedef FInterpMatrixKernelR MatrixKernelClass;
typedef FChebSymKernel<ChebCell,ContainerClass,MatrixKernelClass,7> ChebKernel;
handle->engine = new FInterEngine<ChebCell,ChebKernel>(TreeHeight,BoxWidth,BoxCenter);
handle->engine = new FInterEngine<ChebCell,ChebKernel>(TreeHeight,BoxWidth,BoxCenter, KernelType);
break;
case 2:
//TODO typedefs
......@@ -36,7 +37,7 @@ extern "C" scalfmm_handle scalfmm_init(int TreeHeight,double BoxWidth,double* Bo
typedef FInterpMatrixKernelR MatrixKernelClass;
typedef FUnifKernel<UnifCell,ContainerClass,MatrixKernelClass,7> UnifKernel;
handle->engine = new FInterEngine<UnifCell,UnifKernel>(TreeHeight,BoxWidth,BoxCenter);
handle->engine = new FInterEngine<UnifCell,UnifKernel>(TreeHeight,BoxWidth,BoxCenter,KernelType);
break;
default:
......
// ===================================================================================
// Copyright ScalFmm 2014 I
// 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".
// ===================================================================================
/**
* @file This file contains a class that inherits from FScalFMMEngine,
* and will implement the API functions for a user defined kernel.
*/
#ifndef FUSERKERNELENGINE_HPP
#define FUSERKERNELENGINE_HPP
#include "FScalFMMEngine.hpp"
/**
* @brief CoreCell : Cell used to store User datas
*/
class CoreCell : public FBasicCell {
// Mutable in order to work with the API
mutable void* userData;
public:
CoreCell() : userData(nullptr) {
}
~CoreCell(){
}
/**
* @brief setContainer store the ptr to the user data inside our
* struct
*/
void setContainer(void* inContainer) const {
userData = inContainer;
}
/**
* @brief getContainer : return the user datas (in order to give
* it back to the user defined kernel function)
*/
void* getContainer() const {
return userData;
}
};
/**
* This class simply call the function pointers from Scalfmm_Kernel_Descriptor.
* If not pointer is set the calls are skipped.
* The userData is given at any calls.
*/
template< class CellClass, class ContainerClass>
class CoreKernel : public FAbstractKernels<CellClass,ContainerClass> {
Scalfmm_Kernel_Descriptor kernel;
void* userData;
public:
CoreKernel(Scalfmm_Kernel_Descriptor inKernel, void* inUserData) : kernel(inKernel) , userData(inUserData){
}
/** Default destructor */
virtual ~CoreKernel(){
}
/** Do nothing */
virtual void P2M(CellClass* const cell, const ContainerClass* const container) {
if(kernel.p2m) kernel.p2m(cell->getContainer(), container->getNbParticles(), container->getIndexes().data(), userData);
}
/** Do nothing */
virtual void M2M(CellClass* const FRestrict cell, const CellClass*const FRestrict *const FRestrict children, const int level) {
if(kernel.m2m){
for(int idx = 0 ; idx < 8 ; ++idx){
if( children[idx] ){
kernel.m2m(level, cell->getContainer(), idx, children[idx]->getContainer(), userData);
}
}
}
}
/** 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);
}
}
}
}
/** Do nothing */
virtual void L2L(const CellClass* const FRestrict cell, CellClass* FRestrict *const FRestrict children, const int level) {
if(kernel.l2l){
for(int idx = 0 ; idx < 8 ; ++idx){
if( children[idx] ){
kernel.l2l(level, cell->getContainer(), idx, children[idx]->getContainer(), userData);
}
}
}
}
/** Do nothing */
virtual void L2P(const CellClass* const cell, ContainerClass* const container){
if(kernel.l2p) kernel.l2p(cell->getContainer(), container->getNbParticles(), container->getIndexes().data(), userData);
}
/** Do nothing */
virtual void P2P(const FTreeCoordinate& ,
ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict /*sources*/,
ContainerClass* const neighbors[27], const int ){
if(kernel.p2pinner) kernel.p2pinner(targets->getNbParticles(), targets->getIndexes().data(), 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);
}
}
}
}
/** Do nothing */
virtual void P2PRemote(const FTreeCoordinate& ,
ContainerClass* const FRestrict , const ContainerClass* const FRestrict ,
ContainerClass* const [27], const int ){
}
};
class FUserKernelEngine : public FScalFMMEngine{
private:
//Typedefs :
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf<ContainerClass> LeafClass;
typedef FOctree<CoreCell,ContainerClass,LeafClass> OctreeClass;
typedef CoreKernel<CoreCell,ContainerClass> CoreKernelClass;
typedef FP2PLeafInterface<OctreeClass> LeafInterface;
//For arranger classes
typedef FOctreeArranger<OctreeClass,ContainerClass,LeafInterface> ArrangerClass;
typedef FArrangerPeriodic<OctreeClass,ContainerClass,LeafInterface> ArrangerClassPeriodic;
//Attributes
OctreeClass * octree;
CoreKernelClass * kernel;
ArrangerClass * arranger;
public:
FUserKernelEngine(int TreeHeight, double BoxWidth , double * BoxCenter, scalfmm_kernel_type KernelType) :
octree(nullptr), kernel(nullptr), arranger(nullptr){
octree = new OctreeClass(TreeHeight,FMath::Min(3,TreeHeight-1),BoxWidth,FPoint(BoxCenter));
kernelType = KernelType;
//Kernel is not set now because the user must provide a
//Scalfmm_Kernel_descriptor
//kernel = new CoreKernelClass();
}
void user_kernel_config( Scalfmm_Kernel_Descriptor userKernel, void * userDatas){
kernel = new CoreKernelClass(userKernel,userDatas);
}
void tree_insert_particles( int NbPositions, double * arrayX, double * arrayY, double * arrayZ){
for(int idPart = 0; idPart<NbPositions ; ++idPart){
octree->insert(FPoint(arrayX[idPart],arrayY[idPart],arrayZ[idPart]),idPart);
}
nbPart += NbPositions;
}
void tree_insert_particles_xyz( int NbPositions, double * XYZ){
for(int idPart = 0; idPart<NbPositions ; ++idPart){
octree->insert(FPoint(&XYZ[3*idPart]),idPart);
}
nbPart += NbPositions;
}
/*
* Call the user allocator on userDatas member field of each cell
*/
void init_cell(Callback_init_cell user_cell_initializer){
if(user_cell_initializer){
double boxwidth = octree->getBoxWidth();
FPoint BoxCenter = octree->getBoxCenter();
double boxCorner[3];
boxCorner[0] = BoxCenter.getX() - boxwidth/2.0;
boxCorner[1] = BoxCenter.getY() - boxwidth/2.0;
boxCorner[2] = BoxCenter.getZ() - boxwidth/2.0;
//apply user function on each cell
octree->forEachCellWithLevel([&](CoreCell * currCell,const int currLevel){
FTreeCoordinate currCoord = currCell->getCoordinate();
int arrayCoord[3] = {currCoord.getX(),currCoord.getY(),currCoord.getZ()};
MortonIndex currMorton = currCoord.getMortonIndex(currLevel);
double position[3];
position[0] = boxCorner[0] + currCoord.getX()*boxwidth/double(1<<currLevel);
position[1] = boxCorner[1] + currCoord.getY()*boxwidth/double(1<<currLevel);
position[2] = boxCorner[2] + currCoord.getZ()*boxwidth/double(1<<currLevel);
currCell->setContainer(user_cell_initializer(currLevel,currMorton,arrayCoord,position));
});
}
}
void free_cell(Callback_free_cell user_cell_deallocator){
octree->forEachCell([&](CoreCell * currCell){
user_cell_deallocator(currCell->getContainer());
});
}
void execute_fmm(){
//FAssertLF(kernel,"No kernel set, please use scalfmm_user_kernel_config before calling the execute routine ... Exiting \n");
switch(Algorithm){
case 0:
{
typedef FFmmAlgorithm<OctreeClass,CoreCell,ContainerClass,CoreKernelClass,LeafInterface> AlgoClassSeq;
AlgoClassSeq algoSeq(octree,kernel);
algoSeq.execute();
break;
}
case 1:
{
typedef FFmmAlgorithmThread<OctreeClass,CoreCell,ContainerClass,CoreKernelClass,LeafClass> AlgoClassThread;
AlgoClassThread algoThread(octree,kernel);
algoThread.execute();
break;
}
case 2:
{
typedef FFmmAlgorithmPeriodic<OctreeClass,CoreCell,ContainerClass,CoreKernelClass,LeafClass> AlgoClassPeriodic;
AlgoClassPeriodic algoPeriod(octree,2);
algoPeriod.setKernel(kernel);
algoPeriod.execute();
break;
}
default :
std::cout<< "No algorithm found (probably for strange reasons) : "<< Algorithm <<" exiting" << std::endl;
}
}
};
#endif //FUSERKERNELENGINE_HPP
// ===================================================================================
// Logiciel initial: ScalFmm Version 0.5
// Co-auteurs : Olivier Coulaud, Bérenger Bramas.
// Propriétaires : INRIA.
// Copyright © 2011-2012, diffusé sous les termes et conditions d’une licence propriétaire.
// Initial software: ScalFmm Version 0.5
// Co-authors: Olivier Coulaud, Bérenger Bramas.
// Owners: INRIA.
// Copyright © 2011-2012, spread under the terms and conditions of a proprietary license.
// ===================================================================================
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
/**
* In this file we implement a example kernel which is simply printing information
* about the cells and particles.
* It is recommanded to compile and execute this code in order to understand the behavior
* of the application.
* We mark some part with "JUST-PUT-HERE:" to give instruction to user to create its own kernel.
**/
// Include the FMM API (should be in a different folder for you)
#include "../Src/CScalfmmApi.h"
// Uncomment the next line to avoid the print in the kernel and verbosity
// #define NOT_TOO_MUCH_VERBOSE
#ifdef NOT_TOO_MUCH_VERBOSE
#define VerbosePrint(X...)
#else
#define VerbosePrint(X...) printf(X)
#endif
///////////////////////////////////////////////////////////////////////////
/// Cell struct and functions
///////////////////////////////////////////////////////////////////////////
// This represent a cell
struct MyCellDescriptor{
int level;
long long mortonIndex;
int coord[3];
double position[3];
// JUST-PUT-HERE:
// You local and multipole arrays
};
// This is our function that init a cell (struct MyCellDescriptor)
void* my_Callback_init_cell(int level, long long mortonIndex, int* coord, double* position){
VerbosePrint("\tAllocating cell for level %d, morton index %lld, coord %d/%d/%d\n", level, mortonIndex, coord[0], coord[1], coord[2]);
struct MyCellDescriptor* cellData = (struct MyCellDescriptor*)malloc(sizeof(struct MyCellDescriptor));
memset(cellData, 0, sizeof(struct MyCellDescriptor));
cellData->level = level;
cellData->mortonIndex = mortonIndex;
memcpy(cellData->coord, coord, sizeof(int)*3);
memcpy(cellData->position, position, sizeof(double)*3);
// JUST-PUT-HERE:
// Fill your structure
return cellData;
}
// To dealloc a cell
void my_Callback_free_cell(void* cellData){
free(cellData);
}
///////////////////////////////////////////////////////////////////////////
/// Kernel struct and functions
///////////////////////////////////////////////////////////////////////////
// This is the data passed to our kernel
struct MyData {
// We simply strore the number of call the and particle indexes
double* insertedPositions;
int countP2M;
int countM2M;
int countM2L;
int countL2L;
int countL2P;
int countP2PInner;
int countP2P;
// JUST-PUT-HERE:
// everything your kernel will need for its computation
// pre-computed matrices, etc....
};
// Our P2M
void my_Callback_P2M(void* cellData, int nbParticlesInLeaf, const int* particleIndexes, void* userData){
struct MyData* my_data = (struct MyData*)userData;
my_data->countP2M += 1;
struct MyCellDescriptor* my_cell = (struct MyCellDescriptor*) cellData;
VerbosePrint("Cell morton %lld is doing P2M with %d particles :\n", my_cell->mortonIndex, nbParticlesInLeaf);
int idxPart;
for(idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
double* position = &my_data->insertedPositions[particleIndexes[idxPart]*3];
VerbosePrint("\t particle idx %d position %e/%e%e\n", particleIndexes[idxPart],
position[0], position[1], position[2]);
// JUST-PUT-HERE:
// Your real P2M computation
}
}
void my_Callback_M2M(int level, void* cellData, int childPosition, void* childData, void* userData){
struct MyData* my_data = (struct MyData*)userData;
my_data->countM2M += 1;
struct MyCellDescriptor* my_cell = (struct MyCellDescriptor*) cellData;
struct MyCellDescriptor* my_child = (struct MyCellDescriptor*) childData;
int childFullPosition[3];
Scalfmm_utils_parentChildPosition(childPosition, childFullPosition);
VerbosePrint("Doing a M2M at level %d, between cells %lld and %lld (position %d/%d/%d)\n",
level, my_cell->mortonIndex, my_child->mortonIndex,
childFullPosition[0], childFullPosition[1], childFullPosition[2]);
// JUST-PUT-HERE: your real M2M computation
}
void my_Callback_M2L(int level, void* cellData, int srcPosition, void* srcData, void* userData){
struct MyData* my_data = (struct MyData*)userData;
my_data->countM2L += 1;
struct MyCellDescriptor* my_cell = (struct MyCellDescriptor*) cellData;
struct MyCellDescriptor* my_src_cell = (struct MyCellDescriptor*) srcData;
int interactionFullPosition[3];
Scalfmm_utils_interactionPosition(srcPosition, interactionFullPosition);
VerbosePrint("Doing a M2L at level %d, between cells %lld and %lld (position %d/%d/%d)\n",
level, my_cell->mortonIndex, my_src_cell->mortonIndex,
interactionFullPosition[0], interactionFullPosition[1], interactionFullPosition[2]);
// JUST-PUT-HERE: Your M2L
}
void my_Callback_L2L(int level, void* cellData, int childPosition, void* childData, void* userData){
struct MyData* my_data = (struct MyData*)userData;
my_data->countL2L += 1;
struct MyCellDescriptor* my_cell = (struct MyCellDescriptor*) cellData;
struct MyCellDescriptor* my_child = (struct MyCellDescriptor*) childData;
int childFullPosition[3];
Scalfmm_utils_parentChildPosition(childPosition, childFullPosition);
VerbosePrint("Doing a L2L at level %d, between cells %lld and %lld (position %d/%d/%d)\n",
level, my_cell->mortonIndex, my_child->mortonIndex,
childFullPosition[0], childFullPosition[1], childFullPosition[2]);
// JUST-PUT-HERE: Your L2L
}
void my_Callback_L2P(void* cellData, int nbParticlesInLeaf, const int* particleIndexes, void* userData){
struct MyData* my_data = (struct MyData*)userData;
my_data->countL2P += 1;
struct MyCellDescriptor* my_cell = (struct MyCellDescriptor*) cellData;
VerbosePrint("Cell morton %lld is doing L2P with %d particles :\n", my_cell->mortonIndex, nbParticlesInLeaf);
int idxPart;
for(idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
double* position = &my_data->insertedPositions[particleIndexes[idxPart]*3];
VerbosePrint("\t particle idx %d position %e/%e%e\n", particleIndexes[idxPart],
position[0], position[1], position[2]);
// JUST-PUT-HERE: Your L2P
}
}
void my_Callback_P2P(int nbParticlesInLeaf, const int* particleIndexes, int nbParticlesInSrc, const int* particleIndexesSrc, void* userData){
struct MyData* my_data = (struct MyData*)userData;
my_data->countP2P += 1;
VerbosePrint("Doing P2P between two leaves of %d particles and %d particles :\n", nbParticlesInLeaf, nbParticlesInSrc);
int idxPart;
for(idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
double* position = &my_data->insertedPositions[particleIndexes[idxPart]*3];
VerbosePrint("\t Target >> particle idx %d position %e/%e%e\n", particleIndexes[idxPart],
position[0], position[1], position[2]);
}
for(idxPart = 0 ; idxPart < nbParticlesInSrc ; ++idxPart){
double* position = &my_data->insertedPositions[particleIndexesSrc[idxPart]*3];
VerbosePrint("\t Target >> Src idx %d position %e/%e%e\n", particleIndexesSrc[idxPart],
position[0], position[1], position[2]);
}
// JUST-PUT-HERE:
// Put one loop into the other to make all particles from source
// interacting with the target particles
}
void my_Callback_P2PInner(int nbParticlesInLeaf, const int* particleIndexes, void* userData){
struct MyData* my_data = (struct MyData*)userData;
my_data->countP2PInner += 1;
VerbosePrint("Doing P2P inside a leaf of %d particles :\n", nbParticlesInLeaf);
int idxPart;
for(idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
double* position = &my_data->insertedPositions[particleIndexes[idxPart]*3];
VerbosePrint("\t particle idx %d position %e/%e%e\n", particleIndexes[idxPart],
position[0], position[1], position[2]);
// JUST-PUT-HERE: another loop to have all particles interacting with
// the other from the same leaf
}
}
///////////////////////////////////////////////////////////////////////////
/// Main
///////////////////////////////////////////////////////////////////////////
// Simply create particles and try the kernels
int main(int argc, char ** argv){
// The properties of our tree
int treeHeight = 4;
double boxWidth = 1.0;
double boxCenter[3];
boxCenter[0] = boxCenter[1] = boxCenter[2] = 0.0;
// Create random particles
int nbParticles = 10;
int particleIndexes[nbParticles];
double particleXYZ[nbParticles*3];
{
printf("Creating Particles:\n");
int idxPart;
for(idxPart = 0 ; idxPart < nbParticles ;