6#define BOOST_TEST_DYN_LINK
9#include <boost/test/unit_test.hpp>
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"
18const double a = 2.8e-4;
32template<
typename DX,
typename DY>
39 Fitz(DX & m_ddx, DY & m_ddy) :
45 const double t)
const {
48 auto u{FD::getV<4>(temp)};
49 auto v{FD::getV<5>(temp)};
54 dxdt.data.get<0>() = ddx(u) + ddy(u) + (1.0);
55 dxdt.data.get<1>() = ddx(v) + ddy(v) + (2.0);
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;
75 if (u_max < dxdt.data.get<0>().value(key))
76 u_max = dxdt.data.get<0>().value(key);
78 if (v_max < dxdt.data.get<1>().value(key))
79 v_max = dxdt.data.get<1>().value(key);
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>();
112BOOST_AUTO_TEST_SUITE(odeInt_grid_tests)
114BOOST_AUTO_TEST_CASE(odeint_grid_test_exponential) {
116 size_t edgeSemiSize{40};
117 const size_t sz[dim] = {edgeSemiSize,edgeSemiSize};
121 spacing[0] = 1.0 / (sz[0] - 1);
122 spacing[1] = 1.0 / (sz[1] - 1);
124 BOOST_TEST_MESSAGE(
"Test: exponential");
125 BOOST_TEST_MESSAGE(
"Init grid_dist_id ...");
129 auto it{
grid.getDomainIterator()};
130 while (it.isNext()) {
132 grid.template get<0>(key) = std::exp(0);
133 grid.template get<1>(key) = std::exp(0.4);
138 auto Init{FD::getV<0>(
grid)};
139 auto Sol{FD::getV<1>(
grid)};
140 auto OdeSol{FD::getV<2>(
grid)};
143 x0.data.get<0>() = Init;
147 const double dt{0.1};
152 size_t steps{boost::numeric::odeint::integrate_const(rk4,Exponential,x0,0.0,tf,dt)};
154 OdeSol = x0.data.get<0>();
157 auto it2{
grid.getDomainIterator()};
159 while (it2.isNext()) {
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));
167 std::cout << worst << std::endl;
168 BOOST_REQUIRE(worst < 1e-6);
171 x0.data.get<0>() = Init;
173 rk4.do_step(Exponential,x0,t,dt);
174 OdeSol = x0.data.get<0>();
178 OdeSol = x0.data.get<0>();
181 auto it3{
grid.getDomainIterator()};
183 while (it3.isNext()) {
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));
190 std::cout << worst2 << std::endl;
191 BOOST_REQUIRE(worst2 < 1e-6);
192 BOOST_REQUIRE_EQUAL(worst,worst2);
197BOOST_AUTO_TEST_CASE(odeint_grid_test_STRUCT_exponential) {
199 size_t edgeSemiSize{40};
200 const size_t sz[dim] = {edgeSemiSize,edgeSemiSize};
204 spacing[0] = 1.0 / (sz[0] - 1);
205 spacing[1] = 1.0 / (sz[1] - 1);
207 BOOST_TEST_MESSAGE(
"Test: exponential");
208 BOOST_TEST_MESSAGE(
"Init grid_dist_id ...");
210 grid_dist_id<dim,double,aggregate<double,double,double,double,double,double>>
grid{sz,box,ghost,bc};
212 auto it{
grid.getDomainIterator()};
213 while (it.isNext()) {
216 grid.get<0>(key) = std::exp(0);
217 grid.template get<0>(key) = std::exp(0.0);
218 grid.template get<1>(key) = std::exp(0.4);
219 grid.template get<2>(key) = std::exp(0.0);
220 grid.template get<3>(key) = std::exp(0.8);
225 auto Init1{FD::getV<0>(
grid)};
226 auto Sol1{FD::getV<1>(
grid)};
228 auto Init2{FD::getV<2>(
grid)};
229 auto Sol2{FD::getV<3>(
grid)};
231 auto OdeSol1{FD::getV<4>(
grid)};
232 auto OdeSol2{FD::getV<5>(
grid)};
235 x0.data.get<0>() = Init1;
236 x0.data.get<1>() = Init2;
237 x0.data.get<2>() = Init1;
241 const double dt{0.1};
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);
250 OdeSol1 = x0.data.get<0>();
251 OdeSol2 = x0.data.get<1>();
254 auto it2{
grid.getDomainIterator()};
257 while (it2.isNext()) {
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));
266 std::cout << worst <<
" " << worst2 << std::endl;
267 BOOST_REQUIRE(worst < 1e-6);
268 BOOST_REQUIRE(worst2 < 1e-6);
272 x0.data.get<0>() = Init1;
273 x0.data.get<1>() = Init2;
274 x0.data.get<2>() = Init1;
276 boost::numeric::odeint::runge_kutta4<state_type_3ode,double,state_type_3ode,double,boost::numeric::odeint::vector_space_algebra_ofp> rk4;
278 rk4.do_step(Exponential_struct_ofp2,x0,t,dt);
282 OdeSol1 = x0.data.get<0>();
283 OdeSol2 = x0.data.get<1>();
286 auto it3{
grid.getDomainIterator()};
289 while (it3.isNext()) {
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));
298 std::cout << worst3 <<
" " << worst4 << std::endl;
299 BOOST_REQUIRE(worst3 < 1e-6);
300 BOOST_REQUIRE(worst4 < 5e-5);
306BOOST_AUTO_TEST_CASE(odeint_grid_test2_exponential) {
308 size_t edgeSemiSize{40};
309 const size_t sz[dim] = {edgeSemiSize,edgeSemiSize};
313 spacing[0] = 1.0 / (sz[0] - 1);
314 spacing[1] = 1.0 / (sz[1] - 1);
316 BOOST_TEST_MESSAGE(
"Test: exponential");
317 BOOST_TEST_MESSAGE(
"Init grid_dist_id ...");
323 const double dt{0.1};
325 auto it{
grid.getDomainIterator()};
326 while (it.isNext()) {
328 grid.template get<0>(key) = std::exp(t);
329 grid.template get<1>(key) = std::exp(tf);
334 auto Init{FD::getV<0>(
grid)};
335 auto Sol{FD::getV<1>(
grid)};
336 auto OdeSol{FD::getV<2>(
grid)};
339 x0.data.get<0>() = Init;
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;
344 integrate_adaptive(stepper_type(),Exponential,x0,t,tf,dt);
345 OdeSol = x0.data.get<0>();
348 auto it2{
grid.getDomainIterator()};
350 while (it2.isNext()) {
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));
358 std::cout << worst << std::endl;
359 BOOST_REQUIRE(worst < 1e-6);
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);
368 OdeSol = x0.data.get<0>();
371 auto it3{
grid.getDomainIterator()};
373 while (it3.isNext()) {
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));
380 std::cout << worst2 << std::endl;
381 BOOST_REQUIRE(worst2 < 1e-6);
495BOOST_AUTO_TEST_CASE(dcpse_op_react_diff_test) {
497 size_t edgeSemiSize{5};
498 const size_t sz[dim] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
502 spacing[0] = 1.0 / (sz[0]);
503 spacing[1] = 1.0 / (sz[1]);
506 BOOST_TEST_MESSAGE(
"Test: reaction diffusion");
507 BOOST_TEST_MESSAGE(
"Init grid_dist_id ...");
508 double sigma2 = spacing[0] * spacing[1] / (2 * 4);
511 grid_dist_id<dim,double,aggregate<double,double,double,double,double,double>> domain{sz,box,ghost,bc};
513 auto it{domain.getDomainIterator()};
514 while (it.isNext()) {
516 domain.get<0>(key) = 0.0;
517 domain.get<1>(key) = 0.0;
518 domain.get<2>(key) = 0.0;
519 domain.get<3>(key) = 0.0;
521 auto gkey = it.getGKey(key);
523 if (gkey.get(0)==sz[0] / 2 && gkey.get(1) == sz[1]/2)
525 domain.get<0>(key) = 1.0;
526 domain.get<1>(key) = 1.0;
530 domain.ghost_get<0>();
535 gridGlobal=(
void *) & domain;
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)};
542 Fitz<
decltype(ddx),
decltype(ddy)> system(ddx,ddy);
545 x0.data.get<0>() = u;
546 x0.data.get<1>() = v;
553 typedef boost::numeric::odeint::runge_kutta4<state_type_2ode,double,state_type_2ode,double,boost::numeric::odeint::vector_space_algebra_ofp> stepper_type;
555 integrate_adaptive(stepper_type(),system,x0,t,tf,dt);
556 fu = x0.data.get<0>();
557 fv = x0.data.get<1>();
559 domain.ghost_get<2,3>();
560 u = ddx(fu) + ddy(fu);
561 v = ddx(fv) + ddy(fv);
563 auto it2{domain.getDomainIterator()};
565 if (create_vcluster().rank() == 0)
568 while (it2.isNext()) {
570 BOOST_REQUIRE_CLOSE(domain.get<0>(p),-1.0,1);
576BOOST_AUTO_TEST_SUITE_END()
This class represent an N-dimensional box.
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