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

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

Added energy run-up calculations and breaking criteria to Madsen run-up

File size: 7.7 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
15import os
16import sys
17from os.path import join
18import csv
19import glob
20import numpy
21import run_up
22
23# Change path here!!!
24filepath = r'/nas/gemd/georisk_models/inundation/data/new_south_wales/batemans_bay_tsunami_scenario_2009/anuga/boundaries'
25index = 2872
26filebasename = 'sts_gauge_' + str(index) +'.csv'
27#filename = join(filepath, filebasename)
28
29filepathlist = []
30filenamelist = []
31
32
33# Find all relevant gauges files and put into list
34for 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   
43max_energy_list = []
44max_energy_inst_list = []
45run_up_list = []
46energy_run_up_list = []
47max_stage_list = []
48counter = 0
49
50###########################################################
51# Loop over gauges files
52###########################################################   
53for 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
77   
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   
178counter = 0
179total_max_energy = 0
180total_max_inst_energy = 0
181
182
183
184
185
186# Print results for each file and calculate maximum value across all events
187outFilename= filepath +'/wave_energy_summary.csv'
188outFile = file(outFilename,'w')
189outFile.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')
190for 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
212print '\n File with maximum crest-integrated energy is:', filepathlist[max_energy_index]
213print 'Energy is', max_energy_list[max_energy_index], 'J.s/m^2'
214print '\n File with maximum instantaneous energy is:', filepathlist[max_energy_inst_index]
215print 'Energy is', max_energy_inst_list[max_energy_inst_index], 'J/m^2'
216
217
Note: See TracBrowser for help on using the repository browser.