Energy Estimates and Weak Boundary Procedures for LAM Marco - - PowerPoint PPT Presentation

Energy Estimates and Weak Boundary Procedures for LAM

Marco Kupiainen1,2

1Department of Computational Mathematics

University of Link¨

• ping

2SMHI, Rossby Centre

marco.kupiainen@smhi.se

ASM 2013 Reykjavik, Iceland

Outline

Introduction/Motivation for this work Motivating examples Model problem Why do we impose Boundary Conditions? Short Description of Boundary methods Energy estimate for Davies-relaxation Spectral Radius of operators Computational results Conclusion Conclusion II

Introduction/Motivation for this work

◮ Boundary conditions are set using ad hoc intuition (Davies

relaxation)

◮ Different effects in e.g. Arctic and Europe ◮ ”The apparent success of spectral nudging for one-way

nesting is at least partly an artifact of very bad procedures for windowing and blending”1

1John P. Boyd ”Limited-Area Fourier Spectral Models and Data Analysis

Schemes: Windows, Fourier Extension, Davies Relaxation and All That”, Mon.

• Weat. Rev. Vol. 133 2005

Motivating example

◮ Europe has ”standing waves” ◮ Similar problems over Arctic (mitigated with Spectral

Nudging)

/home/sm_marku/test/Cordex/fc198909010000qq

2 4 6 8 10 12 14 16 18 20

◮ Analyze boundary procedures!

Model problem

Simplification for purpose of illustration, without loosing generality!

◮

ut + 3

i=1 ∇

Fi( u) = 0 N-S/Euler/primitive equations

Model problem

Simplification for purpose of illustration, without loosing generality!

◮

ut + 3

i=1 ∇

Fi( u) = 0 N-S/Euler/primitive equations

◮

⇓ 1D system

Model problem

Simplification for purpose of illustration, without loosing generality!

◮

ut + 3

i=1 ∇

Fi( u) = 0 N-S/Euler/primitive equations

◮

⇓ 1D system

◮

ut + f ( u)x = 0

Model problem

Simplification for purpose of illustration, without loosing generality!

◮

ut + 3

i=1 ∇

Fi( u) = 0 N-S/Euler/primitive equations

◮

⇓ 1D system

◮

ut + f ( u)x = 0

◮

⇓ Semi-linear system

Model problem

Simplification for purpose of illustration, without loosing generality!

◮

ut + 3

i=1 ∇

Fi( u) = 0 N-S/Euler/primitive equations

◮

⇓ 1D system

◮

ut + f ( u)x = 0

◮

⇓ Semi-linear system

◮

ut + A( u) ux = 0

Model problem

Simplification for purpose of illustration, without loosing generality!

◮

ut + 3

i=1 ∇

Fi( u) = 0 N-S/Euler/primitive equations

◮

⇓ 1D system

◮

ut + f ( u)x = 0

◮

⇓ Semi-linear system

◮

ut + A( u) ux = 0

◮

⇓ Linear system

Model problem

Simplification for purpose of illustration, without loosing generality!

◮

ut + 3

i=1 ∇

Fi( u) = 0 N-S/Euler/primitive equations

◮

⇓ 1D system

◮

ut + f ( u)x = 0

◮

⇓ Semi-linear system

◮

ut + A( u) ux = 0

◮

⇓ Linear system

◮

ut + A ux = 0, e.g.A =

• 1

1

Model problem

Simplification for purpose of illustration, without loosing generality!

◮

ut + 3

i=1 ∇

Fi( u) = 0 N-S/Euler/primitive equations

◮

⇓ 1D system

◮

ut + f ( u)x = 0

◮

⇓ Semi-linear system

◮

ut + A( u) ux = 0

◮

⇓ Linear system

◮

ut + A ux = 0, e.g.A =

• 1

1

• ◮

⇓ A = TΛT −1 ⇒ ut + TΛT −1 ux = 0 Scalar PDE

Model problem

Simplification for purpose of illustration, without loosing generality!

◮

ut + 3

i=1 ∇

Fi( u) = 0 N-S/Euler/primitive equations

◮

⇓ 1D system

◮

ut + f ( u)x = 0

◮

⇓ Semi-linear system

◮

ut + A( u) ux = 0

◮

⇓ Linear system

◮

ut + A ux = 0, e.g.A =

• 1

1

• ◮

⇓ A = TΛT −1 ⇒ ut + TΛT −1 ux = 0 Scalar PDE

◮ ut + ux = 0

Why do we impose Boundary Conditions?

(Energy point of view)

ut + ux = 0, x ∈ [0, 1]

Why do we impose Boundary Conditions?

(Energy point of view)

ut + ux = 0, x ∈ [0, 1]

◮ Multiply with u and integrate!

Why do we impose Boundary Conditions?

(Energy point of view)

ut + ux = 0, x ∈ [0, 1]

◮ Multiply with u and integrate! ◮ 2

1

0 uutdx + 2

1

0 uux = 0

