Commit ca884022 authored by BRAMAS Berenger's avatar BRAMAS Berenger
Browse files

Clean rotation kernel is OK, now the efficient rotation kernel is in progress

parent 3ae6306e
......@@ -26,26 +26,14 @@
* This kernels is a complete rotation based kernel with spherical
* harmonic.
*
* This kernel is still not optimized for a pedagogic purpose.
* Please with refer to the optimized version to run real simulation.
*
* Here there is no real precomputation no optimization.
*
* The fallowing code has to be insert inside this class to print
* the details of a vector during the computation if needed.
* @code
* for(int l = 0 ; l <= P ; ++l ){
* for(int m = 0 ; m <= l ; ++m ){
* printf("3 - [%d][%d] %f i%f\t", l, m, target_u[atLm(l,m)].getReal(), target_u[atLm(l,m)].getImag());
* }
* printf("\n");
* }
* @endcode
* Here is the optimizated kernel, please refer to FRotationOriginalKernel
* to see the non optimized easy to understand kernel.
*/
template< class ParticleClass, class CellClass, class ContainerClass, int P>
class FRotationKernel : public FAbstractKernels<ParticleClass,CellClass,ContainerClass> {
//< Size of the data array computed using a suite relation
static const int SizeArray = ((P+2)*(P+1))/2;
static const int P2 = P*2;
///////////////////////////////////////////////////////
// Object attributes
......@@ -61,6 +49,17 @@ class FRotationKernel : public FAbstractKernels<ParticleClass,CellClass,Containe
FSmartPointer<FSpherical> childrenPosition; //< The transfers between level
FSmartPointer<FSpherical> interactionsPosition; //< The transfers at a same level
FReal factorials[P2+1];
FSmartPointer<FReal[P+1]> M2MTranslationCoef;
FSmartPointer<FReal[343][P+1]> M2LTranslationCoef;
FSmartPointer<FReal[P+1]> L2LTranslationCoef;
FComplexe rotationExpMinusImPhi[8][SizeArray];
FComplexe rotationExpImPhi[8][SizeArray];
FSmartPointer<FComplexe[343][SizeArray]> M2LAroundZCoef;
///////////////////////////////////////////////////////
// d_lmk
///////////////////////////////////////////////////////
......@@ -79,20 +78,183 @@ class FRotationKernel : public FAbstractKernels<ParticleClass,CellClass,Containe
return T(0);
}
/** Return the factorial of a number */
FReal fact(const int a){
if(a<0) printf("Error factorial negative!! a=%d\n",a);
FReal result = 1;
for(int i = 1 ; i <= a ; ++i){
result *= FReal(i);
/** Compute the factorial from 0 to P*2 */
void precomputeFactorials(){
factorials[0] = 1;
FReal fidx = 1;
for(int idx = 1 ; idx <= P2 ; ++idx, ++fidx){
factorials[idx] = fidx * factorials[idx-1];
}
}
void precomputeTranslationCoef(){
M2MTranslationCoef = new FReal[treeHeight-1][P+1];
L2LTranslationCoef = new FReal[treeHeight-1][P+1];
FReal widthAtLevel = boxWidth/4;
for( int idxLevel = 0 ; idxLevel < treeHeight - 1 ; ++idxLevel){
FReal b = FMath::Sqrt(widthAtLevel*widthAtLevel*3);
FReal bPowIdx = 1.0;
FReal minus_1_pow_idx = 1.0;
for(int idx = 0 ; idx <= P ; ++idx){
// coef m2m = (-b)^j/j!
M2MTranslationCoef[idxLevel][idx] = minus_1_pow_idx * bPowIdx / factorials[idx];
// coef l2l = b^j/j!
L2LTranslationCoef[idxLevel][idx] = bPowIdx / factorials[idx];
// increase
bPowIdx *= b;
minus_1_pow_idx = -minus_1_pow_idx;
}
widthAtLevel /= 2;
}
M2LTranslationCoef = new FReal[treeHeight][343][P+1];
for( int idxLevel = 0 ; idxLevel < treeHeight ; ++idxLevel){
for(int idxInteraction = 0 ; idxInteraction < 343 ; ++idxInteraction){
int x = (idxInteraction/7*7)-3;
int y = ((idxInteraction - (x+3)*7*7)/7) - 3;
int z = idxInteraction - ((x+3)*7 + (y+3))*7;//idxInteraction%7;
if( x < -1 || 1 < x || y < -1 || 1 < y || z < -1 || 1 < 1 ){
const FReal b = getSphericalInteraction(idxLevel, idxInteraction).getR();
FReal bPowIdx1 = b;
for(int idx = 0 ; idx <= P ; ++idx){
// factorials[j+l] / FMath::pow(b,j+l+1)
M2LTranslationCoef[idxLevel][idxInteraction][idx] = factorials[idx] / bPowIdx1;
bPowIdx1 *= b;
}
}
}
}
}
///////////////////////////////////////////////////////
// Position precomputation
///////////////////////////////////////////////////////
/** This function precompute the position between cells */
void preComputePosition(){
// Compute the parent/children relation for M2M L2L
childrenPosition = new FSpherical[8 * (treeHeight-1)];
{
FReal subBoxWidth = widthAtLeafLevelDiv2;
// For each
for(int idxLevel = treeHeight-2 ; idxLevel > 0 ; --idxLevel){
for(int idxChild = 0 ; idxChild < 8 ; ++idxChild ){
// coord from child to parent
const FReal x = FReal((idxChild&4)? -subBoxWidth : subBoxWidth);
const FReal y = FReal((idxChild&2)? -subBoxWidth : subBoxWidth);
const FReal z = FReal((idxChild&1)? -subBoxWidth : subBoxWidth);
const FPoint relativePosition( x , y , z );
childrenPosition[(idxLevel-1)*8 + idxChild] = FSpherical(relativePosition);
}
subBoxWidth *= FReal(2.0);
}
}
// Compute the interaction relations for M2L
interactionsPosition = new FSpherical[343 * (treeHeight-1)];
{
FReal boxWidthAtLevel = widthAtLeafLevel;
for(int idxLevel = treeHeight-1 ; idxLevel > 0 ; --idxLevel){
for(int idxX = -3 ; idxX <= 3 ; ++idxX ){
for(int idxY = -3 ; idxY <= 3 ; ++idxY ){
for(int idxZ = -3 ; idxZ <= 3 ; ++idxZ ){
if( idxX != 0 || idxY != 0 || idxZ != 0 ){
const FPoint relativePosition( -FReal(idxX)*boxWidthAtLevel,
-FReal(idxY)*boxWidthAtLevel,
-FReal(idxZ)*boxWidthAtLevel);
const int position = ((( (idxX+3) * 7) + (idxY+3))) * 7 + idxZ + 3;
interactionsPosition[(idxLevel-1)*343 + position] = FSpherical(relativePosition);
}
}
}
}
boxWidthAtLevel *= FReal(2.0);
}
}
}
/** This function rotate a multipole vector by an angle azimuth phi
* The formula used is present in several paper, but we refer to
* Implementation of rotation-based operators for Fast Multipole Method in X10
* At page 5 .1
* \f[
* O_{l,m}( \alpha, \beta + \phi ) = e^{-i \phi m} O_{l,m}( \alpha, \beta )
* \f]
* The computation is simply a multiplication per a complex number \f$ e^{-i \phi m} \f$
* Phi should be in [0,2pi]
*/
/** This function rotate a local vector by an angle azimuth phi
* The formula used is present in several paper, but we refer to
* Implementation of rotation-based operators for Fast Multipole Method in X10
* At page 5 .1
* \f[
* M_{l,m}( \alpha, \beta + \phi ) = e^{i \phi m} M_{l,m}( \alpha, \beta )
* \f]
* The computation is simply a multiplication per a complex number \f$ e^{i \phi m} \f$
* Phi should be in [0,2pi]
*/
void precomputeRotationVectors(){
const int index_P0 = atLm(P,0);
for(int idxChild = 0 ; idxChild < 8 ; ++idxChild){
const FReal x = FReal((idxChild&4)? -boxWidth : boxWidth);
const FReal y = FReal((idxChild&2)? -boxWidth : boxWidth);
const FReal z = FReal((idxChild&1)? -boxWidth : boxWidth);
const FPoint relativePosition( x , y , z );
FSpherical sph(relativePosition);
// compute the last part with l == P
{
int index_lm = index_P0;
for(int m = 0 ; m <= P ; ++m, ++index_lm ){
const FReal mphi = (sph.getAzimuth() + FMath::FPiDiv2) * FReal(m);
// O_{l,m}( \alpha, \beta + \phi ) = e^{-i \phi m} O_{l,m}( \alpha, \beta )
rotationExpMinusImPhi[idxChild][index_lm].setRealImag(FMath::Cos(-mphi), FMath::Sin(-mphi));
// M_{l,m}( \alpha, \beta + \phi ) = e^{i \phi m} M_{l,m}( \alpha, \beta )
rotationExpImPhi[idxChild][index_lm].setRealImag(FMath::Cos(mphi), FMath::Sin(mphi));
}
}
// Then copy
{
int index_lm = 0;
for(int l = 0 ; l < P ; ++l){
FMemUtils::copyall(rotationExpMinusImPhi[idxChild] + index_lm,
rotationExpMinusImPhi[idxChild] + index_P0,
l + 1);
FMemUtils::copyall(rotationExpImPhi[idxChild] + index_lm,
rotationExpImPhi[idxChild] + index_P0,
l + 1);
index_lm += l + 1;
}
}
}
return result;
}
/** To get the right spherical object from level and child position */
FSpherical getSphericalChild(const int idxLevel, const int position) const {
return childrenPosition[(idxLevel-1)*8 + position];
}
/** To get the right spherical object from level and interaction position */
FSpherical getSphericalInteraction(const int idxLevel, const int position) const {
return interactionsPosition[(idxLevel-1)*343 + position];
}
/** Return the position of a leaf from its tree coordinate */
FPoint getLeafCenter(const FTreeCoordinate coordinate) const {
return FPoint(
FReal(coordinate.getX()) * widthAtLeafLevel + widthAtLeafLevelDiv2 + boxCorner.getX(),
FReal(coordinate.getY()) * widthAtLeafLevel + widthAtLeafLevelDiv2 + boxCorner.getY(),
FReal(coordinate.getZ()) * widthAtLeafLevel + widthAtLeafLevelDiv2 + boxCorner.getZ());
}
/** Return the combine of a paire of number */
FReal combin(const int& a, const int& b){
if(a-b<0) printf("Error combi negative!! a=%d b=%d\n",a,b);
return fact(a) / (fact(b)*fact(a-b));
return factorials[a] / (factorials[b]*factorials[a-b]);
}
/** To access the rotation matrix value with an analytical computation on the fly
......@@ -115,7 +277,7 @@ class FRotationKernel : public FAbstractKernels<ParticleClass,CellClass,Containe
for(int n = FMath::Max(-(m+k),0) ; n <= FMath::Min(l-m,l-k) ; ++n){
sum += FMath::pow(FReal(-1.0), l-m-n) * combin(l-k, n) * combin(l+k, l-m-n) * FMath::pow(FReal(1.0)+cosTheta,n) * FMath::pow(FReal(1.0)-cosTheta, l-m-n);
}
return (FReal(1.0)/FMath::pow(FReal(2.0),l)) * FMath::Sqrt((fact(l-m)*fact(l+m))/(fact(l-k)*fact(l+k)))
return (FReal(1.0)/FMath::pow(FReal(2.0),l)) * FMath::Sqrt((factorials[l-m]*factorials[l+m])/(factorials[l-k]*factorials[l+k]))
* FMath::pow(FReal(1.0) + FReal(Sgn(k))*cosTheta, FMath::Abs(k)) * FMath::pow(sinTheta, m - FMath::Abs(k)) * sum;
}
// P^{l}_{m,k} = (-1)^{m+k} P^{l}_{k,m} , l \ge 0, -l \leq m \le 0, |m| \leq k \leq l
......@@ -153,22 +315,35 @@ class FRotationKernel : public FAbstractKernels<ParticleClass,CellClass,Containe
*/
void computeLegendre(FReal legendre[], const FReal inCosTheta, const FReal inSinTheta){
const FReal invSinTheta = -inSinTheta;
legendre[atLm(0,0)] = 1.0; // P_0,0(1) = 1
legendre[atLm(1,0)] = inCosTheta; // P_1,0 = cos(theta)
legendre[atLm(1,1)] = invSinTheta; // P_1,1 = -sin(theta)
legendre[0] = 1.0; // P_0,0(1) = 1
legendre[1] = inCosTheta; // P_1,0 = cos(theta)
legendre[2] = invSinTheta; // P_1,1 = -sin(theta)
// work with pointers
FReal* FRestrict legendre_l1_m1 = legendre; // P{l-2,m} starts with P_{0,0}
FReal* FRestrict legendre_l1_m = legendre + 1; // P{l-1,m} starts with P_{1,0}
FReal* FRestrict legendre_lm = legendre + 3; // P{l,m} starts with P_{2,0}
// Compute using recurrence
for(int l = 2; l <= P ; ++l ){
FReal l2_minus_1 = 3; // 2 * l - 1
FReal fl = FReal(2.0);// To get 'l' as a float
for(int l = 2; l <= P ; ++l, ++fl ){
FReal lm_minus_1 = fl - FReal(1.0); // l + m - 1
FReal l_minus_m = fl; // l - m
for( int m = 0; m < l - 1 ; ++m ){
// P_{l,m} = \frac{(2l-1) cos( \theta ) P_{l-1,m} - (l+m-1) P_{l-2,m}x}{(l-k)}
legendre[atLm(l,m)] = (FReal(2*l-1) * inCosTheta * legendre[atLm(l-1,m)] - FReal( l + m - 1 ) * legendre[atLm(l-2,m)] )
/ FReal( l - m );
// P_{l,m} = \frac{(2l-1) cos( \theta ) P_{l-1,m} - (l+m-1) P_{l-2,m}x}{(l-m)}
*(legendre_lm++) = (l2_minus_1 * inCosTheta * (*legendre_l1_m++) - (lm_minus_1++) * (*legendre_l1_m1++) )
/ (l_minus_m--);
}
// P_{l,l-1} = (2l-1) cos( \theta ) P_{l-1,l-1}
legendre[atLm(l,l-1)] = FReal(2*l-1) * inCosTheta * legendre[atLm(l-1,l-1)];
*(legendre_lm++) = l2_minus_1 * inCosTheta * (*legendre_l1_m);
// P_{l,l} = (2l-1) sin( \theta ) P_{l-1,l-1}
legendre[atLm(l,l)] = FReal(2*l-1) * invSinTheta * legendre[atLm(l-1,l-1)];
*(legendre_lm++) = l2_minus_1 * invSinTheta * (*legendre_l1_m);
// goto P_{l-1,0}
++legendre_l1_m;
l2_minus_1 += FReal(2.0); // 2 * l - 1 => progress by two
}
}
......@@ -176,49 +351,6 @@ class FRotationKernel : public FAbstractKernels<ParticleClass,CellClass,Containe
// Real rotation
///////////////////////////////////////////////////////
/** This function rotate a multipole vector by an angle azimuth phi
* The formula used is present in several paper, but we refer to
* Implementation of rotation-based operators for Fast Multipole Method in X10
* At page 5 .1
* \f[
* O_{l,m}( \alpha, \beta + \phi ) = e^{-i \phi m} O_{l,m}( \alpha, \beta )
* \f]
* The computation is simply a multiplication per a complex number \f$ e^{-i \phi m} \f$
* Phi should be in [0,2pi]
*/
void rotateMultipoleAroundZ(FComplexe vec[], const FReal phi){
FComplexe cell_rotate[SizeArray];
for(int l = 0 ; l <= P ; ++l){
for(int m = 0 ; m <= l ; ++m ){
// O_{l,m}( \alpha, \beta + \phi ) = e^{-i \phi m} O_{l,m}( \alpha, \beta )
const FComplexe exp_minus_imphi(FMath::Cos(-phi * FReal(m)), FMath::Sin(-phi * FReal(m)));
cell_rotate[atLm(l,m)].equalMul(exp_minus_imphi , vec[atLm(l,m)]);
}
}
FMemUtils::copyall(vec,cell_rotate,SizeArray);
}
/** This function rotate a local vector by an angle azimuth phi
* The formula used is present in several paper, but we refer to
* Implementation of rotation-based operators for Fast Multipole Method in X10
* At page 5 .1
* \f[
* M_{l,m}( \alpha, \beta + \phi ) = e^{i \phi m} M_{l,m}( \alpha, \beta )
* \f]
* The computation is simply a multiplication per a complex number \f$ e^{i \phi m} \f$
* Phi should be in [0,2pi]
*/
void rotateTaylorAroundZ(FComplexe vec[], const FReal phi){
FComplexe cell_rotate[SizeArray];
for(int l = 0 ; l <= P ; ++l){
for(int m = 0 ; m <= l ; ++m ){
// M_{l,m}( \alpha, \beta + \phi ) = e^{i \phi m} M_{l,m}( \alpha, \beta )
const FComplexe exp_imphi(FMath::Cos(phi * FReal(m)), FMath::Sin(phi * FReal(m)));
cell_rotate[atLm(l,m)].equalMul(exp_imphi , vec[atLm(l,m)]);
}
}
FMemUtils::copyall(vec,cell_rotate,SizeArray);
}
/** This function rotate a multipole vector by an angle inclination \theta
* The formula used is present in several paper, but we refer to
......@@ -243,14 +375,14 @@ class FRotationKernel : public FAbstractKernels<ParticleClass,CellClass,Containe
// O_{l,-m} = \bar{ O_{l,m} } (-1)^m
const FReal d_lmk = d_lmk_analytical(FMath::Cos(theta),FMath::Sin(theta),l,m,k);
// \sqrt{ \frac{(l-k)!(l+k)!}{(l-|m|)!(l+|m|)!} }
const FReal factor = FMath::Sqrt((fact(l-k)*fact(l+k))/(fact(l-abs(m))*fact(l+abs(m))));
const FReal factor = FMath::Sqrt((factorials[l-k]*factorials[l+k])/(factorials[l-abs(m)]*factorials[l+abs(m)]));
w_lkm_real += FMath::pow(FReal(-1.0),-k) * factor * d_lmk * vec[atLm(l,-k)].getReal(); // k<0 => Conjugate * -1^k
w_lkm_imag -= FMath::pow(FReal(-1.0),-k) * factor * d_lmk * vec[atLm(l,-k)].getImag(); // k<0 => Conjugate * -1^k
}
for(int k = 0 ; k <= l ; ++k){
const FReal d_lmk = d_lmk_analytical(FMath::Cos(theta),FMath::Sin(theta),l,m,k);
// \sqrt{ \frac{(l-k)!(l+k)!}{(l-|m|)!(l+|m|)!} }
const FReal factor = FMath::Sqrt((fact(l-k)*fact(l+k))/(fact(l-abs(m))*fact(l+abs(m))));
const FReal factor = FMath::Sqrt((factorials[l-k]*factorials[l+k])/(factorials[l-abs(m)]*factorials[l+abs(m)]));
w_lkm_real += factor * d_lmk * vec[atLm(l,k)].getReal();
w_lkm_imag += factor * d_lmk * vec[atLm(l,k)].getImag();
}
......@@ -283,14 +415,14 @@ class FRotationKernel : public FAbstractKernels<ParticleClass,CellClass,Containe
// M_{l,-m} = \bar{ M_{l,m} } (-1)^m
const FReal d_lmk = d_lmk_analytical(FMath::Cos(theta),FMath::Sin(theta),l,m,k);
// \sqrt{ \frac{(l-|m|)!(l+|m|)!}{(l-k)!(l+k)!} }
const FReal factor = FMath::Sqrt((fact(l-abs(m))*fact(l+abs(m)))/(fact(l-k)*fact(l+k)));
const FReal factor = FMath::Sqrt((factorials[l-abs(m)]*factorials[l+abs(m)])/(factorials[l-k]*factorials[l+k]));
w_lkm_real += FMath::pow(FReal(-1.0),-k) * factor * d_lmk * vec[atLm(l,-k)].getReal(); // k<0 => Conjugate * -1^k
w_lkm_imag -= FMath::pow(FReal(-1.0),-k) * factor * d_lmk * vec[atLm(l,-k)].getImag(); // k<0 => Conjugate * -1^k
}
for(int k = 0 ; k <= l ; ++k){
const FReal d_lmk = d_lmk_analytical(FMath::Cos(theta),FMath::Sin(theta),l,m,k);
// \sqrt{ \frac{(l-|m|)!(l+|m|)!}{(l-k)!(l+k)!} }
const FReal factor = FMath::Sqrt((fact(l-abs(m))*fact(l+abs(m)))/(fact(l-k)*fact(l+k)));
const FReal factor = FMath::Sqrt((factorials[l-abs(m)]*factorials[l+abs(m)])/(factorials[l-k]*factorials[l+k]));
w_lkm_real += factor * d_lmk * vec[atLm(l,k)].getReal();
w_lkm_imag += factor * d_lmk * vec[atLm(l,k)].getImag();
}
......@@ -308,7 +440,7 @@ class FRotationKernel : public FAbstractKernels<ParticleClass,CellClass,Containe
* Rotation are not commutative so we have to do it in the right order
*/
void rotateMultipole(FComplexe vec[], const FReal azimuth, const FReal inclination){
rotateMultipoleAroundZ(vec,(FMath::FPiDiv2 + azimuth));
//rotateMultipoleAroundZ(vec,(FMath::FPiDiv2 + azimuth));
rotateMultipoleAroundY(vec,inclination);
}
/** This function rotate a multipole vector by a angles inclination & azimuth
......@@ -320,7 +452,7 @@ class FRotationKernel : public FAbstractKernels<ParticleClass,CellClass,Containe
*/
void deRotateMultipole(FComplexe vec[], const FReal azimuth, const FReal inclination){
rotateMultipoleAroundY(vec,-inclination);
rotateMultipoleAroundZ(vec,-(FMath::FPiDiv2 + azimuth));
//rotateMultipoleAroundZ(vec,-(FMath::FPiDiv2 + azimuth));
}
/** This function rotate a local vector by a angles inclination & azimuth
......@@ -331,7 +463,7 @@ class FRotationKernel : public FAbstractKernels<ParticleClass,CellClass,Containe
* Rotation are not commutative so we have to do it in the right order
*/
void rotateTaylor(FComplexe vec[], const FReal azimuth, const FReal inclination){
rotateTaylorAroundZ(vec,(FMath::FPiDiv2 + azimuth));
//rotateTaylorAroundZ(vec,(FMath::FPiDiv2 + azimuth));
rotateTaylorAroundY(vec,inclination);
}
/** This function rotate a local vector by a angles inclination & azimuth
......@@ -343,72 +475,13 @@ class FRotationKernel : public FAbstractKernels<ParticleClass,CellClass,Containe
*/
void deRotateTaylor(FComplexe vec[], const FReal azimuth, const FReal inclination){
rotateTaylorAroundY(vec,-inclination);
rotateTaylorAroundZ(vec,-(FMath::FPiDiv2 + azimuth));
//rotateTaylorAroundZ(vec,-(FMath::FPiDiv2 + azimuth));
}
///////////////////////////////////////////////////////
// Position precomputation
///////////////////////////////////////////////////////
/** This function precompute the position between cells */
void preComputePosition(){
// Compute the parent/children relation for M2M L2L
childrenPosition = new FSpherical[8 * (treeHeight-1)];
{
FReal subBoxWidth = widthAtLeafLevelDiv2;
// For each
for(int idxLevel = treeHeight-2 ; idxLevel > 0 ; --idxLevel){
for(int idxChild = 0 ; idxChild < 8 ; ++idxChild ){
// coord from child to parent
const FReal x = FReal((idxChild&4)? -subBoxWidth : subBoxWidth);
const FReal y = FReal((idxChild&2)? -subBoxWidth : subBoxWidth);
const FReal z = FReal((idxChild&1)? -subBoxWidth : subBoxWidth);
const FPoint relativePosition( x , y , z );
childrenPosition[(idxLevel-1)*8 + idxChild] = FSpherical(relativePosition);
}
subBoxWidth *= FReal(2.0);
}
static void ComplexeArrayMulEqual(FComplexe dest[], const FComplexe src[], const int sizeOfArray){
for(int idx = 0 ; idx < sizeOfArray ; ++idx) {
dest[idx] *= src[idx];
}
// Compute the interaction relations for M2L
interactionsPosition = new FSpherical[343 * (treeHeight-1)];
{
FReal boxWidthAtLevel = widthAtLeafLevel;
for(int idxLevel = treeHeight-1 ; idxLevel > 0 ; --idxLevel){
for(int idxX = -3 ; idxX <= 3 ; ++idxX ){
for(int idxY = -3 ; idxY <= 3 ; ++idxY ){
for(int idxZ = -3 ; idxZ <= 3 ; ++idxZ ){
if( idxX != 0 || idxY != 0 || idxZ != 0 ){
const FPoint relativePosition( -FReal(idxX)*boxWidthAtLevel,
-FReal(idxY)*boxWidthAtLevel,
-FReal(idxZ)*boxWidthAtLevel);
const int position = ((( (idxX+3) * 7) + (idxY+3))) * 7 + idxZ + 3;
interactionsPosition[(idxLevel-1)*343 + position] = FSpherical(relativePosition);
}
}
}
}
boxWidthAtLevel *= FReal(2.0);
}
}
}
/** To get the right spherical object from level and child position */
FSpherical getSphericalChild(const int idxLevel, const int position) const {
return childrenPosition[(idxLevel-1)*8 + position];
}
/** To get the right spherical object from level and interaction position */
FSpherical getSphericalInteraction(const int idxLevel, const int position) const {
return interactionsPosition[(idxLevel-1)*343 + position];
}
/** Return the position of a leaf from its tree coordinate */
FPoint getLeafCenter(const FTreeCoordinate coordinate) const {
return FPoint(
FReal(coordinate.getX()) * widthAtLeafLevel + widthAtLeafLevelDiv2 + boxCorner.getX(),
FReal(coordinate.getY()) * widthAtLeafLevel + widthAtLeafLevelDiv2 + boxCorner.getY(),
FReal(coordinate.getZ()) * widthAtLeafLevel + widthAtLeafLevelDiv2 + boxCorner.getZ());
}
public:
......@@ -422,7 +495,10 @@ public:
boxCorner(inBoxCenter.getX()-(inBoxWidth/2),inBoxCenter.getY()-(inBoxWidth/2),inBoxCenter.getZ()-(inBoxWidth/2))
{
// simply does the precomputation
precomputeFactorials();
preComputePosition();
precomputeTranslationCoef();
precomputeRotationVectors();
}
/** Default destructor */
......@@ -464,13 +540,16 @@ public:
computeLegendre(legendre, sph.getCosTheta(), sph.getSinTheta());
// w{l,m}(q,a) = q a^l/(l+|m|)! P{l,m}(cos(alpha)) exp(-i m Beta)
FReal q_aPowL = q; // To consutrct q*a^l continously
int index_l_m = 0; // To construct the index of (l,m) continously
for(int l = 0 ; l <= P ; ++l ){
for(int m = 0 ; m <= l ; ++m ){
const FReal magnitude = q * (FMath::pow( a , l )/fact(l+m))
* legendre[atLm(l,m)];
w[atLm(l,m)].incReal(magnitude * FMath::Cos(FReal(m) * sph.getPhi() + i_pow_m[m & 0x3]));
w[atLm(l,m)].incImag(magnitude * FMath::Sin(FReal(m) * sph.getPhi() + i_pow_m[m & 0x3]));
FReal fm = 0.0; // To have "m" has a float
for(int m = 0 ; m <= l ; ++m, ++index_l_m, ++fm){
const FReal magnitude = q_aPowL * legendre[index_l_m] / factorials[l+m];
w[index_l_m].incReal(magnitude * FMath::Cos(fm * sph.getPhi() + i_pow_m[m & 0x3]));
w[index_l_m].incImag(magnitude * FMath::Sin(fm * sph.getPhi() + i_pow_m[m & 0x3]));
}
q_aPowL *= a;
}
// Goto next particle
......@@ -490,6 +569,8 @@ public:
* and finaly rotate back.
*/
void M2M(CellClass* const FRestrict inPole, const CellClass*const FRestrict *const FRestrict inChildren, const int inLevel) {
// Get the translation coef for this level (same for all chidl)
const FReal*const coef = M2MTranslationCoef[inLevel];
// A buffer to copy the source w allocated once
FComplexe source_w[SizeArray];
// For all children
......@@ -501,27 +582,32 @@ public:
// rotate it forward
const FSpherical sph = getSphericalChild(inLevel, idxChild);
ComplexeArrayMulEqual(source_w,rotationExpMinusImPhi[idxChild],SizeArray);
rotateMultipole(source_w, sph.getAzimuth(), sph.getInclination());
const FReal b = -sph.getR();
//const FReal b = -sph.getR();
// Translate it
FComplexe target_w[SizeArray];
int index_lm = 0;
for(int l = 0 ; l <= P ; ++l ){
for(int m = 0 ; m <= l ; ++m ){
for(int m = 0 ; m <= l ; ++m, ++index_lm ){
// w{l,m}(a+b) = sum(j=m:l, b^(l-j)/(l-j)! w{j,m}(a)
FReal w_lm_real = 0.0;
FReal w_lm_imag = 0.0;
for(int j = m ; j <= l ; ++j ){
const FReal coef = FMath::pow(b,l-j) / fact(l-j);
w_lm_real += coef * source_w[atLm(j,m)].getReal();
w_lm_imag += coef * source_w[atLm(j,m)].getImag();
int index_jm = atLm(m,m); // get atLm(l,m)
int index_l_minus_j = l-m; // get l-j continously
for(int j = m ; j <= l ; ++j, --index_l_minus_j, index_jm += j ){
//const coef = (b^l-j) / (l-j)!;
w_lm_real += coef[index_l_minus_j] * source_w[index_jm].getReal();
w_lm_imag += coef[index_l_minus_j] * source_w[index_jm].getImag();
}
target_w[atLm(l,m)].setRealImag(w_lm_real,w_lm_imag);
target_w[index_lm].setRealImag(w_lm_real,w_lm_imag);
}
}
// Rotate it back
deRotateMultipole(target_w, sph.getAzimuth(), sph.getInclination());
ComplexeArrayMulEqual(target_w,rotationExpImPhi[idxChild],SizeArray);
// Sum the result
FMemUtils::addall( inPole->getMultipole(), target_w, SizeArray);
......@@ -534,7 +620,7 @@ public:
* Implementation of rotation-based operators for Fast Multipole Method in X10
* At page 5 .1 as the operator B
* \f[
* M_{l,m}(a-b') = \sum_{j=|m|}^{\infty}{ \frac{ (j+l)! } { b^{j+l+1} } O_{j,-m}(a) } , \textrm{j bounded by P}
* M_{l,m}(a-b') = \sum_{j=|m|}^{\infty}{ \frac{ (j+l)! } { b^{j+l+1} } O_{j,-m}(a) } , \textrm{j bounded by P-l}
* \f]
* As describe in the paper, when need first to rotate the SH
* then transfer using the formula
......@@ -547,6 +633,7 @@ public:
for(int idxNeigh = 0 ; idxNeigh < 343 ; ++idxNeigh){
// if interaction exits
if(inInteractions[idxNeigh]){
const FReal*const coef = M2LTranslationCoef[inLevel][idxNeigh];