| [7547] | 1 | """ |
|---|
| [7549] | 2 | Program to calcualte the energy density in offshore tsunami time series. |
|---|
| 3 | Energy is integrated over the time of the each wave crest. |
|---|
| 4 | Also calculates time integrated momentum quantities from the boundary time series. |
|---|
| 5 | Used to investiage different methods of estimating offshore tsunami hazard, |
|---|
| 6 | rather than just wave height. |
|---|
| [7547] | 7 | |
|---|
| [7549] | 8 | Energy calculations from Johnson 1997. 'A modern introduction to the mathematical |
|---|
| 9 | theory of water waves' p23. |
|---|
| 10 | |
|---|
| 11 | Creator: Jonathan Griffin, Geoscience Australia |
|---|
| [7547] | 12 | Created: 12 October 2009 |
|---|
| 13 | """ |
|---|
| 14 | |
|---|
| 15 | |
|---|
| [7593] | 16 | def wave_energy(filepath, filebasename): |
|---|
| 17 | |
|---|
| 18 | filepathlist = [] |
|---|
| 19 | filenamelist = [] |
|---|
| [7547] | 20 | |
|---|
| 21 | |
|---|
| [7593] | 22 | # Find all relevant gauges files and put into list |
|---|
| 23 | for root, dirs, files in os.walk(filepath): |
|---|
| [7549] | 24 | |
|---|
| [7593] | 25 | for filenames in files: |
|---|
| 26 | if filenames == filebasename: |
|---|
| 27 | filepathname = join(root, filenames) |
|---|
| 28 | filepathlist.append(filepathname) |
|---|
| 29 | filenamelist.append(filenames) |
|---|
| 30 | |
|---|
| 31 | # Initialise lists for storing maximum values |
|---|
| 32 | max_energy_list = [] |
|---|
| 33 | max_energy_inst_list = [] |
|---|
| 34 | run_up_list = [] |
|---|
| 35 | run_up_list_surf_sim = [] |
|---|
| 36 | energy_run_up_list = [] |
|---|
| 37 | max_stage_list = [] |
|---|
| 38 | counter = 0 |
|---|
| [7549] | 39 | |
|---|
| [7593] | 40 | ########################################################### |
|---|
| 41 | # Loop over gauges files |
|---|
| 42 | ########################################################### |
|---|
| 43 | for filename in filepathlist: |
|---|
| 44 | time = [] |
|---|
| 45 | stage = [] |
|---|
| 46 | x_momentum = [] |
|---|
| 47 | y_momentum = [] |
|---|
| 48 | momentum = [] |
|---|
| 49 | print 'reading file ', filename |
|---|
| 50 | csvreader = csv.reader(open(filename)) |
|---|
| 51 | header = csvreader.next() |
|---|
| 52 | for line in csvreader: |
|---|
| 53 | time.append(float(line[0])) |
|---|
| 54 | stage.append(float(line[1])) |
|---|
| 55 | x_momentum.append(float(line[2])) |
|---|
| 56 | y_momentum.append(float(line[3])) |
|---|
| 57 | # Calculate momentum norm |
|---|
| 58 | |
|---|
| 59 | time_diff = time[1]-time[0] |
|---|
| [7547] | 60 | |
|---|
| [7593] | 61 | sign = 0 |
|---|
| 62 | energy_sum = 0 |
|---|
| 63 | max_energy = 0 |
|---|
| 64 | max_energy_inst = 0 |
|---|
| 65 | max_stage = 0 |
|---|
| 66 | temp_max_stage = 0 |
|---|
| 67 | max_dist = 0 |
|---|
| 68 | |
|---|
| [7549] | 69 | ############################################################### |
|---|
| [7593] | 70 | # Loop through time series and add values (integrate over time) |
|---|
| 71 | # Sum set to zero when stage height changes sign |
|---|
| [7549] | 72 | ############################################################### |
|---|
| [7593] | 73 | for i in range(len(time)): |
|---|
| [7547] | 74 | |
|---|
| [7593] | 75 | ############################################################### |
|---|
| 76 | # Energy calculations - note that quantities from gauge file are |
|---|
| 77 | # already depth integrated. |
|---|
| 78 | ############################################################### |
|---|
| 79 | # Calculate velocities from 'momentum' terms (Note not true momentum, |
|---|
| 80 | # rather defined as u*h. So we devide by the depth to get u). |
|---|
| 81 | # roh = density = 1000 kg.m-3 |
|---|
| 82 | x_vel = x_momentum[i]/(100+stage[i]) |
|---|
| 83 | y_vel = y_momentum[i]/(100+stage[i]) |
|---|
| 84 | velocity = numpy.sqrt(x_vel*x_vel+y_vel*y_vel) |
|---|
| 85 | # Kinetic energy KE = 1/2 roh h u^2 |
|---|
| 86 | KE = (1./2.)*(numpy.power(x_vel,2)+numpy.power(y_vel,2))*(1000)*(100+stage[i]) |
|---|
| 87 | # Potential energy PE = roh g w (w = stage) |
|---|
| 88 | PE = 9.81*(stage[i])*1000 |
|---|
| 89 | energy = KE+PE |
|---|
| [7549] | 90 | |
|---|
| [7593] | 91 | # Maximum instantaneous energy density |
|---|
| 92 | if energy > max_energy_inst: |
|---|
| 93 | max_energy_inst = energy |
|---|
| 94 | |
|---|
| 95 | ############################################################### |
|---|
| 96 | # If sign of wave changes (goes from trough to crest or vice-versa) |
|---|
| 97 | # multiply energy sum by timestep and reset sum value to zero |
|---|
| 98 | ############################################################### |
|---|
| [7547] | 99 | |
|---|
| [7593] | 100 | if stage[i] > 0 and sign != 1: |
|---|
| 101 | sign = 1 |
|---|
| 102 | temp_max_energy = energy_sum*time_diff |
|---|
| 103 | if temp_max_energy > max_energy: |
|---|
| 104 | max_energy = temp_max_energy |
|---|
| 105 | max_energy_time = time[i] - start_time |
|---|
| 106 | if temp_max_stage > max_stage: |
|---|
| 107 | max_stage = temp_max_stage |
|---|
| 108 | max_time = time[i] - start_time |
|---|
| 109 | energy_sum = energy |
|---|
| 110 | start_time = time[i] |
|---|
| 111 | #max_stage_temp = 0 |
|---|
| 112 | |
|---|
| 113 | elif stage[i] < 0 and sign != -1: |
|---|
| 114 | sign = -1 |
|---|
| 115 | temp_max_energy = energy_sum*time_diff |
|---|
| 116 | if temp_max_energy > max_energy: |
|---|
| 117 | max_energy = temp_max_energy |
|---|
| 118 | max_energy_time = time[i] - start_time |
|---|
| 119 | if temp_max_stage > max_stage: |
|---|
| 120 | max_stage = temp_max_stage |
|---|
| 121 | max_time = time[i] - start_time |
|---|
| 122 | energy_sum = energy |
|---|
| 123 | start_time = time[i] |
|---|
| 124 | |
|---|
| 125 | elif stage[i] == 0 and sign != 0: |
|---|
| 126 | sign = 0 |
|---|
| 127 | if temp_max_stage > max_stage: |
|---|
| 128 | max_stage = temp_max_stage |
|---|
| 129 | max_time = time[i] - start_time |
|---|
| 130 | energy_sum = energy |
|---|
| 131 | start_time = time[i] |
|---|
| 132 | |
|---|
| 133 | else: |
|---|
| 134 | energy_sum += energy |
|---|
| 135 | temp_max_stage = max(temp_max_stage,stage[i]) |
|---|
| 136 | |
|---|
| 137 | # Add maximum values to list |
|---|
| 138 | max_energy_inst_list.append(max_energy_inst) |
|---|
| 139 | max_energy_list.append(max_energy) |
|---|
| 140 | |
|---|
| 141 | # get wave period from max time |
|---|
| 142 | ## locator = time.index(max_time) |
|---|
| 143 | ## time_count = 0 |
|---|
| 144 | ## for i in range(locator,len(time)): |
|---|
| [7547] | 145 | |
|---|
| [7549] | 146 | |
|---|
| [7593] | 147 | # Account for both halves of the wave |
|---|
| 148 | max_time = max_time*2 |
|---|
| 149 | print 'max stage',max_stage |
|---|
| 150 | print 'max time',max_time |
|---|
| 151 | ##################################################################### |
|---|
| 152 | # call run-up module |
|---|
| 153 | ##################################################################### |
|---|
| 154 | Run_up = run_up.Madsen(100,0.017,max_stage,max_time) |
|---|
| 155 | print 'Madsen runup', Run_up |
|---|
| 156 | run_up_list.append(Run_up) |
|---|
| 157 | |
|---|
| 158 | Run_up_surf_sim = run_up.surf_sim(100,0.03,max_stage,max_time) |
|---|
| 159 | print 'Madsen surf similarity runup', Run_up_surf_sim |
|---|
| 160 | run_up_list_surf_sim.append(Run_up_surf_sim) |
|---|
| [7584] | 161 | |
|---|
| [7593] | 162 | energy_run_up = run_up.energy_conserv(max_energy,1.0,0.8,9.81) |
|---|
| 163 | energy_run_up_list.append(energy_run_up) |
|---|
| 164 | print 'energy run up', energy_run_up |
|---|
| 165 | max_stage_list.append(max_stage) |
|---|
| [7584] | 166 | |
|---|
| [7593] | 167 | |
|---|
| 168 | counter += 1 |
|---|
| 169 | |
|---|
| [7555] | 170 | |
|---|
| [7549] | 171 | |
|---|
| [7593] | 172 | |
|---|
| 173 | counter = 0 |
|---|
| 174 | total_max_energy = 0 |
|---|
| 175 | total_max_inst_energy = 0 |
|---|
| [7555] | 176 | |
|---|
| [7547] | 177 | |
|---|
| [7555] | 178 | |
|---|
| 179 | |
|---|
| 180 | |
|---|
| [7593] | 181 | # Print results for each file and calculate maximum value across all events |
|---|
| 182 | outFilename= filepath +'/wave_energy_summary.csv' |
|---|
| 183 | outFile = file(outFilename,'w') |
|---|
| 184 | outFile.write('filename, max crest-integrated energy (J.s/m^2), max instantatneous energy J/m^2, Run up (Madsen), Energy run up , max stage (m)\n') |
|---|
| 185 | for filename in filepathlist: |
|---|
| 186 | outFile.write(filename+',') |
|---|
| 187 | outFile.write(str(max_energy_list[counter])+',') |
|---|
| 188 | outFile.write(str(max_energy_inst_list[counter])+',') |
|---|
| 189 | outFile.write(str(run_up_list[counter])+',') |
|---|
| 190 | outFile.write(str(energy_run_up_list[counter])+',') |
|---|
| 191 | outFile.write(str(max_stage_list[counter])+',') |
|---|
| 192 | outFile.write(str(run_up_list_surf_sim[counter])+'\n') |
|---|
| 193 | print filename |
|---|
| 194 | print 'max crest-integrated energy:', max_energy_list[counter],'J/m^2' |
|---|
| 195 | print 'max instantatneous energy:', max_energy_inst_list[counter],'J/s*m^2' |
|---|
| 196 | print 'Run_up (Madsen 2009):', run_up_list[counter], 'm' |
|---|
| 197 | print 'Run_up surf similarity (Madsen 2009):', run_up_list_surf_sim[counter], 'm' |
|---|
| 198 | print 'Run_up (Energy 2009):', energy_run_up_list[counter], 'm' |
|---|
| 199 | |
|---|
| 200 | if max_energy_list[counter] > total_max_energy: |
|---|
| 201 | total_max_energy = max_energy_list[counter] |
|---|
| 202 | max_energy_index = counter |
|---|
| 203 | if max_energy_inst_list[counter] > total_max_inst_energy: |
|---|
| 204 | total_max_inst_energy = max_energy_inst_list[counter] |
|---|
| 205 | max_energy_inst_index = counter |
|---|
| [7555] | 206 | |
|---|
| [7593] | 207 | counter+=1 |
|---|
| [7549] | 208 | |
|---|
| [7593] | 209 | print '\n File with maximum crest-integrated energy is:', filepathlist[max_energy_index] |
|---|
| 210 | print 'Energy is', max_energy_list[max_energy_index], 'J.s/m^2' |
|---|
| 211 | print '\n File with maximum instantaneous energy is:', filepathlist[max_energy_inst_index] |
|---|
| 212 | print 'Energy is', max_energy_inst_list[max_energy_inst_index], 'J/m^2' |
|---|
| [7547] | 213 | |
|---|
| [7555] | 214 | |
|---|
| [7593] | 215 | ############################################################################# |
|---|
| 216 | |
|---|
| 217 | ############################################################################# |
|---|
| 218 | import os |
|---|
| 219 | import sys |
|---|
| 220 | from os.path import join |
|---|
| 221 | import csv |
|---|
| 222 | import glob |
|---|
| 223 | import numpy |
|---|
| 224 | import run_up |
|---|
| 225 | if __name__ == '__main__': |
|---|
| 226 | |
|---|
| 227 | # Change path here!!! |
|---|
| 228 | filepath = r'/nas/gemd/georisk_models/inundation/data/new_south_wales/batemans_bay_tsunami_scenario_2009/anuga/boundaries' |
|---|
| 229 | index = 2872 |
|---|
| 230 | filebasename = 'sts_gauge_' + str(index) +'.csv' |
|---|
| 231 | wave_energy(filepath, filebasename) |
|---|
