Transcript Document

-SELFE Users Trainng Course (2) SELFE numerical formulation

Joseph Zhang, CMOP, OHSU

1

Numerical methods for pde

  1!

(

f

2!

(  ) 2  

f n

!

(  )

n

f

(

n

 1) (

n

 1)!

(  )

n

 1 2        System of algebraic equations Truncation error consistency, and order of convergence (Lax) Stability Higher-order approximation Upwind

c j n

 1 

c j n

D

t

 

u c n j

 D

x c j n

 1 

c j n

 1 

c n j

u

j n

c j n

 1  1 j-1

u

j Finite element Residual Weighting function Galerkin FE Least square j+1 D

x

  0 

L u

ˆ  ,

w

  

wd

  0

w

 

w

  Basis/shape function (satisfied at finite number of points) Positive definite operator  positive definite matrix Symmetric mass matrix

N x

Strong vs. weak formulation

    

Shape function

Used to approximate the unknown function Although usually used within each individual elements, shape function is global not local!

u

i N

  1

u i

i

, 

i

(

x j

)  

ij

Must be sufficiently smooth to allow integration by part in weak formulation Mapping between local and global coordinates Assembly of global matrix 3 

L

 

i N

  1

u

 

i i

j d

 

M

m

 1  

m L

 

N

i

 1

u

 

i i

j d

 

m

0,

j

 1,...,

N

j Conformal Non-conformal 1 1 j-1 j j+1 j-1 j j+1

x N x N

Semi-implicit Eulerian-Lagrangian Finite Element (SELFE)

         A formal semi-implicit finite-element framework  A key step is to decouple the continuity and momentum equations through the bottom boundary layer Unstructured grid in the horizontal Hybrid SZ in the vertical Eulerian-Lagrangian method (ELM) for momentum advection Volume conservation is not enforced numerically (but good) Numerical efficiency: semi-implicit time stepping; ELM   Mild stability constraint: horizontal viscosity, baroclinic pressure Moderately more expensive than ELCIRC Convergence rate: 2 nd order for uniform grid  Optional higher-order numerics for momentum and transport Treats inundation/drying in a natural way (validated for inundation benchmarks) Bed deformation scheme incorporated (e.g., tsunami) 4

Horizontal grid

(

i,j

) 5

k z

Vertical grid (1)

Pure S zone

N z

s zone

h c

S zone

h s

SZ zone

z

s =0

h s S

-levels s =-1

k z

 1     

C h

 2 1

Z

-levels

s

)

z

  (1  s ) 

h c

s  (

h

 

b

) sinh(  s sinh 

f

)  

b

tanh

f

s

c

 1/ 2) 2 tanh( 

f

 tanh( 

f

/ 2) / 2) s 0) (0  

b

 1; 0  

f

 20) 6

S

-coordinates

 f =10 -3 ,  b =1  f =10,  b =0  f =10,  b =1

h c

 f =10,  b =0.7

h s

7

z

= 

h s

Vertical grid (2)

8

Vertical grid (3)

9 • 3D computationa unit: uneven prisms •

S

-coordinates (Song and Haidvogel 1994) are ideal for shallow region • •

Z

-coordinates are necessary to “stabilize” the deep region s -coordinates are used where

S

-coordinates are invalid • Equations are not transformed into

S

- or s -coordinates, but solved in their original (

Z

) form • Interpolation mode: along • Hydrostatic consistency:

Z

  or

S

(in pure

S

   region) • pressure Jacobian with higher-order integration •

Z

-method (preferred)

P

ˆ

n

k

+1

k

+1

S

ˆ

k z

2

s

2

z

1

s

1

Semi-implicit scheme

10 Continuity Momentum b.c.

D

u

Dt

n

 1  D

t

n

  

h

u

n

 1

dz

 )   

h

u

n dz

 0

u

n

 1 

u

* D

t

f

n

 

n

 1 

g

(1     

n n

u

z

u

n

 1

n

 1 

z

 

τ

n

 1 , at

z w

n

u

b n

 1 , at

z

 

n

;  

h

,  

n

+  

z

  

n

n

u

n

 1 

z

  

C D

|

u

b n

|

u

n

 1 D

t

u

* ;

u

*

t t

) dispersion diffusion   Implicit treatment of divergence and pressure gradient terms by-passes most severe CFL condition ELM: takes advantage of both Lagrangian and Eulerian methods   Grid is fixed in time, and time step is not limited by Courant number condition Advections are evaluated by following a particle that starts at certain point at time t and ends right at a pre-given point at time t+ D t.

D

x

   The process of finding the starting point of the path (foot of characteristic line) is called backtracking, which is done by integrating d x/ dt =u 3 backward in time.

To better capture the particle movement, the backward integration is often carried out in small sub-time steps ( D t / N).

 Simple backward Euler method as an option   5th-order embedded R-K method as an alternative Interpolation-ELM does not conserve mass; integration ELM does Numerical diffusion vs. dispersion Other considerations • Implicitness factor 0.5 ≤  ≤1 to ensure stability • Explicit treatment: Coriolis, baroclinicity, horizontal viscosity, part of bottom velocity

t

+ D

t

    Characteristic line

t

 1

  

Finite-element formulation (1)

A Galerkin weighted residual statement for the continuity equation:   

i

n

 1 D 

t

n d

       

i

U

n

 1

d

   

v

i U

ˆ

n n

 1

d

(1   )     

i

i

U

n d

    

i n U d v

+

U

  

u

dz

h

Vertical integration of the momentum equation:

v

 

v

i U n n

 1

d

0, (

i

 1,...,

N p

; 0.5

1)  

U

n

 1 

G

n

H n n

n

 1  

n

D

t

u

n b

 1 

C D

|

u

b n

|;

G

: explicit terms

u

b

: bottom velocity Momentum equation applied to the bottom boundary layer:

u

b n

 1  D

t

u

*

b

f

b n

 

n

 1 

g

(1  

n

+  

z

  

n

u

n

 1 

z

  , at

z

 

b

h

via  

u

z

  0 ln( 

b

/

z

0 )

C

1/ 2

D

|

u

b

|

u

b

, (

z

0

z

b

h

) 

b

u

b

11 z =  h

Finite-element formulation (cont’d)

Vertically integrated velocity becomes:

U

n

 1 

n

 

n

H n

 

n

D

t n

n

 1 ,    Finally, substituting this eq. back to continuity eq. we get one equation for elevations alone: I 1 I 2 I 3   

i n

 1 

g

 2 D 2 ˆ

n

 

i n

 1  

d

 

g

 2 D

t

2  

v

i n

 

n

 1 

n d v

t

 

v

i U

ˆ

n n

 1

d

 

v I n

12

I n

   

i n

 )

t

I 4 

i

U

n

 

i n

d

 )

t

  

i n U d n

I 5

t

 

v

i

I 6

n d

v

Notes: • When node i is on boundaries where essential b.c. is imposed, equations are needed, and so I 2 • When node i the last term is known  and I 6 only the first term is truly unknown!

 is known and so no need not be evaluated there is on boundaries where natural b.c. is imposed, the velocity is known and so

Matrices: I

1

I

1 

N i j

 

j

(  

i

'

n

 1 

g

 2 D ˆ  ˆ

i

'

i

’ 

j N i

 1

l

3   1 

n

 1   1   12

A j

g

 2 D

t

4

A j

is local index of node I in  j .

n

 1 )

d

j

   reaches its max. when

i

’=

l

, and so does I 1 .

i

3 

i

'  0 so diagonal is dominant if the averaged friction-reduced depth ≥0 (can be relaxed) Mass matrix is positive definite and symmetric!

(

i’,1

)

i

i

> ' i  j (

i’,2

) 13

Matrices: I

3

and I

5

I

3 

v

  

i U

ˆ

n n

 1

d

 

v

j

 

ij

i U

ˆ

n n

 1

d

ij

  

j L

4

ij N z

 1 

bs

D

z

 1 (

u

ˆ n

n

 1 )  1  (

u

ˆ n

n

 1 ) Flather b.c. (overbar denotes mean)

u

ˆ n

n

 1 

u

n

I

3   

j L ij

  (

U

n 2 )

ij

n

 1   ) 

gh ij

6 2( 

i n

 1  

i

)  ( 

j n

 1  

j

)   j  ij

n

i j Similarly

I

5  

j

 

ij

 ˆ ˆ

i

'

U d n n

 

ij

 

j L

4

ij N z

 1 

bs

D

z

 1 (

u

ˆ n

n

)  1  (

u

ˆ n

n

) 14

Matrices: I

4 Has most complex form 15

I

3      

i n

'

j

  

j

 )

t

i

U

n

t

i d

 

j

   

j

h

u

*

dzd

  D

j t

 

j

F

1

d

 

d

 

N i

 

j

  

A j

12

l

3  1 

n

(1   ) (1  )

tA j

 ˆ

i

'

U

n j j A t j

  

F

2 

τ

w

  (

u

*

b

 D

tf

2

b

 D

tf

1

b

) 

g

(1   )  

i

G

'

j

    

n

   We have decomposed vector f into two parts: one for averaging and the other for integration

f f

1

f

2 Most complex part is the baroclinicity

f

c

 

g

 0 

z

    at elements/sides

Baroclinicity: pressure Jacobian method

• No longer used in the newest version • In    s S layers   

z

s  s      s  s

z

is evaluated using cubic spline in the vertical In Z layers, the derivative can be evaluated directly The vertical integration is done using higher-order rule (Song and Haidvogel) Advantage: no special treatment is needed near surface and bottom Disadvantage: pressure gradient errors 16

Baroclinicity:

Z

-coordinate method (preferred)

Interpolate density along horizontal planes  

z

Advantage: alleviate pressure gradient errors (Fortunato and Baptista 1996) Disadvantage: near surface or bottom Method: compute derivatives along two direction at node then convert them back to x,y .

i and

ELM transport

• Density is defined at nodes and whole levels • Density gradient calculated at side centers and whole levels • Either constant extrapolation (shallow) or more conservative method (deep) is used near bottom • Cubic spline is used in all interpolation of  (and necessary S,T)  ( ) • Solve 2 eqs for the density gradient ( 

x

' ) (

x

2 

x

1 )  ( 

y

' ) (

y

2 

y

1 )   ' 2   1 ' ( 

x

' ) (

x

4 

x

3 )  ( 

y

' ) (

y

4 

y

3 )   4 '   3 ' Horizontal boundary or dry interface: no flux b.c.

3 ( j,k ) 2 ( j,k ) 1 17 4

Baroclinicity: upwind/TVD transport

• Density is defined at prism centers (half levels) • Density gradient calculated at prism centers first (for continuity eq.); the values at face centers are averaged from those at prism centers (for momentum eq.) • Either constant extrapolation (shallow) or more conservative method (deep) is used near bottom to avoid spurious flow • Cubic spline is used in all interpolation of averaging for sides for 3 pairs; e.g.

 ( )  except for the • 3 eqs for the density gradient vs. 2 unknowns – averaging needed 1 (  ' ) (

x

1 

x

0 )  (  ' ) (

y

1 

y

0 )  1 ' 0 (  ' ) (

x

2 

x

0 )  (  ' ) (

y

2 

y

0 )   2 '   0 ' • degenerate cases • Vertical integration using trapzoidal rule i 0 2 18 3

Horizontal viscosity

Evaluate it at 3 side centers, and integrate using linear quadrature 

u

)       

u

n d

d

 

l

 1

A l



l

3

i

 1

s

   

u

n

 The normal derivative is computed using average between adjacent elements   

u

n

  

l

  

u

n

 

l

Use linear shape function within each element and chain rule to calculate the derivatives 

u

n

 

u

n

s  s side j 19

Momentum equation

  

b

h

l

(

l

 ( u,v ) solved at side centers and whole levels • Easy of imposing b.c.

• Staggering for stability Galerkin FE in the vertical

k b

  

u

n

 1  1,  D

t

 

z

   

N z

) 

u

n

 1 

z

     

dz

  

b

 

h

l

u

*

n

 1  D

t

 

f

 

n

 1 

g

(1

g

Solution at bottom from the b.c. there (u=0 for  ≠0)

D

u

Dt

 

f

 +  

z

  

u

z

 

u

z

u

z

τ

w

 

u

b

, at

z

  

h

 

n

  

dz

( A ( 

z N z

 D

t

τ

n w l

)     1 D

z l

6  1 + A ( (2

u

N z l n

 1 

l

) 

u

l n

 1  1 D

z l

 1 6 ) (2  

g

l l

  1/ 2 D

t

g

l

 1 )

u

l n

  1 1 D 

z l

u

 1

l n

 1  A (

k b

    A (

l

) D

z l

6

k b

(2

g

l n

 1

l

)     D

z l

6

g

l n

  1 1 ) (2

u

l n

 1

l

 

b

1,

N z

) 

u

l n

 1  1 )  

l

 1/ 2 D

t u l n

 1  D

z l u l n

  1 1     

b

 1 D

t

n

u

k b

 1  1  The matrix is tri-diagonal z = h  b 20 N z k b +2 k b +1 k b

Velocity at nodes

Needed for ELM (interpolation) • Method 1: averaging around ball (most diffusive) • Method 2: use linear shape function (conformal or non conformal); optional averaging 3

u

1 

u

II

u

III

u

I

II • Filtering of modes • Needs velocity b.c.

1 III I • Linear interpolation at foot of char. Line • Conformal • Non-conformal – discontinuity along each side 2

i

21

Shaprio filter

Few ocean models are truly free from spurious numerical modes (myriads of processes) A simple filter to reduce sub-grid scale (unphysical) oscillations while leaving the physical signals intact ( a =0.5 optimal) 0 

u

0  a 4  

i

4   1

u

i

 4

u

0   , (0 0.5) No filtering for boundary sides – need b.c. there especially for incoming flow Pre-processor to check the geometry (ipre=1 with 1 processor) Violations usually occur at boundary (fort.11) 1 2 0 3 4 22 internal boundary extreme

ELM with Kriging

• Best linear unbiased estimator for a random function • “Exact” interpolator • Works well on unstructured grid f (x)  a a 1  2

x

 a 3

y

i N

  1 

i K

(|

x

x

i

|) Drift function Fluctuation K is called generalized covariance function Minimizing the variance of the fluctuation we get

i N

  1 

i

 0

i N

  1 

i x i

 0

i N

  1 

i y i

 0

f x y i

) 

f i

  2

h h

3 x 23 2-tier Kriging cloud

Vertical velocity

 Vertical velocity using Finite Volume

S

ˆ

k

 1 (

u k n

  1 1

n k x

 1 

v k n

  1 1

n k y

 1 

w n

 1  1

n k z

 1 )  ˆ

k m

3   1

q

ˆ

j k

u

s

n

j

(

q

ˆ

n

 1 

q

ˆ

n

 1

k n

 1

n k x

v k n

 1

n k y

w n

 1

n k z

)   1 ) / 2  0, (

k

k b

,...,

N z

 1) 24

w n

 1  1 (

k

 

k b k

 1 1 ˆ

k

 1 ,...,

N z

   1) 

m

3   1

s

(

q

ˆ

n

 1 

q

ˆ

n

 1  1 ) / 2 

S

ˆ

k

 1 (

u k n

  1 1

n k x

 1 

v k n

  1 1

n k y

 1 )  ˆ

k k n

 1

n k x

v k n

 1

n k y

w n

 1

n k z

 ) , Bottom b.c.

u k n

 1

n k x

v k n

 1

n k y

w n

 1

n k z

 

b

t

, (

k

k b

)

P

ˆ

n

k

+1

k

+1

S

ˆ

k

Inundation algorithms

  Algorithm 1  Update wet and dry elements, sides, and nodes at the end of each time step based on the newly computed elevations Algorithm 2 (accurate inundation for wetting and drying with sufficient grid resolution)     Wet  Dry  add elements & extrapolate velocity remove elements Iterate until the new interface is found at step n+1 Extrapolate elevation at final interface (smoothing effects) dry 25 A  n B wet Extrapolation of surface

Transport: ELM

Galerkin FE applied to both nodes and sides • no horizontal diffusion   

h

l T n

 1  D

t T

*

dz

   

h

l

   

z

   

T

z n

 1    Apply mass lumping reduces FE to FD

n

 

Dc Dt

   

c

 

z c

 

z

n

    

c

z

 

Q c

,

z

0,

z

    ( 

h

c

),

c

  

h

26 D

z l

 1/ 2

T l n

 1  

l

 1/ 2 D

t n

 1

T l

 1 

T l n

 1 D

z l

 1  

l

 1/ 2 D

t T l n

 1 

T l n

 1  1 D

z l

 D

z l

 1/ 2 (

T l

* 

Q l n

D

t

), (

l

 

b

1, D

z l

2

T l n

 1  

N

 1/ 2 D

t T l n

 1 

T l n

 1  1 D

z l

D

z l

2

T l n

 1  

l

 1/ 2 D

t T l n

 1  1 

T l n

 1 D

z l

 1 

T

ˆ

t

D

z l

2 (

T l

* 

Q l n

D

t

),

l

N z

D

z l

2 (

T l

* 

Q l n

D

t

),

l

k b

ELM: mass conservation not guaranteed • Linear with element splitting in the horizontal; cubic spline in the vertical • Quadratic in 3D (works only in estuary due to dispersion) • Impose bounds from surrounding nodes to alleviate dispersion • Kriging (coming up) • Vertical interpolation: placement of T and  • Upwind bias for quadratic interpolation Char. line ,

N z

 1)

T l

x

l l

+1

l

-1

l

-2

Transport: upwind (1)

T

t

  

T

 

z T

n

(

u

T

)   

z

  

T

z

    ( 

h

T

) 

T

ˆ,

z

   0,

z

 

h

Finite volume discretization in a prism ( i , k ): mass conservation

V i T i m

 1 D

t

 '

T i m

 5 

j



V u S T j j j

* 

V Q i m

D

t

' 3 

j

 1 ( 

h

)

j T j m

 

ij T i m S j

, (

k

k b

 1,  D

i

'     ,

N z

)

T m

 1  1 D

z

T m

 1  1/ 2    1

T m

 1 D

z

T

 1/ 2

m

 1  1     S j S S k T i,k k 1 m≠n, D t’ ≠ T i = T i,k; D t; u j is outward normal velocity 27 Focus on the advection first

T i m

 1 

T i m

 D

t V i n

' 

j



V j

* From continuity equation

Q j

u S j j

Q

 |

Q

|   |

Q

|

j j j

 

j

S + : outflow faces; S  : inflow faces Upwinding

T j

*  

T T i j

,

j

,

j

 

S

S

Transport: upwind (2)

Drop “m” for brevity

T i m

 1 (1  a 

T i

)

T i

 D

t

'

V i n

    a

T j

  |

j i

   |

j j

      1 Max. principle is guaranteed if D

t

'

V i

  |

Q j

|   

T i

 D

t

'

V i

  |

j j

 S j S S T i,k k 1 1  D

t

'

V i

  |

Q j

|

t

'

V i

  |

Q j

| Courant number condition (3D case) Global time step for all prisms Implicit treatment of vertical advective fluxes to bypass vertical Courant number restriction in shallow depths

V i T i m

 1 D

t

 '

T i m

k

2   1

u S T k k k m

 1 *  

j

3   1

u S T V Q j j j

*

i n

ˆ

n

 D

i

'    

T m

 1  1 D

z

T m

 1  1/ 2    1

T m

 1 D

z

T m

 1  1  1/ 2     D

t

'

j

3   1 ( 

h

)

j T j m

 

ij T i m S j

V Q i m

, (

k

k b

 1, ,

N z

) (Modified Courant number condition) Mass conservation • Budget – vertical and horizontal sums • movement of free surface: violation of volume conservation k 28

Solar radiation

Q

ˆ  1  0

C p

H

z

total solar radiation 

H

(0) Re  / 1    2 Budget   0 

H

(0)  0

C p

The source term in the upwind scheme (F.D.)

i i m

 

V i

0

C p H m b

 0

H m

 D

z H m

 1 

i m

H m

 1 )  0

C p

