Changeset 5751

Timestamp:

Sep 9, 2008, 2:14:12 PM (17 years ago)

Author:

kristy

Message:

 

Location:

anuga_work/production/geraldton

Files:

• build_geraldton.py (modified) (1 diff)
• get_gauges.py (modified) (2 diffs)
• project.py (modified) (3 diffs)
• run_geraldton.py (modified) (5 diffs)

anuga_work/production/geraldton/build_geraldton.py

@@ -91 +91 @@
 print'finished reading in files'
 
+G_off1_clip = G_off1.clip(project.poly_CBD_1km, verbose=True)
+
 print'add all geospatial objects'
-G = G1 + G2 + G3 + G_off + G_off1 + G_off2 + G_off3
+G = G1 + G2 + G3 + G_off + G_off1_clip + G_off2 + G_off3
 
 print'clip combined geospatial object by bounding polygon'

anuga_work/production/geraldton/get_gauges.py

@@ -14 +14 @@
         compress,zeros,fabs,take
     fid = NetCDFFile(filename+'.sts', 'r')    #Open existing file for read
+    permutation = fid.variables['permutation'][:]
     x = fid.variables['x'][:]+fid.xllcorner   #x-coordinates of vertices
     y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
@@ -31 +32 @@
     fid_max.write(s)    
     for j in range(len(x)):
-        locx=int(x[j])
-        locy=int(y[j])
+        index= permutation[j]
         stage = quantities['stage'][:,j]
         xmomentum = quantities['xmomentum'][:,j]
         ymomentum = quantities['ymomentum'][:,j]
 
-        s = '%.6f, %.6f, %.6f\n' %(x[j], y[j], max(stage))
+        s = '%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], max(stage))
         fid_max.write(s)
         
-        fid_sts = open(project.boundaries_dir+sep+basename+'_'+str(locx)+'_'+str(locy)+'.csv', 'w')
+        fid_sts = open(project.boundaries_dir+sep+basename+'_'+ str(index)+'.csv', 'w')
         s = 'time, stage, xmomentum, ymomentum \n'
         fid_sts.write(s)

anuga_work/production/geraldton/project.py

@@ -152 +152 @@
 
 poly_neg20_pos20 = read_polygon(polygons_dir+'neg20_pos20.csv')
-res_neg20_pos20 = 20000*res_factor
+res_neg20_pos20 = 25000*res_factor
 
 poly_CBD_1km= read_polygon(polygons_dir+'CBD_1km.csv')
@@ -161 +161 @@
 
 poly_wallabi = read_polygon(polygons_dir+'island_wallabi_poly.csv')
-res_wallabi = 1000*res_factor
+res_wallabi = 5000*res_factor
 
 poly_dingiville = read_polygon(polygons_dir+'island_dingiville_poly.csv')
-res_dingiville = 1000*res_factor
+res_dingiville = 5000*res_factor
 
 poly_pelsaert = read_polygon(polygons_dir+'island_pelsaert_poly.csv')
-res_pelsaert = 1000*res_factor
+res_pelsaert = 5000*res_factor
 
 
@@ -179 +179 @@
 print 'min number triangles', trigs_min
 
-#For no input boundary file
+##For no input boundary file
+#boundary_tags={'back': [2,3,4,5], 'side': [0,1,6], 'ocean': [7,8,9]}
 
-boundary_tags={'back': [2,3,4,5], 'side': [0,1,6], 'ocean': [7,8,9]}
 poly_mainland=read_polygon(topographies_in_dir +'initial_condition.csv')
 

anuga_work/production/geraldton/run_geraldton.py

@@ -78 +78 @@
     # 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_boundary.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+2,N+3], 'side': [N,N+1,N+4],'ocean': range(N)}
-
-    bounding_polygon = project.poly_all
-    
-
-    
+
+    # 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_boundary_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+2,N+3], 'side': [N,N+1,N+4],'ocean': range(N)}
+
+    print 'boundary tags',boundary_tags
+       
     #--------------------------------------------------------------------------
     # Create the triangular mesh based on overall clipping polygon with a
@@ -111 +109 @@
 
         create_mesh_from_regions(bounding_polygon,
-                             boundary_tags=project.boundary_tags,
+                             boundary_tags=boundary_tags,
                              maximum_triangle_area=project.res_poly_all,
                              interior_regions=project.interior_regions,
@@ -208 +206 @@
    
     boundary_urs_out=project.boundaries_dir_name
+
+    Br = Reflective_boundary(domain)
+    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
     
     print 'Available boundary tags', domain.get_boundary_tags()
@@ -213 +214 @@
                    domain, mean_stage= project.tide,
                    time_thinning=1,
+                   default_boundary=Bd,
                    use_cache=True,
                    verbose = True,
                    boundary_polygon=bounding_polygon)
 
-    
-    Br = Reflective_boundary(domain)
-    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
-
-##    Bf = Field_boundary(kwargs['boundary_file'],
-##                domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'], 
-##                use_cache=False, verbose=True)
-
     domain.set_boundary({'back': Br,
                          'side': Bd,
-                         'ocean': Bd}) # change baxk to Bf when running properly
+                         'ocean': Bf}) 
 
     kwargs['input_start_time']=domain.starttime
@@ -290 +284 @@
     run_model(**kwargs)
       
-    if myid==0:
-        export_model(**kwargs)
-    #barrier
+       #barrier