source: trunk/anuga_work/anuga_cuda/src/anuga_HMPP/gravity.c @ 9329

Last change on this file since 9329 was 9017, checked in by steve, 12 years ago

Adding in Zhe (John) Weng's anuga_cuda code as obtained from googlecode ​https://code.google.com/p/anuga-cuda
File size: 16.3 KB
Line 
1
2// OpenHMPP ANUGA grvity function
3//
4// Zhe Weng 2013
5
6//#define ISOLATING_TEST
7
8#ifndef ISOLATING_TEST
9#include "hmpp_fun.h"
10#else // ISOLATING_TEST
11#define DATA_TYPE double
12#define TOLERANCE 1e-11
13#ifdef USING_CPP
14#include <iostream>
15#include <cstdio>
16#include <cstdlib>
17#include <cmath>
18using namespace std;
19
20#else // USING_CPP
21#include <stdio.h>
22#include <math.h>
23#include <stdlib.h>
24#endif
25
26int errs_x =0, errs_y =0;
27#endif
28
29
30int check_tolerance( DATA_TYPE a, DATA_TYPE b)
31{
32    if (fabs(a) <= 0 && fabs(a-b) > TOLERANCE)
33        return 1;
34    else if(fabs(a) < 10 && fabs(a-b) > TOLERANCE*10)
35        return 1;
36    else if(fabs(a) < 100 && fabs(a-b) > TOLERANCE*100)
37        return 1;
38    else if(fabs(a) < 1000 && fabs(a-b) > TOLERANCE*1000)
39        return 1;
40    else if(fabs(a) < 10000 && fabs(a-b) > TOLERANCE*10000)
41        return 1;
42    else return 0;
43}
44       
45
46#ifdef USING_LOCAL_DIRECTIVES
47#pragma hmpp gravity codelet, target=CUDA args[*].transfer=atcall
48#endif
49void gravity_wb( 
50        int n, int n3, int n6, 
51        DATA_TYPE xmom_explicit_update[n], 
52        DATA_TYPE ymom_explicit_update[n], 
53
54        DATA_TYPE stage_vertex_values[n3],
55        DATA_TYPE stage_edge_values[n3],
56        DATA_TYPE stage_centroid_values[n],
57
58        DATA_TYPE bed_edge_values[n3],
59        DATA_TYPE bed_centroid_values[n],
60
61        DATA_TYPE vertex_coordinates[n6],
62
63        DATA_TYPE normals[n6],
64        DATA_TYPE areas[n],
65        DATA_TYPE edgelengths[n3],
66
67        DATA_TYPE g )
68{
69    int k;
70    int k3, k6;
71    double w0, w1, w2,
72           x0, y0, x1, y1, x2, y2,
73           avg_h;
74
75    double wx, wy, det, hh0, hh1,hh2;
76    double sidex_0, sidex_1, sidex_2;
77    double sidey_0, sidey_1, sidey_2;
78    double area;
79
80
81
82    for( k = 0; k < n; ++k ){
83
84        k3 = k*3;
85        w0 = stage_vertex_values[k3];
86        w1 = stage_vertex_values[k3 + 1];
87        w2 = stage_vertex_values[k3 + 2];
88
89        k6 = k*6;
90        x0 = vertex_coordinates[k6 + 0];
91        y0 = vertex_coordinates[k6 + 1];
92        x1 = vertex_coordinates[k6 + 2];
93        y1 = vertex_coordinates[k6 + 3];
94        x2 = vertex_coordinates[k6 + 4];
95        y2 = vertex_coordinates[k6 + 5];
96
97
98        det = (y2 - y0)*(x1 - x0) - (y1 - y0)*(x2 - x0);
99
100        wx = (y2 -y0)*(w1 - w0) - (y1 - y0)*(w2 -w0);
101        wx /= det;
102
103        wy = (x1 - x0)*(w2 - w0) - (x2 - x0)*(w1 -w0);
104        wy /= det;
105
106        avg_h = stage_centroid_values[k] - bed_centroid_values[k];
107
108        xmom_explicit_update[k] += -g * wx * avg_h;
109        ymom_explicit_update[k] += -g * wy * avg_h;
110
111
112
113
114        hh0 = stage_edge_values[k3] - bed_edge_values[k3];
115        hh1 = stage_edge_values[k3+1] - bed_edge_values[k3+1];
116        hh2 = stage_edge_values[k3+2] - bed_edge_values[k3+2];
117
118        area = areas[k];
119       
120        sidex_0 = hh0*hh0*edgelengths[k3] * normals[k6];
121        sidex_1 = hh1*hh1*edgelengths[k3+1] * normals[k6+2];
122        sidex_2 = hh2*hh2*edgelengths[k3+2] * normals[k6+4];
123        xmom_explicit_update[k] += 0.5*g*(sidex_0 + sidex_1 + sidex_2)/area;
124
125
126        sidey_0 = hh0*hh0*edgelengths[k3] * normals[k6+1];
127        sidey_1 = hh1*hh1*edgelengths[k3+1] * normals[k6+3];
128        sidey_2 = hh2*hh2*edgelengths[k3+2] * normals[k6+5];
129        ymom_explicit_update[k] += 0.5*g*(sidey_0 + sidey_1 + sidey_2)/area;
130    }
131}
132
133
134
135void gravity_wb_orig(
136        DATA_TYPE * xmom_explicit_update, 
137        DATA_TYPE * ymom_explicit_update, 
138        DATA_TYPE * stage_vertex_values, 
139        DATA_TYPE * stage_edge_values, 
140        DATA_TYPE * stage_centroid_values, 
141        DATA_TYPE * bed_edge_values, 
142        DATA_TYPE * bed_centroid_values, 
143        DATA_TYPE * vertex_coordinates, 
144        DATA_TYPE * normals, 
145        DATA_TYPE * areas, 
146        DATA_TYPE * edgelengths,
147        DATA_TYPE * test_xe,
148        DATA_TYPE * test_ye,
149        int N,
150        DATA_TYPE g
151        )
152{
153    int i, k, k3, k6;
154    // For testing purpose
155    int errs_x=0, errs_y=0;
156
157    DATA_TYPE w0, w1, w2, 
158           x0, y0, x1, y1, x2, y2,
159           avg_h;
160
161    DATA_TYPE wx, wy, det,
162           hh[3];
163    DATA_TYPE sidex, sidey, area, n0, n1, fact;
164
165    for (k = 0; k < N; k++)
166    {
167        k3 = 3*k;
168        w0 = stage_vertex_values[k3 + 0];
169        w1 = stage_vertex_values[k3 + 1];
170        w2 = stage_vertex_values[k3 + 2];
171
172
173        k6 = 6*k;
174        x0 = vertex_coordinates[k6 + 0];
175        y0 = vertex_coordinates[k6 + 1];
176        x1 = vertex_coordinates[k6 + 2];
177        y1 = vertex_coordinates[k6 + 3];
178        x2 = vertex_coordinates[k6 + 4];
179        y2 = vertex_coordinates[k6 + 5];
180
181
182        //_gradient(x0, y0, x1, y1, x2, y2, w0, w1, w2, &wx, &wy);
183
184        det = (y2 - y0)*(x1 - x0) - (y1 - y0)*(x2 - x0);
185
186        wx = (y2 -y0)*(w1 - w0) - (y1 - y0)*(w2 -w0);
187        wx /= det;
188
189        wy = (x1 - x0)*(w2 - w0) - (x2 - x0)*(w1 -w0);
190        wy /= det;
191
192
193        avg_h = stage_centroid_values[k] - bed_centroid_values[k];
194
195        test_xe[k] += -g * wx * avg_h;
196        test_ye[k] += -g * wy * avg_h;
197
198
199
200        hh[0] = stage_edge_values[k3] - bed_edge_values[k3];
201        hh[1] = stage_edge_values[k3+1] - bed_edge_values[k3+1];
202        hh[2] = stage_edge_values[k3+2] - bed_edge_values[k3+2];
203
204        sidex = 0.0;
205        sidey = 0.0;
206
207        for ( i = 0 ; i < 3 ; i++ )
208        {
209            n0 = normals[k6 + 2*i];
210            n1 = normals[k6 + 2*i + 1];
211
212            fact =  -0.5 * g * hh[i] * hh[i] * edgelengths[k3 + i];
213
214            sidex += fact*n0;
215            sidey += fact*n1;
216        }
217
218        area = areas[k];
219        test_xe[k] += -sidex / area;
220        test_ye[k] += -sidey / area;
221
222        // For testing purpose
223        //test_xe[k] = hh[0];
224        //test_ye[k] = hh[2];
225        if ( check_tolerance(xmom_explicit_update[k], test_xe[k]) )
226        {
227            errs_x += 1;
228            if (errs_x <= 10)
229                printf("   Errors on xe: %d %.9f %.9f\n",
230                        k, xmom_explicit_update[k], test_xe[k]);
231        }
232        if ( check_tolerance(ymom_explicit_update[k], test_ye[k]) )
233        {
234            errs_y += 1;
235            if (errs_y <= 10)
236                printf("   Errors on ye: %d %f %f\n", 
237                        k, ymom_explicit_update[k], test_ye[k]);
238        }
239
240
241    }
242}
243
244
245
246void gravity_call(
247        int n, int n3, int n6, 
248        DATA_TYPE xmom_explicit_update[n], 
249        DATA_TYPE ymom_explicit_update[n], 
250
251        DATA_TYPE stage_vertex_values[n3],
252        DATA_TYPE stage_edge_values[n3],
253        DATA_TYPE stage_centroid_values[n],
254
255        DATA_TYPE bed_edge_values[n3],
256        DATA_TYPE bed_centroid_values[n],
257
258        DATA_TYPE vertex_coordinates[n6],
259
260        DATA_TYPE normals[n6],
261        DATA_TYPE areas[n],
262        DATA_TYPE edgelengths[n3],
263
264        DATA_TYPE g )
265{
266    test_call();
267    /*
268    #pragma hmpp gravity callsite
269    gravity_wb( n, n3, n6,
270        xmom_explicit_update,
271        ymom_explicit_update,
272       
273        stage_vertex_values,
274        stage_edge_values,
275        stage_centroid_values,
276       
277        bed_edge_values,
278        bed_centroid_values,
279       
280        vertex_coordinates,
281       
282        normals,
283        areas,
284        edgelengths,
285       
286        g);
287    */
288}
289
290
291
292#ifdef ISOLATING_TEST
293int main( int argc, char* argv[] ){
294    int n;   /* vector length */
295    DATA_TYPE g = 9.8;
296    int i;
297   
298    DATA_TYPE * xe, //xmom_explicit_update,
299          * ye; //ymom_explicit_update;
300    DATA_TYPE * test_xe, //xmom_explicit_update,
301          * test_ye; //ymom_explicit_update;
302
303
304    DATA_TYPE * sv, //stage_vertex_values,
305          * se, //stage_edge_values,
306          * sc; //stage_centroid_values;
307
308    DATA_TYPE * be, //bed_edge_values,
309          * bc; //bed_centroid_values;
310
311    DATA_TYPE * vc, //vertex_coordinates,
312          * normals, 
313          * a, //areas,
314          * el; //edgelengths;
315
316    char file_name[] = "merimbula.dat";
317    freopen(file_name, "r", stdin);
318
319    scanf("%d\n %lf\n", &n, &g);
320    printf("\n\nThe number of elements is %d\n", n);
321   
322    xe =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE));
323    //assert( xe != NULL);
324    ye =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE));
325    //assert( ye != NULL);
326    test_xe =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE));
327    //assert( test_xe != NULL);
328    test_ye =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE));
329    //assert( test_ye != NULL);
330
331
332    sv =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE) *3);
333    //assert( sv != NULL);
334    se =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE) *3);
335    //assert( se != NULL);
336    sc =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE) );
337    //assert( sc != NULL);
338   
339    be =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE) *3);
340    //assert( be != NULL);
341    bc =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE) );
342    //assert( bc != NULL);
343
344    vc =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE) *6);
345    //assert( vc != NULL);
346    normals=(DATA_TYPE*)malloc( n*sizeof(DATA_TYPE) *6);
347    //assert( normals != NULL);
348    a =     (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE) );
349    //assert( a != NULL);
350    el =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE) *3);
351    //assert( el != NULL);
352
353
354#ifdef USING_DOUBLE
355    for( i = 0; i < n; ++i ){
356        scanf("%lf %lf %lf\n", sv+i*3, sv+i*3 +1, sv+i*3 +2);
357        scanf("%lf %lf %lf\n", se+i*3, se+i*3 +1, se+i*3 +2);
358        scanf("%lf\n", sc+i);
359
360        scanf("%lf %lf %lf\n",be+i*3, be+i*3 +1, be+i*3 +2);
361        scanf("%lf\n", bc+i);
362
363        scanf("%lf %lf\n", vc+i*6, vc+i*6 +1);
364        scanf("%lf %lf\n", vc+i*6+2, vc+i*6 +3);
365        scanf("%lf %lf\n", vc+i*6+4, vc+i*6 +5);
366
367        scanf("%lf\n", xe+i);
368        scanf("%lf\n", ye+i);
369
370        scanf("%lf %lf %lf %lf %lf %lf\n", normals+i*6, normals+i*6+1, normals+i*6+2,normals+i*6+3, normals+i*6+4, normals+i*6+5);
371
372        scanf("%lf\n", a+i);
373
374        scanf("%lf %lf %lf\n", el+i*3, el+i*3 +1, el+i*3 +2);
375    }
376#else
377    for( i = 0; i < n; ++i ){
378        scanf("%f %f %f\n", sv+i*3, sv+i*3 +1, sv+i*3 +2);
379        scanf("%f %f %f\n", se+i*3, se+i*3 +1, se+i*3 +2);
380        scanf("%f\n", sc+i);
381
382        scanf("%f %f %f\n",be+i*3, be+i*3 +1, be+i*3 +2);
383        scanf("%f\n", bc+i);
384
385        scanf("%f %f\n", vc+i*6, vc+i*6 +1);
386        scanf("%f %f\n", vc+i*6+2, vc+i*6 +3);
387        scanf("%f %f\n", vc+i*6+4, vc+i*6 +5);
388
389        scanf("%f\n", xe+i);
390        scanf("%f\n", ye+i);
391
392        scanf("%f %f %f %f %f %f\n", normals+i*6, normals+i*6+1, normals+i*6+2,normals+i*6+3, normals+i*6+4, normals+i*6+5);
393
394        scanf("%f\n", a+i);
395
396        scanf("%f %f %f\n", el+i*3, el+i*3 +1, el+i*3 +2);
397    }
398#endif
399
400#ifndef ON_XE
401    memcpy(test_xe, xe, n*sizeof(DATA_TYPE));
402    memcpy(test_ye, ye, n*sizeof(DATA_TYPE));
403#else
404    for (i=0; i < n; i++)
405    {
406        test_xe[i] = xe[i];
407        test_ye[i] = ye[i];
408    }
409#endif
410    errs_x = 0;
411    errs_y = 0;
412    for( i = 0; i < n; ++i ){
413        if (test_xe[i] != xe[i])
414            errs_x +=1;
415        if (test_ye[i] != ye[i])
416            errs_y +=1;
417    }
418    printf("Verifying initial input: %d %d \n", errs_x, errs_y);
419
420
421    /* compute on the GPU */
422    printf(" --> Entering GPU .... \n");
423   
424    #pragma hmpp gravity callsite
425    gravity_wb( n, n*3, n*6, xe, ye, sv, se, sc, be, bc, vc, normals, a, el, g);
426   
427    /* compare results */
428    printf(" --> Entering CPU .... \n");
429    errs_x = 0;
430    errs_y = 0;
431    gravity_wb_orig( xe, ye, sv, se, sc, be, bc, vc, normals, a, el, 
432            test_xe, test_ye, n, g);
433    //gravity_wb( n, n*3, n*6, test_xe, test_ye, sv, se, sc, be, bc, vc, normals, a, el, g);
434
435   
436   
437    printf( "%d  %derrors found\n", errs_x, errs_y );
438    return errs_x + errs_y;
439}
440#endif
441
442
443
444#ifndef CALL_TEST
445void test_call(){
446    int n;   /* vector length */
447    DATA_TYPE g = 9.8;
448    int i;
449    int errs_x =0, errs_y =0;
450   
451    DATA_TYPE * xe, //xmom_explicit_update,
452          * ye; //ymom_explicit_update;
453    DATA_TYPE * test_xe, //xmom_explicit_update,
454          * test_ye; //ymom_explicit_update;
455
456
457    DATA_TYPE * sv, //stage_vertex_values,
458          * se, //stage_edge_values,
459          * sc; //stage_centroid_values;
460
461    DATA_TYPE * be, //bed_edge_values,
462          * bc; //bed_centroid_values;
463
464    DATA_TYPE * vc, //vertex_coordinates,
465          * normals, 
466          * a, //areas,
467          * el; //edgelengths;
468
469    char file_name[] = "merimbula.dat";
470    freopen(file_name, "r", stdin);
471
472    scanf("%d\n %lf\n", &n, &g);
473    printf("\n\nThe number of elements is %d\n", n);
474   
475    xe =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE));
476    //assert( xe != NULL);
477    ye =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE));
478    //assert( ye != NULL);
479    test_xe =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE));
480    //assert( test_xe != NULL);
481    test_ye =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE));
482    //assert( test_ye != NULL);
483
484
485    sv =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE) *3);
486    //assert( sv != NULL);
487    se =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE) *3);
488    //assert( se != NULL);
489    sc =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE) );
490    //assert( sc != NULL);
491   
492    be =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE) *3);
493    //assert( be != NULL);
494    bc =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE) );
495    //assert( bc != NULL);
496
497    vc =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE) *6);
498    //assert( vc != NULL);
499    normals=(DATA_TYPE*)malloc( n*sizeof(DATA_TYPE) *6);
500    //assert( normals != NULL);
501    a =     (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE) );
502    //assert( a != NULL);
503    el =    (DATA_TYPE*)malloc( n*sizeof(DATA_TYPE) *3);
504    //assert( el != NULL);
505
506
507#ifdef USING_DOUBLE
508    for( i = 0; i < n; ++i ){
509        scanf("%lf %lf %lf\n", sv+i*3, sv+i*3 +1, sv+i*3 +2);
510        scanf("%lf %lf %lf\n", se+i*3, se+i*3 +1, se+i*3 +2);
511        scanf("%lf\n", sc+i);
512
513        scanf("%lf %lf %lf\n",be+i*3, be+i*3 +1, be+i*3 +2);
514        scanf("%lf\n", bc+i);
515
516        scanf("%lf %lf\n", vc+i*6, vc+i*6 +1);
517        scanf("%lf %lf\n", vc+i*6+2, vc+i*6 +3);
518        scanf("%lf %lf\n", vc+i*6+4, vc+i*6 +5);
519
520        scanf("%lf\n", xe+i);
521        scanf("%lf\n", ye+i);
522
523        scanf("%lf %lf %lf %lf %lf %lf\n", normals+i*6, normals+i*6+1, normals+i*6+2,normals+i*6+3, normals+i*6+4, normals+i*6+5);
524
525        scanf("%lf\n", a+i);
526
527        scanf("%lf %lf %lf\n", el+i*3, el+i*3 +1, el+i*3 +2);
528    }
529#else
530    for( i = 0; i < n; ++i ){
531        scanf("%f %f %f\n", sv+i*3, sv+i*3 +1, sv+i*3 +2);
532        scanf("%f %f %f\n", se+i*3, se+i*3 +1, se+i*3 +2);
533        scanf("%f\n", sc+i);
534
535        scanf("%f %f %f\n",be+i*3, be+i*3 +1, be+i*3 +2);
536        scanf("%f\n", bc+i);
537
538        scanf("%f %f\n", vc+i*6, vc+i*6 +1);
539        scanf("%f %f\n", vc+i*6+2, vc+i*6 +3);
540        scanf("%f %f\n", vc+i*6+4, vc+i*6 +5);
541
542        scanf("%f\n", xe+i);
543        scanf("%f\n", ye+i);
544
545        scanf("%f %f %f %f %f %f\n", normals+i*6, normals+i*6+1, normals+i*6+2,normals+i*6+3, normals+i*6+4, normals+i*6+5);
546
547        scanf("%f\n", a+i);
548
549        scanf("%f %f %f\n", el+i*3, el+i*3 +1, el+i*3 +2);
550    }
551#endif
552
553#ifndef ON_XE
554    memcpy(test_xe, xe, n*sizeof(DATA_TYPE));
555    memcpy(test_ye, ye, n*sizeof(DATA_TYPE));
556#else
557    for (i=0; i < n; i++)
558    {
559        test_xe[i] = xe[i];
560        test_ye[i] = ye[i];
561    }
562#endif
563    for( i = 0; i < n; ++i ){
564        if (test_xe[i] != xe[i])
565            errs_x +=1;
566        if (test_ye[i] != ye[i])
567            errs_y +=1;
568    }
569    printf("Verifying initial input: %d %d \n", errs_x, errs_y);
570
571
572    /* compute on the GPU */
573    printf(" --> Entering GPU .... \n");
574   
575    #pragma hmpp gravity callsite
576    gravity_wb( n, n*3, n*6, xe, ye, sv, se, sc, be, bc, vc, normals, a, el, g);
577   
578    /* compare results */
579    printf(" --> Entering CPU .... \n");
580    errs_x = 0;
581    errs_y = 0;
582    gravity_wb_orig( xe, ye, sv, se, sc, be, bc, vc, normals, a, el, 
583            test_xe, test_ye, n, g);
584    //gravity_wb( n, n*3, n*6, test_xe, test_ye, sv, se, sc, be, bc, vc, normals, a, el, g);
585
586   
587   
588    printf( "%d  %derrors found\n", errs_x, errs_y );
589}
590#else
591#pragma hmpp test_call codelet, target=CUDA args[*].transfer=atcall
592void double_plus(int n, double a[n], double b[n], double c[n])
593{   
594    int i;
595    for (i=0; i<n; i++)
596        c[i] = (a[n] + b[n]) * 2;
597}
598
599
600void test_call()
601{
602    double * a, *b, *c, *c_test;
603    int n = 32, i;
604
605    a = (double *) malloc(n * sizeof(double));
606    b = (double *) malloc(n * sizeof(double));
607    c = (double *) malloc(n * sizeof(double));
608    c_test = (double *) malloc(n * sizeof(double));
609
610    #pragma hmpp test_call callsite
611    double_plus(n, a, b, c);
612
613    double_plus(n, a, b, c_test);
614
615    for (i=0; i < n ; i++)
616    {
617        if (c[i] != c_test[i])
618            printf("  -> %d, %lf %lf %lf %lf\n", i, a[i], b[i], c[i], c_test[i]);
619    }
620
621}
622#endif
Note: See TracBrowser for help on using the repository browser.