Changeset 7593


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

Added technical paper on wave energy

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

Legend:

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

    r7584 r7593  
    2020
    2121    R_break = gravity*(power(slope,2.0))/(power(freq,2.0))
    22     run_up = np.min(run_up_madsen, R_break)
     22    run_up = min(run_up_madsen, R_break)
    2323    return run_up
    2424
    25 #def surf_sim(depth, slope, height, wavelength, gravity = 9.81):
    26 # Add Madsen's surf similarity here
     25def surf_sim(depth, slope, height, period, gravity = 9.81):
     26    wavelength = period*sqrt(gravity*depth)
     27    print 'wavelength', wavelength
     28    surf_sim_param = slope/(sqrt(height/wavelength))
     29    print 'surf_sim_param', surf_sim_param
     30    #print 'height',height
     31##    print 'depth', depth
     32    a = 2*(power(pi,(0.75)))
     33    HA = height/depth
     34    #print HA
     35    b = power(HA,-0.25)
     36    c =power(surf_sim_param,(-0.5))
     37##    print 'a', a
     38##    print 'b' ,b
     39##    print 'c', c
     40    R = height*a*b*c
     41    #R = height*2*(power(pi,(3/4)))*(power((height/depth),(-1/4)))*(power(surf_sim_param,(-1/2)))
     42    print 'R',R
     43    R_break = height*(1.0/pi)*power(surf_sim_param,2.0)
     44    print 'R_break',R_break
     45    run_up = min(R, R_break)
     46##    print 'run_up',run_up
     47    return run_up
     48   
     49
    2750
    2851
     
    4366
    4467import numpy as np
    45 from numpy import tan, sin, cos, power,pi
     68from numpy import sin,tan,cos,power,pi, sqrt
    4669if __name__ == '__main__':
    4770    sys.exit(main())
  • misc/tools/event_selection/wave_energy/wave_energy.py

    r7584 r7593  
    1313"""
    1414
     15
     16def wave_energy(filepath, filebasename):
     17   
     18    filepathlist = []
     19    filenamelist = []
     20
     21
     22    # Find all relevant gauges files and put into list
     23    for root, dirs, files in os.walk(filepath):
     24
     25        for filenames in files:
     26            if filenames == filebasename:
     27                filepathname = join(root, filenames)
     28                filepathlist.append(filepathname)
     29                filenamelist.append(filenames)
     30               
     31    # Initialise lists for storing maximum values   
     32    max_energy_list = []
     33    max_energy_inst_list = []
     34    run_up_list = []
     35    run_up_list_surf_sim = []
     36    energy_run_up_list = []
     37    max_stage_list = []
     38    counter = 0
     39
     40    ###########################################################
     41    # Loop over gauges files
     42    ###########################################################   
     43    for filename in filepathlist:
     44        time = []
     45        stage = []
     46        x_momentum = []
     47        y_momentum = []
     48        momentum = []
     49        print 'reading file ', filename
     50        csvreader = csv.reader(open(filename))
     51        header = csvreader.next()
     52        for line in csvreader: 
     53            time.append(float(line[0]))
     54            stage.append(float(line[1]))
     55            x_momentum.append(float(line[2]))
     56            y_momentum.append(float(line[3]))
     57            # Calculate momentum norm
     58           
     59        time_diff = time[1]-time[0]
     60
     61        sign = 0
     62        energy_sum = 0
     63        max_energy = 0
     64        max_energy_inst = 0
     65        max_stage = 0
     66        temp_max_stage = 0
     67        max_dist = 0
     68       
     69        ###############################################################
     70        # Loop through time series and add values (integrate over time)
     71        # Sum set to zero when stage height changes sign
     72        ###############################################################
     73        for i in range(len(time)):
     74
     75            ###############################################################
     76            # Energy calculations - note that quantities from gauge file are
     77            # already depth integrated.
     78            ###############################################################
     79            # Calculate velocities from 'momentum' terms (Note not true momentum,
     80            # rather defined as u*h. So we devide by the depth to get u).
     81            # roh = density = 1000 kg.m-3
     82            x_vel = x_momentum[i]/(100+stage[i])
     83            y_vel = y_momentum[i]/(100+stage[i])
     84            velocity = numpy.sqrt(x_vel*x_vel+y_vel*y_vel)
     85            # Kinetic energy KE = 1/2 roh h u^2
     86            KE = (1./2.)*(numpy.power(x_vel,2)+numpy.power(y_vel,2))*(1000)*(100+stage[i])
     87            # Potential energy PE = roh g w (w = stage)
     88            PE = 9.81*(stage[i])*1000
     89            energy = KE+PE
     90
     91            # Maximum instantaneous energy density
     92            if energy > max_energy_inst:
     93                max_energy_inst = energy
     94
     95            ###############################################################
     96            # If sign of wave changes (goes from trough to crest or vice-versa)
     97            # multiply energy sum by timestep and reset sum value to zero
     98            ###############################################################
     99           
     100            if stage[i] > 0 and sign != 1:
     101                sign = 1
     102                temp_max_energy = energy_sum*time_diff
     103                if temp_max_energy > max_energy:
     104                    max_energy = temp_max_energy
     105                    max_energy_time = time[i] - start_time
     106                if temp_max_stage > max_stage:
     107                    max_stage = temp_max_stage
     108                    max_time = time[i] - start_time   
     109                energy_sum = energy
     110                start_time = time[i]
     111                #max_stage_temp = 0
     112               
     113            elif stage[i] < 0 and sign != -1:
     114                sign = -1
     115                temp_max_energy = energy_sum*time_diff
     116                if temp_max_energy > max_energy:
     117                    max_energy = temp_max_energy
     118                    max_energy_time = time[i] - start_time
     119                if temp_max_stage > max_stage:
     120                    max_stage = temp_max_stage
     121                    max_time = time[i] - start_time                 
     122                energy_sum = energy
     123                start_time = time[i]
     124               
     125            elif stage[i] == 0 and sign != 0:
     126                sign = 0
     127                if temp_max_stage > max_stage:
     128                    max_stage = temp_max_stage
     129                    max_time = time[i] - start_time
     130                energy_sum = energy
     131                start_time = time[i]
     132                         
     133            else:
     134                energy_sum += energy
     135                temp_max_stage = max(temp_max_stage,stage[i])
     136               
     137        # Add maximum values to list       
     138        max_energy_inst_list.append(max_energy_inst)
     139        max_energy_list.append(max_energy)
     140
     141        # get wave period from max time
     142    ##    locator = time.index(max_time)
     143    ##    time_count = 0
     144    ##    for i in range(locator,len(time)):
     145           
     146
     147        # Account for both halves of the wave
     148        max_time = max_time*2
     149        print 'max stage',max_stage
     150        print 'max time',max_time
     151        #####################################################################
     152        # call run-up module
     153        #####################################################################
     154        Run_up = run_up.Madsen(100,0.017,max_stage,max_time)
     155        print 'Madsen runup', Run_up
     156        run_up_list.append(Run_up)
     157
     158        Run_up_surf_sim = run_up.surf_sim(100,0.03,max_stage,max_time)
     159        print 'Madsen surf similarity runup', Run_up_surf_sim
     160        run_up_list_surf_sim.append(Run_up_surf_sim)
     161       
     162        energy_run_up  = run_up.energy_conserv(max_energy,1.0,0.8,9.81)
     163        energy_run_up_list.append(energy_run_up)
     164        print 'energy run up', energy_run_up
     165        max_stage_list.append(max_stage)
     166
     167       
     168        counter += 1
     169       
     170
     171
     172       
     173    counter = 0
     174    total_max_energy = 0
     175    total_max_inst_energy = 0
     176
     177
     178
     179
     180
     181    # Print results for each file and calculate maximum value across all events
     182    outFilename= filepath +'/wave_energy_summary.csv'
     183    outFile = file(outFilename,'w')
     184    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')
     185    for filename in filepathlist:
     186        outFile.write(filename+',')
     187        outFile.write(str(max_energy_list[counter])+',')
     188        outFile.write(str(max_energy_inst_list[counter])+',')
     189        outFile.write(str(run_up_list[counter])+',')
     190        outFile.write(str(energy_run_up_list[counter])+',')
     191        outFile.write(str(max_stage_list[counter])+',')
     192        outFile.write(str(run_up_list_surf_sim[counter])+'\n')
     193        print filename
     194        print 'max crest-integrated energy:', max_energy_list[counter],'J/m^2'
     195        print 'max instantatneous energy:', max_energy_inst_list[counter],'J/s*m^2'
     196        print 'Run_up (Madsen 2009):', run_up_list[counter], 'm'
     197        print 'Run_up surf similarity (Madsen 2009):', run_up_list_surf_sim[counter], 'm'
     198        print 'Run_up (Energy 2009):', energy_run_up_list[counter], 'm'
     199       
     200        if max_energy_list[counter] > total_max_energy:
     201            total_max_energy = max_energy_list[counter]
     202            max_energy_index = counter
     203        if max_energy_inst_list[counter] > total_max_inst_energy:
     204            total_max_inst_energy = max_energy_inst_list[counter]
     205            max_energy_inst_index = counter
     206
     207        counter+=1
     208
     209    print '\n File with maximum crest-integrated energy is:', filepathlist[max_energy_index]
     210    print 'Energy is', max_energy_list[max_energy_index], 'J.s/m^2'
     211    print '\n File with maximum instantaneous energy is:', filepathlist[max_energy_inst_index]
     212    print 'Energy is', max_energy_inst_list[max_energy_inst_index], 'J/m^2'
     213
     214
     215#############################################################################
     216   
     217#############################################################################
    15218import os
    16219import sys
     
    20223import numpy
    21224import run_up
    22 
    23 # Change path here!!!
    24 filepath = r'/nas/gemd/georisk_models/inundation/data/new_south_wales/batemans_bay_tsunami_scenario_2009/anuga/boundaries'
    25 index = 2872
    26 filebasename = 'sts_gauge_' + str(index) +'.csv'
    27 #filename = join(filepath, filebasename)
    28 
    29 filepathlist = []
    30 filenamelist = []
    31 
    32 
    33 # Find all relevant gauges files and put into list
    34 for root, dirs, files in os.walk(filepath):
    35 
    36     for filenames in files:
    37         if filenames == filebasename:
    38             filepathname = join(root, filenames)
    39             filepathlist.append(filepathname)
    40             filenamelist.append(filenames)
    41            
    42 # Initialise lists for storing maximum values   
    43 max_energy_list = []
    44 max_energy_inst_list = []
    45 run_up_list = []
    46 energy_run_up_list = []
    47 max_stage_list = []
    48 counter = 0
    49 
    50 ###########################################################
    51 # Loop over gauges files
    52 ###########################################################   
    53 for filename in filepathlist:
    54     time = []
    55     stage = []
    56     x_momentum = []
    57     y_momentum = []
    58     momentum = []
    59     print 'reading file ', filename
    60     csvreader = csv.reader(open(filename))
    61     header = csvreader.next()
    62     for line in csvreader: 
    63         time.append(float(line[0]))
    64         stage.append(float(line[1]))
    65         x_momentum.append(float(line[2]))
    66         y_momentum.append(float(line[3]))
    67         # Calculate momentum norm
    68        
    69     time_diff = time[1]-time[0]
    70 
    71     sign = 0
    72     energy_sum = 0
    73     max_energy = 0
    74     max_energy_inst = 0
    75     max_stage = 0
    76     temp_max_stage = 0
     225if __name__ == '__main__':
    77226   
    78     ###############################################################
    79     # Loop through time series and add values (integrate over time)
    80     # Sum set to zero when stage height changes sign
    81     ###############################################################
    82     for i in range(len(time)):
    83 
    84         ###############################################################
    85         # Energy calculations - note that quantities from gauge file are
    86         # already depth integrated.
    87         ###############################################################
    88         # Calculate velocities from 'momentum' terms (Note not true momentum,
    89         # rather defined as u*h. So we devide by the depth to get u).
    90         # roh = density = 1000 kg.m-3
    91         x_vel = x_momentum[i]/(100+stage[i])
    92         y_vel = y_momentum[i]/(100+stage[i])
    93         velocity = numpy.sqrt(x_vel*x_vel+y_vel*y_vel)
    94         # Kinetic energy KE = 1/2 roh h u^2
    95         KE = (1./2.)*(numpy.power(x_vel,2)+numpy.power(y_vel,2))*(1000)*(100+stage[i])
    96         # Potential energy PE = roh g w (w = stage)
    97         PE = 9.81*(stage[i])*1000
    98         energy = KE+PE
    99 
    100         # Maximum instantaneous energy density
    101         if energy > max_energy_inst:
    102             max_energy_inst = energy
    103 
    104         ###############################################################
    105         # If sign of wave changes (goes from trough to crest or vice-versa)
    106         # multiply energy sum by timestep and reset sum value to zero
    107         ###############################################################
    108        
    109         if stage[i] > 0 and sign != 1:
    110             sign = 1
    111             temp_max_energy = energy_sum*time_diff
    112             if temp_max_energy > max_energy:
    113                 max_energy = temp_max_energy
    114                 max_energy_time = time[i] - start_time
    115             if temp_max_stage > max_stage:
    116                 max_stage = temp_max_stage
    117                 max_time = time[i] - start_time   
    118             energy_sum = energy
    119             start_time = time[i]
    120             max_stage_temp = 0
    121            
    122         elif stage[i] < 0 and sign != -1:
    123             sign = -1
    124             temp_max_energy = energy_sum*time_diff
    125             if temp_max_energy > max_energy:
    126                 max_energy = temp_max_energy
    127                 max_energy_time = time[i] - start_time
    128             if temp_max_stage > max_stage:
    129                 max_stage = temp_max_stage
    130                 max_time = time[i] - start_time                 
    131             energy_sum = energy
    132             start_time = time[i]
    133            
    134         elif stage[i] == 0 and sign != 0:
    135             sign = 0
    136             if temp_max_stage > max_stage:
    137                 max_stage = temp_max_stage
    138                 max_time = time[i] - start_time
    139             energy_sum = energy
    140             start_time = time[i]
    141                      
    142         else:
    143             energy_sum += energy
    144             temp_max_stage = max(temp_max_stage,stage[i])
    145            
    146     # Add maximum values to list       
    147     max_energy_inst_list.append(max_energy_inst)
    148     max_energy_list.append(max_energy)
    149 
    150     # get wave period from max time
    151 ##    locator = time.index(max_time)
    152 ##    time_count = 0
    153 ##    for i in range(locator,len(time)):
    154        
    155 
    156     # Account for both halves of the wave
    157     max_time = max_time*2
    158     print 'max stage',max_stage
    159     print 'max time',max_time
    160     #####################################################################
    161     # call run-up module
    162     #####################################################################
    163     Run_up = run_up.Madsen(100,0.017,max_stage,max_time)
    164     print 'Madsen runup', Run_up
    165     run_up_list.append(Run_up)
    166    
    167     energy_run_up  = run_up.energy_conserv(max_energy,1.0,0.8,9.81)
    168     energy_run_up_list.append(energy_run_up)
    169     print 'energy run up', energy_run_up
    170     max_stage_list.append(max_stage)
    171 
    172    
    173     counter += 1
    174    
    175 
    176 
    177    
    178 counter = 0
    179 total_max_energy = 0
    180 total_max_inst_energy = 0
    181 
    182 
    183 
    184 
    185 
    186 # Print results for each file and calculate maximum value across all events
    187 outFilename= filepath +'/wave_energy_summary.csv'
    188 outFile = file(outFilename,'w')
    189 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')
    190 for filename in filepathlist:
    191     outFile.write(filename+',')
    192     outFile.write(str(max_energy_list[counter])+',')
    193     outFile.write(str(max_energy_inst_list[counter])+',')
    194     outFile.write(str(run_up_list[counter])+',')
    195     outFile.write(str(energy_run_up_list[counter])+',')
    196     outFile.write(str(max_stage_list[counter])+'\n')
    197     print filename
    198     print 'max crest-integrated energy:', max_energy_list[counter],'J/m^2'
    199     print 'max instantatneous energy:', max_energy_inst_list[counter],'J/s*m^2'
    200     print 'Run_up (Madsen 2009):', run_up_list[counter], 'm'
    201     print 'Run_up (Energy 2009):', energy_run_up_list[counter], 'm'
    202    
    203     if max_energy_list[counter] > total_max_energy:
    204         total_max_energy = max_energy_list[counter]
    205         max_energy_index = counter
    206     if max_energy_inst_list[counter] > total_max_inst_energy:
    207         total_max_inst_energy = max_energy_inst_list[counter]
    208         max_energy_inst_index = counter
    209 
    210     counter+=1
    211 
    212 print '\n File with maximum crest-integrated energy is:', filepathlist[max_energy_index]
    213 print 'Energy is', max_energy_list[max_energy_index], 'J.s/m^2'
    214 print '\n File with maximum instantaneous energy is:', filepathlist[max_energy_inst_index]
    215 print 'Energy is', max_energy_inst_list[max_energy_inst_index], 'J/m^2'
    216 
    217 
     227    # Change path here!!!
     228    filepath = r'/nas/gemd/georisk_models/inundation/data/new_south_wales/batemans_bay_tsunami_scenario_2009/anuga/boundaries'
     229    index = 2872
     230    filebasename = 'sts_gauge_' + str(index) +'.csv'
     231    wave_energy(filepath, filebasename)
Note: See TracChangeset for help on using the changeset viewer.