Godunov Solver for Euler Equations​

This program solves the Euler equations in 1, 2, or 3 dimensions using Godunov's first-order upwind method. (The attachment linked in this header contains source code and makefile for all three of the programs on this page.) Input is taken from a file:

NAMELIST/ProblemData/ starttime,finaltime,frames,xL,xR,xcells,yL,yR,ycells,zL,zR,zcells,BC,IC,CFL,cov,gamma,tol,maxiter,RiemannSolver,Q,alpha,beta,advectionConst,Eq,nDim,num_threads,lthreading
        
! Make sure problem.data exists. If it does, read in. If not, output error and quit.
inquire(file='problem.data', exist=ex)
if (ex) then
 OPEN(UNIT=30, FILE='problem.data', STATUS="OLD")
 READ(30,NML=ProblemData)                               ! read into variables in namelist
 CLOSE(30)
else
 print *, 'Missing problem.data. Exiting.'
 stop           ! exit if missing
end if

The domain is initialized based on selected test (options depend on number of dimensions):

  select case (nDim)
   
   ! 1D case
   case (1)
    ! Allocate domain arrays for conserved variables and intercell flux
    allocate(ProblemRegion%U(0:xcells+1,NrVar),ProblemRegion%flux(0:xcells,NrVar))
    
    select case (IC)
     case (1)
      WL(irho) = 1.0
      WR(irho) = 0.125
      WL(iu) = 0.75
      WR(iu) = 0.0
      WL(ip) = 1.0
      WR(ip) = 0.1
      x0(ix) = 0.3
     case (2)
      WL(irho) = 1.0
      WR(irho) = 1.0
      WL(iu) = -2.0
      WR(iu) = 2.0
      WL(ip) = 0.4
      WR(ip) = 0.4
      x0(ix) = 0.5
     case (3)
      WL(irho) = 1.0
      WR(irho) = 1.0
      WL(iu) = 0.0
      WR(iu) = 0.0
      WL(ip) = 1000.0
      WR(ip) = 0.01
      x0(ix) = 0.5
     case (4)
      WL(irho) = 5.99924
      WR(irho) = 5.99242
      WL(iu) = 19.975
      WR(iu) = -6.19633
      WL(ip) = 460.894
      WR(ip) = 46.0950
      x0(ix) = 0.4
     case (5)
      WL(irho) = 1.0
      WR(irho) = 1.0
      WL(iu) = -19.59745
      WR(iu) = -19.59745
      WL(ip) = 1000.0
      WR(ip) = 0.01
      x0(ix) = 0.8
     case (6)
      WL(irho) = 1.4
      WR(irho) = 1.0
      WL(iu) = 0.0
      WR(iu) = 0.0
      WL(ip) = 1.0
      WR(ip) = 1.0
      x0(ix) = 0.5
     case (7)
      WL(irho) = 1.4
      WR(irho) = 1.0
      WL(iu) = 0.1
      WR(iu) = 0.1
      WL(ip) = 1.0
      WR(ip) = 1.0
      x0(ix) = 0.5
    end select

    ! Initialize domain using set initial conditions (conserved variables)
    do i = 1, xcells
     x(ix) = (i-0.5)*dx(ix) + xL
     if (x(ix) <= x0(ix)) then
      ProblemRegion%U(i,:) = PrimToCons(WL)

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
PrimToCons(irho) = W(irho)
PrimToCons(imom(ix)) = W(iu(ix))*W(irho)
PrimToCons(iE) = W(ip)*(1. - cov*W(irho))/(gamma - 1.) + 0.5*W(irho)*sum(W(iu)**2)
! Tangential directins
if (nDim > 1) PrimToCons(imom(iy)) = W(iu(iy))*W(irho)
if (nDim > 2) PrimToCons(imom(iz)) = W(iu(iz))*W(irho)
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

     else
      ProblemRegion%U(i,:) = PrimToCons(WR)
     end if
    end do
   
   ! 2D case
   case (2)
    allocate(ProblemRegion2D%U(0:xcells+1,0:ycells+1,NrVar),ProblemRegion%U(0:max(xcells,ycells)+1,NrVar),ProblemRegion%flux(0:max(xcells,ycells),NrVar))
  
    select case (IC)
     case (explosion)
      ! "Left" is inside circular boundary, "Right" is outside circular boundary
      WL(irho) = 1.0
      WL(iu(ix)) = 0.0
      WL(iu(iy)) = 0.0
      WL(ip) = 1.0
      WR(irho) = 0.125
      WR(iu(ix)) = 0.0
      WR(iu(iy)) = 0.0
      WR(ip) = 0.1
      x0 = (/ 0, 0 /)
      radius = 0.4
      do i = 1, xcells
       x(ix) = (i-0.5)*dx(ix) + xL
       do j = 1, ycells
        x(iy) = (j-0.5)*dx(iy) + yL
        if (sum((x-x0)**2) < radius**2) then
         ProblemRegion2D%U(i,j,:) = PrimToCons(WL)
         print *, ProblemRegion2D%U(i,j,:)
        else
         ProblemRegion2D%U(i,j,:) = PrimToCons(WR)
         print *, ProblemRegion2D%U(i,j,:)
        end if
       end do
      end do
    end select
   
   ! 3D case  
   case (3)
    allocate(ProblemRegion3D%U(0:xcells+1,0:ycells+1,0:zcells+1,NrVar),ProblemRegion%U(0:max(xcells,ycells,zcells)+1,NrVar),ProblemRegion%flux(0:max(xcells,ycells,zcells),NrVar))
    
    select case (IC)
     case (explosion)
      ! "Left" is inside spherical boundary, "Right" is outside spherical boundary
      WL(irho) = 1.0
      WL(iu(ix)) = 0.0
      WL(iu(iy)) = 0.0
      WL(iu(iz)) = 0.0
      WL(ip) = 1.0
      WR(irho) = 0.125
      WR(iu(ix)) = 0.0
      WR(iu(iy)) = 0.0
      WR(iu(iz)) = 0.0
      WR(ip) = 0.1
      x0 = (/ 0, 0, 0 /)
      radius = 0.4
      do i = 1, xcells
       x(ix) = (i-0.5)*dx(ix) + xL
       do j = 1, ycells
        x(iy) = (j-0.5)*dx(iy) + yL
        do k = 1, zcells
         x(iz) = (k-0.5)*dx(iz) + zL
         if (sum((x-x0)**2) < radius**2) then
          ProblemRegion3D%U(i,j,k,:) = PrimToCons(WL)
         else
          ProblemRegion3D%U(i,j,k,:) = PrimToCons(WR)
         end if
        end do
       end do
      end do
    end select
  
  end select

