# The Godunov Method for Euler Equations, using EXACT Riemann Solver

## Big Picture

- Define a piece-wise function of initial data

- This produces local Riemann Problems (LRPs), centered on intercell boundaries.

- Solve these local RP's using the Exact Riemann Solver, and use solution to compute fluxes

- Use these fluxes to update cells.

Here is a picture of the situation:

In this picture, WL and WR denote the arrays of primitive fluid variables (namely, rho, u, P) designated as the values of those variables at the adjacent cell centers (here at cell=i-1 and i), the asterisk indicates a star region that is bound by two outgoing arrows, the arrows indicate the wave pattern generated in the LRP, F is the numerical flux (described below), and the conservative formula is given by the Godunov scheme, also presented below.

If we were to blow the left star region up (on the intercell boundary between cell i-1 and i), what we would have is 3 non-linear waves (contact, shock(s) and/or rarefaction(s)) being generated at this boundary IF the left and right states (WL and WR) were different. If these states are the same, then NO non-linear waves would be generated. We would like to first find pstar, the value of pressure between the left and right waves, in the region so-called the "star region". Once this is known, we can determine the types of waves generated. Since they are propagating away from the intercell boundary, the exact values of the fluid variables along the intercell boundary depend on their relative speeds. We then would sample the solution to the LRP, along this boundary only. Once we have found the solution at this boundary, given by

we compute the flux at this boundary. These fluxes are the fluxes of the Euler equations, in conservation form. That is,

where gamma is the ratio of specific heats.

These fluxes are then fed into the Godunov scheme to update each (conserved) variable for each cell. In conservative form, the Godunov method is given by:

where

is the numerical flux (equal to the physical flux in the case of the Godunov method for the Euler equations) on the left boundary, and

is the flux along the right intercell boundary.

## The method

### Order of accuracy

This method is first order accurate in both spatial and temporal derivatives. Since the error, E, is proportional to the time step, h, to the nth order of accuracy:

For the Godunov scheme, we have that

### Finding Pstar

Pstar is the solution for pressure in the "star region", which is set up at the intercell boundaries by the left and right cells involved in the local Riemann problems. This region is where the contact discontinuity lies, flanked on either side by a given combination of 2 shock or rarefaction waves. The first step at each intercell boundary is to solve for pstar. Pstar is found by Newton's iterative scheme, with the routine as it was in the Exact Riemann solver described before.

### Sampling the solution

Once Pstar is found to within a desired degree of accuracy, ustar is given immediately by a simple formula in Toro, chapter 4. These quantities are used to determine the combination of waves present along the intercell boundary. Once these are determined, we look for the speeds of the wave, relative to the grid speed, dx/dt = 0, which corresponds to the local intercell boudary. Depending on where the waves are, relative to this interface, the algorithm prescribes a density, velocity, and pressure for the boundary. These quantities are then used to update the cells, as described above.

### Constraining the time step

The time-step is constrained by the CFL condition, which uses the max wave speed on the grid. The constraint used in the God. code here allows any given wave on the grid to travel a full delta_x in a time-step. In allowing this, we achieve a more efficient time-marching scheme, but with the assumption that wave interactions do not lead to wave accelerations. The CFL condition, rearranged for dt, is given by:

dt = CFL*dx/smax

Where CFL is a coefficient that ranges from 0<CFL<1. A CFL close to 1, produces the largest time step, but it should be used cautiously during times when errors can be produced, for which a smaller time step should be used. In the God. code described here, the CFL was reduced by 80% for the first 5 time steps of the code.

Dx is given by the domain length divided by the number of cells.

Smax is found by looping over the grid and checking for the max speed. Smax is computed in the subroutine for calculating dt.

### Boundary Conditions

