Changeset 4570

Timestamp:

Jun 28, 2007, 5:54:00 PM (18 years ago)

Author:

ole

Message:

Setup shark bay initial runs

Location:

anuga_work/production/shark_bay_2007

Files:

• build_shark_bay.py (modified) (5 diffs)
• project.py (modified) (7 diffs)
• run_shark_bay.py (modified) (12 diffs)

anuga_work/production/shark_bay_2007/build_shark_bay.py

@@ -28 +28 @@
 from anuga.shallow_water import File_boundary
 from anuga.shallow_water import Reflective_boundary
-from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf
+from anuga.shallow_water.data_manager import dem2pts
 from anuga.pmesh.mesh_interface import create_mesh_from_regions
 from anuga.geospatial_data.geospatial_data import Geospatial_data
 from anuga.shallow_water.data_manager import start_screen_catcher
 from anuga.shallow_water.data_manager import copy_code_files
-
+from anuga.shallow_water.data_manager import urs_ungridded2sww
 from anuga_parallel.parallel_abstraction import get_processor_name
 
@@ -46 +47 @@
 
 copy_code_files(project.output_build_time_dir,__file__, 
-               dirname(project.__file__)+sep+ project.__name__+'.py' )
+               dirname(project.__file__)+sep+ project.__name__+'.py' )#
 
 start_screen_catcher(project.output_build_time_dir)
@@ -67 +68 @@
 geospatial_data = None
 # create DEMs from asc data
+
+# FIXME (Ole): Clip data by interior regions if possible
+
 print 'creating geospatial data objects from asc data (via dem and pts formats)'
 for filename in project.ascii_grid_filenames:
@@ -74 +78 @@
     dem2pts(filename, use_cache=True, verbose=True)
 
-    geospatial_data += Geospatial_data(file_name = filename + '.pts', verbose=True)
+    geospatial_data += Geospatial_data(file_name=filename + '.pts',
+                                       verbose=True)
 
 
 print 'creating geospatial data objects from txt data'
 for filename in project.point_filenames:
-    geospatial_data += Geospatial_data(file_name = filename + '.txt', verbose=True)
+    geospatial_data += Geospatial_data(file_name=filename + '.txt',
+                                       verbose=True)
 
 
 print 'clip combined geospatial object by bounding polygon'
-G = geospatial_data.clip(project.poly_all)
+G = geospatial_data.clip(project.bounding_polygon)
 
 
@@ -104 +110 @@
 #-------------------------------------------------------------------------
 print 'converting boundary conditions to sww format'
-#boundaries_in_dir_name = project.boundaries_in_dir_name
 
-#print 'minlat=project.north_boundary, maxlat=project.south_boundary',project.north_boundary, project.south_boundary
-#print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary
-       
-from anuga.shallow_water.data_manager import urs_ungridded2sww
+print 'boundary_dir', project.boundary_dir
+print 'project.boundary_name', project.boundary_name
 
-print 'boundaries_dir', project.boundaries_dir
-print 'project.boundaries_name', project.boundaries_name
-
-urs_ungridded2sww(project.boundaries_name,
-                  verbose=True, mint=5000, maxt=35000, zscale=1)
+urs_ungridded2sww(project.boundary_name, verbose=True);
+                  #verbose=True, mint=5000, maxt=35000, zscale=1)
 
 

anuga_work/production/shark_bay_2007/project.py

@@ -32 +32 @@
 
 #tide = -3.9
-#tide = 0.0
-tide = 3.6
+tide = 0.0
+#tide = 3.6
 
 #Maybe will try to make project a class to allow these parameters to be passed in.
@@ -40 +40 @@
 finaltime=25000
 starttime=3600
-setup='final'
+setup='basic'
 source='shark_bay'
 
@@ -67 +67 @@
                    'clipped_bathymetry_final', 'coast_points_final']
 
-#final topo name
+# final topo name
 combined_name ='shark_bay_combined_elevation'
 combined_small_name = 'shark_bay_combined_elevation_small'
@@ -74 +74 @@
 anuga_dir = join(home, state, scenario, 'anuga')
 
-topographies_in_dir = join(home, state, scenario, 'elevation_final', 'test_area')
+topographies_in_dir = join(home, state, scenario,
+                           'elevation_final', 'test_area')
+
 topographies_dir = join(anuga_dir, 'topographies') # Output dir for ANUGA
 
@@ -87 +89 @@
 
 
-meshes_dir = join(anuga_dir, 'meshes')
-meshes_dir_name = join(meshes_dir, scenario_name)
+mesh_dir = join(anuga_dir, 'meshes')
+mesh_name = join(mesh_dir, scenario_name)
 
 
@@ -97 +99 @@
 # locations for boundary conditions
 if source=='shark_bay':
-    boundaries_file_name = 'shark_bay_3867_18052007'
-    boundaries_dir = join(anuga_dir, 'boundaries', 'shark_bay', '1_10000')
+    boundary_file_name = 'shark_bay_3867_18052007'
+    boundary_dir = join(anuga_dir, 'boundaries', 'shark_bay', '1_10000')
 
-boundaries_name = join(boundaries_dir, boundaries_file_name)
+boundary_name = join(boundary_dir, boundary_file_name)
 
 #output locations
-output_dir = join(anuga_dir, 'outputs')
+output_dir = join(anuga_dir, 'outputs')+sep
 output_build_time_dir = output_dir+build_time+sep
 output_run_time_dir = output_dir +run_time+dir_comment+sep
@@ -147 +149 @@
 
 #from anuga.coordinate_transforms.redfearn import redfearn
-poly_all = read_polygon(join(polygons_dir, 'boundary_extent.csv'))
-res_poly_all = 250000*res_factor
+bounding_polygon = read_polygon(join(polygons_dir, 'boundary_extent.csv'))
+res_bounding_polygon = 250000*res_factor
 
 #Interior regions
 poly_inundated_area = read_polygon(join(polygons_dir, 'inundated_area.csv'))
-res_inundated_area = 50000*res_factor
+res_inundated_area = 1000*res_factor
 
 poly_bay_area = read_polygon(join(polygons_dir, 'bay_area.csv'))                                   
-res_bay_area = 100000*res_factor                                   
+res_bay_area = 5000*res_factor                                   
 
-interior_regions = [[poly_bay_area, res_bay_area], [poly_inundated_area, res_inundated_area]]
+interior_regions = [[poly_bay_area, res_bay_area],
+                    [poly_inundated_area, res_inundated_area]]
 
-trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
+trigs_min = number_mesh_triangles(interior_regions,
+                                  bounding_polygon,
+                                  res_bounding_polygon)
 
 print 'min number triangles', trigs_min
 
+onshore_polygon = read_polygon(join(topographies_in_dir,
+                                    'initial_condition.txt'))
 

anuga_work/production/shark_bay_2007/run_shark_bay.py

@@ -32 +32 @@
 
 from anuga.pmesh.mesh_interface import create_mesh_from_regions
-from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files,store_parameters
+from anuga.shallow_water.data_manager import start_screen_catcher
+from anuga.shallow_water.data_manager import copy_code_files
+from anuga.shallow_water.data_manager import store_parameters
 from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
 from anuga_parallel.parallel_abstraction import get_processor_name
 from anuga.caching import myhash
+
 # Application specific imports
-import project_urs                 # Definition of file names and polygons
+import project                 # Definition of file names and polygons
 
 def run_model(**kwargs):
@@ -54 +57 @@
     #copy script must be before screen_catcher
     print 'tide',tide
-    kwargs['est_num_trigs']=project_urs.trigs_min
+    kwargs['est_num_trigs']=project.trigs_min
     kwargs['num_cpu']=numprocs
-    kwargs['host']=project_urs.host
-    kwargs['res_factor']=project_urs.res_factor
-    kwargs['starttime']=project_urs.starttime
-    kwargs['yieldstep']=project_urs.yieldstep
-    kwargs['finaltime']=project_urs.finaltime
-    
-    kwargs['output_dir']=project_urs.output_run_time_dir
-    kwargs['bathy_file']=project_urs.combined_dir_name + '.txt'
-#    kwargs['bathy_file']=project_urs.combined_small_dir_name + '.pts'
-    kwargs['boundary_file']=project_urs.boundaries_dir_name + '.sww'
+    kwargs['host']=project.host
+    kwargs['res_factor']=project.res_factor
+    kwargs['starttime']=project.starttime
+    kwargs['yieldstep']=project.yieldstep
+    kwargs['finaltime']=project.finaltime
+    
+    kwargs['output_dir']=project.output_run_time_dir
+    kwargs['bathy_file']=project.combined_dir_name + '.txt'
+#    kwargs['bathy_file']=project.combined_small_dir_name + '.pts'
+    kwargs['boundary_file']=project.boundary_name + '.sww'
 #    kwargs['Completed']=''
 
@@ -71 +74 @@
     if myid == 0:
         copy_code_files(kwargs['output_dir'],__file__, 
-                 dirname(project_urs.__file__)+sep+ project_urs.__name__+'.py' )
+                 dirname(project.__file__)+sep+ project.__name__+'.py' )
 
         store_parameters(**kwargs)
@@ -79 +82 @@
     start_screen_catcher(kwargs['output_dir'], myid, numprocs)
 
-    print "Processor Name:",get_processor_name()
+    print 'Processor Name:', get_processor_name()
 
     # creates copy of code in output dir
-    print 'min triangles', project_urs.trigs_min,
+    print 'min triangles', project.trigs_min,
     print 'Note: This is generally about 20% less than the final amount'
 
@@ -88 +91 @@
     # Create the triangular mesh based on overall clipping polygon with a
     # tagged
-    # boundary and interior regions defined in project_urs.py along with
+    # boundary and interior regions defined in project.py along with
     # resolutions (maximal area of per triangle) for each polygon
     #--------------------------------------------------------------------------
@@ -97 +100 @@
     
         print 'start create mesh from regions'
-
-        create_mesh_from_regions(project_urs.poly_all,
-                             boundary_tags={'back': [3, 4, 5, 6], 'side': [2, 7],
-                                            'ocean': [0, 1]},
-                             maximum_triangle_area=project_urs.res_poly_all,
-                             interior_regions=project_urs.interior_regions,
-                             filename=project_urs.meshes_dir_name+'.msh',
-                             use_cache=False,
-                             verbose=True)
-    barrier()
+        
+        # FIXME (Ole): What if tags are wrong?
+        create_mesh_from_regions(project.bounding_polygon,
+                                 boundary_tags={'back': [1, 2, 3, 4, 5],
+                                                'side': [0, 6],
+                                                'ocean': [7, 8, 9, 10, 11]},
+                                 maximum_triangle_area=project.res_bounding_polygon,
+                                 interior_regions=project.interior_regions,
+                                 filename=project.mesh_name+'.msh',
+                                 use_cache=False,
+                                 verbose=True)
 
     #-------------------------------------------------------------------------
@@ -113 +117 @@
     print 'Setup computational domain'
 
-    #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
+    #domain = cache(Domain, (mesh_name), {'use_cache':True, 'verbose':True}, verbose=True)
     #above don't work
-    domain = Domain(project_urs.meshes_dir_name+'.msh', use_cache=False, verbose=True)
+    domain = Domain(project.mesh_name+'.msh',
+                    use_cache=False, verbose=True)
        
     print domain.statistics()
@@ -127 +132 @@
     if myid == 0:
 
-        print 'Setup initial conditions'
-
+        print 'Create onshore polygon'
         from polygon import Polygon_function
         #following sets the stage/water to be offcoast only
-        IC = Polygon_function( [(project_urs.poly_mainland, -1.0)], default = tide,
-                                 geo_reference = domain.geo_reference)
+        IC = Polygon_function([(project.onshore_polygon, -1.0)],
+                              default = tide,
+                              geo_reference = domain.geo_reference)
+
+        print 'set initial condition'
         domain.set_quantity('stage', IC)
         domain.set_quantity('friction', friction) 
-        
-        print 'Start Set quantity'
-
         domain.set_quantity('elevation', 
                             filename = kwargs['bathy_file'],
@@ -143 +147 @@
                             verbose = True,
                             alpha = alpha)
-        print 'Finished Set quantity'
+        
+
+    #------------------------------------------------------
+    # Distribute domain to implement parallelism !!!
+    #------------------------------------------------------
     barrier()
-
-    #------------------------------------------------------
-    # Distribute domain to implement parallelism !!!
-    #------------------------------------------------------
-
     if numprocs > 1:
         domain=distribute(domain)
@@ -167 +170 @@
     domain.beta_h = 0
     #domain.limit2007 = 1
+    
 
     #-------------------------------------------------------------------------
@@ -173 +177 @@
     print 'Available boundary tags', domain.get_boundary_tags()
     print 'domain id', id(domain)
-    #print 'Reading Boundary file',project_urs.boundaries_dir_namea + '.sww'
+    #print 'Reading Boundary file',project.boundaries_dir_namea + '.sww'
 
     Bf = Field_boundary(kwargs['boundary_file'],
-                    domain, time_thinning=time_thinning, mean_stage=tide, 
-                    use_cache=True, verbose=True)
+                        domain,
+                        time_thinning=time_thinning,
+                        mean_stage=tide, 
+                        use_cache=True,
+                        verbose=True)
                     
     kwargs['input_start_time']=domain.starttime
@@ -227 +234 @@
 if __name__ == "__main__":
 
-    run_model(file_name=project_urs.home+'detail.csv', aa_scenario_name=project_urs.scenario_name,
-              ab_time=project_urs.time, res_factor= project_urs.res_factor, tide=project_urs.tide, user=project_urs.user,
-              alpha = project_urs.alpha, friction=project_urs.friction,
-              time_thinning = project_urs.time_thinning,
-              dir_comment=project_urs.dir_comment)
-
-
+    run_model(file_name=project.home+'detail.csv', aa_scenario_name=project.scenario_name,
+              ab_time=project.time, res_factor= project.res_factor, tide=project.tide, user=project.user,
+              alpha = project.alpha, friction=project.friction,
+              time_thinning = project.time_thinning,
+              dir_comment=project.dir_comment)
+
+