[7906] | 1 | import os |
---|
| 2 | from math import sqrt |
---|
| 3 | from sww_domain_shv import * |
---|
| 4 | from Numeric import Float |
---|
| 5 | from numpy import zeros |
---|
| 6 | from sf_parameters import * |
---|
| 7 | |
---|
| 8 | print "PRINT IT OUT!" |
---|
| 9 | N = int(N) # number of cells |
---|
| 10 | print "number of cells=",N |
---|
| 11 | boundary = {(0,0):'left', (N-1,1): 'right'} |
---|
| 12 | domain = Domain(points,boundary) |
---|
| 13 | domain.order = 2 |
---|
| 14 | domain.set_timestepping_method('rk2') |
---|
| 15 | domain.cfl = 1.0 |
---|
| 16 | domain.limiter = "minmod" |
---|
| 17 | |
---|
| 18 | def stage(x): |
---|
| 19 | y=zeros(len(x), Float) |
---|
| 20 | for i in range(len(x)): |
---|
| 21 | if x[i] < 11.666: |
---|
| 22 | y[i] = 0.4125 |
---|
| 23 | else: |
---|
| 24 | y[i] = 0.28 #0.33 |
---|
| 25 | return y |
---|
| 26 | domain.set_quantity('stage',stage) |
---|
| 27 | |
---|
| 28 | def elevation(x): |
---|
| 29 | z_b = zeros(len(x),Float) |
---|
| 30 | for i in range(len(x)): |
---|
| 31 | if (x[i] >= 8.0) & (x[i] <= 12.0): |
---|
| 32 | z_b[i] = 0.2 - 0.05*(x[i]-10.0)**2.0 |
---|
| 33 | else: |
---|
| 34 | z_b[i] = 0.0 |
---|
| 35 | return z_b |
---|
| 36 | domain.set_quantity('elevation',elevation) |
---|
| 37 | |
---|
| 38 | ### ================ Define the boundary function ========================= |
---|
| 39 | # ['stage', 'xmomentum', 'elevation', 'height', 'velocity'] |
---|
| 40 | D_left = Dirichlet_boundary([0.4125, 0.18, 0.0, 0.4125, 0.18/0.4125]) |
---|
| 41 | D_right = Dirichlet_boundary([0.33, 0.18, 0.0, 0.33, 0.18/0.33]) |
---|
| 42 | domain.set_boundary({'left':D_left,'right':D_right}) |
---|
| 43 | ### ================ End of the definition of boundary function =========== |
---|
| 44 | |
---|
| 45 | X=domain.vertices |
---|
| 46 | C=domain.centroids |
---|
| 47 | import time |
---|
| 48 | yieldstep=finaltime=10.0 |
---|
| 49 | t0=time.time() |
---|
| 50 | i=1 |
---|
| 51 | |
---|
| 52 | while finaltime < 200.1: |
---|
| 53 | for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime): |
---|
| 54 | domain.write_time() |
---|
| 55 | #if t>0.0: |
---|
| 56 | N = float(N) |
---|
| 57 | StageC = domain.quantities['stage'].centroid_values |
---|
| 58 | XmomC = domain.quantities['xmomentum'].centroid_values |
---|
| 59 | VelC = domain.quantities['velocity'].centroid_values |
---|
| 60 | |
---|
| 61 | X = domain.vertices |
---|
| 62 | StageQ = domain.quantities['stage'].vertex_values |
---|
| 63 | XmomQ = domain.quantities['xmomentum'].vertex_values |
---|
| 64 | VelQ = domain.quantities['velocity'].vertex_values |
---|
| 65 | BedQ = domain.quantities['elevation'].vertex_values |
---|
| 66 | |
---|
| 67 | SD = domain.shock_detector |
---|
| 68 | |
---|
| 69 | #from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot |
---|
| 70 | import matplotlib |
---|
| 71 | matplotlib.use('Agg') |
---|
| 72 | import matplotlib.pyplot as plt |
---|
| 73 | from matplotlib.pylab import hold |
---|
| 74 | |
---|
| 75 | hold(False) |
---|
| 76 | plt.figure(1) |
---|
| 77 | |
---|
| 78 | plt.subplot(311) |
---|
| 79 | plt.plot(X,StageQ, X,BedQ) |
---|
| 80 | #plot1.set_ylim([-1,11]) |
---|
| 81 | #plot1.set_xlim([0.0,2000.0]) |
---|
| 82 | #plt.legend(('Numerical solution', 'Bed elevation'), |
---|
| 83 | # 'upper left', shadow=False) |
---|
| 84 | #xlabel('Position') |
---|
| 85 | plt.ylabel('Stage') |
---|
| 86 | |
---|
| 87 | plt.subplot(312) |
---|
| 88 | plt.plot(points,SD) |
---|
| 89 | #plot2.set_xlim([0.0,2000.0]) |
---|
| 90 | #plt.xlabel('Position') |
---|
| 91 | plt.ylabel('Smoothness indicator') |
---|
| 92 | |
---|
| 93 | plt.subplot(313) |
---|
| 94 | plt.plot(X,XmomQ) |
---|
| 95 | #plot2.set_xlim([0.0,2000.0]) |
---|
| 96 | plt.xlabel('Position') |
---|
| 97 | plt.ylabel('Discharge') |
---|
| 98 | |
---|
| 99 | #print 'That took %.2f seconds'%(time.time()-t0) |
---|
| 100 | print "domain.time=", domain.time |
---|
| 101 | auxtime = int(round(domain.time)) |
---|
| 102 | print "auxtime=", auxtime |
---|
| 103 | filename = "%s%02i%s%i" %("shv_nonuniformtime", i, "_", auxtime) |
---|
| 104 | plt.savefig(filename) |
---|
| 105 | finaltime = finaltime + 10.0 |
---|
| 106 | i = i + 1 |
---|
| 107 | |
---|