Changeset 3928

Timestamp:

Nov 7, 2006, 1:13:03 PM (18 years ago)

Author:

ole

Message:

Parallel domains now store only full triangles in sww files.
Still need to remove ghost nodes.

Location:

anuga_core/source

Files:

• anuga/abstract_2d_finite_volumes/domain.py (modified) (6 diffs)
• anuga/abstract_2d_finite_volumes/general_mesh.py (modified) (11 diffs)
• anuga/abstract_2d_finite_volumes/neighbour_mesh.py (modified) (10 diffs)
• anuga/abstract_2d_finite_volumes/quantity.py (modified) (6 diffs)
• anuga/abstract_2d_finite_volumes/test_general_mesh.py (modified) (1 diff)
• anuga/advection/advection.py (modified) (2 diffs)
• anuga/shallow_water/data_manager.py (modified) (3 diffs)
• anuga/shallow_water/shallow_water_domain.py (modified) (16 diffs)
• anuga/shallow_water/test_shallow_water_domain.py (modified) (5 diffs)
• anuga_parallel/parallel_api.py (modified) (1 diff)
• anuga_parallel/parallel_shallow_water.py (modified) (2 diffs)
• anuga_parallel/test_parallel_sw_runup.py (modified) (2 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

@@ -47 +47 @@
                  processor=0,
                  numproc=1,
-                 number_of_full_nodes=0,
-                 number_of_full_triangles=0):
+                 number_of_full_nodes=None,
+                 number_of_full_triangles=None):
 
 
@@ -87 +87 @@
 
         # Initialise underlying mesh structure
