New Dynamical Core for AGCM

rawsulkyInternet και Εφαρμογές Web

11 Δεκ 2013 (πριν από 3 χρόνια και 8 μήνες)

83 εμφανίσεις

New Dynamical Core for AGCM
using Icosahedral
-
hexagonal Grids






Rashmi Mittal, H.C. Upadhyaya, F. Hourdin,

Om P Sharma

Indian Institute of Technology Delhi

Present Day Emphasis

Incorporation

of

high
-
resolution

observations


in

numerical

models

Data

structures

suiting

to

parallel

computing


Provision

for

local
-
refinement

/


adaptive

meshing

in

models

Increasing

the

resolution

is

key

to

improving

numerical

weather

prediction

by

models

High
-
resolution cloud resolving earth system models

And what we need ?

Mathematical Modelling Framework


Dynamical Models



Fluid Dynamics on a rotating sphere


Basic equations of motion


Advection of hydrometeors, trace
constituents


Discrete representation of system


Parameterizations of atmospheric
and land
-
surface processes

PHYSICS

DYNAMICS

Latitude
-
Longitude Grid


Problems

These grids are not
ISOTROPIC

The POLE PROBLEM: the
convergence of meridians
at poles demands a small
time step in order to
satisfy the CFL
requirement for
computational stability.

Well Proven good
solutions to pole problems
but poor SCALABILITY


Platonic solids
-

Regular grid structures

Source: Wikipedia




Best way to avoid the
POLE PROBLEM is to
use quasi
-
uniform grids.


Geodesic grids:
Tessellations of the
sphere generated from
platonic solids.


Platonic solids can be
enclosed in a

sphere

Source: NCAR



Source: Wikipedia




No Filtering

:
Quasi
-
uniform grid is
free of pole problems.


Most
Isotropic



Scalability
: Run on the JES ~ 2 km
global resolution (NICAM)


Absolutely simple
local
-
refinement

and
adaptive meshing



Why finite element/volume Icosahedral model ?

Traditional

finite

difference

methods

are

not

applicable

to

triangular/hexagonal

meshes
.

Spectral

models

require

global

communication

which

is

inefficient

on

MPP

with

distributed

memory
.

Icosahedral
-
Hexagonal grids
more attractive & could take
care of some problems :

History on Icosahedral Model


Sadourny et al.
1968:

non divergent barotropic vorticity equation

Masuda and Ohnishi
1986:

shallow water equations


Heikes and Randall 1995: shallow
-
water equations


Stuhne and Peltier 1999: shallow water equations.

Thuburn, 1999: shallow
-
water model

Ringler, Heikes and Randall 2000: atmospheric general circulation

H. Tomita et al. 2001 : Shallow Water Model

Ringler and Randall 2002: conserving …shallow
-
water equations

Majewski et al. 2002: operational global gridpoint model

Bonaventura and Ringler 2005: shallow water models


Lipscomb and Ringler 2005: transport scheme

Heikes, Konor, Randall, 2006: atmospheric model


Carfora 2007: Semi
-
Lagrangian advection

Satoh et.al 2007 :Nonhydrostatic Atmospheric Model



Transport Model



Shallow water model






Dry Dynamical Core



Dynamics + Physics ( without orography)
---

Aqua Planet ( Testing Phase)



Dynamics with full physics with orography
(Next Step)

Steps to be followed

Transport Equation



A

conservative,

positive

definite,

monotonic

scheme

based


on

flux

form

of

transport

equations





Evaluation



Upstream

scheme


Centered

scheme





Test

Cases



solid
-
body

rotation




deformational

flow


Model Equations

Integration on the volume of a

polyhedron with
N

faces leads

to the following formulation:

m

total air mass inside polygon

Fn

tracer inward fluxes for face
n

U
n

mass and inward fluxes for face
n

q

