Changeset 3854


Ignore:
Timestamp:
Oct 25, 2006, 4:54:28 PM (18 years ago)
Author:
ole
Message:

Added new mesh interface into original okushiri island validation
and did more comparisons.

Location:
anuga_validation/okushiri_2005
Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed

  • anuga_validation/okushiri_2005/create_okushiri.py

    r3850 r3854  
    11"""Create mesh and time boundary for the Okushiri island validation
    22"""
     3
    34
    45from Numeric import array, zeros, Float, allclose
    56
    67from anuga.pmesh.mesh import *
     8from anuga.pmesh.mesh_interface import create_mesh_from_regions
    79from anuga.coordinate_transforms.geo_reference import Geo_reference
    810
    … …  
    8486
    8587
    86     # Create Mesh
     88    #--------------------------------------------------------------------------
     89    # Create the triangular mesh based on overall clipping polygon with a
     90    # tagged
     91    # boundary and interior regions defined in project.py along with
     92    # resolutions (maximal area of per triangle) for each polygon
     93    #--------------------------------------------------------------------------
    8794
    88     #Basic geometry
    89    
     95    base_resolution = 10
     96
     97    # Basic geometry
    9098    xleft   = 0
    9199    xright  = 5.448
    … …  
    93101    ytop    = 3.402
    94102
    95     #Outline
     103    # Outline
    96104    point_sw = [xleft, ybottom]
    97105    point_se = [xright, ybottom]
    … …  
    99107    point_ne = [xright, ytop]
    100108
    101     #Midway points (left)
     109    # Midway points (left)
    102110    point_mtop = [xleft + (xright-xleft)/3+1, ytop]
    103111    point_mbottom = [xleft + (xright-xleft)/3+1, ybottom]
    104112
    105 
    106     geo = Geo_reference(xllcorner = xleft,
    107                         yllcorner = ybottom)
    108    
    109                        
    110     print "***********************"
    111     print "geo ref", geo
    112     print "***********************"
    113    
    114     m = Mesh(geo_reference=geo)
    115 
    116     #Boundary
    117     dict = {}
    118     dict['points'] = [point_se,   #se
    119                       point_ne,
    120                       point_mtop,
    121                       point_nw,
    122                       point_sw,
    123                       point_mbottom]
    124 
    125    
    126     dict['segments'] = [[0,1], [1,2], [2,3], [3,4],
    127                         [4,5], [5,0],  #The outer border
    128                         [2,5]]         #Separator
    129    
    130     dict['segment_tags'] = ['wall',
    131                             'wall',
    132                             'wall',
    133                             'wave',
    134                             'wall',
    135                             'wall',
    136                             '']           #Interior
    137 
    138        
    139     m.addVertsSegs(dict)
    140 
    141 
    142 
    143 
    144 
    145113    #Localised refined area for gulleys
    146     dict = {}
    147114    xl = 4.8
    148115    xr = 5.3
    … …  
    153120    p2 = [xr, yt]
    154121    p3 = [xl, yt]
    155    
    156     dict['points'] = [p0, p1, p2, p3]
    157     dict['segments'] = [[0,1], [1,2], [2,3], [3,0]]
    158     dict['segment_tags'] = ['', '', '', '']
    159     m.addVertsSegs(dict)   
     122    gulleys = [p0, p1, p2, p3]
    160123
    161124    #Island area
    … …  
    164127    island_2 = [xleft + (xright-xleft)/2+0.3, ybottom + 2*(ytop-ybottom)/3-0.3]
    165128    island_3 = [xleft + (xright-xleft)/2+0.3, ybottom + (ytop-ybottom)/3+0.3]
    166     island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3] 
     129    island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3]
    167130    island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.2]
    168131    island_6 = [xl-.01, yb]  #OK
    169132    island_7 = [xl-.01, yt]  #OK     
     133
     134    island = [island_0, island_1, island_2,
     135              island_3, island_4, island_5,
     136              island_6, island_7]
     137
     138
    170139   
    171140
    172     dict['points'] = [island_0, island_1, island_2,
    173                       island_3, island_4, island_5,
    174                       #p0, p3]                     
    175                       island_6, island_7]
     141    bounding_polygon = [point_se,
     142                        point_ne,
     143                        point_mtop,
     144                        point_nw,
     145                        point_sw,
     146                        point_mbottom]
     147   
    176148
    177                      
    178     dict['segments'] = [[0,1], [1,2], [2,3], [3,4],
    179                         [4,5], [5,6], [6,7], [7,0]]
    180                        
    181     dict['segment_tags'] = ['', '', '', '', '', '', '', '']
     149    interior_regions = [[gulleys, 0.00002*base_resolution],
     150                        [island, 0.00007*base_resolution]]
     151   
    182152
    183153
    184     m.addVertsSegs(dict)   
     154   
    185155
    186    
    187 #
    188    
    189     base_resolution = 1
     156    print 'start create mesh from regions'
    190157
    191     ocean = m.addRegionEN(xleft+.1, ybottom+.1)
    192     ocean.setMaxArea(0.1*base_resolution)
     158    meshname = project.mesh_filename + '.msh'
     159    m = create_mesh_from_regions(bounding_polygon,
     160                                 boundary_tags={'wall': [0, 1, 2, 4, 5],
     161                                                'wave': [3]},
     162                                 maximum_triangle_area=0.1*base_resolution,
     163                                 interior_regions=interior_regions,
     164                                 filename=project.mesh_filename,
     165                                 use_cache=True,
     166                                 verbose=True)
    193167
    194     mid = m.addRegionEN(point_mbottom[0]+.1, ybottom+.1)
    195     mid.setMaxArea(0.0005*base_resolution)
    196  
    197 
    198     inner = m.addRegionEN(island_3[0]+.1, island_3[1]+.1)
    199     inner.setMaxArea(0.00007*base_resolution)
    200 
    201 
    202     gulleys = m.addRegionEN(p0[0]+0.1, p0[1]+0.1)
    203     gulleys.setMaxArea(0.00002*base_resolution)
    204 
    205 
    206     # From r 1709 11 August 2005
    207     #base_resolution = 100
    208     #
    209     #ocean = m.addRegionEN(xleft+.1, ybottom+.1)
    210     #ocean.setMaxArea(0.01*base_resolution)
    211     #
    212     #mid = m.addRegionEN(point_mbottom[0]+.1, ybottom+.1)
    213     #mid.setMaxArea(0.001*base_resolution)
    214 
    215     #inner = m.addRegionEN(island_3[0]+.1, island_3[1]+.1)
    216     #inner.setMaxArea(0.0001*base_resolution)
    217 
    218 
    219     #gulleys = m.addRegionEN(p0[0]+0.1, p0[1]+0.1)
    220     #gulleys.setMaxArea(0.00001*base_resolution)       
    221    
    222    
    223     m.generateMesh('pzq28.0za1000000a')
    224 
    225     import project
    226     m.export_mesh_file(project.mesh_filename)
    227    

  • anuga_validation/okushiri_2005/extract_timeseries.py

    r3850 r3854  
    2121    plotting = True
    2222
     23plotting = False
    2324
    2425#-------------------------
    … …  
    4243                        'ch7': 1.121816242516758e-04,
    4344                        'ch9': 1.249543278366778e-04}
     45
     46# old limiters
     47expected_covariances = {'Boundary': 5.288601162783020386e-05,
     48                        'ch5': 1.167001054284431472e-04,
     49                        'ch7': 1.121474766904651861e-04,
     50                        'ch9': 1.249244820847215335e-04}
    4451
    4552
    … …  
    8188# Read and interpolate model output
    8289#--------------------------------------------------
    83 #f = file_function('okushiri_new_limiters.sww',
     90#f = file_function('okushiri_new_limiters.sww',   #The best so far
    8491#f = file_function('okushiri_as2005_with_mxspd=0.1.sww',
    8592f = file_function(project.output_filename,
    … …  
    111118        print 'Result is better than expected by: %.18e'\
    112119              %(res-expected_covariances[name])
    113         print 'Result = %.18e' %res
    114120        print 'Expect = %.18e' %expected_covariances[name]   
    115121    elif res > expected_covariances[name]+eps:
    116         print 'FAIL: %.18e' %res
     122        print 'FAIL: Result is worse than expected by: %.18e'\
     123                                %(res-expected_covariances[name])
     124        print 'Expect = %.18e' %expected_covariances[name]
    117125    else:
    118126        pass

  • anuga_validation/okushiri_2005/project.py

    r3850 r3854  
    1515
    1616# Model output
    17 output_filename = 'okushiri_old_limiters.sww'
     17#output_filename = 'okushiri_old_limiters.sww'
     18output_filename = 'okushiri_nl_new_mesh.sww'
    1819
    1920

  • anuga_validation/okushiri_2005/run_okushiri.py

    r3850 r3854  
    5353domain.set_name(project.output_filename)  # Name of output sww file
    5454domain.set_default_order(2)               # Apply second order scheme
    55 domain.set_all_limiters(0.9)              # Maximal second order scheme (old lim)
     55#domain.set_all_limiters(0.9)              # Maximal second order scheme (old lim)
    5656domain.set_minimum_storable_height(0.001) # Don't store w < 0.001m
    5757domain.set_maximum_allowed_speed(0.1)     # Allow a little runoff (0.1 is OK)
Note: See TracChangeset for help on using the changeset viewer.