Commit b754185c authored by BRAMAS Berenger's avatar BRAMAS Berenger

add multi rhs for cheb kernel

parent 9af6d563
......@@ -27,7 +27,7 @@
#include "./FChebInterpolator.hpp"
class FTreeCoordinate;
template <KERNEL_FUNCCTION_IDENTIFIER Identifier> struct DirectInteactionComputer;
template <KERNEL_FUNCCTION_IDENTIFIER Identifier, int NVALS> struct DirectInteactionComputer;
/**
* @author Matthias Messner(matthias.messner@inria.fr)
......@@ -42,7 +42,7 @@ template <KERNEL_FUNCCTION_IDENTIFIER Identifier> struct DirectInteactionCompute
* @tparam MatrixKernelClass Type of matrix kernel function
* @tparam ORDER Chebyshev interpolation order
*/
template < class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER>
template < class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER, int NVALS = 1>
class FAbstractChebKernel : public FAbstractKernels< CellClass, ContainerClass>
{
protected:
......@@ -132,22 +132,26 @@ public:
ContainerClass* const NeighborSourceParticles[27],
const int /* size */)
{
DirectInteactionComputer<MatrixKernelClass::Identifier>::P2P(TargetParticles,NeighborSourceParticles);
DirectInteactionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles);
}
void P2PRemote(const FTreeCoordinate& /*inPosition*/,
ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
ContainerClass* const inNeighbors[27], const int /*inSize*/){
DirectInteactionComputer<MatrixKernelClass::Identifier>::P2PRemote(inTargets,inNeighbors,27);
DirectInteactionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27);
}
};
///////////////////////////////////////////////////////
// P2P Wrappers
///////////////////////////////////////////////////////
/*! Specialization for Laplace potential */
template <>
struct DirectInteactionComputer<ONE_OVER_R>
struct DirectInteactionComputer<ONE_OVER_R, 1>
{
template <typename ContainerClass>
static void P2P( ContainerClass* const FRestrict TargetParticles,
......@@ -166,7 +170,7 @@ struct DirectInteactionComputer<ONE_OVER_R>
/*! Specialization for Leonard-Jones potential */
template <>
struct DirectInteactionComputer<LEONARD_JONES_POTENTIAL>
struct DirectInteactionComputer<LEONARD_JONES_POTENTIAL, 1>
{
template <typename ContainerClass>
static void P2P( ContainerClass* const FRestrict TargetParticles,
......@@ -182,6 +186,54 @@ struct DirectInteactionComputer<LEONARD_JONES_POTENTIAL>
}
};
///////////////////////////////////////////////////////
// In case of multi right hand side
///////////////////////////////////////////////////////
template <int NVALS>
struct DirectInteactionComputer<ONE_OVER_R, NVALS>
{
template <typename ContainerClass>
static void P2P( ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FP2P::FullMutual(TargetParticles,NeighborSourceParticles,14);
}
}
template <typename ContainerClass>
static void P2PRemote( ContainerClass* const FRestrict inTargets,
ContainerClass* const inNeighbors[27],
const int inSize){
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FP2P::FullRemote(inTargets,inNeighbors,inSize);
}
}
};
/*! Specialization for Leonard-Jones potential */
template <int NVALS>
struct DirectInteactionComputer<LEONARD_JONES_POTENTIAL, NVALS>
{
template <typename ContainerClass>
static void P2P( ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FP2P::FullMutualLJ(TargetParticles,NeighborSourceParticles,14);
}
}
template <typename ContainerClass>
static void P2PRemote( ContainerClass* const FRestrict inTargets,
ContainerClass* const inNeighbors[27],
const int inSize){
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FP2P::FullRemoteLJ(inTargets,inNeighbors,inSize);
}
}
};
#endif //FCHEBKERNELS_HPP
......
......@@ -28,31 +28,41 @@
* Please read the license
*
* This class defines a cell used in the Chebyshev based FMM.
* @param NVALS is the number of right hand side.
*/
template <int ORDER>
class FChebCell : public FExtendMortonIndex,
public FExtendCoordinate
template <int ORDER, int NVALS = 1>
class FChebCell : public FExtendMortonIndex, public FExtendCoordinate
{
FReal multipole_exp[TensorTraits<ORDER>::nnodes * 2]; //< Multipole expansion
FReal local_exp[TensorTraits<ORDER>::nnodes * 2]; //< Local expansion
static const int VectorSize = TensorTraits<ORDER>::nnodes * 2;
FReal multipole_exp[NVALS * VectorSize]; //< Multipole expansion
FReal local_exp[NVALS * VectorSize]; //< Local expansion
public:
~FChebCell() {}
/** Get Multipole */
const FReal* getMultipole() const
{ return this->multipole_exp; }
const FReal* getMultipole(const int inRhs) const
{ return this->multipole_exp + inRhs*VectorSize;
}
/** Get Local */
const FReal* getLocal() const
{ return this->local_exp; }
const FReal* getLocal(const int inRhs) const{
return this->local_exp + inRhs*VectorSize;
}
/** Get Multipole */
FReal* getMultipole()
{ return this->multipole_exp; }
FReal* getMultipole(const int inRhs){
return this->multipole_exp + inRhs*VectorSize;
}
/** Get Local */
FReal* getLocal()
{ return this->local_exp; }
FReal* getLocal(const int inRhs){
return this->local_exp + inRhs*VectorSize;
}
/** To get the leading dim of a vec */
int getVectorSize() const{
return VectorSize;
}
};
......
This diff is collapsed.
......@@ -40,15 +40,15 @@ class FTreeCoordinate;
* @tparam MatrixKernelClass Type of matrix kernel function
* @tparam ORDER Chebyshev interpolation order
*/
template < class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER>
template < class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER, int NVALS = 1>
class FChebKernel
: public FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER>
: public FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
{
// private types
typedef FChebM2LHandler<ORDER,MatrixKernelClass> M2LHandlerClass;
// using from
typedef FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER>
typedef FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
AbstractBaseClass;
/// Needed for M2L operator
......@@ -64,7 +64,7 @@ public:
const FPoint& inBoxCenter,
const FReal inBoxWidth,
const FReal Epsilon)
: FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER>(inTreeHeight,
: FAbstractChebKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,
inBoxCenter,
inBoxWidth),
M2LHandler(new M2LHandlerClass(Epsilon))
......@@ -78,32 +78,33 @@ public:
void P2M(CellClass* const LeafCell,
const ContainerClass* const SourceParticles)
{
// 1) apply Sy
const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter,
AbstractBaseClass::BoxWidthLeaf,
LeafCell->getMultipole(),
SourceParticles);
// 2) apply B
M2LHandler->applyB(LeafCell->getMultipole(),
LeafCell->getMultipole() + AbstractBaseClass::nnodes);
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply Sy
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getMultipole(idxRhs), SourceParticles);
// 2) apply B
M2LHandler->applyB(LeafCell->getMultipole(idxRhs), LeafCell->getMultipole(idxRhs) + AbstractBaseClass::nnodes);
}
}
void M2M(CellClass* const FRestrict ParentCell,
const CellClass*const FRestrict *const FRestrict ChildCells,
const int TreeLevel)
const int /*TreeLevel*/)
{
// 1) apply Sy
FBlas::scal(AbstractBaseClass::nnodes*2, FReal(0.), ParentCell->getMultipole());
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex)
if (ChildCells[ChildIndex])
AbstractBaseClass::Interpolator->applyM2M(ChildIndex,
ChildCells[ChildIndex]->getMultipole(),
ParentCell->getMultipole());
// 2) apply B
M2LHandler->applyB(ParentCell->getMultipole(),
ParentCell->getMultipole() + AbstractBaseClass::nnodes);
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply Sy
FBlas::scal(AbstractBaseClass::nnodes*2, FReal(0.), ParentCell->getMultipole(idxRhs));
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
if (ChildCells[ChildIndex]){
AbstractBaseClass::Interpolator->applyM2M(ChildIndex, ChildCells[ChildIndex]->getMultipole(idxRhs),
ParentCell->getMultipole(idxRhs));
}
}
// 2) apply B
M2LHandler->applyB(ParentCell->getMultipole(idxRhs), ParentCell->getMultipole(idxRhs) + AbstractBaseClass::nnodes);
}
}
......@@ -127,16 +128,19 @@ public:
void M2L(CellClass* const FRestrict TargetCell,
const CellClass* SourceCells[343],
const int NumSourceCells,
const int /*NumSourceCells*/,
const int TreeLevel)
{
FReal *const CompressedLocalExpansion = TargetCell->getLocal() + AbstractBaseClass::nnodes;
const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
for (int idx=0; idx<343; ++idx)
if (SourceCells[idx])
M2LHandler->applyC(idx, CellWidth,
SourceCells[idx]->getMultipole() + AbstractBaseClass::nnodes,
CompressedLocalExpansion);
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FReal *const CompressedLocalExpansion = TargetCell->getLocal(idxRhs) + AbstractBaseClass::nnodes;
const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
for (int idx=0; idx<343; ++idx){
if (SourceCells[idx]){
M2LHandler->applyC(idx, CellWidth, SourceCells[idx]->getMultipole(idxRhs) + AbstractBaseClass::nnodes,
CompressedLocalExpansion);
}
}
}
}
// void M2L(CellClass* const FRestrict TargetCell,
......@@ -158,44 +162,45 @@ public:
void L2L(const CellClass* const FRestrict ParentCell,
CellClass* FRestrict *const FRestrict ChildCells,
const int TreeLevel)
const int /*TreeLevel*/)
{
// 1) apply U
M2LHandler->applyU(ParentCell->getLocal() + AbstractBaseClass::nnodes,
const_cast<CellClass*>(ParentCell)->getLocal());
// 2) apply Sx
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex)
if (ChildCells[ChildIndex])
AbstractBaseClass::Interpolator->applyL2L(ChildIndex,
ParentCell->getLocal(),
ChildCells[ChildIndex]->getLocal());
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply U
M2LHandler->applyU(ParentCell->getLocal(idxRhs) + AbstractBaseClass::nnodes,
const_cast<CellClass*>(ParentCell)->getLocal(idxRhs));
// 2) apply Sx
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
if (ChildCells[ChildIndex]){
AbstractBaseClass::Interpolator->applyL2L(ChildIndex, ParentCell->getLocal(idxRhs), ChildCells[ChildIndex]->getLocal(idxRhs));
}
}
}
}
void L2P(const CellClass* const LeafCell,
ContainerClass* const TargetParticles)
{
// 1) apply U
M2LHandler->applyU(LeafCell->getLocal() + AbstractBaseClass::nnodes,
const_cast<CellClass*>(LeafCell)->getLocal());
const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
//// 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(),
TargetParticles);
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply U
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);
}
}
};
......
This diff is collapsed.
This diff is collapsed.
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