Changeset 7452

Timestamp:

Sep 1, 2009, 10:29:19 AM (15 years ago)

Author:

steve

Message:

Committing latest parallel code

Location:

anuga_core/source

Files:

• anuga/interface.py (modified) (2 diffs)
• anuga/shallow_water/shallow_water_domain.py (modified) (9 diffs)
• anuga/shallow_water/test_shallow_water_domain.py (modified) (1 diff)
• anuga/visualiser/offline.py (modified) (3 diffs)
• anuga/visualiser/realtime.py (modified) (3 diffs)
• anuga_parallel/test_parallel_sw_flow.py (modified) (10 diffs)
• pymetis/pymetis/metis.c (modified) (3 diffs)
• pymetis/test_metis.py (modified) (1 diff)

anuga_core/source/anuga/interface.py

@@ -18 +18 @@
 
 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Time_boundary
+from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Transmissive_boundary
+
 from anuga.abstract_2d_finite_volumes.util import file_function
+
+from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
 
 from anuga.shallow_water.data_manager import export_grid, create_sts_boundary
@@ -26 +30 @@
 from anuga.utilities.polygon import Polygon_function
 
+from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
+
+
+#-----------------------------
+# rectangular domains
+#-----------------------------
+def rectangular_cross_domain(*args, **kwargs):
+    points, vertices, boundary = rectangular_cross(*args, **kwargs)
+    return Domain(points, vertices, boundary)
+    
+#----------------------------
+# Create domain from file
+#----------------------------
+def create_domain_from_file(file):
+    return pmesh_to_domain_instance(file,Domain)
 
 #---------------------------

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -557 +557 @@
     # @note 'polyline' may contain multiple sections allowing complex shapes.
     # @note Assume absolute UTM coordinates.
