source: inundation/pyvolution/twolevels.py @ 1740

Last change on this file since 1740 was 1014, checked in by ole, 19 years ago
File size: 1.8 KB
Line 
1"""Example of shallow water wave equation.
2
3Two levels - testing that momentum is not genearted in the upward direction
4
5"""
6
7######################
8# Module imports
9#
10from os import sep, path
11from mesh_factory import rectangular
12from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
13     Constant_height
14from Numeric import array
15from util import Polygon_function, read_polygon
16
17
18#Create basic mesh
19N = 2
20points, vertices, boundary = rectangular(3, N, 20, 10)
21
22#Create shallow water domain
23domain = Domain(points, vertices, boundary)
24domain.store = True
25domain.set_name('two_levels')
26print "Output being written to " + domain.get_datadir() + sep + \
27              domain.filename + "_size%d." %len(domain) + domain.format
28       
29
30domain.default_order=2
31domain.smooth = False
32domain.visualise = True
33
34#PLAY WITH THIS [0;1]:
35#
36# beta_h == 0.0 reveals the problem
37# beta_h > 0.2 alleviates it
38domain.beta_h = 0.2   
39
40#IC
41domain.set_quantity('friction', 0.0)
42domain.set_quantity('stage', 3.0)
43
44def twostage(x, y): 
45    z = 0*x     
46    for i in range(len(x)):             
47        if x[i]<10:     
48            z[i] = 4.
49
50    return z                   
51#Set elevation 
52domain.set_quantity('elevation', twostage)
53
54
55######################
56# Boundary conditions
57Br = Reflective_boundary(domain)
58
59domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
60domain.check_integrity()
61
62
63#print domain.quantities['elevation'].vertex_values   
64######################
65#Evolution
66for t in domain.evolve(yieldstep = 0.01, finaltime = 0.05):
67    domain.write_time()
68    print domain.quantities['xmomentum'].edge_values   
69    #print domain.quantities['stage'].centroid_values   
70
71#print domain.quantities['stage'].vertex_values
72#print
73#print domain.quantities['xmomentum'].vertex_values
74#print
75#print domain.quantities['ymomentum'].vertex_values
76   
77print 'Done'   
78
Note: See TracBrowser for help on using the repository browser.