Changeset 6259


Ignore:
Timestamp:
Feb 3, 2009, 9:06:40 AM (15 years ago)
Author:
ole
Message:

Work on simplifying model script

File:
1 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/production/busselton/standardised_version/run_busselton.py

    r6256 r6259  
    3434
    3535# Related major packages
    36 from anuga.shallow_water import Domain
    37 from anuga.shallow_water import Dirichlet_boundary
    38 from anuga.shallow_water import File_boundary
    39 from anuga.shallow_water import Reflective_boundary
    40 from anuga.shallow_water import Field_boundary
    41 from Numeric import allclose
    42 from anuga.shallow_water.data_manager import export_grid, create_sts_boundary
    43 from anuga.pmesh.mesh_interface import create_mesh_from_regions
    44 from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
    45 from anuga_parallel.parallel_abstraction import get_processor_name
    46 from anuga.caching import myhash
    47 #from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
    48 from anuga.fit_interpolate.benchmark_least_squares import mem_usage
    49 from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
    50 from anuga.geospatial_data.geospatial_data import find_optimal_smoothing_parameter
    51 from polygon import Polygon_function
     36from anuga.interface import create_domain_from_regions
     37from anuga.interface import Dirichlet_boundary
     38from anuga.interface import Reflective_boundary
     39from anuga.interface import Field_boundary
     40from anuga.interface import create_sts_boundary
     41from anuga.interface import csv2building_polygons
     42
     43from anuga.shallow_water.data_manager import start_screen_catcher
     44from anuga.shallow_water.data_manager import copy_code_files
     45from anuga.utilities.polygon import read_polygon, Polygon_function
    5246   
    5347# Application specific imports
    5448import project  # Definition of file names and polygons
    55 numprocs = 1
    56 myid = 0
    5749
    58 def run_model(**kwargs):
    5950   
    60     #-----------------------------------------------------------------------
    61     # Copy scripts to time stamped output directory and capture screen
    62     # output to file
    63     #-----------------------------------------------------------------------
    64     print "Processor Name:",get_processor_name()
     51#-----------------------------------------------------------------------
     52# Copy scripts to time stamped output directory and capture screen
     53# output to file
     54#-----------------------------------------------------------------------
    6555
    66     #copy script must be before screen_catcher
    67 
    68     print 'output_dir',kwargs['output_dir']
    69    
    70     copy_code_files(kwargs['output_dir'],__file__,
    71              dirname(project.__file__)+sep+ project.__name__+'.py' )
    72 
    73     store_parameters(**kwargs)
    74 
    75     start_screen_catcher(kwargs['output_dir'], myid, numprocs)
    76 
    77     print "Processor Name:",get_processor_name()
    78    
    79     #-----------------------------------------------------------------------
    80     # Domain definitions
    81     #-----------------------------------------------------------------------
    82 
    83     # Read in boundary from ordered sts file
    84     urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir_event, project.scenario_name))
    85 
    86     # Reading the landward defined points, this incorporates the original clipping
    87     # polygon minus the 100m contour
    88     landward_bounding_polygon = read_polygon(project.landward_dir)
    89 
    90     # Combine sts polyline with landward points
    91     bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
    92    
    93     # counting segments
    94     N = len(urs_bounding_polygon)-1
    95 
    96     # boundary tags refer to project.landward 4 points equals 5 segments start at N
    97     boundary_tags={'back': [N+1,N+2,N+3,N+4, N+5], 'side': [N,N+6], 'ocean': range(N)}
    98 
    99     #--------------------------------------------------------------------------
    100     # Create the triangular mesh based on overall clipping polygon with a tagged
    101     # boundary and interior regions defined in project.py along with
    102     # resolutions (maximal area of per triangle) for each polygon
    103     #--------------------------------------------------------------------------
    104 
    105     # IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
    106     # causes problems with the ability to cache set quantity which takes alot of times
    107        
    108 
    109     # FIXME(Ole): Introduce create_domain_from_regions: Simpler and caches well
    110     print 'Start create mesh from regions'
    111 
    112     create_mesh_from_regions(bounding_polygon,
    113                          boundary_tags=boundary_tags,
    114                          maximum_triangle_area=project.res_poly_all,
    115                          interior_regions=project.interior_regions,
    116                          filename=project.meshes_dir_name,
    117                          use_cache=True,
    118                          verbose=True)
    119    
    120     #-------------------------------------------------------------------------
    121     # Setup computational domain
    122     #-------------------------------------------------------------------------
    123     print 'Setup computational domain'
    124 
    125     domain = Domain(project.meshes_dir_name, use_cache=False, verbose=True)
    126     print domain.statistics()
     56#copy script must be before screen_catcher
     57#copy_code_files(project.output_run_time_dir, __file__,
     58#         dirname(project.__file__)+sep+ project.__name__+'.py' )
     59#start_screen_catcher(project.output_run_time_dir, myid, numprocs)
    12760
    12861
    129     #-------------------------------------------------------------------------
    130     # Setup initial conditions
    131     #-------------------------------------------------------------------------
    132     print 'Setup initial conditions'
     62#--------------------------------------------------------------------------
     63# Create the computational domain based on overall clipping polygon with
     64# a tagged boundary and interior regions defined in project.py along with
     65# resolutions (maximal area of per triangle) for each polygon
     66#--------------------------------------------------------------------------
     67print 'Create computational domain'
    13368
    134     # sets the initial stage in the offcoast region only
    135     IC = Polygon_function( [(project.poly_mainland, 0),
    136                             (project.poly_marina, 0)],
    137                            default = project.tide,
    138                            geo_reference = domain.geo_reference)
    139     domain.set_quantity('stage', IC)
     69# Read in boundary from ordered sts file
     70urs_bounding_polygon = create_sts_boundary(project.urs_boundary_name)
    14071
    141     domain.set_quantity('friction', project.friction)
    142    
    143     print 'Start set quantity', kwargs['elevation_file']
     72# Reading the landward defined points, this incorporates the original clipping
     73# polygon minus the 100m contour
     74landward_bounding_polygon = read_polygon(project.landward_dir)
    14475
    145     domain.set_quantity('elevation',
    146                         filename = kwargs['elevation_file'],
    147                         use_cache = False,
    148                         verbose = True,
    149                         alpha = project.alpha)
    150     print 'Finished Set quantity'
     76# Combine sts polyline with landward points
     77bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
    15178
    152 ##   #------------------------------------------------------
    153 ##    # Distribute domain to implement parallelism !!!
    154 ##    #------------------------------------------------------
    155 ##
    156 ##    if numprocs > 1:
    157 ##        domain=distribute(domain)
     79# Number of boundary segments
     80N = len(urs_bounding_polygon)-1
    15881
    159     #------------------------------------------------------
    160     # Set domain parameters
    161     #------------------------------------------------------
    162     domain.set_name(project.scenario_name)
    163     domain.set_datadir(kwargs['output_dir'])    # FIXME: Use project.output_run_time_dir
    164     domain.set_default_order(2)                 # Apply second order scheme
    165     domain.set_minimum_storable_height(0.01)    # Don't store anything less than 1cm
    166     #domain.set_store_vertices_uniquely(False)
    167     #domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
    168     #domain.tight_slope_limiters = 1
     82# Boundary tags refer to project.landward 4 points equals 5 segments start at N
     83boundary_tags={'back': [N+1,N+2,N+3,N+4, N+5],
     84               'side': [N,N+6],
     85               'ocean': range(N)}
    16986
    170     #-------------------------------------------------------------------------
    171     # Setup boundary conditions
    172     #-------------------------------------------------------------------------
    173     print 'Available boundary tags', domain.get_boundary_tags()
    174    
    175     boundary_urs_out=project.boundaries_dir_event + sep + project.scenario_name
    176 
    177     Br = Reflective_boundary(domain)
    178     Bd = Dirichlet_boundary(project.tide,0,0])
    179    
    180     print 'Available boundary tags', domain.get_boundary_tags()
    181     Bf = Field_boundary(boundary_urs_out+'.sts',
    182                         domain, mean_stage= project.tide,
    183                         time_thinning=1,
    184                         default_boundary=Bd,
    185                         use_cache=True,
    186                         verbose = True,
    187                         boundary_polygon=bounding_polygon)
    188 
    189     domain.set_boundary({'back': Br,
    190                          'side': Bd,
    191                          'ocean': Bf})
     87# Build mesh and domain
     88domain = create_domain_from_regions(bounding_polygon,
     89                                    boundary_tags=boundary_tags,
     90                                    maximum_triangle_area=project.res_poly_all,
     91                                    interior_regions=project.interior_regions,
     92                                    mesh_filename=project.meshes_dir_name,
     93                                    use_cache=True,
     94                                    verbose=True)
     95print domain.statistics()
    19296
    19397
    194     print 'Finish set boundary'
     98domain.set_name(project.scenario_name)
     99domain.set_datadir(project.output_run_time_dir)
     100domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm
    195101
    196     #----------------------------------------------------------------------------
    197     # Evolve system through time
    198     #--------------------------------------------------------------------
    199     t0 = time.time()
     102#-------------------------------------------------------------------------
     103# Setup initial conditions
     104#-------------------------------------------------------------------------
     105print 'Setup initial conditions'
    200106
    201     for t in domain.evolve(yieldstep=project.yieldstep,
    202                            finaltime=project.finaltime,
    203                            skip_initial_step=False):
    204         domain.write_time()
    205         domain.write_boundary_statistics(tags = 'ocean')
     107# Set the initial stage in the offcoast region only
     108IC = Polygon_function([(project.poly_mainland, 0),
     109                       (project.poly_marina, 0)],
     110                      default=project.tide,
     111                      geo_reference=domain.geo_reference)
     112domain.set_quantity('stage', IC, use_cache=True, verbose=True)
    206113
    207     # these outputs should be checked with the resultant inundation map
    208     x, y = domain.get_maximum_inundation_location()
    209     q = domain.get_maximum_inundation_elevation()
    210     print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
     114domain.set_quantity('friction', project.friction)
    211115
    212     print 'Simulation took %.2f seconds' %(time.time()-t0)
     116domain.set_quantity('elevation',
     117                    filename=project.combined_dir_name+'.pts',
     118                    use_cache=True,
     119                    verbose=True,
     120                    alpha = project.alpha)
    213121
    214122
    215      
     123#-------------------------------------------------------------------------
     124# Setup boundary conditions
     125#-------------------------------------------------------------------------
    216126
     127print 'Set boundary - available tags:', domain.get_boundary_tags()
     128
     129boundary_urs_out=project.boundaries_dir_event + sep + project.scenario_name
     130
     131Br = Reflective_boundary(domain)
     132Bd = Dirichlet_boundary([project.tide,0,0])
     133
     134
     135Bf = Field_boundary(boundary_urs_out+'.sts',
     136                    domain, mean_stage=project.tide,
     137                    time_thinning=1,
     138                    default_boundary=Bd,
     139                    boundary_polygon=bounding_polygon,                   
     140                    use_cache=True,
     141                    verbose=True)
     142
     143
     144domain.set_boundary({'back': Br,
     145                     'side': Bd,
     146                     'ocean': Bf})
     147
     148#-------------------------------------------------------------------------
     149# Evolve system through time
     150#-------------------------------------------------------------------------
     151t0 = time.time()
     152for t in domain.evolve(yieldstep=project.yieldstep,
     153                       finaltime=project.finaltime,
     154                       skip_initial_step=False):
     155    print domain.timestepping_statistics()
     156    print domain.boundary_statistics(tags='ocean')
     157
     158print 'Simulation took %.2f seconds' %(time.time()-t0)
    217159
    218160   
    219161   
    220 #-------------------------------------------------------------
    221 if __name__ == "__main__":
    222    
    223     kwargs={}
    224     kwargs['file_name']=project.dir_comment
    225     kwargs['output_dir']=project.output_run_time_dir
    226     kwargs['elevation_file']=project.combined_dir_name+'.pts'
    227      
    228     run_model(**kwargs)
     162
    229163     
    230164   
Note: See TracChangeset for help on using the changeset viewer.