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