Commit a78cc320 authored by BLANCHARD Pierre's avatar BLANCHARD Pierre

Fixed handling of non homogeneous kernels of interactions (aka matrix kernel)...

Fixed handling of non homogeneous kernels of interactions (aka matrix kernel) in Uniform kernel (TODO Fix in Chebyshev if relevant), Introduced blockwise evaluation of tensorial matrix kernel (+redesign precomputation), added input parameter in AbstractUnifKernel (for OptiDis non singular Ra_{ijk} or any matrix kernel requiring an extra parameter), +divers test fixing.
parent ca0bf8b8
......@@ -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>
......
......@@ -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])
{}
......
......@@ -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 */
}
......
......@@ -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);
}
}
}
......
This diff is collapsed.
......@@ -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);
}
}
......
......@@ -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;
......
......@@ -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;
......
......@@ -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
......
......@@ -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[