main.cpp 6.99 KB
Newer Older
Pietro Incardona's avatar
Pietro Incardona committed
1
2
3
4
5
6
7
8
9
10
11
12
13
/*
 * ### WIKI 1 ###
 *
 * ## Simple example
 *
 * In this example we show 1D PSE derivative function approximation
 *
 * ### WIKI END ###
 *
 */

#include "Vector/vector_dist.hpp"
#include "Decomposition/CartDecomposition.hpp"
14
#include "Kernels.hpp"
Pietro Incardona's avatar
Pietro Incardona committed
15
16
17
#include "data_type/aggregate.hpp"
#include <cmath>

18

Pietro Incardona's avatar
Pietro Incardona committed
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
/*
 * ### WIKI 2 ###
 *
 * Here we define the function x*e^(-x*x) and its
 * second derivative in analytic form
 *
 * 2x*(2*x-3)*e^(-x^2)
 *
 */

double f_xex2(double x)
{
	return x*exp(-x*x);
}

double f_xex2(Point<1,double> & x)
{
	return x.get(0)*exp(-x.get(0)*x.get(0));
}

double Lapf_xex2(Point<1,double> & x)
{
	return 2.0*x.get(0)*(2.0*x.get(0)*x.get(0) - 3.0)*exp(-x.get(0)*x.get(0));
}

/*
 *
 * ### WIKI END ###
 *
 */

