source: anuga_validation/analytical solutions/carrier_wave_runup.py @ 7626

Last change on this file since 7626 was 7626, checked in by steve, 14 years ago

Cleaning up directory

File size: 4.3 KB
RevLine 
[2229]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    if domain.get_vertex_coordinate(id,(face+1)%3)[0] < -200.0+1.0e-10:
52        domain.boundary[(id,face)] = 'left'
53
54
55
56
57# Provide file name for storing output
58domain.store = True
59domain.format = 'sww'
[3846]60domain.set_name('run-up_second_order_again')
[2229]61
62print "Number of triangles = ", len(domain)
63
64#Reduction operation for get_vertex_values
[3514]65from anuga.pyvolution.util import mean
[2229]66domain.reduction = mean
67#domain.reduction = min  #Looks better near steep slopes
68
69#Define the boundary condition
70def stage_setup(x,t_star1):
71    vh = 0
72    alpha = 0.1
73    eta = 0.1
74    a = 1.5*sqrt(1.+0.9*eta)
75    l_0 = 200.
76    ii = complex(0,1)
77    g = 9.81
78    v_0 = sqrt(g*l_0*alpha)
79    v1 = 0.
80
81    sigma_max = 100.
82    sigma_min = -100.
83    for j in range (1,50):
84        sigma0 = (sigma_max+sigma_min)/2.
85        lambda_prime = 2./a*(t_star1/sqrt(l_0/alpha/g)+v1)
86        sigma_prime = sigma0/a
87        const = (1.-ii*lambda_prime)**2+sigma_prime**2
88
89        v1 = 8.*eta/a*imag(1./const**(3./2.) \
90            -3./4.*(1.-ii*lambda_prime)/const**(5./2.))
91
92        x1 = -v1**2/2.-a**2*sigma_prime**2/16.+eta \
93            *real(1.-2.*(5./4.-ii*lambda_prime) \
94            /const**(3./2.)+3./2.*(1.-ii*lambda_prime)**2 \
95            /const**(5./2.))
96
97        neta1 = x1 + a*a*sigma_prime**2/16.
98
99        v_star1 = v1*v_0
100        x_star1 = x1*l_0
101        neta_star1 = neta1*alpha*l_0
102        stage = neta_star1
103        z = stage - x_star1*alpha
104        uh = z*v_star1
105
106        if x_star1-x > 0:
107            sigma_max = sigma0
108        else:
109            sigma_min = sigma0
110
111        if abs(abs(sigma0)-100.) < 10:
112
113#     solution does not converge because bed is dry
114            stage = 0.
115            uh = 0.
116            z = 0.
117
118    return [stage, uh, vh]
119
120def boundary_stage(t):
121    x = -200
122    return stage_setup(x,t)
123
124
125######################
126#Initial condition
127print 'Initial condition'
128t_star1 = 0.0
129slope = -0.1
130
131#Set bed-elevation and friction(None)
132def x_slope(x,y):
133    n = x.shape[0]
134    z = 0*x
135    for i in range(n):
136        z[i] = -slope*x[i]
137    return z
138
139domain.set_quantity('elevation', x_slope)
140
141#Set the water depth
142print 'Initial water depth'
143def stage(x,y):
144    z = x_slope(x,y)
145    n = x.shape[0]
146    w = 0*x
147    for i in range(n):
148        w[i], uh, vh = stage_setup(x[i],t_star1)
149        h = w[i] - z[i]
150        if h < 0:
151            h = 0
152            w[i] = z[i]
153    return w
154
155domain.set_quantity('stage', stage)
156
157
158#####################
159#Set up boundary conditions
160Br = Reflective_boundary(domain)
161Bw = Time_boundary(domain, boundary_stage)
162domain.set_boundary({'left': Bw, 'external': Br})
163
164
165
166#domain.visualise = True
167
168######################
169#Evolution
170import time
171t0 = time.time()
172for t in domain.evolve(yieldstep = 1., finaltime = 100):
173    domain.write_time()
174    print boundary_stage(domain.time)
175
176print 'That took %.2f seconds' %(time.time()-t0)
177
Note: See TracBrowser for help on using the repository browser.