| 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 |
|---|
| 5 | |
|---|
| 6 | Creator: Jonathan Griffin |
|---|
| 7 | Created: 12 October 2009 |
|---|
| 8 | """ |
|---|
| 9 | |
|---|
| 10 | import os |
|---|
| 11 | import sys |
|---|
| 12 | from os.path import join |
|---|
| 13 | import csv |
|---|
| 14 | import glob |
|---|
| 15 | ##import matplotlib |
|---|
| 16 | ##matplotlib.use('Agg') |
|---|
| 17 | ##import pylab |
|---|
| 18 | ##from pylab import * |
|---|
| 19 | import numpy |
|---|
| 20 | |
|---|
| 21 | filepath = r'/nas/gemd/georisk_models/inundation/data/new_south_wales/batemans_bay_tsunami_scenario_2009/anuga/boundaries' |
|---|
| 22 | index = 2872 |
|---|
| 23 | filebasename = 'sts_gauge_' + str(index) +'.csv' |
|---|
| 24 | filename = join(filepath, filebasename) |
|---|
| 25 | |
|---|
| 26 | filepathlist = [] |
|---|
| 27 | filenamelist = [] |
|---|
| 28 | |
|---|
| 29 | for root, dirs, files in os.walk(filepath): |
|---|
| 30 | ## print 'root', root |
|---|
| 31 | ## print 'dirs', dirs |
|---|
| 32 | ## print 'files', files |
|---|
| 33 | for filenames in files: |
|---|
| 34 | if filenames == filebasename: |
|---|
| 35 | filepathname = join(root, filenames) |
|---|
| 36 | filepathlist.append(filepathname) |
|---|
| 37 | filenamelist.append(filenames) |
|---|
| 38 | |
|---|
| 39 | max_energy_list = [] |
|---|
| 40 | max_momentum_list = [] |
|---|
| 41 | max_energy_inst_list = [] |
|---|
| 42 | |
|---|
| 43 | counter = 0 |
|---|
| 44 | for filename in filepathlist: |
|---|
| 45 | time = [] |
|---|
| 46 | stage = [] |
|---|
| 47 | x_momentum = [] |
|---|
| 48 | y_momentum = [] |
|---|
| 49 | momentum = [] |
|---|
| 50 | print 'reading file ', filename |
|---|
| 51 | csvreader = csv.reader(open(filename)) |
|---|
| 52 | header = csvreader.next() |
|---|
| 53 | for line in csvreader: |
|---|
| 54 | time.append(float(line[0])) |
|---|
| 55 | stage.append(float(line[1])) |
|---|
| 56 | x_momentum.append(float(line[2])) |
|---|
| 57 | y_momentum.append(float(line[3])) |
|---|
| 58 | # Calculate momentum norm |
|---|
| 59 | momentum.append(numpy.sqrt(numpy.power(float(line[2]),2)+numpy.power(float(line[3]),2))) |
|---|
| 60 | |
|---|
| 61 | time_diff = time[1]-time[0] |
|---|
| 62 | |
|---|
| 63 | sign = 0 |
|---|
| 64 | max_momentum = 0 |
|---|
| 65 | momentum_sum = 0 |
|---|
| 66 | energy_sum = 0 |
|---|
| 67 | max_energy = 0 |
|---|
| 68 | max_energy_inst = 0 |
|---|
| 69 | |
|---|
| 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]) |
|---|
| 76 | energy = KE+PE |
|---|
| 77 | |
|---|
| 78 | if energy > max_energy_inst: |
|---|
| 79 | max_energy_inst = energy |
|---|
| 80 | |
|---|
| 81 | ## print 'max_energy',max_energy |
|---|
| 82 | ## print 'energy_sum', energy_sum |
|---|
| 83 | if stage[i] > 0 and sign != 1: |
|---|
| 84 | sign = 1 |
|---|
| 85 | temp_max_momentum = momentum_sum*time_diff |
|---|
| 86 | temp_max_energy = energy_sum*time_diff |
|---|
| 87 | if temp_max_momentum > max_momentum: |
|---|
| 88 | max_momentum = temp_max_momentum |
|---|
| 89 | if temp_max_energy > max_energy: |
|---|
| 90 | max_energy = temp_max_energy |
|---|
| 91 | energy_sum = energy |
|---|
| 92 | momentum_sum = momentum[i] |
|---|
| 93 | |
|---|
| 94 | start_time = time[i] |
|---|
| 95 | |
|---|
| 96 | elif stage[i] < 0 and sign != -1: |
|---|
| 97 | sign = -1 |
|---|
| 98 | temp_max_momentum = momentum_sum*time_diff |
|---|
| 99 | temp_max_energy = energy_sum*time_diff |
|---|
| 100 | if temp_max_momentum > max_momentum: |
|---|
| 101 | max_momentum = temp_max_momentum |
|---|
| 102 | if temp_max_energy > max_energy: |
|---|
| 103 | max_energy = temp_max_energy |
|---|
| 104 | energy_sum = energy |
|---|
| 105 | momentum_sum = momentum[i] |
|---|
| 106 | # print momentum_sum |
|---|
| 107 | start_time = time[i] |
|---|
| 108 | |
|---|
| 109 | elif stage[i] == 0 and sign != 0: |
|---|
| 110 | sign = 0 |
|---|
| 111 | momentum_sum = momentum[i] |
|---|
| 112 | energy_sum = energy |
|---|
| 113 | |
|---|
| 114 | start_time = time[i] |
|---|
| 115 | |
|---|
| 116 | else: |
|---|
| 117 | momentum_sum += momentum[i] |
|---|
| 118 | energy_sum += energy |
|---|
| 119 | |
|---|
| 120 | ## print 'KE',KE |
|---|
| 121 | ## print 'PE',PE |
|---|
| 122 | |
|---|
| 123 | max_energy_inst_list.append(max_energy_inst) |
|---|
| 124 | 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 |
|---|
| 128 | counter += 1 |
|---|
| 129 | counter = 0 |
|---|
| 130 | total_max_energy = 0 |
|---|
| 131 | total_max_momentum = 0 |
|---|
| 132 | total_max_inst_energy = 0 |
|---|
| 133 | |
|---|
| 134 | for filename in filepathlist: |
|---|
| 135 | 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] |
|---|
| 139 | |
|---|
| 140 | if max_energy_list[counter] > total_max_energy: |
|---|
| 141 | total_max_energy = max_energy_list[counter] |
|---|
| 142 | max_energy_index = counter |
|---|
| 143 | if max_energy_inst_list[counter] > total_max_inst_energy: |
|---|
| 144 | total_max_inst_energy = max_energy_inst_list[counter] |
|---|
| 145 | 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 |
|---|
| 149 | counter+=1 |
|---|
| 150 | |
|---|
| 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' |
|---|
