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
e4b8d3b4
Commit
e4b8d3b4
authored
Mar 17, 2015
by
BRAMAS Berenger
Browse files
Merge branch 'master' of
git+ssh://scm.gforge.inria.fr//gitroot//scalfmm/scalfmm
parents
e120810c
a4869a85
Changes
5
Hide whitespace changes
Inline
Side-by-side
Src/Adaptive/FAdaptUnifKernel.hpp
View file @
e4b8d3b4
...
...
@@ -30,6 +30,7 @@ class FTreeCoordinate;
// ==== CMAKE =====
// @FUSE_BLAS
// @FUSE_FFT
// ================
// for verbosity only!!!
...
...
@@ -221,7 +222,7 @@ public:
// Target cell: local
const
FReal
localCellWidth
(
KernelBaseClass
::
BoxWidth
/
FReal
(
FMath
::
pow
(
2.0
,
localLevel
)));
const
FPoint
localCellCenter
(
KernelBaseClass
::
getCellCenter
(
local
->
getCoordinate
(),
localLevel
));
std
::
cout
<<
" call P2L localLevel "
<<
localLevel
<<
" localCellCenter "
<<
localCellCenter
<<
std
::
endl
;
//
std::cout << " call P2L localLevel "<< localLevel << " localCellCenter "<< localCellCenter <<std::endl;
// interpolation points of target (X) cell
FPoint
X
[
nnodes
];
FUnifTensor
<
order
>::
setRoots
(
localCellCenter
,
localCellWidth
,
X
);
...
...
@@ -458,7 +459,7 @@ public:
for
(
int
idxContainer
=
0
;
idxContainer
<
nbContainers
;
++
idxContainer
){
counterParticles
+=
particles
[
idxContainer
]
->
getNbParticles
();
}
std
::
cout
<<
" Part("
<<
counterParticles
<<
") "
;
//
std::cout << " Part("<<counterParticles<< ") ";
return
counterParticles
>
this
->
sminM
;
}
};
...
...
Src/Core/FCoreCommon.hpp
View file @
e4b8d3b4
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, B
é
renger Bramas, Matthias Messner
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, B
e
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.
//
...
...
@@ -24,16 +24,16 @@
* To chose which operation has to be performed.
*/
enum
FFmmOperations
{
FFmmP2P
=
(
1
<<
0
),
FFmmP2M
=
(
1
<<
1
),
FFmmP2P
=
(
1
<<
0
),
FFmmP2M
=
(
1
<<
1
),
FFmmM2M
=
(
1
<<
2
),
FFmmM2L
=
(
1
<<
3
),
FFmmL2L
=
(
1
<<
4
),
FFmmL2P
=
(
1
<<
5
),
FFmmM2L
=
(
1
<<
3
),
FFmmL2L
=
(
1
<<
4
),
FFmmL2P
=
(
1
<<
5
),
//
FFmmNearField
=
FFmmP2P
,
FFmmFarField
=
(
FFmmP2M
|
FFmmM2M
|
FFmmM2L
|
FFmmL2L
|
FFmmL2P
),
//
FFmmNearAndFarFields
=
(
FFmmNearField
|
FFmmFarField
)
};
...
...
@@ -52,12 +52,15 @@ protected:
int
nbLevelsInTree
;
void
setNbLevelsInTree
(
const
int
inNbLevelsInTree
){
nbLevelsInTree
=
inNbLevelsInTree
;
nbLevelsInTree
=
inNbLevelsInTree
;
lowerWorkingLevel
=
nbLevelsInTree
;
}
void
validateLevels
()
const
{
std
::
cout
<<
"upperWorkingLevel: "
<<
FAbstractAlgorithm
::
upperWorkingLevel
<<
std
::
endl
<<
"lowerWorkingLevel "
<<
FAbstractAlgorithm
::
lowerWorkingLevel
<<
std
::
endl
;
FAssertLF
(
FAbstractAlgorithm
::
upperWorkingLevel
<=
FAbstractAlgorithm
::
lowerWorkingLevel
);
std
::
cout
<<
"End assert 1"
<<
std
::
endl
;
FAssertLF
(
2
<=
FAbstractAlgorithm
::
upperWorkingLevel
);
}
...
...
Tests/noDist/testAdaptiveUnifFMM.cpp
View file @
e4b8d3b4
...
...
@@ -15,7 +15,6 @@
// ===================================================================================
// ==== CMAKE =====
// @FUSE_BLAS
// @FUSE_FFT
// ================
// Keep in private GIT
...
...
@@ -30,9 +29,6 @@
#include
"Utils/FTic.hpp"
#include
"Containers/FOctree.hpp"
//#include "Containers/FVector.hpp"
//#include "Components/FSimpleLeaf.hpp"
#include
"Utils/FPoint.hpp"
...
...
@@ -102,15 +98,12 @@ int main(int argc, char ** argv){
// const unsigned int NbThreads = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, 1);
//
// accuracy
const
unsigned
int
P
=
5
;
const
unsigned
int
P
=
7
;
const
int
sminM
=
FParameters
::
getValue
(
argc
,
argv
,
LocalOptionMinMultipoleThreshod
.
options
,
P
*
P
*
P
);
const
int
sminL
=
FParameters
::
getValue
(
argc
,
argv
,
LocalOptionMinLocalThreshod
.
options
,
P
*
P
*
P
);
// typedef FTestCell CellClass;
// typedef FAdaptiveTestKernel< CellClass, ContainerClass > KernelClass;
//
typedef
FUnifCell
<
P
>
CellClass
;
typedef
FP2PParticleContainerIndexed
<>
ContainerClass
;
typedef
FSimpleIndexedLeaf
<
ContainerClass
>
LeafClass
;
...
...
@@ -191,8 +184,9 @@ int main(int argc, char ** argv){
std
::
cout
<<
"Working on particles ..."
<<
std
::
endl
;
counter
.
tic
();
const
MatrixKernelClass
MatrixKernel
;
KernelWrapperClass
kernels
(
TreeHeight
,
loader
.
getBoxWidth
(),
loader
.
getCenterOfBox
(),
&
MatrixKernel
,
sminM
,
sminL
);
// FTestKernels FBasicKernels
FmmClass
algo
(
&
tree
,
&
kernels
);
//FFmmAlgorithm FFmmAlgorithmThread
KernelWrapperClass
*
kernels
=
new
KernelWrapperClass
(
TreeHeight
,
loader
.
getBoxWidth
(),
loader
.
getCenterOfBox
(),
&
MatrixKernel
,
sminM
,
sminL
);
// FTestKernels FBasicKernels
FmmClass
algo
(
&
tree
,
kernels
);
//FFmmAlgorithm FFmmAlgorithmThread
// For debug purpose
// Set Global id
...
...
@@ -200,7 +194,7 @@ int main(int argc, char ** argv){
long
int
idCell
=
setGlobalID
(
tree
);
//
algo
.
execute
();
//
counter
.
tac
();
std
::
cout
<<
"Done "
<<
"(@Algorithm = "
<<
counter
.
elapsed
()
<<
" s)."
<<
std
::
endl
;
//
...
...
UTests/utest
Chebyshev
MultiRhs.cpp
→
UTests/utest
Interpolation
MultiRhs.cpp
View file @
e4b8d3b4
...
...
@@ -17,6 +17,7 @@
// ==== CMAKE =====
// @FUSE_BLAS
// ==============
#include
<array>
#include
"ScalFmmConfig.h"
#include
"Utils/FGlobal.hpp"
...
...
@@ -35,29 +36,33 @@
#include
"Kernels/Chebyshev/FChebCell.hpp"
#include
"Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include
"Kernels/Chebyshev/FChebKernel.hpp"
#include
"Kernels/Chebyshev/FChebSymKernel.hpp"
#include
"Kernels/Uniform/FUnifCell.hpp"
#include
"Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include
"Kernels/Uniform/FUnifKernel.hpp"
#include
"Kernels/P2P/FP2PParticleContainerIndexed.hpp"
/*
In this test we compare the spherical FMM results and the direct results.
*/
*/
/** the test class
*
*/
class
Test
ChebyshevDirect
:
public
FUTester
<
Test
ChebyshevDirect
>
{
class
Test
InterpolationKernel
:
public
FUTester
<
Test
InterpolationKernel
>
{
///////////////////////////////////////////////////////////
// The tests!
///////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////
// The tests!
///////////////////////////////////////////////////////////
template
<
class
CellClass
,
class
ContainerClass
,
class
KernelClass
,
class
MatrixKernelClass
,
class
LeafClass
,
class
OctreeClass
,
class
FmmClass
,
const
int
NVals
>
void
RunTest
()
{
// Warning in make test the exec dir it Build/UTests
// Load particles
template
<
class
CellClass
,
class
ContainerClass
,
class
KernelClass
,
class
MatrixKernelClass
,
class
LeafClass
,
class
OctreeClass
,
class
FmmClass
,
const
int
NVals
>
void
RunTest
()
{
// Warning in make test the exec dir it Build/UTests
// Load particles
//
// Load particles
//
...
...
@@ -72,101 +77,113 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
std
::
string
filename
(
SCALFMMDataPath
+
parFile
);
//
FFmaGenericLoader
loader
(
filename
);
if
(
!
loader
.
isOpen
()){
Print
(
"Cannot open particles file."
);
uassert
(
false
);
return
;
}
Print
(
"Number of particles:"
);
Print
(
loader
.
getNumberOfParticles
());
const
int
NbLevels
=
4
;
const
int
SizeSubLevels
=
2
;
// Create Matrix Kernel
const
MatrixKernelClass
MatrixKernel
;
// FUKernelTester is only designed to work with 1/R, i.e. matrix kernel ctor takes no argument.
//
if
(
!
loader
.
isOpen
()){
Print
(
"Cannot open particles file."
);
uassert
(
false
);
return
;
}
Print
(
"Number of particles:"
);
Print
(
loader
.
getNumberOfParticles
());
const
int
NbLevels
=
4
;
const
int
SizeSubLevels
=
2
;
// std::cout << "\nInterpolation FMM (ORDER="<< ORDER << ") ... " << std::endl;
// Create Matrix Kernel
const
MatrixKernelClass
MatrixKernel
;
// FUKernelTester is only designed to work with 1/R, i.e. matrix kernel ctor takes no argument.
//
FSize
nbParticles
=
loader
.
getNumberOfParticles
()
;
FmaRWParticle
<
8
,
8
>*
const
particles
=
new
FmaRWParticle
<
8
,
8
>
[
nbParticles
];
loader
.
fillParticle
(
particles
,
nbParticles
);
//
//
// Create octree
OctreeClass
tree
(
NbLevels
,
SizeSubLevels
,
loader
.
getBoxWidth
(),
loader
.
getCenterOfBox
());
// Insert particle in the tree
//
for
(
int
idxPart
=
0
;
idxPart
<
loader
.
getNumberOfParticles
()
;
++
idxPart
){
tree
.
insert
(
particles
[
idxPart
].
getPosition
()
,
idxPart
,
particles
[
idxPart
].
getPhysicalValue
()
,
0.0
,
0.0
,
0.0
);
// For each particles we associate Nvals charge ( q,0,0,0)
//
/* for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
double q = particles[idxPart].getPhysicalValue();
tree.insert(particles[idxPart].getPosition() , idxPart, q,q);//,0.0,0.0,0.0);
}*/
for
(
int
idxPart
=
0
;
idxPart
<
nbParticles
;
++
idxPart
){
// // Convert FReal[NVALS] to std::array<FReal,NVALS>
std
::
array
<
FReal
,
(
1
+
4
*
1
)
*
NVals
>
physicalState
;
for
(
int
idxVals
=
0
;
idxVals
<
NVals
;
++
idxVals
){
double
q
=
particles
[
idxPart
].
getPhysicalValue
();
physicalState
[
0
*
NVals
+
idxVals
]
=
particles
[
idxPart
].
getPhysicalValue
();
physicalState
[
1
*
NVals
+
idxVals
]
=
0.0
;
physicalState
[
2
*
NVals
+
idxVals
]
=
0.0
;
physicalState
[
3
*
NVals
+
idxVals
]
=
0.0
;
physicalState
[
4
*
NVals
+
idxVals
]
=
0.0
;
}
// put in tree
tree
.
insert
(
particles
[
idxPart
].
getPosition
(),
idxPart
,
physicalState
);
}
// //
// // Create octree
// for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
// FPoint position;
// FReal physicalValue;
// loader.fillParticle(&position,&physicalValue);
// // put in tree
// tree.insert(position, idxPart, physicalValue);
// // get copy
// particles[idxPart].position = position;
// particles[idxPart].physicalValue = physicalValue;
// particles[idxPart].potential = 0.0;
// particles[idxPart].forces[0] = 0.0;
// particles[idxPart].forces[1] = 0.0;
// particles[idxPart].forces[2] = 0.0;
// }
// Run FMM
Print
(
"Fmm..."
);
KernelClass
kernels
(
NbLevels
,
loader
.
getBoxWidth
(),
loader
.
getCenterOfBox
(),
&
MatrixKernel
);
FmmClass
algo
(
&
tree
,
&
kernels
);
algo
.
execute
();
//
FReal
energy
=
0.0
,
energyD
=
0.0
;
// Run direct computation
Print
(
"Direct..."
);
for
(
int
idx
=
0
;
idx
<
loader
.
getNumberOfParticles
()
;
++
idx
){
*
(
particles
[
idx
].
setPotential
())
*=
NVals
;
particles
[
idx
].
getForces
()[
0
]
*=
NVals
;
particles
[
idx
].
getForces
()[
1
]
*=
NVals
;
particles
[
idx
].
getForces
()[
2
]
*=
NVals
;
energyD
+=
particles
[
idx
].
getPotential
()
*
particles
[
idx
].
getPhysicalValue
()
;
}
//
// Compare
Print
(
"Compute Diff..."
);
FMath
::
FAccurater
potentialDiff
;
FMath
::
FAccurater
fx
,
fy
,
fz
;
{
// Check that each particle has been summed with all other
//
// Run FMM
Print
(
"Fmm..."
);
KernelClass
kernels
(
NbLevels
,
loader
.
getBoxWidth
(),
loader
.
getCenterOfBox
(),
&
MatrixKernel
);
FmmClass
algo
(
&
tree
,
&
kernels
);
algo
.
execute
();
//
FReal
energy
=
0.0
,
energyD
=
0.0
;
for
(
int
idx
=
0
;
idx
<
loader
.
getNumberOfParticles
()
;
++
idx
){
energyD
+=
particles
[
idx
].
getPotential
()
*
particles
[
idx
].
getPhysicalValue
()
;
}
//
// Compare
Print
(
"Compute Diff..."
);
FMath
::
FAccurater
potentialDiff
[
NVals
];
FMath
::
FAccurater
fx
,
fy
,
fz
;
{
tree
.
forEachLeaf
([
&
](
LeafClass
*
leaf
){
const
FReal
*
const
potentials
=
leaf
->
getTargets
()
->
getPotentials
();
const
FReal
*
const
physicalValues
=
leaf
->
getTargets
()
->
getPhysicalValues
();
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
FVector
<
int
>&
indexes
=
leaf
->
getTargets
()
->
getIndexes
();
for
(
int
idxPart
=
0
;
idxPart
<
nbParticlesInLeaf
;
++
idxPart
){
const
int
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
]);
energy
+=
potentials
[
idxPart
]
*
physicalValues
[
idxPart
];
//
for
(
int
idxVals
=
0
;
idxVals
<
NVals
;
++
idxVals
){
const
FReal
*
const
physicalValues
=
leaf
->
getTargets
()
->
getPhysicalValues
(
idxVals
);
const
FReal
*
const
potentials
=
leaf
->
getTargets
()
->
getPotentials
(
idxVals
);
const
FReal
*
const
forcesX
=
leaf
->
getTargets
()
->
getForcesX
(
idxVals
);
const
FReal
*
const
forcesY
=
leaf
->
getTargets
()
->
getForcesY
(
idxVals
);
const
FReal
*
const
forcesZ
=
leaf
->
getTargets
()
->
getForcesZ
(
idxVals
);
const
int
nbParticlesInLeaf
=
leaf
->
getTargets
()
->
getNbParticles
();
const
FVector
<
int
>&
indexes
=
leaf
->
getTargets
()
->
getIndexes
();
for
(
int
idxPart
=
0
;
idxPart
<
nbParticlesInLeaf
;
++
idxPart
){
const
int
indexPartOrig
=
indexes
[
idxPart
];
// std::cout << " index "<< indexPartOrig << " " << particles[indexPartOrig].getPotential() << " " << potentials[idxPart] << std::endl;
potentialDiff
[
idxVals
].
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
];
}
}
});
}
}
delete
[]
particles
;
// Print for information
energy
/=
NVals
;
// Print for information
double
errorPotRL2
=
0.0
,
errorPotRMS
=
0.0
;
Print
(
"Potential diff is = "
);
printf
(
" Pot L2Norm %e
\n
"
,
potentialDiff
.
getL2Norm
());
printf
(
" Pot RL2Norm %e
\n
"
,
potentialDiff
.
getRelativeL2Norm
());
printf
(
" Pot RMSError %e
\n
"
,
potentialDiff
.
getRMSError
());
for
(
int
idxVals
=
0
;
idxVals
<
NVals
;
++
idxVals
){
printf
(
" Charge: %d
\n
"
,
idxVals
);
printf
(
" Pot L2Norm %e
\n
"
,
potentialDiff
[
idxVals
].
getL2Norm
());
printf
(
" Pot RL2Norm %e
\n
"
,
potentialDiff
[
idxVals
].
getRelativeL2Norm
());
printf
(
" Pot RMSError %e
\n
"
,
potentialDiff
[
idxVals
].
getRMSError
());
errorPotRL2
=
std
::
max
(
errorPotRL2
,
potentialDiff
[
idxVals
].
getRelativeL2Norm
());
errorPotRMS
=
std
::
max
(
errorPotRMS
,
potentialDiff
[
idxVals
].
getRMSError
());
}
Print
(
"Fx diff is = "
);
printf
(
" Fx L2Norm %e
\n
"
,
fx
.
getL2Norm
());
printf
(
" Fx RL2Norm %e
\n
"
,
fx
.
getRelativeL2Norm
());
...
...
@@ -185,14 +202,14 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
printf
(
" Energy FMM = %.12e
\n
"
,
FMath
::
Abs
(
energy
));
printf
(
" Energy DIRECT = %.12e
\n
"
,
FMath
::
Abs
(
energyD
));
// Assert
const
FReal
MaximumDiffPotential
=
FReal
(
9e-3
);
const
FReal
MaximumDiffForces
=
FReal
(
9e-2
);
// Assert
const
FReal
MaximumDiffPotential
=
FReal
(
9e-3
);
const
FReal
MaximumDiffForces
=
FReal
(
9e-2
);
Print
(
"Test1 - Error Relative L2 norm Potential "
);
uassert
(
potentialDiff
.
getRelativeL2Norm
()
<
MaximumDiffPotential
);
//1
uassert
(
errorPotRL2
<
MaximumDiffPotential
);
//1
Print
(
"Test2 - Error RMS L2 norm Potential "
);
uassert
(
potentialDiff
.
getRMSError
()
<
MaximumDiffPotential
);
//2
uassert
(
errorPotRMS
<
MaximumDiffPotential
);
//2
Print
(
"Test3 - Error Relative L2 norm FX "
);
uassert
(
fx
.
getRelativeL2Norm
()
<
MaximumDiffForces
);
//3
Print
(
"Test4 - Error RMS L2 norm FX "
);
...
...
@@ -211,92 +228,96 @@ class TestChebyshevDirect : public FUTester<TestChebyshevDirect> {
uassert
(
FMath
::
Abs
(
energy
-
energyD
)
/
energyD
<
MaximumDiffPotential
);
//10 Total Energy
// Compute multipole local rhs diff
FMath
::
FAccurater
localDiff
;
FMath
::
FAccurater
multiPoleDiff
;
tree
.
forEachCell
([
&
](
CellClass
*
cell
){
for
(
int
idxRhs
=
1
;
idxRhs
<
NVals
;
++
idxRhs
){
localDiff
.
add
(
cell
->
getLocal
(
0
),
cell
->
getLocal
(
idxRhs
),
cell
->
getVectorSize
());
multiPoleDiff
.
add
(
cell
->
getMultipole
(
0
),
cell
->
getMultipole
(
idxRhs
),
cell
->
getVectorSize
());
}
});
Print
(
"Local diff is = "
);
Print
(
localDiff
.
getL2Norm
());
Print
(
localDiff
.
getInfNorm
());
Print
(
"Multipole diff is = "
);
Print
(
multiPoleDiff
.
getL2Norm
());
Print
(
multiPoleDiff
.
getInfNorm
());
uassert
(
localDiff
.
getL2Norm
()
<
1e-10
);
uassert
(
localDiff
.
getInfNorm
()
<
1e-10
);
uassert
(
multiPoleDiff
.
getL2Norm
()
<
1e-10
);
uassert
(
multiPoleDiff
.
getInfNorm
()
<
1e-10
);
}
/** If memstas is running print the memory used */
void
PostTest
()
{
if
(
FMemStats
::
controler
.
isUsed
()
){
std
::
cout
<<
"Memory used at the end "
<<
FMemStats
::
controler
.
getCurrentAllocated
()
<<
" Bytes ("
<<
FMemStats
::
controler
.
getCurrentAllocatedMB
()
<<
"MB)
\n
"
;
std
::
cout
<<
"Max memory used "
<<
FMemStats
::
controler
.
getMaxAllocated
()
<<
" Bytes ("
<<
FMemStats
::
controler
.
getMaxAllocatedMB
()
<<
"MB)
\n
"
;
std
::
cout
<<
"Total memory used "
<<
FMemStats
::
controler
.
getTotalAllocated
()
<<
" Bytes ("
<<
FMemStats
::
controler
.
getTotalAllocatedMB
()
<<
"MB)
\n
"
;
}
}
///////////////////////////////////////////////////////////
// Set the tests!
///////////////////////////////////////////////////////////
/** TestChebKernel */
void
TestChebKernel
(){
const
int
NVals
=
4
;
const
unsigned
int
ORDER
=
6
;
typedef
FP2PParticleContainerIndexed
<>
ContainerClass
;
typedef
FSimpleLeaf
<
ContainerClass
>
LeafClass
;
typedef
FInterpMatrixKernelR
MatrixKernelClass
;
typedef
FChebCell
<
ORDER
,
1
,
1
,
NVals
>
CellClass
;
typedef
FOctree
<
CellClass
,
ContainerClass
,
LeafClass
>
OctreeClass
;
typedef
FChebKernel
<
CellClass
,
ContainerClass
,
MatrixKernelClass
,
ORDER
,
NVals
>
KernelClass
;
typedef
FFmmAlgorithm
<
OctreeClass
,
CellClass
,
ContainerClass
,
KernelClass
,
LeafClass
>
FmmClass
;
// run test
RunTest
<
CellClass
,
ContainerClass
,
KernelClass
,
MatrixKernelClass
,
LeafClass
,
OctreeClass
,
FmmClass
,
NVals
>
();
}
/** TestChebSymKernel */
void
TestChebSymKernel
(){
const
int
NVals
=
4
;
const
unsigned
int
ORDER
=
6
;
typedef
FP2PParticleContainerIndexed
<>
ContainerClass
;
typedef
FSimpleLeaf
<
ContainerClass
>
LeafClass
;
typedef
FInterpMatrixKernelR
MatrixKernelClass
;
typedef
FChebCell
<
ORDER
,
1
,
1
,
NVals
>
CellClass
;
typedef
FOctree
<
CellClass
,
ContainerClass
,
LeafClass
>
OctreeClass
;
typedef
FChebSymKernel
<
CellClass
,
ContainerClass
,
MatrixKernelClass
,
ORDER
,
NVals
>
KernelClass
;
typedef
FFmmAlgorithm
<
OctreeClass
,
CellClass
,
ContainerClass
,
KernelClass
,
LeafClass
>
FmmClass
;
// run test
RunTest
<
CellClass
,
ContainerClass
,
KernelClass
,
MatrixKernelClass
,
LeafClass
,
OctreeClass
,
FmmClass
,
NVals
>
();
}
///////////////////////////////////////////////////////////
// Set the tests!
///////////////////////////////////////////////////////////
/** set test */
void
SetTests
(){
AddTest
(
&
TestChebyshevDirect
::
TestChebKernel
,
"Test Chebyshev Kernel with one big SVD"
);
AddTest
(
&
TestChebyshevDirect
::
TestChebSymKernel
,
"Test Chebyshev Kernel with 16 small SVDs and symmetries"
);
}
// Compute multipole local rhs diff
FMath
::
FAccurater
localDiff
;
FMath
::
FAccurater
multiPoleDiff
;
tree
.
forEachCell
([
&
](
CellClass
*
cell
){
for
(
int
idxRhs
=
1
;
idxRhs
<
NVals
;
++
idxRhs
){
localDiff
.
add
(
cell
->
getLocal
(
0
),
cell
->
getLocal
(
idxRhs
),
cell
->
getVectorSize
());
multiPoleDiff
.
add
(
cell
->
getMultipole
(
0
),
cell
->
getMultipole
(
idxRhs
),
cell
->
getVectorSize
());
}
});
Print
(
"Local diff is = "
);
Print
(
localDiff
.
getL2Norm
());
Print
(
localDiff
.
getInfNorm
());
Print
(
"Multipole diff is = "
);
Print
(
multiPoleDiff
.
getL2Norm
());
Print
(
multiPoleDiff
.
getInfNorm
());
uassert
(
localDiff
.
getL2Norm
()
<
1e-10
);
uassert
(
localDiff
.
getInfNorm
()
<
1e-10
);
uassert
(
multiPoleDiff
.
getL2Norm
()
<
1e-10
);
uassert
(
multiPoleDiff
.
getInfNorm
()
<
1e-10
);
}
/** If memstas is running print the memory used */
void
PostTest
()
{
if
(
FMemStats
::
controler
.
isUsed
()
){
std
::
cout
<<
"Memory used at the end "
<<
FMemStats
::
controler
.
getCurrentAllocated
()
<<
" Bytes ("
<<
FMemStats
::
controler
.
getCurrentAllocatedMB
()
<<
"MB)
\n
"
;
std
::
cout
<<
"Max memory used "
<<
FMemStats
::
controler
.
getMaxAllocated
()
<<
" Bytes ("
<<
FMemStats
::
controler
.
getMaxAllocatedMB
()
<<
"MB)
\n
"
;
std
::
cout
<<
"Total memory used "
<<
FMemStats
::
controler
.
getTotalAllocated
()
<<
" Bytes ("
<<
FMemStats
::
controler
.
getTotalAllocatedMB
()
<<
"MB)
\n
"
;
}
}
///////////////////////////////////////////////////////////
// Set the tests!
///////////////////////////////////////////////////////////
/** TestUnifKernel */
void
TestUnifKernel
(){
const
int
NVals
=
1
;
const
unsigned
int
ORDER
=
6
;
// run test
typedef
FInterpMatrixKernelR
MatrixKernelClass
;
typedef
FP2PParticleContainerIndexed
<
1
,
1
,
NVals
>
ContainerClass
;
typedef
FSimpleLeaf
<
ContainerClass
>
LeafClass
;
typedef
FUnifCell
<
ORDER
,
1
,
1
,
NVals
>
CellClass
;
typedef
FOctree
<
CellClass
,
ContainerClass
,
LeafClass
>
OctreeClass
;
typedef
FUnifKernel
<
CellClass
,
ContainerClass
,
MatrixKernelClass
,
ORDER
,
NVals
>
KernelClass
;
typedef
FFmmAlgorithm
<
OctreeClass
,
CellClass
,
ContainerClass
,
KernelClass
,
LeafClass
>
FmmClass
;
RunTest
<
CellClass
,
ContainerClass
,
KernelClass
,
MatrixKernelClass
,
LeafClass
,
OctreeClass
,
FmmClass
,
NVals
>
();
}
/** TestChebSymKernel */
void
TestChebSymKernel
(){
const
int
NVals
=
2
;
const
unsigned
int
ORDER
=
6
;
typedef
FP2PParticleContainerIndexed
<
1
,
1
,
NVals
>
ContainerClass
;
typedef
FSimpleLeaf
<
ContainerClass
>
LeafClass
;
typedef
FInterpMatrixKernelR
MatrixKernelClass
;
typedef
FChebCell
<
ORDER
,
1
,
1
,
NVals
>
CellClass
;
typedef
FOctree
<
CellClass
,
ContainerClass
,
LeafClass
>
OctreeClass
;
typedef
FChebSymKernel
<
CellClass
,
ContainerClass
,
MatrixKernelClass
,
ORDER
,
NVals
>
KernelClass
;
typedef
FFmmAlgorithm
<
OctreeClass
,
CellClass
,
ContainerClass
,
KernelClass
,
LeafClass
>
FmmClass
;
// run test
RunTest
<
CellClass
,
ContainerClass
,
KernelClass
,
MatrixKernelClass
,
LeafClass
,
OctreeClass
,
FmmClass
,
NVals
>
();
}
///////////////////////////////////////////////////////////
// Set the tests!
///////////////////////////////////////////////////////////
/** set test */
void
SetTests
(){
AddTest
(
&
TestInterpolationKernel
::
TestUnifKernel
,
"Test Lagrange/Uniform grid FMM"
);
AddTest
(
&
TestInterpolationKernel
::
TestChebSymKernel
,
"Test Symmetric Chebyshev Kernel with 16 small SVDs and symmetries"
);
}
};
// You must do this
TestClass
(
Test
ChebyshevDirect
)
TestClass
(
Test
InterpolationKernel
)
...
...
UTests/utestSphericalDirect.cpp
View file @
e4b8d3b4
...
...
@@ -166,8 +166,10 @@ class TestSphericalDirect : public FUTester<TestSphericalDirect> {
printf
(
" Energy DIRECT = %.12e
\n
"
,
FMath
::
Abs
(
energyD
));
// Assert
const
FReal
MaximumDiffPotential
=
FReal
(
9e-3
);
const
FReal
MaximumDiffForces
=
FReal
(
9e-2
);
double
epsilon
=
1.0
/
FMath
::
pow2
(
ORDER
);
const
FReal
MaximumDiffPotential
=
FReal
(
epsilon
);
const
FReal
MaximumDiffForces
=
FReal
(
10
*
epsilon
);
printf
(
" Criteria error - Epsilon %e
\n
"
,
epsilon
);
Print
(
"Test1 - Error Relative L2 norm Potential "
);
uassert
(
potentialDiff
.
getRelativeL2Norm
()
<
MaximumDiffPotential
);
//1
...
...
@@ -328,12 +330,12 @@ class TestSphericalDirect : public FUTester<TestSphericalDirect> {
OctreeClass
,
FmmClass
,
24
>
(
true
);
RunTest
<
CellClass
,
ContainerClass
,
KernelClass
,
LeafClass
,
OctreeClass
,
FmmClass
,
26
>
(
true
);
RunTest
<
CellClass
,
ContainerClass
,
KernelClass
,
LeafClass
,
OctreeClass
,
FmmClass
,
28
>
(
true
);
RunTest
<
CellClass
,
ContainerClass
,
KernelClass
,
LeafClass
,
OctreeClass
,
FmmClass
,
30
>
(
true
);
RunTest
<
CellClass
,
ContainerClass
,
KernelClass
,
LeafClass
,
OctreeClass
,
FmmClass
,
32
>
(
true
);
//
RunTest< CellClass, ContainerClass, KernelClass, LeafClass,
//
OctreeClass, FmmClass, 28>(true);
//
RunTest< CellClass, ContainerClass, KernelClass, LeafClass,
//
OctreeClass, FmmClass, 30>(true);
//
RunTest< CellClass, ContainerClass, KernelClass, LeafClass,