A Meshless Approach to Spectral Wave Modeling Adrean Webb The - - PowerPoint PPT Presentation

SLIDE 1

A Meshless Approach to Spectral Wave Modeling

Adrean Webb

The University of Tokyo Graduate School of Frontier Sciences

August 7, 2017

In Collaboration with:

Natasha Flyer (NCAR), Baylor Fox-Kemper (Brown University)

SLIDE 2

Outline

• 1. Scientific background
• 2. Coupled wave component
• 3. Overview of meshless approach
• 4. Kinematic prototype
• 5. Discussion

2 Adrean Webb / 39

SLIDE 3

Scales of interest

arbitrary energy scale tides seiches

• wind-generated waves

trans-tidal waves surges

• - ---- ----"'-------

__

ts_t_m_a_m_i_s ___ _ }nfra-gravity waves

capillary swell wind sea waves

10-6 10-5

J0-4 10-3

IQ-2

10-1

10° frequency (Hz) 10+1 24 h

3h

15min 100 s

lOs

1 s

period

0.1 s

Figure 1.1 Frequencies and periods of the vertical motions of the ocean surface (after Munk, 1950).

Figure: Illustration of wave spectra from different types of ocean surface waves (Holthuijsen, 2007).

3 Adrean Webb / 39

SLIDE 4

Spectral wave model approach

The statistical properties on intermediate scales are assumed to vary smoothly with location and the generation, propagation, and evolution of wave spectra are modeled deterministically.

Figure: The random sea of each gridded region in (a) is Fourier decomposed in (b). The statistical differences between neighboring gridded regions are assumed to be small enough such that advection and evolution of wave energy can be modeled by a PDE (Holthuijsen, 2007)

4 Adrean Webb / 39

SLIDE 5

Sample spectral model output

5 Adrean Webb / 39

Figure: Global surface Stokes drift on a 1 deg lat x 1.25 deg lon grid.

SLIDE 6

Wave action balance equation

The wave action balance equation can be derived on a coupled slab (~ x, ~ k 2 R2) from the governing linear wave equations and is given by @tW + r

~ k Ω · r ~ x W r ~ x Ω · r ~ k W = Sources

where

~ x, ~ k

• =

n g|~ k| tanh h |~ k|H(~ x) i

• 1/2

is the intermediate-water dispersion relation

• Ω

~ x, ~ k

• = ~

U ~ x

• · ~

k + ~ x, ~ k

• is the Doppler-shifted dispersion relation
• W

~ x, ~ k, t

• = g S~

k

~ k; ~ x, t

• /

~ x, ~ k

• is the spectral adiabatic invariant
• “Sources” encompass the non-kinematic physics of the waves

(generation, dissipation, nonlinear interactions, etc.)

6 Adrean Webb / 39

SLIDE 7

Wave action balance equation

The wave action balance equation can be derived on a coupled slab (~ x, ~ k 2 R2) from the governing linear wave equations and is given by @tW + r

~ k Ω · r ~ x W r ~ x Ω · r ~ k W = Sources

where

~ x, ~ k

• =

n g|~ k| tanh h |~ k|H(~ x) i

• 1/2

is the intermediate-water dispersion relation

• Ω

~ x, ~ k

• = ~

U ~ x

• · ~

k + ~ x, ~ k

• is the Doppler-shifted dispersion relation
• W

~ x, ~ k, t

• = g S~

k

~ k; ~ x, t

• /

~ x, ~ k

• is the spectral adiabatic invariant
• “Sources” encompass the non-kinematic physics of the waves

(generation, dissipation, nonlinear interactions, etc.)

7 Adrean Webb / 39

Spatial changes Spatial transport

SLIDE 8

Wave action balance equation

The wave action balance equation can be derived on a coupled slab (~ x, ~ k 2 R2) from the governing linear wave equations and is given by @tW + r

~ k Ω · r ~ x W r ~ x Ω · r ~ k W = Sources

where

~ x, ~ k

• =

n g|~ k| tanh h |~ k|H(~ x) i

• 1/2

is the intermediate-water dispersion relation

• Ω

~ x, ~ k

• = ~

U ~ x

• · ~

k + ~ x, ~ k

• is the Doppler-shifted dispersion relation
• W

~ x, ~ k, t

• = g S~

k

~ k; ~ x, t

• /

~ x, ~ k

• is the spectral adiabatic invariant
• “Sources” encompass the non-kinematic physics of the waves

(generation, dissipation, nonlinear interactions, etc.)

8 Adrean Webb / 39

Spectral changes Spectral transport