int main(int argc, char* argv[])
{
	//
	// ### WIKI 3 ###
	//
	// Some useful parameters. Like
	//
	// * Number of particles
	// * Minimum number of padding particles
	// * The computational domain
	// * The spacing
	// * The mollification length
	// * Second order Laplacian kernel in 1D
	//

	// Number of particles
	const size_t Npart = 1000;

	// Number of padding particles (At least)
	const size_t Npad = 20;

	// The domain
	Box<1,double> box({0.0},{4.0});

	// Calculated spacing
	double spacing = box.getHigh(0) / Npart;

	// Epsilon of the particle kernel
	const double eps = 2*spacing;

	// Laplacian PSE kernel 1 dimension, on double, second order
	Lap<1,double,2> lker(eps);

	//
	// ### WIKI 2 ###
	//
	// Here we Initialize the library and we define Ghost size
	// and non-periodic boundary conditions
	//
	init_global_v_cluster(&argc,&argv);
	Vcluster & v_cl = *global_v_cluster;

    size_t bc[1]={NON_PERIODIC};
	Ghost<1,double> g(12*eps);

	//
	// ### WIKI 3 ###
	//
	// Here we are creating a distributed vector defined by the following parameters
	//
	// we create a set of N+1 particles to have a fully covered domain of particles between 0.0 and 4.0
	// Suppose we have a spacing given by 1.0 you need 4 +1 particles to cover your domain
	//
	vector_dist<1,double, aggregate<double>, CartDecomposition<1,double> > vd(Npart+1,box,bc,g);

	//
	// ### WIKI 4 ###
	//
	// We assign the position to the particles, the scalar property is set to
	// the function x*e^(-x*x) value.
	// Each processor has parts of the particles and fill part of the space, the
	// position is assigned independently from the decomposition.
	//
	// In this case, if there are 1001 particles and 3 processors the in the
	// domain from 0.0 to 4.0
	//
	// * processor 0 place particles from 0.0 to 1.332 (334 particles)
	// * processor 1 place particles from 1.336 to 2.668 (334 particles)
	// * processor 2 place particles from 2.672 to 4.0 (333 particles)
	//

	// It return how many particles the processors with id < rank has in total
	size_t base = vd.init_size_accum(Npart+1);
	auto it2 = vd.getIterator();

	while (it2.isNext())
	{
		auto key = it2.get();

		// set the position of the particles
		vd.template getPos<0>(key)[0] = (key.getKey() + base) * spacing;
		//set the property of the particles
		vd.template getProp<0>(key) = f_xex2((key.getKey() + base) * spacing);

		++it2;
	}

	//
	// ### WIKI 5 ###
	//
	// Once defined the position, we distribute them across the processors
	// following the decomposition, finally we get the ghost part
	//
	vd.map();
	vd.ghost_get<0>();

	//
	// ### WIKI 6 ###
	//
	// near the boundary we have two options, or we use one sided kernels,
	// or we add additional particles, it is required that such particles
	// produces a 2 time differentiable function. In order to obtain such
	// result we extend for x < 0.0 and x > 4.0 with the test function xe^(-x*x).
	//
	// Note that for x < 0.0 such extension is equivalent to mirror the
	// particles changing the sign of the strength
	//
	// \verbatim
	//
	// 0.6  -0.5      0.5  0.6   Strength
	//  +----+----*----*----*-
	//          0.0              Position
	//
	//  with * = real particle
	//       + = mirror particle
	//
	// \endverbatim
	//
	//
	Box<1,double> m_pad({0.0},{0.1});
	Box<1,double> m_pad2({3.9},{4.0});
	double enlarge = 0.1;

	// Create a box
	if (Npad * spacing > 0.1)
	{
		m_pad.setHigh(0,Npad * spacing);
		m_pad2.setLow(0,4.0 - Npad*spacing);
		enlarge = Npad * spacing;
	}

	auto it = vd.getDomainIterator();

	while (it.isNext())
	{
		auto key = it.get();

		// set the position of the particles
		if (m_pad.isInsideNB(vd.template getPos<0>(key)) == true)
		{
			vd.add();
			vd.template getLastPos<0>()[0] = - vd.template getPos<0>(key)[0];
			vd.template getLastProp<0>() = - vd.template getProp<0>(key);
		}

		// set the position of the particles
		if (m_pad2.isInsideNB(vd.template getPos<0>(key)) == true)
		{
			vd.add();
			vd.template getLastPos<0>()[0] = 2.0 * box.getHigh(0) - vd.template getPos<0>(key)[0];
			vd.template getLastProp<0>() = f_xex2(vd.template getLastPos<0>()[0]);
		}

		++it;
	}

	//
	// ### WIKI 6 ###
	//
	// We create a CellList with cell spacing 12 sigma
	//

    // get and construct the Cell list

	Ghost<1,double> gp(enlarge);
    auto cl = vd.getCellList(8*eps,gp);

    // Maximum infinity norm
    double linf = 0.0;

	//
	// ### WIKI 6 ###
	//
    // For each particle get the neighborhood of each particle
    //
    // This cycle is literally the formula from PSE operator approximation
	//
    //
    //
    //

    auto it_p = vd.getDomainIterator();
    while (it_p.isNext())
    {
    	// double PSE integration accumulator
    	double pse = 0;

    	// key
    	vect_dist_key_dx key = it_p.get();

    	// Get the position of the particles
    	Point<1,double> p = vd.template getPos<0>(key);

    	// We are not interested in calculating out the domain
    	// note added padding particle are considered domain particles
    	if (p.get(0) < 0.0 || p.get(0) >= 4.0)
    	{
    		++it_p;
    		continue;
    	}

    	// Get f(x) at the position of the particle
    	double prp_x = vd.template getProp<0>(key);

    	// Get the neighborhood of the particles
    	auto NN = cl.template getNNIterator<NO_CHECK>(cl.getCell(p));
    	while(NN.isNext())
    	{
    		auto nnp = NN.get();

    		// Calculate contribution given by the kernel value at position p,
    		// given by the Near particle
    		if (nnp != key.getKey())
    		{
    			// W(x-y)
    			double ker = lker.value(p,vd.template getPos<0>(nnp));

    			// f(y)
    			double prp_y = vd.template getProp<0>(nnp);

    			// 1.0/(eps)^2 [f(y)-f(x)]*W(x,y)*V_q
    			double prp = 1.0/eps/eps * (prp_y - prp_x) * spacing;
    			pse += prp * ker;
    		}

    		// Next particle
    		++NN;
    	}

    	// Here we calculate the L_infinity norm or the maximum difference
    	// of the analytic solution from the PSE calculated

    	double sol = Lapf_xex2(p);

    	if (fabs(pse - sol) > linf)
    		linf = fabs(pse - sol);

    	++it_p;
    }

	//
	// ### WIKI 7 ###
	//
    // Calculate the maximum infinity norm across processors and
    // print it
    //

    v_cl.max(linf);
    v_cl.execute();

    if (v_cl.getProcessUnitID() == 0)
    	std::cout << "Norm infinity: " << linf << "\n";

	//
	// ### WIKI 8 ###
	//
	// Deinitialize the library
	//
	delete_global_v_cluster();
}