-        Mesh.__init__(self, coordinates, triangles, boundary,
-                      tagged_elements, geo_reference, use_inscribed_circle,
+        Mesh.__init__(self, coordinates, triangles,
+                      boundary=boundary,
+                      tagged_elements=tagged_elements,
+                      geo_reference=geo_reference,
+                      use_inscribed_circle=use_inscribed_circle,
+                      number_of_full_nodes=number_of_full_nodes,
+                      number_of_full_triangles=number_of_full_triangles,
                       verbose=verbose)
 
@@ -136 +141 @@
         self.numproc = numproc
 
-        self.number_of_full_nodes=number_of_full_nodes
-        self.number_of_full_triangles=number_of_full_triangles
-        
 
         # Setup Communication Buffers
@@ -156 +158 @@
         # =1 for full
         # =0 for ghost
-        N=self.number_of_elements
+        N = len(self) #number_of_elements
         self.tri_full_flag = ones(N, Int)
         for i in self.ghost_recv_dict.keys():
@@ -206 +208 @@
         #MH310505 To avoid calculating the flux across each edge twice, keep an integer (boolean) array,
         #to be used during the flux calculation
-        N=self.number_of_elements
+        N = len(self) #number_of_triangles
         self.already_computed_flux = zeros((N, 3), Int)
 
@@ -923 +925 @@
         from Numeric import ones, sum, equal, Float
 
-        N = self.number_of_elements
+        N = len(self) #number_of_triangles
         d = len(self.conserved_quantities)
 

anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

@@ -56 +56 @@
     #input
     def __init__(self, coordinates, triangles,
+                 number_of_full_nodes=None,
+                 number_of_full_triangles=None,                 
                  geo_reference=None,
                  verbose=False):
@@ -65 +67 @@
         origin is a 3-tuple consisting of UTM zone, easting and northing.
         If specified coordinates are assumed to be relative to this origin.
+
+
+        number_of_full_nodes and number_of_full_triangles relate to
+        parallelism when each mesh has an extra layer of ghost points and
+        ghost triangles attached to the end of the two arrays.
+        In this case it is usefull to specify the number of real (called full)
+        nodes and triangles. If omitted they will default to all.
         """
 
@@ -71 +80 @@
         self.triangles = array(triangles,Int)
         self.coordinates = array(coordinates,Float)
-        
+
+        # Register number of elements and nodes 
+        self.number_of_triangles = N = self.triangles.shape[0]
+        self.number_of_nodes = self.coordinates.shape[0]        
+
+        
+
+        if number_of_full_nodes is None:
+            self.number_of_full_nodes=self.number_of_nodes
+        else:
+            assert int(number_of_full_nodes)
+            self.number_of_full_nodes=number_of_full_nodes            
+
+
+        if number_of_full_triangles is None:
+            self.number_of_full_triangles=self.number_of_triangles           
+        else:
+            assert int(number_of_full_triangles)            
+            self.number_of_full_triangles=number_of_full_triangles
+        
+        
+        #print self.number_of_full_nodes, self.number_of_nodes
+        #print self.number_of_full_triangles, self.number_of_triangles
+        
+            
 
         # FIXME: this stores a geo_reference, but when coords are returned
@@ -94 +127 @@
         assert max(max(self.triangles)) <= self.coordinates.shape[0], msg
 
-
-        # Register number of elements (N)
-        self.number_of_elements = N = self.triangles.shape[0]
 
         # FIXME: Maybe move to statistics?
@@ -175 +205 @@
 
     def __len__(self):
-        return self.number_of_elements
+        return self.number_of_triangles
+    
 
     def __repr__(self):
@@ -213 +244 @@
         """
 
+        N = self.number_of_full_nodes
+
         if unique is True:
-            V = self.coordinates
+            V = self.coordinates[:N,:]
             if absolute is True:
                 if not self.geo_reference.is_absolute():
@@ -227 +260 @@
             if not self.geo_reference.is_absolute():
             
-                V0 = self.geo_reference.get_absolute(V[:,0:2])
-                V1 = self.geo_reference.get_absolute(V[:,2:4])
-                V2 = self.geo_reference.get_absolute(V[:,4:6])
+                V0 = self.geo_reference.get_absolute(V[:N,0:2])
+                V1 = self.geo_reference.get_absolute(V[:N,2:4])
+                V2 = self.geo_reference.get_absolute(V[:N,4:6])
 
                 # This does double the memory need 
@@ -262 +295 @@
         #FIXME (Ole) - oh yes they should
 
-        N = self.number_of_elements
+        N = self.number_of_triangles
         vertex_coordinates = zeros((N, 6), Float)
 
@@ -279 +312 @@
         """
 
-        if (indices ==  None):
-            indices = range(len(self))  #len(self)=number of elements
+        N = self.number_of_full_triangles
+
+        if indices is None:
+            #indices = range(len(self))  #len(self)=number of elements
+            indices = range(N)
 
         return  take(self.triangles, indices)
+    
 
     #FIXME - merge these two (get_vertices and get_triangles)
@@ -314 +351 @@
         return unique_verts.keys()
 
+
     def build_vertexlist(self):
         """Build vertexlist index by vertex ids and for each entry (point id)
@@ -327 +365 @@
 
         vertexlist = [None]*len(self.coordinates)
-        for i in range(self.number_of_elements):
+        for i in range(self.number_of_triangles):
 
             a = self.triangles[i, 0]

anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -67 +67 @@
                  tagged_elements=None,
                  geo_reference=None,
+                 number_of_full_nodes=None,
+                 number_of_full_triangles=None,
                  use_inscribed_circle=False,
                  verbose=False):
@@ -80 +82 @@
 
         General_mesh.__init__(self, coordinates, triangles,
-                              geo_reference, verbose=verbose)
+                              number_of_full_nodes=\
+                              number_of_full_nodes,
+                              number_of_full_triangles=\
+                              number_of_full_triangles,
+                              geo_reference=geo_reference,
+                              verbose=verbose)
 
         if verbose: print 'Initialising mesh'         
 
-        N = self.number_of_elements
+        N = len(self) #Number_of_triangles
 
         self.use_inscribed_circle = use_inscribed_circle
@@ -181 +188 @@
     def set_to_inscribed_circle(self,safety_factor = 1):
         #FIXME phase out eventually
-        N = self.number_of_elements
+        N = self.number_of_triangles
         V = self.vertex_coordinates
 
