Outline Lab on multidimensional problems • Acoustics • Riemann solvers R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011
Wave propagation algorithms in 2D Clawpack requires: Normal Riemann solver rpn2.f Solves 1d Riemann problem q t + Aq x = 0 Decomposes ∆ Q = Q ij − Q i − 1 ,j into A + ∆ Q and A − ∆ Q . For q t + Aq x + Bq y = 0 , split using eigenvalues, vectors: A = R Λ R − 1 = ⇒ A − = R Λ − R − 1 , A + = R Λ + R − 1 Input parameter ixy determines if it’s in x or y direction. In latter case splitting is done using B instead of A . This is all that’s required for dimensional splitting. R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011 [FVMHP Chap. 19]
Wave propagation algorithms in 2D Clawpack requires: Normal Riemann solver rpn2.f Solves 1d Riemann problem q t + Aq x = 0 Decomposes ∆ Q = Q ij − Q i − 1 ,j into A + ∆ Q and A − ∆ Q . For q t + Aq x + Bq y = 0 , split using eigenvalues, vectors: A = R Λ R − 1 = ⇒ A − = R Λ − R − 1 , A + = R Λ + R − 1 Input parameter ixy determines if it’s in x or y direction. In latter case splitting is done using B instead of A . This is all that’s required for dimensional splitting. Transverse Riemann solver rpt2.f Decomposes A + ∆ Q into B − A + ∆ Q and B + A + ∆ Q by splitting this vector into eigenvectors of B . (Or splits vector into eigenvectors of A if ixy=2 .) R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011 [FVMHP Chap. 19]
Wave propagation algorithm for q t + Aq x + Bq y = 0 Decompose A = A + + A − and B = B + + B − . For ∆ Q = Q ij − Q i − 1 ,j : R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011 [FVMHP Chap. 19]
Wave propagation algorithm for q t + Aq x + Bq y = 0 Decompose A = A + + A − and B = B + + B − . For ∆ Q = Q ij − Q i − 1 ,j : R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011 [FVMHP Chap. 19]
Wave propagation algorithm for q t + Aq x + Bq y = 0 Decompose A = A + + A − and B = B + + B − . For ∆ Q = Q ij − Q i − 1 ,j : R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011 [FVMHP Chap. 19]
Wave propagation algorithm for q t + Aq x + Bq y = 0 Decompose A = A + + A − and B = B + + B − . For ∆ Q = Q ij − Q i − 1 ,j : R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011 [FVMHP Chap. 19]
Wave propagation algorithm for q t + Aq x + Bq y = 0 Decompose A = A + + A − and B = B + + B − . For ∆ Q = Q ij − Q i − 1 ,j : R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011 [FVMHP Chap. 19]
Acoustics in 2 dimensions p t + K 0 ( u x + v y ) = 0 ρ 0 u t + p x = 0 ρ 0 v t + p y = 0 Note: pressure responds to compression or expansion and so p t is proportional to divergence of velocity. Second and third equations are F = ma . R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011 [FVMHP Chap. 18]
Acoustics in 2 dimensions p t + K 0 ( u x + v y ) = 0 ρ 0 u t + p x = 0 ρ 0 v t + p y = 0 Note: pressure responds to compression or expansion and so p t is proportional to divergence of velocity. Second and third equations are F = ma . Gives hyperbolic system q t + Aq x + Bq y = 0 with p 0 K 0 0 0 0 K 0 , , . u 1 /ρ 0 0 0 0 0 0 q = A = B = v 0 0 0 1 /ρ 0 0 0 R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011 [FVMHP Chap. 18]
Acoustics in 2 dimensions p 0 K 0 0 0 0 K 0 , , . q = u A = 1 /ρ 0 0 0 B = 0 0 0 v 0 0 0 1 /ρ 0 0 0 Plane waves: 0 K 0 cos θ K 0 sin θ . A cos θ + B sin θ = cos θ/ρ 0 0 0 sin θ/ρ 0 0 0 R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011 [FVMHP Chap. 18]
Acoustics in 2 dimensions p 0 K 0 0 0 0 K 0 , , . q = u A = 1 /ρ 0 0 0 B = 0 0 0 v 0 0 0 1 /ρ 0 0 0 Plane waves: 0 K 0 cos θ K 0 sin θ . A cos θ + B sin θ = cos θ/ρ 0 0 0 sin θ/ρ 0 0 0 Eigenvalues: λ 1 = − c 0 , λ 2 = 0 , λ 3 = + c 0 where � c 0 = K 0 /ρ 0 Independent of angle θ . Isotropic: sound propagates at same speed in any direction. R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011 [FVMHP Chap. 18]
Acoustics in 2 dimensions p 0 K 0 0 0 0 K 0 , , . q = u A = 1 /ρ 0 0 0 B = 0 0 0 v 0 0 0 1 /ρ 0 0 0 Plane waves: 0 K 0 cos θ K 0 sin θ . A cos θ + B sin θ = cos θ/ρ 0 0 0 sin θ/ρ 0 0 0 Eigenvalues: λ 1 = − c 0 , λ 2 = 0 , λ 3 = + c 0 where � c 0 = K 0 /ρ 0 Independent of angle θ . Isotropic: sound propagates at same speed in any direction. Note: Zero wave speed for “shear wave” with variation only in velocity in direction ( − sin θ, cos θ ) . (Fig 18.1) R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011 [FVMHP Chap. 18]
Diagonalization 2 dimensions Can we diagonalize system q t + Aq x + Bq y = 0 ? Only if A and B have the same eigenvectors! If A = R Λ R − 1 and B = RMR − 1 , then let w = R − 1 q and w t + Λ w x + Mw y = 0 R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011 [FVMHP Chap. 18]
Diagonalization 2 dimensions Can we diagonalize system q t + Aq x + Bq y = 0 ? Only if A and B have the same eigenvectors! If A = R Λ R − 1 and B = RMR − 1 , then let w = R − 1 q and w t + Λ w x + Mw y = 0 This decouples into scalar advection equations for each component of w : w p t + λ p w p x + µ p w p ⇒ w p ( x, y, t ) = w p ( x − λ p t, y − µ p t, 0) . y = 0 = Note: In this case information propagates only in a finite number of directions ( λ p , µ p ) for p = 1 , . . . , m . R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011 [FVMHP Chap. 18]
Diagonalization 2 dimensions Can we diagonalize system q t + Aq x + Bq y = 0 ? Only if A and B have the same eigenvectors! If A = R Λ R − 1 and B = RMR − 1 , then let w = R − 1 q and w t + Λ w x + Mw y = 0 This decouples into scalar advection equations for each component of w : w p t + λ p w p x + µ p w p ⇒ w p ( x, y, t ) = w p ( x − λ p t, y − µ p t, 0) . y = 0 = Note: In this case information propagates only in a finite number of directions ( λ p , µ p ) for p = 1 , . . . , m . This is not true for most coupled systems, e.g. acoustics. R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011 [FVMHP Chap. 18]
Acoustics in 2 dimensions p t + K 0 ( u x + v y ) = 0 ρ 0 u t + p x = 0 ρ 0 v t + p y = 0 0 K 0 0 − Z 0 0 Z 0 R x = , A = 1 /ρ 0 0 0 1 0 1 0 0 0 0 1 0 Solving q t + Aq x = 0 gives pressure waves in ( p, u ) . x -variations in v are stationary. R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011 [FVMHP Chap. 18]
Acoustics in 2 dimensions p t + K 0 ( u x + v y ) = 0 ρ 0 u t + p x = 0 ρ 0 v t + p y = 0 0 K 0 0 − Z 0 0 Z 0 R x = , A = 1 /ρ 0 0 0 1 0 1 0 0 0 0 1 0 Solving q t + Aq x = 0 gives pressure waves in ( p, u ) . x -variations in v are stationary. 0 0 K 0 − Z 0 0 Z 0 R y = B = 0 0 0 0 1 0 1 /ρ 0 0 0 1 0 1 Solving q t + Bq y = 0 gives pressure waves in ( p, v ) . y -variations in u are stationary. R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011 [FVMHP Chap. 18]
Storing data in aux arrays In Clawpack, q(i,j,m), m=1,...,meqn holds the solution. Often there is spatially varying data that describes the problem: • Edge velocities for advection, • Density ρ 0 ( x, y ) , bulk modulus K 0 ( x, y ) for acoustics, • Topography or bathymetry for shallow water. • Edge lengths, angles, and cell areas for mapped grids, R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011
Storing data in aux arrays In Clawpack, q(i,j,m), m=1,...,meqn holds the solution. Often there is spatially varying data that describes the problem: • Edge velocities for advection, • Density ρ 0 ( x, y ) , bulk modulus K 0 ( x, y ) for acoustics, • Topography or bathymetry for shallow water. • Edge lengths, angles, and cell areas for mapped grids, These can be stored in aux(i,j,m), m=1,2,...,maux . The Fortran function setaux is called every time a new grid is created (when AMR is used). To use this, copy library version (which does nothing) to application directory and modify this file and Makefile . R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011
Using b4stepN.f The setaux function is only called when grids are created. The b4stepN function (in N dimensions) is called before each time step. Can use this for example to: • Change aux arrays for time-dependent velocities, • Print something out every time step (e.g. total mass), To use this, copy library version (which does nothing) to application directory and modify this file and Makefile . See: $CLAW/apps/advection/2d/swirl/b4step2.f $CLAW/apps/advection/2d/swirl/setaux.f $CLAW/apps/advection/2d/swirl/psi.f R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011
Acoustics in heterogeneous media q = ( p, u, v ) T , q t + A ( x, y ) q x + B ( x, y ) q y = 0 , where 0 K ( x, y ) 0 0 0 K ( x, y ) , . A = 1 /ρ ( x, y ) 0 0 B = 0 0 0 0 0 0 1 /ρ ( x, y ) 0 0 Note: Not in conservation form! R.J. LeVeque, University of Washington IPDE 2011, July 5, 2011 [FVMHP Sec. 21.5]
Recommend
More recommend
