source: trunk/anuga_validation/validation_tests/periodic_forcing/wave.py @ 8078

Last change on this file since 8078 was 8078, checked in by habili, 9 years ago

updated code so that wave.py runs

File size: 4.4 KB
Line 
1"""Simple water flow example using ANUGA
2
3Will Powers example of a simple sinusoidal wave which showed diffusive effects of
4thefirst order and standard second order method. Problem resolved if "rk2" timestepping
5and higher beta = 2 limiter used. Also new edge limiter with rk2 resolves problem
6"""
7
8#------------------------------------------------------------------------------
9# Import necessary modules
10#------------------------------------------------------------------------------
11
12import sys
13import anuga
14from math import cos
15from time import localtime, strftime, gmtime
16from os import sep
17
18
19
20#-------------------------------------------------------------------------------
21# Copy scripts to time stamped output directory and capture screen
22# output to file
23#-------------------------------------------------------------------------------
24time = strftime('%Y%m%d_%H%M%S',localtime())
25
26output_dir = 'wave_'+time
27output_file = 'wave'
28
29#copy_code_files(output_dir,__file__)
30#start_screen_catcher(output_dir+sep)
31
32interactive_visualisation = False
33
34#------------------------------------------------------------------------------
35# Setup domain
36#------------------------------------------------------------------------------
37dx = 1000.
38dy = dx
39L = 100000.
40W = 10*dx
41
42# structured mesh
43points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, -W/2))
44
45domain = anuga.Domain(points, vertices, boundary) 
46
47domain.set_name(output_file)               
48domain.set_datadir(output_dir) 
49
50#------------------------------------------------------------------------------
51# Setup Algorithm
52#------------------------------------------------------------------------------
53domain.set_timestepping_method('rk2')
54domain.set_default_order(2)
55
56print domain.get_timestepping_method()
57
58domain.use_edge_limiter = True
59domain.tight_slope_limiters = False
60domain.use_centroid_velocities = False
61
62domain.CFL = 1.0
63
64domain.beta_w      = 1.0
65domain.beta_w_dry  = 0.0
66domain.beta_uh     = 1.0
67domain.beta_uh_dry = 0.0
68domain.beta_vh     = 1.0
69domain.beta_vh_dry = 0.0
70
71
72#------------------------------------------------------------------------------
73# Setup initial conditions
74#------------------------------------------------------------------------------
75domain.set_quantity('elevation',-100.0)
76domain.set_quantity('friction', 0.00)
77domain.set_quantity('stage', 0.0)           
78
79#-----------------------------------------------------------------------------
80# Setup boundary conditions
81#------------------------------------------------------------------------------
82from math import sin, pi, exp
83Br = anuga.Reflective_boundary(domain)      # Solid reflective wall
84Bt = anuga.Transmissive_boundary(domain)    # Continue all values on boundary
85Bd = anuga.Dirichlet_boundary([1,0.,0.]) # Constant boundary values
86amplitude = 1
87Bw = anuga.Time_boundary(domain=domain,     # Time dependent boundary 
88## Sine wave
89                  f=lambda t: [(-amplitude*sin((1./300.)*t*2*pi)), 0.0, 0.0])
90## Sawtooth?
91#                   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])
92## Sharp rise, linear fall
93#                   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])
94#                   f=lambda t: [amplitude*(1.-2.*(pi*(1./720.)*(t-720.))**2)/exp((pi*(1./720.)*(t-720.))**2) , 0.0, 0.0])
95#                   f=lambda t: [(-8.0*sin((1./720.)*t*2*pi))*((t<720.)-0.5*(t<360.)), 0.0, 0.0])
96
97# Associate boundary tags with boundary objects
98domain.set_boundary({'left': Bw, 'right': Bt, 'top': Br, 'bottom': Br})
99
100
101#===============================================================================
102if interactive_visualisation:
103    from anuga.visualiser import RealtimeVisualiser
104    vis = RealtimeVisualiser(domain)
105    vis.render_quantity_height("stage", zScale =10000, dynamic=True)
106    vis.colour_height_quantity('stage', (1.0, 0.5, 0.5))
107    vis.start()
108#===============================================================================
109
110
111#------------------------------------------------------------------------------
112# Evolve system through time
113#------------------------------------------------------------------------------
114
115for t in domain.evolve(yieldstep = 50.0, finaltime = 60*60.):
116    domain.write_time()
117    if interactive_visualisation:
118        vis.update()
119
120if interactive_visualisation:
121    vis.evolveFinished()
122
Note: See TracBrowser for help on using the repository browser.