OpenFPM  5.2.0
Project that contain the implementation of distributed structures
OdeIntegrator_grid_tests.cpp
1 //
2 // Created by foggia on 19th Jan 2022
3 // It's a modification of Abhinav's test, adapted for grids
4 //
5 
6 #define BOOST_TEST_DYN_LINK
7 
8 #include <iostream>
9 #include <boost/test/unit_test.hpp>
10 
11 #include "config.h"
12 #include "Grid/grid_dist_id.hpp"
13 #include "OdeIntegrators/OdeIntegrators.hpp"
14 #include "FiniteDifference/FD_op.hpp"
15 #include "util/util_debug.hpp"
16 #include "util/common.hpp"
17 
18 const double a = 2.8e-4;
19 const double b = 5e-3;
20 const double tau = .1;
21 const double k = .005;
22 const int dim = 2;
23 
24 void *gridGlobal;
26 
27 // State types for systems with different number of ODEs
31 
32 template<typename DX, typename DY>
33 struct Fitz {
34 
35  DX & ddx;
36  DY & ddy;
37 
38  //Constructor
39  Fitz(DX & m_ddx, DY & m_ddy) :
40  ddx(m_ddx),
41  ddy(m_ddy) {}
42 
43  void operator()(const state_type_2ode & x,
44  state_type_2ode & dxdt,
45  const double t) const {
46 
47  grid_type & temp = *(grid_type *) gridGlobal;
48  auto u{FD::getV<4>(temp)};
49  auto v{FD::getV<5>(temp)};
50  u = x.data.get<0>();
51  v = x.data.get<1>();
52 
53  temp.ghost_get<4,5>();
54  dxdt.data.get<0>() = ddx(u) + ddy(u) + (1.0);
55  dxdt.data.get<1>() = ddx(v) + ddy(v) + (2.0);
56 
57  // One point stay fixed
58 
59  auto key2 = temp.getDomainIterator().get();
60 
61  if (create_vcluster().rank() == 0) {
62  dxdt.data.get<0>().value_ref(key2) = 0.0;
63  dxdt.data.get<1>().value_ref(key2) = 0.0;
64  }
65 
66  double u_max{0.0};
67  double v_max{0.0};
68 
69  auto it = temp.getDomainIterator();
70 
71  while (it.isNext())
72  {
73  auto key = it.get();
74 
75  if (u_max < dxdt.data.get<0>().value(key))
76  u_max = dxdt.data.get<0>().value(key);
77 
78  if (v_max < dxdt.data.get<1>().value(key))
79  v_max = dxdt.data.get<1>().value(key);
80 
81  ++it;
82  }
83  }
84 };
85 
86 void Exponential_struct_ofp2(const state_type_3ode & x,
87  state_type_3ode & dxdt,
88  const double t) {
89 
90  // sytem: dx1/dt = x1 --> solution: x1(t) = exp(t)
91  // sytem: dx2/dt = 2*x2 --> solution: x2(t) = exp(2t)
92  dxdt.data.get<0>() = x.data.get<0>();
93  dxdt.data.get<1>() = 2.0 * x.data.get<1>();
94  dxdt.data.get<2>() = x.data.get<0>();
95 }
96 
97 
98 void Exponential(const state_type_1ode & x,
99  state_type_1ode & dxdt,
100  const double t) {
101 
102  // sytem: dx/dt = x --> solution: x(t) = exp(t)
103  dxdt = x;
104 }
105 
106 // void sigmoid(const state_type_1ode & x,
107 // state_type_1ode & dxdt,
108 // const double t) {
109 // dxdt = x * (1.0 - x);
110 // }
111 
112 BOOST_AUTO_TEST_SUITE(odeInt_grid_tests)
113 
114 BOOST_AUTO_TEST_CASE(odeint_grid_test_exponential) {
115 
116  size_t edgeSemiSize{40};
117  const size_t sz[dim] = {edgeSemiSize,edgeSemiSize};
118  Box<dim,double> box{{0.0,0.0}, {1.0,1.0}};
119  periodicity<dim> bc{{NON_PERIODIC,NON_PERIODIC}};
120  double spacing[dim];
121  spacing[0] = 1.0 / (sz[0] - 1);
122  spacing[1] = 1.0 / (sz[1] - 1);
123  Ghost<dim,long int> ghost{2};
124  BOOST_TEST_MESSAGE("Test: exponential");
125  BOOST_TEST_MESSAGE("Init grid_dist_id ...");
126 
128 
129  auto it{grid.getDomainIterator()};
130  while (it.isNext()) {
131  auto key = it.get();
132  grid.template get<0>(key) = std::exp(0); // Initial state
133  grid.template get<1>(key) = std::exp(0.4); // Analytical solution
134  ++it;
135  }
136  grid.ghost_get<0>();
137 
138  auto Init{FD::getV<0>(grid)}; // Initial state
139  auto Sol{FD::getV<1>(grid)}; // Analytical solution
140  auto OdeSol{FD::getV<2>(grid)}; // Numerical solution
141 
142  state_type_1ode x0;
143  x0.data.get<0>() = Init;
144 
145  double t{0.0};
146  double tf{0.4};
147  const double dt{0.1};
148 
149  boost::numeric::odeint::runge_kutta4<state_type_1ode,double,
150  state_type_1ode,double,
151  boost::numeric::odeint::vector_space_algebra_ofp> rk4; // Time integrator
152  size_t steps{boost::numeric::odeint::integrate_const(rk4,Exponential,x0,0.0,tf,dt)};
153 
154  OdeSol = x0.data.get<0>(); // Numerical solution
155 
156  // Error
157  auto it2{grid.getDomainIterator()};
158  double worst{0.0};
159  while (it2.isNext()) {
160  auto p{it2.get()};
161  if (std::fabs(grid.template get<1>(p) - grid.template get<2>(p)) > worst) {
162  worst = std::fabs(grid.template get<1>(p) - grid.template get<2>(p));
163  }
164  ++it2;
165  }
166 
167  std::cout << worst << std::endl;
168  BOOST_REQUIRE(worst < 1e-6);
169 
170  // Another way
171  x0.data.get<0>() = Init;
172  while (t < tf) {
173  rk4.do_step(Exponential,x0,t,dt);
174  OdeSol = x0.data.get<0>();
175  t += dt;
176  }
177 
178  OdeSol = x0.data.get<0>();
179 
180  // Error
181  auto it3{grid.getDomainIterator()};
182  double worst2{0.0};
183  while (it3.isNext()) {
184  auto p{it3.get()};
185  if (std::fabs(grid.template get<1>(p) - grid.template get<2>(p)) > worst2) {
186  worst2 = fabs(grid.template get<1>(p) - grid.template get<2>(p));
187  }
188  ++it3;
189  }
190  std::cout << worst2 << std::endl;
191  BOOST_REQUIRE(worst2 < 1e-6);
192  BOOST_REQUIRE_EQUAL(worst,worst2);
193 }
194 
195 
196 
197 BOOST_AUTO_TEST_CASE(odeint_grid_test_STRUCT_exponential) {
198 
199  size_t edgeSemiSize{40};
200  const size_t sz[dim] = {edgeSemiSize,edgeSemiSize};
201  Box<dim, double> box{{0.0,0.0},{1.0,1.0}};
202  periodicity<dim> bc{{NON_PERIODIC,NON_PERIODIC}};
203  double spacing[dim];
204  spacing[0] = 1.0 / (sz[0] - 1);
205  spacing[1] = 1.0 / (sz[1] - 1);
206  Ghost<dim,long int> ghost{2};
207  BOOST_TEST_MESSAGE("Test: exponential");
208  BOOST_TEST_MESSAGE("Init grid_dist_id ...");
209 
211 
212  auto it{grid.getDomainIterator()};
213  while (it.isNext()) {
214  auto key = it.get();
215 
216  grid.get<0>(key) = std::exp(0);
217  grid.template get<0>(key) = std::exp(0.0); // Initial state 1
218  grid.template get<1>(key) = std::exp(0.4); // Analytical solution 1
219  grid.template get<2>(key) = std::exp(0.0); // Initial state 2
220  grid.template get<3>(key) = std::exp(0.8); // Analytical solution 2
221  ++it;
222  }
223  grid.ghost_get<0>();
224 
225  auto Init1{FD::getV<0>(grid)}; // Initial state 1
226  auto Sol1{FD::getV<1>(grid)}; // Analytical solution 1
227 
228  auto Init2{FD::getV<2>(grid)}; // Initial state 2
229  auto Sol2{FD::getV<3>(grid)}; // Analytical solution 2
230 
231  auto OdeSol1{FD::getV<4>(grid)}; // Numerical solution 1
232  auto OdeSol2{FD::getV<5>(grid)}; // Numerical solution 2
233 
234  state_type_3ode x0;
235  x0.data.get<0>() = Init1;
236  x0.data.get<1>() = Init2;
237  x0.data.get<2>() = Init1;
238 
239  double t{0};
240  double tf{0.4};
241  const double dt{0.1};
242 
243  // size_t steps{boost::numeric::odeint::integrate(Exponential_struct,x0,0.0,tf,dt)};
244  // size_t steps{boost::numeric::odeint::integrate_const(boost::numeric::odeint::runge_kutta4<state_type_3ode,double,state_type_3ode,double,boost::numeric::odeint::vector_space_algebra_ofp>(),
245  // Exponential_struct_ofp,x0,0.0,tf,dt)};
246 
247  typedef boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_cash_karp54<state_type_3ode,double,state_type_3ode,double,boost::numeric::odeint::vector_space_algebra_ofp>> stepper_type;
248  integrate_adaptive(stepper_type(),Exponential_struct_ofp2,x0,t,tf,dt);
249 
250  OdeSol1 = x0.data.get<0>();
251  OdeSol2 = x0.data.get<1>();
252 
253  // Error
254  auto it2{grid.getDomainIterator()};
255  double worst{0.0};
256  double worst2{0.0};
257  while (it2.isNext()) {
258  auto p{it2.get()};
259  if (std::fabs(grid.getProp<1>(p) - grid.getProp<4>(p)) > worst)
260  worst = std::fabs(grid.getProp<1>(p) - grid.getProp<4>(p));
261  if (std::fabs(grid.getProp<3>(p) - grid.getProp<5>(p)) > worst2)
262  worst2 = std::fabs(grid.getProp<3>(p) - grid.getProp<5>(p));
263  ++it2;
264  }
265 
266  std::cout << worst << " " << worst2 << std::endl;
267  BOOST_REQUIRE(worst < 1e-6);
268  BOOST_REQUIRE(worst2 < 1e-6);
269 
270 
271  // A different way
272  x0.data.get<0>() = Init1;
273  x0.data.get<1>() = Init2;
274  x0.data.get<2>() = Init1;
275 
276  boost::numeric::odeint::runge_kutta4<state_type_3ode,double,state_type_3ode,double,boost::numeric::odeint::vector_space_algebra_ofp> rk4;
277  while (t < tf) {
278  rk4.do_step(Exponential_struct_ofp2,x0,t,dt);
279  t+=dt;
280  }
281 
282  OdeSol1 = x0.data.get<0>();
283  OdeSol2 = x0.data.get<1>();
284 
285  // Error
286  auto it3{grid.getDomainIterator()};
287  double worst3{0.0};
288  double worst4{0.0};
289  while (it3.isNext()) {
290  auto p{it3.get()};
291  if (std::fabs(grid.getProp<1>(p) - grid.getProp<4>(p)) > worst3)
292  worst3 = std::fabs(grid.getProp<1>(p) - grid.getProp<4>(p));
293  if (std::fabs(grid.getProp<3>(p) - grid.getProp<5>(p)) > worst4)
294  worst4 = std::fabs(grid.getProp<3>(p) - grid.getProp<5>(p));
295  ++it3;
296  }
297 
298  std::cout << worst3 << " " << worst4 << std::endl;
299  BOOST_REQUIRE(worst3 < 1e-6);
300  BOOST_REQUIRE(worst4 < 5e-5);
301 }
302 
303 
304 
305 
306 BOOST_AUTO_TEST_CASE(odeint_grid_test2_exponential) {
307 
308  size_t edgeSemiSize{40};
309  const size_t sz[dim] = {edgeSemiSize,edgeSemiSize};
310  Box<dim,double> box{{0.0,0.0}, {1.0,1.0}};
311  periodicity<dim> bc{{NON_PERIODIC,NON_PERIODIC}};
312  double spacing[dim];
313  spacing[0] = 1.0 / (sz[0] - 1);
314  spacing[1] = 1.0 / (sz[1] - 1);
315  Ghost<dim,long int> ghost{2};
316  BOOST_TEST_MESSAGE("Test: exponential");
317  BOOST_TEST_MESSAGE("Init grid_dist_id ...");
318 
320 
321  double t{0.0};
322  double tf{0.5};
323  const double dt{0.1};
324 
325  auto it{grid.getDomainIterator()};
326  while (it.isNext()) {
327  auto key = it.get();
328  grid.template get<0>(key) = std::exp(t); // Initial state
329  grid.template get<1>(key) = std::exp(tf); // Analytical solution
330  ++it;
331  }
332  grid.ghost_get<0>();
333 
334  auto Init{FD::getV<0>(grid)}; // Initial state
335  auto Sol{FD::getV<1>(grid)}; // Analytical solution
336  auto OdeSol{FD::getV<2>(grid)}; // Numerical solution
337 
338  state_type_1ode x0;
339  x0.data.get<0>() = Init;
340 
341  typedef boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_cash_karp54<state_type_1ode,double,state_type_1ode,double,boost::numeric::odeint::vector_space_algebra_ofp>> stepper_type;
342 
343 
344  integrate_adaptive(stepper_type(),Exponential,x0,t,tf,dt);
345  OdeSol = x0.data.get<0>(); // Numerical solution
346 
347  // Error
348  auto it2{grid.getDomainIterator()};
349  double worst{0.0};
350  while (it2.isNext()) {
351  auto p{it2.get()};
352  if (std::fabs(grid.template get<1>(p) - grid.template get<2>(p)) > worst) {
353  worst = std::fabs(grid.template get<1>(p) - grid.template get<2>(p));
354  }
355  ++it2;
356  }
357 
358  std::cout << worst << std::endl;
359  BOOST_REQUIRE(worst < 1e-6);
360 
361  // Another way
362  boost::numeric::odeint::runge_kutta4<state_type_1ode,double,state_type_1ode,double,boost::numeric::odeint::vector_space_algebra_ofp> rk4;
363  x0.data.get<0>() = Init;
364  for (size_t i = 0; i < int(tf/dt); ++i, t += dt) {
365  rk4.do_step(Exponential,x0,t,dt);
366  t += dt;
367  }
368  OdeSol = x0.data.get<0>();
369 
370  // Error
371  auto it3{grid.getDomainIterator()};
372  double worst2{0.0};
373  while (it3.isNext()) {
374  auto p{it3.get()};
375  if (std::fabs(grid.template get<1>(p) - grid.template get<2>(p)) > worst2) {
376  worst2 = fabs(grid.template get<1>(p) - grid.template get<2>(p));
377  }
378  ++it3;
379  }
380  std::cout << worst2 << std::endl;
381  BOOST_REQUIRE(worst2 < 1e-6);
382 
383  // Yet another way
384  // x0.data.get<0>() = Init;
385  // integrate(rk4,Exponential,x0,t,tf,dt);
386 
387  // OdeSol = x0.data.get<0>();
388 
389  // // Error
390  // auto it4{grid.getDomainIterator()};
391  // double worst3{0.0};
392  // while (it4.isNext()) {
393  // auto p{it4.get()};
394  // if (std::fabs(grid.template get<1>(p) - grid.template get<2>(p)) > worst3) {
395  // worst3 = fabs(grid.template get<1>(p) - grid.template get<2>(p));
396  // }
397  // ++it4;
398  // }
399  // std::cout << worst3 << std::endl;
400  // BOOST_REQUIRE(worst3 < 1e-6);
401 
402  // BOOST_REQUIRE_EQUAL(worst,worst2);
403  // BOOST_REQUIRE_EQUAL(worst2,worst3);
404 }
405 
406 
407 
408 // BOOST_AUTO_TEST_CASE(odeint_base_test3)
409 // {
410 // size_t edgeSemiSize = 40;
411 // const size_t sz[2] = {edgeSemiSize,edgeSemiSize };
412 // Box<2, double> box({ 0, 0 }, { 1.0, 1.0 });
413 // size_t bc[2] = {NON_PERIODIC,NON_PERIODIC};
414 // double spacing[2];
415 // spacing[0] = 1.0 / (sz[0] - 1);
416 // spacing[1] = 1.0 / (sz[1] - 1);
417 // double rCut = 3.9 * spacing[0];
418 // Ghost<2, double> ghost(rCut);
419 // BOOST_TEST_MESSAGE("Init vector_dist...");
420 
421 // vector_dist<2, double, aggregate<double, double,double>> Particles(0, box, bc, ghost);
422 
423 // double t=0.0,tf=0.5;
424 // const double dt=0.1;
425 
426 // auto it = Particles.getGridIterator(sz);
427 // while (it.isNext())
428 // {
429 // Particles.add();
430 // auto key = it.get();
431 // mem_id k0 = key.get(0);
432 // double xp0 = k0 * spacing[0];
433 // Particles.getLastPos()[0] = xp0;
434 // mem_id k1 = key.get(1);
435 // double yp0 = k1 * spacing[1];
436 // Particles.getLastPos()[1] = yp0;
437 // Particles.getLastProp<0>() = 1.0/(1.0+exp(-t)); // Carefull in putting the constant, f = A*sigmoid does not respect f' = f*(1.0-f) but f*(1.0-f/A), for simplicity I remove the constant
438 // Particles.getLastProp<1>() = 1.0/(1.0+exp(-tf)); // Carefull in putting the constant, f = A*sigmoid does not respect f' = f*(1.0-f) but f*(1.0-f/A), for simplicity I remove the constant
439 // ++it;
440 // }
441 // Particles.map();
442 // Particles.ghost_get<0>();
443 // auto Init = getV<0>(Particles);
444 // auto Sol = getV<1>(Particles);
445 // auto OdeSol = getV<2>(Particles);
446 
447 // state_type x0;
448 // x0=Init;
449 // // The rhs of x' = f(x)
450 // //size_t steps=boost::numeric::odeint::integrate(sigmoid,x0,0.0,tf,dt);
451 // //typedef boost::numeric::odeint::controlled_runge_kutta< boost::numeric::odeint::runge_kutta_cash_karp54< state_type > > stepper_type;
452 // //integrate_adaptive( stepper_type() , sigmoid , x0 , t , tf , dt);
453 // size_t steps=boost::numeric::odeint::integrate_const( boost::numeric::odeint::runge_kutta4< state_type >(),sigmoid,x0,0.0,tf,dt);
454 
455 // OdeSol=x0;
456 // auto it2 = Particles.getDomainIterator();
457 // double worst = 0.0;
458 // while (it2.isNext()) {
459 // auto p = it2.get();
460 // if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst) {
461 // worst = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
462 // }
463 // ++it2;
464 // }
465 
466 // BOOST_REQUIRE(worst < 1e-8);
467 
468 // x0=Init;
469 // boost::numeric::odeint::runge_kutta4< state_type > rk4;
470 // for( size_t i=0 ; i<int(tf/dt) ; ++i,t+=dt )
471 // {
472 // rk4.do_step(sigmoid,x0,t,dt);
473 // t+=dt;
474 // }
475 
476 // OdeSol=x0;
477 // auto it3 = Particles.getDomainIterator();
478 // double worst2 = 0.0;
479 // while (it3.isNext()) {
480 // auto p = it3.get();
481 // if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst2) {
482 // worst2 = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
483 // }
484 // ++it3;
485 // }
486 // //std::cout<<worst2<<std::endl;
487 // BOOST_REQUIRE(worst < 1e-6);
488 // BOOST_REQUIRE_EQUAL(worst,worst2);
489 // }
490 
491 
492 
493 #ifdef HAVE_EIGEN
494 
495 BOOST_AUTO_TEST_CASE(dcpse_op_react_diff_test) {
496 
497  size_t edgeSemiSize{5};
498  const size_t sz[dim] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
499  Box<dim,double> box{{0.0, 0.0},{1.0, 1.0}};
500  periodicity<dim> bc{{PERIODIC,PERIODIC}};
501  double spacing[dim];
502  spacing[0] = 1.0 / (sz[0]);
503  spacing[1] = 1.0 / (sz[1]);
504  Ghost<dim,double> ghost{spacing[0] * 3};
505 
506  BOOST_TEST_MESSAGE("Test: reaction diffusion");
507  BOOST_TEST_MESSAGE("Init grid_dist_id ...");
508  double sigma2 = spacing[0] * spacing[1] / (2 * 4);
509 
510  // properties: u, v, du, dv
512 
513  auto it{domain.getDomainIterator()};
514  while (it.isNext()) {
515  auto key{it.get()};
516  domain.get<0>(key) = 0.0; // u
517  domain.get<1>(key) = 0.0; // v
518  domain.get<2>(key) = 0.0; // du/dt
519  domain.get<3>(key) = 0.0; // dv/dt
520 
521  auto gkey = it.getGKey(key);
522 
523  if (gkey.get(0)==sz[0] / 2 && gkey.get(1) == sz[1]/2)
524  {
525  domain.get<0>(key) = 1.0;
526  domain.get<1>(key) = 1.0;
527  }
528  ++it;
529  }
530  domain.ghost_get<0>();
531 
534 
535  gridGlobal=(void *) & domain;
536 
537  auto u{FD::getV<0>(domain)};
538  auto v{FD::getV<1>(domain)};
539  auto fu{FD::getV<2>(domain)};
540  auto fv{FD::getV<3>(domain)};
541 
542  Fitz<decltype(ddx),decltype(ddy)> system(ddx,ddy);
543  state_type_2ode x0;
544 
545  x0.data.get<0>() = u;
546  x0.data.get<1>() = v;
547 
548  double dt{0.001};
549  double t{0.0};
550  double tf{10.5};
551 
552  //typedef boost::numeric::odeint::controlled_runge_kutta< boost::numeric::odeint::runge_kutta_cash_karp54< state_type_2d_ofp,double,state_type_2d_ofp,double,boost::numeric::odeint::vector_space_algebra_ofp>> stepper_type;
553  typedef boost::numeric::odeint::runge_kutta4<state_type_2ode,double,state_type_2ode,double,boost::numeric::odeint::vector_space_algebra_ofp> stepper_type;
554 
555  integrate_adaptive(stepper_type(),system,x0,t,tf,dt);
556  fu = x0.data.get<0>();
557  fv = x0.data.get<1>();
558 
559  domain.ghost_get<2,3>();
560  u = ddx(fu) + ddy(fu);
561  v = ddx(fv) + ddy(fv);
562 
563  auto it2{domain.getDomainIterator()};
564 
565  if (create_vcluster().rank() == 0)
566  ++it2;
567 
568  while (it2.isNext()) {
569  auto p{it2.get()};
570  BOOST_REQUIRE_CLOSE(domain.get<0>(p),-1.0,1);
571  ++it2;
572  }
573 }
574 #endif
575 
576 BOOST_AUTO_TEST_SUITE_END()
This class represent an N-dimensional box.
Definition: Box.hpp:60
This is a distributed grid.
void ghost_get(size_t opt=0)
It synchronize the ghost parts.
grid_dist_iterator< dim, device_grid, decltype(device_grid::type_of_subiterator()), FREE > getDomainIterator() const
It return an iterator that span the full grid domain (each processor span its local domain)
KeyT const ValueT ValueT OffsetIteratorT OffsetIteratorT int
[in] The number of segments that comprise the sorting data
Boundary conditions.
Definition: common.hpp:22