source: misc/tools/event_selection/wave_energy/wave_energy.py @ 7593

Last change on this file since 7593 was 7593, checked in by jgriffin, 8 years ago

Added technical paper on wave energy

File size: 9.0 KB
Line 
1"""
2Program to calcualte the energy density in offshore tsunami time series.
3Energy is integrated over the time of the each wave crest.
4Also calculates time integrated momentum quantities from the boundary time series.
5Used to investiage different methods of estimating offshore tsunami hazard,
6rather than just wave height.
7
8Energy calculations from Johnson 1997. 'A modern introduction to the mathematical
9theory of water waves' p23.
10
11Creator: Jonathan Griffin, Geoscience Australia
12Created: 12 October 2009
13"""
14
15
16def 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#############################################################################
218import os
219import sys
220from os.path import join
221import csv
222import glob
223import numpy
224import run_up
225if __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)
Note: See TracBrowser for help on using the repository browser.