[615] | 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 | |
---|