Changeset 785

Timestamp:

11/21/18 16:17:24 (6 years ago)

Author:

jisesh

Message:

devel/Python : added bilaplacian dissipation for theta + diagnostic of curl

Location:

codes/icosagcm/devel/Python/test

Files:

• py/Baroclinic_3D_ullrich.py (modified) (6 diffs)
• xml/field_def_simple.xml (modified) (1 diff)
• xml/filedef_simple.xml (modified) (1 diff)

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

@@ -25 +25 @@
 parser.add_argument("--T", type=float, default=3600.,
                     help="Length of time slice in seconds")
-parser.add_argument("--N", type=int, default=24,
+parser.add_argument("--N", type=int, default=48,
                     help="Number of time slices")
 parser.add_argument("--log",type=str,default='DEBUG',choices=('DEBUG', 'INFO', 'WARNING', 'ERROR', 'CRITICAL'),
@@ -141 +141 @@
     gas = thermo.set_vs(v,s)
     return gas, w, z
+
+def grad(mesh, s_ik):
+    llm=s_ik.shape[1]
+    edge_num, left, right = mesh.edge_num, mesh.left, mesh.right
+    grad_ek = prec.zeros((edge_num, llm))
+    for e in range(edge_num):
+        left_e = left[e]-1 # Fortran => Python indexing
+        right_e = right[e]-1 # Fortran => Python indexing
+        grad_ek[e,:] = s_ik[right_e,:]-s_ik[left_e,:]
+    return grad_ek
+
+def div(mesh, u_ek):
+    llm=u_ek.shape[1]
+    primal_num, primal_deg, primal_edge = mesh.primal_num, mesh.primal_deg, mesh.primal_edge
+    primal_ne, le_de, Ai = mesh.primal_ne, mesh.le_de, mesh.Ai
+    div_ik = prec.zeros((primal_num, llm))
+    for i in range(primal_num):
+        div_i = 0.
+        for iedge in range(primal_deg[i]):
+            e = primal_edge[i,iedge]-1 # Fortran => Python indexing
+            div_i = div_i + u_ek[e,:]*primal_ne[i,iedge]*le_de[e]
+        div_ik[i,:] = div_i / Ai[i]
+    return div_ik
+
+def curl(mesh, u_ek):
+    llm=u_ek.shape[1]
+    dual_num, dual_deg, dual_edge, dual_ne, Av = mesh.dual_num, mesh.dual_deg, mesh.dual_edge, mesh.dual_ne, mesh.Av
+    curl_vk = prec.zeros((dual_num, llm))
+    for v in range(dual_num):
+        curl_v = 0.
+        for iedge in range(dual_deg[v]):
+            e = dual_edge[v,iedge]-1 # Fortran => Python indexing
+            curl_v = curl_v + u_ek[e,:]*dual_ne[v,iedge]
+        curl_vk[v,:] = curl_v / Av[v]
+    return curl_vk
 
 class Davies:
@@ -189 +224 @@
 
     g, Lx, Ly = 9.81, 4e7, 6e6
-    nx, ny, llm = 200, 30, 25
+    nx, ny, llm = 200, 30, 30
     dx,dy = Lx/nx,Ly/ny
 
@@ -227 +262 @@
             caldyn_step.mass[:,:], caldyn_step.theta_rhodz[:,:], caldyn_step.u[:,:] = m,S,u
             caldyn_step.geopot[:,:], caldyn_step.W[:,:] = Phi,W
+            kappa = dx**4/43200.
             for i in range(nt):
                 caldyn_step.next()
                 davies.relax(llm, caldyn_step, flow0)
+                m,S = caldyn_step.mass, caldyn_step.theta_rhodz
+                s = S/m
+                laps = mesh.field_mass()
+                bilaps = mesh.field_mass()
+                unst.ker.dynamico_scalar_laplacian(s,laps)
+                unst.ker.dynamico_scalar_laplacian(laps,bilaps)
+#                grads = grad(mesh,s)
+#                laps = div(mesh,grads)
+#                gradlaps = grad(mesh,laps) # bilaplacian
+#                bilaps = div(mesh,gradlaps)                
+                caldyn_step.theta_rhodz[:] = S - dt*kappa*bilaps*m # Euler step
             return (caldyn_step.mass.copy(), caldyn_step.theta_rhodz.copy(), caldyn_step.u.copy(),
                     caldyn_step.geopot.copy(), caldyn_step.W.copy())
@@ -254 +301 @@
         
         for i in range(Nslice):
-            context.update_calendar(i)
+            context.update_calendar(i+1)
 
             # Diagnose quantities of interest from prognostic variables m,S,u,Phi,W
             gas, w, z = diagnose(Phi,S,m,W)
-
+            curl_vk = curl(mesh, u)
             # write to disk
             context.send_field_primal('ps', davies.beta_i)
-            context.send_field_primal('temp', gas0.T)
+            context.send_field_primal('temp', gas.T)
             context.send_field_primal('p', gas.p)
             context.send_field_primal('theta', gas.s)
             context.send_field_primal('uz', w)
             context.send_field_primal('z', z)
+            context.send_field_primal('curl', curl_vk)
 
             print( 'ptop, model top (m) :', unst.getvar('ptop'), Phi.max()/unst.getvar('g'))
@@ -277 +325 @@
             print( 'nanosec per gridpoint per full time step : ', factor*(time2-time1), factor*(elapsed2-elapsed1))
 
-        context.update_calendar(Nslice+1) # make sure XIOS writes last iteration
 logging.info('************DONE************')

codes/icosagcm/devel/Python/test/xml/field_def_simple.xml

@@ -23 +23 @@
         <field id="z" />
         <field id="theta" />
+        <field id="curl" />
       </field_group>
      

codes/icosagcm/devel/Python/test/xml/filedef_simple.xml

@@ -20 +20 @@
       <field field_ref="z" />
       <field field_ref="theta" />
+      <field field_ref="curl" />
       <field field_ref="phis" operation="once" freq_offset="0ts" name="PHIS"   standard_name="surface_geopotential" long_name="Surface geopotential" unit="m2/m2"/>      
     </field_group>