SLIDE 9

Outline

• 1. Scientific background
• 2. Coupled wave component
• 3. Overview of meshless approach
• 4. Kinematic prototype
• 5. Discussion

9 Adrean Webb / 39

SLIDE 10

Scientific Motivation

10 Adrean Webb / 39

Figure: Satellite observations of Langmuir

turbulence after the Deepwater Horizon oil spill (Quickbird)

There is a persistent, shallow mixed layer bias in the Southern Ocean in global climate models (GCMs): Langmuir turbulence missing??

Figure: Mixed layer bias is reduced in NCAR CESM model runs

(Li, Fox-Kemper, Breivik, & Webb, 2017).

Other studies: Belcher et al. 2012; D’Asaro et al., 2014; Fan & Griffies, 2014

SLIDE 11

Drawbacks for climate use:

• Computational expensive
• 4D problem in time with 6-50 x 10^6 unknowns

and nonlinear source terms

• Spatial and spectral singularities at the poles
• Difficult to model polar-ice-free climate scenarios

Third-generation model:

• Creates and evolves wave spectra in a

coupled spatial-spectral domain

• Uses structured grids (lat-lon, polar)
• Includes extensive physics and

parameterizations

11 Adrean Webb / 39

3.2 lat x 4 lon (deg) 25 freq x 24 dir

Figure: Sample spatial (SWH) and spectral output for the coupled NCAR CESM wave component.

NCAR CESM now includes a prognostic wave field

Modified NOAA WAVEWATCH III has been coupled to NCAR CESM (Craig, Li, Vertenstein, Webb)

SLIDE 12

Simple WAVEWATCH III cost analysis

WAVEWATCH III version 3.14

Figure: WAVEWATCH III cost versus spatial resolution using same time step and fixed spectral grid (25f × 24θ).

• Let N be the number of grid cells.

Then running costs are Cost(N) = αN1.07.

• Using source terms approximately

doubles cost since log2[Costsrc] ≈ 1 + log2[Costin] Costsrc ≈ 2 Costin

• Changing grid resolution can have a

dramatic effect on speed Nglobal = Nx × Ny Ncoarse = (Nx/3.2) × (Ny/3.2) Costglobal ≈ 10 Costcoarse

12 Adrean Webb / 39

SLIDE 13

Outline

• 1. Scientific background
• 2. Coupled wave component
• 3. Overview of meshless approach
• 4. Kinematic prototype
• 5. Discussion

13 Adrean Webb / 39

SLIDE 14

RBF and RBF-FD methods

Figure: Example node layout of a bumpy sphere (Fuselier & Wright,

2012) and possible solid body rotation

• ver the surface (Flyer & Wright, 2007)

Radial Basis Functions (RBF):

• Uses a meshless node layout
• Solves advective problems with spectral accuracy
• Allows geometric flexibility and local node

refinement

• No advective singularities
• Smooth solutions require much fewer unknowns

compared with other methods RBF-Generated Finite Differences (RBF-FD):

• Uses local RBF-generated finite difference

formulas to reduce memory count from O(N2) to O(N)

• Well-suited for parallelization (Bollig et al., 2012)

14 Adrean Webb / 39

SLIDE 15

15 Adrean Webb / 39

Solving PDEs with the local RBF-FD method

Local RBF-FD Method:

• Creates a set of FD weights (w) that are exact at

point ~ ↵s for the RBF () centered at each node in the subset      11 · · · 1n 1 . . . ... . . . . . . n1 · · · nn 1 1 · · · 1           w1 . . . wn wn+1      =      (L)s1 . . . (L)sn      where ji = (k~ ↵ ~ ↵ik)|~

↵=~ ↵j

• Local weights are combined into a banded sparse

matrix D such that  

LW(~ ↵1, t) . . . LW(~ ↵N, t)

  ⇡  

D11 · · · D1N . . . ... . . . DN1 · · · DNN

   

W(~ ↵1, t) . . . W(~ ↵N, t)

 

Figure: Sample (a) local differenti- ation weights and (b) global banded matrix (Flyer et al., 2012)

SLIDE 16

Stabilizing the RBF-FD method with hyperviscosity

Hyperviscosity Filter:

• A high-order Laplacian operator is used to stabilize the RBF-FD method

∂tW = −DW + γ∆pW

• The artificial diffusion coefficient is typically in the range of 10−20 to 10−45

Figure: Sample eigenvalue spectrum with the RK4 stability domain (solid line) (a) without hyperviscosity and (b) with a ∆4-type hyperviscosity added. (Flyer et. al., 2012)

