1 | from dynamico import unstructured as unst |
---|
2 | from dynamico import dyn |
---|
3 | from dynamico import time_step |
---|
4 | from dynamico import DCMIP |
---|
5 | from dynamico import meshes |
---|
6 | from dynamico.meshes import Cartesian_mesh as Mesh |
---|
7 | import math as math |
---|
8 | import matplotlib.pyplot as plt |
---|
9 | import numpy as np |
---|
10 | import time |
---|
11 | |
---|
12 | def thermal_bubble_3D(Lx,nx,Ly,ny,llm,ztop=1000., zc=350., |
---|
13 | rc=250, thetac=0.5, x0=0., y0=0.): |
---|
14 | Cpd, Rd, g, p0,theta0, T0 = 1004.5, 287.,9.81, 1e5, 300., 300. |
---|
15 | nqdyn = 1 |
---|
16 | |
---|
17 | Phi = lambda eta : g*ztop*eta |
---|
18 | p=lambda Phi : p0*np.exp(-Phi/(Rd*T0)) |
---|
19 | zz = lambda p: -(Rd*T0*np.log(p/p0))/g |
---|
20 | rr = lambda x,y,p: np.sqrt((x-x0)**2 + (y-y0)**2 + (zz(p)-zc)**2) |
---|
21 | sa = lambda x,y,p: rr(x,y,p) < rc |
---|
22 | deform = lambda x,y,p: (0.5*thetac*(1+np.cos(np.pi*rr(x,y,p)/rc)))*sa(x,y,p) |
---|
23 | temp = lambda p: theta0*(p/p0)**(Rd/Cpd) |
---|
24 | T = lambda x,y,p: deform(x,y,p) + temp(p) |
---|
25 | |
---|
26 | mesh = Mesh(nx,ny,llm,nqdyn,Lx,Ly,0.) |
---|
27 | thermo = dyn.Ideal_perfect(Cpd, Rd, p0, T0) |
---|
28 | |
---|
29 | Phi_il = Phi(mesh.llp1/float(llm)) |
---|
30 | Phi_ik = Phi((mesh.ll+.5)/llm) |
---|
31 | p_ik = p(Phi_ik) |
---|
32 | T_ik = T(mesh.xx, mesh.yy, p_ik) |
---|
33 | |
---|
34 | gas = thermo.set_pT(p_ik,T_ik) |
---|
35 | mass_ik = mesh.field_mass() |
---|
36 | for l in range(llm): |
---|
37 | mass_ik[:,:,l]=(Phi_il[:,:,l+1]-Phi_il[:,:,l])/(g*gas.v[:,:,l]) |
---|
38 | Sik, ujk, Wil = gas.s*mass_ik, mesh.field_u(), mesh.field_w() |
---|
39 | |
---|
40 | print 'ztop (m) = ', Phi_il[0,0,-1]/g, ztop |
---|
41 | ptop = p(g*ztop) |
---|
42 | print 'ptop (Pa) = ', gas.p[0,0,-1], ptop |
---|
43 | params=dyn.Struct() |
---|
44 | params.ptop=ptop |
---|
45 | params.dx=dx |
---|
46 | params.dx_g0=dx/g |
---|
47 | params.g = g |
---|
48 | pbot = p(Phi_il[:,:,0]) |
---|
49 | gas_bot = thermo.set_pT(pbot, temp(pbot)) |
---|
50 | params.pbot = gas_bot.p |
---|
51 | params.rho_bot = 1e6/gas_bot.v |
---|
52 | |
---|
53 | return thermo, mesh, params, (mass_ik,Sik,ujk,Phi_il,Wil), gas |
---|
54 | |
---|
55 | Lx, nx, llm, thetac, T, Nslice, courant = 2000., 100, 50, 30., 5., 10, 2.8 |
---|
56 | #Lx, nx, llm, thetac, T, Nslice, courant = 2000., 50, 25, 1., 50., 10, 2.8 |
---|
57 | #Lx, nx, llm, thetac, T, Nslice, courant = 3000., 75, 25, -30, 5., 10, 2.8 |
---|
58 | |
---|
59 | nqdyn, dx = 1, Lx/nx |
---|
60 | Ly,ny,dy = Lx,nx,dx |
---|
61 | thermo, mesh, params, flow0, gas0 = thermal_bubble_3D(Lx,nx,Ly,ny,llm, thetac=thetac) |
---|
62 | |
---|
63 | caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass |
---|
64 | caldyn = unst.Caldyn_NH(caldyn_thermo,caldyn_eta, mesh,params,thermo,params,params.g) |
---|
65 | |
---|
66 | # compute hybrid coefs from initial distribution of mass |
---|
67 | mass_bl,mass_dak,mass_dbk = meshes.compute_hybrid_coefs(flow0[0]) |
---|
68 | unst.ker.dynamico_init_hybrid(mass_bl,mass_dak,mass_dbk) |
---|
69 | print mass_dbk |
---|
70 | |
---|
71 | dz = flow0[3].max()/(params.g*llm) |
---|
72 | #dt = courant*.5/np.sqrt(gas0.c2.max()*(dx**-2+dy**-2+dz**-2)) |
---|
73 | dt = courant*.5/np.sqrt(gas0.c2.max()*(dx**-2+dy**-2)) |
---|
74 | nt = int(math.ceil(T/dt)) |
---|
75 | dt = T/nt |
---|
76 | print 'Time step : %g s' % dt |
---|
77 | |
---|
78 | scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt) |
---|
79 | |
---|
80 | flow=flow0 |
---|
81 | w=mesh.field_mass() |
---|
82 | z=mesh.field_mass() |
---|
83 | for it in range(Nslice): |
---|
84 | junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.) |
---|
85 | m,S,u,Phi,W = flow |
---|
86 | s=S/m ; # s=.5*(s+abs(s)) |
---|
87 | for l in range(llm): |
---|
88 | w[:,:,l]=.5*params.g*(W[:,:,l+1]+W[:,:,l])/m[:,:,l] |
---|
89 | z[:,:,l]=.5*(Phi[:,:,l+1]+Phi[:,:,l])/params.g |
---|
90 | |
---|
91 | print 'ptop, model top (m) :', unst.getvar('ptop'), Phi.max()/unst.getvar('g') |
---|
92 | jj=ny/2 |
---|
93 | xx,zz,ss,ww = mesh.xx[:,jj,:]/1000., z[:,jj,:]/1000., s[:,jj,:], w[:,jj,:] |
---|
94 | |
---|
95 | f, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(12,4)) |
---|
96 | |
---|
97 | c=ax1.contourf(xx,zz,ss,20) |
---|
98 | ax1.set_xlim((-.5,.5)), ax1.set_xlabel('x (km)') |
---|
99 | ax1.set_ylim((0.,1.)), ax1.set_ylabel('z (km)') |
---|
100 | plt.colorbar(c,ax=ax1) |
---|
101 | ax1.set_title('Specific entropy at t=%g s (J/K/kg)' % (it*T,)) |
---|
102 | # plt.show() |
---|
103 | |
---|
104 | # plt.figure(figsize=(12,5)) |
---|
105 | c=ax2.contourf(xx,zz,ww,20) |
---|
106 | ax2.set_xlim((-.5,.5)), ax2.set_xlabel('x (km)') |
---|
107 | ax2.set_ylim((0.,1.)), ax2.set_ylabel('z (km)') |
---|
108 | plt.colorbar(c,ax=ax2) |
---|
109 | ax2.set_title('Vertical velocity at t=%g s (m/s)' % (it*T,)) |
---|
110 | # plt.tight_layout() |
---|
111 | # plt.show() |
---|
112 | plt.savefig('fig_NH_3D_bubble/%02d.png'%it) |
---|
113 | |
---|
114 | time1, elapsed1 =time.time(), unst.getvar('elapsed') |
---|
115 | flow = scheme.advance(flow,nt) |
---|
116 | time2, elapsed2 =time.time(), unst.getvar('elapsed') |
---|
117 | factor = 1000./(4*nt) |
---|
118 | print 'ms per call to caldyn_hevi : ', factor*(time2-time1), factor*(elapsed2-elapsed1) |
---|
119 | factor = 1e9/(4*nt*nx*ny*llm) |
---|
120 | print 'nanosec per gridpoint per call to caldyn_hevi : ', factor*(time2-time1), factor*(elapsed2-elapsed1) |
---|
121 | |
---|