source: trunk/anuga_validation/validation_tests/periodic_carrier/carrier_wave_runup.py @ 8067

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

Pulling together basic validate tests for anuga

File size: 5.5 KB
Line 
1"""Example of shallow water wave equation analytical solution of the
2one-dimensional Thacker and Greenspan wave run-up treated as a two-dimensional solution.
3
4   Copyright 2004
5   Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray
6   Geoscience Australia
7
8Specific methods pertaining to the 2D shallow water equation
9are imported from shallow_water
10for use with the generic finite volume framework
11
12Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the
13numerical vector named conserved_quantities.
14"""
15
16#-------------------
17# Module imports
18import sys
19from math import sqrt, cos, sin, pi
20from numpy import asarray
21from time import localtime, strftime, gmtime
22
23
24#-------------------
25# Anuga Imports
26#-------------------
27from anuga.shallow_water_balanced.swb_domain import Domain, Transmissive_boundary, Reflective_boundary,\
28     Dirichlet_boundary, Time_boundary
29
30#from anuga.interface import Domain, Transmissive_boundary, Reflective_boundary,\
31#     Dirichlet_boundary
32
33from anuga.utilities.polygon import inside_polygon, is_inside_triangle
34from anuga.abstract_2d_finite_volumes.mesh_factory import strang_mesh
35
36
37#-------------------
38#Convenience functions
39#-------------------
40def imag(a):
41    return a.imag
42
43def real(a):
44    return a.real
45
46
47#--------------------
48# Domain
49# Strang_domain will search through the file and test to see if there are
50# two or three entries. Two entries are for points and three for triangles.
51
52#points, elements = strang_mesh('Run-up.pt')
53points, elements = strang_mesh('strang_7389.pt')
54domain = Domain(points, elements)
55
56print 'extent ',domain.get_extent()
57
58#----------------
59# Order of scheme
60# Good compromise between
61# limiting and CFL
62#---------------
63domain.set_default_order(2)
64domain.set_timestepping_method(2)
65domain.set_beta(0.7)
66domain.set_CFL(0.6)
67
68
69#-------------------
70#Set a default tagging
71#-------------------
72epsilon = 1.0e-12
73
74for id, face in domain.boundary:
75    domain.boundary[(id,face)] = 'external'
76    if domain.get_vertex_coordinate(id,(face+1)%3)[0] < -200.0+1.0e-10:
77        domain.boundary[(id,face)] = 'left'
78
79#---------------------
80# Provide file name for storing output
81#---------------------
82domain.store = True
83domain.format = 'sww'
84
85time = strftime('%Y%m%d_%H%M%S',localtime())
86output_file= 'carrier_wave_runup_'+time
87
88domain.set_name(output_file)
89
90print "Number of triangles = ", len(domain)
91
92#-----------------------
93#Define the boundary condition
94#-----------------------
95def stage_setup(x,t):
96    vh = 0
97    alpha = 0.1
98    eta = 0.1
99    a = 1.5*sqrt(1.+0.9*eta)
100    l_0 = 200.
101    ii = complex(0,1)
102    g = 9.81
103    v_0 = sqrt(g*l_0*alpha)
104    v1 = 0.
105
106    sigma_max = 100.
107    sigma_min = -100.
108    for j in range (1,50):
109        sigma0 = (sigma_max+sigma_min)/2.
110        lambda_prime = 2./a*(t/sqrt(l_0/alpha/g)+v1)
111        sigma_prime = sigma0/a
112        const = (1.-ii*lambda_prime)**2+sigma_prime**2
113
114        v1 = 8.*eta/a*imag(1./const**(3./2.) \
115            -3./4.*(1.-ii*lambda_prime)/const**(5./2.))
116
117        x1 = -v1**2/2.-a**2*sigma_prime**2/16.+eta \
118            *real(1.-2.*(5./4.-ii*lambda_prime) \
119            /const**(3./2.)+3./2.*(1.-ii*lambda_prime)**2 \
120            /const**(5./2.))
121
122        neta1 = x1 + a*a*sigma_prime**2/16.
123
124        v_star1 = v1*v_0
125        x_star1 = x1*l_0
126        neta_star1 = neta1*alpha*l_0
127        stage = neta_star1
128        z = stage - x_star1*alpha
129        uh = z*v_star1
130
131        if x_star1-x > 0:
132            sigma_max = sigma0
133        else:
134            sigma_min = sigma0
135
136        if abs(abs(sigma0)-100.) < 10:
137
138#     solution does not converge because bed is dry
139            stage = 0.
140            uh = 0.
141            z = 0.
142
143    return [stage, uh, vh]
144
145def boundary_stage(t):
146    x = -200
147    return stage_setup(x,t)
148
149
150#---------------------
151#Initial condition
152#---------------------
153print 'Initial condition'
154t_star1 = 0.0
155slope = -0.1
156
157#Set bed-elevation and friction(None)
158def x_slope(x,y):
159    n = x.shape[0]
160    z = 0*x
161    for i in range(n):
162        z[i] = -slope*x[i]
163    return z
164
165domain.set_quantity('elevation', x_slope)
166
167#Set the water depth
168print 'Initial water depth'
169def stage(x,y):
170    z = x_slope(x,y)
171    n = x.shape[0]
172    w = 0*x
173    for i in range(n):
174        w[i], uh, vh = stage_setup(x[i],t_star1)
175        h = w[i] - z[i]
176        if h < 0:
177            h = 0
178            w[i] = z[i]
179    return w
180
181domain.set_quantity('stage', stage)
182
183
184#####################
185#Set up boundary conditions
186Br = Reflective_boundary(domain)
187Bw = Time_boundary(domain, boundary_stage)
188domain.set_boundary({'left': Bw, 'external': Br})
189
190import time
191visualize = True
192if visualize:
193    from anuga.visualiser import RealtimeVisualiser
194    vis = RealtimeVisualiser(domain)
195    vis.render_quantity_height("elevation", zScale=3.0, offset = 0.01, dynamic=False)
196    vis.render_quantity_height("stage", zScale = 3.0, dynamic=True, opacity = 0.6, wireframe=False)
197    #vis.colour_height_quantity('stage', (lambda q:q['stage'], 1.0, 2.0))
198    vis.colour_height_quantity('stage', (0.4, 0.6, 0.4))
199    vis.start()
200    time.sleep(2.0)
201
202#domain.visualise = True
203
204######################
205#Evolution
206
207t0 = time.time()
208
209
210for t in domain.evolve(yieldstep = 1., finaltime = 100):
211    domain.write_time()
212    print boundary_stage(domain.time)
213    if visualize: vis.update()
214   
215if visualize: vis.evolveFinished()   
216
217print 'That took %.2f seconds' %(time.time()-t0)
218
Note: See TracBrowser for help on using the repository browser.