Gene Golub SIAM Summer School 2012 Numerical Methods for Wave - PowerPoint PPT Presentation

Gene Golub SIAM Summer School 2012 Numerical Methods for Wave Propagation Finite Volume Methods Lecture 2 Randall J. LeVeque Applied Mathematics University of Washington R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

Outline This lecture • Finite difference / finite volume methods • Godunov’s method • High resolution methods (limiters) • Two-dimensional methods • Seismic example R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

Upwind method for advection Scalar advection: q t + uq x = 0 , u > 0 As finite difference method: � Q n +1 � � Q n � − Q n i − Q n i i − 1 i + u = 0 ∆ t ∆ x Gives the explicit method: i − u ∆ t Q n +1 = Q n ∆ x ( Q n i − Q n i − 1 ) . i Stable provided CFL condition satisfied: 0 ≤ u ∆ t ∆ x ≤ 1 and first order accurate on smooth data. R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

The CFL Condition (Courant-Friedrichs-Lewy) Domain of dependence: The solution q ( X, T ) depends on the data q ( x, 0) over some set of x values, x ∈ D ( X, T ) . Advection: q ( X, T ) = q ( X − uT, 0) and so D ( X, T ) = { X − uT } . The CFL Condition: A numerical method can be convergent only if its numerical domain of dependence contains the true domain of dependence of the PDE, at least in the limit as ∆ t and ∆ x go to zero. Note: Necessary but not sufficient for stability! R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

Numerical domain of dependence With a 3-point explicit method: On a finer grid with ∆ t/ ∆ x fixed: R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

The CFL Condition For the method to be stable, the numerical domain of dependence must include the true domain of dependence. For advection, the solution is constant along characteristics, q ( x, t ) = q ( x − ut, 0) � � � u ∆ t � ≤ 1 . For a 3-point method, CFL condition requires ∆ x If this is violated: True solution is determined by data at a point x − ut that is ignored by the numerical method, even as the grid is refined. R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

The CFL Condition For the method to be stable, the numerical domain of dependence must include the true domain of dependence. For advection, the solution is constant along characteristics, q ( x, t ) = q ( x − ut, 0) � � � u ∆ t � ≤ 1 . For a 3-point method, CFL condition requires ∆ x If this is violated: True solution is determined by data at a point x − ut that is ignored by the numerical method, even as the grid is refined. R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

The CFL Condition For the method to be stable, the numerical domain of dependence must include the true domain of dependence. For advection, the solution is constant along characteristics, q ( x, t ) = q ( x − ut, 0) � � � u ∆ t � ≤ 1 . For a 3-point method, CFL condition requires ∆ x If this is violated: True solution is determined by data at a point x − ut that is ignored by the numerical method, even as the grid is refined. R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

The CFL Condition For the method to be stable, the numerical domain of dependence must include the true domain of dependence. For advection, the solution is constant along characteristics, q ( x, t ) = q ( x − ut, 0) � � � u ∆ t � ≤ 1 . For a 3-point method, CFL condition requires ∆ x If this is violated: True solution is determined by data at a point x − ut that is ignored by the numerical method, even as the grid is refined. R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

Stencil CFL Condition 0 ≤ u ∆ t ∆ x ≤ 1 − 1 ≤ u ∆ t ∆ x ≤ 0 − 1 ≤ u ∆ t ∆ x ≤ 1 0 ≤ u ∆ t ∆ x ≤ 2 −∞ < u ∆ t ∆ x < ∞ R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

Upwind method for advection Scalar advection: q t + uq x = 0 , u > 0 As finite difference method: � Q n +1 � � Q n � − Q n i − Q n i i − 1 i + u = 0 ∆ t ∆ x Gives the explicit method: i − u ∆ t Q n +1 = Q n ∆ x ( Q n i − Q n i − 1 ) . i Stable provided CFL condition satisfied: 0 ≤ u ∆ t ∆ x ≤ 1 and first order accurate on smooth data. R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

Finite differences vs. finite volumes Finite difference Methods • Pointwise values Q n i ≈ q ( x i , t n ) • Approximate derivatives by finite differences • Assumes smoothness Finite volume Methods � x i +1 / 2 1 • Approximate cell averages: Q n i ≈ q ( x, t n ) dx ∆ x x i − 1 / 2 • Integral form of conservation law, � x i +1 / 2 ∂ q ( x, t ) dx = f ( q ( x i − 1 / 2 , t )) − f ( q ( x i +1 / 2 , t )) ∂t x i − 1 / 2 leads to conservation law q t + f x = 0 but also directly to numerical method. R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

Finite volume method � x i +1 / 2 Q n i ≈ 1 x i − 1 / 2 q ( x, t n ) dx h Integral form: � x i +1 / 2 ∂ q ( x, t ) dx = f ( q ( x i − 1 / 2 , t )) − f ( q ( x i +1 / 2 , t )) ∂t x i − 1 / 2 R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

