Commit 86ca0e6c authored by COULAUD Olivier's avatar COULAUD Olivier

Taylor : Amelioration du L2P mais il reste un bug

parent 8cee40f2
...@@ -607,7 +607,7 @@ public: ...@@ -607,7 +607,7 @@ public:
* where \f$x_c\f$ is the centre of the cell and \f$x_j\f$ the \f$j^{th}\f$ particles and \f$q_j\f$ its charge. * where \f$x_c\f$ is the centre of the cell and \f$x_j\f$ the \f$j^{th}\f$ particles and \f$q_j\f$ its charge.
*/ */
void L2P(const CellClass* const local, void L2P(const CellClass* const local,
ContainerClass* const particles) ContainerClass* const particles)
{ {
FPoint locCenter = getLeafCenter(local->getCoordinate()); FPoint locCenter = getLeafCenter(local->getCoordinate());
...@@ -631,13 +631,16 @@ public: ...@@ -631,13 +631,16 @@ public:
//Iteration over particles //Iteration over particles
for(int i=0 ; i<nbPart ; ++i){ for(int i=0 ; i<nbPart ; ++i){
//
FReal dx = posX[i]-locCenter.getX(); FReal dx = posX[i] - locCenter.getX();
FReal dy = posY[i]-locCenter.getY(); FReal dy = posY[i] - locCenter.getY();
FReal dz = posZ[i]-locCenter.getZ(); FReal dz = posZ[i] - locCenter.getZ();
printf("dx = %f, dy = %f, dz = %f\n",dx,dy,dz); printf("dx = %f, dy = %f, dz = %f\n",dx,dy,dz);
int a=0,b=0,c=0; int a=0,b=0,c=0;
forceX[i] = 0.0 ;
forceY[i] = 0.0 ;
forceZ[i] = 0.0 ;
// //
// Precompute an arrays of Array[i] = dx^(i-1) // Precompute an arrays of Array[i] = dx^(i-1)
arrayDX[0] = 0.0 ; arrayDX[0] = 0.0 ;
...@@ -655,19 +658,24 @@ public: ...@@ -655,19 +658,24 @@ public:
//Iteration over Local array //Iteration over Local array
const FReal * iterLocal = local->getLocal(); const FReal * iterLocal = local->getLocal();
FReal partPhyValue = phyValues[i];
//
FReal locPot = 0.0, locForceX = 0.0, locForceY = 0.0, locForceZ = 0.0 ;
for(int j=0 ; j<SizeVector ; ++j){ for(int j=0 ; j<SizeVector ; ++j){
FReal partPhyValue = phyValues[i];
FReal locForce = iterLocal[j]; FReal locForce = iterLocal[j];
// compute the potential // compute the potential
targetsPotentials[i] += partPhyValue*iterLocal[j]*arrayDX[a+1]*arrayDY[b+1]*arrayDZ[c+1]; locPot += iterLocal[j]*arrayDX[a+1]*arrayDY[b+1]*arrayDZ[c+1];
//Application of forces //Application of forces
forceX[i] -= partPhyValue*FReal(a)*locForce*arrayDX[a]*arrayDY[b+1]*arrayDZ[c+1]; locForceX += FReal(a)*locForce*arrayDX[a]*arrayDY[b+1]*arrayDZ[c+1];
forceY[i] -= partPhyValue*FReal(b)*locForce*arrayDX[a+1]*arrayDY[b]*arrayDZ[c+1]; locForceY += FReal(b)*locForce*arrayDX[a+1]*arrayDY[b]*arrayDZ[c+1];
forceZ[i] -= partPhyValue*FReal(c)*locForce*arrayDX[a+1]*arrayDY[b+1]*arrayDZ[c]; locForceZ += FReal(c)*locForce*arrayDX[a+1]*arrayDY[b+1]*arrayDZ[c];
// //
incPowers(&a,&b,&c); incPowers(&a,&b,&c);
} }
targetsPotentials[i] = partPhyValue*locPot ;
forceX[i] -= partPhyValue*locForceX ;
forceY[i] -= partPhyValue*locForceY ;
forceZ[i] -= partPhyValue*locForceZ ;
} }
} }
......
...@@ -26,7 +26,7 @@ ...@@ -26,7 +26,7 @@
int main(int argc,char* argv[]){ int main(int argc,char* argv[]){
static const int P = 3; static const int P = 3;
static const int order = 1; static const int order = 1;
FPoint rootCenter(FReal(0),FReal(0),FReal(0)); FPoint rootCenter(FReal(0.0),FReal(0.0),FReal(0.0));
FReal boxWidth = FReal(8); FReal boxWidth = FReal(8);
typedef FTaylorCell<P,order> CellClass; typedef FTaylorCell<P,order> CellClass;
...@@ -76,14 +76,21 @@ int main(int argc,char* argv[]){ ...@@ -76,14 +76,21 @@ int main(int argc,char* argv[]){
for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){ for(int idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
potential += potentials[idxPart]; potential += potentials[idxPart];
fx += forcesX[idxPart]; fx += forcesX[idxPart];
fy += forcesY[idxPart]; fy += forcesY[idxPart];
fz += forcesZ[idxPart]; fz += forcesZ[idxPart];
printf("part : %f, fx = %f, fy = %f, fz = %f\n",leaf->getTargets()->getPositions()[0][0],fx,fy,fz); printf("part : %f, fx = %f, fy = %f, fz = %f\n",leaf->getTargets()->getPositions()[0][0],fx,fy,fz);
} }
}); });
printf("potential : %f\n",potential); //
FReal dx = part1Pos.getX() - part2Pos.getX();
FReal dy = part1Pos.getY() - part2Pos.getY();
FReal dz = part1Pos.getZ() - part2Pos.getZ();
//
FReal potTheo = physVal2*physVal1*FMath::Sqrt(1.0 / (dx*dx + dy*dy + dz*dz));
printf("Exact potential : %f Computed potential : %f Error: %f \n",potTheo, potential,potTheo- potential);
} }
return 0; return 0;
......
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