Commit 731ae45d authored by BLANCHARD Pierre's avatar BLANCHARD Pierre

Introduced tensorial physical value, potential and force (through...

Introduced tensorial physical value, potential and force (through templatization of P2PParticleContainer with dimensions NRHS and NLHS, default NRHS=NLHS=1). ChebTensorial kernel does not support this yet (TODO apply same modif on ChebInterpolator).
parent 8d142ebf
This diff is collapsed.
...@@ -18,48 +18,49 @@ ...@@ -18,48 +18,49 @@
#include "../../Components/FBasicParticleContainer.hpp" #include "../../Components/FBasicParticleContainer.hpp"
class FP2PParticleContainer : public FBasicParticleContainer<5> { template<int NRHS = 1, int NLHS = 1>
typedef FBasicParticleContainer<5> Parent; class FP2PParticleContainer : public FBasicParticleContainer<NRHS+4*NLHS> {
typedef FBasicParticleContainer<NRHS+4*NLHS> Parent;
public: public:
FReal* getPhysicalValues(){ FReal* getPhysicalValues(const int idxRhs = 0){
return Parent::getAttribute<0>(); return Parent::getAttribute(0+idxRhs);
} }
const FReal* getPhysicalValues() const { const FReal* getPhysicalValues(const int idxRhs = 0) const {
return Parent::getAttribute<0>(); return Parent::getAttribute(0+idxRhs);
} }
FReal* getPotentials(){ FReal* getPotentials(const int idxLhs = 0){
return Parent::getAttribute<1>(); return Parent::getAttribute(NRHS+idxLhs);
} }
const FReal* getPotentials() const { const FReal* getPotentials(const int idxLhs = 0) const {
return Parent::getAttribute<1>(); return Parent::getAttribute(NRHS+idxLhs);
} }
FReal* getForcesX(){ FReal* getForcesX(const int idxLhs = 0){
return Parent::getAttribute<2>(); return Parent::getAttribute(NRHS+NLHS+idxLhs);
} }
const FReal* getForcesX() const { const FReal* getForcesX(const int idxLhs = 0) const {
return Parent::getAttribute<2>(); return Parent::getAttribute(NRHS+NLHS+idxLhs);
} }
FReal* getForcesY(){ FReal* getForcesY(const int idxLhs = 0){
return Parent::getAttribute<3>(); return Parent::getAttribute(NRHS+2*NLHS+idxLhs);
} }
const FReal* getForcesY() const { const FReal* getForcesY(const int idxLhs = 0) const {
return Parent::getAttribute<3>(); return Parent::getAttribute(NRHS+2*NLHS+idxLhs);
} }
FReal* getForcesZ(){ FReal* getForcesZ(const int idxLhs = 0){
return Parent::getAttribute<4>(); return Parent::getAttribute(NRHS+3*NLHS+idxLhs);
} }
const FReal* getForcesZ() const { const FReal* getForcesZ(const int idxLhs = 0) const {
return Parent::getAttribute<4>(); return Parent::getAttribute(NRHS+3*NLHS+idxLhs);
} }
}; };
......
...@@ -6,8 +6,9 @@ ...@@ -6,8 +6,9 @@
#include "FP2PParticleContainer.hpp" #include "FP2PParticleContainer.hpp"
#include "../../Components/FParticleType.hpp" #include "../../Components/FParticleType.hpp"
class FP2PParticleContainerIndexed : public FP2PParticleContainer { template<int NRHS = 1, int NLHS = 1>
typedef FP2PParticleContainer Parent; class FP2PParticleContainerIndexed : public FP2PParticleContainer<NRHS,NLHS> {
typedef FP2PParticleContainer<NRHS,NLHS> Parent;
FVector<int> indexes; FVector<int> indexes;
......
...@@ -362,14 +362,13 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyP2M(const FPoint& c ...@@ -362,14 +362,13 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyP2M(const FPoint& c
const ContainerClass *const inParticles) const const ContainerClass *const inParticles) const
{ {
// set all multipole expansions to zero // set all multipole expansions to zero
FBlas::setzero(/*nRhs**/nnodes, multipoleExpansion); FBlas::setzero(nRhs*nnodes, multipoleExpansion);
// allocate stuff // allocate stuff
const map_glob_loc map(center, width); const map_glob_loc map(center, width);
FPoint localPosition; FPoint localPosition;
// loop over source particles // loop over source particles
const FReal*const physicalValues = inParticles->getPhysicalValues(); // PB: TODO add multidim PhysVal
const FReal*const positionsX = inParticles->getPositions()[0]; const FReal*const positionsX = inParticles->getPositions()[0];
const FReal*const positionsY = inParticles->getPositions()[1]; const FReal*const positionsY = inParticles->getPositions()[1];
const FReal*const positionsZ = inParticles->getPositions()[2]; const FReal*const positionsZ = inParticles->getPositions()[2];
...@@ -385,24 +384,24 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyP2M(const FPoint& c ...@@ -385,24 +384,24 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyP2M(const FPoint& c
L_of_x[o][2] = BasisType::L(o, localPosition.getZ()); // 3 * ORDER*(ORDER-1) flops L_of_x[o][2] = BasisType::L(o, localPosition.getZ()); // 3 * ORDER*(ORDER-1) flops
} }
// // PB: More optimal version where the sum over Lhs/Rhs is done inside applyL2P/P2M for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){
// // this avoids nLhs/nRhs evaluations of the same interpolating polynomials
// for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){
// read physicalValue // read physicalValue
const FReal weight = physicalValues[idxPart/*+idxRhs*nParticles*/]; //PB: TODO select next compo const FReal*const physicalValues = inParticles->getPhysicalValues(idxRhs);
// compute weight
const FReal weight = physicalValues[idxPart];
// assemble multipole expansions // assemble multipole expansions
for (unsigned int i=0; i<ORDER; ++i) { for (unsigned int i=0; i<ORDER; ++i) {
for (unsigned int j=0; j<ORDER; ++j) { for (unsigned int j=0; j<ORDER; ++j) {
for (unsigned int k=0; k<ORDER; ++k) { for (unsigned int k=0; k<ORDER; ++k) {
const unsigned int idx = /*idxRhs*nnodes +*/ k*ORDER*ORDER + j*ORDER + i; const unsigned int idx = idxRhs*nnodes + k*ORDER*ORDER + j*ORDER + i;
multipoleExpansion[idx] += L_of_x[i][0] * L_of_x[j][1] * L_of_x[k][2] * weight; // 3 * ORDER*ORDER*ORDER flops multipoleExpansion[idx] += L_of_x[i][0] * L_of_x[j][1] * L_of_x[k][2] * weight; // 3 * ORDER*ORDER*ORDER flops
} }
} }
} }
// } // idxRhs } // idxRhs
} // flops: N * (3 * ORDER*ORDER*ORDER + 3 * 3 * ORDER*(ORDER-1)) flops } // flops: N * (3 * ORDER*ORDER*ORDER + 3 * 3 * ORDER*(ORDER-1)) flops
...@@ -423,11 +422,9 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c ...@@ -423,11 +422,9 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c
const map_glob_loc map(center, width); const map_glob_loc map(center, width);
FPoint localPosition; FPoint localPosition;
//const FReal*const physicalValues = inParticles->getPhysicalValues();
const FReal*const positionsX = inParticles->getPositions()[0]; const FReal*const positionsX = inParticles->getPositions()[0];
const FReal*const positionsY = inParticles->getPositions()[1]; const FReal*const positionsY = inParticles->getPositions()[1];
const FReal*const positionsZ = inParticles->getPositions()[2]; const FReal*const positionsZ = inParticles->getPositions()[2];
FReal*const potentials = inParticles->getPotentials();
const unsigned int nParticles = inParticles->getNbParticles(); const unsigned int nParticles = inParticles->getNbParticles();
...@@ -444,9 +441,7 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c ...@@ -444,9 +441,7 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c
L_of_x[o][2] = BasisType::L(o, localPosition.getZ()); // 3 * ORDER*(ORDER-1) flops L_of_x[o][2] = BasisType::L(o, localPosition.getZ()); // 3 * ORDER*(ORDER-1) flops
} }
// // PB: More optimal version where the sum over Lhs/Rhs is done inside applyL2P/P2M for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
// // this avoid nLhs/nRhs evaluations of the same interpolating polynomials
// for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
// interpolate and increment target value // interpolate and increment target value
FReal targetValue=0.; FReal targetValue=0.;
...@@ -454,7 +449,7 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c ...@@ -454,7 +449,7 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c
for (unsigned int l=0; l<ORDER; ++l) { for (unsigned int l=0; l<ORDER; ++l) {
for (unsigned int m=0; m<ORDER; ++m) { for (unsigned int m=0; m<ORDER; ++m) {
for (unsigned int n=0; n<ORDER; ++n) { for (unsigned int n=0; n<ORDER; ++n) {
const unsigned int idx = /*idxLhs*nnodes +*/ n*ORDER*ORDER + m*ORDER + l; const unsigned int idx = idxLhs*nnodes + n*ORDER*ORDER + m*ORDER + l;
targetValue += targetValue +=
L_of_x[l][0] * L_of_x[m][1] * L_of_x[n][2] * localExpansion[idx]; L_of_x[l][0] * L_of_x[m][1] * L_of_x[n][2] * localExpansion[idx];
} // ORDER * 4 flops } // ORDER * 4 flops
...@@ -462,10 +457,12 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c ...@@ -462,10 +457,12 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c
} // ORDER * ORDER * ORDER * 4 flops } // ORDER * ORDER * ORDER * 4 flops
} }
// set potential // get potential
potentials[idxPart/*+idxLhs*nParticles*/] += (targetValue); FReal*const potentials = inParticles->getPotentials(idxLhs);
// add contribution to potential
potentials[idxPart] += (targetValue);
// } // idxLhs } // idxLhs
} // N * (4 * ORDER * ORDER * ORDER + 9 * ORDER*(ORDER-1) ) flops } // N * (4 * ORDER * ORDER * ORDER + 9 * ORDER*(ORDER-1) ) flops
} }
...@@ -495,14 +492,9 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2PGradient(const F ...@@ -495,14 +492,9 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2PGradient(const F
FReal L_of_x[ORDER][3]; FReal L_of_x[ORDER][3];
FReal dL_of_x[ORDER][3]; FReal dL_of_x[ORDER][3];
const FReal*const physicalValues = inParticles->getPhysicalValues();
const FReal*const positionsX = inParticles->getPositions()[0]; const FReal*const positionsX = inParticles->getPositions()[0];
const FReal*const positionsY = inParticles->getPositions()[1]; const FReal*const positionsY = inParticles->getPositions()[1];
const FReal*const positionsZ = inParticles->getPositions()[2]; const FReal*const positionsZ = inParticles->getPositions()[2];
FReal*const forcesX = inParticles->getForcesX();
FReal*const forcesY = inParticles->getForcesY();
FReal*const forcesZ = inParticles->getForcesZ();
//FReal*const potentials = inParticles->getPotentials();
// const unsigned int nParticles = inParticles->getNbParticles(); // const unsigned int nParticles = inParticles->getNbParticles();
...@@ -521,9 +513,7 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2PGradient(const F ...@@ -521,9 +513,7 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2PGradient(const F
dL_of_x[o][2] = BasisType::dL(o, localPosition.getZ()); // TODO verify 3 * ORDER*(ORDER-1) flops dL_of_x[o][2] = BasisType::dL(o, localPosition.getZ()); // TODO verify 3 * ORDER*(ORDER-1) flops
} }
// // PB: More optimal version where the sum over Lhs/Rhs is done inside applyL2P/P2M for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
// // this avoid nLhs/nLhs evaluations of the same interpolating polynomials
// for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
// interpolate and increment forces value // interpolate and increment forces value
FReal forces[3] = {FReal(0.), FReal(0.), FReal(0.)}; FReal forces[3] = {FReal(0.), FReal(0.), FReal(0.)};
...@@ -531,7 +521,7 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2PGradient(const F ...@@ -531,7 +521,7 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2PGradient(const F
for (unsigned int l=0; l<ORDER; ++l) { for (unsigned int l=0; l<ORDER; ++l) {
for (unsigned int m=0; m<ORDER; ++m) { for (unsigned int m=0; m<ORDER; ++m) {
for (unsigned int n=0; n<ORDER; ++n) { for (unsigned int n=0; n<ORDER; ++n) {
const unsigned int idx = /*idxLhs*nnodes +*/ n*ORDER*ORDER + m*ORDER + l; const unsigned int idx = idxLhs*nnodes + n*ORDER*ORDER + m*ORDER + l;
forces[0] += forces[0] +=
dL_of_x[l][0] * L_of_x[m][1] * L_of_x[n][2] * localExpansion[idx]; dL_of_x[l][0] * L_of_x[m][1] * L_of_x[n][2] * localExpansion[idx];
forces[1] += forces[1] +=
...@@ -545,19 +535,23 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2PGradient(const F ...@@ -545,19 +535,23 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2PGradient(const F
// scale forces // scale forces
forces[0] *= jacobian[0];// / nnodes; forces[0] *= jacobian[0];
forces[1] *= jacobian[1];// / nnodes; forces[1] *= jacobian[1];
forces[2] *= jacobian[2];// / nnodes; forces[2] *= jacobian[2];
} }
const FReal*const physicalValues = inParticles->getPhysicalValues(idxLhs);
FReal*const forcesX = inParticles->getForcesX(idxLhs);
FReal*const forcesY = inParticles->getForcesY(idxLhs);
FReal*const forcesZ = inParticles->getForcesZ(idxLhs);
// set computed forces // set computed forces
// const unsigned int idx = idxPart+idxLhs*nParticles forcesX[idxPart] += forces[0] * physicalValues[idxPart];
forcesX[idxPart] += forces[0] * physicalValues[idxPart]; //PB: TODO update next compo forcesY[idxPart] += forces[1] * physicalValues[idxPart];
forcesY[idxPart] += forces[1] * physicalValues[idxPart]; //PB: TODO update next compo forcesZ[idxPart] += forces[2] * physicalValues[idxPart];
forcesZ[idxPart] += forces[2] * physicalValues[idxPart]; //PB: TODO update next compo }// idxLhs
}
// }// idxLhs }
} }
......
...@@ -41,11 +41,8 @@ class FTreeCoordinate; ...@@ -41,11 +41,8 @@ class FTreeCoordinate;
* 1) Handling tensorial kernels (DIM,NRHS,NLHS) and having multiple rhs (NVALS) * 1) Handling tensorial kernels (DIM,NRHS,NLHS) and having multiple rhs (NVALS)
* are considered 2 separate features and are currently combined. * are considered 2 separate features and are currently combined.
* *
* 2) The present tensorial version is the most naive one. All tensorial aspects * 2) When it comes to applying M2L it is NOT much faster to loop over NRHSxNLHS
* are handled in the kernel. A more optimal version would be to consider looping * inside applyM2L (at least for the Lagrange case).
* over nRhs/nLhs inside Interpolator::applyP2M/L2P in order to avoid extra
* evaluation of the interpolating polynomials. When it comes to applying M2L it is
* NOT much faster to loop over NRHSxNLHS inside applyM2L (at least for the Lagrange case)
* 2-bis) The evaluation of the kernel matrix (see M2LHandler) should be done at once * 2-bis) The evaluation of the kernel matrix (see M2LHandler) should be done at once
* instead of compo-by-compo (TODO). On the other hand, the ChebyshevSym tensorial kernel * instead of compo-by-compo (TODO). On the other hand, the ChebyshevSym tensorial kernel
* requires the matrix kernel to be evaluated compo-by-compo since we currently use a scalar ACA. * requires the matrix kernel to be evaluated compo-by-compo since we currently use a scalar ACA.
...@@ -102,15 +99,17 @@ public: ...@@ -102,15 +99,17 @@ public:
const ContainerClass* const SourceParticles) const ContainerClass* const SourceParticles)
{ {
const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate())); const FPoint LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
for(int idxV = 0 ; idxV < NVALS ; ++idxV){ for(int idxV = 0 ; idxV < NVALS ; ++idxV){
// 1) apply Sy
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->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;
// 1) apply Sy
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getMultipole(idxMul), SourceParticles);
// 2) apply Discrete Fourier Transform // 2) apply Discrete Fourier Transform
M2LHandler->applyZeroPaddingAndDFT(LeafCell->getMultipole(idxMul), M2LHandler->applyZeroPaddingAndDFT(LeafCell->getMultipole(idxMul),
LeafCell->getTransformedMultipole(idxMul)); LeafCell->getTransformedMultipole(idxMul));
...@@ -211,15 +210,16 @@ public: ...@@ -211,15 +210,16 @@ public:
M2LHandler->unapplyZeroPaddingAndDFT(LeafCell->getTransformedLocal(idxLoc), M2LHandler->unapplyZeroPaddingAndDFT(LeafCell->getTransformedLocal(idxLoc),
const_cast<CellClass*>(LeafCell)->getLocal(idxLoc)); const_cast<CellClass*>(LeafCell)->getLocal(idxLoc));
// 2.a) apply Sx }
AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getLocal(idxLoc), TargetParticles); // 2.a) apply Sx
AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getLocal(idxV*nLhs), TargetParticles);
// 2.b) apply Px (grad Sx) // 2.b) apply Px (grad Sx)
AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf, AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getLocal(idxLoc), TargetParticles); LeafCell->getLocal(idxV*nLhs), TargetParticles);
}
}// NVALS }// NVALS
} }
......
...@@ -124,7 +124,7 @@ int main(int argc, char* argv[]) ...@@ -124,7 +124,7 @@ int main(int argc, char* argv[])
const unsigned int ORDER = 7; const unsigned int ORDER = 7;
const FReal epsilon = FReal(1e-7); const FReal epsilon = FReal(1e-7);
// typedefs // typedefs
typedef FP2PParticleContainerIndexed ContainerClass; typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass; typedef FSimpleLeaf< ContainerClass > LeafClass;
//typedef FInterpMatrixKernelLJ MatrixKernelClass; //typedef FInterpMatrixKernelLJ MatrixKernelClass;
typedef FInterpMatrixKernelR MatrixKernelClass; typedef FInterpMatrixKernelR MatrixKernelClass;
......
...@@ -58,7 +58,7 @@ int main(int argc, char* argv[]) ...@@ -58,7 +58,7 @@ int main(int argc, char* argv[])
const unsigned int ORDER = 7; const unsigned int ORDER = 7;
const FReal epsilon = FReal(1e-7); const FReal epsilon = FReal(1e-7);
typedef FP2PParticleContainerIndexed ContainerClass; typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass; typedef FSimpleLeaf< ContainerClass > LeafClass;
typedef FInterpMatrixKernelR MatrixKernelClass; typedef FInterpMatrixKernelR MatrixKernelClass;
......
...@@ -66,12 +66,23 @@ int main(int argc, char* argv[]) ...@@ -66,12 +66,23 @@ int main(int argc, char* argv[])
// init timer // init timer
FTic time; FTic time;
// typedefs
// typedef FInterpMatrixKernelLJ MatrixKernelClass;
// typedef FInterpMatrixKernelR MatrixKernelClass;
// typedef FInterpMatrixKernel_R_IJ MatrixKernelClass; // not working with Non-Symmetric variant ! Because of UCB decomposition.
// and not working either with Symmetric variant because not symmetric...
typedef FInterpMatrixKernel_IOR MatrixKernelClass;
// const KERNEL_FUNCTION_IDENTIFIER MK_ID = MatrixKernelClass::Identifier;
const unsigned int NRHS = MatrixKernelClass::NRHS;
const unsigned int NLHS = MatrixKernelClass::NLHS;
// init particles position and physical value // init particles position and physical value
struct TestParticle{ struct TestParticle{
FPoint position; FPoint position;
FReal forces[3]; FReal forces[3][NLHS];
FReal physicalValue; FReal physicalValue[NRHS];
FReal potential; FReal potential[NLHS];
}; };
// open particle file // open particle file
...@@ -85,11 +96,14 @@ int main(int argc, char* argv[]) ...@@ -85,11 +96,14 @@ int main(int argc, char* argv[])
loader.fillParticle(&position,&physicalValue); loader.fillParticle(&position,&physicalValue);
// get copy // get copy
particles[idxPart].position = position; particles[idxPart].position = position;
particles[idxPart].physicalValue = physicalValue; for(unsigned idxRhs = 0; idxRhs<NRHS;++idxRhs)
particles[idxPart].potential = 0.0; particles[idxPart].physicalValue[idxRhs] = physicalValue;
particles[idxPart].forces[0] = 0.0; for(unsigned idxLhs = 0; idxLhs<NLHS;++idxLhs){
particles[idxPart].forces[1] = 0.0; particles[idxPart].potential[idxLhs] = 0.0;
particles[idxPart].forces[2] = 0.0; particles[idxPart].forces[0][idxLhs] = 0.0;
particles[idxPart].forces[1][idxLhs] = 0.0;
particles[idxPart].forces[2][idxLhs] = 0.0;
}
} }
//////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////
...@@ -110,12 +124,12 @@ int main(int argc, char* argv[]) ...@@ -110,12 +124,12 @@ int main(int argc, char* argv[])
// &particles[idxOther].forces[2], &particles[idxOther].potential); // &particles[idxOther].forces[2], &particles[idxOther].potential);
FP2P::MutualParticlesIOR(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(), FP2P::MutualParticlesIOR(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue, particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue,
&particles[idxTarget].forces[0], &particles[idxTarget].forces[1], particles[idxTarget].forces[0], particles[idxTarget].forces[1],
&particles[idxTarget].forces[2], &particles[idxTarget].potential, particles[idxTarget].forces[2], particles[idxTarget].potential,
particles[idxOther].position.getX(), particles[idxOther].position.getY(), particles[idxOther].position.getX(), particles[idxOther].position.getY(),
particles[idxOther].position.getZ(), particles[idxOther].physicalValue, particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
&particles[idxOther].forces[0], &particles[idxOther].forces[1], particles[idxOther].forces[0], particles[idxOther].forces[1],
&particles[idxOther].forces[2], &particles[idxOther].potential); particles[idxOther].forces[2], particles[idxOther].potential);
} }
} }
} }
...@@ -134,18 +148,7 @@ int main(int argc, char* argv[]) ...@@ -134,18 +148,7 @@ int main(int argc, char* argv[])
const FReal epsilon = FReal(1e-8); const FReal epsilon = FReal(1e-8);
// typedefs // typedefs
//typedef FInterpMatrixKernelLJ MatrixKernelClass; typedef FP2PParticleContainerIndexed<NRHS,NLHS> ContainerClass;
// typedef FInterpMatrixKernelR MatrixKernelClass;
// typedef FInterpMatrixKernel_R_IJ MatrixKernelClass; // not working with Non-Symmetric variant ! Because of UCB decomposition.
// and not working either with Symmetric variant because of lack of symmetry...
typedef FInterpMatrixKernel_IOR MatrixKernelClass;
const unsigned int NRHS = MatrixKernelClass::NRHS;
const unsigned int NLHS = MatrixKernelClass::NLHS;
typedef FP2PParticleContainerIndexed ContainerClass;
// const unsigned int NDIM = NRHS + 4*NLHS;
// typedef FP2PTensorialParticleContainerIndexed<NDIM> ContainerClass; // TODO fix a TensorialParticleContainer for easy access to multidim PhysVal
typedef FSimpleLeaf< ContainerClass > LeafClass; typedef FSimpleLeaf< ContainerClass > LeafClass;
typedef FChebCell<ORDER,NRHS,NLHS> CellClass; typedef FChebCell<ORDER,NRHS,NLHS> CellClass;
...@@ -166,7 +169,8 @@ int main(int argc, char* argv[]) ...@@ -166,7 +169,8 @@ int main(int argc, char* argv[])
for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){ for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
// put in tree // put in tree
tree.insert(particles[idxPart].position, idxPart, particles[idxPart].physicalValue); // PB: here we have to know the number of NRHS...
tree.insert(particles[idxPart].position, idxPart, particles[idxPart].physicalValue[0], particles[idxPart].physicalValue[1], particles[idxPart].physicalValue[2]);
} }
time.tac(); time.tac();
...@@ -187,44 +191,68 @@ int main(int argc, char* argv[]) ...@@ -187,44 +191,68 @@ int main(int argc, char* argv[])
{ // ----------------------------------------------------- { // -----------------------------------------------------
std::cout << "\nError computation ... " << std::endl; std::cout << "\nError computation ... " << std::endl;
FMath::FAccurater potentialDiff; FMath::FAccurater potentialDiff[NLHS];
FMath::FAccurater fx, fy, fz; FMath::FAccurater fx[NLHS], fy[NLHS], fz[NLHS];
FReal checkPotential[20000]; FReal checkPhysVal[20000][NRHS];
FReal checkPotential[20000][NLHS];
FReal checkfx[20000][NLHS];
{ // Check that each particle has been summed with all other { // Check that each particle has been summed with all other
tree.forEachLeaf([&](LeafClass* leaf){ tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const potentials = leaf->getTargets()->getPotentials(); for(unsigned idxLhs = 0; idxLhs<NLHS;++idxLhs){
const FReal*const forcesX = leaf->getTargets()->getForcesX(); const FReal*const physVals = leaf->getTargets()->getPhysicalValues(idxLhs);
const FReal*const forcesY = leaf->getTargets()->getForcesY(); const FReal*const potentials = leaf->getTargets()->getPotentials(idxLhs);