source: trunk/anuga_work/development/sudi/sw_1d/shock_detector_shv/sf_RUN.py @ 7906

Last change on this file since 7906 was 7906, checked in by steve, 14 years ago

Adding sudi's code

File size: 3.0 KB
Line 
1import os
2from math import sqrt
3from sww_domain_shv import *
4from Numeric import Float
5from numpy import zeros
6from sf_parameters import *
7
8print "PRINT IT OUT!"
9N = int(N) # number of cells
10print "number of cells=",N
11boundary = {(0,0):'left', (N-1,1): 'right'}
12domain = Domain(points,boundary)
13domain.order = 2
14domain.set_timestepping_method('rk2')
15domain.cfl = 1.0
16domain.limiter = "minmod"
17
18def 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
26domain.set_quantity('stage',stage)
27
28def 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
36domain.set_quantity('elevation',elevation)
37
38### ================ Define the boundary function =========================
39# ['stage', 'xmomentum', 'elevation', 'height', 'velocity']
40D_left = Dirichlet_boundary([0.4125,  0.18,  0.0,  0.4125,  0.18/0.4125])
41D_right = Dirichlet_boundary([0.33,  0.18,  0.0,  0.33,  0.18/0.33])
42domain.set_boundary({'left':D_left,'right':D_right})
43### ================  End of the definition of boundary function ===========
44
45X=domain.vertices
46C=domain.centroids
47import time
48yieldstep=finaltime=10.0
49t0=time.time()
50i=1
51
52while 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
Note: See TracBrowser for help on using the repository browser.