Commit 8b2eab4d authored by COULAUD Olivier's avatar COULAUD Olivier

Add periodic boundary conditions in Exemples ; Check periodic boundary...

Add periodic boundary conditions in Exemples ; Check periodic boundary conditions. ../Tests/Kernels/testSphericalDlpolyAlgorithm.cpp
parent b41c2e5b
Subproject commit 1deddb2781f62dbbf0ee9199f569e49f7346397a
Subproject commit 31198c66ff51f600dde4282faed234a4c0b9cfd9
......@@ -443,7 +443,7 @@ HIDE_SCOPE_NAMES = NO
# will put a list of the files that are included by a file in the documentation
# of that file.
SHOW_INCLUDE_FILES = YES
SHOW_INCLUDE_FILES = NO
# If the FORCE_LOCAL_INCLUDES tag is set to YES then Doxygen
# will list include files with double quotes in the documentation
......
......@@ -16,7 +16,8 @@ INCLUDE_DIRECTORIES(
${SCALFMM_INCLUDES}
)
set(GENERIC_SOURCE_FILES changeFmaFormat generateDistributions statisticsOnOctree DirectComputation CutOffAlgorithm RotationFMM compareAllPoissonKernels)
set(GENERIC_SOURCE_FILES changeFmaFormat compare2Files generateDistributions statisticsOnOctree
DirectComputation CutOffAlgorithm RotationFMM compareAllPoissonKernels)
if(SCALFMM_USE_BLAS)
set(GENERIC_SOURCE_FILES ${GENERIC_SOURCE_FILES} ChebyshevInterpolationFMM ChebyshevInterpolationAdaptiveFMM )
endif(SCALFMM_USE_BLAS)
......@@ -30,11 +31,11 @@ endif(SCALFMM_USE_FFT)
if(SCALFMM_USE_MPI)
set(GENERIC_SOURCE_FILES ${GENERIC_SOURCE_FILES} RotationMPIFMM )
if(SCALFMM_USE_BLAS)
set(GENERIC_SOURCE_FILES ${GENERIC_SOURCE_FILES} ChebyshevInterpolationMPIFMMSplit ChebyshevInterpolationMPIFMM )
set(GENERIC_SOURCE_FILES ${GENERIC_SOURCE_FILES} MPIChebyshevInterpolationFMM ChebyshevInterpolationMPIFMMSplit ChebyshevInterpolationMPIFMM )
endif(SCALFMM_USE_BLAS)
if(SCALFMM_USE_FFT)
# set(GENERIC_SOURCE_FILES ${GENERIC_SOURCE_FILES} )
set(GENERIC_SOURCE_FILES ${GENERIC_SOURCE_FILES} MPIUniformInterpolationFMM)
endif(SCALFMM_USE_FFT)
endif()
......
/*
* genarateDistributions.cpp
*
* Created on: 23 mars 2014
* Author: Olivier Coulaud
*/
#include <iostream>
#include <fstream>
#include <string>
#include <cstdlib>
//
#include "Utils/FGlobal.hpp"
#include "Utils/FPoint.hpp"
#include "Files/FFmaGenericLoader.hpp"
#include "Utils/FParameters.hpp"
#include "Utils/FCompareResults.hpp"
#include "Utils/FParameterNames.hpp"
#include "Adaptive/FBox.hpp"
//
/// \file compare2files.cpp
//!
//! \brief compare2files: Gives the error between
//!
//! Driver to transform a FMA format and/or to build a visualization file<br>
//! For a description of the FMA format see FFmaGenericLoader<br>
//! <b> General arguments:</b>
//! \param -help (-h) to see the parameters available in this driver
//! \param -fmmfile1 name1: first file name to compare (with extension .fma (ascii) or bfma (binary)
//! \param -fmmfile2 name2: second file name to compare (with extension .fma (ascii) or bfma (binary)
//! \param -ewaldfile2 name2 if name2 contains the result done by the ewald method for 1/r kernel.
//!
//! Error code (return of the main)
//! -1 Number of points is different in the two files.
//!
//! \b examples
//!
//!
//! compare2files -file1 unitCubeXYZQ100.fma -file2 unitCubeXYZQ100
using FReal = double;
int main(int argc, char ** argv){
const FParameterNames LocalSortParameter { {"-sort"}, "sort files according to the maximal morton index available"};
FHelpDescribeAndExit(argc, argv,
"Driver to change the format of the input file.",
FParameterDefinitions::InputFileOne, FParameterDefinitions::InputFileTwow,
LocalSortParameter);
const std::string filename1(FParameters::getStr(argc,argv,FParameterDefinitions::InputFileOne.options, "data.fma"));
const std::string filename2(FParameters::getStr(argc,argv,FParameterDefinitions::InputFileTwow.options, "data.fma"));
FFmaGenericLoader<FReal> loader1(filename1);
FFmaGenericLoader<FReal> loader2(filename2);
//
// Allocation
//
FSize nbParticles = loader1.getNumberOfParticles();
const unsigned int nbData = loader1.getNbRecordPerline() ;
if(nbParticles != loader2.getNumberOfParticles()){
std::cerr << "Number of points is different in the two files."<<std::endl ;
return -1 ;
}
if( (nbData != 8) && (loader2.getNbRecordPerline()) ){
std::cerr << "Wrong files only " << std::min(loader2.getNbRecordPerline(),nbData)<< " to read."<<std::endl ;
return -2 ;
}
FmaRWParticle<FReal,8,8>* particles1 = new FmaRWParticle<FReal,8,8>[nbParticles];
FmaRWParticle<FReal,8,8>* particles2 = new FmaRWParticle<FReal,8,8>[nbParticles];
//
loader1.fillParticle(particles1,nbParticles);
loader2.fillParticle(particles2,nbParticles);
if(FParameters::existParameter(argc, argv, LocalSortParameter.options)) {
// define a box, used in the sort
const FBox<FPoint<FReal>> box{loader1.getBoxWidth(),loader1.getCenterOfBox()};
std::cout << "Sort needed !! " << std::endl;
sortArrayWithMortonIndex(box, nbParticles, particles1) ;
sortArrayWithMortonIndex(box, nbParticles, particles2) ;
}
const FSize error = compareTwoArrays<FReal, FmaRWParticle<FReal,8,8>* >("TAG", nbParticles, particles1, particles2);
//
delete[] particles1 ;
delete[] particles2 ;
//
return int(error);
}
......@@ -63,7 +63,8 @@ using KernelClass = FInterpolationKernel<FReal,CellClass,ContainerClass,Matri
#ifdef _OPENMP
#include "Core/FFmmAlgorithmThread.hpp"
//#include "Core/FFmmAlgorithmThread.hpp"
#include "Core/FFmmAlgorithmPeriodic.hpp"
#include "Core/FFmmAlgorithmSectionTask.hpp"
#else
#include "Core/FFmmAlgorithm.hpp"
......@@ -72,7 +73,8 @@ using KernelClass = FInterpolationKernel<FReal,CellClass,ContainerClass,Matri
#ifdef _OPENMP
//using FmmClass = FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
using FmmClass = FFmmAlgorithmSectionTask<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> ;
using FmmClass = FFmmAlgorithmSectionTask<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> ;
using FmmClassPer = FFmmAlgorithmPeriodic<FReal,OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass>;
#else
using FmmClass = FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass>;
#endif
......@@ -89,7 +91,8 @@ int main(int argc, char* argv[])
FParameterDefinitions::OctreeSubHeight,
FParameterDefinitions::InputFile,
FParameterDefinitions::OutputFile,
FParameterDefinitions::NbThreads
FParameterDefinitions::NbThreads,
FParameterDefinitions::PeriodicityNbLevels
);
......@@ -97,13 +100,20 @@ int main(int argc, char* argv[])
const std::string filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, defaultFile.c_str());
const int TreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 5);
const int SubTreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeSubHeight.options, 2);
bool periodicCondition = false ;
if(FParameters::existParameter(argc, argv, FParameterDefinitions::PeriodicityNbLevels.options)){
periodicCondition = true;
}
const unsigned int aboveTree = FParameters::getValue(argc, argv, FParameterDefinitions::PeriodicityNbLevels.options, -2);
#ifdef _OPENMP
const unsigned int NbThreads = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, omp_get_max_threads());
omp_set_num_threads(NbThreads);
const unsigned int NbThreads = periodicCondition ? 1 : FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, omp_get_max_threads());
omp_set_num_threads(NbThreads);
std::cout << "\n>> Using " << NbThreads << " threads.\n" << std::endl;
#else
const int NbThreads = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, 1);
const int NbThreads = 1;
std::cout << "\n>> Sequential version.\n" << std::endl;
#endif
//
......@@ -111,11 +121,15 @@ int main(int argc, char* argv[])
std::string indent(" ");
auto w = std::setw(18);
std::cout << "Parameters" << std::endl << std::left
<< indent << w << "Octree Depth: " << TreeHeight << std::endl
<< indent << w << "SubOctree depth: " << SubTreeHeight << std::endl
<< indent << w << "Input file name: " << filename << std::endl
<< indent << w << "Thread number: " << NbThreads << std::endl
<< std::endl;
<< indent << w << "Octree Depth: " << TreeHeight << std::endl
<< indent << w << "SubOctree depth: " << SubTreeHeight << std::endl;
if(periodicCondition){
std::cout << indent << w << "AboveTree "<< aboveTree <<std::endl;
}
std::cout << indent << w << "Input file name: " << filename << std::endl
<< indent << w << "Thread number: " << NbThreads << std::endl
<< std::endl;
}
//
// init timer
......@@ -158,30 +172,47 @@ int main(int argc, char* argv[])
//
////////////////////////////////////////////////////////////////////
{ // -----------------------------------------------------
// { // -----------------------------------------------------
std::cout << "\n" << interpolationType << " FMM (ORDER= "<< ORDER << ") ... " << std::endl;
const MatrixKernelClass MatrixKernel;
time.tic();
//
std::unique_ptr<KernelClass> kernels(new KernelClass(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel));
std::unique_ptr<KernelClass> kernelsNoPer(new KernelClass(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel));
//
FmmClass algoNoPer(&tree, kernelsNoPer.get());
// periodic FMM algorithm
FmmClassPer algoPer(&tree, aboveTree);
KernelClass kernelsPer(algoPer.extendedTreeHeight(), algoPer.extendedBoxWidth(),
algoPer.extendedBoxCenter(),&MatrixKernel);
algoPer.setKernel(&kernelsPer);
//
FmmClass algo(&tree, kernels.get());
FAbstractAlgorithm * algorithm = nullptr;
FAlgorithmTimers * timer = nullptr;
if(! periodicCondition) { // Non periodic case
algorithm = &algoNoPer ;
timer = &algoNoPer ;
}
else { // Periodic case
algorithm = &algoPer ;
timer = &algoPer ;
}
//
algo.execute(); // Here the call of the FMM algorithm
algorithm->execute(); // Here the call of the FMM algorithm
//
time.tac();
std::cout << "Timers Far Field \n"
<< "P2M " << algo.getTime(FAlgorithmTimers::P2MTimer) << " seconds\n"
<< "M2M " << algo.getTime(FAlgorithmTimers::M2MTimer) << " seconds\n"
<< "M2L " << algo.getTime(FAlgorithmTimers::M2LTimer) << " seconds\n"
<< "L2L " << algo.getTime(FAlgorithmTimers::L2LTimer) << " seconds\n"
<< "P2P and L2P " << algo.getTime(FAlgorithmTimers::NearTimer) << " seconds\n"
<< "P2M " << timer->getTime(FAlgorithmTimers::P2MTimer) << " seconds\n"
<< "M2M " << timer->getTime(FAlgorithmTimers::M2MTimer) << " seconds\n"
<< "M2L " << timer->getTime(FAlgorithmTimers::M2LTimer) << " seconds\n"
<< "L2L " << timer->getTime(FAlgorithmTimers::L2LTimer) << " seconds\n"
<< "P2P and L2P " << timer->getTime(FAlgorithmTimers::NearTimer) << " seconds\n"
<< std::endl;
std::cout << "Done " << "(@Algorithm = " << time.elapsed() << " s) ." << std::endl;
}
//}
// -----------------------------------------------------
//
// Some output
......@@ -195,8 +226,9 @@ int main(int argc, char* argv[])
//
std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl;
std::cout << std::scientific;
std::cout.precision(10) ;
std::cout.precision(15) ;
// std::size_t count =0 , numLeaf =0 ;
FReal TotalPhysicalValue=0.0 ;
tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const posX = leaf->getTargets()->getPositions()[0];
const FReal*const posY = leaf->getTargets()->getPositions()[1];
......@@ -218,62 +250,31 @@ int main(int argc, char* argv[])
<< " Pos "<<posX[idxPart]<<" "<<posY[idxPart]<<" "<<posZ[idxPart]
<< " Forces: " << forcesX[idxPart] << " " << forcesY[idxPart] << " "<< forcesZ[idxPart] <<std::endl;
}
// if(numLeaf==288){
// std::cout << count << " idxPart " << idxPart << potentials[idxPart]<<" "<<
// physicalValues[idxPart] << " energy " << energy << std::endl;
// }
energy += potentials[idxPart]*physicalValues[idxPart] ;
TotalPhysicalValue += physicalValues[idxPart] ;
// ++count ;
}
// std::cout << "numLeaf "<< numLeaf << " " << energy <<" " << count<<std::endl<<std::endl ;
// ++numLeaf ;
});
std::cout <<std::endl<<"Energy: "<< energy<<std::endl;
std::cout <<std::endl<<"aboveRoot: " << aboveTree << " Energy: "<< energy<<" TotalPhysicalValue: " << TotalPhysicalValue<< std::endl;
std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl<<std::endl;
}
//
// -----------------------------------------------------
//
if(FParameters::existParameter(argc, argv, FParameterDefinitions::OutputFile.options)){
std::string name(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options, "output.fma"));
FFmaGenericWriter<FReal> writer(name) ;
//
FSize NbPoints = loader.getNumberOfParticles();
FReal * particles ;
particles = new FReal[8*NbPoints] ;
memset(particles,0,8*NbPoints*sizeof(FReal));
FSize j = 0 ;
tree.forEachLeaf([&](LeafClass* leaf){
//
// Input
const FReal*const posX = leaf->getTargets()->getPositions()[0];
const FReal*const posY = leaf->getTargets()->getPositions()[1];
const FReal*const posZ = leaf->getTargets()->getPositions()[2];
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
//
// Computed data
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 FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
j = 8*indexes[idxPart];
particles[j] = posX[idxPart] ;
particles[j+1] = posY[idxPart] ;
particles[j+2] = posZ[idxPart] ;
particles[j+3] = physicalValues[idxPart] ;
particles[j+4] = potentials[idxPart] ;
particles[j+5] = forcesX[idxPart] ;
particles[j+6] = forcesY[idxPart] ;
particles[j+7] = forcesZ[idxPart] ;
}
});
writer.writeHeader( loader.getCenterOfBox(), loader.getBoxWidth() , NbPoints, sizeof(FReal), 8) ;
writer.writeArrayOfReal(particles, 8 , NbPoints);
writer.writeDataFromOctree(&tree,loader.getNumberOfParticles());
delete[] particles;
//
std::string name1( "output.fma");
//
FFmaGenericWriter<FReal> writer1(name1) ;
writer1.writeDistributionOfParticlesFromOctree(&tree,NbPoints) ;
}
......
......@@ -2,14 +2,10 @@
#ifndef FCORECOMMON_HPP
#define FCORECOMMON_HPP
#include "../Utils/FGlobal.hpp"
#include "../Utils/FAssert.hpp"
#include <vector>
#include "Utils/FGlobal.hpp"
#include "Utils/FAssert.hpp"
#ifdef SCALFMM_USE_EZTRACE
extern "C" {
#include "eztrace.h"
}
#endif
/**
* @brief The FFmmOperations enum
* To chose which operation has to be performed.
......@@ -111,7 +107,7 @@ protected:
int nbLevelsInTree; ///< Height of the tree
void setNbLevelsInTree(const int inNbLevelsInTree){
nbLevelsInTree = inNbLevelsInTree;
nbLevelsInTree = inNbLevelsInTree;
lowerWorkingLevel = nbLevelsInTree;
}
......@@ -154,6 +150,13 @@ public:
validateLevels();
executeCore(operationsToProceed);
}
#ifdef SCALFMM_USE_MPI
/// Build and dill vector of the MortonIndex Distribution at Leaf level
/// p = mpi process id then
/// [mortonLeafDistribution[2*p], mortonLeafDistribution[2*p+1] is the morton index shared by process p
virtual void getAndBuildMortonIndexAtLeaf(std::vector<MortonIndex> & mortonLeafDistribution) {} ;
#endif
};
......
......@@ -5,25 +5,26 @@
#include <omp.h>
#include <algorithm>
//
#include "../Utils/FAssert.hpp"
#include "../Utils/FLog.hpp"
#include "Utils/FGlobal.hpp"
#include "Utils/FAssert.hpp"
#include "Utils/FLog.hpp"
#include "../Utils/FTic.hpp"
#include "Utils/FTic.hpp"
#include "Utils/FAlgorithmTimers.hpp"
#include "../Utils/FGlobal.hpp"
#include "Utils/FGlobal.hpp"
#include "../Containers/FBoolArray.hpp"
#include "../Containers/FOctree.hpp"
#include "../Containers/FLightOctree.hpp"
#include "../Utils/FEnv.hpp"
#include "Utils/FEnv.hpp"
#include "../Containers/FBufferWriter.hpp"
#include "../Containers/FBufferReader.hpp"
#include "../Containers/FVector.hpp"
#include "../Utils/FMpi.hpp"
#include "Utils/FMpi.hpp"
#include <sys/time.h>
#include "FCoreCommon.hpp"
......@@ -55,9 +56,9 @@
template<class OctreeClass, class CellClass, class ContainerClass, class KernelClass, class LeafClass, class P2PExclusionClass = FP2PMiddleExclusion>
class FFmmAlgorithmThreadProc : public FAbstractAlgorithm, public FAlgorithmTimers {
using multipole_t = typename CellClass::multipole_t;
using multipole_t = typename CellClass::multipole_t;
using local_expansion_t = typename CellClass::local_expansion_t;
using symbolic_data_t = CellClass;
using symbolic_data_t = CellClass;
private:
OctreeClass* const tree; ///< The octree to work on
......@@ -104,6 +105,7 @@ private:
return workingIntervalsPerLevel[OctreeHeight * proc + level];
}
/// Does \a procIdx have work at given \a idxLevel
/** i.e. does it hold cells and is responsible of them ? */
bool procHasWorkAtLevel(const int idxLevel , const int idxProc) const {
......@@ -121,11 +123,33 @@ private:
}
public:
/// Get an interval from a process id and tree level
// Interval& getWorkingInterval( int level, int proc){
// return workingIntervalsPerLevel[OctreeHeight * proc + level];
// }
// /// Get an interval from a process id and tree level
// const Interval& getWorkingInterval( int level, int proc) const {
// return workingIntervalsPerLevel[OctreeHeight * proc + level];
// }
/// Get current process interval at given \a level
Interval& getWorkingInterval( int level){
return getWorkingInterval(level, idProcess);
}
/// Build and dill vector of the MortonIndex Distribution at Leaf level
/// p = mpi process id then
/// [mortonLeafDistribution[2*p], mortonLeafDistribution[2*p+1] is the morton index shared by process p
void getAndBuildMortonIndexAtLeaf(std::vector<MortonIndex> & mortonLeafDistribution) final {
for (int p=0 ; p< comm.processCount() ; ++p ){
auto inter = getWorkingInterval(OctreeHeight-1, p );
mortonLeafDistribution[2*p] = inter.leftIndex;
mortonLeafDistribution[2*p+1] = inter.rightIndex;
}
}
const MortonIndex * getPtrOnMortonIndexAtLeaf() {
return &workingIntervalsPerLevel[0].leftIndex ;
}
/// Does the current process has some work at this level ?
bool hasWorkAtLevel( int level){
return idProcess == 0 || (getWorkingInterval(level, idProcess - 1).rightIndex) < (getWorkingInterval(level, idProcess).rightIndex);
......@@ -139,7 +163,8 @@ public:
*/
FFmmAlgorithmThreadProc(const FMpi::FComm& inComm, OctreeClass* const inTree,
KernelClass* const inKernels,
const int inUserChunkSize = 10, const int inLeafLevelSeperationCriteria = 1) :
const int inUserChunkSize = 10,
const int inLeafLevelSeperationCriteria = 1) :
tree(inTree),
kernels(nullptr),
comm(inComm),
......
......@@ -11,11 +11,13 @@
#include <string>
#include <cstdlib>
#include "FAbstractLoader.hpp"
#include "Utils/FGlobal.hpp"
#include "Utils/FAssert.hpp"
#include "FAbstractLoader.hpp"
#include "Utils/FPoint.hpp"
#include "Containers/FVector.hpp"
#include "Containers/FOctree.hpp"
......@@ -36,9 +38,20 @@ class FmaRWParticle {
static_assert(WRITE >= 4, "Cannot create FmaRWParticle with less than 4 as value for WRITE");
/// Data stored
FReal data[WRITE];
FReal data[WRITE]{};
public:
FmaRWParticle() = default;
static const int NBELEMENT = WRITE ;
// FmaRWParticle() = default;
FmaRWParticle() {
memset(data,0,WRITE);
}
FmaRWParticle(const FmaRWParticle& other)
{
std::copy(other.data, other.data + FmaRWParticle::NBELEMENT, data) ;
}
// FmaRWParticle( FmaRWParticle&& o) noexcept : data(std::move(o.data)) { }
/// Get a FPoint<FReal> from the position
FPoint<FReal> getPosition() const{
......@@ -100,6 +113,7 @@ public:
FAssertLF(WRITE>7,"Cannot set Forces[] with WRITE<=7");
return &data[5];
}
/// Set the forces from three values
void setForces(const FReal& inFX, const FReal& inFY, const FReal &inFZ){
FAssertLF(WRITE>7,"Cannot set Forces[] with WRITE<=7");
......@@ -821,7 +835,7 @@ public:
}
//
for(FSize idxPart = 0 , j = 0; idxPart < nbParticlesInLeaf ; ++idxPart, j+=4){
particles[j] = posX[idxPart] ;
particles[j] = posX[idxPart] ;
particles[j+1] = posY[idxPart] ;
particles[j+2] = posZ[idxPart] ;
particles[j+3] = physicalValues[idxPart] ;
......@@ -832,6 +846,62 @@ public:
});
}
/**
* Write all particles (position, physical value, potential and force) inside the octree
*
* @param tree Octree that contains the particles in the leaves
* @param NbPoints number of particles
* @param fileName name of the output file
*
* example
* \code
* OctreeClass tree(TreeHeight, SubTreeHeight, BoxWidth, CenterOfBox);
* ...
* FFmaGenericWriter<FReal> writer(filenameOut) ;
* Fwriter.writeDataFromOctree(&tree, nbParticles);
* \endcode
*/
template <class Toctree>
void writeDataFromOctree( Toctree *tree, const FSize NbPoints){
FReal * particles ;
particles = new FReal[8*NbPoints] ;
memset(particles,0,8*NbPoints*sizeof(FReal));
FSize j = 0 ;
tree->forEachLeaf([&](typename Toctree::LeafClass_T* leaf){
//
// Input
const FReal*const posX = leaf->getTargets()->getPositions()[0];
const FReal*const posY = leaf->getTargets()->getPositions()[1];
const FReal*const posZ = leaf->getTargets()->getPositions()[2];
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
//
// Computed data
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 FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
j = 8*indexes[idxPart];
particles[j] = posX[idxPart] ;
particles[j+1] = posY[idxPart] ;
particles[j+2] = posZ[idxPart] ;
particles[j+3] = physicalValues[idxPart] ;
particles[j+4] = potentials[idxPart] ;
particles[j+5] = forcesX[idxPart] ;
particles[j+6] = forcesY[idxPart] ;
particles[j+7] = forcesZ[idxPart] ;
}
});
this->writeHeader( tree->getBoxCenter(), tree->getBoxWidth() , NbPoints, sizeof(FReal), 8) ;
this->writeArrayOfReal(particles, 8 , NbPoints);
delete[] particles;
}
private:
void writerAscciHeader( const FPoint<FReal> &centerOfBox,const FReal &boxWidth,
......
......@@ -183,7 +183,7 @@ public:
*
*/
template <class OCTREECLASS>
void writeDistributionOfParticlesFromOctree( OCTREECLASS &myOctree, const FSize& nbParticles, const std::vector<MortonIndex> &mortonLeafDistribution){
void writeDistributionOfParticlesFromGroupedOctree( OCTREECLASS &myOctree, const FSize& nbParticles, const std::vector<MortonIndex> &mortonLeafDistribution){
//
// Write the header
int sizeType = 0,ierr = 0 ;
......@@ -274,6 +274,111 @@ public:
}
/**
* Write all for all particles the position, physical values, potential and forces
*
* @param myOctree the octree
* @param nbParticlesnumber of particles
* @param nbLocalParticles number of local particles (on the MPI processus
* @param mortonLeafDistribution the morton distribution of the leaves (this is a vector of size 2* the number of MPI processes
*
*/
template <class OCTREECLASS>
void writeDistributionOfParticlesFromOctree( OCTREECLASS &myOctree,
const FSize& nbParticles,
const FSize& nbLocalParticles,
const std::vector<MortonIndex> &mortonLeafDistribution){
//
// Write the header
int sizeType = 0,ierr = 0 ;
FReal tt{} ;
MPI_Datatype mpiFSize_t = _parallelManager->GetType(nbParticles) ;
MPI_Datatype mpiFReal_t = _parallelManager->GetType(tt) ;
MPI_Type_size(mpiFReal_t, &sizeType) ;
int myRank = _parallelManager->global().processId() ;
_headerSize = 0 ;
//
unsigned int typeFReal[2] = {sizeof(FReal) , static_cast<unsigned int>(_nbDataTowritePerRecord)};
if(myRank==0){
ierr =MPI_File_write_at(_mpiFile, 0, &typeFReal, 2, MPI_INT, MPI_STATUS_IGNORE);
}
MPI_Type_size(MPI_INT, &sizeType) ;
_headerSize += sizeType*2 ;
if(myRank==0){
ierr = MPI_File_write_at(_mpiFile, _headerSize, &nbParticles, 1, mpiFSize_t, MPI_STATUS_IGNORE);
}
MPI_Type_size(mpiFSize_t, &sizeType) ;
_headerSize += sizeType*1 ;
auto centerOfBox =myOctree.getBoxCenter() ;
FReal boxSim[4] = {myOctree.getBoxWidth()*0.5 , centerOfBox.getX() , centerOfBox.getX() , centerOfBox.getX() } ;
if(myRank==0){
ierr =MPI_File_write_at(_mpiFile, _headerSize, &boxSim[0], 4, mpiFReal_t, MPI_STATUS_IGNORE);
}
if(ierr >0){
std::cerr << "Error during the construction of the header in FMpiFmaGenericWriter::writeDistributionOfParticlesFromOctree"<<std::endl;
}
MPI_Type_size(mpiFReal_t, &sizeType) ;
_headerSize += sizeType*4 ;
//
// Construct the local number of particles on my process
FSize maxPartLeaf =0, nn =0;
MortonIndex starIndex = mortonLeafDistribution[2*myRank], endIndex = mortonLeafDistribution[2*myRank+1];
// myOctree.template forEachCellLeaf<typename OCTREECLASS::LeafClass_T >(
myOctree.forEachCellLeaf(
[&](typename OCTREECLASS::CellClassType* cell,
typename OCTREECLASS::LeafClass_T * leaf)
{
if (! (cell->getMortonIndex() < starIndex || cell->getMortonIndex() > endIndex)) {
auto n = leaf->getSrc()->getNbParticles();
maxPartLeaf = std::max(maxPartLeaf,n);
nn += n ;
}
}
);
std::cout << " nn " << nn << " should be " << nbLocalParticles << std::endl;
std::vector<FReal> particles(maxPartLeaf*_nbDataTowritePerRecord);
// Build the offset for eaxh processes
FSize before=0; // Number of particles before me (rank < myrank)
MPI_Scan(&nbLocalParticles,&before,1,mpiFSize_t,MPI_SUM,_parallelManager->global().getComm());
before -= nbLocalParticles ;
MPI_Offset offset = _headerSize + sizeType*_nbDataTowritePerRecord*before;
//
// Write particles in file
myOctree.forEachCellLeaf(
[&](typename OCTREECLASS::CellClassType* cell,
typename OCTREECLASS::LeafClass_T * leaf )
{
if (! (cell->getMortonIndex() < starIndex || cell->getMortonIndex() > endIndex)) {
const FSize nbPartsInLeaf = leaf->getTargets()->getNbParticles();
const FReal*const posX = leaf->getTargets()->getPositions()[0];
const FReal*const posY = leaf->getTargets()->getPositions()[1];
const FReal*const posZ = leaf->getTargets()->getPositions()[2];
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FReal*const forceX = leaf->getTargets()->getForcesX();
const FReal*const forceY = leaf->getTargets()->getForcesY();
const FReal*const forceZ = leaf->getTargets()->getForcesZ();
const FReal*const potential = leaf->getTargets()->getPotentials();
for (int i=0, k=0 ; i < nbPartsInLeaf ;++i,k+=_nbDataTowritePerRecord ) {
particles[k] = posX[i]; particles[k+1] = posY[i]; particles[k+2] = posZ[i];
particles[k+3] = physicalValues[i]; particles[k+4] = potential[i];
particles[k+5] = forceX[i]; particles[k+6] = forceY[i]; particles[k+7] = forceZ[i];
}