Changeset 4429


Ignore:
Timestamp:
May 15, 2007, 4:13:58 PM (17 years ago)
Author:
nick
Message:

update broome

Location:
anuga_work/production/broome_2006
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/production/broome_2006/build_broome.py

    r4366 r4429  
    8383dem2pts(offshore_dir_name1, use_cache=True, verbose=True)
    8484dem2pts(offshore_dir_name2, use_cache=True, verbose=True)
    85 '''
     85
    8686print'create Geospatial data1 objects from topographies',onshore_dir_name + '.pts'
    8787G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
    … …  
    110110print'export',project.combined_dir_name+'_unclipped' + '.txt'
    111111G.export_points_file(project.combined_dir_name+'_unclipped' + '.txt')
    112 '''
     112
    113113print'split'
    114114G_small, G_other = G_clipped.split(.10,verbose=True)
    … …  
    120120G_small.export_points_file(project.combined_small_dir_name + '.txt')
    121121#G_clipped.export_points_file(project.combined_dir_name + '.xya')
    122 '''
    123122
    124 '''
     123
     124
    125125print'project.combined_dir_name + 1.xya',project.combined_dir_name + '1.xya'
    126126G_all=Geospatial_data(file_name = project.combined_dir_name + '1.xya')
    … …  
    132132G_all_2.export_points_file(project.combined_dir_name+'_other1' + '.xya')
    133133
    134 
     134'''
    135135#-------------------------------------------------------------------------
    136136# Convert URS to SWW file for boundary conditions
    137137#-------------------------------------------------------------------------
    138 '''
     138print 'starting to create boundary conditions'
     139boundaries_in_dir_name = project.boundaries_in_dir_name
     140
     141from anuga.shallow_water.data_manager import urs2sww, urs_ungridded2sww
     142
     143print 'boundaries_in_dir_name',boundaries_in_dir_name
     144print 'project.boundaries_dir_name',project.boundaries_dir_name
     145
     146urs_ungridded2sww(boundaries_in_dir_name, project.boundaries_dir_name,
     147                  verbose=True, mint=4000, maxt=35000, zscale=1)
    139148
    140149
    … …  
    143152
    144153
    145 

  • anuga_work/production/broome_2006/export_results.py

    r4347 r4429  
    55from os import sep
    66
    7 time_dir = '20061113_040953'
     7time_dir = '20070423_040235_run' #
     8cellsize = 25
     9timestep = None
     10directory = project.output_dir
     11#name = directory + time_dir + sep + project.scenario_name
     12name = directory+time_dir+sep+project.scenario_name
    813
    9 directory = project.outputdir
    10 name = directory + time_dir +sep + 'source'
     14is_parallel = True
     15if is_parallel == True: nodes = 4
     16print 'output dir:', name
    1117
    12 print 'output dir:', name
    13 which_var = 4
     18which_var = 2
    1419if which_var == 0:  # Stage
    1520    outname = name + '_stage'
    … …  
    1722
    1823if which_var == 1:  # Absolute Momentum
    19     outname = name + '_momentum'
    20     quantityname = '(xmomentum**2 + ymomentum**2)**0.5'  #Absolute momentum
     24    outname = name + '_momentum_i1'
     25    quantityname = '(xmomentum**2 + ymomentum**2)**0.5' 
    2126
    2227if which_var == 2:  # Depth
    2328    outname = name + '_depth'
    24     quantityname = 'stage-elevation'  #Depth
     29    quantityname = 'stage-elevation' 
    2530
    2631if which_var == 3:  # Speed
    27     outname = name + '_speed'
     32    outname = name + '_speed_i0'
    2833    quantityname = '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6/(stage-elevation))'  #Speed
    2934
    … …  
    3237    quantityname = 'elevation'  #Elevation
    3338
    34 print 'start sww2dem'
    3539
    36 sww2dem(name, basename_out = outname,
    37             quantity = quantityname,
    38             cellsize = 25,      # would prefer this at 25
    39             # define region for viz purposes
    40             easting_min = project.eastingmin,
    41             easting_max = project.eastingmax,
    42             northing_min = project.northingmin,
    43             northing_max = project.northingmax,       
    44             reduction = max, #this is because we want max quantityname
    45             verbose = True,
    46             format = 'asc')
    4740
     41if is_parallel == True:
     42#    print 'is_parallel',is_parallel
     43    for i in range(0,nodes):
     44        namei = name + '_P%d_%d' %(i,nodes)
     45        outnamei = outname + '_P%d_%d' %(i,nodes)
     46        print 'start sww2dem for sww file %d' %(i)
     47        sww2dem(namei, basename_out = outnamei,
     48                    quantity = quantityname,
     49                    timestep = timestep,
     50                    cellsize = cellsize,     
     51                    easting_min = project.e_min_area,
     52                    easting_max = project.e_max_area,
     53                    northing_min = project.n_min_area,
     54                    northing_max = project.n_max_area,       
     55                    reduction = max,
     56                    verbose = True,
     57                    format = 'asc')
     58else:
     59    print 'start sww2dem'
     60    sww2dem(name, basename_out = outname,
     61                quantity = quantityname,
     62                timestep = timestep,
     63                cellsize = cellsize,     
     64                easting_min = project.e_min_area,
     65                easting_max = project.e_max_area,
     66                northing_min = project.n_min_area,
     67                northing_max = project.n_max_area,       
     68                reduction = max,
     69                verbose = True,
     70                format = 'asc')

  • anuga_work/production/broome_2006/project.py

    r4372 r4429  
    2222
    2323#time stuff
    24 time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
     24time = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
     25#time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
    2526gtime = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
    2627build_time = time+'_build'
    … …  
    2930
    3031#tide = -5.3
    31 #tide = 0
    32 tide = 4.9
     32tide = 0
     33#tide = 4.9
    3334
    3435#Making assumptions about the location of scenario data
    … …  
    9394
    9495
    95 boundaries_source = '????'
     96boundaries_source = 'broome_3854_17042007'
    9697#boundaries locations
    97 boundaries_in_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+boundaries_source+sep
    98 boundaries_in_dir_name = boundaries_in_dir + scenario_name
     98boundaries_in_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep
     99#boundaries_in_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+boundaries_source+sep
     100boundaries_in_dir_name = boundaries_in_dir + boundaries_source
    99101boundaries_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep
    100 boundaries_dir_name = boundaries_dir + scenario_name
     102boundaries_dir_name = boundaries_dir + boundaries_source
    101103#boundaries_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+build_time+sep
    102104#boundaries_time_dir_name = boundaries_time_dir + boundaries_name  #Used by post processing
    … …  
    126128poly_all = read_polygon(polygons_dir+'extent_small.csv')
    127129#poly_all = read_polygon(polygons_dir+'extent.csv')
    128 res_poly_all = 500000
    129 #res_poly_all = 150000
     130res_factor = 2
     131#res_poly_all = 500000
     132res_poly_all = 150000*res_factor
    130133
    131134###############################
    … …  
    134137
    135138poly_0 = read_polygon(polygons_dir+'neg20_coast_contour_pts.csv')
    136 res_0 = 100000
    137 #res_0 = 20000
     139#res_0 = 100000
     140res_0 = 20000*res_factor
    138141
    139142poly_1 = read_polygon(polygons_dir+'broome_north_coast_inside_extent.csv')
    140 res_1 = 50000
    141 #res_1 = 5000
     143#res_1 = 50000
     144res_1 = 5000*res_factor
    142145
    143146poly_2 = read_polygon(polygons_dir+'broome_south_coast_inside_extent.csv')
    144 res_2 = 50000
    145 #res_2 = 5000
     147#res_2 = 50000
     148res_2 = 5000*res_factor
    146149
    147150poly_3 = read_polygon(polygons_dir+'Broome_town_pts.csv')
    148 res_3 = 20000
    149 #res_3 = 2000
     151#res_3 = 20000
     152res_3 = 2000*res_factor
    150153
    151154poly_4 = read_polygon(polygons_dir+'Broome_inner_town_pts.csv')
    152 res_4 = 5000
    153 #res_4 = 500
     155#res_4 = 5000
     156res_4 = 500*res_factor
    154157#assert zone == refzone
    155158
    … …  
    161164print 'min number triangles', trigs_min
    162165
     166poly_mainland = read_polygon(polygons_dir+'Initial_Condition.csv')
     167
    163168###################################################################
    164169# Clipping regions for export to asc and regions for clipping data
    … …  
    166171
    167172# exporting asc grid
    168 eastingmin = 406215.87
    169 eastingmax = 440208.78
    170 northingmin = 7983427.73
    171 northingmax = 8032834.52
     173e_min_area = 412000.0
     174e_max_area = 423000.0
     175n_min_area = 8007000.0
     176n_max_area = 8022000.0
    172177
    173178

  • anuga_work/production/broome_2006/run_broome.py

    r4372 r4429  
    8181                         interior_regions=project.interior_regions,
    8282                         filename=meshes_dir_name,
    83                          use_cache=True,
     83                         use_cache=False,
    8484                         verbose=True)
    8585
    … …  
    9292print 'Setup computational domain'
    9393#domain = Domain(meshes_time_dir_name, use_cache=True, verbose=True)
    94 domain = Domain(meshes_dir_name, use_cache=True, verbose=True)
     94domain = Domain(meshes_dir_name, use_cache=False, verbose=True)
    9595print domain.statistics()
    96 '''
    97 print 'starting to create boundary conditions'
    98 boundaries_in_dir_name = project.boundaries_in_dir_name
    99 
    100 from anuga.shallow_water.data_manager import urs2sww, ferret2sww
    101 
    102 # put above distribute
    103 print 'boundary file is: ',boundaries_dir_name
    104 from caching import cache
    105 if myid == 0:
    106     cache(ferret2sww,
    107           (boundaries_in_dir_name,
    108     #       boundaries_time_dir_name),
    109            boundaries_dir_name),
    110           {'verbose': True,
    111            'minlat': project.south_boundary,
    112            'maxlat': project.north_boundary,
    113            'minlon': project.west_boundary,
    114            'maxlon': project.east_boundary,
    115 #           'minlat': project.south,
    116 #           'maxlat': project.north,
    117 #           'minlon': project.west,
    118 #           'maxlon': project.east,
    119            'mint': 0, 'maxt': 35100,
    120            'origin': domain.geo_reference.get_origin(),
    121            'mean_stage': project.tide,
    122 #           'zscale': 1,                 #Enhance tsunami
    123            'fail_on_NaN': False},
    124            verbose = True,
    125            )
    126 barrier()
    127 '''
    12896
    12997#-------------------------------------------------------------------------
    … …  
    134102    print 'Setup initial conditions'
    135103
    136     domain.set_quantity('stage', tide)
    137     domain.set_quantity('friction', 0.01)
     104    from polygon import Polygon_function
     105    #following sets the stage/water to be offcoast only
     106    IC = Polygon_function( [(project.poly_mainland, -1.0)], default = tide,
     107                             geo_reference = domain.geo_reference)
     108    domain.set_quantity('stage', IC)
    138109    #combined_time_dir_name = project.topographies_dir+build_time+project.combined_name
    139110    print 'Start Set quantity'
    140111    print 'project.combined_dir_name_unclipped1',project.combined_dir_name_unclipped1+'.txt'
    141112    domain.set_quantity('elevation',
     113#                    filename = project.combined_dir_name+'.txt',
    142114                    filename = project.combined_dir_name_unclipped1+'.txt',
    143115#                    filename = project.combined_smaller_name_dir+'.xya',
    … …  
    175147print 'domain id', id(domain)
    176148#Bf = File_boundary(boundaries_dir_name + '.sww',
    177 #                  domain, time_thinning=5, use_cache=True, verbose=True)
     149#                  domain, time_thinning=12, use_cache=True, verbose=True)
     150Bf = Field_boundary(boundaries_dir_name + '.sww',
     151                  domain, time_thinning=1, mean_stage=tide,
     152                  use_cache=True, verbose=True)
    178153
    179154print 'finished reading boundary file'
    … …  
    185160domain.set_boundary({'back': Br,
    186161                     'side': Bd,
    187                      'ocean': Bd})
     162                     'ocean': Bf})
     163#                     'ocean': Bd})
    188164print'finish set boundary'
    189165
    … …  
    194170t0 = time.time()
    195171
    196 for t in domain.evolve(yieldstep = 60, finaltime = 10000):
     172for t in domain.evolve(yieldstep = 60, finaltime = 30000):
    197173    domain.write_time()
    198174    domain.write_boundary_statistics(tags = 'ocean')
    199     if allclose(t, 120):
    200         domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bo})
     175#    if allclose(t, 120):
     176#        domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bo})
    201177
    202     if allclose(t, 1020):
    203         domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bd})
     178#    if allclose(t, 1020):
     179#        domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bd})
    204180
    205181   
Note: See TracChangeset for help on using the changeset viewer.