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