Changeset 2350 for production

Timestamp:

Feb 7, 2006, 4:56:20 PM (19 years ago)

Author:

sexton

Message:

tidying up scripts for Sydney tsunami scenario

Location:

production/sydney_2006

Files:

• project.py (modified) (3 diffs)
• run_sydney_smf.py (modified) (4 diffs)

production/sydney_2006/project.py

@@ -44 +44 @@
 #nminvix = 6283000
 
-basename = 'slump'
+basename = 'slump2'
 
 if sys.platform == 'win32':
-    home = '..\..\..\..'     #Sandpit's parent dir
+    home = '..\..\..\..\..'     #Sandpit's parent dir
 else:    
     home = expanduser('~')
@@ -69 +69 @@
 
 #Georeferencing
-from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees
+#from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees
+from coordinate_transforms.redfearn import degminsec2decimal_degrees
 
 #Origin of existing dem (FIXME: Temporary measure)
@@ -77 +78 @@
 x0_origin = 314036.58727982 #revised 100m and 25m data
 y0_origin = 6224951.2960092
-mesh_origin = (refzone, x0_origin, y0_origin)  # input from Neil's data
+#mesh_origin = (refzone, x0_origin, y0_origin)  # input from Neil's data
 
 # define clipping polygon

production/sydney_2006/run_sydney_smf.py

@@ -68 +68 @@
                                  project.combineddemname + '.pts')
                                  
-# this allows the user to switch between different clipping polygons
-# print 'Which total zone are you interested in?'
-# mytest = int(raw_input('0 = all, 1 = harbour and 2 = botany bay     '))
-
-mytest = 0
-
 # Create Triangular Mesh
 # Overall clipping polygon and interior regions defined in project.py
 from pmesh.create_mesh import create_mesh_from_regions
 
-if mytest == 0:
-    # for whole region
-    interior_res = 5000 
-    interior_regions = [[project.harbour_polygon_2, interior_res],
-                        [project.botanybay_polygon_2, interior_res]] # maximal area of per triangle
+# for whole region
+interior_res = 5000 # maximal area of per triangle
+interior_regions = [[project.harbour_polygon_2, interior_res],
+                    [project.botanybay_polygon_2, interior_res]] 
 
-    m = cache(create_mesh_from_regions,
-              project.diffpolygonall,
-              {'boundary_tags': {'bottom': [0], 
-                                 'right1': [1], 'right0': [2],
-                                 'right2': [3], 'top': [4], 'left1': [5],
-                                 'left2': [6], 'left3': [7]},
-               'resolution': 100000,
-               'filename': meshname,           
-               'interior_regions': interior_regions},
-              #evaluate=True,    
-              verbose = True)
-    #import sys; sys.exit()
-    
-    #Add elevation data to the mesh - this is in place of set_quantity
-    mesh_elevname = 'elevation.msh'   
-    cache(fit_to_mesh_file,(meshname,
-                     #project.finedemname + '.pts',
-                     #project.coarsedemname + '.pts',
-                     project.combineddemname + '.pts',
-                     mesh_elevname,
-                     DEFAULT_ALPHA),
-            {'verbose': True, 'expand_search': True, 'precrop': True},
-          dependencies = [meshname],
-          #evaluate = True,   
-          verbose = False)
-    meshname = mesh_elevname
-    
+m = cache(create_mesh_from_regions,
+            project.diffpolygonall,
+            {'boundary_tags': {'bottom': [0], 
+                            'right1': [1], 'right0': [2],
+                            'right2': [3], 'top': [4], 'left1': [5],
+                            'left2': [6], 'left3': [7]},
+            'resolution': 100000,
+            'filename': meshname,           
+            'interior_regions': interior_regions},
+            #evaluate=True,    
+            verbose = True)
+
+#Add elevation data to the mesh - this is in place of set_quantity
+#mesh_elevname = 'elevation.msh'   
+#cache(fit_to_mesh_file,(meshname, project.combineddemname + '.pts',
+#                        mesh_elevname, DEFAULT_ALPHA),
+#      {'verbose': True, 'expand_search': True, 'precrop': True},
+#      dependencies = [meshname],
+#      #evaluate = True,   
+#      verbose = False)
+#meshname = mesh_elevname
+
 # Setup domain
 domain = cache(pmesh_to_domain_instance, (meshname, Domain),
                dependencies = [meshname],                     
-               verbose = True)               
+               verbose = True)
+
+# This section replaced with fit_to_mesh_file above
+domain.set_quantity('elevation',
+                     filename = project.combineddemname + '.pts',
+#                     filename = project.finedemname + '.pts',
+                     use_cache = True,
+                     verbose = True)
+
+          
         
 print 'Number of triangles = ', len(domain)
@@ -122 +119 @@
 domain.set_datadir(project.outputdir)
 domain.store = True
-#domain.quantities_to_be_stored = ['stage']
 domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
 
-#print 'Do you want to run the slump scenario?'
-#q2 = int(raw_input('0 for yes, 1 for no     '))
-
-q2 = 0
-
-if q2 == 0:
-    # slump parameters
-    # as advised by Adrian, this can be as much as 15000-30000 metres
-    # len = 4500.0
-    len = 30000.0
-    #thk = 670 #thk = 760.0
-    dep = 400.0
-    thk = 0.44*dep #thk = 120.0
-    wid = 15000.0 # for slump scenario, wid=len
-    #dep = 400 #1000.0 #d = 400 doesn't make sense with thk = 670!
-    rad = 3330
-    dp = 0.23
-    th = 6
-    alpha = 0.0
-    x0 = project.x0
-    y0 = project.y0
-    # slide parameters
-    # len = 10000.0
-    # dep = 1200.0
-    # thk = 200.0
-    # th = 12
-    # Slump or Slide Initial conditions
-    t = slump_tsunami(length=len, depth=dep, slope=th, thickness=thk, \
-                      radius=rad, dphi=dp, x0=x0, y0=y0, alpha=alpha)
-    #t = slide_tsunami(length=len, depth=dep, slope=th, thickness=thk, \
-    #                  x0=x0, y0=y0, alpha=alpha)
-    domain.set_quantity('stage', t, verbose = True)
-
-if q2 == 1:
-    domain.set_quantity('stage', tide)
-
-
-# These quantities are the same regardless of the initial condition
+# slump parameters, wid=len
+len = 30000.0
+dep = 400.0
+#thk = 120.0 for scenario 1 and 176 for scenario 2 rang 0.2-0.44 of depth
+thk = 0.44*dep 
+rad = 3330
+dp = 0.23
+th = 6
+alpha = 0.0
+x0 = project.x0 # slump origin
+y0 = project.y0
+# Setup Initial conditions
+t = slump_tsunami(length=len, depth=dep, slope=th, thickness=thk, \
+                  radius=rad, dphi=dp, x0=x0, y0=y0, alpha=alpha)
+domain.set_quantity('stage', t)
 domain.set_quantity('friction', 0)
-
-    
+   
 # Setup Boundary Conditions
 print domain.get_boundary_tags()
@@ -176 +148 @@
                    f=lambda t: [(6<t<606)*6, 0, 0])
 
-# for slump scenario
-if q2 == 0:
-    domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
-                                 'right2': Br, 'top': Br, 'left1': Br,
-                                 'left2': Br, 'left3': Br} )
-
-# for square wave scenario
-if q2 == 1:
-    domain.set_boundary( {'top': Br, 'bottom': Br, 'right': Bw, 'left': Br} )
+domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
+                      'right2': Br, 'top': Br, 'left1': Br,
+                      'left2': Br, 'left3': Br} )
 
 # Evolve
@@ -190 +156 @@
 t0 = time.time()
 
-for t in domain.evolve(yieldstep = 120, finaltime = 18000): #120 10800 = 3hrs, 18000 = 5hrs
+for t in domain.evolve(yieldstep = 120, finaltime = 18000): 
     domain.write_time()
-    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage')     
+    domain.write_boundary_statistics(tags = 'bottom')     
     
 print 'That took %.2f seconds' %(time.time()-t0)
-
-#############################################################
-# testing system on a smaller mesh
-    #m = create_mesh_from_regions(project.harbour_polygon,
-    #                             boundary_tags={'bottom': [0], 'top': [2],
-    #                                            'right': [1], 'left': [3]},
-    #                             resolution=100000, filename=meshname)
-    
-# This section replaced with fit_to_mesh_file above
-#domain.set_quantity('elevation',
-#                     #filename = project.combineddemname + '.pts',
-#                     filename = project.finedemname + '.pts',
-                     #filename = project.datadir + 'bathy_land25_orig.pts',
-#                     use_cache = True,
-#                     verbose = True)
-
-
-# for switching between overall clipping regions - NOT USED
-
-if mytest == 1:
-    # for harbour region
-    south = project.hsouth 
-    north = project.hnorth
-    west = project.hwest
-    east = project.heast
-
-    #interior_regions = [[project.harbour_polygon, 25000]] # maximal area of per triangle
-    # meshname = project.meshname + 'harbour.msh'                   
-    m = cache(create_mesh_from_regions,
-              project.polygon_h,
-              {'boundary_tags': {'bottom': [0], 'top': [2],
-                                 'right': [1], 'left': [3]},
-               'resolution': 50000,
-               'filename': meshname},           
-     #          'interior_regions': interior_regions},
-              evaluate=True, 
-              verbose = True)
-
-if mytest == 2:
-    # for botany bay region
-    south = project.bsouth 
-    north = project.bnorth
-    west = project.bwest
-    east = project.beast
-
-    #interior_regions = [[project.botanybay_polygon, 25000]] # maximal area of per triangle
-    # meshname = project.meshname + 'bay.msh'                   
-    m = cache(create_mesh_from_regions,
-              project.polygon_bb,
-              {'boundary_tags': {'bottom': [0], 'top': [2],
-                                 'right': [1], 'left': [3]},
-               'resolution': 50000,
-               'filename': meshname},
-     #          'interior_regions': interior_regions},
-              evaluate=True, 
-              verbose = True)
-