[PPT] - So Solving g the S N Ne Neutron on Transpor ort Equation on Us PowerPoint Presentation

SLIDE 1

So Solving g the SN Ne Neutron

n Transpor
rt Equation
n

Us Using H High O Order L Lax-Fr Friedrich chs W WENO Fa Fast Sweeping Methods

Dean Wang, Tseelmaa Byambaakhuu The Ohio State University Sebastian Schunert Idaho National Laboratory Zeyun Wu Virginia Commonwealth University M&C 2019, Portland, Oregon, USA August 25-29, 2019

SLIDE 2

Outline

Background and motivation
Robustness
High order
Efficiency
LF-WENO methods
Theory (Wang 2019)
Numerical properties (Wang 2019, NSE)
Diffusion limit
Conclusion

2

SLIDE 3

Numerical methods for SN

Finite difference sweeping methods
SD: 1st-order upwind; positivity preserving
DD: 2nd-order; not positivity preserving
SC: weighted DD; 2nd-order; positivity preserving; less accurate than

DD for diffusive problems

Short characteristic methods
SC: constant source
LC: linear source & linear incoming flux; positivity preserving?
QC : Quadratic source & quadratic incoming flux; can be made to be

positivity preserving (on-going work)

Galerkin methods: LD, FEM, DFEM
High-order
FEM or DFEM can be very robust with stabilization; however may

not as efficient as finite difference sweeping methods.

3

SLIDE 4

Motivation

A sweeping based numerical method is more

accurate than DD, and much more robust as well.

A challenging task…
Chen et al. in 2013 proposed Lax-Friedrichs fast

sweeping methods for steady-state hyperbolic conservation laws.

A perfect framework for the SN transport equation!

4

SLIDE 5

SN in 2-D conservative form

𝑔 𝜔 # + 𝑕 𝜔 & + Σ(𝜔 = 𝑡 𝜔, 𝑦, 𝑧

𝑔 𝜔 = 𝜈𝜔 , 𝑕 𝜔 = 𝜃𝜔 , 𝑡 𝜔, 𝑦, 𝑧 = 01

2 𝜚 𝑦, 𝑧 + 4 2 𝑅 𝑦, 𝑧

where,

5

SLIDE 6

Finite difference discretization

𝑗, 𝑘

𝑗 − 1 2 , 𝑘 𝑗 + 1 2 , 𝑘 𝑗, 𝑘 + 1 2 𝑗, 𝑘 − 1 2

; 𝑔<=4

>,? − ;

𝑔<@4

>,?

∆𝑦 + B 𝑕<,?=4

>

− B 𝑕<,?@4

>

