riemann.f

********************************************************************************
* This file is part of python_finite_volume_solver
* Copyright (C) 2017 Kenneth Wood (kw25@st-andrews.ac.uk)
*                    Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
*
* python_finite_volume_solver is free software: you can redistribute it and/or
* modify it under the terms of the GNU Affero General Public License as
* published by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* python_finite_volume_solver is distributed in the hope that it will be useful,
* but WITOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with python_finite_volume_solver. If not, see
* <http://www.gnu.org/licenses/>.
********************************************************************************

********************************************************************************
* @file riemansolver.f
*
* @brief Stanalone Riemann solver library: Fortran 77 version.
*
* This Riemann solver is the original exact Riemann solver as it is presented in
* Toro, E., Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd
* edition (Springer, 2009), chapter 4.
*
* It was carefully copied (by hand) from that book to actual Fortran 77 code by
* Kenneth Wood, and slightly modified to be used as a Riemann solver routine in
* a 1-3D finite volume code.
*
* Additional changes (including this documentation) were made by Bert
* Vandenbroucke to make the Python, C++ and Fortran version of the Riemann
* solver more homogeneous.
*
* @author Kenneth Wood (kw25@st-andrews.ac.uk)
* @author Bert Vandenbroucke (bv7@st-andrews.ac.uk)
********************************************************************************

********************************************************************************
* @brief Solve the Riemann problem with the given left and right state.
*
* @param d1 Left state density.
* @param u1 Left state fluid velocity.
* @param p1 Left state pressure.
* @param d2 Right state density.
* @param u2 Right state fluid velocity.
* @param p2 Right state pressure.
* @param gg Adiabatic index of the gas.
* @param ds Variable to store the resulting density in.
* @param us Variable to store the resulting fluid velocity in.
* @param ps Variable to store the resulting pressure in.
* @param uflag Flag variable used to signal if the solution was sampled from the
* left state (-1) or the right state (1). A vacuum state (0) is currently not
* supported by this version of the Riemann solver.
* @param x_over_t Point (in velocity space) where we want to sample the
* solution. In a finite volume scheme you usually want to set this value to 0.
********************************************************************************
      subroutine riemann(d1, u1, p1, d2, u2, p2, gg, ds,
     &                   us, ps, uflag, x_over_t)

      implicit none

C     Declaration of variables:
      integer I, CELLS, uflag
      REAL(8) GAMMA, G1, G2, G3, G4, G5, G6, G7, G8,
     &     DL, UL, PL, CL, DR, UR, PR, CR,
     &     DIAPH, DOMLEN, DS, DX, PM, MPA, PS, S,
     &     TIMEOUT, UM, US, XPOS

      real(8) gg, d1, u1, p1, d2, u2, p2
      real(8) x_over_t, dxdt

      COMMON /GAMMAS/ GAMMA, G1, G2, G3, G4, G5, G6, G7, G8
      COMMON /STATES/ DL, UL, PL, CL, DR, UR, PR, CR

      dxdt = x_over_t

C     Compute gamma related constants

      GAMMA = gg
      DL = d1
      UL = u1
      PL = p1
      DR = d2
      UR = u2
      PR = p2

      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

      MPA=1.0  ! Not sure about this????

C     Compute sound speeds

      CL = SQRT(GAMMA*PL/DL)
      CR = SQRT(GAMMA*PR/DR)

C     The pressure positivity condition is tested for

      IF(G4*(CL+CR).LE.(UR-UL))THEN
C     The initial data is such that vacuum is generated. Program stopped
        WRITE(6,*) 'Vacuum generated by data. Program stopped.'
        STOP
      ENDIF

C     Exact solution for pressure and velocity in star region is found

      CALL STARPU(PM, UM, MPA)
      if(pm.ne.pm) then
        print*,pm,pl,dl,ul,pr,dr,ur
        stop
      endif
      if(UM .gt.dxdt) then
        uflag=-1 ! "left state"
      else
        uflag=1 ! "right state
      endif


      CALL SAMPLE(PM, UM, dxdt, DS, US, PS)

      RETURN
      END

********************************************************************************
* @brief Get the pressure and velocity in the middle state.
*
* @param p Variable to store the resulting pressure in.
* @param u Variable to store the resulting fluid velocity in.
* @param mpa Mach number? Not used.
********************************************************************************
      SUBROUTINE STARPU(P, U, MPA)

      IMPLICIT NONE

C     Purpose: to compute the solution for pressure and velocity in star region

      INTEGER I, NRITER

      REAL(8) DL, UL, PL, CL, DR, UR, PR, CR,
     &     CHANGE, FL, FLD, FR, FRD, P, POLD, PSTART,
     &     TOLPRE, U, UDIFF, MPA

      COMMON/STATES/ DL, UL, PL, CL, DR, UR, PR, CR
      DATA TOLPRE, NRITER/1.0d-06, 20/

