| 1 | """ |
|---|
| 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. |
|---|
| 7 | |
|---|
| 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 |
|---|
| 12 | Created: 12 October 2009 |
|---|
| 13 | """ |
|---|
| 14 | |
|---|
| 15 | import os |
|---|
| 16 | import sys |
|---|
| 17 | from os.path import join |
|---|
| 18 | import csv |
|---|
| 19 | import glob |
|---|
| 20 | import numpy |
|---|
| 21 | |
|---|
| 22 | # Change path here!!! |
|---|
| 23 | filepath = r'/nas/gemd/georisk_models/inundation/data/new_south_wales/batemans_bay_tsunami_scenario_2009/anuga/boundaries' |
|---|
| 24 | index = 2872 |
|---|
| 25 | filebasename = 'sts_gauge_' + str(index) +'.csv' |
|---|
| 26 | filename = join(filepath, filebasename) |
|---|
| 27 | |
|---|
| 28 | filepathlist = [] |
|---|
| 29 | filenamelist = [] |
|---|
| 30 | |
|---|
| 31 | |
|---|
| 32 | # Find all relevant gauges files and put into list |
|---|
| 33 | for root, dirs, files in os.walk(filepath): |
|---|
| 34 | |
|---|
| 35 | for filenames in files: |
|---|
| 36 | if filenames == filebasename: |
|---|
| 37 | filepathname = join(root, filenames) |
|---|
| 38 | filepathlist.append(filepathname) |
|---|
| 39 | filenamelist.append(filenames) |
|---|
| 40 | |
|---|
| 41 | # Initialise lists for storing maximum values |
|---|
| 42 | max_energy_list = [] |
|---|
| 43 | max_energy_inst_list = [] |
|---|
| 44 | counter = 0 |
|---|
| 45 | |
|---|
| 46 | ########################################################### |
|---|
| 47 | # Loop over gauges files |
|---|
| 48 | ########################################################### |
|---|
| 49 | for filename in filepathlist: |
|---|
| 50 | time = [] |
|---|
| 51 | stage = [] |
|---|
| 52 | x_momentum = [] |
|---|
| 53 | y_momentum = [] |
|---|
| 54 | momentum = [] |
|---|
| 55 | print 'reading file ', filename |
|---|
| 56 | csvreader = csv.reader(open(filename)) |
|---|
| 57 | header = csvreader.next() |
|---|
| 58 | for line in csvreader: |
|---|
| 59 | time.append(float(line[0])) |
|---|
| 60 | stage.append(float(line[1])) |
|---|
| 61 | x_momentum.append(float(line[2])) |
|---|
| 62 | y_momentum.append(float(line[3])) |
|---|
| 63 | # Calculate momentum norm |
|---|
| 64 | |
|---|
| 65 | time_diff = time[1]-time[0] |
|---|
| 66 | |
|---|
| 67 | sign = 0 |
|---|
| 68 | energy_sum = 0 |
|---|
| 69 | max_energy = 0 |
|---|
| 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)): |
|---|
| 77 | |
|---|
| 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 |
|---|
| 91 | energy = KE+PE |
|---|
| 92 | |
|---|
| 93 | # Maximum instantaneous energy density |
|---|
| 94 | if energy > max_energy_inst: |
|---|
| 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 | ############################################################### |
|---|
| 101 | |
|---|
| 102 | if stage[i] > 0 and sign != 1: |
|---|
| 103 | sign = 1 |
|---|
| 104 | temp_max_energy = energy_sum*time_diff |
|---|
| 105 | if temp_max_energy > max_energy: |
|---|
| 106 | max_energy = temp_max_energy |
|---|
| 107 | |
|---|
| 108 | energy_sum = energy |
|---|
| 109 | start_time = time[i] |
|---|
| 110 | |
|---|
| 111 | elif stage[i] < 0 and sign != -1: |
|---|
| 112 | sign = -1 |
|---|
| 113 | temp_max_energy = energy_sum*time_diff |
|---|
| 114 | if temp_max_energy > max_energy: |
|---|
| 115 | max_energy = temp_max_energy |
|---|
| 116 | |
|---|
| 117 | energy_sum = energy |
|---|
| 118 | start_time = time[i] |
|---|
| 119 | |
|---|
| 120 | elif stage[i] == 0 and sign != 0: |
|---|
| 121 | sign = 0 |
|---|
| 122 | energy_sum = energy |
|---|
| 123 | start_time = time[i] |
|---|
| 124 | |
|---|
| 125 | else: |
|---|
| 126 | energy_sum += energy |
|---|
| 127 | |
|---|
| 128 | # Add maximum values to list |
|---|
| 129 | max_energy_inst_list.append(max_energy_inst) |
|---|
| 130 | max_energy_list.append( max_energy) |
|---|
| 131 | |
|---|
| 132 | counter += 1 |
|---|
| 133 | |
|---|
| 134 | |
|---|
| 135 | counter = 0 |
|---|
| 136 | total_max_energy = 0 |
|---|
| 137 | total_max_inst_energy = 0 |
|---|
| 138 | |
|---|
| 139 | # Print results for each file and calculate maximum value across all events |
|---|
| 140 | for filename in filepathlist: |
|---|
| 141 | print filename |
|---|
| 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' |
|---|
| 144 | |
|---|
| 145 | if max_energy_list[counter] > total_max_energy: |
|---|
| 146 | total_max_energy = max_energy_list[counter] |
|---|
| 147 | max_energy_index = counter |
|---|
| 148 | if max_energy_inst_list[counter] > total_max_inst_energy: |
|---|
| 149 | total_max_inst_energy = max_energy_inst_list[counter] |
|---|
| 150 | max_energy_inst_index = counter |
|---|
| 151 | |
|---|
| 152 | counter+=1 |
|---|
| 153 | |
|---|
| 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' |
|---|
