Changeset 2292 for production/sydney_2006

Timestamp:

Jan 26, 2006, 10:40:35 AM (19 years ago)

Author:

sexton

Message:

version 25 Jan - slump initial condition, new data

File:

• production/sydney_2006/run_sydney_smf.py (modified) (7 diffs)

production/sydney_2006/run_sydney_smf.py

@@ -21 +21 @@
 from pyvolution.data_manager import convert_dem_from_ascii2netcdf,\
      dem2pts
-from pyvolution.pmesh2domain import pmesh_to_domain_instance 
+from pyvolution.pmesh2domain import pmesh_to_domain_instance
+from pyvolution.combine_pts import combine_rectangular_points_files 
 from caching import cache
 import project
 from math import pi, sin
 from smf import slump_tsunami, slide_tsunami, Double_gaussian
-
-#Data preparation
-#Convert ASC 2 DEM 2 PTS using source data and store result in source data
-demname = project.demname
+from pyvolution.least_squares import fit_to_mesh_file, DEFAULT_ALPHA 
+
+# Data preparation
+# Convert ASC 2 DEM 2 PTS using source data and store result in source data
+# Do for coarse and fine data
+# Fine pts file to be clipped to area of interest
+coarsedemname = project.coarsedemname
+finedemname = project.finedemname
 meshname = project.meshname+'.msh'
 
-cache(convert_dem_from_ascii2netcdf, demname, {'verbose': True},
-      dependencies = [demname + '.asc'],
+# coarse data
+cache(convert_dem_from_ascii2netcdf, coarsedemname, {'verbose': True},
+      dependencies = [coarsedemname + '.asc'],
       verbose = True)
       #evaluate = True)
 
-cache(dem2pts, demname, {'verbose': True},
-      dependencies = [demname + '.dem'],      
-      verbose = True)
-
-#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     '))
+cache(dem2pts, coarsedemname, {'verbose': True},
+      dependencies = [coarsedemname + '.dem'],      
+      verbose = True)
+
+# fine data
+cache(convert_dem_from_ascii2netcdf, finedemname, {'verbose': True},
+      dependencies = [finedemname + '.asc'],
+      verbose = True)
+      #evaluate = True)
+
+# clipping the fine data
+cache(dem2pts, finedemname, {'verbose': True,
+      'easting_min': project.eastingmin,
+      'easting_max': project.eastingmax,
+      'northing_min': project.northingmin,
+      'northing_max': project.northingmax},
+      dependencies = [finedemname + '.dem'],
+      #evaluate = True,
+      verbose = True)
+
+# combining the coarse and fine data
+combine_rectangular_points_files(project.finedemname + '.pts',
+                                 project.coarsedemname + '.pts',
+                                 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
+# 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
-    #south = project.south 
-    #north = project.north
-    #west = project.west
-    #east = project.east
-
-    interior_regions = [[project.harbour_polygon_2, 25000],
-                        [project.botanybay_polygon_2, 25000]] # maximal area of per triangle
-
-    # meshname = project.meshname + 'all.msh'
+    interior_res = 2000 # run at 2000 for final one
+    interior_regions = [[project.harbour_polygon_2, interior_res],
+                        [project.botanybay_polygon_2, interior_res]] # maximal area of per triangle
 
     m = cache(create_mesh_from_regions,
-              #project.polygonall,
               project.diffpolygonall,
-              {'boundary_tags': {'bottom': [0], 'top': [4],
-                                 'right': [1], 'left1': [5],
-                                 'left2': [6], 'left3': [7],
+              {'boundary_tags': {'bottom': [0], 
                                  'right1': [1], 'right0': [2],
-                                 'right2': [3]},
+                                 'right2': [3], 'top': [4], 'left1': [5],
+                                 'left2': [6], 'left3': [7]},
                'resolution': 100000,
                'filename': meshname,           
@@ -76 +96 @@
               verbose = True)
     #import sys; sys.exit()
-
-#Setup domain
-
+    
+    #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
+    
+# Setup domain
 domain = cache(pmesh_to_domain_instance, (meshname, Domain),
                dependencies = [meshname],                     
@@ -86 +119 @@
 print 'The extent is ', domain.get_extent()
 
-#outputfilename = basename + '_slump'
-
 domain.set_name(project.basename)
-#domain.set_name(project.outfilename)
 domain.set_datadir(project.outputdir)
 domain.store = True
-
-#domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
-domain.quantities_to_be_stored = ['stage']
+#domain.quantities_to_be_stored = ['stage']
+domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
 
 #print 'Do you want to run the slump scenario?'
@@ -102 +131 @@
 
 if q2 == 0:
-    #slump parameters
-    #as advised by Adrian, this can be as much as 15000-30000 metres
+    # slump parameters
+    # as advised by Adrian, this can be as much as 15000-30000 metres
     # len = 4500.0
     len = 15000.0
     thk = 670 #thk = 760.0
-    #wid = 4500.0
     wid = 15000.0 # for slump scenario, wid=len
-    dep = 1000.0 #d = 1200
+    #dep = 400 #1000.0 #d = 400 doesn't make sense with thk = 670!
+    dep = 1000.0
     rad = 3330
     dp = 0.23
     th = 6
-
+    alpha = 0.0
+    x0 = project.x0
+    y0 = project.y0
     # slide parameters
     # len = 10000.0
@@ -119 +150 @@
     # thk = 200.0
     # th = 12
-    x0 = project.x0
-    y0 = project.y0
-    alpha = 0.0
-    ######################
-    # Slump Initial conditions
+    # 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)
+    domain.set_quantity('stage', t, verbose = True)
 
 if q2 == 1:
@@ -136 +163 @@
 # These quantities are the same regardless of the initial condition
 domain.set_quantity('friction', 0)
-domain.set_quantity('elevation',
-                     filename = demname + '.pts',
-                     use_cache = True,
-                     verbose = True)
+
     
 # Setup Boundary Conditions
@@ -151 +175 @@
                    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} )
 
-# for slump scenario
-if q2 == 0:
-    domain.set_boundary( {'top': Br, 'bottom': Br, 'right0': Br,
-                          'left1': Br, 'left2': Br, 'left3': Br,
-                          'right1': Br, 'right2': Br} )
-
 # Evolve
 import time
 t0 = time.time()
 
-for t in domain.evolve(yieldstep = 100, finaltime = 10400):
+for t in domain.evolve(yieldstep = 120, finaltime = 10800): #120 10800 = 3hrs
     domain.write_time()
     domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage')     
-#    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage')
     
 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