OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
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
18const double a = 2.8e-4;
19const double b = 5e-3;
20const double tau = .1;
21const double k = .005;
22const int dim = 2;
23
24void *gridGlobal;
26
27// State types for systems with different number of ODEs
31
32template<typename DX, typename DY>
33struct 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
86void 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
98void 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
112BOOST_AUTO_TEST_SUITE(odeInt_grid_tests)
113
114BOOST_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
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,
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
197BOOST_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
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
306BOOST_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
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
495BOOST_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);
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
576BOOST_AUTO_TEST_SUITE_END()
This class represent an N-dimensional box.
Definition Box.hpp:61
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