source: inundation/pyvolution/cornell_room.py @ 1869

Last change on this file since 1869 was 773, checked in by ole, 19 years ago

Changed quantity name 'level' to 'stage'

File size: 3.0 KB
Line 
1"""Example of shallow water wave equation.
2
3This is called Netherlands because it shows a dam with a gap in it and
4stylised housed behind it and below the water surface.
5
6"""
7
8######################
9# Module imports
10#
11from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
12     Transmissive_boundary, Time_boundary, Constant_height
13
14from mesh_factory import from_polyfile, rectangular
15from Numeric import array
16from math import sqrt
17from least_squares import Interpolation
18   
19
20print 'Creating domain'
21#data_points, _, data_values = from_polyfile('cornell_room_medres')
22#points, triangles, values = from_polyfile('hires2')
23data_points, _, data_values = from_polyfile('hires2')
24
25
26#Regrid onto numerically stable mesh
27#
28#Compute regular mesh based on resolution and extent of data
29data_points = array(data_points)
30pmax = max(data_points)
31pmin = min(data_points)
32
33M = len(data_points)
34
35N = int(0.8*sqrt(M))
36
37#print N
38
39mesh_points, vertices, boundary = rectangular(N, N, 
40                                              len1=pmax[0]-pmin[0],
41                                              len2=pmax[1]-pmin[1],
42                                              origin = pmin)
43                                             
44
45#Compute smooth surface on new mesh based on values from old (regrid)
46print 'Interp'
47interp = Interpolation(mesh_points, vertices, data_points, alpha=10)
48mesh_values = interp.fit(data_values)
49print 'Len mesh values', len(mesh_values)
50print 'Len mesh points', len(mesh_points)
51
52   
53#Create shallow water domain
54print 'Creating domain'
55domain = Domain(mesh_points, vertices) #, boundary)
56
57domain.check_integrity()
58domain.default_order = 2
59domain.smooth = True
60domain.reduction = min  #Looks a lot better on top of steep slopes
61
62print "Number of triangles = ", len(domain)
63
64domain.visualise = False
65domain.checkpoint = False
66domain.store = True    #Store for visualisation purposes
67domain.format = 'sww'   #Native netcdf visualisation format
68import sys, os
69root, ext = os.path.splitext(sys.argv[0])
70if domain.smooth is True:
71    s = 'smooth'
72else:
73    s = 'nonsmooth'       
74domain.filename = root + '_' + s
75
76#Set bed-slope and friction
77manning = 0.0
78
79print 'Field values'
80domain.set_quantity('elevation', mesh_values)
81domain.set_quantity('friction', manning)
82
83
84######################
85# Boundary conditions
86#
87print 'Boundaries'
88Br = Reflective_boundary(domain)
89domain.set_boundary({'exterior': Br})
90
91                 
92
93######################
94#Initial condition
95#
96print 'Initial condition'
97
98#Define water height as a lump in one corner
99def height(x, y):
100    from Numeric import zeros, Float
101   
102    N = len(x)
103    assert N == len(y)   
104
105    xmin = min(x); xmax = max(x)
106    ymin = min(y); ymax = max(y)
107
108    xrange = xmax - xmin
109    yrange = ymax - ymin   
110
111    z = zeros(N, Float)
112    for i in range(N):
113        if x[i] <= xmin + 0.25*xrange and y[i] <= ymin + 0.25*yrange:
114            z[i] = 300
115
116    return z
117
118domain.set_quantity('stage', height)
119
120E = domain.quantities['elevation'].vertex_values
121L = domain.quantities['stage'].vertex_values
122domain.set_quantity('stage', E+L)
123
124#Evolve
125for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
126    domain.write_time()
127   
128
129   
Note: See TracBrowser for help on using the repository browser.