OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
OdeIntegratores_base_tests.cpp
1 //
2 // Created by Abhinav Singh on 01.12.20.
3 //
4 #include "config.h"
5 #include <type_traits>
6 #include <cstring>
7 #include "util/common.hpp"
8 
9 #define BOOST_TEST_DYN_LINK
10 
11 #include "util/util_debug.hpp"
12 #include <boost/test/unit_test.hpp>
13 #include <iostream>
14 #include "Operators/Vector/vector_dist_operators.hpp"
15 #include "Vector/vector_dist_subset.hpp"
16 #include "Decomposition/Distribution/SpaceDistribution.hpp"
17 #include "OdeIntegrators/OdeIntegrators.hpp"
18 #include "DCPSE/DCPSE_op/DCPSE_op.hpp"
19 #include "OdeIntegrators/boost_vector_algebra_ofp.hpp"
20 
22 const double a = 2.8e-4;
23 const double b = 5e-3;
24 const double tau = .1;
25 const double k = .005;
26 
27 void *vectorGlobal;
29 
30 template<typename op_type>
31 struct Fitz
32 {
33  op_type &Lap;
34 
35 
36  //Constructor
37  Fitz(op_type &Lap):Lap(Lap)
38  {}
39 
40  void operator()( const state_type_2d_ofp &x , state_type_2d_ofp &dxdt , const double t ) const
41  {
42  vector_type &Temp= *(vector_type *) vectorGlobal;
43  auto u=getV<4>(Temp);
44  auto v=getV<5>(Temp);
45  u=x.data.get<0>();
46  v=x.data.get<1>();
47 
48  Temp.ghost_get<4,5>();
49  dxdt.data.get<0>() = Lap(u) + (1.0);
50  dxdt.data.get<1>() = Lap(v) + (2.0);
51 
52  // One point stay fixed
53  if (create_vcluster().rank() == 0)
54  {
55  dxdt.data.get<0>().getVector().get<0>(0) = 0.0;
56  dxdt.data.get<1>().getVector().get<0>(0) = 0.0;
57  }
58 
59  double u_max = 0.0;
60  double v_max = 0.0;
61 
62  for (int i = 0 ; i < dxdt.data.get<0>().getVector().size(); i++)
63  {
64  if (u_max < dxdt.data.get<0>().getVector().get<0>(i))
65  {u_max = dxdt.data.get<0>().getVector().get<0>(i);}
66 
67  if (v_max < dxdt.data.get<1>().getVector().get<0>(i))
68  {v_max = dxdt.data.get<1>().getVector().get<0>(i);}
69  }
70  }
71 };
72 
73 void Exponential_struct_ofp( const state_type_2d_ofp &x , state_type_2d_ofp &dxdt , const double t )
74 {
75  dxdt.data.get<0>() = x.data.get<0>();
76  dxdt.data.get<1>() = 2.0*x.data.get<1>();
77 }
78 void Exponential_struct_ofp2( const state_type_3d_ofp &x , state_type_3d_ofp &dxdt , const double t )
79 {
80  dxdt.data.get<0>() = x.data.get<0>();
81  dxdt.data.get<1>() = 2.0*x.data.get<1>();
82 }
83 
84 
85 void Exponential( const state_type &x , state_type &dxdt , const double t )
86 {
87  dxdt = x;
88 }
89 void sigmoid( const state_type &x , state_type &dxdt , const double t )
90 {
91  dxdt = x*(1.0-x);
92 }
93 
94 BOOST_AUTO_TEST_SUITE(odeInt_BASE_tests)
95 
96 BOOST_AUTO_TEST_CASE(odeint_base_test1)
97 {
98  size_t edgeSemiSize = 40;
99  const size_t sz[2] = {edgeSemiSize,edgeSemiSize };
100  Box<2, double> box({ 0, 0 }, { 1.0, 1.0 });
101  size_t bc[2] = { NON_PERIODIC, NON_PERIODIC };
102  double spacing[2];
103  spacing[0] = 1.0 / (sz[0] - 1);
104  spacing[1] = 1.0 / (sz[1] - 1);
105  double rCut = 3.9 * spacing[0];
106  Ghost<2, double> ghost(rCut);
107  BOOST_TEST_MESSAGE("Init vector_dist...");
108 
109  vector_dist<2, double, aggregate<double, double,double>> Particles(0, box, bc, ghost);
110 
111  auto it = Particles.getGridIterator(sz);
112  while (it.isNext())
113  {
114  Particles.add();
115  auto key = it.get();
116  mem_id k0 = key.get(0);
117  double xp0 = k0 * spacing[0];
118  Particles.getLastPos()[0] = xp0;
119  mem_id k1 = key.get(1);
120  double yp0 = k1 * spacing[1];
121  Particles.getLastPos()[1] = yp0;
122  Particles.getLastProp<0>() = xp0*yp0*exp(0);
123  Particles.getLastProp<1>() = xp0*yp0*exp(0.4);
124  ++it;
125  }
126 
127  Particles.map();
128  Particles.ghost_get<0>();
129  auto Init = getV<0>(Particles);
130  auto Sol = getV<1>(Particles);
131  auto OdeSol = getV<2>(Particles);
132 
133  state_type x0;
134  x0=Init;
135  // The rhs of x' = f(x)
136 
137  double t=0,tf=0.4;
138  const double dt=0.1;
139 
140  //This doesnt work Why?
141  //size_t steps=boost::numeric::odeint::integrate(Exponential,x0,0.0,tf,dt);
142 
143  size_t steps=boost::numeric::odeint::integrate_const( boost::numeric::odeint::runge_kutta4< state_type >(),Exponential,x0,0.0,tf,dt);
144 
145  OdeSol=x0;
146  auto it2 = Particles.getDomainIterator();
147  double worst = 0.0;
148  while (it2.isNext()) {
149  auto p = it2.get();
150  if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst) {
151  worst = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
152  }
153  ++it2;
154  }
155 
156  //std::cout<<worst<<std::endl;
157  BOOST_REQUIRE(worst < 1e-6);
158 
159  x0=Init;
160  boost::numeric::odeint::runge_kutta4< state_type > rk4;
161  while (t<tf)
162  {
163  rk4.do_step(Exponential,x0,t,dt);
164  OdeSol=x0;
165  t+=dt;
166  }
167 
168  OdeSol=x0;
169  auto it3 = Particles.getDomainIterator();
170  double worst2 = 0.0;
171  while (it3.isNext()) {
172  auto p = it3.get();
173  if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst2) {
174  worst2 = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
175  }
176  ++it3;
177  }
178  //std::cout<<worst2<<std::endl;
179  BOOST_REQUIRE(worst < 1e-6);
180  BOOST_REQUIRE_EQUAL(worst,worst2);
181 }
182 
183 BOOST_AUTO_TEST_CASE(odeint_base_test_STRUCT_ofp)
184 {
185  size_t edgeSemiSize = 40;
186  const size_t sz[2] = {edgeSemiSize,edgeSemiSize };
187  Box<2, double> box({ 0, 0 }, { 1.0, 1.0 });
188  size_t bc[2] = { NON_PERIODIC, NON_PERIODIC };
189  double spacing[2];
190  spacing[0] = 1.0 / (sz[0] - 1);
191  spacing[1] = 1.0 / (sz[1] - 1);
192  double rCut = 3.9 * spacing[0];
193  Ghost<2, double> ghost(rCut);
194  BOOST_TEST_MESSAGE("Init vector_dist...");
195 
197 
198  auto it = Particles.getGridIterator(sz);
199  while (it.isNext())
200  {
201  Particles.add();
202  auto key = it.get();
203  mem_id k0 = key.get(0);
204  double xp0 = k0 * spacing[0];
205  Particles.getLastPos()[0] = xp0;
206  mem_id k1 = key.get(1);
207  double yp0 = k1 * spacing[1];
208  Particles.getLastPos()[1] = yp0;
209  Particles.getLastProp<0>() = xp0*yp0*exp(0);
210  Particles.getLastProp<1>() = xp0*yp0*exp(0.4);
211  Particles.getLastProp<2>() = xp0*yp0*exp(0);
212  Particles.getLastProp<3>() = xp0*yp0*exp(0.8);
213  ++it;
214  }
215  Particles.map();
216  Particles.ghost_get<0>();
217  auto Init1 = getV<0>(Particles);
218  auto Sol1 = getV<1>(Particles);
219  auto Init2 = getV<2>(Particles);
220  auto Sol2 = getV<3>(Particles);
221  auto OdeSol1 = getV<4>(Particles);
222  auto OdeSol2 = getV<5>(Particles);
223 
224  state_type_3d_ofp x0;//(Init1,Init2);
225  x0.data.get<0>()=Init1;
226  x0.data.get<1>()=Init2;
227  x0.data.get<2>()=Init1;
228 
229  // The rhs of x' = f(x)
230  double t=0,tf=0.4;
231  const double dt=0.1;
232 
233  //size_t steps=boost::numeric::odeint::integrate(Exponential_struct,x0,0.0,tf,dt);
234 
235  //size_t steps=boost::numeric::odeint::integrate_const(
236  // boost::numeric::odeint::runge_kutta4< state_type_struct_ofp<vector_type>,double,state_type_struct_ofp<vector_type>,double,boost::numeric::odeint::vector_space_algebra_ofp > ()
237  // ,Exponential_struct_ofp,x0,0.0,tf,dt);
238 
239  //typedef boost::numeric::odeint::controlled_runge_kutta< boost::numeric::odeint::runge_kutta_cash_karp54< state_type_struct_ofp<vector_type>,double,state_type_struct_ofp<vector_type>,double,boost::numeric::odeint::vector_space_algebra_ofp>> stepper_type;
240  typedef boost::numeric::odeint::controlled_runge_kutta< boost::numeric::odeint::runge_kutta_cash_karp54< state_type_3d_ofp,double,state_type_3d_ofp,double,boost::numeric::odeint::vector_space_algebra_ofp>> stepper_type;
241  integrate_adaptive( stepper_type() , Exponential_struct_ofp2 , x0 , t , tf , dt);
242 
243  OdeSol1=x0.data.get<0>();
244  OdeSol2=x0.data.get<1>();
245  auto it2 = Particles.getDomainIterator();
246  double worst = 0.0;
247  double worst2 = 0.0;
248  while (it2.isNext()) {
249  auto p = it2.get();
250  if (fabs(Particles.getProp<1>(p) - Particles.getProp<4>(p)) > worst) {
251  worst = fabs(Particles.getProp<1>(p) - Particles.getProp<4>(p));
252  }
253  if (fabs(Particles.getProp<3>(p) - Particles.getProp<5>(p)) > worst2) {
254  worst2 = fabs(Particles.getProp<3>(p) - Particles.getProp<5>(p));
255  }
256  ++it2;
257  }
258 
259  BOOST_REQUIRE(worst < 1e-6);
260  BOOST_REQUIRE(worst2 < 1e-6);
261 
262 
263  x0.data.get<0>()=Init1;
264  x0.data.get<1>()=Init2;
265  x0.data.get<2>()=Init1;
266 
267  //x0=Init;
268  boost::numeric::odeint::runge_kutta4< state_type_3d_ofp, double,state_type_3d_ofp,double,boost::numeric::odeint::vector_space_algebra_ofp> rk4;
269  while (t<tf)
270  {
271  rk4.do_step(Exponential_struct_ofp2,x0,t,dt);
272  t+=dt;
273  }
274 
275  OdeSol1=x0.data.get<0>();
276  OdeSol2=x0.data.get<1>();
277  auto it3 = Particles.getDomainIterator();
278  double worst3 = 0.0;
279  double worst4 = 0.0;
280  while (it3.isNext()) {
281  auto p = it3.get();
282  if (fabs(Particles.getProp<1>(p) - Particles.getProp<4>(p)) > worst3) {
283  worst3 = fabs(Particles.getProp<1>(p) - Particles.getProp<4>(p));
284  }
285  if (fabs(Particles.getProp<3>(p) - Particles.getProp<5>(p)) > worst4) {
286  worst4 = fabs(Particles.getProp<3>(p) - Particles.getProp<5>(p));
287  }
288  ++it3;
289  }
290 
291  BOOST_REQUIRE(worst3 < 1e-6);
292  BOOST_REQUIRE(worst4 < 5e-5);
293 }
294 
295 
296 BOOST_AUTO_TEST_CASE(odeint_base_test2)
297 {
298  size_t edgeSemiSize = 40;
299  const size_t sz[2] = {edgeSemiSize,edgeSemiSize };
300  Box<2, double> box({ 0, 0 }, { 1.0, 1.0 });
301  size_t bc[2] = { NON_PERIODIC, NON_PERIODIC };
302  double spacing[2];
303  spacing[0] = 1.0 / (sz[0] - 1);
304  spacing[1] = 1.0 / (sz[1] - 1);
305  double rCut = 3.9 * spacing[0];
306  Ghost<2, double> ghost(rCut);
307  BOOST_TEST_MESSAGE("Init vector_dist...");
308 
309  vector_dist<2, double, aggregate<double, double,double>> Particles(0, box, bc, ghost);
310 
311  double t=0.0,tf=0.5;
312  const double dt=0.1;
313 
314  auto it = Particles.getGridIterator(sz);
315  while (it.isNext())
316  {
317  Particles.add();
318  auto key = it.get();
319  mem_id k0 = key.get(0);
320  double xp0 = k0 * spacing[0];
321  Particles.getLastPos()[0] = xp0;
322  mem_id k1 = key.get(1);
323  double yp0 = k1 * spacing[1];
324  Particles.getLastPos()[1] = yp0;
325  Particles.getLastProp<0>() = xp0*yp0*exp(t);
326  Particles.getLastProp<1>() = xp0*yp0*exp(tf);
327  ++it;
328  }
329  Particles.map();
330  Particles.ghost_get<0>();
331  auto Init = getV<0>(Particles);
332  auto Sol = getV<1>(Particles);
333  auto OdeSol = getV<2>(Particles);
334 
335  state_type x0;
336  x0=Init;
337  // The rhs of x' = f(x)
338 
339 
340 
341 
342  //This doesnt work Why?
343  //size_t steps=boost::numeric::odeint::integrate(Exponential,x0,0.0,tf,dt);
344  typedef boost::numeric::odeint::controlled_runge_kutta< boost::numeric::odeint::runge_kutta_cash_karp54< state_type > > stepper_type;
345  integrate_adaptive( stepper_type() , Exponential , x0 , t , tf , dt);
346  //size_t steps=boost::numeric::odeint::integrate_const( boost::numeric::odeint::runge_kutta4< state_type >(),Exponential,x0,0.0,tf,dt);
347 
348  OdeSol=x0;
349  auto it2 = Particles.getDomainIterator();
350  double worst = 0.0;
351  while (it2.isNext()) {
352  auto p = it2.get();
353  if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst) {
354  worst = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
355  }
356  ++it2;
357  }
358 
359  BOOST_REQUIRE(worst < 1e-6);
360 
361  x0=Init;
362  boost::numeric::odeint::runge_kutta4< state_type > rk4;
363  for( size_t i=0 ; i<int(tf/dt) ; ++i,t+=dt )
364  {
365  rk4.do_step(Exponential,x0,t,dt);
366  t+=dt;
367  }
368 
369  OdeSol=x0;
370  auto it3 = Particles.getDomainIterator();
371  double worst2 = 0.0;
372  while (it3.isNext()) {
373  auto p = it3.get();
374  if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst2) {
375  worst2 = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
376  }
377  ++it3;
378  }
379  //std::cout<<worst2<<std::endl;
380  BOOST_REQUIRE(worst < 1e-6);
381  BOOST_REQUIRE(worst2 < 1e-6);
382 }
383 
384 BOOST_AUTO_TEST_CASE(odeint_base_test3)
385 {
386  size_t edgeSemiSize = 40;
387  const size_t sz[2] = {edgeSemiSize,edgeSemiSize };
388  Box<2, double> box({ 0, 0 }, { 1.0, 1.0 });
389  size_t bc[2] = { NON_PERIODIC, NON_PERIODIC };
390  double spacing[2];
391  spacing[0] = 1.0 / (sz[0] - 1);
392  spacing[1] = 1.0 / (sz[1] - 1);
393  double rCut = 3.9 * spacing[0];
394  Ghost<2, double> ghost(rCut);
395  BOOST_TEST_MESSAGE("Init vector_dist...");
396 
397  vector_dist<2, double, aggregate<double, double,double>> Particles(0, box, bc, ghost);
398 
399  double t=0.0,tf=0.5;
400  const double dt=0.1;
401 
402  auto it = Particles.getGridIterator(sz);
403  while (it.isNext())
404  {
405  Particles.add();
406  auto key = it.get();
407  mem_id k0 = key.get(0);
408  double xp0 = k0 * spacing[0];
409  Particles.getLastPos()[0] = xp0;
410  mem_id k1 = key.get(1);
411  double yp0 = k1 * spacing[1];
412  Particles.getLastPos()[1] = yp0;
413  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
414  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
415  ++it;
416  }
417  Particles.map();
418  Particles.ghost_get<0>();
419  auto Init = getV<0>(Particles);
420  auto Sol = getV<1>(Particles);
421  auto OdeSol = getV<2>(Particles);
422 
423  state_type x0;
424  x0=Init;
425  // The rhs of x' = f(x)
426  //size_t steps=boost::numeric::odeint::integrate(sigmoid,x0,0.0,tf,dt);
427  //typedef boost::numeric::odeint::controlled_runge_kutta< boost::numeric::odeint::runge_kutta_cash_karp54< state_type > > stepper_type;
428  //integrate_adaptive( stepper_type() , sigmoid , x0 , t , tf , dt);
429  size_t steps=boost::numeric::odeint::integrate_const( boost::numeric::odeint::runge_kutta4< state_type >(),sigmoid,x0,0.0,tf,dt);
430 
431  OdeSol=x0;
432  auto it2 = Particles.getDomainIterator();
433  double worst = 0.0;
434  while (it2.isNext()) {
435  auto p = it2.get();
436  if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst) {
437  worst = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
438  }
439  ++it2;
440  }
441 
442  BOOST_REQUIRE(worst < 1e-8);
443 
444  x0=Init;
445  boost::numeric::odeint::runge_kutta4< state_type > rk4;
446  for( size_t i=0 ; i<int(tf/dt) ; ++i,t+=dt )
447  {
448  rk4.do_step(sigmoid,x0,t,dt);
449  t+=dt;
450  }
451 
452  OdeSol=x0;
453  auto it3 = Particles.getDomainIterator();
454  double worst2 = 0.0;
455  while (it3.isNext()) {
456  auto p = it3.get();
457  if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst2) {
458  worst2 = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
459  }
460  ++it3;
461  }
462  //std::cout<<worst2<<std::endl;
463  BOOST_REQUIRE(worst < 1e-6);
464  BOOST_REQUIRE_EQUAL(worst,worst2);
465 }
466 
467 
468 #ifdef HAVE_EIGEN
469 
470 BOOST_AUTO_TEST_CASE(dcpse_op_react_diff_test) {
471  size_t edgeSemiSize = 5;
472  const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
473  Box<2, double> box({0, 0}, {1.0, 1.0});
474  size_t bc[2] = {PERIODIC, PERIODIC};
475  double spacing[2];
476  spacing[0] = 1.0 / (sz[0]);
477  spacing[1] = 1.0 / (sz[1]);
478  Ghost<2, double> ghost(spacing[0] * 3);
479  double rCut = 3.0 * spacing[0];
480  BOOST_TEST_MESSAGE("Init vector_dist...");
481  double sigma2 = spacing[0] * spacing[1] / (2 * 4);
482 
483  //properties: u, v, du, dv
485  bc,
486  ghost);
487 
488  //Init_DCPSE(domain)
489  BOOST_TEST_MESSAGE("Init domain...");
490 
491  std::srand(std::time(nullptr));
492  const double max_rand = 2147483647;
493 
494  auto it = domain.getGridIterator(sz);
495  size_t pointId = 0;
496  size_t counter = 0;
497  double minNormOne = 999;
498  while (it.isNext())
499  {
500  domain.add();
501  auto key = it.get();
502  mem_id k0 = key.get(0);
503  double x = k0 * it.getSpacing(0);
504  domain.getLastPos()[0] = x;//+ gaussian(rng);
505  mem_id k1 = key.get(1);
506  double y = k1 * it.getSpacing(1);
507  domain.getLastPos()[1] = y;//+gaussian(rng);
508  // Here fill the function value
509  domain.template getLastProp<0>() = 0.0;
510  domain.template getLastProp<1>() = 0.0;
511  domain.template getLastProp<2>() = 0;
512  domain.template getLastProp<3>() = 0;
513  if (x==0.5 && y==0.5){
514  domain.template getLastProp<0>() = 1.0;
515  domain.template getLastProp<1>() = 1.0;
516  }
517  ++counter;
518  ++it;
519  }
520  BOOST_TEST_MESSAGE("Sync domain across processors...");
521 
522  domain.map();
523  domain.ghost_get<0>();
524 
525  //Derivative_x Dx(domain, 2, rCut);
526  //Derivative_y Dy(domain, 2, rCut);
527  //Gradient Grad(domain, 2, rCut);
528  vectorGlobal=(void *) &domain;
529 
530  Laplacian Lap(domain, 2, rCut);
531 
532  auto u = getV<0>(domain);
533  auto v = getV<1>(domain);
534  auto fu = getV<2>(domain);
535  auto fv = getV<3>(domain);
536 
539  x0.data.get<0>()=u;
540  x0.data.get<1>()=v;
541 
542  double dt = 0.001;
543  double t = 0.0;
544  double tf = 10.5;
545  //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;
546  typedef boost::numeric::odeint::runge_kutta4< state_type_2d_ofp,double,state_type_2d_ofp,double,boost::numeric::odeint::vector_space_algebra_ofp> stepper_type;
547 
548  integrate_adaptive( stepper_type() , System , x0 , t , tf , dt);
549  fu=x0.data.get<0>();
550  fv=x0.data.get<1>();
551 
552  domain.template ghost_get<2,3>();
553  u = Lap(fu);
554  v = Lap(fv);
555 
556  auto it2 = domain.getDomainIterator();
557 
558  if (create_vcluster().rank() == 0)
559  {++it2;}
560 
561  while (it2.isNext())
562  {
563  auto p = it2.get();
564 
565  BOOST_REQUIRE_CLOSE(domain.template getProp<0>(p),-1.0,1);
566 
567  ++it2;
568  }
569 }
570 #endif
571 
572 BOOST_AUTO_TEST_SUITE_END()
System of equations.
Definition: System.hpp:23
A 3d Odeint and Openfpm compatible structure.
Definition: Ghost.hpp:39
A 2d Odeint and Openfpm compatible structure.
KeyT const ValueT ValueT OffsetIteratorT OffsetIteratorT int
[in] The number of segments that comprise the sorting data
Laplacian second order on h (spacing)
Definition: Laplacian.hpp:22
void ghost_get(size_t opt=WITH_POSITION)
It synchronize the properties and position of the ghost particles.
This class represent an N-dimensional box.
Definition: Box.hpp:60
Main class that encapsulate a vector properties operand to be used for expressions construction Tempo...
Distributed vector.