Changeset 5755 for anuga_work/production/carnarvon/run_carnarvon.py

Timestamp:

Sep 9, 2008, 2:18:01 PM (15 years ago)

Author:

kristy

Message:

localised all python scripts

File:

• anuga_work/production/carnarvon/run_carnarvon.py (modified) (13 diffs)

anuga_work/production/carnarvon/run_carnarvon.py

@@ -17 +17 @@
 # Standard modules
 from os import sep
+import os
 from os.path import dirname, basename
 from os import mkdir, access, F_OK
@@ -30 +31 @@
 from anuga.shallow_water import Field_boundary
 from Numeric import allclose
-from anuga.shallow_water.data_manager import export_grid
+from anuga.shallow_water.data_manager import export_grid, create_sts_boundary
 
 from anuga.pmesh.mesh_interface import create_mesh_from_regions
 from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
-from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
+#from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
 from anuga_parallel.parallel_abstraction import get_processor_name
 from anuga.caching import myhash
 from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
 from anuga.fit_interpolate.benchmark_least_squares import mem_usage 
+from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
+from anuga.geospatial_data.geospatial_data import find_optimal_smoothing_parameter
 
 # Application specific imports
 import project                 # Definition of file names and polygons
+numprocs = 1
+myid = 0
 
 def run_model(**kwargs):
@@ -62 +67 @@
         store_parameters(**kwargs)
 
-    barrier()
+   # barrier()
 
     start_screen_catcher(kwargs['output_dir'], myid, numprocs)
@@ -68 +73 @@
     print "Processor Name:",get_processor_name()
 
-    # filenames
-#    meshes_dir_name = project.meshes_dir_name+'.msh'
-
-    # creates copy of code in output dir
-    print 'min triangles', project.trigs_min,
-    print 'Note: This is generally about 20% less than the final amount'
-
+    
+    #-----------------------------------------------------------------------
+    # Domain definitions
+    #-----------------------------------------------------------------------
+
+    # Read in boundary from ordered sts file
+    urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir,project.scenario_name))
+
+    # Reading the landward defined points, this incorporates the original clipping
+    # polygon minus the 100m contour
+    landward_bounding_polygon = read_polygon(project.boundaries_dir+'landward_bounding_polygon.txt')
+
+    # Combine sts polyline with landward points
+    bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
+    
+    # counting segments
+    N = len(urs_bounding_polygon)-1
+    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
@@ -88 +106 @@
         print 'start create mesh from regions'
 
-        create_mesh_from_regions(project.poly_all,
-                             boundary_tags=project.boundary_tags,
+        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+'.msh',
-                             use_cache=False,
+                             use_cache=True,
                              verbose=True)
-    barrier()
-    
+   # barrier()
+
+##        covariance_value,alpha = find_optimal_smoothing_parameter (data_file= kwargs['elevation_file'],
+##                                alpha_list=[0.001, 0.01, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5],
+##                                mesh_file = project.meshes_dir_name+'.msh')
+##        print 'optimal alpha', covariance_value,alpha       
 
     #-------------------------------------------------------------------------
@@ -103 +125 @@
     print 'Setup computational domain'
 
-    #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
-    #above don't work
     domain = Domain(project.meshes_dir_name+'.msh', use_cache=False, verbose=True)
     print 'memory usage before del domain',mem_usage()
@@ -123 +143 @@
         from polygon import Polygon_function
         #following sets the stage/water to be offcoast only
-        IC = Polygon_function( [(project.poly_mainland, -1.0)], default = kwargs['tide'],
+        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('stage',kwargs['tide'] )
         domain.set_quantity('friction', kwargs['friction']) 
         
-        print 'Start Set quantity',kwargs['bathy_file']
+        print 'Start Set quantity',kwargs['elevation_file']
 
         domain.set_quantity('elevation', 
-                            filename = kwargs['bathy_file'],
-                            use_cache = True,
+                            filename = kwargs['elevation_file'],
+                            use_cache = False,
                             verbose = True,
                             alpha = kwargs['alpha'])
         print 'Finished Set quantity'
-    barrier()
-
+    #barrier()
+
+##    #------------------------------------------------------
+##    # Create x,y,z file of mesh vertex!!!
+##    #------------------------------------------------------
+##        coord = domain.get_vertex_coordinates()
+##        depth = domain.get_quantity('elevation')
+##        
+##        # Write vertex coordinates to file
+##        filename=project.vertex_filename
+##        fid=open(filename,'w')
+##        fid.write('x (m), y (m), z(m)\n')
+##        for i in range(len(coord)):
+##            pt=coord[i]
+##            x=pt[0]
+##            y=pt[1]
+##            z=depth[i]
+##            fid.write('%.6f,%.6f,%.6f\n' %(x, y, z))
+##
 
     #------------------------------------------------------
@@ -157 +194 @@
     domain.set_store_vertices_uniquely(False)
     domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
-    domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
+    domain.tight_slope_limiters = 1
     print 'domain id', id(domain)
-
 
     #-------------------------------------------------------------------------
@@ -166 +202 @@
     print 'Available boundary tags', domain.get_boundary_tags()
     print 'domain id', id(domain)
-    #print 'Reading Boundary file',project.boundaries_dir_namea + '.sww'
-
-    if project.source != 'other':
-        Bf = Field_boundary(kwargs['boundary_file'],
-                    domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'], 
-                    use_cache=True, verbose=True)
-                    
-    kwargs['input_start_time']=domain.starttime
-
-    print 'finished reading boundary file'
+   
+    boundary_urs_out=project.boundaries_dir_name
 
     Br = Reflective_boundary(domain)
     Bd = Dirichlet_boundary([kwargs['tide'],0,0])
-    Bo = Dirichlet_boundary([kwargs['tide']+5.0,0,0])
-
-
-    print'set_boundary'
+    
+    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,
+                   default_boundary=Bd,
+                   use_cache=True,
+                   verbose = True,
+                   boundary_polygon=bounding_polygon)
 
     domain.set_boundary({'back': Br,
                          'side': Bd,
                          'ocean': Bf}) 
+
+    kwargs['input_start_time']=domain.starttime
+
     print'finish set boundary'
 
@@ -194 +230 @@
     t0 = time.time()
 
-    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['starttime']): 
+    for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime']
+                       ,skip_initial_step = False): 
         domain.write_time()
-        domain.write_boundary_statistics(tags = 'ocean')     
-
-        if allclose(t, 120):
-            domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bo})
-
-        if allclose(t, 720):
-            domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bd})
-
-#    for t in domain.evolve(yieldstep = kwargs['yieldstep'], finaltime = kwargs['midtime']
-#                       ,skip_initial_step = True): 
-#        domain.write_time()
-#        domain.write_boundary_statistics(tags = 'ocean')     
-#    
-#    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime']
-#                       ,skip_initial_step = True): 
-#        domain.write_time()
-#        domain.write_boundary_statistics(tags = 'ocean')   
+        domain.write_boundary_statistics(tags = 'ocean')
 
     x, y = domain.get_maximum_inundation_location()
@@ -226 +247 @@
     if myid==0:
         store_parameters(**kwargs)
-    barrier
+   # barrier
     
     print 'memory usage before del domain1',mem_usage()
     
-def export_model(**kwargs):    
-    #store_parameters(**kwargs)
-    
-#    print 'memory usage before del domain',mem_usage()
-    #del domain
-    print 'memory usage after del domain',mem_usage()
-    
-    swwfile = kwargs['output_dir']+kwargs['aa_scenario_name']
-    print'swwfile',swwfile
-    
-    export_grid(swwfile, extra_name_out = 'town',
-            quantities = ['speed','depth','elevation','stage'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
-            #quantities = ['speed','depth'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
-            timestep = None,
-            reduction = max,
-            cellsize = kwargs['export_cellsize'],
-            NODATA_value = -9999,
-            easting_min = project.eastingmin,
-            easting_max = project.eastingmax,
-            northing_min = project.northingmin,
-            northing_max = project.northingmax,
-            verbose = False,
-            origin = None,
-            datum = 'GDA94',
-            format = 'asc')
-            
-    inundation_damage(swwfile+'.sww', project.buildings_filename, 
-                      project.buildings_filename_out)
     
 #-------------------------------------------------------------
@@ -269 +262 @@
     kwargs['starttime']=project.starttime
     kwargs['yieldstep']=project.yieldstep
-    kwargs['midtime']=project.midtime
     kwargs['finaltime']=project.finaltime
     
     kwargs['output_dir']=project.output_run_time_dir
-    kwargs['bathy_file']=project.combined_dir_name
-#    kwargs['bathy_file']=project.combined_small_dir_name + '.pts'
+    kwargs['elevation_file']=project.combined_dir_name+'.txt'
     kwargs['boundary_file']=project.boundaries_in_dir_name + '.sww'
     kwargs['file_name']=project.home+'detail.csv'
@@ -291 +282 @@
     run_model(**kwargs)
       
-    if myid==0:
-        export_model(**kwargs)
-    barrier
+    #barrier