OpenFPM  5.2.0
Project that contain the implementation of distributed structures
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 
12 template<typename T>
13 void 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 
24 template<typename T>
25 void check_small_or_close_abs(T value, T expected, T absTolerance)
26 {
27  BOOST_CHECK_SMALL(value - expected, absTolerance);
28 }
29 
30 
31 BOOST_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 = 3.1 * 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  auto verletList = domain.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
98  Dcpse<2, decltype(verletList), vector_dist<2, double, aggregate<double, double, double>> > dcpse(domain, verletList, Point<2, unsigned int>({1, 0}), 2, rCut, support_options::RADIUS);
99  BOOST_TEST_MESSAGE("DCPSE compute diff operator...");
100  dcpse.template computeDifferentialOperator<0, 1>(domain);
101 
102  //domain.write("OUT_TEST");
103 
104  // Now check against the validation values
105  BOOST_TEST_MESSAGE("Validating against ground truth...");
106 // const double TOL = 1e-6;
107  const double avgSpacing = spacing[0] + spacing[1];
108  const double TOL = 2 * avgSpacing * avgSpacing;
109  auto itVal = domain.getDomainIterator();
110  bool check = true;
111  double computedValue;
112  double validationValue;
113  while (itVal.isNext())
114  {
115  auto key = itVal.get();
116  computedValue = domain.template getProp<1>(key);
117  validationValue = domain.template getProp<2>(key);
118  bool locCheck = (fabs(computedValue - validationValue) < TOL);
119  check = check && locCheck;
120  if (!locCheck)
121  {
122  std::cout << "FAILED CHECK :: pos=" << Point<2, double>(domain.getPos(key)).toString()
123  << ", computedValue=" << computedValue
124  << ", validationValue=" << validationValue
125  << ", difference=" << fabs(computedValue - validationValue)
126  << ", tolerance=" << TOL
127  << std::endl;
128 // break;
129  }
130  ++itVal;
131  }
132  BOOST_REQUIRE(check);
133  }
134 
135  BOOST_AUTO_TEST_CASE(Dcpse_2D_perturbed_test)
136  {
137  int rank;
138  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
139 
140  BOOST_TEST_MESSAGE("Init vars...");
141 
142  // Here build some easy domain and get some points around a given one
143  size_t edgeSemiSize = 20;
144  const size_t sz[2] = {2 * edgeSemiSize, 2 * edgeSemiSize};
145  Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
146  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
147  double spacing[2];
148  spacing[0] = 1.0 / (sz[0] - 1);
149  spacing[1] = 1.0 / (sz[1] - 1);
150 
151  double sigma2 = spacing[0] * spacing[1] / (4);
152  //std::cout
153  // << "sigma2 = " << sigma2
154  // << ", spacing[0] = " << spacing[0]
155  // << std::endl;
156  double rCut = 2 * (spacing[0] + sqrt(sigma2));
157  Ghost<2, double> ghost(rCut);
158 
159 
160  BOOST_TEST_MESSAGE("Init vector_dist...");
161  vector_dist<2, double, aggregate<double, double, double>> domain(0, box, bc, ghost);
162 
163  if (rank == 0)
164  {
165  BOOST_TEST_MESSAGE("Init domain...");
166 // std::random_device rd{};
167 // std::mt19937 rng{rd()};
168  std::mt19937 rng{6666666};
169 
170  std::normal_distribution<> gaussian{0, sigma2};
171 
172  auto it = domain.getGridIterator(sz);
173  size_t pointId = 0;
174  size_t counter = 0;
175  double minNormOne = 999;
176  while (it.isNext())
177  {
178  domain.add();
179  auto key = it.get();
180  mem_id k0 = key.get(0);
181  double x = k0 * spacing[0];
182  domain.getLastPos()[0] = x + gaussian(rng);
183  mem_id k1 = key.get(1);
184  double y = k1 * spacing[1];
185  domain.getLastPos()[1] = y + gaussian(rng);
186  // Here fill the function value
187  domain.template getLastProp<0>() = sin(domain.getLastPos()[0]);
188 // domain.template getLastProp<0>() = x * x;
189 // domain.template getLastProp<0>() = x;
190  // Here fill the validation value for Df/Dx
191  domain.template getLastProp<2>() = cos(domain.getLastPos()[0]);
192 // domain.template getLastProp<2>() = 2 * x;
193 // domain.template getLastProp<2>() = 1;
194 
195  ++counter;
196  ++it;
197  }
198  BOOST_TEST_MESSAGE("Sync domain across processors...");
199 
200  }
201  domain.map(); // Send particles to the right processors
202  domain.ghost_get();
203 
204  // Setup finished, actual test below...
205 // std::cout << "rCut = " << rCut << std::endl;
206  BOOST_TEST_MESSAGE("DCPSE init & compute coefficients...");
207  auto verletList = domain.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
208  Dcpse<2, decltype(verletList), vector_dist<2, double, aggregate<double, double, double>> > dcpse(domain, verletList, Point<2, unsigned int>({1, 0}), 2, rCut, support_options::RADIUS);
209  BOOST_TEST_MESSAGE("DCPSE compute diff operator...");
210  dcpse.template computeDifferentialOperator<0, 1>(domain);
211 
212  // Now check against the validation values
213  BOOST_TEST_MESSAGE("Validating against ground truth...");
214 // const double TOL = 1e-6;
215  const double avgSpacing = spacing[0] + spacing[1] + 2 * sqrt(sigma2);
216  const double TOL = 2 * avgSpacing * avgSpacing;
217  auto itVal = domain.getDomainIterator();
218  bool check = true;
219  double computedValue;
220  double validationValue;
221  while (itVal.isNext())
222  {
223  auto key = itVal.get();
224  computedValue = domain.template getProp<1>(key);
225  validationValue = domain.template getProp<2>(key);
226  bool locCheck = (fabs(computedValue - validationValue) < TOL);
227  check = check && locCheck;
228  if (!locCheck)
229  {
230  std::cout << "FAILED CHECK :: pos=" << Point<2, double>(domain.getPos(key)).toString()
231  << ", computedValue=" << computedValue
232  << ", validationValue=" << validationValue
233  << ", difference=" << fabs(computedValue - validationValue)
234  << ", tolerance=" << TOL
235  << std::endl;
236 // break;
237  }
238  ++itVal;
239  }
240  BOOST_REQUIRE(check);
241  }
242 
243  BOOST_AUTO_TEST_CASE(Dcpse_3D_test)
244  {
245  int rank;
246  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
247 
248  BOOST_TEST_MESSAGE("Init vars...");
249 
250  // Here build some easy domain and get some points around a given one
251  size_t edgeSemiSize = 20;
252  const unsigned int DIM = 3;
253  const size_t sz[DIM] = {2 * edgeSemiSize, 2 * edgeSemiSize, 2 * edgeSemiSize};
254  Box<DIM, double> box({0, 0, 0}, {2 * M_PI, 2 * M_PI, 2 * M_PI});
255  size_t bc[DIM] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
256  double spacing[DIM];
257  spacing[0] = 1.0 / (sz[0] - 1);
258  spacing[1] = 1.0 / (sz[1] - 1);
259  spacing[2] = 1.0 / (sz[2] - 1);
260  Ghost<DIM, double> ghost(0.1);
261 
262  double rCut = 3.1 * spacing[0];
263 
264  BOOST_TEST_MESSAGE("Init vector_dist...");
265  vector_dist<DIM, double, aggregate<double, double, double>> domain(0, box, bc, ghost);
266 
267  if (rank == 0)
268  {
269  BOOST_TEST_MESSAGE("Init domain...");
270  auto it = domain.getGridIterator(sz);
271  size_t pointId = 0;
272  size_t counter = 0;
273  double minNormOne = 999;
274  while (it.isNext())
275  {
276  domain.add();
277  auto key = it.get();
278  mem_id k0 = key.get(0);
279  double x = k0 * spacing[0];
280  domain.getLastPos()[0] = x;
281  mem_id k1 = key.get(1);
282  double y = k1 * spacing[1];
283  domain.getLastPos()[1] = y;
284  mem_id k2 = key.get(2);
285  double z = k2 * spacing[2];
286  domain.getLastPos()[2] = z;
287  // Here fill the function value
288  domain.template getLastProp<0>() = sin(z);
289 // domain.template getLastProp<0>() = x * x;
290 // domain.template getLastProp<0>() = x;
291  // Here fill the validation value for Df/Dx
292  domain.template getLastProp<2>() = cos(z);
293 // domain.template getLastProp<2>() = 2 * x;
294 // domain.template getLastProp<2>() = 1;
295 
296  ++counter;
297  ++it;
298  }
299  BOOST_TEST_MESSAGE("Sync domain across processors...");
300  }
301  domain.map(); // Send particles to the right processors
302 
303  // Setup finished, actual test below...
304  BOOST_TEST_MESSAGE("DCPSE init & compute coefficients...");
305  auto verletList = domain.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
306  Dcpse<DIM, decltype(verletList), vector_dist<DIM, double, aggregate<double, double, double>>> dcpse(domain, verletList, Point<DIM, unsigned int>({0, 0, 1}), 2, rCut, support_options::RADIUS);
307  BOOST_TEST_MESSAGE("DCPSE compute diff operator...");
308  dcpse.template computeDifferentialOperator<0, 1>(domain);
309 
310  // Now check against the validation values
311  BOOST_TEST_MESSAGE("Validating against ground truth...");
312 // const double TOL = 1e-6;
313  const double avgSpacing = spacing[0] + spacing[1] + spacing[2];
314  const double TOL = avgSpacing * avgSpacing;
315  auto itVal = domain.getDomainIterator();
316  bool check = true;
317  double computedValue;
318  double validationValue;
319  while (itVal.isNext())
320  {
321  auto key = itVal.get();
322  computedValue = domain.template getProp<1>(key);
323  validationValue = domain.template getProp<2>(key);
324  check = check && (fabs(computedValue - validationValue) < TOL);
325  if (!check)
326  {
327  std::cout << "FAILED CHECK :: pos=" << Point<DIM, double>(domain.getPos(key)).toString()
328  << ", computedValue=" << computedValue
329  << ", validationValue=" << validationValue
330  << ", tolerance=" << TOL
331  << std::endl;
332  break;
333  }
334  ++itVal;
335  }
336  BOOST_REQUIRE(check);
337  }
338 
339 #endif // HAVE_EIGEN
340 
341 BOOST_AUTO_TEST_SUITE_END()
Point::toString
std::string toString() const
Return the string with the point coordinate.
Definition: Point.hpp:413
Box
This class represent an N-dimensional box.
Definition: Box.hpp:59
vector_dist
Distributed vector.
Definition: vector_dist.hpp:170
Point
This class implement the point shape in an N-dimensional space.
Definition: Point.hpp:27
Ghost
Definition: Ghost.hpp:39