Changes between Version 17 and Version 18 of u/erica/ExactRiemannSolver


Ignore:
Timestamp:
01/04/13 23:40:39 (13 years ago)
Author:
Erica Kaminski
Comment:

—

Legend:

Unmodified
Added
Removed
Modified

  • u/erica/ExactRiemannSolver

    v17 v18  
    1010
    1111=== Exact Riemann Solver ===
    12 PROGRAM ExactRiemannSolver
     12Attached here is the Exact RS program I wrote. It reproduces the results for the Toro tests, as can be seen in the next section. The program is broken up into the following routines:
    1313
    14   !-------------------------------------------------------------------------------------------------
    15   ! This program finds the solution of p*, u*, and rho*_L and rho*_R given initial data states wL
    16   ! and wR of the Riemann problem
    17   ! Author: Erica Kaminski, 12.17.12
    18   !-------------------------------------------------------------------------------------------------
    19 
    20   IMPLICIT NONE
    21 
    22   REAL :: DomLen, TimeOut, DisPos
    23   INTEGER :: cells, i
    24   REAL :: rho_L, p_L, u_L                                             !Left data vars               
    25   REAL :: rho_R, p_R, u_R                                             !Right data vars
    26   REAL, dimension(3) :: wL                                            !Left data state
    27   REAL, dimension(3) :: wR                                            !Right data state
    28   REAL :: Cs_L, Cs_R                                                  !Sound speed Left and Right states
    29   REAL, PARAMETER :: gamma = 1.4                                      ! Ratio specific heats
    30   REAL :: p_Star, u_Star                                              ! Output of program
    31   REAL :: mpa                                                         ! Normalization constant for P
    32   REAL :: u_diff, dx, xpos, s
    33   REAL :: rhos, us, ps
    34 
    35   NAMELIST /InitialData/ DomLen, TimeOut, DisPos, cells, rho_L, p_L, u_L, rho_R, p_R, u_R, mpa
    36 
    37   !--------------------------------------------------------------------------------------------------
    38   ! Main body of program here
    39   !--------------------------------------------------------------------------------------------------
    40 
    41 
    42   OPEN(UNIT = 3, File = "initial.data", status = 'old')              !Initialize data
    43   READ(3, NML = InitialData) 
    44   CLOSE(3)
    45 
    46   wL = (/rho_L, p_L, u_L/)                                           !Populate w arrays
    47   wR = (/rho_R, p_R, u_R/)
    48 
    49   Cs_L = Sqrt(Gamma*p_L/rho_L)
    50   Cs_R = Sqrt(Gamma*p_R/rho_R)
    51 
    52   U_DIFF = u_R - u_L
    53 
    54   !PRESSURE POSITIVITY CONDITION ---
    55 
    56   IF (2.0*Cs_L/(gamma-1.0) + 2.0*Cs_R/(gamma-1.0).LE.U_DIFF) THEN
    57 
    58      PRINT *, 'PRESSURE POSITIVITY CONDITION VIOLATED!'
    59      PRINT *, 'VACUUM STATE IS CREATED BY INITIAL CONDITIONS'
    60      PRINT *, 'THIS RIEMANN SOLVER IS NOT CORRECT METHOD FOR SOLUTION'
    61      PRINT *, 'PROGRAM STOPPED'
    62      STOP
    63 
    64   END IF
    65 
    66   CALL Find_Pstar(p_Star, u_Star)
    67   !check that p_Star = p from subprpgram with a print statement - I did, it does
    68 
    69   OPEN(UNIT=2, FILE='exact.out', STATUS='OLD', POSITION='APPEND')
    70 
    71   dx = DomLen/REAL(cells)
    72   WRITE(2,*)'xpos     rho     u     P     iE'
    73   DO 10 i = 1, cells
    74      xpos = dx*(REAL(i)-0.5) !cell-centered
    75      s = (xpos - DisPos)/TimeOut !characteristic speed
    76      CALL SampleExact(p_Star, u_Star, s, rhos, us, ps)
    77      WRITE(2, 20)xpos, rhos, us, ps/mpa, ps/rhos/(gamma-1.0)/mpa !the last may be internal energy...?   
    78 10   Continue
    79 
    80      CLOSE(2)
    81 
    82 20   FORMAT(5(F14.6,2X))
    83 
    84 
    85 
    86    CONTAINS
    87 
    88      SUBROUTINE SampleExact(p_Star, u_Star, s, rhos, us, ps)
    89        REAL,INTENT(IN) :: p_Star, u_Star, s
    90        REAL,INTENT(OUT) ::rhos, us, ps   !sampled density, velocity, pressure at point (x,t)
    91        REAL :: SHL, STL !Speed of left rarefaction head and tail
    92        REAL :: SHR, STR !Speed of right rarefaction head and tail
    93        REAL :: CsL_Star, CsR_Star !Sound speed behind left and right rarefactions
    94        REAL :: SL, SR   !Left and right shock speeds
    95        REAL :: G1, G2, G3, G4, G5, G6, ppl, ppr
    96 
    97        CsL_Star = Cs_L*(p_Star/p_L)**(G3/G2)
    98        CsR_Star = Cs_R*(p_Star/p_R)**(G3/G2)
    99 
    100        G1 = gamma + 1.0
    101        G2 = 2.0*gamma
    102        G3 = gamma - 1.0
    103        ppl = p_Star/p_L
    104        ppr = p_Star/p_R
    105        G4 = G3/G1
    106        G5 = G1/G2
    107        G6 = G3/G2
    108 
    109        IF(s.LE.u_Star)THEN
    110           ! (x,t) is to the left of contact
    111           IF(p_L.GE.p_Star)THEN
    112              !Left wave is a rarefaction, based on initial data
    113              !3 possible left data states then - WL, WL_Fan, W_Star - next we determine which
    114 !             PRINT*,'IT"S A RAREFACTION'
    115 
    116              !Are these negative quantities:?
    117              SHL = u_L - Cs_L
    118              STL = u_Star - CsL_Star
    119 
    120              IF(S.LE.SHL)THEN
    121 !                PRINT *, 'LEFT DATA STATE'
    122                 !Left data state
    123                 rhos = rho_L
    124                 us = u_L
    125                 ps = p_L
    126              ELSE
    127 
    128                 !                IF((SHL.GE.S).AND.(S.LE.STL))THEN
    129                 IF((SHL.LE.S).AND.(S.LE.STL)) THEN
    130 !                   PRINT*, 'iNSIDE RARE WAVE'
    131                    !Inside rarefaction fan
    132                    rhos = rho_L*(2/G1 + G3/(G1*Cs_L)*(u_L - s))**(2/G3)
    133                    us = 2/G1*(Cs_L + (G3/2)*u_L + s)
    134                    ps = p_L*( 2/G1 + G3/(G1*Cs_L)*(u_L - s) )**(G2/G3)
    135                 ELSE
    136                    !In left star region
    137 !                   PRINT *, 'IN STAR REGION'
    138                    rhos = rho_L*ppl**(1.0/gamma) !follows from isentropic law
    139                    us = u_Star
    140                    ps = p_Star
    141                 END IF
    142 
    143              END IF
    144 
    145           ELSE
    146              !Left wave is a shock
    147              !Two states possible - pre-shock, post-shock
    148 !             PRINT*,'IT"S A SHOCK !  ! !'
    149              SL = (u_L - Cs_L)*SQRT(G1/G2*ppl + G3/G2)
    150              IF(s.LE.SL)THEN
    151                 !Data state is left state/pre-shock
    152                 rhos = rho_L
    153                 us = u_L
    154                 ps = p_L
    155              ELSE
    156                 !Data state is post-shock
    157                 rhos = rho_L*((ppl + G3/G1)/((G3/G1)*ppl + 1.0))
    158                 us = u_Star
    159                 ps = p_Star
    160              END IF
    161           END IF
    162 
    163        ELSE
    164 
    165           !(x,t) is on right of contact
    166           IF(p_Star.LE.p_R)THEN
    167              !Rarefaction, and 3 wave states possible: wR, WRFan, W*R
    168              SHR = Cs_R + u_R
    169              STR = CsR_Star + u_Star
    170              IF(S.LE.STR)THEN
    171                 !In r star region
    172                 rhos = rho_R*ppr**(1.0/gamma) !follows from isentropic law -- CHECK?
    173                 us = u_Star
    174                 ps = p_Star
    175              ELSE
    176                 IF((STR.LE.S).AND.(S.LE.SHR))THEN
    177                    !Inside right rarefaction fan
    178                    rhos = rho_R*(2/G1 - G3/(G1*Cs_R)*(u_R - s))**(2/G3)
    179                    us = 2/G1*(-Cs_R + (G3/2)*u_R + s)
    180                    ps = p_R*( 2/G1 - G3/(G1*Cs_R)*(u_R - s) )**(G2/G3)
    181                    !Check these expressions
    182                 ELSE
    183                    !In right data state
    184                    rhos = rho_R
    185                    us = u_R
    186                    ps = p_R
    187                 END IF
    188              END IF
    189           ELSE
    190              !Right wave is a shock
    191              SR = (u_R + Cs_R)*SQRT(G1/G2*ppr + G3/G2) 
    192              IF(s.LE.SR)THEN
    193                 !In post-shock region
    194                 rhos = rho_R*((ppr + G3/G1)/((G3/G1)*ppr + 1.0))
    195                 us = u_Star
    196                 ps = p_Star
    197              ELSE
    198                 !In pre-shock region
    199                 rhos = rho_R
    200                 us = u_R
    201                 ps = p_R
    202              END IF
    203           END IF
    204        END IF
    205 
    206 
    207      END SUBROUTINE SampleExact
    208 
    209 
    210      SUBROUTINE Find_Pstar(p, u)
    211        !Does mpa NEED to be read into this subroutine, if it is already initialized in main body of program?
    212        !-----------------------------------------------------------------------------------------------
    213        ! Implements Newton's iterative method to solve for p_star
    214        !-----------------------------------------------------------------------------------------------
    215        REAL, INTENT(OUT) :: p, u
    216        REAL :: p0 ! Initial estimate of P_star
    217        REAL, PARAMETER :: tol = 0.0000001
    218        REAL :: change ! Relative change in pressure due to iteration
    219        REAL :: pOld, pNew! Cs_L, Cs_R, rho_L, rho_R, p_L, p_R, u_L, u_R -- these should be all accessible as they are global vars
    220        INTEGER :: i, NmrItr
    221        REAL :: f_p, f_pD, FL, FLD, FR, FRD ! pressure function in Toro
    222        REAL :: Delta
    223        NAMELIST /Itr/ NmrItr
    224 
    225        !statement labels are not global:
    226        OPEN(UNIT = 3, FILE = 'initial.data')
    227        READ(3, NML = Itr)
    228        CLOSE(3)
    229 
    230 
    231        !Get estimate for P_star
    232        CALL P_estimate(p0)
    233 
    234 
    235        pOld = p0  ! Initialize pOld
    236 
    237        !Check that the vars are global:
    238 !       PRINT *, 'u_L, p_L, rho_L, cs_L =', u_L, p_L, rho_L, cs_L
    239 
    240        OPEN(UNIT = 2, FILE = 'exact.out', STATUS = 'UNKNOWN')
    241        WRITE(2,*) '--------------------------------------------------------------------------'
    242        WRITE(2,*) '         Iteration number                      Change                     '
    243        WRITE(2,*) '--------------------------------------------------------------------------'
    244 
    245        !Need define Cs_L here too or main program enough?
    246        DO 10 i = 1,NmrItr
    247           CALL PressFunct(FL, FLD, pOld, rho_L, p_L, Cs_L)
    248           CALL PressFunct(FR, FRD, pOld, rho_R, p_R, Cs_R)
    249           f_p = FL + FR + U_diff
    250           f_pD = FLD + FRD
    251           Delta = f_p/f_pD
    252           pNew = pOld - Delta ! This is corrected p.
    253           change = 2.0*ABS((pNew - pOld)/(pNew + pOld))
    254           WRITE(2,*) i, change
    255           IF (change.LE.tol) GOTO 20
    256           IF (pNew.LT.0.0) pNew = tol
    257           pOld = pNew
    258 10        Continue
    259 
    260           WRITE(2,*) 'Newton Rhapson Method Diverged!'
    261 
    262 20        WRITE(2,*) 'Pressure in Star Region and U in star region'
    263           p = pNew/MPA
    264           u = (1.0/2.0)*(u_L + u_R) + (1.0/2.0)*(fR - fL)
    265           WRITE(2,*) 'p_Star = ', p
    266           WRITE(2,*) 'u_Star = ', u
    267 
    268           CLOSE(2)
    269 
    270         END SUBROUTINE Find_Pstar
    271 
    272         SUBROUTINE P_estimate(pEst)
    273           !Employs the subroutine from Toro for estimating initial P_star
    274 
    275           REAL,INTENT(OUT) :: pEst
    276           REAL :: G1, G2, G3, G4, G5, G6, G7, G8
    277           REAL :: CUP, GEL, GER, pMax, pMin, PPV, PQ, PTL, PTR, QMAX, QUSER, UM
    278 
    279           G1 = (gamma - 1.0)/(2.0*gamma)
    280           G2 = (gamma + 1.0)/(2.0*gamma)
    281           G3 = 2.0*gamma/(gamma - 1.0)
    282           G4 = 2.0/(gamma - 1.0)
    283           G5 = 2.0/(gamma + 1.0)
    284           G6 = (gamma - 1.0)/(gamma + 1.0)
    285           G7 = (gamma - 1.0)/2.0
    286           G8 = (gamma - 1.0)
    287 
    288           QUSER = 2.0
    289 
    290           ! Guess P using the PVRS solver
    291 
    292           CUP = 0.25*(rho_L + rho_R)*(Cs_L + Cs_R)
    293           PPV = 0.5*(p_L + p_R) + 0.5*(u_L - u_R)*CUP
    294           PPV = MAX(0.0, PPV)
    295           PMIN = MIN(p_L, p_R)
    296           PMAX = MAX(p_L, p_R)
    297           QMAX = PMAX/PMIN
    298 
    299           IF (QMAX.LE.QUSER.AND.(PMIN.LE.PPV.AND.PPV.LE.PMAX)) THEN
    300              !SELECT PVRS SOLVER
    301              pEst = PPV
    302 
    303           ELSE
    304 
    305              IF(PPV.LT.PMIN) THEN
    306                 !SELECT 2-RAREFACTION RS
    307                 PQ = (p_L/p_R)**G1
    308                 UM = (PQ*u_L/Cs_L + u_R/Cs_R + G4*(PQ - 1.0))/(PQ/Cs_L + 1.0/Cs_R)
    309                 PTL = 1.0 + G7*(u_L - UM)/Cs_L
    310                 PTR = 1.0 + G7*(UM - u_R)/Cs_R
    311                 pEst = 0.5*(p_L*PTL**G3 + p_R*PTR**G3)
    312 
    313              ELSE
    314                 !Select the 2-shock RS with PVRS as estimate
    315 
    316                 GEL = SQRT((G5/rho_L)/(G6*p_L + PPV))
    317                 GER = SQRT((G5/rho_R)/(G6*p_R + PPV))
    318                 pEst = (GEL*p_L + GER*p_R - (u_R - u_L))/(GEL + GER)
    319 
    320              END IF
    321 
    322           END IF
    323 
    324 
    325         END SUBROUTINE P_estimate
    326 
    327         SUBROUTINE PressFunct(f, fD, pItr, rho, p, Cs)
    328 
    329 
    330           REAL, INTENT(OUT) :: f, fd ! Pressure function and its derivative in Toro
    331           REAL, INTENT(IN) :: pItr  ! Pressure correction as computed by Find_Pstar
    332           REAL, INTENT(IN) :: rho, p, Cs !Initial data
    333           REAL :: A, B  ! Constants in pressure function
    334 
    335 
    336           A =  2.0/((gamma+1.0)*rho)
    337           B =  ((gamma-1.0)/(gamma+1.0))*p
    338 
    339 
    340           IF(pItr.LT.p)THEN !Rarefaction wave
    341 
    342              f = (2.0*Cs/(gamma - 1.0))*((pItr/p)**((gamma-1.0)/(2.0*gamma)) - 1.0) 
    343              fd = (1.0/(rho*Cs))*(pItr/p)**(-(gamma+1)/(2*gamma))
    344 
    345           ELSE
    346 
    347              !Shock
    348              f = (pItr - p)*SQRT(A/(pItr+B))
    349              fd = Sqrt(A /(B + pItr))*(1 - (pItr - p)/(2*(B+pItr)))
    350 
    351           END IF
    352 
    353         END SUBROUTINE PressFunct
    354 
    355       END PROGRAM ExactRiemannSolver
    356 
     141.
    35715
    35816=== Results ===
     17
     18These tests were taken from Toro, chapter 4. They are on a domain of length = 1 computational unit, with 1000 zones. The final time in each of them varies, and were chosen to match those in Toro, as were the initial values for the left and right data states. These results match Toro's exactly, except for the slight mismatch in solutions at the rarefaction-tail, contact-discontinuity boundary.
    35919
    360201. [[Image(SodTest2ndAttempt.png, 50%)]]
    … …  
    36525
    366264. [[Image(LSRR.png, 50%)]]
     27
     28=== Questions ===
     29
     30Is the pressure condition for determining types of waves sufficient?