C     Guessed value PSTART is computed

      CALL GUESSP(PSTART)

      POLD = PSTART
      UDIFF = UR - UL

c      WRITE(6,*)'Iteration number Change '

      DO 10 I = 1, NRITER

         CALL PREFUN(FL, FLD, POLD, DL, PL, CL)
         CALL PREFUN(FR, FRD, POLD, DR, PR, CR)
         P = POLD - (FL + FR + UDIFF)/(FLD + FRD)
         CHANGE = 2.0*ABS((P - POLD)/(P + POLD))
c         WRITE(6, 30)I, CHANGE
         IF(CHANGE.LE.TOLPRE)GOTO 20
         IF(P.LT.0.0)P = TOLPRE
         POLD = P

 10   CONTINUE

c      WRITE(6,*)'Divergence in Newton-Raphson iteration'

 20   CONTINUE

C     Compute velocity in star region

      U = 0.5*(UL + UR + FR - FL)

c      WRITE(6,*)'Pressure	Velocity'
c      WRITE(6,40)P/MPA, U

 30   FORMAT(5X, I5,15X, F12.7)
 40   FORMAT(2(F14.6, 5X))

      RETURN
      END

********************************************************************************
* @brief Get an initial guess for the pressure in the middle state, as a
* starting point for the Newton-Raphson iteration.
*
* @param pm Variable to store the resulting pressure guess in.
********************************************************************************
      SUBROUTINE GUESSP(PM)

C     Purpose: to provide a guess value for pressure PM in the Star Region.
C     The choice is made according to adaptive Riemann solver using the PVRS,
C     TRRS, and TSRS approximate Riemann solvers

      IMPLICIT NONE

      REAL(8) DL, UL, PL, CL, DR, UR, PR, CR,
     &     GAMMA, G1, G2, G3, G4, G5, G6, G7, G8,
     &     CUP, GEL, GER, PM, PMAX, PMIN, PPV, PQ,
     &     PTL, PTR, QMAX, QUSER, UM

      COMMON /GAMMAS/ GAMMA, G1, G2, G3, G4, G5, G6, G7, G8
      COMMON /STATES/ DL, UL, PL, CL, DR, UR, PR, CR

      QUSER = 2.0

C     Compute guess pressure for PVRS Riemann solver

      CUP = 0.25*(DL + DR)*(CL + CR)
      PPV = 0.5*(PL + PR) + 0.5*(UL - UR)*CUP
      PPV = AMAX1(0.0, PPV)
      PMIN = AMIN1(PL, PR)
      PMAX = AMAX1(PL, PR)
      QMAX = PMAX/PMIN

      IF(QMAX.LE.QUSER.AND.
     &   (PMIN.LE.PPV.AND.PPV.LE.PMAX))THEN

C        Select PVRS Riemann solver

         PM = PPV

      ELSE
         IF(PPV.LT.PMIN)THEN

C           Select Two-Rarefaction Riemann solver

            PQ = (PL/PR)**G1
            UM = (PQ*UL/CL + UR/CR +
     &            G4*(PQ - 1.0))/(PQ/CL + 1.0/CR)
            PTL = 1.0 + G7*(UL - UM)/CL
            PTR = 1.0 + G7*(UM - UR)/CR
            PM = 0.5*(PL*PTL**G3 + PR*PTR**G3)
         ELSE

c           Select Two=Shock Riemann solver with PVRS as estimate

            GEL = SQRT((G5/DL)/(G6*PL + PPV))
            GER = SQRT((G5/DR)/(G6*PR + PPV))
            PM = (GEL*PL + GER*PR - (UR - UL))/(GEL + GER)
         ENDIF
      ENDIF

      RETURN
      END

********************************************************************************
* @brief Evaluate the pressure functions for the given pressure guess and the
* given left or right state variables.
*
* @param f Variable to store the pressure function value in.
* @param fd Variable to store the value of the derivative of the presure
* function in.
* @param p Current pressure guess.
* @param dk Left or right state density.
* @param pk Left or right state pressure.
* @param ck Left or right state sound speed.
********************************************************************************
      SUBROUTINE PREFUN(F, FD, P, DK, PK, CK)

C     Purpose: to evaluate the pressure functions FL and FR in the exact Riemann
C     solver

      IMPLICIT NONE

      REAL(8) AK, BK, CK, DK, F, FD, P, PK, PRAT, QRT,
     &     GAMMA, G1, G2, G3, G4, G5, G6, G7, G8

      COMMON /GAMMAS/ GAMMA, G1, G2, G3, G4, G5, G6, G7, G8

      IF(P.LE.PK) THEN

