Changeset 3940


Ignore:
Timestamp:
Nov 8, 2006, 4:04:55 PM (17 years ago)
Author:
nick
Message:

update to dampier

Location:
anuga_work/production/dampier_2006
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/production/dampier_2006/build_dampier.py

    r3905 r3940  
    3232from anuga.pmesh.mesh_interface import create_mesh_from_regions
    3333from anuga.geospatial_data.geospatial_data import *
    34 from anuga.abstract_2d_finite_volumes.util import Screen_Catcher
     34from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
    3535
    3636# Application specific imports
    … …  
    4242#------------------------------------------------------------------------------
    4343
    44 if access(project.output_build_time_dir,F_OK) == 0:
    45     mkdir (project.output_build_time_dir)
    46 copy (dirname(project.__file__) +sep+ project.__name__+'.py',
    47       project.output_build_time_dir + project.__name__+'.py')     #copies project.py
    48 copy (__file__, project.output_build_time_dir + basename(__file__))
    49 print 'files '+ project.__name__+'.py and '+ basename(__file__)+' copied to '+ project.output_build_time_dir       #copies this file
    50 #import sys; sys.exit()
    51 
    52 
    53 #normal screen output is stored in
    54 screen_output_name = project.output_build_time_dir + "screen_output.txt"
    55 screen_error_name = project.output_build_time_dir + "screen_error.txt"
    56 
    57 #used to catch screen output to file
    58 sys.stdout = Screen_Catcher(screen_output_name)
    59 sys.stderr = Screen_Catcher(screen_error_name)
     44copy_code_files(project.output_build_time_dir,__file__,
     45               dirname(project.__file__)+sep+ project.__name__+'.py' )
     46
     47start_screen_catcher(project.output_build_time_dir)
    6048
    6149print 'USER:    ', project.user
    … …  
    125113G_off13 = Geospatial_data(file_name = project.offshore_dir_name13 + '.xya')
    126114G_off14 = Geospatial_data(file_name = project.offshore_dir_name14 + '.xya')
    127 
     115'''
     116
     117
     118print'clip nw', project.clip_poly_nw
     119print'clip e', project.clip_poly_e
     120
     121print'reading combined_dir_name'
     122G_offshore_data = Geospatial_data(file_name = project.topographies_dir+'dampier_combined_elevation_final.pts')
     123
     124print'reading offshore_dir_name_old'
     125G_offshore_old = Geospatial_data(file_name = project.offshore_dir_name_old + '.pts')
     126
     127'''
     128print 'G_offshore_data_old_nw',
     129G_nw_name = project.topographies_dir + 'nw_old_data'
     130G_offshore_nw = G_offshore_old.clip(project.clip_poly_nw)
     131G_offshore_nw.export_points_file(G_nw_name + '.pts')
     132G_offshore_nw.export_points_file(G_nw_name + '.xya')
     133G_nw = array(G_offshore_nw.get_data_points())
     134print'shape of arr nw data', G_nw.shape
     135print' max and min of array 0',max(G_nw[:,0]),min(G_nw[:,0])
     136print' max and min of array 1',max(G_nw[:,1]),min(G_nw[:,1])
     137
     138print 'G_offshore_data_old_e'#, G_offshore_old_nw.get_data_points()
     139G_e_name = project.topographies_dir+'e_old_data'
     140G_offshore_e = G_offshore_old.clip(project.clip_poly_e)
     141G_offshore_e.export_points_file(G_e_name + '.pts')
     142G_offshore_e.export_points_file(G_e_name + '.xya')
     143G_e = array(G_offshore_e.get_data_points())
     144print'shape of arr e data', G_e.shape
     145print' max and min of array 0',max(G_e[:,0]),min(G_e[:,0])
     146print' max and min of array 1',max(G_e[:,1]),min(G_e[:,1])
     147'''
     148
     149print 'G_offshore_data_mid_e'#, G_offshore_old_nw.get_data_points()
     150G_mid_e_name = project.topographies_dir+'mid_e_old_data'
     151G_offshore_mid_e = G_offshore_old.clip(project.clip_poly_mid_e)
     152G_offshore_mid_e.export_points_file(G_mid_e_name + '.pts')
     153G_offshore_mid_e.export_points_file(G_mid_e_name + '.xya')
     154G_mid_e = array(G_offshore_mid_e.get_data_points())
     155print'shape of arr e data', G_mid_e.shape
     156print' max and min of array 0',max(G_mid_e[:,0]),min(G_mid_e[:,0])
     157print' max and min of array 1',max(G_mid_e[:,1]),min(G_mid_e[:,1])
     158
     159print 'G_offshore_data_mid_w'#, G_offshore_old_nw.get_data_points()
     160G_mid_w_name = project.topographies_dir+'mid_w_old_data'
     161G_offshore_mid_w = G_offshore_old.clip(project.clip_poly_mid_w)
     162G_offshore_mid_w.export_points_file(G_mid_w_name + '.pts')
     163G_offshore_mid_w.export_points_file(G_mid_w_name + '.xya')
     164G_mid_w = array(G_offshore_mid_w.get_data_points())
     165print'shape of arr e data', G_mid_w.shape
     166print' max and min of array 0',max(G_mid_w[:,0]),min(G_mid_w[:,0])
     167print' max and min of array 1',max(G_mid_w[:,1]),min(G_mid_w[:,1])
     168
     169
     170print 'G_offshore_data_old_e'#, G_offshore_old_e.get_data_points()
    128171print'add all geospatial objects'
    129 G = G1 + G2 + G3 + G_off + G_off1 + G_off2 + G_off3 + G_off4 + G_off5 \
    130     + G_off6 + G_off7 + G_off8 + G_off9 + G_off10 + G_off11 + G_off12 \
    131     + G_off13 + G_off14
    132 
    133 print'clip combined geospatial object by bounding polygon'
    134 G.clip(project.bounding_polygon)
     172#G = G1 + G2 + G3 + G_off + G_off1 + G_off2 + G_off3 + G_off4 + G_off5 \
     173#    + G_off6 + G_off7 + G_off8 + G_off9 + G_off10 + G_off11 + G_off12 \
     174#    + G_off13 + G_off14
     175
     176G = G_offshore_data + G_offshore_mid_w + G_offshore_mid_e
     177print'shape of arr G data', G.get_data_points().shape
     178
     179
     180#print'clip combined geospatial object by bounding polygon'
     181G_clipped = G.clip(project.bounding_polygon)
    135182#FIXME: add a clip function to pts
     183print'shape of clipped data', G_clipped.get_data_points().shape
    136184
    137185print'export combined DEM file'
    138186if access(project.topographies_time_dir,F_OK) == 0:
    139187    mkdir (project.topographies_time_dir)
    140 G.export_points_file(project.combined_time_dir_name + '.pts')
     188G_clipped.export_points_file(project.combined_time_dir_final_name + '.pts')
     189G_clipped.export_points_file(project.combined_time_dir_final_name + '.xya')
    141190'''
    142191#-------------------------------------------------------------------------
    … …  
    156205if access(project.boundaries_time_dir,F_OK) == 0:
    157206    mkdir (project.boundaries_time_dir)
    158 '''
    159 urs2sww(boundaries_in_dir_name,basename_out= project.boundaries_time_dir_name,
    160         minlat=project.south_boundary, maxlat=project.north_boundary,
    161         minlon= project.west_boundary, maxlon=project.east_boundary,
    162         mint=0, maxt= 35000,
    163         verbose='true')
    164 
    165 
    166 urs2sww(boundaries_in_dir_name, basename_out= project.boundaries_time_dir_name,
    167 #        minlat=project.south, maxlat=project.north,
    168 #        minlon= project.west, maxlon=project.east,
    169         minlat=project.south_boundary, maxlat=project.north_boundary,
    170         minlon= project.east_boundary, maxlon=project.west_boundary,
    171         mint=0, maxt= 35000,
    172         verbose='true')
    173 '''
     207
    174208from caching import cache
    175209cache(urs2sww,
    … …  
    193227       )
    194228#       dependencies = source_dir + project.boundary_basename + '.sww')
    195    
    196 
    197 
    198 
    199 
    200 
    201 
    202 
    203 
     229'''   
     230
     231
     232
     233
     234
     235
     236
     237

  • anuga_work/production/dampier_2006/project.py

    r3937 r3940  
    3333
    3434#mesh_name = 'elevation50m'
    35 boundaries_name = 'dampier1'
     35boundaries_name = 'dampier'
    3636boundaries_source = 'mag_9_corrected'
    3737#boundaries_source = 'test'
    … …  
    6060offshore_name14 = 'XYDM83'
    6161
     62offshore_old = 'elevation50m'
     63
    6264combined_name ='dampier_combined_elevation'
    63 
    64 
     65combined_final_name ='dampier_combined_elevation_final'
     66
     67gauge_name = 'dampier_gauges.csv'
    6568
    6669#Derive subdirectories and filenames
    … …  
    8386#print 'ctime: ', cctime
    8487
     88output_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep
    8589output_build_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep+build_time+sep
    8690output_run_time_dir = home+sep+state+sep+scenario_datas_name+sep+'anuga'+sep+'outputs'+sep+run_time+sep
    … …  
    9498
    9599
    96 #gauge_filename = gaugedir + 'gauge_location_broome.csv'
     100gauges_dir_name = gauges_dir + gauge_name
    97101
    98102onshore_dir_name = topographies_dir + onshore_name
    … …  
    115119offshore_dir_name14 = topographies_dir + offshore_name14
    116120
     121offshore_dir_name_old = topographies_dir + offshore_old
     122
     123
     124
    117125#output dir
    118126combined_dir_name = topographies_dir + combined_name
    119127combined_time_dir_name = topographies_time_dir + combined_name
     128combined_time_dir_final_name = topographies_time_dir + combined_final_name
     129
    120130
    121131meshes_dir_name = meshes_dir + scenario_name
    … …  
    125135boundaries_in_dir_name = boundaries_in_dir + boundaries_name
    126136boundaries_time_dir_name = boundaries_time_dir + boundaries_name  #Used by post processing
    127 
     137boundaries_dir_name = boundaries_dir + boundaries_name
    128138
    129139
    … …  
    164174
    165175bounding_polygon, zone =\
    166                   convert_from_latlon_to_utm([p0, p1, p2, p3, p4,p5, p6, p7, p8])
     176                  convert_from_latlon_to_utm([p0, p1, p2, p3, p4, p5, p6, p7, p8])
    167177#bounding_polygon, zone =\
    168178#                  convert_from_latlon_to_utm([p1, p2, p3, p4, p5, p6, p7])
    … …  
    192202
    193203poly_pipeline = read_polygon(polygons_dir+'pipeline.csv')
     204
     205clip_poly_e = read_polygon(polygons_dir+'gap_e.csv')
     206
     207clip_poly_nw = read_polygon(polygons_dir+'gap_nw.csv')
     208
     209clip_poly_mid_w = read_polygon(polygons_dir+'gap_mid_w.cvs')
     210
     211clip_poly_mid_e = read_polygon(polygons_dir+'gap_mid_e.cvs')
    194212
    195213#Interior regions

  • anuga_work/production/dampier_2006/run_dampier.py

    r3905 r3940  
    3333
    3434from anuga.geospatial_data.geospatial_data import *
    35 from anuga.abstract_2d_finite_volumes.util import Screen_Catcher
     35from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
    3636from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
    3737
    … …  
    4444#------------------------------------------------------------------------------
    4545
     46
    4647# filenames
    4748
    4849#build_time = '20061029_231935_build_tide_24'
    49 build_time = '20061030_165746_build_tide_24'
     50#build_time = '20061030_165746_build_tide_24'
     51#build_time = '20061102_215532_build_plus_old_data'
     52#build_time = '20061103_055258_build'
     53build_time = '20061107_063840_build'
    5054#build_time = '20061025_153643_build_basic'
     55bound_time = '20061102_221245_build'
    5156
    5257boundaries_name = project.boundaries_name
    53 meshes_time_dir_name = project.meshes_time_dir_name+'.msh'
     58#meshes_time_dir_name = project.meshes_time_dir_name+'.msh'
     59meshes_dir_name = project.meshes_dir_name+'.msh'
    5460#source_dir = project.boundarydir
    5561#boundaries_time_dir_name = project.boundaries_dir + build_time + sep + boundaries_name
    5662boundaries_time_dir_name = project.boundaries_time_dir_name
     63boundaries_dir_name = project.boundaries_dir_name
    5764tide = project.tide
    5865
    59 # # creates copy of code in output dir if dir doesn't exist
     66# creates copy of code in output dir
    6067if myid == 0:
    61     if access(project.output_run_time_dir,F_OK) == 0:
    62         print 'project.output_run_time_dir',dirname(project.output_run_time_dir)
    63         mkdir (project.output_run_time_dir,0777)
    64     copy (dirname(project.__file__) +sep+ project.__name__+'.py',
    65              project.output_run_time_dir + project.__name__+'.py')
    66     copy (__file__, project.output_run_time_dir + basename(__file__))
    67     print 'project.output_run_time_dir',project.output_run_time_dir
     68    copy_code_files(project.output_run_time_dir,__file__,
     69                 dirname(project.__file__)+sep+ project.__name__+'.py' )
    6870barrier()
    69 #normal screen output is stored in
    70 screen_output_name = project.output_run_time_dir + "screen_output_%d_%d.txt" %(myid,numprocs)
    71 screen_error_name = project.output_run_time_dir + "screen_error_%d_%d.txt" %(myid,numprocs)
    72 
    73 #used to catch screen output to file
    74 sys.stdout = Screen_Catcher(screen_output_name)
    75 sys.stderr = Screen_Catcher(screen_error_name)
     71start_screen_catcher(project.output_run_time_dir, myid, numprocs)
    7672
    7773print 'USER: ', project.user
    78 
     74#sys.exit()
    7975#--------------------------------------------------------------------------
    8076# Create the triangular mesh based on overall clipping polygon with a
    … …  
    8581
    8682if myid == 0:
    87     if access(project.meshes_time_dir,F_OK) == 0:
    88         mkdir(project.meshes_time_dir)
     83#    if access(project.meshes_time_dir,F_OK) == 0:
     84#        mkdir(project.meshes_time_dir)
     85    if access(project.meshes_dir,F_OK) == 0:
     86        mkdir(project.meshes_dir)
    8987    print 'start create mesh from regions'
    9088    interior_regions = [#[project.karratha_polygon, 25000],
    9189#                    [project.cipma_polygon, 1000],
    92                     [project.poly_pipeline, 1000],
    93                     [project.poly_facility, 10000]]   
     90                    [project.poly_pipeline, 5000],
     91                    [project.poly_facility, 500]]   
    9492#    meshes_dir_name = project.meshes_dir_name + '.msh'
    9593
    … …  
    9997                         maximum_triangle_area=100000,
    10098                         interior_regions=interior_regions,
    101                          filename=meshes_time_dir_name,
     99#                         filename=meshes_time_dir_name,
     100                         filename=meshes_dir_name,
    102101                         use_cache=True,
    103102                         verbose=True)
    … …  
    110109#-------------------------------------------------------------------------
    111110print 'Setup computational domain'
    112 domain = Domain(meshes_time_dir_name, use_cache=True, verbose=True)
     111#domain = Domain(meshes_time_dir_name, use_cache=True, verbose=True)
     112domain = Domain(meshes_dir_name, use_cache=True, verbose=True)
    113113print domain.statistics()
    114114
    … …  
    128128#import sys; sys.exit()
    129129
    130 if access(project.boundaries_time_dir,F_OK) == 0:
    131     mkdir (project.boundaries_time_dir)
     130#if access(project.boundaries_time_dir,F_OK) == 0:
     131#    mkdir (project.boundaries_time_dir)
    132132# put above distribute
     133print 'boundary file is: ',boundaries_dir_name
    133134from caching import cache
    134 cache(ferret2sww,
    135       (boundaries_in_dir_name,
    136        boundaries_time_dir_name),
    137       {'verbose': True,
    138        'minlat': project.south_boundary,
    139        'maxlat': project.north_boundary,
    140        'minlon': project.west_boundary,
    141        'maxlon': project.east_boundary,
    142 #       'minlat': project.south,
    143 #       'maxlat': project.north,
    144 #       'minlon': project.west,
    145 #       'maxlon': project.east,
    146        'mint': 0, 'maxt': 35000,
    147        'origin': domain.geo_reference.get_origin(),
    148        'mean_stage': project.tide,
    149 #       'zscale': 1,                 #Enhance tsunami
    150        'fail_on_NaN': False},
    151        verbose = True,
    152        )
    153 
     135if myid == 0:
     136    cache(urs2sww,
     137          (boundaries_in_dir_name,
     138    #       boundaries_time_dir_name),
     139           boundaries_dir_name),
     140          {'verbose': True,
     141           'minlat': project.south_boundary,
     142           'maxlat': project.north_boundary,
     143           'minlon': project.west_boundary,
     144           'maxlon': project.east_boundary,
     145#           'minlat': project.south,
     146#           'maxlat': project.north,
     147#           'minlon': project.west,
     148#           'maxlon': project.east,
     149           'mint': 0, 'maxt': 35000,
     150           'origin': domain.geo_reference.get_origin(),
     151           'mean_stage': project.tide,
     152#           'zscale': 1,                 #Enhance tsunami
     153           'fail_on_NaN': False},
     154           verbose = True,
     155           )
     156barrier()
    154157
    155158
    … …  
    165168
    166169domain.set_quantity('elevation',
    167                     filename = project.topographies_dir + build_time + sep + project.combined_name + '.pts',
     170                    filename = project.topographies_dir + build_time + sep + project.combined_final_name + '.pts',
    168171                    use_cache = True,
    169172                    verbose = True,
    … …  
    187190domain.set_name(project.scenario_name)
    188191domain.set_datadir(project.output_run_time_dir)
    189 domain.set_default_order(2) #associate to the spatial order of the triangle
     192domain.set_default_order(2) # Apply second order scheme
    190193domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
     194domain.set_store_vertices_uniquely(False)
     195domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
     196domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
    191197
    192198#-------------------------------------------------------------------------
    … …  
    198204#boundariesname = project.boundaries_dir + '20061101_003322_build'+sep+boundaries_name
    199205#print'boundariesname',boundariesname
    200 Bf = File_boundary(boundaries_time_dir_name + '.sww',
     206#Bf = File_boundary(boundaries_time_dir_name + '.sww',
    201207
    202208#Bf = File_boundary(boundariesname + '.sww',
    203                    domain, time_thinning=20, use_cache=True, verbose=True)
     209
     210print 'domain id', id(domain)
     211Bf = File_boundary(boundaries_dir_name + '.sww',
     212                  domain, time_thinning=5, use_cache=True, verbose=True)
    204213
    205214print 'finished reading boundary file'
    … …  
    214223print'finish set boundary'
    215224
    216 
    217225#----------------------------------------------------------------------------
    218226# Evolve system through time
    … …  
    221229t0 = time.time()
    222230
    223 for t in domain.evolve(yieldstep = 60, finaltime = 28800):
    224 #for t in domain.evolve(yieldstep = 120, finaltime = 28800):
     231#for t in domain.evolve(yieldstep = 60, finaltime = 34000):
     232#    domain.write_time()
     233#    domain.write_boundary_statistics(tags = 'ocean')
     234
     235for t in domain.evolve(yieldstep = 120, finaltime = 9000):
     236#for t in domain.evolve(yieldstep = 120, finaltime = 130):
    225237    domain.write_time()
    226 #    domain.write_boundary_statistics(tags = 'ocean')     
     238    domain.write_boundary_statistics(tags = 'ocean')
     239    print 'time: ',time.time()
     240     
     241for t in domain.evolve(yieldstep = 60, finaltime = 28800
     242                       ,skip_initial_step = True):
     243    domain.write_time()
     244    domain.write_boundary_statistics(tags = 'ocean')   
     245
     246for t in domain.evolve(yieldstep = 120, finaltime = 34800
     247                       ,skip_initial_step = True):
     248    domain.write_time()
     249    domain.write_boundary_statistics(tags = 'ocean')   
    227250   
    228251print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracChangeset for help on using the changeset viewer.