wiki:u/adebrech/code

Version 9 (modified by adebrech, 8 years ago) ( diff )

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

with initial conditions

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 WL and WR 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 p0 (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 () for either the left or right regions are given by

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

and the energy and EOS of an ideal gas

where is the velocity of the fluid in the frame of the shock, the function ,

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

and the Generalized Riemann Invariant

where is the speed of sound in the region, can be used to construct

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

and the sound speed is

which reduce to the ideal gas case when .

Newton-Raphson

The Riemann problem can then be solved by finding the value such that

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 and , we want to change so that ; the Taylor expansion is given by

Dropping second order and higher terms, since we have , we know , and solving for ,

and

The iteration procedure is terminated at the term such that

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

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

Attachments (10)

Note: See TracWiki for help on using the wiki.