1#include "RCCE.h" 2#include <stdio.h> 3#include <stdlib.h> 4#include "applu_share.h" 5#include "applu_macros.h" 6 7void pintgr() { 8 9 int i, j, k, ind1, ind2; 10 int ibeg, ifin, ifin1, jbeg, jfin, jfin1; 11 int iglob, iglob1, iglob2, jglob, jglob1, jglob2; 12 double phi1[(isiz2+2)*(isiz3+2)], phi2[(isiz2+2)*(isiz3+2)], 13 frc1, frc2, frc3, dummy; 14 15//c--------------------------------------------------------------------- 16//c set up the sub-domains for integeration in each processor 17//c--------------------------------------------------------------------- 18 ibeg = nx + 1; 19 ifin = 0; 20 iglob1 = ipt + 1; 21 iglob2 = ipt + nx; 22 if ((iglob1 >= ii1) && (iglob2 < ii2+nx)) ibeg = 1; 23 if ((iglob1 > ii1-nx) && (iglob2 <= ii2) ) ifin = nx; 24 if ((ii1 >= iglob1) && (ii1 <= iglob2) ) ibeg = ii1 - ipt; 25 if ((ii2 >= iglob1) && (ii2 <= iglob2) ) ifin = ii2 - ipt; 26 jbeg = ny + 1; 27 jfin = 0; 28 jglob1 = jpt + 1; 29 jglob2 = jpt + ny; 30 if ((jglob1 >= ji1) && (jglob2 < ji2+ny)) jbeg = 1; 31 if ((jglob1 > ji1-ny) && (jglob2 <= ji2) ) jfin = ny; 32 if ((ji1 >= jglob1) && (ji1 <= jglob2) ) jbeg = ji1 - jpt; 33 if ((ji2 >= jglob1) && (ji2 <= jglob2) ) jfin = ji2 - jpt; 34 ifin1 = ifin; 35 jfin1 = jfin; 36 if (ipt + ifin1 == ii2) ifin1 = ifin -1; 37 if (jpt + jfin1 == ji2) jfin1 = jfin -1; 38 39//c--------------------------------------------------------------------- 40//c initialize 41//c--------------------------------------------------------------------- 42 for (i=0; i<=isiz2+1; i++) { 43 for (k=0; k<=isiz3+1; k++) { 44 phi1(i,k) = 0.0; 45 phi2(i,k) = 0.0; 46 } 47 } 48 49 for (j= jbeg; j<=jfin; j++) { 50 jglob = jpt + j; 51 for (i=ibeg; i<=ifin; i++) { 52 iglob = ipt + i; 53 54 k = ki1; 55 56 phi1(i,j) = c2*( u(5,i,j,k) 57 - 0.50 * ( u(2,i,j,k) * u(2,i,j,k) 58 + u(3,i,j,k) * u(3,i,j,k) 59 + u(4,i,j,k) * u(4,i,j,k) ) 60 / u(1,i,j,k) ); 61 62 k = ki2; 63 64 phi2(i,j) = c2*( u(5,i,j,k) 65 - 0.50 * ( u(2,i,j,k) * u(2,i,j,k) 66 + u(3,i,j,k) * u(3,i,j,k) 67 + u(4,i,j,k) * u(4,i,j,k) ) 68 / u(1,i,j,k) ); 69 } 70 } 71 72//c--------------------------------------------------------------------- 73//c communicate in i and j directions 74//c--------------------------------------------------------------------- 75 exchange_4(phi1, phi2, ibeg, ifin1, jbeg, jfin1); 76 77 frc1 = 0.0; 78 79 for (j= jbeg; j<=jfin1; j++) { 80 for (i=ibeg; i<= ifin1; i++) { 81 frc1 = frc1 + ( phi1(i,j) 82 + phi1(i+1,j) 83 + phi1(i,j+1) 84 + phi1(i+1,j+1) 85 + phi2(i,j) 86 + phi2(i+1,j) 87 + phi2(i,j+1) 88 + phi2(i+1,j+1) ); 89 } 90 } 91 92//c--------------------------------------------------------------------- 93//c compute the global sum of individual contributions to frc1 94//c--------------------------------------------------------------------- 95 dummy = frc1; 96 RCCE_allreduce((char*)(&dummy), (char*)(&frc1), 1, RCCE_DOUBLE, RCCE_SUM, RCCE_COMM_WORLD); 97 98 frc1 = dxi * deta * frc1; 99 100//c--------------------------------------------------------------------- 101//c initialize 102//c--------------------------------------------------------------------- 103 for (i=0; i<=isiz2+1; i++) { 104 for (k= 0; k<=isiz3+1; k++) { 105 phi1(i,k) = 0.0; 106 phi2(i,k) = 0.0; 107 } 108 } 109 jglob = jpt + jbeg; 110 ind1 = 0; 111 if (jglob == ji1) { 112 ind1 = 1; 113 for (k= ki1; k<= ki2; k++) { 114 for (i=ibeg; i<= ifin; i++) { 115 iglob = ipt + i; 116 phi1(i,k) = c2*( u(5,i,jbeg,k) 117 - 0.50 * ( u(2,i,jbeg,k) * u(2,i,jbeg,k) 118 + u(3,i,jbeg,k) * u(3,i,jbeg,k) 119 + u(4,i,jbeg,k) *u(4,i,jbeg,k) ) 120 / u(1,i,jbeg,k) ); 121 } 122 } 123 } 124 125 jglob = jpt + jfin; 126 ind2 = 0; 127 if (jglob == ji2) { 128 ind2 = 1; 129 for (k= ki1; k<= ki2; k++) { 130 for (i=ibeg; i<= ifin; i++) { 131 iglob = ipt + i; 132 phi2(i,k) = c2*( u(5,i,jfin,k) 133 - 0.50 * ( u(2,i,jfin,k) * u(2,i,jfin,k) 134 + u(3,i,jfin,k) * u(3,i,jfin,k) 135 + u(4,i,jfin,k) * u(4,i,jfin,k) ) 136 / u(1,i,jfin,k) ); 137 } 138 } 139 } 140 141//c--------------------------------------------------------------------- 142//c communicate in i direction 143//c--------------------------------------------------------------------- 144 if (ind1 == 1) exchange_5(phi1, ibeg, ifin1); 145 if (ind2 == 1) exchange_5(phi2, ibeg, ifin1); 146 147 frc2 = 0.0; 148 for (k= ki1; k<= ki2-1; k++) { 149 for (i=ibeg; i<= ifin1; i++) { 150 frc2 = frc2 + ( phi1(i,k) 151 + phi1(i+1,k) 152 + phi1(i,k+1) 153 + phi1(i+1,k+1) 154 + phi2(i,k) 155 + phi2(i+1,k) 156 + phi2(i,k+1) 157 + phi2(i+1,k+1) ); 158 } 159 } 160 161//c--------------------------------------------------------------------- 162//c compute the global sum of individual contributions to frc2 163//c--------------------------------------------------------------------- 164 dummy = frc2; 165 RCCE_allreduce((char*)(&dummy), (char*)(&frc2), 1, RCCE_DOUBLE, RCCE_SUM, RCCE_COMM_WORLD); 166 167 frc2 = dxi * dzeta * frc2; 168 169//c--------------------------------------------------------------------- 170//c initialize 171//c--------------------------------------------------------------------- 172 for (i=0; i<=isiz2+1; i++) { 173 for (k= 0; k<=isiz3+1; k++) { 174 phi1(i,k) = 0.0; 175 phi2(i,k) = 0.0; 176 } 177 } 178 iglob = ipt + ibeg; 179 ind1 = 0; 180 if (iglob == ii1) { 181 ind1 = 1; 182 for (k= ki1; k<= ki2; k++) { 183 for (j= jbeg; j<= jfin; j++) { 184 jglob = jpt + j; 185 phi1(j,k) = c2*( u(5,ibeg,j,k) 186 - 0.50 * ( u(2,ibeg,j,k) * u(2,ibeg,j,k) 187 + u(3,ibeg,j,k) * u(3,ibeg,j,k) 188 + u(4,ibeg,j,k) * u(4,ibeg,j,k) ) 189 / u(1,ibeg,j,k) ); 190 } 191 } 192 } 193 194 iglob = ipt + ifin; 195 ind2 = 0; 196 if (iglob == ii2) { 197 ind2 = 1; 198 for (k= ki1; k<= ki2; k++) { 199 for (j= jbeg; j<= jfin; j++) { 200 jglob = jpt + j; 201 phi2(j,k) = c2*( u(5,ifin,j,k) 202 - 0.50 * ( u(2,ifin,j,k) * u(2,ifin,j,k) 203 + u(3,ifin,j,k) * u(3,ifin,j,k) 204 + u(4,ifin,j,k) * u(4,ifin,j,k) ) 205 / u(1,ifin,j,k) ); 206 } 207 } 208 } 209 210//c--------------------------------------------------------------------- 211//c communicate in j direction 212//c--------------------------------------------------------------------- 213 if (ind1 == 1) exchange_6(phi1, jbeg, jfin1); 214 if (ind2 == 1) exchange_6(phi2, jbeg, jfin1); 215 216 frc3 = 0.0; 217 218 for (k= ki1; k<= ki2-1; k++) { 219 for (j= jbeg; j<= jfin1; j++) { 220 frc3 = frc3 + ( phi1(j,k) 221 + phi1(j+1,k) 222 + phi1(j,k+1) 223 + phi1(j+1,k+1) 224 + phi2(j,k) 225 + phi2(j+1,k) 226 + phi2(j,k+1) 227 + phi2(j+1,k+1) ); 228 } 229 } 230 231//c--------------------------------------------------------------------- 232//c compute the global sum of individual contributions to frc3 233//c--------------------------------------------------------------------- 234 dummy = frc3; 235 RCCE_allreduce((char*)(&dummy), (char*)(&frc3), 1, RCCE_DOUBLE, RCCE_SUM, RCCE_COMM_WORLD); 236 237 frc3 = deta * dzeta * frc3; 238 frc = 0.25 * ( frc1 + frc2 + frc3 ); 239 240 return; 241} 242