1 | import matplotlib.pyplot as plt |
---|
2 | import time as time |
---|
3 | import math as math |
---|
4 | import numpy as np |
---|
5 | from numpy import cos, sin, exp, log |
---|
6 | from math import pi |
---|
7 | import dyn as dyn |
---|
8 | |
---|
9 | Linf = lambda x : abs(x).max() |
---|
10 | |
---|
11 | def 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 | |
---|
76 | def 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 | |
---|
109 | def 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 | |
---|
161 | def 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 | |
---|