source: anuga_core/source/anuga/examples/pressure_force.py @ 3565

Last change on this file since 3565 was 3565, checked in by ole, 18 years ago

Moved obsolete code and examples to their appropriate locations

File size: 2.3 KB
Line 
1"""Example of shallow water wave equation.
2
3This example applies a dummy forcing function mimicking a
4static cyclone
5"""
6
7######################
8# Module imports
9#
10from shallow_water import Transmissive_boundary, Reflective_boundary,\
11     Dirichlet_boundary, Time_boundary
12from shallow_water import Domain, Constant_height
13from mesh_factory import rectangular
14
15from math import sin, pi
16
17
18#######################
19# Domain
20#
21N = 50
22elements, triangles, boundary = rectangular(N, N)
23domain = Domain(elements, triangles, boundary)
24
25domain.store = True
26domain.filename = 'pressure'
27
28#######################
29#Bed-slope and friction
30def x_slope(x, y):
31    #return -x/3
32    return 0.0*x
33
34domain.set_quantity('elevation', x_slope)
35domain.set_quantity('friction', 0.1)
36
37
38#######################
39#Forcing terms
40def pressure(x,y,t):
41    r=0.1
42    x0 = (1+sin(t*pi/5))/2
43    y0 = (1+sin(t*pi/5))/2
44    if x >= x0-r and x <= x0+r and y >= y0-r and y <= y0+r:
45        return 1000*(1+sin(t*pi))
46    else:
47        return 800
48
49def cyclone(domain):
50    from anuga.config import rho_w
51    from anuga.pyvolution.util import gradient
52
53    xmom = domain.quantities['xmomentum'].explicit_update
54    ymom = domain.quantities['ymomentum'].explicit_update
55
56    Stage = domain.quantities['stage']
57    Elevation = domain.quantities['elevation']   
58    h = Stage.vertex_values - Elevation.vertex_values
59    x = domain.get_vertex_coordinates()
60
61    for k in range(len(domain)):
62        avg_h = sum( h[k,:] )/3
63       
64        #Compute bed slope
65        x0, y0, x1, y1, x2, y2 = x[k,:]   
66
67        p0 = pressure(x0,y0,t)
68        p1 = pressure(x1,y1,t)
69        p2 = pressure(x2,y2,t)       
70       
71        px, py = gradient(x0, y0, x1, y1, x2, y2, p0, p1, p2)
72
73        #Update momentum
74        xmom[k] += -px*avg_h/rho_w
75        ymom[k] += -py*avg_h/rho_w       
76
77
78domain.forcing_terms.append(cyclone)
79
80######################
81# Boundary conditions
82Br = Reflective_boundary(domain)
83
84domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
85
86
87######################
88#Initial condition
89h = 0.05
90domain.set_quantity('stage', Constant_height(x_slope, h))
91
92
93######################
94#Evolution
95for t in domain.evolve(yieldstep = 0.05, finaltime = 7):
96    domain.write_time()
97
98
99   
Note: See TracBrowser for help on using the repository browser.