Commit 7a0d325c authored by BLANCHARD Pierre's avatar BLANCHARD Pierre

Implemented a tensorial version for ChebyshevSym kernel (currently in separate...

Implemented a tensorial version for ChebyshevSym kernel (currently in separate files for clarity/safety, TODO finalize and merge), cleaned up tensorial Uniform kernel.
parent 25e7da9d
......@@ -46,7 +46,8 @@ SET(my_include_dirs "."
"Arranger" "Core" "Utils"
"Kernels/Chebyshev" "Components" "Extensions"
"Containers" "Files" "Kernels/Spherical"
"Kernels/P2P" "Kernels/Taylor" "Kernels/ROtation")
"Kernels/P2P" "Kernels/Taylor" "Kernels/ROtation"
"Kernels/Uniform" "Kernels/Interpolation")
FOREACH(my_dir ${my_include_dirs})
file(GLOB
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -973,7 +973,6 @@ public:
(*sourcePotential) += ( potentialCoef * targetPhysicalValue );
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////
......@@ -1033,6 +1032,7 @@ public:
// potentials
FReal potentialCoef;
// potentialCoef = r[0]*r[0] * inv_distance_pow3; //PB: test
if(i==j)
potentialCoef = inv_distance - ri2 * inv_distance_pow3;
else
......@@ -1101,6 +1101,7 @@ public:
// potentials
FReal potentialCoef;
// potentialCoef = r[0]*r[0] * inv_distance_pow3; //PB: test
if(i==j)
potentialCoef = inv_distance - ri2 * inv_distance_pow3;
else
......@@ -1195,6 +1196,7 @@ public:
// potentials
FReal potentialCoef;
// potentialCoef = r[0]*r[0] * inv_distance_pow3; //PB: test
if(i==j)
potentialCoef = inv_distance - ri2 * inv_distance_pow3;
else
......@@ -1283,6 +1285,7 @@ public:
// potentials
FReal potentialCoef;
// potentialCoef = r[0]*r[0] * inv_distance_pow3; //PB: test
if(i==j)
potentialCoef = inv_distance - ri2 * inv_distance_pow3;
else
......@@ -1388,20 +1391,20 @@ public:
// potentials
FReal potentialCoef;
if(i==j)
//if(i==j)
potentialCoef = inv_distance;
else
potentialCoef = FReal(0.);
//else
// potentialCoef = FReal(0.);
// forces
FReal coef[3];
for(unsigned int k = 0 ; k < 3 ; ++k){
if(i==j){
//if(i==j){
coef[k] = + r[k] * inv_distance_pow3 * (targetsPhysicalValues[idxTarget/* + i*nbParticles*/] * sourcesPhysicalValues[idxSource/* + j*nbParticles*/]);
}
else{
coef[k] = FReal(0.);
}
//}
//else{
// coef[k] = FReal(0.);
//}
}
......@@ -1442,20 +1445,20 @@ public:
// potentials
FReal potentialCoef;
if(i==j)
//if(i==j)
potentialCoef = inv_distance;
else
potentialCoef = FReal(0.);
//else
// potentialCoef = FReal(0.);
// forces
FReal coef[3];
for(unsigned int k = 0 ; k < 3 ; ++k){
if(i==j){
//if(i==j){
coef[k] = + r[k] * inv_distance_pow3 * (targetsPhysicalValues[idxTarget/* + i*nbParticles*/] * targetsPhysicalValues[idxSource/* + j*nbParticles*/]);
}
else{
coef[k] = FReal(0.);
}
//}
//else{
// coef[k] = FReal(0.);
//}
}
......@@ -1523,20 +1526,20 @@ public:
// potentials
FReal potentialCoef;
if(i==j)
//if(i==j)
potentialCoef = inv_distance ;
else
potentialCoef = FReal(0.);
//else
// potentialCoef = FReal(0.);
// forces
FReal coef[3];
for(unsigned int k = 0 ; k < 3 ; ++k){
if(i==j){
//if(i==j){
coef[k] = + r[k] * inv_distance_pow3 * (targetsPhysicalValues[idxTarget/* + i*nbParticles*/] * sourcesPhysicalValues[idxSource/* + j*nbParticles*/]);
}
else{
coef[k] = FReal(0.);
}
//}
//else{
// coef[k] = FReal(0.);
//}
}
......@@ -1598,20 +1601,20 @@ public:
// potentials
FReal potentialCoef;
if(i==j)
//if(i==j)
potentialCoef = inv_distance;
else
potentialCoef = FReal(0.);
//else
// potentialCoef = FReal(0.);
// forces
FReal coef[3];
for(unsigned int k = 0 ; k < 3 ; ++k){
if(i==j){
//if(i==j){
coef[k] = + r[k] * inv_distance_pow3 * (targetPhysicalValue * sourcePhysicalValue);
}
else{
coef[k] = FReal(0.);
}
//}
//else{
// coef[k] = FReal(0.);
//}
}
......
......@@ -362,7 +362,7 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyP2M(const FPoint& c
const ContainerClass *const inParticles) const
{
// set all multipole expansions to zero
FBlas::setzero(nRhs*nnodes, multipoleExpansion);
FBlas::setzero(/*nRhs**/nnodes, multipoleExpansion);
// allocate stuff
const map_glob_loc map(center, width);
......@@ -386,7 +386,7 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyP2M(const FPoint& c
}
// // PB: More optimal version where the sum over Lhs/Rhs is done inside applyL2P/P2M
// // this avoid nLhs/nLhs evaluations of the same interpolating polynomials
// // this avoids nLhs/nRhs evaluations of the same interpolating polynomials
// for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs){
// read physicalValue
......@@ -406,11 +406,6 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyP2M(const FPoint& c
} // flops: N * (3 * ORDER*ORDER*ORDER + 3 * 3 * ORDER*(ORDER-1)) flops
// std::cout<< "multipoleExpansion="<<std::endl;
// for(int idxRhs = 0 ; idxRhs < nRhs ; ++idxRhs)
// for (unsigned int i=0; i<nnodes; ++i)
// std::cout<< multipoleExpansion[idxRhs*nnodes+i] << ", ";
// std::cout<<std::endl;
}
......@@ -428,12 +423,6 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c
const map_glob_loc map(center, width);
FPoint localPosition;
// std::cout<< "localExpansion="<<std::endl;
// for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs)
// for (unsigned int i=0; i<nnodes; ++i)
// std::cout<< localExpansion[idxLhs*nnodes+i] << ", ";
// std::cout<<std::endl;
//const FReal*const physicalValues = inParticles->getPhysicalValues();
const FReal*const positionsX = inParticles->getPositions()[0];
const FReal*const positionsY = inParticles->getPositions()[1];
......@@ -450,13 +439,13 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c
// evaluate Lagrange polynomial at local position
FReal L_of_x[ORDER][3];
for (unsigned int o=0; o<ORDER; ++o) {
L_of_x[o][0] = BasisType::L(o, localPosition.getX()); // 3 * ORDER*(ORDER-1) flops PB: TODO confirm
L_of_x[o][0] = BasisType::L(o, localPosition.getX()); // 3 * ORDER*(ORDER-1) flops
L_of_x[o][1] = BasisType::L(o, localPosition.getY()); // 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
// // this avoid nLhs/nLhs evaluations of the same interpolating polynomials
// // this avoid nLhs/nRhs evaluations of the same interpolating polynomials
// for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
// interpolate and increment target value
......@@ -474,7 +463,7 @@ inline void FUnifInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c
}
// set potential
potentials[idxPart/*+idxLhs*nParticles*/] += (targetValue); //PB: TODO update next compo
potentials[idxPart/*+idxLhs*nParticles*/] += (targetValue);
// } // idxLhs
......
......@@ -56,11 +56,11 @@ class FUnifTensorialM2LHandler : FNoCopyable
nnodes = TensorTraits<ORDER>::nnodes,
ninteractions = 316, // 7^3 - 3^3 (max num cells in far-field)
rc = (2*ORDER-1)*(2*ORDER-1)*(2*ORDER-1),
dim = MatrixKernelClass::Dim,
dim = MatrixKernelClass::DIM,
nRhs = MatrixKernelClass::NRHS,
nLhs = MatrixKernelClass::NLHS};
/*const*/ MatrixKernelClass MatrixKernel;
const MatrixKernelClass MatrixKernel;
FComplexe *FC[dim];
......@@ -86,7 +86,7 @@ class FUnifTensorialM2LHandler : FNoCopyable
public:
FUnifTensorialM2LHandler()
: MatrixKernel(), /*FC(NULL),*/
: MatrixKernel(), // PB: used only to getPosition and getScalFac in applyM2L
opt_rc(rc/2+1),
Dft(rc) // initialize Discrete Fourier Transformator
{
......@@ -273,10 +273,7 @@ FUnifTensorialM2LHandler<ORDER, MatrixKernelClass>::Compute(FComplexe* &FC,
// set roots of target cell (X)
FUnifTensor<order>::setRoots(FPoint(0.,0.,0.), FReal(2.), X);
// init matrix kernel
/*const*/ MatrixKernelClass MatrixKernel;
// update index of compo to be computed
MatrixKernel.updateIndex(d);
const MatrixKernelClass MatrixKernel(d);
// allocate memory and compute 316 m2l operators
FReal *_C;
......
......@@ -30,7 +30,7 @@
#include "../../Src/Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "../../Src/Kernels/Chebyshev/FChebCell.hpp"
#include "../../Src/Kernels/Chebyshev/FChebTensorialKernel.hpp"
#include "../../Src/Kernels/Chebyshev/FChebSymTensorialKernel.hpp"
#include "../../Src/Components/FSimpleLeaf.hpp"
#include "../../Src/Kernels/P2P/FP2PParticleContainerIndexed.hpp"
......@@ -100,14 +100,22 @@ int main(int argc, char* argv[])
{
for(int idxTarget = 0 ; idxTarget < loader.getNumberOfParticles() ; ++idxTarget){
for(int idxOther = idxTarget + 1 ; idxOther < loader.getNumberOfParticles() ; ++idxOther){
FP2P::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue,
&particles[idxTarget].forces[0], &particles[idxTarget].forces[1],
&particles[idxTarget].forces[2], &particles[idxTarget].potential,
particles[idxOther].position.getX(), particles[idxOther].position.getY(),
particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
&particles[idxOther].forces[0], &particles[idxOther].forces[1],
&particles[idxOther].forces[2], &particles[idxOther].potential);
// FP2P::MutualParticles(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
// particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue,
// &particles[idxTarget].forces[0], &particles[idxTarget].forces[1],
// &particles[idxTarget].forces[2], &particles[idxTarget].potential,
// particles[idxOther].position.getX(), particles[idxOther].position.getY(),
// particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
// &particles[idxOther].forces[0], &particles[idxOther].forces[1],
// &particles[idxOther].forces[2], &particles[idxOther].potential);
FP2P::MutualParticlesIOR(particles[idxTarget].position.getX(), particles[idxTarget].position.getY(),
particles[idxTarget].position.getZ(), particles[idxTarget].physicalValue,
&particles[idxTarget].forces[0], &particles[idxTarget].forces[1],
&particles[idxTarget].forces[2], &particles[idxTarget].potential,
particles[idxOther].position.getX(), particles[idxOther].position.getY(),
particles[idxOther].position.getZ(), particles[idxOther].physicalValue,
&particles[idxOther].forces[0], &particles[idxOther].forces[1],
&particles[idxOther].forces[2], &particles[idxOther].potential);
}
}
}
......@@ -122,18 +130,18 @@ int main(int argc, char* argv[])
{ // begin Chebyshev kernel
// accuracy
const unsigned int ORDER = 3;
const FReal epsilon = FReal(1e-7);
const unsigned int ORDER = 4;
const FReal epsilon = FReal(1e-8);
// typedefs
//typedef FInterpMatrixKernelLJ MatrixKernelClass;
typedef FInterpMatrixKernelR 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 of lack of symmetry...
typedef FInterpMatrixKernel_IOR MatrixKernelClass;
const unsigned int NRHS = MatrixKernelClass::NRHS;
const unsigned int NLHS = MatrixKernelClass::NLHS;
std::cout << "NRHS=" << NRHS << std::endl;
std::cout << "NLHS=" << NLHS << std::endl;
typedef FP2PParticleContainerIndexed ContainerClass;
// const unsigned int NDIM = NRHS + 4*NLHS;
......@@ -142,8 +150,7 @@ int main(int argc, char* argv[])
typedef FSimpleLeaf< ContainerClass > LeafClass;
typedef FChebCell<ORDER,NRHS,NLHS> CellClass;
typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
typedef FChebTensorialKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass; // PB: TODO erase when symmetric version done
//typedef FChebTensorialSymKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass; // PB: TODO
typedef FChebSymTensorialKernel<CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass; // PB: tensorial kernel needs to be symmetric !!
typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
// typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
......@@ -168,7 +175,7 @@ int main(int argc, char* argv[])
} // -----------------------------------------------------
{ // -----------------------------------------------------
std::cout << "\nLagrange/Uniform grid FMM (ORDER="<< ORDER << ") ... " << std::endl;
std::cout << "\nChebyshev FMM (ORDER="<< ORDER << ") ... " << std::endl;
time.tic();
KernelClass kernels(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(), epsilon);
FmmClass algorithm(&tree, &kernels);
......
......@@ -148,8 +148,6 @@ int main(int argc, char* argv[])
const unsigned int NRHS = MatrixKernelClass::NRHS;
const unsigned int NLHS = MatrixKernelClass::NLHS;
std::cout << "NRHS=" << NRHS << std::endl;
std::cout << "NLHS=" << NLHS << std::endl;
typedef FP2PParticleContainerIndexed ContainerClass;
// const unsigned int NDIM = NRHS + 4*NLHS;
......
......@@ -68,8 +68,7 @@ int main(int, char **){
std::cout << " direct computation.\n" << std::endl;
//////////////////////////////////////////////////////////////
MatrixKernelClass MatrixKernel;
RIJKMatrixKernelClass RIJKMatrixKernel; // Only used for kernel R_IJ
const unsigned int dim = MatrixKernelClass::DIM;
const unsigned int nrhs = MatrixKernelClass::NRHS;
const unsigned int nlhs = MatrixKernelClass::NLHS;
......@@ -123,6 +122,8 @@ int main(int, char **){
const unsigned int nnodes = TensorTraits<ORDER>::nnodes;
typedef FUnifInterpolator<ORDER,MatrixKernelClass> InterpolatorClass;
InterpolatorClass S;
MatrixKernelClass MatrixKernel;
RIJKMatrixKernelClass RIJKMatrixKernel;
std::cout << "\nCompute interactions approximatively, interpolation order = " << ORDER << " ..." << std::endl;
......@@ -153,9 +154,7 @@ int main(int, char **){
unsigned idxK = idxLhs*3+idxRhs; // or counter
unsigned int d = MatrixKernel.getPosition(idxK);
MatrixKernel.updateIndex(d);
F[i+idxLhs*nnodes] += MatrixKernel.evaluate(rootsX[i], rootsY[j]) * W[j+idxRhs*nnodes];
F[i+idxLhs*nnodes] += MatrixKernelClass(d).evaluate(rootsX[i], rootsY[j]) * W[j+idxRhs*nnodes];
}
}
......@@ -172,14 +171,12 @@ int main(int, char **){
////////////////////////////////////////////////////////////////////////////
// Store M2L in K and apply K
const unsigned int dim=MatrixKernelClass::Dim;
FReal K[dim*nnodes*nnodes]; // local expansion
for (unsigned int i=0; i<nnodes; ++i) {
for (unsigned int j=0; j<nnodes; ++j){
for (unsigned int d=0; d<dim; ++d){
MatrixKernel.updateIndex(d);
K[d*nnodes*nnodes+i*nnodes+j] = MatrixKernel.evaluate(rootsX[i], rootsY[j]);
K[d*nnodes*nnodes+i*nnodes+j] = MatrixKernelClass(d).evaluate(rootsX[i], rootsY[j]);
}
}
......@@ -231,10 +228,9 @@ int main(int, char **){
for(unsigned int n=0; n<2*ORDER-1; ++n){
for (unsigned int d=0; d<dim; ++d){
MatrixKernel.updateIndex(d);
C[d*rc + ido]
= MatrixKernel.evaluate(rootsX[node_ids_pairs[ido][0]],
= MatrixKernelClass(d).evaluate(rootsX[node_ids_pairs[ido][0]],
rootsY[node_ids_pairs[ido][1]]);
}
......@@ -560,8 +556,7 @@ int main(int, char **){
for (unsigned int i=0; i<nlhs; ++i) // sum all compo
for (unsigned int j=0; j<nrhs; ++j){
unsigned int d = MatrixKernel.getPosition(i*nrhs+j);
MatrixKernel.updateIndex(d);
const FReal rij = MatrixKernel.evaluate(x, y);
const FReal rij = MatrixKernelClass(d).evaluate(x, y);
// potential
p[counter] += rij * wy;
// force
......@@ -569,8 +564,7 @@ int main(int, char **){
for (unsigned int k=0; k<3; ++k){
//std::cout << "i,j,k,=" << i << ","<< j << ","<< k << std::endl;
unsigned int dk = RIJKMatrixKernel.getPosition(i*3*3+j*3+k);
RIJKMatrixKernel.updateIndex(dk);
force[k] = RIJKMatrixKernel.evaluate(x, y);
force[k] = RIJKMatrixKernelClass(dk).evaluate(x, y);
f[counter*3 + k] += force[k] * wx * wy;
}
}
......
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