Changeset 756 for codes/icosagcm

Timestamp:

10/09/18 14:28:53 (6 years ago)

Author:

dubos

Message:

devel/Python : moved DYNAMICO_Format to meshes.py + cleanup

Location:

codes/icosagcm/devel/Python

Files:

• dynamico/meshes.py (modified) (4 diffs)
• test/py/RSW_2D_mesh.py (modified) (6 diffs)

codes/icosagcm/devel/Python/dynamico/meshes.py

@@ -190 +190 @@
     def vcomp(self,u): return u[:,range(1,2*self.nx,2),:]
 
-#---------------------- MPAS fully unstructured mesh ------------------------
+#---------------------- DYNAMICO format for fully unstructured mesh and curvilinear meshes ----------------------
 
 class Mesh_Format(Base_class):
@@ -203 +203 @@
             return ret
     
+class DYNAMICO_Format(Mesh_Format):
+    """ Helps class Unstructured_Mesh to read a DYNAMICO mesh file."""
+    start_index=1 # indexing starts at 1 as in Fortran
+    def __init__(self,gridfilename):
+        self.nc = cdf.Dataset(gridfilename, "r")
+        print('#NC4: Opening file:',gridfilename)
+    def get_pdim(self,comm,name): return parallel.PDim(self.nc.dimensions[name], comm)
+    def getvar(self,name): return self.nc.variables[name][:]
+    def get_pvar(self,pdim,name):
+        # provides parallel access to a NetCDF variable, with name translation
+        # first deal with special cases
+        if name == 'edge_deg': return parallel.CstPArray1D(pdim, np.int32, 2)
+        # general case
+        return parallel.PArray(pdim, self.nc.variables[name])
+    def normalize(self, mesh, radius): pass
+
+#---------------------- MPAS fully unstructured mesh ------------------------
+
 class MPAS_Format(Mesh_Format):
