| 1 | /* | 
| 2 |  * sgrid_dist_id_unit_tests.cpp | 
| 3 |  * | 
| 4 |  *  Created on: Nov 18, 2017 | 
| 5 |  *      Author: i-bird | 
| 6 |  */ | 
| 7 |  | 
| 8 |  | 
| 9 | #define BOOST_TEST_DYN_LINK | 
| 10 | #include <boost/test/unit_test.hpp> | 
| 11 | #include "Grid/grid_dist_id.hpp" | 
| 12 | #include "Point_test.hpp" | 
| 13 |  | 
| 14 |  | 
| 15 | const int x = 0; | 
| 16 | const int y = 1; | 
| 17 | const int z = 2; | 
| 18 |  | 
| 19 | BOOST_AUTO_TEST_SUITE( sgrid_dist_id_test ) | 
| 20 |  | 
| 21 |  | 
| 22 |  | 
| 23 | BOOST_AUTO_TEST_CASE (sgrid_dist_id_soa ) | 
| 24 | { | 
| 25 | 	periodicity<3> bc = {PERIODIC, PERIODIC, PERIODIC}; | 
| 26 |  | 
| 27 | 	// Domain | 
| 28 | 	Box<3,double> domain({-0.3,-0.3,-0.3},{1.0,1.0,1.0}); | 
| 29 |  | 
| 30 | 	// grid size | 
| 31 | 	size_t sz[3]; | 
| 32 | 	sz[0] = 1024; | 
| 33 | 	sz[1] = 1024; | 
| 34 | 	sz[2] = 1024; | 
| 35 |  | 
| 36 | 	// Ghost | 
| 37 | 	Ghost<3,double> g(0.01); | 
| 38 |  | 
| 39 | 	sgrid_dist_soa<3,double,Point_test<float>> sg(sz,domain,g,bc); | 
| 40 |  | 
| 41 | 	// create a grid iterator over a bilion point | 
| 42 |  | 
| 43 | 	auto it = sg.getGridIterator(); | 
| 44 |  | 
| 45 | 	while(it.isNext()) | 
| 46 | 	{ | 
| 47 | 		auto gkey = it.get(); | 
| 48 | 		auto key = it.get_dist(); | 
| 49 |  | 
| 50 | 		size_t sx = gkey.get(0) - 512; | 
| 51 | 		size_t sy = gkey.get(1) - 512; | 
| 52 | 		size_t sz = gkey.get(2) - 512; | 
| 53 |  | 
| 54 | 		if (sx*sx + sy*sy + sz*sz < 128*128) | 
| 55 | 		{ | 
| 56 | 			sg.template insert<0>(key) = 1.0; | 
| 57 | 		} | 
| 58 |  | 
| 59 | 		++it; | 
| 60 | 	} | 
| 61 |  | 
| 62 | 	bool match = true; | 
| 63 | 	auto it2 = sg.getGridIterator(); | 
| 64 |  | 
| 65 | 	while(it2.isNext()) | 
| 66 | 	{ | 
| 67 | 		auto gkey = it2.get(); | 
| 68 | 		auto key = it2.get_dist(); | 
| 69 |  | 
| 70 | 		size_t sx = gkey.get(0) - 512; | 
| 71 | 		size_t sy = gkey.get(1) - 512; | 
| 72 | 		size_t sz = gkey.get(2) - 512; | 
| 73 |  | 
| 74 | 		if (sx*sx + sy*sy + sz*sz < 128*128) | 
| 75 | 		{ | 
| 76 | 			match &= (sg.template get<0>(key) == 1.0); | 
| 77 | 		} | 
| 78 |  | 
| 79 | 		++it2; | 
| 80 | 	} | 
| 81 |  | 
| 82 | 	auto & gr = sg.getGridInfo(); | 
| 83 |  | 
| 84 | 	auto it3 = sg.getDomainIterator(); | 
| 85 |  | 
| 86 | 	while (it3.isNext()) | 
| 87 | 	{ | 
| 88 | 		auto key = it3.get(); | 
| 89 | 		auto gkey = it3.getGKey(key); | 
| 90 |  | 
| 91 | 		sg.template insert<0>(key) = gkey.get(0)*gkey.get(0) + gkey.get(1)*gkey.get(1) + gkey.get(2)*gkey.get(2); | 
| 92 |  | 
| 93 | 		++it3; | 
| 94 | 	} | 
| 95 |  | 
| 96 | 	sg.ghost_get<0>(); | 
| 97 | 	// now we check the stencil | 
| 98 |  | 
| 99 | 	bool good = true; | 
| 100 | 	auto it4 = sg.getDomainIterator(); | 
| 101 |  | 
| 102 | 	while (it4.isNext()) | 
| 103 | 	{ | 
| 104 | 		auto key = it4.get(); | 
| 105 | 		auto gkey = it4.getGKey(key); | 
| 106 |  | 
| 107 | 		size_t sx = gkey.get(0) - 512; | 
| 108 | 		size_t sy = gkey.get(1) - 512; | 
| 109 | 		size_t sz = gkey.get(2) - 512; | 
| 110 |  | 
| 111 | 		double lap; | 
| 112 |  | 
| 113 | 		lap = sg.template get<0>(key.move(x,1)) + sg.template get<0>(key.move(x,-1)) + | 
| 114 | 			  sg.template get<0>(key.move(y,1)) + sg.template get<0>(key.move(y,-1)) + | 
| 115 | 			  sg.template get<0>(key.move(z,1)) + sg.template get<0>(key.move(z,-1)) - | 
| 116 | 			  6.0*sg.template get<0>(key); | 
| 117 |  | 
| 118 | 		good &= (lap == 6.0); | 
| 119 |  | 
| 120 | 		++it4; | 
| 121 | 	} | 
| 122 |  | 
| 123 | 	BOOST_REQUIRE_EQUAL(match,true); | 
| 124 | } | 
| 125 |  | 
| 126 | BOOST_AUTO_TEST_CASE( sgrid_dist_id_basic_test_2D) | 
| 127 | { | 
| 128 | 	periodicity<2> bc = {NON_PERIODIC, NON_PERIODIC}; | 
| 129 |  | 
| 130 | 	// Domain | 
| 131 | 	Box<2,double> domain({-0.3,-0.3},{1.0,1.0}); | 
| 132 |  | 
| 133 | 	// grid size | 
| 134 | 	size_t sz[2]; | 
| 135 | 	sz[0] = 1024; | 
| 136 | 	sz[1] = 1024; | 
| 137 |  | 
| 138 | 	// Ghost | 
| 139 | 	Ghost<2,double> g(0.01); | 
| 140 |  | 
| 141 | 	sgrid_dist_id<2,double,Point_test<double>> sg(sz,domain,g,bc); | 
| 142 |  | 
| 143 | 	// create a grid iterator | 
| 144 |  | 
| 145 | 	auto it = sg.getGridIterator(); | 
| 146 |  | 
| 147 | 	while(it.isNext()) | 
| 148 | 	{ | 
| 149 | 		auto gkey = it.get(); | 
| 150 | 		auto key = it.get_dist(); | 
| 151 |  | 
| 152 |  | 
| 153 | 		long int sx = gkey.get(0) - 512; | 
| 154 | 		long int sy = gkey.get(1) - 512; | 
| 155 |  | 
| 156 | 		if (sx*sx + sy*sy < 128*128) | 
| 157 | 		{ | 
| 158 | 			sg.template insert<0>(key) = 1.0; | 
| 159 | 		} | 
| 160 |  | 
| 161 | 		++it; | 
| 162 | 	} | 
| 163 |  | 
| 164 | 	bool match = true; | 
| 165 | 	auto it2 = sg.getGridIterator(); | 
| 166 |  | 
| 167 | 	while(it2.isNext()) | 
| 168 | 	{ | 
| 169 | 		auto gkey = it2.get(); | 
| 170 | 		auto key = it2.get_dist(); | 
| 171 |  | 
| 172 | 		long int sx = gkey.get(0) - 512; | 
| 173 | 		long int sy = gkey.get(1) - 512; | 
| 174 |  | 
| 175 | 		if (sx*sx + sy*sy < 128*128) | 
| 176 | 		{ | 
| 177 | 			match &= (sg.template get<0>(key) == 1.0); | 
| 178 | 		} | 
| 179 |  | 
| 180 | 		++it2; | 
| 181 | 	} | 
| 182 |  | 
| 183 | 	BOOST_REQUIRE_EQUAL(match,true); | 
| 184 |  | 
| 185 | 	auto & gr = sg.getGridInfo(); | 
| 186 |  | 
| 187 | 	auto it3 = sg.getDomainIterator(); | 
| 188 |  | 
| 189 | 	while (it3.isNext()) | 
| 190 | 	{ | 
| 191 | 		auto key = it3.get(); | 
| 192 | 		auto gkey = it3.getGKey(key); | 
| 193 |  | 
| 194 | 		sg.template insert<0>(key) = gkey.get(0)*gkey.get(0) + gkey.get(1)*gkey.get(1); | 
| 195 |  | 
| 196 | 		++it3; | 
| 197 | 	} | 
| 198 |  | 
| 199 | 	sg.ghost_get<0>(); | 
| 200 |  | 
| 201 | 	// now we check the stencil | 
| 202 |  | 
| 203 | 	bool good = true; | 
| 204 | 	auto it4 = sg.getDomainIterator(); | 
| 205 |  | 
| 206 | 	while (it4.isNext()) | 
| 207 | 	{ | 
| 208 | 		auto key = it4.get(); | 
| 209 | 		auto gkey = it4.getGKey(key); | 
| 210 |  | 
| 211 | 		double lap; | 
| 212 |  | 
| 213 | 		// Here we check that all point of the stencil are inside*/ | 
| 214 |  | 
| 215 | 		long int sx = gkey.get(0) - 512; | 
| 216 | 		long int sy = gkey.get(1) - 512; | 
| 217 |  | 
| 218 | 		if (sx*sx + sy*sy < 126*126) | 
| 219 | 		{ | 
| 220 | 			lap = sg.template get<0>(key.move(x,1)) + sg.template get<0>(key.move(x,-1)) + | 
| 221 | 				  sg.template get<0>(key.move(y,1)) + sg.template get<0>(key.move(y,-1)) - | 
| 222 | 				  4.0*sg.template get<0>(key); | 
| 223 |  | 
| 224 | 			good &= (lap == 4.0); | 
| 225 | 		} | 
| 226 |  | 
| 227 | 		++it4; | 
| 228 | 	} | 
| 229 |  | 
| 230 | 	BOOST_REQUIRE_EQUAL(good,true); | 
| 231 | } | 
| 232 |  | 
| 233 | BOOST_AUTO_TEST_CASE( sgrid_dist_id_basic_test) | 
| 234 | { | 
| 235 | 	periodicity<3> bc = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC}; | 
| 236 |  | 
| 237 | 	// Domain | 
| 238 | 	Box<3,double> domain({-0.3,-0.3,-0.3},{1.0,1.0,1.0}); | 
| 239 |  | 
| 240 | 	// grid size | 
| 241 | 	size_t sz[3]; | 
| 242 | 	sz[0] = 1024; | 
| 243 | 	sz[1] = 1024; | 
| 244 | 	sz[2] = 1024; | 
| 245 |  | 
| 246 | 	// Ghost | 
| 247 | 	Ghost<3,double> g(0.01); | 
| 248 |  | 
| 249 | 	sgrid_dist_id<3,double,Point_test<float>> sg(sz,domain,g,bc); | 
| 250 |  | 
| 251 | 	// create a grid iterator over a bilion point | 
| 252 |  | 
| 253 | 	auto it = sg.getGridIterator(); | 
| 254 |  | 
| 255 | 	while(it.isNext()) | 
| 256 | 	{ | 
| 257 | 		auto gkey = it.get(); | 
| 258 | 		auto key = it.get_dist(); | 
| 259 |  | 
| 260 | 		size_t sx = gkey.get(0) - 512; | 
| 261 | 		size_t sy = gkey.get(1) - 512; | 
| 262 | 		size_t sz = gkey.get(2) - 512; | 
| 263 |  | 
| 264 | 		if (sx*sx + sy*sy + sz*sz < 128*128) | 
| 265 | 		{ | 
| 266 | 			sg.template insert<0>(key) = 1.0; | 
| 267 | 		} | 
| 268 |  | 
| 269 | 		++it; | 
| 270 | 	} | 
| 271 |  | 
| 272 | 	bool match = true; | 
| 273 | 	auto it2 = sg.getGridIterator(); | 
| 274 |  | 
| 275 | 	while(it2.isNext()) | 
| 276 | 	{ | 
| 277 | 		auto gkey = it2.get(); | 
| 278 | 		auto key = it2.get_dist(); | 
| 279 |  | 
| 280 | 		size_t sx = gkey.get(0) - 512; | 
| 281 | 		size_t sy = gkey.get(1) - 512; | 
| 282 | 		size_t sz = gkey.get(2) - 512; | 
| 283 |  | 
| 284 | 		if (sx*sx + sy*sy + sz*sz < 128*128) | 
| 285 | 		{ | 
| 286 | 			match &= (sg.template get<0>(key) == 1.0); | 
| 287 | 		} | 
| 288 |  | 
| 289 | 		++it2; | 
| 290 | 	} | 
| 291 |  | 
| 292 | 	auto & gr = sg.getGridInfo(); | 
| 293 |  | 
| 294 | 	auto it3 = sg.getDomainIterator(); | 
| 295 |  | 
| 296 | 	while (it3.isNext()) | 
| 297 | 	{ | 
| 298 | 		auto key = it3.get(); | 
| 299 | 		auto gkey = it3.getGKey(key); | 
| 300 |  | 
| 301 | 		sg.template insert<0>(key) = gkey.get(0)*gkey.get(0) + gkey.get(1)*gkey.get(1) + gkey.get(2)*gkey.get(2); | 
| 302 |  | 
| 303 | 		++it3; | 
| 304 | 	} | 
| 305 |  | 
| 306 | 	sg.ghost_get<0>(); | 
| 307 | 	// now we check the stencil | 
| 308 |  | 
| 309 | 	bool good = true; | 
| 310 | 	auto it4 = sg.getDomainIterator(); | 
| 311 |  | 
| 312 | 	while (it4.isNext()) | 
| 313 | 	{ | 
| 314 | 		auto key = it4.get(); | 
| 315 | 		auto gkey = it4.getGKey(key); | 
| 316 |  | 
| 317 | 		size_t sx = gkey.get(0) - 512; | 
| 318 | 		size_t sy = gkey.get(1) - 512; | 
| 319 | 		size_t sz = gkey.get(2) - 512; | 
| 320 |  | 
| 321 | 		if (sx*sx + sy*sy + sz*sz < 125*125) | 
| 322 | 		{ | 
| 323 | 			double lap; | 
| 324 |  | 
| 325 | 			lap = sg.template get<0>(key.move(x,1)) + sg.template get<0>(key.move(x,-1)) + | 
| 326 | 				  sg.template get<0>(key.move(y,1)) + sg.template get<0>(key.move(y,-1)) + | 
| 327 | 				  sg.template get<0>(key.move(z,1)) + sg.template get<0>(key.move(z,-1)) - | 
| 328 | 				  6.0*sg.template get<0>(key); | 
| 329 |  | 
| 330 | 			good &= (lap == 6.0); | 
| 331 | 		} | 
| 332 |  | 
| 333 | 		++it4; | 
| 334 | 	} | 
| 335 |  | 
| 336 | 	BOOST_REQUIRE_EQUAL(match,true); | 
| 337 | } | 
| 338 |  | 
| 339 |  | 
| 340 | BOOST_AUTO_TEST_CASE( sparse_grid_fast_stencil_vectorized_simplified_conv2) | 
| 341 | { | 
| 342 | 	constexpr int U = 0; | 
| 343 | 	constexpr int V = 1; | 
| 344 |  | 
| 345 | 	constexpr int U_next = 2; | 
| 346 | 	constexpr int V_next = 3; | 
| 347 |  | 
| 348 | 	constexpr int x = 0; | 
| 349 | 	constexpr int y = 1; | 
| 350 | 	constexpr int z = 2; | 
| 351 |  | 
| 352 |     Box<3,double> domain({0.0,0.0,0.0},{2.5,2.5,2.5}); | 
| 353 |  | 
| 354 |     // grid size | 
| 355 |     size_t sz[3] = {32,32,32}; | 
| 356 |  | 
| 357 |     // Define periodicity of the grid | 
| 358 |     periodicity<3> bc = {PERIODIC,PERIODIC,PERIODIC}; | 
| 359 |  | 
| 360 |     // Ghost in grid unit | 
| 361 |     Ghost<3,long int> g(1); | 
| 362 |  | 
| 363 |     // deltaT | 
| 364 |     double deltaT = 1; | 
| 365 |  | 
| 366 |     // Diffusion constant for specie U | 
| 367 |     double du = 2*1e-5; | 
| 368 |  | 
| 369 |     // Diffusion constant for specie V | 
| 370 |     double dv = 1*1e-5; | 
| 371 |  | 
| 372 |     // Number of timesteps | 
| 373 |     size_t timeSteps = 5000; | 
| 374 |  | 
| 375 |     // K and F (Physical constant in the equation) | 
| 376 |     double K = 0.053; | 
| 377 |     double F = 0.014; | 
| 378 |  | 
| 379 |     sgrid_dist_soa<3, double, aggregate<double,double,double,double>> grid(sz,domain,g,bc); | 
| 380 |  | 
| 381 |     auto it = grid.getGridIterator(); | 
| 382 |  | 
| 383 |     while (it.isNext()) | 
| 384 |     { | 
| 385 |             // Get the local grid key | 
| 386 |             auto key = it.get_dist(); | 
| 387 |  | 
| 388 |             // Old values U and V | 
| 389 |             grid.template insert<U>(key) = 1.0; | 
| 390 |             grid.template insert<V>(key) = 0.0; | 
| 391 |  | 
| 392 |             // Old values U and V | 
| 393 |             grid.template insert<U_next>(key) = 0.0; | 
| 394 |             grid.template insert<V_next>(key) = 0.0; | 
| 395 |  | 
| 396 |             ++it; | 
| 397 |     } | 
| 398 |  | 
| 399 |     long int x_start = grid.size(0)*1.55f/domain.getHigh(0); | 
| 400 |     long int y_start = grid.size(1)*1.55f/domain.getHigh(1); | 
| 401 |     long int z_start = grid.size(1)*1.55f/domain.getHigh(2); | 
| 402 |  | 
| 403 |     long int x_stop = grid.size(0)*1.85f/domain.getHigh(0); | 
| 404 |     long int y_stop = grid.size(1)*1.85f/domain.getHigh(1); | 
| 405 |     long int z_stop = grid.size(1)*1.85f/domain.getHigh(2); | 
| 406 |  | 
| 407 |     grid_key_dx<3> start({x_start,y_start,z_start}); | 
| 408 |     grid_key_dx<3> stop ({x_stop,y_stop,z_stop}); | 
| 409 |     auto it_init = grid.getGridIterator(start,stop); | 
| 410 |  | 
| 411 |     while (it_init.isNext()) | 
| 412 |     { | 
| 413 |             auto key = it_init.get_dist(); | 
| 414 |  | 
| 415 |             grid.template insert<U>(key) = 0.5 + (((double)std::rand())/RAND_MAX -0.5)/10.0; | 
| 416 |             grid.template insert<V>(key) = 0.25 + (((double)std::rand())/RAND_MAX -0.5)/20.0; | 
| 417 |  | 
| 418 |             ++it_init; | 
| 419 |     } | 
| 420 |  | 
| 421 |     // spacing of the grid on x and y | 
| 422 |     double spacing[3] = {grid.spacing(0),grid.spacing(1),grid.spacing(2)}; | 
| 423 |     // sync the ghost | 
| 424 |     size_t count = 0; | 
| 425 |     grid.template ghost_get<U,V>(); | 
| 426 |  | 
| 427 |     // because we assume that spacing[x] == spacing[y] we use formula 2 | 
| 428 |     // and we calculate the prefactor of Eq 2 | 
| 429 |     double uFactor = deltaT * du/(spacing[x]*spacing[x]); | 
| 430 |     double vFactor = deltaT * dv/(spacing[x]*spacing[x]); | 
| 431 |  | 
| 432 |     int stencil[6][3] = {{1,0,0},{-1,0,0},{0,-1,0},{0,1,0},{0,0,-1},{0,0,1}}; | 
| 433 |  | 
| 434 |  | 
| 435 |      //! \cond [stencil get and use] \endcond | 
| 436 |  | 
| 437 |  | 
| 438 |     auto func = [uFactor,vFactor,deltaT,F,K](Vc::double_v & u_out,Vc::double_v & v_out, | 
| 439 |                                                                 Vc::double_v (& u)[7],Vc::double_v (& v)[7], | 
| 440 |                                                                 unsigned char * mask){ | 
| 441 |  | 
| 442 | 																													 u_out = u[0] + uFactor *(u[1] + u[2] + | 
| 443 | 																																								  u[3] + u[4] + | 
| 444 | 																																								  u[5] + u[6] - 6.0*u[0]) - deltaT * u[0]*v[0]*v[0] | 
| 445 | 																																								- deltaT * F * (u[0] - 1.0); | 
| 446 |  | 
| 447 | 																													 v_out = v[0] + vFactor *(v[1] + v[2] + | 
| 448 | 																																								  v[3] + v[4] + | 
| 449 | 																																								  v[5] + v[6] - 6.0*v[0]) + deltaT * u[0]*v[0]*v[0] | 
| 450 | 																																								- deltaT * (F+K) * v[0]; | 
| 451 |                                                                                      }; | 
| 452 |  | 
| 453 |     grid.conv2<U,V,U_next,V_next,1>(stencil,{0,0,0},{(long int)sz[0]-1,(long int)sz[1]-1,(long int)sz[2]-1},func); | 
| 454 |     grid.conv2<U,V,U_next,V_next,1>(stencil,{0,0,0},{(long int)sz[0]-1,(long int)sz[1]-1,(long int)sz[2]-1},func); | 
| 455 |  | 
| 456 |     bool match = true; | 
| 457 |  | 
| 458 |     { | 
| 459 | 		auto it = grid.getDomainIterator(); | 
| 460 |  | 
| 461 | 		double max_U = 0.0; | 
| 462 | 		double max_V = 0.0; | 
| 463 | 		grid_dist_key_dx<3> k_max; | 
| 464 | 		while (it.isNext()) | 
| 465 | 		{ | 
| 466 | 			// center point | 
| 467 | 			auto Cp = it.get(); | 
| 468 |  | 
| 469 | 			// plus,minus X,Y,Z | 
| 470 | 			auto mx = Cp.move(0,-1); | 
| 471 | 			auto px = Cp.move(0,+1); | 
| 472 | 			auto my = Cp.move(1,-1); | 
| 473 | 			auto py = Cp.move(1,1); | 
| 474 | 			auto mz = Cp.move(2,-1); | 
| 475 | 			auto pz = Cp.move(2,1); | 
| 476 |  | 
| 477 | 			// update based on Eq 2 | 
| 478 | 			if ( fabs(grid.get<U>(Cp) + uFactor * ( | 
| 479 | 																	grid.get<U>(mz) + | 
| 480 | 																	grid.get<U>(pz) + | 
| 481 | 																	grid.get<U>(my) + | 
| 482 | 																	grid.get<U>(py) + | 
| 483 | 																	grid.get<U>(mx) + | 
| 484 | 																	grid.get<U>(px) - | 
| 485 | 																	6.0*grid.get<U>(Cp)) + | 
| 486 | 																	- deltaT * grid.get<U>(Cp) * grid.get<V>(Cp) * grid.get<V>(Cp) + | 
| 487 | 																	- deltaT * F * (grid.get<U>(Cp) - 1.0) - grid.get<U_next>(Cp)) > 0.000000001 ) | 
| 488 | 			{ | 
| 489 | 				match = false; | 
| 490 | 				break; | 
| 491 | 			} | 
| 492 |  | 
| 493 | 			// update based on Eq 2 | 
| 494 | 			if ( fabs(grid.get<V>(Cp) + vFactor * ( | 
| 495 | 																	grid.get<V>(mz) + | 
| 496 | 																	grid.get<V>(pz) + | 
| 497 | 																	grid.get<V>(my) + | 
| 498 | 																	grid.get<V>(py) + | 
| 499 | 																	grid.get<V>(mx) + | 
| 500 | 																	grid.get<V>(px) - | 
| 501 | 																	6*grid.get<V>(Cp)) + | 
| 502 | 																	deltaT * grid.get<U>(Cp) * grid.get<V>(Cp) * grid.get<V>(Cp) + | 
| 503 | 																	- deltaT * (F+K) * grid.get<V>(Cp) - grid.get<V_next>(Cp)) > 0.000000001 ) | 
| 504 | 			{ | 
| 505 | 				match = false; | 
| 506 | 				break; | 
| 507 | 			} | 
| 508 |  | 
| 509 | 			++it; | 
| 510 | 		} | 
| 511 |     } | 
| 512 |  | 
| 513 |     BOOST_REQUIRE_EQUAL(match,true); | 
| 514 | } | 
| 515 |  | 
| 516 |  | 
| 517 | BOOST_AUTO_TEST_CASE( sparse_grid_fast_stencil_vectorized_simplified_conv2_crossing) | 
| 518 | { | 
| 519 | 	constexpr int U = 0; | 
| 520 | 	constexpr int V = 1; | 
| 521 |  | 
| 522 | 	constexpr int U_next = 2; | 
| 523 | 	constexpr int V_next = 3; | 
| 524 |  | 
| 525 | 	constexpr int x = 0; | 
| 526 | 	constexpr int y = 1; | 
| 527 | 	constexpr int z = 2; | 
| 528 |  | 
| 529 |     Box<3,double> domain({0.0,0.0,0.0},{2.5,2.5,2.5}); | 
| 530 |  | 
| 531 |     // grid size | 
| 532 |     size_t sz[3] = {32,32,32}; | 
| 533 |  | 
| 534 |     // Define periodicity of the grid | 
| 535 |     periodicity<3> bc = {PERIODIC,PERIODIC,PERIODIC}; | 
| 536 |  | 
| 537 |     // Ghost in grid unit | 
| 538 |     Ghost<3,long int> g(1); | 
| 539 |  | 
| 540 |     // deltaT | 
| 541 |     double deltaT = 1; | 
| 542 |  | 
| 543 |     // Diffusion constant for specie U | 
| 544 |     double du = 2*1e-5; | 
| 545 |  | 
| 546 |     // Diffusion constant for specie V | 
| 547 |     double dv = 1*1e-5; | 
| 548 |  | 
| 549 |     // Number of timesteps | 
| 550 |     size_t timeSteps = 5000; | 
| 551 |  | 
| 552 |     // K and F (Physical constant in the equation) | 
| 553 |     double K = 0.053; | 
| 554 |     double F = 0.014; | 
| 555 |  | 
| 556 |     sgrid_dist_soa<3, double, aggregate<double,double,double,double>> grid(sz,domain,g,bc); | 
| 557 |  | 
| 558 |     auto it = grid.getGridIterator(); | 
| 559 |  | 
| 560 |     while (it.isNext()) | 
| 561 |     { | 
| 562 |             // Get the local grid key | 
| 563 |             auto key = it.get_dist(); | 
| 564 |  | 
| 565 |             // Old values U and V | 
| 566 |             grid.template insert<U>(key) = 1.0; | 
| 567 |             grid.template insert<V>(key) = 0.0; | 
| 568 |  | 
| 569 |             // Old values U and V | 
| 570 |             grid.template insert<U_next>(key) = 0.0; | 
| 571 |             grid.template insert<V_next>(key) = 0.0; | 
| 572 |  | 
| 573 |             ++it; | 
| 574 |     } | 
| 575 |  | 
| 576 |     long int x_start = grid.size(0)*1.55f/domain.getHigh(0); | 
| 577 |     long int y_start = grid.size(1)*1.55f/domain.getHigh(1); | 
| 578 |     long int z_start = grid.size(1)*1.55f/domain.getHigh(2); | 
| 579 |  | 
| 580 |     long int x_stop = grid.size(0)*1.85f/domain.getHigh(0); | 
| 581 |     long int y_stop = grid.size(1)*1.85f/domain.getHigh(1); | 
| 582 |     long int z_stop = grid.size(1)*1.85f/domain.getHigh(2); | 
| 583 |  | 
| 584 |     grid_key_dx<3> start({x_start,y_start,z_start}); | 
| 585 |     grid_key_dx<3> stop ({x_stop,y_stop,z_stop}); | 
| 586 |     auto it_init = grid.getGridIterator(start,stop); | 
| 587 |  | 
| 588 |     while (it_init.isNext()) | 
| 589 |     { | 
| 590 |             auto key = it_init.get_dist(); | 
| 591 |  | 
| 592 |             grid.template insert<U>(key) = 0.5 + (((double)std::rand())/RAND_MAX -0.5)/10.0; | 
| 593 |             grid.template insert<V>(key) = 0.25 + (((double)std::rand())/RAND_MAX -0.5)/20.0; | 
| 594 |  | 
| 595 |             ++it_init; | 
| 596 |     } | 
| 597 |  | 
| 598 |     // spacing of the grid on x and y | 
| 599 |     double spacing[3] = {grid.spacing(0),grid.spacing(1),grid.spacing(2)}; | 
| 600 |     // sync the ghost | 
| 601 |     size_t count = 0; | 
| 602 |     grid.template ghost_get<U,V>(); | 
| 603 |  | 
| 604 |     // because we assume that spacing[x] == spacing[y] we use formula 2 | 
| 605 |     // and we calculate the prefactor of Eq 2 | 
| 606 |     double uFactor = deltaT * du/(spacing[x]*spacing[x]); | 
| 607 |     double vFactor = deltaT * dv/(spacing[x]*spacing[x]); | 
| 608 |  | 
| 609 |  | 
| 610 |      //! \cond [stencil get and use] \endcond | 
| 611 |  | 
| 612 |  | 
| 613 |     auto func = [uFactor,vFactor,deltaT,F,K](Vc::double_v & u_out,Vc::double_v & v_out, | 
| 614 |     															Vc::double_v & u,Vc::double_v & v, | 
| 615 |                                                                 cross_stencil_v & us,cross_stencil_v & vs, | 
| 616 |                                                                 unsigned char * mask){ | 
| 617 |  | 
| 618 | 																														 u_out = u + uFactor *(us.xm + us.xp + | 
| 619 | 																																 	           us.ym + us.yp + | 
| 620 | 																																			   us.zm + us.zp - 6.0*u) - deltaT * u*v*v | 
| 621 | 																																									- deltaT * F * (u - 1.0); | 
| 622 |  | 
| 623 | 																														 v_out = v + vFactor *(vs.xm + vs.xp + | 
| 624 | 																																	  	  	   vs.ym + vs.yp + | 
| 625 | 																																			   vs.zm + vs.zp - 6.0*v) + deltaT * u*v*v | 
| 626 | 																																									- deltaT * (F+K) * v; | 
| 627 |                                                                                      }; | 
| 628 |  | 
| 629 |     grid.conv_cross2<U,V,U_next,V_next,1>({0,0,0},{(long int)sz[0]-1,(long int)sz[1]-1,(long int)sz[2]-1},func); | 
| 630 |     grid.conv_cross2<U,V,U_next,V_next,1>({0,0,0},{(long int)sz[0]-1,(long int)sz[1]-1,(long int)sz[2]-1},func); | 
| 631 |  | 
| 632 |     bool match = true; | 
| 633 |  | 
| 634 |     { | 
| 635 | 		auto it = grid.getDomainIterator(); | 
| 636 |  | 
| 637 | 		double max_U = 0.0; | 
| 638 | 		double max_V = 0.0; | 
| 639 | 		grid_dist_key_dx<3> k_max; | 
| 640 | 		while (it.isNext()) | 
| 641 | 		{ | 
| 642 | 			// center point | 
| 643 | 			auto Cp = it.get(); | 
| 644 |  | 
| 645 | 			// plus,minus X,Y,Z | 
| 646 | 			auto mx = Cp.move(0,-1); | 
| 647 | 			auto px = Cp.move(0,+1); | 
| 648 | 			auto my = Cp.move(1,-1); | 
| 649 | 			auto py = Cp.move(1,1); | 
| 650 | 			auto mz = Cp.move(2,-1); | 
| 651 | 			auto pz = Cp.move(2,1); | 
| 652 |  | 
| 653 | 			// update based on Eq 2 | 
| 654 | 			if ( fabs(grid.get<U>(Cp) + uFactor * ( | 
| 655 | 																	grid.get<U>(mz) + | 
| 656 | 																	grid.get<U>(pz) + | 
| 657 | 																	grid.get<U>(my) + | 
| 658 | 																	grid.get<U>(py) + | 
| 659 | 																	grid.get<U>(mx) + | 
| 660 | 																	grid.get<U>(px) - | 
| 661 | 																	6.0*grid.get<U>(Cp)) + | 
| 662 | 																	- deltaT * grid.get<U>(Cp) * grid.get<V>(Cp) * grid.get<V>(Cp) + | 
| 663 | 																	- deltaT * F * (grid.get<U>(Cp) - 1.0) - grid.get<U_next>(Cp)) > 0.000000001 ) | 
| 664 | 			{ | 
| 665 | 				match = false; | 
| 666 | 				break; | 
| 667 | 			} | 
| 668 |  | 
| 669 | 			// update based on Eq 2 | 
| 670 | 			if ( fabs(grid.get<V>(Cp) + vFactor * ( | 
| 671 | 																	grid.get<V>(mz) + | 
| 672 | 																	grid.get<V>(pz) + | 
| 673 | 																	grid.get<V>(my) + | 
| 674 | 																	grid.get<V>(py) + | 
| 675 | 																	grid.get<V>(mx) + | 
| 676 | 																	grid.get<V>(px) - | 
| 677 | 																	6*grid.get<V>(Cp)) + | 
| 678 | 																	deltaT * grid.get<U>(Cp) * grid.get<V>(Cp) * grid.get<V>(Cp) + | 
| 679 | 																	- deltaT * (F+K) * grid.get<V>(Cp) - grid.get<V_next>(Cp)) > 0.000000001 ) | 
| 680 | 			{ | 
| 681 | 				match = false; | 
| 682 | 				break; | 
| 683 | 			} | 
| 684 |  | 
| 685 | 			++it; | 
| 686 | 		} | 
| 687 |     } | 
| 688 |  | 
| 689 |     BOOST_REQUIRE_EQUAL(match,true); | 
| 690 | } | 
| 691 |  | 
| 692 | BOOST_AUTO_TEST_CASE (sgrid_dist_id_soa_write ) | 
| 693 | { | 
| 694 | 	periodicity<3> bc = {PERIODIC, PERIODIC, PERIODIC}; | 
| 695 |  | 
| 696 | 	auto & v_cl = create_vcluster<>(); | 
| 697 |  | 
| 698 | 	if (v_cl.size() > 16) | 
| 699 | 	{return;} | 
| 700 |  | 
| 701 | 	// Domain | 
| 702 | 	Box<3,double> domain({-0.3,-0.3,-0.3},{1.0,1.0,1.0}); | 
| 703 |  | 
| 704 | 	// grid size | 
| 705 | 	size_t sz[3]; | 
| 706 | 	sz[0] = 256; | 
| 707 | 	sz[1] = 256; | 
| 708 | 	sz[2] = 256; | 
| 709 |  | 
| 710 | 	// Ghost | 
| 711 | 	Ghost<3,long int> g(1); | 
| 712 |  | 
| 713 | 	sgrid_dist_soa<3,double,aggregate<double,double[3]>> sg1(sz,domain,g,bc); | 
| 714 | 	sgrid_dist_id<3,double,aggregate<double,double[3]>> sg2(sg1.getDecomposition(),sz,g); | 
| 715 |  | 
| 716 | 	// create a grid iterator over a bilion point | 
| 717 |  | 
| 718 | 	auto it = sg1.getGridIterator(); | 
| 719 |  | 
| 720 | 	while(it.isNext()) | 
| 721 | 	{ | 
| 722 | 		auto gkey = it.get(); | 
| 723 | 		auto key = it.get_dist(); | 
| 724 |  | 
| 725 | 		size_t sx = gkey.get(0) - 128; | 
| 726 | 		size_t sy = gkey.get(1) - 128; | 
| 727 | 		size_t sz = gkey.get(2) - 128; | 
| 728 |  | 
| 729 | 		if (sx*sx + sy*sy + sz*sz < 32*32) | 
| 730 | 		{ | 
| 731 | 			sg1.template insert<0>(key) = 1.0; | 
| 732 | 			sg1.template insert<1>(key)[0] = gkey.get(0); | 
| 733 | 			sg1.template insert<1>(key)[1] = gkey.get(1); | 
| 734 | 			sg1.template insert<1>(key)[2] = gkey.get(2); | 
| 735 |  | 
| 736 | 			sg2.template insert<0>(key) = 1.0; | 
| 737 | 			sg2.template insert<1>(key)[0] = gkey.get(0); | 
| 738 | 			sg2.template insert<1>(key)[1] = gkey.get(1); | 
| 739 | 			sg2.template insert<1>(key)[2] = gkey.get(2); | 
| 740 | 		} | 
| 741 |  | 
| 742 | 		++it; | 
| 743 | 	} | 
| 744 |  | 
| 745 | 	sg1.write("sg1_test" ); | 
| 746 | 	sg2.write("sg2_test" ); | 
| 747 |  | 
| 748 | 	bool test = compare("sg1_test_"  + std::to_string(v_cl.rank()) + ".vtk" ,"sg2_test_"  + std::to_string(v_cl.rank()) + ".vtk" ); | 
| 749 | 	BOOST_REQUIRE_EQUAL(true,test); | 
| 750 | } | 
| 751 |  | 
| 752 | BOOST_AUTO_TEST_SUITE_END() | 
| 753 |  |