source: anuga_validation/automated_validation_tests/okushiri_tank_validation/create_okushiri.py @ 3825

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

Work on validation

File size: 4.4 KB
Line 
1"""Create mesh and time boundary for the Okushiri island validation
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
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    print 'Preparing time boundary from %s' %filename
21    from Numeric import array
22
23    fid = open(filename)
24
25    #Skip first line
26    line = fid.readline()
27
28    #Read remaining lines
29    lines = fid.readlines()
30    fid.close()
31
32
33    N = len(lines)
34    T = zeros(N, Float)  #Time
35    Q = zeros(N, Float)  #Values
36
37    for i, line in enumerate(lines):
38        fields = line.split()
39
40        T[i] = float(fields[0])
41        Q[i] = float(fields[1])
42
43
44    #Create tms file
45    from Scientific.IO.NetCDF import NetCDFFile
46
47    outfile = filename[:-4] + '.tms'
48    print 'Writing to', outfile
49    fid = NetCDFFile(outfile, 'w')
50
51    fid.institution = 'Geoscience Australia'
52    fid.description = 'Input wave for Benchmark 2'
53    fid.starttime = 0.0
54    fid.createDimension('number_of_timesteps', len(T))
55    fid.createVariable('time', Float, ('number_of_timesteps',))
56    fid.variables['time'][:] = T
57
58    fid.createVariable('stage', Float, ('number_of_timesteps',))
59    fid.variables['stage'][:] = Q[:]
60
61    fid.createVariable('xmomentum', Float, ('number_of_timesteps',))
62    fid.variables['xmomentum'][:] = 0.0
63
64    fid.createVariable('ymomentum', Float, ('number_of_timesteps',))
65    fid.variables['ymomentum'][:] = 0.0
66
67    fid.close()
68
69
70#-------------------------------------------------------------
71if __name__ == "__main__":
72
73    # Prepare time boundary
74    prepare_timeboundary(project.boundary_filename)
75
76
77    #--------------------------------------------------------------------------
78    # Create the triangular mesh based on overall clipping polygon with a
79    # tagged
80    # boundary and interior regions defined in project.py along with
81    # resolutions (maximal area of per triangle) for each polygon
82    #--------------------------------------------------------------------------
83
84    base_resolution = 10
85
86    # Basic geometry
87    xleft   = 0
88    xright  = 5.448
89    ybottom = 0
90    ytop    = 3.402
91
92    # Outline
93    point_sw = [xleft, ybottom]
94    point_se = [xright, ybottom]
95    point_nw = [xleft, ytop]   
96    point_ne = [xright, ytop]
97
98    # Midway points (left)
99    point_mtop = [xleft + (xright-xleft)/3+1, ytop]
100    point_mbottom = [xleft + (xright-xleft)/3+1, ybottom]
101
102    #Localised refined area for gulleys
103    xl = 4.8
104    xr = 5.3
105    yb = 1.6
106    yt = 2.3
107    p0 = [xl, yb]
108    p1 = [xr, yb]
109    p2 = [xr, yt]
110    p3 = [xl, yt]
111    gulleys = [p0, p1, p2, p3]
112
113    #Island area
114    island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5]
115    island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3]
116    island_2 = [xleft + (xright-xleft)/2+0.3, ybottom + 2*(ytop-ybottom)/3-0.3]
117    island_3 = [xleft + (xright-xleft)/2+0.3, ybottom + (ytop-ybottom)/3+0.3]
118    island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3]
119    island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.2]
120    island_6 = [xl-.01, yb]  #OK
121    island_7 = [xl-.01, yt]  #OK     
122
123    island = [island_0, island_1, island_2,
124              island_3, island_4, island_5,
125              island_6, island_7]
126
127
128   
129
130    bounding_polygon = [point_se,
131                        point_ne,
132                        point_mtop,
133                        point_nw,
134                        point_sw,
135                        point_mbottom]
136   
137
138    interior_regions = [[gulleys, 0.00002*base_resolution],
139                        [island, 0.00007*base_resolution]]
140   
141
142
143   
144
145    print 'start create mesh from regions'
146
147    meshname = project.mesh_filename + '.msh'
148    m = create_mesh_from_regions(bounding_polygon,
149                                 boundary_tags={'wall': [0, 1, 2, 4, 5],
150                                                'wave': [3]},
151                                 maximum_triangle_area=0.1*base_resolution,
152                                 interior_regions=interior_regions,
153                                 filename=project.mesh_filename,
154                                 use_cache=True,
155                                 verbose=True)
156
Note: See TracBrowser for help on using the repository browser.