Ignore:
Timestamp:
Feb 14, 2006, 5:19:29 PM (19 years ago)
Author:
ole
Message:

Pushed caching into conversion functions and beautified the Sydney scenario

File:
1 edited

Legend:

Unmodified
Added
Removed

  • production/sydney_2006/run_sydney_smf.py

    r2400 r2407  
    66
    77The scenario is defined by a triangular mesh created from project.polygon,
    8 the elevation data and boundary data obtained from a tsunami simulation done with MOST.
     8the elevation data and a simulated submarine landslide.
    99
    10 Ole Nielsen, GA - 2005 and Adrian Hitchman and Jane Sexton, GA - 2006
     10Ole Nielsen and Duncan Gray, GA - 2005 and Adrian Hitchman and Jane Sexton, GA - 2006
    1111"""
    1212
    13 tide = 0       #Australian Height Datum (mean sea level)
    1413
     14#-------------------------------------------------------------------------------# Import necessary modules
     15#-------------------------------------------------------------------------------
     16
     17# Standard modules
    1518import os
    1619import time
    1720
     21# Related major packages
     22from pyvolution.shallow_water import Domain, Reflective_boundary
     23from pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
     24from pyvolution.combine_pts import combine_rectangular_points_files
     25from pyvolution.pmesh2domain import pmesh_to_domain_instance
    1826
    19 from pyvolution.shallow_water import Domain, Reflective_boundary, File_boundary,\
    20      Dirichlet_boundary, Time_boundary, Transmissive_boundary
    21 from pyvolution.data_manager import convert_dem_from_ascii2netcdf,\
    22      dem2pts
    23 from pyvolution.pmesh2domain import pmesh_to_domain_instance
    24 from pyvolution.combine_pts import combine_rectangular_points_files
    25 from caching import cache
    26 import project
    27 from math import pi, sin
    28 from smf import slump_tsunami, slide_tsunami, Double_gaussian
    29 from pyvolution.least_squares import fit_to_mesh_file, DEFAULT_ALPHA
     27# Application specific imports
     28import project                 # Definition of file names and polygons
     29from smf import slump_tsunami  # Function for submarine mudslide
    3030
    31 # Data preparation
     31
     32#-------------------------------------------------------------------------------
     33# Preparation of topographic data
     34#
    3235# Convert ASC 2 DEM 2 PTS using source data and store result in source data
    3336# Do for coarse and fine data
    3437# Fine pts file to be clipped to area of interest
     38#-------------------------------------------------------------------------------
     39
     40# filenames
    3541coarsedemname = project.coarsedemname
    3642finedemname = project.finedemname
    … …  
    3844
    3945# coarse data
    40 cache(convert_dem_from_ascii2netcdf, coarsedemname, {'verbose': True},
    41       dependencies = [coarsedemname + '.asc'],
    42       verbose = True)
    43       #evaluate = True)
     46convert_dem_from_ascii2netcdf(coarsedemname, use_cache=True, verbose=True)
     47dem2pts(coarsedemname, use_cache=True, verbose=True)
    4448
    45 cache(dem2pts, coarsedemname, {'verbose': True},
    46       dependencies = [coarsedemname + '.dem'],     
    47       verbose = True)
     49# fine data (clipping the points file to smaller area)
     50convert_dem_from_ascii2netcdf(finedemname, use_cache=True, verbose=True)
     51dem2pts(finedemname,
     52        easting_min=project.eastingmin,
     53        easting_max=project.eastingmax,
     54        northing_min=project.northingmin,
     55        northing_max= project.northingmax,
     56        use_cache=True,
     57        verbose=True)
    4858
    49 # fine data
    50 cache(convert_dem_from_ascii2netcdf, finedemname, {'verbose': True},
    51       dependencies = [finedemname + '.asc'],
    52       verbose = True)
    53       #evaluate = True)
    54 
    55 # clipping the fine data
    56 cache(dem2pts, finedemname, {'verbose': True,
    57       'easting_min': project.eastingmin,
    58       'easting_max': project.eastingmax,
    59       'northing_min': project.northingmin,
    60       'northing_max': project.northingmax},
    61       dependencies = [finedemname + '.dem'],
    62       #evaluate = True,
    63       verbose = True)
    6459
    6560# combining the coarse and fine data
    … …  
    6762                                 project.coarsedemname + '.pts',
    6863                                 project.combineddemname + '.pts')
    69                                  
    70 # Create Triangular Mesh
    71 # Overall clipping polygon and interior regions defined in project.py
     64
     65
     66#-------------------------------------------------------------------------------                                 
     67# Create the triangular mesh based on overall clipping polygon with a tagged
     68# boundary and interior regions defined in project.py along with
     69# resolutions (maximal area of per triangle) for each polygon
     70#-------------------------------------------------------------------------------                                 
     71
    7272from pmesh.create_mesh import create_mesh_from_regions
    7373
    74 # for whole region
    75 interior_res = 5000 # maximal area of per triangle
     74interior_res = 5000
    7675interior_regions = [[project.harbour_polygon_2, interior_res],
    77                     [project.botanybay_polygon_2, interior_res]] 
     76                    [project.botanybay_polygon_2, interior_res]]
    7877
    79 m = cache(create_mesh_from_regions,
    80             project.diffpolygonall,
    81             {'boundary_tags': {'bottom': [0],
    82                             'right1': [1], 'right0': [2],
    83                             'right2': [3], 'top': [4], 'left1': [5],
    84                             'left2': [6], 'left3': [7]},
    85             'resolution': 100000,
    86             'filename': meshname,           
    87             'interior_regions': interior_regions},
    88             #evaluate=True,   
    89             verbose = True)
    9078
    91 # Setup domain
    92 domain = cache(pmesh_to_domain_instance, (meshname, Domain),
    93                dependencies = [meshname],                     
    94                verbose = True)
     79#FIXME: Fix caching of this one once the mesh_interface is ready
     80from caching import cache
     81_ = cache(create_mesh_from_regions,
     82          project.diffpolygonall,
     83          {'boundary_tags': {'bottom': [0],
     84                             'right1': [1], 'right0': [2],
     85                             'right2': [3], 'top': [4], 'left1': [5],
     86                             'left2': [6], 'left3': [7]},
     87           'resolution': 100000,
     88           'filename': meshname,           
     89           'interior_regions': interior_regions,
     90           'UTM': True,
     91           'refzone': project.refzone},
     92          verbose = True)
    9593
    96 # Bring in elevation data
    97 domain.set_quantity('elevation',
    98                      filename = project.combineddemname + '.pts',
    99                      use_cache = True,
    100                      verbose = True)
    101  
     94
     95#-------------------------------------------------------------------------------                                 
     96# Setup computational domain
     97#-------------------------------------------------------------------------------                                 
     98
     99domain = pmesh_to_domain_instance(meshname, Domain,
     100                                  use_cache = True,
     101                                  verbose = True)
     102
    102103print 'Number of triangles = ', len(domain)
    103104print 'The extent is ', domain.get_extent()
    … …  
    105106domain.set_name(project.basename)
    106107domain.set_datadir(project.outputdir)
    107 domain.store = True
    108 domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
     108domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
    109109
    110 # Setup Initial conditions
    111 t = slump_tsunami(length=30000.0,
    112                   depth=400.0,
    113                   slope=6.0,
    114                   thickness=176.0,
    115                   radius=3330,
    116                   dphi=0.23,
    117                   x0=project.slump_origin[0],
    118                   y0=project.slump_origin[1],
    119                   alpha=0.0,
    120                   domain=domain)
    121 domain.set_quantity('stage', t)
     110
     111#-------------------------------------------------------------------------------
     112# Set up scenario (tsunami_source is a callable object used with set_quantity)
     113#-------------------------------------------------------------------------------
     114
     115tsunami_source = slump_tsunami(length=30000.0,
     116                               depth=400.0,
     117                               slope=6.0,
     118                               thickness=176.0,
     119                               radius=3330,
     120                               dphi=0.23,
     121                               x0=project.slump_origin[0],
     122                               y0=project.slump_origin[1],
     123                               alpha=0.0,
     124                               domain=domain)
     125
     126
     127#-------------------------------------------------------------------------------                                 
     128# Setup initial conditions
     129#-------------------------------------------------------------------------------
     130
     131domain.set_quantity('stage', tsunami_source)
    122132domain.set_quantity('friction', 0)
    123    
    124 # Setup Boundary Conditions
    125 print domain.get_boundary_tags()
     133domain.set_quantity('elevation',
     134                    filename = project.combineddemname + '.pts',
     135                    use_cache = True,
     136                    verbose = True)
     137
     138
     139#-------------------------------------------------------------------------------                                 
     140# Setup boundary conditions (all reflective)
     141#-------------------------------------------------------------------------------
     142
     143print 'Available boundary tags', domain.get_boundary_tags()
    126144
    127145Br = Reflective_boundary(domain)
    128 Bt = Transmissive_boundary(domain)
    129 Bd = Dirichlet_boundary([0,0,0])
    130 # 10 min square wave starting at 1 min, 6m high
    131 Bw = Time_boundary(domain=domain,
    132                    f=lambda t: [(6<t<606)*6, 0, 0])
    133 
    134146domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
    135147                      'right2': Br, 'top': Br, 'left1': Br,
    136148                      'left2': Br, 'left3': Br} )
    137149
    138 # Evolve
     150
     151#-------------------------------------------------------------------------------                                 
     152# Evolve system through time
     153#-------------------------------------------------------------------------------
     154
    139155import time
    140156t0 = time.time()
Note: See TracChangeset for help on using the changeset viewer.