Changeset 3883


Ignore:
Timestamp:
Oct 30, 2006, 11:41:12 AM (16 years ago)
Author:
ole
Message:

More cleanup of okushiri prior to release

Location:
anuga_validation/okushiri_2005
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • anuga_validation/okushiri_2005/create_okushiri.py

    r3878 r3883  
    88from anuga.pmesh.mesh_interface import create_mesh_from_regions
    99from anuga.coordinate_transforms.geo_reference import Geo_reference
     10from anuga.geospatial_data import Geospatial_data
    1011
    1112import project
    1213
     14def prepare_bathymetry(filename):
     15    """Convert benchmark 2 bathymetry to NetCDF pts file.
     16    This is a 'throw-away' code taylor made for files like
     17    'Benchmark_2_bathymetry.txt' from the LWRU2 benchmark
     18    """
     19
     20    print 'Creating', filename
     21   
     22    # Read the ascii (.txt) version of this file,
     23    # make it comma separated and invert the bathymetry
     24    # (Below mean sea level should be negative)
     25    infile = open(filename[:-4] + '.txt')
     26
     27    points = []
     28    attribute = []
     29    for line in infile.readlines()[1:]: #Skip first line (the header)
     30        fields = line.strip().split()
     31
     32        x = float(fields[0])
     33        y = float(fields[1])
     34        z = -float(fields[2]) # Bathymetry is inverted in original file
     35       
     36        points.append([x,y])
     37        attribute.append(z)
     38    infile.close()
     39
     40    # Convert to geospatial data and store as NetCDF
     41    G = Geospatial_data(data_points=points,
     42                        attributes=attribute)
     43    G.export_points_file(filename)
     44   
     45
    1346
    1447def prepare_timeboundary(filename):
    15     """Converting benchmark 2 time series to NetCDF tms file.
     48    """Convert benchmark 2 time series to NetCDF tms file.
    1649    This is a 'throw-away' code taylor made for files like
    1750    'Benchmark_2_input.txt' from the LWRU2 benchmark
    1851    """
    1952
    20     textversion = filename[:-4] + '.txt'
    21 
    22     print 'Preparing time boundary from %s' %textversion
     53    from Scientific.IO.NetCDF import NetCDFFile
    2354    from Numeric import array
    2455
    25     fid = open(textversion)
    2656
    27     #Skip first line
     57    print 'Creating', filename
     58
     59    # Read the ascii (.txt) version of this file
     60    fid = open(filename[:-4] + '.txt')
     61
     62    # Skip first line
    2863    line = fid.readline()
    2964
    30     #Read remaining lines
     65    # Read remaining lines
    3166    lines = fid.readlines()
    3267    fid.close()
     
    4479
    4580
    46     #Create tms file
    47     from Scientific.IO.NetCDF import NetCDFFile
     81    # Create tms NetCDF file
    4882
    49     print 'Writing to', filename
    5083    fid = NetCDFFile(filename, 'w')
    51 
    5284    fid.institution = 'Geoscience Australia'
    5385    fid.description = 'Input wave for Benchmark 2'
     
    73105
    74106
    75     #Preparing points
    76     #(Only when using original Benchmark_2_Bathymetry.txt file)
    77     # FIXME: Replace by using geospatial data
    78     #
    79     #from anuga.abstract_2d_finite_volumes.data_manager import xya2pts
    80     #xya2pts(project.bathymetry_filename, verbose = True,
    81     #        z_func = lambda z: -z)
    82 
     107    # Prepare bathymetry
     108    prepare_bathymetry(project.bathymetry_filename)
    83109
    84110    # Prepare time boundary
     
    94120
    95121
    96     base_resolution = 1
     122    base_resolution = 1 # Use this to coarsen or refine entire mesh.
    97123
    98     # Basic geometry
     124    # Basic geometry and bounding polygon
    99125    xleft   = 0
    100126    xright  = 5.448
     
    102128    ytop    = 3.402
    103129
    104     # Outline
    105130    point_sw = [xleft, ybottom]
    106131    point_se = [xright, ybottom]
     
    108133    point_ne = [xright, ytop]
    109134
    110     # Midway points (left)
    111     point_mtop = [xleft + (xright-xleft)/3+1, ytop]
    112     point_mbottom = [xleft + (xright-xleft)/3+1, ybottom]
     135    bounding_polygon = [point_se,
     136                        point_ne,
     137                        point_nw,
     138                        point_sw]
     139   
    113140
    114     #Localised refined area for gulleys
     141    # Localised refined area for gulleys
    115142    xl = 4.8
    116143    xr = 5.3
     
    121148    p2 = [xr, yt]
    122149    p3 = [xl, yt]
     150   
    123151    gulleys = [p0, p1, p2, p3]
     152   
    124153
    125     #Island area
     154    # Island area and drawdown region
    126155    island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5]
    127156    island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3]
     
    132161    island_6 = [xl-.01, yb]  # Keep right edge just off the gulleys
    133162    island_7 = [xl-.01, yt]
    134 
    135    
    136 
     163 
    137164    island = [island_0, island_1, island_2,
    138165              island_3, island_4, island_5,
     
    140167
    141168
    142     #Midway points just inside boundary and associated rectangular region
    143     point_mtop0 = [xleft + (xright-xleft)/3+1, ytop-0.02]
    144     point_mbottom0 = [xleft + (xright-xleft)/3+1, ybottom+0.02]
    145     point_se0 = [xright-0.02, ybottom+0.02]
    146     point_ne0 = [xright-0.02, ytop-0.02]   
     169    # Region spanning half right hand side of domain just inside boundary
     170    rhs_nw = [xleft + (xright-xleft)/3+1, ytop-0.02]
     171    rhs_sw = [xleft + (xright-xleft)/3+1, ybottom+0.02]
     172    rhs_se = [xright-0.02, ybottom+0.02]
     173    rhs_ne = [xright-0.02, ytop-0.02]   
    147174
    148     mid = [point_mtop0, point_ne0, point_se0, point_mbottom0]
     175    rhs_region = [rhs_nw, rhs_ne, rhs_se, rhs_sw]
    149176
    150177
    151     # Bounding polygon and interior regions
    152     bounding_polygon = [point_se,
    153                         point_ne,
    154                         point_mtop,
    155                         point_nw,
    156                         point_sw,
    157                         point_mbottom]
    158178   
    159 
    160     interior_regions = [[mid, 0.0005],
     179    # Interior regions and creation of mesh
     180    interior_regions = [[rhs_region, 0.0005],
    161181                        [gulleys, 0.00002*base_resolution],
    162182                        [island, 0.00007*base_resolution]]
    163183
    164    
    165 
    166     print 'start create mesh from regions'
    167184    meshname = project.mesh_filename + '.msh'
    168185    m = create_mesh_from_regions(bounding_polygon,
    169                                  boundary_tags={'wall': [0, 1, 2, 4, 5],
    170                                                 'wave': [3]},
     186                                 boundary_tags={'wall': [0, 1, 3],
     187                                                'wave': [2]},     
    171188                                 maximum_triangle_area=0.1*base_resolution,
    172189                                 interior_regions=interior_regions,
  • anuga_validation/okushiri_2005/extract_timeseries.py

    r3878 r3883  
    5656
    5757
    58 expected_covariances = {'Boundary': 5.288392008865989237e-05,
    59                         'ch5': 1.166748190444680592e-04,
    60                         'ch7': 1.121816242516757850e-04,
    61                         'ch9': 1.249543278366777640e-04}
     58#expected_covariances = {'Boundary': 5.288392008865989237e-05,
     59#                        'ch5': 1.166748190444680592e-04,
     60#                        'ch7': 1.121816242516757850e-04,
     61#                        'ch9': 1.249543278366777640e-04}
     62#
     63#expected_differences = {'Boundary':  8.373150808730501615e-04,
     64#                        'ch5': 3.425914311580337875e-03,
     65#                        'ch7': 2.802327594773105189e-03,
     66#                        'ch9': 3.198733498646373370e-03}
    6267
    63 expected_differences = {'Boundary':  8.373150808730501615e-04,
    64                         'ch5': 3.425914311580337875e-03,
    65                         'ch7': 2.802327594773105189e-03,
    66                         'ch9': 3.198733498646373370e-03}
     68
     69expected_covariances = {'Boundary': 5.278365381306500051e-05,
     70                        'ch5': 1.168250251701857399e-04,
     71                        'ch7': 1.125146958894896555e-04,
     72                        'ch9': 1.248634476898675049e-04}
     73
     74expected_differences = {'Boundary': 8.502104901511024380e-04,
     75                        'ch5': 3.442302747303998579e-03,
     76                        'ch7': 2.802987591027187534e-03,
     77                        'ch9': 3.191368834008508938e-03}
     78
     79
    6780
    6881
  • anuga_validation/okushiri_2005/project.py

    r3882 r3883  
    11"""Common filenames for Okushiri Island validation
     2Formats are given as ANUGA native netCDF where applicable.
     3
    24"""
    35
     
    810validation_filename = 'output_ch5-7-9.txt'
    911
    10 # A simple conversion from txt file to the more compact NetCDF format
     12# Digital Elevation Model
    1113bathymetry_filename = 'Benchmark_2_Bathymetry.pts'
    1214
     
    1517
    1618# Model output
    17 #output_filename = 'okushiri_old_limiters.sww'
    18 #output_filename = 'okushiri_new_limiters.sww'
    19 #output_filename = 'okushiri_nl_coarse10_mesh.sww'
    20 #output_filename = 'okushiri_nl_coarse100_mesh.sww'
    21 #output_filename = 'okushiri_original_mesh.sww'
    22 #utput_filename = 'okushiri_new_mesh.sww'
    23 output_filename = 'okushiri_new_mesh1.sww'
     19output_filename = 'okushiri.sww'
    2420
     21
     22
  • anuga_validation/okushiri_2005/run_okushiri.py

    r3882 r3883  
    3838print domain.statistics()
    3939
     40
    4041#-------------------------
    4142# Initial Conditions
     
    4950                    use_cache=True)
    5051
     52
    5153#-------------------------
    52 # Distribute domain
     54# Distribute domain if run in parallel
    5355#-------------------------
    54 domain = distribute(domain)
     56if numprocs > 1:
     57    domain = distribute(domain)
    5558
    56 # Parameters
     59
     60#-------------------------
     61# Set simulation parameters
     62#-------------------------
    5763domain.set_name(project.output_filename)  # Name of output sww file
    5864domain.set_default_order(2)               # Apply second order scheme
Note: See TracChangeset for help on using the changeset viewer.