After the domain is initialized, the first frame is written to file by

 ! Euler solver with Godunov method
 case (GodunovEuler) 
  if (nDim == 1) then
   write (40,*) 'x rho u p'       ! write column headers
   do i = 1, xcells
    x(ix) = xL + (i-0.5)*dx(ix)    ! Data is located at center of each cell
    U = ProblemRegion%U(i,:)
    write (40,*) x, ConsToPrim(U)

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
ConsToPrim(irho) = U(irho)
ConsToPrim(iu(ix)) = U(imom(ix))/U(irho)
ConsToPrim(ip) = (U(iE) - (sum(U(imom)**2))/(2*U(irho)))*(gamma - 1.)/(1. - cov*U(irho))
! Tangential directions
if (nDim > 1) ConsToPrim(iu(iy)) = U(imom(iy))/U(irho)
if (nDim > 2) ConsToPrim(iu(iz)) = U(imom(iz))/U(irho)
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

   end do
  else if (nDim == 2) then
   write (40,*) 'x y rho u v p'
   do i = 1, xcells
    x(ix) = xL + (i-0.5)*dx(ix)    ! Data is located at center of each cell
    do j = 1, ycells
     x(iy) = yL + (j-0.5)*dx(iy)    ! Data is located at center of each cell
     U = ProblemRegion2D%U(i,j,:)
     write (40,*) x, ConsToPrim(U)
    end do
   end do
  else if (nDim == 3) then
   write (40,*) 'x y z rho u v w p'
   do i = 1, xcells
    x(1) = xL + (i-0.5)*dx(1)    ! Data is located at center of each cell
    do j = 1, ycells
     x(2) = yL + (j-0.5)*dx(2)    ! Data is located at center of each cell
     do k = 1, zcells
      x(3) = zL + (k-0.5)*dx(3)    ! Data is located at center of each cell
      U = ProblemRegion3D%U(i,j,k,:)
      write (40,*) x, ConsToPrim(U)
     end do
    end do
   end do
  end if

The time-stepping procedure is then entered. The \Delta t for the current step is calculated by determining the maximum speed in any direction on the grid:

  if (nDim == 1) then
   do i = 1, xcells
    U = ProblemRegion%U(i,:)
    cell = region(ConsToPrim(U))
    if (cell%c .ne. cell%c) then
     print *, 'NaN in sound speed at time ', time, ' in cell ', i,' in TimeStep.'
     stop
    end if
    if (abs(cell%W(iu(ix))) + cell%c > maxSpeed) maxSpeed = abs(cell%W(iu(ix))) + cell%c
   end do
  ! 2D
  else if (nDim == 2) then
   do i = 1, xcells
    do j = 1, ycells
     U = ProblemRegion2D%U(i,j,:)
     cell = region(ConsToPrim(U))
     if (cell%c .ne. cell%c) then
      print *, 'NaN in sound speed at time ', time, ' in cell ', i,' in TimeStep.'
      stop
     end if
     if (abs(cell%W(iu(ix))) + cell%c > maxSpeed) maxSpeed = abs(cell%W(iu(ix))) + cell%c
     if (abs(cell%W(iu(iy))) + cell%c > maxSpeed) maxSpeed = abs(cell%W(iu(iy))) + cell%c
    end do
   end do
  ! 3D
  else if (nDim == 3) then
   do i = 1, xcells
    do j = 1, ycells
     do k = 1, zcells
      U = ProblemRegion3D%U(i,j,k,:)
      cell = region(ConsToPrim(U))
      if (cell%c .ne. cell%c) then
       print *, 'NaN in sound speed at time ', time, ' in cell ', i,' in TimeStep.'
       stop
      end if
      if (abs(cell%W(iu(ix))) + cell%c > maxSpeed) maxSpeed = abs(cell%W(iu(ix))) + cell%c
      if (abs(cell%W(iu(iy))) + cell%c > maxSpeed) maxSpeed = abs(cell%W(iu(iy))) + cell%c
      if (abs(cell%W(iu(iz))) + cell%c > maxSpeed) maxSpeed = abs(cell%W(iu(iz))) + cell%c
     end do
    end do
   end do
  end if

The time-stepping procedure then proceeds in a slightly different manner, depending upon the number of dimensions in the problem. In the 1D case, boundary conditions are applied:

select case (BC)
 case(periodic) ! Wraps around
  ProblemRegion%U(0,1:NrVar) = ProblemRegion%U(cells,1:NrVar)
  ProblemRegion%U(cells+1,1:NrVar) = ProblemRegion%U(1,1:NrVar)
 case(transparent)      ! Wave flows out
  ProblemRegion%U(0,1:NrVar) = ProblemRegion%U(1,1:NrVar)
  ProblemRegion%U(cells+1,1:NrVar) = ProblemRegion%U(cells,1:NrVar)
 case(reflective)       ! Wave bounces off hard boundary
  ProblemRegion%U(cells+1,irho) = ProblemRegion%U(cells,irho)
  ProblemRegion%U(cells+1,imom) = ProblemRegion%U(cells,imom)
  ProblemRegion%U(cells+1,imom(ix)) = -ProblemRegion%U(cells,imom(ix)) ! Momentum normal to boundary is reversed
  ProblemRegion%U(cells+1,iE) = ProblemRegion%U(cells,iE)
  ProblemRegion%U(0,irho) = ProblemRegion%U(1,irho)
  ProblemRegion%U(0,imom) = ProblemRegion%U(1,imom)
  ProblemRegion%U(0,imom(ix)) = -ProblemRegion%U(1,imom(ix))
  ProblemRegion%U(0,iE) = ProblemRegion%U(1,iE)
end select

and the flux across cell boundaries is then calculated. Depending on the method selected, the flux may be estimated directly, or the fluid state may be approximated, which can then be used to calculate the flux. Summaries of the various methods here can be found near the end of this page.

! Loop through domain using equation set in problem.data
  do i = 0, cells
   
   U = ProblemRegion%U(i,:)
   L = region(ConsToPrim(U))
   if (L%c .ne. L%c) then
    print *, 'NaN in sound speed at time ', time, ' in cell ', i,' in CalcFlux.'
    stop
   end if
   U = ProblemRegion%U(i+1,:)
   R = region(ConsToPrim(U))
   
   ! Some solvers approximate flux rather than fluid state
   select case (RiemannSolver)
    case (HLLC)
     call WaveSpeed(L,R,SL,SR,S0)

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
p0 = max(0d0, PrimitiveVariablePressure(L,R))
! Approximate shock or rarefaction speeds
if (p0 > L%W(ip)) then
 SL = L%W(iu(ix)) - L%c*sqrt(1 + (gamma + 1.)*(p0/L%W(ip) - 1.)/(2.*gamma))
else
 SL = L%W(iu(ix)) - L%c
end if
if (p0 > R%W(ip)) then
 SR = R%W(iu(ix)) + R%c*sqrt(1 + (gamma + 1.)*(p0/R%W(ip) - 1.)/(2.*gamma))
else
 SR = R%W(iu(ix)) + R%c
end if
if (RiemannSolver == HLLC) S0 = (R%W(ip) - L%W(ip) + L%W(irho)*L%W(iu(ix))*(SL - L%W(iu(ix))) - R%W(irho)*R%W(iu(ix))*(SR - R%W(iu(ix))))/(L%W(irho)*(SL - L%W(iu(ix))) - R%W(irho)*(SR - R%W(iu(ix))))
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

     ProblemRegion%flux(i,:) = HLLCFlux(L,R,S0,SL,SR)

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
! Approximate with flux from left Riemann state
if (0 <= SL) then
 HLLCflux = EulerFlux(L%W)
