Wave-propagation based
generalized Lagrangian method for
hyperbolic balance laws
Application to inviscid compressible flow
Keh-Ming Shyue
Department of Mathematics National Taiwan University
Taiwan
Objective
Describe simple Lagrangian-like moving grid approach for numerical resolution of nonlinear hyperbolic balance laws of the form
∂
∂tq (~x, t) +
N
X
j=1
∂
∂xj fj (q, ~x) = ψ (q, ~x)
with discontinuous initial data in general N ≥ 1 rectangular or non-rectangular geometry
~x = (x1, x2, . . . , xN): spatial vector, t: time q ∈ lRm: vector of m state quantities
fj ∈ lRm: flux vector, j = 1, 2, . . . , N, ψ ∈ lRm: source term Model equation is hyperbolic if PNj=1d αj (∂fj/∂q) is
diagonalizable with real eigenvalues, αj ∈ lR
Outline
Mathematical model for general balance laws Eulerian formulation
Generalized Lagrangian formulation
Example to single component compressible flow Wave-propagation based finite volume methods
Generalized Riemann problem & approximate solver Flux-based wave decomposition method
Sample numerical examples
Extension to compressible two-phase flow Future work
Mathematical Model
To begin with, consider a general non-rectangular domain Ω (N = 2 shown below) & introduce coordinate change
(~x, t) 7→ (~ξ, τ ) via
ξ = (ξ~ 1, ξ2, . . . , ξN) , ξj = ξj(~x, t), τ = t,
that maps a physical domain Ω to a logical one Ωˆ
−1 0 1
−1.5
−1
−0.5 0
−1 −0.5 0 0.5
−2
−1.5
−1
−0.5
x1 x2
ξ1
ξ2
Ω −→ Ωˆ
mapping
ξ1 = ξ1(x1, x2) ξ2 = ξ2(x1, x2)
logical domain physical domain
Mathematical Model (Cont.)
To derive hyperbolic balance laws in this generalized
coordinate (~ξ, τ ), using chain rule of partial differentiation, derivatives in physical space become
∂
∂t = ∂
∂τ +
N
X
i=1
∂ξi
∂t
∂
∂ξi, ∂
∂xj =
N
X
i=1
∂ξi
∂xj
∂
∂ξi for j = 1, 2, . . . , N ,
yielding the equation
∂q
∂τ +
N
X
j=1
∂ξj
∂t
∂q
∂ξj +
N
X
i=1
∂ξi
∂xj
∂fj
∂ξi
!
= ψ(q)
Note this is not in divergence form, and hence is not conservative, in case the source term ψ is ignored.
Mathematical Model (Cont.)
To obtain a strong conservation-law form as
∂ ˜q
∂τ +
N
X
j=1
∂ ˜fj
∂ξj = ˜ψ
for some q˜, f˜j, & ψ˜, we first multiply J = det
∂~ξ/∂~x−1
to the aforementioned non-conservative equations, and have
J ∂q
∂τ +
N
X
j=1
J ∂ξj
∂t
∂q
∂ξj +
N
X
i=1
∂ξi
∂xj
∂fj
∂ξi
!
= Jψ(q)
Then use differentiation by parts, u dv = d(uv) − v du, yielding
∂ ˜q
∂τ +
N
X
j=1
∂ ˜fj
∂ξj = ˜ψ + G with q = Jq˜ , f˜j = J
q∂ξ∂tj + PN
k=1 fk ∂x∂ξj
k
, ψ = Jψ,˜ & G (see next)
Mathematical Model (Cont.)
Here we have
G = q
∂J
∂τ +
N
X
j=1
∂
∂ξj
J ∂ξj
∂t
+
N
X
j=1
fj
" N X
k=1
∂
∂ξk
J ∂ξk
∂xj
#
Note with the use of basic grid-metric relations (see next), it is known that
∂J
∂τ +
N
X
j=1
∂
∂ξj
J ∂ξj
∂t
= 0 (geometrical conservation law)
N
X
k=1
∂
∂ξk
J ∂ξk
∂xj
= 0 ∀ j = 1, 2, . . . , N (compatibility condition)
and hence G = 0
Basic Grid-Metric Relations
Assume existence of inverse transformation
t = τ, xj = xj(~ξ, t) for j = 1, 2, . . . , N ,
To find basic grid-metric relations between different coordinates, employ elementary differential rule
∂(τ, ~ξ)
∂(t, ~x) = ∂(t, ~x)
∂(τ, ~ξ)
−1
yielding in N = 3 case, for example, as
1 0 0 0
∂tξ1 ∂x1ξ1 ∂x2ξ1 ∂x3ξ1
∂tξ2 ∂x1ξ2 ∂x2ξ2 ∂x3ξ2
∂tξ3 ∂x ξ3 ∂x ξ3 ∂x ξ3
= 1 J
J 0 0 0
J01 J11 J21 J31 J02 J12 J22 J32 J03 J13 J23 J33
Grid-Metric Relations (Cont.)
Here
J =
∂(x1, x2, x3)
∂(ξ1, ξ2, ξ3)
= det ∂(x1, x2, x3)
∂(ξ1, ξ2, ξ3)
, J11 =
∂(x2, x3)
∂(ξ2, ξ3)
, J21 =
∂(x1, x3)
∂(ξ3, ξ2)
, J31 =
∂(x1, x2)
∂(ξ2, ξ3)
, J12 =
∂(x2, x3)
∂(ξ3, ξ1)
, J22 =
∂(x1, x3)
∂(ξ1, ξ3)
, J32 =
∂(x1, x2)
∂(ξ3, ξ1)
, J13 =
∂(x2, x3)
∂(ξ1, ξ2)
, J23 =
∂(x1, x3)
∂(ξ2, ξ1)
, J33 =
∂(x1, x2)
∂(ξ1, ξ2)
,
J0j = −
Nd
X
i=1
Jij∂τxi, j = 1, 2, 3,
and so grid-metric relations between different coordinates
∇ξ = (∂ ξ , ∇ ξ ) = (∂ ξ , ∂ ξ , ∂ ξ , ∂ ξ ) = 1
(J , J , J , J )
Grid-Metric Relations (Cont.)
Note in two dimensions N = 2, we have
∂ξ1
∂t , ∂ξ1
∂x1 , ∂ξ1
∂x2
= 1 J
−∂x1
∂τ
∂x2
∂ξ2 + ∂x2
∂τ
∂x1
∂ξ2 , ∂x2
∂ξ2 , −∂x1
∂ξ2
∂ξ2
∂t , ∂ξ2
∂x1 , ∂ξ2
∂x2
= 1 J
∂x1
∂τ
∂x2
∂ξ1 − ∂x2
∂τ
∂x1
∂ξ1 , −∂x2
∂ξ1 , ∂x1
∂ξ1
J = ∂x1
∂ξ1
∂x2
∂ξ2 − ∂x1
∂ξ2
∂x2
∂ξ1
Thus to have G = 0 fulfilled, grid-metrics should obey
∂J
∂τ + ∂
∂ξ1
J ∂ξ1
∂t
+ ∂
∂ξ2
J ∂ξ2
∂t
= 0
∂
∂ξ1
J ∂ξ1
∂x1
+ ∂
∂ξ2
J ∂ξ2
∂x1
= ∂
∂ξ1
∂x2
∂ξ2
+ ∂
∂ξ2
−∂x2
∂ξ1
= 0
∂
∂ξ
J ∂ξ1
∂x
+ ∂
∂ξ
J ∂ξ2
∂x
= ∂
∂ξ
−∂x1
∂ξ
+ ∂
∂ξ
∂x1
∂ξ
= 0
Mathematical Model (Cont.)
As an example, with gravity effect included, Euler equations for single component compressible gas flow take
Cartesian coordinate case
∂
∂t
ρ ρui
E
+
N
X
j=1
∂
∂xj
ρuj
ρuiuj + pδij Euj + puj
=
0
−ρ∂x∂φ
i
−ρ~u · ∇φ
, i = 1, . . . , N
Generalized coordinate case
∂
∂τ
ρJ ρJui
JE
+
N
X
j=1
∂
∂ξj J
ρUj
ρuiUj + p∂ξ∂xj
i
EUj + pUj − p∂ξ∂tj
=
0
−ρJ ∂x∂φ
i
−ρJ~u · ∇φ
ρ: density, p = p(ρ, e): pressure , e: internal energy
E = ρe + ρPN
j=1 u2j/2: total energy, φ: gravitational potential
Mathematical Model (Cont.)
Note that to complete the model, we must
1. make clear the transformation (~x, t) 7→ (~ξ, τ ) initially
Depending on how complex the geometry is, this can be done by various numerical means
Mathematical Model (Cont.)
Note that to complete the model, we must
1. make clear the transformation (~x, t) 7→ (~ξ, τ ) initially
Depending on how complex the geometry is, this can be done by various numerical means
2. choose a moving grid strategy for ∂τ~x
Mathematical Model (Cont.)
Note that to complete the model, we must
1. make clear the transformation (~x, t) 7→ (~ξ, τ ) initially
Depending on how complex the geometry is, this can be done by various numerical means
2. choose a moving grid strategy for ∂τ~x
When ∂τ~x = 0 or ∂τ~x = ~ub(t) (rigid-body motion)
∂tξ~ & ∇~xξ~ are time-independent; no need to have more additional condition
Mathematical Model (Cont.)
Note that to complete the model, we must
1. make clear the transformation (~x, t) 7→ (~ξ, τ ) initially
Depending on how complex the geometry is, this can be done by various numerical means
2. choose a moving grid strategy for ∂τ~x
When ∂τ~x = 0 or ∂τ~x = ~ub(t) (rigid-body motion)
∂tξ~ & ∇~xξ~ are time-independent; no need to have more additional condition
While ∂τ~x 6= 0 (flow-dependent motion) (see next)
∂tξ~ & ∇~xξ~ would be time-dependent; require
additional conditions to determine ∇ξ~~x (N2 of them in total) over time (see below)
Lagrangian-Like Moving Grid
For compressible flow, to improve numerical resolution of interfaces (material or slip lines), it is popular to take ∂τ~x as
Lagrangian case: ∂τ~x = ~u (flow velocity)
Lagrangian-like case: ∂τ~x = h0~u (pseudo velocity) h0 ∈ [0, 1] (fixed piecewise const.)
Unified coordinate case: ∂τ~x = h~u
h ∈ [0, 1] but is determined from a PDE constraint arising from such as grid-angle or grid-Jacobian preserving condition
ALE-like case: ∂τ~x = ~U (arbitrary velocity)
Lagrangian-Like Moving Grid
For compressible flow, to improve numerical resolution of interfaces (material or slip lines), it is popular to take ∂τ~x as
Lagrangian case: ∂τ~x = ~u (flow velocity)
Lagrangian-like case: ∂τ~x = h0~u (pseudo velocity) h0 ∈ [0, 1] (fixed piecewise const.)
Unified coordinate case: ∂τ~x = h~u
h ∈ [0, 1] but is determined from a PDE constraint arising from such as grid-angle or grid-Jacobian preserving condition
ALE-like case: ∂τ~x = ~U (arbitrary velocity)
Here we will focus on the simple Lagrangian-like case
Unified Coordinate System
Consider N = 2 case, for example, and use simplified
notation ~x = (x, y), ξ = (ξ, η)~ . At given time instance, free parameter h can be chosen based on
Grid-angle preserving condition (Hui et al. JCP 1999)
∂
∂τ cos−1 ∇ξ
|∇ξ| · ∇η
|∇η|
= ∂
∂τ cos−1
−yηxη − yξxξ qyξ2 + yη2q
x2ξ + x2η
= · · ·
= Ahξ + Bhη + Ch = 0 (1st order PDE )
with
A = q
x2η + yη2 (vxξ − uyξ) , B = q
x2ξ + yξ2 (uyη − vxη) C = q
x2ξ + yξ2 (uηyη − vηxη) − q
x2η + yη2 (uξyξ − vξxξ)
Unified Coordinate System
Consider N = 2 case, for example, and use simplified
notation ~x = (x, y), ξ = (ξ, η)~ . Or alternatively, based on Grid-Jacobian preserving condition
∂J
∂τ = ∂
∂τ (xξyη − xηyξ)
= xξτ yη + xξ yητ − xητ yξ − xη yξτ
= · · ·
= Ahξ + Bhη + Ch = 0 (1st order PDE )
with
A = uyη − vxη, B = vxξ − uyξ, C = uξyη + vηxξ − uηyξ − vξxη
Lagrangian-Like Grid (Cont.)
Now with the temporal motion of the coordinate system governed by ∂τ~x = h0~u. We should impose conditions on grid metrics ∂tξ~ & ∇~xξ~ to have the fulfillment of geometrical conservation law
∂J
∂τ +
N
X
j=1
∂
∂ξj
J ∂ξj
∂t
= 0
Here we are interested in an approach that is based on the compatibility condition of ∂τ∂ξjxi & ∂ξj∂τxi, i.e.,
∂
∂τ
∂xi
∂ξj
+ ∂
∂ξj
−∂xi
∂τ
= 0 for i, j = 1, 2, . . . , N.
for unknowns ∂xi/∂ξj, yielding easy computation of J & ∇ξj
Mathematical Model (Cont.)
In summary, our Lagrangina-like model system for single component compressible flow problems consists of
Physical balance laws
∂
∂τ
ρJ ρJui
JE
+
N
X
j=1
∂
∂ξj J
ρUj
ρuiUj + p∂ξ∂xj
i
EUj + pUj − p∂ξ∂tj
=
0
−ρJ ∂x∂φ
i
−ρJ~u · ∇φ
Geometrical conservation laws
∂
∂τ
∂xi
∂ξj
+ ∂
∂ξj
−∂xi
∂τ
= 0 for i, j = 1, 2, . . . , N.
Moving grid condition ∂τ~x = h0~u & pressure law p(ρ, e)
Axisymmetric Compressible Flow
Physical balance laws (ξ1: axisymmetric direction)
∂
∂τ
ρJ ρJui
JE
+
2
X
j=1
∂
∂ξj J
ρUj
ρuiUj + p∂ξ∂xj
i
EUj + pUj − p∂ξ∂tj
=
−x1ρJu1
−x1ρJuiuj − ρJ ∂x∂φ
i
−x1J(E + p)u1 − ρJ~u · ∇φ
for i = 1, 2
Geometrical conservation laws
∂
∂τ
∂xi
∂ξj
+ ∂
∂ξj
−∂xi
∂τ
= 0 for i, j = 1, 2.
Moving grid condition ∂τ~x = h0~u & pressure law p(ρ, e)
Mathematical Models: Remarks
For single component compressible flow model mentioned above, it is known that under some thermodynamic stability conditions
when h0 = 0 (Eulerian case), the model is hyperbolic when h0 = 1 (Lagrangian case), the model is weakly hyperbolic (do not possess complete eigenvectors) when h0 ∈ (0, 1) (Lagrangian-like case), the model is hyperbolic
Mathematical Models: Remarks
For single component compressible flow model mentioned above, it is known that under some thermodynamic stability conditions
when h0 = 0 (Eulerian case), the model is hyperbolic when h0 = 1 (Lagrangian case), the model is weakly hyperbolic (do not possess complete eigenvectors) when h0 ∈ (0, 1) (Lagrangian-like case), the model is hyperbolic
If a prescribed velocity ~ub for a rigid body motion is included in the formulation i.e., with ∂τ~x = h0~u + ~ub, we should be able to use the model to solve some moving body problems as well.
Review of Previous Work
The work presented here is related to
W.H. Hui et al. (JCP 1999, 2001): Unified coordinated system for Euler equations
W.H. Hui (Comm. Phys. Sci. 2007): Unified coordinate system in CFD
C. Jin & K. Xu (JCP 2007): Moving grid gas-kinetic method for viscous flow
P. Jia et al. (Computers and Fluids 2006) Unified
coordinated system for compressible milti-material flow Z. Chen et al. (Int J. Numer. Meth Fluids 2007): Wave speed based moving coordinates for compressible flow equations
Numerical Methods
Employ finite volume formulation of numerical solution
Qnijk ≈ 1
∆ξ1∆ξ2∆ξ3 Z
Cijk
q(ξ1, ξ2, ξ3, τn) dV
that gives approximate value of cell average of solution q over cell Cijk at time τn (sample case in 2D shown below)
i − 1 i − 1
i
i j
j
j + 1 j + 1
Cij Cˆij
ξ1 ξ2
mapping
∆ξ1
∆ξ2 logical domain physical domain
←−
x1 = x1(ξ1, ξ2) x2 = x2(ξ1, ξ2) x1
x2
Numerical Methods (Cont.)
In three dimensions N = 3, equations to be solved take
∂
∂τ q ~ξ, τ + XN
j=1
∂
∂ξj fj
q, ~ξ
= ψ
q, ~ξ
A simple dimensional-splitting method based on f-wave approach of LeVeque et al. is used for approximation, i.e.,
Solve one-dimensional generalized Riemann problem (defined below) at each cell interfaces
Use resulting jumps of fluxes (decomposed into each wave family) of Riemann solution to update cell
averages
Introduce limited jumps of fluxes to achieve high
Numerical Methods (cont.)
Basic steps of a dimensional-splitting scheme ξ1-sweeps: solve
∂q
∂τ + f1 ∂
∂ξ, q, ∇~ξ
= 0 updating Qnijk to Q∗ijk
ξ2-sweeps: solve
∂q
∂τ + f2
∂
∂ξ2, q, ∇~ξ
= 0 updating Q∗ijk to Q∗∗ijk
ξ3-sweeps: solve
∂q
∂τ + f3
∂
∂ξ3 , q, ∇~ξ
= 0 updating Q∗∗ijk to Qn+1ijk
Numerical Methods (cont.)
Consider ξ1-sweeps, for example, First order update is
Q∗ijk = Qnijk − ∆τ
∆ξ1
h A+1 ∆Qn
i−1/2,jk + A−1 ∆Qn
i+1/2,jk
i
with the fluctuations
(A+1 ∆Q)ni−1/2,jk = X
m:(λ1,m)ni−1/2,jk>0
(Z1,m)ni−1/2,jk
and
(A−1 ∆Q)ni+1/2,jk = X
m:(λ1,m)ni+1/2,jk<0
(Z1,m)ni+1/2,jk
(λ1,m)nι−1/2,jk & (Z1,m)nι−1/2,jk are in turn wave speed and f-waves
Numerical Methods (cont.)
High resolution correction is
Q∗ijk := Q∗ijk − ∆τ
∆ξ1
˜F1n
i+1/2,jk − ˜F1n
i−1/2,jk
with ( ˜F1)ni−1/2,jk = 1 2
mw
X
m=1
sign(λ1,m)
1 − ∆τ
∆ξ1 |λ1,m|
Z˜1,m
n
i−1/2,jk
Z˜ι,m is a limited value of Zι,m
It is clear that this method belongs to a class of upwind schemes, and is stable when the typical CFL (Courant-Friedrichs-Lewy) condition:
ν = ∆τ maxm (λ1,m, λ2,m, λ3,m)
min (∆ξ1, ∆ξ2, ∆ξ3) ≤ 1,
Generalized Riemann Problem
Generalized Riemann problem for our model equations at cell interface ξi−1/2 consists of the equation
∂qi−1,jk
∂τ + f1
∂
∂ξ1 , qi−1,jk
= 0 if ξ1 < (ξ1)i−1/2,
∂qijk
∂τ + f1
∂
∂ξ1 , qijk
= 0 if ξ1 > (ξ1)i−1/2,
together with piecewise constant initial data
q(ξ1, 0) =
Qni−1,jk for ξ < ξi−1/2 Qnijk for ξ > ξi−1/2
qijk = q|(∂ξ2~x, ∂ξ3~x)ijk & f1(∂ξ1, qijk) = f1(∂ξ1, q)|(∂ξ2~x, ∂ξ3~x)ijk
Generalized Riemann Problem
Generalized Riemann problem at time τ = 0
∂τqL + f1 (∂ξ1, qL) = 0 ∂τqR + f1 (∂ξ1, qR) = 0 ξ1 τ
Qni−1,jk Qnijk
Generalized Riemann Problem
Exact generalized Riemann solution: basic structure
∂ q + f (∂ , q ) = 0 ∂ q + f (∂ , q ) = 0 ξ1 τ
Qni−1,jk Qnijk
Generalized Riemann Problem
Shock-only approximate Riemann solution: basic structure
∂τqL + f1 (∂ξ1, qL) = 0 ∂τqR + f1 (∂ξ1, qR) = 0 ξ1 τ
Qni−1,jk Qnijk
λ0
λ1 λ2
λ3 qmL− qmL+
qmR
Z1 = fL(qmL− ) − fL(Qni−1,jk) Z2 = fR(qmR) − fR(qmL+ )
Z3 = fR(Qnijk) − fR(qmR)
Shock-only Riemann Solver
Rotate velocity vector in Riemann data normal to each cell interface Find midstate velocity υm and pressure pm by solving
φ(pm) = υmR(pm) − υmL(pm) = 0
derived from Rankine-Hugoniot relation iteratively, where
υmL(p) = υL − p − pL
ML(p), υmR(p) = υR + p − pR MR(p)
Propagation speed of each moving discontinuity is determined by
(λ1,1)i−1/2,jk =
(1 − h0)υm − ML(pm) ρmL(pm)
∇X~ ξ1
i−1/2,jk
(λ1,2)i−1/2,jk = (1 − h0)υm
∇X~ξ1
i−1/2,jk
(λ ) =
(1 − h )υ + MR(pm)
∇ ξ
Lax’s Riemann Problem
Ideal gas EOS with γ = 1.4 h0 = 0 Eulerian result
h0 = 0.99 Lagrangian-like result
sharper resolution for contact discontinuity
0 0.5 1
0 0.5 1 1.5 2
x
ρ
Exact h0=0 h0=0.99
0 0.5 1
0 0.5 1 1.5 2
x
u
0 0.5 1
0 1 2 3 4
x
p
Lax’s Riemann Problem
Physical grid coordinates at selected times
Each little dashed line gives a cell-center location of the proposed Lagrange-like grid system
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0 0.05 0.1 0.15
x
time
Woodward-Colella’s Problem
Ideal gas EOS with γ = 1.4 h0 = 0 Eulerian result
h0 = 0.99 Lagrangian-like result
sharper resolution for contact discontinuity
0 0.5 1
0 2 4 6 8
ρ
0 0.5 1
−10 0 10 20
0 0.5 1
0 100 200 300 400
u p
t = 0.016 t = 0.016
t = 0.016
Woodward-Colella’s Problem
h0 = 0 Eulerian result
h0 = 0.99 Lagrangian-like result
sharper resolution for contact discontinuity
0 0.5 1
0 5 10 15 20
ρ
Fine grid h0=0 h0=0.99
0 0.5 1
−5 0 5 10 15
0 0.5 1
0 200 400 600
u p
t = 0.032 t = 0.032
t = 0.032
Woodward-Colella’s Problem
h0 = 0 Eulerian result
h0 = 0.99 Lagrangian-like result
sharper resolution for contact discontinuity
0 0.5 1
0 2 4 6 8
ρ
0 0.5 1
−5 0 5 10 15
0 0.5 1
0 200 400 600
u p
x x
x
t = 0.038 t = 0.038
t = 0.038
Woodward-Colella’s Problem
Physical grid coordinates at selected times
Each little dashed line gives a cell-center location of the proposed Lagrange-like grid system
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0 0.01 0.02 0.03 0.04
x
time
2 D Riemann Problem
With initial 4-shock wave pattern
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
ρ u v p
=
1.5
0 0 1.5
0.532 1.206
0 0.3
0.138 1.206 1.206 0.029
0.532 0 1.206
0.3
2 D Riemann Problem
With initial 4-shock wave pattern Lagrangian-like result
Occurrence of simple Mach reflection
0.2 0.4 0.6 0.8 0.2
0.4 0.6 0.8
0.2 0.4 0.6 0.8 0.2
0.4 0.6 0.8
0.2 0.4 0.6 0.8 0.2
0.4 0.6 0.8
Density Pressure Physical grid
2 D Riemann Problem
With initial 4-shock wave pattern Eulerian result
Poor resolution around simple Mach reflection
0.2 0.4 0.6 0.8 0.2
0.4 0.6 0.8
0.2 0.4 0.6 0.8 0.2
0.4 0.6 0.8
0.2 0.4 0.6 0.8 0.2
0.4 0.6 0.8
Density Pressure Physical grid
Radially Symmetric Problem
0 0.2 0.4
0 0.1 0.2 0.3 0.4 0.5
0 0.2 0.4
0 0.1 0.2 0.3 0.4 0.5
0 0.2 0.4
0 0.1 0.2 0.3 0.4
Density Pressure 0.5 Physical grid a) h0 = 0.99
0.1 0.2 0.3 0.4 0.5
0.1 0.2 0.3 0.4 0.5
0.1 0.2 0.3 0.4
Density Pressure 0.5 Physical grid b) h0 = 0
Radially Symmetric Prob. (Cont.)
0 0.1 0.2 0.3 0.4 0.5
0.2 0.4 0.6 0.8 1 1.2
one−d h0=0 h0=0.99
0 0.1 0.2 0.3 0.4 0.5
0 0.1 0.2 0.3 0.4
0 0.1 0.2 0.3 0.4 0.5
0 0.5 1 1.5
0 0.1 0.2 0.3 0.4 0.5
0.5 1 1.5 2
r(m) r(m)
ρ(Mg/m3 ) ¯u(km/s)
p(GPa) J
