Changeset 1742 for inundation/validation


Ignore:
Timestamp:
Aug 23, 2005, 2:56:21 PM (20 years ago)
Author:
ole
Message:

Recovered validation stuff from 18 August + wrapping up

Location:
inundation/validation/LWRU2
Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed

  • inundation/validation/LWRU2/create_mesh.py

    r1734 r1742  
    1414
    1515    #Basic geometry
    16 
     16   
    1717    xleft   = 0
    1818    xright  = 5.448
    … …  
    2020    ytop    = 3.402
    2121
    22     #Outline
     22    #Outline 
    2323    point_sw = [xleft, ybottom]
    2424    point_se = [xright, ybottom]
    25     point_nw = [xleft, ytop]
     25    point_nw = [xleft, ytop]   
    2626    point_ne = [xright, ytop]
    2727
    2828    #Midway points (left)
    29     point_ml1 = [xleft + (xright-xleft)/3+1, ytop]
    30     point_ml2 = [xleft + (xright-xleft)/3+1, ybottom]
    31 
    32     #Midway points (right)
    33     point_mr1 = [xleft + 2*(xright-xleft)/3+0.5, ytop]
    34     point_mr2 = [xleft + 2*(xright-xleft)/3+0.5, ybottom]
     29    point_mtop = [xleft + (xright-xleft)/3+1, ytop]
     30    point_mbottom = [xleft + (xright-xleft)/3+1, ybottom]
    3531
    3632
    3733    geo = Geo_reference(xllcorner = xleft,
    3834                        yllcorner = ybottom)
    39 
    40 
     35   
     36                       
    4137    print "***********************"
    4238    print "geo ref", geo
    4339    print "***********************"
    44 
     40   
    4541    m = Mesh(geo_reference=geo)
    4642
    … …  
    4945    dict['points'] = [point_se,   #se
    5046                      point_ne,
    51                       point_mr1,
    52                       point_ml1,
     47                      point_mtop,
    5348                      point_nw,
    5449                      point_sw,
    55                       point_ml2,
    56                       point_mr2]
     50                      point_mbottom]
    5751
    58 
     52   
    5953    dict['segments'] = [[0,1], [1,2], [2,3], [3,4],
    60                         [4,5], [5,6], [6,7], [7,0], #The outer border
    61                         [2,7], [3,6]]               #Separators
    62 
     54                        [4,5], [5,0], #The outer border
     55                        [2,5]]         #Separator
     56   
    6357    dict['segment_tags'] = ['wall',
    64                             'wall',
    6558                            'wall',
    6659                            'wall',
    … …  
    6861                            'wall',
    6962                            'wall',
    70                             'wall',
    71                             '',           #Interior
    7263                            '']           #Interior
    7364
     65       
     66    m.addVertsSegs(dict)
    7467
    75     m.addVertsSegs(dict)
     68
    7669
    7770
    … …  
    8780    p2 = [xr, yt]
    8881    p3 = [xl, yt]
     82   
     83    dict['points'] = [p0, p1, p2, p3]
     84    dict['segments'] = [[0,1], [1,2], [2,3], [3,0]]
     85    dict['segment_tags'] = ['', '', '', '']
     86    m.addVertsSegs(dict)   
    8987
    90     dict['points'] = [p0, p1, p2, p3]
     88    #Island area
     89    island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5]
     90    island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3]
     91    island_2 = [xleft + (xright-xleft)/2+0.3, ybottom + 2*(ytop-ybottom)/3-0.3]
     92    island_3 = [xleft + (xright-xleft)/2+0.3, ybottom + (ytop-ybottom)/3+0.3]
     93    island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3]
     94    island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.2]
     95    island_6 = [xl-.01, yb]  #OK
     96    island_7 = [xl-.01, yt]  #OK     
     97   
     98
     99    dict['points'] = [island_0, island_1, island_2,
     100                      island_3, island_4, island_5,
     101                      #p0, p3]                     
     102                      island_6, island_7]
     103
     104                     
     105    dict['segments'] = [[0,1], [1,2], [2,3], [3,4],
     106                        [4,5], [5,6], [6,7], [7,0]]
     107                       
     108    dict['segment_tags'] = ['', '', '', '', '', '', '', '']
    91109
    92110
    93     dict['segments'] = [[0,1], [1,2], [2,3], [3,0]]
    94     dict['segment_tags'] = ['', '', '', '']
    95     m.addVertsSegs(dict)
     111    m.addVertsSegs(dict)   
    96112
     113   
     114#
     115   
     116    base_resolution = 1
    97117
    98     base_resolution = 100
     118    ocean = m.addRegionEN(xleft+.1, ybottom+.1)
     119    ocean.setMaxArea(0.1*base_resolution)
    99120
    100     ocean = m.addRegionEN(xleft+1, ybottom+1)
    101     ocean.setMaxArea(0.01*base_resolution)
     121    mid = m.addRegionEN(point_mbottom[0]+.1, ybottom+.1)
     122    mid.setMaxArea(0.0005*base_resolution)
     123   
    102124
    103     mid = m.addRegionEN(point_ml2[0]+1, ybottom+1)
    104     mid.setMaxArea(0.001*base_resolution)
    105 
    106 
    107     inner = m.addRegionEN(point_mr2[0]+1, ybottom+1)
    108     inner.setMaxArea(0.0001*base_resolution)
     125    inner = m.addRegionEN(island_3[0]+.1, island_3[1]+.1)
     126    inner.setMaxArea(0.00007*base_resolution)
    109127
    110128
    111129    gulleys = m.addRegionEN(p0[0]+0.1, p0[1]+0.1)
    112     gulleys.setMaxArea(0.00001*base_resolution)
    113 
    114 
     130    gulleys.setMaxArea(0.00002*base_resolution)   
     131   
     132   
    115133    m.generateMesh('pzq28.0za1000000a')
    116134
    117135    import filenames
    118136    m.export_mesh_file(filenames.mesh_filename)
     137   

  • inundation/validation/LWRU2/extract_timeseries.py

    r1736 r1742  
    99gauges = [[0.000, 1.696]]  #Boundary
    1010gauges += [[4.521, 1.196],  [4.521, 1.696],  [4.521, 2.196]] #Ch 5-7-9
    11 #gauges = [[4.521, 1.196],  [4.521, 1.696],  [4.521, 2.196]] #Ch 5-7-9
    1211
    1312gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9']
    14 #gauge_names = ['ch5', 'ch7', 'ch9']
    1513
    1614

  • inundation/validation/LWRU2/lwru2.py

    r1736 r1742  
    7777from Numeric import array, zeros, Float, allclose
    7878import filenames
     79from caching import cache
    7980
    80 picklefile = 'domain.pck'
    81 from cPickle import dump, load
     81#Preparing time boundary
     82prepare_timeboundary(filenames.boundary_filename)
    8283
    83 N = 100
    84 
    85 try:
    86     #raise Exception
    87     fid = open(picklefile)
    88     print 'Read pickled domain'
    89     domain = load(fid)
    90 
    91 except:
    92 
    93     #Preparing time boundary
    94     prepare_timeboundary(filenames.boundary_filename)
    95 
    96     #Preparing points
    97     from pyvolution.data_manager import xya2pts
    98     xya2pts(filenames.bathymetry_filename, verbose = True,
    99             z_func = lambda z: -z)
     84#Preparing points
     85from pyvolution.data_manager import xya2pts
     86xya2pts(filenames.bathymetry_filename, verbose = True,
     87        z_func = lambda z: -z)
    10088
    10189
    102     #######################
    103     # Domain
     90#######################
     91# Domain
    10492   
    105     if read_mesh is True:
    106         print 'Creating domain from', filenames.mesh_filename
    107         domain = pmesh_to_domain_instance(filenames.mesh_filename, Domain)   
    108     else:       
    109         print 'Creating regular domain'
    110         points, vertices, boundary = rectangular(N, N/5*3,
    111                                                  len1=5.448, len2=3.402)
    112         domain = Domain(points, vertices, boundary)
     93if read_mesh is True:
     94    print 'Creating domain from', filenames.mesh_filename
     95    domain = pmesh_to_domain_instance(filenames.mesh_filename, Domain)   
     96else:   
     97    print 'Creating regular mesh'
     98    N = 100
     99    points, vertices, boundary = rectangular(N, N/5*3,
     100                                             len1=5.448, len2=3.402)
     101    print 'Creating domain'
     102   
     103    #domain = Domain(points, vertices, boundary)
     104    domain = cache(Domain, (points, vertices, boundary))
    113105
    114106
    115107
    116     print "Number of triangles = ", len(domain)
    117     print 'The extent is ', domain.get_extent()
     108print "Number of triangles = ", len(domain)
     109print 'The extent is ', domain.get_extent()
    118110
    119111
    120     domain.check_integrity()
    121     domain.default_order = 2
     112#domain.check_integrity()
     113domain.default_order = 2
    122114
    123     print "Number of triangles = ", len(domain)
    124     domain.store = True    #Store for visualisation purposes
     115print "Number of triangles = ", len(domain)
     116domain.store = True    #Store for visualisation purposes
    125117
    126     import sys, os
    127     base = os.path.basename(sys.argv[0])
    128     domain.filename, _ = os.path.splitext(base)
     118import sys, os
     119base = os.path.basename(sys.argv[0])
     120domain.filename, _ = os.path.splitext(base)
    129121
    130     #LS code to be included in set_quantity
    131     from pyvolution import util, least_squares
    132     points, attributes = util.read_xya(filenames.bathymetry_filename[:-4] + '.pts')
     122#LS code to be included in set_quantity
     123print 'Reading points'
     124from pyvolution import util, least_squares
     125points, attributes = util.read_xya(filenames.bathymetry_filename[:-4] + '.pts')
     126elevation = attributes['elevation']
    133127
    134     #Fit attributes to mesh
    135     vertex_attributes = least_squares.fit_to_mesh(domain.coordinates,
    136                                                   domain.triangles,
    137                                                   points,
    138                                                   attributes['elevation'],
    139                                                   alpha = 0.0001,
    140                                                   verbose=True)
     128#Fit attributes to mesh
     129#vertex_attributes = least_squares.fit_to_mesh(domain.coordinates,
     130#                                              domain.triangles,
     131#                                              points,
     132#                                              attributes['elevation'],
     133#                                              alpha = 0.0001,
     134#                                              verbose=True)
    141135
    142     print 'Initial values'
    143     domain.set_quantity('elevation', vertex_attributes)
    144     domain.set_quantity('friction', 0.0)
    145     domain.set_quantity('stage', 0.0)
     136#vertex_attributes = cache(least_squares.fit_to_mesh,
     137#                          (domain.coordinates, domain.triangles,
     138#                           points, elevation) ,
     139#                          {'alpha': 0.0001, 'verbose': True})
    146140
     141print 'Initial values'
     142#domain.set_quantity('elevation', vertex_attributes)
     143domain.set_quantity('elevation', points=points, attributes=attributes['elevation'], alpha=0.0001)
     144domain.set_quantity('friction', 0.0)
     145domain.set_quantity('stage', 0.0)
    147146
    148     fid = open('domain.pck', 'w')
    149     dump(domain, fid)
    150     fid.close()
    151147
    152148
    … …  
    166162
    167163#Set boundary conditions
    168 
    169164if read_mesh is True:
    170165    domain.set_boundary({'wave': Bts, 'wall': Br})     
Note: See TracChangeset for help on using the changeset viewer.