Changeset 800 for codes/icosagcm/devel/Python/test

Timestamp:

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

Author:

jisesh

Message:

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

File:

• codes/icosagcm/devel/Python/test/py/Baroclinic_3D_ullrich.py (modified) (6 diffs)

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)