""" Program to calcualte the energy density in offshore tsunami time series. Energy is integrated over the time of the each wave crest. Also calculates time integrated momentum quantities from the boundary time series. Used to investiage different methods of estimating offshore tsunami hazard, rather than just wave height. Energy calculations from Johnson 1997. 'A modern introduction to the mathematical theory of water waves' p23. Creator: Jonathan Griffin, Geoscience Australia Created: 12 October 2009 """ def wave_energy(filepath, filebasename): filepathlist = [] filenamelist = [] # Find all relevant gauges files and put into list for root, dirs, files in os.walk(filepath): for filenames in files: if filenames == filebasename: filepathname = join(root, filenames) filepathlist.append(filepathname) filenamelist.append(filenames) # Initialise lists for storing maximum values max_energy_list = [] max_energy_inst_list = [] run_up_list = [] run_up_list_surf_sim = [] energy_run_up_list = [] max_stage_list = [] counter = 0 ########################################################### # Loop over gauges files ########################################################### for filename in filepathlist: time = [] stage = [] x_momentum = [] y_momentum = [] momentum = [] print 'reading file ', filename csvreader = csv.reader(open(filename)) header = csvreader.next() for line in csvreader: time.append(float(line[0])) stage.append(float(line[1])) x_momentum.append(float(line[2])) y_momentum.append(float(line[3])) # Calculate momentum norm time_diff = time[1]-time[0] sign = 0 energy_sum = 0 max_energy = 0 max_energy_inst = 0 max_stage = 0 temp_max_stage = 0 max_dist = 0 ############################################################### # Loop through time series and add values (integrate over time) # Sum set to zero when stage height changes sign ############################################################### for i in range(len(time)): ############################################################### # Energy calculations - note that quantities from gauge file are # already depth integrated. ############################################################### # Calculate velocities from 'momentum' terms (Note not true momentum, # rather defined as u*h. So we devide by the depth to get u). # roh = density = 1000 kg.m-3 x_vel = x_momentum[i]/(100+stage[i]) y_vel = y_momentum[i]/(100+stage[i]) velocity = numpy.sqrt(x_vel*x_vel+y_vel*y_vel) # Kinetic energy KE = 1/2 roh h u^2 KE = (1./2.)*(numpy.power(x_vel,2)+numpy.power(y_vel,2))*(1000)*(100+stage[i]) # Potential energy PE = roh g w (w = stage) PE = 9.81*(stage[i])*1000 energy = KE+PE # Maximum instantaneous energy density if energy > max_energy_inst: max_energy_inst = energy ############################################################### # If sign of wave changes (goes from trough to crest or vice-versa) # multiply energy sum by timestep and reset sum value to zero ############################################################### if stage[i] > 0 and sign != 1: sign = 1 temp_max_energy = energy_sum*time_diff if temp_max_energy > max_energy: max_energy = temp_max_energy max_energy_time = time[i] - start_time if temp_max_stage > max_stage: max_stage = temp_max_stage max_time = time[i] - start_time energy_sum = energy start_time = time[i] #max_stage_temp = 0 elif stage[i] < 0 and sign != -1: sign = -1 temp_max_energy = energy_sum*time_diff if temp_max_energy > max_energy: max_energy = temp_max_energy max_energy_time = time[i] - start_time if temp_max_stage > max_stage: max_stage = temp_max_stage max_time = time[i] - start_time energy_sum = energy start_time = time[i] elif stage[i] == 0 and sign != 0: sign = 0 if temp_max_stage > max_stage: max_stage = temp_max_stage max_time = time[i] - start_time energy_sum = energy start_time = time[i] else: energy_sum += energy temp_max_stage = max(temp_max_stage,stage[i]) # Add maximum values to list max_energy_inst_list.append(max_energy_inst) max_energy_list.append(max_energy) # get wave period from max time ## locator = time.index(max_time) ## time_count = 0 ## for i in range(locator,len(time)): # Account for both halves of the wave max_time = max_time*2 print 'max stage',max_stage print 'max time',max_time ##################################################################### # call run-up module ##################################################################### Run_up = run_up.Madsen(100,0.017,max_stage,max_time) print 'Madsen runup', Run_up run_up_list.append(Run_up) Run_up_surf_sim = run_up.surf_sim(100,0.03,max_stage,max_time) print 'Madsen surf similarity runup', Run_up_surf_sim run_up_list_surf_sim.append(Run_up_surf_sim) energy_run_up = run_up.energy_conserv(max_energy,1.0,0.8,9.81) energy_run_up_list.append(energy_run_up) print 'energy run up', energy_run_up max_stage_list.append(max_stage) counter += 1 counter = 0 total_max_energy = 0 total_max_inst_energy = 0 # Print results for each file and calculate maximum value across all events outFilename= filepath +'/wave_energy_summary.csv' outFile = file(outFilename,'w') 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') for filename in filepathlist: outFile.write(filename+',') outFile.write(str(max_energy_list[counter])+',') outFile.write(str(max_energy_inst_list[counter])+',') outFile.write(str(run_up_list[counter])+',') outFile.write(str(energy_run_up_list[counter])+',') outFile.write(str(max_stage_list[counter])+',') outFile.write(str(run_up_list_surf_sim[counter])+'\n') print filename print 'max crest-integrated energy:', max_energy_list[counter],'J/m^2' print 'max instantatneous energy:', max_energy_inst_list[counter],'J/s*m^2' print 'Run_up (Madsen 2009):', run_up_list[counter], 'm' print 'Run_up surf similarity (Madsen 2009):', run_up_list_surf_sim[counter], 'm' print 'Run_up (Energy 2009):', energy_run_up_list[counter], 'm' if max_energy_list[counter] > total_max_energy: total_max_energy = max_energy_list[counter] max_energy_index = counter if max_energy_inst_list[counter] > total_max_inst_energy: total_max_inst_energy = max_energy_inst_list[counter] max_energy_inst_index = counter counter+=1 print '\n File with maximum crest-integrated energy is:', filepathlist[max_energy_index] print 'Energy is', max_energy_list[max_energy_index], 'J.s/m^2' print '\n File with maximum instantaneous energy is:', filepathlist[max_energy_inst_index] print 'Energy is', max_energy_inst_list[max_energy_inst_index], 'J/m^2' ############################################################################# ############################################################################# import os import sys from os.path import join import csv import glob import numpy import run_up if __name__ == '__main__': # Change path here!!! filepath = r'/nas/gemd/georisk_models/inundation/data/new_south_wales/batemans_bay_tsunami_scenario_2009/anuga/boundaries' index = 2872 filebasename = 'sts_gauge_' + str(index) +'.csv' wave_energy(filepath, filebasename)