Changeset 800 for codes

Timestamp:

01/16/19 15:27:39 (5 years ago)

Author:

jisesh

Message:

devel/unstructured : introduced mapping from reference domain to physical domain

Location:

codes/icosagcm/devel/Python

Files:

• dynamico/maps.py (added)
• dynamico/meshes.py (modified) (9 diffs)
• test/py/Baroclinic_3D_ullrich.py (modified) (6 diffs)

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

@@ -315 +315 @@
         # general case
         return parallel.PArray(pdim, self.nc.variables[name])
-    def normalize(self, mesh, radius): pass
+    def normalize(self, mesh): pass
 
 #---------------------- MPAS fully unstructured mesh ------------------------
@@ -364 +364 @@
         name=self.translate[name]
         return parallel.PArray(pdim, self.nc.variables[name])
-    def normalize(self, mesh, radius):
+    def normalize(self, mesh):
         (trisk_deg, trisk, Ai,
          de, le, wee, Av, Riv2) = (mesh.trisk_deg, mesh.trisk, mesh.Ai, 
                                        mesh.de, mesh.le, mesh.wee, mesh.Av, mesh.Riv2)
-        edge_num, dual_num, r2 = de.size, Av.size, radius**2
+        edge_num, dual_num = de.size, Av.size
         # fix normalization of wee and Riv2 weights
         for edge1 in range(edge_num):
@@ -381 +381 @@
             if abs(r-1.)>1e-6 : print 'error with Riv2 at vertex ', ivertex
         # scale lengths and areas
-        mesh.le, mesh.de, mesh.Av, mesh.Ai = radius*le, radius*de, r2*Av, r2*Ai
+        # this is now done by apply_map_factor
+        # mesh.le, mesh.de, mesh.Av, mesh.Ai = radius*le, radius*de, r2*Av, r2*Ai
         
 def compute_ne(num,deg,edges,left,right):
@@ -421 +422 @@
     def plot_e(self,data, **kwargs):
         self.plot(self.triedge,data, **kwargs)
