Changeset 2105

Timestamp:

Dec 1, 2005, 6:04:45 PM (19 years ago)

Author:

steve

Message:

Updating merimbula code

Location:

inundation

Files:

• Merimbula/meri0.xya (added)
• Merimbula/prepare.py (added)
• Merimbula/run_least_squares_merimbula.py (modified) (1 diff)
• Merimbula/run_merimbula_lake.py (modified) (2 diffs)
• parallel/build_submesh.py (modified) (18 diffs)
• pymetis/metis-4.0/Makefile.in (modified) (1 diff)
• pyvolution/data_manager.py (modified) (77 diffs)
• pyvolution/netherlands.py (modified) (4 diffs)
• pyvolution/view_tsh.bat (modified) (1 diff)
• zeus/Merimbula.zpi (modified) (1 diff)
• zeus/Validation.zpi (modified) (1 diff)
• zeus/anuga-workspace.zwi (modified) (1 diff)
• zeus/parallel.zpi (modified) (1 diff)
• zeus/pyvolution.zpi (modified) (2 diffs)

inundation/Merimbula/run_least_squares_merimbula.py

@@ -8 +8 @@
 import least_squares
 
-alpha = 0.5
+alpha = 5
 mesh_file = 'merimbula_10785.tsh'
 point_file = 'meri0.xya'

inundation/Merimbula/run_merimbula_lake.py

@@ -29 +29 @@
 import time
 
+
 #-------
 # Domain
@@ -50 +51 @@
 #---------------------------------------
 #   Tidal cycle recorded at Eden as open
-filename = 'Eden_tide_Sept03.dat'
+filename = 'Eden_tide_Sept03.tms'
 print 'Open sea boundary condition from ',filename
 Bf = File_boundary(filename, domain)

inundation/parallel/build_submesh.py

@@ -51 +51 @@
     submesh = {}
     tsubnodes = concatenate((reshape(arrayrange(nnodes),(nnodes,1)), nodes), 1)
-    
+
     # loop over processors
 
@@ -63 +63 @@
 
         # find the boundary edges on processor p
-        
+
         subboundary = {}
         for k in boundary:
@@ -69 +69 @@
                 subboundary[k]=boundary[k]
         boundary_list.append(subboundary)
-        
+
         # find nodes in processor p
 
@@ -77 +77 @@
             nodemap[t[1]]=1
             nodemap[t[2]]=1
-                   
+
         node_list.append(take(tsubnodes,nonzero(nodemap)))
 
@@ -127 +127 @@
     ncoord = len(mesh.coordinates)
     ntriangles = len(mesh.triangles)
-    
+
     # find the first layer of boundary triangles
 
@@ -166 +166 @@
     nodemap = zeros(ncoord, 'i')
     fullnodes = submesh["full_nodes"][p]
-    
+
     subtriangles = []
     for i in range(len(trianglemap)):
@@ -178 +178 @@
     tsubtriangles = concatenate((trilist, mesh.triangles), 1)
     subtriangles = take(tsubtriangles, nonzero(trianglemap))
-    
+
     # keep a record of the triangle vertices, if they are not already there
 
@@ -188 +188 @@
     tsubnodes = concatenate((nodelist, mesh.coordinates), 1)
     subnodes = take(tsubnodes, nonzero(nodemap))
-     
+
     # clean up before exiting
 
@@ -223 +223 @@
 
     ghost_commun = zeros((len(subtri), 2), Int)
-    
+
     for i in range(len(subtri)):
         global_no = subtri[i][0]
@@ -241 +241 @@
 
         ghost_commun[i] = [global_no, neigh]
-    
+
     return ghost_commun
 
@@ -343 +343 @@
 
         # build the ghost boundary layer
-        
+
         [subnodes, subtri] = ghost_layer(submesh, mesh, p, tupper, tlower)
         ghost_triangles.append(subtri)
         ghost_nodes.append(subnodes)
-    
+
         # build the communication pattern for the ghost nodes
 
@@ -392 +392 @@
 
 def submesh_quantities(submesh, quantities, triangles_per_proc):
-    
+
     nproc = len(triangles_per_proc)
 
@@ -398 +398 @@
 
     # build an empty dictionary to hold the quantites
-    
+
     submesh["full_quan"] = {}
     submesh["ghost_quan"] = {}
@@ -406 +406 @@
 
     # loop trough the subdomains
-    
+
     for p in range(nproc):
         upper =   lower+triangles_per_proc[p]
 
         # find the global ID of the ghost triangles
