Changeset 774

Timestamp:

11/14/18 16:54:54 (5 years ago)

Author:

jisesh

Message:

devel/Python : Added wind perturbation and logging in Baroclinic_3D_ullrich.py

File:

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

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

@@ -12 +12 @@
 import time
 import argparse
+import logging
 from numpy import pi, log, exp, sin, cos
 
@@ -22 +23 @@
 parser.add_argument("--mpi_nj", type=int, default=1,
                     help="number of y processors")
-parser.add_argument("--T", type=float, default=5.,
+parser.add_argument("--T", type=float, default=3600.,
                     help="Length of time slice in seconds")
-parser.add_argument("--Davies_N1", type=int, default=5)
-parser.add_argument("--Davies_N2", type=int, default=5)
+parser.add_argument("--N", type=int, default=24,
+                    help="Number of time slices")
+parser.add_argument("--log",type=str,default='DEBUG',choices=('DEBUG', 'INFO', 'WARNING', 'ERROR', 'CRITICAL'),
+                    help="Set verbosity level of logging")
+parser.add_argument("--Davies_N1", type=int, default=3)
+parser.add_argument("--Davies_N2", type=int, default=3)
 
 
@@ -42 +47 @@
     Cpd   = 1004.5
     p0    = 1e5
+    up    = 1. # amplitude (m/s)
+    xc,yc,lp = 0.,Ly/2.,600000.
 
     omega  = 7.292e-5  # Angular velocity of the Earth (s^-1)    
@@ -55 +62 @@
 
     def Phi_xyeta(y,eta): return T0*g/lap*(1-eta**(Rd*lap/g)) + Phi_xy(y)*log(eta)*exp(-((log(eta)/b)**2))
-    def ulon(x,y,eta): return -u0*(sin(pi*y/Ly)**2)*log(eta)*(eta**(-log(eta)/b/b))
+    def ulon(x,y,eta):
+        u = -u0*(sin(pi*y/Ly)**2)*log(eta)*(eta**(-log(eta)/(b*b))) 
+#        u = u + up*exp(-(((x-xc)**2+(y-yc)**2)/lp**2))
+        return  u
+
     def tmean(eta) : return T0*eta**(Rd*lap/g)
     def T(y,eta) : return tmean(eta)+(Phi_xy(y)/Rd)*(((2/(b*b))*(log(eta))**2)-1)*exp(-((0.5*log(eta))**2)) 
@@ -61 +72 @@
 
     def eta(alpha) : return (1-(lap*ztop*alpha/(T0)))**(g/(Rd*lap)) # roughly equispaced levels
+    def coriolis(x,y): return f0+beta0*y
 
     filename = 'cart_%03d_%03d.nc'%(nx,ny)
-    print 'Reading Cartesian mesh ...'
-    def coriolis(x,y): return f0+beta0*(y+.5*Ly)
+    logging.info('Reading Cartesian mesh ...')
     meshfile = meshes.DYNAMICO_Format(filename)
     pmesh = meshes.Unstructured_PMesh(comm,meshfile)
@@ -74 +85 @@
     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, 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('----------------')
-    print 'ztop(ptop) according to Eq. 7:', T0/lap*(1.-(ptop/p0)**(Rd*lap/g)) 
+    x_ik, alpha_ik = np.meshgrid(mesh.lon_i,       alpha_k, indexing='ij')
+    y_ik, alpha_ik = np.meshgrid(mesh.lat_i+.5*Ly, 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+.5*Ly, 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+.5*Ly, alpha_k, indexing='ij')
+
+    print('ztop(ptop) according to Eq. 7:', T0/lap*(1.-(ptop/p0)**(Rd*lap/g))) 
     print(np.shape(alpha_k),np.shape(alpha_l))
     thermo = dyn.Ideal_perfect(Cpd, Rd, p0, T0)
@@ -100 +110 @@
     for l in range(llm): 
         mass_ik[:,l]=(Phi_il[:,l+1]-Phi_il[:,l])/(g*gas.v[:,l])
-    Sik, ujk, Wil  = gas.s*mass_ik, mesh.field_u(), mesh.field_w()
+    Sik, Wil  = gas.s*mass_ik, mesh.field_w()
     
     u_ek = mesh.ucov3D(ulon(x_ek, y_ek, eta_ek), 0.*eta_ek)
 
-    print 'ztop (m) = ', Phi_il[0,-1]/g, ztop
+    print('ztop (m) = ', Phi_il[0,-1]/g, ztop)
     ptop = p(eta(1.))
-    print 'ptop (Pa) = ', gas.p[0,-1], ptop
+    print( 'ptop (Pa) = ', gas.p[0,-1], ptop)
 
     params=dyn.Struct()
@@ -116 +126 @@
     # define parameters for lower BC
     pbot = p(eta_il[:,0])
-    print 'min p, T :', pbot.min(), tmean(pbot/p0)
+    print( 'min p, T :', pbot.min(), tmean(pbot/p0))
     gas_bot = thermo.set_pT(pbot, tmean(pbot/p0))
     params.pbot = gas_bot.p
     params.rho_bot = 1e6/gas_bot.v
     
-    return thermo, mesh, params, prec.asnum([mass_ik,Sik,ujk,Phi_il,Wil]), gas
+    return thermo, mesh, params, prec.asnum([mass_ik,Sik,u_ek,Phi_il,Wil]), gas
 
 def diagnose(Phi,S,m,W):
@@ -137 +147 @@
         self.beta_i = self.mask(x_i,y_i)
         self.beta_e = self.mask(x_e,y_e)
+        self.x_e,self.y_e = x_e,y_e
     def mask0(self,c,c0,delta): # 1D building block for Davies relaxation
         N1, N2 = self.N1, self.N2
-        N3=N1+N2
         m = np.zeros(c.size)
-
+        c1,c2 = c0-N1*delta, c0-(N1+N2)*delta
         for i in range(c.size):
             ci=c[i]
-            m[i] = (1.+np.cos((ci-c0+N3*delta)*np.pi/(N2*delta)))/2.0
-            if ci < c0-N3*delta : m[i]=1.
-            if ci > c0-N1*delta : m[i]=0.
+            m[i] = (1.+np.cos((ci-c2)*np.pi/(N2*delta)))/2.0
+            if ci < c2 : m[i]=1.
+            if ci > c1 : m[i]=0.
         return m
     def relax(self, llm, step, flow):
@@ -159 +169 @@
             step.W[:,l]           = beta_i*step.W[:,l]           + (1.-beta_i)*W[:,l]
 
-
 class myDavies(Davies):
     def mask(self,x,y):
-#        return self.mask0(y,Ly,dy)
-        return self.mask0(y,-.5*Ly,dy)*self.mask0(-y,-.5*Ly,dy)
+        logging.debug('davies dy = %f' % dy)
+        return self.mask0(y,.5*Ly,dy)*self.mask0(-y,.5*Ly,dy)
+
+
         
 with xios.Client() as client: # setup XIOS which creates the DYNAMICO communicator
     comm = client.comm
     mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size()
-    print '%d/%d starting'%(mpi_rank,mpi_size)
+
+    #-------------Logging verbosity configuration---------------------------
+    myloglevel = args.log
+    myloglevel = getattr(logging, myloglevel.upper())
+    logging.basicConfig(filename='out.log',filemode='w',level=myloglevel, 
+        format='%(processName)s:%(name)s:%(filename)s:%(module)s:%(funcName)s:%(lineno)d:%(levelname)s:%(message)s' )
+
+    logging.debug('%d/%d starting'%(mpi_rank,mpi_size))
 
     g, Lx, Ly = 9.81, 4e7, 6e6
     nx, ny, llm = 200, 30, 25
-    dx,dy=Lx/nx,Ly/ny
+    dx,dy = Lx/nx,Ly/ny
 
     unst.setvar('g',g)
@@ -179 +197 @@
 
     mass_bl,mass_dak,mass_dbk = meshes.compute_hybrid_coefs(flow0[0])
-    print 'Type of mass_bl, mass_dak, mass_dbk : ', [x.dtype for x in mass_bl, mass_dak, mass_dbk]
+    print( 'Type of mass_bl, mass_dak, mass_dbk : ', [x.dtype for x in mass_bl, mass_dak, mass_dbk])
     unst.ker.dynamico_init_hybrid(mass_bl,mass_dak,mass_dbk)
 
-    T=3600.
-    dt=360.
+    T  = args.T
+    dt = 360.
     dz = flow0[3].max()/(params.g*llm)
     nt = int(math.ceil(T/dt))
     dt = T/nt
-    print 'Time step : %d x %g = %g s' % (nt,dt,nt*dt)
+    logging.info( 'Time step : %d x %g = %g s' % (nt,dt,nt*dt))
     
     caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass
@@ -201 +219 @@
         caldyn_step = unst.caldyn_step_NH(mesh,scheme,1, caldyn_thermo,caldyn_eta,thermo,params,params.g)
         davies = myDavies(args.Davies_N1, args.Davies_N2, mesh.lon_i, mesh.lat_i, mesh.lon_e,mesh.lat_e)
+        print('mask min/max :', davies.beta_i.min(), davies.beta_i.max() )
+        logging.debug('mask min/max :')
+        logging.debug(davies.beta_i.min())
+        logging.debug(davies.beta_i.max())
         def next_flow(m,S,u,Phi,W):
             # junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.)
@@ -211 +233 @@
                     caldyn_step.geopot.copy(), caldyn_step.W.copy())
 
