Skip to content

GitLab

S
ScalFMM
Commit 8ca59f73 authored by BLANCHARD Pierre's avatar BLANCHARD Pierre
Browse files
Options

Update tensorial FMM

parents 596ca731 a78cc320
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 795 additions and 580 deletions
Src/Kernels/Chebyshev/FChebM2LHandler.hpp
View file @ 8ca59f73
......@@ -51,6 +51,16 @@ unsigned int Compress(const FReal epsilon, const unsigned int ninteractions,
* size \f$\ell^3\times r\f$, and \f$316\f$ \f$C_t\f$, each of size \f$r\times
* r\f$.
*
* PB: FChebM2LHandler does not seem to support non_homogeneous kernels!
* In fact nothing appears to handle this here (i.e. adapt scaling and storage
* to MatrixKernelClass::Type). Given the relatively important cost of the
* Chebyshev variant, it is probably a choice not to have implemented this
* feature here but instead in the ChebyshevSym variant. But what if the
* kernel is non homogeneous and non symmetric (e.g. Dislocations)...
*
* TODO Specialize class (see UnifM2LHandler) OR prevent from using this
* class with non homogeneous kernels ?!
*
* @tparam ORDER interpolation order \f$\ell\f$
*/
template <int ORDER, class MatrixKernelClass>
......
Src/Kernels/Interpolation/FInterpMatrixKernel.hpp
View file @ 8ca59f73
......@@ -55,7 +55,7 @@ enum KERNEL_FUNCTION_TYPE {HOMOGENEOUS, NON_HOMOGENEOUS};
*
*/
struct FInterpAbstractMatrixKernel : FNoCopyable
{
{
virtual ~FInterpAbstractMatrixKernel(){} // to remove warning
virtual FReal evaluate(const FPoint&, const FPoint&) const = 0;
// I need both functions because required arguments are not always given
......@@ -68,13 +68,14 @@ struct FInterpAbstractMatrixKernel : FNoCopyable
/// One over r
struct FInterpMatrixKernelR : FInterpAbstractMatrixKernel
{
static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS;
static const KERNEL_FUNCTION_TYPE Type = /*NON_*/HOMOGENEOUS;
static const KERNEL_FUNCTION_IDENTIFIER Identifier = ONE_OVER_R;
static const unsigned int DIM = 1; //PB: dimension of kernel
static const unsigned int NIDX = 1; //PB: number of indices
static const unsigned int NRHS = 1;
static const unsigned int NLHS = 1;
FInterpMatrixKernelR(const unsigned int = 0) {}
FInterpMatrixKernelR(const double = 0.0, const unsigned int = 0) {}
// returns position in reduced storage
int getPosition(const unsigned int) const
......@@ -98,6 +99,19 @@ struct FInterpMatrixKernelR : FInterpAbstractMatrixKernel
{
return FReal(2.) / CellWidth;
}
// FReal getScaleFactor(const FReal, const int) const
// {
// // return 1 because non homogeneous kernel functions cannot be scaled!!!
// return FReal(1.);
// }
//
// FReal getScaleFactor(const FReal) const
// {
// // return 1 because non homogeneous kernel functions cannot be scaled!!!
// return FReal(1.);
// }
};
......@@ -108,10 +122,11 @@ struct FInterpMatrixKernelRR : FInterpAbstractMatrixKernel
static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS;
static const KERNEL_FUNCTION_IDENTIFIER Identifier = ONE_OVER_R_SQUARED;
static const unsigned int DIM = 1; //PB: dimension of kernel
static const unsigned int NIDX = 1; //PB: number of indices
static const unsigned int NRHS = 1;
static const unsigned int NLHS = 1;
FInterpMatrixKernelRR(const unsigned int = 0) {}
FInterpMatrixKernelRR(const double = 0.0, const unsigned int = 0) {}
// returns position in reduced storage
int getPosition(const unsigned int) const
......@@ -145,10 +160,11 @@ struct FInterpMatrixKernelLJ : FInterpAbstractMatrixKernel
static const KERNEL_FUNCTION_TYPE Type = NON_HOMOGENEOUS;
static const KERNEL_FUNCTION_IDENTIFIER Identifier = LENNARD_JONES_POTENTIAL;
static const unsigned int DIM = 1; //PB: dimension of kernel
static const unsigned int NIDX = 1; //PB: number of indices
static const unsigned int NRHS = 1;
static const unsigned int NLHS = 1;
FInterpMatrixKernelLJ(const unsigned int = 0) {}
FInterpMatrixKernelLJ(const double = 0.0, const unsigned int = 0) {}
// returns position in reduced storage
int getPosition(const unsigned int) const
......@@ -193,7 +209,8 @@ struct FInterpMatrixKernel_IOR : FInterpAbstractMatrixKernel
static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS;
static const KERNEL_FUNCTION_IDENTIFIER Identifier = ID_OVER_R;
static const unsigned int DIM = 6; //PB: dimension of kernel
const unsigned int indexTab[12]={0,0,0,1,1,2,
static const unsigned int NIDX = 2; //PB: number of indices
static constexpr unsigned int indexTab[12]={0,0,0,1,1,2,
0,1,2,1,2,2};
static const unsigned int NRHS = 3;
static const unsigned int NLHS = 3;
......@@ -205,7 +222,7 @@ struct FInterpMatrixKernel_IOR : FInterpAbstractMatrixKernel
const unsigned int _i,_j;
FInterpMatrixKernel_IOR(const unsigned int d = 0)
FInterpMatrixKernel_IOR(const double = 0.0, const unsigned int d = 0)
: _i(indexTab[d]), _j(indexTab[d+DIM])
{}
......@@ -225,6 +242,22 @@ struct FInterpMatrixKernel_IOR : FInterpAbstractMatrixKernel
}
void evaluateBlock(const FPoint& x, const FPoint& y, FReal* block) const
{
const FPoint xy(x-y);
const FReal one_over_r = FReal(1.)/xy.norm();
for(unsigned int d=0;d<DIM;++d){
// unsigned int i = indexTab[d];
// unsigned int j = indexTab[d+DIM];
// if(i==j)
block[d] = one_over_r;
// else
// block[d] = 0.0;
}
}
FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
{
const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
......@@ -246,19 +279,20 @@ struct FInterpMatrixKernel_R_IJ : FInterpAbstractMatrixKernel
static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS;
static const KERNEL_FUNCTION_IDENTIFIER Identifier = R_IJ;
static const unsigned int DIM = 6; //PB: dimension of kernel
const unsigned int indexTab[12]={0,0,0,1,1,2,
0,1,2,1,2,2};
static const unsigned int NIDX = 2; //PB: number of indices
static constexpr unsigned int indexTab[DIM*NIDX]={0,0,0,1,1,2,
0,1,2,1,2,2};
static const unsigned int NRHS = 3;
static const unsigned int NLHS = 3;
// store positions in sym tensor
// store positions in sym tensor (when looping over NRHSxNLHS)
const unsigned int applyTab[9]={0,1,2,
1,3,4,
2,4,5};
const unsigned int _i,_j;
FInterpMatrixKernel_R_IJ(const unsigned int d = 0)
FInterpMatrixKernel_R_IJ(const double = 0.0, const unsigned int d = 0)
: _i(indexTab[d]), _j(indexTab[d+DIM])
{}
......@@ -283,8 +317,6 @@ struct FInterpMatrixKernel_R_IJ : FInterpAbstractMatrixKernel
else if(_j==2) rj=xy.getZ();
else throw std::runtime_error("Update j!");
// return xy.getX() * xy.getX() * one_over_r3; //PB: test r0^2/R^3
if(_i==_j)
return one_over_r - ri * ri * one_over_r3;
else
......@@ -292,6 +324,24 @@ struct FInterpMatrixKernel_R_IJ : FInterpAbstractMatrixKernel
}
void evaluateBlock(const FPoint& x, const FPoint& y, FReal* block) const
{
const FPoint xy(x-y);
const FReal one_over_r = FReal(1.)/xy.norm();
const FReal one_over_r3 = one_over_r*one_over_r*one_over_r;
const double r[3] = {xy.getX(),xy.getY(),xy.getZ()};
for(unsigned int d=0;d<DIM;++d){
unsigned int i = indexTab[d];
unsigned int j = indexTab[d+DIM];
if(i==j)
block[d] = one_over_r - r[i] * r[i] * one_over_r3;
else
block[d] = - r[i] * r[j] * one_over_r3;
}
}
FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
{
const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
......@@ -307,27 +357,27 @@ struct FInterpMatrixKernel_R_IJ : FInterpAbstractMatrixKernel
};
/// R_{,ijk}
struct FInterpMatrixKernel_R_IJK : FInterpAbstractMatrixKernel
{
static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS;
static const KERNEL_FUNCTION_IDENTIFIER Identifier = R_IJK;
static const unsigned int DIM = 10; //PB: dimension of kernel
const unsigned int indexTab[30]={0,0,0,1,1,1,2,2,2,0,
0,1,2,0,1,2,0,1,2,1,
0,1,2,0,1,2,0,1,2,2};
static const unsigned int NIDX = 3; //PB: number of indices
static constexpr unsigned int indexTab[DIM*NIDX]={0,0,0,1,1,1,2,2,2,0,
0,1,2,0,1,2,0,1,2,1,
0,1,2,0,1,2,0,1,2,2};
static const unsigned int NRHS = 3;
static const unsigned int NLHS = 3;
// store positions in sym tensor
// store positions in sym tensor (when looping over NRHSxNLHS)
const unsigned int applyTab[27]={0,3,6,3,1,9,6,9,2,
3,1,9,1,4,7,9,7,5,
6,9,2,9,7,5,2,5,8};
const unsigned int _i,_j,_k;
FInterpMatrixKernel_R_IJK(const unsigned int d = 0)
FInterpMatrixKernel_R_IJK(const double = 0.0, const unsigned int d = 0)
: _i(indexTab[d]), _j(indexTab[d+DIM]), _k(indexTab[d+2*DIM])
{}
......
Src/Kernels/Uniform/FAbstractUnifKernel.hpp
View file @ 8ca59f73
......@@ -59,6 +59,8 @@ protected:
const FReal BoxWidth;
/// Width of a leaf cell box
const FReal BoxWidthLeaf;
/// Parameter to pass to matrix kernel (material specific or anything)
const double MatParam;
/**
* Compute center of leaf cell from its tree coordinate.
......@@ -79,14 +81,16 @@ public:
* runtime_error is thrown if the required file is not valid).
*/
FAbstractUnifKernel(const int inTreeHeight,
const FReal inBoxWidth,
const FPoint& inBoxCenter)
const FReal inBoxWidth,
const FPoint& inBoxCenter,
const double inMatParam = 0.0)
: Interpolator(new InterpolatorClass()),
MatrixKernel(new MatrixKernelClass()),
MatrixKernel(new MatrixKernelClass(inMatParam)),
TreeHeight(inTreeHeight),
BoxCorner(inBoxCenter - inBoxWidth / FReal(2.)),
BoxWidth(inBoxWidth),
BoxWidthLeaf(BoxWidth / FReal(FMath::pow(2, inTreeHeight - 1)))
BoxWidthLeaf(BoxWidth / FReal(FMath::pow(2, inTreeHeight - 1))),
MatParam(inMatParam)
{
/* empty */
}
......
Src/Kernels/Uniform/FUnifKernel.hpp
View file @ 8ca59f73
......@@ -45,7 +45,7 @@ class FUnifKernel
: public FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
{
// private types
typedef FUnifM2LHandler<ORDER,MatrixKernelClass> M2LHandlerClass;
typedef FUnifM2LHandler<ORDER,MatrixKernelClass::Type> M2LHandlerClass;
// using from
typedef FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
......@@ -61,17 +61,13 @@ public:
* runtime_error is thrown if the required file is not valid).
*/
FUnifKernel(const int inTreeHeight,
const FReal inBoxWidth,
const FPoint& inBoxCenter)
: FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,
inBoxWidth,
inBoxCenter),
M2LHandler(new M2LHandlerClass())
{
// read precomputed compressed m2l operators from binary file
//M2LHandler->ReadFromBinaryFileAndSet(); // PB: TODO?
M2LHandler->ComputeAndSet();
}
const FReal inBoxWidth,
const FPoint& inBoxCenter)
: FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter),
M2LHandler(new M2LHandlerClass(AbstractBaseClass::MatrixKernel.getPtr(),
inTreeHeight,
inBoxWidth))// PB: for non homogeneous case
{ }
void P2M(CellClass* const LeafCell,
......@@ -133,16 +129,17 @@ public:
const int /*NumSourceCells*/,
const int TreeLevel)
{
const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
const FReal scale(AbstractBaseClass::MatrixKernel.getPtr()->getScaleFactor(CellWidth));
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FComplexe *const TransformedLocalExpansion = TargetCell->getTransformedLocal(idxRhs);
const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
for (int idx=0; idx<343; ++idx){
if (SourceCells[idx]){
M2LHandler->applyFC(idx, CellWidth, SourceCells[idx]->getTransformedMultipole(idxRhs),
M2LHandler->applyFC(idx, TreeLevel, scale,
SourceCells[idx]->getTransformedMultipole(idxRhs),
TransformedLocalExpansion);
}
}
}
......
Src/Kernels/Uniform/FUnifM2LHandler.hpp
View file @ 8ca59f73
This diff is collapsed.
Src/Kernels/Uniform/FUnifTensorialKernel.hpp
View file @ 8ca59f73
......@@ -38,14 +38,15 @@ class FTreeCoordinate;
*
* PB: 3 IMPORTANT remarks !!!
*
* 1) Handling tensorial kernels (DIM,NRHS,NLHS) and having multiple rhs (NVALS)
* are considered 2 separate features and are currently combined.
* 1) Handling tensorial kernels (DIM,NRHS,NLHS) and having multiple rhs
* (NVALS) are considered 2 distinct features and are currently combined.
*
* 2) 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
* 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.
* 2) 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) During precomputation the tensorial matrix kernels are evaluated
* blockwise, but this is not always possible.
* In fact, in the ChebyshevSym variant the matrix kernel needs to be
* evaluated compo-by-compo since we currently use a scalar ACA.
*
* 3) We currently use multiple 1D FFT instead of multidim FFT.
* TODO switch to multidim if relevant in considered range of size
......@@ -66,7 +67,7 @@ class FUnifTensorialKernel
protected://PB: for OptiDis
// private types
typedef FUnifTensorialM2LHandler<ORDER,MatrixKernelClass> M2LHandlerClass;
typedef FUnifTensorialM2LHandler<ORDER,MatrixKernelClass,MatrixKernelClass::Type> M2LHandlerClass;
// using from
typedef FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
......@@ -83,16 +84,13 @@ public:
*/
FUnifTensorialKernel(const int inTreeHeight,
const FReal inBoxWidth,
const FPoint& inBoxCenter)
: FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,
inBoxWidth,
inBoxCenter),
M2LHandler(new M2LHandlerClass())
{
// read precomputed compressed m2l operators from binary file
//M2LHandler->ReadFromBinaryFileAndSet(); // PB: TODO?
M2LHandler->ComputeAndSet();
}
const FPoint& inBoxCenter,
const double inMatParam = 0.0)
: FAbstractUnifKernel< CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter,inMatParam),
M2LHandler(new M2LHandlerClass(AbstractBaseClass::MatrixKernel.getPtr(),
inTreeHeight,
inBoxWidth))// PB: for non homogeneous case
{ }
void P2M(CellClass* const LeafCell,
......@@ -150,6 +148,7 @@ public:
const int TreeLevel)
{
const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
const FReal scale(AbstractBaseClass::MatrixKernel.getPtr()->getScaleFactor(CellWidth));
for(int idxV = 0 ; idxV < NVALS ; ++idxV){
for (int idxLhs=0; idxLhs < nLhs; ++idxLhs){
......@@ -163,12 +162,15 @@ public:
int idxMul = idxV*nRhs + idxRhs;
// update kernel index such that: x_i = K_{ij}y_j
int idxK = idxLhs*nRhs + idxRhs;
// get index in matrix kernel
unsigned int d
= AbstractBaseClass::MatrixKernel.getPtr()->getPosition(idxK);
for (int idx=0; idx<343; ++idx){
if (SourceCells[idx]){
M2LHandler->applyFC(idx, CellWidth,
M2LHandler->applyFC(idx, TreeLevel, scale, d,
SourceCells[idx]->getTransformedMultipole(idxMul),
TransformedLocalExpansion,idxK);
TransformedLocalExpansion);
}
}
......
Tests/Kernels/testUnifAlgorithm.cpp
View file @ 8ca59f73
......@@ -122,12 +122,12 @@ int main(int argc, char* argv[])
{ // begin Lagrange kernel
// accuracy
const unsigned int ORDER = 3;
const unsigned int ORDER = 5;
// typedefs
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleLeaf< ContainerClass > LeafClass;
//typedef FInterpMatrixKernelLJ MatrixKernelClass;
// typedef FInterpMatrixKernelLJ MatrixKernelClass;
typedef FInterpMatrixKernelR MatrixKernelClass;
typedef FUnifCell<ORDER> CellClass;
typedef FOctree<CellClass,ContainerClass,LeafClass> OctreeClass;
......
Tests/Kernels/testUnifTensorialAlgorithm.cpp
View file @ 8ca59f73
......@@ -68,8 +68,8 @@ int main(int argc, char* argv[])
// typedefs
typedef FInterpMatrixKernel_R_IJ MatrixKernelClass;
// typedef FInterpMatrixKernel_IOR MatrixKernelClass;
// typedef FInterpMatrixKernel_R_IJ MatrixKernelClass;
typedef FInterpMatrixKernel_IOR MatrixKernelClass;
const KERNEL_FUNCTION_IDENTIFIER MK_ID = MatrixKernelClass::Identifier;
const unsigned int NRHS = MatrixKernelClass::NRHS;
......
Tests/Utils/testUnifInterpolator.cpp
View file @ 8ca59f73
......@@ -509,11 +509,21 @@ int main(int, char **){
counter++;
}
// std::cout << "Check Potential, forceX, forceY, forceZ " << std::endl;
// for(int idxPart = 0 ; idxPart < 20 ; ++idxPart){
// std::cout << approx_p[idxPart] << ", "<< p[idxPart] << "|| ";
// std::cout << approx_f[idxPart] << ", "<< f[idxPart] << "|| ";
// std::cout << approx_f[idxPart+M] << ", "<< f[idxPart+M] << "|| ";
// std::cout << approx_f[idxPart+2*M] << ", "<< f[idxPart+2*M] << "|| ";
// std::cout << std::endl;
// }
// std::cout << std::endl;
std::cout << "\nPotential error:" << std::endl;
std::cout << "Relative error = " << FMath::FAccurater( p, approx_p, M) << std::endl;
std::cout << "Absolute error = " << FMath::FAccurater( p, approx_p, M) << std::endl;
std::cout << "\nForce error:" << std::endl;
std::cout << "Relative L2 error = " << FMath::FAccurater( f, approx_f, M*3) << std::endl;
std::cout << "Absolute L2 error = " << FMath::FAccurater( f, approx_f, M*3) << std::endl;
std::cout << std::endl;
// free memory
......
Tests/Utils/testUnifTensorialInterpolator.cpp
View file @ 8ca59f73
......@@ -55,12 +55,15 @@
*/
int main(int, char **){
typedef FP2PParticleContainer<> ContainerClass;
typedef FSimpleLeaf<ContainerClass> LeafClass;
// typedef FInterpMatrixKernelR MatrixKernelClass;
typedef FInterpMatrixKernel_R_IJ MatrixKernelClass;
typedef FInterpMatrixKernel_R_IJK RIJKMatrixKernelClass; // PB: force computation
const unsigned int dim = MatrixKernelClass::DIM;
const unsigned int nrhs = MatrixKernelClass::NRHS;
const unsigned int nlhs = MatrixKernelClass::NLHS;
typedef FP2PParticleContainer<nrhs,nlhs> ContainerClass;
typedef FSimpleLeaf<ContainerClass> LeafClass;
///////////////////////What we do/////////////////////////////
std::cout << "\nTask: Compute interactions between source particles in leaf Y and target\n";
......@@ -68,10 +71,6 @@ int main(int, char **){
std::cout << " direct computation.\n" << std::endl;
//////////////////////////////////////////////////////////////
const unsigned int dim = MatrixKernelClass::DIM;
const unsigned int nrhs = MatrixKernelClass::NRHS;
const unsigned int nlhs = MatrixKernelClass::NLHS;
const FReal FRandMax = FReal(RAND_MAX);
FTic time;
......@@ -86,12 +85,13 @@ int main(int, char **){
std::cout << "Fill the leaf X of width " << width
<< " centered at cx=" << cx << " with M=" << M << " target particles" << std::endl;
{
for(unsigned long i=0; i<M; ++i){
FReal x = (FReal(rand())/FRandMax - FReal(.5)) * width + cx.getX();
FReal y = (FReal(rand())/FRandMax - FReal(.5)) * width + cx.getY();
FReal z = (FReal(rand())/FRandMax - FReal(.5)) * width + cx.getZ();
// TODO: find a way to push vectorial attribute of arbitrary dim (nlhs)
X.push(FPoint(x, y, z), FReal(rand())/FRandMax);
// PB: need to know the actual value of NLHS
X.push(FPoint(x, y, z), FReal(rand())/FRandMax, FReal(rand())/FRandMax, FReal(rand())/FRandMax);
}
}
......@@ -109,8 +109,8 @@ int main(int, char **){
FReal x = (FReal(rand())/FRandMax - FReal(.5)) * width + cy.getX();
FReal y = (FReal(rand())/FRandMax - FReal(.5)) * width + cy.getY();
FReal z = (FReal(rand())/FRandMax - FReal(.5)) * width + cy.getZ();
// TODO: find a way to push vectorial attribute of arbitrary dim (nrhs)
Y.push(FPoint(x, y, z), FReal(rand())/FRandMax);
// PB: need to know the actual value of NRHS
Y.push(FPoint(x, y, z), FReal(rand())/FRandMax, FReal(rand())/FRandMax, FReal(rand())/FRandMax);
}
}
......@@ -118,7 +118,7 @@ int main(int, char **){
////////////////////////////////////////////////////////////////////
// approximative computation
const unsigned int ORDER = 4;
const unsigned int ORDER = 5;
const unsigned int nnodes = TensorTraits<ORDER>::nnodes;
typedef FUnifInterpolator<ORDER,MatrixKernelClass> InterpolatorClass;
InterpolatorClass S;
......@@ -132,8 +132,7 @@ int main(int, char **){
// Anterpolate: W_n = \sum_j^N S(y_j,\bar y_n) * w_j
FReal W[nrhs*nnodes]; // multipole expansion
// tensorial case interpolate same Y for each component
for (unsigned int idRhs=0; idRhs<nrhs; ++idRhs)
S.applyP2M(cy, width, W + idRhs*nnodes, Y.getSrc()); // the multipole expansions are set to 0 in S.applyP2M
S.applyP2M(cy, width, W, Y.getSrc()); // the multipole expansions are set to 0 in S.applyP2M
std::cout << "took " << time.tacAndElapsed() << "s" << std::endl;
std::cout << "M2L ... " << std::flush;
......@@ -499,15 +498,13 @@ int main(int, char **){
std::cout << "L2P (potential) ... " << std::flush;
time.tic();
// Interpolate p_i = \sum_m^L S(x_i,\bar x_m) * F_m
for (unsigned int idLhs=0; idLhs<nlhs; ++idLhs)
S.applyL2P(cx, width, F+idLhs*nnodes, X.getTargets());
S.applyL2P(cx, width, F, X.getTargets());
std::cout << "took " << time.tacAndElapsed() << "s" << std::endl;
std::cout << "L2P (forces) ... " << std::flush;
time.tic();
// Interpolate f_i = \sum_m^L P(x_i,\bar x_m) * F_m
for (unsigned int idLhs=0; idLhs<nlhs; ++idLhs)
S.applyL2PGradient(cx, width, F+idLhs*nnodes, X.getTargets());
S.applyL2PGradient(cx, width, F, X.getTargets());
std::cout << "took " << time.tacAndElapsed() << "s" << std::endl;
////////////////////////////////////////////////////////////////////
......@@ -517,13 +514,21 @@ int main(int, char **){
std::cout << "Compute interactions directly ..." << std::endl;
time.tic();
FReal* approx_f = new FReal [M * 3];
FReal* f = new FReal [M * 3];
FBlas::setzero(M*3, f);
FReal** approx_f = new FReal* [nlhs];
FReal** f = new FReal* [nlhs];
for (unsigned int i=0; i<nrhs; ++i){
approx_f[i] = new FReal [M * 3];
f[i] = new FReal [M * 3];
FBlas::setzero(M*3, f[i]);
}
FReal* approx_p = new FReal[M];
FReal* p = new FReal[M];
FBlas::setzero(M, p);
FReal** approx_p = new FReal* [nlhs];
FReal** p = new FReal* [nlhs];
for (unsigned int i=0; i<nrhs; ++i){
approx_p[i] = new FReal [M];
p[i] = new FReal [M];
FBlas::setzero(M, p[i]);
}
{ // start direct computation
unsigned int counter = 0;
......@@ -532,14 +537,18 @@ int main(int, char **){
const FPoint x = FPoint(X.getSrc()->getPositions()[0][idxPartX],
X.getSrc()->getPositions()[1][idxPartX],
X.getSrc()->getPositions()[2][idxPartX]);
const FReal wx = X.getSrc()->getPhysicalValues()[idxPartX];
const FReal wx[nrhs] = {X.getSrc()->getPhysicalValues(0)[idxPartX],
X.getSrc()->getPhysicalValues(1)[idxPartX],
X.getSrc()->getPhysicalValues(2)[idxPartX]};
for(int idxPartY = 0 ; idxPartY < Y.getSrc()->getNbParticles() ; ++idxPartY){
const FPoint y = FPoint(Y.getSrc()->getPositions()[0][idxPartY],
Y.getSrc()->getPositions()[1][idxPartY],
Y.getSrc()->getPositions()[2][idxPartY]);
const FReal wy = Y.getSrc()->getPhysicalValues()[idxPartY];
const FReal wy[nrhs] = {Y.getSrc()->getPhysicalValues(0)[idxPartY],
Y.getSrc()->getPhysicalValues(1)[idxPartY],
Y.getSrc()->getPhysicalValues(2)[idxPartY]};
// // 1/R
// const FReal one_over_r = MatrixKernel.evaluate(x, y);
//
......@@ -558,14 +567,14 @@ int main(int, char **){
unsigned int d = MatrixKernel.getPosition(i*nrhs+j);
const FReal rij = MatrixKernelClass(d).evaluate(x, y);
// potential
p[counter] += rij * wy;
p[i][counter] += rij * wy[j];
// force
FReal force[3];
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);
force[k] = RIJKMatrixKernelClass(dk).evaluate(x, y);
f[counter*3 + k] += force[k] * wx * wy;
f[i][counter*3 + k] += force[k] * wx[j] * wy[j];
}
}
......@@ -582,22 +591,41 @@ int main(int, char **){
////////////////////////////////////////////////////////////////////
unsigned int counter = 0;
for(int idxPartX = 0 ; idxPartX < X.getSrc()->getNbParticles() ; ++idxPartX){
approx_p[counter] = X.getSrc()->getPotentials()[idxPartX];
const FPoint force = FPoint(X.getSrc()->getForcesX()[idxPartX],
X.getSrc()->getForcesY()[idxPartX],
X.getSrc()->getForcesZ()[idxPartX]);
approx_f[counter*3 + 0] = force.getX();
approx_f[counter*3 + 1] = force.getY();
approx_f[counter*3 + 2] = force.getZ();
for (unsigned int i=0; i<nlhs; ++i){
approx_p[i][counter] = X.getSrc()->getPotentials(i)[idxPartX];
const FPoint force = FPoint(X.getSrc()->getForcesX(i)[idxPartX],
X.getSrc()->getForcesY(i)[idxPartX],
X.getSrc()->getForcesZ(i)[idxPartX]);
approx_f[i][counter*3 + 0] = force.getX();
approx_f[i][counter*3 + 1] = force.getY();
approx_f[i][counter*3 + 2] = force.getZ();
}
counter++;
}
std::cout << "\nPotential error:" << std::endl;
std::cout << "Relative error = " << FMath::FAccurater( p, approx_p, M) << std::endl;
// std::cout << "Check Potential, forceX, forceY, forceZ " << std::endl;
// for (unsigned int i=0; i<nlhs; ++i){
// std::cout<< "idxLhs="<< i << std::endl;
// for(int idxPart = 0 ; idxPart < 20 ; ++idxPart){
// std::cout << approx_p[i][idxPart] << ", "<< p[i][idxPart] << "|| ";
// std::cout << approx_f[i][idxPart] << ", "<< f[i][idxPart] << "|| ";
// std::cout << approx_f[i][idxPart+M] << ", "<< f[i][idxPart+M] << "|| ";
// std::cout << approx_f[i][idxPart+2*M] << ", "<< f[i][idxPart+2*M] << "|| ";
// std::cout << std::endl;
// }
// std::cout << std::endl;
// }
// std::cout << std::endl;
std::cout << "\nForce error:" << std::endl;
std::cout << "Relative L2 error = " << FMath::FAccurater( f, approx_f, M*3) << std::endl;
std::cout << "\nPotential error:" << std::endl;
std::cout << "Relative error X = " << FMath::FAccurater( p[0], approx_p[0], M) << std::endl;
std::cout << "Relative error Y = " << FMath::FAccurater( p[1], approx_p[1], M) << std::endl;
std::cout << "Relative error Z = " << FMath::FAccurater( p[2], approx_p[2], M) << std::endl;
std::cout << "\nForce error (Something is wrong here => TOFIX):" << std::endl;
std::cout << "Relative L2 error X = " << FMath::FAccurater( f[0], approx_f[0], M*3) << std::endl;