Commit 54c9f4fb authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille
Browse files
parents c8ea6c26 f213bce7
......@@ -45,261 +45,272 @@
template <unsigned NbAttributesPerParticle, class AttributeClass = FReal >
class FBasicParticleContainer : public FAbstractParticleContainer, public FAbstractSerializable {
protected:
/** The number of particles in the container */
int nbParticles;
/** 3 pointers to 3 arrays of real to store the position */
FReal* positions[3];
/** The attributes requested by the user */
AttributeClass* attributes[NbAttributesPerParticle];
/** The allocated memory */
int allocatedParticles;
/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
/** Ending call for pushing the attributes */
template<int index>
void addParticleValue(const int /*insertPosition*/){
}
/** Filling call for each attributes values */
template<int index, typename... Args>
void addParticleValue(const int insertPosition, const AttributeClass value, Args... args){
// Compile test to ensure indexing
static_assert(index < NbAttributesPerParticle, "Index to get attributes is out of scope.");
// insert the value
attributes[index][insertPosition] = value;
// Continue for reamining values
addParticleValue<index+1>( insertPosition, args...);
}
/** The number of particles in the container */
int nbParticles;
/** 3 pointers to 3 arrays of real to store the position */
FReal* positions[3];
/** The attributes requested by the user */
AttributeClass* attributes[NbAttributesPerParticle];
/** The allocated memory */
int allocatedParticles;
/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
/** Ending call for pushing the attributes */
template<int index>
void addParticleValue(const int /*insertPosition*/){
}
/** Filling call for each attributes values */
template<int index, typename... Args>
void addParticleValue(const int insertPosition, const AttributeClass value, Args... args){
// Compile test to ensure indexing
static_assert(index < NbAttributesPerParticle, "Index to get attributes is out of scope.");
// insert the value
attributes[index][insertPosition] = value;
// Continue for reamining values
addParticleValue<index+1>( insertPosition, args...);
}
public:
/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
FBasicParticleContainer(const FBasicParticleContainer&) = delete;
FBasicParticleContainer& operator=(const FBasicParticleContainer&) = delete;
FBasicParticleContainer(const FBasicParticleContainer&) = delete;
FBasicParticleContainer& operator=(const FBasicParticleContainer&) = delete;
/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
/** Basic contructor */
FBasicParticleContainer() : nbParticles(0), allocatedParticles(0){
memset(positions, 0, sizeof(positions[0]) * 3);
memset(attributes, 0, sizeof(attributes[0]) * NbAttributesPerParticle);
}
/** Basic contructor */
FBasicParticleContainer() : nbParticles(0), allocatedParticles(0){
memset(positions, 0, sizeof(positions[0]) * 3);
memset(attributes, 0, sizeof(attributes[0]) * NbAttributesPerParticle);
}
/** Simply dalloc the memory using first pointer
/** Simply dalloc the memory using first pointer
*/
~FBasicParticleContainer(){
FAlignedMemory::Dealloc32BAligned(positions[0]);
}
/**
~FBasicParticleContainer(){
FAlignedMemory::Dealloc32BAligned(positions[0]);
}
/**
* @brief getNbParticles
* @return the number of particles
*/
int getNbParticles() const{
return nbParticles;
}
/**
int getNbParticles() const{
return nbParticles;
}
/**
* @brief reset the number of particles
* @warning Only the number of particles is set to 0, the particles are still here.
*/
void resetNumberOfParticles()
{
nbParticles = 0 ;
}
/**
void resetNumberOfParticles()
{
nbParticles = 0 ;
}
/**
* @brief getPositions
* @return a FReal*[3] to get access to the positions
*/
const FReal*const* getPositions() const {
return positions;
}
const FReal*const* getPositions() const {
return positions;
}
/**
/**
* @brief getWPositions
* @return get the position in write mode
*/
FReal* const* getWPositions() {
return positions;
}
FReal* const* getWPositions() {
return positions;
}
/**
/**
* @brief getAttribute
* @param index
* @return the attribute at index index
*/
AttributeClass* getAttribute(const int index) {
return attributes[index];
}
AttributeClass* getAttribute(const int index) {
return attributes[index];
}
/**
/**
* @brief getAttribute
* @param index
* @return
*/
const AttributeClass* getAttribute(const int index) const {
return attributes[index];
}
const AttributeClass* getAttribute(const int index) const {
return attributes[index];
}
/**
/**
* Get the attribute with a forcing compile optimization
*/
template <int index>
AttributeClass* getAttribute() {
static_assert(index < NbAttributesPerParticle, "Index to get attributes is out of scope.");
return attributes[index];
}
template <int index>
AttributeClass* getAttribute() {
static_assert(index < NbAttributesPerParticle, "Index to get attributes is out of scope.");
return attributes[index];
}
/**
/**
* Get the attribute with a forcing compile optimization
*/
template <int index>
const AttributeClass* getAttribute() const {
static_assert(index < NbAttributesPerParticle, "Index to get attributes is out of scope.");
return attributes[index];
}
template <int index>
const AttributeClass* getAttribute() const {
static_assert(index < NbAttributesPerParticle, "Index to get attributes is out of scope.");
return attributes[index];
}
/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
/**
/**
* Push called bu FSimpleLeaf
* Should have a particle position fallowed by attributes
*/
template<typename... Args>
void push(const FPoint& inParticlePosition, Args... args){
// enought space?
if( nbParticles == allocatedParticles ){
// allocate memory
const int moduloParticlesNumber = (32/sizeof(FReal)); // We want to be rounded to 32B
allocatedParticles = (FMath::Max(10,int(FReal(nbParticles+1)*1.5)) + moduloParticlesNumber - 1) & ~(moduloParticlesNumber-1);
// init with 0
const size_t allocatedBytes = (sizeof(FReal)*3 + sizeof(AttributeClass)*NbAttributesPerParticle)*allocatedParticles;
FReal* newData = reinterpret_cast<FReal*>(FAlignedMemory::Allocate32BAligned(allocatedBytes));
memset( newData, 0, allocatedBytes);
// copy memory
const char*const toDelete = reinterpret_cast<const char*>(positions[0]);
for(int idx = 0 ; idx < 3 ; ++idx){
memcpy(newData + (allocatedParticles * idx), positions[idx], sizeof(FReal) * nbParticles);
positions[idx] = newData + (allocatedParticles * idx);
}
// copy attributes
AttributeClass* startAddress = reinterpret_cast<AttributeClass*>(positions[2] + allocatedParticles);
for(unsigned idx = 0 ; idx < NbAttributesPerParticle ; ++idx){
memcpy(startAddress + (allocatedParticles * idx), attributes[idx], sizeof(AttributeClass) * nbParticles);
attributes[idx] = startAddress + (idx * allocatedParticles);
}
// delete old
FAlignedMemory::Dealloc32BAligned(toDelete);
template<typename... Args>
void push(const FPoint& inParticlePosition, Args... args){
// enought space?
if( nbParticles == allocatedParticles ){
// allocate memory
const int moduloParticlesNumber = (32/sizeof(FReal)); // We want to be rounded to 32B
allocatedParticles = (FMath::Max(10,int(FReal(nbParticles+1)*1.5)) + moduloParticlesNumber - 1) & ~(moduloParticlesNumber-1);
// init with 0
const size_t allocatedBytes = (sizeof(FReal)*3 + sizeof(AttributeClass)*NbAttributesPerParticle)*allocatedParticles;
FReal* newData = reinterpret_cast<FReal*>(FAlignedMemory::Allocate32BAligned(allocatedBytes));
memset( newData, 0, allocatedBytes);
// copy memory
const char*const toDelete = reinterpret_cast<const char*>(positions[0]);
for(int idx = 0 ; idx < 3 ; ++idx){
memcpy(newData + (allocatedParticles * idx), positions[idx], sizeof(FReal) * nbParticles);
positions[idx] = newData + (allocatedParticles * idx);
}
// copy attributes
AttributeClass* startAddress = reinterpret_cast<AttributeClass*>(positions[2] + allocatedParticles);
for(unsigned idx = 0 ; idx < NbAttributesPerParticle ; ++idx){
memcpy(startAddress + (allocatedParticles * idx), attributes[idx], sizeof(AttributeClass) * nbParticles);
attributes[idx] = startAddress + (idx * allocatedParticles);
}
// delete old
FAlignedMemory::Dealloc32BAligned(toDelete);
}
// insert particle data
positions[0][nbParticles] = inParticlePosition.getX();
positions[1][nbParticles] = inParticlePosition.getY();
positions[2][nbParticles] = inParticlePosition.getZ();
// insert attribute data
addParticleValue<0>( nbParticles, args...);
nbParticles += 1;
}
// insert particle data
positions[0][nbParticles] = inParticlePosition.getX();
positions[1][nbParticles] = inParticlePosition.getY();
positions[2][nbParticles] = inParticlePosition.getZ();
// insert attribute data
addParticleValue<0>( nbParticles, args...);
nbParticles += 1;
}
/**
/**
* Push called usually by FTypedLeaf with the isTarget flag in addition
*/
template<typename... Args>
void push(const FPoint& inParticlePosition, const FParticleType /*particleType*/, Args... args){
push(inParticlePosition, args...);
}
template<typename... Args>
void push(const FPoint& inParticlePosition, const FParticleType /*particleType*/, Args... args){
push(inParticlePosition, args...);
}
/** set nb particles to 0 */
void clear(){
nbParticles = 0;
}
/** set nb particles to 0 */
void clear(){
nbParticles = 0;
}
/** to enable rearranging
/** to enable rearranging
* indexesToRemove must be sorted
* it removes all the particles at position indexesToRemove
*/
void removeParticles(const int indexesToRemove[], const int nbParticlesToRemove){
int offset = 1;
int idxIndexes = 1;
int idxIns = indexesToRemove[0] + 1;
for( ; idxIns < nbParticles && idxIndexes < nbParticlesToRemove ; ++idxIns){
if( idxIns == indexesToRemove[idxIndexes] ){
idxIndexes += 1;
offset += 1;
}
else{
for(int idx = 0 ; idx < 3 ; ++idx){
positions[idx][idxIns-offset] = positions[idx][idxIns];
}
for(unsigned idx = 0 ; idx < NbAttributesPerParticle ; ++idx){
attributes[idx][idxIns-offset] = attributes[idx][idxIns];
}
}
}
for( ; idxIns < nbParticles ; ++idxIns){
for(int idx = 0 ; idx < 3 ; ++idx){
positions[idx][idxIns-offset] = positions[idx][idxIns];
}
for(unsigned idx = 0 ; idx < NbAttributesPerParticle ; ++idx){
attributes[idx][idxIns-offset] = attributes[idx][idxIns];
}
void removeParticles(const int indexesToRemove[], const int nbParticlesToRemove){
int offset = 1;
int idxIndexes = 1;
int idxIns = indexesToRemove[0] + 1;
for( ; idxIns < nbParticles && idxIndexes < nbParticlesToRemove ; ++idxIns){
if( idxIns == indexesToRemove[idxIndexes] ){
idxIndexes += 1;
offset += 1;
}
else{
for(int idx = 0 ; idx < 3 ; ++idx){
positions[idx][idxIns-offset] = positions[idx][idxIns];
}
for(unsigned idx = 0 ; idx < NbAttributesPerParticle ; ++idx){
attributes[idx][idxIns-offset] = attributes[idx][idxIns];
}
}
}
for( ; idxIns < nbParticles ; ++idxIns){
for(int idx = 0 ; idx < 3 ; ++idx){
positions[idx][idxIns-offset] = positions[idx][idxIns];
}
for(unsigned idx = 0 ; idx < NbAttributesPerParticle ; ++idx){
attributes[idx][idxIns-offset] = attributes[idx][idxIns];
}
}
nbParticles -= nbParticlesToRemove;
}
nbParticles -= nbParticlesToRemove;
}
/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
/** The size to send a leaf */
int getSavedSize() const{
return int(sizeof(nbParticles) + nbParticles * (3 * sizeof(FReal) + NbAttributesPerParticle * sizeof(AttributeClass)));
}
/** Save the current cell in a buffer */
template <class BufferWriterClass>
void save(BufferWriterClass& buffer) const{
buffer << nbParticles;
for(int idx = 0 ; idx < 3 ; ++idx){
buffer.write(positions[idx], nbParticles);
}
for(unsigned idx = 0 ; idx < NbAttributesPerParticle ; ++idx){
buffer.write(attributes[idx], nbParticles);
/////////////////////////////////////////////////////
/////////////////////////////////////////////////////
/** The size to send a leaf */
int getSavedSize() const{
return int(sizeof(nbParticles) + nbParticles * (3 * sizeof(FReal) + NbAttributesPerParticle * sizeof(AttributeClass)));
}
}
/** Restore the current cell from a buffer */
template <class BufferReaderClass>
void restore(BufferReaderClass& buffer){
buffer >> nbParticles;
if( nbParticles >= allocatedParticles ){
// allocate memory
const int moduloParticlesNumber = (32/sizeof(FReal)); // We want to be rounded to 32B
allocatedParticles = (FMath::Max(10,int(FReal(nbParticles+1)*1.5)) + moduloParticlesNumber - 1) & ~(moduloParticlesNumber-1);
// init with 0
const size_t allocatedBytes = (sizeof(FReal)*3 + sizeof(AttributeClass)*NbAttributesPerParticle)*allocatedParticles;
FReal* newData = reinterpret_cast<FReal*>(FAlignedMemory::Allocate32BAligned(allocatedBytes));
memset( newData, 0, allocatedBytes);
FAlignedMemory::Dealloc32BAligned(positions[0]);
for(int idx = 0 ; idx < 3 ; ++idx){
positions[idx] = newData + (allocatedParticles * idx);
}
AttributeClass* startAddress = reinterpret_cast<AttributeClass*>(positions[2] + allocatedParticles);
for(unsigned idx = 0 ; idx < NbAttributesPerParticle ; ++idx){
attributes[idx] = startAddress + (idx * allocatedParticles);
}
/** Save the current cell in a buffer */
template <class BufferWriterClass>
void save(BufferWriterClass& buffer) const{
buffer << nbParticles;
for(int idx = 0 ; idx < 3 ; ++idx){
buffer.write(positions[idx], nbParticles);
}
for(unsigned idx = 0 ; idx < NbAttributesPerParticle ; ++idx){
buffer.write(attributes[idx], nbParticles);
}
}
for(int idx = 0 ; idx < 3 ; ++idx){
buffer.fillArray(positions[idx], nbParticles);
/** Restore the current cell from a buffer */
template <class BufferReaderClass>
void restore(BufferReaderClass& buffer){
buffer >> nbParticles;
if( nbParticles >= allocatedParticles ){
// allocate memory
const int moduloParticlesNumber = (32/sizeof(FReal)); // We want to be rounded to 32B
allocatedParticles = (FMath::Max(10,int(FReal(nbParticles+1)*1.5)) + moduloParticlesNumber - 1) & ~(moduloParticlesNumber-1);
// init with 0
const size_t allocatedBytes = (sizeof(FReal)*3 + sizeof(AttributeClass)*NbAttributesPerParticle)*allocatedParticles;
FReal* newData = reinterpret_cast<FReal*>(FAlignedMemory::Allocate32BAligned(allocatedBytes));
memset( newData, 0, allocatedBytes);
FAlignedMemory::Dealloc32BAligned(positions[0]);
for(int idx = 0 ; idx < 3 ; ++idx){
positions[idx] = newData + (allocatedParticles * idx);
}
AttributeClass* startAddress = reinterpret_cast<AttributeClass*>(positions[2] + allocatedParticles);
for(unsigned idx = 0 ; idx < NbAttributesPerParticle ; ++idx){
attributes[idx] = startAddress + (idx * allocatedParticles);
}
}
for(int idx = 0 ; idx < 3 ; ++idx){
buffer.fillArray(positions[idx], nbParticles);
}
for(unsigned idx = 0 ; idx < NbAttributesPerParticle ; ++idx){
buffer.fillArray(attributes[idx], nbParticles);
}
}
for(unsigned idx = 0 ; idx < NbAttributesPerParticle ; ++idx){
buffer.fillArray(attributes[idx], nbParticles);
/** Reset the attributes to zeros */
void resetToInitialState(){
for(unsigned idx = 0 ; idx < NbAttributesPerParticle ; ++idx){
memset(attributes[idx], 0, sizeof(AttributeClass) * allocatedParticles);
}
}
}
/** Reset the attributes to zeros */
void resetToInitialState(const int idxAttribute){
memset(attributes[idxAttribute], 0, sizeof(AttributeClass) * allocatedParticles);
}
};
......
......@@ -62,6 +62,12 @@ public:
const FReal* getForcesZ(const int idxLhs = 0) const {
return Parent::getAttribute(NRHS+3*NLHS+idxLhs);
}
void resetForcesAndPotential(){
for(int idx = 0 ; idx < 4*NLHS ; ++idx){
Parent::resetToInitialState(idx + NRHS);
}
}
};
#endif // FP2PPARTICLECONTAINER_HPP
......@@ -59,6 +59,10 @@
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
//Classical Spherical kernel
#include "../../Src/Kernels/Spherical/FSphericalCell.hpp"
#include "../../Src/Kernels/Spherical/FSphericalKernel.hpp"
//Rotation kernel
#include "Kernels/Rotation/FRotationKernel.hpp"
#include "Kernels/Rotation/FRotationCell.hpp"
......@@ -442,7 +446,81 @@ int main(int argc, char* argv[])
} // end Lagrange/Uniform Grid kernel
#endif
//
// Spherical approximation
//
{
//const static int P = 10;
typedef FSphericalCell CellClass;
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass;
typedef FOctree< CellClass, ContainerClass , LeafClass > OctreeClass;
typedef FSphericalKernel< CellClass, ContainerClass > KernelClass;
typedef FFmmAlgorithmThread<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
std::cout << "\nFFmaSpherical FMM ... P: " << DevP<< std::endl;
{ // -----------------------------------------------------
time.tic();
for(int idxPart = 0 ; idxPart < nbParticles; ++idxPart){
// put in tree
tree.insert(particles[idxPart].position, idxPart, particles[idxPart].physicalValue);
}
time.tac();
std::cout << "(FFmaSpherical @Inserting Particles = "<< time.elapsed() << " s)." << std::endl;
} // -----------------------------------------------------
// -----------------------------------------------------
time.tic();
CellClass::Init(DevP);
KernelClass *kernels = new KernelClass(DevP, TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
FmmClass algorithm(&tree, kernels);
algorithm.execute();
time.tac();
std::cout << "(FFmaSpherical @Algorithm = " << time.elapsed() << " s)." << std::endl;
// -----------------------------------------------------
FReal energy = 0.0;
FMath::FAccurater potentialDiff;
FMath::FAccurater fx, fy, fz;
{ // Check that each particle has been summed with all other
tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FReal*const potentials = leaf->getTargets()->getPotentials();
const FReal*const forcesX = leaf->getTargets()->getForcesX();
const FReal*const forcesY = leaf->getTargets()->getForcesY();
const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
const int nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
const FVector<int>& indexes = leaf->getTargets()->getIndexes();
for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
const int indexPartOrig = indexes[idxPart];
potentialDiff.add(particles[indexPartOrig].potential,potentials[idxPart]);
fx.add(particles[indexPartOrig].forces[0],forcesX[idxPart]);
fy.add(particles[indexPartOrig].forces[1],forcesY[idxPart]);
fz.add(particles[indexPartOrig].forces[2],forcesZ[idxPart]);
energy += potentials[idxPart]*physicalValues[idxPart] ;
}
});
}
// Print for information
std::cout << "FFmaSpherical Energy " << FMath::Abs(energy-energyD) /energyD << std::endl;
std::cout << "FFmaSpherical Potential " << potentialDiff << std::endl;
std::cout << "FFmaSpherical Fx " << fx << std::endl;
std::cout << "FFmaSpherical Fy " << fy << std::endl;
std::cout << "FFmaSpherical Fz " << fz << std::endl;
}
//
// Spherical Rotation
//
{
const static int P = 18;
typedef FRotationCell<P> CellClass;
......
......@@ -37,7 +37,7 @@ MESSAGE(STATUS "CMAKE_CXX_COMPILER_ID ${CMAKE_CXX_COMPILER_ID}")
ELSE()
SET(GCOV_CXX_FLAGS "")
SET(GCOV_LINKER_FLAGS "")
MESSAGE(WARNING "Could not find gcov - No Coverage option for Tests")
MESSAGE(STATUS "Could not find gcov in your path - No Coverage option for Tests")
ENDIF()