Changeset 6178

Timestamp:

Jan 15, 2009, 5:45:46 PM (15 years ago)

Author:

ole

Message:

Clean up of Patong script and added new function create_domain_from_regions

Files:

• anuga_core/source/anuga/abstract_2d_finite_volumes/util.py (modified) (1 diff)
• anuga_core/source/anuga/pmesh/mesh_interface.py (modified) (1 diff)
• anuga_core/source/anuga/shallow_water/__init__.py (modified) (1 diff)
• anuga_core/source/anuga/shallow_water/shallow_water_domain.py (modified) (2 diffs)
• anuga_work/production/patong/project.py (modified) (1 diff)
• anuga_work/production/patong/run_patong.py (modified) (3 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

@@ -348 +348 @@
     
     if time_limit is not None:
+        #if verbose is True:
+        #    print '****** Time limit', time_limit
+        #    print '****** Start time', starttime
+        #    print '****** Time in ', time[0], time[-1]
+        
         # Adjust given time limit to given start time
         time_limit = time_limit - starttime
+
 
         # Find limit point

anuga_core/source/anuga/pmesh/mesh_interface.py

@@ -90 +90 @@
               'minimum_triangle_angle': minimum_triangle_angle,
               'fail_if_polygons_outside': fail_if_polygons_outside,
-              'verbose': verbose}   # FIXME (Ole): Should be bypassed one day
-                                    # What should be bypassed? Verbose?
-    
-    #print 'kwargs', kwargs
+              'verbose': verbose}   # FIXME (Ole): Should be bypassed one day. See ticket:14
 
     # Call underlying engine with or without caching

anuga_core/source/anuga/shallow_water/__init__.py

@@ -8 +8 @@
 # Make selected classes available directly
 from shallow_water_domain import Domain,\
-    Transmissive_boundary, Reflective_boundary,\
-    Dirichlet_boundary, Time_boundary, File_boundary,\
-    Transmissive_momentum_set_stage_boundary,\
-    Dirichlet_discharge_boundary,\
-    Field_boundary
+     create_domain_from_regions,\
+     Transmissive_boundary, Reflective_boundary,\
+     Dirichlet_boundary, Time_boundary, File_boundary,\
+     Transmissive_momentum_set_stage_boundary,\
+     Dirichlet_discharge_boundary,\
+     Field_boundary
 
 

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -94 +94 @@
 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
      import Transmissive_boundary
-
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
 from anuga.utilities.numerical_tools import gradient, mean, ensure_numeric
 from anuga.geospatial_data.geospatial_data import ensure_geospatial
