Commit e6db870b authored by BRAMAS Berenger's avatar BRAMAS Berenger

Update the split of particles and cell component for the starpu version

parents 204130fd e65ccfaf
......@@ -365,19 +365,6 @@ if (MORSE_DISTRIB_DIR OR EXISTS "${CMAKE_SOURCE_DIR}/CMakeModules/morse/")
if(CMAKE_CXX_COMPILER_ID STREQUAL "Intel")
execute_process(COMMAND ${CMAKE_CXX_COMPILER} -dumpversion OUTPUT_VARIABLE INTEL_VERSION)
message( STATUS "Intel: ${INTEL_VERSION}")
if(APPLE)
string(FIND ${CMAKE_CXX_COMPILER} "icl++" NAME)
if( ${NAME} GREATER 0)
if( ${INTEL_VERSION} EQUAL 15. OR ${INTEL_VERSION} GREATER 15.0.0)
message( STATUS " Intel compiler is icl++ ( version >= 15.0.0)")
set(INTEL_ICL_COMPILER "ON")
else()
message(FATAL_ERROR " Intel compiler should be icl++ ( version >= 15.0.0)")
endif()
else()
message(FATAL_ERROR " Intel compiler should be icl++ ( version >= 15.0.0)")
endif()
endif()
endif()
#
##################################################################
......
......@@ -424,7 +424,7 @@ if (BLA_VENDOR MATCHES "Intel*" OR BLA_VENDOR STREQUAL "All")
set(OMP_LIB "${OMP_iomp5_LIBRARY}")
endif()
if (NOT WIN32)
if (UNIX AND NOT WIN32)
set(LM "-lm")
set(BLAS_COMPILER_FLAGS "")
if (CMAKE_C_COMPILER_ID STREQUAL "Intel" AND NOT BLA_VENDOR STREQUAL "Intel10_64lp_seq")
......
......@@ -435,6 +435,9 @@ if(FFTW_LIBRARIES)
list(APPEND REQUIRED_LDFLAGS "-Wl,--no-as-needed")
endif()
endif()
if(UNIX OR WIN32)
list(APPEND REQUIRED_LIBS "-lm")
endif()
# set required libraries for link
set(CMAKE_REQUIRED_INCLUDES "${REQUIRED_INCDIRS}")
......
......@@ -332,7 +332,7 @@ if(BLAS_FOUND)
#intel lapack
if (BLA_VENDOR MATCHES "Intel" OR BLA_VENDOR STREQUAL "All")
if (NOT WIN32)
if (UNIX AND NOT WIN32)
set(LM "-lm")
endif ()
if (_LANGUAGES_ MATCHES C OR _LANGUAGES_ MATCHES CXX)
......
......@@ -309,7 +309,9 @@ if (LAPACK_FOUND)
list(APPEND REQUIRED_LIBS "-lifcore")
endif()
# m
list(APPEND REQUIRED_LIBS "-lm")
if(UNIX OR WIN32)
list(APPEND REQUIRED_LIBS "-lm")
endif()
# set required libraries for link
set(CMAKE_REQUIRED_INCLUDES "${REQUIRED_INCDIRS}")
......
......@@ -360,7 +360,7 @@ if (BLA_VENDOR STREQUAL "Generic" OR
endif ()
#intel scalapack
if (BLA_VENDOR MATCHES "Intel" OR BLA_VENDOR STREQUAL "All")
if (NOT WIN32)
if (UNIX AND NOT WIN32)
set(LM "-lm")
endif ()
if (_LANGUAGES_ MATCHES C OR _LANGUAGES_ MATCHES CXX)
......
@article{fong09a,
author = {Fong, W and Darve, E},
doi = {DOI: 10.1016/j.jcp.2009.08.031},
file = {:Users/coulaud/Recherche/Bibliographie/Mendeley/Fong, Darve - 2009 - The black-box fast multipole method.pdf:pdf},
issn = {0021-9991},
journal = {Journal of Computational Physics},
keywords = {Chebyshev polynomials,Fast multipole method,Interpolation,Singular value decomposition},
mendeley-groups = {FMM},
number = {23},
pages = {8712--8725},
title = {{The black-box fast multipole method}},
url = {http://www.sciencedirect.com/science/article/pii/S0021999109004665},
volume = {228},
year = {2009}
}
@inproceedings{Warren1993,
address = {New York, NY, USA},
author = {Warren, M S and Salmon, J K},
booktitle = {Proceedings of the 1993 ACM/IEEE conference on Supercomputing},
doi = {http://doi.acm.org/10.1145/169627.169640},
file = {:Users/coulaud/Recherche/Bibliographie/Mendeley/Warren, Salmon - 1993 - A hashed Oct-Tree N-body algorithm.pdf:pdf},
isbn = {0-8186-4340-4},
mendeley-groups = {FMM},
pages = {12--21},
publisher = {ACM},
series = {Supercomputing '93},
title = {{A hashed Oct-Tree N-body algorithm}},
url = {http://doi.acm.org/10.1145/169627.169640},
year = {1993}
}
@Book{mason2003chebyshev,
author = {J C Mason and D C Handscomb},
title = {Chebyshev polynomials},
publisher = {Chapman \& Hall/CRC},
year = {2003}
}
This diff is collapsed.
\section{Splitting the outer product in the interpolation polynomial}
\label{sec:splitouter}
We introduce the notation of the interpolation polynomial based on the
notation from \cite{mason2003chebyshev}, the interpolation of a function
$p(x)$ reads as
\begin{equation}
\label{eq:masondef}
p(x) \sim \sum_{n=0}^{\ell-1} \gamma_n c_n T_n(x)
\end{equation}
with $\gamma_0=1$ and $\gamma_n=2$ for $n\ge1$ and with
\begin{equation}
\label{eq:chebcoeff}
c_n = \frac{1}{\ell} \sum_{k=1}^\ell p(\bar x_n) T_k(\bar x_n)
\end{equation}
After inserting Eqn.~(\ref{eq:chebcoeff}) into Eqn.~(\ref{eq:masondef}) we end
up with
\begin{equation}
\label{eq:outprod}
p(x) \sim \frac{1}{\ell} \sum_{n=1}^\ell p(\bar x_n) \sum_{k=0}^{\ell-1}
\gamma_k T_k(\bar x_n) T_k(x)
\end{equation}
In the discrete case, i.e., we have $N$ scattered points $\{x_i\}_{i=1}^N
\subset [-1,1]$, we write Eqn.~\eqref{eq:outprod} in matrix notation as
\begin{equation}
\label{eq:interpolmatnot}
\Mat{p} \sim \Mat{T}_x \Mat{\bar T}_x^\top \Mat{\bar p}
\end{equation}
with the matrices $\Mat{T}_x \in \mathbb{R}^{N\times\ell}$ with
$(\Mat{T}_x)_{ik} = T_{k-1}(x_i)$ and $\Mat{\bar T}_x \in
\mathbb{R}^{\ell\times\ell}$ with $(\Mat{\bar T}_x)_{nk} =
\nicefrac{\gamma_{k-1}}{\ell} \, T_{k-1}(\bar x_n)$. Using this notation the
interpolated kernel function $K(x,y)$ reads as
\begin{equation}
\label{eq:kinterpolmatnot}
\Mat{K} \sim \Mat{T}_x \Mat{\bar T}_x^\top \Mat{\bar K} (\Mat{T}_y
\Mat{\bar T}_y^\top)^\top
\end{equation}
With the low-rank representation $\Mat{\bar K} \sim \Mat{UV}^T$ we compute
$\Mat{\bar U} = \Mat{\bar T}_x^\top \Mat{U}$ and $\Mat{\bar V} = \Mat{\bar
T}_y^\top \Mat{V}$ and write Eqn.~\eqref{eq:kinterpolmatnot} as
\begin{equation}
\label{eq:kinterpolmatnot1}
\Mat{K} \sim \Mat{T}_x \Mat{\bar U} \Mat{\bar V}^\top \Mat{T}_y^\top
\end{equation}
\paragraph{Can we still use symmetries?} Recall the approach we proposed in
Sec.~\ref{sec:m2l}. It exploits the symmetries in the arrangement of the
far-field interactions. As we see from Eqn.~\eqref{eq:permut} in that case we
need permutation matrices $\Mat{P}_t$ for the $t$-th far-field interaction and
Eqn.~\eqref{eq:kinterpolmatnot} becomes
\begin{equation}
\label{eq:symsplit}
\Mat{K}_t \sim \Mat{T}_x \Mat{\bar T}_x^\top (\Mat{P}_t \Mat{\bar K}_{p(t)}
\Mat{P}_t^\top) (\Mat{T}_y \Mat{\bar T}_y^\top)^\top.
\end{equation}
With $\Mat{\bar K}_{p(t)} \sim \Mat{U}_{p(t)} \Mat{V}_{p(t)}^\top$ the
matrices in the outer product in Eqn.~\eqref{eq:kinterpolmatnot1} become
$\Mat{\bar U}_t = \Mat{\bar T}_x^\top \Mat{P}_t \Mat{U}_{p(t)}$ and $\Mat{\bar
V}_t = \Mat{\bar T}_y^\top \Mat{P}_t \Mat{V}_{p(t)}$. Thus, it is not
possible anymore to use the fact that due to symmetries we can express all
$316$ far-field interactions by permutations of $16$ only.
%%% Local Variables:
%%% mode: latex
%%% TeX-master: "main"
%%% End:
\documentclass{minimal}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{tikz}
\usetikzlibrary{trees}
\tikzstyle{coordaxis}=[draw=red!50!black, thick, ->]
\tikzstyle{mapping}=[draw=blue!50!black, thick, ->]
\tikzstyle{connect}=[draw=green!50!black, thick, ->]
\newcommand{\sR}{\mathbb{R}}
\begin{document}
% grow=down,
% every node/.style={draw, circle, thin},
% edge from parent/.style={-latex, thick, draw}
\begin{tikzpicture}[scale=0.6]
\begin{tikzpicture}
[level distance=4cm,
level 1/.style={sibling distance=5cm},
level 2/.style={sibling distance=1.4cm},
level 3/.style={sibling distance=0.4cm},
every node/.style={draw, circle, thin}
]
\node (P) {Root}
child {node (C1) {Q}
child {node (C1A) {T}
child {node (C1A1) {T}}
child {node (C1A2) {U}}
child {node (C1A2) {T}}
child {node (C1A4) {U}}
}
child {node (C1B) {U}
child {node (C1B1) {T}}
child {node (C1B2) {U}}
child {node (C1B2) {T}}
child {node (C1B4) {U}}
}
child {node (C1C) {T}}
child {node (C1D) {U}}
}
child {node (C2) {R}
child {node (C2A) {T}}
child {node (C2B) {U}}
child {node (C2C) {T}}
child {node (C2D) {U}}
}
child {node (C3) {S}}
% child {node (C4) {$C^1_4$}
child {node (C4) {$C^1_4$}
child {node (T) {}
child {node (TT) {}
child {node (TTT) {}
}
}
}
child {node (U) {}}
};
\draw[fill=blue!20] (0,8) rectangle +(8,8);
\draw[fill=blue!20] (8,0) rectangle +(8,8);
\draw[fill=blue!20] (0,0) rectangle +(4,4);
\draw[fill=blue!20] (12,12) rectangle +(4,4);
\path (P) -- coordinate[midway] (PQ) (C1);
\path (P) -- coordinate[midway] (PR) (C2);
\draw[fill=brown!50] (4,0) rectangle +(4,4);
\draw[fill=brown!50] (0,4) rectangle +(4,4);
\draw[fill=brown!50] (4,4) rectangle +(4,4);
\draw[fill=brown!50] (8,8) rectangle +(4,4);
\draw[fill=brown!50] (8,12) rectangle +(4,4);
\draw[fill=brown!50] (12,8) rectangle +(4,4);
%Draw the initial octree
\draw[step = 8] (0,0) grid (16,16);
\draw[step = 4] (8,8) grid(16,16);
\draw[step = 4] (0,8) grid(8,16);
\draw[step = 2] (0,4) grid(4,8);
\draw[step = 2] (4,4) grid (8,8);
\draw[step = 2] (4,0) grid (8,4);
\draw[step = 1] (2,4) grid (4,8);
\draw[step = 1] (0,6) grid (2,8);
\draw[step = 0.5] (2,6) grid (3,7);
\draw (PQ) to[bend right=22] (PR);
\draw[step = 2] (8,8) grid(12,16);
\draw[step = 2] (12,8) grid (16,12);
\draw[step = 1] (12,8) grid (14,12);
\draw[step = 1] (14,8) grid (16,10);
\draw[step = 0.5] (13,9) grid (14,10);
\end{tikzpicture}
%sibling distance=2cm,
\begin{tikzpicture}[level distance=1.5cm, grow=down,
every node/.style={draw, circle, thin},
......@@ -84,7 +82,7 @@ level 1/.style={sibling distance=5cm},
}
}
}
% child {node (C2) {}}
child {node (C2) {}}
% child {node (C3) {}}
child {node [fill=blue, text=white] (C4) {$j$}
child {node (T) {}
......@@ -103,4 +101,22 @@ level 1/.style={sibling distance=5cm},
%\draw (PQ) to[bend right=22] (PR);
\end{tikzpicture}
\begin{tikzpicture}
[level distance=30mm, %level/.style={sibling distance=36mm/#1}]
every node/.style={draw, circle, thin},
% edge from parent/.style={-latex, thick, draw},
level 1/.style={sibling distance=12cm},
level 2/.style={sibling distance=30mm},
level 3/.style={sibling distance=7mm},
scale=0.5, rotate=90 ]
\coordinate
child foreach \x in {0,1,2,3} {node{}
{child foreach \y in {0,1,2,3} {node{}
{child foreach \z in {0,1,2,3} {node{}}} }}};
\end{tikzpicture}
\begin{tikzpicture}
%\node {root} child [red] foreach \name in {1,2} {node {\name}}
\end{tikzpicture}
\end{document}
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
......@@ -8,7 +8,7 @@
template<int NRHS = 1, int NLHS = 1>
class FP2PGroupParticleContainer : public FGroupAttachedLeaf<NRHS, 4*NLHS, FReal> {
typedef FGroupAttachedLeaf<NRHS, 4*NLHS, FReal> Parent;
typedef FGroupAttachedLeaf<NVALS*NRHS, NVALS*4*NLHS, FReal> Parent;
public:
FP2PGroupParticleContainer(){}
......@@ -18,45 +18,101 @@ public:
}
FReal* getPhysicalValues(const int idxRhs = 0){
return Parent::getAttribute(0+idxRhs);
FReal* getPhysicalValues(const int idxVals = 0, const int idxRhs = 0){
return Parent::getAttribute((0+idxRhs)*NVALS+idxVals);
}
const FReal* getPhysicalValues(const int idxVals = 0, const int idxRhs = 0) const {
return Parent::getAttribute((0+idxRhs)*NVALS+idxVals);
}
FReal* getPhysicalValuesArray(const int idxVals = 0, const int idxRhs = 0){
return Parent::getRawData() + ((0+idxRhs)*NVALS+idxVals)*Parent::getLeadingRawData();
}
const FReal* getPhysicalValuesArray(const int idxVals = 0, const int idxRhs = 0) const {
return Parent::getRawData() + ((0+idxRhs)*NVALS+idxVals)*Parent::getLeadingRawData();
}
int getLeadingDimension(){
return Parent::getLeadingRawData();
}
FReal* getPotentials(const int idxVals = 0, const int idxLhs = 0){
return Parent::getAttribute((NRHS+idxLhs)*NVALS+idxVals);
}
const FReal* getPotentials(const int idxVals = 0, const int idxLhs = 0) const {
return Parent::getAttribute((NRHS+idxLhs)*NVALS+idxVals);
}
FReal* getPotentialsArray(const int idxVals = 0, const int idxLhs = 0){
return Parent::getRawData() + ((NRHS+idxLhs)*NVALS+idxVals)*Parent::getLeadingRawData();
}
const FReal* getPotentialsArray(const int idxVals = 0, const int idxLhs = 0) const {
return Parent::getRawData() + ((NRHS+idxLhs)*NVALS+idxVals)*Parent::getLeadingRawData();
}
FReal* getForcesX(const int idxVals = 0, const int idxLhs = 0){
return Parent::getAttribute((NRHS+NLHS+idxLhs)*NVALS+idxVals);
}
const FReal* getForcesX(const int idxVals = 0, const int idxLhs = 0) const {
return Parent::getAttribute((NRHS+NLHS+idxLhs)*NVALS+idxVals);
}
const FReal* getPhysicalValues(const int idxRhs = 0) const {
return Parent::getAttribute(0+idxRhs);
FReal* getForcesXArray(const int idxVals = 0, const int idxLhs = 0){
return Parent::getRawData() + ((NRHS+NLHS+idxLhs)*NVALS+idxVals)*Parent::getLeadingRawData();
}
FReal* getPotentials(const int idxLhs = 0){
return Parent::getAttribute(NRHS+idxLhs);
const FReal* getForcesXArray(const int idxVals = 0, const int idxLhs = 0) const {
return Parent::getRawData() + ((NRHS+NLHS+idxLhs)*NVALS+idxVals)*Parent::getLeadingRawData();
}
const FReal* getPotentials(const int idxLhs = 0) const {
return Parent::getAttribute(NRHS+idxLhs);
FReal* getForcesY(const int idxVals = 0, const int idxLhs = 0){
return Parent::getAttribute((NRHS+2*NLHS+idxLhs)*NVALS+idxVals);
}
FReal* getForcesX(const int idxLhs = 0){
return Parent::getAttribute(NRHS+NLHS+idxLhs);
const FReal* getForcesY(const int idxVals = 0, const int idxLhs = 0) const {
return Parent::getAttribute((NRHS+2*NLHS+idxLhs)*NVALS+idxVals);
}
const FReal* getForcesX(const int idxLhs = 0) const {
return Parent::getAttribute(NRHS+NLHS+idxLhs);
FReal* getForcesYArray(const int idxVals = 0, const int idxLhs = 0){
return Parent::getRawData() + ((NRHS+2*NLHS+idxLhs)*NVALS+idxVals)*Parent::getLeadingRawData();
}
FReal* getForcesY(const int idxLhs = 0){
return Parent::getAttribute(NRHS+2*NLHS+idxLhs);
const FReal* getForcesYArray(const int idxVals = 0, const int idxLhs = 0) const {
return Parent::getRawData() + ((NRHS+2*NLHS+idxLhs)*NVALS+idxVals)*Parent::getLeadingRawData();
}
const FReal* getForcesY(const int idxLhs = 0) const {
return Parent::getAttribute(NRHS+2*NLHS+idxLhs);
FReal* getForcesZ(const int idxVals = 0, const int idxLhs = 0){
return Parent::getAttribute((NRHS+3*NLHS+idxLhs)*NVALS+idxVals);
}
FReal* getForcesZ(const int idxLhs = 0){
return Parent::getAttribute(NRHS+3*NLHS+idxLhs);
const FReal* getForcesZ(const int idxVals = 0, const int idxLhs = 0) const {
return Parent::getAttribute((NRHS+3*NLHS+idxLhs)*NVALS+idxVals);
}
const FReal* getForcesZ(const int idxLhs = 0) const {
return Parent::getAttribute(NRHS+3*NLHS+idxLhs);
FReal* getForcesZArray(const int idxVals = 0, const int idxLhs = 0){
return Parent::getRawData() + ((NRHS+3*NLHS+idxLhs)*NVALS+idxVals)*Parent::getLeadingRawData();
}
const FReal* getForcesZArray(const int idxVals = 0, const int idxLhs = 0) const {
return Parent::getRawData() + ((NRHS+3*NLHS+idxLhs)*NVALS+idxVals)*Parent::getLeadingRawData();
}
void resetForcesAndPotential(){
for(int idx = 0 ; idx < 4*NLHS*NVALS ; ++idx){
Parent::resetToInitialState(idx + NRHS*NVALS);
}
}
int getNVALS() const {
return NVALS;
}
};
#endif // FP2PGROUPPARTICLECONTAINER_HPP
......@@ -45,7 +45,7 @@ class FAbstractChebKernel : public FAbstractKernels< CellClass, ContainerClass>
{
protected:
enum {nnodes = TensorTraits<ORDER>::nnodes};
typedef FChebInterpolator<ORDER,MatrixKernelClass> InterpolatorClass;
typedef FChebInterpolator<ORDER,MatrixKernelClass,NVALS> InterpolatorClass;
/// Needed for P2M, M2M, L2L and L2P operators
const FSmartPointer<InterpolatorClass,FSmartPointerMemory> Interpolator;
......
This diff is collapsed.
......@@ -77,19 +77,22 @@ public:
}
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,
const ContainerClass* const SourceParticles)
{
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){
// 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);
}
......@@ -102,7 +105,6 @@ public:
{
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),
......@@ -192,22 +194,23 @@ public:
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);
}
//// 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
......
......@@ -180,10 +180,8 @@ public:
{
// apply Sy
const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getMultipole(idxRhs), SourceParticles);
}
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getMultipole(0), SourceParticles);
}
......@@ -445,7 +443,6 @@ public:
ContainerClass* const TargetParticles)
{
const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// // a) apply Sx
// AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter,
// AbstractBaseClass::BoxWidthLeaf,
......@@ -457,10 +454,10 @@ public:
// LeafCell->getLocal(),
// TargetParticles);
// c) apply Sx and Px (grad Sx)
AbstractBaseClass::Interpolator->applyL2PTotal(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getLocal(idxRhs), TargetParticles);
}
// 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
......
......@@ -16,6 +16,7 @@
// ==== CMAKE =====
// @FUSE_FFT
// @FUSE_BLAS
// ================
// Keep in private GIT
// @SCALFMM_PRIVATE
......
......@@ -16,6 +16,7 @@
// ==== CMAKE =====
// @FUSE_BLAS
// @FUSE_FFT
// ==============
#include <array>
......@@ -37,6 +38,7 @@
#include "Kernels/Chebyshev/FChebCell.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Chebyshev/FChebSymKernel.hpp"
#include "Kernels/Chebyshev/FChebKernel.hpp"
#include "Kernels/Uniform/FUnifCell.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
......@@ -101,33 +103,20 @@ class TestInterpolationKernel : public FUTester<TestInterpolationKernel> {
// Create octree
OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
// Insert particle in the tree
// For each particles we associate Nvals charge ( q,0,0,0)
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
double q = particles[idxPart].getPhysicalValue();
if(NVals == 1){
tree.insert(particles[idxPart].getPosition() , idxPart, q);//,0.0,0.0,0.0);
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){
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){
tree.insert(particles[idxPart].getPosition() , idxPart, q, q);//,0.0,0.0,0.0);
}
else{
FAssertLF(0, "NVALS should be <= 2");
}
}
// 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);
// }
// put in tree
tree.insert(particles[idxPart].getPosition(), idxPart, physicalState);
}
// Run FMM
Print("Fmm...");
......@@ -273,7 +262,7 @@ class TestInterpolationKernel : public FUTester<TestInterpolationKernel> {
/** TestUnifKernel */
void TestUnifKernel(){
const int NVals = 2;
const int NVals = 3;
const unsigned int ORDER = 6 ;
// run test
typedef FInterpMatrixKernelR MatrixKernelClass;
......@@ -291,7 +280,7 @@ class TestInterpolationKernel : public FUTester<TestInterpolationKernel> {
/** TestChebSymKernel */
void TestChebSymKernel(){
const int NVals = 2;
const int NVals = 3;
const unsigned int ORDER = 6;
typedef FP2PParticleContainerIndexed<1,1,NVals> ContainerClass;
typedef FSimpleLeaf<ContainerClass> LeafClass;
......@@ -304,7 +293,20 @@ class TestInterpolationKernel : public FUTester<TestInterpolationKernel> {
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