Outline Notes: Godunovs method for acoustics Riemann solvers in - PDF document

Outline Notes: • Godunov’s method for acoustics • Riemann solvers in Clawpack • Acoustics in heterogeneous media • CFL Condition • Intro to Python plotting tools Reading: Chapters 4 – 6 www.clawpack.org/users Clawpack documentation www.clawpack.org/users/plotting.html Plotting hints See also Python and plotting sections of AMath 583 notes at http://www.amath.washington.edu/~rjl/uwamath583s11 R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 Godunov (upwind) on acoustics Notes: Q n +1 i t n +1 t n Q n i q n ( x, t n ) = Q n Data at time t n : ˜ for x i − 1 / 2 < x < x i +1 / 2 i Solving Riemann problems for small ∆ t gives solution:  Q ∗ if x i − 1 / 2 − c ∆ t < x < x i − 1 / 2 + c ∆ t, i − 1 / 2   q n ( x, t n +1 ) = Q n ˜ if x i − 1 / 2 + c ∆ t < x < x i +1 / 2 − c ∆ t, i Q ∗ if x i +1 / 2 − c ∆ t < x < x i +1 / 2 + c ∆ t,   i +1 / 2 So computing cell average gives: 1 � � Q n +1 i − 1 / 2 + (∆ x − 2 c ∆ t ) Q n = c ∆ tQ ∗ i + c ∆ tQ ∗ . i i +1 / 2 ∆ x R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 3.8, 4.12] R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 3.8, 4.12] Godunov (upwind) on acoustics Notes: 1 � � Q n +1 i − 1 / 2 + (∆ x − 2 c ∆ t ) Q n = c ∆ tQ ∗ i + c ∆ tQ ∗ . i i +1 / 2 ∆ x Solve Riemann problems: i − 1 / 2 r 1 + α 2 Q n i − Q n i − 1 = ∆ Q i − 1 / 2 = W 1 i − 1 / 2 + W 2 i − 1 / 2 = α 1 i − 1 / 2 r 2 , i +1 / 2 r 1 + α 2 Q n i +1 − Q n i = ∆ Q i +1 / 2 = W 1 i +1 / 2 + W 2 i +1 / 2 = α 1 i +1 / 2 r 2 , The intermediate states are: i − 1 / 2 = Q n i − W 2 i +1 / 2 = Q n i + W 1 Q ∗ i − 1 / 2 , Q ∗ i +1 / 2 , So, 1 � � Q n +1 c ∆ t ( Q n i − W 2 i − 1 / 2 ) + (∆ x − 2 c ∆ t ) Q n i + c ∆ t ( Q n i + W 1 = i +1 / 2 ) i ∆ x i − c ∆ t i − 1 / 2 + c ∆ t = Q n ∆ x W 2 ∆ x W 1 i +1 / 2 . R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 3.8, 4.12] R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 3.8, 4.12]

Godunov (upwind) on acoustics Notes: Solve Riemann problems: i − 1 / 2 r 1 + α 2 Q n i − Q n i − 1 = ∆ Q i − 1 / 2 = W 1 i − 1 / 2 + W 2 i − 1 / 2 = α 1 i − 1 / 2 r 2 , i +1 / 2 r 1 + α 2 Q n i +1 − Q n i = ∆ Q i +1 / 2 = W 1 i +1 / 2 + W 2 i +1 / 2 = α 1 i +1 / 2 r 2 , The waves are determined by solving for α from Rα = ∆ Q : � − Z � − 1 � 0 K � Z � R − 1 = 1 Z � A = , R = , . 1 /ρ 0 1 1 2 Z 1 Z So � ∆ p � − Z � Z � � � = α 1 + α 2 ∆ Q = ∆ u 1 1 with α 1 = 1 α 2 = 1 2 Z ( − ∆ p + Z ∆ u ) , 2 Z (∆ p + Z ∆ u ) . R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 4.12] R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 4.12] CLAWPACK Riemann solver Notes: The hyperbolic problem is specified by the Riemann solver • Input: Values of q in each grid cell • Output: Solution to Riemann problem at each interface. • Waves W p ∈ l R m , p = 1 , 2 , . . . , M w • Speeds s p ∈ l R , p = 1 , 2 , . . . , M w , R m • Fluctuations A − ∆ Q, A + ∆ Q ∈ l Note: Number of waves M w often equal to m (length of q ), but could be different (e.g. HLL solver has 2 waves). Fluctuations: A − ∆ Q = Contribution to cell average to left , A + ∆ Q = Contribution to cell average to right For conservation law, A − ∆ Q + A + ∆ Q = f ( Q r ) − f ( Q l ) R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Chap. 5] R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Chap. 5] CLAWPACK Riemann solver Notes: Inputs to rp1 subroutine: Value of q at left edge of i th cell, ql(i,1:m) = Value of q at right edge of i th cell, qr(i,1:m) = Warning: The Riemann problem at the interface between cells i − 1 and i has left state qr(i-1,:) and right state ql(i,:) . rp1 is normally called with ql = qr = q , but designed to allow other methods: R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Chap. 5] R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Chap. 5]

