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 ) |
---|