[547] | 1 | """Example of shallow water wave equation analytical solution of the |
---|
| 2 | one-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 | |
---|
| 8 | Specific methods pertaining to the 2D shallow water equation |
---|
| 9 | are imported from shallow_water |
---|
| 10 | for use with the generic finite volume framework |
---|
| 11 | |
---|
| 12 | Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the |
---|
| 13 | numerical vector named conserved_quantities. |
---|
| 14 | """ |
---|
| 15 | |
---|
| 16 | ###################### |
---|
| 17 | # Module imports |
---|
| 18 | import sys |
---|
| 19 | from os import sep |
---|
| 20 | sys.path.append('..'+sep+'pyvolution') |
---|
| 21 | |
---|
[550] | 22 | from shallow_water import Domain, Reflective_boundary, Time_boundary |
---|
[547] | 23 | from math import sqrt, cos, sin, pi |
---|
[550] | 24 | from mesh_factory import strang_mesh |
---|
[547] | 25 | |
---|
[550] | 26 | #Convenience functions |
---|
[547] | 27 | def imag(a): |
---|
| 28 | return a.imag |
---|
| 29 | |
---|
| 30 | def 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 | |
---|
[550] | 39 | points, elements = strang_mesh('run-up.pt') |
---|
| 40 | domain = Domain(points, elements) |
---|
[547] | 41 | |
---|
| 42 | domain.default_order = 2 |
---|
| 43 | domain.smooth = True |
---|
| 44 | |
---|
| 45 | #Set a default tagging |
---|
| 46 | epsilon = 1.0e-12 |
---|
| 47 | |
---|
| 48 | for id, face in domain.boundary: |
---|
| 49 | domain.boundary[(id,face)] = 'external' |
---|
| 50 | |
---|
| 51 | domain.boundary[(0,0)] = 'left' |
---|
| 52 | |
---|
| 53 | |
---|
| 54 | # Provide file name for storing output |
---|
[550] | 55 | domain.store = True |
---|
| 56 | domain.format = 'sww' |
---|
| 57 | domain.filename = 'run-up_second_order' |
---|
[547] | 58 | |
---|
[550] | 59 | print "Number of triangles = ", len(domain) |
---|
[547] | 60 | |
---|
| 61 | #Reduction operation for get_vertex_values |
---|
[550] | 62 | from util import mean |
---|
[547] | 63 | domain.reduction = mean |
---|
| 64 | #domain.reduction = min #Looks better near steep slopes |
---|
| 65 | |
---|
| 66 | #Define the boundary condition |
---|
| 67 | def stage_setup(x,t_star1): |
---|
| 68 | vh = 0 |
---|
| 69 | alpha = 0.1 |
---|
| 70 | eta = 0.1 |
---|
| 71 | a = 1.5*sqrt(1.+0.9*eta) |
---|
| 72 | l_0 = 200. |
---|
| 73 | ii = complex(0,1) |
---|
| 74 | g = 9.81 |
---|
| 75 | v_0 = sqrt(g*l_0*alpha) |
---|
| 76 | v1 = 0. |
---|
| 77 | |
---|
| 78 | sigma_max = 100. |
---|
| 79 | sigma_min = -100. |
---|
| 80 | for j in range (1,50): |
---|
| 81 | sigma0 = (sigma_max+sigma_min)/2. |
---|
| 82 | lambda_prime = 2./a*(t_star1/sqrt(l_0/alpha/g)+v1) |
---|
| 83 | sigma_prime = sigma0/a |
---|
| 84 | const = (1.-ii*lambda_prime)**2+sigma_prime**2 |
---|
| 85 | |
---|
| 86 | v1 = 8.*eta/a*imag(1./const**(3./2.) \ |
---|
| 87 | -3./4.*(1.-ii*lambda_prime)/const**(5./2.)) |
---|
| 88 | |
---|
| 89 | x1 = -v1**2/2.-a**2*sigma_prime**2/16.+eta \ |
---|
| 90 | *real(1.-2.*(5./4.-ii*lambda_prime) \ |
---|
| 91 | /const**(3./2.)+3./2.*(1.-ii*lambda_prime)**2 \ |
---|
| 92 | /const**(5./2.)) |
---|
| 93 | |
---|
| 94 | neta1 = x1 + a*a*sigma_prime**2/16. |
---|
| 95 | |
---|
| 96 | v_star1 = v1*v_0 |
---|
| 97 | x_star1 = x1*l_0 |
---|
| 98 | neta_star1 = neta1*alpha*l_0 |
---|
| 99 | stage = neta_star1 |
---|
| 100 | z = stage - x_star1*alpha |
---|
| 101 | uh = z*v_star1 |
---|
| 102 | |
---|
| 103 | if x_star1-x > 0: |
---|
| 104 | sigma_max = sigma0 |
---|
| 105 | else: |
---|
| 106 | sigma_min = sigma0 |
---|
| 107 | |
---|
| 108 | if abs(abs(sigma0)-100.) < 10: |
---|
| 109 | |
---|
| 110 | # solution does not converge because bed is dry |
---|
| 111 | stage = 0. |
---|
| 112 | uh = 0. |
---|
| 113 | z = 0. |
---|
| 114 | |
---|
| 115 | return [stage, uh, vh] |
---|
| 116 | |
---|
[774] | 117 | def boundary_stage(t): |
---|
[547] | 118 | x = -200 |
---|
| 119 | return stage_setup(x,t) |
---|
| 120 | |
---|
| 121 | |
---|
| 122 | ###################### |
---|
| 123 | #Initial condition |
---|
| 124 | print 'Initial condition' |
---|
| 125 | t_star1 = 0.0 |
---|
| 126 | slope = -0.1 |
---|
| 127 | |
---|
| 128 | #Set bed-elevation and friction(None) |
---|
| 129 | def x_slope(x,y): |
---|
| 130 | n = x.shape[0] |
---|
| 131 | z = 0*x |
---|
| 132 | for i in range(n): |
---|
| 133 | z[i] = -slope*x[i] |
---|
| 134 | return z |
---|
| 135 | |
---|
[550] | 136 | domain.set_quantity('elevation', x_slope) |
---|
[547] | 137 | |
---|
| 138 | #Set the water depth |
---|
[774] | 139 | def stage(x,y): |
---|
[547] | 140 | z = x_slope(x,y) |
---|
| 141 | n = x.shape[0] |
---|
| 142 | w = 0*x |
---|
| 143 | for i in range(n): |
---|
| 144 | w[i], uh, vh = stage_setup(x[i],t_star1) |
---|
| 145 | h = w[i] - z[i] |
---|
| 146 | if h < 0: |
---|
| 147 | h = 0 |
---|
| 148 | w[i] = z[i] |
---|
| 149 | return w |
---|
| 150 | |
---|
[774] | 151 | domain.set_quantity('stage', stage) |
---|
[547] | 152 | |
---|
[550] | 153 | |
---|
[547] | 154 | ##################### |
---|
| 155 | #Set up boundary conditions |
---|
[550] | 156 | Br = Reflective_boundary(domain) |
---|
[774] | 157 | Bw = Time_boundary(domain, boundary_stage) |
---|
[547] | 158 | domain.set_boundary({'left': Bw, 'external': Br}) |
---|
| 159 | |
---|
| 160 | |
---|
| 161 | ###################### |
---|
| 162 | #Evolution |
---|
| 163 | import time |
---|
| 164 | t0 = time.time() |
---|
[550] | 165 | for t in domain.evolve(yieldstep = 1., finaltime = 100): |
---|
[547] | 166 | domain.write_time() |
---|
[774] | 167 | print boundary_stage(domain.time) |
---|
[547] | 168 | |
---|
| 169 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
| 170 | |
---|
| 171 | |
---|