Commit 6e0dfd84 authored by BLANCHARD Pierre's avatar BLANCHARD Pierre
Browse files

Fixed multiRhs utest for ChebFMM.

parent 94a96644
...@@ -45,7 +45,7 @@ class FAbstractChebKernel : public FAbstractKernels< CellClass, ContainerClass> ...@@ -45,7 +45,7 @@ class FAbstractChebKernel : public FAbstractKernels< CellClass, ContainerClass>
{ {
protected: protected:
enum {nnodes = TensorTraits<ORDER>::nnodes}; enum {nnodes = TensorTraits<ORDER>::nnodes};
typedef FChebInterpolator<ORDER,MatrixKernelClass> InterpolatorClass; typedef FChebInterpolator<ORDER,MatrixKernelClass,NVALS> InterpolatorClass;
/// Needed for P2M, M2M, L2L and L2P operators /// Needed for P2M, M2M, L2L and L2P operators
const FSmartPointer<InterpolatorClass,FSmartPointerMemory> Interpolator; const FSmartPointer<InterpolatorClass,FSmartPointerMemory> Interpolator;
......
This diff is collapsed.
...@@ -77,19 +77,22 @@ public: ...@@ -77,19 +77,22 @@ public:
} }
FChebKernel(const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter, const MatrixKernelClass *const inMatrixKernel) FChebKernel(const int inTreeHeight, const FReal inBoxWidth, const FPoint& inBoxCenter, const MatrixKernelClass *const inMatrixKernel)
: FChebKernel(inTreeHeight, inBoxWidth,inBoxCenter,inMatrixKernel,FMath::pow(10.0,static_cast<FReal>(-ORDER))) : FChebKernel(inTreeHeight, inBoxWidth,inBoxCenter,inMatrixKernel,FReal(0.)/*FMath::pow(10.0,static_cast<FReal>(-ORDER))*/)
{ {
} }
void P2M(CellClass* const LeafCell, void P2M(CellClass* const LeafCell,
const ContainerClass* const SourceParticles) const ContainerClass* const SourceParticles)
{ {
const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate())); const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
// 1) apply Sy
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getMultipole(0), SourceParticles);
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply Sy
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getMultipole(idxRhs), SourceParticles);
// 2) apply B // 2) apply B
M2LHandler->applyB(LeafCell->getMultipole(idxRhs), LeafCell->getMultipole(idxRhs) + AbstractBaseClass::nnodes); M2LHandler->applyB(LeafCell->getMultipole(idxRhs), LeafCell->getMultipole(idxRhs) + AbstractBaseClass::nnodes);
} }
...@@ -102,7 +105,6 @@ public: ...@@ -102,7 +105,6 @@ public:
{ {
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply Sy // 1) apply Sy
FBlas::scal(AbstractBaseClass::nnodes*2, FReal(0.), ParentCell->getMultipole(idxRhs));
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){ for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
if (ChildCells[ChildIndex]){ if (ChildCells[ChildIndex]){
AbstractBaseClass::Interpolator->applyM2M(ChildIndex, ChildCells[ChildIndex]->getMultipole(idxRhs), AbstractBaseClass::Interpolator->applyM2M(ChildIndex, ChildCells[ChildIndex]->getMultipole(idxRhs),
...@@ -192,22 +194,23 @@ public: ...@@ -192,22 +194,23 @@ public:
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply U // 1) apply U
M2LHandler->applyU(LeafCell->getLocal(idxRhs) + AbstractBaseClass::nnodes, const_cast<CellClass*>(LeafCell)->getLocal(idxRhs)); M2LHandler->applyU(LeafCell->getLocal(idxRhs) + AbstractBaseClass::nnodes, const_cast<CellClass*>(LeafCell)->getLocal(idxRhs));
//// 2.a) apply Sx
//AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter,
// AbstractBaseClass::BoxWidthLeaf,
// LeafCell->getLocal(),
// TargetParticles);
//// 2.b) apply Px (grad Sx)
//AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter,
// AbstractBaseClass::BoxWidthLeaf,
// LeafCell->getLocal(),
// TargetParticles);
// 2.c) apply Sx and Px (grad Sx)
AbstractBaseClass::Interpolator->applyL2PTotal(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getLocal(idxRhs), TargetParticles);
} }
//// 2.a) apply Sx
//AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter,
// AbstractBaseClass::BoxWidthLeaf,
// LeafCell->getLocal(0),
// TargetParticles);
//// 2.b) apply Px (grad Sx)
//AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter,
// AbstractBaseClass::BoxWidthLeaf,
// LeafCell->getLocal(0),
// TargetParticles);
// 2.c) apply Sx and Px (grad Sx)
AbstractBaseClass::Interpolator->applyL2PTotal(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getLocal(0), TargetParticles);
} }
void P2P(const FTreeCoordinate& /* LeafCellCoordinate */, // needed for periodic boundary conditions void P2P(const FTreeCoordinate& /* LeafCellCoordinate */, // needed for periodic boundary conditions
......
...@@ -180,10 +180,8 @@ public: ...@@ -180,10 +180,8 @@ public:
{ {
// apply Sy // apply Sy
const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate())); const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf, LeafCell->getMultipole(0), SourceParticles);
LeafCell->getMultipole(idxRhs), SourceParticles);
}
} }
...@@ -445,7 +443,6 @@ public: ...@@ -445,7 +443,6 @@ public:
ContainerClass* const TargetParticles) ContainerClass* const TargetParticles)
{ {
const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate())); const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// // a) apply Sx // // a) apply Sx
// AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter, // AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter,
// AbstractBaseClass::BoxWidthLeaf, // AbstractBaseClass::BoxWidthLeaf,
...@@ -457,10 +454,10 @@ public: ...@@ -457,10 +454,10 @@ public:
// LeafCell->getLocal(), // LeafCell->getLocal(),
// TargetParticles); // TargetParticles);
// c) apply Sx and Px (grad Sx) // c) apply Sx and Px (grad Sx)
AbstractBaseClass::Interpolator->applyL2PTotal(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf, AbstractBaseClass::Interpolator->applyL2PTotal(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getLocal(idxRhs), TargetParticles); LeafCell->getLocal(0), TargetParticles);
}
} }
void P2P(const FTreeCoordinate& /* LeafCellCoordinate */, // needed for periodic boundary conditions void P2P(const FTreeCoordinate& /* LeafCellCoordinate */, // needed for periodic boundary conditions
......
...@@ -38,6 +38,7 @@ ...@@ -38,6 +38,7 @@
#include "Kernels/Chebyshev/FChebCell.hpp" #include "Kernels/Chebyshev/FChebCell.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp" #include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Chebyshev/FChebSymKernel.hpp" #include "Kernels/Chebyshev/FChebSymKernel.hpp"
#include "Kernels/Chebyshev/FChebKernel.hpp"
#include "Kernels/Uniform/FUnifCell.hpp" #include "Kernels/Uniform/FUnifCell.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp" #include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
...@@ -102,35 +103,20 @@ class TestInterpolationKernel : public FUTester<TestInterpolationKernel> { ...@@ -102,35 +103,20 @@ class TestInterpolationKernel : public FUTester<TestInterpolationKernel> {
// Create octree // Create octree
OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox()); OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
// Insert particle in the tree // Insert particle in the tree
// For each particles we associate Nvals charge ( q,0,0,0) for(int idxPart = 0 ; idxPart < nbParticles ; ++idxPart){
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){ // Convert FReal[NVALS] to std::array<FReal,NVALS>
double q = particles[idxPart].getPhysicalValue(); std::array<FReal, (1+4*1)*NVals> physicalState;
if(NVals == 1){ for(int idxVals = 0 ; idxVals < NVals ; ++idxVals){
tree.insert(particles[idxPart].getPosition() , idxPart, q);//,0.0,0.0,0.0); 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;
} }
else if(NVals == 2){ // put in tree
tree.insert(particles[idxPart].getPosition() , idxPart, q, q);//,0.0,0.0,0.0); tree.insert(particles[idxPart].getPosition(), idxPart, physicalState);
} }
else if(NVals == 3){
tree.insert(particles[idxPart].getPosition() , idxPart, q, q,q);//,0.0,0.0,0.0);
} else{
FAssertLF(0, "NVALS should be <= 3");
}
}
// 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]= q;
// 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);
// }
// Run FMM // Run FMM
Print("Fmm..."); Print("Fmm...");
...@@ -307,7 +293,20 @@ class TestInterpolationKernel : public FUTester<TestInterpolationKernel> { ...@@ -307,7 +293,20 @@ class TestInterpolationKernel : public FUTester<TestInterpolationKernel> {
RunTest<CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass, NVals>(); RunTest<CellClass,ContainerClass,KernelClass,MatrixKernelClass,LeafClass,OctreeClass,FmmClass, NVals>();
} }
/** TestChebKernel */
void TestChebKernel(){
const int NVals = 3;
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 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>();
}
/////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////
// Set the tests! // Set the tests!
...@@ -318,6 +317,8 @@ class TestInterpolationKernel : public FUTester<TestInterpolationKernel> { ...@@ -318,6 +317,8 @@ class TestInterpolationKernel : public FUTester<TestInterpolationKernel> {
AddTest(&TestInterpolationKernel::TestUnifKernel,"Test Lagrange/Uniform grid FMM"); AddTest(&TestInterpolationKernel::TestUnifKernel,"Test Lagrange/Uniform grid FMM");
AddTest(&TestInterpolationKernel::TestChebSymKernel,"Test Symmetric Chebyshev Kernel with 16 small SVDs and symmetries"); AddTest(&TestInterpolationKernel::TestChebSymKernel,"Test Symmetric Chebyshev Kernel with 16 small SVDs and symmetries");
AddTest(&TestInterpolationKernel::TestChebKernel,"Test Chebyshev Kernel with 1 large SVD");
} }
}; };
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment