Changeset 4153


Ignore:
Timestamp:
Jan 9, 2007, 3:55:58 PM (17 years ago)
Author:
nick
Message:

update to broome

Location:
anuga_work/production/broome_2006
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/production/broome_2006/project.py

    r3992 r4153  
    77import sys
    88from time import localtime, strftime, gmtime
    9 from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
     9from anuga.utilities.polygon import read_polygon, plot_polygons, is_inside_polygon, number_mesh_triangles
    1010#from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
     11from anuga.utilities.system_tools import get_user_name
    1112
    12 if sys.platform == 'win32':
    13     home = getenv('INUNDATIONHOME')
    14     user = getenv('USERPROFILE')
     13# file and system info
     14#---------------------------------
     15codename = 'project.py'
    1516
    16 else:   
    17     home = getenv('INUNDATIONHOME', sep+'d'+sep+'xrd'+sep+'gem'+sep+'2'+sep+'ramp'+sep+'risk_assessment_methods_project'+sep+'inundation')     
    18     user = getenv('LOGNAME')
    19     print 'USER:', user
     17home = getenv('INUNDATIONHOME') #Sandpit's parent dir   
     18user = get_user_name()
    2019
    2120# INUNDATIONHOME is the inundation directory, not the data directory.
    2221home += sep +'data'
    2322
     23#time stuff
     24time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
     25gtime = strftime('%Y%m%d_%H%M%S',gmtime()) #gets time for new dir
     26build_time = time+'_build'
     27run_time = time+'_run'
     28print 'gtime: ', gtime
     29
     30tide = 0
     31
    2432#Making assumptions about the location of scenario data
    2533state = 'western_australia'
    26 scenario_dir_name = 'broome_tsunami_scenario_2006'
     34scenario_name = 'broome'
     35scenario = 'broome_tsunami_scenario_2006'
    2736
    2837# onshore data provided by WA DLI
    2938onshore_name = 'dted2_z51' # original
    3039
    31 # AHO + DPI data
     40#island
     41#island_name = 'rott_dli_ext' # original
     42
     43# offshore
     44coast_name = 'coastline'
    3245offshore_name1 = 'XY100011610'
    3346offshore_name2 = 'XY100011611'
    … …  
    5366offshore_name22 = 'XYWADPI'
    5467
    55 # developed by NM&I
    56 coast_name = 'coastline'
     68#final topo name
     69combined_name ='perth_combined_elevation'
     70combined_smaller_name = 'perth_combined_elevation_smaller'
    5771
    58 # TIN model to fill in data gap
    59 bathy_interp = 'nearest_neighbour' #'interpolate'
    6072
    61 boundary_basename = 'SU-AU' # Mw ?
     73topographies_in_dir = home+sep+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep
     74topographies_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'topographies'+sep
     75topographies_time_dir = topographies_dir+build_time+sep
    6276
    63 #swollen/ all data output
    64 basename = 'source'
    65 codename = 'project.py'
     77#input topo file location
     78onshore_in_dir_name = topographies_in_dir + onshore_name
     79island_in_dir_name = topographies_in_dir + island_name
     80island_in_dir_name1 = topographies_in_dir + island_name1
     81island_in_dir_name2 = topographies_in_dir + island_name2
     82island_in_dir_name3 = topographies_in_dir + island_name3
    6683
    67 #Derive subdirectories and filenames
    68 local_time = strftime('%Y%m%d_%H%M%S',gmtime())
    69 meshdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep
    70 datadir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep
    71 gaugedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep
    72 polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
    73 boundarydir = home+sep+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
    74 outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'outputs'+sep
    75 outputtimedir = outputdir + local_time + sep
    76 polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
     84coast_in_dir_name = topographies_in_dir + coast_name
     85offshore_in_dir_name = topographies_in_dir + offshore_name
     86offshore1_in_dir_name = topographies_in_dir + offshore1_name
    7787
    78 gauge_filename = gaugedir + 'broome_gauges.csv'
    79 community_filename = gaugedir + 'CHINS_v2.csv'
    80 community_broome = gaugedir + 'community_broome.csv'
    81 codedir = getcwd()+sep                           
    82 codedirname = codedir + 'project.py'
    83 meshname = outputtimedir + 'mesh_' + basename
     88onshore_dir_name = topographies_dir + onshore_name
     89island_dir_name = topographies_dir + island_name
     90island_dir_name1 = topographies_dir + island_name1
     91island_dir_name2 = topographies_dir + island_name2
     92island_dir_name3 = topographies_dir + island_name3
    8493
    85 # Necessary if using point datasets, rather than grid
    86 onshore_dem_name = datadir + onshore_name
    87 offshore_interp_dem_name = datadir + bathy_interp
    88 offshore_dem_name1 = datadir + offshore_name1
    89 offshore_dem_name2 = datadir + offshore_name2
    90 offshore_dem_name3 = datadir + offshore_name3
    91 offshore_dem_name4 = datadir + offshore_name4
    92 offshore_dem_name5 = datadir + offshore_name5
    93 offshore_dem_name6 = datadir + offshore_name6
    94 offshore_dem_name7 = datadir + offshore_name7
    95 offshore_dem_name8 = datadir + offshore_name8
    96 offshore_dem_name9 = datadir + offshore_name9
    97 offshore_dem_name10 = datadir + offshore_name10
    98 offshore_dem_name11 = datadir + offshore_name11
    99 offshore_dem_name12 = datadir + offshore_name12
    100 offshore_dem_name13 = datadir + offshore_name13
    101 offshore_dem_name14 = datadir + offshore_name14
    102 offshore_dem_name15 = datadir + offshore_name15
    103 offshore_dem_name16 = datadir + offshore_name16
    104 offshore_dem_name17 = datadir + offshore_name17
    105 offshore_dem_name18 = datadir + offshore_name18
    106 offshore_dem_name19 = datadir + offshore_name19
    107 offshore_dem_name20 = datadir + offshore_name20
    108 offshore_dem_name21 = datadir + offshore_name21
    109 offshore_dem_name22 = datadir + offshore_name22
     94coast_dir_name = topographies_dir + coast_name
     95offshore_dir_name = topographies_dir + offshore_name
    11096
    111 coast_dem_name = datadir + coast_name
     97#final topo files
     98combined_dir_name = topographies_dir + combined_name
     99combined_time_dir_name = topographies_time_dir + combined_name
     100combined_smaller_name_dir = topographies_dir + combined_smaller_name
     101#combined_time_dir_final_name = topographies_time_dir + combined_final_name
    112102
    113 combined_dem_name   = datadir + 'broome_combined_elevation'
     103meshes_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'meshes'+sep
     104meshes_dir_name = meshes_dir + scenario_name
     105
     106polygons_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'polygons'+sep
     107tide_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'tide_data'+sep
     108
     109
     110boundaries_source = '????'
     111#boundaries locations
     112boundaries_in_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+boundaries_source+sep
     113boundaries_in_dir_name = boundaries_in_dir + scenario_name
     114boundaries_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep
     115boundaries_dir_name = boundaries_dir + scenario_name
     116#boundaries_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+build_time+sep
     117#boundaries_time_dir_name = boundaries_time_dir + boundaries_name  #Used by post processing
     118
     119#output locations
     120output_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep
     121output_build_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+build_time+sep
     122output_run_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+run_time+sep
     123output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
     124
     125#gauges
     126gauge_name = 'broome_gauges.csv'
     127gauges_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'gauges'+sep
     128gauges_dir_name = gauges_dir + gauge_name
     129
     130community_filename = gauges_dir + 'CHINS_v2.csv'
     131community_broome = gauges_dir + 'community_broome.csv'
     132
     133
    114134
    115135###############################
    116136# Domain definitions
    117137###############################
     138from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
    118139
    119 # bounding box for clipping MOST/URS output (much bigger than study area)
    120 ##south = degminsec2decimal_degrees(-19,0,0)
    121 ##north = degminsec2decimal_degrees(-17,15,0)
    122 ##west  = degminsec2decimal_degrees(120,0,0)
    123 ##east  = degminsec2decimal_degrees(122,30,0)
    124 ##
    125 ##d0 = [south, west]
    126 ##d1 = [south, east]
    127 ##d2 = [north, east]
    128 ##d3 = [north, west]
    129 ##poly_bc, zone = convert_points_from_latlon_to_utm([d0, d1, d2, d3])
    130 ##refzone = zone
     140poly_all = read_polygon(polygons_dir+'extent.csv')
     141res_poly_all = 100000
    131142
    132 # bounding polygon for study area
    133 polyAll = read_polygon(polygondir+'extent.csv')
     143###############################
     144# Interior region definitions
     145###############################
    134146
    135 # plot bounding polygon and make sure BC info surrounds it
    136 #plot_polygons([polyAll, poly_bc],'boundingpoly',verbose=False)
    137 print 'Area of bounding polygon', polygon_area(polyAll)/1000000.0
     147poly_1 = read_polygon(polygons_dir+'Broome_Local_Polygon_update.csv')
     148res_1 = 20000
     149
     150poly_2 = read_polygon(polygons_dir+'Broome_Close2_update.csv')
     151res_2 = 1000
     152
     153poly_3 = read_polygon(polygons_dir+'Broome_Coast_update2.csv')
     154res_3 = 1000
     155#assert zone == refzone
     156
     157interior_regions = [[poly_1,res_1],[poly_2,res_2]
     158                     ,[poly_3,res_3]]
     159
     160trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
     161
     162print 'min number triangles', trigs_min
    138163
    139164###################################################################
    … …  
    142167
    143168# exporting asc grid
    144 eastingmin = 412036.43
    145 eastingmax = 427184.37
    146 northingmin = 8007984
    147 northingmax = 8029830
    148 
    149        
    150        
     169eastingmin = 406215.87
     170eastingmax = 440208.78
     171northingmin = 7983427.73
     172northingmax = 8032834.52
    151173
    152174
    153 ###############################
    154 # Interior region definitions
    155 ###############################
    156175
    157 # broome digitized polygons
    158 poly_broome1 = read_polygon(polygondir+'Broome_Local_Polygon_update.csv')
    159 poly_broome2 = read_polygon(polygondir+'Broome_Close2_update.csv')
    160 poly_broome3 = read_polygon(polygondir+'Broome_Coast_update2.csv')
    161 #poly_broome4 = read_polygon(polygondir+'Cable_Beach_revised.csv')
    162 
    163 #plot_polygons([polyAll,poly_broome1,poly_broome2,poly_broome3],'boundingpoly2',verbose=False)
    164 print 'Area of local polygon', polygon_area(poly_broome1)/1000000.0
    165 print 'Area of close polygon', polygon_area(poly_broome2)/1000000.0
    166 print 'Area of coastal polygon', polygon_area(poly_broome3)/1000000.0
    167 #print 'Area of cable beach polygon', polygon_area(poly_broome4)/1000000.0
    168 
    169 for i in poly_broome3:
    170     v = is_inside_polygon(i,poly_broome1, verbose=False)
    171     if v == False: print v
    172 
    173 def number_mesh_triangles(interior_regions, bounding_poly, remainder_res):
    174     from anuga.utilities.polygon import polygon_area
    175    
    176     # TO DO check if any of the regions fall inside one another
    177     no_triangles = 0.0
    178     area = polygon_area(bounding_poly)
    179     for i,j in interior_regions:
    180         this_area = polygon_area(i)
    181         no_triangles += this_area/j
    182         area -= this_area
    183         print j, this_area/1000000., area/1000000.
    184     no_triangles += area/remainder_res
    185     return int(no_triangles/0.7)

  • anuga_work/production/broome_2006/run_broome.py

    r4063 r4153  
    1 """Script for running a tsunami inundation scenario for Broome, WA, Australia.
     1"""Script for running tsunami inundation scenario for Dampier, WA, Australia.
    22
    33Source data such as elevation and boundary data is assumed to be available in
    44directories specified by project.py
    5 The output sww file is stored in project.outputtimedir
     5The output sww file is stored in project.output_run_time_dir
    66
    77The scenario is defined by a triangular mesh created from project.polygon,
    8 the elevation data and a tsunami wave generated by MOST.
    9 
    10 Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006
     8the elevation data and a simulated submarine landslide.
     9
     10Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006
    1111"""
    1212
    13 #-------------------------------------------------------------------------------
     13#------------------------------------------------------------------------------
    1414# Import necessary modules
    15 #-------------------------------------------------------------------------------
     15#------------------------------------------------------------------------------
    1616
    1717# Standard modules
    18 import os
     18from os import sep, mkdir, access, F_OK
     19from os.path import dirname, basename
     20from shutil import copy
    1921import time
    20 from shutil import copy
    21 from os.path import dirname, basename
    22 from os import mkdir, access, F_OK, sep
    2322import sys
    2423
    2524# Related major packages
    26 from anuga.shallow_water import Domain, Reflective_boundary, \
    27                             Dirichlet_boundary, Time_boundary, File_boundary
    28 from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
    29 from anuga.abstract_2d_finite_volumes.combine_pts import combine_rectangular_points_files
     25from anuga.shallow_water import Domain
     26from anuga.shallow_water import Dirichlet_boundary, File_boundary, Reflective_boundary
     27from Numeric import allclose
     28from anuga.pmesh.mesh_interface import create_mesh_from_regions
    3029from anuga.geospatial_data.geospatial_data import *
    3130from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
     31from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
     32from anuga.utilities.polygon import plot_polygons, polygon_area
    3233
    3334# Application specific imports
    3435import project                 # Definition of file names and polygons
    3536
    36 #-------------------------------------------------------------------------------
     37#------------------------------------------------------------------------------
    3738# Copy scripts to time stamped output directory and capture screen
    3839# output to file
    39 #-------------------------------------------------------------------------------
    40 
    41 # creates copy of code in output dir
    42 copy_code_files(project.outputtimedir,__file__,dirname(project.__file__)+sep+ project.__name__+'.py' )
    43 myid = 0
    44 numprocs = 1
    45 #start_screen_catcher(project.outputtimedir, myid, numprocs)
    46 
    47 print 'USER:    ', project.user
    48 
    49 #-------------------------------------------------------------------------------
    50 # Preparation of topographic data
    51 #
    52 # Convert ASC 2 DEM 2 PTS using source data and store result in source data
    53 #-------------------------------------------------------------------------------
     40#------------------------------------------------------------------------------
     41
     42start_screen_catcher(project.output_run_time_dir, myid, numprocs)
    5443
    5544# filenames
    56 onshore_dem_name = project.onshore_dem_name
    57 offshore_interp_dem_name = project.offshore_interp_dem_name
    58 coast_points = project.coast_dem_name
    59 meshname = project.meshname+'.msh'
    60 
    61 # creates DEM from asc data
    62 convert_dem_from_ascii2netcdf(onshore_dem_name, use_cache=True, verbose=True)
    63 
    64 #creates pts file for onshore DEM
    65 dem2pts(onshore_dem_name, use_cache=True, verbose=True)
    66 
    67 # creates DEM from asc data
    68 convert_dem_from_ascii2netcdf(offshore_interp_dem_name, use_cache=True, verbose=True)
    69 
    70 #creates pts file for offshore interpolated DEM
    71 dem2pts(offshore_interp_dem_name, use_cache=True, verbose=True)
    72 
    73 print 'create offshore'
    74 G1 = Geospatial_data(file_name = project.offshore_dem_name1 + '.xya')+\
    75      Geospatial_data(file_name = project.offshore_dem_name2 + '.xya')+\
    76      Geospatial_data(file_name = project.offshore_dem_name3 + '.xya')+\
    77      Geospatial_data(file_name = project.offshore_dem_name4 + '.xya')+\
    78      Geospatial_data(file_name = project.offshore_dem_name5 + '.xya')+\
    79      Geospatial_data(file_name = project.offshore_dem_name6 + '.xya')+\
    80      Geospatial_data(file_name = project.offshore_dem_name7 + '.xya')+\
    81      Geospatial_data(file_name = project.offshore_dem_name8 + '.xya')+\
    82      Geospatial_data(file_name = project.offshore_dem_name9 + '.xya')+\
    83      Geospatial_data(file_name = project.offshore_dem_name10 + '.xya')+\
    84      Geospatial_data(file_name = project.offshore_dem_name11 + '.xya')+\
    85      Geospatial_data(file_name = project.offshore_dem_name12 + '.xya')+\
    86      Geospatial_data(file_name = project.offshore_dem_name13 + '.xya')+\
    87      Geospatial_data(file_name = project.offshore_dem_name14 + '.xya')+\
    88      Geospatial_data(file_name = project.offshore_dem_name15 + '.xya')+\
    89      Geospatial_data(file_name = project.offshore_dem_name16 + '.xya')+\
    90      Geospatial_data(file_name = project.offshore_dem_name17 + '.xya')+\
    91      Geospatial_data(file_name = project.offshore_dem_name18 + '.xya')+\
    92      Geospatial_data(file_name = project.offshore_dem_name19 + '.xya')+\
    93      Geospatial_data(file_name = project.offshore_dem_name20 + '.xya')+\
    94      Geospatial_data(file_name = project.offshore_dem_name21 + '.xya')+\
    95      Geospatial_data(file_name = project.offshore_dem_name22 + '.xya')+\
    96      Geospatial_data(file_name = project.offshore_interp_dem_name + '.pts')
    97 print 'create onshore'
    98 G2 = Geospatial_data(file_name = project.onshore_dem_name + '.pts')
    99 print 'create coast'
    100 G3 = Geospatial_data(file_name = project.coast_dem_name + '.xya')
    101 print 'add'
    102 G = G1 + G2 + G3
    103 print 'export points'
    104 G.export_points_file(project.combined_dem_name + '.pts')
    105 G.export_points_file(project.combined_dem_name + '.xya')
    106 
    107 #----------------------------------------------------------------------------
    108 # Create the triangular mesh based on overall clipping polygon with a tagged
     45boundaries_name = project.scenario
     46meshes_dir_name = project.meshes_dir_name+'.msh'
     47boundaries_dir_name = project.boundaries_dir_name
     48
     49tide = project.tide
     50
     51# creates copy of code in output dir
     52if myid == 0:
     53    copy_code_files(project.output_run_time_dir,__file__,
     54                 dirname(project.__file__)+sep+ project.__name__+'.py' )
     55barrier()
     56
     57print 'USER: ', project.user
     58print 'min triangles', project.trigs_min,
     59print 'Note: This is generally about 20% less than the final amount'
     60
     61#--------------------------------------------------------------------------
     62# Create the triangular mesh based on overall clipping polygon with a
     63# tagged
    10964# boundary and interior regions defined in project.py along with
    11065# resolutions (maximal area of per triangle) for each polygon
    111 #-------------------------------------------------------------------------------
    112 
    113 from anuga.pmesh.mesh_interface import create_mesh_from_regions
    114 remainder_res = 500000
    115 local_res = 25000
    116 broome_res = 5000
    117 coast_res = 500
    118 interior_regions = [[project.poly_broome1, local_res],
    119                     [project.poly_broome2, broome_res],
    120                     [project.poly_broome3, coast_res]]
    121 
    122 from project import number_mesh_triangles
    123 print 'estimate of number of triangles', \
    124       number_mesh_triangles(interior_regions, project.polyAll, remainder_res)
    125 
     66#--------------------------------------------------------------------------
     67
     68if myid == 0:     
     69
     70    print 'start create mesh from regions'
     71    create_mesh_from_regions(project.poly_all,
     72                         boundary_tags={'back': [1, 2], 'side': [0,3],
     73                                        'ocean': [4, 5, 6]},
     74                         maximum_triangle_area=project.res_poly_all,
     75                         interior_regions=project.interior_regions,
     76                         filename=meshes_dir_name,
     77                         use_cache=True,
     78                         verbose=True)
     79
     80# to sync all processors are ready
     81barrier()
     82
     83#-------------------------------------------------------------------------
     84# Setup computational domain
     85#-------------------------------------------------------------------------
     86print 'Setup computational domain'
     87#domain = Domain(meshes_time_dir_name, use_cache=True, verbose=True)
     88domain = Domain(meshes_dir_name, use_cache=True, verbose=True)
     89print domain.statistics()
     90'''
     91print 'starting to create boundary conditions'
     92boundaries_in_dir_name = project.boundaries_in_dir_name
     93
     94from anuga.shallow_water.data_manager import urs2sww, ferret2sww
     95
     96# put above distribute
     97print 'boundary file is: ',boundaries_dir_name
    12698from caching import cache
    127 _ = cache(create_mesh_from_regions,
    128           project.polyAll,
    129            {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2],
    130                               'e3': [3], 'e4':[4], 'e5': [5],
    131                               'e6': [6]},
    132            'maximum_triangle_area': remainder_res,
    133            'filename': meshname,
    134            'interior_regions': interior_regions},
    135           verbose = True, evaluate=False)
    136 print 'created mesh'
    137 
    138 #-------------------------------------------------------------------------------                                 
    139 # Setup computational domain
    140 #-------------------------------------------------------------------------------                                 
    141 domain = Domain(meshname, use_cache = True, verbose = True)
    142 
    143 print 'Number of triangles = ', len(domain)
    144 print 'The extent is ', domain.get_extent()
    145 print domain.statistics()
    146 
    147 domain.set_name(project.basename)
    148 domain.set_datadir(project.outputtimedir)
    149 domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
    150 domain.set_minimum_storable_height(0.01)
    151 
    152 #-------------------------------------------------------------------------------                                 
     99if myid == 0:
     100    cache(ferret2sww,
     101          (boundaries_in_dir_name,
     102    #       boundaries_time_dir_name),
     103           boundaries_dir_name),
     104          {'verbose': True,
     105           'minlat': project.south_boundary,
     106           'maxlat': project.north_boundary,
     107           'minlon': project.west_boundary,
     108           'maxlon': project.east_boundary,
     109#           'minlat': project.south,
     110#           'maxlat': project.north,
     111#           'minlon': project.west,
     112#           'maxlon': project.east,
     113           'mint': 0, 'maxt': 35100,
     114           'origin': domain.geo_reference.get_origin(),
     115           'mean_stage': project.tide,
     116#           'zscale': 1,                 #Enhance tsunami
     117           'fail_on_NaN': False},
     118           verbose = True,
     119           )
     120barrier()
     121'''
     122
     123#-------------------------------------------------------------------------
    153124# Setup initial conditions
    154 #-------------------------------------------------------------------------------
    155 
    156 tide = 0.0
    157 domain.set_quantity('stage', tide)
    158 domain.set_quantity('friction', 0.0)
    159 domain.set_quantity('elevation',
    160                     filename = project.combined_dem_name + '.pts',
     125#-------------------------------------------------------------------------
     126if myid == 0:
     127
     128    print 'Setup initial conditions'
     129
     130    domain.set_quantity('stage', tide)
     131    domain.set_quantity('friction', 0.01)
     132    #combined_time_dir_name = project.topographies_dir+build_time+project.combined_name
     133    print 'Start Set quantity'
     134
     135    domain.set_quantity('elevation',
     136                    filename = project.combined_dir_name+'.txt',
     137#                    filename = project.combined_smaller_name_dir+'.xya',
    161138                    use_cache = True,
    162139                    verbose = True,
    163                     alpha = 0.1
    164                     )
    165 
    166 #-------------------------------------------------------------------------------                                 
    167 # Setup boundary conditions
    168 #-------------------------------------------------------------------------------
    169 '''
    170 print 'start urs2sww'
    171 print '', project.boundary_basename
    172 from anuga.shallow_water.data_manager import urs2sww
    173 
    174 south = project.south
    175 north = project.north
    176 west  = project.west
    177 east  = project.east
    178 
    179 #note only need to do when an SWW file for the MOST boundary doesn't exist
    180 cache(urs2sww,
    181       (source_dir + project.boundary_basename,
    182        source_dir + project.boundary_basename),
    183       {'verbose': True,
    184        'minlat': south,
    185        'maxlat': north,
    186        'minlon': west,
    187        'maxlon': east,
    188        #'origin': domain.geo_reference.get_origin(),
    189        'mean_stage': tide,
    190        'zscale': 1,                 #Enhance tsunami
    191        'fail_on_NaN': False,
    192        'inverted_bathymetry': True},
    193       #evaluate = True,
    194        verbose = True,
    195        dependencies = source_dir + project.boundary_basename + '.sww')
    196 
    197 '''
     140                    alpha = 0.1)
     141    print 'Finished Set quantity'
     142barrier()
     143
     144#------------------------------------------------------
     145# Distribute domain to implement parallelism !!!
     146#------------------------------------------------------
     147
     148if numprocs > 1:
     149    domain=distribute(domain)
     150
     151#------------------------------------------------------
     152# Set domain parameters
     153#------------------------------------------------------
     154
     155domain.set_name(project.scenario_name)
     156domain.set_datadir(project.output_run_time_dir)
     157domain.set_default_order(2) # Apply second order scheme
     158domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
     159domain.set_store_vertices_uniquely(False)
     160domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
     161domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
     162
     163#-------------------------------------------------------------------------
     164# Setup boundary conditions
     165#-------------------------------------------------------------------------
    198166print 'Available boundary tags', domain.get_boundary_tags()
    199167
    200 #Bf = File_boundary(source_dir + project.boundary_basename + '.sww',
    201 #                    domain, verbose = True)
     168print 'Reading Boundary file'
     169print 'domain id', id(domain)
     170#Bf = File_boundary(boundaries_dir_name + '.sww',
     171#                  domain, time_thinning=5, use_cache=True, verbose=True)
     172
     173print 'finished reading boundary file'
    202174Br = Reflective_boundary(domain)
    203175Bd = Dirichlet_boundary([tide,0,0])
    204 
    205 # 7 min square wave starting at 1 min, 6m high
    206 Bw = Time_boundary(domain = domain,
    207                    f=lambda t: [(60<t<480)*10, 0, 0])
    208 
    209 domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd,
    210                       'e5': Bd,  'e6': Bd} )
    211 
    212 
    213 #-------------------------------------------------------------------------------                                 
     176Bo = Dirichlet_boundary([tide+5.0,0,0])
     177
     178print'set_boundary'
     179domain.set_boundary({'back': Br,
     180                     'side': Bd,
     181                     'ocean': Bd})
     182print'finish set boundary'
     183
     184#----------------------------------------------------------------------------
    214185# Evolve system through time
    215 #-------------------------------------------------------------------------------
    216 import time
     186#----------------------------------------------------------------------------
     187
    217188t0 = time.time()
    218189
    219 for t in domain.evolve(yieldstep = 240, finaltime = 480):
     190for t in domain.evolve(yieldstep = 60, finaltime = 10000):
    220191    domain.write_time()
    221     domain.write_boundary_statistics(tags = 'e14')     
     192    domain.write_boundary_statistics(tags = 'ocean')
     193    if allclose(t, 120):
     194        domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bo})
     195
     196    if allclose(t, 1020):
     197        domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bd})
     198
    222199   
    223200print 'That took %.2f seconds' %(time.time()-t0)
    224201
    225 print 'finished'
Note: See TracChangeset for help on using the changeset viewer.