∆𝑧 + Σ(𝜔<,? = 𝑡 𝜔<,?, 𝑦<, 𝑧? ; 𝑔<±D

E,? and B

𝑕<,?±D

E are numerical fluxes

Where (1)

6

SLIDE 7

High order WENO fluxes

For 2𝐿 − 1 – th order WENO scheme, the 𝐿 numerical fluxes are computed as ; 𝑔

<=D

E,?

L

= ∑NOP

Q@4 𝑑LN𝑔 <@L=N,? ,

𝑠 = 0, … , 𝐿 − 1 , which corresponds to 𝐿 different stencils: 𝑇L 𝑗 = 𝑦<@L, 𝑧? , … , 𝑦<@L=Q@4, 𝑧? , 𝑠 = 0, … , 𝐿 − 1. Each of these numerical fluxes is 𝑙– th order accurate. The 2𝐿 − 1 – th order WENO flux is a superposition of all these K numerical fluxes ; 𝑔<=4

>,? = X NOP Q@4

𝑥N ; 𝑔

<=4 >,? N

The nonlinear weights 𝑥N satisfy 𝑥N ≥ 0, ∑NOP

Q@4 𝑥N = 1, and are defined as

𝑥N =

[\ ∑\]^

_`D [\,

𝛽N =

b\ c=d\ .

(2) (3) (4)

7

SLIDE 8

Third order WENO (WENO3)

And the linear weights are given by Smoothness indicators are given by

where 𝜐P = 𝑏 ∗ max 𝑏𝑐𝑡 Σ(<=4,? − Σ(<,? , 𝑏𝑐𝑡 Σn<=4,? − Σn<,? ∆𝑦 𝜐4 = 𝑐 ∗ max 𝑏𝑐𝑡 Σ(<,? − Σ(<@4,? , 𝑏𝑐𝑡 Σn<,? − Σn<@4,? ∆𝑦

𝑒P = >

p ,

𝑒4 = 4

p

𝛾P = 𝜐P 𝑔

<=4,? − 𝑔 <,? > ,

𝛾4 = 𝜐4 𝑔

<,? − 𝑔 <@4,? >

; 𝑔

<=D

E,?

P

= 4

> 𝑔 <,? + 4 > 𝑔 <=4,? ,

; 𝑔

<=D

E,?

4

= − 4

> 𝑔 <@4,? + p > 𝑔 <,?

For 𝐿 = 2 , the 2nd-order accurate numerical fluxes for 𝜈 > 0 are given as

(5) (6) (7)

8

SLIDE 9

Lax-Friedrichs sweeping framework

; ; 𝑔<=D

E,? = ;

𝑔<=D

E,? +

vw >

𝜔<=4,? − 𝜔<,? , 𝑗 = 1, … , 𝑂# y B 𝑕<,?=D

E = B

𝑕<,?=D

E +

vz >

𝜔<,?=4 − 𝜔<,? , 𝑘 = 1, … , 𝑂&

Define Lax–Friedrichs fluxes:

; 𝑔<=D

E,? = ;

; 𝑔<=D

E,? −

vw >

𝜔<=4,? − 𝜔<,? B 𝑕<,?=D

E = y

B 𝑕<,?=D

E −

vz >

𝜔<,?=4 − 𝜔<,?

Then we have

(8𝑏) (8𝑐)

9

SLIDE 10

LF-WENO

; ; 𝑔<=4

>,? − 𝜏𝜈

2 𝜔<=4,? − 𝜔<,? − ; ; 𝑔<@4

>,? + 𝜏𝜈

2 𝜔<,? − 𝜔<@4,? ∆𝑦 + y B 𝑕<,?=4

>

− 𝜏𝜃 2 𝜔<,?=4 − 𝜔<,? − y B 𝑕<,?@4

>

+ 𝜏𝜃 2 𝜔<,? − 𝜔<,?@4 ∆𝑧 + Σ(𝜔<,? = 𝑡 𝜔<,?, 𝑦<, 𝑧? 𝜔<,? =

n }~,•,#~,&• ∆#@ ; ; €~•D

E,•@ ;

; €~`D

E,•@‚ƒ E }~•D,•=}~`D,•

@ y B „~,••D

E

@ y B „~,•`D

E

@‚…

