testChebSxUCBSy.cpp 7.78 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
// ===================================================================================
// Ce LOGICIEL "ScalFmm" est couvert par le copyright Inria 20xx-2012.
// Inria détient tous les droits de propriété sur le LOGICIEL, et souhaite que
// la communauté scientifique l'utilise afin de le tester et de l'évaluer.
// Inria donne gracieusement le droit d'utiliser ce LOGICIEL. Toute utilisation
// dans un but lucratif ou à des fins commerciales est interdite sauf autorisation
// expresse et préalable d'Inria.
// Toute utilisation hors des limites précisées ci-dessus et réalisée sans l'accord
// expresse préalable d'Inria constituerait donc le délit de contrefaçon.
// Le LOGICIEL étant un produit en cours de développement, Inria ne saurait assurer
// aucune responsabilité et notamment en aucune manière et en aucun cas, être tenu
// de répondre d'éventuels dommages directs ou indirects subits par l'utilisateur.
// Tout utilisateur du LOGICIEL s'engage à communiquer à Inria ses remarques
// relatives à l'usage du LOGICIEL
// ===================================================================================

// ==== CMAKE =====
// @FUSE_BLAS
// ================
20 21 22 23 24 25 26

#include <iostream>

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

27 28
#include "../../Src/Utils/FTic.hpp"
#include "../../Src/Utils/FMath.hpp"
29
#include "../../Src/Utils/FParameters.hpp"
30

31
#include "../../Src/Containers/FVector.hpp"
32

33
#include "../../Src/Utils/FAssert.hpp"
34
#include "../../Src/Utils/FPoint.hpp"
35

36

37 38
#include "../../Src/Kernels/Chebyshev/FChebInterpolator.hpp"
#include "../../Src/Kernels/Chebyshev/FChebMatrixKernel.hpp"
39

40 41
#include "../../Src/Components/FSimpleLeaf.hpp"
#include "../../Src/Kernels/P2P/FP2PParticleContainer.hpp"
42

43 44
void applyM2M(FReal *const S,	FReal *const w, const unsigned int n,	FReal *const W, const unsigned int N)
{ FBlas::gemtva(n, N, FReal(1.), S,	w, W); }
45 46 47 48 49 50 51 52


/**
* In this file we show how to use octree
*/

