Commit aafafbd1 authored by Quentin Khan's avatar Quentin Khan

Update Taylor Kernel interface to use the new data organisation

 - Divide the FTaylorCell 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 4b333838
......@@ -4,13 +4,13 @@
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.
//
// abiding by the rules of distribution of free software.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
#ifndef FTAYLORCELL_HPP
......@@ -32,53 +32,68 @@
template < class FReal, int P, int order>
class FTaylorCell : public FBasicCell, public FAbstractSendable {
protected:
//Size of Multipole Vector
static const int MultipoleSize = ((P+1)*(P+2)*(P+3))*order/6;
//Size of Local Vector
static const int LocalSize = ((P+1)*(P+2)*(P+3))*order/6;
//Multipole vector
FReal multipole_exp[MultipoleSize];
//Local vector
FReal local_exp[LocalSize];
template<std::size_t S, class Tag>
struct expansion_impl {
enum {Size = S};
FReal exp[Size];
FReal* get() noexcept {
return exp;
}
const FReal* get() const noexcept {
return exp;
}
void reset() noexcept {
for(int idx = 0; idx < Size; ++idx) {
exp[idx].setRealImag(FReal(0.0), FReal(0.0));
}
}
int getSize() const noexcept {
return Size;
}
template<class BufferWriterClass>
void serialize(BufferWriterClass& buffer) const {
buffer.write(exp, Size);
}
template<class BufferReaderClass>
void deserialize(BufferReaderClass& buffer) {
buffer.fillArray(exp, Size);
}
};
public:
/**
*Default Constructor
*/
FTaylorCell(){
FMemUtils::memset(multipole_exp,0,MultipoleSize*sizeof(FReal(0)));
FMemUtils::memset(local_exp,0,LocalSize*sizeof(FReal(0)));
}
//Get multipole Vector for setting
FReal * getMultipole(void)
{
return multipole_exp;
}
using multipole_t = expansion_impl<((P+1)*(P+2)*(P+3))*order/6, class multipole_tag>;
using local_expansion_t = expansion_impl<((P+1)*(P+2)*(P+3))*order/6, class local_expansion_tag>;
//Get multipole Vector for reading
const FReal * getMultipole(void) const
{
return multipole_exp;
}
protected:
//Get local Vector
FReal * getLocal(void)
{
return local_exp;
}
multipole_t m_data;
local_expansion_t l_data;
//Get local Vector for reading
const FReal * getLocal(void) const
{
return local_exp;
public:
const multipole_t& getMultipoleData() const noexcept {
return m_data;
}
multipole_t& getMultipoleData() {
return m_data;
}
const local_expansion_t& getLocalExpansionData() const noexcept {
return l_data;
}
local_expansion_t& getLocalExpansionData() {
return l_data;
}
/** Make it like the begining */
void resetToInitialState(){
FMemUtils::memset(multipole_exp,0,MultipoleSize*sizeof(FReal(0)));
FMemUtils::memset(local_exp,0,LocalSize*sizeof(FReal(0)));
m_data.reset();
l_data.reset();
}
......@@ -87,28 +102,28 @@ public:
///////////////////////////////////////////////////////
template <class BufferWriterClass>
void serializeUp(BufferWriterClass& buffer) const{
buffer.write(multipole_exp, MultipoleSize);
m_data.serialize(buffer);
}
template <class BufferReaderClass>
void deserializeUp(BufferReaderClass& buffer){
buffer.fillArray(multipole_exp, MultipoleSize);
m_data.deserialize(buffer);
}
template <class BufferWriterClass>
void serializeDown(BufferWriterClass& buffer) const{
buffer.write(local_exp, LocalSize);
l_data.serialize(buffer);
}
template <class BufferReaderClass>
void deserializeDown(BufferReaderClass& buffer){
buffer.fillArray(local_exp, LocalSize);
l_data.deserialize(buffer);
}
FSize getSavedSizeUp() const {
return ((FSize) sizeof(FReal) * (MultipoleSize));
return ((FSize) sizeof(FReal)) * (multipole_t::Size);
}
FSize getSavedSizeDown() const {
return ((FSize) sizeof(FReal) * (LocalSize));
return ((FSize) sizeof(FReal)) * (local_expansion_t::Size);
}
///////////////////////////////////////////////////////
......@@ -117,18 +132,18 @@ public:
template <class BufferWriterClass>
void save(BufferWriterClass& buffer) const{
FBasicCell::save(buffer);
buffer.write(multipole_exp, MultipoleSize);
buffer.write(local_exp, LocalSize);
m_data.serialize();
l_data.serialize();
}
template <class BufferReaderClass>
void restore(BufferReaderClass& buffer){
FBasicCell::restore(buffer);
buffer.fillArray(multipole_exp, MultipoleSize);
buffer.fillArray(local_exp, LocalSize);
m_data.deserialize();
l_data.deserialize();
}
FSize getSavedSize() const {
return FSize((sizeof(FReal) * (MultipoleSize + LocalSize)
return FSize((sizeof(FReal) * (multipole_t::Size + local_expansion_t::Size)
+ FBasicCell::getSavedSize()));
}
};
......
......@@ -4,13 +4,13 @@
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.
//
// abiding by the rules of distribution of free software.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
#ifndef FTAYLORKERNEL_HPP
......@@ -92,13 +92,18 @@ private:
FReal(coordinate.getZ()) * widthAtLeafLevel + widthAtLeafLevelDiv2 + boxCorner.getZ());
}
template<class SymbolicData>
FPoint<FReal> getLeafCenter(const SymbolicData* symb) const {
return getLeafCenter(symb->getCoordinate());
}
/**
* @brief Return the position of the center of a cell from its tree
* coordinate
* @param FTreeCoordinate
* @param inLevel the current level of Cell
*/
FPoint<FReal> getCellCenter(const FTreeCoordinate coordinate, int inLevel)
* @brief Return the position of the center of a cell from its tree
* coordinate
* @param FTreeCoordinate
* @param inLevel the current level of Cell
*/
FPoint<FReal> getCellCenter(const FTreeCoordinate coordinate, int inLevel) const
{
//Set the boxes width needed
FReal widthAtCurrentLevel = widthAtLeafLevel*FReal(1 << (treeHeight-(inLevel+1)));
......@@ -118,6 +123,15 @@ private:
return cCenter;
}
/**
* @brief Return the position of the center of a cell from its symbolic data
* @param symb The cell symbolic data
*/
template<class SymbolicData>
FPoint<FReal> getCellCenter(const SymbolicData* symb) const {
return getCellCenter(symb->getCoordinate(), static_cast<int>(symb->getLevel()));
}
/**
* @brief Incrementation of powers in Taylor expansion
......@@ -542,13 +556,15 @@ public:
* \f]
* where \f$x_c\f$ is the centre of the cell and \f$x_j\f$ the \f$j^{th}\f$ particles and \f$q_j\f$ its charge and \f$N\f$ the particle number.
*/
void P2M(CellClass* const pole,
const ContainerClass* const particles) override
template<class SymbolicData>
void P2M(typename CellClass::multipole_t* const LeafMultipole,
const SymbolicData* const LeafSymb,
const ContainerClass* const particles)
{
//Copying cell center position once and for all
const FPoint<FReal>& cellCenter = getLeafCenter(pole->getCoordinate());
const FPoint<FReal>& cellCenter = getLeafCenter(LeafSymb->getCoordinate());
FReal * FRestrict multipole = pole->getMultipole();
FReal * FRestrict multipole = LeafMultipole->get();
FMemUtils::memset(multipole,0,SizeVector*sizeof(FReal(0.0)));
FReal multipole2[SizeVector] ;
......@@ -606,30 +622,32 @@ public:
*
*
*/
void M2M(CellClass* const FRestrict pole,
const CellClass*const FRestrict *const FRestrict child,
const int inLevel)override
template<class SymbolicData>
void M2M(typename CellClass::multipole_t* const FRestrict ParentMultipole,
const SymbolicData* const ParentSymb,
const typename CellClass::multipole_t*const FRestrict *const FRestrict ChildMultipoles,
const SymbolicData* const ChildSymbs[])
{
//Powers of expansions
int a=0,b=0,c=0;
//Indexes of powers
int idx_a,idx_b,idx_c;
//Distance from current child to parent
//Distance from current child To parent
FReal dx = 0.0;
FReal dy = 0.0;
FReal dz = 0.0;
//Center point of parent cell
const FPoint<FReal>& cellCenter = getCellCenter(pole->getCoordinate(),inLevel);
FReal * FRestrict mult = pole->getMultipole();
const FPoint<FReal>& cellCenter = getCellCenter(ParentSymb);
FReal * FRestrict mult = ParentMultipole->get();
//Iteration over the eight children
int idxChild;
FReal coeff;
for(idxChild=0 ; idxChild<8 ; ++idxChild)
{
if(child[idxChild]){
const FPoint<FReal>& childCenter = getCellCenter(child[idxChild]->getCoordinate(),inLevel+1);
const FReal * FRestrict multChild = child[idxChild]->getMultipole();
if(ChildMultipoles[idxChild]){
const FPoint<FReal>& childCenter = getCellCenter(ChildSymbs[idxChild]);
const FReal * FRestrict multChild = ChildMultipoles[idxChild]->get();
//Set the distance between centers of cells
dx = cellCenter.getX() - childCenter.getX();
......@@ -708,24 +726,30 @@ public:
}
/**
*@brief Convert the multipole expansion into local expansion The
* operator do not use symmetries.
*
* Formula : \f[ L_{\mathbf{n}}^{c} = \frac{|n|!}{n! n!}
* \sum_{\mathbf{k}=0}^{p} \left [ M_\mathbf{k}^c \,
* \Psi_{\mathbf{,n+k}}( \mathbf{x}_c^{target})\right ] \f]
* and \f[ \Psi_{\mathbf{,i}}^{c}(\mathbf{x}) =
* \frac{\partial^i}{\partial x^i} \frac{1}{|x-x_c^{src}|} = \frac{\partial^{i_1}}{\partial x_1^{i_1}} \frac{\partial^{i_2}}{\partial x_2^{i_2}} \frac{\partial^{i_3}}{\partial x_3^{i_3}} \frac{1}{|x-x_c^{src}|}\f]
*
* Where \f$x_c^{src}\f$ is the centre of the cell where the
* multiplole are considered,\f$x_c^{target}\f$ is the centre of the
* current cell. The cell where we compute the local expansion.
*
*/
void M2L(CellClass* const FRestrict inLocal, const CellClass* distantNeighbors[],
const int neighborPositions[], const int inSize, const int inLevel) override {
* @brief Convert the multipole expansion into local expansion The
* operator does not use symetries.
*
* Formula : \f[ L_{\mathbf{n}}^{c} = \frac{|n|!}{n! n!}
* \sum_{\mathbf{k}=0}^{p} \left [ M_\mathbf{k}^c \,
* \Psi_{\mathbf{,n+k}}( \mathbf{x}_c^{target})\right ] \f]
* and \f[ \Psi_{\mathbf{,i}}^{c}(\mathbf{x}) =
* \frac{\partial^i}{\partial x^i} \frac{1}{|x-x_c^{src}|} = \frac{\partial^{i_1}}{\partial x_1^{i_1}} \frac{\partial^{i_2}}{\partial x_2^{i_2}} \frac{\partial^{i_3}}{\partial x_3^{i_3}} \frac{1}{|x-x_c^{src}|}\f]
*
* Where \f$x_c^{src}\f$ is the centre of the cell where the
* multiplole are considered,\f$x_c^{target}\f$ is the centre of the
* current cell. The cell where we compute the local expansion.
*
*/
template<class SymbolicData>
void M2L(typename CellClass::local_expansion_t* const FRestrict TargetLocalExp,
const SymbolicData* const TargetSymb,
const typename CellClass::multipole_t* const FRestrict distantNeighbors[],
const SymbolicData* const FRestrict /*SourceSymbs*/[],
const int neighborPositions[],
const int inSize)
{
//Local expansion to be filled
FReal * FRestrict iterLocal = inLocal->getLocal();
FReal * FRestrict iterLocal = TargetLocalExp->get();
//index for distant neighbors
FReal tmp1;
int * ind;
......@@ -736,9 +760,9 @@ public:
const int idxNeigh = neighborPositions[idxExistingNeigh];
//Get the preComputed values for the derivative
psiVector = M2LpreComputedDerivatives[inLevel][idxNeigh];
psiVector = M2LpreComputedDerivatives[static_cast<int>(TargetSymb->getLevel())][idxNeigh];
//Multipole to be read
multipole = distantNeighbors[idxExistingNeigh]->getMultipole();
multipole = distantNeighbors[idxExistingNeigh]->get();
//Iterating over local array : n
for(int i=0 ; i<SizeVector ; ++i)
......@@ -760,27 +784,28 @@ public:
/**
*@brief Translate the local expansion of parent cell to child cell
*
* One need to translate the local expansion on a father cell
* centered in \f$\mathbf{x}_p\f$ to its eight daughters centered in
* \f$\mathbf{x}_p\f$ .
*
* Local expansion for the daughters will be :
* \f[ \sum_{\mathbf{k}=0}^{|k|<P} L_k * (\mathbf{x}-\mathbf{x}_f)^{\mathbf{k}} \f]
* with :
*
*\f[ L_{n} = \sum_{\mathbf{k}=\mathbf{n}}^{|k|<P} C_{\mathbf{k}}^{\mathbf{n}} * (\mathbf{x}_f-\mathbf{x}_p)^{\mathbf{k}-\mathbf{n}} \f]
*/
void L2L(const CellClass* const FRestrict fatherCell,
CellClass* FRestrict * const FRestrict childCell,
const int inLevel) override
*@brief Translate the local expansion of parent cell to child cell
*
* One need to translate the local expansion on a father cell
* centered in \f$\mathbf{x}_p\f$ to its eight daughters centered in
* \f$\mathbf{x}_p\f$ .
*
* Local expansion for the daughters will be :
* \f[ \sum_{\mathbf{k}=0}^{|k|<P} L_k * (\mathbf{x}-\mathbf{x}_f)^{\mathbf{k}} \f]
* with :
*
*\f[ L_{n} = \sum_{\mathbf{k}=\mathbf{n}}^{|k|<P} C_{\mathbf{k}}^{\mathbf{n}} * (\mathbf{x}_f-\mathbf{x}_p)^{\mathbf{k}-\mathbf{n}} \f]
*/
template<class SymbolicData>
void L2L(const typename CellClass::local_expansion_t * const FRestrict ParentExpansion,
const SymbolicData* const ParentSymb,
typename CellClass::local_expansion_t * FRestrict * const FRestrict ChildExpansions,
const SymbolicData* const ChildSymbs[])
{
FPoint<FReal> locCenter = getCellCenter(fatherCell->getCoordinate(),inLevel);
FPoint<FReal> locCenter = getCellCenter(ParentSymb);
// Get father local expansion
const FReal* FRestrict fatherExpansion = fatherCell->getLocal() ;
const FReal* FRestrict fatherExpansion = ParentExpansion->get();
int idxFatherLoc; //index of Father local expansion to be read.
FReal dx, dy, dz, coeff;
int ap, bp, cp, af, bf, cf; //Indexes of expansion for father and child.
......@@ -788,10 +813,10 @@ public:
// For all children
for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
if(childCell[idxChild]){ // if child exists
if(ChildExpansions[idxChild]){ // if child exists
FReal* FRestrict childExpansion = childCell[idxChild]->getLocal() ;
const FPoint<FReal>& childCenter = getCellCenter(childCell[idxChild]->getCoordinate(),inLevel+1);
FReal* FRestrict childExpansion = ChildExpansions[idxChild]->get() ;
const FPoint<FReal>& childCenter = getCellCenter(ChildSymbs[idxChild]);
//Set the distance between centers of cells
// Child - father
......@@ -866,26 +891,28 @@ public:
* where \f$x_c\f$ is the centre of the local cell and \f$x_j\f$ the
* \f$j^{th}\f$ particles and \f$q_j\f$ its charge.
*/
void L2P(const CellClass* const local,
ContainerClass* const particles) override
template<class SymbolicData>
void L2P(const typename CellClass::local_expansion_t* const LeafExpansion,
const SymbolicData* const LeafSymb,
ContainerClass* const TargetParticles)
{
FPoint<FReal> locCenter = getLeafCenter(local->getCoordinate());
FPoint<FReal> locCenter = getLeafCenter(LeafSymb);
//Iterator over particles
FSize nbPart = particles->getNbParticles();
FSize nbPart = TargetParticles->getNbParticles();
const FReal * iterLocal = local->getLocal();
const FReal * const * positions = particles->getPositions();
const FReal * iterLocal = LeafExpansion->get();
const FReal * const * positions = TargetParticles->getPositions();
const FReal * posX = positions[0];
const FReal * posY = positions[1];
const FReal * posZ = positions[2];
FReal * forceX = particles->getForcesX();
FReal * forceY = particles->getForcesY();
FReal * forceZ = particles->getForcesZ();
FReal * forceX = TargetParticles->getForcesX();
FReal * forceY = TargetParticles->getForcesY();
FReal * forceZ = TargetParticles->getForcesZ();
FReal * targetsPotentials = particles->getPotentials();
FReal * targetsPotentials = TargetParticles->getPotentials();
FReal * phyValues = particles->getPhysicalValues();
FReal * phyValues = TargetParticles->getPhysicalValues();
//Iteration over particles
for(FSize i=0 ; i<nbPart ; ++i){
......@@ -993,6 +1020,18 @@ public:
FP2PRT<FReal>::template FullRemote<ContainerClass>(inTargets,inNeighbors,inSize);
}
// template<class SymbolicData>
// void P2L(typename CellClass::local_expansion_t* const local,
// const SymbolicData * const symb,
// const ContainerClass* const particles,
// const SymbolicData * const /*source_symb*/)
// {}
// template<class SymbolicData>
// void M2P(const typename CellClass::multipole_t* const pole,
// const SymbolicData* const symb,
// ContainerClass* const particles,
// const SymbolicData * const /*target_symb*/)
// {}
};
#endif
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