@@ -117 +117 @@
 
 
-#---------------------
+
+#---------------------------
+# Create domain from regions
+#---------------------------
+
+def create_domain_from_regions(bounding_polygon,
+                               boundary_tags,
+                               maximum_triangle_area=None,
+                               mesh_filename=None,
+                               interior_regions=None,
+                               interior_holes=None,
+                               poly_geo_reference=None,
+                               mesh_geo_reference=None,
+                               minimum_triangle_angle=28.0,
+                               fail_if_polygons_outside=True,
+                               use_cache=False,
+                               verbose=True):
+    """Create domain from bounding polygons and resolutions.
+
+    bounding_polygon is a list of points in Eastings and Northings,
+    relative to the poly_geo_reference.
+
+    Boundary tags is a dictionary of symbolic tags. For every tag there
+    is a list of indices referring to segments associated with that tag.
+    If a segment is omitted it will be assigned the default tag ''.
+
+    maximum_triangle_area is the maximal area per triangle
+    for the bounding polygon, excluding the  interior regions.
+
+    Interior_regions is a list of tuples consisting of (polygon,
+    resolution) for each region to be separately refined. Do not have
+    polygon lines cross or be on-top of each other.  Also do not have
+    polygon close to each other.
+    
+    NOTE: If a interior_region is outside the bounding_polygon it should 
+    throw an error
+    
+    Interior_holes is a list of ploygons for each hole.
+
+    This function does not allow segments to share points - use underlying
+    pmesh functionality for that
+
+    poly_geo_reference is the geo_reference of the bounding polygon and
+    the interior polygons.
+    If none, assume absolute.  Please pass one though, since absolute
+    references have a zone.
+    
+    mesh_geo_reference is the geo_reference of the mesh to be created.
+    If none is given one will be automatically generated.  It was use
+    the lower left hand corner of  bounding_polygon (absolute)
+    as the x and y values for the geo_ref.
+    
+    Returns the shallow water domain instance
+
+    Note, interior regions should be fully nested, as overlaps may cause
+    unintended resolutions. 
+
+    fail_if_polygons_outside: If True (the default) Exception in thrown
+    where interior polygons fall outside bounding polygon. If False, these
+    will be ignored and execution continued.
+        
+    
+    """
+
+
+    # Build arguments and keyword arguments for use with caching or apply.
+    args = (bounding_polygon,
+            boundary_tags)
+    
+    kwargs = {'maximum_triangle_area': maximum_triangle_area,
+              'mesh_filename': mesh_filename,
+              'interior_regions': interior_regions,
+              'interior_holes': interior_holes,
+              'poly_geo_reference': poly_geo_reference,
+              'mesh_geo_reference': mesh_geo_reference,
+              'minimum_triangle_angle': minimum_triangle_angle,
+              'fail_if_polygons_outside': fail_if_polygons_outside,
+              'verbose': verbose} #FIXME (Ole): See ticket:14
+
+    # Call underlying engine with or without caching
+    if use_cache is True:
+        try:
+            from anuga.caching import cache
+        except:
+            msg = 'Caching was requested, but caching module'+\
+                  'could not be imported'
+            raise msg
+
+
+        domain = cache(_create_domain_from_regions,
+                       args, kwargs,
+                       verbose=verbose,
+                       compression=False)
+    else:
+        domain = apply(_create_domain_from_regions,
+                       args, kwargs)
+
+    return domain
+
+        
+def _create_domain_from_regions(bounding_polygon,
+                                boundary_tags,
+                                maximum_triangle_area=None,
+                                mesh_filename=None,                           
+                                interior_regions=None,
+                                interior_holes=None,
+                                poly_geo_reference=None,
+                                mesh_geo_reference=None,
+                                minimum_triangle_angle=28.0,
+                                fail_if_polygons_outside=True,
+                                verbose=True):
+    """_create_domain_from_regions - internal function.
+
+    See create_domain_from_regions for documentation.
+    """
+
+    create_mesh_from_regions(bounding_polygon,
+                             boundary_tags,
+                             maximum_triangle_area=maximum_triangle_area,
+                             interior_regions=interior_regions,
+                             filename=mesh_filename,
+                             interior_holes=interior_holes,
+                             poly_geo_reference=poly_geo_reference,
+                             mesh_geo_reference=mesh_geo_reference,
+                             minimum_triangle_angle=minimum_triangle_angle,
+                             fail_if_polygons_outside=fail_if_polygons_outside,
+                             use_cache=False,
+                             verbose=verbose)
+
+    domain = Domain(mesh_filename, use_cache=False, verbose=verbose)
+
+
+    return domain
+    
+
+
+
+
+
+#
 # Shallow water domain
 #---------------------

anuga_work/production/patong/project.py

@@ -48 +48 @@
 finaltime=15000         # final time for simulation 15000
 
-setup='final'  # Final can be replaced with trial or basic.
+setup='trial'  # Final can be replaced with trial or basic.
                # Either will result in a coarser mesh that will allow a
                # faster, but less accurate, simulation.

anuga_work/production/patong/run_patong.py

@@ -23 +23 @@
 from os import sep
 import os
-from os.path import dirname, basename
-from os import mkdir, access, F_OK
-from shutil import copy
+from os.path import dirname
 import time
