Changeset 642 for codes/icosagcm/devel/Python/test/py/slice_GW_hydro.py

Timestamp:

12/19/17 15:26:51 (7 years ago)

Author:

dubos

Message:

devel/unstructured : bubble test case with Fortran time stepping

File:

• codes/icosagcm/devel/Python/test/py/slice_GW_hydro.py (modified) (3 diffs)

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

@@ -15 +15 @@
     caldyn_thermo = unst.thermo_entropy
     g = mesh.dx/metric.dx_g0
-    caldyn = unst.Caldyn_HPE(caldyn_thermo,caldyn_eta, mesh,metric,thermo,BC,g)
 
     Sik = m0ik*gas0.s
@@ -25 +24 @@
         S0ik = m0ik*gas0.s
 
-    u0=mesh.field_u() 
+    u0=mesh.field_u()
     u0[:,range(0,2*nx,2),:] = 100.*mesh.dx # Doppler shift 100 m/s
     flow0=(m0ik,S0ik,u0)
@@ -35 +34 @@
     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)
 
-    flow=flow0
+    m,S,u=flow0
+
+    if False: # time stepping in Python
+        caldyn = unst.Caldyn_HPE(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):
+            m,S,u = scheme.advance((m,S,u),nt)
+            return m,S,u,caldyn.geopot
+    else: # time stepping in Fortran
+        scheme = time_step.ARK2(None, dt, a32=0.7)
+        caldyn_step = unst.caldyn_step_HPE(mesh,scheme,nt, caldyn_thermo,caldyn_eta,thermo,BC,g)
+        def next_flow(m,S,u):
+            caldyn_step.mass[:,:], caldyn_step.theta_rhodz[:,:], caldyn_step.u[:,:] = m,S,u
+            caldyn_step.next()
+            return caldyn_step.mass, caldyn_step.theta_rhodz, caldyn_step.u, caldyn_step.geopot
+        
     for i in range(10):
         time1=time.time()
-        flow = scheme.advance(flow,nt)
+        m,S,u,geopot = next_flow(m,S,u)
         time2=time.time()
-        m,S,u = flow
-        print 'ms per time step : ', 1000*(time2-time1)/nt
-        print 'ptop, model top (m) :', unst.getvar('ptop'), caldyn.geopot.max()/unst.getvar('g')
-        junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.)
-        zz=caldyn.geopot[:,0:llm]/metric.g0/1000
+        print 'ms per time step : ', 1000*(time2-time1)/nt, 1000*unst.getvar('elapsed')/(nt*(i+1))
+        print 'ptop, model top (m) :', unst.getvar('ptop'), geopot.max()/unst.getvar('g')
+        #        junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.)
+        zz=geopot[:,0:llm]/metric.g0/1000
         plt.figure(figsize=(12,3))
         ucomp = mesh.ucomp(u)