-    m,S,u,Phi,W=flow0
+    m,S,u,Phi,W = flow0
     if caldyn_thermo == unst.thermo_theta:
         s=S/m
@@ -225 +247 @@
 
     #T, nslice, dt = 3600., 1, 3600.
-    Nslice=24
+    Nslice = args.N
 
     with xios.Context_Curvilinear(mesh,1, 24*3600) as context:
@@ -238 +260 @@
 
             # write to disk
+            context.send_field_primal('ps', davies.beta_i)
             context.send_field_primal('temp', gas0.T)
             context.send_field_primal('p', gas.p)
@@ -243 +266 @@
             context.send_field_primal('uz', w)
 
-            print 'ptop, model top (m) :', unst.getvar('ptop'), Phi.max()/unst.getvar('g')
-            #if args.mpi_ni*args.mpi_nj==1: plot()
+            print( 'ptop, model top (m) :', unst.getvar('ptop'), Phi.max()/unst.getvar('g'))
 
             time1, elapsed1 =time.time(), unst.getvar('elapsed')
@@ -250 +272 @@
             time2, elapsed2 =time.time(), unst.getvar('elapsed')
             factor = 1000./nt
-            print 'ms per full time step : ', factor*(time2-time1), factor*(elapsed2-elapsed1)
+            print( 'ms per full time step : ', factor*(time2-time1), factor*(elapsed2-elapsed1))
             factor = 1e9/(4*nt*nx*ny*llm)
-            print 'nanosec per gridpoint per full time step : ', factor*(time2-time1), factor*(elapsed2-elapsed1)
+            print( 'nanosec per gridpoint per full time step : ', factor*(time2-time1), factor*(elapsed2-elapsed1))
 
         context.update_calendar(Nslice+1) # make sure XIOS writes last iteration
-print('************DONE************')
+logging.info('************DONE************')