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