Changeset 5607

Timestamp:

Aug 5, 2008, 1:17:37 PM (16 years ago)

Author:

kristy

Message:

 

Location:

anuga_work/production

Files:

• busselton/project.py (modified) (1 diff)
• perth/get_gauges.py (modified) (4 diffs)
• perth/get_timeseries.py (modified) (1 diff)
• perth/project.py (modified) (2 diffs)
• perth/run_perth.py (modified) (5 diffs)

anuga_work/production/busselton/project.py

@@ -35 +35 @@
 scenario = 'busselton_tsunami_scenario'
 
-tide = 0.6
+tide = 0.0 #0.6
 
 alpha = 0.1 

anuga_work/production/perth/get_gauges.py

@@ -10 +10 @@
 from Scientific.IO.NetCDF import NetCDFFile
 
-def get_sts_gauge_data(filename,polygon,verbose=False):
+def get_sts_gauge_data(filename,verbose=False):
     from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
         compress,zeros,fabs,take
@@ -17 +17 @@
     y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
     points=transpose(asarray([x.tolist(),y.tolist()]))
-    time=fid.variables['time'][:]
+    time=fid.variables['time'][:]+fid.starttime
     elevation=fid.variables['elevation'][:]
-    indices=inside_polygon(points, polygon, closed=True, verbose=False)
-    inside_points=take(points,indices,0)
-    elevation=take(elevation,indices,0)
-    if verbose:
-        inside_points_array=ensure_numeric(inside_points)#for plotting only
-        #plot(points[:,0],points[:,1],'o')
-        #plot(inside_points_array[:,0][::5],inside_points_array[:,1][::5],'o')
-        polygon=ensure_numeric(polygon)
-        #plot(polygon[:,0],polygon[:,1])
-        #show()
-
-    
+       
     basename='sts_gauge'
     quantity_names=['stage','xmomentum','ymomentum']
@@ -36 +25 @@
     for i, name in enumerate(quantity_names):
         quantities[name] = fid.variables[name][:]
-##        if inside_points is not None:
-##            quantities[name] = take(quantities[name],indices,1)
-##        else:
-##            msg = 'No gauges found inside polygon'
-##            raise msg
-##            
-    for j in range(len(x)-1):
+
+    maxname = 'max_sts_stage.csv'
+    fid_max = open(project.boundaries_dir+sep+maxname,'w')
+    s = 'x, y, max_stage \n'
+    fid_max.write(s)    
+    for j in range(len(x)):
         locx=int(x[j])
         locy=int(y[j])
+        stage = quantities['stage'][:,j]
+        xmomentum = quantities['xmomentum'][:,j]
+        ymomentum = quantities['ymomentum'][:,j]
+
+        s = '%.6f, %.6f, %.6f\n' %(x[j], y[j], max(stage))
+        fid_max.write(s)
+        
         fid_sts = open(project.boundaries_dir+sep+basename+'_'+str(locx)+'_'+str(locy)+'.csv', 'w')
         s = 'time, stage, xmomentum, ymomentum \n'
         fid_sts.write(s)
-        stage = quantities['stage'][:,j]
-        xmomentum = quantities['xmomentum'][:,j]
-        ymomentum = quantities['ymomentum'][:,j]
+        
         for k in range(len(time)-1):
-            s = '%.2f, %.2f, %.2f, %.2f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k])
+            s = '%.6f, %.6f, %.6f, %.6f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k])
             fid_sts.write(s)
 
@@ -59 +52 @@
     fid.close()
         
-    return inside_points,quantities,elevation,time
+    return quantities,elevation,time
 
-polygon=project.poly_all
-points,quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir,project.scenario_name),polygon,verbose=False)
-print quantities['stage'][0,:]
-print len(points), len(elevation), len(quantities['stage'][0,:])
-assert len(points)==len(elevation)==len(quantities['stage'][0,:])
+quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir,project.scenario_name),verbose=False)
+
+print len(elevation), len(quantities['stage'][0,:])