-        
+
         global_id = []
         M = len(submesh["ghost_triangles"][p])
@@ -419 +419 @@
         # use the global ID to extract the quantites information from
         # the full domain
-        
+
         for k in quantities:
             submesh["full_quan"][k].append(quantities[k][lower:upper])
@@ -456 +456 @@
 
     submeshf = submesh_full(nodes, triangles, edges, triangles_per_proc)
-        
+
     # add any extra ghost boundary layer information
 
@@ -463 +463 @@
     # order the quantities information to be the same as the triangle
     # information
-    
+
     submesh = submesh_quantities(submeshg, quantities, triangles_per_proc)
-    
+
     return submesh
 
@@ -498 +498 @@
         submesh_cell["full_quan"][k] = submesh["full_quan"][k][0]
         submesh_cell["ghost_quan"][k] = submesh["ghost_quan"][k][0]
-        
+
     return submesh_cell
 
-

inundation/pymetis/metis-4.0/Makefile.in

@@ -1 +1 @@
 
 # Which compiler to use
-CC = cc
+CC = gcc
 
 # What optimization level to use

inundation/pyvolution/data_manager.py

@@ -53 +53 @@
      searchsorted, zeros, allclose, around
 
-    
+
 from coordinate_transforms.geo_reference import Geo_reference
 
@@ -255 +255 @@
 
             fid.order = domain.default_order
-            
-            if hasattr(domain, 'texture'):                  
+
+            if hasattr(domain, 'texture'):
                 fid.texture = domain.texture
-            #else:                          
-            #    fid.texture = 'None'       
+            #else:
+            #    fid.texture = 'None'
 
             #Reference point
@@ -283 +283 @@
             if domain.geo_reference is not None:
                 domain.geo_reference.write_NetCDF(fid)
-                
+
             #FIXME: Backwards compatibility
             fid.createVariable('z', self.precision, ('number_of_points',))
@@ -943 +943 @@
             #easting_min=None, easting_max=None,
             #northing_min=None, northing_max=None,
-            stride = 1,           
+            stride = 1,
             attribute_name = 'elevation',
             z_func = None):
@@ -954 +954 @@
          0.00000e+00  1.40000e-02  1.3535000e-01
 
-    
+
 
     Convert to NetCDF pts format which is
@@ -966 +966 @@
     """
 
-    import os 
+    import os
     #from Scientific.IO.NetCDF import NetCDFFile
     from Numeric import Float, arrayrange, concatenate
@@ -984 +984 @@
 
         if i % stride != 0: continue
-        
+
         fields = line.split()
 
         try:
             assert len(fields) == 3
-        except:    
-            print 'WARNING: Line %d doesn\'t have 3 elements: %s' %(i, line) 
+        except:
+            print 'WARNING: Line %d doesn\'t have 3 elements: %s' %(i, line)
 
         try:
             x = float( fields[0] )
             y = float( fields[1] )
-            z = float( fields[2] )            
+            z = float( fields[2] )
         except:
             continue
@@ -1002 +1002 @@
 
         if callable(z_func):
-            attribute.append(z_func(z))            
-        else:    
+            attribute.append(z_func(z))
+        else:
             attribute.append(z)
 
@@ -1022 +1022 @@
 
     from Numeric import Float
-    
+
     if attribute_name is None:
-        attribute_name = 'attribute'        
-    
+        attribute_name = 'attribute'
+
 
     from Scientific.IO.NetCDF import NetCDFFile
-    
+
     # NetCDF file definition
     outfile = NetCDFFile(ptsname, 'w')
@@ -1037 +1037 @@
     outfile.description = 'NetCDF pts format for compact and '\
                           'portable storage of spatial point data'
-    
+
 
     #Georeferencing
     from coordinate_transforms.geo_reference import Geo_reference
     Geo_reference().write_NetCDF(outfile)
-    
+
 
     outfile.createDimension('number_of_points', len(points))
@@ -1061 +1061 @@
 
     outfile.close()
-    
+
 
 def dem2pts(basename_in, basename_out=None, verbose=False,
@@ -1178 +1178 @@
     for i in range(nrows):
         if verbose and i%((nrows+10)/10)==0:
-            print 'Processing row %d of %d' %(i, nrows)        
-        
+            print 'Processing row %d of %d' %(i, nrows)
+
         lower_index = global_index
         tpoints = zeros((ncols_in_bounding_box, 2), Float)
@@ -1212 +1212 @@
     """Return block of surface lines for each cross section
     Starts with SURFACE LINE,
