[5076] | 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 |
---|
[5350] | 17 | |
---|
| 18 | from anuga.utilities.numerical_tools import ensure_numeric |
---|
[5076] | 19 | |
---|
[5350] | 20 | from interp import interp |
---|
[5076] | 21 | |
---|
[5392] | 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): |
---|
[5076] | 31 | """Convert benchmark 2 time series to NetCDF tms file. |
---|
[5392] | 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'. |
---|
[5370] | 34 | |
---|
[5076] | 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]) |
---|
[5392] | 59 | Q[i] = float(fields[1]) |
---|
| 60 | depth = Q[i] - offshore_bed_elevation |
---|
[5076] | 61 | X[i] = float(fields[2]) * depth |
---|
[5350] | 62 | try: |
---|
| 63 | Y[i] = float(fields[3]) * depth |
---|
| 64 | except: |
---|
| 65 | pass |
---|
[5076] | 66 | |
---|
| 67 | |
---|
[5370] | 68 | |
---|
[5076] | 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 | |
---|
[5350] | 89 | |
---|
[5395] | 90 | def combine_velocity_depth(velocity_file, depth_file, out_file, metadata_dic): |
---|
[5350] | 91 | """ |
---|
| 92 | |
---|
[5370] | 93 | Convert the raw velocity and depth values, which have values at |
---|
[5350] | 94 | different times to a csv file, with values at the same time, with |
---|
| 95 | SI units. |
---|
| 96 | |
---|
[5370] | 97 | Set the depth values to be at the same times as the velocity values. |
---|
| 98 | |
---|
[5350] | 99 | The format for the velocity file is; |
---|
[5395] | 100 | [time, sec], [x-velocity +ve is towards the wave generator, cm/sec], |
---|
[5350] | 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... |
---|
[5395] | 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 | |
---|
[5350] | 113 | """ |
---|
[5370] | 114 | missing = 1e+20 |
---|
| 115 | |
---|
[5350] | 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 |
---|
[5370] | 124 | x_velocities = zeros(n_velocity, Float) # |
---|
| 125 | y_velocities = zeros(n_velocity, Float) # |
---|
[5350] | 126 | |
---|
| 127 | for i, line in enumerate(lines): |
---|
[5370] | 128 | fields = line.split() #(',') |
---|
[5350] | 129 | |
---|
| 130 | vtimes[i] = float(fields[0]) |
---|
[5370] | 131 | x_velocities[i] = float(fields[1]) |
---|
| 132 | y_velocities[i] = float(fields[2]) |
---|
[5350] | 133 | |
---|
| 134 | # Read the depth file |
---|
| 135 | dfid = open(depth_file) |
---|
| 136 | lines = dfid.readlines() |
---|
| 137 | dfid.close() |
---|
| 138 | |
---|
[5370] | 139 | |
---|
[5350] | 140 | n_depth = len(lines) |
---|
[5370] | 141 | n_sensors = len(lines[0].split()) |
---|
[5350] | 142 | dtimes = zeros(n_depth, Float) #Time |
---|
| 143 | depths = zeros(n_depth, Float) # |
---|
[5370] | 144 | sensors = zeros((n_depth,n_sensors), Float) |
---|
[5350] | 145 | |
---|
| 146 | for i, line in enumerate(lines): |
---|
[5370] | 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) |
---|
[5350] | 158 | |
---|
[5370] | 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 |
---|
[5350] | 166 | depths_at_vtimes = depths_at_vtimes/1000.00 # convert from mm to m |
---|
[5370] | 167 | missing=missing/1000.00 # Do to missing what is done to depths_at_vtimes |
---|
| 168 | x_velocities = ensure_numeric(x_velocities) |
---|
[5395] | 169 | # Swap axis around convert cm/sec to m/sec |
---|
| 170 | x_velocities = x_velocities * -0.01 |
---|
[5370] | 171 | y_velocities = ensure_numeric(y_velocities) |
---|
[5395] | 172 | # Swap axis around convert cm/sec to m/sec |
---|
| 173 | y_velocities = y_velocities * -0.01 |
---|
[5350] | 174 | |
---|
| 175 | fid = open(out_file,'w') |
---|
| 176 | |
---|
| 177 | assert len(depths_at_vtimes) == len(vtimes) |
---|
[5370] | 178 | start_time = None |
---|
| 179 | #start_time = 0.0 |
---|
[5350] | 180 | #for vtime, depth_at_vtime, velocity in map(vtimes, depths_at_vtimes, |
---|
| 181 | # velocities): |
---|
| 182 | for i in xrange(len(vtimes)): |
---|
[5370] | 183 | if not depths_at_vtimes[i] == missing: |
---|
| 184 | # Make the times start at zero. |
---|
| 185 | if start_time is None: |
---|
[5395] | 186 | start_time = vtimes[i] |
---|
| 187 | final_time = vtimes[i]-start_time |
---|
| 188 | fid.write(str(final_time) \ |
---|
[5370] | 189 | + ',' + str(depths_at_vtimes[i]) \ |
---|
| 190 | + ',' + str(x_velocities[i]) \ |
---|
| 191 | + ',' + str(y_velocities[i])+'\n') |
---|
[5350] | 192 | |
---|
| 193 | fid.close() |
---|
| 194 | |
---|
[5370] | 195 | # Since there is a new time reference save the depth info using this |
---|
| 196 | # new reference. |
---|
[5395] | 197 | fid = open(depth_file[:-4] + '_exp_depth.csv','w') |
---|
[5370] | 198 | sensors[:,0] -= start_time |
---|
| 199 | #print "depth_file", depth_file |
---|
| 200 | #print "start_time", start_time |
---|
[5395] | 201 | |
---|
[5455] | 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 |
---|
[5395] | 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. |
---|
[5455] | 215 | for j, bed_elevation in map(None, |
---|
| 216 | xrange(2,max_j), |
---|
| 217 | bed_elevation_list): |
---|
[5395] | 218 | # depth, m |
---|
| 219 | gauge = sensors[i,j]/1000 - bed_elevation # Convert from mm to m |
---|
| 220 | fid.write(',' + str(gauge)) |
---|
[5370] | 221 | fid.write('\n') |
---|
| 222 | fid.close() |
---|
[5455] | 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 | |
---|
[5395] | 246 | return final_time |
---|
[5370] | 247 | |
---|
[5392] | 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') |
---|
[5395] | 255 | #out_depth_file = join(output_dir, scenario_id+'_exp_depth.csv') |
---|
[5350] | 256 | |
---|
[5395] | 257 | final_time = combine_velocity_depth(velocity_file, depth_file, out_file, |
---|
| 258 | metadata_dic) |
---|
[5392] | 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. |
---|
[5395] | 262 | print "Warning: The z origin seems incorrect." |
---|
| 263 | tsm_file = out_file[:-4] + '.tsm' |
---|
[5392] | 264 | csv2tms(tsm_file, metadata_dic['xleft'][1]) |
---|
[5395] | 265 | return final_time |
---|
[5370] | 266 | |
---|
[5350] | 267 | #------------------------------------------------------------------- |
---|
[5076] | 268 | if __name__ == "__main__": |
---|
[5392] | 269 | pass |
---|