1 | from dynamico.meshes import Cartesian_mesh as Mesh |
---|
2 | from dynamico import unstructured as unst |
---|
3 | from dynamico import dyn |
---|
4 | from dynamico import time_step |
---|
5 | from dynamico import DCMIP |
---|
6 | |
---|
7 | import math as math |
---|
8 | import matplotlib.pyplot as plt |
---|
9 | import numpy as np |
---|
10 | import time |
---|
11 | |
---|
12 | unst.setvar('nb_threads', 1) |
---|
13 | |
---|
14 | def run(mesh,metric,thermo,BC,caldyn_eta, m0ik,gas0): |
---|
15 | caldyn_thermo = unst.thermo_entropy |
---|
16 | g = mesh.dx/metric.dx_g0 |
---|
17 | caldyn = unst.Caldyn_HPE(caldyn_thermo,caldyn_eta, mesh,metric,thermo,BC,g) |
---|
18 | |
---|
19 | Sik = m0ik*gas0.s |
---|
20 | m0ik = m0ik/dx # NB : different definitions of mass field in DYNAMICO / Slice |
---|
21 | if caldyn_thermo == unst.thermo_theta: |
---|
22 | theta0 = thermo.T0*np.exp(gas0.s/thermo.Cpd) |
---|
23 | S0ik = m0ik*theta0 |
---|
24 | else: |
---|
25 | S0ik = m0ik*gas0.s |
---|
26 | |
---|
27 | u0=mesh.field_u() |
---|
28 | u0[:,range(0,2*nx,2),:] = 100.*mesh.dx # Doppler shift 100 m/s |
---|
29 | flow0=(m0ik,S0ik,u0) |
---|
30 | |
---|
31 | #T, courant = 60., 1. |
---|
32 | T, courant = 60., 2.9 |
---|
33 | dt = courant*.5*dx/np.sqrt(gas0.c2.max()) |
---|
34 | nt=int(T/dt)+1 |
---|
35 | dt=T/nt |
---|
36 | print 'nt,dt,dx', nt,dt,dx |
---|
37 | scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt, a32=0.7) |
---|
38 | #scheme = time_step.ARK3(caldyn.bwd_fast_slow, dt) |
---|
39 | |
---|
40 | flow=flow0 |
---|
41 | for i in range(10): |
---|
42 | time1=time.time() |
---|
43 | flow = scheme.advance(flow,nt) |
---|
44 | time2=time.time() |
---|
45 | m,S,u = flow |
---|
46 | print 'ms per time step : ', 1000*(time2-time1)/nt |
---|
47 | print 'ptop, model top (m) :', unst.getvar('ptop'), caldyn.geopot.max()/unst.getvar('g') |
---|
48 | junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.) |
---|
49 | zz=caldyn.geopot[:,0:llm]/metric.g0/1000 |
---|
50 | plt.figure(figsize=(12,3)) |
---|
51 | ucomp = mesh.ucomp(u) |
---|
52 | plt.pcolor(xx,zz,ucomp[0,:,:]/dx) |
---|
53 | plt.ylim(0,10.) |
---|
54 | plt.colorbar() |
---|
55 | plt.title("u(t = %f s)" % ((i+1)*dt*nt)) |
---|
56 | plt.savefig('fig_slice_GW_hydro/%02d.png'%i) |
---|
57 | plt.close() |
---|
58 | |
---|
59 | Lx, nx, llm = 3e5, 300, 20 |
---|
60 | nqdyn, ny, dx = 1, 1, Lx/nx |
---|
61 | mesh = Mesh(nx,ny,llm,nqdyn,nx*dx,ny*dx,0.) |
---|
62 | xx,ll = mesh.xx[:,0,:]/1000, mesh.ll[:,0,:] |
---|
63 | metric, thermo, BC, m0ik, gas0, Phi0_il, u0_jk = DCMIP.DCMIP31(Lx, nx, llm) |
---|
64 | |
---|
65 | run(mesh,metric,thermo,BC,unst.eta_lag, m0ik,gas0) |
---|