Changeset 656 for codes/icosagcm/devel

Timestamp:

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

Author:

dubos

Message:

devel/unstructured : some improvement of test cases

Location:

codes/icosagcm/devel/Python/test/py

Files:

• NH_3D_bubble.py (modified) (4 diffs)
• slice_GW_NH.py (modified) (4 diffs)
• slice_GW_hydro.py (modified) (1 diff)

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

@@ -54 +54 @@
 
 #Lx, nx, llm, thetac, T, Nslice, courant = 2000., 100, 50, 30., 5., 10, 2.8
-Lx, nx, llm, thetac, T, Nslice, courant = 2000., 8, 80, 30, 5., 10, 2.8
+Lx, nx, llm, thetac, T, Nslice, courant = 2000., 20, 79, 30, 5., 10, 2.8
 #Lx, nx, llm, thetac, T, Nslice, courant = 3000., 75, 25, -30, 5., 10, 2.8
 
@@ -72 +72 @@
 print 'Time step : %d x %g s' % (nt,dt)
 
+#caldyn_thermo, caldyn_eta = unst.thermo_theta, unst.eta_mass
 caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass
 #caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_lag
@@ -95 +96 @@
     
 m,S,u,Phi,W=flow0
+if caldyn_thermo == unst.thermo_theta:
+    s=S/m
+    theta = thermo.T0*np.exp(s/thermo.Cpd)
+    S=m*theta
+    title_format = 'Potential temperature at t=%g s (K)'
+else:
+    title_format = 'Specific entropy at t=%g s (J/K/kg)'
+
 w=mesh.field_mass()
 z=mesh.field_mass()
@@ -114 +123 @@
     ax1.set_ylim((0.,1.)), ax1.set_ylabel('z (km)')        
     plt.colorbar(c,ax=ax1)
-    ax1.set_title('Specific entropy at t=%g s (J/K/kg)' % (it*T,))
+    ax1.set_title(title_format % (it*T,))
+
     #    plt.show()
     

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.)

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

@@ -67 +67 @@
         plt.close()
 
-Lx, nx, llm = 3e5, 300, 20
+Lx, nx, llm = 3e5, 300, 80
 nqdyn, ny, dx = 1, 1, Lx/nx
 mesh = Mesh(nx,ny,llm,nqdyn,nx*dx,ny*dx,0.)