Changeset 8978

Timestamp:

Sep 17, 2013, 7:08:50 PM (12 years ago)

Author:

steve

Message:

Changing anuga init to allow simple interface

Location:

trunk/anuga_core/source

Files:

• anuga/__init__.py (modified) (1 diff)
• anuga/abstract_2d_finite_volumes/quantity.py (modified) (3 diffs)
• anuga/abstract_2d_finite_volumes/test_quantity.py (modified) (3 diffs)
• anuga/operators/rate_operators.py (modified) (5 diffs)
• anuga/operators/test_rate_operators.py (modified) (1 diff)
• anuga_parallel/__init__.py (modified) (3 diffs)

trunk/anuga_core/source/anuga/__init__.py

@@ -387 +387 @@
 from anuga.utilities.model_tools import Create_culvert_bridge_Operator
 
+from anuga_parallel.parallel_api import distribute
+from anuga_parallel.parallel_api import myid, numprocs, get_processor_name
+from anuga_parallel.parallel_api import send, receive
+from anuga_parallel.parallel_api import pypar_available, barrier, finalize
+
+if pypar_available:
+    #from anuga_parallel.parallel_meshes import parallel_rectangle
+    #from anuga_parallel.parallel_shallow_water import Parallel_domain as Parallel_shallow_water_domain
+    #from anuga_parallel.parallel_advection     import Parallel_domain as Parallel_advection_domain
+    from anuga_parallel.parallel_operator_factory import Inlet_operator, Boyd_box_operator
+
+
+
+

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -1091 +1091 @@
         
         
-    def set_values_from_asc_file(self,
+        
+        
+    def set_values_from_utm_grid_file(self,
                              filename,
                              location='vertices',
@@ -1143 +1145 @@
     
         root = filename[:-4]
-    
-#         # Read Meta data
-#         if verbose: log.critical('Reading METADATA from %s' % (root + '.prj'))
-#     
-#         metadatafile = open(root + '.prj')
-#         metalines = metadatafile.readlines()
-#         metadatafile.close()
-#     
-#         L = metalines[0].strip().split()
-#         assert L[0].strip().lower() == 'projection'
-#         projection = L[1].strip()                   #TEXT
-#     
-#         L = metalines[1].strip().split()
-#         assert L[0].strip().lower() == 'zone'
-#         zone = int(L[1].strip())
-#     
-#         L = metalines[2].strip().split()
-#         assert L[0].strip().lower() == 'datum'
-#         datum = L[1].strip()                        #TEXT
-#     
-#         L = metalines[3].strip().split()
-#         assert L[0].strip().lower() == 'zunits'     #IGNORE
-#         zunits = L[1].strip()                       #TEXT
-#     
-#         L = metalines[4].strip().split()
-#         assert L[0].strip().lower() == 'units'
-#         units = L[1].strip()                        #TEXT
-#     
-#         L = metalines[5].strip().split()
-#         assert L[0].strip().lower() == 'spheroid'   #IGNORE
-#         spheroid = L[1].strip()                     #TEXT
-#     
-#         L = metalines[6].strip().split()
-#         assert L[0].strip().lower() == 'xshift'
-#         false_easting = float(L[1].strip())
-#     
-#         L = metalines[7].strip().split()
-#         assert L[0].strip().lower() == 'yshift'
-#         false_northing = float(L[1].strip())
     
     
@@ -1309 +1272 @@
                 # Brute force
                 self.vertex_values[indices] = values.reshape((-1,3))
-                
+            # Cleanup centroid values
+            self.interpolate()
+            
+    def set_values_from_lat_long_grid_file(self,
+                             filename,
+                             location='vertices',
+                             indices=None,
+                             verbose=False):
+        
+        """Read Digital Elevation model from the following ASCII format (.asc or .grd)
+    
+        Example:
+        ncols         3121
+        nrows         1800
+        xllcorner     722000
+        yllcorner     5893000
+        cellsize      25
+        NODATA_value  -9999
+        138.3698 137.4194 136.5062 135.5558 ..........
+    
+    
+        An accompanying file with same basename but extension .prj must exist
+        and is used to fix the UTM zone, datum, false northings and eastings.
+    
+        The prj format is assumed to be as
+    
+        Projection    UTM
+        Zone          56
+        Datum         WGS84
+        Zunits        NO
+        Units         METERS
+        Spheroid      WGS84
+        Xshift        0.0000000000
+        Yshift        10000000.0000000000
+        Parameters
+        """
+
+
+
+
+        msg = 'Filename must be a text string'
+        assert isinstance(filename, basestring), msg
+        
+        msg = 'Extension should be .grd or asc'
+        assert os.path.splitext(filename)[1] in ['.grd', '.asc'], msg
+        
+
+        msg = "set_values_from_grd_file is only defined for location='vertices' or 'centroids'"
+        assert location in ['vertices', 'centroids'], msg
+
+    
+        root = filename[:-4]
+    
+    
+        #Read DEM data
+        datafile = open(filename)
+    
+        if verbose: log.critical('Reading data from %s' % (filename))
+        
+        lines = datafile.readlines()
+        datafile.close()
+    
+    
+        if verbose: log.critical('Got %d lines' % len(lines))
+    
+        # Parse the line data
+        ncols = int(lines[0].split()[1].strip())
+        nrows = int(lines[1].split()[1].strip())
+
+    
+        # Do cellsize (line 4) before line 2 and 3
+        cellsize = float(lines[4].split()[1].strip())
+    
+        # Checks suggested by Joaquim Luis
+        # Our internal representation of xllcorner
+        # and yllcorner is non-standard.
+        xref = lines[2].split()
+        if xref[0].strip() == 'xllcorner':
+            xllcorner = float(xref[1].strip()) # + 0.5*cellsize # Correct offset
+        elif xref[0].strip() == 'xllcenter':
+            xllcorner = float(xref[1].strip())
+        else:
+            msg = 'Unknown keyword: %s' % xref[0].strip()
+            raise Exception, msg
+    
+        yref = lines[3].split()
+        if yref[0].strip() == 'yllcorner':
+            yllcorner = float(yref[1].strip()) # + 0.5*cellsize # Correct offset
+        elif yref[0].strip() == 'yllcenter':
+            yllcorner = float(yref[1].strip())
+        else:
+            msg = 'Unknown keyword: %s' % yref[0].strip()
+            raise Exception, msg
+    
+        NODATA_value = int(float(lines[5].split()[1].strip()))
+    
+        assert len(lines) == nrows + 6
+    
+  
+        #Store data
+        import numpy
+    
+        datafile = open(filename)
+        Z = numpy.loadtxt(datafile, skiprows=6)
+        datafile.close()
+        
+        #print Z.shape
+        #print Z
+        
+        # For raster data we need to a flip and transpose
+        Z = numpy.flipud(Z)
+
+        # Transpose z to have y coordinates along the first axis and x coordinates
+        # along the second axis
+        Z = Z.transpose()
+    
+        x = num.linspace(xllcorner, xllcorner+cellsize*(ncols-1), ncols)
+        y = num.linspace(yllcorner, yllcorner+cellsize*(nrows-1), nrows)
+        
+        
+        if location == 'centroids':
+            points = self.domain.centroid_coordinates
+        else:
+            points = self.domain.vertex_coordinates
+            
+        from anuga.geospatial_data.geospatial_data import Geospatial_data,  ensure_absolute
+        points = ensure_absolute(points, geo_reference=self.domain.geo_reference)
+            
+#         print numpy.max(points[:,0])
+#         print numpy.min(points[:,0])
+#         print numpy.max(points[:,1])
+#         print numpy.min(points[:,1])
+#         
+#         print numpy.max(x)
+#         print numpy.min(x)
+#         print numpy.max(y)
+#         print numpy.min(y)
+        
+        
+        #print x.shape, x
+        #print y.shape, y
+        
+        
+        
+        from  anuga.fit_interpolate.interpolate2d import interpolate2d 
+        
+        #print points 
+        values = interpolate2d(x, y, Z, points, mode='linear', bounds_error=False)
+        
+        #print values
+
+        # Call underlying method using array values
+        if verbose:
+            log.critical('Applying fitted data to quantity')
+            
+        
+        if location == 'centroids':
+            if indices is None:
+                msg = 'Number of values must match number of elements'
+                #assert values.shape[0] == N, msg
+
+                self.centroid_values[:] = values
+            else:
+                msg = 'Number of values must match number of indices'
+                assert values.shape[0] == indices.shape[0], msg
+
+                # Brute force
+                self.centroid_values[indices] = values
+        else:
+            if indices is None:
+                msg = 'Number of values must match number of elements'
+                #assert values.shape[0] == N, msg
+
+                #print values.shape
+                #print self.vertex_values.shape 
+                self.vertex_values[:] = values.reshape((-1,3))
+            else:
+                msg = 'Number of values must match number of indices'
+                assert values.shape[0] == indices.shape[0], msg
+
+                # Brute force
+                self.vertex_values[indices] = values.reshape((-1,3))
+            # Cleanup centroid values
+            self.interpolate()
+                            
     def get_extremum_index(self, mode=None, indices=None):
         """Return index for maximum or minimum value of quantity (on centroids)

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

@@ -1387 +1387 @@
 
 
-    def test_set_values_from_asc_vertices(self):
-        
-        
-        
-
+    def test_set_values_from_utm_grid_file(self):
+        
+        
 
         x0 = 0.0
@@ -1467 +1465 @@
         #print quantity.centroid_values 
 
-        quantity.set_values_from_asc_file(txt_file,
+        quantity.set_values_from_utm_grid_file(txt_file,
                              location='vertices',
                              indices=None,
@@ -1489 +1487 @@
         #print quantity.vertex_values
         #print quantity.centroid_values 
-        quantity.set_values_from_asc_file(txt_file,
+        quantity.set_values_from_utm_grid_file(txt_file,
                              location='centroids',
                              indices=None,

trunk/anuga_core/source/anuga/operators/rate_operators.py

@@ -17 +17 @@
 
 
+from anuga import Quantity
 from anuga.operators.base_operator import Operator
 from anuga.operators.region import Region
@@ -103 +104 @@
 
             rate = self.get_spatial_rate(x,y,t)
+        elif self.rate_type == 'quantity':
+            if indices is None:
+                rate  = self.rate.centroid_values
+            else:
+                rate = self.rate.centroid_values[indices]
         else:
             rate = self.get_non_spatial_rate(t)
@@ -184 +190 @@
 
     def set_rate(self, rate):
-        """Set rate (function)
+        """Set rate 
         Can change rate while running
-        """
-
-        from anuga.utilities.function_utils import determine_function_type
-
-        # Possible types are 'scalar', 't', 'x,y' and 'x,y,t'
-        self.rate_type = determine_function_type(rate)
+        Can be a scalar, or a function of t or x,y or x,y,t or a quantity
+        """
+
+        # Test if rate is a quantity
+        if isinstance(rate, Quantity):
+            self.rate_type = 'quantity'
+        else:
+            # Possible types are 'scalar', 't', 'x,y' and 'x,y,t'
+            from anuga.utilities.function_utils import determine_function_type
+            self.rate_type = determine_function_type(rate)
+
 
         self.rate = rate
+        
 
         if self.rate_type == 'scalar':
+            self.rate_callable = False
+            self.rate_spatial = False
+        elif self.rate_type == 'quantity':
             self.rate_callable = False
             self.rate_spatial = False
@@ -243 +258 @@
                 fid = self.full_indices
                 return num.sum(self.areas[fid]*rate[fid])*self.factor
+            elif self.rate_type == 'quantity':
+                rate = self.rate.centroid_values # rate is a quantity
+                fid = self.full_indices
+                return num.sum(self.areas[fid]*rate[fid])*self.factor
             else:
                 rate = self.get_non_spatial_rate() # rate is a scalar
@@ -251 +270 @@
                 rate = self.get_spatial_rate() # rate is an array
                 return num.sum(self.areas*rate)*self.factor
+            elif self.rate_type == 'quantity':
+                rate = self.rate.centroid_values # rate is a quantity
+                return num.sum(self.areas*rate)*self.factor                
             else:
                 rate = self.get_non_spatial_rate() # rate is a scalar

trunk/anuga_core/source/anuga/operators/test_rate_operators.py

@@ -635 +635 @@
         assert num.allclose(Q_ex, Q)
 
+
+    def test_rate_operator_rate_quantity(self):
+        from anuga.config import rho_a, rho_w, eta_w
+        from math import pi, cos, sin
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #             bac,     bce,     ecf,     dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+
+        #Flat surface with 1m of water
+        domain.set_quantity('elevation', 0.0)
+        domain.set_quantity('stage', 1.0)
+        domain.set_quantity('friction', 0.0)
+
+        Br = Reflective_boundary(domain)
+        domain.set_boundary({'exterior': Br})
+
+        verbose = False
+
+        if verbose:
+            print domain.quantities['elevation'].centroid_values
+            print domain.quantities['stage'].centroid_values
+            print domain.quantities['xmomentum'].centroid_values
+            print domain.quantities['ymomentum'].centroid_values
+
+        # Apply operator to these triangles
+        indices = [0,1,3]
+        factor = 10.0
+
+
+        from anuga import Quantity
+        rate_Q = Quantity(domain)
+        rate_Q.set_values(1.0)
+
+        operator = Rate_operator(domain, rate=rate_Q, factor=factor, \
+                                 indices=indices)
+
+
+        # Apply Operator
+        domain.timestep = 2.0
+        operator()
+        rate = rate_Q.centroid_values[indices]
+        t = operator.get_time()
+        Q = operator.get_Q()
+
+        rate = rate*factor
+        Q_ex = num.sum(domain.areas[indices]*rate)
+        d = operator.get_timestep()*rate + 1
+
+
+        #print "d"
+        #print d
+        #print Q_ex
+        #print Q
+        stage_ex = num.array([ 1.0,  1.0,   1.0,  1.0])
+        stage_ex[indices] = d
+        
+        verbose = False
+        
+        if verbose:
+            print domain.quantities['elevation'].centroid_values
+            print domain.quantities['stage'].centroid_values
+            print domain.quantities['xmomentum'].centroid_values
+            print domain.quantities['ymomentum'].centroid_values
+
+        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
+        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
+        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
+        assert num.allclose(Q_ex, Q)
+
+
     def test_rate_operator_functions_empty_indices(self):
         from anuga.config import rho_a, rho_w, eta_w

trunk/anuga_core/source/anuga_parallel/__init__.py

@@ -5 +5 @@
 imported from this module
 """
+
+# Lets import the standard anuga interface
+#from anuga import *
+
 
 from parallel_api import distribute
@@ -15 +19 @@
     from parallel_shallow_water import Parallel_domain as Parallel_shallow_water_domain
     from parallel_advection     import Parallel_domain as Parallel_advection_domain
+    from parallel_operator_factory import Inlet_operator, Boyd_box_operator
 else:
     from anuga import rectangular_cross as parallel_rectangle
@@ -24 +29 @@
 
 
-#Added by Petar Milevski 10/09/2013
-from anuga_parallel import distribute, myid, numprocs, finalize
-from anuga_parallel.parallel_operator_factory import Inlet_operator, Boyd_box_operator