@@ -225 +232 @@
         #to left hand side edge (facing neighbouring triangle)
 
-        N = self.number_of_elements
+        N = len(self) #Number_of_triangles
         neighbourdict = {}
         for i in range(N):
@@ -289 +296 @@
         """
 
-        N = self.number_of_elements
+        N = len(self) #Number of triangles
         for i in range(N):
             #Find all neighbouring volumes that are not boundaries
@@ -314 +321 @@
         if boundary is None:
             boundary = {}
-            for vol_id in range(self.number_of_elements):
+            for vol_id in range(len(self)):
                 for edge_id in range(0, 3):
                     if self.neighbours[vol_id, edge_id] < 0:
@@ -330 +337 @@
 
             #Check that all boundary segments are assigned a tag
-            for vol_id in range(self.number_of_elements):
+            for vol_id in range(len(self)):
                 for edge_id in range(0, 3):
                     if self.neighbours[vol_id, edge_id] < 0:
@@ -374 +381 @@
 
                 msg = 'Not all elements exist. '
-                assert max(tagged_elements[tag]) < self.number_of_elements, msg
+                assert max(tagged_elements[tag]) < len(self), msg
         #print "tagged_elements", tagged_elements
         self.tagged_elements = tagged_elements
@@ -584 +591 @@
         from anuga.utilities.numerical_tools import anglediff
 
-        N = self.number_of_elements
+        N = len(self)
         #Get x,y coordinates for all vertices for all triangles
         V = self.get_vertex_coordinates()
@@ -732 +739 @@
         str =  '------------------------------------------------\n'
         str += 'Mesh statistics:\n'
-        str += '  Number of triangles = %d\n' %self.number_of_elements
+        str += '  Number of triangles = %d\n' %len(self)
         str += '  Extent [m]:\n'
         str += '    x in [%f, %f]\n' %(min(x), max(x))

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -29 +29 @@
 
         if vertex_values is None:
-            N = domain.number_of_elements
+            N = len(domain) # number_of_elements
             self.vertex_values = zeros((N, 3), Float)
         else:
@@ -42 +42 @@
                   %N
             msg += 'number of elements in specified domain (%d).'\
-                   %domain.number_of_elements
-
-            assert N == domain.number_of_elements, msg
+                   %len(domain)
+
+            assert N == len(domain), msg
 
         self.domain = domain
@@ -1135 +1135 @@
         """
         integral = 0
-        for k in range(self.domain.number_of_elements):
+        for k in range(len(self.domain)):
             area = self.domain.areas[k]
             qc = self.centroid_values[k]
@@ -1165 +1165 @@
         #flux calculations and forcing functions
 
-        N = domain.number_of_elements
+        N = len(domain) # number_of_triangles
         self.explicit_update = zeros(N, Float )
         self.semi_implicit_update = zeros(N, Float )
@@ -1293 +1293 @@
 
     #Check each triangle
-    for k in range(quantity.domain.number_of_elements):
+    for k in range(len(quantity.domain)):
         #Centroid coordinates
         x, y = quantity.domain.centroid_coordinates[k]
@@ -1389 +1389 @@
     from Numeric import zeros, Float
 
-    N = quantity.domain.number_of_elements
+    N = len(quantity.domain)
 
     beta_w = quantity.domain.beta_w

anuga_core/source/anuga/abstract_2d_finite_volumes/test_general_mesh.py

@@ -28 +28 @@
         domain = General_mesh(points, vertices)
 
-        print domain.__class__
 
         assert allclose(domain.get_vertex_coordinates(unique=True), domain.coordinates)

anuga_core/source/anuga/advection/advection.py

@@ -169 +169 @@
         from anuga.config import max_timestep
 
-        N = self.number_of_elements
+        N = len(self)
 
         neighbours = self.neighbours
@@ -260 +260 @@
         from weave import converters
 
