""" 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 anuga.utilities.numerical_tools import ensure_numeric from interp import interp def prepare_time_boundary(filename): """Convert benchmark 2 time series to NetCDF tms file. """ 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] = depth = float(fields[1]) 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): """ 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, m/sec], [y-velocity], [z-velocity] The format for the pressure file is [time, sec], [mm above SWL for sensor A], many other sensors... """ missing = 1e+20 # Read 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 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) # 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) x_velocities = x_velocities * -1.0 # Swap axis around y_velocities = ensure_numeric(y_velocities) y_velocities = y_velocities * -1.0 # Swap axis around 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] fid.write(str(vtimes[i]-start_time) \ + ',' + str(depths_at_vtimes[i]) \ + ',' + str(x_velocities[i]) \ + ',' + str(y_velocities[i])+'\n') fid.close() # Since there is a new time reference save the depth info using this # new reference. fid = open(depth_file[:-4] + '_new_time'+depth_file[-4:],'w') sensors[:,0] -= start_time #print "depth_file", depth_file #print "start_time", start_time for i in xrange(len(dtimes)): fid.write(str(sensors[i,0])) for j in xrange(1,len(sensors[0])): fid.write(' ' + str(sensors[i,j])) fid.write('\n') fid.close() #------------------------------------------------------------------- if __name__ == "__main__": 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') # Test 1 Run 5 combine_velocity_depth(join(raw_data_dir,'T1R5velfilt.txt'), join(raw_data_dir,'T1R5pressfilt.txt'), join(Hinwood_dir, 'T1R5_boundary.csv')) # Create the tsm file prepare_time_boundary(join(Hinwood_dir, 'T1R5_boundary.tsm')) # Test 2 Run 7 combine_velocity_depth(join(raw_data_dir,'T2R7velfilt.txt'), join(raw_data_dir,'T2R7pressfilt.txt'), join(Hinwood_dir, 'T2R7_boundary.csv')) # Create the tsm file prepare_time_boundary(join(Hinwood_dir, 'T2R7_boundary.tsm'))