CS205b/CME306Lecture 11Supplementary Reading: Osher and Fedkiw, §14.2, §14.3.2-4, §14.4-5; Leveque §4.1, §12.9-101 Discrete Conservation FormTo ensure that shocks and other steep gradients are captured by our scheme (i.e. they move atthe right speed even if they are unresolved) we must write the equation in a d iscrete conservationform. That is, a form in which the rate of change of conserved quantities is equal to a difference offluxes. This form guarantees that we discretely conserve the total amount of the states of φ (e.g.mass, momentum and energy) that are present, analogously with the integral form given byddtZΩφ dV +Z∂Ωf(φ) · dA =ZΩs(φ) dV.More importantly, this can be shown to imply that steep gradients or jumps in the discrete profilespropagate at the physically correct speeds.Usually, conservation form is derived for control volume methods, that is, methods that evolvecell average values in time rather than nodal values. In this approach, a grid node xiis assu medto be the center of a grid cell (xi−1/2, xi+1/2), which is taken as the control volume. We integratethe conservation law across this control volume to obtainZxi+1/2xi−1/2φt+ f(φ)xdx =ddtZxi+1/2xi−1/2φ dx + f(φi+1/2) − f(φi−1/2) = 0.If we letˆφidenote the total quantity of φ in the ithgrid cell, i.e.ˆφi=Zxi+1/2xi−1/2φ dxthen we can write this as(ˆφi)t+ f(φi+1/2) − f(φi−1/2) = 0. (1)We will refer to values computed at the xias grid point or cell center values, and values computedat the xi±1/2as half grid point, cell wall, or flux values. We also define the cell average value of φin the grid cell i asφi=1∆xˆφi=1∆xZxi+1/2xi−1/2φ dx.Equation (1) has the desired conservation form in that the rate of change of the cell average is adifference of fluxes. The difficulty with this formulation is that it requires transforming between cell1averages of φ (which are directly evolved in time by the scheme) and cell wall values of φ (which mustbe reconstructed) to evaluate the needed fluxes. We would like to avoid reconstructing pointwisevalues of φ from the cell average values. The distinction between cell average and midpoint valuescan be ignored for schemes whose accuracy is no h igher than second order, since the cell averageand the midpoint value differ by only O(∆x2). This can be seen if we write φ in terms of its Taylorseries expansion about the point xiφ(x) = φ(xi) + (x − xi)φ′(xi) +(x − xi)22φ′′(xi) + · · · .Then,φi=1(xi+1/2− xi−1/2)Zxi+1/2xi−1/2φ(x) dx=1∆xZxi+1/2xi−1/2φ(xi) + (x − xi)φ′(xi) +(x − xi)22φ′′(xi) + · · ·dx=1∆x"φ(xi)∆x +(x − xi)22xi+1/2xi−1/2φ′(xi) +(x − xi)36xi+1/2xi−1/2φ′′(xi) + · · ·#= φ(xi) +∆x224φ′′(xi) + O(∆x4).We also assume that we have a uniform grid, so thatxi+1/2− xi−1/2= ∆xi= ∆x.For i ∈ {1, . . . , m}, we h ave(ˆφi)t+ f(ui+1/2) − f(ui−1/2) = 0.Summing over i, the fluxes cancel except for the ones on either side of the domain, so we getmXi=1(ˆφi)t+ f(um+1/2) − f(u1/2) = 0,or, equivalently,mXi=1(φi∆xi)t+ f(um+1/2) − f(u1/2) = 0.When using the weak form of the conservation law we evolve cell average values of φ in time,but requ ir e pointwise values of φ at the half grid cells in order to evaluate the flux functions. Asnoted above, if we only wanted a s econd order accurate scheme, we could simply approximate thecell average value with the value of φ at the cell center. However, we would like to use the pointwisevalues of φ while still getting better than second order accuracy. To achieve this, we replace thephysical flux function w ith a numerical flux function. We define the numerical flux function Fsuch thatf(u)x=F (x + ∆x/2) − F (x − ∆x/2)∆x(2)2We call F the numerical flux since we require it in our numerical scheme, and also to distinguishit from the closely related “physical flux”, f (u). It is n ot obvious that the numerical flux functionexists, but from relationship (2) one can solve for its Taylor exp ansion to obtainF = f (u) −∆x224f(u)xx+7∆x45760f(u)xxxx− · · · .In summary, we start with the conservation lawut+ f(u)x= 0.Integrating over a grid cell, we have(φi∆x)t+ f(ui+1/2) − f(ui−1/2) = 0.Replacingφiwith the pointwise value uiwe make an O(∆x2) error(ui∆x)t+ f(ui+1/2) − f(ui−1/2) = O(∆x2).Introducing the numerical flux function instead of the physical flux fu nction eliminates the error(ui)t+Fxi+1/2− Fxi−1/2∆x= 0.This is the d esired conservation form.2 Constructing the Numerical Flux FunctionWe define the numerical flux function through the relationf(ui)x=Fxi+1/2− Fxi−1/2∆x.To obtain a convenient algorithm for computing this numerical flux function, we define h(x) im-plicitly through the following equationf(φ(x)) =1∆xZx+∆x/2x−∆x/2h(y) dyand note that taking a derivative on both sides of this equation yieldsf(φ(x))x=h(x + ∆x/2) − h(x − ∆x/2)∆xwhich shows that h is identical to th e numerical flux function at the cell walls. That is Fi±1/2=h(xi±1/2) for all i. We calculate h by finding its primitiveH(x) =Zxx−1/2h(y)dy3using polynomial interpolation, and then take a derivative to get h. We build a divided differencetable to construct H.zeroth order D0i+1/2H at cell wallsfirst order D1iH at cell centerssecond orderD2i+1/2H at cell wallsthird order D3iH at cell centers.........That is, the even divided differences of H are at the cell walls, and the odddivided differences of H are at the cell centers. Since we are actually interested in determining h,we do not need the zeroth order divided differences of H as they will drop out when we differentiateto obtain h. Therefore, we can ignore the zeroth level of the divided difference table and constructit starting at the first level. The first level is given byD1iH =Hxi+1/2− Hxi−1/2∆x= f(ui) = D0if.This is becauseH(xi+1/2) =Zxi+1/2x−1/2h(y)dy=iXj=0 Zxj+1/2xj−1/2h(y)dy!= ∆xiXj=0f(φ(xj)).And similarly,H(xi−1/2) = ∆xi−1Xj=0f(φ(xj)).so thatH(xi+1/2) − H(xi−1/2) = ∆xf (φ(xi)).The higher divided differences areD2i+1/2H =f(φ(xi+1)) − f(φ(xi))2∆x=12D1i+1/2fD3iH =13D2ifcontinuing in that manner.According to the rules of polynomial interpolation, we can take any path along the divideddifference table to construct H, although they do not all give goo d results. ENO reconstructioncomprises two important