Changeset 5652


Ignore:
Timestamp:
Aug 13, 2008, 4:40:45 PM (15 years ago)
Author:
kristy
Message:

updated all Geraldton files to reflex the new boundary file (copied from Perth) and new bathymetry

Location:
anuga_work/production/geraldton
Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/production/geraldton/build_geraldton.py

    r5150 r5652  
    2727from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts,start_screen_catcher, copy_code_files
    2828from anuga.geospatial_data.geospatial_data import *
     29from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files
    2930
    3031# Application specific imports
    … …  
    3637#------------------------------------------------------------------------------
    3738
    38 print 'time stamp: ',project.build_time
    39 
    4039copy_code_files(project.output_build_time_dir,__file__,
    4140               dirname(project.__file__)+sep+ project.__name__+'.py' )
    … …  
    4342start_screen_catcher(project.output_build_time_dir)
    4443
    45 print 'time stamp: ',project.build_time
    4644print 'USER: ', project.user, project.host
    47 
    4845
    4946#-------------------------------------------------------------------------------
    … …  
    6057island_in_dir_name = project.island_in_dir_name
    6158offshore_in_dir_name = project.offshore_in_dir_name
     59offshore_in_dir_name1 = project.offshore_in_dir_name1
     60offshore_in_dir_name2 = project.offshore_in_dir_name2
     61offshore_in_dir_name3 = project.offshore_in_dir_name3
    6262
    6363onshore_dir_name = project.onshore_dir_name
    … …  
    6565island_dir_name = project.island_dir_name
    6666offshore_dir_name = project.offshore_dir_name
     67offshore_dir_name1 = project.offshore_dir_name1
     68offshore_dir_name2 = project.offshore_dir_name2
     69offshore_dir_name3 = project.offshore_dir_name3
    6770
    6871# creates DEM from asc data
    … …  
    7275
    7376#creates pts file for onshore DEM
    74 print "creates pts file for onshore DEM"
    75 dem2pts(onshore_dir_name,
    76         use_cache=True,
    77         verbose=True)
     77dem2pts(onshore_dir_name, use_cache=True, verbose=True)
    7878
    7979#creates pts file for island DEM
    … …  
    8282print'create Geospatial data1 objects from topographies'
    8383G1 = Geospatial_data(file_name = onshore_dir_name + '.pts',verbose=True)
    84 print'finished reading in %s.pts' %onshore_dir_name
    8584G2 = Geospatial_data(file_name = coast_in_dir_name + '.txt',verbose=True)
    8685G3 = Geospatial_data(file_name = island_dir_name + '.pts',verbose=True)
    87 G4 = Geospatial_data(file_name = offshore_in_dir_name + '.txt',verbose=True)
     86print'create Geospatial data1 objects from bathymetry'
     87G_off = Geospatial_data(file_name = offshore_in_dir_name + '.txt',verbose=True)
     88G_off1 = Geospatial_data(file_name = offshore_in_dir_name1 + '.txt',verbose=True)
     89G_off2 = Geospatial_data(file_name = offshore_in_dir_name2 + '.txt',verbose=True)
     90G_off3 = Geospatial_data(file_name = offshore_in_dir_name3 + '.txt',verbose=True)
    8891print'finished reading in files'
    8992
    90 #G_onshore_clip = G1.clip(project.poly_all, verbose=True)
    91 #print 'finished clip'
    92 #G_onshore_clip_clipout=G_onshore_clip.clip_outside(project.poly_cbd, verbose=True)
    93 #print 'finished clipout'
    94 ## reduce resolution by 10 time (1/10) of the original file. This was determined
    95 ## as the resolution for this region is going to be 20000 and the asc grid cellsize is
    96 ## 10m there (100). There is 2 20000m triangles in a 200m x 200m grid sqrt of 40000 is about 200,
    97 ## a 200m x 200m area would contain 400 points at a 10m grid spacing and we only need 4. So 4/400 is 1/100 or 0.01
    98 ## 1/100 reduction is possible but a 1/10 is enough to help with file size and loading
    99 #G_onshore_clip_clipout_reduced=G_onshore_clip_clipout.split(0.1, verbose=True)
    100 #
    101 #G_cbd = G1.clip(projec.poly_cbd,verbose=True)
    102 ##G_cbd_reduced = G_cbd.split(0.25, verbose=True)
    103 #
    104 #G_bathy = G4.clip(project.poly_all)
    105 #G_island = G3.clip(project.poly_all)
    106 #G_coast = G2.clip(project.poly_all)
    107 #
    108 #print'add all geospatial objects'
    109 #G = G_onshore_clip_clipout_reduced + G_cbd_reduced + G_bathy + G_island + G_coast
    110 
    111 G = G1 + G2 + G3 + G4
    112 
    113 G_clip = G.clip(project.poly_all, verbose=True)
     93print'add all geospatial objects'
     94G = G1 + G2 + G3 + G_off + G_off1 + G_off2 + G_off3 + G_off4
    11495
    11596print'clip combined geospatial object by bounding polygon'
    116 #G_clipped = G.clip(project.bounding_polygon)
    117 #FIXME: add a clip function to pts
    118 #print'shape of clipped data', G_clipped.get_data_points().shape
     97G_clip = G.clip(project.poly_all, verbose=True)
    11998
    12099print'export combined DEM file'
    121100if access(project.topographies_dir,F_OK) == 0:
    122101    mkdir (project.topographies_dir)
    123 G_clip.export_points_file(project.combined_dir_name + '.txt')
    124 #G_clipped.export_points_file(project.combined_dir_name + '.xya')
    125 G_small,G_other = G_clip.split(0.01, verbose=True)
    126 G_small.export_points_file(project.combined_dir_name_small + '.txt')
     102G_clip.export_points_file(project.combined_dir_name + '.pts')
    127103
    128 '''
    129 print'project.combined_dir_name + 1.xya',project.combined_dir_name + '1.xya'
    130 G_all=Geospatial_data(file_name = project.combined_dir_name + '1.xya')
    131 print'split'
    132 G_all_1, G_all_2 = G_all.split(.10)
    133 print'export 1'
    134 G_all_1.export_points_file(project.combined_dir_name+'_small1' + '.xya')
    135 print'export 2'
    136 G_all_2.export_points_file(project.combined_dir_name+'_other1' + '.xya')
    137 
    138 
    139 #-------------------------------------------------------------------------
    140 # Convert URS to SWW file for boundary conditions
    141 #-------------------------------------------------------------------------
    142 print 'starting to create boundary conditions'
    143 boundaries_in_dir_name = project.boundaries_in_dir_name
    144 
    145 from anuga.shallow_water.data_manager import urs2sww
    146 
    147 print 'minlat=project.north_boundary, maxlat=project.south_boundary',project.north_boundary, project.south_boundary
    148 print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary
    149 
    150 #import sys; sys.exit()
    151 
    152 #if access(project.boundaries_dir,F_OK) == 0:
    153 #    mkdir (project.boundaries_dir)
    154 
    155 from caching import cache
    156 cache(urs2sww,
    157       (boundaries_in_dir_name,
    158        project.boundaries_dir_name1),
    159       {'verbose': True,
    160        'minlat': project.south_boundary,
    161        'maxlat': project.north_boundary,
    162        'minlon': project.west_boundary,
    163        'maxlon': project.east_boundary,
    164 #       'minlat': project.south,
    165 #       'maxlat': project.north,
    166 #       'minlon': project.west,
    167 #       'maxlon': project.east,
    168        'mint': 0, 'maxt': 40000,
    169 #       'origin': domain.geo_reference.get_origin(),
    170        'mean_stage': project.tide,
    171 #       'zscale': 1,                 #Enhance tsunami
    172        'fail_on_NaN': False},
    173        verbose = True,
    174        )
    175 #       dependencies = source_dir + project.boundary_basename + '.sww')
    176 
    177 '''
     104import sys
     105sys.exit()
    178106
    179107
    … …  
    182110
    183111
    184 

  • anuga_work/production/geraldton/project.py

    r5381 r5652  
    1010#from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
    1111from anuga.utilities.system_tools import get_user_name, get_host_name
     12from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
     13from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
    1214
    1315# file and system info
    … …  
    1517
    1618home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent dir   
     19muxhome = getenv('MUXHOME')
    1720user = get_user_name()
    1821host = get_host_name()
     22
    1923# INUNDATIONHOME is the inundation directory, not the data directory.
    20 
    21 #needed when running using mpirun, mpirun doesn't inherit umask from .bashrc
    22 umask(002)
    2324
    2425#time stuff
    2526time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
     27gtime = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
    2628build_time = time+'_build'
    2729run_time = time+'_run'
    28 
    29 tide = 0.75 #??? must check!!!
     30print 'gtime: ', gtime
    3031
    3132#Making assumptions about the location of scenario data
    … …  
    3435scenario = 'geraldton_tsunami_scenario'
    3536
    36 #Maybe will try to make project a class to allow these parameters to be passed in.
     37
     38tide = 0.75 #??? must check!!!
     39
    3740alpha = 0.1
    3841friction=0.01
    39 starttime=10000
    40 midtime=21600
    41 #finaltime=25000
    42 finaltime=10000
    43 export_cellsize=50
     42starttime=0
     43finaltime=80000
     44export_cellsize=25
    4445setup='final'
    45 source='other'
     46source='polyline'
    4647
    4748
    … …  
    6263    yieldstep=60
    6364
    64 dir_comment='_'+setup+'_'+str(tide)+'_'+str(source)+'_'+str(user)
     65dir_comment='_'+setup+'_'+str(tide)+'_'+str(source)+'_'+ 'alpha' +str(alpha)+'_'+str(user)
    6566
    6667
    … …  
    7071
    7172# AHO + DPI data
    72 coast_name = 'coastline'
    73 offshore_name = 'geraldton_bathy'
     73coast_name = 'XYcoastline_KVP'
     74offshore_name = 'Geraldton_bathymetry'
     75offshore_name1 = 'DPI_Data'
     76offshore_name2 = 'grid250'
     77offshore_name3 = 'Top_bathymetry'
     78
    7479
    7580#final topo name
    … …  
    8287topographies_dir = home+state+sep+scenario+sep+'anuga'+sep+'topographies'+sep
    8388
    84 #input topo file location
     89##input topo file location
    8590onshore_in_dir_name = topographies_in_dir + onshore_name
    8691island_in_dir_name = topographies_in_dir + island_name
    … …  
    8893coast_in_dir_name = topographies_in_dir + coast_name
    8994offshore_in_dir_name = topographies_in_dir + offshore_name
     95offshore_in_dir_name1 = topographies_in_dir + offshore_name1
     96offshore_in_dir_name2 = topographies_in_dir + offshore_name2
     97offshore_in_dir_name3 = topographies_in_dir + offshore_name3
    9098
    9199onshore_dir_name = topographies_dir + onshore_name
    … …  
    93101coast_dir_name = topographies_dir + coast_name
    94102offshore_dir_name = topographies_dir + offshore_name
     103offshore_dir_name1 = topographies_dir + offshore_name1
     104offshore_dir_name2 = topographies_dir + offshore_name2
     105offshore_dir_name3 = topographies_dir + offshore_name3
    95106
    96107#final topo files
    … …  
    104115tide_dir = anuga_dir+'tide_data'+sep
    105116
    106 ##if source =='dampier':
    107 ##    boundaries_name = 'broome_3854_17042007' #Dampier gun
    108 ##    boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'dampier'+sep+'1_10000'+sep
    109 ##
    110 ##if source=='onslow':
    111 ##    boundaries_name = 'broome_3859_16052007' #onslow_hedland_broome gun
    112 ##    boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'onslow_hedland_broome'+sep+'1_10000'+sep
    113 ##   
    114 ##if source=='exmouth':
    115 ##    boundaries_name = 'broome_3103_18052007' #exmouth gun
    116 ##    boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'exmouth'+sep+'1_10000'+sep
    117 ##
    118 ##if source=='other':
    119 ##    boundaries_name = 'other' #exmouth gun
    120 ##    boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'exmouth'+sep+'1_10000'+sep
    121 ##
    122 ##
    123 ###boundaries locations
    124 ##boundaries_in_dir_name = boundaries_in_dir + boundaries_name
    125 boundaries_in_dir_name = anuga_dir + sep + 'boundaries' + sep + scenario_name + '49'
     117#boundaries_source = '1'
     118
     119#boundaries locations
    126120boundaries_dir = anuga_dir+'boundaries'+sep
    127121boundaries_dir_name = boundaries_dir + scenario_name
     122boundaries_dir_mux = muxhome
    128123
    129124#output locations
    130125output_dir = anuga_dir+'outputs'+sep
    131 output_build_time_dir = output_dir+build_time+dir_comment+sep
    132 output_run_time_dir = output_dir +run_time+dir_comment+sep
    133 #output_run_time_dir = output_dir+sep+run_time+sep
     126output_build_time_dir = anuga_dir+'outputs'+sep+build_time+dir_comment+sep
     127output_run_time_dir = anuga_dir+'outputs'+sep+run_time+dir_comment+sep
    134128output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
     129
     130vertex_filename = output_run_time_dir + 'mesh_vertex.csv'
    135131
    136132#gauges
    … …  
    147143
    148144###############################
    149 # Domain definitions
     145# Interior region definitions
    150146###############################
    151147
    152 from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
     148# bounding polygon for study area
     149poly_all = read_polygon(polygons_dir+'poly_all.csv')
     150res_poly_all = 100000*res_factor
    153151
     152poly_neg20_pos20 = read_polygon(polygons_dir+'neg20_pos20.csv')
     153res_neg20_pos20 = 20000*res_factor
     154
     155poly_CBD_1km= read_polygon(polygons_dir+'CBD_1km.csv')
     156res_CBD_1km = 800*res_factor
     157
     158poly_cbd = read_polygon(polygons_dir+'CBD_500m.csv')
     159res_cbd = 500*res_factor
     160
     161poly_wallabi = read_polygon(polygons_dir+'island_wallabi_poly.csv')
     162res_wallabi = 1000*res_factor
     163
     164poly_dingiville = read_polygon(polygons_dir+'island_dingiville_poly.csv')
     165res_dingiville = 1000*res_factor
     166
     167poly_pelsaert = read_polygon(polygons_dir+'island_pelsaert_poly.csv')
     168res_pelsaert = 1000*res_factor
     169
     170
     171#plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3],figname='boundingpoly2',verbose=False)
     172
     173interior_regions = [[poly_neg20_pos20,res_neg20_pos20],[poly_CBD_1km,res_CBD_1km],
     174                    [poly_cbd,res_cbd],[poly_wallabi,res_wallabi],
     175                    [poly_dingiville,res_dingiville],[poly_pelsaert, res_pelsaert]]
     176
     177trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
     178print 'min number triangles', trigs_min
     179
     180poly_mainland=read_polygon(topographies_in_dir +'initial_condition.csv')
    154181
    155182###################################################################
    … …  
    157184###################################################################
    158185
    159 # exporting asc grid
    160 eastingmin = 340000
    161 eastingmax = 350000
    162 northingmin = 6273400
    163 northingmax = 6277700
    164186
    165 ###############################
    166 # Interior region definitions
    167 ###############################
     187# Geraldton CBD extract ascii grid
     188xminCBD = 262700
     189xmaxCBD = 267600
     190yminCBD = 6811500
     191ymaxCBD = 6816400
    168192
    169 #digitized polygons
    170 # bounding polygon for study area
    171 poly_all = read_polygon(polygons_dir+'poly_all_z50.csv')
    172 res_poly_all = 100000*res_factor
    173 
    174 poly_neg20_pos20 = read_polygon(polygons_dir+'neg20_pos20_points.csv')
    175 res_neg20_pos20 = 20000*res_factor
    176 
    177 poly_neg10_pos10= read_polygon(polygons_dir+'neg10_pos10_points.csv')
    178 res_neg10_pos10 = 2000*res_factor
    179 
    180 poly_cbd = read_polygon(polygons_dir+'cbd_500m.csv')
    181 res_cbd = 500*res_factor
    182 
    183 poly_wallabi = read_polygon(polygons_dir+'island_wallabi_poly_points.csv')
    184 res_wallabi = 20000*res_factor
    185 
    186 poly_dingiville = read_polygon(polygons_dir+'island_dingiville_poly_points.csv')
    187 res_dingiville = 20000*res_factor
    188 
    189 poly_pelsaert = read_polygon(polygons_dir+'island_pelsaert_poly_points.csv')
    190 res_pelsaert = 20000*res_factor
    191 
    192 
    193 #plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3],figname='boundingpoly2',verbose=False)
    194 
    195 interior_regions = [[poly_neg20_pos20,res_neg20_pos20],[poly_neg10_pos10,res_neg10_pos10],
    196                     [poly_cbd,res_cbd],[poly_wallabi,res_wallabi],
    197                     [poly_dingiville,res_dingiville],[poly_pelsaert, res_pelsaert]]
    198 
    199 boundary_tags={'back': [2, 3, 4, 5],
    200                'side': [1, 6], 'ocean': [0, 7, 8, 9, 10]}
    201 
    202 trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
    203 
    204 poly_mainland=read_polygon(polygons_dir+'initial_condition.csv')
    205 
    206 print 'min number triangles', trigs_min
    207 
    208 

  • anuga_work/production/geraldton/run_geraldton.py

    r5442 r5652  
    1717# Standard modules
    1818from os import sep
     19import os
    1920from os.path import dirname, basename
    2021from os import mkdir, access, F_OK
    … …  
    3031from anuga.shallow_water import Field_boundary
    3132from Numeric import allclose
    32 from anuga.shallow_water.data_manager import export_grid
     33from anuga.shallow_water.data_manager import export_grid, create_sts_boundary
    3334
    3435from anuga.pmesh.mesh_interface import create_mesh_from_regions
    3536from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
    36 from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
     37#from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
    3738from anuga_parallel.parallel_abstraction import get_processor_name
    3839from anuga.caching import myhash
    3940from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
    4041from anuga.fit_interpolate.benchmark_least_squares import mem_usage
     42from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
     43from anuga.geospatial_data.geospatial_data import find_optimal_smoothing_parameter
     44from Scientific.IO.NetCDF import NetCDFFile
    4145
    4246# Application specific imports
    4347import project                 # Definition of file names and polygons
     48numprocs = 1
     49myid = 0
    4450
    4551def run_model(**kwargs):
    … …  
    6268        store_parameters(**kwargs)
    6369
    64     barrier()
     70   # barrier()
    6571
    6672    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
    … …  
    6874    print "Processor Name:",get_processor_name()
    6975
    70     # filenames
    71 #    meshes_dir_name = project.meshes_dir_name+'.msh'
    72 
    73     # creates copy of code in output dir
    74     print 'min triangles', project.trigs_min,
    75     print 'Note: This is generally about 20% less than the final amount'
    76 
     76   
     77    #-----------------------------------------------------------------------
     78    # Domain definitions
     79    #-----------------------------------------------------------------------
     80
     81    # Read in boundary from ordered sts file
     82    urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir,project.scenario_name))
     83
     84    # Reading the landward defined points, this incorporates the original clipping
     85    # polygon minus the 100m contour
     86    landward_bounding_polygon = read_polygon(project.polygons_dir+'landward_boundary.txt')
     87
     88    # Combine sts polyline with landward points
     89    bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
     90   
     91    # counting segments
     92    N = len(urs_bounding_polygon)-1
     93    boundary_tags={'back': [N+2,N+3], 'side': [N,N+1,N+4],'ocean': range(N)}
     94
     95   
    7796    #--------------------------------------------------------------------------
    7897    # Create the triangular mesh based on overall clipping polygon with a
    … …  
    88107        print 'start create mesh from regions'
    89108
    90         create_mesh_from_regions(project.poly_all,
    91                              boundary_tags=project.boundary_tags,
     109        create_mesh_from_regions(bounding_polygon,
     110                             boundary_tags=boundary_tags,
    92111                             maximum_triangle_area=project.res_poly_all,
    93112                             interior_regions=project.interior_regions,
    94113                             filename=project.meshes_dir_name+'.msh',
    95                              use_cache=False,
     114                             use_cache=True,
    96115                             verbose=True)
    97     barrier()
    98    
     116   # barrier()
     117
     118##        covariance_value,alpha = find_optimal_smoothing_parameter (data_file= kwargs['elevation_file'],
     119##                                alpha_list=[0.001, 0.01, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5],
     120##                                mesh_file = project.meshes_dir_name+'.msh')
     121##        print 'optimal alpha', covariance_value,alpha       
    99122
    100123    #-------------------------------------------------------------------------
    … …  
    103126    print 'Setup computational domain'
    104127
    105     #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
    106     #above don't work
    107128    domain = Domain(project.meshes_dir_name+'.msh', use_cache=False, verbose=True)
    108129    print 'memory usage before del domain',mem_usage()
    … …  
    123144        from polygon import Polygon_function
    124145        #following sets the stage/water to be offcoast only
    125         IC = Polygon_function( [(project.poly_mainland, -1.0)], default = kwargs['tide'],
     146        IC = Polygon_function( [(project.poly_mainland, 0)], default = kwargs['tide'],
    126147                                 geo_reference = domain.geo_reference)
    127148        domain.set_quantity('stage', IC)
    128 #        domain.set_quantity('stage', kwargs['tide'])
     149        #domain.set_quantity('stage',kwargs['tide'] )
    129150        domain.set_quantity('friction', kwargs['friction'])
    130151       
    131         print 'Start Set quantity',kwargs['bathy_file']
     152        print 'Start Set quantity',kwargs['elevation_file']
    132153
    133154        domain.set_quantity('elevation',
    134                             filename = kwargs['bathy_file'],
    135                             use_cache = True,
     155                            filename = kwargs['elevation_file'],
     156                            use_cache = False,
    136157                            verbose = True,
    137158                            alpha = kwargs['alpha'])
    138159        print 'Finished Set quantity'
    139     barrier()
    140 
     160    #barrier()
     161
     162##    #------------------------------------------------------
     163##    # Create x,y,z file of mesh vertex!!!
     164##    #------------------------------------------------------
     165##        coord = domain.get_vertex_coordinates()
     166##        depth = domain.get_quantity('elevation')
     167##       
     168##        # Write vertex coordinates to file
     169##        filename=project.vertex_filename
     170##        fid=open(filename,'w')
     171##        fid.write('x (m), y (m), z(m)\n')
     172##        for i in range(len(coord)):
     173##            pt=coord[i]
     174##            x=pt[0]
     175##            y=pt[1]
     176##            z=depth[i]
     177##            fid.write('%.6f,%.6f,%.6f\n' %(x, y, z))
     178##
    141179
    142180    #------------------------------------------------------
    … …  
    157195    domain.set_store_vertices_uniquely(False)
    158196    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
     197    domain.tight_slope_limiters = 1
    159198    print 'domain id', id(domain)
    160 
    161199
    162200    #-------------------------------------------------------------------------
    … …  
    165203    print 'Available boundary tags', domain.get_boundary_tags()
    166204    print 'domain id', id(domain)
    167     #print 'Reading Boundary file',project.boundaries_dir_namea + '.sww'
    168     print'set_boundary'
    169 
     205   
     206    boundary_urs_out=project.boundaries_dir_name
     207   
     208    print 'Available boundary tags', domain.get_boundary_tags()
     209    Bf = Field_boundary(boundary_urs_out+'.sts',  # Change from file_boundary
     210                   domain, mean_stage= project.tide,
     211                   time_thinning=1,
     212                   use_cache=True,
     213                   verbose = True,
     214                   boundary_polygon=bounding_polygon)
     215
     216   
    170217    Br = Reflective_boundary(domain)
    171218    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
    172     Bo = Dirichlet_boundary([kwargs['tide']+5.0,0,0])
    173 
    174     if project.source != 'other':
    175         Bf = Field_boundary(kwargs['boundary_file'],
    176                     domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'],
    177                     use_cache=True, verbose=True)
    178         print 'finished reading boundary file'
    179         domain.set_boundary({'back': Br,
    180                              'side': Bd,
    181                              'ocean': Bf})
    182     else:
    183         print 'set ocean'               
    184         domain.set_boundary({'back': Br,
    185                              'side': Bd,
    186                              'ocean': Bd})
    187 
    188 #    kwargs['input_start_time']=domain.starttime
    189 
    190 
     219
     220    fid = NetCDFFile(boundary_urs_out+'.sts', 'r')    #Open existing file for read
     221    sts_time=fid.variables['time'][:]+fid.starttime
     222    tmin=min(sts_time)
     223    tmax=max(sts_time)
     224    fid.close()
     225
     226    print 'Boundary end time ', tmax-tmin
     227   
     228##    Bf = Field_boundary(kwargs['boundary_file'],
     229##                domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'],
     230##                use_cache=False, verbose=True)
     231
     232    domain.set_boundary({'back': Br,
     233                         'side': Bd,
     234                         'ocean': Bf})
     235
     236    kwargs['input_start_time']=domain.starttime
    191237
    192238    print'finish set boundary'
    … …  
    197243    t0 = time.time()
    198244
    199     for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime']):
     245    for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime']
     246                       ,skip_initial_step = False):
    200247        domain.write_time()
    201         domain.write_boundary_statistics(tags = 'ocean')     
    202 
    203         if allclose(t, 240):
    204             domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bo})
    205 
    206         if allclose(t, 1440):
    207             domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bd})
    208 
    209 #    for t in domain.evolve(yieldstep = kwargs['yieldstep'], finaltime = kwargs['midtime']
    210 #                       ,skip_initial_step = True):
    211 #        domain.write_time()
    212 #        domain.write_boundary_statistics(tags = 'ocean')     
    213 #   
    214 #    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime']
    215 #                       ,skip_initial_step = True):
    216 #        domain.write_time()
    217 #        domain.write_boundary_statistics(tags = 'ocean')   
    218 
     248        domain.write_boundary_statistics(tags = 'ocean')
     249
     250        if t >= tmax-tmin:
     251            print 'changed to tide boundary condition at ocean'
     252            domain.set_boundary({'ocean': Bd})
     253           
    219254    x, y = domain.get_maximum_inundation_location()
    220255    q = domain.get_maximum_inundation_elevation()
    … …  
    229264    if myid==0:
    230265        store_parameters(**kwargs)
    231     barrier
     266   # barrier
    232267   
    233268    print 'memory usage before del domain1',mem_usage()
    234269   
    235 def export_model(**kwargs):   
    236     #store_parameters(**kwargs)
    237    
    238 #    print 'memory usage before del domain',mem_usage()
    239     #del domain
    240     print 'memory usage after del domain',mem_usage()
    241    
    242     swwfile = kwargs['output_dir']+kwargs['aa_scenario_name']
    243     print'swwfile',swwfile
    244    
    245     export_grid(swwfile, extra_name_out = 'town',
    246             quantities = ['speed','depth','elevation','stage'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
    247             #quantities = ['speed','depth'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
    248             timestep = None,
    249             reduction = max,
    250             cellsize = kwargs['export_cellsize'],
    251             NODATA_value = -9999,
    252             easting_min = project.eastingmin,
    253             easting_max = project.eastingmax,
    254             northing_min = project.northingmin,
    255             northing_max = project.northingmax,
    256             verbose = False,
    257             origin = None,
    258             datum = 'GDA94',
    259             format = 'asc')
    260            
    261     inundation_damage(swwfile+'.sww', project.buildings_filename,
    262                       project.buildings_filename_out)
    263270   
    264271#-------------------------------------------------------------
    … …  
    272279    kwargs['starttime']=project.starttime
    273280    kwargs['yieldstep']=project.yieldstep
    274     kwargs['midtime']=project.midtime
    275281    kwargs['finaltime']=project.finaltime
    276282   
    277283    kwargs['output_dir']=project.output_run_time_dir
    278 #    kwargs['bathy_file']=project.combined_dir_name+'.txt'
    279     kwargs['bathy_file']=project.combined_dir_name_small + '.txt'
     284    kwargs['elevation_file']=project.combined_dir_name+'.pts'
    280285    kwargs['boundary_file']=project.boundaries_in_dir_name + '.sww'
    281286    kwargs['file_name']=project.home+'detail.csv'
    … …  
    296301    if myid==0:
    297302        export_model(**kwargs)
    298     barrier
     303    #barrier
Note: See TracChangeset for help on using the changeset viewer.