+    def apply_map(self, mapping, le, de, Ai, Av):
+        """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
+        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)
+        self.lon_e, self.lat_e = mapping.map(self.lon_e, self.lat_e)
     
 class Unstructured_Mesh(Abstract_Mesh):
@@ -452 +463 @@
                          'Ai', 'Av', 'fv', 'le', 'de',
                          'trisk_deg', 'trisk', 'wee', 'Riv2', 'dual_ne')
-        gridfile.normalize(self, radius)
+        gridfile.normalize(self)
         max_primal_deg, max_dual_deg, max_trisk_deg = [ x.shape[1] for x in primal_edge, dual_edge, trisk]
         unst.init_mesh(llm, nqdyn, edge_num, primal_num,dual_num,
@@ -542 +553 @@
 
 class Local_Mesh(Abstract_Mesh):
-    def __init__(self, pmesh, llm, nqdyn, radius, f):
+    def __init__(self, pmesh, llm, nqdyn, mapping):
         comm, dim_primal, primal_owner, dim_edge, edge_owner, dim_dual, dual_owner = pmesh.members(
             'comm', 'dim_primal', 'primal_owner', 'dim_edge', 'edge_owner', 'dim_dual', 'dual_owner' )
@@ -554 +565 @@
         get_own_cells, get_all_cells = dim_primal.getter(cells_C0), dim_primal.getter(cells_C1)
         get_all_duals, get_all_edges = dim_dual.getter(vertices_V1), dim_edge.getter(edges_E2)
-        
-        self.lon_i, self.lat_i, self.Ai = pmesh.get(get_all_cells, 'lon_i', 'lat_i', 'Ai')
-        self.lon_v, self.lat_v, self.Av = pmesh.get(get_all_duals, 'lon_v', 'lat_v', 'Av')
-        self.fv = f(self.lon_v,self.lat_v) # Coriolis parameter
-        self.le, self.de, self.lon_e, self.lat_e, self.angle_e = pmesh.get(
+
+        # 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))
 
@@ -567 +582 @@
             self.primal_i, self.primal_j = pmesh.get(get_all_cells, 'primal_i', 'primal_j')
             self.dual_i, self.dual_j = pmesh.get(get_all_duals, 'dual_i', 'dual_j')
-
 
         # construct local stencils from global stencils
@@ -624 +638 @@
 
         # normalize and propagate info to Fortran side
-        pmesh.gridfile.normalize(self,radius)
+        pmesh.gridfile.normalize(self)
         unst.init_mesh(llm, nqdyn, self.edge_num, self.primal_num, self.dual_num,
                   max_trisk_deg, max_primal_deg, max_dual_deg,

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

@@ -34 +34 @@
 from dynamico import DCMIP
 from dynamico import meshes
+from dynamico import maps
 from dynamico import xios
 from dynamico import precision as prec
@@ -55 +56 @@
     return pmesh
 
-def baroclinic_3D(pmesh, Lx,Ly,llm,ztop):
+def baroclinic_3D(pmesh, dx,Lx,Ly,llm,ztop):
     Rd    = 287.0      # Gas constant for dryy air (j kg^-1 K^-1)
     T0    = 288.0      # Reference temperature (K) 
@@ -90 +91 @@
 
     def eta(alpha) : return (1.-(lap*ztop*alpha/T0))**(g/(Rd*lap)) # roughly equispaced levels
-    def coriolis(x,y): return f0+beta0*y # here y is 0 at domain center
-
-    nqdyn, radius = 1, None
-    mesh = meshes.Local_Mesh(pmesh, llm, nqdyn, radius, coriolis)
+    
+    betaplane = maps.BetaPlaneMap(dx, .5*Lx, .5*Ly, f0, beta0, .5*Ly)
+    nqdyn = 1
+    mesh = meshes.Local_Mesh(pmesh, llm, nqdyn, betaplane)
 
     alpha_k = (np.arange(llm)  +.5)/llm
     alpha_l = (np.arange(llm+1)+ 0.)/llm
-    x_ik, alpha_ik = np.meshgrid(mesh.lon_i,       alpha_k, indexing='ij')
-    y_ik, alpha_ik = np.meshgrid(mesh.lat_i+y0, alpha_k, indexing='ij')
-    x_il, alpha_il = np.meshgrid(mesh.lon_i,       alpha_l, indexing='ij')
-    y_il, alpha_il = np.meshgrid(mesh.lat_i+y0, alpha_l, indexing='ij')
-    x_ek, alpha_ek = np.meshgrid(mesh.lon_e,       alpha_k, indexing='ij')
-    y_ek, alpha_ek = np.meshgrid(mesh.lat_e+y0, alpha_k, indexing='ij')
+    x_ik, alpha_ik = np.meshgrid(mesh.lon_i, alpha_k, indexing='ij')
+    y_ik, alpha_ik = np.meshgrid(mesh.lat_i, alpha_k, indexing='ij')
+    x_il, alpha_il = np.meshgrid(mesh.lon_i, alpha_l, indexing='ij')
+    y_il, alpha_il = np.meshgrid(mesh.lat_i, alpha_l, indexing='ij')
+    x_ek, alpha_ek = np.meshgrid(mesh.lon_e, alpha_k, indexing='ij')
+    y_ek, alpha_ek = np.meshgrid(mesh.lat_e, alpha_k, indexing='ij')
 
     print(np.shape(alpha_k),np.shape(alpha_l))
@@ -166 +167 @@
 class myDavies(Davies):
     def mask(self,x,y):
-        logging.debug('davies dy = %f'% dy)
-        return self.mask0(y,.5*Ly,dy)*self.mask0(-y,.5*Ly,dy)
+        logging.debug('davies dy = %f'% dx)
+        return self.mask0(y,Ly,dx)*self.mask0(-y,0.,dx)
         
 with xios.Client() as client: # setup XIOS which creates the DYNAMICO communicator
@@ -178 +179 @@
     nx, ny, llm = args.nx, args.ny, args.llm
     dx = Lx/nx
-    dy = dx
-    
-    if True: # physical domain excludes relaxation zone
-        Dy = Ly + dy*2*(args.Davies_N1+args.Davies_N2)
-        print('Dy=%f, Ly = %f'%(Dy,Ly))
-    else: # physical domain includes relaxation zone
-        Dy=Ly
-    #    print('WARNING: Davies_N1 not equal to Davies_N2; average value will be used')
-
+    
     unst.setvar('g',g)
     
     pmesh = create_pmesh(nx,ny)
-    thermo, mesh, params, flow0, gas0 = baroclinic_3D(pmesh, Lx,Ly,llm, args.ztop)
+    thermo, mesh, params, flow0, gas0 = baroclinic_3D(pmesh, dx,Lx,Ly,llm, args.ztop)
 
     mass_bl,mass_dak,mass_dbk = meshes.compute_hybrid_coefs(flow0[0])
@@ -281 +274 @@
 #            context.send_field_primal('p', drhodz)
 
-            context.send_field_primal('theta', gas.s)
+            context.send_field_primal('theta', thermo.T0*exp(gas.s/thermo.Cpd))
+
             context.send_field_primal('uz', w)
             context.send_field_primal('z', z)