1#include "applu_share.h" 2#include "applu_macros.h" 3 4void erhs(){ 5 6// compute the right hand side based on exact solution 7 8//c--------------------------------------------------------------------- 9//c local variables 10//c--------------------------------------------------------------------- 11 int i, j, k, m; 12 int iglob, jglob; 13 int iex; 14 int L1, L2; 15 int ist1, iend1; 16 int jst1, jend1; 17 double dsspm; 18 double xi, eta, zeta; 19 double q; 20 double u21, u31, u41; 21 double tmp; 22 double u21i, u31i, u41i, u51i; 23 double u21j, u31j, u41j, u51j; 24 double u21k, u31k, u41k, u51k; 25 double u21im1, u31im1, u41im1, u51im1; 26 double u21jm1, u31jm1, u41jm1, u51jm1; 27 double u21km1, u31km1, u41km1, u51km1; 28 29 dsspm = dssp; 30 31 for (k=1; k<=nz; k++) 32 for (j=1; j<=ny; j++) 33 for (i=1; i<=nx; i++) 34 for (m=1; m<=5; m++) 35 frct( m, i, j, k ) = 0.0; 36 37 for (k=1; k<=nz; k++) { 38 zeta = ( (double)(k-1) ) / ( nz - 1 ); 39 for (j=1; j<=ny; j++) { 40 jglob = jpt + j; 41 eta = ( (double)(jglob-1) ) / ( ny0 - 1 ); 42 for (i=1; i<=nx; i++) { 43 iglob = ipt + i; 44 xi = ( (double)(iglob-1) ) / ( nx0 - 1 ); 45 for (m=1; m<=5; m++) { 46 rsd(m,i,j,k) = ce(m,1) 47 + ce(m,2) * xi 48 + ce(m,3) * eta 49 + ce(m,4) * zeta 50 + ce(m,5) * xi * xi 51 + ce(m,6) * eta * eta 52 + ce(m,7) * zeta * zeta 53 + ce(m,8) * xi * xi * xi 54 + ce(m,9) * eta * eta * eta 55 + ce(m,10) * zeta * zeta * zeta 56 + ce(m,11) * xi * xi * xi * xi 57 + ce(m,12) * eta * eta * eta * eta 58 + ce(m,13) * zeta * zeta * zeta * zeta; 59 } 60 } 61 } 62 } 63 64//c--------------------------------------------------------------------- 65//c xi-direction flux differences 66//c--------------------------------------------------------------------- 67//c 68//c iex = flag : iex = 0 north/south communication 69//c : iex = 1 east/west communication 70//c 71//c--------------------------------------------------------------------- 72 iex = 0; 73 74//c--------------------------------------------------------------------- 75//c communicate and receive/send two rows of data 76//c--------------------------------------------------------------------- 77// 78 exchange_3(rsd,iex); 79 80 L1 = 0; 81 if (north == -1) L1 = 1; 82 L2 = nx + 1; 83 if (south == -1) L2 = nx; 84 85 ist1 = 1; 86 iend1 = nx; 87 if (north == -1) ist1 = 4; 88 if (south == -1) iend1 = nx - 3; 89 90 for (k=2; k<=nz-1; k++) 91 for (j=jst; j<=jend; j++) 92 for (i=L1; i<=L2; i++) { 93 flux(1,i,j,k) = rsd(2,i,j,k); 94 u21 = rsd(2,i,j,k) / rsd(1,i,j,k); 95 q = 0.50 * ( rsd(2,i,j,k) * rsd(2,i,j,k) 96 + rsd(3,i,j,k) * rsd(3,i,j,k) 97 + rsd(4,i,j,k) * rsd(4,i,j,k) ) 98 / rsd(1,i,j,k); 99 flux(2,i,j,k) = rsd(2,i,j,k) * u21 + c2 * 100 ( rsd(5,i,j,k) - q ); 101 flux(3,i,j,k) = rsd(3,i,j,k) * u21; 102 flux(4,i,j,k) = rsd(4,i,j,k) * u21; 103 flux(5,i,j,k) = ( c1 * rsd(5,i,j,k) - c2 * q ) * u21; 104 } 105 106 for (k=2; k<=nz-1; k++) 107 for (j=jst; j<=jend; j++) { 108 for (i=ist; i<=iend; i++) 109 for (m=1; m<=5; m++) 110 frct(m,i,j,k) = frct(m,i,j,k) 111 - tx2 * ( flux(m,i+1,j,k) - flux(m,i-1,j,k) ); 112 for (i=ist; i<=L2; i++) { 113 tmp = 1.0 / rsd(1,i,j,k); 114 115 u21i = tmp * rsd(2,i,j,k); 116 u31i = tmp * rsd(3,i,j,k); 117 u41i = tmp * rsd(4,i,j,k); 118 u51i = tmp * rsd(5,i,j,k); 119 120 tmp = 1.0 / rsd(1,i-1,j,k); 121 122 u21im1 = tmp * rsd(2,i-1,j,k); 123 u31im1 = tmp * rsd(3,i-1,j,k); 124 u41im1 = tmp * rsd(4,i-1,j,k); 125 u51im1 = tmp * rsd(5,i-1,j,k); 126 127 flux(2,i,j,k) = (4.0/3.0) * tx3 * 128 ( u21i - u21im1 ); 129 flux(3,i,j,k) = tx3 * ( u31i - u31im1 ); 130 flux(4,i,j,k) = tx3 * ( u41i - u41im1 ); 131 flux(5,i,j,k) = 0.50 * ( 1.0 - c1*c5 ) 132 * tx3 * ( ( u21i*u21i + u31i*u31i + u41i*u41i ) 133 - ( u21im1*u21im1 + u31im1*u31im1 + u41im1*u41im1 ) ) 134 + (1.0/6.0) 135 * tx3 * ( u21i*u21i - u21im1*u21im1 ) 136 + c1 * c5 * tx3 * ( u51i - u51im1 ); 137 } 138 139 for (i=ist; i<=iend; i++) { 140 frct(1,i,j,k) = frct(1,i,j,k) 141 + dx1 * tx1 * ( rsd(1,i-1,j,k) 142 - 2.0 * rsd(1,i,j,k) 143 + rsd(1,i+1,j,k) ); 144 frct(2,i,j,k) = frct(2,i,j,k) 145 + tx3 * c3 * c4 * ( flux(2,i+1,j,k) - flux(2,i,j,k) ) 146 + dx2 * tx1 * ( rsd(2,i-1,j,k) 147 - 2.0 * rsd(2,i,j,k) 148 + rsd(2,i+1,j,k) ); 149 frct(3,i,j,k) = frct(3,i,j,k) 150 + tx3 * c3 * c4 * ( flux(3,i+1,j,k) - flux(3,i,j,k) ) 151 + dx3 * tx1 * ( rsd(3,i-1,j,k) 152 - 2.0 * rsd(3,i,j,k) 153 + rsd(3,i+1,j,k) ); 154 frct(4,i,j,k) = frct(4,i,j,k) 155 + tx3 * c3 * c4 * ( flux(4,i+1,j,k) - flux(4,i,j,k) ) 156 + dx4 * tx1 * ( rsd(4,i-1,j,k) 157 - 2.0 * rsd(4,i,j,k) 158 + rsd(4,i+1,j,k) ); 159 frct(5,i,j,k) = frct(5,i,j,k) 160 + tx3 * c3 * c4 * ( flux(5,i+1,j,k) - flux(5,i,j,k) ) 161 + dx5 * tx1 * ( rsd(5,i-1,j,k) 162 - 2.0 * rsd(5,i,j,k) 163 + rsd(5,i+1,j,k) ); 164 } 165 166//c--------------------------------------------------------------------- 167//c Fourth-order dissipation 168//c--------------------------------------------------------------------- 169 if (north == -1) { 170 for (m=1; m<=5; m++) { 171 frct(m,2,j,k) = frct(m,2,j,k) 172 - dsspm * ( + 5.0 * rsd(m,2,j,k) 173 - 4.0 * rsd(m,3,j,k) 174 + rsd(m,4,j,k) ); 175 frct(m,3,j,k) = frct(m,3,j,k) 176 - dsspm * ( - 4.0 * rsd(m,2,j,k) 177 + 6.0 * rsd(m,3,j,k) 178 - 4.0 * rsd(m,4,j,k) 179 + rsd(m,5,j,k) ); 180 } 181 } 182 183 for (i=ist1; i<=iend1; i++) 184 for (m=1; m<=5; m++) 185 frct(m,i,j,k) = frct(m,i,j,k) 186 - dsspm * ( rsd(m,i-2,j,k) 187 - 4.0 * rsd(m,i-1,j,k) 188 + 6.0 * rsd(m,i,j,k) 189 - 4.0 * rsd(m,i+1,j,k) 190 + rsd(m,i+2,j,k) ); 191 192 if (south == -1) { 193 for (m=1; m<=5; m++) { 194 frct(m,nx-2,j,k) = frct(m,nx-2,j,k) 195 - dsspm * ( rsd(m,nx-4,j,k) 196 - 4.0 * rsd(m,nx-3,j,k) 197 + 6.0 * rsd(m,nx-2,j,k) 198 - 4.0 * rsd(m,nx-1,j,k) ); 199 frct(m,nx-1,j,k) = frct(m,nx-1,j,k) 200 - dsspm * ( rsd(m,nx-3,j,k) 201 - 4.0 * rsd(m,nx-2,j,k) 202 + 5.0 * rsd(m,nx-1,j,k) ); 203 } 204 } 205 206 } 207 208//c--------------------------------------------------------------------- 209//c eta-direction flux differences 210//c--------------------------------------------------------------------- 211//c 212//c iex = flag : iex = 0 north/south communication 213//c : iex = 1 east/west communication 214//c 215//c--------------------------------------------------------------------- 216 iex = 1; 217 218//c--------------------------------------------------------------------- 219//c communicate and receive/send two rows of data 220//c--------------------------------------------------------------------- 221 222 exchange_3(rsd,iex); 223 224 L1 = 0; 225 if (west == -1) L1 = 1; 226 L2 = ny + 1; 227 if (east == -1) L2 = ny; 228 229 jst1 = 1; 230 jend1 = ny; 231 if (west == -1) jst1 = 4; 232 if (east == -1) jend1 = ny - 3; 233 234 for (k=2; k<=nz-1; k++) 235 for (j=L1; j<=L2; j++) 236 for (i=ist; i<=iend; i++) { 237 flux(1,i,j,k) = rsd(3,i,j,k); 238 u31 = rsd(3,i,j,k) / rsd(1,i,j,k); 239 q = 0.50 * ( rsd(2,i,j,k) * rsd(2,i,j,k) 240 + rsd(3,i,j,k) * rsd(3,i,j,k) 241 + rsd(4,i,j,k) * rsd(4,i,j,k) ) 242 / rsd(1,i,j,k); 243 flux(2,i,j,k) = rsd(2,i,j,k) * u31 ; 244 flux(3,i,j,k) = rsd(3,i,j,k) * u31 + c2 * 245 ( rsd(5,i,j,k) - q ); 246 flux(4,i,j,k) = rsd(4,i,j,k) * u31; 247 flux(5,i,j,k) = ( c1 * rsd(5,i,j,k) - c2 * q ) * u31; 248 } 249 250 for (k=2; k<=nz-1; k++) { 251 for (i=ist; i<=iend; i++) 252 for (j=jst; j<=jend; j++) 253 for (m=1; m<=5; m++) 254 frct(m,i,j,k) = frct(m,i,j,k) 255 - ty2 * ( flux(m,i,j+1,k) - flux(m,i,j-1,k) ); 256 257 for (j=jst; j<=L2; j++) 258 for (i=ist; i<=iend; i++) { 259 tmp = 1.0 / rsd(1,i,j,k); 260 261 u21j = tmp * rsd(2,i,j,k); 262 u31j = tmp * rsd(3,i,j,k); 263 u41j = tmp * rsd(4,i,j,k); 264 u51j = tmp * rsd(5,i,j,k); 265 266 tmp = 1.0 / rsd(1,i,j-1,k); 267 268 u21jm1 = tmp * rsd(2,i,j-1,k); 269 u31jm1 = tmp * rsd(3,i,j-1,k); 270 u41jm1 = tmp * rsd(4,i,j-1,k); 271 u51jm1 = tmp * rsd(5,i,j-1,k); 272 273 flux(2,i,j,k) = ty3 * ( u21j - u21jm1 ); 274 flux(3,i,j,k) = (4.0/3.0) * ty3 * 275 ( u31j - u31jm1 ); 276 flux(4,i,j,k) = ty3 * ( u41j - u41jm1 ); 277 flux(5,i,j,k) = 0.50 * ( 1.0 - c1*c5 ) 278 * ty3 * ( ( u21j*u21j + u31j*u31j + u41j*u41j ) 279 - ( u21jm1*u21jm1 + u31jm1*u31jm1 + u41jm1*u41jm1 ) ) 280 + (1.0/6.0) 281 * ty3 * ( u31j*u31j - u31jm1*u31jm1 ) 282 + c1 * c5 * ty3 * ( u51j - u51jm1 ); 283 } 284 285 for (j=jst; j<=jend; j++) 286 for (i=ist; i<=iend; i++) { 287 frct(1,i,j,k) = frct(1,i,j,k) 288 + dy1 * ty1 * ( rsd(1,i,j-1,k) 289 - 2.0 * rsd(1,i,j,k) 290 + rsd(1,i,j+1,k) ); 291 frct(2,i,j,k) = frct(2,i,j,k) 292 + ty3 * c3 * c4 * ( flux(2,i,j+1,k) - flux(2,i,j,k) ) 293 + dy2 * ty1 * ( rsd(2,i,j-1,k) 294 - 2.0 * rsd(2,i,j,k) 295 + rsd(2,i,j+1,k) ); 296 frct(3,i,j,k) = frct(3,i,j,k) 297 + ty3 * c3 * c4 * ( flux(3,i,j+1,k) - flux(3,i,j,k) ) 298 + dy3 * ty1 * ( rsd(3,i,j-1,k) 299 - 2.0 * rsd(3,i,j,k) 300 + rsd(3,i,j+1,k) ); 301 frct(4,i,j,k) = frct(4,i,j,k) 302 + ty3 * c3 * c4 * ( flux(4,i,j+1,k) - flux(4,i,j,k) ) 303 + dy4 * ty1 * ( rsd(4,i,j-1,k) 304 - 2.0 * rsd(4,i,j,k) 305 + rsd(4,i,j+1,k) ); 306 frct(5,i,j,k) = frct(5,i,j,k) 307 + ty3 * c3 * c4 * ( flux(5,i,j+1,k) - flux(5,i,j,k) ) 308 + dy5 * ty1 * ( rsd(5,i,j-1,k) 309 - 2.0 * rsd(5,i,j,k) 310 + rsd(5,i,j+1,k) ); 311 } 312 313//c--------------------------------------------------------------------- 314//c fourth-order dissipation 315//c--------------------------------------------------------------------- 316 if (west == -1) { 317 for (i=ist; i<=iend; i++) 318 for (m=1; m<=5; m++) { 319 frct(m,i,2,k) = frct(m,i,2,k) 320 - dsspm * ( + 5.0 * rsd(m,i,2,k) 321 - 4.0 * rsd(m,i,3,k) 322 + rsd(m,i,4,k) ); 323 frct(m,i,3,k) = frct(m,i,3,k) 324 - dsspm * ( - 4.0 * rsd(m,i,2,k) 325 + 6.0 * rsd(m,i,3,k) 326 - 4.0 * rsd(m,i,4,k) 327 + rsd(m,i,5,k) ); 328 } 329 } 330 331 for (j=jst1; j<=jend1; j++) 332 for (i=ist; i<=iend; i++) 333 for (m=1; m<=5; m++) 334 frct(m,i,j,k) = frct(m,i,j,k) 335 - dsspm * ( rsd(m,i,j-2,k) 336 - 4.0 * rsd(m,i,j-1,k) 337 + 6.0 * rsd(m,i,j,k) 338 - 4.0 * rsd(m,i,j+1,k) 339 + rsd(m,i,j+2,k) ); 340 341 if (east == -1) { 342 for (i=ist; i<=iend; i++) 343 for (m=1; m<=5; m++) { 344 frct(m,i,ny-2,k) = frct(m,i,ny-2,k) 345 - dsspm * ( rsd(m,i,ny-4,k) 346 - 4.0 * rsd(m,i,ny-3,k) 347 + 6.0 * rsd(m,i,ny-2,k) 348 - 4.0 * rsd(m,i,ny-1,k) ); 349 frct(m,i,ny-1,k) = frct(m,i,ny-1,k) 350 - dsspm * ( rsd(m,i,ny-3,k) 351 - 4.0 * rsd(m,i,ny-2,k) 352 + 5.0 * rsd(m,i,ny-1,k) ); 353 } 354 } 355 356 } 357 358//c--------------------------------------------------------------------- 359//c zeta-direction flux differences 360//c--------------------------------------------------------------------- 361 for (k=1; k<=nz; k++) 362 for (j=jst; j<=jend; j++) 363 for (i=ist; i<=iend; i++) { 364 flux(1,i,j,k) = rsd(4,i,j,k); 365 u41 = rsd(4,i,j,k) / rsd(1,i,j,k); 366 q = 0.50 * ( rsd(2,i,j,k) * rsd(2,i,j,k) 367 + rsd(3,i,j,k) * rsd(3,i,j,k) 368 + rsd(4,i,j,k) * rsd(4,i,j,k) ) 369 / rsd(1,i,j,k); 370 flux(2,i,j,k) = rsd(2,i,j,k) * u41 ; 371 flux(3,i,j,k) = rsd(3,i,j,k) * u41 ; 372 flux(4,i,j,k) = rsd(4,i,j,k) * u41 + c2 * 373 ( rsd(5,i,j,k) - q ); 374 flux(5,i,j,k) = ( c1 * rsd(5,i,j,k) - c2 * q ) * u41; 375 } 376 377 for (k=2; k<=nz-1; k++) 378 for (j=jst; j<=jend; j++) 379 for (i=ist; i<=iend; i++) 380 for (m=1; m<=5; m++) 381 frct(m,i,j,k) = frct(m,i,j,k) 382 - tz2 * ( flux(m,i,j,k+1) - flux(m,i,j,k-1) ); 383 384 for (k=2; k<=nz; k++) 385 for (j=jst; j<=jend; j++) 386 for (i=ist; i<=iend; i++) { 387 tmp = 1.0 / rsd(1,i,j,k); 388 389 u21k = tmp * rsd(2,i,j,k); 390 u31k = tmp * rsd(3,i,j,k); 391 u41k = tmp * rsd(4,i,j,k); 392 u51k = tmp * rsd(5,i,j,k); 393 394 tmp = 1.0 / rsd(1,i,j,k-1); 395 396 u21km1 = tmp * rsd(2,i,j,k-1); 397 u31km1 = tmp * rsd(3,i,j,k-1); 398 u41km1 = tmp * rsd(4,i,j,k-1); 399 u51km1 = tmp * rsd(5,i,j,k-1); 400 401 flux(2,i,j,k) = tz3 * ( u21k - u21km1 ); 402 flux(3,i,j,k) = tz3 * ( u31k - u31km1 ); 403 flux(4,i,j,k) = (4.0/3.0) * tz3 * ( u41k 404 - u41km1 ); 405 flux(5,i,j,k) = 0.50 * ( 1.0 - c1*c5 ) 406 * tz3 * ( ( u21k*u21k + u31k*u31k + u41k*u41k ) 407 - ( u21km1*u21km1 + u31km1*u31km1 + u41km1*u41km1 ) ) 408 + (1.0/6.0) 409 * tz3 * ( u41k*u41k - u41km1*u41km1 ) 410 + c1 * c5 * tz3 * ( u51k - u51km1 ); 411 } 412 413 for (k=2; k<=nz-1; k++) 414 for (j=jst; j<=jend; j++) 415 for (i=ist; i<=iend; i++) { 416 frct(1,i,j,k) = frct(1,i,j,k) 417 + dz1 * tz1 * ( rsd(1,i,j,k+1) 418 - 2.0 * rsd(1,i,j,k) 419 + rsd(1,i,j,k-1) ); 420 frct(2,i,j,k) = frct(2,i,j,k) 421 + tz3 * c3 * c4 * ( flux(2,i,j,k+1) - flux(2,i,j,k) ) 422 + dz2 * tz1 * ( rsd(2,i,j,k+1) 423 - 2.0 * rsd(2,i,j,k) 424 + rsd(2,i,j,k-1) ); 425 frct(3,i,j,k) = frct(3,i,j,k) 426 + tz3 * c3 * c4 * ( flux(3,i,j,k+1) - flux(3,i,j,k) ) 427 + dz3 * tz1 * ( rsd(3,i,j,k+1) 428 - 2.0 * rsd(3,i,j,k) 429 + rsd(3,i,j,k-1) ); 430 frct(4,i,j,k) = frct(4,i,j,k) 431 + tz3 * c3 * c4 * ( flux(4,i,j,k+1) - flux(4,i,j,k) ) 432 + dz4 * tz1 * ( rsd(4,i,j,k+1) 433 - 2.0 * rsd(4,i,j,k) 434 + rsd(4,i,j,k-1) ); 435 frct(5,i,j,k) = frct(5,i,j,k) 436 + tz3 * c3 * c4 * ( flux(5,i,j,k+1) - flux(5,i,j,k) ) 437 + dz5 * tz1 * ( rsd(5,i,j,k+1) 438 - 2.0 * rsd(5,i,j,k) 439 + rsd(5,i,j,k-1) ); 440 } 441 442//c--------------------------------------------------------------------- 443//c fourth-order dissipation 444//c--------------------------------------------------------------------- 445 for (j=jst; j<=jend; j++) 446 for (i=ist; i<=iend; i++) 447 for (m=1; m<=5; m++) { 448 frct(m,i,j,2) = frct(m,i,j,2) 449 - dsspm * ( + 5.0 * rsd(m,i,j,2) 450 - 4.0 * rsd(m,i,j,3) 451 + rsd(m,i,j,4) ); 452 frct(m,i,j,3) = frct(m,i,j,3) 453 - dsspm * (- 4.0 * rsd(m,i,j,2) 454 + 6.0 * rsd(m,i,j,3) 455 - 4.0 * rsd(m,i,j,4) 456 + rsd(m,i,j,5) ); 457 } 458 459 for (k=4; k<=nz-3; k++) 460 for (j=jst; j<=jend; j++) 461 for (i=ist; i<=iend; i++) 462 for (m=1; m<=5; m++) 463 frct(m,i,j,k) = frct(m,i,j,k) 464 - dsspm * ( rsd(m,i,j,k-2) 465 - 4.0 * rsd(m,i,j,k-1) 466 + 6.0 * rsd(m,i,j,k) 467 - 4.0 * rsd(m,i,j,k+1) 468 + rsd(m,i,j,k+2) ); 469 470 for (j=jst; j<=jend; j++) 471 for (i=ist; i<=iend; i++) 472 for (m=1; m<=5; m++) { 473 frct(m,i,j,nz-2) = frct(m,i,j,nz-2) 474 - dsspm * ( rsd(m,i,j,nz-4) 475 - 4.0 * rsd(m,i,j,nz-3) 476 + 6.0 * rsd(m,i,j,nz-2) 477 - 4.0 * rsd(m,i,j,nz-1) ); 478 frct(m,i,j,nz-1) = frct(m,i,j,nz-1) 479 - dsspm * ( rsd(m,i,j,nz-3) 480 - 4.0 * rsd(m,i,j,nz-2) 481 + 5.0 * rsd(m,i,j,nz-1) ); 482 } 483 484 return; 485} 486 487