CLAWPACK Riemann solver Notes: Outputs from rp1 subroutine: for system of m equations with mw ranging from 1 to M w = # of waves Speed of wave # mw in i th Riemann solution, s(i,mw) = Jump across wave # mw , wave(i,1:m,mw) = Left-going fluctuation, updates Q i − 1 amdq(i,1:m) = Right-going fluctuation, updates Q i apdq(i,1:m) = R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Chap. 5] R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Chap. 5] Clawpack acoustics examples Notes: Constant coefficient acoustics: $CLAW/apps/acoustics/1d/example2/ Riemann solver: rp1.f R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 4.12] R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 4.12] Coupled advection–acoustics Notes: Flow in pipe with constant background velocity ¯ u . φ ( x, t ) = concentration of advected tracer u ( x, t ) , p ( x, t ) = acoustic velocity / pressure perturbation Equations include advection at velocity ¯ u : p t + up x ¯ + Ku x = 0 u t + (1 /ρ ) p x + uu x ¯ = 0 φ t + ¯ uφ x = 0 This is a linear system q t + Aq x = 0 with     p ¯ u K 0  ,  . q = u A = 1 /ρ ¯ u 0   φ 0 0 ¯ u R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 3.10] R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 3.10]

Coupled advection–acoustics Notes:  p   u ¯ K 0   ,  . q = u A = 1 /ρ ¯ u 0   φ 0 0 ¯ u λ 1 = u − c, λ 2 = u λ 3 = u + c, eigenvalues:       − Z 0 Z r 1 = r 2 = r 3 =  ,  ,  , eigenvectors: 1 0 1    0 1 0 κ/ρ, Z = ρc = √ ρκ . � where c =     − Z 0 Z − 1 Z 0 R − 1 = 1  ,  . R = 1 0 1 0 0 1   2 Z 0 1 0 1 Z 0 R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 3.10] R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 3.10] Coupled advection–acoustics Notes: Wave structure of solution in the x – t plane With no advection: R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 3.10] R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 3.10] Coupled advection–acoustics Notes: Wave structure of solution in the x – t plane Subsonic case ( | u 0 | < c ): R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 3.10] R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 3.10]

Coupled advection–acoustics Notes: Wave structure of solution in the x – t plane Supersonic case ( | u 0 | > c ): R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 3.10] R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 3.10] Wave propagation in heterogeneous medium Notes: Linear system q t + A ( x ) q x = 0 . For acoustics: � � 0 K ( x ) A = . 1 /ρ ( x ) 0 λ 1 = − c ( x ) , λ 2 = + c ( x ) , eigenvalues: � where c ( x ) = κ ( x ) /ρ ( x ) = local speed of sound. � − Z ( x ) � Z ( x ) � � r 1 ( x ) = r 2 ( x ) = eigenvectors: , 1 1 where Z ( x ) = ρc = √ ρκ = impedance. � − Z ( x ) � − 1 � 1 � Z ( x ) Z ( x ) R − 1 ( x ) = R ( x ) = , . 1 1 1 Z ( x ) 2 Z ( x ) Cannot diagonalize unless Z ( x ) is constant. R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 9.6] R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 9.6] Wave propagation in heterogeneous medium Notes: Multiply system q t + A ( x ) q x = 0 by R − 1 ( x ) on left to obtain R − 1 ( x ) q t + R − 1 ( x ) A ( x ) R ( x ) R − 1 ( x ) q x = 0 or ( R − 1 ( x ) q ) t + Λ( x ) � ( R − 1 ( x ) q ) x − R − 1 � x ( x ) q = 0 Let w ( x, t ) = R − 1 ( x ) q ( x, t ) (characteristic variable). There is a coupling term on the right: w t + Λ( x ) w x = Λ( x ) R − 1 x ( x ) R ( x ) w ⇒ reflections (unless R − 1 = x ( x ) ≡ 0 ). R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 9.7, 9.8] R.J. LeVeque, University of Washington IPDE 2011, June 27, 2011 [FVMHP Sec. 9.7, 9.8]