+    """ Helps class Unstructured_Mesh to read a MPAS mesh file."""
     start_index=1 # indexing starts at 1 as in Fortran
     translate= {
-        'primal_num':'nCells',
-        'edge_num':'nEdges',
-        'dual_num':'nVertices',
+        'primal_num':'nCells', 'edge_num':'nEdges', 'dual_num':'nVertices',
         'primal_deg':'nEdgesOnCell', 'trisk_deg':'nEdgesOnEdge',
-        'primal_edge':'edgesOnCell', 'left_right':'cellsOnEdge',
-        'dual_edge':'edgesOnVertex', 'down_up':'verticesOnEdge',
+        'primal_edge':'edgesOnCell', 'dual_edge':'edgesOnVertex', 
         'primal_vertex':'verticesOnCell', 'dual_vertex':'cellsOnVertex',
-        'trisk':'edgesOnEdge', 'down_up':'verticesOnEdge', 'left_right':'cellsOnEdge',
+        'trisk':'edgesOnEdge', 'edge_down_up':'verticesOnEdge', 'edge_left_right':'cellsOnEdge',
         'le':'dvEdge', 'de':'dcEdge', 'Ai':'areaCell', 'Av':'areaTriangle',
         'lat_i':'latCell','lon_i':'lonCell','lat_v':'latVertex','lon_v':'lonVertex',
@@ -319 +335 @@
         # get indices for stencils 
         # primal <-> edges
-        primal_edge, left_right = getvars('primal_edge','left_right')
+        primal_edge, left_right = getvars('primal_edge','edge_left_right')
         left,right = left_right[:,0], left_right[:,1]
         primal_ne = compute_ne(primal_num,primal_deg,primal_edge,left,right)
         # dual <-> edges
-        dual_edge, down_up = getvars('dual_edge','down_up')
+        dual_edge, down_up = getvars('dual_edge','edge_down_up')
         down,up = down_up[:,0], down_up[:,1]
         dual_ne = compute_ne(dual_num,dual_deg,dual_edge,up,down)
@@ -357 +373 @@
             dim_primal, 'primal_deg', 'primal_vertex', 'primal_edge', 'lon_i', 'lat_i', 'Ai' )
         edge_deg, trisk_deg, left_right, down_up, trisk, le, de, wee, lon_e, lat_e, angle_e = get_pvars(
-             dim_edge, 'edge_deg', 'trisk_deg', 'left_right', 'down_up', 
+             dim_edge, 'edge_deg', 'trisk_deg', 'edge_left_right', 'edge_down_up', 
              'trisk', 'le', 'de', 'wee', 'lon_e', 'lat_e', 'angle_e')
         dual_deg, dual_vertex, dual_edge, Riv2, lon_v, lat_v, Av = get_pvars(

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

@@ -2 +2 @@
 from dynamico import unstructured as unst
 from dynamico import time_step
+from dynamico import precision as prec
 
 import math as math
@@ -7 +8 @@
 import numpy as np
 import netCDF4 as cdf
-
-class DYNAMICO_Format(meshes.Mesh_Format):
-    """ Netcdf input file handling.
-    Creates a dictionary of names corresponding to original data file."""
-
-    start_index=1 # indexing starts at 1 as in Fortran
-    translate = { 'down_up':'edge_down_up', 'left_right':'edge_left_right' } #new name is edge_left_right
-
-    def __init__(self,gridfilename):
-        try:
-            self.nc = cdf.Dataset(gridfilename, "r")
-            print('#NC4: Opening file:',gridfilename)
-        except RuntimeError:
-            print """ Unable to open grid file %s, maybe you forgot to download it ?
-            To do so, go to the 'Python/' dir and execute './get_MPAS_grids.sh'.
-            """ % gridfile
-            raise
-    def get_pdim(self,comm,name):
-        if name in self.translate: name = self.translate[name]
-        return parallel.PDim(self.nc.dimensions[name], comm)
-    def getvar(self,name):
-        if name in self.translate: name = self.translate[name]
-        return self.nc.variables[name][:]
-    def get_pvar(self,pdim,name):
-        # provides parallel access to a NetCDF variable, with name translation
-        # first deal with special cases
-        if name == 'edge_deg': return parallel.CstPArray1D(pdim, np.int32, 2)
-        # general case
-        if name in self.translate: name = self.translate[name]
-        return parallel.PArray(pdim, self.nc.variables[name])
-    def normalize(self, mesh, radius): pass
 
 #--------------------------------- Main program -----------------------------
@@ -45 +15 @@
 dx,dy=Lx/nx,Ly/ny
 
-filename, llm, nqdyn, g, f, radius = 'cart128.nc', 1, 1, 1., 1., None
+filename, llm, nqdyn, g, f, radius = 'cart_128.nc', 1, 1, 1., 1., None
 unst.setvar('g',g)
 
 print 'Reading Cartesian mesh ...'
-meshfile = DYNAMICO_Format(filename)
-#print("#NC4: check point: 0")
-#print(meshfile)
+meshfile = meshes.DYNAMICO_Format(filename)
 
 def coriolis(lon,lat):
@@ -58 +26 @@
 mesh=meshes.Unstructured_Mesh(meshfile, llm, nqdyn, radius, coriolis)
 caldyn = unst.Caldyn_RSW(mesh)
-#print(caldyn.__dict__.keys())
 
 xx = (mesh.lat_i-nx/2.)*dx
@@ -67 +34 @@
 h0 = 1+0.1*(np.exp(-2.*(x1*x1+yy*yy))+
             np.exp(-2.*(x2*x2+yy*yy)))
-flow0=meshes.asnum([h0,u0])
-
-#print('h0 : ', h0[1234]-1.)
-#print('h0 : ', h0[2345]-1.)
+flow0=prec.asnum([h0,u0])
 
 cfl = .8
@@ -91 +55 @@
 for i in range(10):
     h,u=flow
-    flow, fast, slow = caldyn.bwd_fast_slow(flow, unst.zero_num)
+    flow, fast, slow = caldyn.bwd_fast_slow(flow, prec.zero)
     junk, du_fast = fast
     dh, du_slow = slow