[3644] | 1 | """Create mesh and time boundary for the Okushiri island validation |
---|
[2229] | 2 | """ |
---|
| 3 | |
---|
[3535] | 4 | from anuga.pmesh.mesh import * |
---|
[3854] | 5 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
[3535] | 6 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
[3883] | 7 | from anuga.geospatial_data import Geospatial_data |
---|
[7307] | 8 | from anuga.config import netcdf_float |
---|
[2229] | 9 | |
---|
[7276] | 10 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 11 | |
---|
[3647] | 12 | import project |
---|
| 13 | |
---|
[3883] | 14 | def 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 | """ |
---|
[3647] | 19 | |
---|
[3883] | 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 | |
---|
[3647] | 47 | def prepare_timeboundary(filename): |
---|
[3883] | 48 | """Convert benchmark 2 time series to NetCDF tms file. |
---|
[3647] | 49 | This is a 'throw-away' code taylor made for files like |
---|
| 50 | 'Benchmark_2_input.txt' from the LWRU2 benchmark |
---|
| 51 | """ |
---|
| 52 | |
---|
[3883] | 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 |
---|
[3647] | 59 | line = fid.readline() |
---|
| 60 | |
---|
[3883] | 61 | # Read remaining lines |
---|
[3647] | 62 | lines = fid.readlines() |
---|
| 63 | fid.close() |
---|
| 64 | |
---|
| 65 | |
---|
| 66 | N = len(lines) |
---|
[7276] | 67 | T = num.zeros(N, num.float) #Time |
---|
| 68 | Q = num.zeros(N, num.float) #Values |
---|
[3647] | 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 | |
---|
[3883] | 77 | # Create tms NetCDF file |
---|
[3647] | 78 | |
---|
[3850] | 79 | fid = NetCDFFile(filename, 'w') |
---|
[3647] | 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)) |
---|
[7276] | 84 | fid.createVariable('time', netcdf_float, ('number_of_timesteps',)) |
---|
[3647] | 85 | fid.variables['time'][:] = T |
---|
| 86 | |
---|
[7276] | 87 | fid.createVariable('stage', netcdf_float, ('number_of_timesteps',)) |
---|
[3647] | 88 | fid.variables['stage'][:] = Q[:] |
---|
| 89 | |
---|
[7276] | 90 | fid.createVariable('xmomentum', netcdf_float, ('number_of_timesteps',)) |
---|
[3647] | 91 | fid.variables['xmomentum'][:] = 0.0 |
---|
| 92 | |
---|
[7276] | 93 | fid.createVariable('ymomentum', netcdf_float, ('number_of_timesteps',)) |
---|
[3647] | 94 | fid.variables['ymomentum'][:] = 0.0 |
---|
| 95 | |
---|
| 96 | fid.close() |
---|
| 97 | |
---|
| 98 | |
---|
[2229] | 99 | #------------------------------------------------------------- |
---|
| 100 | if __name__ == "__main__": |
---|
| 101 | |
---|
[3850] | 102 | |
---|
[3883] | 103 | # Prepare bathymetry |
---|
| 104 | prepare_bathymetry(project.bathymetry_filename) |
---|
[2229] | 105 | |
---|
[3647] | 106 | # Prepare time boundary |
---|
| 107 | prepare_timeboundary(project.boundary_filename) |
---|
| 108 | |
---|
| 109 | |
---|
[3854] | 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 | #-------------------------------------------------------------------------- |
---|
[3647] | 116 | |
---|
[3861] | 117 | |
---|
[3883] | 118 | base_resolution = 1 # Use this to coarsen or refine entire mesh. |
---|
[3854] | 119 | |
---|
[3883] | 120 | # Basic geometry and bounding polygon |
---|
[2229] | 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 | |
---|
[3883] | 131 | bounding_polygon = [point_se, |
---|
| 132 | point_ne, |
---|
| 133 | point_nw, |
---|
| 134 | point_sw] |
---|
| 135 | |
---|
[2229] | 136 | |
---|
[3883] | 137 | # Localised refined area for gulleys |
---|
[2229] | 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] |
---|
[3883] | 146 | |
---|
[3854] | 147 | gulleys = [p0, p1, p2, p3] |
---|
[3883] | 148 | |
---|
[2229] | 149 | |
---|
[3913] | 150 | # Island area and drawdown region (original) |
---|
[3915] | 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 |
---|
[3913] | 159 | |
---|
| 160 | |
---|
| 161 | # Island area and drawdown region |
---|
[3915] | 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] |
---|
[3883] | 170 | |
---|
[3854] | 171 | island = [island_0, island_1, island_2, |
---|
| 172 | island_3, island_4, island_5, |
---|
| 173 | island_6, island_7] |
---|
[2229] | 174 | |
---|
| 175 | |
---|
[3913] | 176 | # Region spanning half right hand side of domain just inside boundary (org) |
---|
[3915] | 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] |
---|
[3913] | 181 | |
---|
[3883] | 182 | # Region spanning half right hand side of domain just inside boundary |
---|
[3915] | 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] |
---|
[2229] | 187 | |
---|
[3883] | 188 | rhs_region = [rhs_nw, rhs_ne, rhs_se, rhs_sw] |
---|
[3878] | 189 | |
---|
[3854] | 190 | |
---|
[3883] | 191 | # Interior regions and creation of mesh |
---|
| 192 | interior_regions = [[rhs_region, 0.0005], |
---|
[3915] | 193 | [island, 0.0002*base_resolution], |
---|
| 194 | [gulleys, 0.00002*base_resolution]] |
---|
[2229] | 195 | |
---|
[3854] | 196 | meshname = project.mesh_filename + '.msh' |
---|
| 197 | m = create_mesh_from_regions(bounding_polygon, |
---|
[3883] | 198 | boundary_tags={'wall': [0, 1, 3], |
---|
| 199 | 'wave': [2]}, |
---|
[3854] | 200 | maximum_triangle_area=0.1*base_resolution, |
---|
| 201 | interior_regions=interior_regions, |
---|
| 202 | filename=project.mesh_filename, |
---|
| 203 | verbose=True) |
---|
[2229] | 204 | |
---|