Why do we impose Boundary Conditions?

(Energy point of view)

ut + ux = 0, x ∈ [0, 1]

◮ Multiply with u and integrate! ◮ 2

1

0 uutdx + 2

1

0 uux = 0 ◮ Use integration by parts!

Why do we impose Boundary Conditions?

(Energy point of view)

ut + ux = 0, x ∈ [0, 1]

◮ Multiply with u and integrate! ◮ 2

1

0 uutdx + 2

1

0 uux = 0 ◮ Use integration by parts! ◮ d dt u2 + [u2]1 0 = 0,

u2 = 1

0 u2dx

Why do we impose Boundary Conditions?

(Energy point of view)

ut + ux = 0, x ∈ [0, 1]

◮ Multiply with u and integrate! ◮ 2

1

0 uutdx + 2

1

0 uux = 0 ◮ Use integration by parts! ◮ d dt u2 + [u2]1 0 = 0,

u2 = 1

0 u2dx ◮ d dt u2 =

(u(0, t))2

• Faster than exponential growth

− (u(1, t))2

• Decay

Why do we impose Boundary Conditions?

(Energy point of view)

ut + ux = 0, x ∈ [0, 1]

◮ Multiply with u and integrate! ◮ 2

1

0 uutdx + 2

1

0 uux = 0 ◮ Use integration by parts! ◮ d dt u2 + [u2]1 0 = 0,

u2 = 1

0 u2dx ◮ d dt u2 =

(u(0, t))2

• Faster than exponential growth

− (u(1, t))2

• Decay

◮ We must set conditions to bound the energy for u(0, t) with

data: g(t) < C for some constant C ∈ R

Discrete Boundary procedures

CKD : Ut + P−1QU = 0 U(·, t) = (I − W )U(·, t) + WG(t) W = diag(tanh)

◮ G is data! ◮ U = (U0, U1, . . . UN)T is the solution. ◮

CKD = Classic K˚ allberg-Davies Relaxation

Discrete Boundary procedures

WKD : Ut + P−1QU = P−1W (G − U) W = diag(tanh)

◮

WKD = Weak K˚ allberg-Davies Relaxation, proven stable

Discrete Boundary procedures

SAT : Ut + P−1QU = P−1E0(G − U) E0 = diag(1, 0, . . . , 0)

◮

SAT = Simultaneous Approximation Term (Carpenter et. al.), proven stable

Discrete Boundary procedures

CKD : Ut + P−1QU = 0 U(·, t) = (I − W )U(·, t) + WG(t) W = diag(tanh) WKD : Ut + P−1QU = P−1W (G − U) W = diag(tanh) SAT : Ut + P−1QU = P−1E0(G − U) E0 = diag(1, 0, . . . , 0)

◮

CKD = Classic K˚ allberg-Davies Relaxation

◮

WKD = Weak K˚ allberg-Davies Relaxation, proven stable

◮

SAT = Simultaneous Approximation Term (Carpenter et. al.), proven stable

Energy estimate for Davies-relaxation

Strong Davies relaxation Ut + P−1QU = 0 U = U + Wtanh(G − U)

◮ No energy estimate! (For

stability(?) proof use GKS theory DIFFICULT!) Weak Davies relaxation Ut + P−1QU = τP−1Wtanh(G − U)

◮ Discrete Energy Method to

prove stability

◮ Multiply with UTP ◮ Use summation-by-parts

property of difference operator

◮ In fact there are counterexamples shoving instability for strong

methods.2

2”High Order Difference Methods for Time Dependent PDE”, Bertil

Gustafsson ISBN 978-3-540-74992-9, Springer Verlag 2008

Energy estimate for Davies-relaxation cont.

d dt U2

P = U2 1(1 − 2w1) + 2w1U1G1 − U2 N(1 + 2wN) + 2wNUNGN

−

N−1

• i=2

(2wiU2

i − 2wiUiGi)

(1) Since w1 ≥ 1

2 and wi ≥ 0 ∴ Davies Relaxation is proven

energy-stable!

Energy estimate for Davies-relaxation cont.

d dt U2

P = U2 1(1 − 2w1) + 2w1U1G1 − U2 N(1 + 2wN) + 2wNUNGN

−

N−1

• i=2

(2wiU2

i − 2wiUiGi)

(1) Since w1 ≥ 1

2 and wi ≥ 0 ∴ Davies Relaxation is proven

energy-stable!

◮ Stable, BUT imposing conditions where not needed!

Energy estimate for Davies-relaxation cont.

d dt U2

P = U2 1(1 − 2w1) + 2w1U1G1 − U2 N(1 + 2wN) + 2wNUNGN

−

N−1

• i=2

(2wiU2

i − 2wiUiGi)

(1) Since w1 ≥ 1

2 and wi ≥ 0 ∴ Davies Relaxation is proven

energy-stable!

◮ Stable, BUT imposing conditions where not needed! ◮ Solution quality is not affected if Ui − Gi is ”small”.