-    Ends with END CROSS-SECTION  
+    Ends with END CROSS-SECTION
     """
 
@@ -1218 +1218 @@
 
     reading_surface = False
-    for i, line in enumerate(lines): 
-        
+    for i, line in enumerate(lines):
+
         if len(line.strip()) == 0:    #Ignore blanks
-            continue                       
+            continue
 
         if lines[i].strip().startswith('SURFACE LINE'):
@@ -1236 +1236 @@
             easting = float(fields[0])
             northing = float(fields[1])
-            elevation = float(fields[2])                        
-            points.append([easting, northing, elevation])                
+            elevation = float(fields[2])
+            points.append([easting, northing, elevation])
 
 
@@ -1246 +1246 @@
                               verbose=False):
     """Read HEC-RAS Elevation datal from the following ASCII format (.sdf)
-    
+
     Example:
-    
+
 
 # RAS export file created on Mon 15Aug2005 11:42
@@ -1262 +1262 @@
   PROJECTION ZONE: 50
   DATUM: AGD66
-  VERTICAL DATUM: 
-  NUMBER OF REACHES:  19 
-  NUMBER OF CROSS-SECTIONS:  14206 
+  VERTICAL DATUM:
+  NUMBER OF REACHES:  19
+  NUMBER OF CROSS-SECTIONS:  14206
 END HEADER:
 
@@ -1273 +1273 @@
     STREAM ID:Southern-Wungong
     REACH ID:Southern-Wungong
-    STATION:19040.* 
+    STATION:19040.*
     CUT LINE:
-      405548.671603161 , 6438142.7594925 
-      405734.536092045 , 6438326.10404912 
-      405745.130459356 , 6438331.48627354 
-      405813.89633823 , 6438368.6272789 
+      405548.671603161 , 6438142.7594925
+      405734.536092045 , 6438326.10404912
+      405745.130459356 , 6438331.48627354
+      405813.89633823 , 6438368.6272789
     SURFACE LINE:
      405548.67,   6438142.76,   35.37
@@ -1290 +1290 @@
      405566.99,   6438160.83,   35.37
      ...
-   END CROSS-SECTION  
+   END CROSS-SECTION
 
     Convert to NetCDF pts format which is
@@ -1322 +1322 @@
     assert lines[i].strip().upper() == 'BEGIN HEADER:'
     i += 1
-    
+
     assert lines[i].strip().upper().startswith('UNITS:')
     units = lines[i].strip().split()[1]
-    i += 1    
-
-    assert lines[i].strip().upper().startswith('DTM TYPE:')    
     i += 1
 
-    assert lines[i].strip().upper().startswith('DTM:')    
-    i += 1            
-
-    assert lines[i].strip().upper().startswith('STREAM')    
+    assert lines[i].strip().upper().startswith('DTM TYPE:')
     i += 1
 
-    assert lines[i].strip().upper().startswith('CROSS')    
+    assert lines[i].strip().upper().startswith('DTM:')
+    i += 1
+
+    assert lines[i].strip().upper().startswith('STREAM')
+    i += 1
+
+    assert lines[i].strip().upper().startswith('CROSS')
     i += 1
 
@@ -1354 +1354 @@
     i += 1
 
-    assert lines[i].strip().upper().startswith('NUMBER OF REACHES:') 
+    assert lines[i].strip().upper().startswith('NUMBER OF REACHES:')
     i += 1
 
     assert lines[i].strip().upper().startswith('NUMBER OF CROSS-SECTIONS:')
     number_of_cross_sections = int(lines[i].strip().split(':')[1])
-    i += 1                        
+    i += 1
 
 
@@ -1368 +1368 @@
         for k, entry in enumerate(entries):
             points.append(entry[:2])
-            elevation.append(entry[2])            
+            elevation.append(entry[2])
 
 
     msg = 'Actual #number_of_cross_sections == %d, Reported as %d'\
-          %(j+1, number_of_cross_sections)          
+          %(j+1, number_of_cross_sections)
     assert j+1 == number_of_cross_sections, msg
-    
+
     #Get output file
     if basename_out == None:
@@ -1390 +1390 @@
     outfile.description = 'NetCDF pts format for compact and portable ' +\
                           'storage of spatial point data derived from HEC-RAS'
-    
+
     #Georeferencing
     outfile.zone = zone
@@ -1429 +1429 @@
             northing_min = None,
             northing_max = None,
-            expand_search = False, #To avoid intractable situations (This will be fixed when least_squares gets redesigned)           
+            expand_search = False, #To avoid intractable situations (This will be fixed when least_squares gets redesigned)
             verbose = False,
             origin = None,
             datum = 'WGS84',
             format = 'ers'):
