Changeset 5084


Ignore:
Timestamp:
Feb 27, 2008, 10:02:55 AM (16 years ago)
Author:
ole
Message:

Simplified debug of transmissive_momentum_set_stage

Location:
anuga_work/debug/bad_time_step
Files:
3 edited
2 moved

Legend:

Unmodified
Added
Removed
  • anuga_work/debug/bad_time_step/create_mesh.py

    r5073 r5084  
    1414   
    1515    """
     16   
    1617    #Basic geometry
    17    
    1818    xright  = 19.0
    1919    ybottom = 0
  • anuga_work/debug/bad_time_step/project.py

    r5075 r5084  
    1212from anuga.utilities.system_tools import get_user_name
    1313
    14 boundary_file = 'boundary_mellow.tsm'
     14from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
     15     import File_boundary_time
     16boundary_file = 'boundary.tms'
     17#boundary_file = 'boundary_mellow.tms'
    1518
    1619class Project:   
  • anuga_work/debug/bad_time_step/run_dam.py

    r5079 r5084  
    66
    77Duncan Gray, GA - 2007
    8 
    9 
    10 
    118"""
    129
     
    1613#----------------------------------------------------------------------------
    1714
    18 # Standard modules
    19 import 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  #, basename
     15# Related major packages
     16from anuga.shallow_water import Domain
     17from anuga.shallow_water import Reflective_boundary
     18from anuga.shallow_water import Dirichlet_boundary
     19from anuga.shallow_water import Time_boundary
     20from anuga.shallow_water import File_boundary
     21from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
    2522
    2623
    27 # 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
    32 from anuga.fit_interpolate.interpolate import interpolate_sww2csv
    33 from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, \
    34      copy_code_files, file_function
    35 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
    36      import File_boundary_time
     24from anuga.abstract_2d_finite_volumes.util import file_function
     25#from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
     26from anuga.pmesh.mesh_interface import create_mesh_from_regions
     27
     28
    3729
    3830# Scenario specific imports
    39 import project                 # Definition of file names and polygons
    4031import create_mesh
     32
    4133
    4234def elevation_function(x,y):
     
    5244    return z
    5345
    54 def main(friction=0.01, outputdir_name=None, is_trial_run=False):
    55     basename = 'zz' + str(friction)
    56     if is_trial_run is True:
    57         outputdir_name += '_test'
    58         yieldstep = 0.1
    59         finaltime = 15.
    60         maximum_triangle_area=0.01
    61         #maximum_triangle_area=0.1 #OK
    62     else:
    63         yieldstep = 0.02
    64         finaltime = 15.1
    65         maximum_triangle_area=0.0001
    66        
    67         maximum_triangle_area=0.001
    68        
    69     pro_instance = project.Project([],
    70                                    outputdir_name=outputdir_name,
    71                                    home='.')
    72     print 'The output dir is', pro_instance.outputdir
    73     copy_code_files(pro_instance.outputdir,__file__,
    74                     dirname(project.__file__) \
    75                     + sep + project.__name__+'.py')
    76     copy (pro_instance.codedir + 'run_dam.py',
    77           pro_instance.outputdir + 'run_dam.py')
    78     copy (pro_instance.codedir + 'create_mesh.py',
    79           pro_instance.outputdir + 'create_mesh.py')
     46
     47#-------------------------------------------------------------------------
     48# Model parameters
     49#-------------------------------------------------------------------------
     50friction=0.01
     51maximum_triangle_area=0.01
     52#maximum_triangle_area=0.1 #OK
     53
     54basename = 'debug_boundary'
     55mesh_filename = basename + '.msh'
     56
     57boundary_file = 'boundary.tms'
     58#boundary_file = 'boundary_mellow.tms'
     59
     60
     61
     62#-------------------------------------------------------------------------
     63# Create the triangular mesh
     64#-------------------------------------------------------------------------
     65
     66# Original unstructured mesh ()
     67#create_mesh.generate(mesh_filename,
     68#                     maximum_triangle_area=maximum_triangle_area)
     69
     70
     71# Simple mesh with same dimensions
     72
     73xright  = 19.0
     74ybottom = 0
     75ytop    = 0.45
     76xleft = 0.0
     77#xslope = 4.0
     78
     79#Outline
     80point_sw = [xleft, ybottom]
     81point_se = [xright, ybottom]
     82point_nw = [xleft, ytop]   
     83point_ne = [xright, ytop]
     84
     85
     86create_mesh_from_regions([point_sw,
     87                          point_se,
     88                          point_ne,
     89                          point_nw],
     90                         {'wall': [0,1,2], 'wave': [3]},
     91                         maximum_triangle_area=maximum_triangle_area,
     92                         filename=mesh_filename)
     93
     94
     95#TODO
     96#points, elements, boundary = rectangular_cross(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
     97
     98
     99
     100#-------------------------------------------------------------------------
     101# Setup computational domain
     102#-------------------------------------------------------------------------
     103domain = Domain(mesh_filename, use_cache = False, verbose = True)
     104print domain.statistics()
     105
     106domain.set_name(basename)
     107domain.set_datadir('.')
     108
     109
     110#-------------------------------------------------------------------------
     111# Setup initial conditions
     112#-------------------------------------------------------------------------
     113
     114domain.set_quantity('stage', 0.06)
     115domain.set_quantity('friction', friction)
     116domain.set_quantity('elevation', elevation_function)
     117
     118
     119print 'Available boundary tags', domain.get_boundary_tags()
     120
     121# Create boundary function from timeseries provided in file
     122function = file_function(boundary_file, domain, verbose=True)
     123
     124
     125Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
     126Br = Reflective_boundary(domain)
     127domain.set_boundary( {'wall': Br, 'wave': Bts} )
     128
     129
     130#-------------------------------------------------------------------------
     131# Evolve system through time
     132#-------------------------------------------------------------------------
     133
     134for t in domain.evolve(0.1, 15):
     135    print domain.timestepping_statistics(track_speeds=False)
     136    print domain.boundary_statistics(tags=['wave'])
    80137   
    81     mesh_filename = pro_instance.meshdir + basename + '.msh'
     138print 'finished'
    82139
    83     #--------------------------------------------------------------------------
    84     # Copy scripts to output directory and capture screen
    85     # output to file
    86     #--------------------------------------------------------------------------
    87 
    88     # creates copy of code in output dir
    89     if is_trial_run is False:
    90         start_screen_catcher(pro_instance.outputdir, rank, pypar.size())
    91 
    92     print 'USER:    ', pro_instance.user
    93 
    94    
    95     #-------------------------------------------------------------------------
    96     # Create the triangular mesh
    97     #-------------------------------------------------------------------------
    98 
    99     # this creates the mesh
    100     #gate_position = 12.0
    101     create_mesh.generate(mesh_filename,
    102                          maximum_triangle_area=maximum_triangle_area)
    103 
    104     head,tail = path.split(mesh_filename)
    105     copy (mesh_filename,
    106           pro_instance.outputdir + tail )
    107     #-------------------------------------------------------------------------
    108     # Setup computational domain
    109     #-------------------------------------------------------------------------
    110     domain = Domain(mesh_filename, use_cache = False, verbose = True)
    111    
    112 
    113     print 'Number of triangles = ', len(domain)
    114     print 'The extent is ', domain.get_extent()
    115     print domain.statistics()
    116 
    117    
    118     domain.set_name(basename)
    119     domain.set_datadir(pro_instance.outputdir)
    120     domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
    121     domain.set_minimum_storable_height(0.001)
    122     #domain.set_store_vertices_uniquely(True)  # for writting to sww
    123 
    124     #-------------------------------------------------------------------------
    125     # Setup initial conditions
    126     #-------------------------------------------------------------------------
    127 
    128     domain.set_quantity('stage', 0.06)
    129     domain.set_quantity('friction', friction)
    130     domain.set_quantity('elevation', elevation_function)
    131 
    132    
    133     print 'Available boundary tags', domain.get_boundary_tags()
    134 
    135     # Create boundary function from timeseries provided in file
    136     function = file_function(project.boundary_file, domain, verbose=True)
    137 
    138     #print dir(function)
    139     #print function.time, function.quantities_range
    140     #import sys; sys.exit()
    141    
    142     Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
    143 
    144     Br = Reflective_boundary(domain)
    145     #Bd = Dirichlet_boundary([0.3,0,0])
    146     #Bts = Time_boundary(domain, function)
    147     domain.set_boundary( {'wall': Br, 'wave': Bts} )
    148 
    149     #-------------------------------------------------------------------------
    150     # Evolve system through time
    151     #-------------------------------------------------------------------------
    152     t0 = time.time()
    153 
    154     for t in domain.evolve(yieldstep, finaltime):
    155         print domain.timestepping_statistics(track_speeds=False)
    156         print domain.boundary_statistics(tags=['wave'])
    157        
    158     print 'That took %.2f seconds' %(time.time()-t0)
    159     print 'finished'
    160 
    161     points = [[2.8,0.225],  #-1.8m from SWL
    162               [5.1,0.225],  #0.5m from SWL
    163               [6.6,0.225],  #2m from SWL
    164               [6.95,0.255], #2.35m from SWL
    165               [7.6,0.255],  #3m from SWL
    166               [8.2,0.255],  #3.5m from SWL
    167               [9.2,0.255]  #4.5m from SWL
    168               ]
    169 
    170 #     points = [[gate_position - 0.65,0.2],
    171 #               [gate_position - 0.55,0.2],
    172 #               [gate_position - 0.45,0.2],
    173 #               [gate_position - 0.35,0.2],
    174 #               [gate_position - 0.25,0.2]
    175 #               ]
    176 
    177     #-------------------------------------------------------------------------
    178     # Calculate gauge info
    179     #-------------------------------------------------------------------------
    180 
    181     if False:
    182         interpolate_sww2csv(pro_instance.outputdir + basename +".sww",
    183                             points,
    184                             pro_instance.outputdir + "depth_manning_"+str(friction)+".csv",
    185                             pro_instance.outputdir + "velocity_x.csv",
    186                             pro_instance.outputdir + "velocity_y.csv")
    187  
    188     return pro_instance
    189 
    190 #-------------------------------------------------------------
    191 if __name__ == "__main__":
    192     main( is_trial_run = True,
    193          outputdir_name='Hinwood_draft')
Note: See TracChangeset for help on using the changeset viewer.