Commit f2b4213b authored by messner's avatar messner

added force computation to Chebyshev FMM


git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/scalfmm/scalfmm/trunk@406 2616d619-271b-44dc-8df4-d4a8f33a7222
parent b80ce6cc
......@@ -38,7 +38,7 @@ class FChebInterpolator : FNoCopyable
* S which is precomputed and reused for all M2M and L2L operations, ie for
* all non leaf inter/anterpolations.
*/
void initChildParentInterpolator()
void initM2MandL2L()
{
F3DPosition ParentRoots[nnodes], ChildRoots[nnodes];
const FReal ParentWidth(2.);
......@@ -82,7 +82,7 @@ public:
// initialize interpolation operator for non M2M and L2L (non leaf
// operations)
this -> initChildParentInterpolator();
this -> initM2MandL2L();
}
......@@ -139,8 +139,8 @@ public:
/**
* The anterpolation corresponds to the P2M and M2M operation. It is indeed
* the transposed interpolation.
* Particle to moment: application of \f$S_\ell(y,\bar y_n)\f$
* (anterpolation, it is the transposed interpolation)
*/
template <class ContainerClass>
void applyP2M(const F3DPosition& center,
......@@ -153,9 +153,6 @@ public:
F3DPosition localPosition;
FReal T_of_x[ORDER][3];
// std::cout << "ClusterCenter = " << center
// << ", width = " << width << std::endl;
// set all multipole expansions to zero
for (unsigned int n=0; n<nnodes; ++n) multipoleExpansion[n] = FReal(0.);
......@@ -165,9 +162,6 @@ public:
// map global position to [-1,1]
map(iter.data().getPosition(), localPosition);
// std::cout << "\tParticlePosition = " << iter.data().getPosition()
// << ", local Position = " << localPosition << std::endl;
// get source value
const FReal sourceValue = iter.data().getPhysicalValue();
......@@ -200,7 +194,7 @@ public:
/**
* The interpolation corresponds to the L2L and L2P operation.
* Local to particle operation: application of \f$S_\ell(x,\bar x_m)\f$ (interpolation)
*/
template <class ContainerClass>
void applyL2P(const F3DPosition& center,
......@@ -248,6 +242,79 @@ public:
}
}
/**
* Local to particle operation: application of \f$\nabla_x S_\ell(x,\bar x_m)\f$ (interpolation)
*/
template <class ContainerClass>
void applyL2PGradient(const F3DPosition& center,
const FReal width,
const FReal *const localExpansion,
ContainerClass *const localParticles) const
{
// setup local to global mapping
const map_glob_loc map(center, width);
F3DPosition Jacobian;
map.computeJacobian(Jacobian);
const FReal jacobian[3] = {Jacobian.getX(), Jacobian.getY(), Jacobian.getZ()};
F3DPosition localPosition;
FReal T_of_x[ORDER][3];
FReal U_of_x[ORDER][3];
typename ContainerClass::BasicIterator iter(*localParticles);
while(iter.hasNotFinished()){
// map global position to [-1,1]
map(iter.data().getPosition(), localPosition);
// get target value
FReal forces[3] = {iter.data().getForces().getX(),
iter.data().getForces().getY(),
iter.data().getForces().getZ()};
// evaluate chebyshev polynomials of source particle
for (unsigned int o=1; o<ORDER; ++o) {
// T_o(x_i)
T_of_x[o][0] = BasisType::T(o, localPosition.getX());
T_of_x[o][1] = BasisType::T(o, localPosition.getY());
T_of_x[o][2] = BasisType::T(o, localPosition.getZ());
// T_o(x_i)
U_of_x[o][0] = BasisType::U(o, localPosition.getX());
U_of_x[o][1] = BasisType::U(o, localPosition.getY());
U_of_x[o][2] = BasisType::U(o, localPosition.getZ());
}
// apply P
for (unsigned int n=0; n<nnodes; ++n) {
for (unsigned int i=0; i<3; ++i) {
FReal P = FReal(1.);
for (unsigned int d=0; d<3; ++d) {
const unsigned int j = node_ids[n][d];
FReal P_d;
if (d==i) {
P_d = 0.;
for (unsigned int o=1; o<ORDER; ++o)
P_d += 2. / ORDER * U_of_x[o][d] * T_of_roots[o][j] * jacobian[d];
} else {
P_d = 1. / ORDER;
for (unsigned int o=1; o<ORDER; ++o)
P_d += 2. / ORDER * T_of_x[o][d] * T_of_roots[o][j];
}
P *= P_d;
}
// the minus sign comes due to the \f$- \nabla_x K(x,y) = \nabla_y K(x,y) = \bar K(x,y)\f$
forces[i] -= P * localExpansion[n];
}
}
iter.data().setForces(forces[0] * iter.data().getPhysicalValue(),
forces[1] * iter.data().getPhysicalValue(),
forces[2] * iter.data().getPhysicalValue());
// increment iterator
iter.gotoNext();
}
}
void applyM2M(const unsigned int ChildIndex,
const FReal *const ChildExpansion,
......
......@@ -197,12 +197,17 @@ public:
M2LHandler->applyU(LeafCell->getLocal() + nnodes,
const_cast<CellClass *const>(LeafCell)->getLocal());
// 2) apply Sx
// 2.a) apply Sx
const F3DPosition LeafCellCenter(getLeafCellCenter(LeafCell->getCoordinate()));
Interpolator->applyL2P(LeafCellCenter,
BoxWidthLeaf,
LeafCell->getLocal(),
TargetParticles);
// 2.b) apply Px (grad Sx)
Interpolator->applyL2PGradient(LeafCellCenter,
BoxWidthLeaf,
LeafCell->getLocal(),
TargetParticles);
}
......@@ -216,15 +221,23 @@ public:
typename ContainerClass::BasicIterator iTargets(*TargetParticles);
while (iTargets.hasNotFinished()) {
ParticleClass& Target = iTargets.data();
const FReal wt = Target.getPhysicalValue();
{ // loop: source particles (target leaf cell == source leaf cell)
typename ContainerClass::ConstBasicIterator iSources(*SourceParticles);
while (iSources.hasNotFinished()) {
const ParticleClass& Source = iSources.data();
// only if target and source are not identical
if (&Target != &Source)
Target.incPotential(MatrixKernel->evaluate(Target.getPosition(), Source.getPosition())
* Source.getPhysicalValue());
if (&Target != &Source) {
const FReal one_over_r = MatrixKernel->evaluate(Target.getPosition(), Source.getPosition());
const FReal ws = Source.getPhysicalValue();
// potential
Target.incPotential(one_over_r * ws);
F3DPosition force(Target.getPosition() - Source.getPosition());
force *= ((ws*wt) * (one_over_r*one_over_r*one_over_r));
// force
Target.incForces(force.getX(), force.getY(), force.getZ());
}
// progress sources
iSources.gotoNext();
}
......@@ -237,8 +250,14 @@ public:
while (iSources.hasNotFinished()) {
const ParticleClass& Source = iSources.data();
// target and source cannot be identical
Target.incPotential(MatrixKernel->evaluate(Target.getPosition(), Source.getPosition())
* Source.getPhysicalValue());
const FReal one_over_r = MatrixKernel->evaluate(Target.getPosition(), Source.getPosition());
const FReal ws = Source.getPhysicalValue();
// potential
Target.incPotential(one_over_r * ws);
F3DPosition force(Target.getPosition() - Source.getPosition());
force *= ((ws*wt) * (one_over_r*one_over_r*one_over_r));
// force
Target.incForces(force.getX(), force.getY(), force.getZ());
// progress sources
iSources.gotoNext();
}
......
......@@ -7,6 +7,7 @@
#include "../Extensions/FExtendPosition.hpp"
#include "../Extensions/FExtendPhysicalValue.hpp"
#include "../Extensions/FExtendPotential.hpp"
#include "../Extensions/FExtendForces.hpp"
/**
* @author Matthias Messner (matthias.matthias@inria.fr)
......@@ -18,7 +19,8 @@
*/
class FChebParticle : public FExtendPosition,
public FExtendPhysicalValue,
public FExtendPotential
public FExtendPotential,
public FExtendForces
{
public:
~FChebParticle() {}
......
......@@ -74,7 +74,7 @@ const FReal computeINFnorm(unsigned int N, FReal *const u, FReal *const v)
// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
const unsigned int ORDER = 3;
const unsigned int ORDER = 4;
const FReal epsilon = FParameters::getValue(argc, argv, "-eps", FReal(1e-3));
const long NbPart = FParameters::getValue(argc, argv, "-num", 100000);
const unsigned int TreeHeight = FParameters::getValue(argc, argv, "-h", 5);
......@@ -91,16 +91,16 @@ int main(int argc, char* argv[])
typedef FChebParticle ParticleClass;
typedef FVector<FChebParticle> ContainerClass;
typedef FChebLeaf<ParticleClass,ContainerClass> LeafClass;
typedef FChebMatrixKernelRR MatrixKernelClass;
typedef FChebMatrixKernelR MatrixKernelClass;
typedef FChebCell<ORDER> CellClass;
typedef FOctree<ParticleClass,CellClass,ContainerClass,LeafClass> OctreeClass;
typedef FChebKernels<ParticleClass,CellClass,ContainerClass,MatrixKernelClass,ORDER> KernelClass;
//typedef FFmmAlgorithm<OctreeClass,ParticleClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
typedef FFmmAlgorithmThread<OctreeClass,ParticleClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
typedef FFmmAlgorithm<OctreeClass,ParticleClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
//typedef FFmmAlgorithmThread<OctreeClass,ParticleClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
// What we do //////////////////////////////////////////////////////
std::cout << ">> Testing the Chebyshev interpolation base FMM algorithm.\n";
const F3DPosition BoxCenter(.5*Width, .5*Width, .5*Width);
OctreeClass tree(TreeHeight, SubTreeHeight, Width, BoxCenter);
......@@ -135,8 +135,12 @@ int main(int argc, char* argv[])
const ContainerClass *const Targets = iLeafs.getCurrentListTargets();
const unsigned int NumTargets = Targets->getSize();
FReal* Potential = new FReal [NumTargets];
FBlas::scal(NumTargets, FReal(0.), Potential);
FBlas::setzero(NumTargets, Potential);
FReal* Force = new FReal [NumTargets * 3];
FBlas::setzero(NumTargets * 3, Force);
std::cout << "\nDirect computation of " << NumTargets << " target particles ..." << std::endl;
const MatrixKernelClass MatrixKernel;
......@@ -145,12 +149,23 @@ int main(int argc, char* argv[])
unsigned int counter = 0;
ContainerClass::ConstBasicIterator iTarget(*Targets);
while(iTarget.hasNotFinished()) {
const FReal wt = iTarget.data().getPhysicalValue();
ContainerClass::ConstBasicIterator iSource(*Sources);
while(iSource.hasNotFinished()) {
if (&iTarget.data() != &iSource.data())
Potential[counter] += MatrixKernel.evaluate(iTarget.data().getPosition(),
iSource.data().getPosition())
* iSource.data().getPhysicalValue();
if (&iTarget.data() != &iSource.data()) {
const FReal one_over_r = MatrixKernel.evaluate(iTarget.data().getPosition(),
iSource.data().getPosition());
const FReal ws = iSource.data().getPhysicalValue();
// potential
Potential[counter] += one_over_r * ws;
// force
F3DPosition force(iTarget.data().getPosition() - iSource.data().getPosition());
force *= ((ws*wt) * (one_over_r*one_over_r*one_over_r));
// force
Force[counter*3 + 0] += force.getX();
Force[counter*3 + 1] += force.getY();
Force[counter*3 + 2] += force.getZ();
}
iSource.gotoNext();
}
counter++;
......@@ -160,26 +175,37 @@ int main(int argc, char* argv[])
FReal* ApproxPotential = new FReal [NumTargets];
FReal* ApproxForce = new FReal [NumTargets * 3];
unsigned int counter = 0;
ContainerClass::ConstBasicIterator iTarget(*Targets);
while(iTarget.hasNotFinished()) {
ApproxPotential[counter] = iTarget.data().getPotential();
//std::cout << Potential[counter] << " - " << ApproxPotential[counter] << " = "
// << Potential[counter]-ApproxPotential[counter] << "\t rel error = "
// << (Potential[counter]-ApproxPotential[counter]) / Potential[counter]
// << std::endl;
ApproxPotential[counter] = iTarget.data().getPotential();
ApproxForce[counter*3 + 0] = iTarget.data().getForces().getX();
ApproxForce[counter*3 + 1] = iTarget.data().getForces().getY();
ApproxForce[counter*3 + 2] = iTarget.data().getForces().getZ();
counter++;
iTarget.gotoNext();
}
std::cout << "\nRelative L2 error = " << computeL2norm( NumTargets, Potential, ApproxPotential)
std::cout << "\nPotential error:" << std::endl;
std::cout << "Relative L2 error = " << computeL2norm( NumTargets, Potential, ApproxPotential)
<< std::endl;
std::cout << "Relative Lmax error = " << computeINFnorm(NumTargets, Potential, ApproxPotential)
<< std::endl;
std::cout << "\nForce error:" << std::endl;
std::cout << "Relative L2 error = " << computeL2norm( NumTargets*3, Force, ApproxForce)
<< std::endl;
std::cout << "Relative Lmax error = " << computeINFnorm(NumTargets*3, Force, ApproxForce)
<< std::endl;
std::cout << "Relative Lmax error = " << computeINFnorm(NumTargets, Potential, ApproxPotential)
<< "\n" << std::endl;
std::cout << std::endl;
// free memory
delete [] Potential;
delete [] ApproxPotential;
delete [] Force;
delete [] ApproxForce;
/*
......
......@@ -49,6 +49,7 @@ const FReal computeINFnorm(unsigned int N, FReal *const u, FReal *const v)
/**
* In this file we show how to use octree
*/
......@@ -73,15 +74,14 @@ int main(int, char **){
// Leaf size
FReal width = 1.;
FReal width = 3.723;
////////////////////////////////////////////////////////////////////
LeafClass X;
F3DPosition cx(0., 0., 0.);
const long M = 10000;
std::cout << "Fill the leaf X of width " << width
<< " centered at cx=[" << cx.getX() << "," << cx.getY() << "," << cx.getZ()
<< "] with M=" << M << " target particles" << std::endl;
<< " centered at cx=" << cx << " with M=" << M << " target particles" << std::endl;
{
FChebParticle particle;
for(long i=0; i<M; ++i){
......@@ -89,6 +89,7 @@ int main(int, char **){
FReal y = (FReal(rand())/FRandMax - FReal(.5)) * width + cx.getY();
FReal z = (FReal(rand())/FRandMax - FReal(.5)) * width + cx.getZ();
particle.setPosition(x, y, z);
particle.setPhysicalValue(FReal(rand())/FRandMax);
X.push(particle);
}
}
......@@ -96,11 +97,10 @@ int main(int, char **){
////////////////////////////////////////////////////////////////////
LeafClass Y;
F3DPosition cy(2., 0., 0.);
F3DPosition cy(FReal(2.)*width, 0., 0.);
const long N = 10000;
std::cout << "Fill the leaf Y of width" << width
<< " centered at cy=[" << cy.getX() << "," << cy.getY() << "," << cy.getZ()
<< "] with N=" << N << " target particles" << std::endl;
std::cout << "Fill the leaf Y of width " << width
<< " centered at cy=" << cy << " with N=" << N << " target particles" << std::endl;
{
FChebParticle particle;
for(long i=0; i<N; ++i){
......@@ -117,12 +117,13 @@ int main(int, char **){
////////////////////////////////////////////////////////////////////
// approximative computation
const unsigned int ORDER = 5;
const unsigned int ORDER = 8;
const unsigned int nnodes = TensorTraits<ORDER>::nnodes;
typedef FChebInterpolator<ORDER> InterpolatorClass;
InterpolatorClass S;
std::cout << "\nCompute interactions approximatively, interpolation order = " << ORDER << " ..." << std::endl;
std::cout << "\nCompute interactions approximatively, interpolation order = " << ORDER << " ..."
<< std::endl;
time.tic();
// Anterpolate: W_n = \sum_j^N S(y_j,\bar y_n) * w_j
......@@ -141,9 +142,13 @@ int main(int, char **){
F[i] += MatrixKernel.evaluate(rootsX[i], rootsY[j]) * W[j];
}
// Interpolate f_i = \sum_m^L S(x_i,\bar x_m) * F_m
// Interpolate p_i = \sum_m^L S(x_i,\bar x_m) * F_m
S.applyL2P(cx, width, F, X.getTargets());
// Interpolate f_i = \sum_m^L P(x_i,\bar x_m) * F_m
S.applyL2PGradient(cx, width, F, X.getTargets());
time.tac();
std::cout << "Done in " << time.elapsed() << "sec." << std::endl;
// -----------------------------------------------------
......@@ -153,22 +158,44 @@ int main(int, char **){
std::cout << "Compute interactions directly ..." << std::endl;
time.tic();
FReal* approx_f = new FReal[M];
FReal* f = new FReal[M];
for (unsigned int i=0; i<M; ++i) f[i] = FReal(0.);
ContainerClass::ConstBasicIterator iterY(*(Y.getSrc()));
while(iterY.hasNotFinished()){
const F3DPosition& y = iterY.data().getPosition();
const FReal w = iterY.data().getPhysicalValue();
FReal* approx_f = new FReal [M * 3];
FReal* f = new FReal [M * 3];
FBlas::setzero(M*3, f);
FReal* approx_p = new FReal[M];
FReal* p = new FReal[M];
FBlas::setzero(M, p);
{ // start direct computation
unsigned int counter = 0;
ContainerClass::ConstBasicIterator iterX(*(X.getSrc()));
while(iterX.hasNotFinished()){
const F3DPosition& x = iterX.data().getPosition();
f[counter++] += MatrixKernel.evaluate(x,y) * w;
const FReal wx = iterX.data().getPhysicalValue();
ContainerClass::ConstBasicIterator iterY(*(Y.getSrc()));
while(iterY.hasNotFinished()){
const F3DPosition& y = iterY.data().getPosition();
const FReal wy = iterY.data().getPhysicalValue();
const FReal one_over_r = MatrixKernel.evaluate(x, y);
// potential
p[counter] += one_over_r * wy;
// force
F3DPosition force(x - y);
force *= one_over_r*one_over_r*one_over_r;
f[counter*3 + 0] += force.getX() * wx * wy;
f[counter*3 + 1] += force.getY() * wx * wy;
f[counter*3 + 2] += force.getZ() * wx * wy;
iterY.gotoNext();
}
counter++;
iterX.gotoNext();
}
iterY.gotoNext();
}
} // end direct computation
time.tac();
std::cout << "Done in " << time.elapsed() << "sec." << std::endl;
......@@ -177,14 +204,28 @@ int main(int, char **){
ContainerClass::ConstBasicIterator iterX(*(X.getSrc()));
unsigned int counter = 0;
while(iterX.hasNotFinished()) {
approx_f[counter++] = iterX.data().getPotential();
approx_p[counter] = iterX.data().getPotential();
const F3DPosition& force = iterX.data().getForces();
approx_f[counter*3 + 0] = force.getX();
approx_f[counter*3 + 1] = force.getY();
approx_f[counter*3 + 2] = force.getZ();
counter++;
iterX.gotoNext();
}
std::cout << "\nRelative L2 error = " << computeL2norm( M, f, approx_f) << std::endl;
std::cout << "Relative Lmax error = " << computeINFnorm(M, f, approx_f) << "\n" << std::endl;
std::cout << "\nPotential error:" << std::endl;
std::cout << "Relative L2 error = " << computeL2norm( M, p, approx_p) << std::endl;
std::cout << "Relative Lmax error = " << computeINFnorm(M, p, approx_p) << std::endl;
std::cout << "\nForce error:" << std::endl;
std::cout << "Relative L2 error = " << computeL2norm( M*3, f, approx_f) << std::endl;
std::cout << "Relative Lmax error = " << computeINFnorm(M*3, f, approx_f) << std::endl;
std::cout << std::endl;
// free memory
delete [] approx_p;
delete [] p;
delete [] approx_f;
delete [] f;
......
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