source: inundation/pyvolution/wind_example_variable.py @ 1884

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

Changed quantity name 'level' to 'stage'

File size: 2.2 KB
Line 
1"""Example of shallow water wave equation.
2
3Flat bed with rotational wind stress
4
5"""
6
7######################
8# Module imports
9#
10from mesh_factory import rectangular
11from shallow_water import Domain, Dirichlet_boundary, Wind_stress
12
13#Create basic mesh
14N = 20
15length = 200
16points, vertices, boundary = rectangular(N, N, length, length)
17
18#Create shallow water domain
19domain = Domain(points, vertices, boundary)
20domain.smooth = True
21domain.visualise = False
22domain.store = True
23domain.default_order=2
24domain.set_name('wind_rotation')
25
26#Set initial conditions
27h = 1.0
28domain.set_quantity('elevation', 0.0)
29domain.set_quantity('stage', h)
30domain.set_quantity('friction', 0.01)
31
32
33#Variable windfield implemented using functions
34def speed(t,x,y):
35    """Large speeds halfway between center and edges
36    Low speeds at center and edges
37    """
38
39    from math import pi
40    from Numeric import sqrt, exp, cos
41
42    c = (length/2, length/2)
43    r = sqrt((x - c[0])**2 + (y - c[1])**2)/length
44    factor = exp( -(r-0.15)**2 )
45
46    #return 9000 * factor
47    return 4000 * factor * (cos(t*2*pi/150) + 2)
48
49
50def phi(t,x,y):
51    """Rotating field
52    """
53   
54    from math import pi
55    from Numeric import sqrt, exp, cos, arctan2, choose, less
56
57    c = (length/2, length/2)
58    xx = (x - c[0])/length
59    yy = (y - c[1])/length
60    angle = arctan2(yy,xx)
61
62    #Take normal direction (but reverse direction every 50 seconds)
63    #if sin(t*2*pi/100) < 0:
64    #    sgn = -1
65    #else:
66    #    sgn = 1
67    #angle += sgn*pi/2
68    angle -= pi/2   
69
70    #Convert to degrees and return
71    return angle/pi*180
72
73   
74domain.forcing_terms.append( Wind_stress(speed, phi) )
75
76
77#Add lateral wind gusts bearing 25 deg
78def gust(t,x,y): 
79    from math import sin, pi
80    from Numeric import zeros, ones, Float
81
82    N = len(x)
83
84    tt = sin(2*pi*t/200)
85
86    if tt > 0.9:
87        return 6000*tt*ones(N, Float)
88    else:
89        return zeros(N, Float)
90   
91domain.forcing_terms.append(Wind_stress(gust, 25))
92
93
94######################
95# Boundary conditions
96#Br = Reflective_boundary(domain)
97Bd = Dirichlet_boundary([h, 0, 0])
98domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
99
100######################
101#Evolution
102for t in domain.evolve(yieldstep = 0.5, finaltime = 1000):
103    domain.write_time()
104
105print 'Done'   
106
Note: See TracBrowser for help on using the repository browser.