Changeset 5607


Ignore:
Timestamp:
Aug 5, 2008, 1:17:37 PM (16 years ago)
Author:
kristy
Message:
 
Location:
anuga_work/production
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/production/busselton/project.py

    r5575 r5607  
    3535scenario = 'busselton_tsunami_scenario'
    3636
    37 tide = 0.6
     37tide = 0.0 #0.6
    3838
    3939alpha = 0.1
  • anuga_work/production/perth/get_gauges.py

    r5592 r5607  
    1010from Scientific.IO.NetCDF import NetCDFFile
    1111
    12 def get_sts_gauge_data(filename,polygon,verbose=False):
     12def get_sts_gauge_data(filename,verbose=False):
    1313    from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
    1414        compress,zeros,fabs,take
     
    1717    y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
    1818    points=transpose(asarray([x.tolist(),y.tolist()]))
    19     time=fid.variables['time'][:]
     19    time=fid.variables['time'][:]+fid.starttime
    2020    elevation=fid.variables['elevation'][:]
    21     indices=inside_polygon(points, polygon, closed=True, verbose=False)
    22     inside_points=take(points,indices,0)
    23     elevation=take(elevation,indices,0)
    24     if verbose:
    25         inside_points_array=ensure_numeric(inside_points)#for plotting only
    26         #plot(points[:,0],points[:,1],'o')
    27         #plot(inside_points_array[:,0][::5],inside_points_array[:,1][::5],'o')
    28         polygon=ensure_numeric(polygon)
    29         #plot(polygon[:,0],polygon[:,1])
    30         #show()
    31 
    32    
     21       
    3322    basename='sts_gauge'
    3423    quantity_names=['stage','xmomentum','ymomentum']
     
    3625    for i, name in enumerate(quantity_names):
    3726        quantities[name] = fid.variables[name][:]
    38 ##        if inside_points is not None:
    39 ##            quantities[name] = take(quantities[name],indices,1)
    40 ##        else:
    41 ##            msg = 'No gauges found inside polygon'
    42 ##            raise msg
    43 ##           
    44     for j in range(len(x)-1):
     27
     28    maxname = 'max_sts_stage.csv'
     29    fid_max = open(project.boundaries_dir+sep+maxname,'w')
     30    s = 'x, y, max_stage \n'
     31    fid_max.write(s)   
     32    for j in range(len(x)):
    4533        locx=int(x[j])
    4634        locy=int(y[j])
     35        stage = quantities['stage'][:,j]
     36        xmomentum = quantities['xmomentum'][:,j]
     37        ymomentum = quantities['ymomentum'][:,j]
     38
     39        s = '%.6f, %.6f, %.6f\n' %(x[j], y[j], max(stage))
     40        fid_max.write(s)
     41       
    4742        fid_sts = open(project.boundaries_dir+sep+basename+'_'+str(locx)+'_'+str(locy)+'.csv', 'w')
    4843        s = 'time, stage, xmomentum, ymomentum \n'
    4944        fid_sts.write(s)
    50         stage = quantities['stage'][:,j]
    51         xmomentum = quantities['xmomentum'][:,j]
    52         ymomentum = quantities['ymomentum'][:,j]
     45       
    5346        for k in range(len(time)-1):
    54             s = '%.2f, %.2f, %.2f, %.2f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k])
     47            s = '%.6f, %.6f, %.6f, %.6f\n' %(time[k], stage[k], xmomentum[k], ymomentum[k])
    5548            fid_sts.write(s)
    5649
     
    5952    fid.close()
    6053       
    61     return inside_points,quantities,elevation,time
     54    return quantities,elevation,time
    6255
    63 polygon=project.poly_all
    64 points,quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir,project.scenario_name),polygon,verbose=False)
    65 print quantities['stage'][0,:]
    66 print len(points), len(elevation), len(quantities['stage'][0,:])
    67 assert len(points)==len(elevation)==len(quantities['stage'][0,:])
     56quantities,elevation,time=get_sts_gauge_data(os.path.join(project.boundaries_dir,project.scenario_name),verbose=False)
     57
     58print len(elevation), len(quantities['stage'][0,:])
  • anuga_work/production/perth/get_timeseries.py

    r5592 r5607  
    1818#timestamp='20080724_121200_run_trial_0.6_polyline_alpha0.1_kvanputt'
    1919#timestamp='20080724_161830_run_final_0.6_polyline_alpha0.1_kvanputt'
    20 timestamp='20080725_173911_run_final_0.6_polyline_alpha0.1_kvanputt'
     20timestamp='20080804_044523_run_final_0.0_polyline_alpha0.1_kvanputt'
    2121
    2222filename=project.output_dir+timestamp+sep+project.scenario_name+'.sww'
  • anuga_work/production/perth/project.py

    r5592 r5607  
    3737scenario = 'perth_tsunami_scenario'
    3838
     39
    3940tide = 0.0 #0.6
    4041
    41 alpha = 0.15
     42alpha = 0.1
    4243friction=0.01
    4344starttime=0
     
    151152output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
    152153
     154
     155vertex_filename = output_run_time_dir + 'mesh_vertex.csv'
     156
    153157#gauges
    154158gauge_name = 'perth.csv'
  • anuga_work/production/perth/run_perth.py

    r5578 r5607  
    4141from anuga.fit_interpolate.benchmark_least_squares import mem_usage
    4242from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
     43from anuga.geospatial_data.geospatial_data import find_optimal_smoothing_parameter
    4344
    4445# Application specific imports
     
    113114                             verbose=True)
    114115   # barrier()
    115    
     116
     117        covariance_value,alpha = find_optimal_smoothing_parameter (data_file= kwargs['elevation_file'],
     118                                alpha_list=[0.001, 0.01, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5],
     119                                mesh_file = project.meshes_dir_name+'.msh')
     120        print 'optimal alpha', covariance_value,alpha       
     121
    116122
    117123    #-------------------------------------------------------------------------
     
    154160    #barrier()
    155161
     162##    #------------------------------------------------------
     163##    # Create x,y,z file of mesh vertex!!!
     164##    #------------------------------------------------------
     165##        coord = domain.get_vertex_coordinates()
     166##        depth = domain.get_quantity('elevation')
     167##       
     168##        # Write vertex coordinates to file
     169##        filename=project.vertex_filename
     170##        fid=open(filename,'w')
     171##        fid.write('x (m), y (m), z(m)\n')
     172##        for i in range(len(coord)):
     173##            pt=coord[i]
     174##            x=pt[0]
     175##            y=pt[1]
     176##            z=depth[i]
     177##            fid.write('%.6f,%.6f,%.6f\n' %(x, y, z))
     178##
    156179
    157180    #------------------------------------------------------
     
    186209   
    187210    print 'Available boundary tags', domain.get_boundary_tags()
    188     Bf = File_boundary(boundary_urs_out+'.sts',
    189                    domain, time_thinning=1,
     211    Bf = Field_boundary(boundary_urs_out+'.sts',  # Change from file_boundary
     212                   domain, mean_stage= project.tide,
     213                   time_thinning=1,
    190214                   use_cache=True,
    191                    verbose = True,
    192                   boundary_polygon=bounding_polygon)
     215                   verbose = True)
     216                  # boundary_polygon=bounding_polygon)
    193217
    194218    Br = Reflective_boundary(domain)
    195219    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
    196220
     221    print dir(Bf)
    197222    print 'start reading boundary file'
     223   
    198224
    199225##    Bf = Field_boundary(kwargs['boundary_file'],
     
    214240    t0 = time.time()
    215241
     242    using_sts_boundary = True
     243   
    216244    for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime']
    217245                       ,skip_initial_step = False):
    218246        domain.write_time()
    219         domain.write_boundary_statistics(tags = 'ocean')   
     247        domain.write_boundary_statistics(tags = 'ocean')
     248        if using_sts_boundary is True:
     249            if Bf.values >= 99:
     250                domain.set_boundary({'ocean': Bd})
     251                using_sts_boundary = False
    220252
    221253    x, y = domain.get_maximum_inundation_location()
Note: See TracChangeset for help on using the changeset viewer.