-        N = self.number_of_elements
+        N = len(self)
 
 

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -213 +213 @@
         #Create filename
         self.filename = create_filename(domain.get_datadir(),
-            domain.get_name(), extension)
+                                        domain.get_name(), extension)
 
         #print 'F', self.filename
         self.timestep = 0
-        self.number_of_volumes = len(domain)
         self.domain = domain
+        
+
+
+        # Exclude ghosts in case this is a parallel domain
+        self.number_of_nodes = domain.number_of_full_nodes        
+        self.number_of_volumes = domain.number_of_full_triangles
+        #self.number_of_volumes = len(domain)        
+
+
 
 
@@ -277 +285 @@
             if hasattr(domain, 'texture'):
                 fid.texture = domain.texture
-        #else:
+            #else:
             #    fid.texture = 'None'
 
@@ -291 +299 @@
             if domain.smooth is True:
                 fid.createDimension('number_of_points', len(domain.vertexlist))
+                #fid.createDimension('number_of_points', self.number_of_nodes)
+
+                # FIXME(Ole): This will cause sww files for paralle domains to
+                # have ghost nodes stored (but not used by triangles).
+                # To clean this up, we have to change get_vertex_values and friends in
+                # quantity.py (but I can't be bothered right now)
             else:
                 fid.createDimension('number_of_points', 3*self.number_of_volumes)

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -112 +112 @@
                  processor=0,
                  numproc=1,
-                 number_of_full_nodes=0,
-                 number_of_full_triangles=0):
+                 number_of_full_nodes=None,
+                 number_of_full_triangles=None):
 
 
@@ -746 +746 @@
     import sys
 
-    N = domain.number_of_elements
+    N = len(domain) # number_of_triangles
 
     #Shortcuts
@@ -820 +820 @@
     import sys
 
-    N = domain.number_of_elements
+    N = len(domain) # number_of_triangles
 
     #Shortcuts
@@ -848 +848 @@
     import sys
 
-    N = domain.number_of_elements
+    N = len(domain) # number_of_triangles
 
     #Shortcuts
@@ -960 +960 @@
 
     #Update
-    #FIXME: Modify accroditg to c-version - or discard altogether.
-    for k in range(domain.number_of_elements):
+    #FIXME: Modify according to c-version - or discard altogether.
+    for k in range(len(domain)):
 
         if hc[k] < domain.minimum_allowed_height:
@@ -999 +999 @@
     """
 
-    N = domain.number_of_elements
+    N = len(domain)
     beta_h = domain.beta_h
 
@@ -1069 +1069 @@
     """
 
-    N = domain.number_of_elements
+    N = len(domain) # number_of_triangles
     beta_h = domain.beta_h
 
@@ -1114 +1114 @@
     hvbar = h_limiter(domain)
 
-    for k in range(domain.number_of_elements):
+    for k in range(len(domain)):
         #Compute maximal variation in bed elevation
         #  This quantitiy is
@@ -1402 +1402 @@
     g = domain.g
 
-    for k in range(domain.number_of_elements):
+    for k in range(len(domain)):
         avg_h = sum( h[k,:] )/3
 
@@ -1454 +1454 @@
     ymom_update = domain.quantities['ymomentum'].semi_implicit_update
 
-    N = domain.number_of_elements
+    N = len(domain)
     eps = domain.minimum_allowed_height
     g = domain.g
@@ -1489 +1489 @@
     ymom_update = ymom.semi_implicit_update
 
-    N = domain.number_of_elements
+    N = len(domain)
     eps = domain.minimum_allowed_height
     g = domain.g
@@ -1516 +1516 @@
     ymom_update = ymom.explicit_update
 
-    N = domain.number_of_elements
+    N = len(domain)
     eps = domain.minimum_allowed_height
     g = domain.g
@@ -1543 +1543 @@
     ymom_update = domain.quantities['ymomentum'].semi_implicit_update
 
-    N = domain.number_of_elements
+    N = len(domain) # number_of_triangles
     eps = domain.minimum_allowed_height
     g = domain.g #Not necessary? Why was this added?
@@ -1681 +1681 @@
         ymom_update = domain.quantities['ymomentum'].explicit_update
 