16 Adrean Webb / 39

SLIDE 17

Outline

• 1. Scientific background
• 2. Coupled wave component
• 3. Overview of meshless approach
• 4. Kinematic prototype
• 5. Discussion

17 Adrean Webb / 39

SLIDE 18

Coupled model domains

Spatial: Directional-Frequency: WAVEWATCH III: RBF-FD WAVE:

18 Adrean Webb / 39

SLIDE 19

RBF-FD-W model (full-4D)

RBF-FD-W Details:

• (r) = e(εr)2 for ", r 2 R⇤

+

• ~

↵ = (~ x, ~ k) = (x, y, z, kx, ky, kz) m ~ ↵ = (~ x, ~ f ) = (x, y, z, ✓x, ✓y, f )

• Polar ice caps for latitudes above

and below ±75

• For constant depth with no

currents, the wave action balance equation is given below:

Figure: Sample coupled spatial and directional-frequency do- main ∂tW + ˆ cg (kz ) p x2 + y2 n ykx xzky

• ∂x +
• xkx yzky
• ∂y + ky
• x2 + y2

∂z +

• zkx ky
• ∂kx +

⇣ zk2

x

⌘ ∂ky

• W

= Sources

19 Adrean Webb / 39

SLIDE 20

20 Adrean Webb / 39

RBF-FD-W stencil selection (full-3D)

Convergence rates are determined by both the total number of nodes and local stencil size

(a) Spatial stencil (b) Directional stencil

Figure: Sample spatial and directional stencils for spatial node (1, 0, 0) and dominant direction π/6. The combined stencil is 17~

x × 9~ k = 153~ ↵.

SLIDE 21

Tuning: Spectral domain

21 Adrean Webb / 39

Spectral domain tuning details:

• Stencil sizes: 5, 7, 9, 11
• Global nodes: 24, 36, 48, 60
• Initial bell widths: 60 - 180 deg

Initial conditions: Idealized spectra:

• n

Equation Normalization W0(θ) = exp[−(9θ/2d)2] (2√πd/9) erf[9π/2d] W0(θ) = cos2[πθ/d] d/2 for |θ| < d/2; 0 otherwise Type Ab Gaussian bell G Cosine bell C

SLIDE 22

Tuning: Spectral domain

22 Adrean Webb / 39

• 180
• 90

90 180 θ 0.2 0.4 0.6 0.8 1 1.2 Wave action dθ = 120◦ 0.2 0.4 0.6 0.8 1 a = ∆t/(∆θ/vθ) 10-6 10-4 10-2 100 Normalized ℓ2 error dθ = 120◦, ε = 3, Nθ = 120

nθ = 5 nθ = 7 nθ = 9 nθ = 11

1 2 3 4 5 ε 10-5 10-4 10-3 10-2 10-1 100 Normalized ℓ2 error a = 0.1, dθ = 120◦, Nθ = 36

nθ = 5 nθ = 7 nθ = 9 nθ = 11

101 102 Nθ 10-8 10-6 10-4 10-2 100 Normalized ℓ2 error a = 0.1, dθ = 120◦, ε = 3

nθ = 5 nθ = 7 nθ = 9 nθ = 11

1 2 3 4 5 ε 10-5 10-4 10-3 10-2 10-1 100 Normalized ℓ2 error a = 0.1, dθ = 120◦, nθ = 9

Nθ = 24 Nθ = 36 Nθ = 48 Nθ = 60

5 10 15 20 b = dθ/∆θ 10-3 10-2 10-1 100 Normalized ℓ2 error a = 0.1, ε = 3, nθ = 9

∆θ = 15◦ (Nθ = 24) ∆θ = 10◦ (Nθ = 36) ∆θ = 7.5◦ (Nθ = 48) ∆θ = 6◦ (Nθ = 60)

Figure: Spectral advection test with (a) Gaussian bell initial condition and (b-f) relative errors after 1/4 revolution.

