source: codes/icosagcm/devel/Python/test/py/DCMIP2008c5.py @ 680

Last change on this file since 680 was 680, checked in by dubos, 6 years ago

devel/unstructured : moved mesh partitioning code into meshes.py
File size: 8.0 KB
Line 
1from dynamico import unstructured as unst
2from dynamico import dyn
3from dynamico import time_step
4from dynamico import DCMIP
5from dynamico import meshes
6import math as math
7import matplotlib.pyplot as plt
8import numpy as np
9import time
10
11#------------------------ initial condition -------------------------
12
13# Parameters for the simulation
14g              = 9.80616   # gravitational acceleration in meters per second squared
15omega          = 7.29211e-5  # Earth's angular velocity in radians per second
16f0             = 2.0*omega   # Coriolis parameter
17u_0            = 20.0        # velocity in meters per second
18T_0            = 288.0       # temperature in Kelvin
19d2             = 1.5e6**2    # square of half width of Gaussian mountain profile in meters
20h_0            = 2.0e3       # mountain height in meters
21lon_c          = np.pi/2.0   # mountain peak longitudinal location in radians
22lat_c          = np.pi/6.0   # mountain peak latitudinal location in radians
23radius         = 6.371229e6  # mean radius of the Earth in meters
24ref_press      = 100145.6    # reference pressure (mean surface pressure) in Pascals
25ref_surf_press = 930.0e2     # South Pole surface pressure in Pascals
26Rd             = 287.04      # ideal gas constant for dry air in joules per kilogram Kelvin
27Cpd            = 1004.64     # specific heat at constant pressure in joules per kilogram Kelvin
28kappa          = Rd/Cpd      # kappa=R_d/c_p
29N_freq         = np.sqrt(g**2/(Cpd*T_0)) # Brunt-Vaisala buoyancy frequency
30N2, g2kappa = N_freq**2, g*g*kappa
31
32def DCMIP2008c5(grid,llm):
33    def Phis(lon,lat): 
34        rgrc = radius*np.arccos(np.sin(lat_c)*np.sin(lat)+np.cos(lat_c)*np.cos(lat)*np.cos(lon-lon_c))
35        return g*h_0*np.exp(-rgrc**2/d2)
36    def ps(lon, lat, Phis): 
37        return ref_surf_press * np.exp(
38            - radius*N2*u_0/(2.0*g2kappa)*(u_0/radius+f0)*(np.sin(lat)**2-1.)
39            - N2/(g2kappa)*Phis )
40    def ulon(lat): return u_0*np.cos(lat)
41    def ulat(lat): return 0.*lat
42    def f(lon,lat): return f0*np.sin(lat) # Coriolis parameter
43
44    nqdyn, preff, Treff = 1, 1e5, 300.
45    thermo = dyn.Ideal_perfect(Cpd, Rd, preff, Treff)
46    meshfile = meshes.MPAS_Format('grids/x1.%d.grid.nc'%grid)
47    mesh = meshes.Unstructured_Mesh(meshfile, llm, nqdyn, radius, f)
48    lev,levp1 = np.arange(llm),np.arange(llm+1)
49    lon_i, lat_i, lon_e, lat_e =  mesh.lon_i, mesh.lat_i, mesh.lon_e, mesh.lat_e
50    lat_ik,k_i = np.meshgrid(mesh.lat_i,lev, indexing='ij')
51    lon_ik,k_i = np.meshgrid(mesh.lon_i,lev, indexing='ij')
52    lat_il,l_i = np.meshgrid(mesh.lat_i,levp1, indexing='ij')           
53    lon_il,l_i = np.meshgrid(mesh.lon_i,levp1, indexing='ij')
54    lat_ek,k_e = np.meshgrid(mesh.lat_e,lev, indexing='ij')
55
56    Phis_i = Phis(lon_i, lat_i)
57    ps_i = ps(lon_i, lat_i, Phis_i)
58
59    if llm==18:
60        ap_l=[0.00251499, 0.00710361, 0.01904260, 0.04607560, 0.08181860,
61              0.07869805, 0.07463175, 0.06955308, 0.06339061, 0.05621774, 0.04815296,
62              0.03949230, 0.03058456, 0.02193336, 0.01403670, 0.007458598, 0.002646866,
63              0.0, 0.0 ]
64        bp_l=[0.0, 0.0, 0.0, 0.0, 0.0, 0.03756984, 0.08652625, 0.1476709, 0.221864,
65              0.308222, 0.4053179, 0.509588, 0.6168328, 0.7209891, 0.816061, 0.8952581, 
66              0.953189, 0.985056, 1.0 ]
67    if llm==26:
68        ap_l=[0.002194067, 0.004895209, 0.009882418, 0.01805201, 0.02983724, 0.04462334, 0.06160587, 
69               0.07851243, 0.07731271, 0.07590131, 0.07424086, 0.07228744, 0.06998933, 0.06728574, 0.06410509, 
70               0.06036322, 0.05596111, 0.05078225, 0.04468960, 0.03752191, 0.02908949, 0.02084739, 0.01334443, 
71               0.00708499, 0.00252136, 0.0, 0.0 ]
72        bp_l=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.01505309, 0.03276228, 0.05359622, 
73                0.07810627, 0.1069411, 0.1408637, 0.1807720, 0.2277220, 0.2829562, 0.3479364, 0.4243822, 
74                0.5143168, 0.6201202, 0.7235355, 0.8176768, 0.8962153, 0.9534761, 0.9851122, 1.0 ]
75
76    ap_l, bp_l = ref_press*np.asarray(ap_l[-1::-1]), bp_l[-1::-1] # revserse indices
77    ptop = ap_l[-1] # pressure BC
78   
79    print ptop, ap_l, bp_l
80    ps_il,ap_il = np.meshgrid(ps_i,ap_l, indexing='ij')
81    ps_il,bp_il = np.meshgrid(ps_i,bp_l, indexing='ij')
82    hybrid_coefs = meshes.mass_coefs_from_pressure_coefs(g, ap_il, bp_il)
83    pb_il = ap_il + bp_il*ps_il
84    mass_ik, pb_ik = mesh.field_mass(), mesh.field_mass()
85    for l in range(llm):
86        pb_ik[:,l]=.5*(pb_il[:,l]+pb_il[:,l+1])
87        mass_ik[:,l]=(pb_il[:,l]-pb_il[:,l+1])/g
88    Tb_ik = T_0 + 0.*pb_ik
89    gas = thermo.set_pT(pb_ik,Tb_ik)
90    Sik  = gas.s*mass_ik
91# start at hydrostatic geopotential
92    Phi_il = mesh.field_w()
93    Phi_il[:,0]=Phis_i
94    for l in range(llm):
95        Phi_il[:,l+1] = Phi_il[:,l] + g*mass_ik[:,l]*gas.v[:,l]
96
97    ujk, Wil = mesh.ucov3D(ulon(lat_ek),ulat(lat_ek)), mesh.field_w()
98   
99    print 'ztop (m) = ', Phi_il[0,-1]/g
100    print 'ptop (Pa) = ', gas.p[0,-1], ptop
101    dx=mesh.de.min()
102    params=dyn.Struct()
103    params.g, params.ptop = g, ptop
104    params.dx, params.dx_g0 = dx, dx/g
105    params.pbot, params.rho_bot = 1e5+0.*mesh.lat_i, 1e6+0.*mesh.lat_i
106    return thermo, mesh, hybrid_coefs, params, (mass_ik,Sik,ujk,Phi_il,Wil), gas
107
108
109#------------------------ main program -------------------------
110
111grid, llm = 40962, 26
112T, Nslice, courant = 14400, 24, 3.0
113#caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_lag
114caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass
115thermo, mesh, hybrid_coefs, params, flow0, gas0 = DCMIP2008c5(grid,llm)
116llm, dx = mesh.llm, params.dx
117print 'llm, nb_hex, dx =', llm, mesh.Ai.size, dx
118if caldyn_eta == unst.eta_lag:
119    print 'Lagrangian coordinate.'
120else:
121    print 'Mass-based coordinate.'
122    unst.ker.dynamico_init_hybrid(*hybrid_coefs)
123   
124dt = courant*.5*dx/np.sqrt(gas0.c2.max())
125
126nt = int(math.ceil(T/dt))
127dt = T/nt
128print 'Time step : %g s' % dt
129
130mesh.plot_e(mesh.le/mesh.de) ; plt.title('le/de')
131plt.savefig('fig_DCMIP2008c5/le_de.png') ;plt.close()
132
133mesh.plot_i(mesh.Ai) ; plt.title('Ai')
134plt.savefig('fig_DCMIP2008c5/Ai.png'); plt.close()
135
136
137if False: # time stepping in Python
138    caldyn = unst.Caldyn_HPE(caldyn_thermo,caldyn_eta, mesh,thermo,params,params.g)
139    scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt, a32=0.7)
140    def next_flow(m,S,u): 
141        m,S,u = scheme.advance((m,S,u),nt)
142        return m,S,u, caldyn.geopot
143else: # time stepping in Fortran
144    scheme = time_step.ARK2(None, dt, a32=0.7)
145    caldyn_step = unst.caldyn_step_HPE(mesh,scheme,nt, caldyn_thermo,caldyn_eta, thermo,params,params.g)
146    def next_flow(m,S,u):
147        caldyn_step.mass[:,:], caldyn_step.theta_rhodz[:,:], caldyn_step.u[:,:] = m,S,u
148        caldyn_step.next()
149        return (caldyn_step.mass.copy(), caldyn_step.theta_rhodz.copy(), 
150                caldyn_step.u.copy(), caldyn_step.geopot.copy())
151
152def plots(it):
153    s=S/m
154    for l in range(llm):
155        z[:,l]=.5*(Phi[:,l+1]+Phi[:,l])/params.g
156        vol[:,l]=(Phi[:,l+1]-Phi[:,l])/params.g/m[:,l] # specific volume
157    gas = thermo.set_vs(vol, s)
158    s=.5*(s+abs(s))
159    t = (it*T)/3600
160    print( 'ptop, model top (m) :', unst.getvar('ptop'), Phi.max()/unst.getvar('g') )
161    mesh.plot_i(gas.T[:,llm/2])
162    plt.title('T at t=%dh'%(t))
163    plt.savefig('fig_DCMIP2008c5/T%02d.png'%it)
164    plt.close()
165
166    mesh.plot_i(m[:,llm/2])
167    plt.title('mass at t=%dh'%(t))
168    plt.savefig('fig_DCMIP2008c5/m%02d.png'%it)
169    plt.close()
170
171    mesh.plot_i(Phi[:,0])
172    plt.title('Surface geopotential at t=%dh'%(t))
173    plt.savefig('fig_DCMIP2008c5/Phis%02d.png'%it)
174    plt.close()
175
176z, vol = mesh.field_mass(), mesh.field_mass()
177m,S,u,Phi,W = flow0
178caldyn_step.geopot[:,0]=Phi[:,0]
179plots(0)
180
181for it in range(Nslice):
182    unst.setvar('debug_hevi_solver',False)
183    time1, elapsed1 =time.time(), unst.getvar('elapsed')
184    m,S,u,Phi = next_flow(m,S,u)
185    time2, elapsed2 = time.time(), unst.getvar('elapsed')
186    factor = 1000./nt
187    print 'ms per full time step : ', factor*(time2-time1), factor*(elapsed2-elapsed1)
188    factor = 1e9/(4*nt*m.size)
189    print 'nanosec per gridpoint per call to caldyn_hevi : ', factor*(time2-time1), factor*(elapsed2-elapsed1)
190    plots(it+1)
Note: See TracBrowser for help on using the repository browser.