-        N = domain.number_of_elements
+        N = len(domain) # number_of_triangles
         t = domain.time
 
@@ -1767 +1767 @@
     #Computed linear combination between constant stages and and
     #stages parallel to the bed elevation.
-    for k in range(domain.number_of_elements):
+    for k in range(len(domain)):
         #Compute maximal variation in bed elevation
         #  This quantitiy is
@@ -2034 +2034 @@
     from Numeric import zeros, Float
 
-    N = domain.number_of_elements
+    N = len(domain) # number_of_triangles
 
     #Shortcuts

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -1311 +1311 @@
         const = eta_w*rho_a/rho_w
 
-        N = domain.number_of_elements
+        N = len(domain) # number_of_triangles
 
         xc = domain.get_centroid_coordinates()
@@ -1427 +1427 @@
         const = eta_w*rho_a/rho_w
 
-        N = domain.number_of_elements
+        N = len(domain) # number_of_triangles
 
         t = domain.time
@@ -1922 +1922 @@
         #Assert that quantities are conserved
         from Numeric import sum
-        for k in range(mesh.number_of_elements):
+        for k in range(len(mesh)):
             assert allclose (ref_centroid_values[k],
                              sum(stage.vertex_values[k,:])/3)
@@ -1944 +1944 @@
 
         #Assert that all vertex quantities have changed
-        for k in range(mesh.number_of_elements):
+        for k in range(len(mesh)):
             #print ref_vertex_values[k,:], stage.vertex_values[k,:]
             assert not allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])
         #and assert that quantities are still conserved
         from Numeric import sum
-        for k in range(mesh.number_of_elements):
+        for k in range(len(mesh)):
             assert allclose (ref_centroid_values[k],
                              sum(stage.vertex_values[k,:])/3)
@@ -2212 +2212 @@
 
             #print t, volume, initial_volume
-
-
-            #if not allclose (volume, initial_volume):
-            #    print 't==4.05'
-            #    for k in range(domain.number_of_elements):
-            #       pass
-            #       print domain.quantities['stage'].centroid_values[k] -\
-            #             ref_centroid_values[k]
 
             assert allclose (volume, initial_volume)

anuga_core/source/anuga_parallel/parallel_api.py

@@ -171 +171 @@
 
     for p in range(numprocs):
+        N = len(submesh['ghost_nodes'][p])                
         M = len(submesh['ghost_triangles'][p])
-        print 'There are %d ghost triangles on proc %d' %(M, p)
-        N = len(submesh['ghost_nodes'][p])
-        print 'There are %d ghost nodes on proc %d' %(N, p)
+        print 'There are %d ghost nodes and %d ghost triangles on proc %d'\
+              %(N, M, p)
 
 

anuga_core/source/anuga_parallel/parallel_shallow_water.py

@@ -34 +34 @@
                  full_send_dict=None,
                  ghost_recv_dict=None,
-                 number_of_full_nodes=0,
-                 number_of_full_triangles=0):
+                 number_of_full_nodes=None,
+                 number_of_full_triangles=None):
 
         Domain.__init__(self,
@@ -48 +48 @@
                         number_of_full_triangles=number_of_full_triangles)
 
-        N = self.number_of_elements
+        N = len(self) # number_of_triangles
 
 #        self.processor = pypar.rank()

anuga_core/source/anuga_parallel/test_parallel_sw_runup.py

@@ -96 +96 @@
         local_interpolation_points.append(i)
 
-# Hack
+# Hack before we excluded ghosts.
 if numprocs == 2:
     if myid == 0:
         del local_interpolation_points[0]                
         #local_interpolation_points = [1,2,3]
-
 if numprocs == 3:
     if myid == 1:
         del local_interpolation_points[0]
-
-
 if numprocs == 4:
     if myid == 0:
@@ -122 +119 @@
 
 for t in domain.evolve(yieldstep = 0.1, finaltime = 5.0):
-    #domain.write_time()
-
+    domain.write_time()
     
     # Record time series at known points