[617] | 1 | import math as math |
---|
| 2 | import matplotlib.pyplot as plt |
---|
| 3 | import numpy as np |
---|
| 4 | import dynamico.dyn as dyn |
---|
| 5 | import dynamico.time_step as time_step |
---|
| 6 | import dynamico.DCMIP as DCMIP |
---|
| 7 | import dynamico.advect as advect |
---|
| 8 | |
---|
| 9 | dx, Nx, llm, tau_relax = 500., 400, 60, 25. |
---|
| 10 | #dx, Nx, llm, tau_relax = 250., 800, 120, 25. |
---|
| 11 | T1, N1, N2 = 30., 6, 20 # check energy every 60s, plot every 10min, run 10h |
---|
| 12 | #dx, Nx, llm, tau_relax = 20000., 200, 200, 1e3 |
---|
| 13 | #T1, N1, N2 = 1800., 4, 6 # check energy every 60s, plot every 10min, run 2h |
---|
| 14 | |
---|
| 15 | Lx=Nx*dx |
---|
| 16 | metric, thermo, BC, m0ik, gas0, Phi0_il, u0_jk = DCMIP.DCMIP21( |
---|
| 17 | Lx, Nx, llm, h0=250., ztop=3e4, dd=4000, xi=3200) |
---|
| 18 | # Lx, Nx, llm, h0=250., ztop=3e4, dd=10*dx, xi=8*dx) |
---|
| 19 | Sik = m0ik*gas0.s |
---|
| 20 | start =(m0ik, Sik, u0_jk, Phi0_il, 0*Phi0_il) |
---|
| 21 | scalars = (Sik, m0ik*metric.mk_int(Phi0_il)) |
---|
| 22 | |
---|
| 23 | relax = lambda Phi : (1./tau_relax)*(1+np.tanh((Phi/9.81-23e3)/3e3)) |
---|
| 24 | |
---|
| 25 | def sponge(Phi_il, vjk, dujk_dt, tau): |
---|
| 26 | tau_rel = 1./relax(Phi_il[:,0:-1]) |
---|
| 27 | dujk = (vjk-u0_jk + tau*dujk_dt)/(tau+tau_rel) |
---|
| 28 | ujk = u0_jk + tau_rel*dujk |
---|
| 29 | dujk_dt = dujk_dt - dujk |
---|
| 30 | return ujk, dujk_dt |
---|
| 31 | |
---|
| 32 | BC.sponge=sponge # will be called by fast_bwd |
---|
| 33 | BC.rigid_top = False |
---|
| 34 | BC.rigid_bottom = False |
---|
| 35 | rho_bot = BC.rho_bot |
---|
| 36 | BC.rho_bot = 1e6*rho_bot |
---|
| 37 | |
---|
| 38 | def go(model,courant,scheme,start): |
---|
| 39 | def plotter(metric, flow, scalars, t): |
---|
| 40 | mik, gas_ik, ujk, Phi_il, Wil = model.diagnose(*flow) |
---|
| 41 | mil = metric.ml_ext(mik) |
---|
| 42 | wil = metric.g0*Wil/mil |
---|
| 43 | xil = metric.xil/1000 |
---|
| 44 | xjl = metric.xjl/1000 |
---|
| 45 | zil = Phi_il/metric.g0/1000 |
---|
| 46 | zjl = metric.mj(zil) |
---|
| 47 | plt.figure(figsize=(12,5)) |
---|
| 48 | # plt.figure(figsize=(6,5)) |
---|
| 49 | # plt.pcolormesh(xjl, zjl, gas_ik.T-gas0.T) |
---|
| 50 | # plt.contourf(xil, zil, wil, np.arange(-.2,.2,.02)) |
---|
| 51 | plt.contourf(xil, zil, wil, np.arange(-1.,1.,.1)) |
---|
| 52 | plt.plot(xil,zil[:,llm/10]) |
---|
| 53 | plt.plot(xil,zil[:,llm/5]) |
---|
| 54 | plt.plot(xil,zil[:,llm/4]) |
---|
| 55 | plt.plot(xil,zil[:,llm/2]) |
---|
| 56 | plt.xlabel('x (km)') |
---|
| 57 | plt.ylabel('z (km)') |
---|
| 58 | # plt.ylim((-1,30.)), plt.xlim((-dx/20.,dx/10.)) |
---|
| 59 | plt.ylim((0.,10.)), plt.xlim((-10.,10.)) |
---|
| 60 | plt.colorbar() |
---|
| 61 | plt.title("t = %f s" % (t)) |
---|
| 62 | plt.savefig("fig_slice_orographic/%04d"%t) |
---|
| 63 | plt.close() |
---|
| 64 | transport = advect.SplitTransport(model.metric, advect.HorizVanLeer, advect.VertVanLeer) |
---|
| 65 | replace = lambda fd,fv : fv |
---|
| 66 | DCMIP.run_with_scalars(model,scheme,transport,replace, plotter, courant, start,scalars, T1, N1, N2) |
---|
| 67 | |
---|
| 68 | model = dyn.Nonhydro(thermo,metric,BC) |
---|
| 69 | solver = dyn.MassBasedNH(model) |
---|
| 70 | #solver = dyn.LagrangianNH(model) |
---|
| 71 | |
---|
| 72 | def bwd_fast_slow(flow,tau): |
---|
| 73 | flow,fast,slow=solver.bwd_fast_slow(flow,tau) |
---|
| 74 | m,S,u,Phi,W,X,Y = flow |
---|
| 75 | dm,dS,du,dPhi,dW,dX,dY = slow |
---|
| 76 | dW = dW - relax(Phi)*W # damp vertical velocity |
---|
| 77 | slow = dm,dS,du,dPhi,dW,dX,dY |
---|
| 78 | return flow,fast,slow |
---|
| 79 | |
---|
| 80 | time_scheme = lambda dt : time_step.ARK2(bwd_fast_slow, dt, a32=0.7) |
---|
| 81 | go(model, 3.0, time_scheme, start ) |
---|