Changeset 5755


Ignore:
Timestamp:
Sep 9, 2008, 2:18:01 PM (17 years ago)
Author:
kristy
Message:

localised all python scripts

Location:
anuga_work/production/carnarvon
Files:
3 added
3 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/production/carnarvon/build_carnarvon.py

    r5005 r5755  
    2727from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
    2828from anuga.geospatial_data.geospatial_data import *
    29 from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
     29from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
    3030
    3131# Application specific imports
    … …  
    3737#------------------------------------------------------------------------------
    3838
    39 start_screen_catcher(project.output_build_time_dir)
    40 
    41 print 'time stamp: ',project.gtime
    42 print 'USER: ', project.user
    43 
    44 
    4539copy_code_files(project.output_build_time_dir,__file__,
    4640               dirname(project.__file__)+sep+ project.__name__+'.py' )
    4741
     42start_screen_catcher(project.output_build_time_dir)
    4843
     44print 'USER: ', project.user
    4945
    5046#-------------------------------------------------------------------------------
    … …  
    5551# Fine pts file to be clipped to area of interest
    5652#-------------------------------------------------------------------------------
     53print"project.poly_all",project.poly_all
     54print"project.combined_dir_name",project.combined_dir_name
    5755
    5856# topography directory filenames
    5957onshore_in_dir_name = project.onshore_in_dir_name
    6058coast_in_dir_name = project.coast_in_dir_name
    61 #island_in_dir_name = project.island_in_dir_name
    6259offshore_in_dir_name = project.offshore_in_dir_name
     60offshore_in_dir_name1 = project.offshore_in_dir_name1
     61offshore_in_dir_name2 = project.offshore_in_dir_name2
    6362
    6463onshore_dir_name = project.onshore_dir_name
    6564coast_dir_name = project.coast_dir_name
    66 #island_dir_name = project.island_dir_name
     65
    6766offshore_dir_name = project.offshore_dir_name
     67offshore_dir_name1 = project.offshore_dir_name1
     68offshore_dir_name2 = project.offshore_dir_name2
     69
    6870
    6971# creates DEM from asc data
    7072print "creates DEMs from ascii data"
    7173convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
    72 #convert_dem_from_ascii2netcdf(island_in_dir_name, basename_out=island_dir_name, use_cache=True, verbose=True)
    7374
    7475#creates pts file for onshore DEM
    … …  
    8687
    8788print'create Geospatial data1 objects from topographies'
    88 G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
    89 G2 = Geospatial_data(file_name = coast_in_dir_name + '.xya')
    90 #G3 = Geospatial_data(file_name = island_dir_name + '.pts')
    91 G_off = Geospatial_data(file_name = offshore_in_dir_name + '.xya')
     89G1 = Geospatial_data(file_name = onshore_dir_name + '.pts',verbose=True)
     90G2 = Geospatial_data(file_name = coast_in_dir_name + '.txt',verbose=True)
    9291
     92G_off = Geospatial_data(file_name = offshore_in_dir_name + '.txt',verbose=True)
     93G_off1 = Geospatial_data(file_name = offshore_in_dir_name1 + '.txt',verbose=True)
     94G_off2 = Geospatial_data(file_name = offshore_in_dir_name2 + '.txt',verbose=True)
    9395print'add all geospatial objects'
    94 G = G1 + G2 + G_off
     96G = G1 + G2 + G_off + G_off1 + G_off2
    9597
    9698print'clip combined geospatial object by bounding polygon'
    97 #G_clipped = G.clip(project.bounding_polygon)
    98 #FIXME: add a clip function to pts
    99 #print'shape of clipped data', G_clipped.get_data_points().shape
     99G_clipped = G.clip(project.poly_all)
    100100
    101101print'export combined DEM file'
    102102if access(project.topographies_dir,F_OK) == 0:
    103103    mkdir (project.topographies_dir)
    104 G.export_points_file(project.combined_dir_name + '.xya')
    105 #G_clipped.export_points_file(project.combined_dir_name + '.xya')
     104G_clipped.export_points_file(project.combined_dir_name + '.txt')
    106105
    107 '''
    108 print'project.combined_dir_name + 1.xya',project.combined_dir_name + '1.xya'
    109 G_all=Geospatial_data(file_name = project.combined_dir_name + '1.xya')
    110 print'split'
    111 G_all_1, G_all_2 = G_all.split(.10)
    112 print'export 1'
    113 G_all_1.export_points_file(project.combined_dir_name+'_small1' + '.xya')
    114 print'export 2'
    115 G_all_2.export_points_file(project.combined_dir_name+'_other1' + '.xya')
     106print'project.combined_dir_name + .txt',project.combined_dir_name + '.txt'
    116107
    117 
    118 #-------------------------------------------------------------------------
    119 # Convert URS to SWW file for boundary conditions
    120 #-------------------------------------------------------------------------
    121 print 'starting to create boundary conditions'
    122 boundaries_in_dir_name = project.boundaries_in_dir_name
    123 
    124 from anuga.shallow_water.data_manager import urs2sww
    125 
    126 print 'minlat=project.north_boundary, maxlat=project.south_boundary',project.north_boundary, project.south_boundary
    127 print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary
    128 
    129 #import sys; sys.exit()
    130 
    131 #if access(project.boundaries_dir,F_OK) == 0:
    132 #    mkdir (project.boundaries_dir)
    133 
    134 from caching import cache
    135 cache(urs2sww,
    136       (boundaries_in_dir_name,
    137        project.boundaries_dir_name1),
    138       {'verbose': True,
    139        'minlat': project.south_boundary,
    140        'maxlat': project.north_boundary,
    141        'minlon': project.west_boundary,
    142        'maxlon': project.east_boundary,
    143 #       'minlat': project.south,
    144 #       'maxlat': project.north,
    145 #       'minlon': project.west,
    146 #       'maxlon': project.east,
    147        'mint': 0, 'maxt': 40000,
    148 #       'origin': domain.geo_reference.get_origin(),
    149        'mean_stage': project.tide,
    150 #       'zscale': 1,                 #Enhance tsunami
    151        'fail_on_NaN': False},
    152        verbose = True,
    153        )
    154 #       dependencies = source_dir + project.boundary_basename + '.sww')
    155 
    156 '''
    157 
    158 
    159 
    160 
    161 
    162 
    163 
     108###-------------------------------------------------------------------------
     109### Convert URS to SWW file for boundary conditions
     110###-------------------------------------------------------------------------
     111##print 'starting to create boundary conditions'
     112##from anuga.shallow_water.data_manager import urs2sww, urs_ungridded2sww
     113##
     114##boundaries_in_dir_name = project.boundaries_in_dir_name
     115##print 'boundaries_in_dir_name',project.boundaries_in_dir_name
     116##
     117##
     118###import sys; sys.exit()
     119##
     120##urs_ungridded2sww(project.boundaries_in_dir_name, project.boundaries_in_dir_name,
     121##                  verbose=True, mint=4000, maxt=80000, zscale=1)
     122##

  • anuga_work/production/carnarvon/project.py

    r5381 r5755  
    33"""
    44
    5 from os import sep, environ, getenv, getcwd ,umask
     5from os import sep, environ, getenv, getcwd
    66from os.path import expanduser
    77import sys
    … …  
    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
    1416#---------------------------------
     17#codename = 'project.py'
    1518
    16 home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent dir   
     19home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent diruser = get_user_name()
     20muxhome = getenv('MUXHOME')
    1721user = get_user_name()
    1822host = get_host_name()
     23
    1924# 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)
    2325
    2426#time stuff
    2527time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
     28gtime = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
    2629build_time = time+'_build'
    2730run_time = time+'_run'
    28 
    29 tide = 0.6
     31print 'gtime: ', gtime
    3032
    3133#Making assumptions about the location of scenario data
    … …  
    3436scenario = 'carnarvon_tsunami_scenario'
    3537
    36 #Maybe will try to make project a class to allow these parameters to be passed in.
     38
     39tide = 0.0 #1.0
     40
    3741alpha = 0.1
    3842friction=0.01
    39 starttime=10000
    40 midtime=21600
    41 finaltime=25000
    42 export_cellsize=50
     43starttime=0
     44finaltime=1000
     45export_cellsize=25
    4346setup='final'
    44 source='other'
    45 
     47source='polyline'
    4648
    4749if setup =='trial':
    … …  
    6163    yieldstep=60
    6264
    63 dir_comment='_'+setup+'_'+str(tide)+'_'+str(source)+'_'+str(user)
     65dir_comment='_'+setup+'_'+str(tide)+'_'+str(source)+'_'+ 'alpha' +str(alpha)+'_'+str(user)
     66
     67# onshore data provided by WA DLI
     68onshore_name = 'bathy250_clipland' # original
     69
     70# AHO + DPI data + colin French coastline
     71coast_name = 'DPI_coastlineP'
     72offshore_name = 'Shark_Bay_Clip'
     73offshore_name1 = 'XYAHD_Clip'
     74offshore_name2 = 'DPI_data'
    6475
    6576
    66 # onshore data provided by WA DLI
    67 onshore_name = 'DLI_orthophoto_DEM' # original
    68 #islands
    69 
    70 # AHO + DPI data
    71 coast_name = '100k_coastline'
    72 offshore_name = 'Bathymetry'
    73 
    7477#final topo name
    75 combined_name ='carnarvon_combined_elevation.pts'
    76 combined_smaller_name = 'carnarvon_combined_elevation_smaller'
     78combined_name = scenario_name+'_combined_elevation'
     79combined_smaller_name = scenario_name+'_combined_elevation_smaller'
    7780
    7881anuga_dir = home+state+sep+scenario+sep+'anuga'+sep
    7982
    80 topographies_in_dir = home+sep+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep
    81 topographies_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'topographies'+sep
     83topographies_in_dir = home+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep
     84topographies_dir = anuga_dir+'topographies'+sep
    8285
    83 #input topo file location
     86# input topo file location
    8487onshore_in_dir_name = topographies_in_dir + onshore_name
    85 #island_in_dir_name = topographies_in_dir + island_name
    8688
    8789coast_in_dir_name = topographies_in_dir + coast_name
    8890offshore_in_dir_name = topographies_in_dir + offshore_name
     91offshore_in_dir_name1 = topographies_in_dir + offshore_name1
     92offshore_in_dir_name2 = topographies_in_dir + offshore_name2
    8993
    9094onshore_dir_name = topographies_dir + onshore_name
    91 #island_dir_name = topographies_dir + island_name
     95
    9296coast_dir_name = topographies_dir + coast_name
    9397offshore_dir_name = topographies_dir + offshore_name
     98offshore_dir_name1 = topographies_dir + offshore_name1
     99offshore_dir_name2 = topographies_dir + offshore_name2
    94100
    95101#final topo files
    96102combined_dir_name = topographies_dir + combined_name
    97 combined_smaller_dir_name = topographies_dir + combined_smaller_name
     103#combined_time_dir_name = topographies_time_dir + combined_name
     104combined_smaller_name_dir = topographies_dir + combined_smaller_name
     105#combined_time_dir_final_name = topographies_time_dir + combined_final_name
    98106
    99107meshes_dir = anuga_dir+'meshes'+sep
    … …  
    103111tide_dir = anuga_dir+'tide_data'+sep
    104112
    105 if source =='dampier':
    106     boundaries_name = 'broome_3854_17042007' #Dampier gun
    107     boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'dampier'+sep+'1_10000'+sep
    108113
    109 if source=='onslow':
    110     boundaries_name = 'broome_3859_16052007' #onslow_hedland_broome gun
    111     boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'onslow_hedland_broome'+sep+'1_10000'+sep
     114#boundaries_source = '1'
    112115   
    113 if source=='exmouth':
    114     boundaries_name = 'broome_3103_18052007' #exmouth gun
    115     boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'exmouth'+sep+'1_10000'+sep
     116if source=='polyline':
     117    boundaries_name = 'perth_3103_28052008' #polyline gun
     118    boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'polyline'+sep+'1_10000'+sep
    116119
    117 if source=='other':
    118     boundaries_name = 'other' #exmouth gun
    119     boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'exmouth'+sep+'1_10000'+sep
     120if source=='test':
     121    boundaries_name = 'other' #polyline
     122    boundaries_in_dir = anuga_dir+'boundaries'+sep
    120123
    121124
    122125#boundaries locations
    123 boundaries_in_dir_name = boundaries_in_dir + boundaries_name
     126boundaries_in_dir_name = boundaries_in_dir + boundaries_name 
    124127boundaries_dir = anuga_dir+'boundaries'+sep
    125 boundaries_dir_name = boundaries_dir + scenario_name
     128boundaries_dir_name = boundaries_dir + scenario_name # what it creates???
     129boundaries_dir_mux = muxhome
    126130
    127131#output locations
    128132output_dir = anuga_dir+'outputs'+sep
    129 output_build_time_dir = output_dir+build_time+'_'+source+sep
    130 #output_run_time_dir = output_dir +run_time+dir_comment+sep
    131 output_run_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+run_time+sep
     133output_build_time_dir = anuga_dir+'outputs'+sep+build_time+dir_comment+sep
     134output_run_time_dir = anuga_dir+'outputs'+sep+run_time+dir_comment+sep
    132135output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
    133136
     137vertex_filename = output_run_time_dir + 'mesh_vertex.csv'
     138
    134139#gauges
    135 gauge_name = '???.csv'
     140gauge_name = 'carnarvon.csv'
     141
    136142gauges_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'gauges'+sep
     143beach_gauges = gauges_dir + 'beach_gauges.csv'
    137144gauges_dir_name = gauges_dir + gauge_name
    138145
    139 buildings_filename = gauges_dir + 'carnarvon_res_Project.csv'
    140 buildings_filename_out = 'carnarvon_res_Project_modified.csv'
     146##buildings_filename = gauges_dir + 'Perth_resA.csv'
     147##buildings_filename_out = 'Perth_res_Project_modified.csv'
    141148
    142 community_filename = gauges_dir +''
    143 community_broome = gauges_dir + ''
     149###############################
     150# Interior region definitions
     151###############################
     152
     153#Initial bounding polygon for data clipping
     154poly_all = read_polygon(polygons_dir+'poly_all.csv')
     155res_poly_all = 100000*res_factor
     156
     157#Polygon designed
     158poly_internal_20 = read_polygon(polygons_dir+'Carnarvon20m.csv')
     159res_internal_20 = 25000*res_factor
     160
     161#Polygon designed to
     162poly_internal_5 = read_polygon(polygons_dir+'Carnarvon5m.csv')
     163res_internal_5 = 500*res_factor
     164
     165#Polygon designed to
     166poly_internal_10 = read_polygon(polygons_dir+'Carnarvon10m.csv')
     167res_internal_10 = 1500*res_factor
     168
     169#Polygon designed to incorporate
     170poly_island_20 = read_polygon(polygons_dir+'Island20m.csv')
     171res_island_20 = 50000*res_factor
    144172
    145173
    146 ###############################
    147 # Domain definitions
    148 ###############################
     174interior_regions = [[poly_internal_20,res_internal_20],[poly_internal_5,res_internal_5]
     175                     ,[poly_island_20,res_island_20],[poly_internal_10,res_internal_10]]
    149176
    150 from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
     177   
     178trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
     179print 'min number triangles', trigs_min
     180   
    151181
    152 # bounding polygon for study area
    153 poly_all = read_polygon(polygons_dir+'poly_all.csv')
    154 res_poly_all = 100000*res_factor
     182poly_mainland = read_polygon(polygons_dir+'initial_condition.csv')
    155183
    156184###################################################################
    … …  
    158186###################################################################
    159187
    160 # exporting asc grid
    161 eastingmin = 340000
    162 eastingmax = 350000
    163 northingmin = 6273400
    164 northingmax = 6277700
    165188
    166 ###############################
    167 # Interior region definitions
    168 ###############################
    169 
    170 #digitized polygons
    171 poly_carnarvon1 = read_polygon(polygons_dir+'neg20_pos10_polygon.csv')
    172 res_carnarvon1 = 10000*res_factor
    173 poly_carnarvon2 = read_polygon(polygons_dir+'neg5_pos5_poly_.csv')
    174 res_carnarvon2 = 500*res_factor
    175 
    176 #plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3],figname='boundingpoly2',verbose=False)
    177 
    178 interior_regions = [[poly_carnarvon1,res_carnarvon1],[poly_carnarvon2,res_carnarvon2]]
    179 
    180 boundary_tags={'back': [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 14],
    181                'side': [10,13], 'ocean': [11, 12]}
    182 
    183 trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
    184 
    185 poly_mainland=read_polygon(polygons_dir+'initial_condition_south.csv')
    186 
    187 print 'min number triangles', trigs_min
    188 
    189 

  • anuga_work/production/carnarvon/run_carnarvon.py

    r5442 r5755  
    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
    4144
    4245# Application specific imports
    4346import project                 # Definition of file names and polygons
     47numprocs = 1
     48myid = 0
    4449
    4550def run_model(**kwargs):
    … …  
    6267        store_parameters(**kwargs)
    6368
    64     barrier()
     69   # barrier()
    6570
    6671    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
    … …  
    6873    print "Processor Name:",get_processor_name()
    6974
    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 
     75   
     76    #-----------------------------------------------------------------------
     77    # Domain definitions
     78    #-----------------------------------------------------------------------
     79
     80    # Read in boundary from ordered sts file
     81    urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir,project.scenario_name))
     82
     83    # Reading the landward defined points, this incorporates the original clipping
     84    # polygon minus the 100m contour
     85    landward_bounding_polygon = read_polygon(project.boundaries_dir+'landward_bounding_polygon.txt')
     86
     87    # Combine sts polyline with landward points
     88    bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
     89   
     90    # counting segments
     91    N = len(urs_bounding_polygon)-1
     92    boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4],'ocean': range(N)}
     93
     94   
    7795    #--------------------------------------------------------------------------
    7896    # Create the triangular mesh based on overall clipping polygon with a
    … …  
    88106        print 'start create mesh from regions'
    89107
    90         create_mesh_from_regions(project.poly_all,
    91                              boundary_tags=project.boundary_tags,
     108        create_mesh_from_regions(bounding_polygon,
     109                             boundary_tags=boundary_tags,
    92110                             maximum_triangle_area=project.res_poly_all,
    93111                             interior_regions=project.interior_regions,
    94112                             filename=project.meshes_dir_name+'.msh',
    95                              use_cache=False,
     113                             use_cache=True,
    96114                             verbose=True)
    97     barrier()
    98    
     115   # barrier()
     116
     117##        covariance_value,alpha = find_optimal_smoothing_parameter (data_file= kwargs['elevation_file'],
     118##                                alpha_list=[0.001, 0.01, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5],
     119##                                mesh_file = project.meshes_dir_name+'.msh')
     120##        print 'optimal alpha', covariance_value,alpha       
    99121
    100122    #-------------------------------------------------------------------------
    … …  
    103125    print 'Setup computational domain'
    104126
    105     #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
    106     #above don't work
    107127    domain = Domain(project.meshes_dir_name+'.msh', use_cache=False, verbose=True)
    108128    print 'memory usage before del domain',mem_usage()
    … …  
    123143        from polygon import Polygon_function
    124144        #following sets the stage/water to be offcoast only
    125         IC = Polygon_function( [(project.poly_mainland, -1.0)], default = kwargs['tide'],
     145        IC = Polygon_function( [(project.poly_mainland, 0)], default = kwargs['tide'],
    126146                                 geo_reference = domain.geo_reference)
    127147        domain.set_quantity('stage', IC)
    128 #        domain.set_quantity('stage', kwargs['tide'])
     148        #domain.set_quantity('stage',kwargs['tide'] )
    129149        domain.set_quantity('friction', kwargs['friction'])
    130150       
    131         print 'Start Set quantity',kwargs['bathy_file']
     151        print 'Start Set quantity',kwargs['elevation_file']
    132152
    133153        domain.set_quantity('elevation',
    134                             filename = kwargs['bathy_file'],
    135                             use_cache = True,
     154                            filename = kwargs['elevation_file'],
     155                            use_cache = False,
    136156                            verbose = True,
    137157                            alpha = kwargs['alpha'])
    138158        print 'Finished Set quantity'
    139     barrier()
    140 
     159    #barrier()
     160
     161##    #------------------------------------------------------
     162##    # Create x,y,z file of mesh vertex!!!
     163##    #------------------------------------------------------
     164##        coord = domain.get_vertex_coordinates()
     165##        depth = domain.get_quantity('elevation')
     166##       
     167##        # Write vertex coordinates to file
     168##        filename=project.vertex_filename
     169##        fid=open(filename,'w')
     170##        fid.write('x (m), y (m), z(m)\n')
     171##        for i in range(len(coord)):
     172##            pt=coord[i]
     173##            x=pt[0]
     174##            y=pt[1]
     175##            z=depth[i]
     176##            fid.write('%.6f,%.6f,%.6f\n' %(x, y, z))
     177##
    141178
    142179    #------------------------------------------------------
    … …  
    157194    domain.set_store_vertices_uniquely(False)
    158195    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
    159     domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
     196    domain.tight_slope_limiters = 1
    160197    print 'domain id', id(domain)
    161 
    162198
    163199    #-------------------------------------------------------------------------
    … …  
    166202    print 'Available boundary tags', domain.get_boundary_tags()
    167203    print 'domain id', id(domain)
    168     #print 'Reading Boundary file',project.boundaries_dir_namea + '.sww'
    169 
    170     if project.source != 'other':
    171         Bf = Field_boundary(kwargs['boundary_file'],
    172                     domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'],
    173                     use_cache=True, verbose=True)
    174                    
    175     kwargs['input_start_time']=domain.starttime
    176 
    177     print 'finished reading boundary file'
     204   
     205    boundary_urs_out=project.boundaries_dir_name
    178206
    179207    Br = Reflective_boundary(domain)
    180208    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
    181     Bo = Dirichlet_boundary([kwargs['tide']+5.0,0,0])
    182 
    183 
    184     print'set_boundary'
     209   
     210    print 'Available boundary tags', domain.get_boundary_tags()
     211    Bf = Field_boundary(boundary_urs_out+'.sts',  # Change from file_boundary
     212                   domain, mean_stage= project.tide,
     213                   time_thinning=1,
     214                   default_boundary=Bd,
     215                   use_cache=True,
     216                   verbose = True,
     217                   boundary_polygon=bounding_polygon)
    185218
    186219    domain.set_boundary({'back': Br,
    187220                         'side': Bd,
    188221                         'ocean': Bf})
     222
     223    kwargs['input_start_time']=domain.starttime
     224
    189225    print'finish set boundary'
    190226
    … …  
    194230    t0 = time.time()
    195231
    196     for t in domain.evolve(yieldstep = 240, finaltime = kwargs['starttime']):
     232    for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime']
     233                       ,skip_initial_step = False):
    197234        domain.write_time()
    198         domain.write_boundary_statistics(tags = 'ocean')     
    199 
    200         if allclose(t, 120):
    201             domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bo})
    202 
    203         if allclose(t, 720):
    204             domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bd})
    205 
    206 #    for t in domain.evolve(yieldstep = kwargs['yieldstep'], finaltime = kwargs['midtime']
    207 #                       ,skip_initial_step = True):
    208 #        domain.write_time()
    209 #        domain.write_boundary_statistics(tags = 'ocean')     
    210 #   
    211 #    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime']
    212 #                       ,skip_initial_step = True):
    213 #        domain.write_time()
    214 #        domain.write_boundary_statistics(tags = 'ocean')   
     235        domain.write_boundary_statistics(tags = 'ocean')
    215236
    216237    x, y = domain.get_maximum_inundation_location()
    … …  
    226247    if myid==0:
    227248        store_parameters(**kwargs)
    228     barrier
     249   # barrier
    229250   
    230251    print 'memory usage before del domain1',mem_usage()
    231252   
    232 def export_model(**kwargs):   
    233     #store_parameters(**kwargs)
    234    
    235 #    print 'memory usage before del domain',mem_usage()
    236     #del domain
    237     print 'memory usage after del domain',mem_usage()
    238    
    239     swwfile = kwargs['output_dir']+kwargs['aa_scenario_name']
    240     print'swwfile',swwfile
    241    
    242     export_grid(swwfile, extra_name_out = 'town',
    243             quantities = ['speed','depth','elevation','stage'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
    244             #quantities = ['speed','depth'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
    245             timestep = None,
    246             reduction = max,
    247             cellsize = kwargs['export_cellsize'],
    248             NODATA_value = -9999,
    249             easting_min = project.eastingmin,
    250             easting_max = project.eastingmax,
    251             northing_min = project.northingmin,
    252             northing_max = project.northingmax,
    253             verbose = False,
    254             origin = None,
    255             datum = 'GDA94',
    256             format = 'asc')
    257            
    258     inundation_damage(swwfile+'.sww', project.buildings_filename,
    259                       project.buildings_filename_out)
    260253   
    261254#-------------------------------------------------------------
    … …  
    269262    kwargs['starttime']=project.starttime
    270263    kwargs['yieldstep']=project.yieldstep
    271     kwargs['midtime']=project.midtime
    272264    kwargs['finaltime']=project.finaltime
    273265   
    274266    kwargs['output_dir']=project.output_run_time_dir
    275     kwargs['bathy_file']=project.combined_dir_name
    276 #    kwargs['bathy_file']=project.combined_small_dir_name + '.pts'
     267    kwargs['elevation_file']=project.combined_dir_name+'.txt'
    277268    kwargs['boundary_file']=project.boundaries_in_dir_name + '.sww'
    278269    kwargs['file_name']=project.home+'detail.csv'
    … …  
    291282    run_model(**kwargs)
    292283     
    293     if myid==0:
    294         export_model(**kwargs)
    295     barrier
     284    #barrier
Note: See TracChangeset for help on using the changeset viewer.