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

Timestamp:

01/04/13 23:28:52 (12 years ago)

Author:

Erica Kaminski

Comment:

—

u/erica/ExactRiemannSolver

@@ -7 +7 @@
 Given the above properties of the Riemann problem, an algebraic expression can be derived which gives p in the central "star-region" (denoted by '*'). The overall structure of the Riemann problem then is to solve this algebraic equation for p*. Once p* is known, u* follows immediately. The remaining rho*_L and rho*_R follow from expressions valid for the specific L- or R- wave present.
 
-Specifically, this means that the code will determine at every point (x,t), that point's relative position to the different waves present. Once the position is determined, the solution is given by analytic expressions for the following five possibilities: pre- or post- shock, ahead of rarefaction head, behind rarefaction tail, or within the rarefaction fan. The position of each sampling point (x,t) is determined by a characteristic speed in the grid given by s = dx/t, where dx is the distance from the initial discontinuity to the sampling point (x,t), and t is the simulation end time. This s for every point on the grid can be compared to the known present waves, as their speeds are known exactly. In this way, the relative position of the sampling point with respect to the waves on the grid is determined.
-
-The program thus will consist of 4 parts: 1. Code to find p* using an iterative process, 2. Comparing p* to p_L and p_R to identify the types of nonlinear waves present in the solution, 3. Find remaining quantities in star region, where u* is given by simple expression in Toro (eqn. 4.9), and wave-specific expressions give rho*_L and rho*_R. 
+Specifically, this means that the code will determine at every point (x,t), that point's relative position to the different waves present. Once the position is determined, the solution is given by analytic expressions for the following five possibilities: pre- or post- shock, ahead of rarefaction head, behind rarefaction tail, or within the rarefaction fan. The position of each sampling point (x,t) is determined by a characteristic speed in the grid given by s = dx/t, where dx is the distance from the initial discontinuity to the sampling point (x,t), and t is the simulation end time. This s for every point on the grid can be compared to the known present waves, as their speeds are known exactly. In this way, the relative position of the sampling point with respect to the waves on the grid is determined. 
 
 === Exact Riemann Solver ===