`2

SLIDE 23

Tuning: Spatial domain

23 Adrean Webb / 39

Spatial domain tuning details:

• Stencil sizes: 17, 31, 50
• Global nodes: 2500, 3600,

4900

• Initial bell widths: 45 - 60 deg
• n

Equation Normalization W0(θ) = exp[−(9θ/2d)2] (2√πd/9) erf[9π/2d] W0(θ) = cos2[πθ/d] d/2 for |θ| < d/2; 0 otherwise

Initial conditions:

Type Ab Gaussian bell G Cosine bell C

• 1
• 0.5
• 1

0.5 z 1

• 0.5

1 x 0.5 y 0.5

• 0.5

1

• 1
• 1
• 0.5
• 1

0.5 z 1

• 0.5

1 x 0.5 y 0.5

• 0.5

1

• 1

Width = 30 deg Width = 60 deg

SLIDE 24

Tuning: Spatial domain

24 Adrean Webb / 39

0.2 0.4 0.6 0.8 1 a⃗

x = ∆t/(∆x/cg)

10-4 10-3 10-2 10-1 100 Normalized ℓ2 error d⃗

x = 60◦, ε = 3, N⃗ x = 3600

n⃗

x = 17

n⃗

x = 31

n⃗

x = 50

1 2 3 4 5 ε 10-4 10-3 10-2 10-1 100 Normalized ℓ2 error a⃗

x = 0.1, d⃗ x = 60◦, N⃗ x = 3600

n⃗

x = 17

n⃗

x = 31

n⃗

x = 50

1 2 3 4 5 ε 10-3 10-2 10-1 100 Normalized ℓ2 error a⃗

x = 0.1, d⃗ x = 60◦, n⃗ x = 17

N⃗

x = 2500

N⃗

x = 3600

N⃗

x = 4900

103 104 N⃗

x

10-4 10-3 10-2 10-1 100 Normalized ℓ2 error a⃗

x = 0.1, d⃗ x = 60◦, ε = 3

n⃗

x = 17

n⃗

x = 31

n⃗

x = 50

• 1
• 0.5
• 1

0.5 z 1

• 0.5

1 x 0.5 y 0.5

• 0.5

1

• 1

Figure: Spatial advection test along equator with (a) Gaussian bell initial condition and (b-f) relative errors after 1/2

revolution. `2

5 10 15 20 25 30 35 40 b = d⃗

x/∆⃗

x 10-4 10-3 10-2 10-1 100 Normalized ℓ2 error a = 0.1, ε = 3, n⃗

x = 17

∆⃗ x = 7.2◦ (N⃗

x = 2500)

∆⃗ x = 6◦ (N⃗

x = 3600)

∆⃗ x = 5.1429◦ (N⃗

x = 4900)

SLIDE 25

25 Adrean Webb / 39

Advection in the coupled domains

Test initialized with a spatial and directional Gaussian bell and a 30 dominant direction 3600~

x × 36~ k global nodes with 17~ x × 9~ k stencil Figure: (a) Node set with advection path (red), peak Gaussian bell edge (blue), and ice edges (black); (b) initial spatial Gaussian profile; (c) directional spread about dominant direction

SLIDE 26

26 Adrean Webb / 39

Spectral advection in coupled domain

3600~

x × 36~ k global nodes with 17~ x × 9~ k stencil

Figure: Select directional node values initialized with a spatial and directional Gaussian bell (60 width).

SLIDE 27

Initial direction error comparison

Figure: Relative `2

2 errors after 1/2 revolution for select initial directions.

4.36 % 1.41 % 13.93 % 3.16 %

SLIDE 28

Global node error comparison

Figure: Relative `2

2 errors after 1/2 revolution for select global node sets.

4.36 % 5.39 % 21.77 % 21.49 %

SLIDE 29

Moment advection in the coupled domain

3600~

x × 36~ k global nodes with 17~ x × 9~ k stencil

Figure: Integrated wave action (Hm2

0/16) initialized with a spatial and directional Gaussian bell (45 and 120 width)

and 20 dominant direction

29 Adrean Webb / 39

SLIDE 30

Coupled model domains

Spatial: Directional-Frequency: WAVE- WATCH III: RBF-FD WAVE:

30 Adrean Webb / 39

SLIDE 31

WAVEWATCH III comparison

N~

x = 4320,

revolution = 1/2, dominant direction = −30 Solution WW3

Figure: Selected directional grid cell value initialized with a spatial Gaussian (57.2 width) and directional cosine- 20-power (64 width) after 1/2 revolution. The significant wave height is 2.5 m.

31 Adrean Webb / 39

SLIDE 32

WAVEWATCH III comparison

N~

x = 44064,

revolution = 1/2, dominant direction = −30 Solution WW3

Figure: Selected directional grid cell value initialized with a spatial Gaussian (57.2 width) and directional cosine- 20-power (64 width) after 1/2 revolution. The significant wave height is 2.5 m.

32 Adrean Webb / 39

SLIDE 33

WAVEWATCH III error

Relative total `2

2 errors over all directions

WW3 WW3

Figure: Relative total `2

2 errors after 1/2 revolution for different WW3 grid resolutions.

33 Adrean Webb / 39