-import sys
 
 # Related major packages
@@ -35 +32 @@
 from anuga.shallow_water import Reflective_boundary
 from anuga.shallow_water import Field_boundary
-from Numeric import allclose
 from anuga.shallow_water.data_manager import export_grid, create_sts_boundary, csv2building_polygons
-from anuga.pmesh.mesh_interface import create_mesh_from_regions
+from anuga.shallow_water import create_domain_from_regions
 from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
 from anuga.caching import myhash
@@ -51 +47 @@
 myid = 0
 
-def run_model(**kwargs):
-    
-    #------------------------------------------------------------------------------
-    # Copy scripts to time stamped output directory and capture screen
-    # output to file
-    #------------------------------------------------------------------------------
-
-    #copy script must be before screen_catcher
-
-    print 'output_dir',kwargs['output_dir']
-    
-    copy_code_files(kwargs['output_dir'],__file__, 
-             dirname(project.__file__)+sep+ project.__name__+'.py' )
-
-    store_parameters(**kwargs)
-
-    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
 
     
-    #-----------------------------------------------------------------------
-    # Domain definitions
-    #-----------------------------------------------------------------------
+#------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#------------------------------------------------------------------------------
 
-    # Read in boundary from ordered sts file
-    urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir,project.scenario_name))
+#copy script must be before screen_catcher
 
-    # Reading the landward defined points, this incorporates the original clipping
-    # polygon minus the 100m contour
-    landward_bounding_polygon = read_polygon(project.landward_dir)
+output_dir = project.output_run_time_dir
+print 'output_dir', output_dir
 
-    # Combine sts polyline with landward points
-    bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
-    
-    # counting segments
-    N = len(urs_bounding_polygon)-1
+#copy_code_files(output_dir,__file__, 
+#         dirname(project.__file__)+sep+ project.__name__+'.py' )
 
-    # boundary tags refer to project.landward 4 points equals 5 segments start at N
-    boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4], 'ocean': range(N)}
-
-    #--------------------------------------------------------------------------
-    # Create the triangular mesh based on overall clipping polygon with a tagged
-    # boundary and interior regions defined in project.py along with
-    # resolutions (maximal area of per triangle) for each polygon
-    #--------------------------------------------------------------------------
-
-    # IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
-    # causes problems with the ability to cache set quantity which takes alot of times
-        
-    print 'start create mesh from regions'
-
-    create_mesh_from_regions(bounding_polygon,
-                             boundary_tags=boundary_tags,
-                             maximum_triangle_area=project.res_poly_all,
-                             interior_regions=project.interior_regions,
-                             filename=project.meshes_dir_name,
-                             use_cache=True,
-                             verbose=True)
-   
-    #-------------------------------------------------------------------------
-    # Setup computational domain
-    #-------------------------------------------------------------------------
-    print 'Setup computational domain'
-
-    domain = Domain(project.meshes_dir_name, use_cache=False, verbose=True)
-    print 'memory usage before del domain',mem_usage()
-       
-    print domain.statistics()
-    print 'triangles',len(domain)
-    
-    kwargs['act_num_trigs']=len(domain)
+#start_screen_catcher(output_dir, myid, numprocs)
 
 
-    #-------------------------------------------------------------------------
-    # Setup initial conditions
-    #-------------------------------------------------------------------------
-    print 'Setup initial conditions'
+#-----------------------------------------------------------------------
+# Domain definitions
+#-----------------------------------------------------------------------
 
-    # sets the initial stage in the offcoast region only
-    IC = Polygon_function( [(project.poly_mainland, 0)], default = kwargs['tide'],
-                             geo_reference = domain.geo_reference)
-    domain.set_quantity('stage', IC)
-    #domain.set_quantity('stage',kwargs['tide'] )
-    domain.set_quantity('friction', kwargs['friction']) 
-    
-    print 'Start Set quantity',kwargs['elevation_file']
+# Read in boundary from ordered sts file
+urs_boundary_name = os.path.join(project.boundaries_dir,
+                                 project.scenario_name)
+urs_bounding_polygon = create_sts_boundary(urs_boundary_name)
 