Energy estimate for Davies-relaxation cont.

d dt U2

P = U2 1(1 − 2w1) + 2w1U1G1 − U2 N(1 + 2wN) + 2wNUNGN

−

N−1

• i=2

(2wiU2

i − 2wiUiGi)

(1) Since w1 ≥ 1

2 and wi ≥ 0 ∴ Davies Relaxation is proven

energy-stable!

◮ Stable, BUT imposing conditions where not needed! ◮ Solution quality is not affected if Ui − Gi is ”small”. ◮ Usually we impose time-interpolated G i.e. Gi(t) = π6hGi(tn)

Spectral Radius of operators part I

−2 −1.5 −1 −0.5 0.5 1 1.5 2 x 10

−5

−150 −100 −50 50 100 150 Imaginary Real Spectral radius of operator WITHOUT Boundary conditions

(a) Operator without bound- ary conditions

−140 −120 −100 −80 −60 −40 −20 −50 −40 −30 −20 −10 10 20 30 40 50 Spectral radius of WEAK DAVIES RELAXATION Real Imaginary

(b) WKD

−6 −5 −4 −3 −2 −1 −50 −40 −30 −20 −10 10 20 30 40 50 Spectral radius of SAT Real Imaginary

(c) SAT

Spectral Radius of operators part II

Increasing resolution by a factor of 5:

−350 −300 −250 −200 −150 −100 −50 −150 −100 −50 50 100 150 Real Imaginary Spectral radius of WEAK DAVIES RELAXATION

(d) WKD

−8 −7 −6 −5 −4 −3 −2 −1 −150 −100 −50 50 100 150 Real Imaginary Spectral radius of SAT

(e) SAT

ρ(WKD) ≈ 350, ρ(SAT) ≈ 140

Computational results

ut + ux = 0 u(0, t) = G(0, t), x ∈ [0, 1] Exact solution is u(x, t) = sin(2π(x − t − 1 2)) = G(x, t) We use the exact solution as initial data, i.e. assume perfect assimilation

Computational results

ut + ux = 0 u(0, t) = G(0, t), x ∈ [0, 1] Exact solution is u(x, t) = sin(2π(x − t − 1 2)) = G(x, t) We use the exact solution as initial data, i.e. assume perfect assimilation

◮ Use exact G(·, t) as boundary data (show movie)

Computational results

ut + ux = 0 u(0, t) = G(0, t), x ∈ [0, 1] Exact solution is u(x, t) = sin(2π(x − t − 1 2)) = G(x, t) We use the exact solution as initial data, i.e. assume perfect assimilation

◮ πG is linear interpolation in time (show movie)

Computational results

ut + ux = 0 u(0, t) = G(0, t), x ∈ [0, 1] Exact solution is u(x, t) = sin(2π(x − t − 1 2)) = G(x, t) We use the exact solution as initial data, i.e. assume perfect assimilation

◮ πG is constant interpolation in time (closest in time) (show

movie)

Computational results

ut + ux = 0 u(0, t) = G(0, t), x ∈ [0, 1] Exact solution is u(x, t) = sin(2π(x − t − 1 2)) = G(x, t) We use the exact solution as initial data, i.e. assume perfect assimilation

◮ πG is P3-Hermite (no new minima or maxima are introduced)

Computational results

ut + ux = 0 u(0, t) = G(0, t), x ∈ [0, 1] Exact solution is u(x, t) = sin(2π(x − t − 1 2)) = G(x, t) We use the exact solution as initial data, i.e. assume perfect assimilation

◮ πG is P3-Spline (new minima, maxima can be introduced)

Computational results

ut + ux = 0 u(0, t) = G(0, t), x ∈ [0, 1] Exact solution is u(x, t) = sin(2π(x − t − 1 2)) = G(x, t) We use the exact solution as initial data, i.e. assume perfect assimilation

◮ Mismatching data on outflow boundary

Conclusion

◮ Davies relaxation in its weak form is a penalty method ◮ WKD is proven stable and yields similar results with CKD ◮ CKD and WKD impose data on all boundaries, even when not

needed.

◮ SAT is proven stable, imposes data only where needed! ◮ The WKD operator is much more stiff than SAT by a factor

• f ≈20!

◮ When data is not close to the solution on outflow boundaries

(and WKD) yields unsatisfactory results ⇒ horizontal diffusion mitigates this problem

◮ The excessive diffusion can be the reason for the poor results

• ver the Arctic

◮ Non-matching data causes a ”standing wave” on the outflow

boundary with WKD and CKD

◮ SAT yields similar results for exact and almost matching data

(time interpolation error is visible), but outperforms CKD and WKD for non-matching data.

Conclusion II

◮ Theory for penalty-based boundary conditions is considered

mature and ready to be used for NWP and climate

◮ Results are already extended to non-linear multidimensional

systems, but for purpose of illustration a model problem was shown here.

Conclusion II

◮ THANK YOU FOR LISTENING!