source: inundation/ga/storm_surge/analytical solutions/Analytical solution_wave_runup.py @ 1701

Last change on this file since 1701 was 1701, checked in by chris, 19 years ago
File size: 4.0 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 os import sep
20sys.path.append('..'+sep+'pyvolution')
21
22from shallow_water import Domain, Reflective_boundary, Time_boundary
23from math import sqrt, cos, sin, pi
24from mesh_factory import strang_mesh
25
26#Convenience functions
27def imag(a):
28    return a.imag
29
30def real(a):
31    return a.real
32
33
34######################
35# Domain
36# Strang_domain will search through the file and test to see if there are
37# two or three entries. Two entries are for points and three for triangles.
38
39#points, elements = strang_mesh('Run-up.pt')
40points, elements = strang_mesh('strang_7389.pt')
41domain = Domain(points, elements) 
42
43domain.default_order = 2
44domain.smooth = True
45
46#Set a default tagging
47epsilon = 1.0e-12
48
49for id, face in domain.boundary:
50    domain.boundary[(id,face)] = 'external'
51   
52domain.boundary[(0,0)] = 'left'
53
54
55# Provide file name for storing output
56domain.store = True
57domain.format = 'sww'
58domain.filename = 'run-up_second_order_again'
59
60print "Number of triangles = ", len(domain)
61
62#Reduction operation for get_vertex_values             
63from util import mean
64domain.reduction = mean
65#domain.reduction = min  #Looks better near steep slopes
66
67#Define the boundary condition
68def stage_setup(x,t_star1):
69    vh = 0
70    alpha = 0.1
71    eta = 0.1
72    a = 1.5*sqrt(1.+0.9*eta)
73    l_0 = 200.
74    ii = complex(0,1)
75    g = 9.81
76    v_0 = sqrt(g*l_0*alpha)
77    v1 = 0.
78
79    sigma_max = 100.
80    sigma_min = -100.
81    for j in range (1,50):
82        sigma0 = (sigma_max+sigma_min)/2.
83        lambda_prime = 2./a*(t_star1/sqrt(l_0/alpha/g)+v1)
84        sigma_prime = sigma0/a
85        const = (1.-ii*lambda_prime)**2+sigma_prime**2
86         
87        v1 = 8.*eta/a*imag(1./const**(3./2.) \
88            -3./4.*(1.-ii*lambda_prime)/const**(5./2.))
89           
90        x1 = -v1**2/2.-a**2*sigma_prime**2/16.+eta \
91            *real(1.-2.*(5./4.-ii*lambda_prime) \
92            /const**(3./2.)+3./2.*(1.-ii*lambda_prime)**2 \
93            /const**(5./2.))
94
95        neta1 = x1 + a*a*sigma_prime**2/16.
96
97        v_star1 = v1*v_0
98        x_star1 = x1*l_0
99        neta_star1 = neta1*alpha*l_0
100        stage = neta_star1
101        z = stage - x_star1*alpha
102        uh = z*v_star1
103
104        if x_star1-x > 0:
105            sigma_max = sigma0
106        else:
107            sigma_min = sigma0
108
109        if abs(abs(sigma0)-100.) < 10:
110
111#     solution does not converge because bed is dry
112            stage = 0.
113            uh = 0.
114            z = 0.
115
116    return [stage, uh, vh]
117
118def boundary_stage(t):
119    x = -200
120    return stage_setup(x,t)
121
122
123######################
124#Initial condition
125print 'Initial condition'
126t_star1 = 0.0
127slope = -0.1
128
129#Set bed-elevation and friction(None)
130def x_slope(x,y):
131    n = x.shape[0]
132    z = 0*x
133    for i in range(n):
134        z[i] = -slope*x[i]
135    return z
136
137domain.set_quantity('elevation', x_slope)
138
139#Set the water depth
140print 'Initial water depth'
141def stage(x,y):
142    z = x_slope(x,y)
143    n = x.shape[0]
144    w = 0*x
145    for i in range(n):
146        w[i], uh, vh = stage_setup(x[i],t_star1)
147        h = w[i] - z[i]
148        if h < 0:
149            h = 0
150            w[i] = z[i]
151    return w   
152
153domain.set_quantity('stage', stage)
154
155
156#####################
157#Set up boundary conditions
158Br = Reflective_boundary(domain)
159Bw = Time_boundary(domain, boundary_stage)
160domain.set_boundary({'left': Bw, 'external': Br})
161
162
163######################
164#Evolution
165import time
166t0 = time.time()
167for t in domain.evolve(yieldstep = 1., finaltime = 100):
168    domain.write_time()
169    print boundary_stage(domain.time)
170   
171print 'That took %.2f seconds' %(time.time()-t0)
172
173   
Note: See TracBrowser for help on using the repository browser.