-    def new_get_energy_through_cross_section(self, polyline,
+    def get_energy_through_cross_section(self, polyline,
                                          kind='total',
                                          verbose=False):
@@ -588 +588 @@
         cross_section = Cross_section(self, polyline, verbose)
 
-        return cross_section.get_energy_through_cross_section()
+        return cross_section.get_energy_through_cross_section(kind)
 
 
@@ -598 +598 @@
     # @note 'polyline' may contain multiple sections allowing complex shapes.
     # @note Assume absolute UTM coordinates.
-    def get_energy_through_cross_section(self, polyline,
+    def old_get_energy_through_cross_section(self, polyline,
                                          kind='total',
                                          verbose=False):
@@ -2786 +2786 @@
         self.midpoints = ensure_geospatial(self.midpoints, self.domain.geo_reference)
 
+    ##
+    # @brief set verbose mode
+    def set_verbose(self,verbose=True):
+        """Set verbose mode true or flase
+        """
+
+        self.verbose=verbose
 
     ##
@@ -2820 +2827 @@
     ##
     # @brief calculate current energy flow through cross section
-    def get_energy_through_cross_section(self):
+    def get_energy_through_cross_section(self, kind='total'):
         """Obtain average energy head [m] across specified cross section.
 
@@ -2837 +2844 @@
         """
 
+        from anuga.config import g, epsilon, velocity_protection as h0
+        
         # Get interpolated values
         stage = self.domain.get_quantity('stage')
@@ -2843 +2852 @@
         ymomentum = self.domain.get_quantity('ymomentum')
 
-        w = stage.get_values(interpolation_points=midpoints, use_cache=True)
-        z = elevation.get_values(interpolation_points=midpoints, use_cache=True)
-        uh = xmomentum.get_values(interpolation_points=midpoints,
+        w = stage.get_values(interpolation_points=self.midpoints, use_cache=True)
+        z = elevation.get_values(interpolation_points=self.midpoints, use_cache=True)
+        uh = xmomentum.get_values(interpolation_points=self.midpoints,
                                   use_cache=True)
-        vh = ymomentum.get_values(interpolation_points=midpoints,
+        vh = ymomentum.get_values(interpolation_points=self.midpoints,
                                   use_cache=True)
         h = w-z                # Depth
@@ -2853 +2862 @@
         # Compute total length of polyline for use with weighted averages
         total_line_length = 0.0
-        for segment in segments:
+        for segment in self.segments:
             total_line_length += segment.length
 
@@ -2879 +2888 @@
 
             # Add to weighted average
-            weigth = segments[i].length/total_line_length
+            weigth = self.segments[i].length/total_line_length
             average_energy += segment_energy*weigth
 

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -1644 +1644 @@
                                                              verbose=False)
                 assert num.allclose(Et, w + 0.5*u*u/g)
+
+
+    def test_cross_section_class(self):
+        """test_cross_section_class(self):
+
+        Test that the total and specific energy through a cross section can be
+        correctly obtained at run-time from the ANUGA cross section class.
+
+        This test creates a flat bed with a known flow through it, creates a cross
+        section and tests that the correct flow and energies are calculated
+
+        The specifics are
+        e = -1 m
+        u = 2 m/s
+        h = 2 m
+        w = 3 m (width of channel)
+
+        q = u*h*w = 12 m^3/s
+
+        This run tries it with georeferencing and with elevation = -1
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+        from mesh_factory import rectangular
+
+        # Create basic mesh (20m x 3m)
+        width = 3
+        length = 20
+        t_end = 1
+        points, vertices, boundary = rectangular(length, width, length, width)
+
+        # Create shallow water domain
+        domain = Domain(points, vertices, boundary,
+                        geo_reference=Geo_reference(56, 308500, 6189000))
+
+        domain.default_order = 2
+        domain.set_quantities_to_be_stored(None)
+
+        e = -1.0
+        w = 1.0
+        h = w-e
+        u = 2.0
+        uh = u*h
+
+        Br = Reflective_boundary(domain)       # Side walls
+        Bd = Dirichlet_boundary([w, uh, 0])    # 2 m/s across the 3 m inlet:
+
+        # Initial conditions
+        domain.set_quantity('elevation', e)
+        domain.set_quantity('stage', w)
+        domain.set_quantity('xmomentum', uh)
+        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
+
+        # Interpolation points down the middle
+        I = [[0, width/2.],
+             [length/2., width/2.],
+             [length, width/2.]]
+        interpolation_points = domain.geo_reference.get_absolute(I)
+
+        # Shortcuts to quantites
+        stage = domain.get_quantity('stage')
+        xmomentum = domain.get_quantity('xmomentum')
+        ymomentum = domain.get_quantity('ymomentum')
+
+
+        # Create some cross sections
+        cross_sections = []
+        for i in range(5):
+            x = length/2. + i*0.23674563    # Arbitrary
+            polyline = [[x, 0], [x, width]]
+
+            polyline = domain.geo_reference.get_absolute(polyline)
+        
+            cross_sections.append(Cross_section(domain,polyline))
+
+
+ 
+        for t in domain.evolve(yieldstep=0.1, finaltime=t_end):
+            # Check that quantities are they should be in the interior
+            w_t = stage.get_values(interpolation_points)
+            uh_t = xmomentum.get_values(interpolation_points)
+            vh_t = ymomentum.get_values(interpolation_points)
+
+            assert num.allclose(w_t, w)
+            assert num.allclose(uh_t, uh)
+            assert num.allclose(vh_t, 0.0, atol=1.0e-6)
+
+
+            # Check flows and energies through the middle
+            for cross_section in cross_sections:
+                
+                Q = cross_section.get_flow_through_cross_section()
+
+                assert num.allclose(Q, uh*width)
+
+                Es = cross_section.get_energy_through_cross_section(kind='specific')
+
+                assert num.allclose(Es, h + 0.5*u*u/g)
+
+                Et = cross_section.get_energy_through_cross_section(kind='total')
+                
+                assert num.allclose(Et, w + 0.5*u*u/g)
+
+
+
+
+
 
     def test_another_runup_example(self):

anuga_core/source/anuga/visualiser/offline.py

@@ -1 @@
-from Numeric import array, Float, ravel, zeros
+#from Numeric import array, Float, ravel, zeros
+import numpy as num
 from Scientific.IO.NetCDF import NetCDFFile
 from Tkinter import Button, E, Tk, W, Label, StringVar, Scale, HORIZONTAL
@@ -80 +81 @@
         else:
             N_vert = len(fin.variables[quantityName])
-        x = ravel(array(fin.variables['x'], Float))
-        y = ravel(array(fin.variables['y'], Float))
+        x = num.ravel(num.array(fin.variables['x'], num.float))
+        y = num.ravel(num.array(fin.variables['y'], num.float))
         if dynamic is True:
-            q = array(fin.variables[quantityName][frameNumber], Float)
-        else:
-            q = ravel(array(fin.variables[quantityName], Float))
+            q = num.array(fin.variables[quantityName][frameNumber], num.float)
+        else:
+            q = num.ravel(num.array(fin.variables[quantityName], num.float))
 
         q *= self.height_zScales[quantityName]
@@ -123 +124 @@
         for q in filter(lambda n:n != 'x' and n != 'y' and n != 'z' and n != 'time' and n != 'volumes', fin.variables.keys()):
             if len(fin.variables[q].shape) == 1: # Not a time-varying quantity
-                quantities[q] = ravel(array(fin.variables[q], Float))
+                quantities[q] = num.ravel(num.array(fin.variables[q], num.float))
             else: # Time-varying, get the current timestep data
-                quantities[q] = array(fin.variables[q][self.frameNumber], Float)
+                quantities[q] = num.array(fin.variables[q][self.frameNumber], num.float)
         fin.close()
         return quantities

anuga_core/source/anuga/visualiser/realtime.py

@@ -1 @@
-from Numeric import Float, zeros, shape
+#from  import Float, zeros, shape
+import numpy as num
 from Tkinter import Button, E, Tk, W
 from threading import Event
@@ -47 +48 @@
 
         # Also build vert_index - a list of the x & y values of each vertex
-        self.vert_index = zeros((N_vert,2), Float)
+        self.vert_index = num.zeros((N_vert,2), num.float)
         for n in range(N_tri):
             self.vtk_cells.InsertNextCell(3)
@@ -56 +57 @@
     def update_height_quantity(self, quantityName, dynamic=True):
         N_vert = len(self.source.get_vertex_coordinates())
-        qty_index = zeros(N_vert, Float)
+        qty_index = num.zeros(N_vert, num.float)
         triangles = self.source.get_triangles()
         vertex_values, _ = self.source.get_quantity(quantityName).get_vertex_values(xy=False, smooth=False)

anuga_core/source/anuga_parallel/test_parallel_sw_flow.py

@@ -15 +15 @@
 import os
 import sys
-import pypar
+#import pypar
 import numpy as num
-
-from anuga.utilities.numerical_tools import ensure_numeric
-from anuga.utilities.polygon import is_inside_polygon
 
 from anuga.interface import Domain
@@ -30 +27 @@
 
 
-from anuga_parallel.interface import distribute, myid, numprocs, send, receive
+from anuga_parallel.interface import distribute, myid, numprocs, send, receive, barrier
 
 #--------------------------------------------------------------------------
@@ -38 +35 @@
 finaltime = 6.0
 nprocs = 4
-N = 29
-M = 29
-verbose = True
+N = 25
+M = 25
+verbose =  True
 
 #---------------------------------
@@ -51 +48 @@
 # Setup Test
 ##########################################################################
-def evolution_test(parallel=False, G = None):
+def evolution_test(parallel=False, G = None, seq_interpolation_points=None):
 
     #--------------------------------------------------------------------------
@@ -65 +62 @@
     #--------------------------------------------------------------------------
     if parallel:
-        if myid == 0: print 'DISTRIBUTING PARALLEL DOMAIN'
+        if myid == 0 and verbose : print 'DISTRIBUTING PARALLEL DOMAIN'
         domain = distribute(domain, verbose=False)
 
@@ -91 +88 @@
     #------------------------------------------------------------------------------
     # Find which sub_domain in which the interpolation points are located 
-    #------------------------------------------------------------------------------
-    interpolation_points = [[0.4,0.51], [0.6,0.51], [0.8,0.51], [0.9,0.51]]
+    #
+    # Sometimes the interpolation points sit exactly
+    # between to centroids, so in the parallel run we
+    # reset the interpolation points to the centroids
+    # found in the sequential run
+    #------------------------------------------------------------------------------
+    interpolation_points = [[0.4,0.5], [0.6,0.5], [0.8,0.5], [0.9,0.5]]
+
+    if parallel:
+        interpolation_points = seq_interpolation_points
+
+
     gauge_values = []
     tri_ids = []
@@ -108 +115 @@
             tri_ids.append(-2)
 
+
     if verbose: print 'P%d has points = %s' %(myid, tri_ids)
 
-
-
+    if not parallel:
+        c_coord = domain.get_centroid_coordinates()
+        interpolation_points = []
+        for id in tri_ids:
+            if id<1:
+                print 'ERROR: All interpolation points be within the sequential domain!'
+            interpolation_points.append(c_coord[id,:])
+            
     #------------------------------------------------------------------------------
     # Evolve system through time
@@ -152 +166 @@
     assert_(success)
     
-    return G
+    return G, interpolation_points
 
 # Test an nprocs-way run of the shallow water equations
@@ -182 +196 @@
         # array G
         #------------------------------------------
-        pypar.barrier()
+        barrier()
         if myid == 0 and verbose: print 'SEQUENTIAL START'
 
-        G = evolution_test(parallel=False)
+        G , interpolation_points = evolution_test(parallel=False)
         G = num.array(G,num.float)
 
-        pypar.barrier()
+        barrier()
         
         #------------------------------------------
@@ -196 +210 @@
         if myid ==0 and verbose: print 'PARALLEL START'
 
-        evolution_test(parallel=True, G=G)
+        evolution_test(parallel=True, G=G, seq_interpolation_points = interpolation_points)
         
 

anuga_core/source/pymetis/pymetis/metis.c

@@ -1 +1 @@
 #include <Python.h>
-#include <Numeric/arrayobject.h>
+#include <numpy/arrayobject.h>
 
 /* This must be the same as the metis idxtype */
@@ -50 +50 @@
   int edgecut;
   int numflag = 0; // The metis routine requires an int * for numflag.
-  int dims[1]; // PyArray_FromDimsAndData needs an int[] of array sizes.
+  npy_intp dims[1]; // PyArray_SimpleNewFromData needs an npy_intp[] of array sizes.
 
   PyObject * elements;
@@ -104 +104 @@
 
   dims[0] = ne;
-  epart_pyarr = (PyArrayObject *)PyArray_FromDimsAndData(1, dims, PyArray_INT, (char *)epart);
+  epart_pyarr = (PyArrayObject *)PyArray_SimpleNewFromData(1, dims, PyArray_INT, (void *)epart);
+  //epart_pyarr = (PyArrayObject *)PyArray_FromDimsAndData(1, dims, PyArray_INT, (char *)epart);
   dims[0] = nn;
-  npart_pyarr = (PyArrayObject *)PyArray_FromDimsAndData(1, dims, PyArray_INT, (char *)npart);
+  npart_pyarr = (PyArrayObject *)PyArray_SimpleNewFromData(1, dims, PyArray_INT, (void *)npart);
+  //npart_pyarr = (PyArrayObject *)PyArray_FromDimsAndData(1, dims, PyArray_INT, (char *)npart);
 
   

anuga_core/source/pymetis/test_metis.py

@@ -3 +3 @@
 from sys import path
 
-from Numeric import array, allclose
+from numpy import array, allclose
 
 import unittest