source: anuga_validation/convergence_study/convergence_structured.py @ 5306

Last change on this file since 5306 was 5162, checked in by steve, 16 years ago

Updated some methods for quantity. Looks like we can use old
limiting system with larger values of beta.

File size: 4.3 KB
Line 
1"""Simple water flow example using ANUGA
2
3Water driven up a linear slope and time varying boundary,
4similar to a beach environment
5"""
6
7#------------------------------------------------------------------------------
8# Import necessary modules
9#------------------------------------------------------------------------------
10
11import sys
12from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
13from anuga.shallow_water import Domain
14from anuga.shallow_water import Reflective_boundary
15from anuga.shallow_water import Dirichlet_boundary
16from anuga.shallow_water import Time_boundary
17from anuga.shallow_water import Transmissive_boundary
18from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
19#from anuga.geospatial_data.geospatial_data import *
20from math import cos
21from Numeric import zeros, Float
22
23#------------------------------------------------------------------------------
24# Setup computational domain
25#------------------------------------------------------------------------------
26dx = 200.
27dy = dx
28L = 100000.
29W = dx
30
31# structured mesh
32points, vertices, boundary = rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, -W/2))
33
34domain = Domain(points, vertices, boundary) 
35
36domain.set_timestepping_method('rk2')
37domain.set_default_order(2)
38domain.set_name('myexample9')               
39domain.set_datadir('.')                     # Use current directory for output
40
41domain.set_all_limiters(0.9)
42
43print domain.beta_w
44domain.use_old_limiter = False
45domain.CFL = 1.0
46
47#------------------------------------------------------------------------------
48# Setup initial conditions
49#------------------------------------------------------------------------------
50#domain.set_quantity('elevation', topography) # Use function for elevation
51domain.set_quantity('elevation',0.0)
52domain.set_quantity('friction', 0.00)
53
54h0 = 10.0
55h1 = 1.0
56
57def height(x,y):
58    z = zeros(len(x),Float)
59    for i in range(len(x)):
60        if x[i]<=50000.0:
61            z[i] = h0
62        else:
63            z[i] = h1
64    return z
65
66
67domain.set_quantity('stage', height)
68#domain.set_quantity('stage', 0.0)           
69
70#-----------------------------------------------------------------------------
71# Setup boundary conditions
72#------------------------------------------------------------------------------
73from math import sin, pi, exp
74Br = Reflective_boundary(domain)      # Solid reflective wall
75Bt = Transmissive_boundary(domain)    # Continue all values on boundary
76Bd = Dirichlet_boundary([1,0.,0.]) # Constant boundary values
77amplitude = 1
78#Bw = Transmissive_Momentum_Set_Stage_boundary(domain=domain,
79Bw = Time_boundary(domain=domain,     # Time dependent boundary 
80## Sine wave
81#                   f=lambda t: [(-amplitude*sin((1./300.)*t*2*pi)), 0.0, 0.0])
82## Single wave
83                   f=lambda t: [h0, 0.0, 0.0])
84## Sawtooth?
85#                   f=lambda t: [(-8.0*(sin((1./180.)*t*2*pi))+(1./2.)*sin((2./180.)*t*2*pi)+(1./3.)*sin((3./180.)*t*2*pi)), 0.0, 0.0])
86## Sharp rise, linear fall
87#                   f=lambda t: [(5.0*(-((t-0.)/300.)*(t<300.)-cos((t-300.)*2.*pi*(1./240.))*(t>=300. and t<420.)+(1.-(t-420.)/300.)*(t>=420. and t <720.))), 0.0, 0.0])
88#                   f=lambda t: [amplitude*(1.-2.*(pi*(1./720.)*(t-720.))**2)/exp((pi*(1./720.)*(t-720.))**2) , 0.0, 0.0])
89#                   f=lambda t: [(-8.0*sin((1./720.)*t*2*pi))*((t<720.)-0.5*(t<360.)), 0.0, 0.0])
90
91# Associate boundary tags with boundary objects
92domain.set_boundary({'left': Bw, 'right': Bt, 'top': Br, 'bottom': Br})
93
94
95#===============================================================================
96from anuga.visualiser import RealtimeVisualiser
97vis = RealtimeVisualiser(domain)
98#vis.render_quantity_height("elevation", zScale=1, offset = 5.0, dynamic=False)
99vis.render_quantity_height("stage", zScale =10000, dynamic=True)
100#vis.colour_height_quantity('stage', (lambda q:q['stage'], -1.0, 1.0))
101vis.colour_height_quantity('stage', (1.0, 0.5, 0.5))
102vis.start()
103#===============================================================================
104
105
106#------------------------------------------------------------------------------
107# Evolve system through time
108#------------------------------------------------------------------------------
109
110for t in domain.evolve(yieldstep = 50.0, finaltime = 10*40*60.):
111    domain.write_time()
112    vis.update()
113   
114vis.evolveFinished()
115
Note: See TracBrowser for help on using the repository browser.