Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
solverstack
ScalFMM
Commits
e003316f
Commit
e003316f
authored
Sep 08, 2014
by
PIACIBELLO Cyrille
Browse files
Merge branch 'master' of
git+ssh://scm.gforge.inria.fr//gitroot/scalfmm/scalfmm
parents
7b4af394
2397f7f9
Changes
3
Hide whitespace changes
Inline
Side-by-side
Src/Core/FFmmAlgorithmThreadProc.hpp
View file @
e003316f
...
...
@@ -685,10 +685,7 @@ private:
// Do M2L
//////////////////////////////////////////////////////////////////
KernelClass
*
const
myThreadkernels
=
kernels
[
omp_get_thread_num
()];
const
CellClass
*
neighbors
[
343
];
#pragma omp single nowait
#pragma omp single nowait
{
typename
OctreeClass
::
Iterator
octreeIterator
(
tree
);
octreeIterator
.
moveDown
();
...
...
@@ -718,8 +715,11 @@ private:
{
const
int
chunckSize
=
FMath
::
Max
(
1
,
numberOfCells
/
(
omp_get_num_threads
()
*
omp_get_num_threads
()));
for
(
int
idxCell
=
0
;
idxCell
<
numberOfCells
;
idxCell
+=
chunckSize
){
#pragma omp task
#pragma omp task
{
KernelClass
*
const
myThreadkernels
=
kernels
[
omp_get_thread_num
()];
const
CellClass
*
neighbors
[
343
];
const
int
nbCellToCompute
=
FMath
::
Min
(
chunckSize
,
numberOfCells
-
idxCell
);
for
(
int
idxCellToCompute
=
idxCell
;
idxCellToCompute
<
idxCell
+
nbCellToCompute
;
++
idxCellToCompute
){
const
int
counter
=
tree
->
getInteractionNeighbors
(
neighbors
,
iterArray
[
idxCellToCompute
].
getCurrentGlobalCoordinate
(),
idxLevel
);
...
...
@@ -729,10 +729,10 @@ private:
}
}
#pragma omp taskwait
#pragma omp taskwait
for
(
int
idxThread
=
0
;
idxThread
<
omp_get_num_threads
()
;
++
idxThread
){
#pragma omp task
#pragma omp task
{
kernels
[
idxThread
]
->
finishedLevelM2L
(
idxLevel
);
}
...
...
@@ -1271,8 +1271,6 @@ private:
//////////////////////////////////////////////////////////
{
KernelClass
*
myThreadkernels
=
(
kernels
[
omp_get_thread_num
()]);
#pragma omp single nowait
{
FLOG
(
computationCounter
.
tic
());
...
...
@@ -1286,6 +1284,7 @@ private:
const
int
nbLeavesInTask
=
FMath
::
Min
(
endAtThisShape
-
idxLeafs
,
chunckSize
);
#pragma omp task
{
KernelClass
*
myThreadkernels
=
(
kernels
[
omp_get_thread_num
()]);
// There is a maximum of 26 neighbors
ContainerClass
*
neighbors
[
27
];
...
...
Tests/noDist/ChebyshevPeriodic.cpp
0 → 100755
View file @
e003316f
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas
// 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".
// ===================================================================================
// ==== CMAKE =====
// @FUSE_BLAS
// ================
#include
<iostream>
#include
<cstdio>
#include
<cstdlib>
#include
<string>
#include
"ScalFmmConfig.h"
#include
"Utils/FGlobal.hpp"
#include
"Files/FFmaGenericLoader.hpp"
#include
"Kernels/Chebyshev/FChebCell.hpp"
#include
"Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include
"Kernels/Chebyshev/FChebSymKernel.hpp"
#include
"Components/FSimpleLeaf.hpp"
#include
"Kernels/P2P/FP2PParticleContainerIndexed.hpp"
#include
"Utils/FParameters.hpp"
#include
"Containers/FOctree.hpp"
#include
"Core/FFmmAlgorithmPeriodic.hpp"
template
<
class
Output
>
void
Print
(
const
Output
&
value
){
std
::
cout
<<
"--- Output from program : "
<<
value
<<
"
\n
"
;
}
/**
* This program runs the FMM Algorithm with the Chebyshev kernel and compares the results with a direct computation.
*/
/// \file ChebyshevInterpolationFMM.cpp
//!
//! \brief This program runs the FMM Algorithm with the interpolation kernel based on Chebyshev interpolation (1/r kernel)
//! \authors B. Bramas, O. Coulaud
//!
//! This code is a short example to use the Chebyshev Interpolation approach for the 1/r kernel
//!
//!@Algorithm
//! <b> General arguments:</b>
//! \param -help(-h) to see the parameters available in this driver
//! \param -depth The depth of the octree
//! \param -subdepth Specifies the size of the sub octree
//! \param -t The number of threads
//!
//! \param -f name Name of the particles file with extension (.fma or .bfma). The data in file have to be in our FMA format
//!
//
void
usage
()
{
std
::
cout
<<
"Driver for Chebyshev interpolation kernel (1/r kernel)"
<<
std
::
endl
;
std
::
cout
<<
"Options "
<<
std
::
endl
<<
" -help to see the parameters "
<<
std
::
endl
<<
" -depth the depth of the octree "
<<
std
::
endl
<<
" -subdepth specifies the size of the sub octree "
<<
std
::
endl
<<
" -fin name name specifies the name of the particle distribution"
<<
std
::
endl
<<
" -fout name to write the computation in a file"
<<
std
::
endl
<<
" -t n specifies the number of threads used in the computations"
<<
std
::
endl
;
}
// Simply create particles and try the kernels
int
main
(
int
argc
,
char
*
argv
[])
{
const
std
::
string
defaultFile
(
/*SCALFMMDataPath+*/
"../Data/test20k.fma"
);
const
std
::
string
filename
=
FParameters
::
getStr
(
argc
,
argv
,
"-fin"
,
defaultFile
.
c_str
());
const
std
::
string
filenameOut
=
FParameters
::
getStr
(
argc
,
argv
,
"-fout"
,
"resultPer.fma"
);
const
unsigned
int
TreeHeight
=
FParameters
::
getValue
(
argc
,
argv
,
"-depth"
,
5
);
const
unsigned
int
SubTreeHeight
=
FParameters
::
getValue
(
argc
,
argv
,
"-subdepth"
,
2
);
const
unsigned
int
NbThreads
=
FParameters
::
getValue
(
argc
,
argv
,
"-t"
,
1
);
const
int
PeriodicDeep
=
FParameters
::
getValue
(
argc
,
argv
,
"-per"
,
3
);
if
(
FParameters
::
existParameter
(
argc
,
argv
,
"-h"
)
||
FParameters
::
existParameter
(
argc
,
argv
,
"-help"
)){
usage
()
;
exit
(
EXIT_SUCCESS
);
}
#ifdef _OPENMP
omp_set_num_threads
(
NbThreads
);
std
::
cout
<<
"
\n
>> Using "
<<
omp_get_max_threads
()
<<
" threads.
\n
"
<<
std
::
endl
;
#else
std
::
cout
<<
"
\n
>> Sequential version.
\n
"
<<
std
::
endl
;
#endif
//
std
::
cout
<<
"Parameters "
<<
std
::
endl
<<
" Octree Depth "
<<
TreeHeight
<<
std
::
endl
<<
" SubOctree depth "
<<
SubTreeHeight
<<
std
::
endl
<<
" Input file name: "
<<
filename
<<
std
::
endl
<<
" Thread number: "
<<
NbThreads
<<
std
::
endl
<<
std
::
endl
;
//
// init timer
FTic
time
;
// open particle file
////////////////////////////////////////////////////////////////////
//
FFmaGenericLoader
loader
(
filename
);
FSize
nbParticles
=
loader
.
getNumberOfParticles
()
;
FmaRWParticle
<
8
,
8
>*
const
particles
=
new
FmaRWParticle
<
8
,
8
>
[
nbParticles
];
//
////////////////////////////////////////////////////////////////////
// begin Chebyshev kernel
// accuracy
const
unsigned
int
ORDER
=
7
;
// typedefs
typedef
FP2PParticleContainerIndexed
<>
ContainerClass
;
typedef
FSimpleLeaf
<
ContainerClass
>
LeafClass
;
typedef
FChebCell
<
ORDER
>
CellClass
;
typedef
FOctree
<
CellClass
,
ContainerClass
,
LeafClass
>
OctreeClass
;
//
typedef
FInterpMatrixKernelR
MatrixKernelClass
;
const
MatrixKernelClass
MatrixKernel
;
typedef
FChebSymKernel
<
CellClass
,
ContainerClass
,
MatrixKernelClass
,
ORDER
>
KernelClass
;
//
typedef
FFmmAlgorithmPeriodic
<
OctreeClass
,
CellClass
,
ContainerClass
,
KernelClass
,
LeafClass
>
FmmClass
;
// init oct-tree
OctreeClass
tree
(
TreeHeight
,
SubTreeHeight
,
loader
.
getBoxWidth
(),
loader
.
getCenterOfBox
());
{
// -----------------------------------------------------
std
::
cout
<<
"Creating & Inserting "
<<
loader
.
getNumberOfParticles
()
<<
" particles ..."
<<
std
::
endl
;
std
::
cout
<<
"
\t
Height : "
<<
TreeHeight
<<
"
\t
sub-height : "
<<
SubTreeHeight
<<
std
::
endl
;
time
.
tic
();
//
FPoint
position
;
//
loader
.
fillParticle
(
particles
,
nbParticles
);
for
(
int
idxPart
=
0
;
idxPart
<
loader
.
getNumberOfParticles
()
;
++
idxPart
){
tree
.
insert
(
particles
[
idxPart
].
getPosition
(),
idxPart
,
particles
[
idxPart
].
getPhysicalValue
()
);
}
time
.
tac
();
std
::
cout
<<
"Done "
<<
"(@Creating and Inserting Particles = "
<<
time
.
elapsed
()
<<
" s) ."
<<
std
::
endl
;
}
// -----------------------------------------------------
/////////////////////////////////////////////////////////////////////////////////////////////////
{
// -----------------------------------------------------
std
::
cout
<<
"
\n
Chebyshev FMM (ORDER="
<<
ORDER
<<
") ... "
<<
std
::
endl
;
time
.
tic
();
FmmClass
algo
(
&
tree
,
PeriodicDeep
);
const
MatrixKernelClass
MatrixKernel
;
KernelClass
kernels
(
algo
.
extendedTreeHeight
(),
algo
.
extendedBoxWidth
(),
algo
.
extendedBoxCenter
(),
&
MatrixKernel
);
algo
.
setKernel
(
&
kernels
);
algo
.
execute
();
//
time
.
tac
();
std
::
cout
<<
"Done "
<<
"(@Algorithm = "
<<
time
.
elapsed
()
<<
" s) ."
<<
std
::
endl
;
}
// -----------------------------------------------------
//
// Some output
//
//
FmaRWParticle
<
8
,
8
>*
const
particlesOut
=
new
FmaRWParticle
<
8
,
8
>
[
nbParticles
];
{
// -----------------------------------------------------
//
FReal
energy
=
0.0
,
energyD
=
0.0
;
/////////////////////////////////////////////////////////////////////////////////////////////////
// Compute direct energy
/////////////////////////////////////////////////////////////////////////////////////////////////
for
(
int
idx
=
0
;
idx
<
loader
.
getNumberOfParticles
()
;
++
idx
){
energyD
+=
particles
[
idx
].
getPotential
()
*
particles
[
idx
].
getPhysicalValue
()
;
}
//
// Loop over all leaves
//
std
::
cout
<<
std
::
endl
<<
" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "
<<
std
::
endl
;
std
::
cout
<<
std
::
scientific
;
std
::
cout
.
precision
(
10
)
;
/////////////////////////////////////////////////////////////////////////////////////////////////
// Compare
/////////////////////////////////////////////////////////////////////////////////////////////////
FMath
::
FAccurater
potentialDiff
;
FMath
::
FAccurater
fx
,
fy
,
fz
;
tree
.
forEachLeaf
([
&
](
LeafClass
*
leaf
){
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
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
FReal
*
const
physicalValues
=
leaf
->
getTargets
()
->
getPhysicalValues
();
const
FVector
<
int
>&
indexes
=
leaf
->
getTargets
()
->
getIndexes
();
for
(
int
idxPart
=
0
;
idxPart
<
nbParticlesInLeaf
;
++
idxPart
){
const
int
indexPartOrig
=
indexes
[
idxPart
];
//
particlesOut
[
indexPartOrig
].
setPosition
(
posX
[
idxPart
],
posY
[
idxPart
],
posZ
[
idxPart
])
;
particlesOut
[
indexPartOrig
].
setPhysicalValue
(
physicalValues
[
idxPart
])
;
particlesOut
[
indexPartOrig
].
setPotential
(
potentials
[
idxPart
])
;
particlesOut
[
indexPartOrig
].
setForces
(
forcesX
[
idxPart
],
forcesY
[
idxPart
],
forcesZ
[
idxPart
])
;
potentialDiff
.
add
(
particles
[
indexPartOrig
].
getPotential
(),
potentials
[
idxPart
]);
fx
.
add
(
particles
[
indexPartOrig
].
getForces
()[
0
],
forcesX
[
idxPart
]);
fy
.
add
(
particles
[
indexPartOrig
].
getForces
()[
1
],
forcesY
[
idxPart
]);
fz
.
add
(
particles
[
indexPartOrig
].
getForces
()[
2
],
forcesZ
[
idxPart
]);
energy
+=
potentials
[
idxPart
]
*
physicalValues
[
idxPart
];
}
});
energy
*=
0.5
;
std
::
cout
<<
std
::
endl
<<
"Energy: "
<<
energy
<<
std
::
endl
;
std
::
cout
<<
std
::
endl
<<
" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "
<<
std
::
endl
<<
std
::
endl
;
// remove index
Print
(
"Test1 - Error Relative L2 norm Potential "
);
printf
(
" Pot L2Norm %e
\n
"
,
potentialDiff
.
getL2Norm
());
printf
(
" Pot RL2Norm %e
\n
"
,
potentialDiff
.
getRelativeL2Norm
());
printf
(
" Pot RMSError %e
\n
"
,
potentialDiff
.
getRMSError
());
Print
(
"Fx diff is = "
);
printf
(
" Fx L2Norm %e
\n
"
,
fx
.
getL2Norm
());
printf
(
" Fx RL2Norm %e
\n
"
,
fx
.
getRelativeL2Norm
());
printf
(
" Fx RMSError %e
\n
"
,
fx
.
getRMSError
());
Print
(
"Fy diff is = "
);
printf
(
" Fy L2Norm %e
\n
"
,
fy
.
getL2Norm
());
printf
(
" Fy RL2Norm %e
\n
"
,
fy
.
getRelativeL2Norm
());
printf
(
" Fy RMSError %e
\n
"
,
fy
.
getRMSError
());
Print
(
"Fz diff is = "
);
printf
(
" Fz L2Norm %e
\n
"
,
fz
.
getL2Norm
());
printf
(
" Fz RL2Norm %e
\n
"
,
fz
.
getRelativeL2Norm
());
printf
(
" Fz RMSError %e
\n
"
,
fz
.
getRMSError
());
FReal
L2error
=
(
fx
.
getRelativeL2Norm
()
*
fx
.
getRelativeL2Norm
()
+
fy
.
getRelativeL2Norm
()
*
fy
.
getRelativeL2Norm
()
+
fz
.
getRelativeL2Norm
()
*
fz
.
getRelativeL2Norm
()
);
printf
(
" Total L2 Force Error= %e
\n
"
,
FMath
::
Sqrt
(
L2error
))
;
printf
(
" Energy Error = %.12e
\n
"
,
FMath
::
Abs
(
energy
-
energyD
));
printf
(
" Energy FMM = %.12e
\n
"
,
FMath
::
Abs
(
energy
));
printf
(
" Energy DIRECT = %.12e
\n
"
,
FMath
::
Abs
(
energyD
));
}
// -----------------------------------------------------
if
(
FParameters
::
existParameter
(
argc
,
argv
,
"-fout"
)
){
std
::
cout
<<
"Generate "
<<
filenameOut
<<
" for output file"
<<
std
::
endl
;
//
std
::
cout
<<
" numberofParticles: "
<<
nbParticles
<<
" "
<<
sizeof
(
nbParticles
)
<<
std
::
endl
;
std
::
cout
<<
" Box size: "
<<
loader
.
getBoxWidth
()
<<
" "
<<
sizeof
(
loader
.
getBoxWidth
())
<<
std
::
endl
;
//
FFmaGenericWriter
writer
(
filenameOut
)
;
writer
.
writeHeader
(
loader
.
getCenterOfBox
(),
loader
.
getBoxWidth
()
,
nbParticles
,
*
particlesOut
)
;
writer
.
writeArrayOfParticles
(
particlesOut
,
nbParticles
);
}
delete
[]
particlesOut
;
return
0
;
}
Utils/noDist/removeMoment.cpp
0 → 100755
View file @
e003316f
// ===================================================================================
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger 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".
// ===================================================================================
//
//
//
//
// Tests/Release/removeMoment -Ewald2FMM -fin ../Data/forceNacl_2000_dlpolyPer.bin -dlpoly -fout ../Data/forceNacl_2000.fma
#include
<iostream>
#include
<iomanip>
#include
<cstdio>
#include
<cstdlib>
#include
<cmath>
#include
<algorithm>
#include
"ScalFmmConfig.h"
#include
"Utils/FTic.hpp"
#include
"Utils/FParameters.hpp"
#include
"Files/FDlpolyLoader.hpp"
#include
"Files/FFmaGenericLoader.hpp"
/** DLpoly particle is used in the gadget program
* here we try to make the same simulation
*/
struct
MDParticle
{
FPoint
position
;
FReal
forces
[
3
];
FReal
physicalValue
;
FReal
potential
;
int
index
;
};
// Simply create particles and try the kernels
int
main
(
int
argc
,
char
**
argv
){
//
FReal
scaleEnergy
,
scaleForce
;
const
FReal
r4pie0
=
FReal
(
138935.4835
);
///////////////////////What we do/////////////////////////////
if
(
FParameters
::
existParameter
(
argc
,
argv
,
"-help"
)
||
FParameters
::
existParameter
(
argc
,
argv
,
"-h"
)
){
std
::
cout
<<
">> Remove or Add First moment This executable has to be used to compute direct interaction either for periodic or non periodic system.
\n
"
;
std
::
cout
<<
" Options "
<<
std
::
endl
<<
" -fin fileIn.{bin, txt) input file name "
<<
std
::
endl
<<
" -fout fileout.{bfma, fma) output file name "
<<
std
::
endl
<<
" -Ewald2FMM to add the Electric Moment in the potential and the force "
<<
std
::
endl
<<
" -FMM2Ewald to remove the Electric Moment in the potential and the force "
<<
std
::
endl
;
std
::
cout
<<
" -verbose : print index x y z fx fy fy Q and V"
<<
std
::
endl
;
exit
(
-
1
);
}
if
(
FParameters
::
existParameter
(
argc
,
argv
,
"-Ewald2FMM"
)
&&
FParameters
::
existParameter
(
argc
,
argv
,
"-FMM2Ewald"
)
){
std
::
cerr
<<
" Only -Ewald2FMM or FDMM2Ewald have to be set "
<<
std
::
endl
;
exit
(
EXIT_FAILURE
);
}
///
// ---------------------------------------------------
// DL_POLY CONSTANT
// ---------------------------------------------------
scaleEnergy
=
1.0
;
scaleForce
=
1.0
;
// bool scale = true ;
if
(
FParameters
::
existParameter
(
argc
,
argv
,
"-dlpoly"
)){
// scale = false ;
scaleEnergy
=
r4pie0
/
418.4
;
// kcal mol^{-1}
scaleForce
=
-
r4pie0
;
// 10 J mol^{-1} A^{-1}
}
std
::
cout
<<
"Scaling factor "
<<
std
::
endl
<<
"Energy factor "
<<
scaleEnergy
<<
std
::
endl
<<
" Force factor "
<<
scaleForce
<<
std
::
endl
;
//////////////////////////////////////////////////////////////
//
const
std
::
string
filenameIn
(
FParameters
::
getStr
(
argc
,
argv
,
"-fin"
,
"../Data/forceNacl_128_dlpolyPer.bin"
));
const
std
::
string
filenameOut
(
FParameters
::
getStr
(
argc
,
argv
,
"-fout"
,
"result.fma"
));
// file for -saveError option
FTic
counter
;
// -----------------------------------------------------
// LOADER EWALD PARTICLES FROM DLPOLY
// -----------------------------------------------------
std
::
cout
<<
"Opening : "
<<
filenameIn
<<
"
\n
"
;
FDlpolyLoader
*
loader
=
nullptr
;
std
::
string
ext
(
".bin"
);
// open particle file
if
(
filenameIn
.
find
(
ext
)
!=
std
::
string
::
npos
)
{
loader
=
new
FDlpolyBinLoader
(
filenameIn
.
c_str
());
}
else
if
(
filenameIn
.
find
(
".txt"
)
!=
std
::
string
::
npos
)
{
loader
=
new
FDlpolyAsciiLoader
(
filenameIn
.
c_str
());
}
else
{
std
::
cout
<<
"Input file not allowed only .bin or .txt extensions"
<<
std
::
endl
;
std
::
exit
(
EXIT_FAILURE
)
;
}
if
(
!
loader
->
isOpen
()){
std
::
cout
<<
"Loader Error, "
<<
filenameIn
<<
" is missing
\n
"
;
exit
(
EXIT_FAILURE
);
}
// ---------------------------------------------------------------------------------
// Read particles
// ---------------------------------------------------------------------------------
double
boxsize
[
3
]
;
boxsize
[
0
]
=
boxsize
[
1
]
=
boxsize
[
2
]
=
loader
->
getBoxWidth
()
;
FSize
numberofParticles
=
loader
->
getNumberOfParticles
()
;
std
::
cout
<<
"Creating & Inserting "
<<
numberofParticles
<<
" particles ..."
<<
std
::
endl
;
std
::
cout
<<
"
\t
Width : "
<<
loader
->
getBoxWidth
()
<<
"
\t
center x : "
<<
loader
->
getCenterOfBox
().
getX
()
<<
" y : "
<<
loader
->
getCenterOfBox
().
getY
()
<<
" z : "
<<
loader
->
getCenterOfBox
().
getZ
()
<<
std
::
endl
;
counter
.
tic
();
FPoint
electricMoment
(
0.0
,
0.0
,
0.0
)
;
// const --> then shared
MDParticle
*
const
particles
=
new
MDParticle
[
numberofParticles
];
memset
(
particles
,
0
,
sizeof
(
MDParticle
)
*
numberofParticles
)
;
//
double
totalCharge
=
0.0
;
//
for
(
int
idxPart
=
0
;
idxPart
<
numberofParticles
;
++
idxPart
){
//
loader
->
fillParticle
(
&
particles
[
idxPart
].
position
,
particles
[
idxPart
].
forces
,
&
particles
[
idxPart
].
physicalValue
,
&
particles
[
idxPart
].
index
);
//
totalCharge
+=
particles
[
idxPart
].
physicalValue
;
//
electricMoment
.
incX
(
particles
[
idxPart
].
physicalValue
*
particles
[
idxPart
].
position
.
getX
()
);
electricMoment
.
incY
(
particles
[
idxPart
].
physicalValue
*
particles
[
idxPart
].
position
.
getY
()
);
electricMoment
.
incZ
(
particles
[
idxPart
].
physicalValue
*
particles
[
idxPart
].
position
.
getZ
()
);
}
counter
.
tac
();
std
::
cout
<<
std
::
endl
;
std
::
cout
<<
"Total Charge = "
<<
totalCharge
<<
std
::
endl
;
std
::
cout
<<
"Electric Moment = "
<<
electricMoment
<<
std
::
endl
;
std
::
cout
<<
"Electric Moment norm = "
<<
electricMoment
.
norm2
()
<<
std
::
endl
;
std
::
cout
<<
"----------------------------------------------------"
<<
std
::
endl
;
std
::
cout
<<
std
::
endl
;
std
::
cout
<<
"Done "
<<
"(@Reading Ewald file = "
<<
counter
.
elapsed
()
<<
" s)."
<<
std
::
endl
;
//
//
// ---------------------------------------------------------------------------------
// READ DIRECT COMPUTATION
// ----------------------------------------------------------------
// Save direct computation in binary format
// write binary output file
//
// ----------------------------------------------------------------------------
// Remove or Add Electric Moment
//-----------------------------------------------------------------------------
// remove polarization component
//
double
volume
=
boxsize
[
0
]
*
boxsize
[
1
]
*
boxsize
[
2
]
;
double
coeffCorrection
=
2.0
*
FMath
::
FPi
/
volume
/
3.0
;
double
preScaleEnergy
=
1.0
,
postScaleEnergy
=
1.0
,
preScaleForce
=
1.0
,
postScaleForce
=
1.0
;
if
(
FParameters
::
existParameter
(
argc
,
argv
,
"-Ewald2FMM"
)
){
preScaleEnergy
=
1.0
/
scaleEnergy
;
postScaleEnergy
=
1.0
;
preScaleForce
=
1.0
/
scaleForce
;
postScaleForce
=
1.0
;
}
else
if
(
FParameters
::
existParameter
(
argc
,
argv
,
"-FMM2Ewald"
)
){
preScaleEnergy
=
1.0
;
postScaleEnergy
=
scaleEnergy
;
preScaleForce
=
1.0
/
scaleForce
;
postScaleForce
=
scaleForce
;
coeffCorrection
=
-
coeffCorrection
;
}
else
{
std
::
cout
<<
" -Ewald2FMM ou -FMM2Ewald should be set"
<<
std
::
endl
;
exit
(
EXIT_FAILURE
);
}
std
::
cout
<<
"coeffCorrection: "
<<
coeffCorrection
<<
std
::
endl
;
std
::
cout
<<
"preScaleEnergy: "
<<
preScaleEnergy
<<
std
::
endl
;
std
::
cout
<<
"postScaleEnergy: "
<<
postScaleEnergy
<<
std
::
endl
;
std
::
cout
<<
"preScaleForce: "
<<
preScaleForce
<<
std
::
endl
;
std
::
cout
<<
"postScaleForce: "
<<
postScaleForce
<<
std
::
endl
;
//
double
tmp
,
newEnergy
=
0.0
;
for
(
int
idx
=
0
;
idx
<
numberofParticles
;
++
idx
){
// std::cout << " Pos " << particles[idx].position.getX()<< " "<< particles[idx].position.getY()<< " "<< particles[idx].position.getZ()<< std::endl
// << " F ori " << particles[idx].forces[0]<< " "<< particles[idx].forces[1]<< " "<< particles[idx].forces[2]<< std::endl;
tmp
=
particles
[
idx
].
position
.
getX
()
*
electricMoment
.
getX
()
+
particles
[
idx
].
position
.
getY
()
*
electricMoment
.
getY
()
+
particles
[
idx
].
position
.
getZ
()
*
electricMoment
.
getZ
()
;
//
particles
[
idx
].
potential
+=
2.0
*
tmp
*
coeffCorrection
;
//
particles
[
idx
].
forces
[
0
]
=
preScaleForce
*
particles
[
idx
].
forces
[
0
]
+
2.0
*
particles
[
idx
].
physicalValue
*
coeffCorrection
*
electricMoment
.
getX
()
;
particles
[
idx
].
forces
[
1
]
=
preScaleForce
*
particles
[
idx
].
forces
[
1
]
+
2.0
*
particles
[
idx
].
physicalValue
*
coeffCorrection
*
electricMoment
.
getY
()
;
particles
[
idx
].
forces
[
2
]
=
preScaleForce
*
particles
[
idx
].
forces
[
2
]
+
2.0
*
particles
[
idx
].
physicalValue
*
coeffCorrection
*
electricMoment
.
getZ
()
;
//
particles
[
idx
].
forces
[
0
]
*=
postScaleForce
;
particles
[
idx
].
forces
[
1
]
*=
postScaleForce
;
particles
[
idx
].
forces
[
2
]
*=
postScaleForce
;
// std::cout << " F new " << particles[idx].forces[0]<< " "<< particles[idx].forces[1]<< " "<< particles[idx].forces[2]<< std::endl<< std::endl;
//
}
std
::
cout
<<
std
::
endl
<<
" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "
<<
std
::
endl
;
std
::
cout
<<
std
::
scientific
;
std
::
cout
.
precision
(
10
)
;
std
::
cout
<<
" "
<<
coeffCorrection
*
electricMoment
.
norm2
()
<<
std
::
endl
;
std
::
cout
<<
" "
<<
preScaleEnergy
*
loader
->
getEnergy
()
<<
std
::
endl
;
newEnergy
=
preScaleEnergy
*
loader
->
getEnergy
()
+
coeffCorrection
*
electricMoment
.
norm2
()
;
newEnergy
*=
postScaleEnergy
;
//
std
::
cout
.
precision
(
10
)
;
std
::
cout
<<
"Energy EWALD = "
<<
loader
->
getEnergy
()
<<
std
::
endl
;
std
::
cout
<<
"Energy New = "
<<
newEnergy
<<
std
::
endl
;
//
// Save data in FMAFormat
//FSize nbParticles = loader.getNumberOfParticles() ;
FmaRWParticle
<
8
,
8
>*
const
particlesOut
=
new
FmaRWParticle
<
8
,
8
>
[
numberofParticles
];
// remove index
memset
(
particlesOut
,
0
,
sizeof
(
FmaRWParticle
<
8
,
8
>
)
*
numberofParticles
)
;
for
(
int
idx
=
0
;
idx
<
numberofParticles
;
++
idx
){
particlesOut
[
idx
].
setPosition
(
particles
[
idx
].
position
)
;
particlesOut
[
idx
].
setPhysicalValue
(
particles
[
idx
].
physicalValue
)
;
particlesOut
[
idx
].
setPotential
(
particles
[
idx
].
potential
)
;
particlesOut
[
idx
].
setForces
(
particles
[
idx
].
forces
[
0
],
particles
[
idx
].
forces
[
0
],
particles
[
idx
].
forces
[
2
])
;
}
//
// ----------------------------------------------------------------
// Save computation in binary format
//
//
std
::
cout
<<
"Generate "
<<
filenameOut
<<
" for output file"
<<
std
::
endl
;
//
// std::cout << " numberofParticles: " << nbParticles <<" " << sizeof(numberofParticles) <<std::endl;
// std::cout << " denergy: " << newEnergy <<" " << sizeof(newEnergy) <<std::endl;
// std::cout << " Box size: " << loader.getBoxWidth() << " " << sizeof(loader.getBoxWidth())<<std::endl;
//
FFmaGenericWriter
writer
(
filenameOut
)
;
writer
.
writeHeader
(
loader
->
getCenterOfBox
(),
loader
->
getBoxWidth
()
,
numberofParticles
,
*
particlesOut
)
;
writer
.
writeArrayOfParticles
(
particlesOut
,
numberofParticles
);
//
delete
[]
particles
;
return
0
;
}
Write
Preview