! Approximate with integral average of left star region
else if (SL <= 0 .and. 0 <= S0) then
 UL = PrimToCons(L%W)
 U0L(irho) = L%W(irho)*(SL - L%W(iu(ix)))/(SL - S0)
 U0L(imom(ix)) = S0*L%W(irho)*(SL - L%W(iu(ix)))/(SL - S0)
 if (nDim > 1) U0L(imom(iy)) = L%W(iu(iy))*L%W(irho)*(SL - L%W(iu(ix)))/(SL - S0)
 if (nDim > 2) U0L(imom(iz)) = L%W(iu(iz))*L%W(irho)*(SL - L%W(iu(ix)))/(SL - S0)
 U0L(iE) = L%W(irho)*((SL - L%W(iu(ix)))/(SL - S0))*(UL(iE)/L%W(irho) + (S0 - L%W(iu(ix)))*(S0 + L%W(ip)/(L%W(irho)*(SL - L%W(iu(ix))))))
 HLLCflux = EulerFlux(L%W) + SL*(U0L - UL)
! Approximate with integral average of right star region
else if (S0 <= 0 .and. 0 <= SR) then
 UR = PrimToCons(R%W)
 U0R(irho) = R%W(irho)*(SR - R%W(iu(ix)))/(SR - S0)
 U0R(imom(ix)) = S0*R%W(irho)*(SR - R%W(iu(ix)))/(SR - S0)
 if (nDim > 1) U0R(imom(iy)) = R%W(iu(iy))*R%W(irho)*(SR - R%W(iu(ix)))/(SR - S0)
 if (nDim > 2) U0R(imom(iz)) = R%W(iu(iz))*R%W(irho)*(SR - R%W(iu(ix)))/(SR - S0)
 U0R(iE) = R%W(irho)*((SR - R%W(iu(ix)))/(SR - S0))*(UR(iE)/R%W(irho) + (S0 - R%W(iu(ix)))*(S0 + R%W(ip)/(R%W(irho)*(SR - R%W(iu(ix))))))
 HLLCflux = EulerFlux(R%W) + SR*(U0R - UR)
! Approximate with right Riemann state
else
 HLLCflux = EulerFlux(R%W)
end if
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

    case (HLL)
     call WaveSpeed(L,R,SL,SR)
     ProblemRegion%flux(i,:) = HLLFlux(L,R,SL,SR)

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
! Approximate with left Riemann state
if (0 <= SL) then
 HLLflux = EulerFlux(L%W)
! Approximate with integral average of star region
else if (SL <= 0 .and. 0 <= SR) then
 UL = PrimToCons(L%W)
 UR = PrimToCons(R%W)
 HLLflux = (SR*EulerFlux(L%W) - SL*EulerFlux(R%W) + SL*SR*(UR - UL))/(SR - SL)
! Approximate with right Riemann state
else
 HLLflux = EulerFlux(R%W)
end if
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

    ! If solver approximates fluid state, use it to calculate flux
    case default
     ProblemRegion%flux(i,:) = EulerFlux(RiemannSoln(L,R))

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
function EulerRiemannSoln

! Check for vacuum
call check_vacuum(L,R)

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
critdiff = 2.*L%c/(gamma - 1.) + 2.*R%c/(gamma - 1.)      ! Difference between right and left velocities must be strictly less than this for pressure to remain positive

if ((L%W(ip) == 0 .and. L%W(irho) == 0) .or. (R%W(ip) == 0 .and. R%w(irho) == 0)) then  ! one of the initial states is vacuum
 vacuum = 1
 print *, 'One of initial states is vacuum.'
else if (critdiff <= R%W(iu(ix)) - L%W(iu(ix))) then     ! vacuum is created by initial states
 vacuum = 2
 print *, 'Vacuum is created by initial states.'
else    ! no vacuum
 vacuum = 3
end if
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

! Vacuum present, don't need star region - exact solution is cheap
if (vacuum /= 3) then
 EulerRiemannSoln = sample(L,R,0d0)

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
! Calculate shock and rarefaction speeds
SvacL = L%W(iu(ix)) + 2.*L%c/(gamma - 1.)   ! Speed of vacuum front (vacuum on right)
SvacR = R%W(iu(ix)) - 2.*R%c/(gamma - 1.)   ! Speed of vacuum front (vacuum on left)

! Initial state is vacuum
if (vacuum == 1) then
 if (R%W(ip) == 0 .and. R%W(irho) == 0) then !vacuum on the right
  lright = .false.
  R%W(iu(ix)) = 0
  if (location <= L%W(iu(ix)) - L%c) then     ! Left of rarefaction fan
   sampleWithVacuum = L%W
  else if (location >= SvacL) then! Right of rarefaction fan (in vacuum)
   sampleWithVacuum = R%W
  else       ! In rarefaction fan
   call rarefaction(L,R,Wfan,lright,location)

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
! If fluid has covolume constant, iterate to solve for density
if (cov > 0) then
 if (.not. lright) then
  beta = ((((gamma - 1.)*(L%W(iu(ix)) + 2.*L%c*(1. - cov*L%W(irho))/(gamma - 1.) - location))**2.)/(gamma*L%W(ip)*((1. - cov*L%W(irho))/L%W(irho))**gamma)) ! Constant related to Riemann invariant
  Wfan(irho) = NRSolver(rhoguess(L,R),densityF,densityDf,L,R)
  Wfan(ip) = L%W(ip)*(((1. - cov*L%W(irho))/L%W(irho))*(Wfan(irho)/(1. - cov*Wfan(irho)))**gamma)
  c = sqrt((Wfan(ip)*gamma)/(Wfan(irho)*(1. - cov*Wfan(irho)))) ! Speed of sound at current location
  Wfan(iu(ix)) = location + c
  if (nDim > 1) Wfan(iu(iy)) = L%W(iu(iy))
  if (nDim > 2) Wfan(iu(iz)) = L%W(iu(iz))
 else
  beta = (((location - R%W(iu(ix)) + (2.*R%c*(1. - cov*R%W(irho))/(gamma - 1.)))**2.)/(gamma*R%W(ip)*((1. - cov*R%W(irho))/R%W(irho))**gamma))   ! Constant related to Riemann invariant
  Wfan(irho) = NRSolver(rhoguess(L,R),densityF,densityDf,L,R)
  Wfan(ip) = R%W(ip)*(((1. - cov*R%W(irho))/R%W(irho))*(Wfan(irho)/(1. - cov*Wfan(irho)))**gamma)
  c = sqrt((Wfan(ip)*gamma)/(Wfan(irho)*(1. - cov*Wfan(irho)))) ! Speed of sound at current location
  Wfan(iu(ix)) = location - c
  if (nDim > 1) Wfan(iu(iy)) = R%W(iu(iy))
  if (nDim > 2) Wfan(iu(iz)) = R%W(iu(iz))
 end if
