Changeset 869


Ignore:
Timestamp:
05/17/19 13:20:04 (5 years ago)
Author:
dubos
Message:

devel/Python : write cell bounds into mesh file / partition file

Location:
codes/icosagcm/devel/Python
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r838 r869  
    2222INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug 
    2323 
    24 radian=180/math.pi 
     24radian=180/math.pi # convert from radians to degrees 
    2525 
    2626#------------------- Hybrid mass-based coordinate ------------- 
    … …  
    131131        Riv2                           = zeros((nx*ny,4)  ) 
    132132        wee                            = zeros((2*nx*ny,4)) 
     133        primal_bounds_lon, primal_bounds_lat, dual_bounds_lon, dual_bounds_lat = zeros( (nx*ny,4), 4) 
    133134 
    134135        primal_deg,dual_deg                 = indices( nx*ny,   2) 
    … …  
    146147                ij=index(x,y)-1 
    147148                lon_i[ij], lat_i[ij] = x, y 
    148                 put(index(x,y)-1, primal_deg,(primal_edge,primal_vertex,primal_ne), 
    149                     ((indexu(x,y),index(x,y),1),  
    150                     (indexv(x,y),index(x-1,y),1), 
    151                     (indexu(x-1,y),index(x-1,y-1),-1), 
    152                     (indexv(x,y-1),index(x,y-1),-1) )) 
     149                put(ij, primal_deg,(primal_edge,primal_vertex,primal_ne, primal_bounds_lon, primal_bounds_lat), 
     150                    ((indexu(x,y),   index(x,y),      1, x+.5, y+.5),  
     151                     (indexv(x,y),   index(x-1,y),    1, x-.5, y+.5), 
     152                     (indexu(x-1,y), index(x-1,y-1), -1, x-.5, y-.5), 
     153                     (indexv(x,y-1), index(x,y-1),   -1, x+.5, y-.5) )) 
    153154                # dual cells 
    154155                ij=index(x,y)-1 
    155156                lon_v[ij], lat_v[ij] = x+.5, y+.5 
    156                 put(ij,dual_deg,(dual_edge,dual_vertex,dual_ne,Riv2), 
    157                    ((indexv(x+1,y),index(x,y),1,.25), 
    158                     (indexu(x,y+1),index(x+1,y),-1,.25), 
    159                     (indexv(x,y),index(x+1,y+1),-1,.25), 
    160                     (indexu(x,y),index(x,y+1),1,.25) )) 
     157                put(ij,dual_deg,(dual_edge,dual_vertex,dual_ne,Riv2, dual_bounds_lon, dual_bounds_lat), 
     158                    ((indexv(x+1,y),index(x,y),      1, .25, x,    y), 
     159                     (indexu(x,y+1),index(x+1,y),   -1, .25, x+1., y), 
     160                     (indexv(x,y),  index(x+1,y+1), -1, .25, x+1., y+1.), 
     161                     (indexu(x,y),  index(x,y+1),    1, .25, x,    y+1.) )) 
    161162                # edges : 
    162163                # left and right are adjacent primal cells 
    … …  
    193194        self.set_members(locals(), 'primal_num',  'dual_num', 'edge_num', 
    194195                         'primal_deg', 'primal_edge', 'primal_ne', 'primal_vertex', 
     196                         'primal_bounds_lon', 'primal_bounds_lat', 
    195197                         'dual_deg', 'dual_edge', 'dual_ne', 'dual_vertex', 
     198                         'dual_bounds_lon', 'dual_bounds_lat', 
    196199                         'left','right','down','up', 
    197200                         'trisk_deg','trisk','Riv2','wee', 
    … …  
    253256                     [("primal_edge", "i4", self.primal_edge), 
    254257                      ("primal_ne", "i4", self.primal_ne), 
    255                       ("primal_vertex", "i4", self.primal_vertex)] ) 
     258                      ("primal_vertex", "i4", self.primal_vertex), 
     259                      ("primal_bounds_lon", "f8", self.primal_bounds_lon), 
     260                      ("primal_bounds_lat", "f8", self.primal_bounds_lat)] ) 
    256261 
    257262        create_vars( ("dual_cell","dual_edge_or_vertex"),  
    258263                     [("dual_edge", "i4", self.dual_edge), 
    259264                      ("dual_vertex","i4",self.dual_vertex), 
     265                      ("dual_bounds_lon", "f8", self.dual_bounds_lon), 
     266                      ("dual_bounds_lat", "f8", self.dual_bounds_lat), 
    260267                      ("dual_ne","i4",self.dual_ne), 
    261268                      ("Riv2","f8",self.Riv2)] ) 
    … …  
    347354        'lat_e':'latEdge','lon_e':'lonEdge','angle_e':'angleEdge', 
    348355        'wee':'weightsOnEdge','Riv2':'kiteAreasOnVertex', 
    349         'primal_ne':'signOnCell', 'dual_ne':'signOnVertex'} 
     356        'primal_ne':'signOnCell', 'dual_ne':'signOnVertex', 
     357        'primal_bounds_lon':None, 'primal_bounds_lat':None, 
     358        'dual_bounds_lon':None, 'dual_bounds_lat':None } 
    350359    def __init__(self,gridfilename): 
    351360        try: 
    … …  
    377386        # general case 
    378387        name=self.translate[name] 
    379         return parallel.PArray(pdim, self.nc.variables[name]) 
     388        return None if name is None else parallel.PArray(pdim, self.nc.variables[name]) 
    380389    def normalize(self, mesh): 
    381390        (trisk_deg, trisk, Ai, 
    … …  
    397406        # this is now done by apply_map 
    398407        # mesh.le, mesh.de, mesh.Av, mesh.Ai = radius*le, radius*de, r2*Av, r2*Ai 
    399          
     408    def get_bounds(self, mesh): 
     409        def bounds(deg, vertex, lon, lat): 
     410            num, maxdeg = deg.size, deg.max() 
     411            vertex_lon, vertex_lat = np.zeros((num, maxdeg)), np.zeros((num, maxdeg)) 
     412            for cell in range(num): 
     413                for ivertex in range(deg[cell]): 
     414                    vtx = vertex[cell,ivertex] - 1 
     415                    vertex_lon[cell, ivertex] = lon[vtx] 
     416                    vertex_lat[cell, ivertex] = lat[vtx] 
     417            return vertex_lon, vertex_lat 
     418        mesh.primal_bounds_lon, mesh.primal_bounds_lat = bounds(mesh.primal_deg, mesh.primal_vertex, mesh.lon_v, mesh.lat_v) 
     419        mesh.dual_bounds_lon, mesh.dual_bounds_lat = bounds(mesh.dual_deg, mesh.dual_vertex, mesh.lon_i, mesh.lat_i) 
     420 
    400421def compute_ne(num,deg,edges,left,right): 
    401422    ne = np.zeros_like(edges) 
    … …  
    450471        self.lon_v, self.lat_v = mapping.map(lon_v, lat_v) 
    451472        self.lon_e, self.lat_e = mapping.map(lon_e, lat_e) 
     473        self.primal_bounds_lon, self.primal_bounds_lat = mapping.map( self.primal_bounds_lon, self.primal_bounds_lat) 
     474        self.dual_bounds_lon, self.dual_bounds_lat = mapping.map( self.dual_bounds_lon, self.dual_bounds_lat) 
    452475        self.ref_lon_i, self.ref_lat_i = lon_i, lat_i 
    453476        self.ref_lon_e, self.ref_lat_e = lon_e, lat_e 
    … …  
    500523        primal_deg, primal_vertex, primal_edge, primal_ne, lon_i, lat_i, Ai = get_pvars( 
    501524            dim_primal, 'primal_deg', 'primal_vertex', 'primal_edge', 'primal_ne', 'lon_i', 'lat_i', 'Ai' ) 
     525        primal_bounds_lon, primal_bounds_lat = get_pvars(dim_primal, 'primal_bounds_lon', 'primal_bounds_lat') 
    502526        edge_deg, trisk_deg, left_right, down_up, trisk, le, de, wee, lon_e, lat_e, angle_e = get_pvars( 
    503527             dim_edge, 'edge_deg', 'trisk_deg', 'edge_left_right', 'edge_down_up',  
    … …  
    505529        dual_deg, dual_vertex, dual_edge, dual_ne, Riv2, lon_v, lat_v, Av = get_pvars( 
    506530            dim_dual, 'dual_deg', 'dual_vertex', 'dual_edge', 'dual_ne', 'Riv2', 'lon_v', 'lat_v', 'Av') 
     531        dual_bounds_lon, dual_bounds_lat = get_pvars(dim_dual, 'dual_bounds_lon', 'dual_bounds_lat') 
    507532        if gridfile.meshtype == 'curvilinear': 
    508533            self.primal_i, self.primal_j = get_pvars(dim_primal, 'primal_i','primal_j') 
    … …  
    521546        self.edge2edge   = Stencil_glob(trisk_deg, trisk, wee) 
    522547        self.set_members(locals(), 'dim_primal', 'dim_dual', 'dim_edge', 
    523                          'primal_deg', 'primal_vertex', 'primal_edge',  
     548                         'primal_deg', 'primal_vertex', 'primal_edge', 
     549                         'primal_bounds_lon', 'primal_bounds_lat', 
    524550                         'trisk_deg', 'trisk', 'dual_deg', 'dual_edge', 
     551                         'dual_bounds_lon', 'dual_bounds_lat', 
    525552                         'lon_v', 'lat_v', 'Av', 'lon_i', 'lat_i', 'Ai',  
    526553                         'le', 'de', 'lon_e', 'lat_e', 'angle_e') 
    … …  
    640667        # read from mesh file : coordinates, lengths and areas in reference domain 
    641668        lon_i, lat_i, Ai = pmesh.get(get_all_cells, 'lon_i', 'lat_i', 'Ai') 
     669 
     670 
    642671        lon_v, lat_v, Av = pmesh.get(get_all_duals, 'lon_v', 'lat_v', 'Av') 
    643672        le, de, lon_e, lat_e, angle_e = pmesh.get(get_all_edges, 'le', 'de', 'lon_e', 'lat_e', 'angle_e') 
    … …  
    657686                         'le', 'de', 'le_de', 'lon_e', 'lat_e', 'angle_e') 
    658687 
     688        # obtain bounds of primal and dual cells (for XIOS) 
     689        if pmesh.primal_bounds_lon is not None: 
     690            INFO('Reading cell bounds from mesh file.') 
     691            self.primal_bounds_lon, self.primal_bounds_lat = pmesh.get(get_all_cells, 'primal_bounds_lon', 'primal_bounds_lat') 
     692            self.dual_bounds_lon, self.dual_bounds_lat = pmesh.get(get_all_duals, 'dual_bounds_lon', 'dual_bounds_lat') 
     693        else: 
     694            INFO('Mesh file does not contain information abound cell bounds.') 
     695            pmesh.gridfile.get_bounds(self) 
     696 
    659697        # normalize data read from file depending on file format 
    660698        pmesh.gridfile.normalize(self) 
    … …  
    662700        self.apply_map(mapping) 
    663701        # now latitudes and longitudes refer to the physical domain 
     702 
    664703        self.fv = mapping.coriolis(self.lon_v,self.lat_v) # Coriolis parameter 
    665704        # propagate info to Fortran side 

  • codes/icosagcm/devel/Python/test/py/dump_partition.py

    r825 r869  
    107107    def create_vars(self,f,dimname, info): 
    108108        for vname, vtype, vdata in info: 
     109            vdata = np.asarray(vdata) 
     110            print('creating variable %s with shape %s.' %(vname, vdata.shape)) 
    109111            var = f.createVariable(vname,vtype,dimname) 
    110112            var[:] = vdata 
    111113     
    112114    def meshwrite(self,name): 
     115        mesh = self._mesh 
    113116        """ Writing Mesh information to disk for Fortran segment""" 
    114117     
    … …  
    124127     
    125128        #----------------DEFINING DIMENSIONS-------------------- 
     129        max_primal_deg = mesh.primal_deg.max() 
     130        max_dual_deg = mesh.dual_deg.max() 
     131        max_trisk_deg = mesh.trisk_deg.max() 
     132        def crop(dim, *data) : return [ x[:,0:dim] for x in data ] 
     133        mesh.primal_edge, mesh.primal_ne, mesh.primal_vertex = crop(max_primal_deg, mesh.primal_edge, mesh.primal_ne, mesh.primal_vertex) 
     134        mesh.dual_edge, mesh.dual_ne = crop(max_dual_deg, mesh.dual_edge, mesh.dual_ne) 
     135        mesh.trisk, mesh.wee = crop(max_trisk_deg, mesh.trisk, mesh.wee) 
     136 
    126137        for dimname, dimsize in [("primal_cell", self._mesh.primal_deg.size), 
    127138                                    ("dual_cell", self._mesh.dual_deg.size), 
    … …  
    144155     
    145156        self.create_vars(f,("primal_cell","primal_edge_or_vertex"), 
    146                     [("primal_vertex","i4",self._mesh.primal_vertex), 
    147                     ("primal_edge","i4",self._mesh.primal_edge), 
    148                     ("primal_ne","i4",self._mesh.primal_ne) ]) 
    149      
     157                         [("primal_vertex","i4",self._mesh.primal_vertex), 
     158                          ("primal_bounds_lon","f8",self._mesh.primal_bounds_lon), 
     159                          ("primal_bounds_lat","f8",self._mesh.primal_bounds_lat), 
     160                          ("primal_edge","i4",self._mesh.primal_edge), 
     161                          ("primal_ne","i4",self._mesh.primal_ne) ]) 
     162         
    150163        self.create_vars(f,"dual_cell", 
    151164                    [("dual_deg","i4",self._mesh.dual_deg), 
    … …  
    161174        self.create_vars(f,("dual_cell","dual_edge_or_vertex"), 
    162175                    [("dual_edge","i4",self._mesh.dual_edge), 
    163                     ("dual_vertex","i4",self._mesh.dual_vertex), 
    164                     ("dual_ne","i4",self._mesh.dual_ne), 
    165                     ("Riv2","f8",self._mesh.Riv2)]) 
     176                     ("dual_vertex","i4",self._mesh.dual_vertex), 
     177                     ("dual_bounds_lon","f8",self._mesh.dual_bounds_lon), 
     178                     ("dual_bounds_lat","f8",self._mesh.dual_bounds_lat), 
     179                     ("dual_ne","i4",self._mesh.dual_ne), 
     180                     ("Riv2","f8",self._mesh.Riv2)]) 
    166181     
    167182        self.create_vars(f,"edge", 

  • codes/icosagcm/devel/Python/test/py/write_Cartesian_mesh.py

    r825 r869  
    11from dynamico import getargs 
    2 getargs.add("--nx", type=int, 
    3                     default=64, choices=None, 
    4                     help="number of x points") 
    5 getargs.add("--ny", type=int, 
    6                     default=64, choices=None, 
    7                     help="number of y points") 
    8 getargs.add("--Lx", type=float, 
    9                     default=8., choices=None, 
    10                     help="Lx") 
    11 getargs.add("--Ly", type=float, 
    12                     default=8., choices=None, 
    13                     help="Ly") 
     2getargs.add("--nx", type=int, default=64,  
     3                    help="number of points along the x axis") 
     4getargs.add("--ny", type=int, default=0, 
     5                    help="number of points along the y axis, defaults to nx") 
     6getargs.add("--Lx", type=float, default=0.,  
     7                    help="domain size along x, defaults to nx") 
     8getargs.add("--Ly", type=float, default=0.,  
     9                    help="domain size along y, defaults to ny") 
    1410 
    1511args = getargs.parse() 
    … …  
    2218 
    2319nx, ny, Lx, Ly, llm, nqdyn = args.nx, args.ny,args.Lx, args.Ly, 1, 1 
     20if ny==0 : ny=nx 
     21if Lx==0. : Lx=nx 
     22if Ly==0. : Ly=ny 
    2423 
    2524dx,dy=Lx/nx,Ly/ny 
Note: See TracChangeset for help on using the changeset viewer.