source: anuga_core/source/anuga/examples/wind_example_variable.py @ 4026

Last change on this file since 4026 was 4026, checked in by jack, 17 years ago

Moved the vpython visualiser to obsolete_code and cleared out things that call it from other code.

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