| 1 | """ |
|---|
| 2 | Script for running a breaking wave simulation of Jon Hinwoods wave tank. |
|---|
| 3 | Note: this is based on the frinction_ua_flume_2006 structure. |
|---|
| 4 | |
|---|
| 5 | |
|---|
| 6 | Duncan Gray, GA - 2007 |
|---|
| 7 | |
|---|
| 8 | """ |
|---|
| 9 | |
|---|
| 10 | |
|---|
| 11 | #---------------------------------------------------------------------------- |
|---|
| 12 | # Import necessary modules |
|---|
| 13 | #---------------------------------------------------------------------------- |
|---|
| 14 | |
|---|
| 15 | from Scientific.IO.NetCDF import NetCDFFile |
|---|
| 16 | from Numeric import array, zeros, Float |
|---|
| 17 | |
|---|
| 18 | from anuga.utilities.numerical_tools import ensure_numeric |
|---|
| 19 | |
|---|
| 20 | from interp import interp |
|---|
| 21 | |
|---|
| 22 | from os.path import join |
|---|
| 23 | |
|---|
| 24 | # from os import getenv |
|---|
| 25 | # from os.path import join |
|---|
| 26 | # home = getenv('INUNDATIONHOME') |
|---|
| 27 | # Hinwood_dir = join(home,'data','flumes','Hinwood_2008') |
|---|
| 28 | # raw_data_dir = join(Hinwood_dir, 'raw_data') |
|---|
| 29 | |
|---|
| 30 | def csv2tms(filename, offshore_bed_elevation): |
|---|
| 31 | """Convert benchmark 2 time series to NetCDF tms file. |
|---|
| 32 | the filename is the name of the output tms file, eg 'hi.tsm'. |
|---|
| 33 | There must be an equivalent .csv file, eg 'hi.csv'. |
|---|
| 34 | |
|---|
| 35 | """ |
|---|
| 36 | |
|---|
| 37 | |
|---|
| 38 | |
|---|
| 39 | print 'Creating', filename |
|---|
| 40 | |
|---|
| 41 | # Read the ascii (.csv) version of this file |
|---|
| 42 | fid = open(filename[:-4] + '.csv') |
|---|
| 43 | |
|---|
| 44 | # Read remaining lines |
|---|
| 45 | lines = fid.readlines() |
|---|
| 46 | fid.close() |
|---|
| 47 | |
|---|
| 48 | |
|---|
| 49 | N = len(lines) |
|---|
| 50 | T = zeros(N, Float) #Time |
|---|
| 51 | Q = zeros(N, Float) #Stage |
|---|
| 52 | X = zeros(N, Float) #XMomentum |
|---|
| 53 | Y = zeros(N, Float) #YMomentum |
|---|
| 54 | |
|---|
| 55 | for i, line in enumerate(lines): |
|---|
| 56 | fields = line.split(',') |
|---|
| 57 | |
|---|
| 58 | T[i] = float(fields[0]) |
|---|
| 59 | Q[i] = float(fields[1]) |
|---|
| 60 | depth = Q[i] - offshore_bed_elevation |
|---|
| 61 | X[i] = float(fields[2]) * depth |
|---|
| 62 | try: |
|---|
| 63 | Y[i] = float(fields[3]) * depth |
|---|
| 64 | except: |
|---|
| 65 | pass |
|---|
| 66 | |
|---|
| 67 | |
|---|
| 68 | |
|---|
| 69 | # Create tms NetCDF file |
|---|
| 70 | fid = NetCDFFile(filename, 'w') |
|---|
| 71 | fid.institution = 'Geoscience Australia' |
|---|
| 72 | fid.description = 'Input wave for Benchmark 2' |
|---|
| 73 | fid.starttime = 0.0 |
|---|
| 74 | fid.createDimension('number_of_timesteps', len(T)) |
|---|
| 75 | fid.createVariable('time', Float, ('number_of_timesteps',)) |
|---|
| 76 | fid.variables['time'][:] = T |
|---|
| 77 | |
|---|
| 78 | fid.createVariable('stage', Float, ('number_of_timesteps',)) |
|---|
| 79 | fid.variables['stage'][:] = Q[:] |
|---|
| 80 | |
|---|
| 81 | fid.createVariable('xmomentum', Float, ('number_of_timesteps',)) |
|---|
| 82 | fid.variables['xmomentum'][:] = X[:] |
|---|
| 83 | |
|---|
| 84 | fid.createVariable('ymomentum', Float, ('number_of_timesteps',)) |
|---|
| 85 | fid.variables['ymomentum'][:] = Y[:] |
|---|
| 86 | |
|---|
| 87 | fid.close() |
|---|
| 88 | |
|---|
| 89 | |
|---|
| 90 | def combine_velocity_depth(velocity_file, depth_file, out_file, metadata_dic): |
|---|
| 91 | """ |
|---|
| 92 | |
|---|
| 93 | Convert the raw velocity and depth values, which have values at |
|---|
| 94 | different times to a csv file, with values at the same time, with |
|---|
| 95 | SI units. |
|---|
| 96 | |
|---|
| 97 | Set the depth values to be at the same times as the velocity values. |
|---|
| 98 | |
|---|
| 99 | The format for the velocity file is; |
|---|
| 100 | [time, sec], [x-velocity +ve is towards the wave generator, cm/sec], |
|---|
| 101 | [y-velocity], [z-velocity] |
|---|
| 102 | |
|---|
| 103 | The format for the pressure file is |
|---|
| 104 | [time, sec], [mm above SWL for sensor A], many other sensors... |
|---|
| 105 | |
|---|
| 106 | The format of the output file is; |
|---|
| 107 | [time, sec], [m above SWL for sensor A], |
|---|
| 108 | [x-velocity +ve is towards the shore, m/sec], |
|---|
| 109 | [y-velocity +ve is towards the left of the tank, looking from the |
|---|
| 110 | generator to the shore, m/sec] |
|---|
| 111 | |
|---|
| 112 | |
|---|
| 113 | """ |
|---|
| 114 | missing = 1e+20 |
|---|
| 115 | |
|---|
| 116 | # Read velocity file |
|---|
| 117 | vfid = open(velocity_file) |
|---|
| 118 | lines = vfid.readlines() |
|---|
| 119 | vfid.close() |
|---|
| 120 | |
|---|
| 121 | |
|---|
| 122 | n_velocity = len(lines) |
|---|
| 123 | vtimes = zeros(n_velocity, Float) #Time |
|---|
| 124 | x_velocities = zeros(n_velocity, Float) # |
|---|
| 125 | y_velocities = zeros(n_velocity, Float) # |
|---|
| 126 | |
|---|
| 127 | for i, line in enumerate(lines): |
|---|
| 128 | fields = line.split() #(',') |
|---|
| 129 | |
|---|
| 130 | vtimes[i] = float(fields[0]) |
|---|
| 131 | x_velocities[i] = float(fields[1]) |
|---|
| 132 | y_velocities[i] = float(fields[2]) |
|---|
| 133 | |
|---|
| 134 | # Read the depth file |
|---|
| 135 | dfid = open(depth_file) |
|---|
| 136 | lines = dfid.readlines() |
|---|
| 137 | dfid.close() |
|---|
| 138 | |
|---|
| 139 | |
|---|
| 140 | n_depth = len(lines) |
|---|
| 141 | n_sensors = len(lines[0].split()) |
|---|
| 142 | dtimes = zeros(n_depth, Float) #Time |
|---|
| 143 | depths = zeros(n_depth, Float) # |
|---|
| 144 | sensors = zeros((n_depth,n_sensors), Float) |
|---|
| 145 | |
|---|
| 146 | for i, line in enumerate(lines): |
|---|
| 147 | fields = line.split() #(',') |
|---|
| 148 | fields = [float(j) for j in fields] |
|---|
| 149 | dtimes[i] = fields[0] |
|---|
| 150 | depths[i] = fields[1] |
|---|
| 151 | sensors[i] = fields |
|---|
| 152 | #print "dtimes", dtimes |
|---|
| 153 | #print "depths", depths |
|---|
| 154 | #print "vtimes", vtimes |
|---|
| 155 | depths_at_vtimes = interp( depths, dtimes, |
|---|
| 156 | vtimes, missing=missing) |
|---|
| 157 | depths_at_vtimes = ensure_numeric(depths_at_vtimes) |
|---|
| 158 | |
|---|
| 159 | #print "len(dtimes)", len(vtimes) |
|---|
| 160 | #print "len(depths_at_vtimes)", len(depths_at_vtimes) |
|---|
| 161 | # for i in range(len(depths_at_vtimes)): |
|---|
| 162 | # print "i", i |
|---|
| 163 | # print "vtimes[i]", vtimes[i] |
|---|
| 164 | # print "depths_at_vtimes[i]", depths_at_vtimes[i] |
|---|
| 165 | # print "depths_at_vtimes", depths_at_vtimes |
|---|
| 166 | depths_at_vtimes = depths_at_vtimes/1000.00 # convert from mm to m |
|---|
| 167 | missing=missing/1000.00 # Do to missing what is done to depths_at_vtimes |
|---|
| 168 | x_velocities = ensure_numeric(x_velocities) |
|---|
| 169 | # Swap axis around convert cm/sec to m/sec |
|---|
| 170 | x_velocities = x_velocities * -0.01 |
|---|
| 171 | y_velocities = ensure_numeric(y_velocities) |
|---|
| 172 | # Swap axis around convert cm/sec to m/sec |
|---|
| 173 | y_velocities = y_velocities * -0.01 |
|---|
| 174 | |
|---|
| 175 | fid = open(out_file,'w') |
|---|
| 176 | |
|---|
| 177 | assert len(depths_at_vtimes) == len(vtimes) |
|---|
| 178 | start_time = None |
|---|
| 179 | #start_time = 0.0 |
|---|
| 180 | #for vtime, depth_at_vtime, velocity in map(vtimes, depths_at_vtimes, |
|---|
| 181 | # velocities): |
|---|
| 182 | for i in xrange(len(vtimes)): |
|---|
| 183 | if not depths_at_vtimes[i] == missing: |
|---|
| 184 | # Make the times start at zero. |
|---|
| 185 | if start_time is None: |
|---|
| 186 | start_time = vtimes[i] |
|---|
| 187 | final_time = vtimes[i]-start_time |
|---|
| 188 | fid.write(str(final_time) \ |
|---|
| 189 | + ',' + str(depths_at_vtimes[i]) \ |
|---|
| 190 | + ',' + str(x_velocities[i]) \ |
|---|
| 191 | + ',' + str(y_velocities[i])+'\n') |
|---|
| 192 | |
|---|
| 193 | fid.close() |
|---|
| 194 | |
|---|
| 195 | # Since there is a new time reference save the depth info using this |
|---|
| 196 | # new reference. |
|---|
| 197 | fid = open(depth_file[:-4] + '_exp_depth.csv','w') |
|---|
| 198 | sensors[:,0] -= start_time |
|---|
| 199 | #print "depth_file", depth_file |
|---|
| 200 | #print "start_time", start_time |
|---|
| 201 | |
|---|
| 202 | |
|---|
| 203 | bed_elevation_list = metadata_dic['gauge_bed_elevation'] |
|---|
| 204 | # +2 due to Time column and gauge A |
|---|
| 205 | max_j = len(bed_elevation_list)+2 |
|---|
| 206 | # Write a header |
|---|
| 207 | fid.write('Time') |
|---|
| 208 | for gauge_x in metadata_dic['gauge_x']: |
|---|
| 209 | fid.write(',' + str(gauge_x)) |
|---|
| 210 | fid.write('\n') |
|---|
| 211 | for i in xrange(len(dtimes)): |
|---|
| 212 | fid.write(str(sensors[i,0])) # Time |
|---|
| 213 | # Don't write sensor A. |
|---|
| 214 | # It is the boundary condition. |
|---|
| 215 | for j, bed_elevation in map(None, |
|---|
| 216 | xrange(2,max_j), |
|---|
| 217 | bed_elevation_list): |
|---|
| 218 | # depth, m |
|---|
| 219 | gauge = sensors[i,j]/1000 - bed_elevation # Convert from mm to m |
|---|
| 220 | fid.write(',' + str(gauge)) |
|---|
| 221 | fid.write('\n') |
|---|
| 222 | fid.close() |
|---|
| 223 | |
|---|
| 224 | |
|---|
| 225 | # Since there is a new time reference save the stage info using this |
|---|
| 226 | # new reference. |
|---|
| 227 | fid = open(depth_file[:-4] + '_exp_stage.csv','w') |
|---|
| 228 | sensors[:,0] -= start_time |
|---|
| 229 | |
|---|
| 230 | # Write a header |
|---|
| 231 | fid.write('Time') |
|---|
| 232 | for gauge_x in metadata_dic['gauge_x']: |
|---|
| 233 | fid.write(',' + str(gauge_x)) |
|---|
| 234 | fid.write('\n') |
|---|
| 235 | for i in xrange(len(dtimes)): |
|---|
| 236 | fid.write(str(sensors[i,0])) # Time |
|---|
| 237 | # Don't write sensor A. |
|---|
| 238 | # It is the boundary condition. |
|---|
| 239 | for j in xrange(2,max_j): |
|---|
| 240 | # stage, m |
|---|
| 241 | gauge = sensors[i,j]/1000 # Convert from mm to m |
|---|
| 242 | fid.write(',' + str(gauge)) |
|---|
| 243 | fid.write('\n') |
|---|
| 244 | fid.close() |
|---|
| 245 | |
|---|
| 246 | return final_time |
|---|
| 247 | |
|---|
| 248 | def prepare_time_boundary(metadata_dic, raw_data_dir, output_dir): |
|---|
| 249 | """ |
|---|
| 250 | """ |
|---|
| 251 | scenario_id = metadata_dic['scenario_id'] |
|---|
| 252 | velocity_file = join(raw_data_dir,scenario_id+'velfilt.txt') |
|---|
| 253 | depth_file = join(raw_data_dir,scenario_id+'pressfilt.txt') |
|---|
| 254 | out_file = join(output_dir, scenario_id+'_boundary.csv') |
|---|
| 255 | #out_depth_file = join(output_dir, scenario_id+'_exp_depth.csv') |
|---|
| 256 | |
|---|
| 257 | final_time = combine_velocity_depth(velocity_file, depth_file, out_file, |
|---|
| 258 | metadata_dic) |
|---|
| 259 | if metadata_dic['xleft'][1] >= 0.0: |
|---|
| 260 | # This should be a -ve value, since the still water level is the |
|---|
| 261 | # z origin. |
|---|
| 262 | print "Warning: The z origin seems incorrect." |
|---|
| 263 | tsm_file = out_file[:-4] + '.tsm' |
|---|
| 264 | csv2tms(tsm_file, metadata_dic['xleft'][1]) |
|---|
| 265 | return final_time |
|---|
| 266 | |
|---|
| 267 | #------------------------------------------------------------------- |
|---|
| 268 | if __name__ == "__main__": |
|---|
| 269 | pass |
|---|
