Changeset 3190


Ignore:
Timestamp:
Jun 21, 2006, 9:31:57 AM (19 years ago)
Author:
sexton
Message:

MOST and ANUGA comparisons and updates

Files:
4 added
1 deleted
40 edited

Legend:

Unmodified
Added
Removed

  • documentation/experimentation/boundary_ANUGA_MOST/report/MOST_ANUGA.tex

    r3189 r3190  
    226226see for example the output for the Ocean polygon 1 and 2 locations.
    227227
     228\input{50100MOSTcomparison_onslow}
    228229
    229230It is more instructive in this case to compare differences in

  • documentation/user_manual/examples/project.py

    r3150 r3190  
    8383# demo poly which isn't rectangular
    8484j0 = [385000, 6280000]
    85 j1 = [350000, 6270000]
    86 j2 = [325000, 6263000]
    87 j3 = [325000, 6260000]
     85j1 = [360000, 6272500]
     86j2 = [335000, 6272500]
     87j3 = [330000, 6265000]
     88j31 = [325000, 6260000]
    8889j4 = [316000, 6260000]
    8990j5 = [316000, 6247000]
    … …  
    9192j7 = [385000, 6238000]
    9293
    93 demopoly = [j0, j2, j3, j4, j5, j6, j7]
     94demopoly = [j0, j1, j2, j3, j31, j4, j5, j6, j7]
    9495
     96polygonptsfile4 = polygondir + 'poly1'
     97polygonptsfile0 = polygondir + 'poly2'
    9598polygonptsfile1 = polygondir + 'poly3'
    9699polygonptsfile2 = polygondir + 'poly4'
    97100polygonptsfile3 = polygondir + 'poly5'
    98 northern_polygon = read_polygon(polygonptsfile1 + '.csv')
     101northern_polygon = read_polygon(polygonptsfile0 + '.csv')
     102manly_polygon = read_polygon(polygonptsfile1 + '.csv')
    99103harbour_polygon = read_polygon(polygonptsfile2 + '.csv')
    100104southern_polygon = read_polygon(polygonptsfile3 + '.csv')
     105top_polygon = read_polygon(polygonptsfile4 + '.csv')
    101106
    102107#Interior regions - the Harbour

  • documentation/user_manual/examples/runsydney.py

    r3150 r3190  
    5757interior_regions = [[project.northern_polygon, interior_res],
    5858                    [project.harbour_polygon, interior_res],
    59                     [project.southern_polygon, interior_res]]
     59                    [project.southern_polygon, interior_res],
     60                    [project.manly_polygon, interior_res],
     61                    [project.top_polygon, interior_res]]
    6062
    6163
    6264create_mesh_from_regions(project.demopoly,
    6365                         boundary_tags= {'oceannorth': [0],
    64                                          'onshorewest': [1],
    65                                          'harbournorth': [2],
    66                                          'harbourwest': [3],
    67                                          'harboursouth': [4],
    68                                          'oceansouth': [5],
    69                                          'oceaneast': [6]},
     66                                         'onshorenorth1': [1],
     67                                         'onshorenorthwest1': [2],
     68                                         'onshorenorthwest2': [3],
     69                                         'onshorenorth2': [4],
     70                                         'onshorewest': [5],
     71                                         'onshoresouth': [6],
     72                                         'oceansouth': [7],
     73                                         'oceaneast': [8]},
    7074                         maximum_triangle_area=100000,
    7175                         filename=meshname,           
    … …  
    115119                    verbose = True)
    116120
    117 
    118121#------------------------------------------------------------------------------
    119122# Setup boundary conditions (all Dirichlet)
    … …  
    123126
    124127Bd = Dirichlet_boundary([0.0,0.0,0.0])
    125 domain.set_boundary( {'oceannorth': Bd, 'onshorewest': Bd,
    126                       'harbournorth': Bd, 'harbourwest': Bd,
    127                       'harboursouth': Bd, 'oceansouth': Bd,
    128                       'oceaneast': Bd})
    129 
     128domain.set_boundary( {'oceannorth': Bd,
     129                      'onshorenorth1': Bd,
     130                      'onshorenorthwest1': Bd,
     131                      'onshorenorthwest2': Bd,
     132                      'onshorenorth2': Bd,
     133                      'onshorewest': Bd,
     134                      'onshoresouth': Bd,
     135                      'oceansouth': Bd,
     136                      'oceaneast': Bd} )
    130137
    131138#------------------------------------------------------------------------------
    … …  
    135142import time
    136143t0 = time.time()
    137 tend = 7200.0
    138 print 'Simulate to time = %d secs' %tend
    139144
    140 for t in domain.evolve(yieldstep = 60, duration = tend):
     145for t in domain.evolve(yieldstep = 60, finaltime = 7200):
    141146    print domain.timestepping_statistics()
    142147    print domain.boundary_statistics(tags = 'oceaneast')   

  • inundation/pyvolution/smf.py

    r3136 r3190  
    346346        kappa = self.kappa
    347347        kappad = self.kappad
    348         #amin = self.find_min(x0,wa,kappad,kappa,dx,am)
    349         amin = 1.0
     348        amin = self.find_min(x0,wa,kappad,kappa,dx,am)
     349        print 'hello amin', amin
     350        #amin = 1.0
    350351
    351352        #double Gaussian calculation assumes water displacement is oriented
    … …  
    369370                if z[i] > maxz: maxz = z[i]
    370371                if z[i] < minz: minz = z[i]
     372               
    371373            except OverflowError:
    372374                pass
    373        
     375
     376        print 'max z', maxz
     377        print 'min z', minz
     378               
    374379        return z
    375380
    … …  
    407412        count_max = 1000000
    408413        c = 0
     414        am = 1.0
    409415
    410416        while c < count_max and deriv > 0:

  • inundation/pyvolution/util.py

    r3176 r3190  
    548548                   reportname = None,
    549549                   plot_quantity = None,
     550                   surface = None,
    550551                   time_min = None,
    551552                   time_max = None,
    … …  
    596597                    - default will be ['stage', 'speed', 'bearing']
    597598
     599    surface          - if True, then generate solution surface with 3d plot
     600                       and save to current working directory
     601   
    598602    time_min         - beginning of user defined time range for plotting purposes
    599603                        - default will be first available time found in swwfile
    … …  
    631635                        reportname,
    632636                        plot_quantity,
     637                        surface,
    633638                        time_min,
    634639                        time_max,
    … …  
    644649                    reportname = None,
    645650                    plot_quantity = None,
     651                    surface = None,
    646652                    time_min = None,
    647653                    time_max = None,
    … …  
    670676        check_list(plot_quantity)
    671677
     678    if surface is None:
     679        surface = False
     680       
    672681    if title_on is None:
    673682        title_on = True
    … …  
    727736    if verbose: print 'Inputs OK - going to generate figures'
    728737
    729     return generate_figures(plot_quantity, file_loc, report, reportname, leg_label,
    730                             f_list, gauges, locations, elev, production_dirs,
     738    return generate_figures(plot_quantity, file_loc, report, reportname, surface,
     739                            leg_label, f_list, gauges, locations, elev, production_dirs,
    731740                            time_min, time_max, title_on, label_id, verbose)
    732741                         
    … …  
    789798    return bearing
    790799
    791 def generate_figures(plot_quantity, file_loc, report, reportname, leg_label, f_list, gauges,
    792                      locations, elev, production_dirs, time_min, time_max,
     800def generate_figures(plot_quantity, file_loc, report, reportname, surface,
     801                     leg_label, f_list,
     802                     gauges, locations, elev, production_dirs, time_min, time_max,
    793803                     title_on, label_id, verbose):
    794804
    795805    from math import sqrt, atan, degrees
    796     from Numeric import ones, allclose, zeros, Float
     806    from Numeric import ones, allclose, zeros, Float, ravel
    797807    from os import sep, altsep, getcwd, mkdir, access, F_OK, environ
    798808    from pylab import ion, hold, plot, axis, figure, legend, savefig, \
    799809         xlabel, ylabel, title, close, subplot
     810
     811    import pylab as p1
     812    import mpl3d.mplot3d as p3
    800813
    801814    if report == True:   
    … …  
    837850    bearings = zeros((n0,m,p), Float)
    838851    depths = zeros((n0,m,p), Float)
     852    eastings = zeros((n0,m,p), Float)
    839853    min_stages = []
    840854    max_stages = []
    … …  
    842856    max_speeds = []
    843857    c = 0
     858    model_time_plot3d = zeros((n0,m), Float)
     859    stages_plot3d = zeros((n0,m), Float)
     860    eastings_plot3d = zeros((n0,m),Float)
    844861    ##### loop over each swwfile #####
    845862    for j, f in enumerate(f_list):
    … …  
    858875            fid_out.write(s)
    859876            #### generate quantities #######
    860            
    861877            for i, t in enumerate(f.get_time()):
    862878                if time_min <= t <= time_max:
    … …  
    872888                        vel = m / (depth + 1.e-30)
    873889                    bearing = calc_bearing(uh, vh)
    874                     model_time[i,k,j] = t/60.0         
    875                     stages[i,k,j] = w 
     890                    model_time[i,k,j] = t/60.0
     891                    stages[i,k,j] = w
    876892                    elevations[i,k,j] = z
    877893                    xmom[i,k,j] = uh
    … …  
    880896                    speed[i,k,j] = vel
    881897                    bearings[i,k,j] = bearing
    882                     depths[i,k,j] = depth
     898                    depths[i,k,j] = depth
     899                    thisgauge = gauges[k]
     900                    eastings[i,k,j] = thisgauge[0]
    883901                    s = '%.2f, %.2f, %.2f, %.2f\n' %(t, w, m, vel)
    884902                    fid_out.write(s)
    … …  
    886904                    if w < min_stage: min_stage = w
    887905                    if m > max_momentum: max_momentum = m
    888                     if vel > max_speed: max_speed = vel
    889 
     906                    if vel > max_speed: max_speed = vel           
     907           
    890908            max_stages.append(max_stage)
    891909            min_stages.append(min_stage)
    892910            max_momentums.append(max_momentum)
    893911            max_speeds.append(max_speed)       
     912            #### finished generating quantities for each swwfile #####
     913
     914        model_time_plot3d[:,:] = model_time[:,:,j]
     915        stages_plot3d[:,:] = stages[:,:,j]
     916        eastings_plot3d[:,] = eastings[:,:,j]
    894917           
    895             #### finished generating quantities for each swwfile #####
    896            
     918        if surface ==  True:
     919            print 'Printing surface figure'
     920            for i in range(2):
     921                fig = p1.figure(10)
     922                ax = p3.Axes3D(fig)
     923                if len(gauges) > 80:
     924                    ax.plot_surface(eastings,model_time_plot3d,stages_plot3d)
     925                    #ax.plot_surface(eastings[:,:,j],model_time[:,:,j],stages[:,:,j])
     926                else:
     927                    #ax.plot_wireframe(eastings[:,:,j],model_time[:,:,j],stages[:,:,j])
     928                    ax.plot_wireframe(eastings,model_time_plot3d,stages_plot3d)
     929                    #ax.plot3D(ravel(eastings[:,:,j]),ravel(model_time[:,:,j]),ravel(stages[:,:,j]))
     930                ax.set_xlabel('x')
     931                ax.set_ylabel('time')
     932                ax.set_zlabel('stage')
     933                fig.add_axes(ax)
     934                p1.show()
     935                surfacefig = 'solution_surface%s' %leg_label[j]
     936                p1.savefig(surfacefig)
     937                p1.close()
     938       
    897939    #### finished generating quantities for all swwfiles #####
    898            
     940
    899941    stage_axis = axis([time_min/60.0, time_max/60.0, min(min_stages), max(max_stages)*1.1])
    900942    vel_axis = axis([time_min/60.0, time_max/60.0, min(max_speeds), max(max_speeds)*1.1])
    … …  
    9751017                        mkdir (figdir)
    9761018                    latex_file_loc = figdir.replace(sep,altsep)
    977                     # storing files in production directory
    978                     graphname_latex = '%sgauge%s%s' %(latex_file_loc, gaugeloc2, which_quantity)     
    979                     # giving location in latex output file
    980                     graphname_report_input = '%sgauge%s%s' %('..'+altsep+'report_figures'+altsep, gaugeloc2, which_quantity)
     1019                    graphname_latex = '%sgauge%s%s' %(latex_file_loc, gaugeloc2, which_quantity) # storing files in production directory   
     1020                    graphname_report_input = '%sgauge%s%s' %('..'+altsep+'report_figures'+altsep, gaugeloc2, which_quantity) # giving location in latex output file
    9811021                    graphname_report.append(graphname_report_input)
    9821022                   
    … …  
    9911031
    9921032                    if len(label_id) == 1:
    993                         # storing files in production directory
    994                         graphname_latex = '%sgauge%s%s%s' %(latex_file_loc, gaugeloc2, which_quantity, label_id2)     
    995                         # giving location in latex output file
    996                         graphname_report = '%sgauge%s%s%s' %('..'+altsep+'report_figures'+altsep, gaugeloc2, which_quantity, label_id2)
     1033                        graphname_latex = '%sgauge%s%s%s' %(latex_file_loc, gaugeloc2, which_quantity, label_id2) # storing files in production directory 
     1034                        graphname_report = '%sgauge%s%s%s' %('..'+altsep+'report_figures'+altsep, gaugeloc2, which_quantity, label_id2) # giving location in latex output file
    9971035                        s = '\includegraphics[width=0.49\linewidth, height=50mm]{%s%s}' %(graphname_report, '.png')
    9981036                        fid.write(s)

  • production/onslow_2006/MOST_timeseries.py

    r3188 r3190  
    6161
    6262# Start report generation
    63 input_name = project.comparereportdir + sep + 'comparison_onslow.tex'
     63input_name = project.comparereportdir + sep + 'MOST_ANUGA_comparison_onslow.tex'
    6464fid = open(input_name, 'w')
    6565

  • production/onslow_2006/compare_timeseries.py

    r3188 r3190  
    2323# User defined inputs
    2424production_dirs = {'20060515_001733': '100m boundary',
    25                    '20060530_102753': '50m boundary'}
     25                   '20060530_102753': '50m boundary',
     26                   'MOST': 'MOST'}
    2627
    2728gauge_map = 'onslow_boundary_gauges.png'
    … …  
    3637    file_loc = project.outputdir + label_id + sep
    3738    swwfile = file_loc + project.basename + '.sww'
     39    if label_id == 'MOST':
     40        swwfile = project.boundarydir + project.boundary_basename + '.sww'
    3841    swwfiles[swwfile] = label_id
    3942
    … …  
    4245                              production_dirs,
    4346                              report = True,
    44                               reportname = 'latexoutput_compare',
     47                              reportname = 'latexoutput_compare50_100_MOST',
    4548                              surface = False,
    4649                              plot_quantity = plot_quantity,
    … …  
    5760
    5861# Start report generation
    59 input_name = project.comparereportdir + sep + 'comparison_onslow.tex'
     62input_name = project.comparereportdir + sep + '50100MOSTcomparison_onslow.tex'
    6063fid = open(input_name, 'w')
    6164

  • production/onslow_2006/make_report.py

    r3188 r3190  
    2020The report structure is
    2121
     22* Executive Summary
    2223* Introduction
    2324* Modelling Methodology
    … …  
    159160\\begin{document}
    160161  \maketitle
    161  
    162   \\begin{abstract}
    163     \input{abstract}
    164   \end{abstract}
    165  
     162   
    166163  \\tableofcontents
     164
     165    \section{Executive Summary}
     166    \label{sec:execsum}
     167  \input{execsum}
    167168 
    168169  \section{Introduction}

  • production/pt_hedland_2006/make_report.py

    r3094 r3190  
    2020The report structure is
    2121
     22* Executive Summary
    2223* Introduction
     24* Modelling Methodology
    2325* Tsunami scenario
    2426* Inundation model
    … …  
    6870report_title = 'Tsunami impact modelling for the North West shelf: %s' %scenario_name.title()
    6971
    70 production_dirs = {'': '1.5 AHD',
    71                    '': '-1.5 AHD',
     72production_dirs = {'': '3.6 AHD',
     73                   '': '-3.9 AHD',
    7274                   '': '0 AHD'}
    7375
    74 max_maps = {'1.5 AHD': 'HAT_map',
    75             '-1.5 AHD': 'LAT_map',
     76max_maps = {'3.6 AHD': 'HAT_map',
     77            '-3.9 AHD': 'LAT_map',
    7678            '0 AHD': 'MSL_map'}
    77 
    78 damage_maps = {'1.5 AHD': 'HAT_damage',
    79                '-1.5 AHD': 'LAT_damage',
    80                '0 AHD': 'MSL_damage'}
    8179
    8280gauge_map = 'pt_hedland_gauge_map.jpg'
    … …  
    9189    swwfiles[swwfile] = label_id
    9290
    93 texname = sww2timeseries(swwfiles,
    94                          project.gauge_filename,
    95                          #label_id,
    96                          production_dirs,
    97                          report = True,
    98                          plot_quantity = ['stage', 'speed'],
    99                          time_min = None,
    100                          time_max = None,
    101                          title_on = False,
    102                          verbose = True)
     91texname, elev_output = sww2timeseries(swwfiles,
     92                                      project.gauge_filename,
     93                                      production_dirs,
     94                                      report = True,
     95                                      reportname = 'latexoutput',
     96                                      plot_quantity = ['stage', 'speed'],
     97                                      surface = False,
     98                                      time_min = None,
     99                                      time_max = None,
     100                                      title_on = False,
     101                                      verbose = True)
    103102
    104103latex_output.append(texname)
    … …  
    117116% * an introduction must be written in introduction.tex; a basic outline and
    118117%   some of the core inputs are already in place
     118% * outline of the modelling methodology provided in modelling_methodology.tex
    119119% * the tsunami-genic event should be discussed in tsunami_scenario.tex
    120120% * an computational_setup.tex file needs to be written for the particular scenario
    … …  
    157157  \maketitle
    158158 
    159   \\begin{abstract}
    160     \input{abstract}
    161   \end{abstract}
    162  
    163159  \\tableofcontents
     160
     161   \section{Executive Summary}
     162    \label{sec:execsum}
     163  \input{execsum}
    164164 
    165165  \section{Introduction}
    166166    \label{sec:intro}
    167167  \input{introduction}
    168    
     168
     169   \section{Modelling methodology}
     170    \label{sec:methodology}
     171    \input{modelling_methodology}
     172   
    169173  \section{Tsunami scenarios}
    170174    \label{sec:tsunamiscenario}
    171175    \input{tsunami_scenario}
    172176
    173    \section{Inundation Model}
     177   \section{Inundation model}
    174178    \label{sec:anuga}
    175179    \input{anuga}
    … …  
    213217
    214218s  = """
    215 \caption{Point locations used for Pt Hedland study.} 
     219\caption{Point locations used for Onslow study.} 
    216220\label{fig:points}
    217221\end{figure}
    … …  
    235239
    236240s = """
    237    \section{Damage modelling}
    238     \label{sec:damage}
     241   \section{Impact modelling}
     242    \label{sec:impact}
    239243     \input{damage}
    240244"""
    241245fid.write(s)
    242246
    243 for i, name in enumerate(production_dirs.keys()):
    244 
    245     s = '\input{%s} \n \clearpage \n \n' %damage_maps[production_dirs[name]]
    246     fid.write(s)
     247#for i, name in enumerate(production_dirs.keys()):
     248
     249#    s = '\input{%s} \n \clearpage \n \n' %damage_maps[production_dirs[name]]
     250#    fid.write(s)
    247251
    248252s = """

  • production/pt_hedland_2006/project.py

    r3094 r3190  
    108108# region to export (used from export_results.py)
    109109
    110 e_min_area = 659000#633000
    111 e_max_area = 678000#690000
    112 n_min_area = 7746000#7740000
    113 n_max_area = 7757000#7761000
     110e_min_area = 633000
     111e_max_area = 690000
     112n_min_area = 7740000
     113n_max_area = 7761000
    114114
    115115refzone = 50 # confirm with Hamish

  • production/sydney_2006/export_results.py

    r2875 r3190  
    44from pyvolution.ermapper_grids import read_ermapper_grid
    55
    6 name = project.outputname2
     6name = project.outputname
    77
    88#print 'Which variable do you want to export?'

  • production/sydney_2006/find_max.py

    r2878 r3190  
    3535c = 0
    3636direction = 'negative'
     37xstar = x0
     38print 'start of loop', deriv
    3739while c < 100000 and deriv > 0:
    3840    deriv = dfuncdx(x,0.83)
    39     if deriv < 0: xstar = x
     41    if deriv < 0:
     42        xstar = x
     43        print 'hello', deriv, xstar/lam0, c
    4044    if direction == 'positive': x += step
    4145    if direction == 'negative': x -= step
    4246    c += 1
     47   
    4348
    4449print 'location of maximum of surface elevation function: xstar = %f' % (xstar/lam0)
    4550const = 1.0 #a3D = ? #sydney = 86? and kappad=1
    4651
    47 x = arange(-20,25,0.001)
    48 y = arange(-10,10,0.1)
    49 X,Y = meshgrid(x,y)
     52#x = arange(-20,25,0.001)
     53#y = arange(-10,10,0.1)
     54#X,Y = meshgrid(x,y)
    5055
    51 test = 0.0
    52 for xi in x:
    53     testi = func(xi,0.0,0.83)
    54     if direction == 'positive':
    55         if testi > test:
    56             test = testi
    57             xstar = xi
    58     else:
    59         if testi < test:
    60             test = testi
    61             xstar = xi
     56#test = 0.0
     57#for xi in x:
     58#    testi = func(xi,0.0,0.83)
     59#    if direction == 'positive':
     60#        if testi > test:
     61#            test = testi
     62#            xstar = xi
     63#    else:
     64#        if testi < test:
     65#            test = testi
     66#            xstar = xi
    6267   
    63 print 'check: xstar = %f and eta(xstar) = %f' %(xstar/lam0, test)
     68#print 'check: xstar = %f and eta(xstar) = %f' %(xstar/lam0, test)

  • production/sydney_2006/get_timeseries.py

    r2885 r3190  
    55import project
    66from pyvolution.util import file_function
    7 from coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
     7#from coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
    88from pylab import *
    99from matplotlib.ticker import MultipleLocator, FormatStrFormatter
    1010
    11 swwfile = project.outputname + '.sww'
    12 
    13 
     11swwfile = project.outputname + '.sww'
     12#swwfile = project.outputdir + timestampdir + sep + project.outputname + '.sww'
     13# place to store figures
     14#graphloc = project.outputdir + timestampdir + sep
     15graphloc = project.outputdir
    1416
    1517#Time interval to plot
    16 tmin = 13000
    17 tmax = 21000
     18tmin = 2*60
     19tmax = 10*60
    1820
    1921def get_gauges_from_file(filename):
    … …  
    3436        gaugelocation.append(location)
    3537       
    36     #Return gauges and raw data for subsequent storage
    3738    return gauges, lines, gaugelocation
    3839
    … …  
    4041gauges, lines, locations = get_gauges_from_file(project.gauge_filename)
    4142
    42 print 'number of gauges for Benfield:   ', len(gauges)
     43print 'number of gauges:   ', len(gauges)
    4344
    4445#Read model output
    … …  
    5051                  use_cache = True)
    5152
    52 
    53                
     53print 'size f', size(f.quantities['stage'],axis=0), size(f.quantities['stage'],axis=1)
     54
    5455from math import sqrt, atan, degrees
    5556from Numeric import ones
    … …  
    8586        uh = f(t, point_id = k)[2]
    8687        vh = f(t, point_id = k)[3]
    87         myloc = locations[k]
     88        gaugeloc = locations[k]
    8889        depth = w-z
    89 
     90       
    9091        m = sqrt(uh*uh + vh*vh)   #Absolute momentum
    9192        vel = sqrt(uh*uh + vh*vh) / (w-z + 1.e-30) #Absolute velocity
    … …  
    135136             model_time, elevations, '-k')
    136137    #name = 'Gauge_%d: (%.1f, %.1f)' %(k, g[0], g[1])
    137     name = 'Gauge_%d: (%.1f, %.1f) Location: %s' %(k, g[0], g[1], myloc)
     138    name = 'Gauge_%d: (%.1f, %.1f) Location: %s' %(k, g[0], g[1], gaugeloc)
    138139    title(name)
    139140
    … …  
    144145           shadow=True,
    145146           loc='upper right')
    146     #savefig('Gauge_%d_stage' %k)
    147     savefig('Gauge_%s_stage' %myloc)
    148 
    149     # raw_input('Next')
     147    #('Gauge_%d_stage' %k)
     148    savefig('%sGauge_%s_stage' %(graphloc, gaugeloc))
     149    #savefig('Gauge_%s_stage.eps' %gaugeloc)
    150150   
    151151    #Momentum plot
    … …  
    159159    ylabel('sqrt( uh^2 + vh^2 ) [m^2/s]')   
    160160    #savefig('Gauge_%d_momentum' %k)
    161     savefig('Gauge_%s_momentum' %myloc)
     161    savefig('%sGauge_%s_momentum' %(graphloc, gaugeloc))
    162162   
    163     # raw_input('Next')
    164 
    165163    #Bearing plot
    166164    ion()
    … …  
    181179    ylabel(' atan(vh/uh) [degrees from North]')   
    182180    #savefig('Gauge_%d_bearing' %k)
    183     savefig('Gauge_%s_bearing' %myloc)
    184    
    185     # raw_input('Next')
     181    savefig('%sGauge_%s_bearing' %(graphloc, gaugeloc))
    186182
    187183    #Speed plot
    … …  
    195191    ylabel('sqrt( uh^2 + vh^2 ) / depth [m/s]')   
    196192    #savefig('Gauge_%d_speed' %k)
    197     savefig('Gauge_%s_speed' %myloc)
    198    
    199     # raw_input('Next')
    200 
    201     whichone = '_%s' %myloc
     193    savefig('%sGauge_%s_speed' %(graphloc, gaugeloc))
     194
     195    whichone = '_%s' %gaugeloc
    202196    thisfile = project.gaugetimeseries+whichone+'.csv'
    203197    fid = open(thisfile, 'w')

  • production/sydney_2006/plot_integral.py

    r2487 r3190  
    3030hold(False)
    3131plot(t, integral, '.-b')
     32axis([0, max(t), min(integral)-0.1, min(integral) + 1])
    3233title('Checking for water conservation')
    3334xlabel('Time (sec)')

  • production/sydney_2006/process_gauges.py

    r2885 r3190  
    55import project
    66from pyvolution.util import file_function
    7 from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
     7#from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees, redfearn
    88from pylab import *
    99
    … …  
    1111plot = False
    1212   
    13 swwfile = project.outputname + '.sww'
     13swwfile = project.outputname4 + '.sww'
    1414
    1515def get_gauges_from_file(filename):

  • production/sydney_2006/project.py

    r2875 r3190  
    3030nmaxviz = 6283000
    3131
    32 basename = 'slump_315res'
    33 basename2 = 'slump_1000res'
     32basename = 'slump_duncan'
     33basename2 = 'slump_1000res_1min'
    3434basename4 = 'slump_poly_ingo_test'
    3535
    … …  
    6262
    6363gauge_filename = gaugedir + 'sydney_gauges.xya'
     64buildings_filename = gaugedir + 'all_bld_ind.csv'
    6465#gauge_filename = gaugedir + 'sydney_gauges_test.xya'
    6566gauge_outname = gaugedir + 'gauges_max_output.xya'
    … …  
    137138
    138139# clipping used for demo purposes to fit around poly3 and poly4 (Ingo's files)
    139 j0 = [318000, 6249000]
    140 j1 = [387000, 6249000]
    141 j2 = [387000, 6270000]
    142 j3 = [318000, 6270000]
    143 
    144 demopoly = [j0, j1, j2, j3]
     140#j0 = [318000, 6249000]
     141#j1 = [387000, 6249000]
     142#j2 = [387000, 6270000]
     143#j3 = [318000, 6270000]
     144
     145#demopoly = [j0, j1, j2, j3]
     146
     147# second demo poly which isn't rectangular
     148j0 = [385000, 6280000]
     149j1 = [360000, 6272500]
     150j2 = [335000, 6272500]
     151j3 = [330000, 6265000]
     152j31 = [325000, 6260000]
     153j4 = [316000, 6260000]
     154j5 = [316000, 6247000]
     155j6 = [350000, 6247000]
     156j7 = [385000, 6238000]
     157
     158demopoly = [j0, j1, j2, j3, j31, j4, j5, j6, j7]
    145159
    146160# to put chunk back in

  • production/sydney_2006/run_sydney_smf.py

    r2640 r3190  
    2727from pyvolution.quantity import Quantity
    2828from Numeric import allclose
     29from utilities.polygon import inside_polygon
    2930
    3031# Application specific imports
    … …  
    7071
    7172# original issue to Benfield
    72 #interior_res = 5000
    73 #interior_regions = [[project.harbour_polygon_2, interior_res],
    74 #                    [project.botanybay_polygon_2, interior_res]]
     73interior_res = 5000
     74interior_regions = [[project.harbour_polygon_2, interior_res],
     75                    [project.botanybay_polygon_2, interior_res]]
    7576
    7677# used for finer mesh
    77 interior_res1 = 5000
    78 interior_res2 = 315
    79 interior_regions = [[project.newpoly1, interior_res1],
    80                     [project.south1, interior_res1],
    81                     [project.finepolymanly, interior_res2],
    82                     [project.finepolyquay, interior_res2]]
    83 
    84 print 'number of interior regions', len(interior_regions)
     78#interior_res1 = 5000
     79#interior_res2 = 315
     80#interior_regions = [[project.newpoly1, interior_res1],
     81#                    [project.south1, interior_res1],
     82#                    [project.finepolymanly, interior_res2],
     83#                    [project.finepolyquay, interior_res2]]
     84
     85# used for coastal polygons constructed by GIS guys
     86def get_polygon_from_file(filename):
     87    """ Function to read in output from GIS determined polygon
     88    """
     89    fid = open(filename)
     90    lines = fid.readlines()
     91    fid.close()
     92   
     93    polygon = []
     94    for line in lines[1:]:
     95       fields = line.split(',')
     96       x = float(fields[1])
     97       y = float(fields[2])
     98       polygon.append([x, y])
     99       
     100    return polygon
     101
     102num_polygons = 9
     103fileext = '.csv'
     104filename = project.polygonptsfile
     105
     106#interior_res = 1000
     107#interior_regions = []
     108#bounding_polygon = project.diffpolygonall#project.demopoly
     109#count = 0
     110#for p in range(1, num_polygons+1):
     111#    thefilename = filename + str(p) + fileext
     112#    print 'reading in polygon points', thefilename
     113#    interior_polygon = get_polygon_from_file(thefilename)
     114#    interior_regions.append([interior_polygon, interior_res])
     115#    n = len(interior_polygon)
     116#    # check interior polygon falls in bounding polygon
     117#    if len(inside_polygon(interior_polygon, bounding_polygon,
     118#                   closed = True, verbose = False)) <> len(interior_polygon):
     119#        print 'WARNING: interior polygon %d is outside bounding polygon' %(p)
     120#        count += 1
     121#    # check for duplicate points in interior polygon
     122   
     123print 'number of interior polygons: ', len(interior_regions)
     124#if count == 0: print 'interior polygons OK'
    85125
    86126#FIXME: Fix caching of this one once the mesh_interface is ready
    … …  
    89129# original + refined region
    90130_ = cache(create_mesh_from_regions,
    91           project.diffpolygonall,
     131          #project.demopoly,
     132          project.diffpolygonall2,
     133          #{'boundary_tags': {'bottom': [0],
     134          #                   'right1': [1], 'right0': [2],
     135          #                   'right2': [3], 'top': [4], 'left1': [5],
     136          #                   'left2': [6], 'left3': [7]},
    92137          {'boundary_tags': {'bottom': [0],
    93                              'right1': [1], 'right0': [2],
    94                              'right2': [3], 'top': [4], 'left1': [5],
     138                             'bottom1': [1], 'right': [2],
     139                             'top1': [3], 'top': [4], 'left1': [5],
    95140                             'left2': [6], 'left3': [7]},
     141          #{'boundary_tags': {'bottom': [0], 'right': [1],
     142          #                   'top': [2], 'left': [3]},
    96143           'maximum_triangle_area': 100000,
    97144           'filename': meshname,           
    … …  
    126173                               radius=3330,
    127174                               dphi=0.23,
    128                                x0=project.slump_origin[0],
    129                                y0=project.slump_origin[1],
     175                               x0=project.slump_origin2[0],
     176                               y0=project.slump_origin2[1],
    130177                               alpha=0.0,
    131                                domain=domain)
    132 
     178                               domain=domain,
     179                               verbose=True)
    133180
    134181#-------------------------------------------------------------------------------                                 
    … …  
    138185# apply slump after 30 mins, initialise to water level of tide = 0
    139186domain.set_quantity('stage', 0.0)
    140 domain.set_quantity('friction', 0.03)
     187domain.set_quantity('friction', 0.04)
    141188domain.set_quantity('elevation',
    142189                    filename = project.combineddemname + '.pts',
    … …  
    158205#                      'right2': Br, 'top': Br, 'left1': Br,
    159206#                      'left2': Br, 'left3': Br} )
    160 
     207# for new tests 4 April 2006
     208#domain.set_boundary( {'bottom': Br, 'bottom1': Br, 'right': Br,
     209#                      'top1': Br, 'top': Br, 'left1': Br,
     210#                      'left2': Br, 'left3': Br} )
     211                             
    161212# enforce Dirichlet BC - from 30/03/06 Benfield visit
    162 domain.set_boundary( {'bottom': Bd, 'right1': Bd, 'right0': Bd,
    163                       'right2': Bd, 'top': Bd, 'left1': Bd,
     213domain.set_boundary( {'bottom': Bd, 'bottom1': Bd, 'right': Bd,
     214                      'top1': Bd, 'top': Bd, 'left1': Bd,
    164215                      'left2': Bd, 'left3': Bd} )
    165216
    166217# increasingly finer interior regions
    167 #domain.set_boundary( {'bottom': Br, 'right': Br, 'left': Br, 'top': Br} )
     218#domain.set_boundary( {'bottom': Bd, 'right': Bd, 'left': Bd, 'top': Bd} )
    168219
    169220
    … …  
    177228fid = open(thisfile, 'w')
    178229
    179 for t in domain.evolve(yieldstep = 120, finaltime = 18000):
     230# save every 10 secs leading up to slump initiation
     231for t in domain.evolve(yieldstep = 10, finaltime = 60): # 6 steps
    180232    domain.write_time()
    181233    domain.write_boundary_statistics(tags = 'bottom')
    182 
    183     # calculate integral
    184     thisstagestep = domain.get_quantity('stage')
    185     s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
    186     fid.write(s)
    187 
    188     # add slump after 30 mins
    189     if allclose(t, 30*60):
     234    # calculate integral
     235    thisstagestep = domain.get_quantity('stage')
     236    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
     237    fid.write(s)
     238    # add slump after 30 secs
     239    if allclose(t, 30):
    190240        slump = Quantity(domain)
    191241        slump.set_values(tsunami_source)
    192242        domain.set_quantity('stage', slump + thisstagestep)
    193    
     243        #test_stage = domain.get_quantity('stage')
     244        #test_elevation = domain.get_quantity('elevation')
     245        #test_depth = test_stage - test_elevation
     246        #test_max = max(test_depth.get_values())
     247        #print 'testing', test_max
     248
     249import sys; sys.exit()
     250
     251# save every two minutes leading up to interesting period
     252for t in domain.evolve(yieldstep = 120, finaltime = 660, # steps
     253                       skip_initial_step = True):
     254    domain.write_time()
     255    domain.write_boundary_statistics(tags = 'bottom')
     256    # calculate integral
     257    thisstagestep = domain.get_quantity('stage')
     258    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
     259    fid.write(s)
     260
     261
     262# save every thirty secs during interesting period
     263for t in domain.evolve(yieldstep = 60, finaltime = 5000, # steps
     264                       skip_initial_step = True):
     265    domain.write_time()
     266    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage')
     267    # calculate integral
     268    thisstagestep = domain.get_quantity('stage')
     269    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
     270    fid.write(s)
     271
     272
     273import sys; sys.exit()
     274# save every two mins for next 5000 secs
     275for t in domain.evolve(yieldstep = 120, finaltime = 10000, # about 42 steps
     276                       skip_initial_step = True):
     277    domain.write_time()
     278    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage')
     279    # calculate integral
     280    thisstagestep = domain.get_quantity('stage')
     281    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
     282    fid.write(s)
     283   
     284# save every half hour to end of simulation   
     285for t in domain.evolve(yieldstep = 1800, finaltime = 10*60*60, # 14 steps
     286                       skip_initial_step = True):
     287    domain.write_time()
     288    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage'
     289    # calculate integral
     290    thisstagestep = domain.get_quantity('stage')
     291    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
     292    fid.write(s)
    194293   
    195294print 'That took %.2f seconds' %(time.time()-t0)

  • production/sydney_2006/run_sydney_smf_polyingo.py

    r2732 r3190  
    2020
    2121# Related major packages
    22 from pyvolution.shallow_water import Domain, Reflective_boundary, Time_boundary
     22from pyvolution.shallow_water import Domain, Reflective_boundary, Time_boundary, Dirichlet_boundary
    2323from pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
    2424from pyvolution.combine_pts import combine_rectangular_points_files
    … …  
    2626from pyvolution.quantity import Quantity
    2727from geospatial_data import *
     28from utilities.polygon import inside_polygon
     29from Numeric import allclose
    2830
    2931# Application specific imports
    … …  
    8183meshname = project.meshname4+'.msh'
    8284
    83 interior_res1 = 25000
     85interior_res1 = 1000
    8486interior_res2 = 1315
    85 print 'interior resolution', interior_res1, interior_res2
    8687
    8788# Read in files containing polygon points (generated by GIS)
    … …  
    104105    return polygon
    105106
    106 num_polygons = 3
     107num_polygons = 4
    107108fileext = '.csv'
    108109filename = project.polygonptsfile
    109110
    110111interior_regions = []
    111 
     112bounding_polygon = project.diffpolygonall
     113count = 0
    112114for p in range(3, num_polygons+1):
    113115    thefilename = filename + str(p) + fileext
    … …  
    115117    interior_polygon = get_polygon_from_file(thefilename)
    116118    interior_regions.append([interior_polygon, interior_res1])
    117 
     119    n = len(interior_polygon)
     120    # check interior polygon falls in bounding polygon
     121    if len(inside_polygon(interior_polygon, bounding_polygon,
     122                   closed = True, verbose = False)) <> len(interior_polygon):
     123        print 'WARNING: interior polygon %d is outside bounding polygon' %(p)
     124        count += 1
     125    # check for duplicate points in interior polygon
     126   
    118127print 'number of interior polygons: ', len(interior_regions)
    119 
    120 #print 'Number of interior polygons: ', len(interior_regions)
     128if count == 0: print 'interior polygons OK'
    121129
    122130# original
    … …  
    132140#                    [project.finepolyquay, interior_res2]]
    133141
     142
    134143#FIXME: Fix caching of this one once the mesh_interface is ready
    135144from caching import cache
     145
    136146
    137147#_ = cache(create_mesh_from_regions,
    … …  
    141151#                             'right2': [3], 'top': [4], 'left1': [5],
    142152#                             'left2': [6], 'left3': [7]},
    143 #           'maximum_triangle_area': 100000,
     153#           'maximum_triangle_area': 150000,
    144154#           'filename': meshname,           
    145155#           'interior_regions': interior_regions},
    146156#          verbose = True)
    147157
     158# for demo construction
    148159_ = cache(create_mesh_from_regions,
    149           project.diffpolygonall_test,
     160          project.demopoly,
    150161          {'boundary_tags': {'bottom': [0], 'right': [1],
    151162                             'top': [2], 'left': [3]},
    … …  
    171182domain.set_datadir(project.outputdir)
    172183domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
    173 
    174 
    175184
    176185
    … …  
    190199                               domain=domain)
    191200
    192 #print 'testing', tsunami_source.get_integral()
    193 
    194201#-------------------------------------------------------------------------------                                 
    195202# Setup initial conditions
    … …  
    199206domain.set_quantity('stage', 0) #tsunami_source)
    200207domain.set_quantity('friction', 0.0) # supplied by Benfield, initial value 0.03
    201 domain.set_quantity('elevation', alpha = 10.0,#G, alpha = 550.0,
    202                     #filename = project.combineddemname + '.pts',
     208domain.set_quantity('elevation',
    203209                    filename = project.coarsedemname + '.pts',
    204                     use_cache = False,
     210                    use_cache = True,
    205211                    verbose = True)
    206212
    … …  
    231237
    232238Br = Reflective_boundary(domain)
    233 
    234 #domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
    235 #                      'right2': Br, 'top': Br, 'left1': Br,
    236 #                      'left2': Br, 'left3': Br} )
    237 domain.set_boundary( {'bottom': Br, 'right': Br, 'left': Br, 'top': Br} )
     239Bd = Dirichlet_boundary([0, 0, 0])
     240
     241# for demo
     242domain.set_boundary( {'bottom': Bd, 'right': Bd, 'left': Bd, 'top': Bd} )
    238243#domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
    239244#                      'right2': Br, 'top': Br, 'left1': Br,
    … …  
    250255fid = open(thisfile, 'w')
    251256   
    252 # original
    253 for t in domain.evolve(yieldstep = 1, finaltime = 10):
     257for t in domain.evolve(yieldstep = 10, finaltime = 600):
    254258    domain.write_time()
    255259    domain.write_boundary_statistics(tags = 'bottom')
    … …  
    257261    s = '%.2f, %.2f\n' %(t, stagestep.get_integral())
    258262    fid.write(s)
     263    if allclose(t, 30):
     264        slump = Quantity(domain)
     265        slump.set_values(tsunami_source)
     266        domain.set_quantity('stage', slump + stagestep)
     267
     268# save every two minutes leading up to interesting period
     269for t in domain.evolve(yieldstep = 120, finaltime = 660, # steps
     270                       skip_initial_step = True):
     271    domain.write_time()
     272    domain.write_boundary_statistics(tags = 'bottom')
     273    # calculate integral
     274    thisstagestep = domain.get_quantity('stage')
     275    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
     276    fid.write(s)
     277   
     278# save every thirty secs during interesting period
     279for t in domain.evolve(yieldstep = 30, finaltime = 3500, # steps
     280                       skip_initial_step = True):
     281    domain.write_time()
     282    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage')
     283    # calculate integral
     284    thisstagestep = domain.get_quantity('stage')
     285    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
     286    fid.write(s)
     287
     288# save every two mins for next 5000 secs
     289for t in domain.evolve(yieldstep = 120, finaltime = 10000, # about 42 steps
     290                       skip_initial_step = True):
     291    domain.write_time()
     292    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage')
     293    # calculate integral
     294    thisstagestep = domain.get_quantity('stage')
     295    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
     296    fid.write(s)
     297   
     298# save every half hour to end of simulation   
     299for t in domain.evolve(yieldstep = 1800, finaltime = 10*60*60, # 14 steps
     300                       skip_initial_step = True):
     301    domain.write_time()
     302    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage'
     303    # calculate integral
     304    thisstagestep = domain.get_quantity('stage')
     305    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
     306    fid.write(s)
    259307
    260308print 'That took %.2f seconds' %(time.time()-t0)

  • production/sydney_2006/run_timeseries.py

    r2795 r3190  
    1414
    1515# False to generate latex output, True otherwise
    16 title_on = False
     16title_on = True
    1717
    1818# generate figures - see sww2timeseries for documentation

  • production/sydney_2006/smf_volume_conservation.py

    r2673 r3190  
    77from Numeric import ones
    88from pylab import *
    9 import scipy
    10 from scipy.special import erf
     9#import scipy
     10#from scipy.special import erf
    1111
    1212ion()
    … …  
    3939eta2 = func(X,0,1)
    4040
     41#import sys; sys.exit()
    4142#eta = func(X,Y,kappad)
    4243#im1 = imshow(eta, cmap=cm.jet, interpolation='nearest')
Note: See TracChangeset for help on using the changeset viewer.