Hiro Nishikawa
Divergence Formulation of Source Term Hiro Nishikawa National Institute of Aerospace CFD Seminar, December 4, 2012
Give Up or Never Give Up “ There're times when it is good to give up.” “I agree.”
I give up to get creative.
Interesting Schemes Economical high-order schemes Residual-distribution schemes (Roe,VKI, INRIA, etc.) Residual-based compact schemes (Corre and Lerat, JCP2001) Third-order edge-based finite-volume scheme (Katz and Sankaran JCP2011) These schemes contain the target equation (or residual) in the truncation error (TE): E.g., for linear advection, an RD scheme has the following TE,
h TE = (a∂x + b∂y )(a∂x u + b∂y u) + O(h2 ) 2a Leading term vanishes in steady state, and accuracy upgraded to second-order (Residual property).
This talk will focus on the third-order FV scheme for conservation laws with a source term.
Second-Order FV Scheme k
Conservation law: ∂x f + ∂y g = 0
nrjk n!jk
Edge-based finite-volume scheme: 0=−
!
φjk Ajk
j
k∈{kj }
with the upwind flux at edge midpoint:
φjk
1 1 ˆ jk − |λ|(uR − uL ) = (FL + FR ) · n 2 2
with the left and right solution values: 1 1 uL = uj + (∇u)j · ∆ljk , uR = uk − (∇u)k · ∆ljk 2 2
Ajk = |n!jk + nrjk |
ˆ jk = (n!jk + nrjk )/Ajk n
F = (f, g) ∆ljk = (xk − xj , yk − yj ) ˆ jk λ = (∂u f, ∂u g) · n
Second-order accurate with first-order accurate gradients. NASA’s FUN3D, Software Cradle’s SC/Tetra, etc.
Third-Order FV Scheme
(Katz and Sankaran JCP2011)
1. Extrapolate the fluxes:
φjk =
1 1 ˆ jk − |λ|(uR − uL ) (FL + FR ) · n 2 2
Left and right fluxes are computed by 1 FL = Fj + (∇F)j · ∆ljk , 2
1 FR = Fk − (∇F)k · ∆ljk 2
2. Second-order gradients (e.g., LSQ quadratic fit) (∇u)j =
!
k∈{kj }
(uk − uj )
cxjk y cjk
LSQ coefficients
(∇F)j =
∂f ∂u ux
∂f ∂u uy
∂g ∂u ux
∂g ∂u uy
j
The resulting scheme has the truncation error on triangular(tetrahedral) grids: 2
3
T E = (C1 ∂xx + C2 ∂xy + C3 ∂yy )(∂x f + ∂y g)h + O(h )
Third-order scheme on second-order stencil
Conservation Law with Source For a conservation law with a source term:
∂x f + ∂y g = s
Source term includes a time-derivative term.
We add a source term discretization to the third-order scheme: " ! ! 0=− φjk Ajk + s dV k∈{kj }
Vj
s dV = sj Vj
Vj
Source term must be discretized to yield
T E = (C1 ∂xx + C2 ∂xy + C3 ∂yy )(∂x f + ∂y g − s)h2 + O(h3 ) This is critical for extending the third-order scheme to time-dependent problems.
Special formulas exist for regular grids.
Formulas for Regular Grids 5
4
Equilateral-triangular stencil: !
Vj s dV = (6sj + s1 + s2 + s3 + s4 + s5 + s6 ) 12 Vj
Galerkin Discretization
6
3
j
(Katz and Sankaran JCP2011)
2
1 5
Right-triangular stencil: !
Vj s dV = (30sj − s1 − s2 − s3 − s4 − s5 − s6 ) 24 Vj (Nishikawa, JCP2012)
How can I come up with such a formula for irregular grids?
6
1
4
j
2
3
I can’t.
New Problem Can we write the source term in the divergence form? s
s −→ ∂x f + ∂y g
s
If possible, we can write
∂x f + ∂y g = s s
∂x f + ∂y g = ∂x f + ∂y g s
s
s
∂x (f − f ) + ∂y (g − g ) = 0
Then, source term discretization will not be needed.
Nishikawa, JCP2012
Divergence Formulation of Source Term where
s −→ ∂x f s + ∂y g s
1 1 1 2 f = (x − xj )s + (x − xj ) ∂x s + (x − xj )3 ∂xx s 2 4 12 1 1 1 s 2 g = (y − yj )s + (y − yj ) ∂y s + (y − yj )3 ∂yy s 2 4 12 s
(xj, yj) is a point in a computational grid.
Then, ∂x f + ∂y g = s can be written as a single divergence form: ∂x (f − f s ) + ∂y (g − g s ) = 0
Source term discretization is no longer needed. Gradient and Hessian of the source are needed, which can be computed by the quadratic fit.
Nishikawa, JCP2012
Equivalent up to Third-Order The divergence form, s
s
∂x (f − f ) + ∂y (g − g ) = 0
can be expanded as 1 1 3 3 ∂x f + ∂y g = s + (x − xj ) ∂xxx s + (y − yj ) ∂yyy s 12 12
At node j, it is equivalent to the original equation. In the neighborhood, equivalent up to third-order, which is sufficient for third-order scheme.
Nishikawa, JCP2012
One-Component Forms The divergence form is equivalent to the following:
1 1 2 3 f = (x − xj )s + (x − xj ) ∂x s + (x − xj ) ∂xx s 2 6 gs = 0 s
or s
f =0 1 1 2 3 g = (y − yj )s + (y − yj ) ∂y s + (y − yj ) ∂yy s 2 6 s
All are equivalent to one another up to third-order.
Nishikawa, JCP2012
Third-Order Scheme Conservation law: ∂x (f − f s ) + ∂y (g − g s ) = 0 Edge-based finite-volume scheme: ! 0=− (φjk + ψjk )Ajk k∈{kj }
with the central flux for the source flux: ψjk
1 s ˆ jk = (FL + FsR ) · n 2
Fs = (f s , g s )
with the left and right flux values: FsL
=
Fsj
1 + (∇Fs )j · ∆ljk , 2
FsR
=
Fsk
1 − (∇Fs )k · ∆ljk 2
The resulting scheme has the truncation error: T E = (C1 ∂xx + C2 ∂xy + C3 ∂yy )(∂x f + ∂y g − s)h2 + O(h3 )
Third-order achieved without source term discretization.
Nishikawa, JCP2012
One-Component Case 1 1 2 f = (x − xj )s + (x − xj ) ∂x s + (x − xj )3 ∂xx s 2 6 gs = 0 s
Edge-based finite-volume scheme: ! 0=− (φjk + ψjk )Ajk k∈{kj }
0=−
with the central flux for the source flux: 1 s ˆ jk ψjk = (fL + fRs , 0) · n 2
!
φjk Ajk +
k∈{kj }
"
s dV Vj
with the left and right flux values: 1 1 s s s s s fL = fj + (∇f )j · ∆ljk , fR = fk − (∇f s )k · ∆ljk 2 2
Source discretization replaced by a scalar central scheme.
Some details not given in Nishikawa, JCP2012
Source Flux and Flux Gradients Source flux:
1 1 2 f = (x − xj )s + (x − xj ) ∂x s + (x − xj )3 ∂xx s 2 6 s
j
1. Compute the source flux at nodes: fjs
fks
= 0,
(∇f )j =
sj 0
(∇f )k =
sk +
1 6 (xk
Ignored for third-order 3
− xj ) ∂xxx sk
=
fjs
Ignored for third-order
(xk − xj )∂y sk + 12 (xk − xj )2 ∂xy sk + 16 (xk − xj )3 ∂xxy sk
3. Compute the left and right fluxes: fLs
ψkj
1 1 2 = (xk − xj )sk + (xk − xj ) ∂x sk + (xk − xj )3 ∂xx sk 2 6
2. Compute the gradient of the source flux:
k
ψjk
1 + (∇f s )j · ∆ljk , 2
fRs
=
fks
1 − (∇f s )k · ∆ljk 2
4. Compute the central flux for node j: ψjk
1 s ˆ jk = (fL + fRs , 0) · n 2
The other flux needs to be computed separately because ψjk != −ψkj .
Source flux discretization is not conservative (of course).
Nishikawa, JCP2012
Exact Divergence Form If the source term is simple enough, e.g.,
∂x f + ∂y g = cos(x − y) it can be written exactly as s
∂x f + ∂y g = ∂x f + ∂y g f s = sin(x − y),
s
gs = 0
( The choice is not unique. )
Therefore, again, source term discretization is not needed. 1. Second derivatives are not needed. 2. Gradient of the source term is needed. (It can be computed analytically or by a quadratic fit.) 3. This is not possible for time-derivative terms (only discrete values).
Nishikawa, JCP2012
Burgers Equation: Regular Grids ∂x (u2 /2) + ∂y u = cos(x − y)(sin(x − y) − 1)
Exact solution: u(x, y) = sin(x − y)
10
-4
10
-5
Slope 2 Slope 3 Divergence Form Exact Divergence Point (5/4, -1/4)
- Dirichlet boundary condition. - 6 neighbors for quadratic fit. - Time-stepping by RK2 to steady state.
L1 error
- n x n grids: n = 9, 17, 33, 65, 129, 257.
10
-6
10
-7
y
- (5/4, -1/4) indicates the special formula.
1
Second-order with the point discretization.
0 0
h
x
0.05 0.1
1
Nishikawa, JCP2012
Burgers Equation: Irregular Grids ∂x (u2 /2) + ∂y u = cos(x − y)(sin(x − y) − 1)
Exact solution: u(x, y) = sin(x − y)
10
-4
Slope 2 Slope 3 Divergence Form Exact Divergence Point (5/4, -1/4)
- n x n grids: n = 9, 17, 33, 65, 129, 257.
- Time-stepping by RK2 to steady state.
10-6 1
10
-7
10
-8
y
- 10 neighbors for quadratic fit. (to avoid ill-conditioning of LSQ matrix)
L1 error
- Dirichlet boundary condition.
10-5
- (5/4, -1/4) indicates the special formula.
Only the divergence formulation achieved third-order accuracy.
0 0
h
x
0.05 0.1
1
Conclusion Third-order finite-volume scheme made simple for source term by the divergence formulation. Future work: Relation with the formula of Katz (Katz 2012, unpublished); looks similar. Application to unsteady computation (time derivative as a source). Application to other discretization methods. Application to other types of source terms (involving the solution). Hyperbolic
which I might have never been able to even think about if I had not given up.
Give Up or Never Give Up “ There're times when it is good to give up.” “I agree.”
I give up to get creative.
Interesting Schemes Economical high-order schemes Residual-distribution schemes (Roe,VKI, INRIA, etc.) Residual-based compact schemes (Corre and Lerat, JCP2001) Third-order edge-based finite-volume scheme (Katz and Sankaran JCP2011) These schemes contain the target equation (or residual) in the truncation error (TE): E.g., for linear advection, an RD scheme has the following TE,
h TE = (a∂x + b∂y )(a∂x u + b∂y u) + O(h2 ) 2a Leading term vanishes in steady state, and accuracy upgraded to second-order (Residual property).
This talk will focus on the third-order FV scheme for conservation laws with a source term.
Second-Order FV Scheme k
Conservation law: ∂x f + ∂y g = 0
nrjk n!jk
Edge-based finite-volume scheme: 0=−
!
φjk Ajk
j
k∈{kj }
with the upwind flux at edge midpoint:
φjk
1 1 ˆ jk − |λ|(uR − uL ) = (FL + FR ) · n 2 2
with the left and right solution values: 1 1 uL = uj + (∇u)j · ∆ljk , uR = uk − (∇u)k · ∆ljk 2 2
Ajk = |n!jk + nrjk |
ˆ jk = (n!jk + nrjk )/Ajk n
F = (f, g) ∆ljk = (xk − xj , yk − yj ) ˆ jk λ = (∂u f, ∂u g) · n
Second-order accurate with first-order accurate gradients. NASA’s FUN3D, Software Cradle’s SC/Tetra, etc.
Third-Order FV Scheme
(Katz and Sankaran JCP2011)
1. Extrapolate the fluxes:
φjk =
1 1 ˆ jk − |λ|(uR − uL ) (FL + FR ) · n 2 2
Left and right fluxes are computed by 1 FL = Fj + (∇F)j · ∆ljk , 2
1 FR = Fk − (∇F)k · ∆ljk 2
2. Second-order gradients (e.g., LSQ quadratic fit) (∇u)j =
!
k∈{kj }
(uk − uj )
cxjk y cjk
LSQ coefficients
(∇F)j =
∂f ∂u ux
∂f ∂u uy
∂g ∂u ux
∂g ∂u uy
j
The resulting scheme has the truncation error on triangular(tetrahedral) grids: 2
3
T E = (C1 ∂xx + C2 ∂xy + C3 ∂yy )(∂x f + ∂y g)h + O(h )
Third-order scheme on second-order stencil
Conservation Law with Source For a conservation law with a source term:
∂x f + ∂y g = s
Source term includes a time-derivative term.
We add a source term discretization to the third-order scheme: " ! ! 0=− φjk Ajk + s dV k∈{kj }
Vj
s dV = sj Vj
Vj
Source term must be discretized to yield
T E = (C1 ∂xx + C2 ∂xy + C3 ∂yy )(∂x f + ∂y g − s)h2 + O(h3 ) This is critical for extending the third-order scheme to time-dependent problems.
Special formulas exist for regular grids.
Formulas for Regular Grids 5
4
Equilateral-triangular stencil: !
Vj s dV = (6sj + s1 + s2 + s3 + s4 + s5 + s6 ) 12 Vj
Galerkin Discretization
6
3
j
(Katz and Sankaran JCP2011)
2
1 5
Right-triangular stencil: !
Vj s dV = (30sj − s1 − s2 − s3 − s4 − s5 − s6 ) 24 Vj (Nishikawa, JCP2012)
How can I come up with such a formula for irregular grids?
6
1
4
j
2
3
I can’t.
New Problem Can we write the source term in the divergence form? s
s −→ ∂x f + ∂y g
s
If possible, we can write
∂x f + ∂y g = s s
∂x f + ∂y g = ∂x f + ∂y g s
s
s
∂x (f − f ) + ∂y (g − g ) = 0
Then, source term discretization will not be needed.
Nishikawa, JCP2012
Divergence Formulation of Source Term where
s −→ ∂x f s + ∂y g s
1 1 1 2 f = (x − xj )s + (x − xj ) ∂x s + (x − xj )3 ∂xx s 2 4 12 1 1 1 s 2 g = (y − yj )s + (y − yj ) ∂y s + (y − yj )3 ∂yy s 2 4 12 s
(xj, yj) is a point in a computational grid.
Then, ∂x f + ∂y g = s can be written as a single divergence form: ∂x (f − f s ) + ∂y (g − g s ) = 0
Source term discretization is no longer needed. Gradient and Hessian of the source are needed, which can be computed by the quadratic fit.
Nishikawa, JCP2012
Equivalent up to Third-Order The divergence form, s
s
∂x (f − f ) + ∂y (g − g ) = 0
can be expanded as 1 1 3 3 ∂x f + ∂y g = s + (x − xj ) ∂xxx s + (y − yj ) ∂yyy s 12 12
At node j, it is equivalent to the original equation. In the neighborhood, equivalent up to third-order, which is sufficient for third-order scheme.
Nishikawa, JCP2012
One-Component Forms The divergence form is equivalent to the following:
1 1 2 3 f = (x − xj )s + (x − xj ) ∂x s + (x − xj ) ∂xx s 2 6 gs = 0 s
or s
f =0 1 1 2 3 g = (y − yj )s + (y − yj ) ∂y s + (y − yj ) ∂yy s 2 6 s
All are equivalent to one another up to third-order.
Nishikawa, JCP2012
Third-Order Scheme Conservation law: ∂x (f − f s ) + ∂y (g − g s ) = 0 Edge-based finite-volume scheme: ! 0=− (φjk + ψjk )Ajk k∈{kj }
with the central flux for the source flux: ψjk
1 s ˆ jk = (FL + FsR ) · n 2
Fs = (f s , g s )
with the left and right flux values: FsL
=
Fsj
1 + (∇Fs )j · ∆ljk , 2
FsR
=
Fsk
1 − (∇Fs )k · ∆ljk 2
The resulting scheme has the truncation error: T E = (C1 ∂xx + C2 ∂xy + C3 ∂yy )(∂x f + ∂y g − s)h2 + O(h3 )
Third-order achieved without source term discretization.
Nishikawa, JCP2012
One-Component Case 1 1 2 f = (x − xj )s + (x − xj ) ∂x s + (x − xj )3 ∂xx s 2 6 gs = 0 s
Edge-based finite-volume scheme: ! 0=− (φjk + ψjk )Ajk k∈{kj }
0=−
with the central flux for the source flux: 1 s ˆ jk ψjk = (fL + fRs , 0) · n 2
!
φjk Ajk +
k∈{kj }
"
s dV Vj
with the left and right flux values: 1 1 s s s s s fL = fj + (∇f )j · ∆ljk , fR = fk − (∇f s )k · ∆ljk 2 2
Source discretization replaced by a scalar central scheme.
Some details not given in Nishikawa, JCP2012
Source Flux and Flux Gradients Source flux:
1 1 2 f = (x − xj )s + (x − xj ) ∂x s + (x − xj )3 ∂xx s 2 6 s
j
1. Compute the source flux at nodes: fjs
fks
= 0,
(∇f )j =
sj 0
(∇f )k =
sk +
1 6 (xk
Ignored for third-order 3
− xj ) ∂xxx sk
=
fjs
Ignored for third-order
(xk − xj )∂y sk + 12 (xk − xj )2 ∂xy sk + 16 (xk − xj )3 ∂xxy sk
3. Compute the left and right fluxes: fLs
ψkj
1 1 2 = (xk − xj )sk + (xk − xj ) ∂x sk + (xk − xj )3 ∂xx sk 2 6
2. Compute the gradient of the source flux:
k
ψjk
1 + (∇f s )j · ∆ljk , 2
fRs
=
fks
1 − (∇f s )k · ∆ljk 2
4. Compute the central flux for node j: ψjk
1 s ˆ jk = (fL + fRs , 0) · n 2
The other flux needs to be computed separately because ψjk != −ψkj .
Source flux discretization is not conservative (of course).
Nishikawa, JCP2012
Exact Divergence Form If the source term is simple enough, e.g.,
∂x f + ∂y g = cos(x − y) it can be written exactly as s
∂x f + ∂y g = ∂x f + ∂y g f s = sin(x − y),
s
gs = 0
( The choice is not unique. )
Therefore, again, source term discretization is not needed. 1. Second derivatives are not needed. 2. Gradient of the source term is needed. (It can be computed analytically or by a quadratic fit.) 3. This is not possible for time-derivative terms (only discrete values).
Nishikawa, JCP2012
Burgers Equation: Regular Grids ∂x (u2 /2) + ∂y u = cos(x − y)(sin(x − y) − 1)
Exact solution: u(x, y) = sin(x − y)
10
-4
10
-5
Slope 2 Slope 3 Divergence Form Exact Divergence Point (5/4, -1/4)
- Dirichlet boundary condition. - 6 neighbors for quadratic fit. - Time-stepping by RK2 to steady state.
L1 error
- n x n grids: n = 9, 17, 33, 65, 129, 257.
10
-6
10
-7
y
- (5/4, -1/4) indicates the special formula.
1
Second-order with the point discretization.
0 0
h
x
0.05 0.1
1
Nishikawa, JCP2012
Burgers Equation: Irregular Grids ∂x (u2 /2) + ∂y u = cos(x − y)(sin(x − y) − 1)
Exact solution: u(x, y) = sin(x − y)
10
-4
Slope 2 Slope 3 Divergence Form Exact Divergence Point (5/4, -1/4)
- n x n grids: n = 9, 17, 33, 65, 129, 257.
- Time-stepping by RK2 to steady state.
10-6 1
10
-7
10
-8
y
- 10 neighbors for quadratic fit. (to avoid ill-conditioning of LSQ matrix)
L1 error
- Dirichlet boundary condition.
10-5
- (5/4, -1/4) indicates the special formula.
Only the divergence formulation achieved third-order accuracy.
0 0
h
x
0.05 0.1
1
Conclusion Third-order finite-volume scheme made simple for source term by the divergence formulation. Future work: Relation with the formula of Katz (Katz 2012, unpublished); looks similar. Application to unsteady computation (time derivative as a source). Application to other discretization methods. Application to other types of source terms (involving the solution). Hyperbolic
which I might have never been able to even think about if I had not given up.