mean tracer mixing ratio in the same volume






)
5
(
6
1
n
n
U
t
m





)
5
(
6
1
n
n
F
t
qm
Upstream & Centered Scheme

Centered:

n
ˆ
t
ˆ
V

*
1
q
5
q
1
q
2
q
3
q
4
q
0
q
6
q
)
(
5
.
0
1
0
*
1
q
q
q


Upstream:

0
*
1
q
q

0
,
*
1

U
,
1
*
1
q
q

otherwise
Simulation Result

INITIAL

CENTERED

UPSTREAM

NEW
SCHEME

Stepwise Procedure

Step 1 : Compute gradient in each each triangle

Step 2 : Estimate the location of edge centre w.r.t.


cell center (0,0)

Step 4 : Estimate




Step 3 : Compute gradient at cell center





0
.
1
)
(
)
(
0
.
0
0
*
*
*
*
*
*
*
*
0
0
*
*
*
*










q
q
q
y
q
x
if
q
y
q
x
q
q
q
q
q
y
q
x
if
i
i
y
i
i
x
i
i
y
i
i
x
i
i
i
i
y
i
i
x
i

Step 5: Estimate the Slope limiter

2
*
2
*
2
*
i
i
i
y
x
x

2
*
2
*
2
*
i
i
i
y
x
y

2
*
2
*
*
*
i
i
i
i
y
x
y
x

OPEN PROBLEM :

Design of van Leer type scheme on
polygonal meshes

Stepwise Procedure

Step 6: Recomputed Gradient value recursively for all the polygonal edges


Step 7: Compute q estimate at edge

)
)(
1
(
*
2
*
2
*
*
*
*
2
*
2
*
2
*
*
*
i
y
i
i
i
i
i
x
i
i
i
i
x
i
x
q
y
x
y
x
q
y
x
x
q
q







)
)(
1
(
*
2
*
2
*
*
*
*
2
*
2
*
2
*
*
*
i
y
i
i
i
i
i
x
i
i
i
i
y
i
y
q
y
x
y
x
q
y
x
y
q
q







*
*
*
*
0
*
i
y
i
i
x
i
i
q
y
q
x
q
q



Note : Better estimate is the q value at location

2
ˆ
.
dt
n
V
i
This provides very promising monotonic, positive
and conservative scheme for polygonal meshes

Test Case 1 :
Solid
-
Body Rotation


)
sin
cos
sin
cos
(cos
0







u
u


sin
sin
0
u
v











R
r
if
R
r
if
R
r
q
o
0
,
/
cos(
1
)
,
(




r



great circle distance from centre of cosine bell

)
0
,
2
/
3
(
)
,
(
0




o
center of the cosine bell

64
/
7


R
bell radius

The

numerical

experiments

done

at

resolutions
:

level
-
32

(
np

=

10242
),

level
-
48

(
np

=

23042
)

and

level
-
64

(
np

=

40962
),

where

np

is

the

total

number

of

gridpoints

on

the

globe
.

The

time
-
step

is

so

chosen

that

a

complete

revolution

around

the

globe

takes

288

iterations
.


Comparison between different schemes

Global error norms at different resolutions

Simulation Results

Test Case 2:
Deformational Flow


Idealized Cyclogenesis
: Two steady vortices are located in the global
flow field. The flow field is obtained by transforming the coordinate system

such that poles coincide with the vortex centers in rotated system
with a constant angular velocity

)
,
(


)
,
(




)
(






The

normalized

tangential

velocity

of

the

vortex

field

is















dt
d
0


dt
d

)
tanh(
)
(
sec
2
3
3
2





h
V
t



cos















)
sin(
tanh
1
)
,
,
(
t
t
F






Exact

Solution










The

model

is

integrated

on

level
-
32

grid

for

6

time

units

with

different

time

steps

:


16
/
6


t
64
/
6


t
128
/
6


t

Results

Analytical Field

Absolute Error (np = 10242)


Simulated Field

Shallow Water Equations

(two
-
dimensional case)



Involves of differential operators such as:

Gradient, divergence, curl, Laplacian and Jacobian



Grid Geometry is triangular/hexagonal and not

rectangular



Mesh boundaries do not coincide with

meteorological zonal / meridional boundaries



Approximating Spherical hexagons by plane can

cause serious shape distortions at the boundary

Model Equations

.(hV)
-
t
h




)]
h
g(h
[K
V

k
ˆ
t
V
s









Z
h

height of free surface

h
s

surface height

K

kinetic energy

Z

absolute vorticity

Field Distortion : Laplacian

Hât Function on Hexagonal Element

Weak Formulation

An integral form that implicitly contains the differential
equation is called the
weak form


Let

ψ
be the the stream function then the vorticity is:

ψ
2




Choosing the shape function on a hexagonal cell as

S
kl
k
l
S


)
(
1
6



i
k
S
s. t.

We get :

Ω
d
ψ
S
d
S
Ω
Ω





2


Where is hexagonal cell

Weak Formulation Continued …

Now applying the green theorem and integrating by part:











d
S
Sd


Discrete form :








i
i
i
i
T
Area
T
Area
S
3
)
(
)
(
)
(



S=0


S=1

T
1

T
2

T
3

T
4

T
5

T
6


S=0


S=0

S=0


S=0

where gradients on approximated on
linear triangular elements


S=0

(weak form)

Weak Formulation Continued …

k
T
Area
Area
k
k









)]
(
[
2
1
)
(
k
T
Area
q
p
Area
q
p
k
k







)]
(
).
,
[(
2
1
)
(
)
,
.(
k
T
Area
q
p
Area
q
p
k
k
k









)]
(
).
,
[(
2
1
)
(
)
,
(
.
ˆ
GRADIENT :

DIVERGENCE :

CURL:

Analytical vs Numerical : Laplacian

Global Error Norms

)
cos(
)
(
cos
)
,
(
4
2




n
m
a
f
mn

x
f


y
f


f
2

Resolution


M= 1, N= 1


Level
-
24

2.19E
-
03

3.18E
-
03

0.012498

Level
-
48

6.69E
-
04

0.012498

6.52E
-
03

Level
-
96

2.15E
-
04

9.18E
-
04

4.30E
-
03

M= 1, N= 3


Level
-
24

4.67E
-
03

6.38E
-
03

0.014033

Level
-
48

1.32E
-
03

1.85E
-
03

7.57E
-
03

Level
-
96

1.12E
-
03

1.49E
-
03

5.53E
-
03

M= 3, N= 1

Level
-
24

0.012744

0.022382

0.063417

Level
-
48

3.68E
-
03

5.76E
-
03

0.017198

Level
-
96

3.96E
-
04

5.70E
-
04

5.04E
-
03

M= 3, N= 3

Level
-
24

0.021552

0.028671

0.060789

Level
-
48

5.68E
-
03

7.33E
-
03

0.017205

Level
-
96

1.54E
-
03

1.88E
-
03

6.53E
-
03

Time Integration

Matsuno leapfrog Scheme

L

L

L

L

L

L

L

L

L

L

L

L

M

M

M

The

Matsuno

scheme

(M)

is






(Euler

Forward)

















(Backward)






n
n
*
1
n
A
f
t
Δ
A
A





1)*
(n
n
1
A
f
t
Δ
A
A




n
The

Leap

Frog

scheme

(L)

is














n
1
n
1
n
A
f
t
2
Δ
A
A




Test Case 5 : Height field

Global Error Norm

Day

L1 norm

L2 norm

1

5.16E
-
04

5.50E
-
04

2

5.20E
-
04

5.58E
-
04

3

5.48E
-
04

6.29E
-
04

4

5.61E
-
04

6.59E
-
04

5

5.72E
-
04

6.79E
-
04

6

5.98E
-
04

7.05E
-
04

7

6.20E
-
04

7.61E
-
04

8

6.55E
-
04

8.13E
-
04

9

6.97E
-
04

8.83E
-
04

10

7.32E
-
04

9.32E
-
04

11

7.77E
-
04

9.80E
-
04

12

8.26E
-
04

1.08E
-
03

13

9.59E
-
04

1.26E
-
03

14

1.00E
-
03

1.37E
-
03

Test Case 6 : Height and Error field

Test case 6 : Rossby Haurwitz Wave

Day 1

Higher resolution ~ 1 degrees

Test case 6 : Rossby Haurwitz Wave

Day 7

Higher resolution ~ 1 degrees

Test case 6 : Rossby Hurwitz Wave

Day 14

Higher resolution ~ 1 degrees

Global Error Norm

Day

L1 norm

L2 norm

1

9.79E
-
04

1.08E
-
03

2

1.13E
-
03

1.31E
-
03

3

1.41E
-
03

1.64E
-
03

4

1.89E
-
03

2.26E
-
03

5

2.34E
-
03

2.81E
-
03

6

2.79E
-
03

3.37E
-
03

7

3.46E
-
03

4.37E
-
03

8

4.12E
-
03

5.43E
-
03

9

4.55E
-
03

6.07E
-
03

10

5.45E
-
03

6.53E
-
03

11

6.52E
-
03

7.06E
-
03

12

7.52E
-
03

6.96E
-
03

13

8.76E
-
03

7.56E
-
03

14

9.66E
-
03

7.66E
-
03

Baroclinic Model

( three
-
dimensional case )



Description of discrete dynamics of more comprehensive three

dimensional General Circulation Model of LMD, France



Development of new dynamics on Icosadedral
-
hexagonal grid



Testing LMDZ dynamics and new model dynamics w. r. t.

most recent and very popular Jablonskwi
-
Williamson (JW)

baroclinic wave test case

Hydrostatic Governing Equations

0
p
W
Y
V
X
U
p
b
t
p
~
s




















u
~
k
S
p
u
~
m
W
X
π
θ
X
B
ZV
t
u
~













v
~
k
S
p
v
~
m
W
Y
π
θ
Y
B
ZU
t
v
~





















k
θ
k
m
S
θW
p
θV
Y
θU
X
m
t





















k
Q
k
m
S
qW
p
qV
Y
qU
X
q
m
t












0
θ
π
Φ




Continuity Equation


Momentum Equations

Thermodynamic
Equation

Moisture conservation

Hydrostatic Equation

Governing Equations in terrain following pressure coordinate
in hydrostatic limits:

Horizontal Distribution

Vertical Distribution

P
s

P

P

P

P

P = 0


max


max
-
1


2


1

surface

W= 0

W

W

W

W

W= 0


,
,
,
,
,
~
,
~
B
Z
Q
H
v
u

,
,
,
,
,
~
,
~
B
Z
Q
H
v
u

,
,
,
,
,
~
,
~
B
Z
Q
H
v
u

,
,
,
,
,
~
,
~
B
Z
Q
H
v
u
Top of the atmosphere

J
-
W Test Case


LEVEL 48

LMDZ Model versus ICOSA MODEL


ICOSA


LMDZ

LMDZ Model versus ICOSA MODEL


ICOSA


192 x 144


96 x 72


LEVEL 24


LMDZ

LMDZ Model versus ICOSA MODEL


ICOSA

Day 9: Models Comparison

Relative Vorticity

Level 24

Level 48


96 x 72


192 x 144

7

7

9

9

DAY 7 & DAY 9

Relative Vorticity

DAY 7

DAY 9

Inclusion of Dissipation

Day 6

Day 8

Day 10

Beyond 10 days

Day 15

Beyond 10 days

Day 20

Beyond 10 days

Day 25

Beyond 10 days

Day 30

Steady State Case for 30 days period

INCLUSION OF DISSIPATION

Future Research Directions



Incorporation of LMDZ physics to make it a
comprehensive forecasting system on Icosahedral
-
hexagonal grid geometry.



Long term integration with this dynamical core for
thorough evaluation in realistic situations.



Exploring the possibilities of developing numerical
methods on staggered grids.



Inclusion of Adaptive Meshing

THANKS

Meteorological Primitive Equations


Applicable to wide scale of motions; > 1hour, >100km


q
p
S
dq/dt
RT/p
p
Φ/
p
ω/
V
Q/c
κTω/p
dT/dt
F
Φ
V
k
f
/dt
V
d


















0
0
ˆ



(horizontal momentum)

(thermodynamic equation)

(mass continuity)

(hydrostatic approximation)

(water vapor mass continuity)