OpenFPM_pdata  4.1.0 Project that contain the implementation of distributed structures
main.cpp
1 /*
2  * ### WIKI 1 ###
3  *
4  * ## Simple example
5  *
6  * In this example we show The 1D PSE Diffusion equation
7  *
8  * in particular we integrate the following equation
9  *
10  * \$\frac{u}{t} = D \laplacian u \$
11  *
12  * ### WIKI END ###
13  *
14  */
15
16 #include "Vector/vector_dist.hpp"
17 #include "Decomposition/CartDecomposition.hpp"
18 #include "PSE/Kernels.hpp"
19 #include "Plot/util.hpp"
20 #include "Plot/GoogleChart.hpp"
21 #include "data_type/aggregate.hpp"
22 #include <cmath>
23
24
25 /*
26  * ### WIKI 2 ###
27  *
28  * Here we define the function x*e^(-x*x) the initial condition
29  *
30  */
31
32 double f_xex2(double x)
33 {
34  return x*exp(-x*x);
35 }
36
37 double f_xex2(Point<1,double> & x)
38 {
39  return x.get(0)*exp(-x.get(0)*x.get(0));
40 }
41
42
43 /*
44  *
45  * ### WIKI 3 ###
46  *
47  * It calculate the Laplacian of the function in one point using the PSE
48  * Approximation
49  *
50  * \$\frac{1}{\epsilon^{2}} h (u_{q} - u_{p}) \eta_{\epsilon}(x_q - x_p) \$
51  *
52  * \param p point
53  * \param key vector
54  * \param vd Distributed vector
55  * \param eps epsilon of the kernel
56  * \param spacing between particles
57  * \param cl CellList
58  *
59  */
60 template<typename CellL> double calcLap(Point<1,double> p, vect_dist_key_dx key, vector_dist<1,double, aggregate<double,double> > & vd, double eps, double spacing, CellL & cl)
61 {
62  // Laplacian PSE kernel<1,double,2> = 1 dimension, on double, second order
63  Lap_PSE<1,double,2> lker(eps);
64
65  // double PSE integration accumulator
66  double pse = 0;
67
68  // f(x_p)
69  double prp_x = vd.template getProp<0>(key);
70
71  // Get the neighborhood of the particles
72  auto NN = cl.template getNNIterator<NO_CHECK>(cl.getCell(p));
73  while(NN.isNext())
74  {
75  auto nnp = NN.get();
76
77  // Calculate contribution given by the particle q near to p
78  if (nnp != key.getKey())
79  {
80  // W(x_p-x_q) x = position of p, y = position of q
81  double ker = lker.value(p,vd.getPos(nnp));
82
83  // f(x_q)
84  double prp_y = vd.template getProp<0>(nnp);
85
86  // 1.0/(eps)^2 [f(x_q)-f(x_p)]*W(x_p-x_q)*V_q
87  double prp = 1.0/eps/eps * (prp_y - prp_x) * spacing;
88  pse += prp * ker;
89  }
90
91  // Next particle
92  ++NN;
93  }
94
95  // return the calculated laplacian
96  return pse;
97 }
98
99 /*
100  *
101  * ### WIKI 4 ###
102  *
103  * It mirror the particles at the border of the domain, in particular on the left we do the operation
104  * show in figure
105  *
106  * \verbatim
107  *
108  * -0.6 -0.5 0.5 0.6 Strength
109  +----+----*----*----*-
110  0.0 Position
111
112  with * = real particle
113  + = mirror particle
114
115  \endverbatim
116  *
117  * We create particle on the negative part of the domain and we assign a strength opposite to the particle in the
118  * positive part
119  *
120  * On the right instead we prolong the initial condition (Ideally should remain almost fixed)
121  *
122  * \param vd vector of particles
123  * \param key particle
124  * \param m_pad border box(interval) on the left
125  * \param m_pad2 border box(interval) on the right
126  * \param box domain
127  *
128  */
129 inline void mirror(vector_dist<1,double, aggregate<double,double> > & vd, vect_dist_key_dx & key, const Box<1,double> & m_pad, Box<1,double> & m_pad2, const Box<1,double> & box)
130 {
131  // If the particle is inside the box (interval), it mean that is at the left border
132  if (m_pad.isInsideNB(vd.getPos(key)) == true)
133  {
134  // Add a new particle
135  vd.add();
136
137  //Set the added particle position, symmetrically positioned on the negative part of the domain
138  vd.getLastPos()[0] = - vd.getPos(key)[0];
139
140  // Set the property of the particle to be negative
141  vd.template getLastProp<0>() = - vd.template getProp<0>(key);
142  }
143
144  // If the particle is inside the box (interval), it mean that is at the right border
145  if (m_pad2.isInsideNB(vd.getPos(key)) == true)
146  {
147  // Add a new particle
148  vd.add();
149
150  //Set the added particle position, symmetrically positioned around x = 4.0
151  // 4.0 Added
152  // *------|------*
153  //
154  vd.getLastPos()[0] = 2.0 * box.getHigh(0) - vd.getPos(key)[0];
155
156  // no flow boundary condition
157  vd.template getLastProp<0>() = vd.template getProp<0>(key);
158  }
159 }
160
161 /*
162  *
163  * ### WIKI END ###
164  *
165  */
166
167 int main(int argc, char* argv[])
168 {
169  //
170  // ### WIKI 5 ###
171  //
172  // Some useful parameters.
173  //
174
175  // Number of particles
176  const size_t Npart = 1000;
177
178  // A different way to write Npart
179  size_t NpartA[1] = {Npart+1};
180
181  // Number of steps
182  const size_t Nstep = 1000;
183
184  // Delta t
185  const double dt = 10.0 / Nstep;
186
187  // Number of padding particles, particles outside the domain
188  const size_t Npad = 20;
189
190  // The domain
191  Box<1,double> box({0.0},{4.0});
192
193  // Calculated spacing
194  double spacing = box.getHigh(0) / Npart;
195
196  // The mollification length
197  const double eps = 2*spacing;
198
199  // Diffusion constant
200  double D = 1e-4;
201
202  // 2 constants
203  constexpr int U = 0;
204  constexpr int Lap_U = 1;
205
206  //
207  // ### WIKI 6 ###
208  //
209  // Here we Initialize the library and we define Ghost size
210  // and non-periodic boundary conditions
211  //
212  openfpm_init(&argc,&argv);
213  Vcluster<> & v_cl = create_vcluster();
214
215  size_t bc[1]={NON_PERIODIC};
216  Ghost<1,double> g(12*eps);
217
218  //
219  // ### WIKI 7 ###
220  //
221  // Here we are creating a distributed vector defined by the following parameters
222  //
223  // <1,double, aggregate<double,double> > = 1 dimension, space with type double, 2 properties double,double (aggregate is like tuples, a way to encapsulate types). The first property is the field value u, the second is the Laplacian of u
224  // with Npart+1 particles, on the box domain [0.0,4.0], with bc boundary condition NON_PERIODIC, and an extended ghost part of size g
225  //
226  vector_dist<1,double, aggregate<double,double> > vd(Npart+1,box,bc,g);
227
228  //
229  // ### WIKI 8 ###
230  //
231  // We assign the position to the particles in order to be equally spaced in a grid like position, the scalar property is set to
232  // the function x*e^(-x*x) value the initial condition
233  //
234  // P.S. In a Multi processor context getGridIterator adapt based on the underline space division across processors
235  //
236
237  // Create a grid like iterator
238  auto itp = vd.getGridIterator(NpartA);
239
240  // Counter starting from 0
241  size_t p = 0;
242
243  // For each point in the grid
244  while (itp.isNext())
245  {
246  // Get the grid index
247  auto key = itp.get();
248
249  // Assign the particle position
250  vd.getPos(p)[0] = key.get(0) * spacing;
251
252  // Assign the property 0 based on the initial condition
253  vd.template getProp<U>(p) = f_xex2(key.get(0) * spacing);
254
255  // Next particle
256  ++p;
257
258  // Next grid point
259  ++itp;
260  }
261
262  //
263  // ### WIKI 9 ###
264  //
265  // Once defined the position, we distribute them across the processors
266  // following the decomposition, finally we update the ghost part for the vield U
267  //
268  vd.map();
269  vd.ghost_get<U>();
270
271  //
272  // ### WIKI 10 ###
273  //
274  // near the boundary we have two options, or we use one sided kernels,
275  // or we add additional particles, it is required that such particles
276  // produces a 2 time differentiable function. In order to obtain such
277  // result we extend for x < 0.0 and x > 4.0 with the test function xe^(-x*x).
278  //
279  // Note that for x < 0.0 such extension is equivalent to mirror the
280  // particles changing the sign of the strength
281  //
282  // \verbatim
283  //
284  //-0.6 -0.5 0.5 0.6 Strength
285  // +----+----*----*----*-
286  // 0.0 Position
287  //
288  // with * = real particle
289  // + = mirror particle
290  //
291  // \endverbatim
292  //
293  //
294
295  // The particle that we have to mirror are ...
296
297  // ... on the box m_pad for the laft part
298  Box<1,double> m_pad({0.0},{0.1});
299
300  // ... on the box m_pad2 for the right part
301  Box<1,double> m_pad2({3.9},{4.0});
302
303  // It contain how much we enlarge the domain
304  double enlarge = m_pad.getHigh(0) - m_pad.getLow(0);
305
306  // In general 0.1 should be a sufficent padding, but in case it is not we recalculate the left and right boxes
307  if (Npad * spacing > 0.1)
308  {
309  m_pad.setHigh(0,Npad * spacing);
310  m_pad2.setLow(0,4.0 - Npad*spacing);
311  enlarge = Npad * spacing;
312  }
313
314
315  //
316  // ### WIKI 11 ###
317  //
318  // Here we mirror the particles if needed
319  //
320  //
321
322  auto it = vd.getDomainIterator();
323  size_t n_part = vd.size_local();
324
325  while (it.isNext())
326  {
327  auto key = it.get();
328
329  mirror(vd,key,m_pad,m_pad2,box);
330
331  ++it;
332  }
333
334  //
335  // ### WIKI 12 ###
336  //
337  // We create a CellList with cell spacing 8 epsilon, on the part of space domain + ghost + padding
338  // The padding area, can be bigger than the ghost part. In practice there is no reason to do it, but
339  // keeping ghost size and padding area unrelated give us the possibility to show how to create a CellList
340  // on a area bigger than the domain + ghost
341
342  Ghost<1,double> gp(enlarge);
343
344  // Create a Cell list with Cell spaping 8*epsilon, the CellList is created on a domain+ghost+enlarge space
345  auto cl = vd.getCellList(8*eps,gp);
346
347  // Maximum infinity norm
348  double linf = 0.0;
349
350  //
351  // ### WIKI 13 ###
352  //
353  // Euler time integration until t = 10.0
354  //
355  //
356
357  for (double t = 0; t <= 10.0 ; t += dt)
358  {
359
360  // Iterate all the particles, including padding because are real particles
361  auto it_p = vd.getDomainIterator();
362  while (it_p.isNext())
363  {
364  // Get particle p
365  auto p = it_p.get();
366
367  // Get the position of the particle p
368  Point<1,double> x_p = vd.getPos(p);
369
370  // We are not interested in calculating out the domain (So on the padded particles)
371  if (x_p.get(0) < 0.0 || x_p.get(0) >= 4.0)
372  {
373  ++it_p;
374  continue;
375  }
376
377  // Here we calculate the Laplacian of u on position x_p and we store the result on the particle p
378  vd.template getProp<Lap_U>(p) = calcLap(x_p,p,vd,eps,spacing,cl);
379
380  // Next particle
381  ++it_p;
382  }
383
384  // Once we calculated the Laplacian we do not need enymore the padding particle eliminate them
385  // n_part is the original size of the vector without padding particles
386  vd.resize(n_part);
387
388  // Iterate again on all particles
389  auto it_p2 = vd.getDomainIterator();
390  while (it_p2.isNext())
391  {
392  // Get particle p
393  auto p = it_p2.get();
394
395  // Get the Laplacian
396  double pse = vd.template getProp<Lap_U>(p);
397
398  // Integrate with Euler step
399  vd.template getProp<0>(p) += D * pse * dt;
400
401  // mirror again the particles if needed (this for the next step)
402  mirror(vd,p,m_pad,m_pad2,box);
403
404  // Next particle
405  ++it_p2;
406  }
407
408  // Update the ghost part
409  vd.ghost_get<U>();
410  }
411
412  //
413  // ### WIKI 14 ###
414  //
415  // U now contain the solution scattered across processors. Because it is not possible plot on multiple
416  // processors (Google chart does not support it) we collect the solution on one processor
417  //
418
419  struct pos_val
420  {
421  // position of the particle
422  double pos;
423
424  // value of U of the particle
425  double val;
426
427  // usefull to sort the particle by position
428  bool operator<(const pos_val & p) const
429  {
430  return pos < p.pos;
431  }
432  };
433
434  // Resize a vector with the local size of the vector
435  openfpm::vector<pos_val> v(vd.size_local());
436
437  // This will contail all the particles on the master processor
439
440
441  // Copy the particle position and the field U, into the vector v
442  for (size_t i = 0 ; i < vd.size_local() ; i++)
443  {
444  // particle position
445  v.get(i).pos = vd.getPos(i)[0];
446
447  // Field U
448  v.get(i).val = vd.template getProp<U>(i);
449  }
450
451  // from the partial solution vector v we Gather the total solution v_tot on master (processor=0)
452  v_cl.SGather(v,v_tot,0);
453
454  //
455  // ### WIKI 15 ###
456  //
457  // Once we gather the solution on master we plot it. Only processor 0 do this
458  //
459
460  // Only processor = 0
461  if (v_cl.getProcessUnitID() == 0)
462  {
463  // sort the solution by particle position
464  v_tot.sort();
465
466  // y contain 2 vectors the first is the calculated solution
467  // the second contain the initial condition
469
470  // vector with the particle position for plotting
472
473  // We fill the plotting variables
474  for (size_t i = 0 ; i < v_tot.size() ; i++)
475  {
476  // fill x with the particle position
477  x.add(v_tot.get(i).pos);
478
479  // fill y with the value of the field U contained by the particles
480  y.get(0).add(v_tot.get(i).val);
481
482  // fill y with the initial condition
483  y.get(1).add(f_xex2(v_tot.get(i).pos));
484  }
485
486  // Plot with GooglePlot
487  // Google charts options
488  GCoptions options;
489
490  // Assign a name for the Y axis
491  options.yAxis = std::string("U Axis");
492
493  // Assign a name for the X Axis
494  options.xAxis = std::string("X Axis");
495
496  // width of the line
497  options.lineWidth = 1.0;
498
499  // Create a google plot object
500  GoogleChart cg;
501
502  // add a line plot
503  cg.AddLinesGraphT(x,y,options);
504
505  // write the plot file
506  cg.write("PSE_plot.html");
507  }
508
509  //
510  // ### WIKI 16 ###
511  //
512  // Deinitialize the library
513  //
514  openfpm_finalize();
515 }
Derivative second order on h (spacing)
Definition: Derivative.hpp:28
size_t getProcessUnitID()
Get the process unit id.
__device__ __host__ T getLow(int i) const
get the i-coordinate of the low bound interval of the box
Definition: Box.hpp:556
This class implement the point shape in an N-dimensional space.
Definition: Point.hpp:27
Grid key for a distributed grid.
size_t size()
Stub size.
Definition: map_vector.hpp:211
void AddLinesGraphT(openfpm::vector< X > &x, openfpm::vector< Y > &y, const GCoptions &opt)
Add lines graph.
Definition: Ghost.hpp:39
Implementation of the Laplacian kernels for PSE.
Definition: Kernels.hpp:25
Implementation of VCluster class.
Definition: VCluster.hpp:58
__device__ __host__ void setHigh(int i, T val)
set the high interval of the box
Definition: Box.hpp:544
std::string yAxis
Y axis name.
Definition: GoogleChart.hpp:30
Small class to produce graph with Google chart in HTML.
bool SGather(T &send, S &recv, size_t root)
Semantic Gather, gather the data from all processors into one node.
Definition: VCluster.hpp:450
__device__ __host__ const T & get(unsigned int i) const
Get coordinate.
Definition: Point.hpp:172
__device__ __host__ void setLow(int i, T val)
set the low interval of the box
Definition: Box.hpp:533
size_t lineWidth
Width of the line.
Definition: GoogleChart.hpp:56
This class represent an N-dimensional box.
Definition: Box.hpp:60
__device__ __host__ size_t getKey() const
Get the key.
void write(std::string file)
It write the graphs on file in html format using Google charts.
Distributed vector.
Definition: eq.hpp:67
bool isInsideNB(const Point< dim, T > &p) const
Check if the point is inside the region excluding the borders.
Definition: Box.hpp:1095
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
Definition: aggregate.hpp:214
Implementation of 1-D std::vector like structure.
Definition: map_vector.hpp:202
__device__ __host__ T getHigh(int i) const
get the high interval of the box
Definition: Box.hpp:567
std::string xAxis
X axis name.
Definition: GoogleChart.hpp:32
Google chart options.
Definition: GoogleChart.hpp:25