1 | /* SolarSim.c */ |
---|
2 | #include "SolarSim.h" |
---|
3 | |
---|
4 | typedef struct task { |
---|
5 | int i, j; |
---|
6 | } Task, *TaskPtr; |
---|
7 | |
---|
8 | int sock = 0; |
---|
9 | long SIM_STEPS, DT, SAVE; |
---|
10 | |
---|
11 | void |
---|
12 | do_work(AState state) { |
---|
13 | int tid, nthreads; |
---|
14 | long i, j, t, ntasks, sim_step; |
---|
15 | int n = state->particles; |
---|
16 | vector *r = state->positions; |
---|
17 | vector *v = state->velocities; |
---|
18 | vector *a = state->accelerations; |
---|
19 | vector *ac; |
---|
20 | real *m = state->masses; |
---|
21 | TaskPtr tasks; |
---|
22 | |
---|
23 | tasks = malloc( n*(n-1)/2 * sizeof(Task) ); |
---|
24 | t = 0; |
---|
25 | for( i = 0; i < n-1; i++ ) |
---|
26 | for( j = i+1; j < n; j++ ) { |
---|
27 | tasks[t].i = i; |
---|
28 | tasks[t].j = j; |
---|
29 | t++; |
---|
30 | } |
---|
31 | ntasks = t; |
---|
32 | |
---|
33 | #pragma omp parallel shared(nthreads) |
---|
34 | #pragma omp master |
---|
35 | nthreads = omp_get_num_threads(); |
---|
36 | |
---|
37 | ac = malloc(nthreads * n * sizeof(vector)); |
---|
38 | |
---|
39 | #pragma omp parallel shared(r,v,a,ac,m,n,tasks,ntasks,nthreads,DT,sim_step,SIM_STEPS,SAVE) private(i,j,t,tid) default(none) |
---|
40 | sim_step = 0; |
---|
41 | while( sim_step < SIM_STEPS ) { |
---|
42 | if( (SAVE != 0) && (sim_step % SAVE == 0) ) { |
---|
43 | write_state(sock, state); |
---|
44 | } |
---|
45 | |
---|
46 | #pragma omp for schedule(static,n) //shared(ac,n,nthreads) private(i) default(none) nowait |
---|
47 | for( i = 0; i < n*nthreads; i++ ) { |
---|
48 | ac[i].x = 0; |
---|
49 | ac[i].y = 0; |
---|
50 | ac[i].z = 0; |
---|
51 | } |
---|
52 | |
---|
53 | // compute accelerations. each thread uses ac[tid] vector |
---|
54 | #pragma omp for schedule(static) //shared(tasks, ntasks, r, v, ac, m, n) private(i,j,t,tid) default(none) nowait |
---|
55 | for( t = 0; t < ntasks; t++ ) { |
---|
56 | i = tasks[t].i; |
---|
57 | j = tasks[t].j; |
---|
58 | tid = omp_get_thread_num(); |
---|
59 | compute_acc(r+i, r+j, v+i, v+j, ac+tid*n+i, ac+tid*n+j, m[i], m[j]); |
---|
60 | } |
---|
61 | |
---|
62 | // add up the accelerations computed by each thread |
---|
63 | #pragma omp parallel for schedule(static,1) //shared(a,ac,n,nthreads) private(i,j) default(none) nowait |
---|
64 | for( i = 0; i < n; i++ ) { |
---|
65 | v_init(a+i,0,0,0); |
---|
66 | for( j = 0; j < nthreads; j++ ) { |
---|
67 | a[i].x += ac[j*n+i].x; |
---|
68 | a[i].y += ac[j*n+i].y; |
---|
69 | a[i].z += ac[j*n+i].z; |
---|
70 | } |
---|
71 | } |
---|
72 | |
---|
73 | // update positions and speeds |
---|
74 | #pragma omp parallel for schedule(static,1) //shared(r,v,a,n,DT) private(i) default(none) nowait |
---|
75 | for( i = 0; i < n; i++ ) { |
---|
76 | update_pos_vel(r+i, v+i, a+i, DT); |
---|
77 | } |
---|
78 | |
---|
79 | #pragma omp master |
---|
80 | sim_step ++; |
---|
81 | |
---|
82 | #pragma omp barrier |
---|
83 | } |
---|
84 | if( SAVE != 0 ) |
---|
85 | write_state(sock, state); |
---|
86 | } |
---|
87 | |
---|
88 | int |
---|
89 | main( int argc, char **argv ) { |
---|
90 | TState state; |
---|
91 | double before, after; |
---|
92 | int i, n=8; |
---|
93 | |
---|
94 | SIM_STEPS = atol(argv[3]); |
---|
95 | DT = atol(argv[4]); |
---|
96 | SAVE = atol(argv[5]); |
---|
97 | |
---|
98 | if( SAVE != 0 ) { |
---|
99 | sock = open( argv[2], O_WRONLY | O_CREAT | O_TRUNC, 644 ); |
---|
100 | } |
---|
101 | |
---|
102 | // read_initial( &state, argv[1] ); |
---|
103 | rand_data( &state, 100 ); |
---|
104 | |
---|
105 | before = omp_get_wtime(); |
---|
106 | do_work(&state); |
---|
107 | after = omp_get_wtime(); |
---|
108 | |
---|
109 | if( SAVE != 0 ) |
---|
110 | close(sock); |
---|
111 | |
---|
112 | #pragma omp parallel |
---|
113 | #pragma omp master |
---|
114 | printf("%d %lf\n", omp_get_num_threads(), after-before); |
---|
115 | |
---|
116 | return EXIT_SUCCESS; |
---|
117 | } |
---|
118 | |
---|