[631] | 1 | from dynamico.meshes import Cartesian_mesh as Mesh |
---|
[646] | 2 | from dynamico import meshes |
---|
[617] | 3 | from dynamico import unstructured as unst |
---|
| 4 | from dynamico import dyn |
---|
| 5 | from dynamico import time_step |
---|
| 6 | from dynamico import DCMIP |
---|
| 7 | |
---|
| 8 | import math as math |
---|
| 9 | import matplotlib.pyplot as plt |
---|
| 10 | import numpy as np |
---|
| 11 | import time |
---|
| 12 | |
---|
| 13 | unst.setvar('nb_threads', 1) |
---|
| 14 | |
---|
| 15 | def run(mesh,metric,thermo,BC,caldyn_eta,u, m0ik,gas0,Phi0_il, T=300., courant=2.9): |
---|
| 16 | caldyn_thermo = unst.thermo_entropy |
---|
| 17 | g = mesh.dx/metric.dx_g0 |
---|
| 18 | |
---|
[656] | 19 | dt = courant*.5*dx/np.sqrt(gas0.c2.max()) |
---|
| 20 | nt=int(T/dt)+1 |
---|
| 21 | dt=T/nt |
---|
| 22 | print 'nt,dt,dx', nt,dt,dx |
---|
| 23 | |
---|
[617] | 24 | Sik = m0ik*gas0.s |
---|
| 25 | m0ik = m0ik/dx # NB : different definitions of mass field in DYNAMICO / Slice |
---|
| 26 | if caldyn_thermo == unst.thermo_theta: |
---|
| 27 | theta0 = thermo.T0*np.exp(gas0.s/thermo.Cpd) |
---|
| 28 | S0ik = m0ik*theta0 |
---|
| 29 | else: |
---|
| 30 | S0ik = m0ik*gas0.s |
---|
| 31 | |
---|
| 32 | u0=mesh.field_u() |
---|
| 33 | u0[:,range(0,2*nx,2),:] = u*mesh.dx # Doppler shift by u |
---|
| 34 | flow0=(m0ik,S0ik,u0,Phi0_il,0*Phi0_il) |
---|
| 35 | |
---|
[656] | 36 | if False: # time stepping in Python |
---|
| 37 | caldyn = unst.Caldyn_NH(caldyn_thermo,caldyn_eta, mesh,thermo,BC,g) |
---|
| 38 | scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt, a32=0.7) |
---|
| 39 | #scheme = time_step.ARK3(caldyn.bwd_fast_slow, dt) |
---|
| 40 | def next_flow(m,S,u,Phi,W): |
---|
| 41 | m,S,u,Phi,W = scheme.advance((m,S,u,Phi,W),nt) |
---|
| 42 | return m,S,u,Phi,W |
---|
| 43 | else: # time stepping in Fortran |
---|
| 44 | scheme = time_step.ARK2(None, dt, a32=0.7) |
---|
| 45 | caldyn_step = unst.caldyn_step_NH(mesh,scheme,nt, |
---|
| 46 | caldyn_thermo,caldyn_eta,thermo,BC,g) |
---|
| 47 | def next_flow(m,S,u,Phi,W): |
---|
| 48 | caldyn_step.mass[:,:], caldyn_step.theta_rhodz[:,:], caldyn_step.u[:,:] = m,S,u |
---|
| 49 | caldyn_step.geopot[:,:], caldyn_step.W[:,:] = Phi,W |
---|
| 50 | caldyn_step.next() |
---|
| 51 | return (caldyn_step.mass.copy(), caldyn_step.theta_rhodz.copy(), caldyn_step.u.copy(), |
---|
| 52 | caldyn_step.geopot.copy(), caldyn_step.W.copy()) |
---|
[617] | 53 | |
---|
[656] | 54 | m,S,u,Phi,W=flow0 |
---|
| 55 | |
---|
[617] | 56 | for i in range(12): |
---|
| 57 | time1=time.time() |
---|
[656] | 58 | m,S,u,Phi,W = next_flow(m,S,u,Phi,W) |
---|
[617] | 59 | time2=time.time() |
---|
| 60 | print 'ms per time step : ', 1000*(time2-time1)/nt |
---|
| 61 | |
---|
| 62 | w=.5*(W[:,0:llm]+W[:,1:llm+1])/m |
---|
[656] | 63 | # junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.) |
---|
[617] | 64 | # zz=caldyn.geopot[:,0:llm]/metric.g0/1000 |
---|
| 65 | zz=Phi[:,0:llm]/metric.g0/1000 |
---|
| 66 | print 'ptop, model top (m) :', unst.getvar('ptop'), Phi.max()/unst.getvar('g') |
---|
| 67 | |
---|
| 68 | def plt_axes(): |
---|
| 69 | plt.plot(xx,zz[:,llm/4]) |
---|
| 70 | plt.plot(xx,zz[:,llm/2]) |
---|
| 71 | plt.plot(xx,zz[:,3*llm/4]) |
---|
| 72 | plt.ylim(0,ztop/1000.) |
---|
| 73 | plt.colorbar() |
---|
| 74 | |
---|
[656] | 75 | # plt.figure(figsize=(12,3)) |
---|
| 76 | # ucomp = mesh.ucomp(u) |
---|
| 77 | # plt.pcolor(xx,zz,caldyn.pk) |
---|
| 78 | # plt_axes() |
---|
| 79 | # plt.title("T(t = %f s)" % ((i+1)*dt*nt)) |
---|
| 80 | # plt.savefig('fig_slice_GW_NH/T%02d.png'%i) |
---|
| 81 | # plt.close() |
---|
| 82 | # |
---|
[617] | 83 | plt.figure(figsize=(12,3)) |
---|
| 84 | ucomp = mesh.ucomp(u) |
---|
| 85 | plt.pcolor(xx,zz,ucomp[0,:,:]/dx) |
---|
| 86 | plt_axes() |
---|
| 87 | plt.title("u(t = %f s)" % ((i+1)*dt*nt)) |
---|
| 88 | plt.savefig('fig_slice_GW_NH/U%02d.png'%i) |
---|
| 89 | plt.close() |
---|
| 90 | |
---|
| 91 | plt.figure(figsize=(12,3)) |
---|
| 92 | plt.pcolor(xx,zz,w) |
---|
| 93 | plt_axes() |
---|
| 94 | plt.title("w(t = %f s)" % ((i+1)*dt*nt)) |
---|
| 95 | plt.savefig('fig_slice_GW_NH/W%02d.png'%i) |
---|
| 96 | plt.close() |
---|
| 97 | |
---|
| 98 | def deform(x,eta): # x and eta go from 0 to 1 |
---|
| 99 | h0, dd, xi = 2000., 5000., 4000. |
---|
| 100 | cos2=lambda x: np.cos(np.pi*x)**2 |
---|
| 101 | gauss=lambda x: np.exp(-x*x) |
---|
| 102 | xx = x*Lx |
---|
| 103 | phis = h0*gauss(xx/dd)*cos2(xx/xi) |
---|
| 104 | return (phis/ztop)*np.sin(np.pi*eta) |
---|
| 105 | |
---|
| 106 | #Lx, nx, llm = 3e5, 100, 1e4, 30 |
---|
[656] | 107 | Lx, nx, ztop, llm = 3e5, 300, 3e4, 80 |
---|
[617] | 108 | nqdyn, ny, dx = 1, 1, Lx/nx |
---|
[631] | 109 | mesh = Mesh(nx,ny,llm,nqdyn,nx*dx,ny*dx,0.) |
---|
[617] | 110 | xx,ll = mesh.xx[:,0,:]/1000, mesh.ll[:,0,:] |
---|
| 111 | |
---|
[646] | 112 | metric, thermo, BC, m0ik, gas0, Phi0_il, u0_jk = DCMIP.DCMIP31(Lx, nx, llm, deform=deform, dTheta=1.) |
---|
| 113 | #metric, thermo, BC, m0ik, gas0, Phi0_il, u0_jk = DCMIP.DCMIP21(Lx, nx, llm, deform=deform, h0=0.) |
---|
[617] | 114 | |
---|
| 115 | BC.rho_bot = 1e6*BC.rho_bot # large stiffness |
---|
| 116 | |
---|
[646] | 117 | mass_bl,mass_dak,mass_dbk = meshes.compute_hybrid_coefs(m0ik) |
---|
| 118 | unst.ker.dynamico_init_hybrid(mass_bl,mass_dak,mass_dbk) |
---|
[617] | 119 | |
---|
| 120 | run(mesh,metric,thermo,BC,unst.eta_mass,20., m0ik,gas0,Phi0_il, courant=2.) # lower courant number required if h0=2000m |
---|
| 121 | #run(mesh,metric,thermo,BC,unst.eta_lag,20., m0ik,gas0,Phi0_il) |
---|