to ensure total heat is conserved (black body) This can cause overheating in shallow depths: adjust albedo 29

Transport: TVD (1)

Start from discretized advection Now include downwind component '

T j

*   Flux limiter function (depends only on face upwind j ); Downwind and central schemes; The key is to find an appropriate function that does not violate max principle

T i m T j

* '    D

t V i n

2

j

' 

T

j



V jD j

* 

T j

*  upwind 0   

j

 2

T i m

 1 

T i m

 D

t V i

 

j

S

|

Q j

| (

T j T i m

 1 

T i

)  D

t V i

 

j

S Q j

2

j

(

T jD

T j

* ) 

T i m

 D

t V i

'  

j

S

|

Q j

|   1   2

j

  (

T j

T i

)  D

t V i

'  

j

S

|

Q j

| 

j

2 (

T i

T j

) Sum of coefficients = 1; so max principle is satisfied if and only if the coefficient of negative T i is non Put the last term in similar form as the upwind term D

t V i

' 

j

S

 |

Q j

| 2

j

(

T i

T j

)  D

t V i

'  

j

S

 |

Q j

| 

j

2  (

T i

T j

)    

j

S

|

Q j

| (

T j

T i

)  30  D

V i t

'  

j

S

|

Q j

| 

j

2 (

T i

T j

)  D

t

' 

V i

 

j

S

|

Q j

| (

T j

T i

)

Transport: TVD (2)

T i m

 1 

T i T i

 D

V i t

'   |

Q j

 | 1  

j

2    1  D

V i t

'   |

Q j

 | 1  

j

2   (

T j

T i

)  D

t

' 

V i

  |

Q j

| 

j

2 (

T j

T i

)          D

t

'

V i

  |

Q j

 | 1  

j

2   

T j

Max principle is satisfied if (Courant number condition)   1    |

Q j

| 

j

2 

T i

T j

    |

Q j

| 

j

2   |

Q m

T i

| 

T m

T j

 

T i

    

j

2

r j r j

 Since   |

Q m

|

Q j

| | 

T m

T i

T j

T i

  Upwind ratio r 0  

j

 2, 0  

j r j

 2, 0   

N

 The Courant number condition is

t V i

(1 

N

 )   |

Q j

| as compared to upwind

t V i

  |

Q j

|   D

t

'

V i

  1  

j

2 |

Q j

 | 1  

j

2 0      0 j S + i S  m 31

Transport: TVD (3)

TVD scheme • Mass conservative with max principle; second-order accuracy • Fully explicit scheme due to nonlinearity – vertical Courant number most severe • Boundaries and wet/dry interface – revert back to upwind • Other higher-order schemes: MUSCL; WENO Choice of flux limiting function (gradient detector)         max(0, min(1, 2 ), min(2, )) 

r

  Super Bee (compressive) Minmod (diffusive) Osher Van Leer 32 Second order TVD region (Sweby 1984)

Turbulence closure (1)

Dk Dt

  

z

  

k

k

z

  

M

2  

N

2  

D

   

z

    

z

  

c

  2 

c

 

N

2  

k l

     1 (

c

0  0 ) 2  D ,

z

(

c

0  )

p k

Essential b.c.

33 

u

z

 

h

or 

m

(  ,

z

0 D )

n

 

h

,

z

or   

h

or 

wall

  and  

k

    

k

 

z

  

z

 

z

 Natural b.c.

0,

z

 

h

or      0 

n

  0

n

   ,

z l

l

  ,

z

h

 A (

l

N z

)   D

z l

 1 6  2

k l n

 1 

k l n

  1 1  2

k l n

A (

l

k b

)   D

z l

6  2

k l n

 1 

k l n

  1 1  2

k l n

k l n

 1   (  

k

)

l

 1/ 2 D

t k l n

  1 1 D 

z l k

 1

l n

 1    

k l n

 1   (  

k

)

l

 1/ 2 D

t k l n

 1  D

