Changeset 794 for codes

Timestamp:

12/13/18 18:01:37 (6 years ago)

Author:

jisesh

Message:

devel/unstructured : moved Davies relaxation zone outside physical domain

Location:

codes/icosagcm/devel/Python

Files:

• dynamico/LAM.py (modified) (1 diff)
• test/py/Baroclinic_3D_ullrich.py (modified) (5 diffs)

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

@@ -10 +10 @@
         N1, N2 = self.N1, self.N2
         m = zeros(c.size)
-        c1,c2 = c0-N1*delta, c0-(N1+N2)*delta
+        c1 = c0+N2*delta
         for i in range(c.size):
             ci=c[i]
-            m[i] = (1.+cos((ci-c2)*pi/(N2*delta)))/2.0
-            if ci < c2 : m[i]=1.
+            m[i] = (1.+cos((ci-c0)*pi/(N2*delta)))/2.0
+            if ci < c0 : m[i]=1.
             if ci > c1 : m[i]=0.
         return m

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

@@ -36 +36 @@
 # Baroclinic instability test based on Ullrich et al. 2015, QJRMS
 
-
-def baroclinic_3D(Lx,nx,Ly,ny,llm,ztop=25000.):
+def create_pmesh(nx,ny):
+    filename = 'cart_%03d_%03d.nc'%(nx,ny)
+    logging.info('Reading Cartesian mesh ...')
+    meshfile = meshes.DYNAMICO_Format(filename)
+    pmesh = meshes.Unstructured_PMesh(comm,meshfile)
+    pmesh.partition_curvilinear(args.mpi_ni,args.mpi_nj)
+    return pmesh
+
+def baroclinic_3D(pmesh, Lx,Ly,llm,ztop=25000.):
     Rd    = 287.0      # Gas constant for dryy air (j kg^-1 K^-1)
     T0    = 288.0      # Reference temperature (K) 
@@ -75 +82 @@
     def coriolis(x,y): return f0+beta0*y # here y is 0 at domain center
 
-    filename = 'cart_%03d_%03d.nc'%(nx,ny)
-    logging.info('Reading Cartesian mesh ...')
-    meshfile = meshes.DYNAMICO_Format(filename)
-    pmesh = meshes.Unstructured_PMesh(comm,meshfile)
-    pmesh.partition_curvilinear(args.mpi_ni,args.mpi_nj)
     nqdyn, radius = 1, None
-    mesh  = meshes.Local_Mesh(pmesh, llm, nqdyn, radius, coriolis)
+    mesh = meshes.Local_Mesh(pmesh, llm, nqdyn, radius, coriolis)
 
     alpha_k = (np.arange(llm)  +.5)/llm
@@ -158 +160 @@
     g, Lx, Ly = 9.81, 4e7, 6e6
     nx, ny, llm = args.nx, args.ny, args.llm
-    dx,dy = Lx/nx,Ly/ny
+    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)
-
-    thermo, mesh, params, flow0, gas0 = baroclinic_3D(Lx,nx,Ly,ny,llm)
+    
+    pmesh = create_pmesh(nx,ny)
+    thermo, mesh, params, flow0, gas0 = baroclinic_3D(pmesh, Lx,Ly,llm)
 
     mass_bl,mass_dak,mass_dbk = meshes.compute_hybrid_coefs(flow0[0])
@@ -173 +185 @@
     dt = T/nt
     logging.info( 'Time step : %d x %g = %g s' % (nt,dt,nt*dt))
+
     
     caldyn_thermo, caldyn_eta = unst.thermo_entropy, unst.eta_mass
@@ -197 +210 @@
             caldyn_step.mass[:,:], caldyn_step.theta_rhodz[:,:], caldyn_step.u[:,:] = m,S,u
             caldyn_step.geopot[:,:], caldyn_step.W[:,:] = Phi,W
-            kappa = dx**4/432000.
+            #Relaxation zone inside the domain
+            kappa = 3.26e+15 #nonhyd=>#Y3.3#Y3.27#Y3.26#Y3.255#dblchkd;N3.25#N3.2#N3.0#Y3.4#dx**4/432000.
+            #Relaxation zone outside the domain
+            #kappa = 3.6e+15 #nonhyd=>#Y3.3#Y3.27#Y3.26#Y3.255#dblchkd;N3.25#N3.2#N3.0#Y3.4#dx**4/432000.
+            #print('Kappa=',kappa)
             for i in range(nt):
                 caldyn_step.next()