! Otherwise solution is closed-form
else
 ! Rarefaction fan differs depending on whether it's a right or left rarefaction
 if (.not. lright) then
  Wfan(irho) = L%W(irho)*(2/(gamma + 1) + (gamma - 1)*(L%W(iu(ix)) - location)/((gamma + 1)*L%c))**(2/(gamma - 1))
  Wfan(iu(ix)) = (2/(gamma + 1))*(L%c + (gamma - 1)*L%W(iu(ix))/2 + location)
  Wfan(ip) = L%W(ip)*(2/(gamma + 1) + (gamma - 1)*(L%W(iu(ix)) - location)/((gamma + 1)*L%c))**(2*gamma/(gamma - 1))
  ! Tangential directions unchanged
  if (nDim > 1) Wfan(iu(iy)) = L%W(iu(iy))
  if (nDim > 2) Wfan(iu(iz)) = L%W(iu(iz))
 else
  Wfan(irho) = R%W(irho)*(2/(gamma + 1) - (gamma - 1)*(R%W(iu(ix)) - location)/((gamma + 1)*R%c))**(2/(gamma - 1))
  Wfan(iu(ix)) = (2/(gamma + 1))*(-R%c + (gamma - 1)*R%W(iu(ix))/2 + location)
  Wfan(ip) = R%W(ip)*(2/(gamma + 1) - (gamma - 1)*(R%W(iu(ix)) - location)/((gamma + 1)*R%c))**(2*gamma/(gamma - 1))
  ! Tangential directions unchanged
  if (nDim > 1) Wfan(iu(iy)) = R%W(iu(iy))
  if (nDim > 2) Wfan(iu(iz)) = R%W(iu(iz))
 end if
end if
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

   sampleWithVacuum = Wfan
  end if
 else!vacuum on the left
  lright = .true.
  L%W(iu(ix)) = 0
  if (location >= R%W(iu(ix)) - R%c) then     ! Right of rarefaction fan
   sampleWithVacuum = R%W
  else if (location <= SvacR) then! Left of rarefaction fan (in vacuum)
   sampleWithVacuum = L%W
  else
   call rarefaction(L,R,Wfan,lright,location)       ! In rarefaction fan
   sampleWithVacuum = Wfan
  end if
 end if
! Initial regions result in vacuum between
else if (vacuum == 2) then
 if (location <= L%W(iu(ix)) - L%c) then      ! Left of left rarefaction fan
  sampleWithVacuum = L%W
 else if (location <= SvacL) then ! In left rarefaction fan
  lright = .false.
  call rarefaction(L,R,Wfan,lright,location)
  sampleWithVacuum = Wfan
 else if (location >= R%W(iu(ix)) + R%c) then ! Right of right rarefaction fan
  sampleWithVacuum = R%W
 else if (location >= SvacR) then ! In right rarefaction fan
  lright = .true.
  call rarefaction(L,R,Wfan,lright,location)
  sampleWithVacuum = Wfan
 else! In vacuum (between rarefaction fans)
  sampleWithVacuum(irho) = 0
  sampleWithVacuum(iu(ix)) = 0
  sampleWithVacuum(ip) = 0
  if (nDim > 1) sampleWithVacuum(iu(iy)) = 0
  if (nDim > 2) sampleWithVacuum(iu(iz)) = 0
 endif
end if
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

! No vacuum, calculate star region and sample along t axis
else
 ! Star region may be approximated
 call StarRegion(L,R,Lstar,Rstar)

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
! Passively advected - no math, so no approximations
if (nDim > 1) then
 Lstar%W(iu(iy)) = L%W(iu(iy))
 Rstar%W(iu(iy)) = R%W(iu(iy))
end if
if (nDim > 2) then
 Lstar%W(iu(iz)) = L%W(iu(iz))
 Rstar%W(iu(iz)) = R%W(iu(iz))
end if

select case (RiemannSolver)
 ! Exact Riemann solver
 case (exact)
  ! Use Newton-Raphson solver to find pressure in star region
  Lstar%W(ip) = NRSolver(pguess(L,R),covolumeF,covolumeDf,L,R)
  ! Calculate velocity and densities in star regions
  Lstar%W(iu(ix)) = exactVelocity(Lstar%W(ip),L,R)
  Rstar%W(ip) = Lstar%W(ip)
  Rstar%W(iu(ix)) = Rstar%W(iu(ix))
  Lstar%W(irho) = exactDensity(Lstar,L)
  Rstar%W(irho) = exactDensity(Rstar,R)
!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
!            See below for exact Riemann solver summary           !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 ! Two-shock approximation
 case (TSRS)
  p0 = max(PrimitiveVariablePressure(L,R),0d0)
  Lstar%W(ip) = TwoShockPressure(L,R,p0)
  Lstar%W(iu(ix)) = TwoShockVelocity(L,R,Lstar,p0)
  Rstar%W(ip) = Lstar%W(ip)
  Rstar%W(iu(ix)) = Rstar%W(iu(ix))
  Lstar%W(irho) = TwoShockDensity(L,Lstar)
  Rstar%W(irho) = TwoShockDensity(R,Rstar)
 ! Adaptive non-iterative method
 case (ANRS)
  p0 = PrimitiveVariablePressure(L,R)
  pmax = max(L%W(ip),R%W(ip))
  pmin = min(L%W(ip),R%W(ip))
  if (pmax/pmin < Q .and. pmin < p0 .and. pmax > p0) then
   Lstar%W(ip) = p0
   Lstar%W(iu(ix)) = PrimitiveVariableVelocity(L,R)
   Rstar%W(ip) = Lstar%W(ip)
   Rstar%W(iu(ix)) = Rstar%W(iu(ix))
   Lstar%W(irho) = PrimitiveVariableDensity(L,Lstar,R)
   Rstar%W(irho) = PrimitiveVariableDensity(R,Rstar,L)
  else if (p0 < pmin) then
   Lstar%W(ip) = TwoRarefactionPressure(L,R)
   Lstar%W(iu(ix)) = TwoRarefactionVelocity(L,R)
   Rstar%W(ip) = Lstar%W(ip)
   Rstar%W(iu(ix)) = Rstar%W(iu(ix))
   Lstar%W(irho) = TwoRarefactionDensity(L,Lstar)
   Rstar%W(irho) = TwoRarefactionDensity(R,Rstar)
  else
   Lstar%W(ip) = TwoShockPressure(L,R,p0)
   Lstar%W(iu(ix)) = TwoShockVelocity(L,R,Lstar,p0)
   Rstar%W(ip) = Lstar%W(ip)
   Rstar%W(iu(ix)) = Rstar%W(iu(ix))
   Lstar%W(irho) = TwoShockDensity(L,Lstar)
   Rstar%W(irho) = TwoShockDensity(R,Rstar)
  end if
end select
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

 EulerRiemannSoln = sample(L,R,Lstar,Rstar,0d0)

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
! Calculate shock and rarefaction speeds
SL = L%W(iu(ix)) - L%c*sqrt((gamma + 1)*Lstar%W(ip)/(2*gamma*L%W(ip)) + (gamma - 1)/(2*gamma))      ! Left shock speed
SHL = L%W(iu(ix)) - L%c     ! Head speed of left rarefaction
STL = Lstar%W(iu(ix)) - L%c*(Lstar%W(ip)/L%W(ip))**((gamma - 1)/(2*gamma))  ! Tail speed of left rarefaction
SR = R%W(iu(ix)) + R%c*sqrt((gamma + 1)*Rstar%W(ip)/(2*gamma*R%W(ip)) + (gamma - 1)/(2*gamma))      ! right
SHR = R%W(iu(ix)) + R%c     ! right
STR = Rstar%W(iu(ix)) + R%c*(Lstar%W(ip)/L%W(ip))**((gamma - 1)/(2*gamma))  ! right

