Changeset 5370 for anuga_work


Ignore:
Timestamp:
May 27, 2008, 5:39:41 PM (17 years ago)
Author:
duncan
Message:

Hinwood scenario

Location:
anuga_work/development
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/development/Hinwood_2008/create_mesh.py

    r5092 r5370  
    1 """Create mesh for University of Queensland dam break flume
     1"""Create mesh for Jon Hinwoods wave flume
    22"""
    33
    … …  
    66from anuga.coordinate_transforms.geo_reference import Geo_reference
    77
    8 xslope = 5.0 # Distance between the boundary ande the start of the slope
     8xslope = 5.0 # Distance between the boundary and the start of the slope
    99
    1010
    11 def generate(mesh_filename, maximum_triangle_area=0.01):
     11def generate(mesh_filename, slope, maximum_triangle_area=0.01):
    1212    """
    13     Generate mesh for University of Aberbeen dam break flume.
    14     The gate position is the distance from the right wall of the wave tank
    15     to the gate.
     13    Generate mesh for Monash University wave generation flume.
     14   
     15    Origin of x coordinate is the toe of the beach, x measured
     16    positive shorewards.
    1617   
    1718    """
    … …  
    1920    global xslope
    2021   
    21     xright  = 19.0
     22    xleft = slope['xleft'][0]
     23    xtoe = slope['xtoe'][0]
     24    xbeach = slope['xbeach'][0]
     25    xright = slope['xright'][0]
    2226    ybottom = 0
    2327    ytop    = 1.00
    24     xleft = 0.0
    2528    ###xslope = slope
    2629
    … …  
    3235
    3336
    34     # slope seperation (middle)
    35     point_slope_top = [xslope, ytop]
    36     point_slope_bottom = [xslope, ybottom]   
     37    # 1st slope seperation (middle)
     38    point_toe_top = [xtoe, ytop]
     39    point_toe_bottom = [xtoe, ybottom]   
    3740
     41    # 2nd slope seperation (middle)
     42    point_beach_top = [xbeach, ytop]
     43    point_beach_bottom = [xbeach, ybottom]
     44   
    3845    m = Mesh()
    3946
    … …  
    4148    points = [point_sw,   #se
    4249              point_nw,
    43               point_slope_top,
     50              point_toe_top, # 2
     51              point_beach_top, # 3
    4452              point_ne,
    4553              point_se,
    46               point_slope_bottom
     54              point_beach_bottom, # 6
     55              point_toe_bottom # 7
    4756              ]
    4857   
    … …  
    5362        [3,4],
    5463        [4,5],
    55         [5,0],  #The outer border
    56         [2,5]       # slope separator
     64        [5,6],
     65        [6,7],
     66        [7,0],  #The outer border
     67        [2,7],       # toe separator
     68        [3,6]       # _beach separator
    5769        ]   
    5870   
    59     segment_tags = {'wall':[1,2,3,4,5],
     71    segment_tags = {'wall':[1,2,3,4,5,6,7],
    6072                    'wave':[0]} # '':[6]
    6173       
    6274    m.add_points_and_segments(points, segments, segment_tags)
    6375   
    64     dam = m.add_region(xslope - 0.0000001,(ytop - ybottom)/2)
    65     # this is the location of the reservoir region.
    66     dam.setTag("flat")
     76#     dam = m.add_region(xslope - 0.0000001,(ytop - ybottom)/2)
     77#     # this is the location of the reservoir region.
     78#     dam.setTag("flat")
    6779   
    68     slope = m.add_region(xslope + 0.0000001,(ytop - ybottom)/2)
    69     # this is the location of the slope region.
    70     slope.setTag("slope")
     80#     slope = m.add_region(xslope + 0.0000001,(ytop - ybottom)/2)
     81#     # this is the location of the slope region.
     82#     slope.setTag("slope")
    7183   
    7284    m.generate_mesh(maximum_triangle_area=maximum_triangle_area)

  • anuga_work/development/Hinwood_2008/prepare_time_boundary.py

    r5350 r5370  
    2525def prepare_time_boundary(filename):
    2626    """Convert benchmark 2 time series to NetCDF tms file.
    27     This is a 'throw-away' code taylor made for files like
    28     'Benchmark_2_input.txt' from the LWRU2 benchmark
     27   
    2928    """
    3029
    … …  
    3534    # Read the ascii (.csv) version of this file
    3635    fid = open(filename[:-4] + '.csv')
    37 
    38     # Skip first line
    39     line = fid.readline()
    4036
    4137    # Read remaining lines
    … …  
    6258
    6359
     60   
    6461    # Create tms NetCDF file
    65 
    6662    fid = NetCDFFile(filename, 'w')
    6763    fid.institution = 'Geoscience Australia'
    … …  
    8783    """
    8884   
    89     Convert the rawish velocity and depth values, which have values at
     85    Convert the raw velocity and depth values, which have values at
    9086    different times to a csv file, with values at the same time, with
    9187    SI units.
     88
     89    Set the depth values to be at the same times as the velocity values.
    9290
    9391    The format for the velocity file is;
    … …  
    9896    [time, sec], [mm above SWL for sensor A], many other sensors...
    9997    """
    100 
     98    missing = 1e+20
     99   
    101100    # Read velocity file
    102101    vfid = open(velocity_file)
    … …  
    107106    n_velocity = len(lines)
    108107    vtimes = zeros(n_velocity, Float)  #Time
    109     velocities = zeros(n_velocity, Float)  #
     108    x_velocities = zeros(n_velocity, Float)  #
     109    y_velocities = zeros(n_velocity, Float)  #
    110110
    111111    for i, line in enumerate(lines):
    112         fields = line.split(',')
     112        fields = line.split() #(',')
    113113
    114114        vtimes[i] = float(fields[0])
    115         velocities[i] = float(fields[1])
     115        x_velocities[i] = float(fields[1])
     116        y_velocities[i] = float(fields[2])
    116117
    117118    # Read the depth file
    … …  
    120121    dfid.close()
    121122
    122 
     123   
    123124    n_depth = len(lines)
     125    n_sensors = len(lines[0].split())
    124126    dtimes = zeros(n_depth, Float)  #Time
    125127    depths = zeros(n_depth, Float)  #
     128    sensors = zeros((n_depth,n_sensors), Float)
    126129   
    127130    for i, line in enumerate(lines):
    128         fields = line.split(',')
    129 
    130         dtimes[i] = float(fields[0])
    131         depths[i] = float(fields[1])
    132 
    133     depths_at_vtimes = interp(dtimes, depths, vtimes, missing=1e+20)
     131        fields = line.split() #(',')
     132        fields = [float(j) for j in fields]
     133        dtimes[i] = fields[0]
     134        depths[i] = fields[1]
     135        sensors[i] = fields
     136    #print "dtimes", dtimes
     137    #print "depths", depths
     138    #print "vtimes", vtimes
     139    depths_at_vtimes = interp( depths, dtimes,
     140                               vtimes, missing=missing)
    134141    depths_at_vtimes = ensure_numeric(depths_at_vtimes)
     142
     143    #print "len(dtimes)", len(vtimes)
     144    #print "len(depths_at_vtimes)", len(depths_at_vtimes)
     145#     for i in range(len(depths_at_vtimes)):
     146#         print "i", i
     147#         print "vtimes[i]", vtimes[i]
     148#         print "depths_at_vtimes[i]", depths_at_vtimes[i]
     149#     print "depths_at_vtimes",  depths_at_vtimes
    135150    depths_at_vtimes = depths_at_vtimes/1000.00 # convert from mm to m
    136     velocities = ensure_numeric(velocities)
    137     velocities = velocities * -1.0 # Swap axis around
     151    missing=missing/1000.00 # Do to missing what is done to depths_at_vtimes
     152    x_velocities = ensure_numeric(x_velocities)
     153    x_velocities = x_velocities * -1.0 # Swap axis around
     154    y_velocities = ensure_numeric(y_velocities)
     155    y_velocities = y_velocities * -1.0 # Swap axis around
    138156   
    139157    fid = open(out_file,'w')
    140158
    141159    assert len(depths_at_vtimes) == len(vtimes)
    142    
     160    start_time = None
     161    #start_time = 0.0
    143162    #for vtime, depth_at_vtime, velocity in map(vtimes, depths_at_vtimes,
    144163    #                                           velocities):
    145164    for i in xrange(len(vtimes)):
    146         fid.write(str(vtimes[i]) + ',' + str(depths_at_vtimes[i]) \
    147                   + ',' + str(velocities[i])+'\n')
    148 
    149     fid.close()
     165        if not depths_at_vtimes[i] == missing:
     166            # Make the times start at zero.
     167            if start_time is None:
     168                start_time = vtimes[i]           
     169            fid.write(str(vtimes[i]-start_time) \
     170                      + ',' + str(depths_at_vtimes[i]) \
     171                      + ',' + str(x_velocities[i]) \
     172                      + ',' + str(y_velocities[i])+'\n')
     173
     174    fid.close()
     175   
     176    # Since there is a new time reference save the depth info using this
     177    # new reference.
     178    fid = open(depth_file[:-4] + '_new_time'+depth_file[-4:],'w')
     179    sensors[:,0] -= start_time
     180    #print "depth_file", depth_file
     181    #print "start_time", start_time
     182    for i in xrange(len(dtimes)):
     183        fid.write(str(sensors[i,0]))
     184        for j in xrange(1,len(sensors[0])):
     185            fid.write(' ' + str(sensors[i,j]))
     186        fid.write('\n')
     187    fid.close()
     188
    150189   
    151190   
    152191#-------------------------------------------------------------------
    153192if __name__ == "__main__":
    154     combine_velocity_depth('T2R7velfilt.csv','T2R7pressfilt.csv', 'cyeah.csv')
    155 
     193    from os import getenv
     194    from os.path import join
     195    home = getenv('INUNDATIONHOME')
     196    Hinwood_dir = join(home,'data','flumes','Hinwood_2008')
     197    raw_data_dir = join(Hinwood_dir, 'raw_data')
     198
     199    # Test 1 Run 5
     200    combine_velocity_depth(join(raw_data_dir,'T1R5velfilt.txt'),
     201                           join(raw_data_dir,'T1R5pressfilt.txt'),
     202                           join(Hinwood_dir, 'T1R5_boundary.csv'))   
     203    # Create the tsm file
     204    prepare_time_boundary(join(Hinwood_dir, 'T1R5_boundary.tsm'))
     205   
     206    # Test 2 Run 7
     207    combine_velocity_depth(join(raw_data_dir,'T2R7velfilt.txt'),
     208                           join(raw_data_dir,'T2R7pressfilt.txt'),
     209                           join(Hinwood_dir, 'T2R7_boundary.csv'))
     210    # Create the tsm file
     211    prepare_time_boundary(join(Hinwood_dir, 'T2R7_boundary.tsm'))

  • anuga_work/development/Hinwood_2008/project.py

    r5350 r5370  
    3838        self.meshdir = scenariodir+sep+'meshes'+sep
    3939        self.scenariodir = scenariodir+sep
    40         self.boundary_file = self.scenariodir + 'boundary.tsm'
     40        self.boundarydir = self.scenariodir
    4141
    4242        # creates copy of output dir structure, if it doesn't exist

  • anuga_work/development/Hinwood_2008/run_dam.py

    r5350 r5370  
    2222from shutil import copy
    2323from os import path, sep
    24 from os.path import dirname  #, basename
    25 
     24from os.path import dirname, join  #, basename
     25from Numeric import zeros, size, Float
    2626
    2727# Related major packages
    … …  
    4040import create_mesh
    4141from prepare_time_boundary import prepare_time_boundary
    42 
    43 
    44 
    45 def elevation_function(x,y):
    46     from Numeric import zeros, size, Float
    47 
    48     xslope = create_mesh.xslope  #4 ## Bit of a magic Number
    49     print "xslope",xslope
    50     z = zeros(size(x), Float)
    51     for i in range(len(x)):
    52         if x[i] < xslope:
    53             z[i] = 0.0  #WARNING: the code in prepare_time_boundary
    54                         # that calc's momentum assumes this is 0.0
    55         else:
    56             z[i] = (x[i]-xslope)*(1./16.)
    57     return z
    58 
    59 def main(friction=0.01, outputdir_name=None, is_trial_run=False):
     42from interp import interp
     43
     44
     45class Elevation_function:
     46    def __init__(self, slope):
     47        self.xslope_position = [slope['xleft'][0],slope['xtoe'][0],
     48                  slope['xbeach'][0],slope['xright'][0]]
     49        self.yslope_height = [slope['xleft'][1],slope['xtoe'][1],
     50                  slope['xbeach'][1],slope['xright'][1]]
     51       
     52    def __call__(self, x,y):
     53       
     54        z = interp(self.yslope_height, self.xslope_position, x)
     55        return z
     56
     57def main(boundary_file,
     58         boundary_location,
     59         slope,
     60         boundary_path=None,
     61         friction=0.01,
     62         outputdir_name=None,
     63         is_trial_run=False):
    6064
    6165   
    … …  
    8488          pro_instance.outputdir + 'create_mesh.py')
    8589
    86    
     90    # Boundary file manipulation
     91    if boundary_path is None:
     92        boundary_path = pro_instance.boundarydir
     93    boundary_file_path = join(boundary_path, boundary_file)
    8794    # Convert the boundary file, .csv to .tsm
    8895    try:
    89         temp = open(pro_instance.boundary_file)
     96        temp = open(boundary_file_path)
    9097        temp.close()
    9198    except IOError:
    92         prepare_time_boundary(pro_instance.boundary_file)
     99        prepare_time_boundary(boundary_file_path)
    93100   
    94101    mesh_filename = pro_instance.meshdir + basename + '.msh'
    … …  
    110117    # this creates the mesh
    111118    #gate_position = 12.0
    112     create_mesh.generate(mesh_filename,
     119    create_mesh.generate(mesh_filename, slope,
    113120                         maximum_triangle_area=maximum_triangle_area)
    114121
    … …  
    138145
    139146    domain.set_quantity('stage', 0.4)
    140     domain.set_quantity('friction', friction)
     147    domain.set_quantity('friction', friction)
     148    elevation_function = Elevation_function(slope)
    141149    domain.set_quantity('elevation', elevation_function)
    142150
    … …  
    148156    #Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
    149157    try:
    150         function = file_function(pro_instance.boundary_file, domain,
     158        function = file_function(boundary_file_path, domain,
    151159                                 verbose=True)
    152160    except IOError:
    … …  
    156164       
    157165    Br = Reflective_boundary(domain)
    158     #Bd = Dirichlet_boundary([0.3,0,0])
     166    Bd = Dirichlet_boundary([0.3,0,0])
    159167    Bts = Time_boundary(domain, function)
    160168    domain.set_boundary( {'wall': Br, 'wave': Bts} )
     169    domain.set_boundary( {'wall': Br, 'wave': Bd} )
    161170
    162171    #-------------------------------------------------------------------------
    … …  
    165174    t0 = time.time()
    166175
     176    # It seems that ANUGA can't handle a starttime that is >0.
     177    domain.starttime = 1.0
    167178    for t in domain.evolve(yieldstep, finaltime):   
    168179        domain.write_time()
    … …  
    196207#-------------------------------------------------------------
    197208if __name__ == "__main__":
    198     main( is_trial_run = True,
     209    #slopes = [[-4.5,0.0],[0.0,0.0],[1.285,0.090],[16.1,.960]]
     210    slope = {'xleft':[-4.5,0.0],
     211              'xtoe':[0.0,0.0],
     212              'xbeach':[1.285,0.090],
     213              'xright':[16.1,.960]}
     214    main( 'T1R5_boundary.tsm', 8, slope,
     215          is_trial_run = True,
    199216         outputdir_name='Hinwood_low_stage_low_velocity_draft')

  • anuga_work/development/friction_UA_flume_2006/parallel_run_dam.py

    r4443 r5370  
    99
    1010e.g.
    11 miprun -c 4 python parallel_run_dam.py
    12 
    13 or
    14 e.g.
    15 miprun c4-8 python parallel_run_dam.py
     11mpirun -c 4 python parallel_run_dam.py
     12
    1613
    1714
Note: See TracChangeset for help on using the changeset viewer.