Skip to content
GitLab
Commit 8cc1b074 authored by Matthias Messner's avatar Matthias Messner
Browse files

only one direct interaction computation in testCompareKernels; changed vector 

size to 512 in FSphericalBlockBlasKernel
parent cb425050
Hide whitespace changes
Inline Side-by-side
Src/Kernels/Spherical/FSphericalBlockBlasKernel.hpp
View file @ 8cc1b074
......@@ -128,7 +128,7 @@ public:
* @param inBoxWidth the size of the simulation box
* @param inPeriodicLevel the number of level upper to 0 that will be requiried
*/
FSphericalBlockBlasKernel(const int inDevP, const int inTreeHeight, const FReal inBoxWidth, const F3DPosition& inBoxCenter, const int inBlockSize = 90, const int inPeriodicLevel = 0)
FSphericalBlockBlasKernel(const int inDevP, const int inTreeHeight, const FReal inBoxWidth, const F3DPosition& inBoxCenter, const int inBlockSize = 512, const int inPeriodicLevel = 0)
: Parent(inDevP, inTreeHeight, inBoxWidth, inBoxCenter, inPeriodicLevel),
FF_MATRIX_ROW_DIM(Parent::harmonic.getExpSize()), FF_MATRIX_COLUMN_DIM(Parent::harmonic.getNExpSize()),
FF_MATRIX_SIZE(FF_MATRIX_ROW_DIM * FF_MATRIX_COLUMN_DIM),
......
Tests/testCompareKernels.cpp
View file @ 8cc1b074
......@@ -40,6 +40,7 @@
#include "../Src/Kernels/Chebyshev/FChebCell.hpp"
#include "../Src/Kernels/Chebyshev/FChebMatrixKernel.hpp"
#include "../Src/Kernels/Chebyshev/FChebKernel.hpp"
#include "../Src/Kernels/Chebyshev/FChebSymKernel.hpp"
// spherical kernel
#include "../Src/Components/FSimpleLeaf.hpp"
......@@ -84,16 +85,98 @@ FReal computeINFnorm(const unsigned int N, const FReal *const u, const FReal *co
}
template <class OctreeClass, class MatrixKernelClass, class ContainerClass>
class DirectInteractionComputer
{
const MatrixKernelClass MatrixKernel;
public:
DirectInteractionComputer()
: MatrixKernel()
{}
unsigned int operator()(OctreeClass& tree, const unsigned int NumTargetCells, FReal* &p, FReal* &f) const
{
// begin direct computation of first leaf cell
typename OctreeClass::Iterator iLeafs(&tree);
iLeafs.gotoBottomLeft();
const ContainerClass* *const TargetSets = new const ContainerClass* [NumTargetCells];
unsigned int n = 0;
for (unsigned int t=0; t<NumTargetCells; ++t) {
TargetSets[t] = iLeafs.getCurrentListTargets();
n += TargetSets[t]->getSize();
iLeafs.moveRight();
}
p = new FReal [n];
FBlas::setzero(n, p);
f = new FReal [n * 3];
FBlas::setzero(n * 3, f);
std::cout << "\nDirect computation of " << n << " target particles ..." << std::endl;
unsigned int start = 0;
for (unsigned int t=0; t<NumTargetCells; ++t) {
iLeafs.gotoBottomLeft();
// retrieve targets
const ContainerClass *const Targets = TargetSets[t];
do {
const ContainerClass *const Sources = iLeafs.getCurrentListSrc();
unsigned int counter = start;
typename ContainerClass::ConstBasicIterator iTarget(*Targets);
while(iTarget.hasNotFinished()) {
const FReal wt = iTarget.data().getPhysicalValue();
typename ContainerClass::ConstBasicIterator iSource(*Sources);
while(iSource.hasNotFinished()) {
if (&iTarget.data() != &iSource.data()) {
const FReal ws = iSource.data().getPhysicalValue();
const FReal one_over_r = MatrixKernel.evaluate(iTarget.data().getPosition(),
iSource.data().getPosition());
// potential
p[counter] += one_over_r * ws;
// force
F3DPosition force(iTarget.data().getPosition() - iSource.data().getPosition());
force *= ((ws*wt) * (one_over_r*one_over_r*one_over_r));
f[counter*3 + 0] += force.getX();
f[counter*3 + 1] += force.getY();
f[counter*3 + 2] += force.getZ();
}
iSource.gotoNext();
}
counter++;
iTarget.gotoNext();
}
} while(iLeafs.moveRight());
start += Targets->getSize();
}
delete [] TargetSets;
return n;
}
};
// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
// get info from commandline
const char* const filename = FParameters::getStr(argc,argv,"-f", "../Data/test20k.fma");
const char* const filename = FParameters::getStr(argc,argv,"-f", "../Data/test20k.fma");
const unsigned int TreeHeight = FParameters::getValue(argc, argv, "-h", 5);
const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, "-sh", 2);
const unsigned int NbThreads = FParameters::getValue(argc, argv, "-t", 1);
omp_set_num_threads(NbThreads);
std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl;
// init timer
FTic time;
......@@ -107,19 +190,22 @@ int main(int argc, char* argv[])
// // number of particles
// unsigned int NumParticles = 0;
unsigned int nt1 = 0;
FReal* p10; p10 = NULL;
FReal* f10; p10 = NULL;
////////////////////////////////////////////////////////////////////
{ // begin Chebyshef kernel
// accuracy
const unsigned int ORDER = 4;
const FReal epsilon = FReal(1e-5);
const unsigned int ORDER = 7;
const FReal epsilon = FReal(1e-7);
unsigned int nt1 = 0;
//unsigned int nt1 = 0;
FReal* p1; p1 = NULL;
FReal* p10; p10 = NULL;
//FReal* p10; p10 = NULL;
FReal* f1; p1 = NULL;
FReal* f10; p10 = NULL;
//FReal* f10; p10 = NULL;
// typedefs
typedef FChebParticle ParticleClass;
......@@ -128,9 +214,12 @@ int main(int argc, char* argv[])
typedef FChebMatrixKernelR MatrixKernelClass;
typedef FChebCell<ORDER> CellClass;
typedef FOctree<ParticleClass,CellClass,ContainerClass,LeafClass> OctreeClass;
typedef FChebKernel<ParticleClass,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
typedef FFmmAlgorithm<OctreeClass,ParticleClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
//typedef FFmmAlgorithmThread<OctreeClass,ParticleClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
//typedef FChebKernel<ParticleClass,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
typedef FChebSymKernel<ParticleClass,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
//typedef FFmmAlgorithm<OctreeClass,ParticleClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
typedef FFmmAlgorithmThread<OctreeClass,ParticleClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
// open particle file
FFmaScanfLoader<ParticleClass> loader(filename);
......@@ -182,52 +271,12 @@ int main(int argc, char* argv[])
} while(iLeafs.moveRight());
} // -----------------------------------------------------
// compute direct interaction
const unsigned int NumTargetCells = 3;
DirectInteractionComputer<OctreeClass,MatrixKernelClass,ContainerClass> direct;
nt1 = direct(tree, NumTargetCells, p10, f10);
{ // -----------------------------------------------------
// begin direct computation of first leaf cell
OctreeClass::Iterator iLeafs(&tree);
iLeafs.gotoBottomLeft();
// retrieve targets
const ContainerClass *const Targets = iLeafs.getCurrentListTargets();
nt1 = Targets->getSize();
p10 = new FReal [nt1];
FBlas::setzero(nt1, p10);
f10 = new FReal [nt1 * 3];
FBlas::setzero(nt1 * 3, f10);
std::cout << "\nDirect computation of " << nt1 << " target particles ..." << std::endl;
const MatrixKernelClass MatrixKernel;
do {
const ContainerClass *const Sources = iLeafs.getCurrentListSrc();
unsigned int counter = 0;
ContainerClass::ConstBasicIterator iTarget(*Targets);
while(iTarget.hasNotFinished()) {
const FReal wt = iTarget.data().getPhysicalValue();
ContainerClass::ConstBasicIterator iSource(*Sources);
while(iSource.hasNotFinished()) {
if (&iTarget.data() != &iSource.data()) {
const FReal ws = iSource.data().getPhysicalValue();
const FReal one_over_r = MatrixKernel.evaluate(iTarget.data().getPosition(),
iSource.data().getPosition());
// potential
p10[counter] += one_over_r * ws;
// force
F3DPosition force(iTarget.data().getPosition() - iSource.data().getPosition());
force *= ((ws*wt) * (one_over_r*one_over_r*one_over_r));
f10[counter*3 + 0] += force.getX();
f10[counter*3 + 1] += force.getY();
f10[counter*3 + 2] += force.getZ();
}
iSource.gotoNext();
}
counter++;
iTarget.gotoNext();
}
} while(iLeafs.moveRight());
} // -----------------------------------------------------
// for (unsigned int n=0; n<nt1; ++n)
// std::cout << p10[n] << " - " << p1[n] << " = " << p10[n]-p1[n] << std::endl;
......@@ -240,9 +289,7 @@ int main(int argc, char* argv[])
std::cout << "Relative Lmax error = " << computeINFnorm(nt1*3, f10, f1) << std::endl << std::endl;
if (p10!=NULL) delete [] p10;
if (p1 !=NULL) delete [] p1;
if (f10!=NULL) delete [] f10;
if (f1 !=NULL) delete [] f1;
} // end Chebyshev kernel
......@@ -252,7 +299,7 @@ int main(int argc, char* argv[])
{ // begin FFmaBlas kernel
// accuracy
const int DevP = FParameters::getValue(argc, argv, "-p", 8);
const int DevP = FParameters::getValue(argc, argv, "-p", 5);
// typedefs
typedef FSphericalParticle ParticleClass;
......@@ -263,15 +310,18 @@ int main(int argc, char* argv[])
//typedef FSphericalBlasKernel<ParticleClass, CellClass, ContainerClass > KernelClass;
//typedef FSphericalBlockBlasKernel<ParticleClass, CellClass, ContainerClass > KernelClass;
//typedef FSphericalKernel<ParticleClass, CellClass, ContainerClass > KernelClass;
typedef FSphericalBlockBlasKernel<ParticleClass, CellClass, ContainerClass > KernelClass;
typedef FFmmAlgorithm<OctreeClass, ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
//typedef FSphericalBlockBlasKernel<ParticleClass, CellClass, ContainerClass > KernelClass;
typedef FSphericalRotationKernel<ParticleClass, CellClass, ContainerClass > KernelClass;
//typedef FFmmAlgorithm<OctreeClass, ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
typedef FFmmAlgorithmThread<OctreeClass, ParticleClass, CellClass, ContainerClass, KernelClass, LeafClass > FmmClass;
// open particle file
FFmaScanfLoader<ParticleClass> loader(filename);
if(!loader.isOpen()) throw std::runtime_error("Particle file couldn't be opened!");
// init cell class and oct-tree
CellClass::Init(DevP, true);
CellClass::Init(DevP, true); // only for blas
CellClass::Init(DevP, false);
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
{ // -----------------------------------------------------
......@@ -318,73 +368,19 @@ int main(int argc, char* argv[])
} while(leavesIterator.moveRight());
} // -----------------------------------------------------
unsigned int sizeFirstLeaf = 0;
FReal* firstLeafPotential = 0;
FReal* firstLeafForces = 0;
{ // -----------------------------------------------------
// begin direct computation of first leaf cell
OctreeClass::Iterator leavesIterator(&tree);
leavesIterator.gotoBottomLeft();
ContainerClass *const firstLeaf = leavesIterator.getCurrentListTargets();
{ // Compute the first leaf with it self
ContainerClass::BasicIterator targetParticles(*firstLeaf);
while(targetParticles.hasNotFinished()) {
ContainerClass::ConstBasicIterator sourceParticles(*firstLeaf);
while(sourceParticles.hasNotFinished()) {
if (&targetParticles.data() != &sourceParticles.data()){
kernels.directInteraction( &targetParticles.data(), sourceParticles.data());
}
sourceParticles.gotoNext();
}
targetParticles.gotoNext();
}
}
while( leavesIterator.moveRight() ) {
ContainerClass::ConstBasicIterator sourceParticles(*leavesIterator.getCurrentListSrc());
while(sourceParticles.hasNotFinished()) {
ContainerClass::BasicIterator targetParticles(*firstLeaf);
while(targetParticles.hasNotFinished()) {
kernels.directInteraction( &targetParticles.data(), sourceParticles.data());
targetParticles.gotoNext();
}
sourceParticles.gotoNext();
}
}
// Copy data
sizeFirstLeaf = firstLeaf->getSize();
firstLeafPotential = new FReal[sizeFirstLeaf];
FBlas::setzero(sizeFirstLeaf, firstLeafPotential);
firstLeafForces = new FReal[sizeFirstLeaf * 3];
FBlas::setzero(sizeFirstLeaf * 3, firstLeafForces);
unsigned int particlePosition = 0;
ContainerClass::ConstBasicIterator targetParticles(*firstLeaf);
while(targetParticles.hasNotFinished()) {
firstLeafPotential[particlePosition] = targetParticles.data().getPotential();
firstLeafForces[particlePosition*3 + 0] = targetParticles.data().getForces().getX();
firstLeafForces[particlePosition*3 + 1] = targetParticles.data().getForces().getY();
firstLeafForces[particlePosition*3 + 2] = targetParticles.data().getForces().getZ();
targetParticles.gotoNext();
++particlePosition;
}
} // -----------------------------------------------------
FBlas::scal(nt1*3, FReal(-1.), f10);
std::cout << "\nPotential error:" << std::endl;
std::cout << "Relative L2 error = " << computeL2norm( sizeFirstLeaf, firstLeafPotential, fmmParticlesPotential) << std::endl;
std::cout << "Relative Lmax error = " << computeINFnorm(sizeFirstLeaf, firstLeafPotential, fmmParticlesPotential) << std::endl;
std::cout << "Relative L2 error = " << computeL2norm( nt1, p10, fmmParticlesPotential) << std::endl;
std::cout << "Relative Lmax error = " << computeINFnorm(nt1, p10, fmmParticlesPotential) << std::endl;
std::cout << "\nForces error:" << std::endl;
std::cout << "Relative L2 error = " << computeL2norm( sizeFirstLeaf*3, firstLeafForces, fmmParticlesForces) << std::endl;
std::cout << "Relative Lmax error = " << computeINFnorm(sizeFirstLeaf*3, firstLeafForces, fmmParticlesForces) << std::endl;
std::cout << "Relative L2 error = " << computeL2norm( nt1*3, f10, fmmParticlesForces) << std::endl;
std::cout << "Relative Lmax error = " << computeINFnorm(nt1*3, f10, fmmParticlesForces) << std::endl;
delete [] firstLeafPotential;
delete [] fmmParticlesPotential;
delete [] firstLeafForces;
delete [] fmmParticlesForces;
} // end FFmaBlas kernel
......@@ -399,11 +395,9 @@ int main(int argc, char* argv[])
//std::cout << "Relative Lmax error = " << computeINFnorm(NumParticles, p1, p2) << "\n" << std::endl;
// // free memory
// if (p10!=NULL) delete [] p10;
// if (p20!=NULL) delete [] p20;
// if (p1!=NULL) delete [] p1;
// if (p2!=NULL) delete [] p2;
// free memory
if (p10!=NULL) delete [] p10;
if (f10!=NULL) delete [] f10;
return 0;
}
......@@ -437,3 +431,73 @@ particleIterator.gotoNext();
}
} while(octreeIterator.moveRight());
*/
//unsigned int sizeFirstLeaf = 0;
//FReal* firstLeafPotential = 0;
//FReal* firstLeafForces = 0;
//{ // -----------------------------------------------------
// // begin direct computation of first leaf cell
// OctreeClass::Iterator leavesIterator(&tree);
// leavesIterator.gotoBottomLeft();
// ContainerClass *const firstLeaf = leavesIterator.getCurrentListTargets();
//
// { // Compute the first leaf with it self
// ContainerClass::BasicIterator targetParticles(*firstLeaf);
// while(targetParticles.hasNotFinished()) {
// ContainerClass::ConstBasicIterator sourceParticles(*firstLeaf);
// while(sourceParticles.hasNotFinished()) {
// if (&targetParticles.data() != &sourceParticles.data()){
// kernels.directInteraction( &targetParticles.data(), sourceParticles.data());
// }
// sourceParticles.gotoNext();
// }
// targetParticles.gotoNext();
// }
// }
//
// while( leavesIterator.moveRight() ) {
// ContainerClass::ConstBasicIterator sourceParticles(*leavesIterator.getCurrentListSrc());
// while(sourceParticles.hasNotFinished()) {
// ContainerClass::BasicIterator targetParticles(*firstLeaf);
// while(targetParticles.hasNotFinished()) {
// kernels.directInteraction( &targetParticles.data(), sourceParticles.data());
// targetParticles.gotoNext();
// }
// sourceParticles.gotoNext();
// }
// }
//
// // Copy data
// sizeFirstLeaf = firstLeaf->getSize();
//
// firstLeafPotential = new FReal[sizeFirstLeaf];
// FBlas::setzero(sizeFirstLeaf, firstLeafPotential);
//
// firstLeafForces = new FReal[sizeFirstLeaf * 3];
// FBlas::setzero(sizeFirstLeaf * 3, firstLeafForces);
//
// unsigned int particlePosition = 0;
// ContainerClass::ConstBasicIterator targetParticles(*firstLeaf);
// while(targetParticles.hasNotFinished()) {
// firstLeafPotential[particlePosition] = targetParticles.data().getPotential();
// firstLeafForces[particlePosition*3 + 0] = targetParticles.data().getForces().getX();
// firstLeafForces[particlePosition*3 + 1] = targetParticles.data().getForces().getY();
// firstLeafForces[particlePosition*3 + 2] = targetParticles.data().getForces().getZ();
// targetParticles.gotoNext();
// ++particlePosition;
// }
//} // -----------------------------------------------------
//std::cout << "\nPotential error:" << std::endl;
//std::cout << "Relative L2 error = " << computeL2norm( sizeFirstLeaf, firstLeafPotential, fmmParticlesPotential) << std::endl;
//std::cout << "Relative Lmax error = " << computeINFnorm(sizeFirstLeaf, firstLeafPotential, fmmParticlesPotential) << std::endl;
//
//std::cout << "\nForces error:" << std::endl;
//std::cout << "Relative L2 error = " << computeL2norm( sizeFirstLeaf*3, firstLeafForces, fmmParticlesForces) << std::endl;
//std::cout << "Relative Lmax error = " << computeINFnorm(sizeFirstLeaf*3, firstLeafForces, fmmParticlesForces) << std::endl;
Supports Markdown
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