Changeset 7593 for misc/tools/event_selection/wave_energy
- Timestamp:
- Dec 17, 2009, 9:36:34 AM (15 years ago)
- Location:
- misc/tools/event_selection/wave_energy
- Files:
-
- 1 added
- 2 edited
-
run_up.py (modified) (2 diffs)
-
wave_energy.py (modified) (2 diffs)
-
wave_energy_technical_paper.doc (added)
Legend:
- Unmodified
- Added
- Removed
-
misc/tools/event_selection/wave_energy/run_up.py
r7584 r7593 20 20 21 21 R_break = gravity*(power(slope,2.0))/(power(freq,2.0)) 22 run_up = np.min(run_up_madsen, R_break)22 run_up = min(run_up_madsen, R_break) 23 23 return run_up 24 24 25 #def surf_sim(depth, slope, height, wavelength, gravity = 9.81): 26 # Add Madsen's surf similarity here 25 def surf_sim(depth, slope, height, period, gravity = 9.81): 26 wavelength = period*sqrt(gravity*depth) 27 print 'wavelength', wavelength 28 surf_sim_param = slope/(sqrt(height/wavelength)) 29 print 'surf_sim_param', surf_sim_param 30 #print 'height',height 31 ## print 'depth', depth 32 a = 2*(power(pi,(0.75))) 33 HA = height/depth 34 #print HA 35 b = power(HA,-0.25) 36 c =power(surf_sim_param,(-0.5)) 37 ## print 'a', a 38 ## print 'b' ,b 39 ## print 'c', c 40 R = height*a*b*c 41 #R = height*2*(power(pi,(3/4)))*(power((height/depth),(-1/4)))*(power(surf_sim_param,(-1/2))) 42 print 'R',R 43 R_break = height*(1.0/pi)*power(surf_sim_param,2.0) 44 print 'R_break',R_break 45 run_up = min(R, R_break) 46 ## print 'run_up',run_up 47 return run_up 48 49 27 50 28 51 … … 43 66 44 67 import numpy as np 45 from numpy import tan, sin, cos, power,pi68 from numpy import sin,tan,cos,power,pi, sqrt 46 69 if __name__ == '__main__': 47 70 sys.exit(main()) -
misc/tools/event_selection/wave_energy/wave_energy.py
r7584 r7593 13 13 """ 14 14 15 16 def wave_energy(filepath, filebasename): 17 18 filepathlist = [] 19 filenamelist = [] 20 21 22 # Find all relevant gauges files and put into list 23 for root, dirs, files in os.walk(filepath): 24 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 39 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] 60 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 69 ############################################################### 70 # Loop through time series and add values (integrate over time) 71 # Sum set to zero when stage height changes sign 72 ############################################################### 73 for i in range(len(time)): 74 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 90 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 ############################################################### 99 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)): 145 146 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) 161 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) 166 167 168 counter += 1 169 170 171 172 173 counter = 0 174 total_max_energy = 0 175 total_max_inst_energy = 0 176 177 178 179 180 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 206 207 counter+=1 208 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' 213 214 215 ############################################################################# 216 217 ############################################################################# 15 218 import os 16 219 import sys … … 20 223 import numpy 21 224 import run_up 22 23 # Change path here!!! 24 filepath = r'/nas/gemd/georisk_models/inundation/data/new_south_wales/batemans_bay_tsunami_scenario_2009/anuga/boundaries' 25 index = 2872 26 filebasename = 'sts_gauge_' + str(index) +'.csv' 27 #filename = join(filepath, filebasename) 28 29 filepathlist = [] 30 filenamelist = [] 31 32 33 # Find all relevant gauges files and put into list 34 for root, dirs, files in os.walk(filepath): 35 36 for filenames in files: 37 if filenames == filebasename: 38 filepathname = join(root, filenames) 39 filepathlist.append(filepathname) 40 filenamelist.append(filenames) 41 42 # Initialise lists for storing maximum values 43 max_energy_list = [] 44 max_energy_inst_list = [] 45 run_up_list = [] 46 energy_run_up_list = [] 47 max_stage_list = [] 48 counter = 0 49 50 ########################################################### 51 # Loop over gauges files 52 ########################################################### 53 for filename in filepathlist: 54 time = [] 55 stage = [] 56 x_momentum = [] 57 y_momentum = [] 58 momentum = [] 59 print 'reading file ', filename 60 csvreader = csv.reader(open(filename)) 61 header = csvreader.next() 62 for line in csvreader: 63 time.append(float(line[0])) 64 stage.append(float(line[1])) 65 x_momentum.append(float(line[2])) 66 y_momentum.append(float(line[3])) 67 # Calculate momentum norm 68 69 time_diff = time[1]-time[0] 70 71 sign = 0 72 energy_sum = 0 73 max_energy = 0 74 max_energy_inst = 0 75 max_stage = 0 76 temp_max_stage = 0 225 if __name__ == '__main__': 77 226 78 ############################################################### 79 # Loop through time series and add values (integrate over time) 80 # Sum set to zero when stage height changes sign 81 ############################################################### 82 for i in range(len(time)): 83 84 ############################################################### 85 # Energy calculations - note that quantities from gauge file are 86 # already depth integrated. 87 ############################################################### 88 # Calculate velocities from 'momentum' terms (Note not true momentum, 89 # rather defined as u*h. So we devide by the depth to get u). 90 # roh = density = 1000 kg.m-3 91 x_vel = x_momentum[i]/(100+stage[i]) 92 y_vel = y_momentum[i]/(100+stage[i]) 93 velocity = numpy.sqrt(x_vel*x_vel+y_vel*y_vel) 94 # Kinetic energy KE = 1/2 roh h u^2 95 KE = (1./2.)*(numpy.power(x_vel,2)+numpy.power(y_vel,2))*(1000)*(100+stage[i]) 96 # Potential energy PE = roh g w (w = stage) 97 PE = 9.81*(stage[i])*1000 98 energy = KE+PE 99 100 # Maximum instantaneous energy density 101 if energy > max_energy_inst: 102 max_energy_inst = energy 103 104 ############################################################### 105 # If sign of wave changes (goes from trough to crest or vice-versa) 106 # multiply energy sum by timestep and reset sum value to zero 107 ############################################################### 108 109 if stage[i] > 0 and sign != 1: 110 sign = 1 111 temp_max_energy = energy_sum*time_diff 112 if temp_max_energy > max_energy: 113 max_energy = temp_max_energy 114 max_energy_time = time[i] - start_time 115 if temp_max_stage > max_stage: 116 max_stage = temp_max_stage 117 max_time = time[i] - start_time 118 energy_sum = energy 119 start_time = time[i] 120 max_stage_temp = 0 121 122 elif stage[i] < 0 and sign != -1: 123 sign = -1 124 temp_max_energy = energy_sum*time_diff 125 if temp_max_energy > max_energy: 126 max_energy = temp_max_energy 127 max_energy_time = time[i] - start_time 128 if temp_max_stage > max_stage: 129 max_stage = temp_max_stage 130 max_time = time[i] - start_time 131 energy_sum = energy 132 start_time = time[i] 133 134 elif stage[i] == 0 and sign != 0: 135 sign = 0 136 if temp_max_stage > max_stage: 137 max_stage = temp_max_stage 138 max_time = time[i] - start_time 139 energy_sum = energy 140 start_time = time[i] 141 142 else: 143 energy_sum += energy 144 temp_max_stage = max(temp_max_stage,stage[i]) 145 146 # Add maximum values to list 147 max_energy_inst_list.append(max_energy_inst) 148 max_energy_list.append(max_energy) 149 150 # get wave period from max time 151 ## locator = time.index(max_time) 152 ## time_count = 0 153 ## for i in range(locator,len(time)): 154 155 156 # Account for both halves of the wave 157 max_time = max_time*2 158 print 'max stage',max_stage 159 print 'max time',max_time 160 ##################################################################### 161 # call run-up module 162 ##################################################################### 163 Run_up = run_up.Madsen(100,0.017,max_stage,max_time) 164 print 'Madsen runup', Run_up 165 run_up_list.append(Run_up) 166 167 energy_run_up = run_up.energy_conserv(max_energy,1.0,0.8,9.81) 168 energy_run_up_list.append(energy_run_up) 169 print 'energy run up', energy_run_up 170 max_stage_list.append(max_stage) 171 172 173 counter += 1 174 175 176 177 178 counter = 0 179 total_max_energy = 0 180 total_max_inst_energy = 0 181 182 183 184 185 186 # Print results for each file and calculate maximum value across all events 187 outFilename= filepath +'/wave_energy_summary.csv' 188 outFile = file(outFilename,'w') 189 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') 190 for filename in filepathlist: 191 outFile.write(filename+',') 192 outFile.write(str(max_energy_list[counter])+',') 193 outFile.write(str(max_energy_inst_list[counter])+',') 194 outFile.write(str(run_up_list[counter])+',') 195 outFile.write(str(energy_run_up_list[counter])+',') 196 outFile.write(str(max_stage_list[counter])+'\n') 197 print filename 198 print 'max crest-integrated energy:', max_energy_list[counter],'J/m^2' 199 print 'max instantatneous energy:', max_energy_inst_list[counter],'J/s*m^2' 200 print 'Run_up (Madsen 2009):', run_up_list[counter], 'm' 201 print 'Run_up (Energy 2009):', energy_run_up_list[counter], 'm' 202 203 if max_energy_list[counter] > total_max_energy: 204 total_max_energy = max_energy_list[counter] 205 max_energy_index = counter 206 if max_energy_inst_list[counter] > total_max_inst_energy: 207 total_max_inst_energy = max_energy_inst_list[counter] 208 max_energy_inst_index = counter 209 210 counter+=1 211 212 print '\n File with maximum crest-integrated energy is:', filepathlist[max_energy_index] 213 print 'Energy is', max_energy_list[max_energy_index], 'J.s/m^2' 214 print '\n File with maximum instantaneous energy is:', filepathlist[max_energy_inst_index] 215 print 'Energy is', max_energy_inst_list[max_energy_inst_index], 'J/m^2' 216 217 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)
Note: See TracChangeset
for help on using the changeset viewer.
