Changeset 6178


Ignore:
Timestamp:
Jan 15, 2009, 5:45:46 PM (10 years ago)
Author:
ole
Message:

Clean up of Patong script and added new function create_domain_from_regions

Files:
6 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

    r6175 r6178  
    348348   
    349349    if time_limit is not None:
     350        #if verbose is True:
     351        #    print '****** Time limit', time_limit
     352        #    print '****** Start time', starttime
     353        #    print '****** Time in ', time[0], time[-1]
     354       
    350355        # Adjust given time limit to given start time
    351356        time_limit = time_limit - starttime
     357
    352358
    353359        # Find limit point
  • anuga_core/source/anuga/pmesh/mesh_interface.py

    r6177 r6178  
    9090              'minimum_triangle_angle': minimum_triangle_angle,
    9191              'fail_if_polygons_outside': fail_if_polygons_outside,
    92               'verbose': verbose}   # FIXME (Ole): Should be bypassed one day
    93                                     # What should be bypassed? Verbose?
    94    
    95     #print 'kwargs', kwargs
     92              'verbose': verbose}   # FIXME (Ole): Should be bypassed one day. See ticket:14
    9693
    9794    # Call underlying engine with or without caching
  • anuga_core/source/anuga/shallow_water/__init__.py

    r6045 r6178  
    88# Make selected classes available directly
    99from shallow_water_domain import Domain,\
    10     Transmissive_boundary, Reflective_boundary,\
    11     Dirichlet_boundary, Time_boundary, File_boundary,\
    12     Transmissive_momentum_set_stage_boundary,\
    13     Dirichlet_discharge_boundary,\
    14     Field_boundary
     10     create_domain_from_regions,\
     11     Transmissive_boundary, Reflective_boundary,\
     12     Dirichlet_boundary, Time_boundary, File_boundary,\
     13     Transmissive_momentum_set_stage_boundary,\
     14     Dirichlet_discharge_boundary,\
     15     Field_boundary
    1516
    1617
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r6173 r6178  
    9494from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
    9595     import Transmissive_boundary
    96 
     96from anuga.pmesh.mesh_interface import create_mesh_from_regions
    9797from anuga.utilities.numerical_tools import gradient, mean, ensure_numeric
    9898from anuga.geospatial_data.geospatial_data import ensure_geospatial
     
    117117
    118118
    119 #---------------------
     119
     120#---------------------------
     121# Create domain from regions
     122#---------------------------
     123
     124def create_domain_from_regions(bounding_polygon,
     125                               boundary_tags,
     126                               maximum_triangle_area=None,
     127                               mesh_filename=None,
     128                               interior_regions=None,
     129                               interior_holes=None,
     130                               poly_geo_reference=None,
     131                               mesh_geo_reference=None,
     132                               minimum_triangle_angle=28.0,
     133                               fail_if_polygons_outside=True,
     134                               use_cache=False,
     135                               verbose=True):
     136    """Create domain from bounding polygons and resolutions.
     137
     138    bounding_polygon is a list of points in Eastings and Northings,
     139    relative to the poly_geo_reference.
     140
     141    Boundary tags is a dictionary of symbolic tags. For every tag there
     142    is a list of indices referring to segments associated with that tag.
     143    If a segment is omitted it will be assigned the default tag ''.
     144
     145    maximum_triangle_area is the maximal area per triangle
     146    for the bounding polygon, excluding the  interior regions.
     147
     148    Interior_regions is a list of tuples consisting of (polygon,
     149    resolution) for each region to be separately refined. Do not have
     150    polygon lines cross or be on-top of each other.  Also do not have
     151    polygon close to each other.
     152   
     153    NOTE: If a interior_region is outside the bounding_polygon it should
     154    throw an error
     155   
     156    Interior_holes is a list of ploygons for each hole.
     157
     158    This function does not allow segments to share points - use underlying
     159    pmesh functionality for that
     160
     161    poly_geo_reference is the geo_reference of the bounding polygon and
     162    the interior polygons.
     163    If none, assume absolute.  Please pass one though, since absolute
     164    references have a zone.
     165   
     166    mesh_geo_reference is the geo_reference of the mesh to be created.
     167    If none is given one will be automatically generated.  It was use
     168    the lower left hand corner of  bounding_polygon (absolute)
     169    as the x and y values for the geo_ref.
     170   
     171    Returns the shallow water domain instance
     172
     173    Note, interior regions should be fully nested, as overlaps may cause
     174    unintended resolutions.
     175
     176    fail_if_polygons_outside: If True (the default) Exception in thrown
     177    where interior polygons fall outside bounding polygon. If False, these
     178    will be ignored and execution continued.
     179       
     180   
     181    """
     182
     183
     184    # Build arguments and keyword arguments for use with caching or apply.
     185    args = (bounding_polygon,
     186            boundary_tags)
     187   
     188    kwargs = {'maximum_triangle_area': maximum_triangle_area,
     189              'mesh_filename': mesh_filename,
     190              'interior_regions': interior_regions,
     191              'interior_holes': interior_holes,
     192              'poly_geo_reference': poly_geo_reference,
     193              'mesh_geo_reference': mesh_geo_reference,
     194              'minimum_triangle_angle': minimum_triangle_angle,
     195              'fail_if_polygons_outside': fail_if_polygons_outside,
     196              'verbose': verbose} #FIXME (Ole): See ticket:14
     197
     198    # Call underlying engine with or without caching
     199    if use_cache is True:
     200        try:
     201            from anuga.caching import cache
     202        except:
     203            msg = 'Caching was requested, but caching module'+\
     204                  'could not be imported'
     205            raise msg
     206
     207
     208        domain = cache(_create_domain_from_regions,
     209                       args, kwargs,
     210                       verbose=verbose,
     211                       compression=False)
     212    else:
     213        domain = apply(_create_domain_from_regions,
     214                       args, kwargs)
     215
     216    return domain
     217
     218       
     219def _create_domain_from_regions(bounding_polygon,
     220                                boundary_tags,
     221                                maximum_triangle_area=None,
     222                                mesh_filename=None,                           
     223                                interior_regions=None,
     224                                interior_holes=None,
     225                                poly_geo_reference=None,
     226                                mesh_geo_reference=None,
     227                                minimum_triangle_angle=28.0,
     228                                fail_if_polygons_outside=True,
     229                                verbose=True):
     230    """_create_domain_from_regions - internal function.
     231
     232    See create_domain_from_regions for documentation.
     233    """
     234
     235    create_mesh_from_regions(bounding_polygon,
     236                             boundary_tags,
     237                             maximum_triangle_area=maximum_triangle_area,
     238                             interior_regions=interior_regions,
     239                             filename=mesh_filename,
     240                             interior_holes=interior_holes,
     241                             poly_geo_reference=poly_geo_reference,
     242                             mesh_geo_reference=mesh_geo_reference,
     243                             minimum_triangle_angle=minimum_triangle_angle,
     244                             fail_if_polygons_outside=fail_if_polygons_outside,
     245                             use_cache=False,
     246                             verbose=verbose)
     247
     248    domain = Domain(mesh_filename, use_cache=False, verbose=verbose)
     249
     250
     251    return domain
     252   
     253
     254
     255
     256
     257
     258#
    120259# Shallow water domain
    121260#---------------------
  • anuga_work/production/patong/project.py

    r6168 r6178  
    4848finaltime=15000         # final time for simulation 15000
    4949
    50 setup='final'  # Final can be replaced with trial or basic.
     50setup='trial'  # Final can be replaced with trial or basic.
    5151               # Either will result in a coarser mesh that will allow a
    5252               # faster, but less accurate, simulation.
  • anuga_work/production/patong/run_patong.py

    r6175 r6178  
    2323from os import sep
    2424import os
    25 from os.path import dirname, basename
    26 from os import mkdir, access, F_OK
    27 from shutil import copy
     25from os.path import dirname
    2826import time
    29 import sys
    3027
    3128# Related major packages
     
    3532from anuga.shallow_water import Reflective_boundary
    3633from anuga.shallow_water import Field_boundary
    37 from Numeric import allclose
    3834from anuga.shallow_water.data_manager import export_grid, create_sts_boundary, csv2building_polygons
    39 from anuga.pmesh.mesh_interface import create_mesh_from_regions
     35from anuga.shallow_water import create_domain_from_regions
    4036from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
    4137from anuga.caching import myhash
     
    5147myid = 0
    5248
    53 def run_model(**kwargs):
    54    
    55     #------------------------------------------------------------------------------
    56     # Copy scripts to time stamped output directory and capture screen
    57     # output to file
    58     #------------------------------------------------------------------------------
    59 
    60     #copy script must be before screen_catcher
    61 
    62     print 'output_dir',kwargs['output_dir']
    63    
    64     copy_code_files(kwargs['output_dir'],__file__,
    65              dirname(project.__file__)+sep+ project.__name__+'.py' )
    66 
    67     store_parameters(**kwargs)
    68 
    69     start_screen_catcher(kwargs['output_dir'], myid, numprocs)
    7049
    7150   
    72     #-----------------------------------------------------------------------
    73     # Domain definitions
    74     #-----------------------------------------------------------------------
     51#------------------------------------------------------------------------------
     52# Copy scripts to time stamped output directory and capture screen
     53# output to file
     54#------------------------------------------------------------------------------
    7555
    76     # Read in boundary from ordered sts file
    77     urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir,project.scenario_name))
     56#copy script must be before screen_catcher
    7857
    79     # Reading the landward defined points, this incorporates the original clipping
    80     # polygon minus the 100m contour
    81     landward_bounding_polygon = read_polygon(project.landward_dir)
     58output_dir = project.output_run_time_dir
     59print 'output_dir', output_dir
    8260
    83     # Combine sts polyline with landward points
    84     bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
    85    
    86     # counting segments
    87     N = len(urs_bounding_polygon)-1
     61#copy_code_files(output_dir,__file__,
     62#         dirname(project.__file__)+sep+ project.__name__+'.py' )
    8863
    89     # boundary tags refer to project.landward 4 points equals 5 segments start at N
    90     boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4], 'ocean': range(N)}
    91 
    92     #--------------------------------------------------------------------------
    93     # Create the triangular mesh based on overall clipping polygon with a tagged
    94     # boundary and interior regions defined in project.py along with
    95     # resolutions (maximal area of per triangle) for each polygon
    96     #--------------------------------------------------------------------------
    97 
    98     # IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
    99     # causes problems with the ability to cache set quantity which takes alot of times
    100        
    101     print 'start create mesh from regions'
    102 
    103     create_mesh_from_regions(bounding_polygon,
    104                              boundary_tags=boundary_tags,
    105                              maximum_triangle_area=project.res_poly_all,
    106                              interior_regions=project.interior_regions,
    107                              filename=project.meshes_dir_name,
    108                              use_cache=True,
    109                              verbose=True)
    110    
    111     #-------------------------------------------------------------------------
    112     # Setup computational domain
    113     #-------------------------------------------------------------------------
    114     print 'Setup computational domain'
    115 
    116     domain = Domain(project.meshes_dir_name, use_cache=False, verbose=True)
    117     print 'memory usage before del domain',mem_usage()
    118        
    119     print domain.statistics()
    120     print 'triangles',len(domain)
    121    
    122     kwargs['act_num_trigs']=len(domain)
     64#start_screen_catcher(output_dir, myid, numprocs)
    12365
    12466
    125     #-------------------------------------------------------------------------
    126     # Setup initial conditions
    127     #-------------------------------------------------------------------------
    128     print 'Setup initial conditions'
     67#-----------------------------------------------------------------------
     68# Domain definitions
     69#-----------------------------------------------------------------------
    12970
    130     # sets the initial stage in the offcoast region only
    131     IC = Polygon_function( [(project.poly_mainland, 0)], default = kwargs['tide'],
    132                              geo_reference = domain.geo_reference)
    133     domain.set_quantity('stage', IC)
    134     #domain.set_quantity('stage',kwargs['tide'] )
    135     domain.set_quantity('friction', kwargs['friction'])
    136    
    137     print 'Start Set quantity',kwargs['elevation_file']
     71# Read in boundary from ordered sts file
     72urs_boundary_name = os.path.join(project.boundaries_dir,
     73                                 project.scenario_name)
     74urs_bounding_polygon = create_sts_boundary(urs_boundary_name)
    13875
    139     domain.set_quantity('elevation',
    140                         filename = kwargs['elevation_file'],
    141                         use_cache = False,
    142                         verbose = True,
    143                         alpha = kwargs['alpha'])
     76# Reading the landward defined points, this incorporates the original clipping
     77# polygon minus the 100m contour
     78landward_bounding_polygon = read_polygon(project.landward_dir)
    14479
    145     # Add buildings from file
    146     building_polygons, building_heights = csv2building_polygons(project.building_polygon_file)
     80# Combine sts polyline with landward points
     81bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
    14782
    148     L = []
    149     for key in building_polygons:
    150         poly = building_polygons[key]
    151         elev = building_heights[key]
    152         L.append((poly, elev))
    153     domain.add_quantity('elevation', Polygon_function(L, default=0.0))
    154    
    155     print 'Finished Set quantity'
     83# Counting segments
     84N = len(urs_bounding_polygon)-1
    15685
    157 ##   #------------------------------------------------------
    158 ##    # Distribute domain to implement parallelism !!!
    159 ##    #------------------------------------------------------
    160 ##
    161 ##    if numprocs > 1:
    162 ##        domain=distribute(domain)
     86# boundary tags refer to project.landward 4 points equals 5 segments start at N
     87boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4], 'ocean': range(N)}
    16388
    164     #------------------------------------------------------
    165     # Set domain parameters
    166     #------------------------------------------------------
    167     print 'domain id', id(domain)
    168     domain.set_name(kwargs['scenario_name'])
    169     domain.set_datadir(kwargs['output_dir'])
    170     domain.set_default_order(2)                 # Apply second order scheme
    171     domain.set_minimum_storable_height(0.01)    # Don't store anything less than 1cm
    172     domain.set_store_vertices_uniquely(False)
    173     domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
    174     domain.tight_slope_limiters = 1
    175     print 'domain id', id(domain)
     89#--------------------------------------------------------------------------
     90# Create the computational domain based on overall clipping polygon with
     91# a tagged boundary and interior regions defined in project.py along with
     92# resolutions (maximal area of per triangle) for each polygon
     93#--------------------------------------------------------------------------
    17694
    177     #-------------------------------------------------------------------------
    178     # Setup boundary conditions
    179     #-------------------------------------------------------------------------
    180     print 'Available boundary tags', domain.get_boundary_tags()
    181     print 'domain id', id(domain)
    182    
    183     boundary_urs_out=project.boundaries_dir + sep + project.scenario_name
     95domain = create_domain_from_regions(bounding_polygon,
     96                                    boundary_tags=boundary_tags,
     97                                    maximum_triangle_area=project.res_poly_all,
     98                                    interior_regions=project.interior_regions,
     99                                    mesh_filename=project.meshes_dir_name,
     100                                    use_cache=True,
     101                                    verbose=True)
    184102
    185     Br = Reflective_boundary(domain)
    186     Bd = Dirichlet_boundary([kwargs['tide'],0,0])
    187    
    188     print 'Available boundary tags', domain.get_boundary_tags()
    189     Bf = Field_boundary(boundary_urs_out+'.sts',  # Change from file_boundary
    190                    domain, mean_stage= project.tide,
    191                    time_thinning=1,
    192                    time_limit=project.finaltime,
    193                    default_boundary=Bd,
    194                    use_cache=True,
    195                    verbose=True,
    196                    boundary_polygon=bounding_polygon)
     103print 'memory usage after creation of domain', mem_usage()
     104print domain.statistics()
    197105
    198     domain.set_boundary({'back': Br,
    199                          'side': Bd,
    200                          'ocean': Bf})
    201106
    202     kwargs['input_start_time']=domain.starttime
     107#-------------------------------------------------------------------------
     108# Setup initial conditions
     109#-------------------------------------------------------------------------
     110print 'Setup initial conditions'
    203111
    204     print'finish set boundary'
     112# sets the initial stage in the offcoast region only
     113IC = Polygon_function([(project.poly_mainland, 0)],
     114                      default=project.tide,
     115                      geo_reference=domain.geo_reference)
     116domain.set_quantity('stage', IC)
     117domain.set_quantity('friction', project.friction)
     118domain.set_quantity('elevation',
     119                    filename=project.combined_dir_name+'.pts',
     120                    use_cache=True,
     121                    verbose=True,
     122                    alpha=project.alpha)
    205123
    206     #----------------------------------------------------------------------------
    207     # Evolve system through time
    208     #--------------------------------------------------------------------
    209     t0 = time.time()
     124# Add buildings from file
     125print 'Creating building polygons'   
     126building_polygons, building_heights = csv2building_polygons(project.building_polygon_file)
    210127
    211     for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime']
    212                        ,skip_initial_step = False):
    213         domain.write_time()
    214         domain.write_boundary_statistics(tags = 'ocean')
     128L = []
     129for key in building_polygons:
     130    poly = building_polygons[key]
     131    elev = building_heights[key]
     132    L.append((poly, elev))
    215133
    216     # these outputs should be checked with the resultant inundation map
    217     x, y = domain.get_maximum_inundation_location()
    218     q = domain.get_maximum_inundation_elevation()
    219     print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
     134buildings = Polygon_function(L, default=0.0)
     135domain.add_quantity('elevation', buildings)
    220136
    221     print 'Simulation took %.2f seconds' %(time.time()-t0)
    222137
    223     #kwargs 'completed' must be added to write the final parameters to file
    224     kwargs['completed']=str(time.time()-t0)
    225      
    226     store_parameters(**kwargs)
    227      
    228    
    229    
    230 #-------------------------------------------------------------
    231 if __name__ == "__main__":
    232    
    233     kwargs={}
    234     kwargs['file_name']=project.dir_comment
    235     kwargs['finaltime']=project.finaltime
    236     kwargs['output_dir']=project.output_run_time_dir
    237     kwargs['elevation_file']=project.combined_dir_name+'.pts'
    238     kwargs['scenario_name']=project.scenario_name
    239     kwargs['tide']=project.tide
    240     kwargs['alpha'] = project.alpha
    241     kwargs['friction']=project.friction
    242     #kwargs['num_cpu']=numprocs
    243     #kwargs['host']=project.host
    244     #kwargs['starttime']=project.starttime
    245     #kwargs['yieldstep']=project.yieldstep
    246     #kwargs['user']=project.user
    247     #kwargs['time_thinning'] = project.time_thinning
    248      
    249     run_model(**kwargs)
     138#------------------------------------------------------
     139# Distribute domain to implement parallelism !!!
     140#------------------------------------------------------
     141
     142#  if numprocs > 1:
     143#      domain=distribute(domain)
     144
     145#------------------------------------------------------
     146# Set domain parameters
     147#------------------------------------------------------
     148domain.set_name(project.scenario_name)
     149domain.set_datadir(output_dir)
     150domain.set_default_order(2)                 # Apply second order scheme
     151domain.set_minimum_storable_height(0.01)    # Don't store anything less than 1cm
     152
     153#-------------------------------------------------------------------------
     154# Setup boundary conditions
     155#-------------------------------------------------------------------------
     156print 'Set boundary conditions'
     157
     158Br = Reflective_boundary(domain)
     159Bd = Dirichlet_boundary([project.tide,0,0])
     160Bf = Field_boundary(urs_boundary_name+'.sts',
     161                    domain,
     162                    mean_stage= project.tide,
     163                    time_thinning=project.time_thinning,
     164                    default_boundary=Bd,
     165                    use_cache=True,
     166                    verbose=True,
     167                    boundary_polygon=bounding_polygon)
     168
     169domain.set_boundary({'back': Br,
     170                     'side': Bd,
     171                     'ocean': Bf})
     172
     173
     174#----------------------------------------------------------------------------
     175# Evolve system through time
     176#--------------------------------------------------------------------
     177t0 = time.time()
     178
     179for t in domain.evolve(yieldstep=project.yieldstep,
     180                       finaltime=project.finaltime)
     181    print domain.timestepping_statistics()
     182    print domain.boundary_statistics(tags='ocean')
     183
     184print 'Simulation took %.2f seconds' %(time.time()-t0)
    250185     
    251186   
Note: See TracChangeset for help on using the changeset viewer.