source: branches/numpy_anuga_validation/convergence_study/dam_break.py @ 7952

Last change on this file since 7952 was 5442, checked in by ole, 17 years ago

Retired h-limiter and beta_h as per ticket:194.
All unit tests and validation tests pass.

File size: 3.9 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
19from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files
20
21from math import cos
22from Numeric import zeros, Float
23from time import localtime, strftime, gmtime
24
25
26
27#-------------------------------------------------------------------------------
28# Copy scripts to time stamped output directory and capture screen
29# output to file
30#-------------------------------------------------------------------------------
31time = strftime('%Y%m%d_%H%M%S',localtime())
32
33output_dir = 'dam_break_'+time
34output_file = 'dam_break'
35
36copy_code_files(output_dir,__file__)
37#start_screen_catcher(output_dir+'_')
38
39
40#------------------------------------------------------------------------------
41# Setup domain
42#------------------------------------------------------------------------------
43dx = 200.
44dy = dx
45L = 100000.
46W = 10*dx
47
48# structured mesh
49points, vertices, boundary = rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, -W/2))
50
51domain = Domain(points, vertices, boundary) 
52
53domain.set_name(output_file)               
54domain.set_datadir(output_dir) 
55
56#------------------------------------------------------------------------------
57# Setup Algorithm
58#------------------------------------------------------------------------------
59domain.set_timestepping_method('rk2')
60domain.set_default_order(2)
61
62print domain.get_timestepping_method()
63
64domain.use_edge_limiter = True
65domain.tight_slope_limiters = True
66domain.use_centroid_velocities = False
67
68domain.CFL = 1.0
69
70domain.beta_w      = 0.6
71domain.beta_uh     = 0.6
72domain.beta_vh     = 0.6
73
74
75#------------------------------------------------------------------------------
76# Setup initial conditions
77#------------------------------------------------------------------------------
78domain.set_quantity('elevation',0.0)
79domain.set_quantity('friction', 0.0)
80
81h0 = 10.0
82h1 = 0.0
83
84def height(x,y):
85    z = zeros(len(x),Float)
86    for i in range(len(x)):
87        if x[i]<=50000.0:
88            z[i] = h0
89        else:
90            z[i] = h1
91    return z
92
93
94domain.set_quantity('stage', height)
95#-----------------------------------------------------------------------------
96# Setup boundary conditions
97#------------------------------------------------------------------------------
98from math import sin, pi, exp
99Br = Reflective_boundary(domain)      # Solid reflective wall
100Bt = Transmissive_boundary(domain)    # Continue all values on boundary
101Bd = Dirichlet_boundary([1,0.,0.]) # Constant boundary values
102
103
104# Associate boundary tags with boundary objects
105domain.set_boundary({'left': Bt, 'right': Bt, 'top': Br, 'bottom': Br})
106
107
108#===============================================================================
109from anuga.visualiser import RealtimeVisualiser
110vis = RealtimeVisualiser(domain)
111vis.render_quantity_height("stage", zScale =h0*500, dynamic=True)
112vis.colour_height_quantity('stage', (0.0, 0.5, 1.0))
113vis.start()
114#===============================================================================
115
116
117#------------------------------------------------------------------------------
118# Evolve system through time
119#------------------------------------------------------------------------------
120
121for t in domain.evolve(yieldstep = 100.0, finaltime = 60*60.):
122    print domain.timestepping_statistics(track_speeds=True)
123    vis.update()
124   
125vis.evolveFinished()
126
Note: See TracBrowser for help on using the repository browser.