-
-I. Declarations
- 1. Wl=()
- 2. Wr=()
- 3. tol =
-
-II. Main Program Body
- 1. Call P* routine - find and record p*
- 2. Call Identification routine - find types of waves present, record
- 3. Call U*() - reads in vars and returns U*
- 4. Call rho* - reads in vars and returns rho*
- 2. Print to File - u*, rho*, p*
-
-III. Functions
- 1. FL = ()
- 2. FR = ()
- 3. U* = ()
- 4. rho*L = {}
- 5. rho*R = {}
-
-IV. Subroutines
- 1. Find p* using initial estimate
-  a. type of approximation (1)-(5) for different options (on screen options)
-  b. let p0=(approximation)
-  c. Begin iterative process, that ends once | | < tol
-
- 2. Determine types of waves present
-  a. comparing p* to pL and pR
-  b. series of if statements corresponding to flow chart
-  c. left_wave = 1 for shock, 2 for rarefaction
+PROGRAM ExactRiemannSolver
+
+  !-------------------------------------------------------------------------------------------------
+  ! This program finds the solution of p*, u*, and rho*_L and rho*_R given initial data states wL 
+  ! and wR of the Riemann problem
+  ! Author: Erica Kaminski, 12.17.12
+  !-------------------------------------------------------------------------------------------------
+
+  IMPLICIT NONE
+
+  REAL :: DomLen, TimeOut, DisPos
+  INTEGER :: cells, i
+  REAL :: rho_L, p_L, u_L                                             !Left data vars                
+  REAL :: rho_R, p_R, u_R                                             !Right data vars
+  REAL, dimension(3) :: wL                                            !Left data state
+  REAL, dimension(3) :: wR                                            !Right data state
+  REAL :: Cs_L, Cs_R                                                  !Sound speed Left and Right states
+  REAL, PARAMETER :: gamma = 1.4                                      ! Ratio specific heats
+  REAL :: p_Star, u_Star                                              ! Output of program
+  REAL :: mpa                                                         ! Normalization constant for P
+  REAL :: u_diff, dx, xpos, s
+  REAL :: rhos, us, ps
+
+  NAMELIST /InitialData/ DomLen, TimeOut, DisPos, cells, rho_L, p_L, u_L, rho_R, p_R, u_R, mpa
+
+  !--------------------------------------------------------------------------------------------------
+  ! Main body of program here 
+  !--------------------------------------------------------------------------------------------------
+
+
+  OPEN(UNIT = 3, File = "initial.data", status = 'old')              !Initialize data 
+  READ(3, NML = InitialData)  
+  CLOSE(3)
+
+  wL = (/rho_L, p_L, u_L/)                                           !Populate w arrays
+  wR = (/rho_R, p_R, u_R/)
+
+  Cs_L = Sqrt(Gamma*p_L/rho_L)
+  Cs_R = Sqrt(Gamma*p_R/rho_R)
+
+  U_DIFF = u_R - u_L
+
+  !PRESSURE POSITIVITY CONDITION ---
+
+  IF (2.0*Cs_L/(gamma-1.0) + 2.0*Cs_R/(gamma-1.0).LE.U_DIFF) THEN
+
+     PRINT *, 'PRESSURE POSITIVITY CONDITION VIOLATED!'
+     PRINT *, 'VACUUM STATE IS CREATED BY INITIAL CONDITIONS'
+     PRINT *, 'THIS RIEMANN SOLVER IS NOT CORRECT METHOD FOR SOLUTION'
+     PRINT *, 'PROGRAM STOPPED'
+     STOP
+
+  END IF
+
+  CALL Find_Pstar(p_Star, u_Star)
+  !check that p_Star = p from subprpgram with a print statement - I did, it does
+
+  OPEN(UNIT=2, FILE='exact.out', STATUS='OLD', POSITION='APPEND')
+
+  dx = DomLen/REAL(cells)
+  WRITE(2,*)'xpos     rho     u     P     iE'
+  DO 10 i = 1, cells
+     xpos = dx*(REAL(i)-0.5) !cell-centered
+     s = (xpos - DisPos)/TimeOut !characteristic speed
+     CALL SampleExact(p_Star, u_Star, s, rhos, us, ps)
+     WRITE(2, 20)xpos, rhos, us, ps/mpa, ps/rhos/(gamma-1.0)/mpa !the last may be internal energy...?   
+10   Continue
+
+     CLOSE(2)
+
+20   FORMAT(5(F14.6,2X))
+
+
+
+   CONTAINS
+
+     SUBROUTINE SampleExact(p_Star, u_Star, s, rhos, us, ps)
+       REAL,INTENT(IN) :: p_Star, u_Star, s
+       REAL,INTENT(OUT) ::rhos, us, ps   !sampled density, velocity, pressure at point (x,t)
+       REAL :: SHL, STL !Speed of left rarefaction head and tail
+       REAL :: SHR, STR !Speed of right rarefaction head and tail
+       REAL :: CsL_Star, CsR_Star !Sound speed behind left and right rarefactions
+       REAL :: SL, SR   !Left and right shock speeds
+       REAL :: G1, G2, G3, G4, G5, G6, ppl, ppr
+
+       CsL_Star = Cs_L*(p_Star/p_L)**(G3/G2)
+       CsR_Star = Cs_R*(p_Star/p_R)**(G3/G2)
+
+       G1 = gamma + 1.0
+       G2 = 2.0*gamma
+       G3 = gamma - 1.0
+       ppl = p_Star/p_L
+       ppr = p_Star/p_R
+       G4 = G3/G1
+       G5 = G1/G2
+       G6 = G3/G2 
+
+       IF(s.LE.u_Star)THEN
+          ! (x,t) is to the left of contact
+          IF(p_L.GE.p_Star)THEN
+             !Left wave is a rarefaction, based on initial data
+             !3 possible left data states then - WL, WL_Fan, W_Star - next we determine which 
+!             PRINT*,'IT"S A RAREFACTION'
+
+             !Are these negative quantities:?
+             SHL = u_L - Cs_L
+             STL = u_Star - CsL_Star
+
+             IF(S.LE.SHL)THEN
+!                PRINT *, 'LEFT DATA STATE'
+                !Left data state
+                rhos = rho_L
+                us = u_L
+                ps = p_L
+             ELSE 
+
+                !                IF((SHL.GE.S).AND.(S.LE.STL))THEN
+                IF((SHL.LE.S).AND.(S.LE.STL)) THEN
+!                   PRINT*, 'iNSIDE RARE WAVE'
+                   !Inside rarefaction fan
+                   rhos = rho_L*(2/G1 + G3/(G1*Cs_L)*(u_L - s))**(2/G3)
+                   us = 2/G1*(Cs_L + (G3/2)*u_L + s)
+                   ps = p_L*( 2/G1 + G3/(G1*Cs_L)*(u_L - s) )**(G2/G3)
+                ELSE
+                   !In left star region
+!                   PRINT *, 'IN STAR REGION'
+                   rhos = rho_L*ppl**(1.0/gamma) !follows from isentropic law
+                   us = u_Star
+                   ps = p_Star
+                END IF
+
+             END IF
+
+          ELSE
+             !Left wave is a shock
+             !Two states possible - pre-shock, post-shock
+!             PRINT*,'IT"S A SHOCK !  ! !'
+             SL = (u_L - Cs_L)*SQRT(G1/G2*ppl + G3/G2)
+             IF(s.LE.SL)THEN
+                !Data state is left state/pre-shock
+                rhos = rho_L
+                us = u_L
+                ps = p_L
+             ELSE
+                !Data state is post-shock
+                rhos = rho_L*((ppl + G3/G1)/((G3/G1)*ppl + 1.0))
+                us = u_Star
+                ps = p_Star
+             END IF
+          END IF
+
+       ELSE
+
+          !(x,t) is on right of contact
+          IF(p_Star.LE.p_R)THEN
+             !Rarefaction, and 3 wave states possible: wR, WRFan, W*R
+             SHR = Cs_R + u_R
+             STR = CsR_Star + u_Star
+             IF(S.LE.STR)THEN
+                !In r star region
+                rhos = rho_R*ppr**(1.0/gamma) !follows from isentropic law -- CHECK?
+                us = u_Star
+                ps = p_Star
+             ELSE
+                IF((STR.LE.S).AND.(S.LE.SHR))THEN
+                   !Inside right rarefaction fan
+                   rhos = rho_R*(2/G1 - G3/(G1*Cs_R)*(u_R - s))**(2/G3)
+                   us = 2/G1*(-Cs_R + (G3/2)*u_R + s)
+                   ps = p_R*( 2/G1 - G3/(G1*Cs_R)*(u_R - s) )**(G2/G3)
+                   !Check these expressions
+                ELSE
+                   !In right data state
+                   rhos = rho_R
+                   us = u_R
+                   ps = p_R
+                END IF
+             END IF
+          ELSE
+             !Right wave is a shock
+             SR = (u_R + Cs_R)*SQRT(G1/G2*ppr + G3/G2)  
+             IF(s.LE.SR)THEN
+                !In post-shock region
+                rhos = rho_R*((ppr + G3/G1)/((G3/G1)*ppr + 1.0)) 
+                us = u_Star
+                ps = p_Star
+             ELSE 
+                !In pre-shock region
+                rhos = rho_R
+                us = u_R
+                ps = p_R
+             END IF
+          END IF
+       END IF
+
+
+     END SUBROUTINE SampleExact
+
+
+     SUBROUTINE Find_Pstar(p, u)
+       !Does mpa NEED to be read into this subroutine, if it is already initialized in main body of program?
+       !-----------------------------------------------------------------------------------------------
+       ! Implements Newton's iterative method to solve for p_star
+       !-----------------------------------------------------------------------------------------------
+       REAL, INTENT(OUT) :: p, u 
+       REAL :: p0 ! Initial estimate of P_star
+       REAL, PARAMETER :: tol = 0.0000001 
+       REAL :: change ! Relative change in pressure due to iteration
+       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
+       INTEGER :: i, NmrItr
+       REAL :: f_p, f_pD, FL, FLD, FR, FRD ! pressure function in Toro
+       REAL :: Delta
+       NAMELIST /Itr/ NmrItr
+
+       !statement labels are not global:
+       OPEN(UNIT = 3, FILE = 'initial.data')
+       READ(3, NML = Itr)
+       CLOSE(3)
+
+
+       !Get estimate for P_star
+       CALL P_estimate(p0)
+
+
+       pOld = p0  ! Initialize pOld
+
+       !Check that the vars are global:
+!       PRINT *, 'u_L, p_L, rho_L, cs_L =', u_L, p_L, rho_L, cs_L 
+
+       OPEN(UNIT = 2, FILE = 'exact.out', STATUS = 'UNKNOWN')
+       WRITE(2,*) '--------------------------------------------------------------------------'
+       WRITE(2,*) '         Iteration number                      Change                     '
+       WRITE(2,*) '--------------------------------------------------------------------------'
+
+       !Need define Cs_L here too or main program enough?
+       DO 10 i = 1,NmrItr
+          CALL PressFunct(FL, FLD, pOld, rho_L, p_L, Cs_L)
+          CALL PressFunct(FR, FRD, pOld, rho_R, p_R, Cs_R)
+          f_p = FL + FR + U_diff
+          f_pD = FLD + FRD
+          Delta = f_p/f_pD
+          pNew = pOld - Delta ! This is corrected p.
+          change = 2.0*ABS((pNew - pOld)/(pNew + pOld))
+          WRITE(2,*) i, change
+          IF (change.LE.tol) GOTO 20 
+          IF (pNew.LT.0.0) pNew = tol
+          pOld = pNew
+10        Continue 
+
+          WRITE(2,*) 'Newton Rhapson Method Diverged!'
+
+20        WRITE(2,*) 'Pressure in Star Region and U in star region'
+          p = pNew/MPA
+          u = (1.0/2.0)*(u_L + u_R) + (1.0/2.0)*(fR - fL)
+          WRITE(2,*) 'p_Star = ', p
+          WRITE(2,*) 'u_Star = ', u
+
+          CLOSE(2)
+
+        END SUBROUTINE Find_Pstar
+
+        SUBROUTINE P_estimate(pEst)
+          !Employs the subroutine from Toro for estimating initial P_star
+
+          REAL,INTENT(OUT) :: pEst 
+          REAL :: G1, G2, G3, G4, G5, G6, G7, G8
+          REAL :: CUP, GEL, GER, pMax, pMin, PPV, PQ, PTL, PTR, QMAX, QUSER, UM
+
+          G1 = (gamma - 1.0)/(2.0*gamma)
+          G2 = (gamma + 1.0)/(2.0*gamma)
+          G3 = 2.0*gamma/(gamma - 1.0)
+          G4 = 2.0/(gamma - 1.0)
+          G5 = 2.0/(gamma + 1.0)
+          G6 = (gamma - 1.0)/(gamma + 1.0)
+          G7 = (gamma - 1.0)/2.0
+          G8 = (gamma - 1.0)
+
+          QUSER = 2.0
+
+          ! Guess P using the PVRS solver
+
+          CUP = 0.25*(rho_L + rho_R)*(Cs_L + Cs_R)
+          PPV = 0.5*(p_L + p_R) + 0.5*(u_L - u_R)*CUP
+          PPV = MAX(0.0, PPV)
+          PMIN = MIN(p_L, p_R)
+          PMAX = MAX(p_L, p_R)
+          QMAX = PMAX/PMIN
+
+          IF (QMAX.LE.QUSER.AND.(PMIN.LE.PPV.AND.PPV.LE.PMAX)) THEN
+             !SELECT PVRS SOLVER
+             pEst = PPV 
+
+          ELSE
+
+             IF(PPV.LT.PMIN) THEN
+                !SELECT 2-RAREFACTION RS
+                PQ = (p_L/p_R)**G1
+                UM = (PQ*u_L/Cs_L + u_R/Cs_R + G4*(PQ - 1.0))/(PQ/Cs_L + 1.0/Cs_R)
+                PTL = 1.0 + G7*(u_L - UM)/Cs_L
+                PTR = 1.0 + G7*(UM - u_R)/Cs_R
+                pEst = 0.5*(p_L*PTL**G3 + p_R*PTR**G3)
+
+             ELSE
+                !Select the 2-shock RS with PVRS as estimate
+
+                GEL = SQRT((G5/rho_L)/(G6*p_L + PPV))
+                GER = SQRT((G5/rho_R)/(G6*p_R + PPV))
+                pEst = (GEL*p_L + GER*p_R - (u_R - u_L))/(GEL + GER)
+
+             END IF
+
+          END IF
+
+
+        END SUBROUTINE P_estimate
+
+        SUBROUTINE PressFunct(f, fD, pItr, rho, p, Cs)
+
+
+          REAL, INTENT(OUT) :: f, fd ! Pressure function and its derivative in Toro
+          REAL, INTENT(IN) :: pItr  ! Pressure correction as computed by Find_Pstar
+          REAL, INTENT(IN) :: rho, p, Cs !Initial data
+          REAL :: A, B  ! Constants in pressure function
+
+
+          A =  2.0/((gamma+1.0)*rho)
+          B =  ((gamma-1.0)/(gamma+1.0))*p
+
+
+          IF(pItr.LT.p)THEN !Rarefaction wave
+
+             f = (2.0*Cs/(gamma - 1.0))*((pItr/p)**((gamma-1.0)/(2.0*gamma)) - 1.0)  
+             fd = (1.0/(rho*Cs))*(pItr/p)**(-(gamma+1)/(2*gamma))
+
+          ELSE
+
+             !Shock
+             f = (pItr - p)*SQRT(A/(pItr+B)) 
+             fd = Sqrt(A /(B + pItr))*(1 - (pItr - p)/(2*(B+pItr)))
+
+          END IF
+
+        END SUBROUTINE PressFunct
+
+      END PROGRAM ExactRiemannSolver
+
 
 === Results ===