1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC 2 ! C 3 ! Module name: DISCRETIZATION C 4 ! Purpose: A collection of functions for higher-order discretization C 5 ! C 6 ! Author: M. Syamlal Date: 18-MAR-97 C 7 ! Reviewer: Date: C 8 ! C 9 ! C 10 ! Literature/Document References: C 11 ! C 12 ! Variables referenced: C 13 ! Variables modified: C 14 ! C 15 ! Local variables: C 16 ! C 17 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C 18 19 MODULE discretization 20 21 CONTAINS 22 23 DOUBLE PRECISION FUNCTION SUPERBEE (PHI_C) 24 !...Translated by Pacific-Sierra Research VAST-90 2.06G5 12:17:31 12/09/98 25 !...Switches: -xf 26 !----------------------------------------------- 27 ! M o d u l e s 28 !----------------------------------------------- 29 USE param 30 USE param1 31 32 IMPLICIT NONE 33 !----------------------------------------------- 34 ! D u m m y A r g u m e n t s 35 !----------------------------------------------- 36 DOUBLE PRECISION PHI_C 37 !----------------------------------------------- 38 ! L o c a l P a r a m e t e r s 39 !----------------------------------------------- 40 DOUBLE PRECISION, PARAMETER :: TWO = 2.D0 41 !----------------------------------------------- 42 ! L o c a l V a r i a b l e s 43 !----------------------------------------------- 44 DOUBLE PRECISION :: TH 45 !----------------------------------------------- 46 ! 47 IF (PHI_C>=ZERO .AND. PHI_C<ONE) THEN !monotonic region 48 TH = PHI_C/(ONE - PHI_C) 49 SUPERBEE = HALF*MAX(ZERO,MIN(ONE,TWO*TH),MIN(TWO,TH)) 50 ELSE IF (PHI_C == ONE) THEN 51 SUPERBEE = ONE 52 ELSE !first order upwinding 53 SUPERBEE = ZERO 54 ENDIF 55 ! 56 RETURN 57 END FUNCTION SUPERBEE 58 ! 59 ! 60 DOUBLE PRECISION FUNCTION SMART (PHI_C) 61 !...Translated by Pacific-Sierra Research VAST-90 2.06G5 12:17:31 12/09/98 62 !...Switches: -xf 63 !----------------------------------------------- 64 ! M o d u l e s 65 !----------------------------------------------- 66 USE param 67 USE param1 68 IMPLICIT NONE 69 !----------------------------------------------- 70 ! D u m m y A r g u m e n t s 71 !----------------------------------------------- 72 DOUBLE PRECISION PHI_C 73 !----------------------------------------------- 74 ! L o c a l P a r a m e t e r s 75 !----------------------------------------------- 76 !----------------------------------------------- 77 ! L o c a l V a r i a b l e s 78 !----------------------------------------------- 79 DOUBLE PRECISION :: TH 80 !----------------------------------------------- 81 ! 82 IF (PHI_C>=ZERO .AND. PHI_C<ONE) THEN !monotonic region 83 TH = PHI_C/(ONE - PHI_C) 84 SMART = HALF*MAX(ZERO,MIN(4.0D0*TH,0.75D0 + 0.25D0*TH,2.0D0)) 85 ELSE IF (PHI_C == ONE) THEN 86 SMART = ONE 87 ELSE !first order upwinding 88 SMART = ZERO 89 ENDIF 90 ! 91 RETURN 92 END FUNCTION SMART 93 94 DOUBLE PRECISION FUNCTION Chi_SMART (PHI_C, Chi) 95 ! calculate DWF from Chi_SMART scheme 96 !----------------------------------------------- 97 ! M o d u l e s 98 !----------------------------------------------- 99 USE param 100 USE param1 101 IMPLICIT NONE 102 !----------------------------------------------- 103 ! D u m m y A r g u m e n t s 104 !----------------------------------------------- 105 DOUBLE PRECISION PHI_C, Chi 106 107 if(PHI_C < ONE)then 108 Chi_SMART = Chi * (3.D0/8.D0 - PHI_C/4.d0) / (ONE-PHI_C) 109 elseif(Chi == zero)then ! insures that all species equations uses the same chi 110 Chi_SMART = zero 111 else 112 Chi_SMART = ONE 113 endif 114 115 RETURN 116 END FUNCTION Chi_SMART 117 118 DOUBLE PRECISION FUNCTION Chi4SMART (PHI_C, PHIU, PHIC, PHID) 119 ! calculate CHI for SMART scheme 120 !----------------------------------------------- 121 ! M o d u l e s 122 !----------------------------------------------- 123 USE param 124 USE param1 125 IMPLICIT NONE 126 !----------------------------------------------- 127 ! D u m m y A r g u m e n t s 128 !----------------------------------------------- 129 DOUBLE PRECISION PHI_C , PHIU, PHIC, PHID 130 131 IF(PHI_C > ZERO .AND. PHI_C <= (1.D0/6.D0))THEN 132 Chi4SMART = 16.D0 * PHI_C/(3.D0 - 2.D0 * PHI_C) 133 ELSEIF(PHI_C > (1.D0/6.D0) .AND. PHI_C <= (5.D0/6.D0))THEN 134 Chi4SMART = ONE 135 ELSEIF(PHI_C > (5.D0/6.D0) .AND. PHI_C <= ONE)THEN 136 Chi4SMART = 8.D0*(ONE-PHI_C)/(3.D0 - 2.D0 * PHI_C) 137 ELSEIF( (PHIU < 1d-15) .AND. (PHIC < 1d-15) .AND. (PHID < 1d-15) )THEN 138 ! if a species Xg is less than machine precision, do not use its chi. 139 ! this will avoid calculating chi = 0 for a vanishing species. 140 Chi4SMART = LARGE_NUMBER 141 ELSE 142 Chi4SMART = ZERO 143 ENDIF 144 145 RETURN 146 END FUNCTION Chi4SMART 147 ! 148 ! 149 ! 150 DOUBLE PRECISION FUNCTION ULTRA_QUICK (PHI_C, CF) 151 !...Translated by Pacific-Sierra Research VAST-90 2.06G5 12:17:31 12/09/98 152 !...Switches: -xf 153 !----------------------------------------------- 154 ! M o d u l e s 155 !----------------------------------------------- 156 USE param 157 USE param1 158 IMPLICIT NONE 159 !----------------------------------------------- 160 ! D u m m y A r g u m e n t s 161 !----------------------------------------------- 162 DOUBLE PRECISION PHI_C, CF 163 !----------------------------------------------- 164 ! L o c a l P a r a m e t e r s 165 !----------------------------------------------- 166 DOUBLE PRECISION, PARAMETER :: FIVEOSIX = 5.D0/6.D0 167 !----------------------------------------------- 168 ! L o c a l V a r i a b l e s 169 !----------------------------------------------- 170 DOUBLE PRECISION :: TH, OCF 171 !----------------------------------------------- 172 ! 173 OCF = MAX(ONE,ONE/MAX(1.D-2,CF)) 174 IF (PHI_C > ONE) THEN 175 ULTRA_QUICK = HALF 176 ELSE IF (PHI_C > FIVEOSIX) THEN 177 ULTRA_QUICK = ONE 178 ELSE IF (PHI_C > 3.D0/(8.D0*OCF - 6.D0)) THEN 179 TH = PHI_C/(ONE - PHI_C) 180 ULTRA_QUICK = HALF - 0.125D0*(ONE - TH) 181 ELSE IF (PHI_C > ZERO) THEN 182 TH = PHI_C/(ONE - PHI_C) 183 ULTRA_QUICK = (OCF - ONE)*TH 184 ELSE 185 TH = PHI_C/(ONE - PHI_C) 186 ULTRA_QUICK = HALF*TH 187 ENDIF 188 ! 189 RETURN 190 END FUNCTION ULTRA_QUICK 191 ! 192 ! 193 ! 194 DOUBLE PRECISION FUNCTION QUICKEST (PHI_C, CF, ODXC, ODXUC, ODXCD) 195 !...Translated by Pacific-Sierra Research VAST-90 2.06G5 12:17:31 12/09/98 196 !...Switches: -xf 197 !----------------------------------------------- 198 ! M o d u l e s 199 !----------------------------------------------- 200 USE param 201 USE param1 202 IMPLICIT NONE 203 !----------------------------------------------- 204 ! D u m m y A r g u m e n t s 205 !----------------------------------------------- 206 DOUBLE PRECISION PHI_C, CF, ODXC, ODXUC, ODXCD 207 !----------------------------------------------- 208 ! L o c a l P a r a m e t e r s 209 !----------------------------------------------- 210 !----------------------------------------------- 211 ! L o c a l V a r i a b l e s 212 !----------------------------------------------- 213 DOUBLE PRECISION :: FCF, TH 214 !----------------------------------------------- 215 ! 216 IF (PHI_C>ZERO .AND. PHI_C<ONE) THEN !monotonic region 217 FCF = -(ONE - CF*CF)/3.D0 218 TH = PHI_C/(ONE - PHI_C + SMALL_NUMBER) 219 QUICKEST = HALF*(ONE - CF) + FCF*(ODXC/ODXCD - ODXC*ODXUC*TH/ODXCD**2) 220 IF(PHI_C<CF)QUICKEST=MIN(QUICKEST,(ONE/CF-ONE)*PHI_C/(ONE-PHI_C)) 221 QUICKEST = MAX(ZERO,MIN(ONE,QUICKEST)) 222 ELSE !first order upwinding 223 QUICKEST = ZERO 224 ENDIF 225 ! 226 RETURN 227 END FUNCTION QUICKEST 228 ! 229 ! 230 DOUBLE PRECISION FUNCTION MUSCL (PHI_C) 231 !...Translated by Pacific-Sierra Research VAST-90 2.06G5 12:17:31 12/09/98 232 !...Switches: -xf 233 !----------------------------------------------- 234 ! M o d u l e s 235 !----------------------------------------------- 236 USE param 237 USE param1 238 IMPLICIT NONE 239 !----------------------------------------------- 240 ! D u m m y A r g u m e n t s 241 !----------------------------------------------- 242 DOUBLE PRECISION PHI_C 243 !----------------------------------------------- 244 ! L o c a l P a r a m e t e r s 245 !----------------------------------------------- 246 !----------------------------------------------- 247 ! L o c a l V a r i a b l e s 248 !----------------------------------------------- 249 DOUBLE PRECISION :: TH 250 !----------------------------------------------- 251 ! 252 IF (PHI_C>=ZERO .AND. PHI_C<ONE) THEN !monotonic region 253 TH = PHI_C/(ONE - PHI_C) 254 MUSCL = HALF*MAX(ZERO,MIN(2.0D0*TH,(ONE + TH)/2.0D0,2.0D0)) 255 ELSE IF (PHI_C == ONE) THEN 256 MUSCL = ONE 257 ELSE !first order upwinding 258 MUSCL = ZERO 259 ENDIF 260 ! 261 RETURN 262 END FUNCTION MUSCL 263 264 DOUBLE PRECISION FUNCTION Chi_MUSCL (PHI_C, Chi) 265 ! calculate DWF from Chi_MUSCL scheme 266 !----------------------------------------------- 267 ! M o d u l e s 268 !----------------------------------------------- 269 USE param 270 USE param1 271 IMPLICIT NONE 272 !----------------------------------------------- 273 ! D u m m y A r g u m e n t s 274 !----------------------------------------------- 275 DOUBLE PRECISION PHI_C, Chi 276 277 if(PHI_C < ONE)then 278 Chi_MUSCL = Chi /(4.D0 * (ONE - PHI_C)) 279 elseif(Chi == zero)then ! insures that all species equations uses the same chi 280 Chi_MUSCL = zero 281 else 282 Chi_MUSCL = ONE 283 endif 284 285 RETURN 286 END FUNCTION Chi_MUSCL 287 288 DOUBLE PRECISION FUNCTION Chi4MUSCL (PHI_C, PHIU, PHIC, PHID) 289 ! calculate CHI for MUSCL scheme 290 !----------------------------------------------- 291 ! M o d u l e s 292 !----------------------------------------------- 293 USE param 294 USE param1 295 IMPLICIT NONE 296 !----------------------------------------------- 297 ! D u m m y A r g u m e n t s 298 !----------------------------------------------- 299 DOUBLE PRECISION PHI_C, PHIU, PHIC, PHID 300 301 IF(PHI_C > ZERO .AND. PHI_C <= 0.25d0)THEN 302 Chi4MUSCL = 4.D0 * PHI_C 303 ELSEIF(PHI_C > 0.25d0 .AND. PHI_C <= 0.75d0)THEN 304 Chi4MUSCL = ONE 305 ELSEIF(PHI_C > 0.75d0 .AND. PHI_C <= ONE)THEN 306 Chi4MUSCL = 4.D0*(ONE - PHI_C) 307 ELSEIF( (PHIU < 1d-15) .AND. (PHIC < 1d-15) .AND. (PHID < 1d-15) )THEN 308 ! if a species Xg is less than machine precision, do not use its chi. 309 ! this will avoid calculating chi = 0 for a vanishing species. 310 Chi4MUSCL = LARGE_NUMBER 311 ELSE 312 Chi4MUSCL = ZERO 313 ENDIF 314 315 RETURN 316 END FUNCTION Chi4MUSCL 317 ! 318 ! 319 ! 320 DOUBLE PRECISION FUNCTION VANLEER (PHI_C) 321 !...Translated by Pacific-Sierra Research VAST-90 2.06G5 12:17:31 12/09/98 322 !...Switches: -xf 323 !----------------------------------------------- 324 ! M o d u l e s 325 !----------------------------------------------- 326 USE param 327 USE param1 328 IMPLICIT NONE 329 !----------------------------------------------- 330 ! D u m m y A r g u m e n t s 331 !----------------------------------------------- 332 DOUBLE PRECISION PHI_C 333 !----------------------------------------------- 334 ! L o c a l P a r a m e t e r s 335 !----------------------------------------------- 336 !----------------------------------------------- 337 ! L o c a l V a r i a b l e s 338 !----------------------------------------------- 339 !----------------------------------------------- 340 ! 341 IF (PHI_C>=ZERO .AND. PHI_C<=ONE) THEN !monotonic region 342 VANLEER = PHI_C 343 ELSE !first order upwinding 344 VANLEER = ZERO 345 ENDIF 346 ! 347 RETURN 348 END FUNCTION VANLEER 349 ! 350 ! 351 ! 352 DOUBLE PRECISION FUNCTION MINMOD (PHI_C) 353 !...Translated by Pacific-Sierra Research VAST-90 2.06G5 12:17:31 12/09/98 354 !...Switches: -xf 355 !----------------------------------------------- 356 ! M o d u l e s 357 !----------------------------------------------- 358 USE param 359 USE param1 360 IMPLICIT NONE 361 !----------------------------------------------- 362 ! D u m m y A r g u m e n t s 363 !----------------------------------------------- 364 DOUBLE PRECISION PHI_C 365 !----------------------------------------------- 366 ! L o c a l P a r a m e t e r s 367 !----------------------------------------------- 368 !----------------------------------------------- 369 ! L o c a l V a r i a b l e s 370 !----------------------------------------------- 371 !----------------------------------------------- 372 ! 373 IF (PHI_C>=ZERO .AND. PHI_C<=ONE) THEN !monotonic region 374 IF (PHI_C >= HALF) THEN !central differencing 375 MINMOD = HALF 376 ELSE !second order upwinding 377 MINMOD = HALF*PHI_C/(ONE - PHI_C) 378 ENDIF 379 ELSE !first order upwinding 380 MINMOD = ZERO 381 ENDIF 382 ! 383 RETURN 384 END FUNCTION MINMOD 385 ! 386 387 388 389 390 ! Central scheme for MMS cases. 391 DOUBLE PRECISION FUNCTION CENTRAL_SCHEME (PHI_C) 392 !----------------------------------------------- 393 ! M o d u l e s 394 !----------------------------------------------- 395 USE param 396 USE param1 397 398 IMPLICIT NONE 399 !----------------------------------------------- 400 ! D u m m y A r g u m e n t s 401 !----------------------------------------------- 402 DOUBLE PRECISION PHI_C 403 ! 404 CENTRAL_SCHEME = HALF 405 406 RETURN 407 END FUNCTION CENTRAL_SCHEME 408 409 ! 410 DOUBLE PRECISION FUNCTION UMIST (PHI_C) 411 !...Translated by Pacific-Sierra Research VAST-90 2.06G5 12:17:31 12/09/98 412 !...Switches: -xf 413 !----------------------------------------------- 414 ! M o d u l e s 415 !----------------------------------------------- 416 USE param 417 USE param1 418 IMPLICIT NONE 419 !----------------------------------------------- 420 ! D u m m y A r g u m e n t s 421 !----------------------------------------------- 422 DOUBLE PRECISION PHI_C 423 !----------------------------------------------- 424 ! L o c a l P a r a m e t e r s 425 !----------------------------------------------- 426 !----------------------------------------------- 427 ! L o c a l V a r i a b l e s 428 !----------------------------------------------- 429 DOUBLE PRECISION :: TH 430 !----------------------------------------------- 431 ! 432 IF (PHI_C>=ZERO .AND. PHI_C<ONE) THEN !monotonic region 433 TH = PHI_C/(ONE - PHI_C) 434 UMIST = HALF*MAX(ZERO,MIN(2.0D0*TH,0.75D0 + 0.25D0*TH,2.0D0)) 435 ELSE IF (PHI_C == ONE) THEN 436 UMIST = ONE 437 ELSE !first order upwinding 438 UMIST = ZERO 439 ENDIF 440 ! 441 RETURN 442 END FUNCTION UMIST 443 ! 444 !---------------------------------------------------------------------- 445 ! 446 DOUBLE PRECISION FUNCTION XSI (V, DWF) 447 !...Translated by Pacific-Sierra Research VAST-90 2.06G5 12:17:31 12/09/98 448 !...Switches: -xf 449 !----------------------------------------------- 450 ! M o d u l e s 451 !----------------------------------------------- 452 USE param 453 USE param1 454 IMPLICIT NONE 455 !----------------------------------------------- 456 ! D u m m y A r g u m e n t s 457 !----------------------------------------------- 458 DOUBLE PRECISION V, DWF 459 !----------------------------------------------- 460 ! L o c a l P a r a m e t e r s 461 !----------------------------------------------- 462 !----------------------------------------------- 463 ! L o c a l V a r i a b l e s 464 !----------------------------------------------- 465 !----------------------------------------------- 466 ! 467 IF (V >= ZERO) THEN 468 XSI = DWF 469 ELSE 470 XSI = ONE - DWF 471 ENDIF 472 ! 473 RETURN 474 END FUNCTION XSI 475 ! 476 DOUBLE PRECISION FUNCTION PHI_C_OF (PHI_U, PHI_C, PHI_D) 477 !...Translated by Pacific-Sierra Research VAST-90 2.06G5 12:17:31 12/09/98 478 !...Switches: -xf 479 !----------------------------------------------- 480 ! M o d u l e s 481 !----------------------------------------------- 482 USE param 483 USE param1 484 USE toleranc 485 IMPLICIT NONE 486 !----------------------------------------------- 487 ! D u m m y A r g u m e n t s 488 !----------------------------------------------- 489 DOUBLE PRECISION PHI_U, PHI_C, PHI_D 490 !----------------------------------------------- 491 ! L o c a l P a r a m e t e r s 492 !----------------------------------------------- 493 !----------------------------------------------- 494 ! L o c a l V a r i a b l e s 495 !----------------------------------------------- 496 DOUBLE PRECISION :: DEN 497 !----------------------------------------------- 498 ! 499 IF (COMPARE(PHI_D,PHI_U)) THEN 500 PHI_C_OF = ZERO 501 ELSE 502 DEN = PHI_D - PHI_U 503 PHI_C_OF = (PHI_C - PHI_U)/DEN 504 ENDIF 505 ! 506 RETURN 507 END FUNCTION PHI_C_OF 508 ! 509 DOUBLE PRECISION FUNCTION FPFOI_OF (PHI_D, PHI_C, & 510 PHI_U, PHI_UU) 511 !----------------------------------------------- 512 ! M o d u l e s 513 !----------------------------------------------- 514 USE param 515 USE param1 516 IMPLICIT NONE 517 !----------------------------------------------- 518 ! D u m m y A r g u m e n t s 519 !----------------------------------------------- 520 DOUBLE PRECISION PHI_D, PHI_C, PHI_U, PHI_UU 521 !----------------------------------------------- 522 ! L o c a l P a r a m e t e r s 523 !----------------------------------------------- 524 !----------------------------------------------- 525 ! L o c a l V a r i a b l e s 526 !----------------------------------------------- 527 ! 528 FPFOI_OF = PHI_C + (5.0D0/16.0D0)*(PHI_D - PHI_C) + & 529 (1.0D0/4.0D0)*(PHI_C - PHI_U) - & 530 (1.0D0/16.0D0)*(PHI_U - PHI_UU) 531 ! 532 ! LIMIT THE HIGH ORDER INTERPOLATION 533 ! 534 FPFOI_OF = UNIV_LIMITER_OF(FPFOI_OF , PHI_D, PHI_C, PHI_U) 535 RETURN 536 END FUNCTION FPFOI_OF 537 ! 538 ! 539 DOUBLE PRECISION FUNCTION UNIV_LIMITER_OF (PHI_TEMP, PHI_D, PHI_C, & 540 PHI_U) 541 ! 542 !----------------------------------------------- 543 ! M o d u l e s 544 !----------------------------------------------- 545 USE param 546 USE param1 547 USE run 548 ! USE fldvar 549 IMPLICIT NONE 550 !----------------------------------------------- 551 ! D u m m y A r g u m e n t s 552 !----------------------------------------------- 553 DOUBLE PRECISION PHI_TEMP, PHI_D, PHI_C, PHI_U 554 !----------------------------------------------- 555 ! L o c a l P a r a m e t e r s 556 !----------------------------------------------- 557 !----------------------------------------------- 558 ! L o c a l V a r i a b l e s 559 !----------------------------------------------- 560 DOUBLE PRECISION PHI_REF, DEL, CURV 561 !----------------------------------------------- 562 ! E x t e r n a l F u n c t i o n s 563 !----------------------------------------------- 564 !----------------------------------------------- 565 ! 566 ! 567 DEL = PHI_D - PHI_U 568 CURV = PHI_D - 2.0*PHI_C + PHI_U 569 IF (ABS (CURV) >= ABS (DEL)) THEN ! NON-MONOTONIC 570 UNIV_LIMITER_OF = PHI_C 571 ELSE 572 PHI_REF = PHI_U + (PHI_C-PHI_U)/C_FAC 573 IF (DEL > ZERO) THEN 574 IF (PHI_TEMP < PHI_C)UNIV_LIMITER_OF = PHI_C 575 IF (PHI_TEMP > MIN(PHI_REF,PHI_D)) & 576 UNIV_LIMITER_OF = MIN(PHI_REF,PHI_D) 577 IF (PHI_TEMP >= PHI_C .AND. PHI_TEMP <= MIN(PHI_REF,PHI_D)) & 578 UNIV_LIMITER_OF = PHI_TEMP 579 ELSE 580 IF (PHI_TEMP > PHI_C)UNIV_LIMITER_OF = PHI_C 581 IF (PHI_TEMP < MAX(PHI_REF,PHI_D)) & 582 UNIV_LIMITER_OF = MAX(PHI_REF,PHI_D) 583 IF (PHI_TEMP >= MAX(PHI_REF,PHI_D) .AND. PHI_TEMP <= PHI_C) & 584 UNIV_LIMITER_OF = PHI_TEMP 585 ENDIF 586 ENDIF 587 RETURN 588 END FUNCTION UNIV_LIMITER_OF 589 590 END MODULE discretization 591