source: anuga_validation/okushiri_2005/create_okushiri_original.py @ 7795

Last change on this file since 7795 was 7276, checked in by ole, 15 years ago

Merged numpy branch back into the trunk.

In ~/sandpit/anuga/anuga_core/source
svn merge -r 6246:HEAD ../../branches/numpy .

In ~/sandpit/anuga/anuga_validation
svn merge -r 6417:HEAD ../branches/numpy_anuga_validation .

In ~/sandpit/anuga/misc
svn merge -r 6809:HEAD ../branches/numpy_misc .

For all merges, I used numpy version where conflicts existed

The suites test_all.py (in source/anuga) and validate_all.py passed using Python2.5 with numpy on my Ubuntu Linux box.

File size: 5.9 KB
Line 
1"""Create mesh and time boundary for the Okushiri island validation
2"""
3
4import numpy as num
5from Scientific.IO.NetCDF import NetCDFFile
6
7from anuga.pmesh.mesh import *
8from anuga.coordinate_transforms.geo_reference import Geo_reference
9from anuga.config import netcdf_float
10
11import project
12
13
14def prepare_timeboundary(filename):
15    """Converting benchmark 2 time series to NetCDF tms file.
16    This is a 'throw-away' code taylor made for files like
17    'Benchmark_2_input.txt' from the LWRU2 benchmark
18    """
19
20    textversion = filename[:-4] + '.txt'
21
22    print 'Preparing time boundary from %s' %textversion
23
24    fid = open(textversion)
25
26    #Skip first line
27    line = fid.readline()
28
29    #Read remaining lines
30    lines = fid.readlines()
31    fid.close()
32
33
34    N = len(lines)
35    T = num.zeros(N, num.float)  #Time
36    Q = num.zeros(N, num.float)  #Values
37
38    for i, line in enumerate(lines):
39        fields = line.split()
40
41        T[i] = float(fields[0])
42        Q[i] = float(fields[1])
43
44
45    #Create tms file
46    print 'Writing to', filename
47    fid = NetCDFFile(filename, 'w')
48
49    fid.institution = 'Geoscience Australia'
50    fid.description = 'Input wave for Benchmark 2'
51    fid.starttime = 0.0
52    fid.createDimension('number_of_timesteps', len(T))
53    fid.createVariable('time', netcdf_float, ('number_of_timesteps',))
54    fid.variables['time'][:] = T
55
56    fid.createVariable('stage', netcdf_float, ('number_of_timesteps',))
57    fid.variables['stage'][:] = Q[:]
58
59    fid.createVariable('xmomentum', netcdf_float, ('number_of_timesteps',))
60    fid.variables['xmomentum'][:] = 0.0
61
62    fid.createVariable('ymomentum', netcdf_float, ('number_of_timesteps',))
63    fid.variables['ymomentum'][:] = 0.0
64
65    fid.close()
66
67
68#-------------------------------------------------------------
69if __name__ == "__main__":
70
71
72    #Preparing points
73    #(Only when using original Benchmark_2_Bathymetry.txt file)
74    # FIXME: Replace by using geospatial data
75    #
76    #from anuga.abstract_2d_finite_volumes.data_manager import xya2pts
77    #xya2pts(project.bathymetry_filename, verbose = True,
78    #        z_func = lambda z: -z)
79
80
81    # Prepare time boundary
82    prepare_timeboundary(project.boundary_filename)
83
84
85    # Create Mesh
86
87    #Basic geometry
88   
89    xleft   = 0
90    xright  = 5.448
91    ybottom = 0
92    ytop    = 3.402
93
94    #Outline
95    point_sw = [xleft, ybottom]
96    point_se = [xright, ybottom]
97    point_nw = [xleft, ytop]   
98    point_ne = [xright, ytop]
99
100    #Midway points (left)
101    point_mtop = [xleft + (xright-xleft)/3+1, ytop]
102    point_mbottom = [xleft + (xright-xleft)/3+1, ybottom]
103
104
105    geo = Geo_reference(xllcorner = xleft,
106                        yllcorner = ybottom)
107   
108                       
109    print "***********************"
110    print "geo ref", geo
111    print "***********************"
112   
113    m = Mesh(geo_reference=geo)
114
115    #Boundary
116    dict = {}
117    dict['points'] = [point_se,   #se
118                      point_ne,
119                      point_mtop,
120                      point_nw,
121                      point_sw,
122                      point_mbottom]
123
124   
125    dict['segments'] = [[0,1], [1,2], [2,3], [3,4],
126                        [4,5], [5,0],  #The outer border
127                        [2,5]]         #Separator
128   
129    dict['segment_tags'] = ['wall',
130                            'wall',
131                            'wall',
132                            'wave',
133                            'wall',
134                            'wall',
135                            '']           #Interior
136
137       
138    m.addVertsSegs(dict)
139
140
141
142
143
144    #Localised refined area for gulleys
145    dict = {}
146    xl = 4.8
147    xr = 5.3
148    yb = 1.6
149    yt = 2.3
150    p0 = [xl, yb]
151    p1 = [xr, yb]
152    p2 = [xr, yt]
153    p3 = [xl, yt]
154   
155    dict['points'] = [p0, p1, p2, p3]
156    dict['segments'] = [[0,1], [1,2], [2,3], [3,0]]
157    dict['segment_tags'] = ['', '', '', '']
158    m.addVertsSegs(dict)   
159
160    #Island area
161    island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5]
162    island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3]
163    island_2 = [xleft + (xright-xleft)/2+0.3, ybottom + 2*(ytop-ybottom)/3-0.3]
164    island_3 = [xleft + (xright-xleft)/2+0.3, ybottom + (ytop-ybottom)/3+0.3]
165    island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3] 
166    island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.2]
167    island_6 = [xl-.01, yb]  #OK
168    island_7 = [xl-.01, yt]  #OK     
169   
170
171    dict['points'] = [island_0, island_1, island_2,
172                      island_3, island_4, island_5,
173                      #p0, p3]                     
174                      island_6, island_7]
175
176                     
177    dict['segments'] = [[0,1], [1,2], [2,3], [3,4],
178                        [4,5], [5,6], [6,7], [7,0]]
179                       
180    dict['segment_tags'] = ['', '', '', '', '', '', '', '']
181
182
183    m.addVertsSegs(dict)   
184
185   
186#
187   
188    base_resolution = 1
189
190    ocean = m.addRegionEN(xleft+.1, ybottom+.1)
191    ocean.setMaxArea(0.1*base_resolution)
192
193    mid = m.addRegionEN(point_mbottom[0]+.1, ybottom+.1)
194    mid.setMaxArea(0.0005*base_resolution)
195 
196
197    #inner = m.addRegionEN(island_3[0]+.1, island_3[1]+.1)
198    #inner.setMaxArea(0.00007*base_resolution)
199
200    inner = m.addRegionEN(island_3[0]+.1, island_3[1]+.1)
201    inner.setMaxArea(0.0003*base_resolution)   
202
203
204    gulleys = m.addRegionEN(p0[0]+0.1, p0[1]+0.1)
205    gulleys.setMaxArea(0.00003*base_resolution)
206
207
208    # From r 1709 11 August 2005
209    #base_resolution = 100
210    #
211    #ocean = m.addRegionEN(xleft+.1, ybottom+.1)
212    #ocean.setMaxArea(0.01*base_resolution)
213    #
214    #mid = m.addRegionEN(point_mbottom[0]+.1, ybottom+.1)
215    #mid.setMaxArea(0.001*base_resolution)
216
217    #inner = m.addRegionEN(island_3[0]+.1, island_3[1]+.1)
218    #inner.setMaxArea(0.0001*base_resolution)
219
220
221    #gulleys = m.addRegionEN(p0[0]+0.1, p0[1]+0.1)
222    #gulleys.setMaxArea(0.00001*base_resolution)       
223   
224   
225    m.generateMesh('pzq28.0za1000000a')
226
227    m.export_mesh_file(project.mesh_filename)
228   
Note: See TracBrowser for help on using the repository browser.