source: anuga_work/development/Rudy_vandrie_2007/create_example_version1.py @ 4440

Last change on this file since 4440 was 4440, checked in by ole, 17 years ago

Examples prepared for Rudy

File size: 6.5 KB
Line 
1"""Create mesh and time boundary for Hydrograph example
2"""
3
4
5from Numeric import array, zeros, Float, allclose
6
7from anuga.pmesh.mesh import *
8from anuga.pmesh.mesh_interface import create_mesh_from_regions
9from anuga.coordinate_transforms.geo_reference import Geo_reference
10from anuga.geospatial_data import Geospatial_data
11
12import project
13
14
15def prepare_timeboundary(filename, hydrograph, field):
16    """Convert .out file to NetCDF tms file.
17    This is a 'throw-away' code taylor made for this type of file
18
19    hydrograph is an identifier e.g.
20    HYDROGRAPHS_SUB17       
21
22    field is one of the column headers.
23
24    Time will always be read and converted into seconds.
25
26    The file has contents like
27
28    #####START_HYDROGRAPHS_SUB17       
29     Time    Rain  Rainperv    Qtop      Qbot      Qper      Qimp    Qinto_OS   Qout_OS  Stage
30      0.0    0.00    0.00     0.000     0.000     0.000     0.000     0.000     0.000     0.000
31      5.0    2.65    0.00     0.000     0.000     0.000     0.004     0.004     0.004     0.000
32     10.0    2.65    0.00     0.000     0.000     0.000     0.004     0.004     0.004     0.000
33     15.0    2.65    0.00     0.000     0.000     0.000     0.004     0.004     0.004     0.000
34
35   
36    """
37
38    from Scientific.IO.NetCDF import NetCDFFile
39    from Numeric import array
40
41    assert filename[-4:] == '.tms'
42   
43    outfilename = filename
44    infilename = filename[:-4] + '.out'
45
46    print 'Creating', outfilename
47
48    # Read the ascii (.txt) version of this file
49    fid = open(infilename)
50
51    # Read all lines and search for selected hydrograph
52    lines = fid.readlines()
53    fid.close()
54
55    found = False
56    search_string = 'START_'+hydrograph
57    print 'looking for', search_string
58    for i, line in enumerate(lines):
59        if line.startswith('####'):
60            if line.find(search_string) >= 0:
61                print 'Found', line
62                found = True
63                break
64
65    if found is False:
66        msg = 'Did not find', hydrograph
67        raise Exception, msg
68
69    # Select data column
70    headers = lines[i+1]
71    fields = headers.split()
72    index = fields.index(field)
73    print 'Found header #%d: "%s"' %(index, fields[index])
74
75
76    # Read data
77    search_string = 'END_'+hydrograph
78    time = []
79    data = []   
80    for line in lines[i+2:]:
81        if line.startswith('####'):
82            if line.find(search_string) >= 0:
83                break       
84       
85        fields = line.split()
86        time.append(float(fields[0]))
87        data.append(float(fields[index]))
88
89
90    # Convert to NetCDF
91    N = len(time)
92    T = array(time, Float)*60  # Time (seconds)
93    Q = array(data, Float)     # Values (m^3/s)
94
95    inflow_stage = 0           # Assumed stage at inflow
96    boundary_width = 3.402     # Use width to convert flow into x-momentum
97
98    # Create tms NetCDF file
99
100    fid = NetCDFFile(outfilename, 'w')
101    fid.institution = 'Geoscience Australia'
102    fid.description = 'Input wave for Benchmark 2'
103    fid.starttime = 0.0
104    fid.createDimension('number_of_timesteps', len(T))
105    fid.createVariable('time', Float, ('number_of_timesteps',))
106    fid.variables['time'][:] = T
107
108    fid.createVariable('stage', Float, ('number_of_timesteps',))
109    fid.variables['stage'][:] = inflow_stage
110
111    fid.createVariable('xmomentum', Float, ('number_of_timesteps',))
112    fid.variables['xmomentum'][:] = Q/boundary_width
113
114    fid.createVariable('ymomentum', Float, ('number_of_timesteps',))
115    fid.variables['ymomentum'][:] = 0.0
116
117    fid.close()
118
119
120#-------------------------------------------------------------
121if __name__ == "__main__":
122
123
124    # Prepare time boundary
125    prepare_timeboundary(filename=project.boundary_filename,
126                         hydrograph='HYDROGRAPHS_SUB1',
127                         field='Qinto_OS')
128
129
130    #--------------------------------------------------------------------------
131    # Create the triangular mesh based on overall clipping polygon with a
132    # tagged
133    # boundary and interior regions defined in project.py along with
134    # resolutions (maximal area of per triangle) for each polygon
135    #--------------------------------------------------------------------------
136
137
138    base_resolution = 1 # Use this to coarsen or refine entire mesh.
139
140    # Basic geometry and bounding polygon
141    xleft   = 0
142    xright  = 5.448
143    ybottom = 0
144    ytop    = 3.402
145
146    point_sw = [xleft, ybottom]
147    point_se = [xright, ybottom]
148    point_nw = [xleft, ytop]   
149    point_ne = [xright, ytop]
150
151    bounding_polygon = [point_se,
152                        point_ne,
153                        point_nw,
154                        point_sw]
155   
156
157    # Localised refined area for gulleys
158    xl = 4.8
159    xr = 5.3
160    yb = 1.6
161    yt = 2.3
162    p0 = [xl, yb]
163    p1 = [xr, yb]
164    p2 = [xr, yt]
165    p3 = [xl, yt]
166   
167    gulleys = [p0, p1, p2, p3]
168   
169
170
171    # Island area and drawdown region
172    island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5]
173    island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3]
174    island_2 = [xleft + (xright-xleft)/2+0.4, ybottom + 2*(ytop-ybottom)/3-0.3]
175    island_3 = [xleft + (xright-xleft)/2+0.4, ybottom + (ytop-ybottom)/3+0.3]
176    island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3]
177    island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.8]
178    island_6 = [xl-.01, yb]  # Keep right edge just off the gulleys
179    island_7 = [xl-.01, yt]
180 
181    island = [island_0, island_1, island_2,
182              island_3, island_4, island_5,
183              island_6, island_7]
184
185
186
187    # Region spanning half right hand side of domain just inside boundary
188    rhs_nw = [xleft + (xright-xleft)/3+1, ytop-1.4]
189    rhs_sw = [xleft + (xright-xleft)/3+1, ybottom+0.5]
190    rhs_se = [xright-0.1, ybottom+0.2]
191    rhs_ne = [xright-0.1, ytop-0.2]       
192
193    rhs_region = [rhs_nw, rhs_ne, rhs_se, rhs_sw]
194
195   
196    # Interior regions and creation of mesh
197    interior_regions = [[rhs_region, 0.0005],
198                        [island, 0.0003*base_resolution],
199                        [gulleys, 0.00003*base_resolution]]   
200
201    meshname = project.mesh_filename + '.msh'
202    m = create_mesh_from_regions(bounding_polygon,
203                                 boundary_tags={'wall': [0, 1, 3],
204                                                'wave': [2]},     
205                                 maximum_triangle_area=0.1*base_resolution,
206                                 interior_regions=interior_regions,
207                                 filename=project.mesh_filename,
208                                 use_cache=False,
209                                 verbose=True)
210
Note: See TracBrowser for help on using the repository browser.