if (location <= Lstar%W(iu(ix))) then !left side of contact
 lright = .false.
 if (Lstar%W(ip) > L%W(ip)) then !shock on left side
  if (location <= SL) then !left side of shock
   sampleWithoutVacuum = L%W
  else
   sampleWithoutVacuum = Lstar%W
  end if
 else !rarefaction on left side
  if (location <= SHL) then !left of the rarefaction fan
   sampleWithoutVacuum = L%W
  else if (location >= STL) then !right of the rarefaction fan
   sampleWithoutVacuum = Lstar%W
  else !in rarefaction fan
   call rarefaction(L,R,Wfan,lright,location)
   sampleWithoutVacuum = Wfan
  end if
 end if
else !right side of contact
 lright = .true.
 if (Rstar%W(ip) > R%W(ip)) then !shock on right side
  if (location >= SR) then !right side of shock
   sampleWithoutVacuum = R%W
  else
   sampleWithoutVacuum = Rstar%W
  end if
 else !rarefaction on right side
  if (location >= SHR) then !right of the rarefaction fan
   sampleWithoutVacuum = R%W
  else if (location <= STR) then !left of the rarefaction fan
   sampleWithoutVacuum = Rstar%W
  else !in rarefaction fan
   call rarefaction(L,R,Wfan,lright,location)
   sampleWithoutVacuum = Wfan
  end if
 end if
end if
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

 deallocate(Lstar%W,Rstar%W)
end if
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
EulerFlux(irho) = W(irho)*W(iu(ix))
EulerFlux(imom(ix)) = W(irho)*W(iu(ix))**2 + W(ip)
EulerFlux(iE) = W(iu(ix))*((W(ip)*(1. - cov*W(ip)))/(gamma - 1.) + 0.5*W(irho)*W(iu(ix))**2 + W(ip))
if (nDim > 1) EulerFlux(imom(iy)) = W(irho)*W(iu(ix))*W(iu(iy))
if (nDim > 2) EulerFlux(imom(iz)) = W(irho)*W(iu(ix))*W(iu(iz))
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

   end select
  end do

The grid is then updated using Godunov's method to find the new state of the system:

do i = 1, cells
 ProblemRegion%U(i,:) = ProblemRegion%U(i,:) + (dt/dl)*(ProblemRegion%flux(i-1,:) - ProblemRegion%flux(i,:))
 if (any(ProblemRegion%U(i,:) .ne. ProblemRegion%U(i,:))) then
  print *,' NaN at time ', time, 'in cell ',i, ' in UpdateCells.'
  stop
 end if
end do

The time is updated, and if the program is at one of the times requested for output, information is written to file (as above). In 2D, the problem is first reduced to a 1D problem along each strip of cells in the x direction:

   do j = 1, ycells
    ProblemRegion%U(:,:) = ProblemRegion2D%U(:,j,:)

The solution then proceeds as above, and the current strip is updated with the 1D solution:

    ! Update 2D grid from 1D solution
    ProblemRegion2D%U(:,j,:) = ProblemRegion%U(:,:)
   end do

Once the solution in the x direction has been obtained, the problem is reduced to 1D along each strip in the y direction:

do i = 1, xcells
    ProblemRegion%U(:,:) = ProblemRegion2D%U(i,:,:)
    ! x is considered the direction normal to the cells in 1D routines - so switch v to u for y sweep
    ProblemRegion%U(:,imom(ix)) = ProblemRegion2D%U(i,:,imom(iy))
    ProblemRegion%U(:,imom(iy)) = ProblemRegion2D%U(i,:,imom(ix))

The solution proceeds once again as the 1D problem, and the grid is updated with the new solution in the y direction.

    ! Update 2D grid
    ProblemRegion2D%U(i,:,:) = ProblemRegion%U(:,:)
    ! Switch speeds back
    ProblemRegion2D%U(i,:,imom(iy)) = ProblemRegion%U(:,imom(ix))
    ProblemRegion2D%U(i,:,imom(ix)) = ProblemRegion%U(:,imom(iy))
   end do

Once complete, the time is advanced and output may once again occur. In 3D, the procedure is the same as in 2D, incorporating the z direction as appropriate:

   ! For each strip in x direction, copy to 1D region, then solve
   do j = 1, ycells
    do k = 1, zcells
     ProblemRegion%U(:,:) = ProblemRegion3D%U(:,j,k,:)
     
