7 #include "util/common.hpp" 9 #define BOOST_TEST_DYN_LINK 11 #include "util/util_debug.hpp" 12 #include <boost/test/unit_test.hpp> 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" 22 const double a = 2.8e-4;
23 const double b = 5e-3;
24 const double tau = .1;
25 const double k = .005;
30 template<
typename op_type>
49 dxdt.data.get<0>() =
Lap(u) + (1.0);
50 dxdt.data.get<1>() =
Lap(v) + (2.0);
53 if (create_vcluster().rank() == 0)
55 dxdt.data.get<0>().getVector().get<0>(0) = 0.0;
56 dxdt.data.get<1>().getVector().get<0>(0) = 0.0;
62 for (
int i = 0 ; i < dxdt.data.get<0>().getVector().size(); i++)
64 if (u_max < dxdt.data.get<0>().getVector().get<0>(i))
65 {u_max = dxdt.data.get<0>().getVector().get<0>(i);}
67 if (v_max < dxdt.data.get<1>().getVector().get<0>(i))
68 {v_max = dxdt.data.get<1>().getVector().get<0>(i);}
75 dxdt.data.get<0>() = x.data.get<0>();
76 dxdt.data.get<1>() = 2.0*x.data.get<1>();
80 dxdt.data.get<0>() = x.data.get<0>();
81 dxdt.data.get<1>() = 2.0*x.data.get<1>();
94 BOOST_AUTO_TEST_SUITE(odeInt_BASE_tests)
96 BOOST_AUTO_TEST_CASE(odeint_base_test1)
98 size_t edgeSemiSize = 40;
99 const size_t sz[2] = {edgeSemiSize,edgeSemiSize };
101 size_t bc[2] = { NON_PERIODIC, NON_PERIODIC };
103 spacing[0] = 1.0 / (sz[0] - 1);
104 spacing[1] = 1.0 / (sz[1] - 1);
105 double rCut = 3.9 * spacing[0];
107 BOOST_TEST_MESSAGE(
"Init vector_dist...");
111 auto it = Particles.getGridIterator(sz);
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);
128 Particles.ghost_get<0>();
129 auto Init = getV<0>(Particles);
130 auto Sol = getV<1>(Particles);
131 auto OdeSol = getV<2>(Particles);
143 size_t steps=boost::numeric::odeint::integrate_const( boost::numeric::odeint::runge_kutta4< state_type >(),Exponential,x0,0.0,tf,dt);
146 auto it2 = Particles.getDomainIterator();
148 while (it2.isNext()) {
150 if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst) {
151 worst = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
157 BOOST_REQUIRE(worst < 1e-6);
160 boost::numeric::odeint::runge_kutta4< state_type > rk4;
163 rk4.do_step(Exponential,x0,t,dt);
169 auto it3 = Particles.getDomainIterator();
171 while (it3.isNext()) {
173 if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst2) {
174 worst2 = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
179 BOOST_REQUIRE(worst < 1e-6);
180 BOOST_REQUIRE_EQUAL(worst,worst2);
183 BOOST_AUTO_TEST_CASE(odeint_base_test_STRUCT_ofp)
185 size_t edgeSemiSize = 40;
186 const size_t sz[2] = {edgeSemiSize,edgeSemiSize };
188 size_t bc[2] = { NON_PERIODIC, NON_PERIODIC };
190 spacing[0] = 1.0 / (sz[0] - 1);
191 spacing[1] = 1.0 / (sz[1] - 1);
192 double rCut = 3.9 * spacing[0];
194 BOOST_TEST_MESSAGE(
"Init vector_dist...");
196 vector_dist<2, double, aggregate<double,double,double,double,double,double>> Particles(0, box, bc, ghost);
198 auto it = Particles.getGridIterator(sz);
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);
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);
225 x0.data.get<0>()=Init1;
226 x0.data.get<1>()=Init2;
227 x0.data.get<2>()=Init1;
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);
243 OdeSol1=x0.data.get<0>();
244 OdeSol2=x0.data.get<1>();
245 auto it2 = Particles.getDomainIterator();
248 while (it2.isNext()) {
250 if (fabs(Particles.getProp<1>(p) - Particles.getProp<4>(p)) > worst) {
251 worst = fabs(Particles.getProp<1>(p) - Particles.getProp<4>(p));
253 if (fabs(Particles.getProp<3>(p) - Particles.getProp<5>(p)) > worst2) {
254 worst2 = fabs(Particles.getProp<3>(p) - Particles.getProp<5>(p));
259 BOOST_REQUIRE(worst < 1e-6);
260 BOOST_REQUIRE(worst2 < 1e-6);
263 x0.data.get<0>()=Init1;
264 x0.data.get<1>()=Init2;
265 x0.data.get<2>()=Init1;
268 boost::numeric::odeint::runge_kutta4< state_type_3d_ofp, double,state_type_3d_ofp,double,boost::numeric::odeint::vector_space_algebra_ofp> rk4;
271 rk4.do_step(Exponential_struct_ofp2,x0,t,dt);
275 OdeSol1=x0.data.get<0>();
276 OdeSol2=x0.data.get<1>();
277 auto it3 = Particles.getDomainIterator();
280 while (it3.isNext()) {
282 if (fabs(Particles.getProp<1>(p) - Particles.getProp<4>(p)) > worst3) {
283 worst3 = fabs(Particles.getProp<1>(p) - Particles.getProp<4>(p));
285 if (fabs(Particles.getProp<3>(p) - Particles.getProp<5>(p)) > worst4) {
286 worst4 = fabs(Particles.getProp<3>(p) - Particles.getProp<5>(p));
291 BOOST_REQUIRE(worst3 < 1e-6);
292 BOOST_REQUIRE(worst4 < 5e-5);
296 BOOST_AUTO_TEST_CASE(odeint_base_test2)
298 size_t edgeSemiSize = 40;
299 const size_t sz[2] = {edgeSemiSize,edgeSemiSize };
301 size_t bc[2] = { NON_PERIODIC, NON_PERIODIC };
303 spacing[0] = 1.0 / (sz[0] - 1);
304 spacing[1] = 1.0 / (sz[1] - 1);
305 double rCut = 3.9 * spacing[0];
307 BOOST_TEST_MESSAGE(
"Init vector_dist...");
314 auto it = Particles.getGridIterator(sz);
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);
330 Particles.ghost_get<0>();
331 auto Init = getV<0>(Particles);
332 auto Sol = getV<1>(Particles);
333 auto OdeSol = getV<2>(Particles);
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);
349 auto it2 = Particles.getDomainIterator();
351 while (it2.isNext()) {
353 if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst) {
354 worst = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
359 BOOST_REQUIRE(worst < 1e-6);
362 boost::numeric::odeint::runge_kutta4< state_type > rk4;
363 for(
size_t i=0 ; i<
int(tf/dt) ; ++i,t+=dt )
365 rk4.do_step(Exponential,x0,t,dt);
370 auto it3 = Particles.getDomainIterator();
372 while (it3.isNext()) {
374 if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst2) {
375 worst2 = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
380 BOOST_REQUIRE(worst < 1e-6);
381 BOOST_REQUIRE(worst2 < 1e-6);
384 BOOST_AUTO_TEST_CASE(odeint_base_test3)
386 size_t edgeSemiSize = 40;
387 const size_t sz[2] = {edgeSemiSize,edgeSemiSize };
389 size_t bc[2] = { NON_PERIODIC, NON_PERIODIC };
391 spacing[0] = 1.0 / (sz[0] - 1);
392 spacing[1] = 1.0 / (sz[1] - 1);
393 double rCut = 3.9 * spacing[0];
395 BOOST_TEST_MESSAGE(
"Init vector_dist...");
402 auto it = Particles.getGridIterator(sz);
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));
414 Particles.getLastProp<1>() = 1.0/(1.0+exp(-tf));
418 Particles.ghost_get<0>();
419 auto Init = getV<0>(Particles);
420 auto Sol = getV<1>(Particles);
421 auto OdeSol = getV<2>(Particles);
429 size_t steps=boost::numeric::odeint::integrate_const( boost::numeric::odeint::runge_kutta4< state_type >(),sigmoid,x0,0.0,tf,dt);
432 auto it2 = Particles.getDomainIterator();
434 while (it2.isNext()) {
436 if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst) {
437 worst = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
442 BOOST_REQUIRE(worst < 1e-8);
445 boost::numeric::odeint::runge_kutta4< state_type > rk4;
446 for(
size_t i=0 ; i<
int(tf/dt) ; ++i,t+=dt )
448 rk4.do_step(sigmoid,x0,t,dt);
453 auto it3 = Particles.getDomainIterator();
455 while (it3.isNext()) {
457 if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst2) {
458 worst2 = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
463 BOOST_REQUIRE(worst < 1e-6);
464 BOOST_REQUIRE_EQUAL(worst,worst2);
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};
474 size_t bc[2] = {PERIODIC, PERIODIC};
476 spacing[0] = 1.0 / (sz[0]);
477 spacing[1] = 1.0 / (sz[1]);
479 double rCut = 3.0 * spacing[0];
480 BOOST_TEST_MESSAGE(
"Init vector_dist...");
481 double sigma2 = spacing[0] * spacing[1] / (2 * 4);
489 BOOST_TEST_MESSAGE(
"Init domain...");
491 std::srand(std::time(
nullptr));
492 const double max_rand = 2147483647;
494 auto it = domain.getGridIterator(sz);
497 double minNormOne = 999;
502 mem_id k0 = key.get(0);
503 double x = k0 * it.getSpacing(0);
504 domain.getLastPos()[0] = x;
505 mem_id k1 = key.get(1);
506 double y = k1 * it.getSpacing(1);
507 domain.getLastPos()[1] = y;
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;
520 BOOST_TEST_MESSAGE(
"Sync domain across processors...");
523 domain.ghost_get<0>();
528 vectorGlobal=(
void *) &domain;
530 Laplacian
Lap(domain, 2, rCut);
532 auto u = getV<0>(domain);
533 auto v = getV<1>(domain);
534 auto fu = getV<2>(domain);
535 auto fv = getV<3>(domain);
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;
548 integrate_adaptive( stepper_type() ,
System , x0 , t , tf , dt);
552 domain.template ghost_get<2,3>();
556 auto it2 = domain.getDomainIterator();
558 if (create_vcluster().rank() == 0)
565 BOOST_REQUIRE_CLOSE(domain.template getProp<0>(p),-1.0,1);
572 BOOST_AUTO_TEST_SUITE_END()
A 3d Odeint and Openfpm compatible structure.
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)
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.
Main class that encapsulate a vector properties operand to be used for expressions construction Tempo...