""" Script for running a breaking wave simulation of Jon Hinwoods wave tank. Note: this is based on the frinction_ua_flume_2006 structure. Duncan Gray, GA - 2007 """ #---------------------------------------------------------------------------- # Import necessary modules #---------------------------------------------------------------------------- from Scientific.IO.NetCDF import NetCDFFile from Numeric import array, zeros, Float from os.path import join from anuga.utilities.numerical_tools import ensure_numeric from anuga.utilities.interp import interp import project # from os import getenv # from os.path import join # home = getenv('INUNDATIONHOME') # Hinwood_dir = join(home,'data','flumes','Hinwood_2008') # raw_data_dir = join(Hinwood_dir, 'raw_data') def csv2tms(filename, offshore_bed_elevation): """ Convert Hinwood boundary file to NetCDF tms file. the filename is the name of the output tms file, eg 'hi.tsm'. There must be an equivalent .csv file, eg 'hi.csv'. """ print 'Creating', filename # Read the ascii (.csv) version of this file fid = open(filename[:-4] + '.csv') # Read remaining lines lines = fid.readlines() fid.close() N = len(lines) T = zeros(N, Float) #Time Q = zeros(N, Float) #Stage X = zeros(N, Float) #XMomentum Y = zeros(N, Float) #YMomentum for i, line in enumerate(lines): fields = line.split(',') T[i] = float(fields[0]) Q[i] = float(fields[1]) depth = Q[i] - offshore_bed_elevation X[i] = float(fields[2]) * depth try: Y[i] = float(fields[3]) * depth except: pass # Create tms NetCDF file fid = NetCDFFile(filename, 'w') fid.institution = 'Geoscience Australia' fid.description = 'Input wave for Benchmark 2' fid.starttime = 0.0 fid.createDimension('number_of_timesteps', len(T)) fid.createVariable('time', Float, ('number_of_timesteps',)) fid.variables['time'][:] = T fid.createVariable('stage', Float, ('number_of_timesteps',)) fid.variables['stage'][:] = Q[:] fid.createVariable('xmomentum', Float, ('number_of_timesteps',)) fid.variables['xmomentum'][:] = X[:] fid.createVariable('ymomentum', Float, ('number_of_timesteps',)) fid.variables['ymomentum'][:] = Y[:] fid.close() def combine_velocity_depth(velocity_file, depth_file, out_file, metadata_dic): """ Convert the raw velocity and depth values, which have values at different times to a csv file, with values at the same time, with SI units. Set the depth values to be at the same times as the velocity values. The format for the velocity file is; [time, sec], [x-velocity +ve is towards the wave generator, cm/sec], [y-velocity], [z-velocity] The format for the pressure file is [time, sec], [mm above SWL for sensor A], many other sensors... The format of the output file is; [time, sec], [m above SWL for sensor A], [x-velocity +ve is towards the shore, m/sec], [y-velocity +ve is towards the left of the tank, looking from the generator to the shore, m/sec] """ missing = 1e+20 # Read velocity file #print "*********************" #print "velocity_file", velocity_file vfid = open(velocity_file) lines = vfid.readlines() vfid.close() n_velocity = len(lines) vtimes = zeros(n_velocity, Float) #Time x_velocities = zeros(n_velocity, Float) # y_velocities = zeros(n_velocity, Float) # for i, line in enumerate(lines): fields = line.split() #(',') vtimes[i] = float(fields[0]) x_velocities[i] = float(fields[1]) y_velocities[i] = float(fields[2]) # Read the depth file #print "depth_file", depth_file dfid = open(depth_file) lines = dfid.readlines() dfid.close() n_depth = len(lines) n_sensors = len(lines[0].split()) dtimes = zeros(n_depth, Float) #Time depths = zeros(n_depth, Float) # sensors = zeros((n_depth,n_sensors), Float) for i, line in enumerate(lines): fields = line.split() #(',') fields = [float(j) for j in fields] dtimes[i] = fields[0] depths[i] = fields[1] sensors[i] = fields #print "dtimes", dtimes #print "depths", depths #print "vtimes", vtimes depths_at_vtimes = interp( depths, dtimes, vtimes, missing=missing) depths_at_vtimes = ensure_numeric(depths_at_vtimes) #print "len(dtimes)", len(vtimes) #print "len(depths_at_vtimes)", len(depths_at_vtimes) #print "metadata_dic['scenario_id']", metadata_dic['scenario_id'] # for i in range(len(depths_at_vtimes)): # print "i", i # print "vtimes[i]", vtimes[i] # print "depths_at_vtimes[i]", depths_at_vtimes[i] #print "depths_at_vtimes", depths_at_vtimes depths_at_vtimes = depths_at_vtimes/1000.00 # convert from mm to m missing=missing/1000.00 # Do to missing what is done to depths_at_vtimes x_velocities = ensure_numeric(x_velocities) # Swap axis around convert cm/sec to m/sec x_velocities = x_velocities * -0.01 y_velocities = ensure_numeric(y_velocities) # Swap axis around convert cm/sec to m/sec y_velocities = y_velocities * -0.01 fid = open(out_file,'w') assert len(depths_at_vtimes) == len(vtimes) start_time = None #start_time = 0.0 #for vtime, depth_at_vtime, velocity in map(vtimes, depths_at_vtimes, # velocities): for i in xrange(len(vtimes)): if not depths_at_vtimes[i] == missing: # Make the times start at zero. if start_time is None: start_time = vtimes[i] final_time = vtimes[i]-start_time fid.write(str(final_time) \ + ',' + str(depths_at_vtimes[i]) \ + ',' + str(x_velocities[i]) \ + ',' + str(y_velocities[i])+'\n') fid.close() print "The start time for is", start_time if not start_time == metadata_dic['ANUGA_start_time']: raise ValueError, "The calc'ed and recorded start times are different" # Modify the sensor array to reflect the new start time sensors[:,0] -= start_time # Since there is a new time reference save the depth info using this # new reference. fid = open(depth_file[:-4] + '_exp_depth.csv','w') #print "depth_file", depth_file #print "start_time", start_time bed_elevation_list = metadata_dic['gauge_bed_elevation'] # +2 due to Time column and gauge A max_j = len(bed_elevation_list)+2 # Write a header fid.write('Time') for gauge_x in metadata_dic['gauge_x']: fid.write(',' + str(gauge_x)) fid.write('\n') for i in xrange(len(dtimes)): fid.write(str(sensors[i,0])) # Time # Don't write sensor A. # It is the boundary condition. for j, bed_elevation in map(None, xrange(2,max_j), bed_elevation_list): # depth, m gauge = sensors[i,j]/1000 - bed_elevation # Convert from mm to m fid.write(',' + str(gauge)) fid.write('\n') fid.close() # Since there is a new time reference save the stage info using this # new reference. fid = open(depth_file[:-4] + '_exp_stage.csv','w') # Write a header fid.write('Time') for gauge_x in metadata_dic['gauge_x']: fid.write(',' + str(gauge_x)) fid.write('\n') for i in xrange(len(dtimes)): fid.write(str(sensors[i,0])) # Time # Don't write sensor A. # It is the boundary condition. for j in xrange(2,max_j): # stage, m gauge = sensors[i,j]/1000 # Convert from mm to m fid.write(',' + str(gauge)) fid.write('\n') fid.close() return final_time def prepare_time_boundary(metadata_dic, raw_data_dir, output_dir): """ Use this if a project instance has already been created. """ scenario_id = metadata_dic['scenario_id'] velocity_file = join(raw_data_dir,scenario_id+'velfilt.txt') depth_file = join(raw_data_dir,scenario_id+'pressfilt.txt') out_file = join(output_dir, scenario_id+'_boundary.csv') #out_depth_file = join(output_dir, scenario_id+'_exp_depth.csv') final_time = combine_velocity_depth(velocity_file, depth_file, out_file, metadata_dic) #print "metadata_dic['xleft'][1]", metadata_dic['xleft'][1] if metadata_dic['xleft'][1] >= 0.0: # This should be a -ve value, since the still water level is the # z origin. print "Warning: The z origin seems incorrect." tsm_file = out_file[:-4] + '.tsm' csv2tms(tsm_file, metadata_dic['xleft'][1]) return final_time # Don't do this, since run-dam changes the metadata_dic['xleft'][1], # which is used by this function def prepare_time_boundary_for_scenarios_warning(scenarios, outputdir_tag): for run_data in scenarios: id = run_data['scenario_id'] outputdir_name = id + outputdir_tag pro_instance = project.Project(['data','flumes','Hinwood_2008'], outputdir_name=outputdir_name) prepare_time_boundary(run_data, pro_instance.raw_data_dir, pro_instance.boundarydir) #------------------------------------------------------------------- if __name__ == "__main__": import sys; sys.exit() from scenarios import scenarios outputdir_tag = "XXXX" prepare_time_boundary_for_scenarios(scenarios, outputdir_tag)