63.25 % 44.72 %

SLIDE 34

RBF-FD-W versus WAVEWATCH III

Relative total `2

2 errors over all directions

RBF-FD-W WW3

Figure: Comparison of relative total `2

2 errors after 1/2 revolution between RBF-FD-WAVE and WW3.

34 Adrean Webb / 39

44.72 % 4.24 %

SLIDE 35

Outline

• 1. Scientific background
• 2. Coupled wave component
• 3. Overview of meshless approach
• 4. Kinematic prototype
• 5. Discussion

35 Adrean Webb / 39

SLIDE 36

Wave action balance equation

The wave action balance equation can be derived on a coupled slab (~ x, ~ k 2 R2) from the governing linear wave equations and is given by @tW + r

~ k Ω · r ~ x W r ~ x Ω · r ~ k W = Sources

where

~ x, ~ k

• =

n g|~ k| tanh h |~ k|H(~ x) i

• 1/2

is the intermediate-water dispersion relation

• Ω

~ x, ~ k

• = ~

U ~ x

• · ~

k + ~ x, ~ k

• is the Doppler-shifted dispersion relation
• W

~ x, ~ k, t

• = g S~

k

~ k; ~ x, t

• /

~ x, ~ k

• is the spectral adiabatic invariant
• “Sources” encompass the non-kinematic physics of the waves

(generation, dissipation, nonlinear interactions, etc.)

36 Adrean Webb / 39

Spectral changes Spectral transport

SLIDE 37

Spatial refinement

37 Adrean Webb / 39

Figure: Comparison of current and potential spatial refinement tuning.

5 10 15 20 25 30 35 40 b = d⃗

x/∆⃗

x 10-4 10-3 10-2 10-1 100 Normalized ℓ2 error a = 0.1, ε = 3, n⃗

x = 17

∆⃗ x = 7.2◦ (N⃗

x = 2500)

∆⃗ x = 6◦ (N⃗

x = 3600)

∆⃗ x = 5.1429◦ (N⃗

x = 4900)

1 2 3 4 5 ε 10-4 10-3 10-2 10-1 100 Normalized ℓ2 error a⃗

x = 0.1, d⃗ x = 60◦, N⃗ x = 3600

n⃗

x = 17

n⃗

x = 31

n⃗

x = 50

1 2 3 4 5 ε 10-5 10-4 10-3 10-2 10-1 100 Normalized ℓ2 error a = 0.1, dθ = 120◦, Nθ = 36

nθ = 5 nθ = 7 nθ = 9 nθ = 11

3 4 5 6 7 ε 10-4 10-3 10-2 10-1 100 Normalized ℓ2 error a⃗

x = 0.1, d⃗ x = 30◦, N⃗ x = 4900

n⃗

x = 17

n⃗

x = 31

n⃗

x = 50

5 10 15 20 25 30 35 40 b = d⃗

x/∆⃗

x 10-6 10-4 10-2 100 Normalized ℓ2 error a = 0.1, ε = 5, n⃗

x = 31

∆⃗ x = 6◦ (N⃗

x = 3600)

∆⃗ x = 5.1429◦ (N⃗

x = 4900)

∆⃗ x = 4.5◦ (N⃗

x = 6400)

3 4 5 6 7 ε 10-5 10-4 10-3 10-2 10-1 100 Normalized ℓ2 error a = 0.1, dθ = 120◦, Nθ = 60

nθ = 5 nθ = 7 nθ = 9 nθ = 11

Width = 60 deg Width = 30 deg Spectral Spatial Spatial

SLIDE 38

Next steps

Stage 2

• Replace current attenuation filter approach with physical boundaries

Stage 3

• Add source terms and wavenumber shifting (via time-splitting)

Figure: Examples of (a) boundary attenuation filter and (b) variable spatial node density.

38 Adrean Webb / 39

SLIDE 39

Summary and conclusions

• WAVEWATCH III is now coupled to the NCAR Community Earth System Model

(3.2× 4 spatial, 25f × 24θ spectral grid)

• A meshless approach using RBF-FDs is a viable alternative for future spectral

wave models (particularly for use in global climate model simulations)

• A monochromatic kinematic prototype performs well with limited computation

(e.g., runs on a laptop)

• In simple tests the prototype compares well against higher resolution

WAVEWATCH III runs: requires fewer unknowns (∼ 1/12) and is more accurate (∼ 10×)

39 Adrean Webb / 39

SLIDE 40

Thank You!

Figure: A modern reinterpretation of Hokusai’s “The Great Wave” (Murakami, “727”).

40