| 1 | //#define USING_MAIN |
|---|
| 2 | |
|---|
| 3 | #ifndef USING_MAIN |
|---|
| 4 | #include "hmpp_fun.h" |
|---|
| 5 | #else |
|---|
| 6 | #include <math.h> |
|---|
| 7 | #endif |
|---|
| 8 | |
|---|
| 9 | //#ifdef USING_MAIN |
|---|
| 10 | #ifdef USING_LOCAL_DIRECTIVES |
|---|
| 11 | #pragma hmpp protectSW codelet, target=CUDA args[*].transfer=atcall |
|---|
| 12 | #endif |
|---|
| 13 | void protect_sw( |
|---|
| 14 | int N, |
|---|
| 15 | int N3, |
|---|
| 16 | double minimum_allowed_height, |
|---|
| 17 | double maximum_allowed_speed, |
|---|
| 18 | double epsilon, |
|---|
| 19 | |
|---|
| 20 | double wc[N], // stage_centroid_values |
|---|
| 21 | double zc[N], // bed_centroid_values |
|---|
| 22 | double xmomc[N],// xmom_centroid_values |
|---|
| 23 | double ymomc[N] //ymom_centroid_values |
|---|
| 24 | ) |
|---|
| 25 | { |
|---|
| 26 | int k; |
|---|
| 27 | double hc; |
|---|
| 28 | double u, v, reduced_speed; |
|---|
| 29 | |
|---|
| 30 | |
|---|
| 31 | // Protect against initesimal and negative heights |
|---|
| 32 | if (maximum_allowed_speed < epsilon) { |
|---|
| 33 | #pragma hmppcg gridify(k), & |
|---|
| 34 | #pragma hmppcg & private( hc, u, v, reduced_speed), & |
|---|
| 35 | #pragma hmppcg & global( minimum_allowed_height, maximum_allowed_speed, & |
|---|
| 36 | #pragma hmppcg & epsilon, wc, zc, xmomc, ymomc) |
|---|
| 37 | for (k = 0; k < N; k++) { |
|---|
| 38 | hc = wc[k] - zc[k]; |
|---|
| 39 | if (hc < minimum_allowed_height) { |
|---|
| 40 | // Set momentum to zero and ensure h is non negative |
|---|
| 41 | xmomc[k] = 0.0; |
|---|
| 42 | ymomc[k] = 0.0; |
|---|
| 43 | if (hc <= 0.0) wc[k] = zc[k]; |
|---|
| 44 | } |
|---|
| 45 | } |
|---|
| 46 | |
|---|
| 47 | } else { |
|---|
| 48 | // Protect against initesimal and negative heights |
|---|
| 49 | #pragma hmppcg gridify(k), & |
|---|
| 50 | #pragma hmppcg & private( hc, u, v, reduced_speed), & |
|---|
| 51 | #pragma hmppcg & global( minimum_allowed_height, maximum_allowed_speed, & |
|---|
| 52 | #pragma hmppcg & epsilon, wc, zc, xmomc, ymomc) |
|---|
| 53 | for (k = 0; k < N; k++) { |
|---|
| 54 | hc = wc[k] - zc[k]; |
|---|
| 55 | if (hc < minimum_allowed_height) { |
|---|
| 56 | //New code: Adjust momentum to guarantee speeds are physical |
|---|
| 57 | // ensure h is non negative |
|---|
| 58 | if (hc <= 0.0) { |
|---|
| 59 | wc[k] = zc[k]; |
|---|
| 60 | xmomc[k] = 0.0; |
|---|
| 61 | ymomc[k] = 0.0; |
|---|
| 62 | } else { |
|---|
| 63 | //Reduce excessive speeds derived from division by small hc |
|---|
| 64 | //FIXME (Ole): This may be unnecessary with new slope limiters |
|---|
| 65 | //in effect. |
|---|
| 66 | |
|---|
| 67 | u = xmomc[k] / hc; |
|---|
| 68 | if (fabs(u) > maximum_allowed_speed) { |
|---|
| 69 | reduced_speed = maximum_allowed_speed * u / fabs(u); |
|---|
| 70 | xmomc[k] = reduced_speed * hc; |
|---|
| 71 | } |
|---|
| 72 | |
|---|
| 73 | v = ymomc[k] / hc; |
|---|
| 74 | if (fabs(v) > maximum_allowed_speed) { |
|---|
| 75 | reduced_speed = maximum_allowed_speed * v / fabs(v); |
|---|
| 76 | ymomc[k] = reduced_speed * hc; |
|---|
| 77 | } |
|---|
| 78 | } |
|---|
| 79 | } |
|---|
| 80 | } |
|---|
| 81 | } |
|---|
| 82 | } |
|---|
| 83 | |
|---|
| 84 | |
|---|
| 85 | |
|---|
| 86 | #ifdef USING_MAIN |
|---|
| 87 | int main() |
|---|
| 88 | #else |
|---|
| 89 | void test_protect_sw() |
|---|
| 90 | #endif |
|---|
| 91 | { |
|---|
| 92 | int N=4; |
|---|
| 93 | double minimum_allowed_height=0.001, |
|---|
| 94 | maximum_allowed_speed=0, |
|---|
| 95 | epsilon=1e-12; |
|---|
| 96 | |
|---|
| 97 | double wc[]={1,2,3,4}, |
|---|
| 98 | zc[]={1,2,3,4}, |
|---|
| 99 | xmomc[]={1,2,3,4}, |
|---|
| 100 | ymomc[]={1,2,3,4}; |
|---|
| 101 | |
|---|
| 102 | #pragma hmpp protectSW callsite |
|---|
| 103 | protect_sw(N, N*3, minimum_allowed_height, |
|---|
| 104 | maximum_allowed_speed, epsilon, wc, zc, xmomc, ymomc); |
|---|
| 105 | } |
|---|
