source: trunk/anuga_work/shallow_water_balanced_steve/run_wave_shelf.py @ 8288

Last change on this file since 8288 was 8203, checked in by steve, 13 years ago

Now swb_domain is calling a well balanced gravity calculation

File size: 4.6 KB
Line 
1"""Simple water flow example using ANUGA
2
3Oscillatory wave from off shore coming into shore
4"""
5
6#------------------------------------------------------------------------------
7# Import necessary modules
8#------------------------------------------------------------------------------
9
10import sys
11import anuga
12from swb_domain import *
13
14from math import cos
15import numpy as num
16from time import localtime, strftime, gmtime
17from os import sep
18
19
20
21#-------------------------------------------------------------------------------
22# Copy scripts to time stamped output directory and capture screen
23# output to file
24#-------------------------------------------------------------------------------
25time = strftime('%Y%m%d_%H%M%S',localtime())
26
27output_dir = '.'
28output_file = 'data_wave_shelf_'+time
29
30#copy_code_files(output_dir,__file__)
31#start_screen_catcher(output_dir+sep)
32
33interactive_visualisation = True
34
35#------------------------------------------------------------------------------
36# Setup domain
37#------------------------------------------------------------------------------
38dx = 1000.
39dy = dx
40L = 150000.
41W = 10*dx
42
43# structured mesh
44points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, -W/2))
45
46domain = Domain(points, vertices, boundary)
47
48domain.set_name(output_file)               
49domain.set_datadir(output_dir) 
50
51#------------------------------------------------------------------------------
52# Setup Algorithm
53#------------------------------------------------------------------------------
54domain.set_timestepping_method('rk2')
55domain.set_default_order(2)
56
57print domain.get_timestepping_method()
58
59#domain.use_edge_limiter = True
60#domain.tight_slope_limiters = False
61#domain.use_centroid_velocities = False
62
63
64#------------------------------------------------------------------------------
65# Setup initial conditions
66#------------------------------------------------------------------------------
67depth1 = -100.0
68depth2 = -30.00
69def topography(x,y):
70    z = num.zeros_like(x)
71    z[:] = depth1
72
73    N = len(x)
74    for i in range(N):
75        if  40000.0 < x[i] < 60000.0:
76            z[i] = depth1 + (depth2-depth1)*((x[i] - 40000.0)/20000.0)
77
78        if x[i] > 60000.0:
79            z[i] = depth2
80
81        if x[i] > 100000.0:
82            z[i] = depth2*(1.0 - 1.1*(x[i] - 100000.0)/50000)
83    return z
84
85domain.set_quantity('elevation',topography)
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 = anuga.Reflective_boundary(domain)      # Solid reflective wall
94Bt = anuga.Transmissive_boundary(domain)    # Continue all values on boundary
95Bd = anuga.Dirichlet_boundary([1,0.,0.]) # Constant boundary values
96amplitude = 1
97wave_length = 2000.0
98Bw = anuga.Time_boundary(domain=domain,     # Time dependent boundary
99## Sine wave
100                  f=lambda t: [(amplitude*sin((1./wave_length)*t*2*pi)), 0.0, 0.0])
101## Sawtooth?
102#                   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])
103## Sharp rise, linear fall
104#                   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])
105#                   f=lambda t: [amplitude*(1.-2.*(pi*(1./720.)*(t-720.))**2)/exp((pi*(1./720.)*(t-720.))**2) , 0.0, 0.0])
106#                   f=lambda t: [(-8.0*sin((1./720.)*t*2*pi))*((t<720.)-0.5*(t<360.)), 0.0, 0.0])
107
108# Associate boundary tags with boundary objects
109domain.set_boundary({'left': Bw, 'right': Br, 'top': Br, 'bottom': Br})
110
111
112#===============================================================================
113if interactive_visualisation:
114    from anuga.visualiser import RealtimeVisualiser
115    vis = RealtimeVisualiser(domain)
116    vis.render_quantity_height("stage", zScale =10000, dynamic=True)
117    vis.colour_height_quantity('stage', (1.0, 0.5, 0.5))
118    vis.start()
119#===============================================================================
120
121
122#------------------------------------------------------------------------------
123# Evolve system through time
124#------------------------------------------------------------------------------
125
126min = 60.0
127hr = 60*min
128
129for t in domain.evolve(yieldstep = min, finaltime = 5*hr):
130    domain.write_time()
131    if interactive_visualisation:
132        vis.update()
133
134    if domain.get_time()> 2*hr:
135        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
136
137if interactive_visualisation:
138    vis.evolveFinished()
139
Note: See TracBrowser for help on using the repository browser.