Changeset 801 for codes

Timestamp:

01/17/19 10:42:38 (5 years ago)

Author:

dubos

Message:

devel/unstructured : implement reference vs physical mesh for spherical meshes

Location:

codes/icosagcm/devel/Python

Files:

• dynamico/maps.py (modified) (1 diff)
• dynamico/meshes.py (modified) (6 diffs)
• test/py/RSW2_MPAS_W02.py (modified) (2 diffs)

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

@@ -22 +22 @@
         radius : planetary radius
         Omega : rotation rate => f = 2*Omega*sin(lat)"""
+        self.radius, self.Omega = radius, Omega
     def coriolis(self, lon,lat):
        return 2.*self.Omega*np.sin(lat)

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

@@ -381 +381 @@
             if abs(r-1.)>1e-6 : print 'error with Riv2 at vertex ', ivertex
         # scale lengths and areas
-        # this is now done by apply_map_factor
+        # this is now done by apply_map
         # mesh.le, mesh.de, mesh.Av, mesh.Ai = radius*le, radius*de, r2*Av, r2*Ai
         
@@ -422 +422 @@
     def plot_e(self,data, **kwargs):
         self.plot(self.triedge,data, **kwargs)
-    def apply_map(self, mapping, le, de, Ai, Av):
+    def apply_map(self, mapping):
         """Before calling apply_map, lat and lon are coordinates in the reference domain.
         After calling apply_map, lat and lon are coordinates in the physical domain."""
         factor_e = mapping.map_factor(self.lon_e, self.lat_e)
-        factor_i = mapping.map_factor(self.lon_i, self.lat_i)**2
-        factor_v = mapping.map_factor(self.lon_v, self.lat_v)**2
-        self.le, self.de, self.Av, self.Ai = factor_e*le, factor_e*de, Av*factor_v, Ai*factor_i
+        factor2_i = mapping.map_factor(self.lon_i, self.lat_i)**2
+        factor2_v = mapping.map_factor(self.lon_v, self.lat_v)**2
+        self.le, self.de, self.Av, self.Ai = self.le*factor_e, self.de*factor_e, self.Av*factor2_v, self.Ai*factor2_i
         self.lon_i, self.lat_i = mapping.map(self.lon_i, self.lat_i)
         self.lon_v, self.lat_v = mapping.map(self.lon_v, self.lat_v)
@@ -503 +503 @@
         self.edge2edge   = Stencil_glob(trisk_deg, trisk, wee)
         self.set_members(locals(), 'dim_primal', 'dim_dual', 'dim_edge',
-                         'primal_deg', 'primal_vertex', 'primal_edge', 'trisk_deg', 'trisk',
+                         'primal_deg', 'primal_vertex', 'primal_edge', 
+                         'trisk_deg', 'trisk', 'dual_deg', 'dual_edge',
                          'lon_v', 'lat_v', 'Av', 'lon_i', 'lat_i', 'Ai', 
                          'le', 'de', 'lon_e', 'lat_e', 'angle_e')
@@ -566 +567 @@
         get_all_duals, get_all_edges = dim_dual.getter(vertices_V1), dim_edge.getter(edges_E2)
 
-        # read from mesh file coordinates in reference domain
-        self.lon_i, self.lat_i, Ai = pmesh.get(get_all_cells, 'lon_i', 'lat_i', 'Ai')
-        self.lon_v, self.lat_v, Av = pmesh.get(get_all_duals, 'lon_v', 'lat_v', 'Av')
-        le, de, self.lon_e, self.lat_e, self.angle_e = pmesh.get(
-            get_all_edges, 'le', 'de', 'lon_e', 'lat_e', 'angle_e')
-        # map from reference domain to physical domain
-        self.apply_map(mapping, le, de, Ai, Av)
-        # now latitudes and longitudes refer to the physical domain
-        self.fv = mapping.coriolis(self.lon_v,self.lat_v) # Coriolis parameter
         primal_num, dual_num, edge_num = map(len, (cells_C1, vertices_V1, edges_E2))
 
@@ -627 +619 @@
         left, right, down, up = left+1, right+1, down+1, up+1
 
-        # store what we have computed
-        le_de = self.le/self.de
+        # read from mesh file : coordinates, lengths and areas in reference domain
+        lon_i, lat_i, Ai = pmesh.get(get_all_cells, 'lon_i', 'lat_i', 'Ai')
+        lon_v, lat_v, Av = pmesh.get(get_all_duals, 'lon_v', 'lat_v', 'Av')
+        le, de, lon_e, lat_e, angle_e = pmesh.get(get_all_edges, 'le', 'de', 'lon_e', 'lat_e', 'angle_e')
+        le_de = le/de
+
+        # store what we have read/computed
         self.set_members(locals(), 'nqdyn', 'llm', 'primal_num', 'dual_num', 'edge_num',
                          'primal_own_glo', 'primal_own_loc', 'primal_own_deg',
                          'primal_deg', 'primal_vertex', 'displ', 'dual_own_loc', 
-                         'trisk_deg', 'trisk', 'wee', 
+                         'trisk_deg', 'trisk', 'wee',
                          'dual_deg', 'dual_edge', 'dual_ne', 'dual_vertex', 'Riv2',
-                         'left','right','primal_deg','primal_edge','primal_ne','le_de',
-                         'vertices_V1', 'edges_E2')
-
-        # normalize and propagate info to Fortran side
+                         'left','right','primal_deg','primal_edge','primal_ne',
+                         'vertices_V1', 'edges_E2',
+                         'lon_i', 'lat_i', 'Ai',
+                         'lon_v', 'lat_v', 'Av',
+                         'le', 'de', 'le_de', 'lon_e', 'lat_e', 'angle_e')
+
+        # normalize data read from file depending on file format
         pmesh.gridfile.normalize(self)
+        # map from reference domain to physical domain
+        self.apply_map(mapping)
+        # now latitudes and longitudes refer to the physical domain
+        self.fv = mapping.coriolis(self.lon_v,self.lat_v) # Coriolis parameter
+        # propagate info to Fortran side
         unst.init_mesh(llm, nqdyn, self.edge_num, self.primal_num, self.dual_num,
                   max_trisk_deg, max_primal_deg, max_dual_deg,
@@ -644 +649 @@
                   dual_deg,dual_edge,dual_ne,dual_vertex,
                   left,right,down,up,trisk_deg,trisk,
-                  self.Ai,self.Av,self.fv,le_de,Riv2,wee)
+                  self.Ai,self.Av,self.fv,self.le_de,self.Riv2,self.wee)
 
         # setup halo transfers - NB : llm must be set before we call set_dynamico_transfer

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

@@ -11 +11 @@
 from dynamico.meshes import MPAS_Format, Unstructured_PMesh as PMesh, Local_Mesh as Mesh
 from dynamico import time_step
+from dynamico import maps
 print '...Done'
 
@@ -24 +25 @@
 
 print 'Omega, planetary PV', Omega, 2*Omega/gh0
+planet = maps.SphereMap(radius, Omega)
 
-def f(lon,lat): return 2*Omega*np.sin(lat) # Coriolis parameter
 print 'Reading MPAS mesh ...'
 meshfile = MPAS_Format('grids/x1.%d.grid.nc'%grid)
 pmesh = PMesh(comm,meshfile)
 pmesh.partition_metis()
-mesh = Mesh(pmesh, llm, nqdyn, radius, f)
+mesh = Mesh(pmesh, llm, nqdyn, planet)
 print '...Done'
 lon, lat = mesh.lon_i, mesh.lat_i