Commit 71a2140e authored by BLANCHARD Pierre's avatar BLANCHARD Pierre
Browse files

Merge branch 'master' into PB_tensorialFMM

parents 731ae45d 7dcb5cc7
@article{Schmidt1991,
abstract = {The Rokhlin-Greengard fast multipole algorithm for evaluating Coulomb and multipole potentials has been implemented and analyzed in three dimensions. The implementation is presented for bounded charged systems and systems with periodic boundary conditions. The results include timings and error characterizations},
author = {Schmidt, K. E. and Lee, Michael a.},
doi = {10.1007/BF01030008},
file = {:Users/coulaud/Recherche/Bibliographie/Mendeley/Schmidt, Lee - 1991 - Implementing the fast multipole method in three dimensions.pdf:pdf},
issn = {0022-4715},
journal = {Journal of Statistical Physics},
keywords = {Ewald summation,Periodic boundary conditions,fast multipole method,macro summation,many-body problem,n-body,spherical harmonique},
mendeley-groups = {FMM},
month = jun,
number = {5-6},
pages = {1223--1235},
title = {{Implementing the fast multipole method in three dimensions}},
url = {http://link.springer.com/10.1007/BF01030008},
volume = {63},
year = {1991}
}
@article{Heyes1981,
abstract = {The Coulomb potentials and fields in infinite point charge lattices are represented by series expansions in real and reciprocal space using a method similar to that devised by Bertaut. A systematic investigation is made into the relative convergence characteristics of series derived by varying the charge spreading function which is required in the theory. A number of formulas which are already in use are compared with new expressions derived by this method.},
author = {Heyes, D. M.},
doi = {10.1063/1.441285},
file = {:Users/coulaud/Recherche/Bibliographie/Mendeley/Heyes - 1981 - Electrostatic potentials and fields in infinite point charge lattices.pdf:pdf},
issn = {00219606},
journal = {The Journal of Chemical Physics},
keywords = {Electrostatics,Lattice theory},
mendeley-groups = {Chimie},
mendeley-tags = {Electrostatics,Lattice theory},
number = {3},
pages = {1924},
title = {{Electrostatic potentials and fields in infinite point charge lattices}},
url = {http://link.aip.org/link/?JCP/74/1924/1\&Agg=doi},
volume = {74},
year = {1981}
}
@article{DeLeeuw1980,
author = {{De Leeuw}, S W and Perram, J W and Smith, E R},
file = {:Users/coulaud/Recherche/Bibliographie/Mendeley/De Leeuw, Perram, Smith - 1980 - Simulation of electrostatic systems in periodic boundary conditions. I. Lattice sums and dielectric con.pdf:pdf},
journal = {Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences},
keywords = {Periodic Boundary Conditions,dieclectric constant,lattice sums},
mendeley-groups = {Chimie},
mendeley-tags = {Periodic Boundary Conditions,dieclectric constant,lattice sums},
number = {1752},
pages = {27--56},
title = {{Simulation of electrostatic systems in periodic boundary conditions. I. Lattice sums and dielectric constants}},
url = {http://www.jstor.org/stable/2397238 http://rspa.royalsocietypublishing.org/content/373/1752/27.short},
volume = {373},
year = {1980}
}
......@@ -77,19 +77,22 @@
\begin{minipage}{0.4\textwidth}
\begin{flushleft} \large
\emph{Author:}\\
Berenger \textsc{Bramas} % Your name
Berenger \textsc{Bramas} \\% Your name
Olivier \textsc{Coulaud} % Your name
\end{flushleft}
\end{minipage}
~
\begin{minipage}{0.4\textwidth}
\begin{flushright} \large
\emph{Supervisor:} \\
\textsc{} % Supervisor's Name
%\emph{Supervisor:} \\
%\textsc{} % Supervisor's Name
\end{flushright}
\end{minipage}\\[4cm]
{\large \today}\\[3cm] % Date, change the \today to a set date if you want to be precise
{\large \today}\\[1cm] % Date, change the \today to a set date if you want to be precise
{\Large \textbf{Working note}}
\vfill % Fill the rest of the page with whitespace
\end{titlepage}
......@@ -387,6 +390,16 @@ and the force on atom $x_i$
$$
f(x_i) = \frac{q_i }{4 \pi\epsilon_0}\sum_{j=0,i\neq j}^{N}{q_j\frac{x_i-x_j}{\|x_i-x_j\|^3}}
$$
The Ewald method is equivalent to consider a surrounding conductor ($\epsilon=\infty$). The method described in this document and implemented in ScalFMM consider a finite set of the original box surrounded by a vacuum with dielectric constant $\epsilon=1$. The results by these two approaches are different \cite{DeLeeuw1980, Heyes1981} and for a very large sphere surrounding all boxes we have the following equation
\begin{equation}
E_{Ewald} = E_{FMM} -\frac{2\pi}{3V}|\sum_{i=1}^{N}{q_i x_i}|^2.
\end{equation}
By tacking the gradient of this equation, the relationship between FMM forces and Ewald forces is
\begin{equation}
F_{Ewald}(x_k) = -\nabla_k E_{Ewald} = F_{FMM}(x_k) + \frac{4\pi}{3V} q_k \sum_{i=1}^{N}{q_i x_i}.
\end{equation}
\subsubsection{DL\_Poly comparisons}
DL\_POLY\_2 uses the following internal molecular units \\
......@@ -428,12 +441,49 @@ Force &$C_{force}$& -138935.4835 & $10\,J mol^{-1}\, \AA^{-1}$\\
The force unit is the internal DL\_Poly unit.
The first test consists in a small crystal $4\times 4\times 4$ of NaCl. It is composed of 128 atoms The second tests is a larger crystal $10\times 10\times 10$ of NaCl and have X atoms.
The first test consists in a small crystal $4\times 4\times 4$ of NaCl. It is composed of 128 atoms The second tests is a larger crystal $10\times 10\times 10$ of NaCl and have 2000 atoms. The positions and the forces are stored in the \texttt{REVCON} file and the energy in the \texttt{STATIS} file. The Ewald's parameter are chosen such that the Ewald sum precision is $10^{-08}$.
\begin{center}
\begin{tabular}{|l|c|c|c|c|}
\hline
Crystal & Cube size & Cut off & reciprocal lattice vector &Energy \\
& $\AA$ & $\AA$ & & $kcal\, mol^{-1}$ \\
\hline
Small & 16 & 8 & (9,9,9) & $ -2.346403 \,10^{4} $ \\
\hline
Big & 40 & 12 & (15,15,15) & $-3.666255 \,10^{5} $ \\
\hline
\end{tabular}
\end{center}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion}
A new approach has been proposed to compute a periodic FMM.
This approach can be used with any kernel and without the need of developing new operators.
Moreover, it is possible to compute an important repetition size for a small cost.
\bibliographystyle{plain}
\bibliography{fmm}
\newpage
\section*{Appendix 1 }
\begin{center}
\large DL\_poly CONTROL file
\end{center}
\begin{verbatim}
timestep 1.e-15
steps 1
temperature 0.0
ensemble nve
delr 0.5
cutoff 8.0
ewald precision 1.e-8
no vdw
traj 1,1,2
print 1
stats 1
job time 3600
close time 10
finish
\end{verbatim}
\end{document}
......@@ -71,7 +71,8 @@ public :
/** Move the read index to a position */
void seek(const int inIndex){
FAssertLF(inIndex < arrayCapacity, "FMpiBufferReader :: Aborting :: Can't move index because buffer isn't long enough");
FAssertLF(inIndex <= arrayCapacity, "FMpiBufferReader :: Aborting :: Can't move index because buffer isn't long enough ",inIndex," ",arrayCapacity);
currentIndex = inIndex;
}
......
// ===================================================================================
// 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 dune licence propriétaire.
// 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.
// Co-authors: Olivier Coulaud, B��renger Bramas.
// Owners: INRIA.
// Copyright © 2011-2012, spread under the terms and conditions of a proprietardata[1] license.
// Copyright �� 2011-2012, spread under the terms and conditions of a proprietardata[1] license.
// ===================================================================================
#ifndef FPOINT_HPP
#define FPOINT_HPP
......@@ -62,7 +62,7 @@ public:
}
/**
* Assignement operator
* Assignment operator
* @param other the source class to copy
*/
FPoint(const FPoint& other, const FReal addset) {
......@@ -172,14 +172,23 @@ public:
FReal * getDataValue(){
return this->data ;
}
/**
*Compute the distance to the origin
* @return the norm of the Fpoint
*/
FReal norm() const {
return FMath::Sqrt(this->data[0]*this->data[0]+this->data[1]*this->data[1]
+this->data[2]*this->data[2]) ;
}
/**
*Compute the distance to the origin
* @return the norm of the Fpoint
*/
FReal norm() const {
return FMath::Sqrt(this->data[0]*this->data[0]+this->data[1]*this->data[1]
+this->data[2]*this->data[2]) ;
}
/**
*Compute the distance to the origin
* @return the square norm of the Fpoint
*/
FReal norm2() const {
return (this->data[0]*this->data[0]+this->data[1]*this->data[1]
+this->data[2]*this->data[2]) ;
}
/**
* Subtract to all dim the inValue
......@@ -195,7 +204,7 @@ public:
/**
* Affect to all dim the inValue
* @param inValue the value to afect
* @param inValue the value to affect
* @return the current object after being affected
*/
FPoint& operator+=(const FReal inValue){
......
#include "../../Src/Utils/FAssert.hpp"
#include "../../Src/Containers/FTreeCoordinate.hpp"
#include <list>
#include <functional>
template <class CellClass>
class FBlockedTree {
protected:
/**
* @brief The FBlockOfCells class manages the cells in block allocation.
*/
class FBlockOfCells {
/** One header is allocated at the beginning of each block */
struct BlockHeader{
MortonIndex startingIndex;
MortonIndex endingIndex;
int numberOfCellsInBlock;
int blockIndexesTableSize;
};
//< This value is for not used cells
static const int CellIsEmptyFlag = -1;
protected:
//< Pointer to a block memory
unsigned char* memoryBuffer;
//< Pointer to the header inside the block memory
BlockHeader* blockHeader;
//< Pointer to the indexes table inside the block memory
int* blockIndexesTable;
//< Pointer to the cells inside the block memory
CellClass* blockCells;
public:
/**
* @brief FBlockOfCells
* @param inStartingIndex first cell morton index
* @param inEndingIndex last cell morton index + 1
* @param inNumberOfCells total number of cells in the interval (should be <= inEndingIndex-inEndingIndex)
*/
FBlockOfCells(const MortonIndex inStartingIndex, const MortonIndex inEndingIndex, const int inNumberOfCells)
: memoryBuffer(0), blockHeader(0), blockIndexesTable(0), blockCells(0){
// Find the number of cell to allocate in the blocks
const int blockIndexesTableSize = int(inEndingIndex-inStartingIndex);
FAssertLF(inNumberOfCells <= blockIndexesTableSize);
// Total number of bytes in the block
const size_t memoryToAlloc = sizeof(BlockHeader) + (blockIndexesTableSize*sizeof(int)) + (inNumberOfCells*sizeof(CellClass));
// Allocate
memoryBuffer = new unsigned char[memoryToAlloc];
FAssertLF(memoryBuffer);
memset(memoryBuffer, 0, memoryToAlloc);
// Move the pointers to the correct position
blockHeader = reinterpret_cast<BlockHeader*>(memoryBuffer);
blockIndexesTable = reinterpret_cast<int*>(memoryBuffer+sizeof(BlockHeader));
blockCells = reinterpret_cast<CellClass*>(memoryBuffer+sizeof(BlockHeader)+(blockIndexesTableSize*sizeof(int)));
// Init header
blockHeader->startingIndex = inStartingIndex;
blockHeader->endingIndex = inEndingIndex;
blockHeader->numberOfCellsInBlock = inNumberOfCells;
blockHeader->blockIndexesTableSize = blockIndexesTableSize;
// Set all index to not used
for(int idxCellPtr = 0 ; idxCellPtr < blockIndexesTableSize ; ++idxCellPtr){
blockIndexesTable[idxCellPtr] = CellIsEmptyFlag;
}
}
/** Call the destructor of cells and dealloc block memory */
~FBlockOfCells(){
for(int idxCellPtr = 0 ; idxCellPtr < blockHeader->blockIndexesTableSize ; ++idxCellPtr){
if(blockIndexesTable[idxCellPtr] != CellIsEmptyFlag){
(&blockCells[blockIndexesTable[idxCellPtr]])->~CellClass();
}
}
delete[] memoryBuffer;
}
/** The index of the fist cell (set from the constructor) */
MortonIndex getStartingIndex() const {
return blockHeader->startingIndex;
}
/** The index of the last cell + 1 (set from the constructor) */
MortonIndex getEndingIndex() const {
return blockHeader->endingIndex;
}
/** The number of cell (set from the constructor) */
int getNumberOfCellsInBlock() const {
return blockHeader->numberOfCellsInBlock;
}
/** The size of the interval endingIndex-startingIndex (set from the constructor) */
int getSizeOfInterval() const {
return blockHeader->blockIndexesTableSize;
}
/** Return true if inIndex should be located in the current block */
bool isInside(const MortonIndex inIndex) const{
return blockHeader->startingIndex <= inIndex && inIndex < blockHeader->endingIndex;
}
/** Return true if inIndex is located in the current block and is not empty */
bool exists(const MortonIndex inIndex) const {
return isInside(inIndex) && (blockIndexesTable[inIndex-blockHeader->startingIndex] != CellIsEmptyFlag);
}
/** Return the address of the cell if it exists (or NULL) */
CellClass* getCell(const MortonIndex inIndex){
if( exists(inIndex) ) return &blockCells[blockIndexesTable[inIndex-blockHeader->startingIndex]];
else return 0;
}
/** Return the address of the cell if it exists (or NULL) */
const CellClass* getCell(const MortonIndex inIndex) const {
if( exists(inIndex) ) return &blockCells[blockIndexesTable[inIndex-blockHeader->startingIndex]];
else return 0;
}
/** Allocate a new cell by calling its constructor */
template<typename... CellConstructorParams>
void newCell(const MortonIndex inIndex, const int id, CellConstructorParams... args){
FAssertLF(isInside(inIndex));
FAssertLF(!exists(inIndex));
FAssertLF(id < blockHeader->blockIndexesTableSize);
new((void*)&blockCells[id]) CellClass(args...);
blockIndexesTable[inIndex-blockHeader->startingIndex] = id;
}
/** Iterate on each allocated cells */
template<typename... FunctionParams>
void forEachCell(std::function<void(CellClass*)> function, FunctionParams... args){
for(int idxCellPtr = 0 ; idxCellPtr < blockHeader->blockIndexesTableSize ; ++idxCellPtr){
if(blockIndexesTable[idxCellPtr] != CellIsEmptyFlag){
function(&blockCells[blockIndexesTable[idxCellPtr]], args...);
}
}
}
};
//< height of the tree (1 => only the root)
const int treeHeight;
//< max number of cells in a block
const int nbCellsPerBlocks;
//< all the blocks of the tree
std::list<FBlockOfCells*>* cellBlocksPerLevel;
public:
/** This constructor create a blocked octree from a usual octree
* The cell are allocated as in the usual octree (no copy constructor are called!)
* Once allocated each cell receive its morton index and tree coordinate.
* No blocks are allocated at level 0.
*/
template<class OctreeClass>
FBlockedTree(const int inTreeHeight, const int inNbCellsPerBlock, OctreeClass*const inOctreeSrc)
: treeHeight(inTreeHeight), nbCellsPerBlocks(inNbCellsPerBlock), cellBlocksPerLevel(0){
cellBlocksPerLevel = new std::list<FBlockOfCells*>[treeHeight];
// Iterate on the tree and build
typename OctreeClass::Iterator octreeIterator(inOctreeSrc);
octreeIterator.gotoBottomLeft();
// For each level from heigth - 1 to 1
for(int idxLevel = treeHeight-1; idxLevel > 0 ; --idxLevel){
typename OctreeClass::Iterator avoidGotoLeft = octreeIterator;
// For each cell at this level
do {
typename OctreeClass::Iterator blockIteratorInOctree = octreeIterator;
// Move the iterator per nbCellsPerBlocks (or until it cannot move right)
int sizeOfBlock = 1;
while(sizeOfBlock < nbCellsPerBlocks && octreeIterator.moveRight()){
sizeOfBlock += 1;
}
// Create a block with the apropriate parameters
FBlockOfCells*const newBlock = new FBlockOfCells(blockIteratorInOctree.getCurrentGlobalIndex(),
octreeIterator.getCurrentGlobalIndex()+1,
sizeOfBlock);
// Initialize each cell of the block
int cellIdInBlock = 0;
while(cellIdInBlock != sizeOfBlock){
const CellClass*const oldNode = blockIteratorInOctree.getCurrentCell();
newBlock->newCell(oldNode->getMortonIndex(), cellIdInBlock);
CellClass* newNode = newBlock->getCell(oldNode->getMortonIndex());
newNode->setMortonIndex(oldNode->getMortonIndex());
newNode->setCoordinate(oldNode->getCoordinate());
cellIdInBlock += 1;
blockIteratorInOctree.moveRight();
}
// Keep the block
cellBlocksPerLevel[idxLevel].push_back(newBlock);
// If we can move right then add another block
} while(octreeIterator.moveRight());
avoidGotoLeft.moveUp();
octreeIterator = avoidGotoLeft;
}
}
/** This function dealloc the tree by deleting each block */
~FBlockedTree(){
for(int idxLevel = 0 ; idxLevel < treeHeight ; ++idxLevel){
std::list<FBlockOfCells*>& levelBlocks = cellBlocksPerLevel[idxLevel];
for (FBlockOfCells* &block: levelBlocks){
delete block;
}
}
delete[] cellBlocksPerLevel;
}
/////////////////////////////////////////////////////////
// Lambda function to apply to all member
/////////////////////////////////////////////////////////
/**
* @brief forEachLeaf iterate on the leaf and apply the function
* @param function
*/
//void forEachLeaf(std::function<void(LeafClass*)> function){
// TODO
//}
/**
* @brief forEachLeaf iterate on the cell and apply the function
* @param function
*/
void forEachCell(std::function<void(CellClass*)> function){
for(int idxLevel = 0 ; idxLevel < treeHeight ; ++idxLevel){
std::list<FBlockOfCells*>& levelBlocks = cellBlocksPerLevel[idxLevel];
for (FBlockOfCells* &block: levelBlocks){
block->forEachCell(function);
}
}
}
/**
* @brief forEachLeaf iterate on the cell and apply the function
* @param function
*/
void forEachCellWithLevel(std::function<void(CellClass*,const int)> function){
for(int idxLevel = 0 ; idxLevel < treeHeight ; ++idxLevel){
std::list<FBlockOfCells*>& levelBlocks = cellBlocksPerLevel[idxLevel];
for (FBlockOfCells* &block: levelBlocks){
block->forEachCell(function, idxLevel);
}
}
}
/**
* @brief forEachLeaf iterate on the cell and apply the function
* @param function
*/
//void forEachCellLeaf(std::function<void(CellClass*,LeafClass*)> function){
// TODO
//}
};
#include "../../Src/Components/FSimpleLeaf.hpp"
#include "../../Src/Containers/FVector.hpp"
#include "../../Src/Containers/FOctree.hpp"
#include "../../Src/Core/FFmmAlgorithm.hpp"
#include "../../Src/Kernels/P2P/FP2PParticleContainer.hpp"
#include "../../Src/Kernels/Rotation/FRotationKernel.hpp"
#include "../../Src/Kernels/Rotation/FRotationCell.hpp"
#include "../../Src/Utils/FMath.hpp"
#include "../../Src/Utils/FMemUtils.hpp"
#include "../../Src/Utils/FParameters.hpp"
#include "../../Src/Core/FFmmAlgorithm.hpp"
#include "../../Src/Core/FFmmAlgorithmThread.hpp"
#include "../../Src/Core/FFmmAlgorithmTask.hpp"
#include "../../Src/Files/FFmaLoader.hpp"
int main(int argc, char* argv[]){
static const int P = 9;
typedef FRotationCell<P> CellClass;
typedef FP2PParticleContainer ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass;
typedef FOctree< CellClass, ContainerClass , LeafClass > OctreeClass;
typedef FRotationKernel< CellClass, ContainerClass , P> KernelClass;
typedef FBlockedTree< CellClass > BlockedOctreeClass;
FTic counter;
const int NbLevels = FParameters::getValue(argc,argv,"-h", 5);
const int SizeSubLevels = FParameters::getValue(argc,argv,"-sh", 3);
const char* const filename = FParameters::getStr(argc,argv,"-f", "../Data/test20k.fma");
FFmaLoader loader(filename);
FAssertLF(loader.isOpen());
OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
FPoint particlePosition;
FReal physicalValue;
loader.fillParticle(&particlePosition,&physicalValue);
tree.insert(particlePosition, physicalValue );
}
std::cout << "Done " << "(@Creating and Inserting Particles = " << counter.tacAndElapsed() << "s)." << std::endl;
const int blockSize = FParameters::getValue(argc,argv,"-bs", 250);
counter.tic();
BlockedOctreeClass blockedTree(NbLevels, blockSize, &tree);
std::cout << "Done " << "(@Converting the tree = " << counter.tacAndElapsed() << "s). Block size is " << blockSize << "." << std::endl;
blockedTree.forEachCell([&](CellClass * cell){
std::cout << "Cell " << cell->getMortonIndex() << std::endl;
});
return 0;
}
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// 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".
// ===================================================================================
#include "FUTester.hpp"
#include "../Src/Utils/FParameters.hpp"
// compile by g++ utestMorton.cpp -o utestMorton.exe
/**
* This file is a unit test for the FTreeCoordinate class
*/
/** this class test the list container */
class TestParameters : public FUTester<TestParameters> {
void Lower(){
uassert(FParameters::toLower('A') == 'a');
uassert(FParameters::toLower('Z') == 'z');
uassert(FParameters::toLower('a') == 'a');
uassert(FParameters::toLower('z') == 'z');
uassert(FParameters::toLower('m') == 'm');
uassert(FParameters::toLower(';') == ';');
}
void CharsEquals(){