testChebSxUCBSy.cpp 8.07 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
#include "../../Src/Kernels/Chebyshev/FChebInterpolator.hpp"
38
#include "../../Src/Kernels/Interpolation/FInterpMatrixKernel.hpp"
39

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

43 44
#include "../../Src/Utils/FParameterNames.hpp"

45 46

template <class FReal>
47 48
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); }
49 50 51 52 53 54 55 56


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

int main(int argc, char* argv[])
{
57 58
    FHelpDescribeAndExit(argc, argv, "Test Chebyshev interation computations.");

59 60 61 62
    typedef double FReal;
    typedef FP2PParticleContainer<FReal> ContainerClass;
    typedef FSimpleLeaf<FReal, ContainerClass> LeafClass;
	typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
63 64 65 66 67 68 69

	///////////////////////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;
	//////////////////////////////////////////////////////////////

BRAMAS Berenger's avatar
BRAMAS Berenger committed
70
    MatrixKernelClass MatrixKernel;
71 72 73 74
	FTic time;

	
	// Leaf size
BRAMAS Berenger's avatar
BRAMAS Berenger committed
75
    FReal width(FReal(drand48()) * FReal(10.));
76 77 78

	////////////////////////////////////////////////////////////////////
	LeafClass X;
79
    FPoint<FReal> cx(0., 0., 0.);
80
	const long M = 1000;
81 82 83
	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;
84
    {
85
		for(long i=0; i<M; ++i){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
86 87 88
            FReal x = (FReal(drand48()) - FReal(.5)) * width + cx.getX();
            FReal y = (FReal(drand48()) - FReal(.5)) * width + cx.getY();
            FReal z = (FReal(drand48()) - FReal(.5)) * width + cx.getZ();
89
            X.push(FPoint<FReal>(x,y,z));
90 91 92 93 94 95
		}
	}


	////////////////////////////////////////////////////////////////////
	LeafClass Y;
96
    FPoint<FReal> cy(FReal(2.)*width, 0., 0.);
97
	const long N = 1000;
98 99 100
	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;
101
    {
102
		for(long i=0; i<N; ++i){
BRAMAS Berenger's avatar
BRAMAS Berenger committed
103 104 105
            FReal x = (FReal(drand48()) - FReal(.5)) * width + cy.getX();
            FReal y = (FReal(drand48()) - FReal(.5)) * width + cy.getY();
            FReal z = (FReal(drand48()) - FReal(.5)) * width + cy.getZ();
106
            Y.push(FPoint<FReal>(x, y, z),FReal(drand48()));
107 108 109 110 111 112
		}
	}


	////////////////////////////////////////////////////////////////////
	// approximative computation
113
	const unsigned int ORDER = 10;
114
	const unsigned int nnodes = TensorTraits<ORDER>::nnodes;
115
	typedef FChebInterpolator<FReal,ORDER,MatrixKernelClass> InterpolatorClass;
116
	InterpolatorClass S;
117
	
118

119 120 121 122
	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)
123 124 125
		W[n] = F[n] = FReal(0.);

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

	// M2L (direct)
131 132 133
    FPoint<FReal> rootsX[nnodes], rootsY[nnodes];
	FChebTensor<FReal,ORDER>::setRoots(cx, width, rootsX);
	FChebTensor<FReal,ORDER>::setRoots(cy, width, rootsY);
134 135 136 137 138 139 140

	{
		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];
		}
141 142
	}

143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158
//	{
//		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];
//							}
//				}
//	}
		

159
	// Interpolate f_i = \sum_m^L S(x_i,\bar x_m) * F_m
160 161
	time.tic();
	//S.applyL2PTotal(cx, width, F, X.getTargets());
162
	S.applyL2PTotal(cx, width, F, X.getTargets());
163
	std::cout << "L2P done in " << time.tacAndElapsed() << "s" << std::endl;
164 165 166 167 168 169 170 171 172 173

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

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

	FReal* approx_f = new FReal[M];
	FReal*        f = new FReal[M];
174 175 176 177 178 179 180 181
    {
        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();

182
        for(FSize idxPart = 0 ; idxPart < Y.getSrc()->getNbParticles() ; ++idxPart){
183
            const FPoint<FReal> y = FPoint<FReal>(positionsX[idxPart],positionsY[idxPart],positionsZ[idxPart]);
184 185 186 187 188 189
            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];

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


	////////////////////////////////////////////////////////////////////
201
    const FReal*const potentials = X.getSrc()->getPotentials();
202
    for(FSize idxPart = 0 ; idxPart < X.getSrc()->getNbParticles() ; ++idxPart){
203
        approx_f[idxPart] = potentials[idxPart];
204 205
	}

206
	
207 208
	//for (unsigned int i=0; i<8; ++i)
	//	std::cout << f[i] << "\t" << approx_f[i] << "\t" << approx_f[i]/f[i] << std::endl;
209 210
	

211 212
    std::cout << "\nRelative L2 error  = " << FMath::FAccurater<FReal>( f, approx_f, M) << std::endl;
    std::cout << "Relative Lmax error = "  << FMath::FAccurater<FReal>( f, approx_f, M) << "\n" << std::endl;
213 214 215 216 217 218 219 220 221 222 223

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


	return 0;
}


// [--END--]