[5535] | 1 | import os |
---|
[5827] | 2 | #import Gnuplot |
---|
[5535] | 3 | from math import sqrt, pi |
---|
[5827] | 4 | from shallow_water_domain import * |
---|
[5535] | 5 | from Numeric import allclose, array, zeros, ones, Float, take, sqrt |
---|
| 6 | from config import g, epsilon |
---|
| 7 | |
---|
| 8 | |
---|
| 9 | #h_0 = 0.48 |
---|
| 10 | #u_0 = 0.625 |
---|
| 11 | #h_1 = 0.48 |
---|
| 12 | #u_1 = 0.625 |
---|
| 13 | |
---|
| 14 | #h_0 = 0.48 |
---|
| 15 | #u_0 = 0.625 |
---|
| 16 | #h_1 = 0.106 |
---|
| 17 | #u_1 = 2.829228 |
---|
| 18 | |
---|
| 19 | h_0 = 0.5 |
---|
| 20 | u_0 = 0.6 |
---|
| 21 | h_1 = 0.5 |
---|
| 22 | u_1 = 0.6 |
---|
| 23 | |
---|
[5827] | 24 | |
---|
[5535] | 25 | #h_c = (q**2/g)**(1/3) |
---|
| 26 | #u_c = sqrt(g*h_c) |
---|
| 27 | |
---|
| 28 | q = 0.3 |
---|
| 29 | g = 9.81 |
---|
| 30 | |
---|
| 31 | b = 4.0 |
---|
| 32 | xmax = 25.0 |
---|
| 33 | z_bmax = 0.2 |
---|
| 34 | Fr_0 = u_0/sqrt(g*h_0) |
---|
| 35 | |
---|
| 36 | def stage(x): |
---|
| 37 | y,u,h = analytical_sol(x) |
---|
| 38 | return y |
---|
| 39 | |
---|
| 40 | |
---|
| 41 | def elevation(x): |
---|
| 42 | b = 4.0 |
---|
| 43 | xmax = 25.0 |
---|
| 44 | z_bmax = 0.2 |
---|
| 45 | z_b = zeros(len(x),Float) |
---|
| 46 | for i in range(len(x)): |
---|
| 47 | if (x[i] >= xmax-b) & (x[i] <= xmax+b): |
---|
| 48 | z_b[i] = z_bmax*(1.0-(x[i]-xmax)**2/b**2) |
---|
| 49 | else: |
---|
| 50 | z_b[i] = 0.0 |
---|
| 51 | return z_b |
---|
| 52 | #return 0.0 |
---|
| 53 | |
---|
| 54 | def xmomentum(x): |
---|
| 55 | return u_0*h_0 |
---|
| 56 | |
---|
| 57 | |
---|
| 58 | def bisection(func,xR,xL,H): |
---|
| 59 | while ((xR - xL) > epsilon): |
---|
| 60 | |
---|
| 61 | #Calculate midpoint of domain |
---|
| 62 | xM = xL + (xR - xL) / 2.0 |
---|
| 63 | |
---|
| 64 | #Find f(xM) |
---|
| 65 | if ((f(xL,H) * f(xM,H)) > 0): |
---|
| 66 | #Throw away left half |
---|
| 67 | xL = xM |
---|
| 68 | else: |
---|
| 69 | #Throw away right half |
---|
| 70 | xR = xM |
---|
| 71 | return xR |
---|
| 72 | |
---|
| 73 | def f(D,H): |
---|
| 74 | return D**3+D**2*(H-1.0-Fr_0**2/2.0)+Fr_0**2/2.0 |
---|
| 75 | |
---|
| 76 | def fprime(D,H): |
---|
| 77 | Fr_0 = u_0/sqrt(g*h_0) |
---|
| 78 | return 3*D**2+2*D*(H-1-Fr_0**2/2.0) |
---|
| 79 | |
---|
| 80 | def analytical_sol(x): |
---|
| 81 | y = zeros(len(x),Float) |
---|
| 82 | u = zeros(len(x),Float) |
---|
| 83 | height = zeros(len(x),Float) |
---|
| 84 | for i in range(len(x)): |
---|
| 85 | if (x[i] >= xmax-b) & (x[i] <= xmax+b): |
---|
| 86 | zb = z_bmax*(1.0-(x[i]-xmax)**2/b**2) |
---|
| 87 | else: |
---|
| 88 | zb = 0.0 |
---|
| 89 | |
---|
| 90 | H = zb/h_0 |
---|
| 91 | y1 = bisection(f,1.0,0.2/0.5,H) |
---|
| 92 | height[i] = y1*h_0 |
---|
| 93 | y[i] = height[i]+zb |
---|
| 94 | u[i] = q/(height[i]) |
---|
| 95 | |
---|
| 96 | |
---|
| 97 | return y,u,height |
---|
| 98 | |
---|
| 99 | Fr_0 = u_0/sqrt(g*h_0) |
---|
| 100 | zb=0.2 |
---|
| 101 | H = zb/h_0 |
---|
| 102 | D = bisection(f,1,0.2/0.5,H) |
---|
| 103 | |
---|
| 104 | L = 50.0 |
---|
[5827] | 105 | N = 400 |
---|
[5535] | 106 | cell_len = L/N |
---|
| 107 | points = zeros(N+1,Float) |
---|
| 108 | for j in range(N+1): |
---|
| 109 | points[j] = j*cell_len |
---|
| 110 | boundary = { (0,0): 'left',(N-1,1): 'right'} |
---|
| 111 | domain = Domain(points,boundary) |
---|
| 112 | D1 = Dirichlet_boundary([h_0,u_0*h_0]) |
---|
| 113 | D2 = Dirichlet_boundary([h_1,u_1*h_1]) |
---|
| 114 | domain.set_boundary({'left':D1,'right':D2}) |
---|
| 115 | domain.set_quantity('elevation',elevation) |
---|
| 116 | domain.set_quantity('stage',stage) |
---|
| 117 | domain.set_quantity('xmomentum',xmomentum) |
---|
| 118 | X = domain.vertices |
---|
| 119 | C = domain.centroids |
---|
| 120 | |
---|
| 121 | Stage = domain.quantities['stage'] |
---|
| 122 | Xmom = domain.quantities['xmomentum'] |
---|
| 123 | |
---|
| 124 | StageQ = Stage.vertex_values |
---|
| 125 | StageC = Stage.centroid_values |
---|
| 126 | XmomQ = Xmom.vertex_values |
---|
| 127 | XmomC = Xmom.centroid_values |
---|
| 128 | ElevationQ = domain.quantities['elevation'].vertex_values |
---|
| 129 | |
---|
| 130 | stage_bdry = Stage.boundary_values |
---|
| 131 | xmom_bdry = Xmom.boundary_values |
---|
| 132 | |
---|
| 133 | w,u,h = analytical_sol(X.flat) |
---|
| 134 | wc, uc, hc = analytical_sol(C) |
---|
| 135 | |
---|
[5827] | 136 | #from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ylim,xlim, rc |
---|
| 137 | #hold(False) |
---|
| 138 | #rc('text', usetex=True) |
---|
[5535] | 139 | h_error = zeros(1,Float) |
---|
| 140 | uh_error = zeros(1,Float) |
---|
| 141 | k = 0 |
---|
[5827] | 142 | yieldstep = 1.0 |
---|
| 143 | finaltime = 1.0 |
---|
[5535] | 144 | domain.limiter = "vanleer" |
---|
[5827] | 145 | domain.order = 2 |
---|
| 146 | domain.set_timestepping_method('rk2') |
---|
| 147 | domain.cfl = 1.0 |
---|
| 148 | print 'THE DOMAIN LIMITER is', domain.limiter |
---|
| 149 | import time |
---|
| 150 | t0=time.time() |
---|
[5535] | 151 | for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): |
---|
| 152 | domain.write_time() |
---|
[5827] | 153 | print "integral", domain.quantities['stage'].get_integral() |
---|
| 154 | if t>0.0: |
---|
| 155 | print 'Test 1' |
---|
| 156 | from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ion |
---|
| 157 | #ion() |
---|
| 158 | hold(False) |
---|
| 159 | #rc('text', usetex=True) |
---|
| 160 | clf() |
---|
| 161 | |
---|
| 162 | plot1 = subplot(211) |
---|
| 163 | print 'Test 2' |
---|
| 164 | plot(X,w,X,StageQ,X,ElevationQ) |
---|
| 165 | print 'Test 3' |
---|
| 166 | xlabel('Position') |
---|
| 167 | ylabel('Stage (m)') |
---|
| 168 | legend(('Analytical Solution', 'Numerical Solution', 'Channel Bed'), |
---|
| 169 | 'upper right', shadow=True) |
---|
| 170 | plot1.set_ylim([-0.1,1.0]) |
---|
| 171 | print 'Test 4' |
---|
| 172 | plot2 = subplot(212) |
---|
| 173 | print 'Test 5' |
---|
| 174 | plot(X,u*h,X,XmomQ)#/(StageQ-ElevationQ)) |
---|
| 175 | print 'Test 6' |
---|
| 176 | xlabel('x (m)') |
---|
| 177 | ylabel(r'X-momentum (m^2/s)') |
---|
| 178 | plot2.set_ylim([0.297,0.301]) |
---|
| 179 | h_error[k] = 1.0/(N)*sum(abs(wc-StageC)) |
---|
| 180 | uh_error[k] = 1.0/(N)*sum(abs(uc*hc-XmomC)) |
---|
| 181 | print "h_error %.10f" %(h_error[k]) |
---|
| 182 | print "uh_error %.10f"% (uh_error[k]) |
---|
| 183 | print "h_max %.10f"%max(wc-StageC) |
---|
| 184 | print "uh max %.10f"%max(uc*hc-XmomC) |
---|
| 185 | #filename = "steady_flow_" |
---|
| 186 | #filename += domain.limiter |
---|
| 187 | #filename += str(400) |
---|
| 188 | #filename += ".ps" |
---|
| 189 | #savefig(filename) |
---|
| 190 | show() |
---|
| 191 | print 'That took %.2f seconds'%(time.time()-t0) |
---|
| 192 | raw_input("Press return key!") |
---|