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

MOST and ANUGA comparisons and updates

Location:
production/sydney_2006
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • 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.