Changeset 1832


Ignore:
Timestamp:
Sep 14, 2005, 6:22:00 PM (19 years ago)
Author:
ole
Message:

First run of karratha with MOST. New mesh

Location:
production/karratha_2005
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • production/karratha_2005/create_mesh.py

    r1830 r1832  
    66#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    77
    8 from pmesh.mesh import *
     8
    99from pyvolution.coordinate_transforms.geo_reference import Geo_reference
    1010from pyvolution.coordinate_transforms.redfearn import redfearn
    1111
    1212
    13 import project
    1413
    15 #-------------------------------------------------------------
    16 if __name__ == "__main__":
     14def create_mesh_boundary(polygon, tags, resolution, filename = None):
     15    """Create mesh from bounding polygon, tags for all segments and resolution.
     16
     17    Polygon is a list of points in latitude and longitudes - decimal degrees
    1718
    1819
     20    Tags is a list of symbolic tags
     21
     22    Resolution is the maximal area
    1923   
    20     #Basic geometry
    2124
    22     south = project.south
    23     north = project.north
    24     west = project.west
    25     east = project.east
     25    FIXME - FUTURE:
     26      Use class Point
     27      Take multiple polygons
    2628
    27     refzone = project.refzone
     29      Accept deg, min, sec, e.g.
     30      [ [(-20,30,55), (116, 20, 00)], ...]
     31     
     32   
     33    """
     34
     35    from pmesh.mesh import *
     36    from pyvolution.coordinate_transforms.redfearn import redfearn
     37   
     38
     39    #Convert to UTM
     40    import project
     41    refzone = project.refzone  #FIXME
    2842
    2943
    30     #SW
    31     zone, easting, northing  = redfearn(south, west)
    32     assert zone == refzone
    33     point_sw = [easting, northing]
     44    points = []
     45    for point in polygon:
     46        zone, easting, northing  = redfearn(point[0], point[1]) #FIXME: Use point.latitude etc
     47        assert zone == refzone
     48        points.append([easting, northing])
    3449
    35     #SE
    36     zone, easting, northing  = redfearn(south, east)
    37     assert zone == refzone
    38     point_se = [easting, northing]   
     50    polygon = points   
    3951
    40     #NW
    41     zone, easting, northing  = redfearn(north, west)
    42     assert zone == refzone
    43     point_nw = [easting, northing]
     52   
     53    mesh_origin = project.mesh_origin
    4454
    45     #NE
    46     zone, easting, northing  = redfearn(north, east)
    47     assert zone == refzone
    48     point_ne = [easting, northing]   
     55    geo = Geo_reference(xllcorner = mesh_origin[1],  #From dem
     56                        yllcorner = mesh_origin[2],
     57                        zone = refzone)
     58   
     59
     60    #
     61    dict = {}
     62    dict['points'] = polygon
     63
     64    #Create segments
     65    #E.g. [[0,1], [1,2], [2,3], [3,0]]
     66    segments = []
     67    N = len(polygon)
     68    for i in range(N):
     69        lo = i
     70        hi = (lo + 1) % N
     71        segments.append( [lo, hi] )
     72       
     73       
     74    dict['segments'] = segments
     75
     76
     77    #Create tags
     78    #E.g. ['wall', 'wall', 'ocean', 'wall']
     79    segment_tags = [0]*N
     80    for key in tags:
     81        indices = tags[key]
     82        for i in indices:
     83            segment_tags[i] = key
    4984   
     85    dict['segment_tags'] = segment_tags
    5086
    51     #geo = Geo_reference(xllcorner = 421544.35127423,  #from dem
    52     #                    yllcorner = 7677669.5257159,
    53     #                    zone = refzone)
    5487
    55     geo = Geo_reference(xllcorner = 420468.31429902,  #From dem
    56                         yllcorner = 7677669.5257159,
    57                         zone = refzone)
    5888
    59     #geo = Geo_reference(xllcorner = point_sw[0]+1000,  #FIXME: BAD WHY
    60     #                    yllcorner = point_sw[1]+1000,
    61     #                    zone = refzone)       
    62    
    63                        
    6489    print "***********************"
    6590    print "geo ref", geo
    … …  
    6792   
    6893    m = Mesh(geo_reference=geo)
    69 
    70     #Boundary
    71     dict = {}
    72     dict['points'] = [point_sw,
    73                       point_se,
    74                       point_ne,
    75                       point_nw]
    76 
    77    
    78     dict['segments'] = [[0,1], [1,2], [2,3], [3,0]]
    79     dict['segment_tags'] = ['wall', 'wall', 'ocean', 'wall']
    8094       
    8195    m.addVertsSegs(dict)
    8296
    8397
    84     base_resolution = 50000
     98    #FIXME: Find centroid automatically
     99    inner = m.addRegionEN(points[0][0]+1, points[0][1]+1)
     100    inner.setMaxArea(resolution)
     101   
     102    m.generateMesh('pzq28.0za1000000a')
     103   
     104    m.export_mesh_file(filename)
     105   
    85106
    86     inner = m.addRegionEN(point_sw[0]+1, point_sw[1]+1)
    87     inner.setMaxArea(base_resolution)
    88107
    89     m.generateMesh('pzq28.0za1000000a')
    90108
    91     m.export_mesh_file(project.meshname + '.msh')
     109
    92110   

  • production/karratha_2005/project.py

    r1830 r1832  
    3232from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees
    3333
     34#Origin of existing dem (FIXME: Temporary measure)
     35#mesh_origin = (50, 421544.35127423, 7677669.5257159)  #250m
     36mesh_origin = (50, 420468.31429902, 7677669.5257159)  #100m
     37
     38
     39#south = degminsec2decimal_degrees(-20,45,0)
     40#north = degminsec2decimal_degrees(-20,15,0)
     41#west = degminsec2decimal_degrees(116,30,0)
     42#east = degminsec2decimal_degrees(117,0,0)
     43
    3444south = degminsec2decimal_degrees(-20,45,0)
    3545north = degminsec2decimal_degrees(-20,15,0)
    36 west = degminsec2decimal_degrees(116,30,0)
    37 east = degminsec2decimal_degrees(117,0,0)
     46west = degminsec2decimal_degrees(116,20,0)
     47east = degminsec2decimal_degrees(117,10,0)
     48
     49p0 = [south, degminsec2decimal_degrees(116,32,0)]
     50p1 = [south, west]
     51p2 = [degminsec2decimal_degrees(-20,23,0), west]
     52p3 = [north, degminsec2decimal_degrees(116,45,0)]
     53p4 = [north, degminsec2decimal_degrees(117,0,0)]
     54p5 = [p2[0], degminsec2decimal_degrees(117,8,0)]
     55p6 = [degminsec2decimal_degrees(-20,30,0), east]
     56p7 = [degminsec2decimal_degrees(-20,38,0), east]
     57p8 = [south, east]
     58   
     59polygon = [p0, p1, p2, p3, p4, p5, p6, p7, p8]
    3860
    3961refzone = 50

  • production/karratha_2005/run_karratha.py

    r1830 r1832  
    77
    88
    9 from pyvolution.shallow_water import Domain, Reflective_boundary, File_boundary, Dirichlet_boundary, Time_boundary
     9from pyvolution.shallow_water import Domain, Reflective_boundary, File_boundary, Dirichlet_boundary, Time_boundary, Transmissive_boundary
    1010from pyvolution.data_manager import convert_dem_from_ascii2netcdf,\
    1111     dem2pts, ferret2sww
    … …  
    4141east = project.east
    4242
    43 #Origin of existing dem (FIXME: Temporary measure)
    44 #mesh_origin = (50, 421544.35127423, 7677669.5257159)  #250m
    45 mesh_origin = (50, 420468.31429902, 7677669.5257159)  #100m
     43
     44#Mesh
     45from create_mesh import create_mesh_boundary
     46
     47m = cache(create_mesh_boundary,
     48          project.polygon,
     49          {'tags': {'back': [7, 8], 'side': [0, 6], 'ocean': [1, 2, 3, 4, 5]},
     50          'resolution': 52000,
     51          'filename': project.meshname + '.msh'},
     52          verbose = True)
     53
     54
     55
     56
     57
     58          #tags = {'back': [7, 8], 'side': [0, 6], 'ocean': [1, 2, 3, 4, 5]},
     59          #resolution = 200000,
     60          #filename = project.meshname + '.msh')
     61
     62
     63
    4664
    4765ferret2sww(source_dir+project.boundary_basename,
    … …  
    5068           minlat=south-1, maxlat=north+1,
    5169           minlon=west-1, maxlon=east+1,
    52            origin = mesh_origin,
     70           origin = project.mesh_origin,
    5371           zscale = 1,
    5472           fail_on_NaN = False,
    … …  
    6886domain.store = True
    6987
    70 print "Number of triangles = ", len(domain)
     88print 'Number of triangles = ', len(domain)
    7189print 'The extent is ', domain.get_extent()
    7290
    7391
    7492#Setup IC
    75 domain.set_quantity('stage', 0)
     93tide = 0
     94domain.set_quantity('stage', tide)
    7695domain.set_quantity('elevation',
    7796                    filename = demname + '.pts',
    … …  
    83102print domain.get_boundary_tags()
    84103
    85 #Bf = File_boundary(project.boundary_basename + '.sww', domain, verbose = True)
     104Bf = File_boundary(project.boundary_basename + '.sww', domain, verbose = True)
     105domain.starttime = 3000  #Obtained from MOST
     106
    86107
    87108Br = Reflective_boundary(domain)
    88 Bd = Dirichlet_boundary([3,0,0])
     109Bt = Transmissive_boundary(domain)
     110Bd = Dirichlet_boundary([tide,0,0])
    89111Bw = Time_boundary(domain=domain,
    90                    f=lambda t: [(60<t<660)*12, 0, 0])
    91 domain.set_boundary( {'wall': Br, 'ocean': Bw} )
     112                   f=lambda t: [(60<t<660)*4, 0, 0])
     113
     114#domain.set_boundary( {'wall': Bd, 'ocean': Bf} )
     115domain.set_boundary( {'back': Br,'side': Bd, 'ocean': Bf} )
    92116
    93117
    94118#Run
    95 for t in domain.evolve(yieldstep = 15, finaltime = 4800):
     119
     120for t in domain.evolve(yieldstep = 60, finaltime = 28610):
    96121    domain.write_time()
    97     domain.write_boundary_statistics(quantities = 'stage')
     122    domain.write_boundary_statistics(tags = 'ocean') #quantities = 'stage')
Note: See TracChangeset for help on using the changeset viewer.