-    
+
     """Read SWW file and convert to Digitial Elevation model format (.asc or .ers)
 
@@ -1465 +1465 @@
     The parameter quantity must be the name of an existing quantity or
     an expression involving existing quantities. The default is
-    'elevation'. 
+    'elevation'.
 
     if timestep (an index) is given, output quantity at that timestep
@@ -1471 +1471 @@
     if reduction is given use that to reduce quantity over all timesteps.
 
-    datum 
-    
+    datum
+
     format can be either 'asc' or 'ers'
     """
@@ -1479 +1479 @@
     from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue
     from Numeric import array2string
-    
+
     from utilities.polygon import inside_polygon, outside_polygon, separate_points_by_polygon
     from util import apply_expression_to_dictionary
-    
+
     msg = 'Format must be either asc or ers'
     assert format.lower() in ['asc', 'ers'], msg
 
-        
+
     false_easting = 500000
     false_northing = 10000000
@@ -1498 +1498 @@
     if basename_out is None:
         basename_out = basename_in + '_%s' %quantity
-  
+
     swwfile = basename_in + '.sww'
     demfile = basename_out + '.' + format
@@ -1590 +1590 @@
     #Post condition: Now q has dimension: number_of_points
     assert len(q.shape) == 1
-    assert q.shape[0] == number_of_points    
-    
+    assert q.shape[0] == number_of_points
+
 
     if verbose:
         print 'Processed values for %s are in [%f, %f]' %(quantity, min(q), max(q))
-    
+
 
     #Create grid and update xll/yll corner and x,y
@@ -1609 +1609 @@
     else:
         xmax = easting_max - xllcorner
-        
-    if northing_min is None:        
+
+    if northing_min is None:
         ymin = min(y)
     else:
         ymin = northing_min - yllcorner
-        
-    if northing_max is None:        
+
+    if northing_max is None:
         ymax = max(y)
     else:
         ymax = northing_max - yllcorner
-        
-        
-    
+
+
+
     if verbose: print 'Creating grid'
     ncols = int((xmax-xmin)/cellsize)+1
@@ -1645 +1645 @@
         if format.lower() == 'asc':
             yg = i*cellsize
-        else:    
+        else:
             #this will flip the order of the y values for ers
-            yg = (nrows-i)*cellsize 
-            
+            yg = (nrows-i)*cellsize
+
         for j in xrange(ncols):
             xg = j*cellsize
@@ -1658 +1658 @@
     #Interpolate
     from least_squares import Interpolation
-    
+
 
     #FIXME: This should be done with precrop = True (?), otherwise it'll
@@ -1669 +1669 @@
                            verbose = verbose)
 
-    
+
 
     #Interpolate using quantity values
@@ -1685 +1685 @@
 
     for i in outside_indices:
-        grid_values[i] = NODATA_value        
+        grid_values[i] = NODATA_value
 
 
@@ -1692 +1692 @@
     if format.lower() == 'ers':
         # setup ERS header information
-        grid_values = reshape(grid_values,(nrows, ncols))    
+        grid_values = reshape(grid_values,(nrows, ncols))
         header = {}
         header['datum'] = '"' + datum + '"'
@@ -1717 +1717 @@
         ermapper_grids.write_ermapper_grid(demfile, grid_values, header)
 
-        fid.close()    
+        fid.close()
     else:
         #Write to Ascii format
-        
+
         #Write prj file
-        prjfile = basename_out + '.prj' 
-        
+        prjfile = basename_out + '.prj'
+
         if verbose: print 'Writing %s' %prjfile
         prjid = open(prjfile, 'w')
@@ -1738 +1738 @@
 
 
-    
+
         if verbose: print 'Writing %s' %demfile
- 
+
         ascid = open(demfile, 'w')
 
@@ -1765 +1765 @@
             ascid.write(s[1:-1] + '\n')
 
-            
+
             #print
             #for j in range(ncols):
@@ -1777 +1777 @@
         fid.close()
 
-#Backwards compatibility        
+#Backwards compatibility
 def sww2asc(basename_in, basename_out = None,
             quantity = None,
@@ -1796 +1796 @@
             datum = 'WGS84',
             format = 'asc')
-            
+
 def sww2ers(basename_in, basename_out = None,
             quantity = None,
@@ -1806 +1806 @@
             datum = 'WGS84'):
     print 'sww2ers will soon be obsoleted - please use sww2dem'
