source: anuga_validation/convergence_study/wave.py @ 5162

Last change on this file since 5162 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.8 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
13from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
14from anuga.shallow_water import Domain
15from anuga.shallow_water import Reflective_boundary
16from anuga.shallow_water import Dirichlet_boundary
17from anuga.shallow_water import Time_boundary
18from anuga.shallow_water import Transmissive_boundary
19from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
20from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files
21
22from math import cos
23from Numeric import zeros, Float
24from time import localtime, strftime, gmtime
25from os import sep
26
27
28
29#-------------------------------------------------------------------------------
30# Copy scripts to time stamped output directory and capture screen
31# output to file
32#-------------------------------------------------------------------------------
33time = strftime('%Y%m%d_%H%M%S',localtime())
34
35output_dir = 'wave_'+time
36output_file = 'wave'
37
38copy_code_files(output_dir,__file__)
39#start_screen_catcher(output_dir+sep)
40
41#------------------------------------------------------------------------------
42# Setup domain
43#------------------------------------------------------------------------------
44dx = 500.
45dy = dx
46L = 100000.
47W = 10*dx
48
49# structured mesh
50points, vertices, boundary = rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, -W/2))
51
52domain = Domain(points, vertices, boundary) 
53
54domain.set_name(output_file)               
55domain.set_datadir(output_dir) 
56
57#------------------------------------------------------------------------------
58# Setup Algorithm
59#------------------------------------------------------------------------------
60#domain.set_timestepping_method('euler')
61
62#domain.set_default_order(2)
63#domain.set_timestepping_method('euler')
64#domain.set_all_limiters(2.0)
65#domain.use_old_limiter = True
66#domain.CFL = 1.0
67
68
69domain.set_default_order(2)
70#domain.set_timestepping_method('euler')
71
72domain.use_old_limiter = True
73
74domain.beta_w      = 2.0
75domain.beta_w_dry  = 0.5
76domain.beta_uh     = 2.0
77domain.beta_uh_dry = 0.5
78domain.beta_vh     = 2.0
79domain.beta_vh_dry = 0.5
80
81
82#------------------------------------------------------------------------------
83# Setup initial conditions
84#------------------------------------------------------------------------------
85domain.set_quantity('elevation',-100.0)
86domain.set_quantity('friction', 0.00)
87domain.set_quantity('stage', 0.0)           
88
89#-----------------------------------------------------------------------------
90# Setup boundary conditions
91#------------------------------------------------------------------------------
92from math import sin, pi, exp
93Br = Reflective_boundary(domain)      # Solid reflective wall
94Bt = Transmissive_boundary(domain)    # Continue all values on boundary
95Bd = Dirichlet_boundary([1,0.,0.]) # Constant boundary values
96amplitude = 1
97Bw = Time_boundary(domain=domain,     # Time dependent boundary 
98## Sine wave
99                  f=lambda t: [(-amplitude*sin((1./300.)*t*2*pi)), 0.0, 0.0])
100## Sawtooth?
101#                   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])
102## Sharp rise, linear fall
103#                   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])
104#                   f=lambda t: [amplitude*(1.-2.*(pi*(1./720.)*(t-720.))**2)/exp((pi*(1./720.)*(t-720.))**2) , 0.0, 0.0])
105#                   f=lambda t: [(-8.0*sin((1./720.)*t*2*pi))*((t<720.)-0.5*(t<360.)), 0.0, 0.0])
106
107# Associate boundary tags with boundary objects
108domain.set_boundary({'left': Bw, 'right': Bt, 'top': Br, 'bottom': Br})
109
110
111#===============================================================================
112from anuga.visualiser import RealtimeVisualiser
113vis = RealtimeVisualiser(domain)
114vis.render_quantity_height("stage", zScale =10000, dynamic=True)
115vis.colour_height_quantity('stage', (1.0, 0.5, 0.5))
116vis.start()
117#===============================================================================
118
119
120#------------------------------------------------------------------------------
121# Evolve system through time
122#------------------------------------------------------------------------------
123
124for t in domain.evolve(yieldstep = 50.0, finaltime = 60*60.):
125    domain.write_time()
126    vis.update()
127   
128vis.evolveFinished()
129
Note: See TracBrowser for help on using the repository browser.