source: anuga_validation/okushiri_2005/create_okushiri.py @ 3647

Last change on this file since 3647 was 3647, checked in by ole, 18 years ago

Rearranging okushiri

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