OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
Dcpse_unit_tests.cpp
1//
2// Created by tommaso on 1/04/19.
3//
4
5#define BOOST_TEST_DYN_LINK
6
7#include <boost/test/unit_test.hpp>
8#include <boost/test/tools/floating_point_comparison.hpp>
9#include <Vector/vector_dist.hpp>
10#include <DCPSE/Dcpse.hpp>
11
12template<typename T>
13void check_small_or_close(T value, T expected, T tolerance)
14{
15 if (fabs(expected) < tolerance)
16 {
17 BOOST_CHECK_SMALL(value, tolerance);
18 } else
19 {
20 BOOST_CHECK_CLOSE(value, expected, tolerance);
21 }
22}
23
24template<typename T>
25void check_small_or_close_abs(T value, T expected, T absTolerance)
26{
27 BOOST_CHECK_SMALL(value - expected, absTolerance);
28}
29
30
31BOOST_AUTO_TEST_SUITE(Dcpse_tests)
32
33// If EIGEN is not present, EMatrix is not available and we don't need to build this test
34#ifdef HAVE_EIGEN
35
36 BOOST_AUTO_TEST_CASE(Dcpse_2D_test)
37 {
38 int rank;
39 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
40
41 BOOST_TEST_MESSAGE("Init vars...");
42
43 // Here build some easy domain and get some points around a given one
44 size_t edgeSemiSize = 20;
45 const size_t sz[2] = {2 * edgeSemiSize, 2 * edgeSemiSize};
46 Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
47 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
48 double spacing[2];
49 spacing[0] = 1.0 / (sz[0] - 1);
50 spacing[1] = 1.0 / (sz[1] - 1);
51 Ghost<2, double> ghost(0.1);
52
53 double rCut = 2 * spacing[0];
54
55 BOOST_TEST_MESSAGE("Init vector_dist...");
56 vector_dist<2, double, aggregate<double, double, double>> domain(0, box, bc, ghost);
57
58 if (rank == 0)
59 {
60 BOOST_TEST_MESSAGE("Init domain...");
61 auto it = domain.getGridIterator(sz);
62 size_t pointId = 0;
63 size_t counter = 0;
64 double minNormOne = 999;
65 while (it.isNext())
66 {
67 domain.add();
68 auto key = it.get();
69 mem_id k0 = key.get(0);
70 double x = k0 * spacing[0];
71 domain.getLastPos()[0] = x;
72 mem_id k1 = key.get(1);
73 double y = k1 * spacing[1];
74 domain.getLastPos()[1] = y;
75 // Here fill the function value
76 domain.template getLastProp<0>() = sin(x);
77// domain.template getLastProp<0>() = x * x;
78// domain.template getLastProp<0>() = x;
79 // Here fill the validation value for Df/Dx
80 domain.template getLastProp<2>() = cos(x);
81// domain.template getLastProp<2>() = 2 * x;
82// domain.template getLastProp<2>() = 1;
83
84 ++counter;
85 ++it;
86 }
87 BOOST_TEST_MESSAGE("Sync domain across processors...");
88 }
89 domain.map(); // Send particles to the right processors
90
91 BOOST_TEST_MESSAGE("Getting ghost...");
92 domain.ghost_get();
93
94 // Setup finished, actual test below...
95// std::cout << "rCut = " << rCut << std::endl;
96 BOOST_TEST_MESSAGE("DCPSE init & compute coefficients...");
97 Dcpse<2, vector_dist<2, double, aggregate<double, double, double>> > dcpse(domain, Point<2, unsigned int>({1, 0}), 2, rCut,2,support_options::N_PARTICLES);
98 BOOST_TEST_MESSAGE("DCPSE compute diff operator...");
99 dcpse.template computeDifferentialOperator<0, 1>(domain);
100
101 //domain.write("OUT_TEST");
102
103 // Now check against the validation values
104 BOOST_TEST_MESSAGE("Validating against ground truth...");
105// const double TOL = 1e-6;
106 const double avgSpacing = spacing[0] + spacing[1];
107 const double TOL = 2 * avgSpacing * avgSpacing;
108 auto itVal = domain.getDomainIterator();
109 bool check = true;
110 double computedValue;
111 double validationValue;
112 while (itVal.isNext())
113 {
114 auto key = itVal.get();
115 computedValue = domain.template getProp<1>(key);
116 validationValue = domain.template getProp<2>(key);
117 bool locCheck = (fabs(computedValue - validationValue) < TOL);
118 check = check && locCheck;
119 if (!locCheck)
120 {
121 std::cout << "FAILED CHECK :: pos=" << Point<2, double>(domain.getPos(key)).toString()
122 << ", computedValue=" << computedValue
123 << ", validationValue=" << validationValue
124 << ", difference=" << fabs(computedValue - validationValue)
125 << ", tolerance=" << TOL
126 << std::endl;
127// break;
128 }
129 ++itVal;
130 }
131 BOOST_REQUIRE(check);
132 }
133
134 BOOST_AUTO_TEST_CASE(Dcpse_2D_perturbed_test)
135 {
136 int rank;
137 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
138
139 BOOST_TEST_MESSAGE("Init vars...");
140
141 // Here build some easy domain and get some points around a given one
142 size_t edgeSemiSize = 20;
143 const size_t sz[2] = {2 * edgeSemiSize, 2 * edgeSemiSize};
144 Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
145 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
146 double spacing[2];
147 spacing[0] = 1.0 / (sz[0] - 1);
148 spacing[1] = 1.0 / (sz[1] - 1);
149
150 double sigma2 = spacing[0] * spacing[1] / (2 * 4);
151 //std::cout
152 // << "sigma2 = " << sigma2
153 // << ", spacing[0] = " << spacing[0]
154 // << std::endl;
155 double rCut = 2 * (spacing[0] + sqrt(sigma2));
156 Ghost<2, double> ghost(rCut);
157
158
159 BOOST_TEST_MESSAGE("Init vector_dist...");
160 vector_dist<2, double, aggregate<double, double, double>> domain(0, box, bc, ghost);
161
162 if (rank == 0)
163 {
164 BOOST_TEST_MESSAGE("Init domain...");
165// std::random_device rd{};
166// std::mt19937 rng{rd()};
167 std::mt19937 rng{6666666};
168
169 std::normal_distribution<> gaussian{0, sigma2};
170
171 auto it = domain.getGridIterator(sz);
172 size_t pointId = 0;
173 size_t counter = 0;
174 double minNormOne = 999;
175 while (it.isNext())
176 {
177 domain.add();
178 auto key = it.get();
179 mem_id k0 = key.get(0);
180 double x = k0 * spacing[0];
181 domain.getLastPos()[0] = x + gaussian(rng);
182 mem_id k1 = key.get(1);
183 double y = k1 * spacing[1];
184 domain.getLastPos()[1] = y + gaussian(rng);
185 // Here fill the function value
186 domain.template getLastProp<0>() = sin(domain.getLastPos()[0]);
187// domain.template getLastProp<0>() = x * x;
188// domain.template getLastProp<0>() = x;
189 // Here fill the validation value for Df/Dx
190 domain.template getLastProp<2>() = cos(domain.getLastPos()[0]);
191// domain.template getLastProp<2>() = 2 * x;
192// domain.template getLastProp<2>() = 1;
193
194 ++counter;
195 ++it;
196 }
197 BOOST_TEST_MESSAGE("Sync domain across processors...");
198
199 }
200 domain.map(); // Send particles to the right processors
201
202 // Setup finished, actual test below...
203// std::cout << "rCut = " << rCut << std::endl;
204 BOOST_TEST_MESSAGE("DCPSE init & compute coefficients...");
205 Dcpse<2, vector_dist<2, double, aggregate<double, double, double>> > dcpse(domain, Point<2, unsigned int>({1, 0}), 2, rCut, 2,support_options::N_PARTICLES);
206 BOOST_TEST_MESSAGE("DCPSE compute diff operator...");
207 dcpse.template computeDifferentialOperator<0, 1>(domain);
208
209 // Now check against the validation values
210 BOOST_TEST_MESSAGE("Validating against ground truth...");
211// const double TOL = 1e-6;
212 const double avgSpacing = spacing[0] + spacing[1] + 2 * sqrt(sigma2);
213 const double TOL = 2 * avgSpacing * avgSpacing;
214 auto itVal = domain.getDomainIterator();
215 bool check = true;
216 double computedValue;
217 double validationValue;
218 while (itVal.isNext())
219 {
220 auto key = itVal.get();
221 computedValue = domain.template getProp<1>(key);
222 validationValue = domain.template getProp<2>(key);
223 bool locCheck = (fabs(computedValue - validationValue) < TOL);
224 check = check && locCheck;
225 if (!locCheck)
226 {
227 std::cout << "FAILED CHECK :: pos=" << Point<2, double>(domain.getPos(key)).toString()
228 << ", computedValue=" << computedValue
229 << ", validationValue=" << validationValue
230 << ", difference=" << fabs(computedValue - validationValue)
231 << ", tolerance=" << TOL
232 << std::endl;
233// break;
234 }
235 ++itVal;
236 }
237 BOOST_REQUIRE(check);
238 }
239
240 BOOST_AUTO_TEST_CASE(Dcpse_3D_test)
241 {
242 int rank;
243 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
244
245 BOOST_TEST_MESSAGE("Init vars...");
246
247 // Here build some easy domain and get some points around a given one
248 size_t edgeSemiSize = 20;
249 const unsigned int DIM = 3;
250 const size_t sz[DIM] = {2 * edgeSemiSize, 2 * edgeSemiSize, 2 * edgeSemiSize};
251 Box<DIM, double> box({0, 0, 0}, {2 * M_PI, 2 * M_PI, 2 * M_PI});
252 size_t bc[DIM] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
253 double spacing[DIM];
254 spacing[0] = 1.0 / (sz[0] - 1);
255 spacing[1] = 1.0 / (sz[1] - 1);
256 spacing[2] = 1.0 / (sz[2] - 1);
257 Ghost<DIM, double> ghost(0.1);
258
259 double rCut = 2 * spacing[0];
260
261 BOOST_TEST_MESSAGE("Init vector_dist...");
262 vector_dist<DIM, double, aggregate<double, double, double>> domain(0, box, bc, ghost);
263
264 if (rank == 0)
265 {
266 BOOST_TEST_MESSAGE("Init domain...");
267 auto it = domain.getGridIterator(sz);
268 size_t pointId = 0;
269 size_t counter = 0;
270 double minNormOne = 999;
271 while (it.isNext())
272 {
273 domain.add();
274 auto key = it.get();
275 mem_id k0 = key.get(0);
276 double x = k0 * spacing[0];
277 domain.getLastPos()[0] = x;
278 mem_id k1 = key.get(1);
279 double y = k1 * spacing[1];
280 domain.getLastPos()[1] = y;
281 mem_id k2 = key.get(2);
282 double z = k2 * spacing[2];
283 domain.getLastPos()[2] = z;
284 // Here fill the function value
285 domain.template getLastProp<0>() = sin(z);
286// domain.template getLastProp<0>() = x * x;
287// domain.template getLastProp<0>() = x;
288 // Here fill the validation value for Df/Dx
289 domain.template getLastProp<2>() = cos(z);
290// domain.template getLastProp<2>() = 2 * x;
291// domain.template getLastProp<2>() = 1;
292
293 ++counter;
294 ++it;
295 }
296 BOOST_TEST_MESSAGE("Sync domain across processors...");
297 }
298 domain.map(); // Send particles to the right processors
299
300 // Setup finished, actual test below...
301 BOOST_TEST_MESSAGE("DCPSE init & compute coefficients...");
302 Dcpse<DIM, vector_dist<DIM, double, aggregate<double, double, double>>> dcpse(domain, Point<DIM, unsigned int>({0, 0, 1}), 2, rCut,2.5,support_options::N_PARTICLES);
303 BOOST_TEST_MESSAGE("DCPSE compute diff operator...");
304 dcpse.template computeDifferentialOperator<0, 1>(domain);
305
306 // Now check against the validation values
307 BOOST_TEST_MESSAGE("Validating against ground truth...");
308// const double TOL = 1e-6;
309 const double avgSpacing = spacing[0] + spacing[1] + spacing[2];
310 const double TOL = avgSpacing * avgSpacing;
311 auto itVal = domain.getDomainIterator();
312 bool check = true;
313 double computedValue;
314 double validationValue;
315 while (itVal.isNext())
316 {
317 auto key = itVal.get();
318 computedValue = domain.template getProp<1>(key);
319 validationValue = domain.template getProp<2>(key);
320 check = check && (fabs(computedValue - validationValue) < TOL);
321 if (!check)
322 {
323 std::cout << "FAILED CHECK :: pos=" << Point<DIM, double>(domain.getPos(key)).toString()
324 << ", computedValue=" << computedValue
325 << ", validationValue=" << validationValue
326 << ", tolerance=" << TOL
327 << std::endl;
328 break;
329 }
330 ++itVal;
331 }
332 BOOST_REQUIRE(check);
333 }
334
335 BOOST_AUTO_TEST_CASE(Dcpse_3D_2_test)
336 {
337 int rank;
338 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
339
340 BOOST_TEST_MESSAGE("Init vars...");
341
342 // Here build some easy domain and get some points around a given one
343 size_t edgeSemiSize = 10;
344 const unsigned int DIM = 3;
345 const size_t sz[DIM] = {2 * edgeSemiSize, 2 * edgeSemiSize, 2 * edgeSemiSize};
346 Box<DIM, double> box({0, 0, 0}, {2 * M_PI, 2 * M_PI, 2 * M_PI});
347 size_t bc[DIM] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
348 double spacing[DIM];
349 spacing[0] = 1.0 / (sz[0] - 1);
350 spacing[1] = 1.0 / (sz[1] - 1);
351 spacing[2] = 1.0 / (sz[2] - 1);
352
353 double rCut = 2 * spacing[0];
354 Ghost<DIM, double> ghost(rCut);
355
356
357 BOOST_TEST_MESSAGE("Init vector_dist...");
358 vector_dist<DIM, double, aggregate<double, double, double>> domain(0, box, bc, ghost);
359
360 if (rank == 0)
361 {
362 BOOST_TEST_MESSAGE("Init domain...");
363 auto it = domain.getGridIterator(sz);
364 size_t pointId = 0;
365 size_t counter = 0;
366 double minNormOne = 999;
367 while (it.isNext())
368 {
369 domain.add();
370 auto key = it.get();
371 mem_id k0 = key.get(0);
372 double x = k0 * spacing[0];
373 domain.getLastPos()[0] = x;
374 mem_id k1 = key.get(1);;
375 double y = k1 * spacing[1];
376 domain.getLastPos()[1] = y;
377 mem_id k2 = key.get(2);
378 double z = k2 * spacing[2];
379 domain.getLastPos()[2] = z;
380 // Here fill the function value
381 domain.template getLastProp<0>() = x * x * sin(z);
382 // Here fill the validation value for Df/Dx
383 domain.template getLastProp<2>() = -2*sin(z);
384
385 ++counter;
386 ++it;
387 }
388 BOOST_TEST_MESSAGE("Sync domain across processors...");
389 }
390 domain.map(); // Send particles to the right processors
391
392 // Setup finished, actual test below...
393 unsigned int r = 2;
394 BOOST_TEST_MESSAGE("DCPSE init & compute coefficients...");
395 Dcpse<DIM, vector_dist<DIM, double, aggregate<double, double, double>> > dcpse(domain, Point<DIM, unsigned int>({2, 0, 2}), r, rCut,2.5,support_options::N_PARTICLES);
396 BOOST_TEST_MESSAGE("DCPSE compute diff operator...");
397 dcpse.template computeDifferentialOperator<0, 1>(domain);
398
399 // Now check against the validation values
400 BOOST_TEST_MESSAGE("Validating against ground truth...");
401// const double TOL = 1e-6;
402 const double avgSpacing = spacing[0] + spacing[1] + spacing[2];
403 const double TOL = pow(avgSpacing, r);
404 auto itVal = domain.getDomainIterator();
405 bool check = true;
406 double computedValue;
407 double validationValue;
408 while (itVal.isNext())
409 {
410 auto key = itVal.get();
411 computedValue = domain.template getProp<1>(key);
412 validationValue = domain.template getProp<2>(key);
413 check = check && (fabs(computedValue - validationValue) < TOL);
414 if (!check)
415 {
416 std::cout << "FAILED CHECK :: pos=" << Point<DIM, double>(domain.getPos(key)).toString()
417 << ", computedValue=" << computedValue
418 << ", validationValue=" << validationValue
419 << ", tolerance=" << TOL
420 << std::endl;
421 break;
422 }
423 ++itVal;
424 }
425 BOOST_REQUIRE(check);
426 }
427
428#endif // HAVE_EIGEN
429
430BOOST_AUTO_TEST_SUITE_END()
Box
This class represent an N-dimensional box.
Definition Box.hpp:61
Point
This class implement the point shape in an N-dimensional space.
Definition Point.hpp:28
Point::toString
std::string toString() const
Return the string with the point coordinate.
Definition Point.hpp:398
vector_dist
Distributed vector.
Definition vector_dist.hpp:305