source: codes/icosagcm/devel/Python/dynamico/DCMIP.py @ 615

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

devel/unstructured : import Python layer from HEAT
File size: 7.0 KB
Line 
1import matplotlib.pyplot as plt
2import time as time
3import math as math
4import numpy as np
5from numpy import cos, sin, exp, log
6from math import pi
7import dyn as dyn
8
9Linf = lambda x : abs(x).max()
10
11def run_with_scalars(model,dyn_scheme,transport,replace, plotter, 
12                     courant,flow0,scalars, T1, N1, N2):
13    metric, BC = model.metric, model.BC
14    mik, gas0, ujk, Phi_il, Wil = model.diagnose(*flow0)
15    ztop = Phi_il.max()/metric.g0
16    Phi_ik = metric.mk_int(Phi_il)
17    vol, jac = metric.volume(Phi_il[:,-1])
18    internal0 = (mik*(gas0.e+Phi_ik)).sum() + BC.ptop*vol.sum()
19
20    dtmax, dt_hydro = model.max_time_step(gas0,mik)
21    dt=courant*dtmax # RK4 stability limit is courant=2.sqrt(2)
22   
23    N0 = 2*int(math.ceil(.5*float(T1)/dt))
24    dt=T1/N0 # adjust time step to fraction of T1
25    step = dyn_scheme(dt)
26    print 'dt (s)', step.dt, 'Horizontal acoustic Courant number', dt/dt_hydro
27   
28    Sik = mik*gas0.s
29    flow0 = mik, Sik, ujk, Phi_il, Wil, 0.,0.
30    elapsed=time.time()
31    flow=step.advance(flow0,10)
32    elapsed = (time.time()-elapsed)/10
33    print "Elapsed time per time step (ms)", elapsed*1e3
34    print "Simulated time / elapsed time  (s/s)", step.dt/elapsed
35       
36    plotter(metric, (mik, Sik, ujk, Phi_il, Wil), scalars, 0.)
37
38    out=[flow0]
39    kinetic=[]
40    internal=[]
41    for i in range(N2):
42        for j in range(N1):
43            for k in range(N0/2):
44                mik_old = mik
45                mik, Sik, ujk, Phi_il, Wil, Fjk, Fil = step.next((mik, Sik, ujk, Phi_il, Wil, 0.,0.))
46                mik, Sik, ujk, Phi_il, Wil, Fjk, Fil = step.next((mik, Sik, ujk, Phi_il, Wil, Fjk, Fil))
47                if type(Fil) is type(Sik):# fix Lagrangian coordinate
48                    # transport scalars using time-integrated mass fluxes Fjk, Fil
49                    mik_new, scalars = transport.next(Fjk,Fil,mik_old,scalars)
50                    Sik = replace(Sik,scalars[0]) # replace (or not) entropy by FV-transported scalar
51                    # print Linf(mik-mik_new)/Linf(mik)
52            mik, gas_ik, ujk, Phi_il, Wil = model.diagnose(mik, Sik, ujk, Phi_il, Wil)
53            Phi_ik = metric.mk_int(Phi_il)
54            vol, jac = metric.volume(Phi_il[:,-1])
55            kinetic.append(model.kinetic(mik, ujk, Wil, Phi_il))
56            internal.append((mik*(gas_ik.e+Phi_ik)).sum() + BC.ptop*vol.sum())
57            print '\r', i*N1+j+1, '/', N1*N2,
58        plotter(metric, (mik, Sik, ujk, Phi_il, Wil), scalars, (i+1)*T1*N1)
59       
60    kinetic=np.array(kinetic)
61    internal=np.array(internal)
62    dE=kinetic+internal-internal0
63    plt.figure(figsize=(12,3))
64    plt.plot(T1*np.arange(N1*N2),dE/kinetic.max())
65    plt.xlabel('time')
66    plt.ylabel('(E(t)-E(0))/Kmax')
67   
68    plt.figure(figsize=(12,3))
69    plt.plot(T1*np.arange(N1*N2),kinetic/kinetic.max())
70    plt.xlabel('time')
71    plt.ylabel('K(t)/Kmax')
72    plt.show()
73
74    return model.diagnose(mik, Sik, ujk, Phi_il, Wil)
75
76def DCMIP21(Lx,Nx,llm, ztop=3e4, h0=250., 
77            dd=5000., xi=4000., Ueq=20.,
78            deform=lambda x,eta:0 ) :
79    """ DCMIP2.1 experiment at equator (phi=0)"""
80    Cpd, Rd, pref, Tref, = 1004.5, 287., 1e5, 273.15
81    g, p0, peq, Teq = 9.81, 1e5, 1e5, 300.
82
83    p=lambda x,Phi : peq*np.exp(-Phi/(Rd*Teq))
84    ptop = p(0,g*ztop)
85    print 'ptop (Pa) =', ptop
86
87    xmesh  = dyn.Periodic_FD(Nx)
88    zmesh  = dyn.Lorenz(llm)
89    mesh   = dyn.Slice(xmesh,zmesh)
90    metric = dyn.Shallow(mesh, Lx/Nx,g)
91    thermo = dyn.Ideal_perfect(Cpd, Rd, pref, Tref)
92
93    T=lambda x,Phi : 0*x+Teq
94    gauss=lambda x: np.exp(-x*x)
95    cos2=lambda x: cos(pi*x)**2
96    Phis = lambda x:g*h0*gauss(x/dd)*cos2(x/xi)
97    Phi = lambda x,eta : g*ztop*eta + (1.-eta)*Phis(x)
98    Phi_deform = lambda x,eta : Phi(x,eta+deform(x/Lx,eta))
99    mik, gas, Phi_il = metric.inicond_NH(thermo, Phi_deform,p,T)
100    ujk = 0*mik + Ueq*(Lx/Nx)
101   
102    print 'ztop (m) = ', Phi_il[0,-1]/g, ztop
103    print 'ptop (Pa) = ', gas.p[0,-1], ptop
104    Phi_bot = Phi_il[:,0:1]
105    gas_bot=thermo.set_pT(p(metric.xil[:,0:1], Phi_bot), Teq)
106    BC = dyn.BoundCond(ptop, Phi_bot, gas_bot.p, 1./gas_bot.v)
107    return metric, thermo, BC, mik, gas, Phi_il, ujk
108
109def DCMIP31(Lx,Nx,llm,deform=lambda x,eta:0, dTheta=1., d2=25e6) :
110    """ DCMIP3.1 initial condition with Phi(x,eta)=g.ztop*(eta + deform(eta,x/L)) """
111    # DCMIP 3.1 at equator (phi=0)
112    Cpd, Rd, pref, Tref, = 1004.5, 287., 1e5, 273.15
113    g, p0, ptop = 9.81, 1e5, 27391.9
114    peq, N2, Teq, ztop = 1e5, 1e-4, 300., 1e4
115
116    G=g**2/(Cpd*N2)
117    inv_G=1./G
118    inv_Teq=1./Teq
119    G_Teq = G/Teq
120
121    p=lambda x,Phi : peq*((1-G_Teq*(1-np.exp(-N2*Phi/g**2)))**(Cpd/Rd))
122    ptop = p(0,g*ztop)
123
124    xmesh  = dyn.Periodic_FD(Nx)
125    zmesh  = dyn.Lorenz(llm)
126    mesh   = dyn.Slice(xmesh,zmesh)
127    metric = dyn.Shallow(mesh, Lx/Nx,g)
128    thermo = dyn.Ideal_perfect(Cpd, Rd, pref, Tref)
129
130    ps = lambda xi: peq+0*xi
131    fac = lambda p : (p/peq)**(Rd/Cpd)
132    T = lambda x,p : fac(p)/(inv_G*(fac(p)-1)+inv_Teq)
133#    T = lambda x,p : 0.*p+Teq
134
135    Phi = lambda x,eta : g*ztop*(eta+deform(x/Lx,eta))
136    m0ik, gas, Phi0_il = metric.inicond_NH(thermo, Phi,p,T)
137    u0_jk = 0*m0ik
138#    u0_jk = 0*m0ik + 100.*Lx/Nx
139   
140    Tb  = gas.T # unperturbed temperature
141
142    print 'ztop (m) = ', Phi0_il[0,-1]/g, ztop
143    print 'ptop (Pa) = ', gas.p[0,-1], ptop
144    Phi_bot = Phi0_il[:,0:1]
145    pbot = p(metric.xil[:,0:1], Phi_bot)
146    gas_bot=thermo.set_pT(pbot, T(metric.xil[:,0:1],pbot))
147    BC = dyn.BoundCond(ptop, Phi_bot, gas_bot.p, 1./gas_bot.v)
148    print BC.rho_bot.shape
149   
150    # now add localized temperature perturbation
151    xik = metric.xik
152    pik = gas.p
153    zik = zmesh.mk_int(Phi0_il)/g
154    shape = d2/(d2+xik*xik)
155    shape = shape*np.sin(math.pi*zik/ztop)
156    shape = shape*((pik/p0)**(Rd/Cpd))
157    gas0 = thermo.set_pT(pik, Tb + dTheta*shape)
158
159    return metric, thermo, BC, m0ik, gas0, Phi0_il, u0_jk
160
161def thermal_bubble(Lx,Nx,llm, ztop=1000., zc=200., rc=250, thetac=0.5, x0=0.0):
162    """Hi"""
163    Cpd, Rd, g, p0,theta0, T0 = 1004.5, 287.,9.81, 1e5, 300., 300.
164
165    Phis = lambda x: 0.0
166    Phi = lambda x,eta : g*ztop*eta + (1.-eta)*Phis(x)
167    p=lambda x, Phi : p0*np.exp(-Phi/(Rd*T0))
168    ptop = p(0,g*ztop)
169    print 'ptop (Pa) =', ptop
170    xmesh  = dyn.Periodic_FD(Nx)
171    zmesh  = dyn.Lorenz(llm)
172    mesh   = dyn.Slice(xmesh,zmesh)
173    metric = dyn.Shallow(mesh, Lx/Nx,g)
174    thermo = dyn.Ideal_perfect(Cpd, Rd, p0, T0)
175   
176    zz = lambda p: -(Rd*T0*np.log(p/p0))/g
177    rr =   lambda x, p: np.sqrt((x-x0)**2 + (zz(p)-zc)**2)
178    sa = lambda x,p: rr(x,p) < rc
179    deform = lambda x,p: (0.5*thetac*(1+cos(pi*rr(x,p)/rc)))*sa(x,p)
180    temp =  lambda x, p: theta0*(p/p0)**(Rd/Cpd)
181    T = lambda x,p: deform(x,p) + temp(x,p)
182       
183    #T = lambda x, Phi: T0
184    mik, gas, Phi_il = metric.inicond_NH(thermo, Phi,p,T)
185    ujk = 0*mik
186       
187    print 'ztop (m) = ', Phi_il[0,-1]/g, ztop
188    print 'ptop (Pa) = ', gas.p[0,-1], ptop
189    Phi_bot = Phi_il[:,0:1]
190    gas_bot=thermo.set_pT(p(metric.xil[:,0:1], Phi_bot), T0) #### T ot T0
191    BC = dyn.BoundCond(ptop, Phi_bot, gas_bot.p, 1./gas_bot.v)
192    return metric, thermo, BC, mik, gas, Phi_il, ujk
193
Note: See TracBrowser for help on using the repository browser.