Ignore:
Timestamp:
Aug 26, 2008, 3:23:31 PM (16 years ago)
Author:
duncan
Message:

Simplifying files

File:
1 edited

Legend:

Unmodified
Added
Removed

  • anuga_validation/Hinwood_2008/run_dam.py

    r5682 r5686  
    22
    33Script for running a breaking wave simulation of Jon Hinwoods wave tank.
    4 Note: this is based on the frinction_ua_flume_2006 structure.
    5 
    64
    75Duncan Gray, GA - 2007
    8 
    9 
    10 
    116"""
    127
    … …  
    1813# Standard modules
    1914import time
    20 from time import localtime, strftime
    21 import sys
    22 from shutil import copy
    23 from os import path, sep
    24 from os.path import dirname, join  #, basename
    25 from Numeric import zeros, size, Float
     15
    2616
    2717# Related major packages
    28 from anuga.shallow_water import Domain, Reflective_boundary, \
    29                             Dirichlet_boundary,  Time_boundary, \
    30                             File_boundary, \
    31                             Transmissive_Momentum_Set_Stage_boundary
     18from anuga.shallow_water import Domain, Reflective_boundary, Time_boundary
    3219from anuga.fit_interpolate.interpolate import interpolate_sww2csv
    33 from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, \
    34       file_function
    35 from anuga.shallow_water.data_manager import copy_code_files
    36 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
    37      import File_boundary_time
     20from anuga.abstract_2d_finite_volumes.util import file_function
     21from anuga.utilities.interp import interp
    3822
    3923# Scenario specific imports
    40 import project                 # Definition of file names and polygons
    4124import create_mesh
    42 from prepare_time_boundary import prepare_time_boundary
    43 from anuga.utilities.interp import interp
    44 
    45 
    46 class Elevation_function:
    47     def __init__(self, slope):
    48         self.xslope_position = [slope['xleft'][0],slope['xtoe'][0],
    49                   slope['xbeach'][0],slope['xright'][0]]
    50         self.yslope_height = [slope['xleft'][1],slope['xtoe'][1],
    51                   slope['xbeach'][1],slope['xright'][1]]
    52        
    53     def __call__(self, x,y):
    54        
    55         z = interp(self.yslope_height, self.xslope_position, x)
    56         return z
    57 
    58 def main(boundary_file,
    59          metadata_dic,
     25from prepare_time_boundary import csv2tms
     26
     27
     28def main(metadata_dic,
    6029         maximum_triangle_area,
    6130         yieldstep,
    6231         boundary_path=None,
    63          friction=0.012,  # planed wood. http://www.lmnoeng.com/manningn.htm
     32         friction=0.,  # planed wood 0.012 http://www.lmnoeng.com/manningn.htm
    6433         outputdir_name=None,
    65          run_type=1,
     34         isTest=False,
    6635         width=1.0,
    6736         use_limits=True,
    68          end_tag = '_limiterD'):
    69 
    70    
    71     basename = 'zz_' + metadata_dic['scenario_id']
    72    
    73     if run_type == 0:
     37         end_tag = ''):
     38
     39    id = metadata_dic['scenario_id']
     40   
     41    if isTest is True:
    7442        yieldstep = 1.0
    7543        finaltime = 15.
    7644        maximum_triangle_area=0.1
    77         outputdir_name += '_test'
    78          
    79     elif run_type == 1:
    80         finaltime = None
    81        
    82     elif run_type == 2:
    83         yieldstep = 0.5
     45        outputdir_name += '_test'         
     46    else:
    8447        finaltime = None
    85         maximum_triangle_area=0.01
    86         outputdir_name += '_test_long_time'
    87        
    88     elif run_type == 3:
    89         yieldstep = 0.1
    90         finaltime = None       
    91         maximum_triangle_area=0.01
    92         #outputdir_name += '_yieldstep_0.1'
    93        
    94     elif run_type == 4:
    95         # this is not a test
    96         # Output will go to a file
    97         # The sww file will be interpolated
    98         yieldstep = 0.01
    99         finaltime = None       
    100         maximum_triangle_area=0.01
    101        
    102     elif run_type == 5:
    103         # this is not a test
    104         # Output will go to a file
    105         # The sww file will be interpolated
    106         yieldstep = 0.01
    107         finaltime = None       
    108         maximum_triangle_area=0.001
    109         #outputdir_name += '_good'
    110        
    111     elif run_type == 6:
    112         # this is not a test
    113         # Output will go to a file
    114         # The sww file will be interpolated
    115         yieldstep = 0.01
    116         finaltime = None       
    117         maximum_triangle_area=0.0001
    118         #outputdir_name += '_good'
    119        
    120     elif run_type == 7:
    121         # this is not a test
    122         # Output will go to a file
    123         # The sww file will be interpolated
    124         yieldstep = 0.01
    125         finaltime = None       
    126         maximum_triangle_area=0.00001
    127         #outputdir_name += '_good'
    128        
    129     elif run_type == 8:
    130         # this is not a test
    131         # Output will go to a file
    132         # The sww file will be interpolated
    133         yieldstep = 0.05
    134         finaltime = None       
    135         maximum_triangle_area=0.00001
    136         #outputdir_name += '_good'
    137     #outputdir_name += '_no_velocity'
    138 
     48       
    13949    outputdir_name = paras2outputdir_tag(mta=maximum_triangle_area,
    14050                        yieldstep=yieldstep,
    … …  
    14353                        friction=friction,
    14454                        end_tag=end_tag,
    145                         outputdir_name=outputdir_name
    146                         )
    147 
     55                        outputdir_name=outputdir_name)
     56    basename = outputdir_name   
    14857    metadata_dic = set_z_origin_to_water_depth(metadata_dic)   
    149        
    150     pro_instance = project.Project(['data','flumes','Hinwood_2008'],
    151                                    outputdir_name=outputdir_name)
    152     print "The output dir is", pro_instance.outputdir
    153     copy_code_files(pro_instance.outputdir,__file__,
    154                     dirname(project.__file__) \
    155                     + sep + project.__name__+'.py')
    156     copy (pro_instance.codedir + 'run_dam.py',
    157           pro_instance.outputdir + 'run_dam.py')
    158     copy (pro_instance.codedir + 'create_mesh.py',
    159           pro_instance.outputdir + 'create_mesh.py')
    160 
    161     boundary_final_time = prepare_time_boundary(metadata_dic,
    162                                        pro_instance.raw_data_dir,
    163                                        pro_instance.boundarydir)
    164     #return pro_instance
     58
    16559    if finaltime is None:
    166         finaltime = boundary_final_time - 0.1 # Edge boundary problems
    167     # Boundary file manipulation
    168     if boundary_path is None:
    169         boundary_path = pro_instance.boundarydir
    170     boundary_file_path = join(boundary_path, boundary_file)
    171    #  # Convert the boundary file, .csv to .tsm
    172 #     try:
    173 #         temp = open(boundary_file_path)
    174 #         temp.close()
    175 #     except IOError:
    176 #         prepare_time_boundary(boundary_file_path)
    177    
    178     mesh_filename = pro_instance.meshdir + basename + '.msh'
    179 
    180     #--------------------------------------------------------------------------
    181     # Copy scripts to output directory and capture screen
    182     # output to file
    183     #--------------------------------------------------------------------------
    184 
    185     # creates copy of code in output dir
    186     if run_type >= 1:
    187         #start_screen_catcher(pro_instance.outputdir, rank, pypar.size())
    188         start_screen_catcher(pro_instance.outputdir)
    189 
    190     print 'USER:    ', pro_instance.user
     60        finaltime = metadata_dic['wave_times'][1] + 0.1
     61       
     62    boundary_file_path = id + '_boundary.tsm'
     63       
     64    # Convert the boundary file, .csv to .tsm
     65    csv2tms(boundary_file_path, metadata_dic['xleft'][1])
     66   
     67    mesh_filename =  basename + '.msh'
     68
    19169    #-------------------------------------------------------------------------
    19270    # Create the triangular mesh
    19371    #-------------------------------------------------------------------------
    19472
    195     # this creates the mesh
    196     #gate_position = 12.0
    19773    create_mesh.generate(mesh_filename, metadata_dic, width=width,
    19874                         maximum_triangle_area=maximum_triangle_area)
    199 
    200     head,tail = path.split(mesh_filename)
    201     copy (mesh_filename,
    202           pro_instance.outputdir + tail )
     75   
    20376    #-------------------------------------------------------------------------
    20477    # Setup computational domain
    20578    #-------------------------------------------------------------------------
    206     domain = Domain(mesh_filename, use_cache = False, verbose = True)
    207    
     79    domain = Domain(mesh_filename, use_cache = False, verbose = True) 
    20880
    20981    print 'Number of triangles = ', len(domain)
    21082    print 'The extent is ', domain.get_extent()
    21183    print domain.statistics()
    212 
    21384   
    21485    domain.set_name(basename)
    215     domain.set_datadir(pro_instance.outputdir)
     86    domain.set_datadir('.')
    21687    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
    21788    domain.set_minimum_storable_height(0.0001)
    … …  
    237108    domain.set_quantity('elevation', elevation_function)
    238109
    239    
    240     print 'Available boundary tags', domain.get_boundary_tags()
    241 
    242     # Create boundary function from timeseries provided in file
    243     #function = file_function(project.boundary_file, domain, verbose=True)
    244     #Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
    245     try:
    246         function = file_function(boundary_file_path, domain,
    247                                  verbose=True)
    248         print "boundary_file_path", boundary_file_path
    249     except IOError:
    250         msg = 'Run prepare_time_boundary.py. File "%s" could not be opened.'\
    251                   %(pro_instance.boundary_file)
    252         raise msg
     110    function = file_function(boundary_file_path, domain, verbose=True)
    253111       
    254112    Br = Reflective_boundary(domain)
    255     Bd = Dirichlet_boundary([0.3,0,0])
    256113    Bts = Time_boundary(domain, function)
    257     #Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
    258114    domain.set_boundary( {'wall': Br, 'wave': Bts} )
    259     #domain.set_boundary( {'wall': Br, 'wave': Bd} )
    260115
    261116    #-------------------------------------------------------------------------
    … …  
    264119    t0 = time.time()
    265120
    266     # It seems that ANUGA can't handle a starttime that is >0.
    267     #domain.starttime = 1.0 #!!! what was this doing?
    268121    for t in domain.evolve(yieldstep, finaltime):   
    269122        domain.write_time()
    … …  
    271124    print 'finished'
    272125
    273     flume_y_middle = 0.0
    274     points = []
    275     for gauge_x in metadata_dic['gauge_x']:
    276         points.append([gauge_x, flume_y_middle])
    277 
    278126
    279127    #-------------------------------------------------------------------------
    … …  
    281129    #-------------------------------------------------------------------------
    282130
    283     if run_type >= 1:
     131    if isTest is not True:
     132        points = []
     133        for gauge_x in metadata_dic['gauge_x']:
     134            points.append([gauge_x, 0.0])
     135
    284136        id = metadata_dic['scenario_id'] + ".csv"
    285         interpolate_sww2csv(pro_instance.outputdir + basename +".sww",
     137        interpolate_sww2csv(basename +".sww",
    286138                            points,
    287                             pro_instance.outputdir + "depth_" + id,
    288                             pro_instance.outputdir + "velocity_x_" + id,
    289                             pro_instance.outputdir + "velocity_y_" + id,
    290                             pro_instance.outputdir + "stage_" + id,
    291                             pro_instance.outputdir + "froude_" + id)
     139                            basename + "_depth_" + id,
     140                            basename + "_velocity_x_" + id,
     141                            basename + "_velocity_y_" + id,
     142                            basename + "_stage_" + id)
    292143 
    293     return pro_instance
    294 
    295 
    296144
    297145def paras2outputdir_tag(mta,
    … …  
    301149                        friction,
    302150                        end_tag,
    303                         outputdir_name=None
    304                         ):
    305 
     151                        outputdir_name=None):
     152    """
     153   
     154    """
    306155    if outputdir_name is None:
    307156        outputdir_name = ''
    … …  
    319168    return outputdir_name
    320169
    321    
    322 def set_z_origin_to_water_depth(seabed_coords):
    323     offset = seabed_coords['offshore_water_depth']
     170
     171class Elevation_function:
     172    """
     173    """
     174    def __init__(self, slope):
     175        self.xslope_position = [slope['xleft'][0],slope['xtoe'][0],
     176                  slope['xbeach'][0],slope['xright'][0]]
     177        self.yslope_height = [slope['xleft'][1],slope['xtoe'][1],
     178                  slope['xbeach'][1],slope['xright'][1]]
     179       
     180    def __call__(self, x,y):       
     181        z = interp(self.yslope_height, self.xslope_position, x)
     182        return z
     183   
     184def set_z_origin_to_water_depth(metadata_dic):
     185    """
     186    """
     187    offset = metadata_dic['offshore_water_depth']
    324188    keys = ['xleft', 'xtoe', 'xbeach', 'xright']
    325189    for x in keys:
    326             seabed_coords[x][1] -= offset
    327     return seabed_coords
     190            metadata_dic[x][1] -= offset
     191    return metadata_dic
     192
    328193#-------------------------------------------------------------
    329194if __name__ == "__main__":
    330195   
    331196    from scenarios import scenarios
    332     from slope import gauges_for_slope
    333     #from plot import plot
    334 
    335     # 1 is fast and dirty
    336     # 4 is 0.01
    337     # 5 is 0.001
    338     # 6 is 0.0001
    339     # 7 is 0.00001
    340     # 8 is 0.00001  yieldstep = 0.05
    341    
    342     run_type = 1
    343     run_type = 0  # for testing
    344    
    345     #for run_data in [scenarios[5]]:
    346     #scenarios = scenarios[3:]
     197
    347198    scenarios = [scenarios[0]]
    348199   
    349     width = 1.0
    350200    width = 0.1
    351     #width = 0.01
    352 
    353     maximum_triangle_area=0.0001
     201    maximum_triangle_area=0.01
    354202    yieldstep = 0.01
     203    yieldstep = 0.5
    355204    friction=0.0
     205    isTest=True
     206    #isTest=False
    356207   
    357208    for run_data in scenarios:
    358         pro_instance = main( run_data['scenario_id'] + '_boundary.tsm'  ,
    359                              run_data,
    360                              maximum_triangle_area=maximum_triangle_area,
    361                              yieldstep=yieldstep,
    362                              width=width,
    363                              run_type=run_type,
    364                              outputdir_name=run_data['scenario_id'],
    365                              use_limits=False,
    366                              friction=friction,
    367                              end_tag='_I')
    368         #gauges_for_slope(pro_instance.outputdir,[run_data])
     209        main(run_data,
     210             maximum_triangle_area=maximum_triangle_area,
     211             yieldstep=yieldstep,
     212             width=width,
     213             isTest=isTest,
     214             outputdir_name=run_data['scenario_id'],
     215             use_limits=False,
     216             friction=friction,
     217             end_tag='')
Note: See TracChangeset for help on using the changeset viewer.