Changeset 656 for codes/icosagcm/devel/Python/test/py/slice_GW_NH.py

Timestamp:

12/30/17 01:06:26 (7 years ago)

Author:

dubos

Message:

devel/unstructured : some improvement of test cases

File:

• codes/icosagcm/devel/Python/test/py/slice_GW_NH.py (modified) (4 diffs)

codes/icosagcm/devel/Python/test/py/slice_GW_NH.py

@@ -16 +16 @@
     caldyn_thermo = unst.thermo_entropy
     g = mesh.dx/metric.dx_g0
-    caldyn = unst.Caldyn_NH(caldyn_thermo,caldyn_eta, mesh,thermo,BC,g)
+
+    dt = courant*.5*dx/np.sqrt(gas0.c2.max())
+    nt=int(T/dt)+1
+    dt=T/nt
+    print 'nt,dt,dx', nt,dt,dx
 
     Sik = m0ik*gas0.s
@@ -30 +34 @@
     flow0=(m0ik,S0ik,u0,Phi0_il,0*Phi0_il)
 
-    dt = courant*.5*dx/np.sqrt(gas0.c2.max())
-    nt=int(T/dt)+1
-    dt=T/nt
-    print 'nt,dt,dx', nt,dt,dx
-    scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt, a32=0.7)
-    #scheme = time_step.ARK3(caldyn.bwd_fast_slow, dt)
+    if False: # time stepping in Python
+        caldyn = unst.Caldyn_NH(caldyn_thermo,caldyn_eta, mesh,thermo,BC,g)
+        scheme = time_step.ARK2(caldyn.bwd_fast_slow, dt, a32=0.7)
+        #scheme = time_step.ARK3(caldyn.bwd_fast_slow, dt)
+        def next_flow(m,S,u,Phi,W):
+            m,S,u,Phi,W = scheme.advance((m,S,u,Phi,W),nt)
+            return m,S,u,Phi,W
+    else: # time stepping in Fortran
+        scheme = time_step.ARK2(None, dt, a32=0.7)
+        caldyn_step = unst.caldyn_step_NH(mesh,scheme,nt, 
+                                          caldyn_thermo,caldyn_eta,thermo,BC,g)
+        def next_flow(m,S,u,Phi,W):
+            caldyn_step.mass[:,:], caldyn_step.theta_rhodz[:,:], caldyn_step.u[:,:] = m,S,u
+            caldyn_step.geopot[:,:], caldyn_step.W[:,:] = Phi,W
+            caldyn_step.next()
+            return (caldyn_step.mass.copy(), caldyn_step.theta_rhodz.copy(), caldyn_step.u.copy(),
+                    caldyn_step.geopot.copy(), caldyn_step.W.copy())
 
-#    nt=1 # FIXME
-    flow=flow0
+    m,S,u,Phi,W=flow0
+
     for i in range(12):
         time1=time.time()
-        flow = scheme.advance(flow,nt)
+        m,S,u,Phi,W = next_flow(m,S,u,Phi,W)
         time2=time.time()
         print 'ms per time step : ', 1000*(time2-time1)/nt
 
-        m,S,u,Phi,W = flow
         w=.5*(W[:,0:llm]+W[:,1:llm+1])/m
-        junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.)
+#        junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.)
 #        zz=caldyn.geopot[:,0:llm]/metric.g0/1000
         zz=Phi[:,0:llm]/metric.g0/1000
@@ -59 +73 @@
             plt.colorbar()
 
-        plt.figure(figsize=(12,3))
-        ucomp = mesh.ucomp(u)
-        plt.pcolor(xx,zz,caldyn.pk)
-        plt_axes()
-        plt.title("T(t = %f s)" % ((i+1)*dt*nt))
-        plt.savefig('fig_slice_GW_NH/T%02d.png'%i)
-        plt.close()
-
+#        plt.figure(figsize=(12,3))
+#        ucomp = mesh.ucomp(u)
+#        plt.pcolor(xx,zz,caldyn.pk)
+#        plt_axes()
+#        plt.title("T(t = %f s)" % ((i+1)*dt*nt))
+#        plt.savefig('fig_slice_GW_NH/T%02d.png'%i)
+#        plt.close()
+#
         plt.figure(figsize=(12,3))
         ucomp = mesh.ucomp(u)
@@ -91 +105 @@
 
 #Lx, nx, llm = 3e5, 100, 1e4, 30
-Lx, nx, ztop, llm = 2e5, 400, 3e4, 60
+Lx, nx, ztop, llm = 3e5, 300, 3e4, 80
 nqdyn, ny, dx = 1, 1, Lx/nx
 mesh = Mesh(nx,ny,llm,nqdyn,nx*dx,ny*dx,0.)