Commit 9eec32f8 authored by Quentin Khan's avatar Quentin Khan
Browse files

Update Lagrange Kernel interface to use the new data organisation

 - Divide the FChebCell inner data layout into two sub-types:
   multipole_t and local_expansion_t. Two class attribute are accessible
   through the `getMultipoledata` and `getLocalExpansionData` methods.

 - Change FMM operators signature to take advantage of the new
   layout. The operators only take as parameters the data they may need
   instead of the whole cell. For instance (simplified):

   void M2M(CellClass* parent, CellClass** children);

   becomes

   void M2M(multipole_t*  parent_m,    symbolic_data_t* parent_s,
            multipole_t** children_ms, symbolic_data_t* children_ss);
parent b26fd3c1
...@@ -58,7 +58,9 @@ public: ...@@ -58,7 +58,9 @@ public:
FComplex<FReal>* getTransformedMultipole(const int inRhs) FComplex<FReal>* getTransformedMultipole(const int inRhs)
{ return this->transformed_multipole_exp + inRhs*TransformedVectorSize; } { return this->transformed_multipole_exp + inRhs*TransformedVectorSize; }
constexpr int getVectorSize() const {
return VectorSize;
}
// to extend FAbstractSendable // to extend FAbstractSendable
template <class BufferWriterClass> template <class BufferWriterClass>
...@@ -91,6 +93,10 @@ public: ...@@ -91,6 +93,10 @@ public:
FComplex<FReal>* getTransformedLocal(const int inRhs) FComplex<FReal>* getTransformedLocal(const int inRhs)
{ return this->transformed_local_exp + inRhs*TransformedVectorSize; } { return this->transformed_local_exp + inRhs*TransformedVectorSize; }
constexpr int getVectorSize() const {
return VectorSize;
}
// to extend FAbstractSendable // to extend FAbstractSendable
template <class BufferWriterClass> template <class BufferWriterClass>
void serialize(BufferWriterClass& buffer) const{ void serialize(BufferWriterClass& buffer) const{
...@@ -105,8 +111,8 @@ public: ...@@ -105,8 +111,8 @@ public:
} }
}; };
multipole_t m_data; multipole_t m_data {};
local_expansion_t l_data; local_expansion_t l_data {};
multipole_t* up = &m_data; multipole_t* up = &m_data;
local_expansion_t* down = &l_data; local_expansion_t* down = &l_data;
......
...@@ -109,10 +109,11 @@ public: ...@@ -109,10 +109,11 @@ public:
} }
template<class SymbolicData>
void M2M(typename CellClass::multipole_t* const FRestrict ParentMultipole, void M2M(typename CellClass::multipole_t* const FRestrict ParentMultipole,
const SymbolicData* const /*ParentSymb*/,
const typename CellClass::multipole_t*const FRestrict *const FRestrict ChildMultipoles, const typename CellClass::multipole_t*const FRestrict *const FRestrict ChildMultipoles,
const int /*TreeLevel*/) const SymbolicData* const /*ChildSymbs*/[])
override
{ {
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply Sy // 1) apply Sy
...@@ -132,13 +133,15 @@ public: ...@@ -132,13 +133,15 @@ public:
template<class SymbolicData>
void M2L(typename CellClass::local_expansion_t * const FRestrict TargetExpansion, void M2L(typename CellClass::local_expansion_t * const FRestrict TargetExpansion,
const SymbolicData* const TargetSymb,
const typename CellClass::multipole_t * const FRestrict SourceMultipoles[], const typename CellClass::multipole_t * const FRestrict SourceMultipoles[],
const int neighborPositions[], const int inSize, const int TreeLevel) const SymbolicData* const FRestrict /*SourceSymbs*/[],
override const int neighborPositions[],
const int inSize)
{ {
const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel))); const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(1 << TargetSymb->getLevel()));
const FReal scale(MatrixKernel->getScaleFactor(CellWidth)); const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
...@@ -146,7 +149,7 @@ public: ...@@ -146,7 +149,7 @@ public:
for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){ for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
const int idxNeigh = neighborPositions[idxExistingNeigh]; const int idxNeigh = neighborPositions[idxExistingNeigh];
M2LHandler.applyFC(idxNeigh, TreeLevel, scale, M2LHandler.applyFC(idxNeigh, static_cast<int>(TargetSymb->getLevel()), scale,
SourceMultipoles[idxExistingNeigh]->getTransformedMultipole(idxRhs), SourceMultipoles[idxExistingNeigh]->getTransformedMultipole(idxRhs),
TransformedLocalExpansion); TransformedLocalExpansion);
} }
...@@ -154,10 +157,11 @@ public: ...@@ -154,10 +157,11 @@ public:
} }
template<class SymbolicData>
void L2L(const typename CellClass::local_expansion_t * const FRestrict ParentExpansion, void L2L(const typename CellClass::local_expansion_t * const FRestrict ParentExpansion,
typename CellClass::local_expansion_t * FRestrict * const FRestrict ChildExpansion, const SymbolicData* const /*ParentSymb*/,
const int /*TreeLevel*/) typename CellClass::local_expansion_t * FRestrict * const FRestrict ChildExpansions,
override const SymbolicData* const /*ChildSymbs*/[])
{ {
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) Apply Inverse Discete Fourier Transform // 1) Apply Inverse Discete Fourier Transform
...@@ -168,9 +172,9 @@ public: ...@@ -168,9 +172,9 @@ public:
const_cast<FReal*>(ParentExpansion->getLocal(idxRhs)),localExp); const_cast<FReal*>(ParentExpansion->getLocal(idxRhs)),localExp);
// 2) apply Sx // 2) apply Sx
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){ for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
if (ChildExpansion[ChildIndex]){ if (ChildExpansions[ChildIndex]){
AbstractBaseClass::Interpolator->applyL2L( AbstractBaseClass::Interpolator->applyL2L(
ChildIndex, localExp, ChildExpansion[ChildIndex]->getLocal(idxRhs)); ChildIndex, localExp, ChildExpansions[ChildIndex]->getLocal(idxRhs));
} }
} }
} }
...@@ -178,7 +182,7 @@ public: ...@@ -178,7 +182,7 @@ public:
template<class SymbolicData> template<class SymbolicData>
void L2P(const typename CellClass::local_expansion_t* const LeafCell, void L2P(const typename CellClass::local_expansion_t* const LeafExpansion,
const SymbolicData* const LeafSymbData, const SymbolicData* const LeafSymbData,
ContainerClass* const TargetParticles) ContainerClass* const TargetParticles)
{ {
...@@ -192,9 +196,9 @@ public: ...@@ -192,9 +196,9 @@ public:
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){ for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) Apply Inverse Discete Fourier Transform // 1) Apply Inverse Discete Fourier Transform
this->M2LHandler.unapplyZeroPaddingAndDFT(LeafCell->getTransformedLocal(idxRhs), this->M2LHandler.unapplyZeroPaddingAndDFT(LeafExpansion->getTransformedLocal(idxRhs),
localExp + idxRhs*AbstractBaseClass::nnodes); localExp + idxRhs*AbstractBaseClass::nnodes);
FBlas::add(AbstractBaseClass::nnodes,const_cast<FReal*>(LeafCell->getLocal(idxRhs)),localExp + idxRhs*AbstractBaseClass::nnodes); FBlas::add(AbstractBaseClass::nnodes,const_cast<FReal*>(LeafExpansion->getLocal(idxRhs)),localExp + idxRhs*AbstractBaseClass::nnodes);
} }
// 2.a) apply Sx // 2.a) apply Sx
AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter, leafBoxWidth, AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter, leafBoxWidth,
......
...@@ -109,10 +109,12 @@ public: ...@@ -109,10 +109,12 @@ public:
{ } { }
void P2M(CellClass* const LeafCell, template<class SymbolicData>
void P2M(typename CellClass::multipole_t* const LeafMultipole,
const SymbolicData* const LeafSymbData,
const ContainerClass* const SourceParticles) const ContainerClass* const SourceParticles)
{ {
const FPoint<FReal> LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate())); const FPoint<FReal> LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafSymbData->getCoordinate()));
const FReal ExtendedLeafCellWidth(AbstractBaseClass::BoxWidthLeaf const FReal ExtendedLeafCellWidth(AbstractBaseClass::BoxWidthLeaf
+ AbstractBaseClass::BoxWidthExtension); + AbstractBaseClass::BoxWidthExtension);
...@@ -120,24 +122,26 @@ public: ...@@ -120,24 +122,26 @@ public:
// 1) apply Sy // 1) apply Sy
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, ExtendedLeafCellWidth, AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, ExtendedLeafCellWidth,
LeafCell->getMultipole(idxV*nRhs), SourceParticles); LeafMultipole->getMultipole(idxV*nRhs), SourceParticles);
for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){ for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){
// update multipole index // update multipole index
int idxMul = idxV*nRhs + idxRhs; int idxMul = idxV*nRhs + idxRhs;
// 2) apply Discrete Fourier Transform // 2) apply Discrete Fourier Transform
M2LHandler.applyZeroPaddingAndDFT(LeafCell->getMultipole(idxMul), M2LHandler.applyZeroPaddingAndDFT(LeafMultipole->getMultipole(idxMul),
LeafCell->getTransformedMultipole(idxMul)); LeafMultipole->getTransformedMultipole(idxMul));
} }
}// NVALS }// NVALS
} }
void M2M(CellClass* const FRestrict ParentCell, template<class SymbolicData>
const CellClass*const FRestrict *const FRestrict ChildCells, void M2M(typename CellClass::multipole_t* const FRestrict ParentMultipole,
const int TreeLevel) const SymbolicData* const ParentSymb,
const typename CellClass::multipole_t*const FRestrict *const FRestrict ChildMultipoles,
const SymbolicData* const /*ChildSymbs*/[])
{ {
for(int idxV = 0 ; idxV < NVALS ; ++idxV){ for(int idxV = 0 ; idxV < NVALS ; ++idxV){
for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){ for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){
...@@ -145,26 +149,33 @@ public: ...@@ -145,26 +149,33 @@ public:
int idxMul = idxV*nRhs + idxRhs; int idxMul = idxV*nRhs + idxRhs;
// 1) apply Sy // 1) apply Sy
FBlas::scal(AbstractBaseClass::nnodes, FReal(0.), ParentCell->getMultipole(idxMul)); FBlas::scal(AbstractBaseClass::nnodes, FReal(0.), ParentMultipole->getMultipole(idxMul));
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){ for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
if (ChildCells[ChildIndex]){ if (ChildMultipoles[ChildIndex]){
AbstractBaseClass::Interpolator->applyM2M(ChildIndex, AbstractBaseClass::Interpolator->applyM2M(
ChildCells[ChildIndex]->getMultipole(idxMul), ChildIndex,
ParentCell->getMultipole(idxMul), ChildMultipoles[ChildIndex]->getMultipole(idxMul),
TreeLevel/*Cell width extension specific*/); ParentMultipole->getMultipole(idxMul),
ParentSymb->getLevel()/*Cell width extension specific*/);
} }
} }
// 2) Apply Discete Fourier Transform // 2) Apply Discete Fourier Transform
M2LHandler.applyZeroPaddingAndDFT(ParentCell->getMultipole(idxMul), M2LHandler.applyZeroPaddingAndDFT(ParentMultipole->getMultipole(idxMul),
ParentCell->getTransformedMultipole(idxMul)); ParentMultipole->getTransformedMultipole(idxMul));
} }
}// NVALS }// NVALS
} }
void M2L(CellClass* const FRestrict TargetCell, const CellClass* SourceCells[], template<class SymbolicData>
const int neighborPositions[], const int inSize, const int TreeLevel) override { void M2L(typename CellClass::local_expansion_t * const FRestrict TargetExpansion,
const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel))); const SymbolicData* const TargetSymb,
const typename CellClass::multipole_t * const FRestrict SourceMultipoles[],
const SymbolicData* const FRestrict /*SourceSymbs*/[],
const int neighborPositions[],
const int inSize)
{
const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(1 << TargetSymb->getLevel()));
const FReal ExtendedCellWidth(CellWidth + AbstractBaseClass::BoxWidthExtension); const FReal ExtendedCellWidth(CellWidth + AbstractBaseClass::BoxWidthExtension);
const FReal scale(MatrixKernel->getScaleFactor(ExtendedCellWidth)); const FReal scale(MatrixKernel->getScaleFactor(ExtendedCellWidth));
...@@ -175,7 +186,7 @@ public: ...@@ -175,7 +186,7 @@ public:
const int idxLoc = idxV*nLhs + idxLhs; const int idxLoc = idxV*nLhs + idxLhs;
// load transformed local expansion // load transformed local expansion
FComplex<FReal> *const TransformedLocalExpansion = TargetCell->getTransformedLocal(idxLoc); FComplex<FReal> *const TransformedLocalExpansion = TargetExpansion->getTransformedLocal(idxLoc);
// update idxRhs // update idxRhs
const int idxRhs = idxLhs % nPV; const int idxRhs = idxLhs % nPV;
...@@ -189,43 +200,48 @@ public: ...@@ -189,43 +200,48 @@ public:
for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){ for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
const int idx = neighborPositions[idxExistingNeigh]; const int idx = neighborPositions[idxExistingNeigh];
M2LHandler.applyFC(idx, TreeLevel, scale, d, M2LHandler.applyFC(idx, static_cast<int>(TargetSymb->getLevel()), scale, d,
SourceCells[idxExistingNeigh]->getTransformedMultipole(idxMul), SourceMultipoles[idxExistingNeigh]->getTransformedMultipole(idxMul),
TransformedLocalExpansion); TransformedLocalExpansion);
} }
}// NLHS=NPOT*NPV }// NLHS=NPOT*NPV
}// NVALS }// NVALS
} }
void L2L(const CellClass* const FRestrict ParentCell, template<class SymbolicData>
CellClass* FRestrict *const FRestrict ChildCells, void L2L(const typename CellClass::local_expansion_t * const FRestrict ParentExpansion,
const int TreeLevel) const SymbolicData* const ParentSymb,
typename CellClass::local_expansion_t * FRestrict * const FRestrict ChildExpansions,
const SymbolicData* const /*ChildSymbs*/[])
{ {
for(int idxV = 0 ; idxV < NVALS ; ++idxV){ for(int idxV = 0 ; idxV < NVALS ; ++idxV){
for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){ for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
int idxLoc = idxV*nLhs + idxLhs; int idxLoc = idxV*nLhs + idxLhs;
// 1) Apply Inverse Discete Fourier Transform // 1) Apply Inverse Discete Fourier Transform
M2LHandler.unapplyZeroPaddingAndDFT(ParentCell->getTransformedLocal(idxLoc), M2LHandler.unapplyZeroPaddingAndDFT(
const_cast<CellClass*>(ParentCell)->getLocal(idxLoc)); ParentExpansion->getTransformedLocal(idxLoc),
const_cast<typename CellClass::local_expansion_t*>(ParentExpansion)->getLocal(idxLoc));
// 2) apply Sx // 2) apply Sx
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){ for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
if (ChildCells[ChildIndex]){ if (ChildExpansions[ChildIndex]){
AbstractBaseClass::Interpolator->applyL2L(ChildIndex, AbstractBaseClass::Interpolator->applyL2L(
ParentCell->getLocal(idxLoc), ChildIndex,
ChildCells[ChildIndex]->getLocal(idxLoc), ParentExpansion->getLocal(idxLoc),
TreeLevel/*Cell width extension specific*/); ChildExpansions[ChildIndex]->getLocal(idxLoc),
static_cast<int>(ParentSymb->getLevel())/*Cell width extension specific*/);
} }
} }
} }
}// NVALS }// NVALS
} }
void L2P(const CellClass* const LeafCell, template<class SymbolicData>
void L2P(const typename CellClass::local_expansion_t* const LeafExpansion,
const SymbolicData* const LeafSymbData,
ContainerClass* const TargetParticles) ContainerClass* const TargetParticles)
{ {
const FPoint<FReal> LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate())); const FPoint<FReal> LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafSymbData->getCoordinate()));
const FReal ExtendedLeafCellWidth(AbstractBaseClass::BoxWidthLeaf const FReal ExtendedLeafCellWidth(AbstractBaseClass::BoxWidthLeaf
+ AbstractBaseClass::BoxWidthExtension); + AbstractBaseClass::BoxWidthExtension);
...@@ -233,18 +249,19 @@ public: ...@@ -233,18 +249,19 @@ public:
for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){ for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
int idxLoc = idxV*nLhs + idxLhs; int idxLoc = idxV*nLhs + idxLhs;
// 1) Apply Inverse Discete Fourier Transform // 1) Apply Inverse Discete Fourier Transform
M2LHandler.unapplyZeroPaddingAndDFT(LeafCell->getTransformedLocal(idxLoc), M2LHandler.unapplyZeroPaddingAndDFT(
const_cast<CellClass*>(LeafCell)->getLocal(idxLoc)); LeafExpansion->getTransformedLocal(idxLoc),
const_cast<typename CellClass::local_expansion_t*>(LeafExpansion)->getLocal(idxLoc));
} }
// 2.a) apply Sx // 2.a) apply Sx
AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter, ExtendedLeafCellWidth, AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter, ExtendedLeafCellWidth,
LeafCell->getLocal(idxV*nLhs), TargetParticles); LeafExpansion->getLocal(idxV*nLhs), TargetParticles);
// 2.b) apply Px (grad Sx) // 2.b) apply Px (grad Sx)
AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter, ExtendedLeafCellWidth, AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter, ExtendedLeafCellWidth,
LeafCell->getLocal(idxV*nLhs), TargetParticles); LeafExpansion->getLocal(idxV*nLhs), TargetParticles);
}// NVALS }// NVALS
} }
......
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