# Changeset 7593

Ignore:
Timestamp:
Dec 17, 2009, 9:36:34 AM (9 years ago)
Message:

Added technical paper on wave energy

Location:
misc/tools/event_selection/wave_energy
Files:
1 added
2 edited

Unmodified
Added
Removed
• ## misc/tools/event_selection/wave_energy/run_up.py

 r7584 R_break = gravity*(power(slope,2.0))/(power(freq,2.0)) run_up = np.min(run_up_madsen, R_break) run_up = min(run_up_madsen, R_break) return run_up #def surf_sim(depth, slope, height, wavelength, gravity = 9.81): # Add Madsen's surf similarity here def surf_sim(depth, slope, height, period, gravity = 9.81): wavelength = period*sqrt(gravity*depth) print 'wavelength', wavelength surf_sim_param = slope/(sqrt(height/wavelength)) print 'surf_sim_param', surf_sim_param #print 'height',height ##    print 'depth', depth a = 2*(power(pi,(0.75))) HA = height/depth #print HA b = power(HA,-0.25) c =power(surf_sim_param,(-0.5)) ##    print 'a', a ##    print 'b' ,b ##    print 'c', c R = height*a*b*c #R = height*2*(power(pi,(3/4)))*(power((height/depth),(-1/4)))*(power(surf_sim_param,(-1/2))) print 'R',R R_break = height*(1.0/pi)*power(surf_sim_param,2.0) print 'R_break',R_break run_up = min(R, R_break) ##    print 'run_up',run_up return run_up import numpy as np from numpy import tan, sin, cos, power,pi from numpy import sin,tan,cos,power,pi, sqrt if __name__ == '__main__': sys.exit(main())
• ## misc/tools/event_selection/wave_energy/wave_energy.py

 r7584 """ 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 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 if __name__ == '__main__': ############################################################### # 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' # 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)
Note: See TracChangeset for help on using the changeset viewer.