Finite volume method � x i +1 / 2 Q n i ≈ 1 x i − 1 / 2 q ( x, t n ) dx h Integral form: � x i +1 / 2 ∂ q ( x, t ) dx = f ( q ( x i − 1 / 2 , t )) − f ( q ( x i +1 / 2 , t )) ∂t x i − 1 / 2 Integrate from t n to t n +1 = ⇒ � � � t n +1 q ( x, t n +1 ) dx = q ( x, t n ) dx + f ( q ( x i − 1 / 2 , t )) − f ( q ( x i +1 / 2 , t )) dt t n R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

Finite volume method � x i +1 / 2 Q n i ≈ 1 x i − 1 / 2 q ( x, t n ) dx h Integral form: � x i +1 / 2 ∂ q ( x, t ) dx = f ( q ( x i − 1 / 2 , t )) − f ( q ( x i +1 / 2 , t )) ∂t x i − 1 / 2 Integrate from t n to t n +1 = ⇒ � � � t n +1 q ( x, t n +1 ) dx = q ( x, t n ) dx + f ( q ( x i − 1 / 2 , t )) − f ( q ( x i +1 / 2 , t )) dt t n i − ∆ t Q n +1 = Q n ∆ x ( F n i +1 / 2 − F n Numerical method: i − 1 / 2 ) i � t n +1 i − 1 / 2 ≈ 1 Numerical flux: F n f ( q ( x i − 1 / 2 , t )) dt . ∆ t t n R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

Upwind method for advection Flux: f ( q ) = uq � t n +1 i − 1 / 2 ≈ 1 F n Numerical flux: f ( q ( x i − 1 / 2 , t )) dt . ∆ t t n If q ( x, t n ) is piecewise constant in each cell, then � uQ n if u > 0 , F n i − 1 i − 1 / 2 = uQ n if u < 0 . i R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

Upwind method for advection Flux: f ( q ) = uq � t n +1 i − 1 / 2 ≈ 1 F n Numerical flux: f ( q ( x i − 1 / 2 , t )) dt . ∆ t t n If q ( x, t n ) is piecewise constant in each cell, then � uQ n if u > 0 , F n i − 1 i − 1 / 2 = uQ n if u < 0 . i This gives the upwind method: i − u ∆ t Q n +1 = Q n ∆ x ( Q n i − Q n i − 1 ) if u > 0 i i − u ∆ t Q n +1 = Q n ∆ x ( Q n i +1 − Q n i ) if u < 0 i R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

Upwind for advection as a finite volume method i − ∆ t Q n +1 = Q n ∆ x ( F n i +1 / 2 − F n i − 1 / 2 ) i Advection equation: f ( q ) = uq � t n +1 F i − 1 / 2 ≈ 1 uq ( x i − 1 / 2 , t ) dt. ∆ t t n First order upwind: F i − 1 / 2 = u + Q n i − 1 + u − Q n i i − ∆ t Q n +1 = Q n ∆ x ( u + ( Q n i − Q n i − 1 ) + u − ( Q n i +1 − Q n i )) . i where u + = max( u, 0) , u − = min( u, 0) . R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

Generalize upwind to a linear system? Consider q t + Aq x = 0 . Eigenvalues are wave speeds. Upwind method if all λ p > 0 : i − ∆ t Q n +1 = Q n ∆ x ( AQ n i − AQ n i − 1 ) i R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

Generalize upwind to a linear system? Consider q t + Aq x = 0 . Eigenvalues are wave speeds. Upwind method if all λ p > 0 : i − ∆ t Q n +1 = Q n ∆ x ( AQ n i − AQ n i − 1 ) i Upwind method if all λ p < 0 : i − ∆ t Q n +1 = Q n ∆ x ( AQ n i +1 − AQ n i ) i R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

Generalize upwind to a linear system? Consider q t + Aq x = 0 . Eigenvalues are wave speeds. Upwind method if all λ p > 0 : i − ∆ t Q n +1 = Q n ∆ x ( AQ n i − AQ n i − 1 ) i Upwind method if all λ p < 0 : i − ∆ t Q n +1 = Q n ∆ x ( AQ n i +1 − AQ n i ) i What if some eigenvalues of each sign? R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012

Generalize upwind to a linear system? Consider q t + Aq x = 0 . Eigenvalues are wave speeds. Upwind method if all λ p > 0 : i − ∆ t Q n +1 = Q n ∆ x ( AQ n i − AQ n i − 1 ) i Upwind method if all λ p < 0 : i − ∆ t Q n +1 = Q n ∆ x ( AQ n i +1 − AQ n i ) i What if some eigenvalues of each sign? Diagonalize and apply scalar upwind to each wave family. Easier ways to accomplish this! R.J. LeVeque, University of Washington Gene Golub SIAM Summer School 2012