Révision 162 NBody/NBodyGL.py
| NBodyGL.py (revision 162) | ||
|---|---|---|
| 158 | 158 |
xf=x0+dt*(v1+(MYFLOAT)2.e0f*(v2+v3)+v4)/(MYFLOAT)6.e0f; |
| 159 | 159 |
vf=v0+dt*(a1+(MYFLOAT)2.e0f*(a2+a3)+a4)/(MYFLOAT)6.e0f; |
| 160 | 160 |
|
| 161 |
return((MYFLOAT8)(xf.s0,xf.s1,xf.s2,1.e0f,vf.s0,vf.s1,vf.s2,1.e0f));
|
|
| 161 |
return((MYFLOAT8)(xf.s0,xf.s1,xf.s2,1.e0f,vf.s0,vf.s1,vf.s2,0.e0f));
|
|
| 162 | 162 |
} |
| 163 | 163 |
|
| 164 | 164 |
// Elements from : http://doswa.com/2009/01/02/fourth-order-runge-kutta-numerical-integration.html |
| ... | ... | |
| 236 | 236 |
__kernel void SplutterPoints(__global MYFLOAT4* clDataX, MYFLOAT box, |
| 237 | 237 |
uint seed_z,uint seed_w) |
| 238 | 238 |
{
|
| 239 |
uint gid=(uint)get_global_id(0);
|
|
| 240 |
uint z=(seed_z+gid)%gid;
|
|
| 239 |
int gid=get_global_id(0);
|
|
| 240 |
uint z=seed_z+gid;
|
|
| 241 | 241 |
uint w=seed_w+gid; |
| 242 |
|
|
| 243 |
// cast to float for sin,cos are NEEDED by Mesa FP64 implementation! |
|
| 242 |
|
|
| 243 |
for (int i=0;i<gid;i++) |
|
| 244 |
{
|
|
| 245 |
private MYFLOAT heat=MWCfp; |
|
| 246 |
} |
|
| 247 |
|
|
| 248 |
// Distribute in sphere |
|
| 244 | 249 |
MYFLOAT radius=MWCfp*box; |
| 245 | 250 |
MYFLOAT theta=MWCfp*PI; |
| 246 | 251 |
MYFLOAT phi=MWCfp*PI*(MYFLOAT)2.e0f; |
| 247 |
MYFLOAT sinTheta=sin((float)theta);
|
|
| 248 |
|
|
| 252 |
// cast to float for sin,cos are NEEDED by Mesa FP64 implementation!
|
|
| 253 |
MYFLOAT sinTheta=sin((float)theta); |
|
| 249 | 254 |
clDataX[gid].s0=radius*sinTheta*cos((float)phi)/2.e0f; |
| 250 | 255 |
clDataX[gid].s1=radius*sinTheta*sin((float)phi)/2.e0f; |
| 251 | 256 |
clDataX[gid].s2=radius*cos((float)theta)/2.e0f; |
| 252 |
clDataX[gid].s3=(MYFLOAT)1.e0f;
|
|
| 257 |
clDataX[gid].s3=(MYFLOAT)0.e0f;
|
|
| 253 | 258 |
barrier(CLK_GLOBAL_MEM_FENCE); |
| 254 | 259 |
|
| 255 |
//MYFLOAT x0=box*(MYFLOAT)(MWCfp-0.5); |
|
| 256 |
//MYFLOAT y0=box*(MYFLOAT)(MWCfp-0.5); |
|
| 257 |
//MYFLOAT z0=box*(MYFLOAT)(MWCfp-0.5); |
|
| 258 |
//clDataX[gid].s0123 = (MYFLOAT4) (x0,y0,z0,1.e0f); |
|
| 260 |
// Distribute in box |
|
| 261 |
// MYFLOAT x0=box*(MYFLOAT)(MWCfp-0.5); |
|
| 262 |
// MYFLOAT y0=box*(MYFLOAT)(MWCfp-0.5); |
|
| 263 |
// MYFLOAT z0=box*(MYFLOAT)(MWCfp-0.5); |
|
| 264 |
// clDataX[gid].s0123 = (MYFLOAT4) (x0,y0,z0,0.e0f); |
|
| 259 | 265 |
} |
| 260 | 266 |
|
| 261 | 267 |
__kernel void SplutterStress(__global MYFLOAT4* clDataX,__global MYFLOAT4* clDataV,__global MYFLOAT4* clCoM, MYFLOAT velocity,uint seed_z,uint seed_w) |
| ... | ... | |
| 279 | 285 |
clDataV[gid].s0=velocity*sinTheta*cos((float)phi); |
| 280 | 286 |
clDataV[gid].s1=velocity*sinTheta*sin((float)phi); |
| 281 | 287 |
clDataV[gid].s2=velocity*cos((float)theta); |
| 282 |
clDataV[gid].s3=(MYFLOAT)1.e0f;
|
|
| 288 |
clDataV[gid].s3=(MYFLOAT)0.e0f;
|
|
| 283 | 289 |
} |
| 284 | 290 |
barrier(CLK_GLOBAL_MEM_FENCE); |
| 285 | 291 |
} |
| ... | ... | |
| 373 | 379 |
|
| 374 | 380 |
def display(*args): |
| 375 | 381 |
|
| 376 |
global MyDataX,clDataX,clDataV,Step,Method,Number
|
|
| 382 |
global MyDataX,MyDataV,clDataX,clDataV,Step,Method,Number,Iterations
|
|
| 377 | 383 |
|
| 378 | 384 |
glClearColor(0.0, 0.0, 0.0, 0.0) |
| 379 | 385 |
glClear(GL_COLOR_BUFFER_BIT) |
| 380 |
glColor3f(1.0,1.0,0.0)
|
|
| 381 |
|
|
| 386 |
glColor3f(1.0,1.0,1.0)
|
|
| 387 |
|
|
| 382 | 388 |
time_start=time.time() |
| 383 | 389 |
if Method=="RungeKutta": |
| 384 | 390 |
CLLaunch=MyRoutines.RungeKutta(queue,(Number,1),None,clDataX,clDataV,Step) |
| ... | ... | |
| 389 | 395 |
else: |
| 390 | 396 |
CLLaunch=MyRoutines.ImplicitEuler(queue,(Number,1),None,clDataX,clDataV,Step) |
| 391 | 397 |
CLLaunch.wait() |
| 392 |
print("Duration of 1 iteration %s" % (time.time()-time_start))
|
|
| 398 |
print("Duration of #%s iteration: %s" % (Iterations,(time.time()-time_start)))
|
|
| 393 | 399 |
|
| 394 | 400 |
cl.enqueue_copy(queue, MyDataX, clDataX) |
| 395 |
print(MyDataX.reshape(Number,4)) |
|
| 401 |
# print(MyDataX.reshape(Number,4)) |
|
| 402 |
MyDataX.reshape(Number,4)[:,3]=1 |
|
| 396 | 403 |
glVertexPointerf(MyDataX.reshape(Number,4)) |
| 404 |
# cl.enqueue_copy(queue, MyDataV, clDataV) |
|
| 405 |
# print(MyDataV.reshape(Number,4)) |
|
| 406 |
# MyDataV.reshape(Number,4)[:,3]=1 |
|
| 407 |
# glVertexPointerf(MyDataV.reshape(Number,4)) |
|
| 397 | 408 |
glEnableClientState(GL_VERTEX_ARRAY) |
| 398 | 409 |
glDrawArrays(GL_POINTS, 0, Number) |
| 399 | 410 |
glDisableClientState(GL_VERTEX_ARRAY) |
| 400 | 411 |
glFlush() |
| 412 |
Iterations+=1 |
|
| 401 | 413 |
glutSwapBuffers() |
| 402 | 414 |
|
| 403 | 415 |
def halt(): |
| 404 | 416 |
pass |
| 405 | 417 |
|
| 406 |
def keyboard(*args): |
|
| 407 |
return() |
|
| 418 |
def keyboard(k, x, y): |
|
| 419 |
global view_rotz |
|
| 420 |
LC_Z = as_8_bit( 'z' ) |
|
| 421 |
UC_Z = as_8_bit( 'Z' ) |
|
| 408 | 422 |
|
| 423 |
if k == LC_Z: |
|
| 424 |
view_rotz += 1.0 |
|
| 425 |
elif k == UC_Z: |
|
| 426 |
view_rotz -= 1.0 |
|
| 427 |
elif ord(k) == 27: # Escape |
|
| 428 |
glutSetOption(GLUT_ACTION_ON_WINDOW_CLOSE,GLUT_ACTION_CONTINUE_EXECUTION) |
|
| 429 |
glutSetOption(GLUT_ACTION_GLUTMAINLOOP_RETURNS,GLUT_ACTION_CONTINUE_EXECUTION) |
|
| 430 |
glutLeaveMainLoop() |
|
| 431 |
return(False) |
|
| 432 |
else: |
|
| 433 |
return |
|
| 434 |
glRotatef(view_rotz, 0.0, 0.0, 1.0) |
|
| 435 |
glutPostRedisplay() |
|
| 436 |
|
|
| 437 |
def special(k, x, y): |
|
| 438 |
global view_rotx, view_roty, view_rotz |
|
| 439 |
|
|
| 440 |
if k == GLUT_KEY_UP: |
|
| 441 |
view_rotx += 1.0 |
|
| 442 |
elif k == GLUT_KEY_DOWN: |
|
| 443 |
view_rotx -= 1.0 |
|
| 444 |
elif k == GLUT_KEY_LEFT: |
|
| 445 |
view_roty += 1.0 |
|
| 446 |
elif k == GLUT_KEY_RIGHT: |
|
| 447 |
view_roty -= 1.0 |
|
| 448 |
else: |
|
| 449 |
return |
|
| 450 |
glRotatef(view_rotx, 1.0, 0.0, 0.0) |
|
| 451 |
glRotatef(view_roty, 0.0, 1.0, 0.0) |
|
| 452 |
glutPostRedisplay() |
|
| 453 |
|
|
| 409 | 454 |
def mouse(button, state, x, y): |
| 410 | 455 |
global angle, delta_angle, halted |
| 411 | 456 |
if button == GLUT_LEFT_BUTTON: |
| ... | ... | |
| 425 | 470 |
glMatrixMode(GL_PROJECTION) |
| 426 | 471 |
glLoadIdentity() |
| 427 | 472 |
glOrtho(-SizeOfBox, SizeOfBox, -SizeOfBox, SizeOfBox, -SizeOfBox, SizeOfBox) |
| 428 |
|
|
| 473 |
|
|
| 429 | 474 |
def reshape(w, h): |
| 430 | 475 |
glViewport(0, 0, w, h) |
| 431 | 476 |
setup_viewport() |
| 432 | 477 |
|
| 433 |
def main(): |
|
| 434 |
|
|
| 435 |
glutInit(sys.argv) |
|
| 436 |
glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB) |
|
| 437 |
glutInitWindowSize(400,400) |
|
| 438 |
glutCreateWindow(b'NBodyGL') |
|
| 439 |
setup_viewport() |
|
| 440 |
# glutReshapeFunc(reshape) |
|
| 441 |
glutDisplayFunc(display) |
|
| 442 |
glutIdleFunc(display) |
|
| 443 |
# glutMouseFunc(mouse) |
|
| 444 |
glutKeyboardFunc(keyboard) |
|
| 445 |
glutMainLoop() |
|
| 446 |
|
|
| 447 | 478 |
if __name__=='__main__': |
| 448 | 479 |
|
| 449 | 480 |
global Number,Step,clDataX,clDataV,MyDataX,MyDataV,Method,SizeOfBox |
| ... | ... | |
| 593 | 624 |
|
| 594 | 625 |
print(Marsaglia[RNG],Computing[ValueType]) |
| 595 | 626 |
# Build all routines used for the computing |
| 596 |
#BuildOptions="-DTRNG=%i -DTYPE=%i" % (Marsaglia[RNG],Computing[ValueType]) |
|
| 597 |
#BuildOptions="-cl-mad-enable -cl-fast-relaxed-math -DTRNG=%i -DTYPE=%i" % (Marsaglia[RNG],Computing[ValueType]) |
|
| 598 | 627 |
BuildOptions="-cl-mad-enable -cl-kernel-arg-info -cl-fast-relaxed-math -cl-std=CL1.2 -DTRNG=%i -DTYPE=%i" % (Marsaglia[RNG],Computing[ValueType]) |
| 599 | 628 |
|
| 600 |
if 'Intel' in PlatForm.name or 'Clover' in PlatForm.name or 'Portable' in PlatForm.name : |
|
| 629 |
if 'Intel' in PlatForm.name or 'Experimental' in PlatForm.name or 'Clover' in PlatForm.name or 'Portable' in PlatForm.name :
|
|
| 601 | 630 |
MyRoutines = cl.Program(ctx, BlobOpenCL).build(options = BuildOptions) |
| 602 | 631 |
else: |
| 603 | 632 |
MyRoutines = cl.Program(ctx, BlobOpenCL).build(options = BuildOptions+" -cl-strict-aliasing") |
| 604 | 633 |
|
| 605 | 634 |
mf = cl.mem_flags |
| 635 |
# Read/Write approach for buffering |
|
| 606 | 636 |
# clDataX = cl.Buffer(ctx, mf.READ_WRITE, MyDataX.nbytes) |
| 607 | 637 |
# clDataV = cl.Buffer(ctx, mf.READ_WRITE, MyDataV.nbytes) |
| 608 | 638 |
# clPotential = cl.Buffer(ctx, mf.READ_WRITE, MyPotential.nbytes) |
| 609 | 639 |
# clKinetic = cl.Buffer(ctx, mf.READ_WRITE, MyKinetic.nbytes) |
| 610 | 640 |
# clCoM = cl.Buffer(ctx, mf.READ_WRITE, MyCoM.nbytes) |
| 611 | 641 |
|
| 642 |
# Write/HostPoniter approach for buffering |
|
| 612 | 643 |
clDataX = cl.Buffer(ctx, mf.WRITE_ONLY|mf.COPY_HOST_PTR,hostbuf=MyDataX) |
| 613 | 644 |
clDataV = cl.Buffer(ctx, mf.WRITE_ONLY|mf.COPY_HOST_PTR,hostbuf=MyDataV) |
| 614 | 645 |
clPotential = cl.Buffer(ctx, mf.WRITE_ONLY|mf.COPY_HOST_PTR,hostbuf=MyPotential) |
| ... | ... | |
| 617 | 648 |
|
| 618 | 649 |
print('All particles superimposed.')
|
| 619 | 650 |
|
| 620 |
print(SizeOfBox.dtype) |
|
| 621 |
|
|
| 622 | 651 |
# Set particles to RNG points |
| 623 | 652 |
if InitialRandom: |
| 624 | 653 |
MyRoutines.SplutterPoints(queue,(Number,1),None,clDataX,SizeOfBox,np.uint32(nprnd(2**32)),np.uint32(nprnd(2**32))) |
| ... | ... | |
| 657 | 686 |
# t1=np.array([[MyDataX[1][0],MyDataX[1][1],MyDataX[1][2]]]) |
| 658 | 687 |
# tL=np.array([[MyDataX[-1][0],MyDataX[-1][1],MyDataX[-1][2]]]) |
| 659 | 688 |
|
| 660 |
time_start=time.time() |
|
| 689 |
wall_time_start=time.time() |
|
| 690 |
|
|
| 661 | 691 |
print("Use the mouse buttons to control some of the dots.")
|
| 662 | 692 |
print("Hit any key to quit.")
|
| 663 |
main() |
|
| 664 |
print("\nDuration on %s for each %s\n" % (Device,(time.time()-time_start)/Iterations))
|
|
| 693 |
global view_rotx,view_roty,view_rotz,Iterations |
|
| 694 |
(view_rotx,view_roty,view_rotz)=(0.0, 0.0, 0.0) |
|
| 695 |
Iterations=0 |
|
| 696 |
glutInit(sys.argv) |
|
| 697 |
glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB) |
|
| 698 |
glutInitWindowSize(512,512) |
|
| 699 |
glutCreateWindow(b'NBodyGL') |
|
| 700 |
setup_viewport() |
|
| 701 |
glutReshapeFunc(reshape) |
|
| 702 |
glutDisplayFunc(display) |
|
| 703 |
glutIdleFunc(display) |
|
| 704 |
# glutMouseFunc(mouse) |
|
| 705 |
glutSpecialFunc(special) |
|
| 706 |
Loop=glutKeyboardFunc(keyboard) |
|
| 707 |
glutMainLoop() |
|
| 665 | 708 |
|
| 709 |
print("\nWall Duration on %s for each %s\n" % (Device,(time.time()-wall_time_start)/Iterations))
|
|
| 710 |
|
|
| 666 | 711 |
MyRoutines.CenterOfMass(queue,(1,1),None,clDataX,clCoM,np.int32(Number)) |
| 667 | 712 |
CLLaunch=MyRoutines.Potential(queue,(Number,1),None,clDataX,clPotential) |
| 668 | 713 |
CLLaunch=MyRoutines.Kinetic(queue,(Number,1),None,clDataV,clKinetic) |
Formats disponibles : Unified diff