Changeset 7549
- Timestamp:
- Oct 16, 2009, 12:04:30 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
misc/tools/event_selection/wave_energy/wave_energy.py
r7547 r7549 1 1 """ 2 Program to integrate the momentum quantities from a boundary time series. 3 Used to develop different methods of estimating offshore tsunami hazard, 4 rather than just wave height 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. 5 7 6 Creator: Jonathan Griffin 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 7 12 Created: 12 October 2009 8 13 """ … … 13 18 import csv 14 19 import glob 15 ##import matplotlib16 ##matplotlib.use('Agg')17 ##import pylab18 ##from pylab import *19 20 import numpy 20 21 22 # Change path here!!! 21 23 filepath = r'/nas/gemd/georisk_models/inundation/data/new_south_wales/batemans_bay_tsunami_scenario_2009/anuga/boundaries' 22 24 index = 2872 … … 27 29 filenamelist = [] 28 30 31 32 # Find all relevant gauges files and put into list 29 33 for root, dirs, files in os.walk(filepath): 30 ## print 'root', root 31 ## print 'dirs', dirs 32 ## print 'files', files 34 33 35 for filenames in files: 34 36 if filenames == filebasename: … … 36 38 filepathlist.append(filepathname) 37 39 filenamelist.append(filenames) 40 41 # Initialise lists for storing maximum values 42 max_energy_list = [] 43 max_energy_inst_list = [] 44 counter = 0 38 45 39 max_energy_list = [] 40 max_momentum_list = [] 41 max_energy_inst_list = [] 42 43 counter = 0 46 ########################################################### 47 # Loop over gauges files 48 ########################################################### 44 49 for filename in filepathlist: 45 50 time = [] … … 57 62 y_momentum.append(float(line[3])) 58 63 # Calculate momentum norm 59 momentum.append(numpy.sqrt(numpy.power(float(line[2]),2)+numpy.power(float(line[3]),2))) 60 64 61 65 time_diff = time[1]-time[0] 62 66 63 67 sign = 0 64 max_momentum = 065 momentum_sum = 066 68 energy_sum = 0 67 69 max_energy = 0 68 70 max_energy_inst = 0 71 72 ############################################################### 73 # Loop through time series and add values (integrate over time) 74 # Sum set to zero when stage height changes sign 75 ############################################################### 76 for i in range(len(time)): 69 77 70 for i in range(len(time)): 71 72 # Kinetic energy 73 KE = (1./2.)*(numpy.power(x_momentum[i],2)+numpy.power(y_momentum[i],2)) 74 # Potential energy 75 PE = 9.81*(stage[i]) 78 ############################################################### 79 # Energy calculations - note that quantities from gauge file are 80 # already depth integrated. 81 ############################################################### 82 # Calculate velocities from 'momentum' terms (Note not true momentum, 83 # rather defined as u*h. So we devide by the depth to get u). 84 # roh = density = 1000 kg.m-3 85 x_vel = x_momentum[i]/(100+stage[i]) 86 y_vel = y_momentum[i]/(100+stage[i]) 87 # Kinetic energy KE = 1/2 roh h u^2 88 KE = (1./2.)*(numpy.power(x_vel,2)+numpy.power(y_vel,2))*(1000)*(100+stage[i]) 89 # Potential energy PE = roh g w (w = stage) 90 PE = 9.81*(stage[i])*1000 76 91 energy = KE+PE 77 92 93 # Maximum instantaneous energy density 78 94 if energy > max_energy_inst: 79 95 max_energy_inst = energy 96 97 ############################################################### 98 # If sign of wave changes (goes from trough to crest or vice-versa) 99 # multiply energy sum by timestep and reset sum value to zero 100 ############################################################### 80 101 81 ## print 'max_energy',max_energy82 ## print 'energy_sum', energy_sum83 102 if stage[i] > 0 and sign != 1: 84 103 sign = 1 85 temp_max_momentum = momentum_sum*time_diff86 104 temp_max_energy = energy_sum*time_diff 87 if temp_max_momentum > max_momentum:88 max_momentum = temp_max_momentum89 105 if temp_max_energy > max_energy: 90 106 max_energy = temp_max_energy 107 91 108 energy_sum = energy 92 momentum_sum = momentum[i]93 94 109 start_time = time[i] 95 110 96 111 elif stage[i] < 0 and sign != -1: 97 112 sign = -1 98 temp_max_momentum = momentum_sum*time_diff99 113 temp_max_energy = energy_sum*time_diff 100 if temp_max_momentum > max_momentum:101 max_momentum = temp_max_momentum102 114 if temp_max_energy > max_energy: 103 115 max_energy = temp_max_energy 116 104 117 energy_sum = energy 105 momentum_sum = momentum[i]106 # print momentum_sum107 118 start_time = time[i] 108 119 109 120 elif stage[i] == 0 and sign != 0: 110 121 sign = 0 111 momentum_sum = momentum[i]112 122 energy_sum = energy 113 114 123 start_time = time[i] 115 124 116 125 else: 117 momentum_sum += momentum[i]118 126 energy_sum += energy 119 120 ## print 'KE',KE121 ## print 'PE',PE122 127 128 # Add maximum values to list 123 129 max_energy_inst_list.append(max_energy_inst) 124 130 max_energy_list.append( max_energy) 125 max_momentum_list.append(max_momentum) 126 ## print 'max_momentum', max_momentum 127 ## print 'max_energy', max_energy 131 128 132 counter += 1 133 134 129 135 counter = 0 130 136 total_max_energy = 0 131 total_max_momentum = 0132 137 total_max_inst_energy = 0 133 138 139 # Print results for each file and calculate maximum value across all events 134 140 for filename in filepathlist: 135 141 print filename 136 print 'max energy', max_energy_list[counter] 137 print 'max instantatneous energy', max_energy_inst_list[counter] 138 print 'max momentum', max_momentum_list[counter] 142 print 'max crest-integrated energy:', max_energy_list[counter],'J.s/m^2' 143 print 'max instantatneous energy:', max_energy_inst_list[counter],'J/m^2' 139 144 140 145 if max_energy_list[counter] > total_max_energy: … … 144 149 total_max_inst_energy = max_energy_inst_list[counter] 145 150 max_energy_inst_index = counter 146 if max_momentum_list[counter] > total_max_momentum: 147 total_max_momentum = max_momentum_list[counter] 148 max_momentum_index = counter 151 149 152 counter+=1 150 153 151 print '\n File with maximum energy integrated over a wave crest is: ', filepathlist[max_energy_index] 152 print 'Energy is ', max_energy_list[max_energy_index], ' J.s/m^2' 153 print '\n File with maximum instantaneous energy is: ', filepathlist[max_energy_inst_index] 154 print 'Energy is ', max_energy_inst_list[max_energy_inst_index], ' J/m^2' 155 print '\n File with maximum momentum is: ', filepathlist[max_momentum_index] 156 print 'Momentum is ', max_energy_list[max_momentum_index], 'kg.m/s' 154 print '\n File with maximum crest-integrated energy is:', filepathlist[max_energy_index] 155 print 'Energy is', max_energy_list[max_energy_index], 'J.s/m^2' 156 print '\n File with maximum instantaneous energy is:', filepathlist[max_energy_inst_index] 157 print 'Energy is', max_energy_inst_list[max_energy_inst_index], 'J/m^2'
Note: See TracChangeset
for help on using the changeset viewer.