-    domain.set_quantity('elevation', 
-                        filename = kwargs['elevation_file'],
-                        use_cache = False,
-                        verbose = True,
-                        alpha = kwargs['alpha'])
+# Reading the landward defined points, this incorporates the original clipping
+# polygon minus the 100m contour
+landward_bounding_polygon = read_polygon(project.landward_dir)
 
-    # Add buildings from file
-    building_polygons, building_heights = csv2building_polygons(project.building_polygon_file)
+# Combine sts polyline with landward points
+bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
 
-    L = []
-    for key in building_polygons:
-        poly = building_polygons[key]
-        elev = building_heights[key]
-        L.append((poly, elev))
-    domain.add_quantity('elevation', Polygon_function(L, default=0.0))
-    
-    print 'Finished Set quantity'
+# Counting segments
+N = len(urs_bounding_polygon)-1
 
-##   #------------------------------------------------------
-##    # Distribute domain to implement parallelism !!!
-##    #------------------------------------------------------
-##
-##    if numprocs > 1:
-##        domain=distribute(domain)
+# boundary tags refer to project.landward 4 points equals 5 segments start at N
+boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4], 'ocean': range(N)}
 
-    #------------------------------------------------------
-    # Set domain parameters
-    #------------------------------------------------------
-    print 'domain id', id(domain)
-    domain.set_name(kwargs['scenario_name'])
-    domain.set_datadir(kwargs['output_dir'])
-    domain.set_default_order(2)                 # Apply second order scheme
-    domain.set_minimum_storable_height(0.01)    # Don't store anything less than 1cm
-    domain.set_store_vertices_uniquely(False)
-    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
-    domain.tight_slope_limiters = 1
-    print 'domain id', id(domain)
+#--------------------------------------------------------------------------
+# Create the computational domain based on overall clipping polygon with
+# a tagged boundary and interior regions defined in project.py along with
+# resolutions (maximal area of per triangle) for each polygon
+#--------------------------------------------------------------------------
 
-    #-------------------------------------------------------------------------
-    # Setup boundary conditions 
-    #-------------------------------------------------------------------------
-    print 'Available boundary tags', domain.get_boundary_tags()
-    print 'domain id', id(domain)
-   
-    boundary_urs_out=project.boundaries_dir + sep + project.scenario_name
+domain = create_domain_from_regions(bounding_polygon,
+                                    boundary_tags=boundary_tags,
+                                    maximum_triangle_area=project.res_poly_all,
+                                    interior_regions=project.interior_regions,
+                                    mesh_filename=project.meshes_dir_name,
+                                    use_cache=True,
+                                    verbose=True)
 
-    Br = Reflective_boundary(domain)
-    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
-    
-    print 'Available boundary tags', domain.get_boundary_tags()
-    Bf = Field_boundary(boundary_urs_out+'.sts',  # Change from file_boundary
-                   domain, mean_stage= project.tide,
-                   time_thinning=1,
-                   time_limit=project.finaltime,
-                   default_boundary=Bd,
-                   use_cache=True,
-                   verbose=True,
-                   boundary_polygon=bounding_polygon)
+print 'memory usage after creation of domain', mem_usage()
+print domain.statistics()
 
-    domain.set_boundary({'back': Br,
-                         'side': Bd,
-                         'ocean': Bf}) 
 
-    kwargs['input_start_time']=domain.starttime
+#-------------------------------------------------------------------------
+# Setup initial conditions
+#-------------------------------------------------------------------------
+print 'Setup initial conditions'
 
-    print'finish set boundary'
+# sets the initial stage in the offcoast region only
+IC = Polygon_function([(project.poly_mainland, 0)],
+                      default=project.tide,
+                      geo_reference=domain.geo_reference)
+domain.set_quantity('stage', IC)
+domain.set_quantity('friction', project.friction) 
+domain.set_quantity('elevation', 
+                    filename=project.combined_dir_name+'.pts',
+                    use_cache=True,
+                    verbose=True,
+                    alpha=project.alpha)
 
-    #----------------------------------------------------------------------------
-    # Evolve system through time
-    #--------------------------------------------------------------------
-    t0 = time.time()
+# Add buildings from file
+print 'Creating building polygons'    
+building_polygons, building_heights = csv2building_polygons(project.building_polygon_file)
 
-    for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime']
-                       ,skip_initial_step = False): 
-        domain.write_time()
-        domain.write_boundary_statistics(tags = 'ocean')
+L = []
+for key in building_polygons:
+    poly = building_polygons[key]
+    elev = building_heights[key]
+    L.append((poly, elev))
 
-    # these outputs should be checked with the resultant inundation map
-    x, y = domain.get_maximum_inundation_location()
-    q = domain.get_maximum_inundation_elevation()
-    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
+buildings = Polygon_function(L, default=0.0)
+domain.add_quantity('elevation', buildings)
 
-    print 'Simulation took %.2f seconds' %(time.time()-t0)
 
-    #kwargs 'completed' must be added to write the final parameters to file
-    kwargs['completed']=str(time.time()-t0)
-     
-    store_parameters(**kwargs)
-     
-    
-    
-#-------------------------------------------------------------
-if __name__ == "__main__":
-    
-    kwargs={}
-    kwargs['file_name']=project.dir_comment
-    kwargs['finaltime']=project.finaltime
-    kwargs['output_dir']=project.output_run_time_dir
-    kwargs['elevation_file']=project.combined_dir_name+'.pts'
-    kwargs['scenario_name']=project.scenario_name
-    kwargs['tide']=project.tide
-    kwargs['alpha'] = project.alpha
-    kwargs['friction']=project.friction
-    #kwargs['num_cpu']=numprocs
-    #kwargs['host']=project.host
-    #kwargs['starttime']=project.starttime
-    #kwargs['yieldstep']=project.yieldstep
-    #kwargs['user']=project.user
-    #kwargs['time_thinning'] = project.time_thinning
-      
-    run_model(**kwargs)
+#------------------------------------------------------
+# Distribute domain to implement parallelism !!!
+#------------------------------------------------------
+
+#  if numprocs > 1:
+#      domain=distribute(domain)
+
+#------------------------------------------------------
+# Set domain parameters
+#------------------------------------------------------
+domain.set_name(project.scenario_name)
+domain.set_datadir(output_dir)
+domain.set_default_order(2)                 # Apply second order scheme
+domain.set_minimum_storable_height(0.01)    # Don't store anything less than 1cm
+
+#-------------------------------------------------------------------------
+# Setup boundary conditions 
+#-------------------------------------------------------------------------
+print 'Set boundary conditions'
+
+Br = Reflective_boundary(domain)
+Bd = Dirichlet_boundary([project.tide,0,0])
+Bf = Field_boundary(urs_boundary_name+'.sts', 
+                    domain,
+                    mean_stage= project.tide,
+                    time_thinning=project.time_thinning,
+                    default_boundary=Bd,
+                    use_cache=True,
+                    verbose=True,
+                    boundary_polygon=bounding_polygon)
+
+domain.set_boundary({'back': Br,
+                     'side': Bd,
+                     'ocean': Bf}) 
+
+
+#----------------------------------------------------------------------------
+# Evolve system through time
+#--------------------------------------------------------------------
+t0 = time.time()
+
+for t in domain.evolve(yieldstep=project.yieldstep,
+                       finaltime=project.finaltime)
+    print domain.timestepping_statistics()
+    print domain.boundary_statistics(tags='ocean')
+
+print 'Simulation took %.2f seconds' %(time.time()-t0)