1 | from dynamico import meshes |
---|
2 | from dynamico import unstructured as unst |
---|
3 | from dynamico import time_step |
---|
4 | from dynamico import precision as prec |
---|
5 | |
---|
6 | import math as math |
---|
7 | import matplotlib.pyplot as plt |
---|
8 | import numpy as np |
---|
9 | import netCDF4 as cdf |
---|
10 | |
---|
11 | #--------------------------------- Main program ----------------------------- |
---|
12 | |
---|
13 | nx,ny,llm,nqdyn=128,128,1,1 |
---|
14 | Lx,Ly,g,f = 8.,8.,1.,1. |
---|
15 | dx,dy=Lx/nx,Ly/ny |
---|
16 | |
---|
17 | filename, llm, nqdyn, g, f, radius = 'cart_%03d_%03d.nc'%(nx,ny), 1, 1, 1., 1., None |
---|
18 | unst.setvar('g',g) |
---|
19 | |
---|
20 | print 'Reading Cartesian mesh ...' |
---|
21 | meshfile = meshes.DYNAMICO_Format(filename) |
---|
22 | |
---|
23 | def coriolis(lon,lat): |
---|
24 | return f+0.*lon |
---|
25 | |
---|
26 | mesh=meshes.Unstructured_Mesh(meshfile, llm, nqdyn, radius, coriolis) |
---|
27 | caldyn = unst.Caldyn_RSW(mesh) |
---|
28 | |
---|
29 | xx,yy = mesh.lon_i, mesh.lat_i |
---|
30 | |
---|
31 | x1,x2,yy = xx-1., xx+1., yy |
---|
32 | u0 = mesh.field_u() # flow initally at rest |
---|
33 | h0 = 1+0.1*(np.exp(-2.*(x1*x1+yy*yy))+ |
---|
34 | np.exp(-2.*(x2*x2+yy*yy))) |
---|
35 | #h0 = 1+0.1*(np.exp(-5.*yy*yy)) |
---|
36 | |
---|
37 | flow0=prec.asnum([h0,u0]) |
---|
38 | |
---|
39 | cfl = .8 |
---|
40 | dt = cfl/math.sqrt((nx/Lx)**2+(ny/Ly)**2) |
---|
41 | |
---|
42 | T=10. |
---|
43 | N=int(T/dt)+1 |
---|
44 | dt=T/N |
---|
45 | print N,dt,Lx/nx |
---|
46 | #scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt, precision=unst.np_num) |
---|
47 | scheme = time_step.RK4(caldyn.bwd_fast_slow, dt, precision=unst.np_num) |
---|
48 | |
---|
49 | flow=flow0 |
---|
50 | |
---|
51 | def minmax(name, x): print('Min/max %s :'%name, x.min(), x.max() ) |
---|
52 | def reshape(data): return data.reshape((nx,ny)) |
---|
53 | |
---|
54 | x, y = map(reshape, (xx,yy) ) |
---|
55 | |
---|
56 | for i in range(10): |
---|
57 | h,u=flow |
---|
58 | flow, fast, slow = caldyn.bwd_fast_slow(flow, prec.zero) |
---|
59 | junk, du_fast = fast |
---|
60 | dh, du_slow = slow |
---|
61 | # minmax('lon',mesh.lon_i) |
---|
62 | # minmax('lat',mesh.lat_i) |
---|
63 | # minmax('x',xx) |
---|
64 | # minmax('y',yy) |
---|
65 | minmax('PV', caldyn.qv-1.) |
---|
66 | # minmax('geopot', caldyn.geopot) |
---|
67 | # minmax('du_fast', du_fast) |
---|
68 | plt.figure(); plt.pcolor(x,y,reshape(caldyn.qv)) |
---|
69 | plt.colorbar(); plt.title('potential vorticity') |
---|
70 | plt.savefig('fig_RSW_2D_mesh/%03d.png'%i) |
---|
71 | plt.close() |
---|
72 | # plt.figure(); plt.pcolor(mesh.x,mesh.y,h) |
---|
73 | # plt.colorbar(); plt.title('h') |
---|
74 | # plt.show() |
---|
75 | # plt.pcolor(x,y,vcomp(u)/dx) |
---|
76 | # plt.colorbar(); plt.title('v') |
---|
77 | # plt.show() |
---|
78 | for j in range(5): |
---|
79 | unst.elapsed=0. |
---|
80 | flow = scheme.advance(flow,N) |
---|
81 | # flow = scheme.advance(flow,5) |
---|
82 | # flops |
---|
83 | # mass flux : 1add+1mul per edge => 4 |
---|
84 | # div U : 4 add per cell => 4 |
---|
85 | # KE : 4*(2add+1mul) per cell => 12 |
---|
86 | # grad KE : 1 add per edge => 2 |
---|
87 | # grad h : 1 add per edge => 2 |
---|
88 | # qv : 4+4+1 add +4mul + 1div per cell => 10 |
---|
89 | # qu : 1add+1mul per edge => 4 |
---|
90 | # TrisK : 4add+4mul+4add+1add per edge => 26 |
---|
91 | # Total : 64 FLOPS/cell |
---|
92 | print i,j, unst.elapsed, 100.*(1.-unst.getvar('elapsed')/unst.elapsed) |
---|
93 | print 'GFlops', 64*(N*nx*ny)/unst.elapsed/1e9 |
---|
94 | unst.setvar('elapsed',0.) |
---|