Changeset 825 for codes/icosagcm/devel/Python/dynamico/meshes.py

Timestamp:

04/30/19 18:14:37 (5 years ago)

Author:

dubos

Message:

devel/Python : moved Fortran bindings and *.pyx to dynamico/dev module + necessary changes to test/py/*.py

File:

• codes/icosagcm/devel/Python/dynamico/meshes.py (modified) (12 diffs)

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

@@ +1 @@
+from __future__ import absolute_import     
+
+from dynamico.dev import unstructured as unst
+
 import time
 import math
@@ -8 +12 @@
 from matplotlib.collections import PatchCollection
 
-import unstructured as unst
-import parallel
-from util import list_stencil, Base_class, inverse_list
-from precision import zeros
-
-import getargs
+from dynamico import parallel
+from dynamico.util import list_stencil, BaseClass, inverse_list
+from dynamico import getargs
+from dynamico.parmetis import partition_mesh
+
 log_master, log_world = getargs.getLogger(__name__)
 INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
@@ -19 +22 @@
 
 radian=180/math.pi
-  
+
 #------------------- Hybrid mass-based coordinate -------------
 
@@ -73 +76 @@
     return [ np.asarray(x) for x in bp_il, mass_dak, mass_dbk ]
 
+#----------------------- Base class for "user" meshes ---------
+
+class BaseMesh(BaseClass):
+    def zeros(self,dims): return np.zeros([n for n in dims if n>1]) # overriden by meshes.dev.DevMesh
+
 #----------------------- Cartesian mesh -----------------------
 
@@ -90 +98 @@
     return z.transpose()
 
-class Cartesian_Mesh(Base_class):
+class Cartesian_Mesh(BaseMesh):
     def __init__(self,nx,ny,llm,nqdyn,Lx,Ly,f):
         dx,dy = Lx/nx, Ly/ny
@@ -273 +281 @@
         f.close()
 
-    def field_theta(self,n=1): return zeros((n,self.nqdyn,self.ny,self.nx,self.llm))
-    def field_mass(self,n=1): return zeros((n,self.ny,self.nx,self.llm))
-    def field_z(self,n=1): return zeros((n,self.ny,self.nx,self.llm))
-    def field_w(self,n=1): return zeros((n,self.ny,self.nx,self.llm+1))
-    def field_u(self,n=1): return zeros((n,self.ny,2*self.nx,self.llm))
-    def field_ps(self,n=1): return zeros((n,self.ny,self.nx))
-    def field_ucomp(self,n=1): return zeros((n,self.ny,self.nx,self.llm))
+    def field_theta(self,n=1): return self.zeros((n,self.nqdyn,self.ny,self.nx,self.llm))
+    def field_mass(self,n=1): return self.zeros((n,self.ny,self.nx,self.llm))
+    def field_z(self,n=1): return self.zeros((n,self.ny,self.nx,self.llm))
+    def field_w(self,n=1): return self.zeros((n,self.ny,self.nx,self.llm+1))
+    def field_u(self,n=1): return self.zeros((n,self.ny,2*self.nx,self.llm))
+    def field_ps(self,n=1): return self.zeros((n,self.ny,self.nx))
+    def field_ucomp(self,n=1): return self.zeros((n,self.ny,self.nx,self.llm))
     def ucomp(self,u): 
         return u[range(0,2*self.nx,2),:] if self.ny==1 else u[:,range(0,2*self.nx,2),:]
@@ -293 +301 @@
 #---------------------- DYNAMICO format for fully unstructured mesh and curvilinear meshes ----------------------
 
-class Mesh_Format(Base_class):
+class Mesh_Format(BaseMesh):
     def getdims(self,*names): return [len(self.nc.dimensions[name]) for name in names]
     def get_pdims(self,comm,*names): return [self.get_pdim(comm,name) for name in names]
@@ -399 +407 @@
     return ne
 
-class Abstract_Mesh(Base_class):
+class Abstract_Mesh(BaseMesh):
     def to_dynamico(self): pass
-    def field_theta(self,n=1): return zeros((n,self.nqdyn,self.primal_num,self.llm))
-    def field_mass(self,n=1):  return zeros((n,self.primal_num,self.llm))
-    def field_z(self,n=1):     return zeros((n,self.dual_num,self.llm))
-    def field_w(self,n=1):     return zeros((n,self.primal_num,self.llm+1))
-    def field_u(self,n=1):     return zeros((n,self.edge_num,self.llm))
-    def field_ps(self,n=1):    return zeros((n,self.primal_num,))
+    def field_theta(self,n=1): return self.zeros((n,self.nqdyn,self.primal_num,self.llm))
+    def field_mass(self,n=1):  return self.zeros((n,self.primal_num,self.llm))
+    def field_z(self,n=1):     return self.zeros((n,self.dual_num,self.llm))
+    def field_w(self,n=1):     return self.zeros((n,self.primal_num,self.llm+1))
+    def field_u(self,n=1):     return self.zeros((n,self.edge_num,self.llm))
+    def field_ps(self,n=1):    return self.zeros((n,self.primal_num,))
     def ucov2D(self, ulon, ulat): 
         return self.de*(ulon*np.cos(self.angle_e)+ulat*np.sin(self.angle_e))
     def ucov3D(self, ulon, ulat):
-        ucov = zeros((self.edge_num,ulon.shape[1]))
+        ucov = self.zeros((self.edge_num,ulon.shape[1]))
         for edge in range(self.edge_num):
             angle=self.angle_e[edge]
@@ -518 +526 @@
     def partition_metis(self):
         INFO('Partitioning with ParMetis...')
-        edge_owner = unst.partition_mesh(self.trisk_deg, self.trisk, self.comm.Get_size())
+        edge_owner = partition_mesh(self.comm, self.trisk_deg, self.trisk)
+#        edge_owner = unst.partition_mesh(self.trisk_deg, self.trisk, self.comm.Get_size())
         self.edge_owner = parallel.LocPArray1D(self.dim_edge, edge_owner)
         primal_owner = partition_from_stencil(self.edge_owner, self.primal_deg, self.primal_edge)
@@ -563 +572 @@
 
 class Local_Mesh(Abstract_Mesh):
+    Halo_Xchange = parallel.Halo_Xchange # overriden in dev.meshes
     def __init__(self, pmesh, llm, nqdyn, mapping):
         comm, dim_primal, primal_owner, dim_edge, edge_owner, dim_dual, dual_owner = pmesh.members(
@@ -620 +630 @@
         # construct own dual mesh for XIOS output
         dual_own_glo = find_my_cells(dual_owner)
-        dual_own_glo = np.asarray(primal_own_glo, dtype=np.int32)
+        dual_own_glo = np.asarray(dual_own_glo, dtype=np.int32)
         dual_own_loc = [dict_V1[i] for i in dual_own_glo]
         
@@ -655 +665 @@
         self.to_dynamico()
         # setup halo transfers - NB : llm must be set before we call set_dynamico_transfer
-        self.com_primal = parallel.Halo_Xchange(
+        self.com_primal = self.Halo_Xchange(
             42, dim_primal, cells_C1, get_all_cells(primal_owner))
-        self.com_edges = parallel.Halo_Xchange(
+        self.com_edges = self.Halo_Xchange(
             73, dim_edge, edges_E2, get_all_edges(edge_owner))
         self.com_primal.set_dynamico_transfer('primal')