z l k l n

  1 1    ( 

b

, , A(

l

N z

) D

t

        D

z l

 1 6

k l n

 1/ 2   D

z l

 1 2

M

2  

M

2  

N

2  

N

2 

l n

 1/ 2

l

n

 1/ 2 (2

k l n

 1 A (

l

k b

) D

t

        D

z l

6

k l n

 1/ 2 D

z l

2  

M

2  

M

2  

N

2  

N

2 

l n

 1/ 2

l

n

 1/ 2 (2

k l n

 1 

k l n

  1 1 )      (

c

0  ) 3 

k

1/ 2

l

 1 

n l

 1/ 2 

k l n

  1 1 )      (

c

0  ) 3 

k

1/ 2

l

 1 

l n

 1/ 2 D

z l

6 D

z l

 1 (2

k l n

 1 6 

k l n

  1 1   )    (2

k l n

 1 

k l n

 1  1   )  

N z

)

Turbulence closure (2)

A (

l

A (

l

N z

)   D

z l

6  1 

k b

)   D

z l

6  2 

l

 2 

l n

 1

n

 1  

l n

  1 1  

l n

  1 1  2 

l n

 2 

l n

 

l n

 1   

l n

 1   (    (   )

l

 1/ 2 D

t

l n

  1 1 D 

z l

  1

l n

 1    )

l

 1/ 2 D

t

l n

 1  

l n

  1 1 D

z l

   A(

l

N z

) D

t

           (

c

6 D 0  )

z l

3 D

k l n

 1/ 2   1

z

2

l c

   1 2

c

k

c

 1  1/ 2 1 

M M

2  2 

l F wall

c

c

 3  

l n

 1/ 2 3 

N

2

N

l

D

z l

6 2

n

 1 

l n

 1/ 2  1/ 2 (2  (2 

l n

 1

l n

 1 

l n

 1/ 2   

l n

 1  1

l n

  1 1 ) )     A (

l

k b

) D

t

            (

c

D

z l

D 6

k l n

 1/ 2 0  ) 3  2

c z

l

 2 

c

k c

 1  1/ 2 1 

M M

2 2   

l F wall c

c

 3  3  

n l

 1/ 2

N N

2 6 2  D

z l l n

l n

 1/ 2  1/ 2 (2  (

l

k l n

 1/ 2 2 

n

 1

l n

 1    

l n

 1  1 )

l n

  1 1 )                    (

l

k b

, ,

N z

) 34

   

GOTM turbulence library

General Ocean Turbulence Model One-dimensional water column model for marine and limnological applications Coupled to a choice of traditional as well as state-of-the-art parameterizations for vertical turbulent mixing (including KPP) Some perplexing problems on some platforms – robustness?

  Bugs page Division by 0?

35 www.gotm.net

   

Numerical stability

Semi-implicitness circumvents CFL (most severe) (Casulli and Cattani 1994) ELM bypasses Courant number condition for advection Implicit transport scheme along the vertical bypasses Courant number condition Explicit terms   Baroclinicity  internal Courant number restriction Horizontal viscosity/diffusivity  D

x

diffusion number condition   Upwind and TVD transport – Courant number condition Coriolis – stable but prone to “modes”  1  D

t

D

x

2  0.5

3 2

u

1

n

 1 

u

2

n

 1 

u

3

n

 1  0 36 1 2

Lateral b.c. and nudging

37  S&T u,v Variable Type 1 (*.th) Time history; uniform across bnd Type 2 (param.in); Const.

Type 3 (param.in); tides Type 4 (*3D.th); Time history Type -1 (param.in); radiation Type -4 (*3D.th); nudging Nudging near bnd Via discharge Via discharge Clamp at i.c.

Tides (uniform across bnd) Nudge to i.c.

Nudge to *3D.th

Flather Nudge to uv3D.th (2 separate relaxation for in & outflow)

Summary

Time stepping Horizontal grid Vertical grid Numerical algorithm Convergence rate for non orthogonal grid Advection (momentum) Advection (transport) Volume conservation Wetting/drying ELCIRC Semi-implicit Orthogonal unstructured Z -coordinates Finite difference / finite volume Divergence ELM ELM/upwind Enforced Yes SELFE Semi-implicit unstructured Hybrid SZ coordinates Finite element / finite volume At least 1 st ELM (optional higher order Kriging) ELM/upwind/TVD Not enforced (but good) Yes 38