Transmissive boundaries were set up for the current code. These boundary conditions (BC's) theoretically do not produce waves at the boundary, and so allow waves to travel right through them. Such BC's are prescribed by setting the values of the fluid variables at the ghost zone cell equal to the adjacent physical cell. Thus, the local Riemann problem has WL = WR at the boundary of the physical and ghost domains, and so the solution is trivial (no wave generated).

### Variable transformations

The Exact Riemann Solver (contained in the separate module below — that includes routines to solve for pstar, ustar, and sampling along intercell boundary) uses PRIMITIVE variables, that is: density (rho), velocity (u), pressure (P). But, the Godunov scheme discretizes the Euler equations in CONSERVATIVE form, and so updating cells using this scheme, requires transforming to the conservative variables (density, momentum (p), energy (E)) updating them using Godunov formula, and then converting back to primitive form for use with the Exact Riemann Solver in the next time step. The transformation from primitive to conservative is as follows:

## Expectations

For initial data that consists of left and right states separated by a discontinuity at position x0, one would expect the generation of non-linear waves to originate at the discontinuity, and spread outward from the discontinuity over time. Outside of the region in causal contact with the generated waves, the initial condition is not affected. Indeed, this is what happens, as can be seen below in the results section.

## The code

Here is an html version of the main program, with encompassing module, and problem.data file.

### Program Structure

Here is an outline for the order in which the different routines are called in the program, and what they do. The Arabic numerals give the order of operations in the main program and the lower case Roman numerals give the order of operations in the time iteration part of the program. The schematic of the do-loop presented in this part of the outline, give the order of operations for the subroutine "Get Fluxes", which loops through the cells on the grid to set up and solve the LRP's. The 2nd set of Roman numerals should actually read: iv, v, vi, but trac is getting confused.

- Call init — Sets up initial conditions. This subroutine (SR) loops through i=1,mx and checks whether the cell center ((i-0.5)*dx) is to left or right of discontinuity (x0), and then assigns appropriate value to that cell. (All initial conditions for the tests in Toro consist of a left and right state separated by a discontinuity.
- Call output — begins the output file — records the following info - position, density, velocity, pressure, and time and timestep. This output file is named output.out.
- Sets time = 0.0 — starts the clock for iteration scheme
- Enters the update cycle, where from n=1 to n=maxitr (n is the counter for the time iteration scheme) it performs the following commands -

- Call Bcond — Boundary condition SR that sets left and right ghost zones to adjacent physical cells (recall boundary conditions are set to transmissive).
- Call CFLcond — Computes timestep using cfl coeff and max speed, where the max speed is also computed in this SR.
- Call GetFluxes — This SR sets up the local Riemann Problem (LRP) and solves it exactly to get rho(x/t), u(x/t), P(x/t) at x/t = 0 (see above explanation of sampling routine). This is the most complicated structure of the program. A schematic of the loop is as follows:

- Call Update — Updates cells using flux computed in previous step
- Call Output — Write quantities to file
- Begin loop again for n+1 time step.

### Standard Out and Program Outfile

You can run the program with the command ** ./program > program.out & **, and the standard out will be printed to file. The following image is of the information the standard out produces:

Basically, it gives the time, the values of q, the solution to the LRP, and the flux at every cell during that time step.

Additionally, the program's output is printed to output.out, and it gives the position and fluid variables for each time, in neat columns that are easy to plot. Such an output.out file for Test 4 is here.

## Debugging

I ran into issues with the original revisions of this code. I started debugging by first adding sequential print statements throughout the code, module, and subroutines, to make sure the routines were called in the correct order, and that values of cells, were being delivered correctly to the various parts of the program. These all checked out. Next I began adding print statements to the iterative scheme, that checked the time step (dt), and the cfl. For the first time step, dt was calculated correctly, and the cfl was also correct. The next thing I checked was that between time = 0 and time = 0 + dt, the values of q were being computed correctly at those cells near the discontinuity where I expected to see non-linear wave formation. I made a spreadsheet of values, and calculated the updated values of q for these cells. Here is a screenshot of that spreadsheet:

Using my standard out, I checked the updated values of q I calculated for the cells adjacent to the discontinuity, and it checked out. Running the code again, after the 8th or 9th time step, the code stopped after the Newton Iteration Scheme failed to converge… But everything was being computed correctly! A look at the time by this 8th or 9th time step hinted to a bug… The time was extremely small, at this rate, the simulation wouldn't get done for hundreds and hundreds of time steps. This was wrong. Checking then the cfl, dt, and max speed for the time steps between 2 and 8, provided the answer. Smax was not changing, but the cfl was getting smaller and smaller, decreasing the time step by many factors. I checked the code I was using to set up the cfl coefficient, and found an error in the way I was doing so. I was rewriting the cfl coeff each time, and then re-using this value when I'd ask for 20% of the cfl coeff during those first 5 time steps (see above cfl condition subroutine explanation). I adjusted the code to not do this any more, and sure enough, the code ran smoothly afterward… No more divergence in iteration scheme for pstar, and the code finished with the correct output. Interesting an error that led to divergence in the iteration scheme was produced with such small time steps…

## Results

To check the hypothesis that the initial condition would change initially at the discontinuity with the generation of non-linear waves, and then move outward with time, I plotted some sequential time steps in density. Here is the initial condition for density:

As time proceeds, you can see that a disturbance in density is generated at the border between the discontinuous left and right states and travels outward with time:

In the last time step here, we begin to see the shape of the density curve shown in the following plot to emerge.

The above plot is that result that Toro published in chapter 6 for Test 4, but without the exact solution given by the Exact Riemann Solver. I simply wanted to check that the Godonuv scheme I wrote was performing correctly, so decided to just post these results for test 4, analyze them, and move on with my studies of numerics. In the next section I will analyse Toro's results to determine the difference in solution between the Exact Riemann Solver and Godunov method, to understand why the differences occur, and which is a better model for differing initial conditions.

Here are the results for Test 4 of velocity and pressure:

They match Toro's results exactly.

## Discussion

### Guessing the Solution

It seems like the types of waves generated in the Riemann problem can be 'guessed' at from the initial conditions, but there is not a hard and strict rule to learn, it's more of a type of intuition to be developed. For simple cases, such as those given in chapter 4 of Toro, you can guess that if you have a left and right state that are flowing into each other, you likely will get out 2 shocks. If they are flowing away from each other, you may get 2 rarefactions. If there is no motion in the left and right states, but the left is at a higher density or pressure, you can imagine it forcing a shock into the right state (a right traveling shock). By the same token, you can imagine this would create a left traveling rarefaction. Instead of using the iterative method to find pstar, you could guess at the value of pstar by using exact solutions if the 2 waves consisted of 2 shocks or 2 rare's. Another approximation could be the geometric mean of the left and right pressures. Once you'd have an estimate for pstar, to further your estimate of the solution, you'd have to consider how far those waves traveled, whether they disturbed the left and right states by the time you are interested, and so on. Obviously this would be quite tricky business, and so it seems reasonable to try to hone your skills for estimating what types of waves may be generated at any given discontinuity, and qualitatively how this may effect the data.

### Test 4's Exact Solution

The initial conditions for test 4 produce a left facing (originating from the left initial state) shock wave, traveling very slowly to the right, a right traveling contact, and a right traveling shock wave. This makes intuitive sense from the initial conditions:

Now the EXACT solution at t=0.035 could be computed from the Exact Riemann Solver alone, and this would be done by looping through the grid at different dx/dt, where dt=0.035 and dx is each cell's position in the loop, to determine the relative position of this 'grid speed' to the shocks and contact generated at the boundary. If the grid speed is ahead of all shocks (far enough to the right, given all waves in this solution are traveling to the right), then the initial state is not yet disturbed, and so the EXACT solution there just is the initial condition there. The same holds true if the grid speed were far enough to the left, such that the entire system of generated waves has not interacted with the initial state. Now the tricky thing is to get the solution inside the wave structure, and honestly, it seems like something you'd ever only want to do with the help of a computer program, looping through tests and conditions to determine the new values of the disturbed initial condition. Here is the plot from Toro of the exact solution of Test 4 using the Exact Riemann Solver:

These plots make qualitative sense, based on the aforementioned concepts. Now… how does this change when we solve these initial conditions with the Godunov scheme?

### How the Godunov Solution Compares

As I have talked about a LOT in this entire page, the method for 'evolving' the initial data is different for the Godunov scheme, but you'd imagine that it should produce roughly the same results as the Exact Riemann Solver, as it is doing the same thing, but just over successive time steps. So how does the method compare? Here is the plot from Toro showing the Godunov method over the Exact method:

It is interesting to see the Godunov method performs quite well for the pressure and velocity, but falls short with density (and hence internal energy, which is a function of density). As I have shown above, we see the development of the structure that is approaching the discontinuous jump in density within the star region of the exact solution. By solving the Riemann problem with the Godunov method over successive time steps, we have smoothed out this region in our discretization of the problem. This must happen when we are inside of the star region *and* there is a discontinuity that occurs in a fluid variable. I can imagine this smoothing happening only within the (global) star region , as outside of that region, the waves have not yet penetrated the surrounding medium. Given the only variable with the potential of not being continuous in the global star region is density, I imagine the only smoothing we see of sharp peaks in the solution is with density, and any quantities that depend on density. This then supports why we see close agreement with pressure and velocity between methods.

### Why Use One Method Over the Other

You can imagine a situation with a much more complicated initial condition than just a flat left and right state separated by a discontintuity, but rather a series of discontinuites. For a system such as this, it seems to be straightforward to evolve the initial condition over time. In contrast, an exact solution for such a system seems non-trivial.

## Outstanding Questions

- Do we see an effect of the waves to spread away from the discontinuity? Yes!
- What is the effect of choosing a gamma=1.4? How would a different gamma affect the outcome of these results?

### Attachments (14)

- GodEulRho.png (19.0 KB ) - added by 11 years ago.
- GodEulP.png (21.1 KB ) - added by 11 years ago.
- GodEulU.png (22.3 KB ) - added by 11 years ago.
- God.2.png (50.2 KB ) - added by 11 years ago.
- God.png (50.2 KB ) - added by 11 years ago.
- LoopSR.png (44.1 KB ) - added by 11 years ago.
- program_out.png (28.7 KB ) - added by 11 years ago.
- God_density_init.png (12.8 KB ) - added by 11 years ago.
- God_density.png (21.5 KB ) - added by 11 years ago.
- GodDebug.png (72.1 KB ) - added by 11 years ago.
- test.png (5.8 KB ) - added by 11 years ago.
- test4.png (8.7 KB ) - added by 11 years ago.
- ExactTest4.png (29.1 KB ) - added by 11 years ago.
- GodTest4.png (37.9 KB ) - added by 11 years ago.

Download all attachments as: .zip

**Note:**See TracWiki for help on using the wiki.