Commit 602599bb authored by COULAUD Olivier's avatar COULAUD Olivier

Add a new example to compare all our implementations of 1/r kernel.

Add flag to print the transfer (M2L) matrix
parent 56fb0060
......@@ -422,21 +422,23 @@ F_{Ewald}(x_k) = -\nabla_k E_{Ewald} = F_{FMM}(x_k) + \frac{4\pi}{3V} q_k \math
\end{equation}
\subsubsection{DL\_Poly comparisons}
DL\_POLY\_2 uses the following internal molecular units \\
\begin{tabular}{|l|c|c|l|}
\subsection{Units}
Consider the following dimensional units
\begin{center}
\begin{tabular}{|l|c|}
\hline
type & Expression & Numerical value & \\
time &$t_0$ & $10^{-12}\;s$& picosecond\\
length & $l_0$ &$10^{-10}\; m $& $\AA$ Angstrom\\
mass & $m_0$ & $ 1.6605402 \; 10^{27}\; kg $& atomic mass unit\\
charge & $q_0$ & $1.60217733 \; 10^{19} Coulombs$& unit of proton charge\\
energy & $E_0 = m_0(l_0/t_0)^2$&$1.6605402 \; 10^{23}\; Joules $ & $10\; J\, mol^{-1}$\\
type & Expression \\
\hline
time &$t_0$ \\
length & $l_0$ \\
mass & $m_0$ \\
charge & $q_0$ \\
energy & $E_0 = m_0(l_0/t_0)^2$\\
\hline
\end{tabular}
\end{center}
In internal variables the energy writes
In adimensional variables the energy writes
$$
U = \frac{q_0^2}{4 \pi\epsilon_0 l_0}\sum_{i=0}^{N}{\sum_{j<i}{\frac{q_i q_j}{\|x_i-x_j\|}}} = C_{energy} \;\sum_{i=0}^{N}{\sum_{j<i}{\frac{q_i q_j}{\|x_i-x_j\|}}}
$$
......@@ -446,9 +448,24 @@ f(x_i) = -\frac{q_0^2}{4 \pi\epsilon_0 l_0^2} q_i \sum_{j=0,i\neq j}^{N}{q_j\fra
= C_{force} \; q_i \sum_{j=0,i\neq j}^{N}{q_j\frac{x_i-x_j}{\|x_i-x_j\|^3}}
$$
The Energy conversion factor is $\displaystyle \gamma_0 = \frac{q_0^2}{4 \pi\epsilon_0 l_0}/E_0 = 138935.4835$. The energy unit is in Joules and if you want $kcal\, mol^{-1}$ unit the the factor becomes $\gamma_0/418.400$.
The Energy conversion factor is $\displaystyle \gamma_0 = \frac{q_0^2}{4 \pi\epsilon_0 l_0}/E_0 $.
\subsubsection{DL\_Poly comparisons}
DL\_POLY\_2 uses the following internal molecular units \\
\begin{tabular}{|l|c|c|l|}
\hline
type & Expression & Numerical value & \\
time &$t_0$ & $10^{-12}\;s$& picosecond\\
length & $l_0$ &$10^{-10}\; m $& $\AA$ Angstrom\\
mass & $m_0$ & $ 1.6605402 \; 10^{27}\; kg $& atomic mass unit\\
charge & $q_0$ & $1.60217733 \; 10^{19} Coulombs$& unit of proton charge\\
energy & $E_0 = m_0(l_0/t_0)^2$&$1.6605402 \; 10^{23}\; Joules $ & $10\; J\, mol^{-1}$\\
\hline
\end{tabular}
The Energy conversion factor is $\displaystyle \gamma_0 = \frac{q_0^2}{4 \pi\epsilon_0 l_0}/E_0 = 138935.4835$. The energy unit is in Joules and if you want $kcal\, mol^{-1}$ unit the the factor becomes $\gamma_0/418.400$.
To compare with the molecular dynamics code, the forces and the energy is multiplied by $C_{force}$ ($C_{energy}$ respectively).
\begin{center}
\begin{tabular}{|l|c|c|c|}
......@@ -476,6 +493,46 @@ test2 & 40 &$10x10x10$ & 2000 & $-3.666255 \,10^{5} $ &$183.312729$\\
\hline
\end{tabular}
\end{center}
\subsubsection{Stamp comparisons}
Stamp uses the International unit system as internal molecular units \\
\begin{tabular}{|l|c|c|l|}
\hline
type & Expression & Numerical value & \\
time &$t_0$ & $\;s$& second \\
length & $l_0$ &$\; m $& $\AA$ meter \\
mass & $m_0$ & $ kg $& atomic mass unit\\
charge & $q_0$ & $C$& Coulombs\\
energy & $E_0 = m_0(l_0/t_0)^2$& $\frac{kg m^2}{s^2}= J$ & $ Joules $\\
\hline
\end{tabular}
\begin{center}
\begin{tabular}{|l|c|c|l|}
\hline
Constant & Numerical value \\
$\epsilon_0$ & $8.854187\; 10^{-12}$\\
$\displaystyle\frac{1}{4\pi\epsilon_0}$& $8.98755179\; 10^{9}$ \\
\hline
\end{tabular}
\end{center}
To compare with the molecular dynamics code, the forces and the energy is multiplied by $C_{force}$ ($C_{energy}$ respectively).
\begin{center}
\begin{tabular}{|l|c|c|c|}
\hline
& Constant & Value & Unit \\
\hline
Energy & $C_{energy}$& $8.98755179\; 10^{9}$ & $J$ \\
Force &$C_{force}$& -$C_{energy}$
& $J m^{-1}$\\
\hline
\end{tabular}
\end{center}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion}
A new approach has been proposed to compute a periodic FMM.
......
......@@ -14,6 +14,7 @@
#include "Utils/FGlobal.hpp"
#include "Utils/FPoint.hpp"
#include "Files/FFmaGenericLoader.hpp"
#include "Files/FDlpolyLoader.hpp"
#include "Utils/FParameters.hpp"
#include "Utils/FGenerateDistribution.hpp"
......@@ -27,6 +28,7 @@
//! <b> General arguments:</b>
//! \param -help (-h) to see the parameters available in this driver
//! \param -fin name: file name to convert (with extension .fma (ascii) or bfma (binary)
//! \param -fdlpoly name file coming from a DLpoly simulation
//!
//! \param -fout name: generic name for files (without extension) and save data
//! with following format in name.fma or name.bfma if -bin is set"
......@@ -38,7 +40,7 @@
//!
//! Transform an ascii file in a binary file
//!
//! changeFormat -filename unitCubeXYZQ100.fma -outfilename unitCubeXYZQ100 -bin
//! changeFormat -fin unitCubeXYZQ100.fma -fout unitCubeXYZQ100 -bin
//
......@@ -48,7 +50,9 @@ void genDistusage() {
<< std::endl;
std::cout << "Options "<< std::endl
<< " -help to see the parameters " << std::endl
<< " Input: only one option is allowed" << std::endl
<< " -fin name: file name to convert (with extension .fma (ascii) or bfma (binary) " <<std::endl
<< " -fdlpoly name: file name to convert with extension (.bin if binary file) " <<std::endl
<< " Output " << std::endl
<< " -fout name: generic name for files (without extension) and save data" <<std::endl
<< " with following format in name.fma or name.bfma if -bin is set" <<std::endl
......@@ -62,109 +66,138 @@ int main(int argc, char ** argv){
genDistusage() ;
exit(-1);
}
const std::string filename(FParameters::getStr(argc,argv,"-fin", "data.fma"));
FSize NbPoints;
FReal * particles = nullptr ;
if (FParameters::existParameter(argc, argv, "-fin")) {
const std::string filename(FParameters::getStr(argc,argv,"-fin", "data.fma"));
FFmaGenericLoader loader(filename);
//
// Allocation
//
FSize NbPoints = loader.getNumberOfParticles();
const unsigned int nbData = loader.getNbRecordPerline() ;
const unsigned int arraySize =nbData*NbPoints;
FReal * particles ;
particles = new FReal[arraySize] ;
std::memset(particles,0,arraySize*sizeof(FReal));
//
// Read Data
int j = 0 ;
for(int idxPart = 0 ; idxPart < NbPoints ;++idxPart, j+=nbData){
// //
loader.fillParticle(&particles[j],nbData);
// std::cout << "idxPart "<< idxPart << " ";
// for (int jj= 0 ; jj<nbData ; ++jj, ++k){
// std::cout << particles[k] << " ";
// }
// std::cout << std::endl;
}
//
/////////////////////////////////////////////////////////////////////////
// Save data
/////////////////////////////////////////////////////////////////////////
//
// Generate file for ScalFMM FMAGenericLoader
//
if(FParameters::existParameter(argc, argv, "-fout")){
std::string name(FParameters::getStr(argc,argv,"-fout", "output"));
std::string ext(".");
if(name.find(ext) !=std::string::npos) {
std::cout << "No file with extension permitted for output name : " << name << std::endl;
exit(-1);
FFmaGenericLoader loader(filename);
//
// Allocation
//
NbPoints = loader.getNumberOfParticles();
const unsigned int nbData = loader.getNbRecordPerline() ;
const unsigned int arraySize =nbData*NbPoints;
//
particles = new FReal[arraySize] ;
std::memset(particles,0,arraySize*sizeof(FReal));
//
// Read Data
int j = 0 ;
for(int idxPart = 0 ; idxPart < NbPoints ;++idxPart, j+=nbData){
// //
loader.fillParticle(&particles[j],nbData);
// std::cout << "idxPart "<< idxPart << " ";
// for (int jj= 0 ; jj<nbData ; ++jj, ++k){
// std::cout << particles[k] << " ";
// }
// std::cout << std::endl;
}
if( FParameters::existParameter(argc, argv, "-bin")){
name += ".bfma";
else if(FParameters::existParameter(argc, argv, "-fdlpoly")){
FDlpolyLoader *loader = nullptr ;
// if(FParameters::existParameter(argc, argv, "-bin")){
// loader = new FDlpolyBinLoader(filenameEwaldIn.c_str());
// }
// else {
// loader = new FDlpolyAsciiLoader(filenameEwaldIn.c_str());
// }
// NbPoints = loader->getNumberOfParticles() ;
// particles = new FReal[arraySize] ;
// std::memset(particles,0,arraySize*sizeof(FReal));
// for(int idxPart = 0 ; idxPart < NbPoints ; ++idxPart){
//
// int index ;
// FPoint P ; FReal t[3];
/// / loader->fillParticle(&P, t, &physicalValue,&index);
// particles[(index-1)*]
//
// totalCharge += physicalValue ;
}
}
else {
name += ".fma";
genDistusage() ;
}
FFmaGenericWriter writer(name) ;
writer.writeHeader( loader.getCenterOfBox(), loader.getBoxWidth() , NbPoints, sizeof(FReal), nbData) ;
writer.writeArrayOfReal(particles, nbData, NbPoints);
}
//
// Generate file for visualization purpose
//
if(FParameters::existParameter(argc, argv, "-visufmt")){
std::string outfilename(FParameters::getStr(argc,argv,"-fout", "output"));
std::string visufile(""), fmt(FParameters::getStr(argc,argv,"-visufmt", "vtp"));
if( fmt == "vtp" ){
visufile = outfilename + ".vtp" ;
}
else if( fmt == "vtk" ){
visufile = outfilename + ".vtk" ;
}
else if( fmt == "cosmo" ){
if(nbData !=4) {
std::cerr << "Cosmos export accept only 4 data per particles. here: "<<nbData<<std::endl;
std::exit(EXIT_FAILURE);
FReal totalPhysical Value ;
//
/////////////////////////////////////////////////////////////////////////
// Save data
/////////////////////////////////////////////////////////////////////////
//
// Generate file for ScalFMM FMAGenericLoader
//
if(FParameters::existParameter(argc, argv, "-fout")){
std::string name(FParameters::getStr(argc,argv,"-fout", "output"));
std::string ext(".");
if(name.find(ext) !=std::string::npos) {
std::cout << "No file with extension permitted for output name : " << name << std::endl;
exit(-1);
}
visufile = outfilename + ".cosmo" ;
}
else {
visufile = outfilename + ".csv" ;
if( FParameters::existParameter(argc, argv, "-bin")){
name += ".bfma";
}
else {
name += ".fma";
}
FFmaGenericWriter writer(name) ;
writer.writeHeader( loader.getCenterOfBox(), loader.getBoxWidth() , NbPoints, sizeof(FReal), nbData) ;
writer.writeArrayOfReal(particles, nbData, NbPoints);
}
std::ofstream file( visufile, std::ofstream::out);
if(!file) {
std::cout << "Cannot open file."<< std::endl;
exit(-1) ;
} //
//
// Export data in cvs format
// Generate file for visualization purpose
//
if( fmt == "vtp" ){
std::cout << "Writes in XML VTP format (visualization) in file "<< visufile <<std::endl ;
if(nbData==4){
exportVTKxml( file, particles, NbPoints) ;
if(FParameters::existParameter(argc, argv, "-visufmt")){
std::string outfilename(FParameters::getStr(argc,argv,"-fout", "output"));
std::string visufile(""), fmt(FParameters::getStr(argc,argv,"-visufmt", "vtp"));
if( fmt == "vtp" ){
visufile = outfilename + ".vtp" ;
}
else if( fmt == "vtk" ){
visufile = outfilename + ".vtk" ;
}
else if( fmt == "cosmo" ){
if(nbData !=4) {
std::cerr << "Cosmos export accept only 4 data per particles. here: "<<nbData<<std::endl;
std::exit(EXIT_FAILURE);
}
visufile = outfilename + ".cosmo" ;
}
else {
exportVTKxml( file, particles, NbPoints,nbData) ;
visufile = outfilename + ".csv" ;
}
std::ofstream file( visufile, std::ofstream::out);
if(!file) {
std::cout << "Cannot open file."<< std::endl;
exit(-1) ;
} //
//
// Export data in cvs format
//
if( fmt == "vtp" ){
std::cout << "Writes in XML VTP format (visualization) in file "<< visufile <<std::endl ;
if(nbData==4){
exportVTKxml( file, particles, NbPoints) ;
}
else {
exportVTKxml( file, particles, NbPoints,nbData) ;
}
}
else if( fmt == "vtk" ){
std::cout << "Writes in VTK format (visualization) in file "<< visufile <<std::endl ;
exportVTK( file, particles, NbPoints,nbData) ;
}
else if( fmt == "cosmo" ){
std::cout << "Writes in COSMO format (visualization) in file "<< visufile <<std::endl ;
exportCOSMOS( file, particles, NbPoints) ;
}
else {
std::cout << "Writes in CVS format (visualization) in file "<<visufile<<std::endl ;
exportCVS( file, particles, NbPoints,nbData) ;
}
}
else if( fmt == "vtk" ){
std::cout << "Writes in VTK format (visualization) in file "<< visufile <<std::endl ;
exportVTK( file, particles, NbPoints,nbData) ;
}
else if( fmt == "cosmo" ){
std::cout << "Writes in COSMO format (visualization) in file "<< visufile <<std::endl ;
exportCOSMOS( file, particles, NbPoints) ;
}
else {
std::cout << "Writes in CVS format (visualization) in file "<<visufile<<std::endl ;
exportCVS( file, particles, NbPoints,nbData) ;
}
}
//
delete particles ;
//
delete particles ;
//
return 1;
}
//
return 1;
}
This diff is collapsed.
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Brenger Bramas
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
......
......@@ -223,7 +223,19 @@ public:
preExpNExp(temporaryMultiSource);
const FReal one[2] = {1.0 , 0.0};
#ifdef DEBUG_SPHERICAL_M2L
//
std::cout << "\n ====== MultipolToLocal MatrixTransfer ====== \n"<<std::endl;
for(int idxRow =0 ;idxRow< FF_MATRIX_ROW_DIM ; ++idxRow ){ // Row
std::cout << "Row="<<idxRow<<" : " ;
for(int idxCol =0 ;idxCol< FF_MATRIX_COLUMN_DIM ; ++idxCol ){ // Col
std::cout << M2L_Outer_transfer[idxCol * FF_MATRIX_ROW_DIM + idxRow] <<" ";
}
std::cout << std::endl ;
}
std::cout << "============================ \n"<<std::endl;
#endif
FBlas::c_gemva(
FF_MATRIX_ROW_DIM,
FF_MATRIX_COLUMN_DIM,
......
......@@ -142,8 +142,9 @@ public:
// work with a local variable
FComplexe L_j_k = local_exp[index_j_k];
// n from 0 to P
for (int n = 0 ; n <= /*devP or*/ Parent::devP-j ; ++n){
// n from 0 to P or do P-j
// for (int n = 0 ; n <= Parent::devP /*or*/ /*Parent::devP-j*/ ; ++n){
for (int n = 0 ; n <= Parent::devP-j ; ++n){ // faster than double height Parent::devP
// O_n^l : here points on the source multipole expansion term of degree n and order |l|
const int index_n = Parent::harmonic.getPreExpRedirJ(n);
......
......@@ -318,7 +318,16 @@ public:
friend inline FPoint operator+(const FPoint& inPosition, const FPoint& inOther){
return FPoint(inPosition.data[0] + inOther.data[0], inPosition.data[1] + inOther.data[1], inPosition.data[2] + inOther.data[2]);
}
/**
* Compare two points
*/
inline bool operator==(const FPoint& pp){
/* do actual comparison */
return this->data[0]==pp.data[0] && this->data[1]==pp.data[1]&& this->data[2]==pp.data[2];
}
inline bool operator!=( const FPoint& rhs)
{return !(*this == rhs);}
/**
* Operator stream FPoint to std::ostream
......
......@@ -14,6 +14,7 @@
// "http://www.gnu.org/licenses".
// ===================================================================================
#define DEBUG_SPHERICAL_M2L
#include <iostream>
#include "../Src/Utils/FGlobal.hpp"
......@@ -47,7 +48,7 @@ class TestSphericalDirect : public FUTester<TestSphericalDirect> {
template < class CellClass, class ContainerClass, class KernelClass, class LeafClass,
class OctreeClass, class FmmClass>
void RunTest(const bool isBlasKernel){
const int DevP = 26;
const int DevP = 5;
// Warning in make test the exec dir it Build/UTests
// Load particles
const int nbParticles = 2;
......@@ -92,16 +93,16 @@ class TestSphericalDirect : public FUTester<TestSphericalDirect> {
/* for(int idxLeafZ = 0 ; idxLeafZ < dimGrid ; ++idxLeafZ)*/{
//std::cout << "Shift : " << idxLeafX << " " << idxLeafY << " " << idxLeafZ << std::endl;
particles[1].position = FPoint(FReal(idxLeafX)*dimLeaf + 3*quarterDimLeaf,
/* particles[1].position = FPoint(FReal(idxLeafX)*dimLeaf + 3*quarterDimLeaf,
FReal(idxLeafY)*dimLeaf + quarterDimLeaf,
FReal(idxLeafZ)*dimLeaf + quarterDimLeaf);
/* particles[1].position = FPoint(FReal(idxLeafX)*dimLeaf + quarterDimLeaf,
2*quarterDimLeaf,
2*quarterDimLeaf);
particles[1].position = FPoint(FReal(idxLeafX)*dimLeaf + 2*quarterDimLeaf,
FReal(idxLeafZ)*dimLeaf + quarterDimLeaf);*/
// particles[1].position = FPoint(FReal(idxLeafX)*dimLeaf + quarterDimLeaf,
// quarterDimLeaf,
// quarterDimLeaf);
particles[1].position = FPoint(FReal(idxLeafX)*dimLeaf + 2*quarterDimLeaf,
quarterDimLeaf,
quarterDimLeaf);
particles[1].position = FPoint(FReal(idxLeafX)*dimLeaf + 2*quarterDimLeaf,
/* particles[1].position = FPoint(FReal(idxLeafX)*dimLeaf + 2*quarterDimLeaf,
quarterDimLeaf,
quarterDimLeaf);*/
......@@ -188,7 +189,7 @@ class TestSphericalDirect : public FUTester<TestSphericalDirect> {
for (int j = 0 ; j <= DevP ; ++j ){
std::cout <<"[" << j << "] ----- level \n";
for (int k=0; k<=j ;++k, ++index_j_k){
std::cout << "[" << k << "] ( " << cell->getLocal()[index_j_k].getReal() << " , " << cell->getLocal()[index_j_k].getImag() << " ) ";
std::cout << "[" << k << "] ( " << cell->getLocal()[index_j_k].getReal() << " , " << cell->getLocal()[index_j_k].getImag() << " ) ";
}
std::cout << "\n";
}
......@@ -262,6 +263,10 @@ class TestSphericalDirect : public FUTester<TestSphericalDirect> {
typedef FOctree< CellClass, ContainerClass , LeafClass > OctreeClass;
typedef FFmmAlgorithm<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
std::cout << std::endl << std::endl << std::endl
<< " $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<<std::endl
<< " $$$$$$$$$ TestSpherical $$$$$$$$$$$$$$$$"<<std::endl
<< " $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<<std::endl;
RunTest< CellClass, ContainerClass, KernelClass, LeafClass,
OctreeClass, FmmClass>(false);
......@@ -311,7 +316,11 @@ class TestSphericalDirect : public FUTester<TestSphericalDirect> {
typedef FOctree< CellClass, ContainerClass , LeafClass > OctreeClass;
typedef FFmmAlgorithm<OctreeClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
//
std::cout << std::endl << std::endl << std::endl
<< " $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<<std::endl
<< " $$$$$$$$$ TestSphericalBlockBlas $$$$$$$$$$$$$$$$"<<std::endl
<< " $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<<std::endl;
RunTest< CellClass, ContainerClass, KernelClass, LeafClass,
OctreeClass, FmmClass>(true);
}
......@@ -326,8 +335,8 @@ class TestSphericalDirect : public FUTester<TestSphericalDirect> {
AddTest(&TestSphericalDirect::TestSpherical,"Test Spherical Kernel");
//AddTest(&TestSphericalDirect::TestRotation,"Test Rotation Spherical Kernel");
#ifdef ScalFMM_USE_BLAS
//AddTest(&TestSphericalDirect::TestSphericalBlas,"Test Spherical Blas Kernel");
//AddTest(&TestSphericalDirect::TestSphericalBlockBlas,"Test Spherical Block Blas Kernel");
// AddTest(&TestSphericalDirect::TestSphericalBlas,"Test Spherical Blas Kernel");
// AddTest(&TestSphericalDirect::TestSphericalBlockBlas,"Test Spherical Block Blas Kernel");
#endif
}
};
......
Improvemnents
Improvements
1) Periodic boundary conditions for cubic box
2)
......@@ -7,7 +7,7 @@ Improvemnents
Bugs
1) Spherical -- check the values of multipole and local coefficients.
use : utestSphericalDirectDebug.cpp
2)
2) Intel compiler : need -fp-model precise -fp-model source -fimf-precision=low for Spherical (?????)
ToDO
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment