Changeset 2640

Timestamp:

Mar 31, 2006, 1:29:48 PM (17 years ago)

Author:

sexton

Message:

Update Sydney example post Benfield visit (30 March 2006) to mimic changes from island.py

File:

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

production/sydney_2006/run_sydney_smf.py

@@ -20 +20 @@
 
 # Related major packages
-from pyvolution.shallow_water import Domain, Reflective_boundary
+from pyvolution.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary
 from pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+#from pyvolution.data_manager_old import convert_dem_from_ascii2netcdf, dem2pts
 from pyvolution.combine_pts import combine_rectangular_points_files 
 from pyvolution.pmesh2domain import pmesh_to_domain_instance
+from pyvolution.quantity import Quantity
+from Numeric import allclose
 
 # Application specific imports
@@ -57 +60 @@
         verbose=True)
 
-
 # combining the coarse and fine data
 combine_rectangular_points_files(project.finedemname + '.pts',
@@ -64 +66 @@
 
 #from pmesh.create_mesh import create_mesh_from_regions
-# new interface
+#new interface
 from pmesh.mesh_interface import create_mesh_from_regions
 
-# original
-interior_res = 5000
-interior_regions = [[project.harbour_polygon_2, interior_res],
-                    [project.botanybay_polygon_2, interior_res]]
+# original issue to Benfield
+#interior_res = 5000
+#interior_regions = [[project.harbour_polygon_2, interior_res],
+#                    [project.botanybay_polygon_2, interior_res]]
+
+# used for finer mesh
+interior_res1 = 5000
+interior_res2 = 315
+interior_regions = [[project.newpoly1, interior_res1],
+                    [project.south1, interior_res1],
+                    [project.finepolymanly, interior_res2],
+                    [project.finepolyquay, interior_res2]]
+
+print 'number of interior regions', len(interior_regions)
 
 #FIXME: Fix caching of this one once the mesh_interface is ready
 from caching import cache
 
+# original + refined region
 _ = cache(create_mesh_from_regions,
           project.diffpolygonall,
@@ -86 +99 @@
           verbose = True)
 
-
-#-----------------------------------------------------------------------------
+#-------------------------------------------------------------------------------                                 
 # Setup computational domain
-#----------------------------------------------------------------------------- 
+#-------------------------------------------------------------------------------                                 
 
 domain = pmesh_to_domain_instance(meshname, Domain,
@@ -97 +109 @@
 print 'Number of triangles = ', len(domain)
 print 'The extent is ', domain.get_extent()
+print domain.statistics()
 
 domain.set_name(project.basename)
@@ -119 +132 @@
 
 
-#------------------------------------------------------------------------------                                 
+#-------------------------------------------------------------------------------                                 
 # Setup initial conditions
 #-------------------------------------------------------------------------------
 
-domain.set_quantity('stage', tsunami_source)
-domain.set_quantity('friction', 0.0) 
+# apply slump after 30 mins, initialise to water level of tide = 0
+domain.set_quantity('stage', 0.0)
+domain.set_quantity('friction', 0.03) 
 domain.set_quantity('elevation',
                     filename = project.combineddemname + '.pts',
@@ -138 +152 @@
 
 Br = Reflective_boundary(domain)
-domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
-                      'right2': Br, 'top': Br, 'left1': Br,
-                      'left2': Br, 'left3': Br} )
+Bd = Dirichlet_boundary([0, 0, 0])
+
+# original + refined regions
+#domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
+#                      'right2': Br, 'top': Br, 'left1': Br,
+#                      'left2': Br, 'left3': Br} )
+
+# enforce Dirichlet BC - from 30/03/06 Benfield visit
+domain.set_boundary( {'bottom': Bd, 'right1': Bd, 'right0': Bd,
+                      'right2': Bd, 'top': Bd, 'left1': Bd,
+                      'left2': Bd, 'left3': Bd} )
+
+# increasingly finer interior regions
+#domain.set_boundary( {'bottom': Br, 'right': Br, 'left': Br, 'top': Br} )
 
 
@@ -149 +174 @@
 import time
 t0 = time.time()
+thisfile = project.integraltimeseries+'.csv'
+fid = open(thisfile, 'w')
 
-for t in domain.evolve(yieldstep = 120, finaltime = 18000): 
+for t in domain.evolve(yieldstep = 120, finaltime = 18000):
     domain.write_time()
-    domain.write_boundary_statistics(tags = 'bottom')     
+    domain.write_boundary_statistics(tags = 'bottom')
+
+    # calculate integral
+    thisstagestep = domain.get_quantity('stage') 
+    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
+    fid.write(s)
+
+    # add slump after 30 mins
+    if allclose(t, 30*60):
+        slump = Quantity(domain)
+        slump.set_values(tsunami_source)
+        domain.set_quantity('stage', slump + thisstagestep)
+    
     
 print 'That took %.2f seconds' %(time.time()-t0)