-    sww2dem(basename_in, 
+    sww2dem(basename_in,
             basename_out = basename_out,
             quantity = quantity,
@@ -1815 +1815 @@
             origin = origin,
             datum = datum,
-            format = 'ers')         
-################################# END COMPATIBILITY ##############          
+            format = 'ers')
+################################# END COMPATIBILITY ##############
 
 
@@ -1973 +1973 @@
         fields = line.split()
         if verbose and i%((n+10)/10)==0:
-            print 'Processing row %d of %d' %(i, nrows)            
+            print 'Processing row %d of %d' %(i, nrows)
 
         elevation[i, :] = array([float(x) for x in fields])
@@ -2030 +2030 @@
     from Numeric import Float, Int, Int32, searchsorted, zeros, array
     from Numeric import allclose, around
-        
+
     precision = Float
 
@@ -2063 +2063 @@
     assert allclose(longitudes, file_u.variables['LON'])
     assert allclose(longitudes, file_v.variables['LON'])
-    assert allclose(longitudes, e_lon)    
+    assert allclose(longitudes, e_lon)
 
 
@@ -2107 +2107 @@
     zname = 'ELEVATION'
 
-    
+
     amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax]
     uspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lon
@@ -2113 +2113 @@
     elevations = file_e.variables[zname][kmin:kmax, lmin:lmax]
 
-#    if latitudes2[0]==latitudes[0] and latitudes2[-1]==latitudes[-1]:
-#        elevations = file_e.variables['ELEVATION'][kmin:kmax, lmin:lmax]
-#    elif latitudes2[0]==latitudes[-1] and latitudes2[-1]==latitudes[0]:
-#        from Numeric import asarray
-#        elevations=elevations.tolist()
-#        elevations.reverse()
-#        elevations=asarray(elevations)
-#    else:
-#        from Numeric import asarray
-#        elevations=elevations.tolist()
-#        elevations.reverse()
-#        elevations=asarray(elevations)
-#        'print hmmm'
+    #    if latitudes2[0]==latitudes[0] and latitudes2[-1]==latitudes[-1]:
+    #        elevations = file_e.variables['ELEVATION'][kmin:kmax, lmin:lmax]
+    #    elif latitudes2[0]==latitudes[-1] and latitudes2[-1]==latitudes[0]:
+    #        from Numeric import asarray
+    #        elevations=elevations.tolist()
+    #        elevations.reverse()
+    #        elevations=asarray(elevations)
+    #    else:
+    #        from Numeric import asarray
+    #        elevations=elevations.tolist()
+    #        elevations.reverse()
+    #        elevations=asarray(elevations)
+    #        'print hmmm'
 
 
@@ -2359 +2359 @@
     outfile.variables['y'][:] = y - yllcorner
     outfile.variables['z'][:] = z             #FIXME HACK
-    outfile.variables['elevation'][:] = z  
+    outfile.variables['elevation'][:] = z
     outfile.variables['time'][:] = times   #Store time relative
     outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
@@ -2384 +2384 @@
                 i += 1
 
-    #outfile.close()                
+    #outfile.close()
 
     #FIXME: Refactor using code from file_function.statistics
@@ -2413 +2413 @@
 
 
-    outfile.close()                
+    outfile.close()
 
 
@@ -2491 +2491 @@
         fields = line.split(',')
         realtime = calendar.timegm(time.strptime(fields[0], time_format))
-        
+
         T[i] = realtime - starttime
 
@@ -2516 +2516 @@
     #FIXME: Use Georef
     fid.starttime = starttime
-    
-    
+
+
     # dimension definitions
     #fid.createDimension('number_of_volumes', self.number_of_volumes)
@@ -2528 +2528 @@
 
     fid.variables['time'][:] = T
-    
+
     for i in range(Q.shape[1]):
         try:
             name = quantity_names[i]
-        except:    
+        except:
             name = 'Attribute%d'%i
-            
+
         fid.createVariable(name, Float, ('number_of_timesteps',))
-        fid.variables[name][:] = Q[:,i]        
+        fid.variables[name][:] = Q[:,i]
 
     fid.close()
@@ -2642 +2642 @@
 
     if verbose: print '    building domain'
-#    From domain.Domain:
-#    domain = Domain(coordinates, volumes,\
-#                    conserved_quantities = conserved_quantities,\
-#                    other_quantities = other_quantities,zone=zone,\
-#                    xllcorner=xllcorner, yllcorner=yllcorner)
-
-#   From shallow_water.Domain:
+    #    From domain.Domain:
+    #    domain = Domain(coordinates, volumes,\
+    #                    conserved_quantities = conserved_quantities,\
+    #                    other_quantities = other_quantities,zone=zone,\
+    #                    xllcorner=xllcorner, yllcorner=yllcorner)
+
+    #   From shallow_water.Domain:
     coordinates=coordinates.tolist()
     volumes=volumes.tolist()