int main(int argc, char* argv[])
{
53 54
    typedef FP2PParticleContainer ContainerClass;
    typedef FSimpleLeaf<ContainerClass> LeafClass;
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
	typedef FChebMatrixKernelR MatrixKernelClass;

	///////////////////////What we do/////////////////////////////
	std::cout << "\nTask: Compute interactions between source particles in leaf Y and target\n";
	std::cout << " particles in leaf X. Compare the fast summation K ~ Sx K Sy' with the\n";
	std::cout << " direct computation.\n" << std::endl;
	//////////////////////////////////////////////////////////////

	MatrixKernelClass MatrixKernel;
	const FReal FRandMax = FReal(RAND_MAX);
	FTic time;

	
	// Leaf size
	FReal width(FReal(rand()) / FRandMax * FReal(10.));

	////////////////////////////////////////////////////////////////////
	LeafClass X;
73
	FPoint cx(0., 0., 0.);
74
	const long M = 1000;
75 76 77
	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;
78
    {
79 80 81
		for(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();
82 83
            FReal z = (FReal(rand())/FRandMax - FReal(.5)) * width + cx.getZ();
            X.push(FPoint(x,y,z));
84 85 86 87 88 89
		}
	}


	////////////////////////////////////////////////////////////////////
	LeafClass Y;
90
	FPoint cy(FReal(2.)*width, 0., 0.);
91
	const long N = 1000;
92 93 94
	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;
95
    {
96 97 98 99
		for(long i=0; i<N; ++i){
			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();
100
            Y.push(FPoint(x, y, z),FReal(rand())/FRandMax);
101 102 103 104 105 106
		}
	}


	////////////////////////////////////////////////////////////////////
	// approximative computation
107
	const unsigned int ORDER = 10;
108 109 110
	const unsigned int nnodes = TensorTraits<ORDER>::nnodes;
	typedef FChebInterpolator<ORDER> InterpolatorClass;
	InterpolatorClass S;
111
	
112

113 114 115 116
	std::cout << "\nCompute interactions approximatively, ORDER = " << ORDER << std::endl;
	FReal W[nnodes]; // multipole expansion
	FReal F[nnodes]; // local expansion
	for (unsigned int n=0; n<nnodes; ++n)
117 118 119
		W[n] = F[n] = FReal(0.);

	// Anterpolate: W_n = \sum_j^N S(y_j,\bar y_n) * w_j
120
	time.tic();
121
	S.applyP2M(cy, width, W, Y.getSrc()); // the multipole expansions are set to 0 in S.applyP2M
122
	std::cout << "P2M done in " << time.tacAndElapsed() << "s" << std::endl;
123 124

	// M2L (direct)
125
	FPoint rootsX[nnodes], rootsY[nnodes];
126 127
	FChebTensor<ORDER>::setRoots(cx, width, rootsX);
	FChebTensor<ORDER>::setRoots(cy, width, rootsY);
128 129 130 131 132 133 134

	{
		for (unsigned int i=0; i<nnodes; ++i) {
			F[i] = FReal(0.);
			for (unsigned int j=0; j<nnodes; ++j)
				F[i] += MatrixKernel.evaluate(rootsX[i], rootsY[j]) * W[j];
		}
135 136
	}

137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152
//	{
//		for (unsigned int ix=0; ix<ORDER; ++ix)
//			for (unsigned int jx=0; jx<ORDER; ++jx)
//				for (unsigned int kx=0; kx<ORDER; ++kx)  {
//					const unsigned int idx = kx*ORDER*ORDER + jx*ORDER + ix;
//					F[idx] = FReal(0.);
//					for (unsigned int iy=0; iy<ORDER; ++iy)
//						for (unsigned int jy=0; jy<ORDER; ++jy)
//							for (unsigned int ky=0; ky<ORDER; ++ky) {
//								const unsigned int idy = ky*ORDER*ORDER + jy*ORDER + iy;
//								F[idx] += MatrixKernel.evaluate(rootsX[idx], rootsY[idy]) * W[idy];
//							}
//				}
//	}
		

153
	// Interpolate f_i = \sum_m^L S(x_i,\bar x_m) * F_m
154 155
	time.tic();
	//S.applyL2PTotal(cx, width, F, X.getTargets());
156
	S.applyL2PTotal(cx, width, F, X.getTargets());
157
	std::cout << "L2P done in " << time.tacAndElapsed() << "s" << std::endl;
158 159 160 161 162 163 164 165 166 167

	// -----------------------------------------------------

	////////////////////////////////////////////////////////////////////
	// direct computation
	std::cout << "Compute interactions directly ..." << std::endl;
	time.tic();

	FReal* approx_f = new FReal[M];
	FReal*        f = new FReal[M];
168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189
    {
        for (unsigned int i=0; i<M; ++i) f[i] = FReal(0.);

        const FReal*const positionsX = Y.getSrc()->getPositions()[0];
        const FReal*const positionsY = Y.getSrc()->getPositions()[1];
        const FReal*const positionsZ = Y.getSrc()->getPositions()[2];
        const FReal*const physicalValues = Y.getSrc()->getPhysicalValues();

        for(int idxPart = 0 ; idxPart < Y.getSrc()->getNbParticles() ; ++idxPart){
            const FPoint y = FPoint(positionsX[idxPart],positionsY[idxPart],positionsZ[idxPart]);
            const FReal        w = physicalValues[idxPart];

            const FReal*const xpositionsX = X.getSrc()->getPositions()[0];
            const FReal*const xpositionsY = X.getSrc()->getPositions()[1];
            const FReal*const xpositionsZ = X.getSrc()->getPositions()[2];

            for(int idxPartX = 0 ; idxPartX < X.getSrc()->getNbParticles() ; ++idxPartX){
                const FPoint x = FPoint(xpositionsX[idxPart],xpositionsY[idxPart],xpositionsZ[idxPart]);
                f[idxPartX] += MatrixKernel.evaluate(x,y) * w;
            }
        }
    }
190 191 192 193 194
	time.tac();
	std::cout << "Done in " << time.elapsed() << "sec." << std::endl;


	////////////////////////////////////////////////////////////////////
195 196 197
    const FReal*const potentials = X.getSrc()->getPotentials();
    for(int idxPart = 0 ; idxPart < X.getSrc()->getNbParticles() ; ++idxPart){
        approx_f[idxPart] = potentials[idxPart];
198 199
	}

200
	
201 202
	//for (unsigned int i=0; i<8; ++i)
	//	std::cout << f[i] << "\t" << approx_f[i] << "\t" << approx_f[i]/f[i] << std::endl;
203 204
	

205 206
    std::cout << "\nRelative L2 error  = " << FMath::FAccurater( f, approx_f, M) << std::endl;
    std::cout << "Relative Lmax error = "  << FMath::FAccurater( f, approx_f, M) << "\n" << std::endl;
207 208 209 210 211 212 213 214 215 216 217

	// free memory
	delete [] approx_f;
	delete [] f;


	return 0;
}


// [--END--]