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"
19#include "DCPSE/DCPSE_op/DCPSE_op.hpp"
21#include "OdeIntegrators/vector_algebra_ofp.hpp"
24const double a = 2.8e-4;
32template<
typename op_type>
51 dxdt.data.get<0>() =
Lap(u) + (1.0);
52 dxdt.data.get<1>() =
Lap(v) + (2.0);
55 if (create_vcluster().rank() == 0)
57 dxdt.data.get<0>().getVector().get<0>(0) = 0.0;
58 dxdt.data.get<1>().getVector().get<0>(0) = 0.0;
64 for (
int i = 0 ; i < dxdt.data.get<0>().getVector().size(); i++)
66 if (u_max < dxdt.data.get<0>().getVector().get<0>(i))
67 {u_max = dxdt.data.get<0>().getVector().get<0>(i);}
69 if (v_max < dxdt.data.get<1>().getVector().get<0>(i))
70 {v_max = dxdt.data.get<1>().getVector().get<0>(i);}
77 dxdt.data.get<0>() = x.data.get<0>();
78 dxdt.data.get<1>() = 2.0*x.data.get<1>();
82 dxdt.data.get<0>() = x.data.get<0>();
83 dxdt.data.get<1>() = 2.0*x.data.get<1>();
100BOOST_AUTO_TEST_SUITE(odeInt_BASE_tests)
102BOOST_AUTO_TEST_CASE(odeint_base_test1)
104 size_t edgeSemiSize = 512;
105 const size_t sz[2] = {edgeSemiSize,edgeSemiSize };
107 size_t bc[2] = { NON_PERIODIC, NON_PERIODIC };
109 spacing[0] = 1.0 / (sz[0] - 1);
110 spacing[1] = 1.0 / (sz[1] - 1);
111 double rCut = 3.9 * spacing[0];
113 BOOST_TEST_MESSAGE(
"Init vector_dist...");
118 const double dt=0.01;
119 auto it = Particles.getGridIterator(sz);
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);
136 Particles.ghost_get<0>();
137 auto Init = getV<0>(Particles);
138 auto Sol = getV<1>(Particles);
139 auto OdeSol = getV<2>(Particles);
151 size_t steps=boost::numeric::odeint::integrate_const( boost::numeric::odeint::runge_kutta4< state_type >(),Exponential,x0,t,tf,dt);
153 std::cout<<
"Time taken by CPU:"<<tt.
getwct()<<std::endl;
155 auto it2 = Particles.getDomainIterator();
157 while (it2.isNext()) {
159 if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst) {
160 worst = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
166 BOOST_REQUIRE(worst < 1e-6);
170 boost::numeric::odeint::runge_kutta4< state_type > rk4;
171 for(
int i=0;i<floor((tf-t0)/dt+0.5);i++)
173 rk4.do_step(Exponential,x0,t,dt);
179 auto it3 = Particles.getDomainIterator();
181 while (it3.isNext()) {
183 if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst2) {
184 worst2 = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
189 BOOST_REQUIRE(worst < 1e-6);
190 BOOST_REQUIRE_EQUAL(worst,worst2);
193BOOST_AUTO_TEST_CASE(odeint_base_test_STRUCT_ofp)
195 size_t edgeSemiSize = 40;
196 const size_t sz[2] = {edgeSemiSize,edgeSemiSize };
198 size_t bc[2] = { NON_PERIODIC, NON_PERIODIC };
200 spacing[0] = 1.0 / (sz[0] - 1);
201 spacing[1] = 1.0 / (sz[1] - 1);
202 double rCut = 3.9 * spacing[0];
204 BOOST_TEST_MESSAGE(
"Init vector_dist...");
206 vector_dist<2, double, aggregate<double,double,double,double,double,double>> Particles(0, box, bc, ghost);
208 auto it = Particles.getGridIterator(sz);
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);
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);
235 x0.data.get<0>()=Init1;
236 x0.data.get<1>()=Init2;
237 x0.data.get<2>()=Init1;
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);
253 OdeSol1=x0.data.get<0>();
254 OdeSol2=x0.data.get<1>();
255 auto it2 = Particles.getDomainIterator();
258 while (it2.isNext()) {
260 if (fabs(Particles.getProp<1>(p) - Particles.getProp<4>(p)) > worst) {
261 worst = fabs(Particles.getProp<1>(p) - Particles.getProp<4>(p));
263 if (fabs(Particles.getProp<3>(p) - Particles.getProp<5>(p)) > worst2) {
264 worst2 = fabs(Particles.getProp<3>(p) - Particles.getProp<5>(p));
269 BOOST_REQUIRE(worst < 1e-6);
270 BOOST_REQUIRE(worst2 < 1e-6);
273 x0.data.get<0>()=Init1;
274 x0.data.get<1>()=Init2;
275 x0.data.get<2>()=Init1;
278 boost::numeric::odeint::runge_kutta4< state_type_3d_ofp, double,state_type_3d_ofp,double,boost::numeric::odeint::vector_space_algebra_ofp> rk4;
281 rk4.do_step(Exponential_struct_ofp2,x0,t,dt);
285 OdeSol1=x0.data.get<0>();
286 OdeSol2=x0.data.get<1>();
287 auto it3 = Particles.getDomainIterator();
290 while (it3.isNext()) {
292 if (fabs(Particles.getProp<1>(p) - Particles.getProp<4>(p)) > worst3) {
293 worst3 = fabs(Particles.getProp<1>(p) - Particles.getProp<4>(p));
295 if (fabs(Particles.getProp<3>(p) - Particles.getProp<5>(p)) > worst4) {
296 worst4 = fabs(Particles.getProp<3>(p) - Particles.getProp<5>(p));
301 BOOST_REQUIRE(worst3 < 1e-6);
302 BOOST_REQUIRE(worst4 < 5e-5);
306BOOST_AUTO_TEST_CASE(odeint_base_test2)
308 size_t edgeSemiSize = 40;
309 const size_t sz[2] = {edgeSemiSize,edgeSemiSize };
311 size_t bc[2] = { NON_PERIODIC, NON_PERIODIC };
313 spacing[0] = 1.0 / (sz[0] - 1);
314 spacing[1] = 1.0 / (sz[1] - 1);
315 double rCut = 3.9 * spacing[0];
317 BOOST_TEST_MESSAGE(
"Init vector_dist...");
324 auto it = Particles.getGridIterator(sz);
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);
340 Particles.ghost_get<0>();
341 auto Init = getV<0>(Particles);
342 auto Sol = getV<1>(Particles);
343 auto OdeSol = getV<2>(Particles);
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);
359 auto it2 = Particles.getDomainIterator();
361 while (it2.isNext()) {
363 if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst) {
364 worst = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
369 BOOST_REQUIRE(worst < 1e-6);
372 boost::numeric::odeint::runge_kutta4< state_type > rk4;
373 for(
size_t i=0 ; i<
int(tf/dt) ; ++i,t+=dt )
375 rk4.do_step(Exponential,x0,t,dt);
380 auto it3 = Particles.getDomainIterator();
382 while (it3.isNext()) {
384 if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst2) {
385 worst2 = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
390 BOOST_REQUIRE(worst < 1e-6);
391 BOOST_REQUIRE(worst2 < 1e-6);
394BOOST_AUTO_TEST_CASE(odeint_base_test3)
396 size_t edgeSemiSize = 40;
397 const size_t sz[2] = {edgeSemiSize,edgeSemiSize };
399 size_t bc[2] = { NON_PERIODIC, NON_PERIODIC };
401 spacing[0] = 1.0 / (sz[0] - 1);
402 spacing[1] = 1.0 / (sz[1] - 1);
403 double rCut = 3.9 * spacing[0];
405 BOOST_TEST_MESSAGE(
"Init vector_dist...");
412 auto it = Particles.getGridIterator(sz);
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));
424 Particles.getLastProp<1>() = 1.0/(1.0+exp(-tf));
428 Particles.ghost_get<0>();
429 auto Init = getV<0>(Particles);
430 auto Sol = getV<1>(Particles);
431 auto OdeSol = getV<2>(Particles);
439 size_t steps=boost::numeric::odeint::integrate_const( boost::numeric::odeint::runge_kutta4< state_type >(),sigmoid,x0,0.0,tf,dt);
442 auto it2 = Particles.getDomainIterator();
444 while (it2.isNext()) {
446 if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst) {
447 worst = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
452 BOOST_REQUIRE(worst < 1e-8);
455 boost::numeric::odeint::runge_kutta4< state_type > rk4;
456 for(
size_t i=0 ; i<
int(tf/dt) ; ++i,t+=dt )
458 rk4.do_step(sigmoid,x0,t,dt);
463 auto it3 = Particles.getDomainIterator();
465 while (it3.isNext()) {
467 if (fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p)) > worst2) {
468 worst2 = fabs(Particles.getProp<1>(p) - Particles.getProp<2>(p));
473 BOOST_REQUIRE(worst < 1e-6);
474 BOOST_REQUIRE_EQUAL(worst,worst2);
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};
482 size_t bc[2] = {PERIODIC, PERIODIC};
484 spacing[0] = 1.0 / (sz[0]);
485 spacing[1] = 1.0 / (sz[1]);
487 double rCut = 3.0 * spacing[0];
488 BOOST_TEST_MESSAGE(
"Init vector_dist...");
489 double sigma2 = spacing[0] * spacing[1] / (2 * 4);
497 BOOST_TEST_MESSAGE(
"Init domain...");
499 std::srand(std::time(
nullptr));
500 const double max_rand = 2147483647;
502 auto it = domain.getGridIterator(sz);
505 double minNormOne = 999;
510 mem_id k0 = key.get(0);
511 double x = k0 * it.getSpacing(0);
512 domain.getLastPos()[0] = x;
513 mem_id k1 = key.get(1);
514 double y = k1 * it.getSpacing(1);
515 domain.getLastPos()[1] = y;
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;
528 BOOST_TEST_MESSAGE(
"Sync domain across processors...");
531 domain.ghost_get<0>();
536 vectorGlobal=(
void *) &domain;
538 Laplacian
Lap(domain, 2, rCut);
540 auto u = getV<0>(domain);
541 auto v = getV<1>(domain);
542 auto fu = getV<2>(domain);
543 auto fv = getV<3>(domain);
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;
556 integrate_adaptive( stepper_type() ,
System , x0 , t , tf , dt);
560 domain.template ghost_get<2,3>();
564 auto it2 = domain.getDomainIterator();
566 if (create_vcluster().rank() == 0)
573 BOOST_REQUIRE_CLOSE(domain.template getProp<0>(p),-1.0,1);
579BOOST_AUTO_TEST_SUITE_END()
This class represent an N-dimensional box.
Laplacian second order on h (spacing)
Class for cpu time benchmarking.
void stop()
Stop the timer.
void start()
Start the timer.
double getwct()
Return the elapsed real time.
Main class that encapsulate a vector properties operand to be used for expressions construction.
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.