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

"""
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
"""

import os
import sys
from os.path import join
import csv
import glob
import numpy
import run_up

# 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'
#filename = join(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 = []
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
    
    ###############################################################
    # 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)
    
    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])+'\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 (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'