@@ -2692 +2692 @@
             X = resize(X,(len(X)/3,3))
         domain.set_quantity(quantity,X)
-#
+    #
     for quantity in conserved_quantities:
         try:
@@ -3118 +3118 @@
 
 
-    
+
 def sww2asc_obsolete(basename_in, basename_out = None,
             quantity = None,
@@ -3268 +3268 @@
     #Interpolate
     from least_squares import Interpolation
-    
+
 
     #FIXME: This should be done with precrop = True, otherwise it'll
@@ -3298 +3298 @@
 
 
-    
+
     if verbose: print 'Writing %s' %ascfile
     NODATA_value = -9999
- 
+
     ascid = open(ascfile, 'w')
 
@@ -3390 +3390 @@
     """
     Produce an sww boundary file, from esri ascii data from CSIRO.
-    
+
     Also convert latitude and longitude to UTM. All coordinates are
     assumed to be given in the GDA94 datum.
-    
+
     assume:
     All files are in esri ascii format
@@ -3410 +3410 @@
     """
     from Scientific.IO.NetCDF import NetCDFFile
-        
+
     from coordinate_transforms.redfearn import redfearn
-    
+
     precision = Float # So if we want to change the precision its done here
-    
+
     # go in to the bath dir and load the only file,
     bath_files = os.listdir(bath_dir)
-    #print "bath_files",bath_files 
+    #print "bath_files",bath_files
 
     #fixme: make more general?
@@ -3424 +3424 @@
     bath_metadata,bath_grid =  _read_asc(bath_dir_file)
     #print "bath_metadata",bath_metadata
-    #print "bath_grid",bath_grid 
+    #print "bath_grid",bath_grid
 
     #Use the date.time of the bath file as a basis for
     #the start time for other files
     base_start = bath_file[-12:]
-    
+
     #go into the elevation dir and load the 000 file
     elevation_dir_file = elevation_dir  + os.sep + elevation_prefix \
                          + base_start
     #elevation_metadata,elevation_grid =  _read_asc(elevation_dir_file)
-    #print "elevation_dir_file",elevation_dir_file 
-    #print "elevation_grid", elevation_grid 
-    
+    #print "elevation_dir_file",elevation_dir_file
+    #print "elevation_grid", elevation_grid
+
     elevation_files = os.listdir(elevation_dir)
     ucur_files = os.listdir(ucur_dir)
@@ -3444 +3444 @@
     # file with the same base name as the bath data
     #print "elevation_files[0]",elevation_files[0]
-    #print "'el' + base_start",'el' + base_start 
+    #print "'el' + base_start",'el' + base_start
     assert elevation_files[0] == 'el' + base_start
-    
+
     #assert bath_metadata == elevation_metadata
 
 
 
-    number_of_latitudes = bath_grid.shape[0] 
+    number_of_latitudes = bath_grid.shape[0]
     number_of_longitudes = bath_grid.shape[1]
     #number_of_times = len(os.listdir(elevation_dir))
@@ -3461 +3461 @@
     latitudes = [bath_metadata['yllcorner']+y*bath_metadata['cellsize'] \
                  for y in range(number_of_latitudes)]
-    
+
      # reverse order of lat, so the fist lat represents the first grid row
     latitudes.reverse()
-    
-    #print "latitudes - before _get_min_max_indexes",latitudes 
+
+    #print "latitudes - before _get_min_max_indexes",latitudes
     kmin, kmax, lmin, lmax = _get_min_max_indexes(latitudes,longitudes,
                                                  minlat=minlat, maxlat=maxlat,
                                                  minlon=minlon, maxlon=maxlon)
-    
+
 
     bath_grid = bath_grid[kmin:kmax,lmin:lmax]
@@ -3480 +3480 @@
     number_of_points = number_of_latitudes*number_of_longitudes
     number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
-    #print "number_of_points",number_of_points 
+    #print "number_of_points",number_of_points
 
     #Work out the times
@@ -3492 +3492 @@
     #print "times", times
     #print "number_of_latitudes", number_of_latitudes
-    #print "number_of_longitudes", number_of_longitudes 
-    #print "number_of_times", number_of_times 
+    #print "number_of_longitudes", number_of_longitudes
+    #print "number_of_times", number_of_times
     #print "latitudes", latitudes
     #print "longitudes", longitudes
@@ -3510 +3510 @@
         print '    t in [%f, %f], len(t) == %d'\
               %(min(times), max(times), len(times))
-    
+
     ######### WRITE THE SWW FILE #############
     # NetCDF file definition
@@ -3523 +3523 @@
     outfile.smoothing = 'Yes'
     outfile.order = 1
-    
+
     #Start time in seconds since the epoch (midnight 1/1/1970)
     outfile.starttime = starttime = times[0]
-    
-    
+
+
     # dimension definitions
     outfile.createDimension('number_of_volumes', number_of_volumes)
@@ -3570 +3570 @@
     #Get zone of 1st point.
     refzone, _, _ = redfearn(latitudes[0],longitudes[0])
-    
+
     vertices = {}
     i = 0
-    for k, lat in enumerate(latitudes):       
-        for l, lon in enumerate(longitudes):  
+    for k, lat in enumerate(latitudes):
+        for l, lon in enumerate(longitudes):
 
             vertices[l,k] = i
@@ -3639 +3639 @@
     xmomentum = outfile.variables['xmomentum']
     ymomentum = outfile.variables['ymomentum']
-    
+
     outfile.variables['time'][:] = times   #Store time relative
-    
+
     if verbose: print 'Converting quantities'
-    n = number_of_times   
+    n = number_of_times
     for j in range(number_of_times):
         # load in files
         elevation_meta, elevation_grid = \
            _read_asc(elevation_dir + os.sep + elevation_files[j])
-                         
+
         _, u_momentum_grid =  _read_asc(ucur_dir + os.sep + ucur_files[j])
         _, v_momentum_grid =  _read_asc(vcur_dir + os.sep + vcur_files[j])
 
-        #print "elevation_grid",elevation_grid 
+        #print "elevation_grid",elevation_grid
         #cut matrix to desired size
         elevation_grid = elevation_grid[kmin:kmax,lmin:lmax]
         u_momentum_grid = u_momentum_grid[kmin:kmax,lmin:lmax]
         v_momentum_grid = v_momentum_grid[kmin:kmax,lmin:lmax]
-        #print "elevation_grid",elevation_grid 
+        #print "elevation_grid",elevation_grid
         # handle missing values
         missing = (elevation_grid == elevation_meta['NODATA_value'])
@@ -3681 +3681 @@
                 i += 1
     outfile.close()
-    
+
 def _get_min_max_indexes(latitudes,longitudes,
                         minlat=None, maxlat=None,
@@ -3709 +3709 @@
         if lat_min_index <0:
             lat_min_index = 0
-            
+
 
     if maxlat == None:
@@ -3717 +3717 @@
         if lat_max_index > largest_lat_index:
             lat_max_index = largest_lat_index
-        
+
     if minlon == None:
         lon_min_index = 0
@@ -3724 +3724 @@
         if lon_min_index <0:
             lon_min_index = 0
-        
+
     if maxlon == None:
         lon_max_index = len(longitudes)
@@ -3730 +3730 @@
         lon_max_index = searchsorted(longitudes, maxlon)
 
-    #Take into account that the latitude list was reversed 
+    #Take into account that the latitude list was reversed
     latitudes.reverse() # Python passes by reference, need to swap back
     lat_min_index, lat_max_index = largest_lat_index - lat_max_index , \
@@ -3736 +3736 @@
     lat_max_index = lat_max_index + 1 # taking into account how slicing works
     lon_max_index = lon_max_index + 1 # taking into account how slicing works
-    
+
     return lat_min_index, lat_max_index, lon_min_index, lon_max_index
-        
-    
+
+
 def _read_asc(filename, verbose=False):
     """Read esri file from the following ASCII format (.asc)
@@ -3769 +3769 @@
     cellsize = float(lines.pop(0).split()[1].strip())
     NODATA_value = float(lines.pop(0).split()[1].strip())
-    
-    assert len(lines) == nrows 
-    
+
+    assert len(lines) == nrows
+
     #Store data
     grid = []
-    
+
     n = len(lines)
     for i, line in enumerate(lines):
@@ -3786 +3786 @@
             'cellsize':cellsize,
             'NODATA_value':NODATA_value}, grid
-    
+

inundation/pyvolution/netherlands.py

@@ -90 +90 @@
 print 'Creating domain'
 #Create basic mesh
-points, vertices, boundary = rectangular(N, N)
+points, elements, boundary = rectangular(N, N)
 
 #Create shallow water domain
-domain = Domain(points, vertices, boundary, use_inscribed_circle=True)
+domain = Domain(points, elements, boundary, use_inscribed_circle=True)
 
 domain.check_integrity()
@@ -107 +107 @@
 
 #Set bed-slope and friction
-inflow_stage = 0.1
+inflow_stage = 0.5
 manning = 0.03
 Z = Weir(inflow_stage)
@@ -155 +155 @@
 #
 print 'Initial condition'
-#domain.set_quantity('stage', Constant_height(Z, 0.0))
-domain.set_quantity('stage', Constant_stage(inflow_stage/2.0))
+domain.set_quantity('stage', Constant_height(Z, 0.0))
+#domain.set_quantity('stage', Constant_stage(inflow_stage/2.0))
 
 #Evolve
@@ -168 +168 @@
 #V.update_quantity('elevation')
 
-for t in domain.evolve(yieldstep = 0.02, finaltime = 15.0):
+for t in domain.evolve(yieldstep = 0.005, finaltime = 15.0):
     domain.write_time()
-    domain.write_boundary()    
-    
+    #domain.write_boundary()
+
     print domain.quantities['stage'].get_values(location='centroids',
                                                 indices=[0])

inundation/pyvolution/view_tsh.bat

@@ -1 @@
-python C:\Home\Projects\anuga\pyvolution\view_tsh.py %1  %2
+python C:\Home\Projects\anuga\inundation\pyvolution\view_tsh.py %1  %2

inundation/zeus/Merimbula.zpi

@@ -73 +73 @@
     <ReleaseProjectSaveAll>Off</ReleaseProjectSaveAll>
     <ReleaseProjectReload>Off</ReleaseProjectReload>
+    <file>..\merimbula\prepare.py</file>
     <file>..\merimbula\run_least_squares_merimbula.py</file>
     <file>..\merimbula\run_merimbula_lake.py</file>

inundation/zeus/Validation.zpi

@@ -73 +73 @@
     <ReleaseProjectSaveAll>Off</ReleaseProjectSaveAll>
     <ReleaseProjectReload>Off</ReleaseProjectReload>
-    <file>..\validation\lwru2\create_mesh.py</file>
-    <file>..\validation\lwru2\extract_timeseries.py</file>
-    <file>..\validation\lwru2\filenames.py</file>
-    <file>..\validation\lwru2\lwru2.py</file>
+    <file>..\validation\completed\lwru2\create_mesh.py</file>
+    <file>..\validation\completed\lwru2\extract_timeseries.py</file>
+    <file>..\validation\completed\lwru2\lwru2.py</file>
+    <file>..\validation\completed\lwru2\project.py</file>
     <folder name="Header Files" />
     <folder name="Resource Files" />

inundation/zeus/anuga-workspace.zwi

@@ -2 +2 @@
 <workspace name="anuga-workspace">
     <mode>Debug</mode>
-    <active>analytic_solutions</active>
+    <active>parallel</active>
     <project name="analytic_solutions">analytic_solutions.zpi</project>
     <project name="euler">euler.zpi</project>

inundation/zeus/parallel.zpi

@@ -86 +86 @@
     <file>..\parallel\run_parallel_mesh.py</file>
     <file>..\parallel\run_parallel_sw_merimbula.py</file>
+    <file>..\parallel\run_parallel_sw_merimbula_metis.py</file>
     <file>..\parallel\run_parallel_sw_rectangle.py</file>
     <file>..\parallel\run_sw_lost_mass.py</file>

inundation/zeus/pyvolution.zpi

@@ -79 +79 @@
     <file>..\pyvolution\check_sww_tsh.py</file>
     <file>..\pyvolution\combine_pts.py</file>
-    <file>..\pyvolution\compile.py</file>
     <file>..\pyvolution\config.py</file>
     <file>..\pyvolution\cornell_room.py</file>
     <file>..\pyvolution\data_manager.py</file>
     <file>..\pyvolution\domain.py</file>
-    <file>..\pyvolution\euler.py</file>
-    <file>..\pyvolution\euler_config.py</file>
-    <file>..\pyvolution\euler_ext.c</file>
     <file>..\pyvolution\flatbed.py</file>
     <file>..\pyvolution\flatbed_compare.py</file>
@@ -123 +119 @@
     <file>..\pyvolution\test_data_manager.py</file>
     <file>..\pyvolution\test_domain.py</file>
-    <file>..\pyvolution\test_euler.py</file>
     <file>..\pyvolution\test_general_mesh.py</file>
     <file>..\pyvolution\test_generic_boundary_conditions.py</file>