[7837] | 1 | from scipy import sin, cos, sqrt, linspace, pi, dot |
---|
| 2 | from Numeric import zeros, Float, array |
---|
| 3 | from gaussPivot import * |
---|
| 4 | from analytical_prescription import * |
---|
| 5 | from parameter import * |
---|
| 6 | import os, time, csv, pprint |
---|
| 7 | from domain_johns import * |
---|
| 8 | from config import g, epsilon |
---|
| 9 | from rootsearch import * |
---|
| 10 | from bisect import * |
---|
| 11 | |
---|
| 12 | #Analytical computations################################################################# |
---|
| 13 | def root_g(a,b,t): |
---|
| 14 | dx = 0.01 |
---|
| 15 | def g(u): |
---|
| 16 | return u + 2.0*A*pi/T*sin(2.0*pi/T*(t+u)) |
---|
| 17 | while 1: |
---|
| 18 | x1,x2 = rootsearch(g,a,b,dx) |
---|
| 19 | if x1 != None: |
---|
| 20 | a = x2 |
---|
| 21 | root = bisect(g,x1,x2,1) |
---|
| 22 | else: |
---|
| 23 | break |
---|
| 24 | return root |
---|
| 25 | def shore(t): |
---|
| 26 | a = -1.0#-0.2#-1.0 |
---|
| 27 | b = 1.0#0.2#1.0 |
---|
| 28 | #dx = 0.01 |
---|
| 29 | u = root_g(a,b,t) |
---|
| 30 | xi = -0.5*u*u + A*cos(2.0*pi/T*(t+u)) |
---|
| 31 | position = 1.0 + xi |
---|
| 32 | return position, u |
---|
| 33 | |
---|
| 34 | |
---|
| 35 | |
---|
| 36 | |
---|
| 37 | #Numerical computations################################################################### |
---|
| 38 | def newtonRaphson2(f,q,tol=1.0e-15): ##1.0e-9 may be too large |
---|
| 39 | for i in range(30): |
---|
| 40 | h = 1.0e-4 ##1.0e-4 may be too large. |
---|
| 41 | n = len(q) |
---|
| 42 | jac = zeros((n,n),Float) |
---|
| 43 | if 1.0+q[0]-x<0.0: |
---|
| 44 | temp1 = 1.0+q[0]-x |
---|
| 45 | q[0] = q[0]-temp1 |
---|
| 46 | q[1] = v |
---|
| 47 | return q |
---|
| 48 | f0 = f(q) |
---|
| 49 | for i in range(n): |
---|
| 50 | temp = q[i] |
---|
| 51 | q[i] = temp + h |
---|
| 52 | f1 = f(q) |
---|
| 53 | q[i] = temp |
---|
| 54 | jac[:,i] = (f1 - f0)/h |
---|
| 55 | if sqrt(dot(f0,f0)/len(q)) < tol: return q |
---|
| 56 | dq = gaussPivot(jac,-f0) |
---|
| 57 | q = q + dq |
---|
| 58 | if sqrt(dot(dq,dq)) < tol*max(max(abs(q)),1.0): return q |
---|
| 59 | print 'Too many iterations' |
---|
| 60 | |
---|
| 61 | def elevation(X): |
---|
| 62 | N = len(X) |
---|
| 63 | z = zeros(N,Float) |
---|
| 64 | for i in range(N): |
---|
| 65 | z[i] = (h_0/L)*X[i] - h_0 |
---|
| 66 | return z |
---|
| 67 | |
---|
| 68 | def height(X): |
---|
| 69 | N = len(X) |
---|
| 70 | z = zeros(N,Float) |
---|
| 71 | for i in range(N): |
---|
| 72 | z[i] = h_0 - (h_0/L)*X[i] |
---|
| 73 | return z |
---|
| 74 | |
---|
| 75 | def velocity(X): |
---|
| 76 | N = len(X) |
---|
| 77 | return zeros(N,Float) |
---|
| 78 | |
---|
| 79 | boundary = { (0,0): 'left',(N-1,1): 'right'} |
---|
| 80 | |
---|
| 81 | domain = Domain(points,boundary) |
---|
| 82 | domain.order = 2 |
---|
| 83 | domain.set_timestepping_method('rk2') |
---|
| 84 | domain.set_CFL(1.0) |
---|
| 85 | domain.beta = 1.0 |
---|
| 86 | domain.set_limiter("minmod") |
---|
| 87 | |
---|
| 88 | |
---|
| 89 | def f_CG(t): |
---|
| 90 | timing = t*sqrt(g*h_0)/L |
---|
| 91 | w, u = prescribe(0.0,timing) |
---|
| 92 | wO = w*h_0 |
---|
| 93 | uO = u*sqrt(g*h_0) |
---|
| 94 | zO = -h_0 |
---|
| 95 | hO = wO - zO |
---|
| 96 | pO = uO * hO |
---|
| 97 | #[ 'stage', 'xmomentum', 'elevation', 'height', 'velocity'] |
---|
| 98 | return [wO, pO, zO, hO, uO] |
---|
| 99 | def f_JOHNS(t): |
---|
| 100 | timing = t*sqrt(g*h_0)/L |
---|
| 101 | w, u = prescribe_at_O_JOHNS(timing) |
---|
| 102 | wO = w*h_0 |
---|
| 103 | uO = u*sqrt(g*h_0) |
---|
| 104 | zO = -h_0 |
---|
| 105 | hO = wO - zO |
---|
| 106 | pO = uO * hO |
---|
| 107 | #[ 'stage', 'xmomentum', 'elevation', 'height', 'velocity'] |
---|
| 108 | return [wO, pO, zO, hO, uO] |
---|
| 109 | |
---|
| 110 | T1 = Time_boundary(domain,f_JOHNS) |
---|
| 111 | D2 = Dirichlet_boundary([50.5, 0.0, 50.5, 0.0, 0.0]) |
---|
| 112 | domain.set_boundary({'left':T1,'right':D2}) |
---|
| 113 | |
---|
| 114 | domain.set_quantity('height',height) |
---|
| 115 | domain.set_quantity('elevation',elevation) |
---|
| 116 | domain.set_quantity('velocity',velocity) |
---|
| 117 | |
---|
| 118 | |
---|
| 119 | Ver = domain.vertices |
---|
| 120 | n_V = len(Ver) |
---|
| 121 | AnalitW_V = zeros((n_V,2), Float) |
---|
| 122 | AnalitP_V = zeros((n_V,2), Float) |
---|
| 123 | AnalitZ_V = zeros((n_V,2), Float) |
---|
| 124 | AnalitH_V = zeros((n_V,2), Float) |
---|
| 125 | AnalitU_V = zeros((n_V,2), Float) |
---|
| 126 | |
---|
| 127 | Cen = domain.centroids |
---|
| 128 | n_C = len(Cen) |
---|
| 129 | AnalitW_C = zeros(n_C, Float) |
---|
| 130 | AnalitP_C = zeros(n_C, Float) |
---|
| 131 | AnalitZ_C = zeros(n_C, Float) |
---|
| 132 | AnalitH_C = zeros(n_C, Float) |
---|
| 133 | AnalitU_C = zeros(n_C, Float) |
---|
| 134 | |
---|
| 135 | waktu = 10.0 #3.0*60.0 |
---|
| 136 | WAKTU = 12690.0 #Note: Tp=15.0*60.0 |
---|
| 137 | yieldstep = finaltime = waktu |
---|
| 138 | t0 = time.time() |
---|
| 139 | counter=1 |
---|
| 140 | |
---|
| 141 | shorelines_numerical_johns = zeros(int(WAKTU/waktu), Float) |
---|
| 142 | shorelines_analytical = zeros(int(WAKTU/waktu), Float) |
---|
| 143 | time_instants = zeros(int(WAKTU/waktu), Float) |
---|
| 144 | print "the initial time_instants=", time_instants |
---|
| 145 | |
---|
| 146 | |
---|
| 147 | while finaltime < WAKTU+0.1: |
---|
| 148 | for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): |
---|
| 149 | domain.write_time() |
---|
| 150 | time_instants[counter-1] = domain.time |
---|
| 151 | |
---|
| 152 | Stage = domain.quantities['stage'] |
---|
| 153 | Momentum = domain.quantities['xmomentum'] |
---|
| 154 | Elevation = domain.quantities['elevation'] |
---|
| 155 | Height = domain.quantities['height'] |
---|
| 156 | Velocity = domain.quantities['velocity'] |
---|
| 157 | |
---|
| 158 | StageV = Stage.vertex_values |
---|
| 159 | MomV = Momentum.vertex_values |
---|
| 160 | ElevationV = Elevation.vertex_values |
---|
| 161 | HeightV = Height.vertex_values |
---|
| 162 | VelV = Velocity.vertex_values |
---|
| 163 | |
---|
| 164 | StageC = Stage.centroid_values |
---|
| 165 | MomC = Momentum.centroid_values |
---|
| 166 | ElevationC = Elevation.centroid_values |
---|
| 167 | HeightC = Height.centroid_values |
---|
| 168 | VelC = Velocity.centroid_values |
---|
| 169 | |
---|
| 170 | table = zeros((len(Ver.flat),6),Float) |
---|
| 171 | for r in range(len(Ver.flat)): |
---|
| 172 | for c in range(6): |
---|
| 173 | if c==0: |
---|
| 174 | table[r][c] = Ver.flat[r] |
---|
| 175 | elif c==1: |
---|
| 176 | table[r][c] = StageV.flat[r] |
---|
| 177 | elif c==2: |
---|
| 178 | table[r][c] = MomV.flat[r] |
---|
| 179 | elif c==3: |
---|
| 180 | table[r][c] = ElevationV.flat[r] |
---|
| 181 | elif c==4: |
---|
| 182 | table[r][c] = HeightV.flat[r] |
---|
| 183 | else: |
---|
| 184 | table[r][c] = VelV.flat[r] |
---|
| 185 | |
---|
| 186 | outname = "%s%04i%s%f%s" %("numerical_johns_", counter, "_", domain.time, ".csv") |
---|
| 187 | outfile = open(outname, 'w') |
---|
| 188 | writer = csv.writer(outfile) |
---|
| 189 | for row in table: |
---|
| 190 | writer.writerow(row) |
---|
| 191 | outfile.close() |
---|
| 192 | |
---|
| 193 | |
---|
| 194 | for s in range(2*n_V): |
---|
| 195 | heiL = HeightV.flat[s] |
---|
| 196 | momR = MomV.flat[s+1] |
---|
| 197 | if heiL >= 1e-6: |
---|
| 198 | if abs(momR)==0.0: #<1e-15: |
---|
| 199 | break |
---|
| 200 | #print "s+1=",s+1 |
---|
| 201 | shorelines_numerical_johns[counter-1] = Ver.flat[s+1] |
---|
| 202 | #print "shorelines_numerical_johns=",shorelines_numerical_johns[counter-1] |
---|
| 203 | |
---|
| 204 | |
---|
| 205 | |
---|
| 206 | #print "Now the ANALYTIC" |
---|
| 207 | pos_shore, vel_shore = shore(domain.time*sqrt(g*h_0)/L) #dimensionless |
---|
| 208 | pos_shore = pos_shore*L |
---|
| 209 | vel_shore = vel_shore*sqrt(g*h_0) |
---|
| 210 | |
---|
| 211 | shorelines_analytical[counter-1] = pos_shore |
---|
| 212 | |
---|
| 213 | #The following is for calculating the error at centroids |
---|
| 214 | for i in range(n_C): |
---|
| 215 | x = Cen[i] |
---|
| 216 | if x < pos_shore: |
---|
| 217 | sta, vel = prescribe(x/L,domain.time*sqrt(g*h_0)/L) #dimensionless |
---|
| 218 | AnalitW_C[i] = sta*h_0 |
---|
| 219 | AnalitU_C[i] = vel #It needs dimensionalisation |
---|
| 220 | AnalitZ_C[i] = (h_0/L)*x - h_0 |
---|
| 221 | AnalitH_C[i] = AnalitW_C[i] - AnalitZ_C[i] |
---|
| 222 | AnalitP_C[i] = AnalitH_C[i]*AnalitU_C[i] #It needs dimensionalisation |
---|
| 223 | else: |
---|
| 224 | AnalitW_C[i] = (h_0/L)*x - h_0 |
---|
| 225 | AnalitU_C[i] = 0.0 |
---|
| 226 | AnalitZ_C[i] = (h_0/L)*x - h_0 |
---|
| 227 | AnalitH_C[i] = 0.0 |
---|
| 228 | AnalitP_C[i] = 0.0 |
---|
| 229 | AnalitU_C = AnalitU_C*sqrt(g*h_0) #This is the dimensionalisation |
---|
| 230 | AnalitP_C = AnalitP_C*sqrt(g*h_0) #This is the dimensionalisation |
---|
| 231 | |
---|
| 232 | error_W = (1.0/n_C)*sum(abs(StageC-AnalitW_C)) |
---|
| 233 | error_P = (1.0/n_C)*sum(abs(MomC-AnalitP_C)) |
---|
| 234 | error_U = (1.0/n_C)*sum(abs(VelC-AnalitU_C)) |
---|
| 235 | table = zeros((1,4),Float) |
---|
| 236 | for r in range(1): |
---|
| 237 | for c in range(4): |
---|
| 238 | if c==0: |
---|
| 239 | table[r][c] = domain.time |
---|
| 240 | elif c==1: |
---|
| 241 | table[r][c] = error_W |
---|
| 242 | elif c==2: |
---|
| 243 | table[r][c] = error_P |
---|
| 244 | else: |
---|
| 245 | table[r][c] = error_U |
---|
| 246 | outname = "%s%04i%s%f%s" %("error_johns_", counter, "_", domain.time, ".csv") |
---|
| 247 | outfile = open(outname, 'w') |
---|
| 248 | writer = csv.writer(outfile) |
---|
| 249 | for row in table: |
---|
| 250 | writer.writerow(row) |
---|
| 251 | outfile.close() |
---|
| 252 | |
---|
| 253 | |
---|
| 254 | #The following is for ploting the quantities at vertex values |
---|
| 255 | for i in range(n_V): |
---|
| 256 | vector_x = Ver[i] |
---|
| 257 | for k in range(2): |
---|
| 258 | x = vector_x[k] |
---|
| 259 | if x < pos_shore: |
---|
| 260 | sta, vel = prescribe(x/L,domain.time*sqrt(g*h_0)/L) #dimensionless |
---|
| 261 | AnalitW_V[i,k] = sta*h_0 |
---|
| 262 | AnalitU_V[i,k] = vel #It needs dimensionalisation |
---|
| 263 | AnalitZ_V[i,k] = (h_0/L)*x - h_0 |
---|
| 264 | AnalitH_V[i,k] = AnalitW_V[i,k] - AnalitZ_V[i,k] |
---|
| 265 | AnalitP_V[i,k] = AnalitH_V[i,k]*AnalitU_V[i,k] #It needs dimensionalisation |
---|
| 266 | else: |
---|
| 267 | AnalitW_V[i,k] = (h_0/L)*x - h_0 |
---|
| 268 | AnalitU_V[i,k] = 0.0 |
---|
| 269 | AnalitZ_V[i,k] = (h_0/L)*x - h_0 |
---|
| 270 | AnalitH_V[i,k] = 0.0 |
---|
| 271 | AnalitP_V[i,k] = 0.0 |
---|
| 272 | AnalitU_V = AnalitU_V*sqrt(g*h_0) #This is the dimensionalisation |
---|
| 273 | AnalitP_V = AnalitP_V*sqrt(g*h_0) #This is the dimensionalisation |
---|
| 274 | |
---|
| 275 | |
---|
| 276 | table = zeros((len(Ver.flat),6),Float) |
---|
| 277 | for r in range(len(Ver.flat)): |
---|
| 278 | for c in range(6): |
---|
| 279 | if c==0: |
---|
| 280 | table[r][c] = Ver.flat[r] |
---|
| 281 | elif c==1: |
---|
| 282 | table[r][c] = AnalitW_V.flat[r] |
---|
| 283 | elif c==2: |
---|
| 284 | table[r][c] = AnalitP_V.flat[r] |
---|
| 285 | elif c==3: |
---|
| 286 | table[r][c] = AnalitZ_V.flat[r] |
---|
| 287 | elif c==4: |
---|
| 288 | table[r][c] = AnalitH_V.flat[r] |
---|
| 289 | else: |
---|
| 290 | table[r][c] = AnalitU_V.flat[r] |
---|
| 291 | |
---|
| 292 | outname = "%s%04i%s%f%s" %("analytical_", counter, "_", domain.time, ".csv") |
---|
| 293 | outfile = open(outname, 'w') |
---|
| 294 | writer = csv.writer(outfile) |
---|
| 295 | for row in table: |
---|
| 296 | writer.writerow(row) |
---|
| 297 | outfile.close() |
---|
| 298 | |
---|
| 299 | #put this |
---|
| 300 | from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ion |
---|
| 301 | hold(False) |
---|
| 302 | clf() |
---|
| 303 | |
---|
| 304 | plot1 = subplot(311) |
---|
| 305 | plot(Ver/1e+4,AnalitW_V,'b-', Ver/1e+4,StageV,'g-', Ver/1e+4,ElevationV,'k-') |
---|
| 306 | #plot(Ver,StageV, Ver,ElevationV) |
---|
| 307 | #plot(Ver/L,StageV/h_0, Ver/L,ElevationV/h_0) |
---|
| 308 | #xlabel('Position') |
---|
| 309 | ylabel('Stage') |
---|
| 310 | #plot1.set_xlim([0.0,1.2]) |
---|
| 311 | plot1.set_ylim([-6.0,6.0])#([-9.0e-3,9.0e-3]) |
---|
| 312 | #legend(('Analytical Solution', 'Numerical Solution', 'Discretized Bed'), |
---|
| 313 | # 'upper right', shadow=False) |
---|
| 314 | |
---|
| 315 | plot2 = subplot(312) |
---|
| 316 | plot(Ver/1e+4,AnalitP_V,'b-', Ver/1e+4,MomV,'g-') |
---|
| 317 | #plot(Ver/L, VelV/sqrt(g*h_0)) |
---|
| 318 | #xlabel('Position') |
---|
| 319 | ylabel('Momentum') |
---|
| 320 | #plot2.set_xlim([0.0,1.2]) |
---|
| 321 | #legend(('Analytical Solution','Numerical Solution'), |
---|
| 322 | # 'upper right', shadow=False) |
---|
| 323 | |
---|
| 324 | plot3 = subplot(313) |
---|
| 325 | plot(Ver/1e+4,AnalitU_V,'b-', Ver/1e+4,VelV,'g-') |
---|
| 326 | #plot(Ver/L, VelV/sqrt(g*h_0)) |
---|
| 327 | xlabel('Position / 10,000') |
---|
| 328 | ylabel('Velocity') |
---|
| 329 | #plot2.set_xlim([0.0,1.2]) |
---|
| 330 | legend(('Analytical Solution','Numerical Solution'), |
---|
| 331 | 'upper center', shadow=False) |
---|
| 332 | |
---|
| 333 | filename = "%s%04i%s%f%s" %("numerical_johns_", counter,"_", domain.time, ".png") |
---|
| 334 | savefig(filename) |
---|
| 335 | #show() |
---|
| 336 | #raw_input("Press ENTER to continue") |
---|
| 337 | #put this |
---|
| 338 | counter = counter+1 |
---|
| 339 | finaltime = finaltime + waktu |
---|
| 340 | |
---|
| 341 | |
---|
| 342 | |
---|
| 343 | table = zeros((int(WAKTU/waktu), 2),Float) |
---|
| 344 | for r in range(int(WAKTU/waktu)): |
---|
| 345 | for c in range(2): |
---|
| 346 | if c==0: |
---|
| 347 | table[r][c] = time_instants[r] |
---|
| 348 | else: |
---|
| 349 | table[r][c] = shorelines_numerical_johns[r] |
---|
| 350 | |
---|
| 351 | outname = "%s" %("shore_numerical_johns.csv") |
---|
| 352 | outfile = open(outname, 'w') |
---|
| 353 | writer = csv.writer(outfile) |
---|
| 354 | for row in table: |
---|
| 355 | writer.writerow(row) |
---|
| 356 | outfile.close() |
---|
| 357 | |
---|
| 358 | |
---|
| 359 | table = zeros((int(WAKTU/waktu), 2),Float) |
---|
| 360 | for r in range(int(WAKTU/waktu)): |
---|
| 361 | for c in range(2): |
---|
| 362 | if c==0: |
---|
| 363 | table[r][c] = time_instants[r] |
---|
| 364 | else: |
---|
| 365 | table[r][c] = shorelines_analytical[r] |
---|
| 366 | |
---|
| 367 | outname = "%s" %("shore_analytical.csv") |
---|
| 368 | outfile = open(outname, 'w') |
---|
| 369 | writer = csv.writer(outfile) |
---|
| 370 | for row in table: |
---|
| 371 | writer.writerow(row) |
---|
| 372 | outfile.close() |
---|
| 373 | |
---|