Changeset 771 for codes

Timestamp:

11/08/18 17:29:34 (5 years ago)

Author:

jisesh

Message:

devel/Python : Implemented Davies relaxation in Baroclinic_3D_ullrich.py

File:

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

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

@@ -18 +18 @@
 parser = argparse.ArgumentParser()
 
-parser.add_argument("--mpi_ni", type=int,
-                    default=64, choices=None,
+parser.add_argument("--mpi_ni", type=int, default=1,
                     help="number of x processors")
-parser.add_argument("--mpi_nj", type=int,
-                    default=64, choices=None,
+parser.add_argument("--mpi_nj", type=int, default=1,
                     help="number of y processors")
 parser.add_argument("--T", type=float, default=5.,
                     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)
+
 
 args = parser.parse_args()
@@ -43 +44 @@
 
     omega  = 7.292e-5  # Angular velocity of the Earth (s^-1)    
-    phi0   = 45.*np.pi/180.0 # Reference latitude North pi/4 (deg)
+    phi0   = 90.*np.pi/180.0 # Reference latitude North pi/4 (deg)
     f0     = 2*omega*np.sin(phi0) 
     beta0  = 2*omega*np.cos(phi0)/a
@@ -63 +64 @@
     filename = 'cart_%03d_%03d.nc'%(nx,ny)
     print 'Reading Cartesian mesh ...'
-    def coriolis(lon,lat): return f0+0.*lon
+    def coriolis(x,y): return f0+beta0*(y+.5*Ly)
     meshfile = meshes.DYNAMICO_Format(filename)
     pmesh = meshes.Unstructured_PMesh(comm,meshfile)
@@ -131 +132 @@
     return gas, w, z
 
+class Davies:
+    def __init__(self,N1,N2,x_i,y_i,x_e,y_e):
+        self.N1, self.N2 = N1, N2
+        self.beta_i = self.mask(x_i,y_i)
+        self.beta_e = self.mask(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)
+
+        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.
+        return m
+    def relax(self, llm, step, flow):
+        beta_i, beta_e = self.beta_i, self.beta_e
+        m,S,u,Phi,W=flow
+        for l in range(llm):
+            step.mass[:,l]        = beta_i*step.mass[:,l]        + (1.-beta_i)*m[:,l]
+            step.theta_rhodz[:,l] = beta_i*step.theta_rhodz[:,l] + (1.-beta_i)*S[:,l]
+            step.u[:,l]           = beta_e*step.u[:,l]           + (1.-beta_e)*u[:,l]
+        for l in range(llm+1):
+            step.geopot[:,l]      = beta_i*step.geopot[:,l]      + (1.-beta_i)*Phi[:,l]
+            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)
+        
 with xios.Client() as client: # setup XIOS which creates the DYNAMICO communicator
     comm = client.comm
@@ -165 +199 @@
     else: # time stepping in Fortran
         scheme = time_step.ARK2(None, dt)
-        caldyn_step = unst.caldyn_step_NH(mesh,scheme,nt, caldyn_thermo,caldyn_eta,
-                                      thermo,params,params.g)
+        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)
         def next_flow(m,S,u,Phi,W):
             # junk,fast,slow = caldyn.bwd_fast_slow(flow, 0.)
             caldyn_step.mass[:,:], caldyn_step.theta_rhodz[:,:], caldyn_step.u[:,:] = m,S,u
             caldyn_step.geopot[:,:], caldyn_step.W[:,:] = Phi,W
-            caldyn_step.next()
+            for i in range(nt):
+                caldyn_step.next()
+                davies.relax(llm, caldyn_step, flow0)
             return (caldyn_step.mass.copy(), caldyn_step.theta_rhodz.copy(), caldyn_step.u.copy(),
-                caldyn_step.geopot.copy(), caldyn_step.W.copy())
+                    caldyn_step.geopot.copy(), caldyn_step.W.copy())
 
     m,S,u,Phi,W=flow0
@@ -186 +222 @@
     w=mesh.field_mass()
     z=mesh.field_mass()
-    xx,yy = mesh.lat_i, mesh.lon_i