[631] | 1 | from dynamico.meshes import Cartesian_mesh as Mesh |
---|
[617] | 2 | from dynamico import unstructured as unst |
---|
| 3 | from dynamico import dyn |
---|
| 4 | from dynamico import time_step |
---|
| 5 | from dynamico import DCMIP |
---|
[754] | 6 | from dynamico.precision import zero, np_num |
---|
| 7 | |
---|
[617] | 8 | import math as math |
---|
| 9 | import matplotlib.pyplot as plt |
---|
| 10 | import numpy as np |
---|
| 11 | import time |
---|
| 12 | |
---|
[754] | 13 | nx,ny,llm,nqdyn=128,192,1,1 |
---|
| 14 | Lx,Ly,g,f = 8.,12.,1.,1. |
---|
[617] | 15 | dx,dy=Lx/nx,Ly/ny |
---|
| 16 | |
---|
| 17 | unst.setvar('g',g) |
---|
[631] | 18 | mesh = Mesh(nx,ny,llm,nqdyn,Lx,Ly,f) |
---|
[617] | 19 | caldyn = unst.Caldyn_RSW(mesh) |
---|
| 20 | |
---|
| 21 | x1,x2,yy = mesh.xx[:,:,0]-1., mesh.xx[:,:,0]+1., mesh.yy[:,:,0] |
---|
[747] | 22 | h0, u0 = mesh.field_mass(), mesh.field_u() # flow initally at rest |
---|
| 23 | h0[:] = 1+0.1*(np.exp(-2.*(x1*x1+yy*yy))+ |
---|
[617] | 24 | np.exp(-2.*(x2*x2+yy*yy))) |
---|
[773] | 25 | #h0[:] = 1+0.1*(np.exp(-2.*yy*yy)) |
---|
| 26 | |
---|
[617] | 27 | flow0=(h0,u0) |
---|
| 28 | |
---|
| 29 | cfl = .8 |
---|
| 30 | dt = cfl/math.sqrt((nx/Lx)**2+(ny/Ly)**2) |
---|
| 31 | |
---|
| 32 | T=10. |
---|
| 33 | N=int(T/dt)+1 |
---|
| 34 | dt=T/N |
---|
| 35 | print N,dt,Lx/nx |
---|
[754] | 36 | scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt, precision=np_num) |
---|
[617] | 37 | |
---|
[747] | 38 | print 'types : yy, u0, h0, u0', yy.dtype, u0.dtype, h0.dtype, u0.dtype |
---|
| 39 | |
---|
[617] | 40 | flow=flow0 |
---|
| 41 | for i in range(10): |
---|
| 42 | h,u=flow |
---|
[747] | 43 | print 'types : h,u', h.dtype, u.dtype |
---|
[754] | 44 | caldyn.bwd_fast_slow(flow, zero) |
---|
[617] | 45 | plt.figure(); plt.pcolor(mesh.x,mesh.y,caldyn.qv) |
---|
| 46 | plt.colorbar(); plt.title('potential vorticity') |
---|
| 47 | plt.savefig('fig_RSW_2D/%02d.png'%i) |
---|
| 48 | # plt.figure(); plt.pcolor(mesh.x,mesh.y,h) |
---|
| 49 | # plt.colorbar(); plt.title('h') |
---|
| 50 | # plt.show() |
---|
| 51 | # plt.pcolor(x,y,vcomp(u)/dx) |
---|
| 52 | # plt.colorbar(); plt.title('v') |
---|
| 53 | # plt.show() |
---|
| 54 | for j in range(5): |
---|
| 55 | unst.elapsed=0. |
---|
| 56 | flow = scheme.advance(flow,N) |
---|
[747] | 57 | h,u=flow |
---|
| 58 | print 'types : h,u', h.dtype, u.dtype |
---|
[617] | 59 | # flops |
---|
| 60 | # mass flux : 1add+1mul per edge => 4 |
---|
| 61 | # div U : 4 add per cell => 4 |
---|
| 62 | # KE : 4*(2add+1mul) per cell => 12 |
---|
| 63 | # grad KE : 1 add per edge => 2 |
---|
| 64 | # grad h : 1 add per edge => 2 |
---|
| 65 | # qv : 4+4+1 add +4mul + 1div per cell => 10 |
---|
| 66 | # qu : 1add+1mul per edge => 4 |
---|
| 67 | # TrisK : 4add+4mul+4add+1add per edge => 26 |
---|
| 68 | # Total : 64 FLOPS/cell |
---|
| 69 | print i,j, unst.elapsed, 100.*(1.-unst.getvar('elapsed')/unst.elapsed) |
---|
| 70 | print 'GFlops', 64*(N*nx*ny)/unst.elapsed/1e9 |
---|
| 71 | unst.setvar('elapsed',0.) |
---|