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