Skip to content

GitLab

S
ScalFMM
Commit 38ec2433 authored by BRAMAS Berenger's avatar BRAMAS Berenger
Browse files
Options

Merge branch 'master' of git+ssh://scm.gforge.inria.fr//gitroot//scalfmm/scalfmm

parents cca87d2d 0b684086
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 2546 additions and 911 deletions
Addons/CKernelApi/Src/CScalfmmApi.h
View file @ 38ec2433
......@@ -90,7 +90,8 @@ typedef enum kernel_type {
typedef enum scalfmm_algorithm_config {
sequential = 0, /* Use the sequential version of Scalfmm*/
multi_thread = 1, /* Use the Multi thread version of Scalfmm*/
periodic = 2 /* Use the periodic version of Scalfmm*/
periodic = 2, /* Use the periodic version of Scalfmm*/
source_target = 3 /* USe the source/target algorithm */
} scalfmm_algorithm;
......@@ -110,7 +111,7 @@ typedef void* scalfmm_handle;
* Every data will be stored in order to crete later through a builder
* what is needed for the simulation
*/
scalfmm_handle scalfmm_init( scalfmm_kernel_type KernelType);
scalfmm_handle scalfmm_init( scalfmm_kernel_type KernelType,scalfmm_algorithm algo);
/////////////////////////////////////////////////////////////////////
......@@ -128,8 +129,9 @@ scalfmm_handle scalfmm_init( scalfmm_kernel_type KernelType);
* @param tree_position int[3] position inside the tree (number of boxes in
* each direction)
* @param spatial_position double[3] center of the cell
* @param inDatas user generic pointer to kernel.
*/
typedef void* (*Callback_init_cell)(int level, long long morton_index, int* tree_position, double* spatial_position);
typedef void* (*Callback_init_cell)(int level, long long morton_index, int* tree_position, double* spatial_position, void * inDatas);
/**
* Function to destroy what have bee initialized by the user (should
......@@ -160,33 +162,67 @@ typedef struct User_Scalfmm_Cell_Descriptor{
void scalfmm_build_tree(scalfmm_handle handle,int TreeHeight,double BoxWidth,double* BoxCenter,Scalfmm_Cell_Descriptor user_cell_descriptor);
/**
* @brief This enum flag is to know if function calling will deal with
* source, target or both
*/
typedef enum particule_type{
SOURCE=0,
TARGET=1,
BOTH=2
}PartType;
/**
* @brief This function insert an array of position into the octree
* @brief This function insert alongside to position an arbitrary
* number of attributes.
* @param Handle scalfmm_handle provided by scalfmm_init.
* @param NbPartToInsert number of particles to be inserted
* @param nbAttributeToInsert number of attribute to insert (this
* number will be > 3 (because we need at least 3 doubles for
* position))
* @param strideForEachAtt How to get each attribute for each particles
* @param rawDatas datas to be read
*
* Example :
struct part{
double[3] position;
double charge;
double test; //not used
double coeff;
};
Then nbAttributeToInsert will be 3+1+1
and strideForEachAtt will be : [0,1,2,3,5]
*/
void scalfmm_tree_abstract_insert(scalfmm_handle Handle, int NbPartToInsert, int nbAttributeToInsert, int * strideForEachAtt,
double* rawDatas);
/**
* @brief This function insert an array of position into the
* octree. THis fonction will insert particules with no SOURCE/TARGET
* type.
* @param Handle scalfmm_handle provided by scalfmm_init.
* @param NbPositions Number of position to be inserted
* @param arrayX Array containing the X coordinate for all the parts, size : NbPositions
* @param arrayY Array containing the Y coordinate for all the parts, size : NbPositions
* @param arrayZ Array containing the Z coordinate for all the parts, size : NbPositions
*
* @param type : type to insert
* The parts will be inserted with their indices inside the
* array. Each index will be unique.
*
* In case several calls are performed to scalfmm_tree_insert_arrays,
* first call with N particles, and then with M particles.
* The second call particles will have index from [N ; N+M].
*
* Default physical values, potential and forces are set to 0.
*/
void scalfmm_tree_insert_particles(scalfmm_handle Handle, int NbPositions, double * arrayX, double * arrayY, double * arrayZ);
void scalfmm_tree_insert_particles(scalfmm_handle Handle, int NbPositions, double * arrayX, double * arrayY, double * arrayZ, PartType type);
/**
* This function is equivalent to scalfmm_tree_insert_particles
* but the given array XYZ should contains a triple value per paticles.
*/
void scalfmm_tree_insert_particles_xyz(scalfmm_handle Handle, int NbPositions, double * XYZ);
void scalfmm_tree_insert_particles_xyz(scalfmm_handle Handle, int NbPositions, double * XYZ, PartType type);
/**
* @brief This function set the physical values of all the particles
......@@ -196,17 +232,17 @@ void scalfmm_tree_insert_particles_xyz(scalfmm_handle Handle, int NbPositions, d
* inserted.
* @param physicalValues Array containing the physical values to be
* associated to each parts.
*
* @param type : type of the particules to be setted.
* The physical values will be stored according to their indices in
* the array. First particle inserted will take value physicalValues[0].
*/
void scalfmm_set_physical_values(scalfmm_handle Handle, int nbPhysicalValues, double * physicalValues);
void scalfmm_set_physical_values(scalfmm_handle Handle, int nbPhysicalValues, double * physicalValues, PartType type);
/**
* @brief get the physical values.
*
* WARNING : the user must allocate (and initialize) the array given
*/
void scalfmm_get_physical_values(scalfmm_handle Handle, int nbPhysicalValues, double * physicalValues);
void scalfmm_get_physical_values(scalfmm_handle Handle, int nbPhysicalValues, double * physicalValues, PartType type);
/**
*
......@@ -215,6 +251,7 @@ void scalfmm_get_physical_values(scalfmm_handle Handle, int nbPhysicalValues, do
* @param idxOfParticles an array of indexes of size nbPhysicalValues to know which particles
* to set the values to.
* @param physicalValues the physical values.
* @param type : type of the particules to be setted.
*
* For example to set the physical values to particles 0 and 1 to values 1.1 and 1.4:
* @code nbPhysicalValues = 2;
......@@ -224,9 +261,9 @@ void scalfmm_get_physical_values(scalfmm_handle Handle, int nbPhysicalValues, do
* Be aware that such approach requiere to find particles in the tree which can have high cost.
*/
void scalfmm_set_physical_values_npart(scalfmm_handle Handle, int nbPhysicalValues,
int* idxOfParticles, double * physicalValues);
int* idxOfParticles, double * physicalValues, PartType type);
void scalfmm_get_physical_values_npart(scalfmm_handle Handle, int nbPhysicalValues,
int* idxOfParticles, double * physicalValues);
int* idxOfParticles, double * physicalValues, PartType type);
/**
......@@ -236,14 +273,17 @@ void scalfmm_get_physical_values_npart(scalfmm_handle Handle, int nbPhysicalValu
* to the number of parts inserted)
* @param forcesToFill array of size nbParts*3, that will contains the
* forces. WARNING : User must allocate the array before call.
*
* @param type : type of the particules to be setted.
* Forces will be stored sequentially, according to the indices in the
* array. (ie fx1,fy1,fz1,fx2,fy2,fz2,fx3 ....)
* @param idxOfParticles : array of indices of the particles
* wanted. If used, then the parts will be given in the order of
* idxOfParticles
*/
void scalfmm_get_forces_xyz(scalfmm_handle Handle, int nbParts, double * forcesToFill);
void scalfmm_get_forces_xyz_npart(scalfmm_handle Handle, int nbParts, int* idxOfParticles, double * forcesToFill);
void scalfmm_get_forces(scalfmm_handle Handle, int nbParts, double * fX, double* fY, double* fZ);
void scalfmm_get_forces_npart(scalfmm_handle Handle, int nbParts, int* idxOfParticles, double * fX, double* fY, double* fZ);
void scalfmm_get_forces_xyz(scalfmm_handle Handle, int nbParts, double * forcesToFill, PartType type);
void scalfmm_get_forces_xyz_npart(scalfmm_handle Handle, int nbParts, int* idxOfParticles, double * forcesToFill, PartType type);
void scalfmm_get_forces(scalfmm_handle Handle, int nbParts, double * fX, double* fY, double* fZ, PartType type);
void scalfmm_get_forces_npart(scalfmm_handle Handle, int nbParts, int* idxOfParticles, double * fX, double* fY, double* fZ, PartType type);
/**
......@@ -257,12 +297,12 @@ void scalfmm_get_forces_npart(scalfmm_handle Handle, int nbParts, int* idxOfPart
* forces . WARNING : User must allocate the array before call.
* @param forcesZ array of size nbParts, that will contains the
* forces . WARNING : User must allocate the array before call.
*
* @param type : type of the particules to be setted.
*/
void scalfmm_set_forces_xyz(scalfmm_handle Handle, int nbParts, double * forcesToFill);
void scalfmm_set_forces_xyz_npart(scalfmm_handle Handle, int nbParts, int* idxOfParticles, double * forcesToFill);
void scalfmm_set_forces(scalfmm_handle Handle, int nbParts, double * fX, double* fY, double* fZ);
void scalfmm_set_forces_npart(scalfmm_handle Handle, int nbParts, int* idxOfParticles, double * fX, double* fY, double* fZ);
void scalfmm_set_forces_xyz(scalfmm_handle Handle, int nbParts, double * forcesToFill, PartType type);
void scalfmm_set_forces_xyz_npart(scalfmm_handle Handle, int nbParts, int* idxOfParticles, double * forcesToFill, PartType type);
void scalfmm_set_forces(scalfmm_handle Handle, int nbParts, double * fX, double* fY, double* fZ, PartType type);
void scalfmm_set_forces_npart(scalfmm_handle Handle, int nbParts, int* idxOfParticles, double * fX, double* fY, double* fZ, PartType type);
/**
......@@ -272,14 +312,15 @@ void scalfmm_set_forces_npart(scalfmm_handle Handle, int nbParts, int* idxOfPart
* to the number of parts inserted)
* @param potentialsToFill array of potentials to be filled. WARNING :
* User must allocate the array before call.
* @param type : type of the particules to be setted.
*
* Potentials will be stored sequentially, according to the indices in the
* array.
*/
void scalfmm_get_potentials(scalfmm_handle Handle, int nbParts, double * potentialsToFill);
void scalfmm_set_potentials(scalfmm_handle Handle, int nbParts, double * potentialsToRead);
void scalfmm_get_potentials_npart(scalfmm_handle Handle, int nbParts, int* idxOfParticles, double * potentialsToFill);
void scalfmm_set_potentials_npart(scalfmm_handle Handle, int nbParts, int* idxOfParticles, double * potentialsToRead);
void scalfmm_get_potentials(scalfmm_handle Handle, int nbParts, double * potentialsToFill, PartType type);
void scalfmm_set_potentials(scalfmm_handle Handle, int nbParts, double * potentialsToRead, PartType type);
void scalfmm_get_potentials_npart(scalfmm_handle Handle, int nbParts, int* idxOfParticles, double * potentialsToFill, PartType type);
void scalfmm_set_potentials_npart(scalfmm_handle Handle, int nbParts, int* idxOfParticles, double * potentialsToRead, PartType type);
/**
......@@ -290,11 +331,12 @@ void scalfmm_set_potentials_npart(scalfmm_handle Handle, int nbParts, int* idxOf
* number of parts inserted)
* @param updatedXYZ array of displacement (ie
* dx1,dy1,dz1,dx2,dy2,dz2,dx3 ...)
* @param type : type of the particules to be setted.
*/
void scalfmm_add_to_positions_xyz(scalfmm_handle Handle, int NbPositions, double * updatedXYZ);
void scalfmm_add_to_positions(scalfmm_handle Handle, int NbPositions, double * X, double * Y , double * Z);
void scalfmm_add_to_positions_xyz_npart(scalfmm_handle Handle, int NbPositions, int* idxOfParticles, double * updatedXYZ);
void scalfmm_add_to_positions_npart(scalfmm_handle Handle, int NbPositions, int* idxOfParticles, double * X, double * Y , double * Z);
void scalfmm_add_to_positions_xyz(scalfmm_handle Handle, int NbPositions, double * updatedXYZ, PartType type);
void scalfmm_add_to_positions(scalfmm_handle Handle, int NbPositions, double * X, double * Y , double * Z, PartType type);
void scalfmm_add_to_positions_xyz_npart(scalfmm_handle Handle, int NbPositions, int* idxOfParticles, double * updatedXYZ, PartType type);
void scalfmm_add_to_positions_npart(scalfmm_handle Handle, int NbPositions, int* idxOfParticles, double * X, double * Y , double * Z, PartType type);
/**
......@@ -305,17 +347,17 @@ void scalfmm_add_to_positions_npart(scalfmm_handle Handle, int NbPositions, int*
* number of parts inserted)
* @param newXYZ array of new positions (ie
* dx1,dy1,dz1,dx2,dy2,dz2,dx3 ...)
*
* @param type : type of the particules to be setted.
* @return Error code, a parts may move out of the simulation
* box. ScalFMM cannot deals with that specific case. Error code :
* 0. Success code 1. Could be an arg in order to be Fortran
* compliant.
*
*/
void scalfmm_set_positions_xyz(scalfmm_handle Handle, int NbPositions, double * updatedXYZ);
void scalfmm_set_positions(scalfmm_handle Handle, int NbPositions, double * X, double * Y , double * Z);
void scalfmm_set_positions_xyz_npart(scalfmm_handle Handle, int NbPositions, int* idxOfParticles, double * updatedXYZ);
void scalfmm_set_positions_npart(scalfmm_handle Handle, int NbPositions, int* idxOfParticles, double * X, double * Y , double * Z);
void scalfmm_set_positions_xyz(scalfmm_handle Handle, int NbPositions, double * updatedXYZ, PartType type);
void scalfmm_set_positions(scalfmm_handle Handle, int NbPositions, double * X, double * Y , double * Z, PartType type);
void scalfmm_set_positions_xyz_npart(scalfmm_handle Handle, int NbPositions, int* idxOfParticles, double * updatedXYZ, PartType type);
void scalfmm_set_positions_npart(scalfmm_handle Handle, int NbPositions, int* idxOfParticles, double * X, double * Y , double * Z, PartType type);
/**
*@brief This is function is to be called after a call modifying some
......@@ -324,10 +366,10 @@ void scalfmm_set_positions_npart(scalfmm_handle Handle, int NbPositions, int* id
*/
void scalfmm_update_tree(scalfmm_handle handle);
void scalfmm_get_positions_xyz(scalfmm_handle Handle, int NbPositions, double * positionsToFill);
void scalfmm_get_positions(scalfmm_handle Handle, int NbPositions, double * X, double * Y , double * Z);
void scalfmm_get_positions_xyz_npart(scalfmm_handle Handle, int NbPositions, int* idxOfParticles, double * positionsToFill);
void scalfmm_get_positions_npart(scalfmm_handle Handle, int NbPositions, int* idxOfParticles, double * X, double * Y , double * Z);
void scalfmm_get_positions_xyz(scalfmm_handle Handle, int NbPositions, double * positionsToFill, PartType type);
void scalfmm_get_positions(scalfmm_handle Handle, int NbPositions, double * X, double * Y , double * Z, PartType type);
void scalfmm_get_positions_xyz_npart(scalfmm_handle Handle, int NbPositions, int* idxOfParticles, double * positionsToFill, PartType type);
void scalfmm_get_positions_npart(scalfmm_handle Handle, int NbPositions, int* idxOfParticles, double * X, double * Y , double * Z, PartType type);
/* /\** */
/* * @brief This function provides a way for the user to define scalfmm */
......@@ -446,7 +488,6 @@ typedef void (*Callback_P2PFull)(FSize nbParticles, const FSize* particleIndexes
*/
typedef void (*Callback_P2PInner)(FSize nbParticles, const FSize* particleIndexes, void* userData);
/**
* @brief Function to be filled by user's method to reset a user's cell
* @param level level of the cell.
......@@ -458,7 +499,6 @@ typedef void (*Callback_P2PInner)(FSize nbParticles, const FSize* particleIndexe
*/
typedef void (*Callback_reset_cell)(int level, long long morton_index, int* tree_position, double* spatial_position, void * userCell);
/**
* @brief Structure containing callbacks to fill in order to define
* user kernel.
......@@ -535,13 +575,30 @@ void scalfmm_dealloc_handle(scalfmm_handle handle, Callback_free_cell cellDestro
*/
void scalfmm_reset_tree(scalfmm_handle handle, Callback_reset_cell cellReseter);
/////////////////////////////////////////////////////////////////////
/////////////// Monitoring functions /////////////////
/////////////////////////////////////////////////////////////////////
/**
* @brief This function shouldn't be there !! display information
* about the octree built versus the octree hibox want.
* @brief Scalfmm has a built in feature to get the time elapsed in
* each operator. This function returns the number of different timers.
* @param Handle scalfmm_handle provided by scalfmm_init.
* @param Rinflu influence radius for each particle previously
* inserted. Tree must be built before calling this function
* @return Number of timers
*/
void scalfmm_hibox_Rinflu_display(scalfmm_handle Handle, FSize nbPart, double * Rinflu);
int scalfmm_get_nb_timers(scalfmm_handle handle);
/**
* @brief Scalfmm has a built in feature to get the time elapsed in
* each operator. This function fill the array with elapsed time for
* each operator.
* @param Handle scalfmm_handle provided by scalfmm_init.
* @param Array of Timers, to be allocated by the user (using
* scalfmm_get_nb_timers)
* Order inside the array : P2M, M2M, M2L, L2L, L2P, P2P, NearField
* (P2P+L2P).
*/
void scalfmm_get_timers(scalfmm_handle handle,double * Timers);
#endif
Addons/CKernelApi/Src/FInterEngine.hpp
View file @ 38ec2433
This diff is collapsed.
Addons/CKernelApi/Src/FScalFMMEngine.hpp
View file @ 38ec2433
This diff is collapsed.
Addons/CKernelApi/Src/FScalfmmApiInit.cpp
View file @ 38ec2433
......@@ -6,48 +6,94 @@ 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();
extern "C" scalfmm_handle scalfmm_init(/*int TreeHeight,double BoxWidth,double* BoxCenter, */scalfmm_kernel_type KernelType,
scalfmm_algorithm algo){
ScalFmmCoreHandle<double> * handle = new ScalFmmCoreHandle<double>();
typedef double FReal;
switch(KernelType){
case 0:
handle->engine = new FUserKernelEngine<FReal>(/*TreeHeight, BoxWidth, BoxCenter, */KernelType);
break;
if(algo == source_target){
case 1:
//TODO typedefs
typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
typedef FChebCell<FReal,7> ChebCell;
switch(KernelType){
case 0:
typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
typedef FTypedLeaf<FReal,ContainerClass> LeafClass;
typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
typedef FChebSymKernel<FReal,ChebCell,ContainerClass,MatrixKernelClass,7> ChebKernel;
handle->engine = new FUserKernelEngine<FReal,LeafClass>(/*TreeHeight, BoxWidth, BoxCenter, */KernelType);
break;
handle->engine = new FInterEngine<FReal,ChebCell,ChebKernel>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
break;
// case 2:
// //TODO typedefs
// typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
// typedef FUnifCell<7> UnifCell;
case 1:
//TODO typedefs
typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
typedef FTypedChebCell<FReal,7> ChebCell;
typedef FTypedLeaf<FReal,ContainerClass> LeafClass;
// typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
// typedef FUnifKernel<UnifCell,ContainerClass,MatrixKernelClass,7> UnifKernel;
typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
typedef FChebSymKernel<FReal,ChebCell,ContainerClass,MatrixKernelClass,7> ChebKernel;
// handle->engine = new FInterEngine<UnifCell,UnifKernel>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
// break;
handle->engine = new FInterEngine<FReal,ChebCell,ChebKernel,LeafClass>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
break;
// case 2:
// //TODO typedefs
// typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
// typedef FUnifCell<7> UnifCell;
default:
std::cout<< "Kernel type unsupported" << std::endl;
exit(0);
break;
// typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
// typedef FUnifKernel<UnifCell,ContainerClass,MatrixKernelClass,7> UnifKernel;
// handle->engine = new FInterEngine<UnifCell,UnifKernel>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
// break;
default:
std::cout<< "Kernel type unsupported" << std::endl;
exit(0);
break;
}
}
else{ //No Source/Targets distinction
switch(KernelType){
case 0:
typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
typedef FSimpleLeaf<FReal,ContainerClass> LeafClass;
handle->engine = new FUserKernelEngine<FReal,LeafClass>(/*TreeHeight, BoxWidth, BoxCenter, */KernelType);
break;
case 1:
//TODO typedefs
typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
typedef FChebCell<FReal,7> ChebCell;
typedef FSimpleLeaf<FReal,ContainerClass> LeafClass;
typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
typedef FChebSymKernel<FReal,ChebCell,ContainerClass,MatrixKernelClass,7> ChebKernel;
handle->engine = new FInterEngine<FReal,ChebCell,ChebKernel,LeafClass>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
break;
// case 2:
// //TODO typedefs
// typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
// typedef FUnifCell<7> UnifCell;
// typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
// typedef FUnifKernel<UnifCell,ContainerClass,MatrixKernelClass,7> UnifKernel;
// handle->engine = new FInterEngine<UnifCell,UnifKernel>(/*TreeHeight,BoxWidth,BoxCenter, */KernelType);
// break;
default:
std::cout<< "Kernel type unsupported" << std::endl;
exit(0);
break;
}
}
return handle;
return handle;
}
extern "C" void scalfmm_dealloc_handle(scalfmm_handle handle, Callback_free_cell userDeallocator){
((ScalFmmCoreHandle *) handle)->engine->intern_dealloc_handle(userDeallocator);
delete ((ScalFmmCoreHandle *) handle)->engine ;
delete (ScalFmmCoreHandle *) handle;
((ScalFmmCoreHandle<double> *) handle)->engine->intern_dealloc_handle(userDeallocator);
delete ((ScalFmmCoreHandle<double> *) handle)->engine ;
delete (ScalFmmCoreHandle<double> *) handle;
}
/**
......
Addons/CKernelApi/Src/FUserKernelEngine.hpp
View file @ 38ec2433
This diff is collapsed.
Addons/CKernelApi/Tests/testChebInterface.c
View file @ 38ec2433
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
//For timing monitoring
#include <time.h>
......@@ -20,7 +20,7 @@
/**
* @brief Wrapper to init internal ChebCell
*/
void* cheb_init_cell(int level, long long morton_index, int* tree_position, double* spatial_position){
void* cheb_init_cell(int level, long long morton_index, int* tree_position, double* spatial_position, void * KernelDatas){
return ChebCellStruct_create(morton_index,tree_position);
}
......@@ -55,6 +55,10 @@ void cheb_p2pFull(FSize nbParticles, const FSize* particleIndexes,
ChebKernel_P2P(nbParticles, particleIndexes, sourceParticleIndexes, sourceNbPart, userData);
}
void cheb_resetCell(int level, long long morton_index, int* tree_position, double* spatial_position, void * userCell){
ChebCell_reset(level,morton_index,tree_position,spatial_position,userCell);
}
/**
* @brief Wrapper on timeval struct
......@@ -158,10 +162,11 @@ int main(int argc, char ** av){
/* fclose(fd); */
}
scalfmm_handle handle = scalfmm_init(user_defined_kernel);
scalfmm_handle handle = scalfmm_init(user_defined_kernel,multi_thread);
//For Reference
scalfmm_handle handle_ref = scalfmm_init(chebyshev);
scalfmm_handle handle_ref = scalfmm_init(chebyshev,multi_thread);
//Struct for user defined kernel
struct User_Scalfmm_Cell_Descriptor cellDescriptor;
......@@ -180,13 +185,8 @@ int main(int argc, char ** av){
scalfmm_build_tree(handle,treeHeight, boxWidth, boxCenter, cellDescriptor);
scalfmm_build_tree(handle_ref,treeHeight, boxWidth, boxCenter, user_descr);
// Insert particles
printf("Inserting particles...\n");
scalfmm_tree_insert_particles_xyz(handle, nbPart, particleXYZ);
scalfmm_tree_insert_particles_xyz(handle_ref, nbPart, particleXYZ);
//Once is the tree built, one must set the kernel before inserting particles
//Set physical values for Cheb_ref
scalfmm_set_physical_values(handle_ref,nbPart,physicalValues);
//Set our callbacks
struct User_Scalfmm_Kernel_Descriptor kernel;
kernel.p2m = cheb_p2m;
......@@ -239,83 +239,125 @@ int main(int argc, char ** av){
//Give ScalFMM the datas before calling fmm (this will set as well the kernel)
scalfmm_user_kernel_config(handle,kernel,&userDatas);
//Set timers
Timer interface_timer,ref_timer;
//Execute FMM
tic(&interface_timer);
scalfmm_execute_fmm(handle/*, kernel, &my_data*/);
tac(&interface_timer);
// Insert particles
printf("Inserting particles...\n");
scalfmm_tree_insert_particles_xyz(handle, nbPart, particleXYZ,BOTH);
scalfmm_tree_insert_particles_xyz(handle_ref, nbPart, particleXYZ,BOTH);
//Reduction on forces & potential arrays
{
FSize idxPart;
for(idThreads=1 ; idThreads<nb_threads ; ++idThreads){
for(idxPart=0 ; idxPart<nbPart ; ++idxPart){
//Everything is stored in first array
forcesToStore[0][3*idxPart+0] += forcesToStore[idThreads][3*idxPart+0];
forcesToStore[0][3*idxPart+1] += forcesToStore[idThreads][3*idxPart+1];
forcesToStore[0][3*idxPart+2] += forcesToStore[idThreads][3*idxPart+2];
potentialToStore[0][idxPart] += potentialToStore[idThreads][idxPart];
}
}
}
printf("User defined Chebyshev done\n");
print_elapsed(&interface_timer);
//Set physical values for Cheb_ref
scalfmm_set_physical_values(handle_ref,nbPart,physicalValues,BOTH);
tic(&ref_timer);
scalfmm_execute_fmm(handle_ref/*, kernel, &my_data*/);
tac(&ref_timer);
printf("Intern Chebyshev done\n");
print_elapsed(&ref_timer);
//Set timers
Timer interface_timer,ref_timer;
int ite=0, max_ite=5;
while(ite<max_ite){
//Execute FMM
tic(&interface_timer);
scalfmm_execute_fmm(handle/*, kernel, &my_data*/);
tac(&interface_timer);
//Reduction on forces & potential arrays
{
FSize idxPart;
for(idThreads=1 ; idThreads<nb_threads ; ++idThreads){
for(idxPart=0 ; idxPart<nbPart ; ++idxPart){
//Everything is stored in first array
forcesToStore[0][3*idxPart+0] += forcesToStore[idThreads][3*idxPart+0];
forcesToStore[0][3*idxPart+1] += forcesToStore[idThreads][3*idxPart+1];
forcesToStore[0][3*idxPart+2] += forcesToStore[idThreads][3*idxPart+2];
potentialToStore[0][idxPart] += potentialToStore[idThreads][idxPart];
}
}
}
//Print time results
print_difference_elapsed(&interface_timer,&ref_timer);
printf("User defined Chebyshev done\n");
print_elapsed(&interface_timer);
tic(&ref_timer);
scalfmm_execute_fmm(handle_ref/*, kernel, &my_data*/);
tac(&ref_timer);
printf("Intern Chebyshev done\n");
print_elapsed(&ref_timer);
//Print time results
print_difference_elapsed(&interface_timer,&ref_timer);
//get back the forces & potentials for ref_cheb execution
double * forcesRef = malloc(sizeof(double)*3*nbPart);
double * potentialsRef = malloc(sizeof(double)*nbPart);
memset(forcesRef,0,sizeof(double)*3*nbPart);
memset(potentialsRef,0,sizeof(double)*nbPart);
scalfmm_get_forces_xyz(handle_ref,nbPart,forcesRef,BOTH);
scalfmm_get_potentials(handle_ref,nbPart,potentialsRef,BOTH);
//scalfmm_print_everything(handle_ref);
{//Comparison part
FSize idxPart;
int nbPartOkay = 0;
for(idxPart=0 ; idxPart<nbPart ; ++idxPart ){
double diffX,diffY,diffZ,diffPot;
diffX = fabs( forcesToStore[0][idxPart*3+0]-forcesRef[