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