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

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

Adding sudi's code

File size: 3.7 KB
Line 
1import os
2from math import sqrt
3from sww_domain_shm import *
4from Numeric import Float
5from numpy import zeros
6#from analytic_dam_sudi import AnalyticDam
7#Note:apply analytical_sol given in debris avalanche solution
8from parameters import *
9
10N = int(N) # number of cells
11print "number of cells=",N
12#analytical_sol=AnalyticDam(h_0,h_1)
13boundary = {(0,0):'left', (N-1,1): 'right'}
14domain = Domain(points,boundary)
15domain.order = 2
16domain.set_timestepping_method('rk2')
17domain.cfl = 1.0
18domain.limiter = "minmod"
19
20def stage(x):
21    y=zeros(len(x), Float)
22    for i in range (len(x)):
23        if x[i]<=L/4.0:
24            y[i]=bed_slope*x[i]#0.0
25        elif x[i]<=3*L/4.0:
26            y[i]=h_1
27        else:
28            y[i]=h_0
29    return y
30domain.set_quantity('stage',stage)
31
32def elevation(x):
33    y=zeros(len(x), Float)
34    for i in range (len(x)):
35        y[i] = bed_slope*x[i]
36    return y
37domain.set_quantity('elevation',elevation)
38
39### ================ Define the boundary function =========================
40#def f_right(t):
41#    z_r = bed_slope*(0.5*L)
42#    h_r = h_0 #+ bed_slope*cell_len
43#    w_r = z_r + h_r
44#    u_r = m*t   
45#    #['stage', 'xmomentum', 'elevation', 'height', 'velocity']
46#    return [w_r,  u_r*h_r,  z_r,  h_r,  u_r]
47
48#T_right = Time_boundary(domain,f_right)
49#D_right = Dirichlet_boundary([bed_slope*(0.5*L)+h_0,(m*domain.time)*h_0,bed_slope*(0.5*L),h_0,m*domain.time])
50#D_left = Dirichlet_boundary([-1.0*bed_slope*(0.5*L),  0.0,  -1.0*bed_slope*(0.5*L),  0.0,  0.0])
51D_right = Dirichlet_boundary([h_0,  0.0,  0.0,  h_0,  0.0])
52D_left = Dirichlet_boundary([0.0,  0.0,  0.0,  0.0,  0.0])
53domain.set_boundary({'left':D_left,'right':D_right})
54### ================  End of the definition of boundary function ===========
55
56X=domain.vertices
57C=domain.centroids
58import time
59yieldstep=finaltime=1.0
60t0=time.time()
61i=1
62
63while finaltime < 1.01:
64    for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
65        domain.write_time()
66    finaltime = finaltime + 10.0 
67if t>0.0:
68    N = float(N)
69    StageC = domain.quantities['stage'].centroid_values
70    XmomC = domain.quantities['xmomentum'].centroid_values
71    VelC = domain.quantities['velocity'].centroid_values
72    #hC, uhC, uC = analytical_sol(C,domain.time)
73    #h_error = sum(abs(hC-StageC))/N
74    #uh_error = sum(abs(uhC-XmomC))/N
75    #u_error = sum(abs(uC-VelC))/N
76    #print "h_error %.10f" %(h_error)
77    #print "uh_error %.10f"%(uh_error)
78    #print "u_error %.10f" %(u_error)
79    #print 'That took %.2f seconds' %(time.time()-t0)
80    X = domain.vertices
81    StageQ = domain.quantities['stage'].vertex_values
82    XmomQ = domain.quantities['xmomentum'].vertex_values
83    VelQ = domain.quantities['velocity'].vertex_values
84    BedQ = domain.quantities['elevation'].vertex_values
85    #h, uh, u = analytical_sol(X.flat, domain.time)
86
87    SD = domain.shock_detector
88
89    from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
90    hold(False)
91    plot1 = subplot(211)
92    #plot(X,h, X,StageQ, X,BedQ) #plot(X,h,'k-',X,StageQ,'k--')
93    plot(X,StageQ, X,BedQ)
94    plot1.set_ylim([-1,11])
95    plot1.set_xlim([0.0,2000.0])
96    legend(('Numerical solution', 'Bed elevation'),
97            'upper left', shadow=False)   
98    #xlabel('Position')
99    ylabel('Stage')
100   
101
102    plot2 = subplot(212)
103    plot(points,SD)
104    plot2.set_xlim([0.0,2000.0])
105    xlabel('Position')
106    ylabel('Smoothness indicator')
107
108    show()
109   
110    #filename = "%s%04i%s" %("dam_", i, ".eps")
111    #savefig(filename)
112    #finaltime = finaltime + 0.25
113    #print "finaltime=", finaltime
114    #i = i + 1
115    #print "The domain.limiter is", domain.limiter
116    #print 'That took %.2f seconds'%(time.time()-t0)
117    #print '============================================================================='
Note: See TracBrowser for help on using the repository browser.