1/*
2 * Kernels_unit_tests.hpp
3 *
4 * Created on: Feb 16, 2016
5 * Author: i-bird
6 */
7
8#ifndef OPENFPM_NUMERICS_SRC_PSE_KERNELS_TEST_UTIL_HPP_
9#define OPENFPM_NUMERICS_SRC_PSE_KERNELS_TEST_UTIL_HPP_
10
11#include "Kernels.hpp"
12#include "Space/Ghost.hpp"
13#include "Vector/vector_dist.hpp"
14#include "data_type/aggregate.hpp"
15#include "Decomposition/CartDecomposition.hpp"
16
17struct PSEError
18{
19 double l2_error;
20 double linf_error;
21};
22
23
24template<typename T> T f_xex2(T x)
25{
26 return x*exp(-x*x);
27}
28
29template<typename T> T f_xex2(Point<1,T> & x)
30{
31 return x.get(0)*exp(-x.get(0)*x.get(0));
32}
33
34template<typename T> T Lapf_xex2(Point<1,T> & x)
35{
36 return 2.0*x.get(0)*(2.0*x.get(0)*x.get(0) - 3.0)*exp(-x.get(0)*x.get(0));
37}
38
39/*
40 * Given the Number of particles, it calculate the error produced by the standard
41 * second order PSE kernel to approximate the laplacian of x*e^{-x^2} in the point
42 * 0.5 (1D)
43 *
44 * The spacing h is calculated as
45 *
46 * \$ h = 1.0 / N_{part} \$
47 *
48 * epsilon of the 1D kernel is calculated as
49 *
50 * \$ \epsilon=overlap * h \$
51 *
52 * this mean that
53 *
54 * \$ overlap = \frac{\epsilon}{h} \$
55 *
56 * While the particles that are considered in the neighborhood of 0.5 are
57 * calculated as
58 *
59 * \$ \N_contrib = 20*eps/spacing \$
60 *
61 * \param N_part Number of particles in the domain
62 * \param overlap overlap parameter
63 *
64 *
65 */
66template<typename T, typename Kernel> void PSE_test(size_t Npart, size_t overlap,struct PSEError & error)
67{
68 // The domain
69 Box<1,T> box({0.0},{4.0});
70
71 // Calculated spacing
72 T spacing = box.getHigh(0) / Npart;
73
74 // Epsilon of the particle kernel
75 const T eps = overlap*spacing;
76
77 // Laplacian PSE kernel 1 dimension, on double, second order
78 Kernel lker(eps);
79
80 // For convergence we need less particles
81 Npart = static_cast<size_t>(20*eps/spacing);
82
83 // Middle particle
84 long int mp = Npart / 2;
85
86 size_t bc[1]={NON_PERIODIC};
87 Ghost<1,T> g(20*eps);
88
89 vector_dist<1,T, aggregate<T> > vd(Npart,box,bc,g);
90
91 auto it2 = vd.getIterator();
92
93 while (it2.isNext())
94 {
95 auto key = it2.get();
96
97 // set the position of the particles
98 vd.getPos(key)[0] = 0.448000 - ((long int)key.getKey() - mp) * spacing;
99 //set the property of the particles
100 T d = vd.getPos(key)[0];
101
102 vd.template getProp<0>(key) = f_xex2(d);
103
104 ++it2;
105 }
106
107 vect_dist_key_dx key;
108 key.setKey(mp);
109
110 // get and construct the Cell list
111 auto cl = vd.getCellList(0.5);
112
113 // Maximum infinity norm
114 double linf = 0.0;
115
116 // PSE integration accumulator
117 T pse = 0;
118
119 // Get the position of the particle
120 Point<1,T> p = vd.getPos(key);
121
122 // Get f(x) at the position of the particle
123 T prp_x = vd.template getProp<0>(key);
124
125 // Get the neighborhood of the particles
126 auto NN = cl.template getNNIterator<NO_CHECK>(cl.getCell(p));
127 while(NN.isNext())
128 {
129 auto nnp = NN.get();
130
131 // Calculate contribution given by the kernel value at position p,
132 // given by the Near particle, exclude itself
133 if (nnp != key.getKey())
134 {
135 Point<1,T> xp = vd.getPos(nnp);
136
137 // W(x-y)
138 T ker = lker.value(p,xp);
139
140 // f(y)
141 T prp_y = vd.template getProp<0>(nnp);
142
143 // 1.0/(eps)^2 [f(y)-f(x)]*W(x,y)*V_q
144 T prp = 1.0/eps/eps * spacing * (prp_y - prp_x);
145 pse += prp * ker;
146 }
147
148 // Next particle
149 ++NN;
150 }
151
152 // Here we calculate the L_infinity norm or the maximum difference
153 // of the analytic solution from the PSE calculated
154
155 T sol = Lapf_xex2(p);
156
157 if (fabs(pse - sol) > linf)
158 linf = static_cast<double>(fabs(pse - sol));
159
160
161 error.linf_error = (double)linf;
162}
163
164
165#endif /* OPENFPM_NUMERICS_SRC_PSE_KERNELS_TEST_UTIL_HPP_ */
166