source: anuga_validation/okushiri_2005/create_okushiri.py @ 7795

Last change on this file since 7795 was 7307, checked in by rwilson, 15 years ago

Added import that fixes 'netcdf_float unknown' error.

File size: 6.5 KB
Line 
1"""Create mesh and time boundary for the Okushiri island validation
2"""
3
4from anuga.pmesh.mesh import *
5from anuga.pmesh.mesh_interface import create_mesh_from_regions
6from anuga.coordinate_transforms.geo_reference import Geo_reference
7from anuga.geospatial_data import Geospatial_data
8from anuga.config import netcdf_float
9
10from Scientific.IO.NetCDF import NetCDFFile
11
12import project
13
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
46
47def prepare_timeboundary(filename):
48    """Convert benchmark 2 time series to NetCDF tms file.
49    This is a 'throw-away' code taylor made for files like
50    'Benchmark_2_input.txt' from the LWRU2 benchmark
51    """
52
53    print 'Creating', filename
54
55    # Read the ascii (.txt) version of this file
56    fid = open(filename[:-4] + '.txt')
57
58    # Skip first line
59    line = fid.readline()
60
61    # Read remaining lines
62    lines = fid.readlines()
63    fid.close()
64
65
66    N = len(lines)
67    T = num.zeros(N, num.float)  #Time
68    Q = num.zeros(N, num.float)  #Values
69
70    for i, line in enumerate(lines):
71        fields = line.split()
72
73        T[i] = float(fields[0])
74        Q[i] = float(fields[1])
75
76
77    # Create tms NetCDF file
78
79    fid = NetCDFFile(filename, 'w')
80    fid.institution = 'Geoscience Australia'
81    fid.description = 'Input wave for Benchmark 2'
82    fid.starttime = 0.0
83    fid.createDimension('number_of_timesteps', len(T))
84    fid.createVariable('time', netcdf_float, ('number_of_timesteps',))
85    fid.variables['time'][:] = T
86
87    fid.createVariable('stage', netcdf_float, ('number_of_timesteps',))
88    fid.variables['stage'][:] = Q[:]
89
90    fid.createVariable('xmomentum', netcdf_float, ('number_of_timesteps',))
91    fid.variables['xmomentum'][:] = 0.0
92
93    fid.createVariable('ymomentum', netcdf_float, ('number_of_timesteps',))
94    fid.variables['ymomentum'][:] = 0.0
95
96    fid.close()
97
98
99#-------------------------------------------------------------
100if __name__ == "__main__":
101
102
103    # Prepare bathymetry
104    prepare_bathymetry(project.bathymetry_filename)
105
106    # Prepare time boundary
107    prepare_timeboundary(project.boundary_filename)
108
109
110    #--------------------------------------------------------------------------
111    # Create the triangular mesh based on overall clipping polygon with a
112    # tagged
113    # boundary and interior regions defined in project.py along with
114    # resolutions (maximal area of per triangle) for each polygon
115    #--------------------------------------------------------------------------
116
117
118    base_resolution = 1 # Use this to coarsen or refine entire mesh.
119
120    # Basic geometry and bounding polygon
121    xleft   = 0
122    xright  = 5.448
123    ybottom = 0
124    ytop    = 3.402
125
126    point_sw = [xleft, ybottom]
127    point_se = [xright, ybottom]
128    point_nw = [xleft, ytop]   
129    point_ne = [xright, ytop]
130
131    bounding_polygon = [point_se,
132                        point_ne,
133                        point_nw,
134                        point_sw]
135   
136
137    # Localised refined area for gulleys
138    xl = 4.8
139    xr = 5.3
140    yb = 1.6
141    yt = 2.3
142    p0 = [xl, yb]
143    p1 = [xr, yb]
144    p2 = [xr, yt]
145    p3 = [xl, yt]
146   
147    gulleys = [p0, p1, p2, p3]
148   
149
150    # Island area and drawdown region (original)
151    #island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5]
152    #island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3]
153    #island_2 = [xleft + (xright-xleft)/2+0.3, ybottom + 2*(ytop-ybottom)/3-0.3]
154    #island_3 = [xleft + (xright-xleft)/2+0.3, ybottom + (ytop-ybottom)/3+0.3]
155    #island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3]
156    #island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.2]
157    #island_6 = [xl-.01, yb]  #OK
158    #island_7 = [xl-.01, yt]  #OK     
159
160
161    # Island area and drawdown region
162    island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5]
163    island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3]
164    island_2 = [xleft + (xright-xleft)/2+0.4, ybottom + 2*(ytop-ybottom)/3-0.3]
165    island_3 = [xleft + (xright-xleft)/2+0.4, ybottom + (ytop-ybottom)/3+0.3]
166    island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3]
167    island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.8]
168    island_6 = [xl-.01, yb]  # Keep right edge just off the gulleys
169    island_7 = [xl-.01, yt]
170 
171    island = [island_0, island_1, island_2,
172              island_3, island_4, island_5,
173              island_6, island_7]
174
175
176    # Region spanning half right hand side of domain just inside boundary (org)
177    #rhs_nw = [xleft + (xright-xleft)/3+1, ytop-0.02]
178    #rhs_sw = [xleft + (xright-xleft)/3+1, ybottom+0.02]
179    #rhs_se = [xright-0.02, ybottom+0.02]
180    #rhs_ne = [xright-0.02, ytop-0.02]
181
182    # Region spanning half right hand side of domain just inside boundary
183    rhs_nw = [xleft + (xright-xleft)/3+1, ytop-1.4]
184    rhs_sw = [xleft + (xright-xleft)/3+1, ybottom+0.5]
185    rhs_se = [xright-0.1, ybottom+0.2]
186    rhs_ne = [xright-0.1, ytop-0.2]       
187
188    rhs_region = [rhs_nw, rhs_ne, rhs_se, rhs_sw]
189
190   
191    # Interior regions and creation of mesh
192    interior_regions = [[rhs_region, 0.0005],
193                        [island, 0.0002*base_resolution],
194                        [gulleys, 0.00002*base_resolution]]   
195
196    meshname = project.mesh_filename + '.msh'
197    m = create_mesh_from_regions(bounding_polygon,
198                                 boundary_tags={'wall': [0, 1, 3],
199                                                'wave': [2]},     
200                                 maximum_triangle_area=0.1*base_resolution,
201                                 interior_regions=interior_regions,
202                                 filename=project.mesh_filename,
203                                 verbose=True)
204
Note: See TracBrowser for help on using the repository browser.