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