C        Rarefaction wave

         PRAT = P/PK
         F = G4*CK*(PRAT**G1 - 1.0)
         FD = (1.0/(DK*CK))*PRAT**(-G2)
      ELSE

C        Shock wave

         AK = G5/DK
         BK = G6*PK
         QRT = SQRT(AK/(BK + P))
         F = (P - PK)*QRT
         FD = (1.0 - 0.5*(P - PK)/(BK + P))*QRT
      ENDIF

      RETURN
      END

********************************************************************************
* @brief Sample the Riemann solution with the given middle state pressure and
* fluid velocity at the given point in velocity space.
*
* @param pm Middle state pressure.
* @param um Middle state fluid velocity.
* @param s Sampling point in velocity space.
* @param d Variable to store the resulting density in.
* @param u Variable to store the resulting fluid velocity in.
* @param p Variable to store the resulting pressure in.
********************************************************************************
      SUBROUTINE SAMPLE(PM, UM, S, D, U, P)

C     Purpose to sample the solution throughout the wave pattern. Pressure PM
C     and velocity UM in the Star Region are known. Sampling is performed in
C     terms of the 'speed' S = X/T. Sampled values are D, U, P

C     Input variables: PM, UM, S, /GAMMAS. /STATES/
C     Output variables: D, U, P

      IMPLICIT NONE

      REAL(8) DL, UL, PL, CL, DR, UR, PR, CR,
     &     GAMMA, G1, G2, G3, G4, G5, G6, G7, G8,
     &     C, CML, CMR, D, P, PM, PML, PMR, S,
     &     SHL, SHR, SL, SR, STL, STR, U, UM

      COMMON /GAMMAS/ GAMMA, G1, G2, G3, G4, G5, G6, G7, G8
      COMMON /STATES/ DL, UL, PL, CL, DR, UR, PR, CR

      IF(S.LE.UM)THEN

C     Sampling point lies to the left of the contact discontinuity

         IF(PM.LE.PL)THEN

C           Left rarefaction

            SHL = UL - CL

            IF(S.LE.SHL)THEN

C              Sampled point is left data state

               D = DL
               U = UL
               P = PL
            ELSE
               CML = CL*(PM/PL)**G1
               STL = UM - CML

               IF(S.GT.STL)THEN

C                 Sampled point is Star Left state

                  D = DL*(PM/PL)**(1.0/GAMMA)
                  U = UM
                  P = PM
               ELSE

C                 Sampled point is inside left fan

                  U = G5*(CL + G7*UL + S)
                  C = G5*(CL + G7*(UL - S))
                  D = DL*(C/CL)**G4
                  P = PL*(C/CL)**G3
               ENDIF
            ENDIF
c         ENDIF
         ELSE

C            Left shock

            PML = PM/PL
            SL = UL - CL*SQRT(G2*PML + G1)

            IF(S.LE.SL)THEN

C              Sampled point is left data state

               D = DL
               U = UL
               P = PL

             ELSE

C              Sampled point is Star Left state

               D = DL*(PML+G6)/(PML*G6 + 1.0)
               U = UM
               P = PM
             ENDIF
           ENDIF
          ELSE

C     Sampling point lies to the right of the contact discontinuity

         IF(PM.GT.PR)THEN

C           Right shock

            PMR = PM/PR
            SR = UR + CR*SQRT(G2*PMR + G1)

            IF(S.GE.SR)THEN

C              Sampled point is right data state

               D = DR
               U = UR
               P = PR
            ELSE

C              Sampled point is Star Right state

               D = DR*(PMR + G6)/(PMR*G6 + 1.0)
               U = UM
               P = PM
            ENDIF
         ELSE

C           Right rarefaction

            SHR = UR + CR

            IF(S.GE.SHR)THEN

C              Sampled point is right data state

               D = DR
               U = UR
               P = PR
            ELSE
               CMR = CR*(PM/PR)**G1
               STR = UM + CMR

               IF(S.LE.STR)THEN

C                 Sampled point is Star Right state

                  D = DR*(PM/PR)**(1.0/GAMMA)
                  U = UM
                  P = PM
               ELSE

C                 Sampled point is inside left fan

                  U=G5*(-CR + G7*UR + S)
                  C = G5*(CR - G7*(UR-S))
                  D = DR*(C/CR)**G4
                  P = PR*(C/CR)**G3
               ENDIF
            ENDIF
         ENDIF
      ENDIF

      RETURN
      END