Changeset 2105


Ignore:
Timestamp:
Dec 1, 2005, 6:04:45 PM (18 years ago)
Author:
steve
Message:

Updating merimbula code

Location:
inundation
Files:
2 added
12 edited

Legend:

Unmodified
Added
Removed
  • inundation/Merimbula/run_least_squares_merimbula.py

    r1292 r2105  
    88import least_squares
    99
    10 alpha = 0.5
     10alpha = 5
    1111mesh_file = 'merimbula_10785.tsh'
    1212point_file = 'meri0.xya'
  • inundation/Merimbula/run_merimbula_lake.py

    r1393 r2105  
    2929import time
    3030
     31
    3132#-------
    3233# Domain
     
    5051#---------------------------------------
    5152#   Tidal cycle recorded at Eden as open
    52 filename = 'Eden_tide_Sept03.dat'
     53filename = 'Eden_tide_Sept03.tms'
    5354print 'Open sea boundary condition from ',filename
    5455Bf = File_boundary(filename, domain)
  • inundation/parallel/build_submesh.py

    r2102 r2105  
    5151    submesh = {}
    5252    tsubnodes = concatenate((reshape(arrayrange(nnodes),(nnodes,1)), nodes), 1)
    53    
     53
    5454    # loop over processors
    5555
     
    6363
    6464        # find the boundary edges on processor p
    65        
     65
    6666        subboundary = {}
    6767        for k in boundary:
     
    6969                subboundary[k]=boundary[k]
    7070        boundary_list.append(subboundary)
    71        
     71
    7272        # find nodes in processor p
    7373
     
    7777            nodemap[t[1]]=1
    7878            nodemap[t[2]]=1
    79                    
     79
    8080        node_list.append(take(tsubnodes,nonzero(nodemap)))
    8181
     
    127127    ncoord = len(mesh.coordinates)
    128128    ntriangles = len(mesh.triangles)
    129    
     129
    130130    # find the first layer of boundary triangles
    131131
     
    166166    nodemap = zeros(ncoord, 'i')
    167167    fullnodes = submesh["full_nodes"][p]
    168    
     168
    169169    subtriangles = []
    170170    for i in range(len(trianglemap)):
     
    178178    tsubtriangles = concatenate((trilist, mesh.triangles), 1)
    179179    subtriangles = take(tsubtriangles, nonzero(trianglemap))
    180    
     180
    181181    # keep a record of the triangle vertices, if they are not already there
    182182
     
    188188    tsubnodes = concatenate((nodelist, mesh.coordinates), 1)
    189189    subnodes = take(tsubnodes, nonzero(nodemap))
    190      
     190
    191191    # clean up before exiting
    192192
     
    223223
    224224    ghost_commun = zeros((len(subtri), 2), Int)
    225    
     225
    226226    for i in range(len(subtri)):
    227227        global_no = subtri[i][0]
     
    241241
    242242        ghost_commun[i] = [global_no, neigh]
    243    
     243
    244244    return ghost_commun
    245245
     
    343343
    344344        # build the ghost boundary layer
    345        
     345
    346346        [subnodes, subtri] = ghost_layer(submesh, mesh, p, tupper, tlower)
    347347        ghost_triangles.append(subtri)
    348348        ghost_nodes.append(subnodes)
    349    
     349
    350350        # build the communication pattern for the ghost nodes
    351351
     
    392392
    393393def submesh_quantities(submesh, quantities, triangles_per_proc):
    394    
     394
    395395    nproc = len(triangles_per_proc)
    396396
     
    398398
    399399    # build an empty dictionary to hold the quantites
    400    
     400
    401401    submesh["full_quan"] = {}
    402402    submesh["ghost_quan"] = {}
     
    406406
    407407    # loop trough the subdomains
    408    
     408
    409409    for p in range(nproc):
    410410        upper =   lower+triangles_per_proc[p]
    411411
    412412        # find the global ID of the ghost triangles
    413        
     413
    414414        global_id = []
    415415        M = len(submesh["ghost_triangles"][p])
     
    419419        # use the global ID to extract the quantites information from
    420420        # the full domain
    421        
     421
    422422        for k in quantities:
    423423            submesh["full_quan"][k].append(quantities[k][lower:upper])
     
    456456
    457457    submeshf = submesh_full(nodes, triangles, edges, triangles_per_proc)
    458        
     458
    459459    # add any extra ghost boundary layer information
    460460
     
    463463    # order the quantities information to be the same as the triangle
    464464    # information
    465    
     465
    466466    submesh = submesh_quantities(submeshg, quantities, triangles_per_proc)
    467    
     467
    468468    return submesh
    469469
     
    498498        submesh_cell["full_quan"][k] = submesh["full_quan"][k][0]
    499499        submesh_cell["ghost_quan"][k] = submesh["ghost_quan"][k][0]
    500        
     500
    501501    return submesh_cell
    502502
    503 
  • inundation/pymetis/metis-4.0/Makefile.in

    r2051 r2105  
    11
    22# Which compiler to use
    3 CC = cc
     3CC = gcc
    44
    55# What optimization level to use
  • inundation/pyvolution/data_manager.py

    r2074 r2105  
    5353     searchsorted, zeros, allclose, around
    5454
    55    
     55
    5656from coordinate_transforms.geo_reference import Geo_reference
    5757
     
    255255
    256256            fid.order = domain.default_order
    257            
    258             if hasattr(domain, 'texture'):                 
     257
     258            if hasattr(domain, 'texture'):
    259259                fid.texture = domain.texture
    260             #else:                         
    261             #    fid.texture = 'None'       
     260            #else:
     261            #    fid.texture = 'None'
    262262
    263263            #Reference point
     
    283283            if domain.geo_reference is not None:
    284284                domain.geo_reference.write_NetCDF(fid)
    285                
     285
    286286            #FIXME: Backwards compatibility
    287287            fid.createVariable('z', self.precision, ('number_of_points',))
     
    943943            #easting_min=None, easting_max=None,
    944944            #northing_min=None, northing_max=None,
    945             stride = 1,           
     945            stride = 1,
    946946            attribute_name = 'elevation',
    947947            z_func = None):
     
    954954         0.00000e+00  1.40000e-02  1.3535000e-01
    955955
    956    
     956
    957957
    958958    Convert to NetCDF pts format which is
     
    966966    """
    967967
    968     import os 
     968    import os
    969969    #from Scientific.IO.NetCDF import NetCDFFile
    970970    from Numeric import Float, arrayrange, concatenate
     
    984984
    985985        if i % stride != 0: continue
    986        
     986
    987987        fields = line.split()
    988988
    989989        try:
    990990            assert len(fields) == 3
    991         except:   
    992             print 'WARNING: Line %d doesn\'t have 3 elements: %s' %(i, line) 
     991        except:
     992            print 'WARNING: Line %d doesn\'t have 3 elements: %s' %(i, line)
    993993
    994994        try:
    995995            x = float( fields[0] )
    996996            y = float( fields[1] )
    997             z = float( fields[2] )           
     997            z = float( fields[2] )
    998998        except:
    999999            continue
     
    10021002
    10031003        if callable(z_func):
    1004             attribute.append(z_func(z))           
    1005         else:   
     1004            attribute.append(z_func(z))
     1005        else:
    10061006            attribute.append(z)
    10071007
     
    10221022
    10231023    from Numeric import Float
    1024    
     1024
    10251025    if attribute_name is None:
    1026         attribute_name = 'attribute'       
    1027    
     1026        attribute_name = 'attribute'
     1027
    10281028
    10291029    from Scientific.IO.NetCDF import NetCDFFile
    1030    
     1030
    10311031    # NetCDF file definition
    10321032    outfile = NetCDFFile(ptsname, 'w')
     
    10371037    outfile.description = 'NetCDF pts format for compact and '\
    10381038                          'portable storage of spatial point data'
    1039    
     1039
    10401040
    10411041    #Georeferencing
    10421042    from coordinate_transforms.geo_reference import Geo_reference
    10431043    Geo_reference().write_NetCDF(outfile)
    1044    
     1044
    10451045
    10461046    outfile.createDimension('number_of_points', len(points))
     
    10611061
    10621062    outfile.close()
    1063    
     1063
    10641064
    10651065def dem2pts(basename_in, basename_out=None, verbose=False,
     
    11781178    for i in range(nrows):
    11791179        if verbose and i%((nrows+10)/10)==0:
    1180             print 'Processing row %d of %d' %(i, nrows)       
    1181        
     1180            print 'Processing row %d of %d' %(i, nrows)
     1181
    11821182        lower_index = global_index
    11831183        tpoints = zeros((ncols_in_bounding_box, 2), Float)
     
    12121212    """Return block of surface lines for each cross section
    12131213    Starts with SURFACE LINE,
    1214     Ends with END CROSS-SECTION 
     1214    Ends with END CROSS-SECTION
    12151215    """
    12161216
     
    12181218
    12191219    reading_surface = False
    1220     for i, line in enumerate(lines): 
    1221        
     1220    for i, line in enumerate(lines):
     1221
    12221222        if len(line.strip()) == 0:    #Ignore blanks
    1223             continue                       
     1223            continue
    12241224
    12251225        if lines[i].strip().startswith('SURFACE LINE'):
     
    12361236            easting = float(fields[0])
    12371237            northing = float(fields[1])
    1238             elevation = float(fields[2])                       
    1239             points.append([easting, northing, elevation])               
     1238            elevation = float(fields[2])
     1239            points.append([easting, northing, elevation])
    12401240
    12411241
     
    12461246                              verbose=False):
    12471247    """Read HEC-RAS Elevation datal from the following ASCII format (.sdf)
    1248    
     1248
    12491249    Example:
    1250    
     1250
    12511251
    12521252# RAS export file created on Mon 15Aug2005 11:42
     
    12621262  PROJECTION ZONE: 50
    12631263  DATUM: AGD66
    1264   VERTICAL DATUM: 
    1265   NUMBER OF REACHES:  19 
    1266   NUMBER OF CROSS-SECTIONS:  14206 
     1264  VERTICAL DATUM:
     1265  NUMBER OF REACHES:  19
     1266  NUMBER OF CROSS-SECTIONS:  14206
    12671267END HEADER:
    12681268
     
    12731273    STREAM ID:Southern-Wungong
    12741274    REACH ID:Southern-Wungong
    1275     STATION:19040.* 
     1275    STATION:19040.*
    12761276    CUT LINE:
    1277       405548.671603161 , 6438142.7594925 
    1278       405734.536092045 , 6438326.10404912 
    1279       405745.130459356 , 6438331.48627354 
    1280       405813.89633823 , 6438368.6272789 
     1277      405548.671603161 , 6438142.7594925
     1278      405734.536092045 , 6438326.10404912
     1279      405745.130459356 , 6438331.48627354
     1280      405813.89633823 , 6438368.6272789
    12811281    SURFACE LINE:
    12821282     405548.67,   6438142.76,   35.37
     
    12901290     405566.99,   6438160.83,   35.37
    12911291     ...
    1292    END CROSS-SECTION 
     1292   END CROSS-SECTION
    12931293
    12941294    Convert to NetCDF pts format which is
     
    13221322    assert lines[i].strip().upper() == 'BEGIN HEADER:'
    13231323    i += 1
    1324    
     1324
    13251325    assert lines[i].strip().upper().startswith('UNITS:')
    13261326    units = lines[i].strip().split()[1]
    1327     i += 1   
    1328 
    1329     assert lines[i].strip().upper().startswith('DTM TYPE:')   
    13301327    i += 1
    13311328
    1332     assert lines[i].strip().upper().startswith('DTM:')   
    1333     i += 1           
    1334 
    1335     assert lines[i].strip().upper().startswith('STREAM')   
     1329    assert lines[i].strip().upper().startswith('DTM TYPE:')
    13361330    i += 1
    13371331
    1338     assert lines[i].strip().upper().startswith('CROSS')   
     1332    assert lines[i].strip().upper().startswith('DTM:')
     1333    i += 1
     1334
     1335    assert lines[i].strip().upper().startswith('STREAM')
     1336    i += 1
     1337
     1338    assert lines[i].strip().upper().startswith('CROSS')
    13391339    i += 1
    13401340
     
    13541354    i += 1
    13551355
    1356     assert lines[i].strip().upper().startswith('NUMBER OF REACHES:') 
     1356    assert lines[i].strip().upper().startswith('NUMBER OF REACHES:')
    13571357    i += 1
    13581358
    13591359    assert lines[i].strip().upper().startswith('NUMBER OF CROSS-SECTIONS:')
    13601360    number_of_cross_sections = int(lines[i].strip().split(':')[1])
    1361     i += 1                       
     1361    i += 1
    13621362
    13631363
     
    13681368        for k, entry in enumerate(entries):
    13691369            points.append(entry[:2])
    1370             elevation.append(entry[2])           
     1370            elevation.append(entry[2])
    13711371
    13721372
    13731373    msg = 'Actual #number_of_cross_sections == %d, Reported as %d'\
    1374           %(j+1, number_of_cross_sections)         
     1374          %(j+1, number_of_cross_sections)
    13751375    assert j+1 == number_of_cross_sections, msg
    1376    
     1376
    13771377    #Get output file
    13781378    if basename_out == None:
     
    13901390    outfile.description = 'NetCDF pts format for compact and portable ' +\
    13911391                          'storage of spatial point data derived from HEC-RAS'
    1392    
     1392
    13931393    #Georeferencing
    13941394    outfile.zone = zone
     
    14291429            northing_min = None,
    14301430            northing_max = None,
    1431             expand_search = False, #To avoid intractable situations (This will be fixed when least_squares gets redesigned)           
     1431            expand_search = False, #To avoid intractable situations (This will be fixed when least_squares gets redesigned)
    14321432            verbose = False,
    14331433            origin = None,
    14341434            datum = 'WGS84',
    14351435            format = 'ers'):
    1436    
     1436
    14371437    """Read SWW file and convert to Digitial Elevation model format (.asc or .ers)
    14381438
     
    14651465    The parameter quantity must be the name of an existing quantity or
    14661466    an expression involving existing quantities. The default is
    1467     'elevation'. 
     1467    'elevation'.
    14681468
    14691469    if timestep (an index) is given, output quantity at that timestep
     
    14711471    if reduction is given use that to reduce quantity over all timesteps.
    14721472
    1473     datum 
    1474    
     1473    datum
     1474
    14751475    format can be either 'asc' or 'ers'
    14761476    """
     
    14791479    from Numeric import array, Float, concatenate, NewAxis, zeros, reshape, sometrue
    14801480    from Numeric import array2string
    1481    
     1481
    14821482    from utilities.polygon import inside_polygon, outside_polygon, separate_points_by_polygon
    14831483    from util import apply_expression_to_dictionary
    1484    
     1484
    14851485    msg = 'Format must be either asc or ers'
    14861486    assert format.lower() in ['asc', 'ers'], msg
    14871487
    1488        
     1488
    14891489    false_easting = 500000
    14901490    false_northing = 10000000
     
    14981498    if basename_out is None:
    14991499        basename_out = basename_in + '_%s' %quantity
    1500  
     1500
    15011501    swwfile = basename_in + '.sww'
    15021502    demfile = basename_out + '.' + format
     
    15901590    #Post condition: Now q has dimension: number_of_points
    15911591    assert len(q.shape) == 1
    1592     assert q.shape[0] == number_of_points   
    1593    
     1592    assert q.shape[0] == number_of_points
     1593
    15941594
    15951595    if verbose:
    15961596        print 'Processed values for %s are in [%f, %f]' %(quantity, min(q), max(q))
    1597    
     1597
    15981598
    15991599    #Create grid and update xll/yll corner and x,y
     
    16091609    else:
    16101610        xmax = easting_max - xllcorner
    1611        
    1612     if northing_min is None:       
     1611
     1612    if northing_min is None:
    16131613        ymin = min(y)
    16141614    else:
    16151615        ymin = northing_min - yllcorner
    1616        
    1617     if northing_max is None:       
     1616
     1617    if northing_max is None:
    16181618        ymax = max(y)
    16191619    else:
    16201620        ymax = northing_max - yllcorner
    1621        
    1622        
    1623    
     1621
     1622
     1623
    16241624    if verbose: print 'Creating grid'
    16251625    ncols = int((xmax-xmin)/cellsize)+1
     
    16451645        if format.lower() == 'asc':
    16461646            yg = i*cellsize
    1647         else:   
     1647        else:
    16481648            #this will flip the order of the y values for ers
    1649             yg = (nrows-i)*cellsize 
    1650            
     1649            yg = (nrows-i)*cellsize
     1650
    16511651        for j in xrange(ncols):
    16521652            xg = j*cellsize
     
    16581658    #Interpolate
    16591659    from least_squares import Interpolation
    1660    
     1660
    16611661
    16621662    #FIXME: This should be done with precrop = True (?), otherwise it'll
     
    16691669                           verbose = verbose)
    16701670
    1671    
     1671
    16721672
    16731673    #Interpolate using quantity values
     
    16851685
    16861686    for i in outside_indices:
    1687         grid_values[i] = NODATA_value       
     1687        grid_values[i] = NODATA_value
    16881688
    16891689
     
    16921692    if format.lower() == 'ers':
    16931693        # setup ERS header information
    1694         grid_values = reshape(grid_values,(nrows, ncols))   
     1694        grid_values = reshape(grid_values,(nrows, ncols))
    16951695        header = {}
    16961696        header['datum'] = '"' + datum + '"'
     
    17171717        ermapper_grids.write_ermapper_grid(demfile, grid_values, header)
    17181718
    1719         fid.close()   
     1719        fid.close()
    17201720    else:
    17211721        #Write to Ascii format
    1722        
     1722
    17231723        #Write prj file
    1724         prjfile = basename_out + '.prj' 
    1725        
     1724        prjfile = basename_out + '.prj'
     1725
    17261726        if verbose: print 'Writing %s' %prjfile
    17271727        prjid = open(prjfile, 'w')
     
    17381738
    17391739
    1740    
     1740
    17411741        if verbose: print 'Writing %s' %demfile
    1742  
     1742
    17431743        ascid = open(demfile, 'w')
    17441744
     
    17651765            ascid.write(s[1:-1] + '\n')
    17661766
    1767            
     1767
    17681768            #print
    17691769            #for j in range(ncols):
     
    17771777        fid.close()
    17781778
    1779 #Backwards compatibility       
     1779#Backwards compatibility
    17801780def sww2asc(basename_in, basename_out = None,
    17811781            quantity = None,
     
    17961796            datum = 'WGS84',
    17971797            format = 'asc')
    1798            
     1798
    17991799def sww2ers(basename_in, basename_out = None,
    18001800            quantity = None,
     
    18061806            datum = 'WGS84'):
    18071807    print 'sww2ers will soon be obsoleted - please use sww2dem'
    1808     sww2dem(basename_in, 
     1808    sww2dem(basename_in,
    18091809            basename_out = basename_out,
    18101810            quantity = quantity,
     
    18151815            origin = origin,
    18161816            datum = datum,
    1817             format = 'ers')         
    1818 ################################# END COMPATIBILITY ##############         
     1817            format = 'ers')
     1818################################# END COMPATIBILITY ##############
    18191819
    18201820
     
    19731973        fields = line.split()
    19741974        if verbose and i%((n+10)/10)==0:
    1975             print 'Processing row %d of %d' %(i, nrows)           
     1975            print 'Processing row %d of %d' %(i, nrows)
    19761976
    19771977        elevation[i, :] = array([float(x) for x in fields])
     
    20302030    from Numeric import Float, Int, Int32, searchsorted, zeros, array
    20312031    from Numeric import allclose, around
    2032        
     2032
    20332033    precision = Float
    20342034
     
    20632063    assert allclose(longitudes, file_u.variables['LON'])
    20642064    assert allclose(longitudes, file_v.variables['LON'])
    2065     assert allclose(longitudes, e_lon)   
     2065    assert allclose(longitudes, e_lon)
    20662066
    20672067
     
    21072107    zname = 'ELEVATION'
    21082108
    2109    
     2109
    21102110    amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax]
    21112111    uspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax] #Lon
     
    21132113    elevations = file_e.variables[zname][kmin:kmax, lmin:lmax]
    21142114
    2115 #    if latitudes2[0]==latitudes[0] and latitudes2[-1]==latitudes[-1]:
    2116 #        elevations = file_e.variables['ELEVATION'][kmin:kmax, lmin:lmax]
    2117 #    elif latitudes2[0]==latitudes[-1] and latitudes2[-1]==latitudes[0]:
    2118 #        from Numeric import asarray
    2119 #        elevations=elevations.tolist()
    2120 #        elevations.reverse()
    2121 #        elevations=asarray(elevations)
    2122 #    else:
    2123 #        from Numeric import asarray
    2124 #        elevations=elevations.tolist()
    2125 #        elevations.reverse()
    2126 #        elevations=asarray(elevations)
    2127 #        'print hmmm'
     2115    #    if latitudes2[0]==latitudes[0] and latitudes2[-1]==latitudes[-1]:
     2116    #        elevations = file_e.variables['ELEVATION'][kmin:kmax, lmin:lmax]
     2117    #    elif latitudes2[0]==latitudes[-1] and latitudes2[-1]==latitudes[0]:
     2118    #        from Numeric import asarray
     2119    #        elevations=elevations.tolist()
     2120    #        elevations.reverse()
     2121    #        elevations=asarray(elevations)
     2122    #    else:
     2123    #        from Numeric import asarray
     2124    #        elevations=elevations.tolist()
     2125    #        elevations.reverse()
     2126    #        elevations=asarray(elevations)
     2127    #        'print hmmm'
    21282128
    21292129
     
    23592359    outfile.variables['y'][:] = y - yllcorner
    23602360    outfile.variables['z'][:] = z             #FIXME HACK
    2361     outfile.variables['elevation'][:] = z 
     2361    outfile.variables['elevation'][:] = z
    23622362    outfile.variables['time'][:] = times   #Store time relative
    23632363    outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
     
    23842384                i += 1
    23852385
    2386     #outfile.close()               
     2386    #outfile.close()
    23872387
    23882388    #FIXME: Refactor using code from file_function.statistics
     
    24132413
    24142414
    2415     outfile.close()               
     2415    outfile.close()
    24162416
    24172417
     
    24912491        fields = line.split(',')
    24922492        realtime = calendar.timegm(time.strptime(fields[0], time_format))
    2493        
     2493
    24942494        T[i] = realtime - starttime
    24952495
     
    25162516    #FIXME: Use Georef
    25172517    fid.starttime = starttime
    2518    
    2519    
     2518
     2519
    25202520    # dimension definitions
    25212521    #fid.createDimension('number_of_volumes', self.number_of_volumes)
     
    25282528
    25292529    fid.variables['time'][:] = T
    2530    
     2530
    25312531    for i in range(Q.shape[1]):
    25322532        try:
    25332533            name = quantity_names[i]
    2534         except:   
     2534        except:
    25352535            name = 'Attribute%d'%i
    2536            
     2536
    25372537        fid.createVariable(name, Float, ('number_of_timesteps',))
    2538         fid.variables[name][:] = Q[:,i]       
     2538        fid.variables[name][:] = Q[:,i]
    25392539
    25402540    fid.close()
     
    26422642
    26432643    if verbose: print '    building domain'
    2644 #    From domain.Domain:
    2645 #    domain = Domain(coordinates, volumes,\
    2646 #                    conserved_quantities = conserved_quantities,\
    2647 #                    other_quantities = other_quantities,zone=zone,\
    2648 #                    xllcorner=xllcorner, yllcorner=yllcorner)
    2649 
    2650 #   From shallow_water.Domain:
     2644    #    From domain.Domain:
     2645    #    domain = Domain(coordinates, volumes,\
     2646    #                    conserved_quantities = conserved_quantities,\
     2647    #                    other_quantities = other_quantities,zone=zone,\
     2648    #                    xllcorner=xllcorner, yllcorner=yllcorner)
     2649
     2650    #   From shallow_water.Domain:
    26512651    coordinates=coordinates.tolist()
    26522652    volumes=volumes.tolist()
     
    26922692            X = resize(X,(len(X)/3,3))
    26932693        domain.set_quantity(quantity,X)
    2694 #
     2694    #
    26952695    for quantity in conserved_quantities:
    26962696        try:
     
    31183118
    31193119
    3120    
     3120
    31213121def sww2asc_obsolete(basename_in, basename_out = None,
    31223122            quantity = None,
     
    32683268    #Interpolate
    32693269    from least_squares import Interpolation
    3270    
     3270
    32713271
    32723272    #FIXME: This should be done with precrop = True, otherwise it'll
     
    32983298
    32993299
    3300    
     3300
    33013301    if verbose: print 'Writing %s' %ascfile
    33023302    NODATA_value = -9999
    3303  
     3303
    33043304    ascid = open(ascfile, 'w')
    33053305
     
    33903390    """
    33913391    Produce an sww boundary file, from esri ascii data from CSIRO.
    3392    
     3392
    33933393    Also convert latitude and longitude to UTM. All coordinates are
    33943394    assumed to be given in the GDA94 datum.
    3395    
     3395
    33963396    assume:
    33973397    All files are in esri ascii format
     
    34103410    """
    34113411    from Scientific.IO.NetCDF import NetCDFFile
    3412        
     3412
    34133413    from coordinate_transforms.redfearn import redfearn
    3414    
     3414
    34153415    precision = Float # So if we want to change the precision its done here
    3416    
     3416
    34173417    # go in to the bath dir and load the only file,
    34183418    bath_files = os.listdir(bath_dir)
    3419     #print "bath_files",bath_files 
     3419    #print "bath_files",bath_files
    34203420
    34213421    #fixme: make more general?
     
    34243424    bath_metadata,bath_grid =  _read_asc(bath_dir_file)
    34253425    #print "bath_metadata",bath_metadata
    3426     #print "bath_grid",bath_grid 
     3426    #print "bath_grid",bath_grid
    34273427
    34283428    #Use the date.time of the bath file as a basis for
    34293429    #the start time for other files
    34303430    base_start = bath_file[-12:]
    3431    
     3431
    34323432    #go into the elevation dir and load the 000 file
    34333433    elevation_dir_file = elevation_dir  + os.sep + elevation_prefix \
    34343434                         + base_start
    34353435    #elevation_metadata,elevation_grid =  _read_asc(elevation_dir_file)
    3436     #print "elevation_dir_file",elevation_dir_file 
    3437     #print "elevation_grid", elevation_grid 
    3438    
     3436    #print "elevation_dir_file",elevation_dir_file
     3437    #print "elevation_grid", elevation_grid
     3438
    34393439    elevation_files = os.listdir(elevation_dir)
    34403440    ucur_files = os.listdir(ucur_dir)
     
    34443444    # file with the same base name as the bath data
    34453445    #print "elevation_files[0]",elevation_files[0]
    3446     #print "'el' + base_start",'el' + base_start 
     3446    #print "'el' + base_start",'el' + base_start
    34473447    assert elevation_files[0] == 'el' + base_start
    3448    
     3448
    34493449    #assert bath_metadata == elevation_metadata
    34503450
    34513451
    34523452
    3453     number_of_latitudes = bath_grid.shape[0] 
     3453    number_of_latitudes = bath_grid.shape[0]
    34543454    number_of_longitudes = bath_grid.shape[1]
    34553455    #number_of_times = len(os.listdir(elevation_dir))
     
    34613461    latitudes = [bath_metadata['yllcorner']+y*bath_metadata['cellsize'] \
    34623462                 for y in range(number_of_latitudes)]
    3463    
     3463
    34643464     # reverse order of lat, so the fist lat represents the first grid row
    34653465    latitudes.reverse()
    3466    
    3467     #print "latitudes - before _get_min_max_indexes",latitudes 
     3466
     3467    #print "latitudes - before _get_min_max_indexes",latitudes
    34683468    kmin, kmax, lmin, lmax = _get_min_max_indexes(latitudes,longitudes,
    34693469                                                 minlat=minlat, maxlat=maxlat,
    34703470                                                 minlon=minlon, maxlon=maxlon)
    3471    
     3471
    34723472
    34733473    bath_grid = bath_grid[kmin:kmax,lmin:lmax]
     
    34803480    number_of_points = number_of_latitudes*number_of_longitudes
    34813481    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
    3482     #print "number_of_points",number_of_points 
     3482    #print "number_of_points",number_of_points
    34833483
    34843484    #Work out the times
     
    34923492    #print "times", times
    34933493    #print "number_of_latitudes", number_of_latitudes
    3494     #print "number_of_longitudes", number_of_longitudes 
    3495     #print "number_of_times", number_of_times 
     3494    #print "number_of_longitudes", number_of_longitudes
     3495    #print "number_of_times", number_of_times
    34963496    #print "latitudes", latitudes
    34973497    #print "longitudes", longitudes
     
    35103510        print '    t in [%f, %f], len(t) == %d'\
    35113511              %(min(times), max(times), len(times))
    3512    
     3512
    35133513    ######### WRITE THE SWW FILE #############
    35143514    # NetCDF file definition
     
    35233523    outfile.smoothing = 'Yes'
    35243524    outfile.order = 1
    3525    
     3525
    35263526    #Start time in seconds since the epoch (midnight 1/1/1970)
    35273527    outfile.starttime = starttime = times[0]
    3528    
    3529    
     3528
     3529
    35303530    # dimension definitions
    35313531    outfile.createDimension('number_of_volumes', number_of_volumes)
     
    35703570    #Get zone of 1st point.
    35713571    refzone, _, _ = redfearn(latitudes[0],longitudes[0])
    3572    
     3572
    35733573    vertices = {}
    35743574    i = 0
    3575     for k, lat in enumerate(latitudes):       
    3576         for l, lon in enumerate(longitudes): 
     3575    for k, lat in enumerate(latitudes):
     3576        for l, lon in enumerate(longitudes):
    35773577
    35783578            vertices[l,k] = i
     
    36393639    xmomentum = outfile.variables['xmomentum']
    36403640    ymomentum = outfile.variables['ymomentum']
    3641    
     3641
    36423642    outfile.variables['time'][:] = times   #Store time relative
    3643    
     3643
    36443644    if verbose: print 'Converting quantities'
    3645     n = number_of_times   
     3645    n = number_of_times
    36463646    for j in range(number_of_times):
    36473647        # load in files
    36483648        elevation_meta, elevation_grid = \
    36493649           _read_asc(elevation_dir + os.sep + elevation_files[j])
    3650                          
     3650
    36513651        _, u_momentum_grid =  _read_asc(ucur_dir + os.sep + ucur_files[j])
    36523652        _, v_momentum_grid =  _read_asc(vcur_dir + os.sep + vcur_files[j])
    36533653
    3654         #print "elevation_grid",elevation_grid 
     3654        #print "elevation_grid",elevation_grid
    36553655        #cut matrix to desired size
    36563656        elevation_grid = elevation_grid[kmin:kmax,lmin:lmax]
    36573657        u_momentum_grid = u_momentum_grid[kmin:kmax,lmin:lmax]
    36583658        v_momentum_grid = v_momentum_grid[kmin:kmax,lmin:lmax]
    3659         #print "elevation_grid",elevation_grid 
     3659        #print "elevation_grid",elevation_grid
    36603660        # handle missing values
    36613661        missing = (elevation_grid == elevation_meta['NODATA_value'])
     
    36813681                i += 1
    36823682    outfile.close()
    3683    
     3683
    36843684def _get_min_max_indexes(latitudes,longitudes,
    36853685                        minlat=None, maxlat=None,
     
    37093709        if lat_min_index <0:
    37103710            lat_min_index = 0
    3711            
     3711
    37123712
    37133713    if maxlat == None:
     
    37173717        if lat_max_index > largest_lat_index:
    37183718            lat_max_index = largest_lat_index
    3719        
     3719
    37203720    if minlon == None:
    37213721        lon_min_index = 0
     
    37243724        if lon_min_index <0:
    37253725            lon_min_index = 0
    3726        
     3726
    37273727    if maxlon == None:
    37283728        lon_max_index = len(longitudes)
     
    37303730        lon_max_index = searchsorted(longitudes, maxlon)
    37313731
    3732     #Take into account that the latitude list was reversed 
     3732    #Take into account that the latitude list was reversed
    37333733    latitudes.reverse() # Python passes by reference, need to swap back
    37343734    lat_min_index, lat_max_index = largest_lat_index - lat_max_index , \
     
    37363736    lat_max_index = lat_max_index + 1 # taking into account how slicing works
    37373737    lon_max_index = lon_max_index + 1 # taking into account how slicing works
    3738    
     3738
    37393739    return lat_min_index, lat_max_index, lon_min_index, lon_max_index
    3740        
    3741    
     3740
     3741
    37423742def _read_asc(filename, verbose=False):
    37433743    """Read esri file from the following ASCII format (.asc)
     
    37693769    cellsize = float(lines.pop(0).split()[1].strip())
    37703770    NODATA_value = float(lines.pop(0).split()[1].strip())
    3771    
    3772     assert len(lines) == nrows 
    3773    
     3771
     3772    assert len(lines) == nrows
     3773
    37743774    #Store data
    37753775    grid = []
    3776    
     3776
    37773777    n = len(lines)
    37783778    for i, line in enumerate(lines):
     
    37863786            'cellsize':cellsize,
    37873787            'NODATA_value':NODATA_value}, grid
    3788    
     3788
  • inundation/pyvolution/netherlands.py

    r1828 r2105  
    9090print 'Creating domain'
    9191#Create basic mesh
    92 points, vertices, boundary = rectangular(N, N)
     92points, elements, boundary = rectangular(N, N)
    9393
    9494#Create shallow water domain
    95 domain = Domain(points, vertices, boundary, use_inscribed_circle=True)
     95domain = Domain(points, elements, boundary, use_inscribed_circle=True)
    9696
    9797domain.check_integrity()
     
    107107
    108108#Set bed-slope and friction
    109 inflow_stage = 0.1
     109inflow_stage = 0.5
    110110manning = 0.03
    111111Z = Weir(inflow_stage)
     
    155155#
    156156print 'Initial condition'
    157 #domain.set_quantity('stage', Constant_height(Z, 0.0))
    158 domain.set_quantity('stage', Constant_stage(inflow_stage/2.0))
     157domain.set_quantity('stage', Constant_height(Z, 0.0))
     158#domain.set_quantity('stage', Constant_stage(inflow_stage/2.0))
    159159
    160160#Evolve
     
    168168#V.update_quantity('elevation')
    169169
    170 for t in domain.evolve(yieldstep = 0.02, finaltime = 15.0):
     170for t in domain.evolve(yieldstep = 0.005, finaltime = 15.0):
    171171    domain.write_time()
    172     domain.write_boundary()   
    173    
     172    #domain.write_boundary()
     173
    174174    print domain.quantities['stage'].get_values(location='centroids',
    175175                                                indices=[0])
  • inundation/pyvolution/view_tsh.bat

    r1293 r2105  
    1 python C:\Home\Projects\anuga\pyvolution\view_tsh.py %1  %2
     1python C:\Home\Projects\anuga\inundation\pyvolution\view_tsh.py %1  %2
  • inundation/zeus/Merimbula.zpi

    r1358 r2105  
    7373    <ReleaseProjectSaveAll>Off</ReleaseProjectSaveAll>
    7474    <ReleaseProjectReload>Off</ReleaseProjectReload>
     75    <file>..\merimbula\prepare.py</file>
    7576    <file>..\merimbula\run_least_squares_merimbula.py</file>
    7677    <file>..\merimbula\run_merimbula_lake.py</file>
  • inundation/zeus/Validation.zpi

    r1735 r2105  
    7373    <ReleaseProjectSaveAll>Off</ReleaseProjectSaveAll>
    7474    <ReleaseProjectReload>Off</ReleaseProjectReload>
    75     <file>..\validation\lwru2\create_mesh.py</file>
    76     <file>..\validation\lwru2\extract_timeseries.py</file>
    77     <file>..\validation\lwru2\filenames.py</file>
    78     <file>..\validation\lwru2\lwru2.py</file>
     75    <file>..\validation\completed\lwru2\create_mesh.py</file>
     76    <file>..\validation\completed\lwru2\extract_timeseries.py</file>
     77    <file>..\validation\completed\lwru2\lwru2.py</file>
     78    <file>..\validation\completed\lwru2\project.py</file>
    7979    <folder name="Header Files" />
    8080    <folder name="Resource Files" />
  • inundation/zeus/anuga-workspace.zwi

    r1995 r2105  
    22<workspace name="anuga-workspace">
    33    <mode>Debug</mode>
    4     <active>analytic_solutions</active>
     4    <active>parallel</active>
    55    <project name="analytic_solutions">analytic_solutions.zpi</project>
    66    <project name="euler">euler.zpi</project>
  • inundation/zeus/parallel.zpi

    r1697 r2105  
    8686    <file>..\parallel\run_parallel_mesh.py</file>
    8787    <file>..\parallel\run_parallel_sw_merimbula.py</file>
     88    <file>..\parallel\run_parallel_sw_merimbula_metis.py</file>
    8889    <file>..\parallel\run_parallel_sw_rectangle.py</file>
    8990    <file>..\parallel\run_sw_lost_mass.py</file>
  • inundation/zeus/pyvolution.zpi

    r1949 r2105  
    7979    <file>..\pyvolution\check_sww_tsh.py</file>
    8080    <file>..\pyvolution\combine_pts.py</file>
    81     <file>..\pyvolution\compile.py</file>
    8281    <file>..\pyvolution\config.py</file>
    8382    <file>..\pyvolution\cornell_room.py</file>
    8483    <file>..\pyvolution\data_manager.py</file>
    8584    <file>..\pyvolution\domain.py</file>
    86     <file>..\pyvolution\euler.py</file>
    87     <file>..\pyvolution\euler_config.py</file>
    88     <file>..\pyvolution\euler_ext.c</file>
    8985    <file>..\pyvolution\flatbed.py</file>
    9086    <file>..\pyvolution\flatbed_compare.py</file>
     
    123119    <file>..\pyvolution\test_data_manager.py</file>
    124120    <file>..\pyvolution\test_domain.py</file>
    125     <file>..\pyvolution\test_euler.py</file>
    126121    <file>..\pyvolution\test_general_mesh.py</file>
    127122    <file>..\pyvolution\test_generic_boundary_conditions.py</file>
Note: See TracChangeset for help on using the changeset viewer.