1 | from dynamico.meshes import Cartesian_mesh as Mesh |
---|
2 | import dynamico.dyn as dyn |
---|
3 | import dynamico.time_step as time_step |
---|
4 | import dynamico.DCMIP as DCMIP |
---|
5 | import dynamico.advect as advect |
---|
6 | from dynamico import unstructured as unst |
---|
7 | |
---|
8 | import math as math |
---|
9 | import numpy as np |
---|
10 | import matplotlib.pyplot as plt |
---|
11 | import matplotlib.animation as manimation |
---|
12 | |
---|
13 | |
---|
14 | def reldiff(a,b): |
---|
15 | return (a-b)/(np.absolute(a).max()+np.absolute(b).max()+1e-20) |
---|
16 | |
---|
17 | def go(model,courant,scheme,start): |
---|
18 | def plotter(metric, flow, scalars, t): |
---|
19 | mik, gas_ik, ujk, Phi_il, Wil = model.diagnose(*flow) |
---|
20 | Q1, Q2 = scalars |
---|
21 | |
---|
22 | mil = metric.ml_ext(mik) |
---|
23 | wil = Wil/mil*metric.g0 |
---|
24 | wik = metric.mk_ext(Wil)/mik*metric.g0 |
---|
25 | xik = metric.xik/1000 |
---|
26 | xjl = metric.xjl/1000 |
---|
27 | xil = metric.xil/1000 |
---|
28 | zil = Phi_il/metric.g0/1000 |
---|
29 | zik = metric.mk_int(zil) |
---|
30 | zjl = metric.mj(zil) |
---|
31 | print 's ',np.max(gas_ik.s),np.min(gas_ik.s) |
---|
32 | print 'ujk ',np.max(ujk),np.min(ujk) |
---|
33 | print 'wil ',np.max(wil),np.min(wil) |
---|
34 | sik = gas_ik.s |
---|
35 | sik = .5*(abs(sik)+sik) |
---|
36 | |
---|
37 | plt.clf() |
---|
38 | plt.contourf(xik, zik, sik, 20 ) |
---|
39 | # plt.contourf(xik, zik, wik, 20 ) |
---|
40 | plt.xlim((-1,1)), plt.xlabel('x (km)') |
---|
41 | plt.ylim((0,ztop/1000.)), plt.ylabel('z (km)') |
---|
42 | plt.colorbar() |
---|
43 | plt.title("t = %f s" % (t)) |
---|
44 | plt.savefig('fig_bubble/bubble.png',bbox_inches='tight') |
---|
45 | writer.grab_frame() |
---|
46 | |
---|
47 | transport = advect.SplitTransport(model.metric, advect.HorizVanLeer, advect.VertVanLeer) |
---|
48 | replace = lambda fd,fv : fv |
---|
49 | return DCMIP.run_with_scalars(model,scheme,transport,replace, plotter, courant, start,scalars, T1, N1, N2) |
---|
50 | |
---|
51 | Lx, ztop, nx, llm = 2000.,1000., 100, 50 |
---|
52 | nqdyn, ny, dx = 1, 1, Lx/nx |
---|
53 | T1, N1, N2 = 5., 1, 200 |
---|
54 | #T1, N1, N2 = 5., 2, 5 |
---|
55 | |
---|
56 | metric, thermo, BC, m0ik, gas0, Phi0_il, u0_jk = DCMIP.thermal_bubble(Lx,nx,llm, ztop=ztop, zc=350.) |
---|
57 | Sik = m0ik*gas0.s |
---|
58 | start = (m0ik, Sik, u0_jk, Phi0_il, 0*Phi0_il) |
---|
59 | scalars = (Sik, m0ik*metric.mk_int(Phi0_il)) |
---|
60 | |
---|
61 | model = dyn.Nonhydro(thermo,metric,BC) |
---|
62 | solver, caldyn_eta = dyn.MassBasedNH(model), unst.eta_mass |
---|
63 | #solver, caldyn_eta = dyn.LagrangianNH(model), unst.eta_lag |
---|
64 | |
---|
65 | rho_bot=BC.rho_bot |
---|
66 | if False: # explicit time-stepping with rigid lid |
---|
67 | def rigid(flow,tau): # impose rigid BCs at top and bottom |
---|
68 | flow, fast, slow = solver.bwd_fast_slow(flow,0.) |
---|
69 | fast[4][:,0]=0. |
---|
70 | fast[4][:,-1]=0. |
---|
71 | return flow,fast,slow |
---|
72 | BC.rho_bot = 100*rho_bot # moderate stiffness |
---|
73 | courant, time_scheme = 2.0, lambda dt : time_step.RK4(rigid, dt) |
---|
74 | else: # HEVI time-stepping with pressure BC |
---|
75 | # Note : parameters BC.rigid_top and BC.rigid_bottom are True by default, but are ignored by the Python code... |
---|
76 | BC.rigid_top=False |
---|
77 | BC.rigid_bottom=False |
---|
78 | BC.rho_bot = 1e6*rho_bot # large stiffness |
---|
79 | # BC.rho_bot = 100*rho_bot # moderate stiffness |
---|
80 | courant, time_scheme = 3.0, lambda dt : time_step.ARK2(solver.bwd_fast_slow, dt, a32=0.7) |
---|
81 | |
---|
82 | FFMpegWriter = manimation.writers['ffmpeg'] |
---|
83 | metadata = dict(title='Aquaplanet', artist='DYNAMICO_LMDZ5', |
---|
84 | comment='Movie support!') |
---|
85 | writer = FFMpegWriter(fps=15, metadata=metadata, bitrate=4000, codec="h263p") |
---|
86 | fig = plt.figure(figsize=(12,5)) |
---|
87 | with writer.saving(fig, "fig_bubble/bubble.avi", 100): |
---|
88 | m0ik, gas0_ik, u0jk, Phi0_il, W0il = go(model, courant, time_scheme, start) |
---|
89 | |
---|
90 | #------------------------------------------------------------------------------------------------------- |
---|
91 | # Now use final state to check that Fortran and Python implementations produce the exact same tendencies |
---|
92 | #------------------------------------------------------------------------------------------------------- |
---|
93 | |
---|
94 | #W0il[:,0]=0. |
---|
95 | #W0il[:,-1]=0. |
---|
96 | |
---|
97 | unst.setvar('nb_threads', 1) |
---|
98 | mesh = Mesh(nx,ny,llm,nqdyn,Lx,ny*dx,0.) |
---|
99 | xx_ik, xx_il, ll = mesh.xx[:,0,:]/1000, mesh.xxp1[:,0,:]/1000, mesh.ll[:,0,:] |
---|
100 | |
---|
101 | mass_bl,mass_dak,mass_dbk=unst.compute_hybrid_coefs(m0ik) |
---|
102 | unst.init_hybrid(mass_bl,mass_dak,mass_dbk) |
---|
103 | |
---|
104 | caldyn_thermo = unst.thermo_entropy |
---|
105 | #caldyn_thermo = unst.thermo_theta |
---|
106 | g = mesh.dx/metric.dx_g0 |
---|
107 | caldyn = unst.Caldyn_NH(caldyn_thermo,caldyn_eta, mesh,metric,thermo,BC,g) |
---|
108 | unst.setvar('rigid',False) |
---|
109 | |
---|
110 | s0ik = gas0_ik.s |
---|
111 | if caldyn_thermo == unst.thermo_theta: |
---|
112 | s0ik = thermo.T0*np.exp(s0ik/thermo.Cpd) |
---|
113 | |
---|
114 | # 2 momentum components in DYNAMICO vs 1 component in slice |
---|
115 | u0jk = u0jk + 100*dx # shift horizontal wind by 100 m/s |
---|
116 | u0_ker = mesh.field_u() |
---|
117 | mesh.set_ucomp(u0_ker,u0jk) |
---|
118 | |
---|
119 | tau=0.1 |
---|
120 | # impose hybdrid distribution of mass (Fortran) |
---|
121 | flow = (m0ik/dx,s0ik*m0ik/dx,u0_ker,Phi0_il,W0il/dx) # NB : different definitions of mass field in DYNAMICO / Python |
---|
122 | flow,fast,slow = caldyn.bwd_fast_slow(flow, 0.) |
---|
123 | # shift Phi from target value Phi_bot |
---|
124 | Phi0_il = Phi0_il+g*100.*np.exp(-10.*xx_il*xx_il) |
---|
125 | m,S,u,Phi,W = flow |
---|
126 | flow = m,S,u,Phi0_il,W |
---|
127 | # compute tendencies (0=Python) |
---|
128 | flow0 = (m*dx,S*dx,u0jk,Phi0_il,W0il, 0., 0.) # NB : different definitions of mass field in DYNAMICO / Python |
---|
129 | flow0, fast0, slow0 = solver.bwd_fast_slow(flow0,tau,True) |
---|
130 | # compute tendencies (Fortran) |
---|
131 | flow,fast,slow = caldyn.bwd_fast_slow(flow, tau) |
---|
132 | |
---|
133 | |
---|
134 | zz_il = Phi0_il/metric.g0/1000. |
---|
135 | zz_ik=.5*(zz_il[:,0:-1]+zz_il[:,1:]) |
---|
136 | print |
---|
137 | print 'ptop, model top (km) :', unst.getvar('ptop'), zz_il.max() |
---|
138 | |
---|
139 | for data0,data,dname in ((flow0,flow,'flow'),(slow0,slow,'slow'),(fast0,fast,'fast')): |
---|
140 | def plot(xx,zz,var,vname): |
---|
141 | title = "%s_%s" % (vname,dname) |
---|
142 | filename = "bubble/bubble_%s.png" % title |
---|
143 | var = 1.*var |
---|
144 | if var.size>1: |
---|
145 | plt.figure(figsize=(12,5)) |
---|
146 | plt.contourf(xx,zz,var,20) |
---|
147 | # plt.pcolor(xx,zz,var) |
---|
148 | plt.xlim((-1.,1.)), plt.xlabel('x (km)') |
---|
149 | plt.ylim((0.,1.)), plt.ylabel('z (km)') |
---|
150 | plt.colorbar() |
---|
151 | plt.title(title) |
---|
152 | plt.savefig(filename,bbox_inches='tight') |
---|
153 | print title, var.min(), var.max() |
---|
154 | if var.max()>1e100: |
---|
155 | raise OverflowError |
---|
156 | |
---|
157 | dm,dS,du,dPhi,dW = data |
---|
158 | du=mesh.ucomp(du)/dx |
---|
159 | du=du[0,:,:] |
---|
160 | |
---|
161 | dm0, dS0, du0, dPhi0, dW0, junk, junk = data0 |
---|
162 | dm0, dS0, du0, dW0 = dm0/dx, dS0/dx, du0/dx, dW0/dx |
---|
163 | |
---|
164 | if dname == 'flow': |
---|
165 | hflux=mesh.ucomp(caldyn.hflux) |
---|
166 | plot(xx_ik, zz_ik, dm,'m') |
---|
167 | plot(xx_ik, zz_ik, du,'u') |
---|
168 | plot(xx_ik, zz_ik, du-du0,'uu') |
---|
169 | plot(xx_ik, zz_ik, dS/m,'s') |
---|
170 | plot(xx_il, zz_il, dW,'w') |
---|
171 | plot(xx_il, zz_il, dW-dW0,'ww') |
---|
172 | else: |
---|
173 | dS=reldiff(dS,dS0) |
---|
174 | dW=dW-dW0 |
---|
175 | dPhi=dPhi-dPhi0 |
---|
176 | du=reldiff(du,du0) |
---|
177 | plot(xx_ik, zz_ik, dS,'dS') |
---|
178 | plot(xx_il, zz_il, dPhi,'dPhi') |
---|
179 | plot(xx_il, zz_il, dW,'dW') |
---|
180 | plot(xx_ik, zz_ik, du,'du') |
---|
181 | |
---|
182 | if unst.getvar('rigid'): |
---|
183 | print 'Fortran BCs : rigid' |
---|
184 | else: |
---|
185 | print 'Fortran BCs : pressure/soft', unst.getvars('ptop','rho_bot','Phi_bot','pbot') |
---|
186 | |
---|
187 | if BC.rigid_top: |
---|
188 | print 'Top BC : rigid' |
---|
189 | else: |
---|
190 | print 'Top BC : pressure', BC.ptop |
---|
191 | |
---|
192 | if BC.rigid_bottom: |
---|
193 | print 'Bottom BC : rigid' |
---|
194 | else: |
---|
195 | print 'Bottom BC : soft', BC.rho_bot[0], BC.Phi_bot[0], BC.pbot[0] |
---|