anuga_work/production/perth/get_timeseries.py

@@ -18 +18 @@
 #timestamp='20080724_121200_run_trial_0.6_polyline_alpha0.1_kvanputt'
 #timestamp='20080724_161830_run_final_0.6_polyline_alpha0.1_kvanputt'
-timestamp='20080725_173911_run_final_0.6_polyline_alpha0.1_kvanputt'
+timestamp='20080804_044523_run_final_0.0_polyline_alpha0.1_kvanputt'
 
 filename=project.output_dir+timestamp+sep+project.scenario_name+'.sww'

anuga_work/production/perth/project.py

@@ -37 +37 @@
 scenario = 'perth_tsunami_scenario'
 
+
 tide = 0.0 #0.6
 
-alpha = 0.15
+alpha = 0.1
 friction=0.01
 starttime=0
@@ -151 +152 @@
 output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
 
+
+vertex_filename = output_run_time_dir + 'mesh_vertex.csv'
+
 #gauges
 gauge_name = 'perth.csv'

anuga_work/production/perth/run_perth.py

@@ -41 +41 @@
 from anuga.fit_interpolate.benchmark_least_squares import mem_usage 
 from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
+from anuga.geospatial_data.geospatial_data import find_optimal_smoothing_parameter
 
 # Application specific imports
@@ -113 +114 @@
                              verbose=True)
    # barrier()
-    
+
+        covariance_value,alpha = find_optimal_smoothing_parameter (data_file= kwargs['elevation_file'],
+                                alpha_list=[0.001, 0.01, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5],
+                                mesh_file = project.meshes_dir_name+'.msh')
+        print 'optimal alpha', covariance_value,alpha       
+
 
     #-------------------------------------------------------------------------
@@ -154 +160 @@
     #barrier()
 
+##    #------------------------------------------------------
+##    # Create x,y,z file of mesh vertex!!!
+##    #------------------------------------------------------
+##        coord = domain.get_vertex_coordinates()
+##        depth = domain.get_quantity('elevation')
+##        
+##        # Write vertex coordinates to file
+##        filename=project.vertex_filename
+##        fid=open(filename,'w')
+##        fid.write('x (m), y (m), z(m)\n')
+##        for i in range(len(coord)):
+##            pt=coord[i]
+##            x=pt[0]
+##            y=pt[1]
+##            z=depth[i]
+##            fid.write('%.6f,%.6f,%.6f\n' %(x, y, z))
+##
 
     #------------------------------------------------------
@@ -186 +209 @@
     
     print 'Available boundary tags', domain.get_boundary_tags()
-    Bf = File_boundary(boundary_urs_out+'.sts', 
-                   domain, time_thinning=1,
+    Bf = Field_boundary(boundary_urs_out+'.sts',  # Change from file_boundary
+                   domain, mean_stage= project.tide,
+                   time_thinning=1,
                    use_cache=True,
-                   verbose = True,
-                   boundary_polygon=bounding_polygon)
+                   verbose = True)
+                  # boundary_polygon=bounding_polygon)
 
     Br = Reflective_boundary(domain)
     Bd = Dirichlet_boundary([kwargs['tide'],0,0])
 
+    print dir(Bf)
     print 'start reading boundary file'
+    
 
 ##    Bf = Field_boundary(kwargs['boundary_file'],
@@ -214 +240 @@
     t0 = time.time()
 
+    using_sts_boundary = True 
+    
     for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime']
                        ,skip_initial_step = False): 
         domain.write_time()
-        domain.write_boundary_statistics(tags = 'ocean')   
+        domain.write_boundary_statistics(tags = 'ocean')
+        if using_sts_boundary is True:
+            if Bf.values >= 99:
+                domain.set_boundary({'ocean': Bd})
+                using_sts_boundary = False
 
     x, y = domain.get_maximum_inundation_location()