E }~,••D=}~,•`D ∆† ∆‡

v w=z ∆†

∆‡

=ˆ‰∆#

(9)

10

SLIDE 11

Computing algorithm

Initialize 𝜔<,? and 𝑇<,?
While 𝑓 > etol
1. for 𝑜 = 1:

Ž 2

% sweeping in angle (𝜈 > 0, 𝜃 > 0) for 𝑗 = 1: 𝑂𝑦 % sweeping in x for 𝑘 = 1: 𝑂𝑧 % sweeping in y

Calculate ;

𝑔

<±D

E,?

N

and B 𝑕

<±D

E,?

N

, 𝑙 = 1,2 % Eq (5)

Calculate 𝛾P , 𝛾4

% Eq (7)

Calculate 𝛽N , 𝑥N, 𝑙 = 1,2

% Eq (4)

Calculate ;

𝑔<±D

E,? and ;

𝑔<±D

E,?

% Eq (5)

Calculate ;

; 𝑔<±D

E,? and y

B 𝑕<±D

E,?

% Eq (8)

Calculate 𝜔<,?

% Eq (9)

Calculate 𝑇<,?
2. for 𝑜 =

Ž 2 + 1: Ž >

% sweeping in angle (𝜈 < 0, 𝜃 > 0)

…

3. for 𝑜 =

Ž > + 1: pŽ 2

% sweeping in angle (𝜈 < 0, 𝜃 < 0) …

4. for 𝑜 =

pŽ 2 + 1: 𝑂

% sweeping in angle (𝜈 > 0, 𝜃 < 0)

11

SLIDE 12

Spatial convergence

12

SLIDE 13

Manufactured solution

𝜔 𝑦, 𝑧, 𝜈N, 𝜃N = 𝑦p𝑧p 2 − 𝑦 p 2 − 𝑧 p

𝑅N 𝑦, 𝑧 = 4 24𝑦> − 48𝑦p + 30𝑦2 − 6𝑦’ 𝑧p 2 − 𝑧 p𝜈N +𝑦p 2 − 𝑦 p 24𝑧> − 48𝑧p + 30𝑧2 − 6𝑧’ 𝜃N − Σ“𝜚 𝑦, 𝑧

13

SLIDE 14

Sweeping convergence rate

Σ( = 1 cm@4 and c =0.6 Σ( = 5 cm@4 and c =0.6

14

SLIDE 15

Computational complexity

* Computational complexity: the number of grid points x the number of iterations

15

SLIDE 16

Positivity?

Note that LF-WENO3 can be rendered to be positivity preserving using the linear scaling limiter proposed by Zhang and Shu (2010).

16

SLIDE 17

Diffusion limit of SN

𝜈• 𝑒 𝑒𝑦 𝜔• + Σ(𝜔• = Σn 2 𝜚 + 𝑅 2 Scaling Σ( →

0‰ — ,

Σ“ → 𝜁Σ“ , 𝑅 → 𝜁𝑅 , We have 𝜔• =

™ > + 𝑃 𝜁 , for 𝜁 → 0

− 𝑒 𝑒𝑦 1 3Σ› 𝑒 𝑒𝑦 𝜚 + Σœ𝜚 = 𝑅 Where 𝜚 satisfies the following diffusion equation

17

SLIDE 18

Diffusion limit – smooth solution

𝑀 = 1, ℎ = 0.1 Σ› =

4 ,

Σ¡ =

4 − 0.8ε,

𝑅 = ε

18

SLIDE 19

Diffusion limit – nonsmooth solution with boundary layer

𝜁 = 0.01

DD LF-WENO3

19

SLIDE 20

Diffusion limit – 2D

L×L = 2×2, h# = h& = 0.2 Σ› =

4 ,

Σ¡ =

4 − 0.8ε,

𝑅 = ε

20

SLIDE 21

A theoretical result on diffusion limit

(Wang 2019, NSE): Δ𝑦 = 𝜁¦ℎ = 𝜁 ⁄

4 Nℎ

Δ𝑦 = ℎ Δ𝑦 = 𝜁 ⁄

4 Nℎ

𝑚 = 0: Thick diffusion limit 𝑚 = 1: Intermidiate diffusion limit Larsen et al. 1987:

21

SLIDE 22

Conclusions

LF-WENO3 is a sweeping scheme based on the Lax–

Friedrichs fluxes with the WENO reconstruction.

It can achieve better accuracy than DD, and more

importantly it possesses good positivity-preserving property.

In addition, LF-WENO3 can achieve almost linear

computational complexity with underrelaxation.

Finally, LF-WENO3 has the diffusion limit of 𝑚 =

1/3, which lies between the thick diffusion regime (𝑚 = 0) and the intermediate regime (𝑚 = 1).

22

SLIDE 23

References

W. Chen, C.-S. Chou, and C.-Y. Kao, “Lax-Friedrichs Fast Sweeping

Methods for Steady State Problems for Hyperbolic Conservation Laws,”

J. Comput. Phys., 234, 452 (2013)
D. Wang, T. Byambaakhuu, "High Order Lax-Friedrichs WENO Fast

Sweeping Methods for the SN Neutron Transport Equation," Nucl. Sci. Eng., 193, 9, 982 (2019).

X. Zhang, C.-W. Shu, “On maximum-principle-satisfying high order

schemes for scalar conservation laws,” J. Comput. Phys., 229 (2010).

X. Zhang, C.-W. Shu, “On positivity-preserving high order discontinuous

Galerkin schemes for compressible Euler equations on rectangular meshes,” J. Comput. Phys., 229 (2010).

D. Wang, "The Asymptotic Diffusion Limit of Numerical Schemes for the

SN Transport Equation," Nucl. Sci. Eng., (2019).

E. W. Larsen, J. E. Morel, and W. F. Miller Jr., “Asymptotic Solutions of

Numerical Transport Problems in Optically Thick, Diffusive Regimes,” J.

Comput. Phys., 69, 283 (1987).

23

SLIDE 24

Thank You!

24