Skip to content
GitLab
Menu
Projects
Groups
Snippets
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
solverstack
ScalFMM
Commits
5d5b6d5d
Commit
5d5b6d5d
authored
Mar 10, 2017
by
COULAUD Olivier
Browse files
Mmove Utils/noDist in Obsolete
parent
ae20c7f1
Changes
10
Expand all
Hide whitespace changes
Inline
Side-by-side
Utils/noDist/ChebyshevInterpolationCmpAlgo.cpp
deleted
100644 → 0
View file @
ae20c7f1
// ===================================================================================
// Copyright ScalFmm 2016 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.
// An extension to the license is given to allow static linking of scalfmm
// inside a proprietary application (no matter its license).
// See the main license file for more details.
//
// 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 "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"
#ifdef _OPENMP
#include "Core/FFmmAlgorithmThread.hpp"
#include "Core/FFmmAlgorithmTask.hpp"
#include "Core/FFmmAlgorithmSectionTask.hpp"
#include "Core/FFmmAlgorithmThreadBalance.hpp"
#else
#include "Core/FFmmAlgorithm.hpp"
#endif
#include "Utils/FParameterNames.hpp"
/**
* This program runs the FMM Algorithm with the Chebyshev kernel and compares the results with a direct computation.
*/
/// \file ChebyshevInterpolationCmpAlgo.cpp
//!
//! \brief This program runs the FMM Algorithm with the interpolation kernel based on Chebyshev interpolation (1/r kernel)
//! \authors O. Coulaud
//!
//! This code is a short example to use the Chebyshev Interpolation approach for the 1/r kernel
//
int
main
(
int
argc
,
char
*
argv
[])
{
const
FParameterNames
LocalOptionAlgo
=
{
{
"-algo"
}
,
" Algorithm to run (basic, balanced, task, sectiontask)
\n
"
};
const
FParameterNames
LocalOptionCmp
=
{
{
"-cmp"
}
,
"Use to check the result with the exact solution given in the input file
\n
"
};
FHelpDescribeAndExit
(
argc
,
argv
,
"Driver for Chebyshev interpolation kernel (1/r kernel)."
,
FParameterDefinitions
::
InputFile
,
FParameterDefinitions
::
OctreeHeight
,
FParameterDefinitions
::
OctreeSubHeight
,
FParameterDefinitions
::
InputFile
,
FParameterDefinitions
::
OctreeSubHeight
,
FParameterDefinitions
::
OutputFile
,
FParameterDefinitions
::
NbThreads
,
LocalOptionAlgo
,
LocalOptionCmp
);
/////////////////////////////////////////////////////////////////////////////////////
// Set parameters
/////////////////////////////////////////////////////////////////////////////////////
const
std
::
string
defaultFile
(
SCALFMMDataPath
+
"Data/unitCubeXYZQ100.bfma"
);
const
std
::
string
filename
=
FParameters
::
getStr
(
argc
,
argv
,
FParameterDefinitions
::
InputFile
.
options
,
defaultFile
.
c_str
());
const
unsigned
int
TreeHeight
=
FParameters
::
getValue
(
argc
,
argv
,
FParameterDefinitions
::
OctreeHeight
.
options
,
5
);
const
unsigned
int
SubTreeHeight
=
FParameters
::
getValue
(
argc
,
argv
,
FParameterDefinitions
::
OctreeSubHeight
.
options
,
2
);
const
unsigned
int
NbThreads
=
FParameters
::
getValue
(
argc
,
argv
,
FParameterDefinitions
::
NbThreads
.
options
,
1
);
#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
;
/////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////
typedef
double
FReal
;
//
////////////////////////////////////////////////////////////////////
// begin Chebyshev kernel
// accuracy
const
unsigned
int
ORDER
=
7
;
// typedefs
typedef
FP2PParticleContainerIndexed
<
FReal
>
ContainerClass
;
typedef
FSimpleLeaf
<
FReal
,
ContainerClass
>
LeafClass
;
typedef
FChebCell
<
FReal
,
ORDER
>
CellClass
;
typedef
FOctree
<
FReal
,
CellClass
,
ContainerClass
,
LeafClass
>
OctreeClass
;
//
typedef
FInterpMatrixKernelR
<
FReal
>
MatrixKernelClass
;
const
MatrixKernelClass
MatrixKernel
;
typedef
FChebSymKernel
<
FReal
,
CellClass
,
ContainerClass
,
MatrixKernelClass
,
ORDER
>
KernelClass
;
//
// Different algorithms
#ifdef _OPENMP
typedef
FFmmAlgorithmThread
<
OctreeClass
,
CellClass
,
ContainerClass
,
KernelClass
,
LeafClass
>
ForFmmClass
;
typedef
FFmmAlgorithmTask
<
OctreeClass
,
CellClass
,
ContainerClass
,
KernelClass
,
LeafClass
>
TaskFmmClass
;
typedef
FFmmAlgorithmSectionTask
<
OctreeClass
,
CellClass
,
ContainerClass
,
KernelClass
,
LeafClass
>
SectionTaskFmmClass
;
typedef
FFmmAlgorithmThreadBalance
<
OctreeClass
,
CellClass
,
ContainerClass
,
KernelClass
,
LeafClass
>
ForBalFmmClass
;
#else
typedef
FFmmAlgorithm
<
OctreeClass
,
CellClass
,
ContainerClass
,
KernelClass
,
LeafClass
>
FmmClass
;
#endif
/////////////////////////////////////////////////////////////////////
// open particle file
////////////////////////////////////////////////////////////////////
FFmaGenericLoader
<
FReal
>
loader
(
filename
);
if
(
loader
.
getNbRecordPerline
()
!=
8
){
std
::
cerr
<<
"File should contain 8 data to read (x,y,z,q,p,fx,fy,fz)
\n
"
;
std
::
exit
(
EXIT_FAILURE
);
}
FSize
nbParticles
=
loader
.
getNumberOfParticles
()
;
FmaRWParticle
<
FReal
,
8
,
8
>
*
const
particles
=
new
FmaRWParticle
<
FReal
,
8
,
8
>
[
nbParticles
];
//
std
::
cout
<<
"Creating & Inserting "
<<
nbParticles
<<
" particles ..."
<<
std
::
endl
;
std
::
cout
<<
"
\t
Height : "
<<
TreeHeight
<<
"
\t
sub-height : "
<<
SubTreeHeight
<<
std
::
endl
;
time
.
tic
();
//
// open particle file
////////////////////////////////////////////////////////////////////
//
loader
.
fillParticle
(
particles
,
nbParticles
);
time
.
tac
();
//
// init oct-tree
OctreeClass
tree
(
TreeHeight
,
SubTreeHeight
,
loader
.
getBoxWidth
(),
loader
.
getCenterOfBox
());
//
FReal
energyD
=
0.0
;
/////////////////////////////////////////////////////////////////////////////////////////////////
// Compute direct energy
/////////////////////////////////////////////////////////////////////////////////////////////////
for
(
int
idx
=
0
;
idx
<
nbParticles
;
++
idx
){
tree
.
insert
(
particles
[
idx
].
getPosition
()
,
idx
,
particles
[
idx
].
getPhysicalValue
()
);
energyD
+=
particles
[
idx
].
getPotential
()
*
particles
[
idx
].
getPhysicalValue
()
;
}
std
::
cout
<<
"Done "
<<
"(@Creating and Inserting Particles = "
<<
time
.
elapsed
()
<<
" s) ."
<<
std
::
endl
;
// -----------------------------------------------------
{
// -----------------------------------------------------
std
::
cout
<<
"
\n
Chebyshev FMM (ORDER="
<<
ORDER
<<
") ... "
<<
std
::
endl
;
time
.
tic
();
//
KernelClass
kernels
(
TreeHeight
,
loader
.
getBoxWidth
(),
loader
.
getCenterOfBox
(),
&
MatrixKernel
);
//
// false : dynamic schedule.
int
inUserChunckSize
=
10
;
// To specify the chunck size in the loops (-1 is static, 0 is N/p^2, otherwise i)
std
::
string
algoStr
=
FParameters
::
getStr
(
argc
,
argv
,
"-algo"
,
"basic"
);
ForFmmClass
algo1
(
&
tree
,
&
kernels
,
inUserChunckSize
);
ForBalFmmClass
algo4
(
&
tree
,
&
kernels
);
TaskFmmClass
algo2
(
&
tree
,
&
kernels
);
SectionTaskFmmClass
algo3
(
&
tree
,
&
kernels
);
FAbstractAlgorithm
*
algo
=
nullptr
;
FAlgorithmTimers
*
timer
=
nullptr
;
if
(
"basic"
==
algoStr
)
{
algo
=
&
algo1
;
timer
=
&
algo1
;
}
else
if
(
"balanced"
==
algoStr
)
{
algo
=
&
algo4
;
timer
=
&
algo4
;
}
else
if
(
"task"
==
algoStr
)
{
algo
=
&
algo2
;
timer
=
&
algo2
;
}
else
if
(
"sectiontask"
==
algoStr
)
{
algo
=
&
algo3
;
timer
=
&
algo3
;
}
else
{
std
::
cout
<<
"Unknown algorithm: "
<<
algoStr
<<
std
::
endl
;
return
1
;
}
std
::
cout
<<
"Algorithm to check: "
<<
algoStr
<<
std
::
endl
;
time
.
tic
();
// ---------------------------------------------
// algo->execute(FFmmNearField); // Here the call of the FMM algorithm
// algo->execute(FFmmFarField); // Here the call of the FMM algorithm
algo
->
execute
();
// Here the call of the FMM algorithm
// ---------------------------------------------
time
.
tac
();
std
::
cout
<<
"Timers Far Field
\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
//
//
// -----------------------------------------------------
FReal
energy
=
0.0
;
//
// Loop over all leaves
//
if
(
FParameters
::
existParameter
(
argc
,
argv
,
"-cmp"
)
){
std
::
cout
<<
std
::
endl
<<
" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "
<<
std
::
endl
;
std
::
cout
<<
std
::
scientific
;
std
::
cout
.
precision
(
10
)
;
printf
(
"Compute Diff..."
);
FMath
::
FAccurater
<
FReal
>
potentialDiff
;
FMath
::
FAccurater
<
FReal
>
fx
,
fy
,
fz
,
f
;
{
// Check that each particle has been summed with all other
tree
.
forEachLeaf
([
&
](
LeafClass
*
leaf
){
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
();
const
FReal
*
const
physicalValues
=
leaf
->
getTargets
()
->
getPhysicalValues
();
const
FVector
<
FSize
>&
indexes
=
leaf
->
getTargets
()
->
getIndexes
();
for
(
FSize
idxPart
=
0
;
idxPart
<
nbParticlesInLeaf
;
++
idxPart
){
const
FSize
indexPartOrig
=
indexes
[
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
]);
f
.
add
(
particles
[
indexPartOrig
].
getForces
()[
0
],
forcesX
[
idxPart
]);
f
.
add
(
particles
[
indexPartOrig
].
getForces
()[
1
],
forcesY
[
idxPart
]);
f
.
add
(
particles
[
indexPartOrig
].
getForces
()[
2
],
forcesZ
[
idxPart
]);
energy
+=
potentials
[
idxPart
]
*
physicalValues
[
idxPart
];
}
});
std
::
cout
<<
energy
<<
" "
<<
energyD
<<
std
::
endl
;
delete
[]
particles
;
f
.
setNbElements
(
nbParticles
);
std
::
cout
<<
"FChebSymKernel Energy "
<<
FMath
::
Abs
(
energy
-
energyD
)
<<
" Relative "
<<
FMath
::
Abs
(
energy
-
energyD
)
/
FMath
::
Abs
(
energyD
)
<<
std
::
endl
;
std
::
cout
<<
"FChebSymKernel Potential "
<<
potentialDiff
<<
std
::
endl
;
std
::
cout
<<
"FChebSymKernel Fx "
<<
fx
<<
std
::
endl
;
std
::
cout
<<
"FChebSymKernel Fy "
<<
fy
<<
std
::
endl
;
std
::
cout
<<
"FChebSymKernel Fz "
<<
fz
<<
std
::
endl
;
std
::
cout
<<
"FChebSymKernel F "
<<
f
<<
std
::
endl
;
std
::
cout
<<
std
::
endl
<<
"Energy: "
<<
energy
<<
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
*
particlesW
;
particlesW
=
new
FReal
[
8
*
NbPoints
]
;
memset
(
particlesW
,
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
];
particlesW
[
j
]
=
posX
[
idxPart
]
;
particlesW
[
j
+
1
]
=
posY
[
idxPart
]
;
particlesW
[
j
+
2
]
=
posZ
[
idxPart
]
;
particlesW
[
j
+
3
]
=
physicalValues
[
idxPart
]
;
particlesW
[
j
+
4
]
=
potentials
[
idxPart
]
;
particlesW
[
j
+
5
]
=
forcesX
[
idxPart
]
;
particlesW
[
j
+
6
]
=
forcesY
[
idxPart
]
;
particlesW
[
j
+
7
]
=
forcesZ
[
idxPart
]
;
}
});
writer
.
writeHeader
(
loader
.
getCenterOfBox
(),
loader
.
getBoxWidth
()
,
NbPoints
,
sizeof
(
FReal
),
8
)
;
writer
.
writeArrayOfReal
(
particlesW
,
8
,
NbPoints
);
delete
[]
particles
;
//
std
::
string
name1
(
"output.fma"
);
//
FFmaGenericWriter
<
FReal
>
writer1
(
name1
)
;
writer1
.
writeDistributionOfParticlesFromOctree
(
&
tree
,
NbPoints
)
;
}
return
0
;
}
Utils/noDist/FmmAlgorithmTsm.cpp
deleted
100644 → 0
View file @
ae20c7f1
// ===================================================================================
// Copyright ScalFmm 2016 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.
// An extension to the license is given to allow static linking of scalfmm
// inside a proprietary application (no matter its license).
// See the main license file for more details.
//
// 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".
// ===================================================================================
//
//
// Utils/Release/FmmAlgorithmTsm -N 2 -h 4
// plot "output.txt" using ($1):(log10($3)) w l,"output.txt" using ($1):(log10($4)) w l,"output.txt" using ($1):(log10($5)) w l,"output.txt" using ($1):(log10($6)) w l
// ==== CMAKE =====
// @FUSE_BLAS
// ================
// Keep in private GIT
// @SCALFMM_PRIVATE
#include <iostream>
#include <fstream>
#include <sstream>
#include <cstdio>
#include <cstdlib>
#include "Containers/FOctree.hpp"
#include "Containers/FVector.hpp"
#include "Utils/FMath.hpp"
#include "Utils/FParameters.hpp"
#include "Utils/FPoint.hpp"
#include "Utils/FParameterNames.hpp"
#include "Utils/FTic.hpp"
#include "Utils/FTemplate.hpp"
#include "Components/FTypedLeaf.hpp"
#include "Extensions/FExtendCellType.hpp"
#include "Core/FFmmAlgorithmTsm.hpp"
#include "Core/FFmmAlgorithmThreadTsm.hpp"
#include "Kernels/Chebyshev/FChebCell.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Chebyshev/FChebSymKernel.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
#include "Files/FRandomLoader.hpp"
template
<
class
FReal
>
struct
sourcePart
{
FPoint
<
FReal
>
position
;
double
physicalValue
;
};
/** This program show an example of use of
* the fmm basic algo
* it also check that each particles is impacted each other particles
*/
template
<
class
FReal
,
int
ORDER
>
class
FTestCellTsm
:
public
FChebCell
<
FReal
,
ORDER
>
,
public
FExtendCellType
{
};
// Simply create particles and try the kernels
struct
TempMainStruct
{
template
<
const
unsigned
int
ORDER
>
static
int
Run
(
int
argc
,
char
**
argv
){
// int main(int argc, char ** argv){
const
FParameterNames
LocalOptionPoints
=
{
{
"-MT"
}
,
"Number of target points to build."
};
FHelpDescribeAndExit
(
argc
,
argv
,
"Test FMM TSM (target source model) algorithm by counting the nb of interactions each particle receive."
,
FParameterDefinitions
::
OctreeHeight
,
FParameterDefinitions
::
OctreeSubHeight
,
FParameterDefinitions
::
NbParticles
,
LocalOptionPoints
);
// const int ORDER=4 ;
typedef
double
FReal
;
typedef
FP2PParticleContainerIndexed
<
FReal
>
ContainerClassTyped
;
typedef
FTestCellTsm
<
FReal
,
ORDER
>
CellClassTyped
;
typedef
FTypedLeaf
<
FReal
,
ContainerClassTyped
>
LeafClassTyped
;
typedef
FOctree
<
FReal
,
CellClassTyped
,
ContainerClassTyped
,
LeafClassTyped
>
OctreeClassTyped
;
// typedef FTestKernels< CellClassTyped, ContainerClassTyped > KernelClassTyped;
typedef
FInterpMatrixKernelR
<
FReal
>
MatrixKernelClass
;
const
MatrixKernelClass
MatrixKernel
;
typedef
FChebSymKernel
<
FReal
,
CellClassTyped
,
ContainerClassTyped
,
MatrixKernelClass
,
ORDER
>
KernelClassTyped
;
typedef
FFmmAlgorithmThreadTsm
<
OctreeClassTyped
,
CellClassTyped
,
ContainerClassTyped
,
KernelClassTyped
,
LeafClassTyped
>
FmmClassTyped
;
///////////////////////What we do/////////////////////////////
std
::
cout
<<
">> This executable has to be used to test the FMM algorithm.
\n
"
;
//////////////////////////////////////////////////////////////
const
int
TreeHeight
=
FParameters
::
getValue
(
argc
,
argv
,
FParameterDefinitions
::
OctreeHeight
.
options
,
5
);
const
int
SubTreeHeight
=
FParameters
::
getValue
(
argc
,
argv
,
FParameterDefinitions
::
OctreeSubHeight
.
options
,
3
);
const
FSize
NbPart
=
FParameters
::
getValue
(
argc
,
argv
,
FParameterDefinitions
::
NbParticles
.
options
,
FSize
(
20
));
const
std
::
string
filename
=
FParameters
::
getStr
(
argc
,
argv
,
FParameterDefinitions
::
OutputFile
.
options
,
"output"
);
const
int
nbTargets
=
FParameters
::
getValue
(
argc
,
argv
,
LocalOptionPoints
.
options
,
100
);
FTic
counter
;
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
//
std
::
cout
<<
"Parameters "
<<
std
::
endl
<<
" Octree Depth "
<<
TreeHeight
<<
std
::
endl
<<
" SubOctree depth "
<<
SubTreeHeight
<<
std
::
endl
<<
std
::
endl
;
//
//
//
double
boxSize
=
1.0
;
FPoint
<
FReal
>
CentreBox
(
0.5
*
boxSize
,
0.5
*
boxSize
,
0.5
*
boxSize
)
;
double
cellSize
=
1.0
/
FMath
::
pow2
(
TreeHeight
);
std
::
cout
<<
" Cell size: "
<<
cellSize
<<
" "
<<
1.0
/
FMath
::
pow2
(
TreeHeight
-
1
)
<<
std
::
endl
;
const
int
dimGrid
=
(
1
<<
(
TreeHeight
-
1
));
const
FReal
dimLeaf
=
(
boxSize
/
FReal
(
dimGrid
));
const
FReal
quarterDimLeaf
=
(
dimLeaf
/
4.0
);
std
::
cout
<<
"dimLeaf: "
<<
dimLeaf
<<
" quarterDimLeaf: "
<<
quarterDimLeaf
<<
" nbGrid "
<<
dimGrid
<<
std
::
endl
;
//////////////////////////////////////////////////////////////////////////////////
//
// Sources are inside the cell number 1
//
FPoint
<
FReal
>
CentreBoxSource
(
0.5
*
dimLeaf
,
0.5
*
dimLeaf
,
0.5
*
dimLeaf
)
;
FRandomLoaderTsm
<
FReal
>
loader
(
NbPart
,
dimLeaf
,
CentreBoxSource
,
1
);
OctreeClassTyped
tree
(
TreeHeight
,
SubTreeHeight
,
boxSize
,
CentreBox
);
//
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
std
::
cout
<<
"Creating "
<<
NbPart
<<
" particles ..."
<<
std
::
endl
;
counter
.
tic
();
FPoint
<
FReal
>
particlePosition
;
double
physicalValue
=
1.0
;
FSize
nbSRC
=
loader
.
getNumberOfParticles
();
sourcePart
<
FReal
>
*
tabSrc
=
new
sourcePart
<
FReal
>
[
nbSRC
];
{
// Insert sources
FParticleType
particleType
,
source
=
FParticleType
::
FParticleTypeSource
;
for
(
FSize
idxPart
=
0
;
idxPart
<
nbSRC
;
++
idxPart
){
loader
.
fillParticle
(
&
particlePosition
,
&
particleType
);
// std::cout << idxPart << " " << particlePosition << " type " << particleType
// << " physicalValue: " << physicalValue<< std::endl;
tree
.
insert
(
particlePosition
,
source
,
idxPart
,
physicalValue
);
tabSrc
[
idxPart
].
position
=
particlePosition
;
tabSrc
[
idxPart
].
physicalValue
=
physicalValue
;
physicalValue
=-
physicalValue
;
}
}
{
//
// Insert target
//
physicalValue
=
1.0
;
double
dx
=
boxSize
/
(
nbTargets
-
1
)
;
//
std
::
cout
<<
" TARGETS "
<<
std
::
endl
;
FPoint
<
FReal
>
particlePosition2
(
-
dx
,
dimLeaf
+
quarterDimLeaf
,
quarterDimLeaf
);
// int nbTargets = 256;
for
(
FSize
idxPart
=
0
;
idxPart
<
nbTargets
;
++
idxPart
){
particlePosition2
.
incX
(
dx
);
std
::
cout
<<
idxPart
<<
" "
<<
particlePosition2
.
getX
()
/
dimLeaf
<<
" "
<<
particlePosition2
<<
" type "
<<
static_cast
<
int
>
(
FParticleType
::
FParticleTypeTarget
)
<<
" "
<<
physicalValue
<<
std
::
endl
;
tree
.
insert
(
particlePosition2
,
FParticleType
::
FParticleTypeTarget
,
idxPart
,
physicalValue
);
}
}
counter
.
tac
();
std
::
cout
<<
"Done "
<<
"(@Creating and Inserting Particles = "
<<
counter
.
elapsed
()
<<
"s)."
<<
std
::
endl
;
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
std
::
cout
<<
"Working on particles ..."
<<
std
::
endl
;
counter
.
tic
();
//// KernelClassTyped kernels;
KernelClassTyped
kernels
(
TreeHeight
,
boxSize
,
CentreBox
,
&
MatrixKernel
);
//
FmmClassTyped
algo
(
&
tree
,
&
kernels
);
algo
.
execute
();
counter
.
tac
();
std
::
cout
<<
"Done "
<<
"(@Algorithm = "
<<
counter
.
elapsed
()
<<
"s)."
<<
std
::
endl
;
std
::
stringstream
ss
;
ss
<<
filename
<<
"-"
<<
ORDER
<<
".txt"
;
std
::
ofstream
out
(
ss
.
str
(),
std
::
ofstream
::
out
);
std
::
ofstream
outMax
(
"outputMax.txt"
,
std
::
fstream
::
out
|
std
::
fstream
::
app
);
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
double
errorMax
=
0.0
,
errorNorm2
=
0.0
,
d2
=
0.0
;
tree
.
forEachLeaf
([
&
](
LeafClassTyped
*
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
FSize
nbParticlesInLeaf
=
leaf
->
getTargets
()
->
getNbParticles
();
//
for
(
FSize
idxPart
=
0
;
idxPart
<
nbParticlesInLeaf
;
++
idxPart
){
double
pot
=
0.0
,
xx
,
yy
,
zz
;
FPoint
<
FReal
>
FF
{}
;
for
(
int
idxSrc
=
0
;
idxSrc
<
nbSRC
;
++
idxSrc
){
xx
=
posX
[
idxPart
]
-
tabSrc
[
idxSrc
].
position
.
getX
();
yy
=
posY
[
idxPart
]
-
tabSrc
[
idxSrc
].
position
.
getY
();
zz
=
posZ
[
idxPart
]
-
tabSrc
[
idxSrc
].
position
.
getZ
();
FPoint
<
FReal
>
dx
(
xx
,
yy
,
zz
)
;
double
dd
=
1.0
/
dx
.
norm
()
;
double
coeff
=
tabSrc
[
idxSrc
].
physicalValue
*
dd
;
pot
+=
coeff
;
dx
*=
coeff
*
dd
*
dd
;
FF
-=
dx
;
}
d2
+=
pot
*
pot
;
double
d
=
FMath
::
Abs
(
potentials
[
idxPart
]
-
pot
);
errorMax
=
FMath
::
Max
(
errorMax
,
d
);
errorNorm2
+=
d
*
d
;
out
<<
posX
[
idxPart
]
<<
" "
<<
potentials
[
idxPart
]
<<
" "
<<
d
<<
" "
<<
FMath
::
Abs
(
forcesX
[
idxPart
]
-
FF
.
getX
())
<<
" "
<<
FMath
::
Abs
(
forcesY
[
idxPart
]
-
FF
.
getY
())
<<
" "
<<
FMath
::
Abs
(
forcesZ
[
idxPart
]
-
FF
.
getZ
())
<<
std
::
endl
;
}
});
outMax
<<
ORDER
<<
" "
<<
errorMax
<<
" "
<<
FMath
::
Sqrt
(
errorNorm2
/
d2
)
<<
std
::
endl
;
out
.
close
();
return
0
;