!!!!!!!!!!!!!!!!!!!!!!!
!      As before      !
!!!!!!!!!!!!!!!!!!!!!!!
     
     ! Update 3D grid
     ProblemRegion3D%U(:,j,k,:) = ProblemRegion%U(:,:)
    end do
   end do

   ! For each strip in y direction, copy to 1D region, then solve
   do i = 1, xcells
    do k = 1, zcells
     ProblemRegion%U(:,:) = ProblemRegion3D%U(i,:,k,:)
     ! x is considered the direction normal to the cells in 1D routines - so switch v to u for y sweep (w doesn't change - still passively advected)
     ProblemRegion%U(:,imom(ix)) = ProblemRegion3D%U(i,:,k,imom(iy))
     ProblemRegion%U(:,imom(iy)) = ProblemRegion3D%U(i,:,k,imom(ix))
     
!!!!!!!!!!!!!!!!!!!!!!!
!      As before      !
!!!!!!!!!!!!!!!!!!!!!!!

     ! Update 3D grid
     ProblemRegion3D%U(i,:,k,:) = ProblemRegion%U(:,:)
     ! switch speeds back
     ProblemRegion3D%U(i,:,k,imom(ix)) = ProblemRegion%U(:,imom(iy))
     ProblemRegion3D%U(i,:,k,imom(iy)) = ProblemRegion%U(:,imom(ix))
    end do
   end do

   ! For each strip in z direction, copy to 1D region, then solve
   do i = 1, xcells
    do j = 1, ycells
     
     ProblemRegion%U(:,:) = ProblemRegion3D%U(i,j,:,:)
     ! x is considered the direction normal to the cells for 1D routines - so switch w to u for y sweep (v doesn't change - still passively advected)
     ProblemRegion%U(:,imom(ix)) = ProblemRegion3D%U(i,j,:,imom(iz))
     ProblemRegion%U(:,imom(iz)) = ProblemRegion3D%U(i,j,:,imom(ix))

!!!!!!!!!!!!!!!!!!!!!!!
!      As before      !
!!!!!!!!!!!!!!!!!!!!!!!
     
     ! Update 3D grid
     ProblemRegion3D%U(i,j,:,:) = ProblemRegion%U(:,:)
     ! switch speeds back
     ProblemRegion3D%U(i,j,:,imom(ix)) = ProblemRegion%U(:,imom(iz))
     ProblemRegion3D%U(i,j,:,imom(iz)) = ProblemRegion%U(:,imom(ix))
     
    end do
   end do

and time is updated and output is done if requested.

Godunov Solver​

This program solves either the inviscid Burgers equation or the linear advection equation with the given boundary or initial conditions. First, input is taken from a file:

NAMELIST/ProblemData/ starttime,finaltime,frames,xL,xR,cells,alpha,beta,BC,Eq,IC,advectionConst,CFL

! Make sure problem.data exists. If it does, read in. If not, output error and quit.
inquire(file='GodunovProblem.data', exist=ex)
if (ex) then
 OPEN(UNIT=30, FILE='GodunovProblem.data', STATUS="OLD")
 READ(30,NML=ProblemData)       ! read into variables in namelist
 CLOSE(30)
else
 print *, 'Missing GodunovProblem.data. Exiting.'
 stop   ! exit if missing
end if

The domain is then initialized according to the selected conditions:

select case (IC)
 case (Gaussian)
  do i = 1, cells
   ProblemRegion%u(i) = alpha*exp(-1.*beta*((i-0.5)*dx+xL)**2)
  end do
 case (SquareWave)
  do i = 1, cells
  x=(i-0.5)*dx + xL
   if (x <= 0.3) then
    ProblemRegion%u(i) = 0
   else if (x >= 0.7) then
    ProblemRegion%u(i) = 0
   else
    ProblemRegion%u(i) = 1
   end if
  end do
 case (BurgersShock)
  do i = 1, cells
  x=(i-1)*dx + xL
   if (x <= 0.5) then
    ProblemRegion%u(i) = -0.5
   else if (x >= 1.) then
    ProblemRegion%u(i) = 0
   else
    ProblemRegion%u(i) = 1
   end if
  end do
end select

The boundary conditions are then applied:

select case (BC)
 case(periodic) ! Wraps around
  ProblemRegion%u(0) = ProblemRegion%u(cells)
  ProblemRegion%u(cells+1) = ProblemRegion%u(1)
 case(transparent)      ! Wave flows out
  ProblemRegion%u(0) = ProblemRegion%u(1)
  ProblemRegion%u(cells+1) = ProblemRegion%u(cells)
end select

After this, the time step is determined using the maximum speed in the domain and the given CFL coefficient.

! Determine maximum speed in domain
maxSpeed = 0.
do i = 0, cells+1
 if (abs(ProblemRegion%u(i)) > maxSpeed) maxSpeed = abs(ProblemRegion%u(i))     ! (See Toro, Ch. 5.6)
end do

! Calculate delta t
dt = CFL*dx/maxSpeed

! If time step would pass the next output frame, recalculate
if (time + dt > nexttime) dt = nexttime - time

The intercell flux is then calculated and the velocity array is updated to the next time.

do i = 0, cells
 uL = ProblemRegion%u(i)
 uR = ProblemRegion%u(i+1)   
 select case (Eq)
  case (advection)
   ProblemRegion%flux(i) = LinearAdvectionFlux(RiemannSoln(uL,uR))

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
function RiemannSoln(uL,uR)

if (uL > uR) then
 S = 0.5*(uL + uR)
 if (S <= 0) then
  RiemannSoln = uR
 else
  RiemannSoln = uL
 end if
else
 if (uL >= 0) then
  RiemannSoln = uL
 else if (uR <= 0) then
  RiemannSoln = uR
 else
  RiemannSoln = 0
 end if
end if
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

  case (inviscidBurgers)
   ProblemRegion%flux(i) = BurgersFlux(RiemannSoln(uL,uR))
 end select
end do

! Loop through cells using previously calculated flux
do i = 1, cells
 ProblemRegion%u(i) = ProblemRegion%u(i) + (dt/dx)*(ProblemRegion%flux(i-1) - ProblemRegion%flux(i))
end do

Finally, at each requested output frame, the current system state is written to a file:

write (40,*) 'x u'       ! write column headers
do i = 1, cells
 x = xL + (i-0.5)*dx    ! Data is located at center of each cell
 write (40,*) x, ProblemRegion%u(i)
end do

Exact Riemann Solver​

The Riemann problem is the solution of a system of hyperbolic conservation laws (in our case the Euler equations). This solver assumes two one-dimensional gases (ideal or with a covolume parameter b) separated at the origin x=0, with initial densities, velocities, and pressures given as input (vacuum is allowed, both as input and as a result of system evolution). It gives a solution to the set of equations

\rho + \rho u = 0$ $\rho u + \rho u^2 + p = 0$ $E + u(E + p) = 0$

with initial conditions

\rho(x,0) = \begin{cases} \rho_L, & x < 0 \\ \rho_R, & x > 0 \end{cases}$ $\rho(x,0)u(x,0) = \begin{cases} \rho_L u_L, & x < 0 \\ \rho_R u_R, & x > 0 \end{cases}$ $E(x,0) = \begin{cases} E_L, & x < 0 \\ E_R, & x > 0 \end{cases}$

The solutions to this IVP at later times are split into four regions - to the left, the initial left state, unperturbed by interaction with the right; between the left and center-left regions, there will either be a rarefaction fan or a shock (assumed ideal), across which the new equilibrium left state holds; the left and right fluids will be separated by a contact discontinuity, across which the new equilibrium right state holds (note that the right and left equilibrium states differ only in their densities); the center-right region will then be separated by a rarefaction fan or (again ideal) shock from the rightmost region, which is still in the initial right state. If vacuum is present, one or more rarefaction fans connect the unperturbed states with the region of vacuum.

The program reads input from a file called problem.data, and stores many of the relevant parameters in global variables.

NAMELIST/ProblemData/ rhoL,uL,pL,rhoR,uR,pR,gamma,tol,starttime,finaltime,frames,xL,xR,res

inquire(file='problem.data', exist=ex)
if (ex) then
 OPEN(UNIT=30, FILE='problem.data', STATUS="OLD")
 READ(30,NML=ProblemData)                               ! read into variables in namelist
 CLOSE(30)
else
 print *, 'Missing problem.data. Exiting.'
 stop           ! exit if missing
end if

The state vectors W_L and W_R are created with overloaded constructors that also calculate useful parameters A and B and the speed of sound c for use later in the program. (Note code between comment lines is excerpted from its actual location.)

L = region(rhoL,uL,pL)
R = region(rhoR,uR,pR)

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
function ConstructRegion(rho,u,p)
        
! Fill state
ConstructRegion%W(irho) = rho
ConstructRegion%W(iu) = u
ConstructRegion%W(ip) = p
! Calculate useful constants
if (.not. (rho == 0 .and. p == 0)) then              ! state is not vacuum
 ConstructRegion%A = 2.*(1. - cov*ConstructRegion%W(irho))/((gamma + 1.)*ConstructRegion%W(irho))
 ConstructRegion%B = (gamma - 1.)*ConstructRegion%W(ip)/(gamma + 1.)
 ConstructRegion%c = sqrt(gamma*ConstructRegion%W(ip)/(ConstructRegion%W(irho)*(1. - cov*ConstructRegion%W(irho))))
else
 ConstructRegion%c = 0          ! Set speed of sound to zero if state is vacuum (to avoid NaN)          !!! could also set A,B to zero, but since they don't represent physical quantities, not much point
end if
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The initial states are then checked for vacuum, and the state of the system (initial vacuum, vacuum created, or no vacuum) is set.

call check_vacuum()

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
check_vacuum()

critdiff = 2.*L%c/(gamma - 1.) + 2.*R%c/(gamma - 1.)      ! Difference between right and left velocities must be strictly less than this for pressure to remain positive
        
if ((L%W(ip) == 0 .and. L%W(irho) == 0) .or. (R%W(ip) == 0 .and. R%w(irho) == 0)) then  ! one of the initial states is vacuum
 vacuum = 1
 print *, 'One of initial states is vacuum.'
else if (critdiff <= R%W(iu) - L%W(iu)) then             ! vacuum is created by initial states
 vacuum = 2
 print *, 'Vacuum is created by initial states.'
else            ! no vacuum
 vacuum = 3
end if
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

If there is no vacuum, the state vectors L and R are then passed to the Newton-Raphson solver, along with a guess at a good value for p₀ (currently not very effective - just uses the average value of the initial states). The solver takes a guess at a root for a function, the function itself, and its derivative, and returns the root of the function (to within specified tolerance). If a negative value is encountered at any point during iteration, the guess is reseeded with the tolerance of the solver (since the two values solved for are density and pressure, negative solutions would be unphysical - would be pretty easy to rework to avoid NaNs if negative values were allowed).

! Use Newton-Raphson solver to find pressure in star region
Lstar%W(ip) = NRSolver(pguess(L,R),covolumeF,covolumeDf)

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨ 
! Newton-Raphson iterative solver. Takes guess of root of equation, equation, and derivative and calculates root to tol.
function NRSolver(guess,f,df)

! do while loop (exits once change is less than tolerance)
do
 guess_old = guess
 guess = guess_old - f(guess_old)/df(guess_old)           ! Newton-Raphson formula
 delta = 2.*abs((guess - guess_old)/(guess + guess_old))  ! Relative change in pressure from previous iteration
 if (guess < 0) guess = tol                               ! If we get a negative value (unphysical for our applications), set guess to tol      !!! could use nan inequality for more general case
 if (delta < tol) then                                    ! Exit upon sufficient accuracy
  NRSolver = guess
  exit
 end if
end do
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The functions f and df (f') for either the left or right regions are given by

f_{reg}(p_*) = \begin{cases} (p_* - p_{reg}) [\frac{A_{reg}}{p_* + B_{reg})}], & p_* > p_{reg} \\ f_{reg}(p_*) = \frac{2 c_{reg}}{\gamma - 1}[(\frac{p_*}{p_{reg}})^{\frac{\gamma - 1}{2 \gamma}} - 1], & p_* \leq p_{reg} \end{cases}$ $f'_{reg}(p_*) = \begin{cases} (\frac{A_{reg}}{B_{reg} + p_*})^{\frac{1}{2}}[1 - \frac{p_* - p_{reg}}{2(B_{reg} + p_*)}], & p_* > p_{reg} \\ \frac{1}{\rho_{reg} c_{reg}}(\frac{p_*}{p_{reg}})^{-\frac{\gamma + 1}{2 \gamma}}, & p_* \leq p_{reg} \end{cases}$

See below for derivation of function f (in the ideal gas case. The derivation is identical for the covolume case, except as noted below).

Once the pressure is solved for, velocity and densities for the star regions can then be calculated:

! Calculate velocity and densities in star regions
Lstar%W(iu) = exactVelocity(Lstar%W(ip))

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
! Calculates velocity in star region
exactVelocity(L,R,pstar)
        
exactVelocity = 0.5*(L%W(iu) + R%W(iu)) + 0.5*(f(pstar,R) - f(pstar,L)) ! Velocity in star region is the same on either side of CD
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Rstar%W = Lstar%W
Lstar%W(irho) = exactDensity(Lstar,L)
Rstar%W(irho) = exactDensity(Rstar,R)

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
! Calculates densities in star regions
exactDensity(starReg,reg)
        
! Density is different for shocks and rarefactions - may differ on either side of CD
if (starReg%W(ip) > reg%W(ip)) then
 exactDensity = reg%W(irho)*((gamma - 1.)/(gamma + 1.) + (starReg%W(ip)/reg%W(ip)))/(((gamma - 1.)/(gamma + 1.))*(starReg%W(ip)/reg%W(ip)) + 1.)
else
 exactDensity = reg%W(irho)*(starReg%W(ip)/reg%W(ip))**(1./gamma)
end if
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The final value of pstar is then printed to console and output is written to files. If vacuum is present, no calculation is required; the output can just be written to files. One file is written for each time state requested. During output for the case with no vacuum, the speeds of shocks and rarefaction boundaries are calculated:

SL = L%W(iu) - L%c*sqrt((gamma + 1)*Lstar%W(ip)/(2*gamma*L%W(ip)) + (gamma - 1)/(2*gamma))      ! Left shock speed
SHL = L%W(iu) - L%c                                                                             ! Head speed of left rarefaction
STL = Lstar%W(iu) - L%c*(Lstar%W(ip)/L%W(ip))**((gamma - 1)/(2*gamma))                          ! Tail speed of left rarefaction
SR = R%W(iu) + R%c*sqrt((gamma + 1)*Rstar%W(ip)/(2*gamma*R%W(ip)) + (gamma - 1)/(2*gamma))      ! right
SHR = R%W(iu) + R%c                                                                             ! right
STR = Rstar%W(iu) + R%c*(Lstar%W(ip)/L%W(ip))**((gamma - 1)/(2*gamma))                          ! right

Then a large conditional is entered to determine which region we are sampling at the given time and position. If we are located in a rarefaction fan, the density, velocity, and pressure are calculated, which requires the solution of an equation for the density using the Newton-Raphson solver if the fluid has a covolume parameter. A similar conditional is used for the vacuum states, but no star region is required.

! End once t > finaltime (last output will be t = finaltime)
do while (t <= finaltime)
 !print fileformat, "out/result",frame,".dat"
 write(filename,fileformat) "out/result",frame,".dat"   ! filename to output - goes to out/resultxxxxx.dat
 open(unit=40,file=filename,iostat=stat,status='replace')
 write (40,*) 'x rho u p'                               ! write column headers
 x = xL
 ! Calculate values for each position - similar to time above
 !!! could parallelize if I changed from while to for loop
 do while (x <= xR)
  if (t == 0) then      !initial state (so we don't divide by zero)
   if (x < 0) write (40,*) x, L%W
   if (x > 0) write (40,*) x, R%W
   if (x == 0) write (40,*) x, 0, 0, 0
  else
   if (x/t <= Lstar%W(iu)) then !left side of contact
   lright = .false.
    if (Lstar%W(ip) > L%W(ip)) then !shock on left side
     if (x/t <= SL) then !left side of shock
      write (40,*) x, L%W
     else
      write (40,*) x, Lstar%W
     end if
    else !rarefaction on left side
     if (x/t <= SHL) then !left of the rarefaction fan
      write (40,*) x, L%W
     else if (x/t >= STL) then !right of the rarefaction fan
      write (40,*) x, Lstar%W
     else !in rarefaction fan
      call rarefaction(L,x,t,Wfan,lright)

!∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨∨
subroutine rarefaction(reg,x,t,Wfan,lright)

! If fluid has covolume constant, iterate to solve for density
if (cov > 0) then
 if (.not. lright) then
  beta = ((((gamma - 1.)*(L%W(iu) + 2.*L%c*(1. - cov*L%W(irho))/(gamma - 1.) - x/t))**2.)/(gamma*L%W(ip)*((1. - cov*L%W(irho))/L%W(irho))**gamma)) ! Constant related to Riemann invariant
  Wfan(irho) = NRSolver(rhoguess(L,R),densityF,densityDf)
  Wfan(ip) = reg%W(ip)*(((1. - cov*reg%W(irho))/reg%W(irho))*(Wfan(irho)/(1. - cov*Wfan(irho)))**gamma)
  c = sqrt((Wfan(ip)*gamma)/(Wfan(irho)*(1. - cov*Wfan(irho))))                                         ! Speed of sound at current location
  Wfan(iu) = x/t + c
 else
  beta = (((x/t - R%W(iu) + (2.*R%c*(1. - cov*R%W(irho))/(gamma - 1.)))**2.)/(gamma*R%W(ip)*((1. - cov*R%W(irho))/R%W(irho))**gamma))   ! Constant related to Riemann invariant
  Wfan(irho) = NRSolver(rhoguess(L,R),densityF,densityDf)
  Wfan(ip) = reg%W(ip)*(((1. - cov*reg%W(irho))/reg%W(irho))*(Wfan(irho)/(1. - cov*Wfan(irho)))**gamma)
  c = sqrt((Wfan(ip)*gamma)/(Wfan(irho)*(1. - cov*Wfan(irho))))                                         ! Speed of sound at current location
  Wfan(iu) = x/t - c
 end if
! Otherwise solution is closed-form
else
 ! Rarefaction fan differs depending on whether it's a right or left rarefaction
 if (.not. lright) then
  Wfan(irho) = reg%W(irho)*(2/(gamma + 1) + (gamma - 1)*(reg%W(iu) - x/t)/((gamma + 1)*reg%c))**(2/(gamma + 1))
  Wfan(iu) = (2/(gamma + 1))*(reg%c + (gamma - 1)*reg%W(iu)/2 + x/t)
  Wfan(ip) = reg%W(ip)*(2/(gamma + 1) + (gamma - 1)*(reg%W(iu) - x/t)/((gamma + 1)*reg%c))**(2*gamma/(gamma + 1))
 else
  Wfan(irho) = reg%W(irho)*(2/(gamma + 1) - (gamma - 1)*(reg%W(iu) - x/t)/((gamma + 1)*reg%c))**(2/(gamma + 1))
  Wfan(iu) = (2/(gamma + 1))*(-reg%c + (gamma - 1)*reg%W(iu)/2 + x/t)
  Wfan(ip) = reg%W(ip)*(2/(gamma + 1) - (gamma - 1)*(reg%W(iu) - x/t)/((gamma + 1)*reg%c))**(2*gamma/(gamma + 1))
 end if
end if
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

      write (40,*) x, Wfan
     end if
    end if
   ! Repeat entire sequence above for right side (lright = .true.)
   end if
  end if
  x = x + dx            ! increment x to next position
 end do
 t = t + dt             ! increment t to next time
 frame = frame + 1      ! update frame counter
end do

These output files can then be processed by VisIt, Mathematica, etc. to make nice pictures.

View ​online or download 1​, 2​, 3​, 4​, 5​ the tests from Toro. The vacuum routines have only been sanity tested, while the covolume routines have been tested once against the case in Table III(b) in Toro, E.F. "A Fast Riemann Solver with Constant Covolume Applied to the Random Choice Method," Int. J. Num. Meth. in Fluids, Vol. 9 (1989).

Deriving the function f and the Newton-Raphson procedure

f

From the shock jump conditions

\rho_L \hat{u_L} = \rho_{*L} \hat{u_*}$ $\rho_L \hat{u_L}^2 + p_L = \rho_{*L} \hat{u_*}^2 + p_*$ $\hat{u_L}(\hat{E_L} + p_L) = \hat{u_*}(\hat{E_{*L}} + p_*)$

and the energy and EOS of an ideal gas

e = \frac{p}{(\gamma - 1) \rho}$ $p = \rho R T$

where \hat{u} = u_L - S_L is the velocity of the fluid in the frame of the shock, the function f_L(p_*) = u_L - u_*,

f_L(p_*) = (p_* - p_L) [\frac{2}{(\gamma + 1) \rho_L (p_* + \frac{\gamma - 1}{\gamma + 1} p_L)}]$

can be derived to connect the left state with the pressure in the central, interaction region in the case of a shock. Similarly, in the case of a rarefaction, instead of the shock jump conditions, the isentropic equation

\frac{p_*}{\rho_*} = \frac{p_L}{\rho_L}$

and the Generalized Riemann Invariant

u_L + \frac{2 c_L}{\gamma - 1} = u_* + \frac{2 c_{*L}}{\gamma - 1}$

where c = \sqrt{\frac{\gamma p}{\rho}} is the speed of sound in the region, can be used to construct

f_L(p_*) = \frac{2 c_L}{\gamma - 1}[(\frac{p_*}{p_L})^{\frac{\gamma - 1}{2 \gamma}} - 1]$

The same analysis can be performed on the right side of the system to yield the same results.

Covolume

In the covolume case, the internal energy becomes

e = \frac{p(1 - b \rho)}{\rho(\gamma - 1)}$

and the sound speed is

c = (\frac{p \gamma}{\rho (1 - b \rho)})^\frac{1}{2}$

which reduce to the ideal gas case when b = 0.

Newton-Raphson

The Riemann problem can then be solved by finding the value p = p_* such that

f_L(p_*) + f_R(p_*) + u_R - u_L = 0$

The function and its first and second derivatives are well-behaved and monotonic, so we can use a simple Newton-Raphson iterative method, derived from a Taylor series expansion of f(p_*):

If our initial guess is p_k and f(p_k) \ne 0, we want to change p_k so that f(p_k + \delta) = f(p_{k+1}) = 0; the Taylor expansion is given by

f(p_k + \delta) = f(p_k) + \delta f'(p_k) + O(\delta^2)$

Dropping second order and higher terms, since we have f(p_k + \delta) = 0, we know f(p_k) + \delta f'(p_k) = 0, and solving for \delta,

\delta = -\frac{f(p_k)}{f'(p_k)}$

and

p_{k + 1} = p_k + \delta = p_k - \frac{f(p_k)}{f'(p_k)}$

The iteration procedure is terminated at the term such that

|\frac{2 (p_{k + 1} - p_k)}{p_{k + 1} + p_k}| < Tol$

where Tol is the (generally small) tolerance passed to the solver.

Source: See Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics