source: trunk/anuga_work/development/sudi/sw_1d/numerical_entropy_production/dam_h_2sides_windowsRUN.py @ 8041

Last change on this file since 8041 was 8041, checked in by mungkasi, 13 years ago

Adding numerical_entropy_production folder containing codes on some smoothness indicators of solutions.

File size: 6.1 KB
Line 
1import os
2from math import sqrt
3from sww_domain_shv import *
4from numpy import zeros
5#from analytic_dam_sudi import AnalyticDam
6#Note:apply analytical_sol given in debris avalanche solution
7from parameters import *
8
9N = int(N) # number of cells
10print "number of cells=",N
11#analytical_sol=AnalyticDam(h_0,h_1)
12boundary = {(0,0):'left', (N-1,1): 'right'}
13domain = Domain(points,boundary)
14domain.order = 2
15domain.set_timestepping_method('rk2')
16domain.cfl = 1.0
17domain.limiter = "minmod"
18
19#print "Check 1, bed_slope=", bed_slope
20
21def stage(x):
22    y=zeros(len(x), float)
23    for i in range (len(x)):
24        if x[i]<=L/4.0:
25            y[i]=bed_slope*x[i]#0.0
26        elif x[i]<=3*L/4.0:
27            y[i]=h_1
28        else:
29            y[i]=h_0
30    return y
31domain.set_quantity('stage',stage)
32
33def elevation(x):
34    y=zeros(len(x), float)
35    for i in range (len(x)):
36        y[i] = bed_slope*x[i]
37    return y
38domain.set_quantity('elevation',elevation)
39
40### ================ Define the boundary function =========================
41#def f_right(t):
42#    z_r = bed_slope*(0.5*L)
43#    h_r = h_0 #+ bed_slope*cell_len
44#    w_r = z_r + h_r
45#    u_r = m*t   
46#    #['stage', 'xmomentum', 'elevation', 'height', 'velocity']
47#    return [w_r,  u_r*h_r,  z_r,  h_r,  u_r]
48
49#T_right = Time_boundary(domain,f_right)
50#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])
51#D_left = Dirichlet_boundary([-1.0*bed_slope*(0.5*L),  0.0,  -1.0*bed_slope*(0.5*L),  0.0,  0.0])
52D_right = Dirichlet_boundary([h_0,  0.0,  0.0,  h_0,  0.0])
53D_left = Dirichlet_boundary([0.0,  0.0,  0.0,  0.0,  0.0])
54domain.set_boundary({'left':D_left,'right':D_right})
55### ================  End of the definition of boundary function ===========
56
57import time
58yieldstep=finaltime=10.0
59t0=time.time()
60i=1
61
62X = domain.vertices.flat
63C = domain.centroids
64
65while finaltime < 10.01:
66    for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
67        domain.write_time()
68    """
69    #finaltime = finaltime + 10.0
70    #if t>0.0:
71    N = float(N)
72    StageC = domain.quantities['stage'].centroid_values
73    XmomC = domain.quantities['xmomentum'].centroid_values
74    VelC = domain.quantities['velocity'].centroid_values
75    #hC, uhC, uC = analytical_sol(C,domain.time)
76    #h_error = sum(abs(hC-StageC))/N
77    #uh_error = sum(abs(uhC-XmomC))/N
78    #u_error = sum(abs(uC-VelC))/N
79    #print "h_error %.10f" %(h_error)
80    #print "uh_error %.10f"%(uh_error)
81    #print "u_error %.10f" %(u_error)
82    #print 'That took %.2f seconds' %(time.time()-t0)
83    X = domain.vertices.flat
84    C = domain.centroids
85    EP = domain.entropy_production
86    StageQ = domain.quantities['stage'].vertex_values.flat
87    XmomQ = domain.quantities['xmomentum'].vertex_values.flat
88    VelQ = domain.quantities['velocity'].vertex_values.flat
89    BedQ = domain.quantities['elevation'].vertex_values.flat
90    #h, uh, u = analytical_sol(X.flat, domain.time)
91
92    CK = domain.local_truncation_error_CK
93    KKP = domain.local_truncation_error_KKP
94
95    #print "Check 2, bed_slope=", bed_slope   
96
97    from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
98    hold(False)
99    plot1 = subplot(411)
100    plot(X,StageQ, X,BedQ) #plot(X,h,'k-',X,StageQ,'k--')
101    plot1.set_ylim([-1,11])
102    #plot1.set_xlim([0.0,2000.0])
103    #legend(('Numerical solution', 'Bed elevation'),
104    #        'upper left', shadow=False)   
105    #xlabel('Position')
106    ylabel('Stage')
107   
108
109    plot2 = subplot(412)
110    plot(points,CK)
111    #plot2.set_xlim([0.0,2000.0])
112    #xlabel('Position')
113    ylabel('CK')
114
115    plot3 = subplot(413)
116    plot(C,KKP)
117    #plot3.set_xlim([0.0,2000.0])
118    #xlabel('Position')
119    ylabel('KKP')
120
121    plot4 = subplot(414)
122    plot(C,EP)
123    #plot4.set_xlim([0.0,2000.0])
124    xlabel('Position')
125    ylabel('NEP')   
126   
127    #show()
128   
129    filename = "%s%03i%s" %("2dam_", i, ".png")
130    savefig(filename)
131    finaltime = finaltime + 0.25
132    #print "finaltime=", finaltime
133    i = i + 1
134    #print "The domain.limiter is", domain.limiter
135    #print 'That took %.2f seconds'%(time.time()-t0)
136    #print '============================================================================='
137    """
138    finaltime = finaltime + 100.0#0.5
139    StageC = domain.quantities['stage'].centroid_values
140    XmomC = domain.quantities['xmomentum'].centroid_values
141    VelC = domain.quantities['velocity'].centroid_values
142    HeiC = domain.quantities['height'].centroid_values
143    Energy = zeros(N,float)
144    for r in range(N):
145        Energy[r] = 0.5*VelC[r]*VelC[r] + g*StageC[r]   
146
147    StageQ = domain.quantities['stage'].vertex_values.flat
148    XmomQ = domain.quantities['xmomentum'].vertex_values.flat
149    VelQ = domain.quantities['velocity'].vertex_values.flat
150    BedQ = domain.quantities['elevation'].vertex_values.flat
151
152    CK = domain.local_truncation_error_CK
153    KKP = domain.local_truncation_error_KKP
154    EP = domain.entropy_production   
155    ECV = domain.entropy_centroid_values
156    from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
157    hold(False)
158   
159    plot1 = subplot(511)
160    plot(X,StageQ, X,BedQ)
161    #plot1.set_ylim([-1,11])
162    #plot1.set_xlim([0.0,2000.0])
163    #legend(('Numerical solution', 'Bed elevation'),
164    #        'center left', shadow=False)   
165    #xlabel('Position')
166    ylabel('Sta')
167   
168
169    plot2 = subplot(512)
170    plot(X,XmomQ)
171    #plot2.set_xlim([0.0,2000.0])
172    #xlabel('Position')
173    ylabel('Mom')
174
175    plot3 = subplot(513)
176    plot(C,ECV)
177    #plot3.set_xlim([0.0,2000.0])
178    #xlabel('Position')
179    ylabel('Ent')
180
181    plot4 = subplot(514)
182    plot(C,EP)
183    #plot4.set_xlim([0.0,2000.0])
184    xlabel('Position')
185    ylabel('NEP')     
186   
187    plot5 = subplot(515)
188    plot(C,Energy)
189    #plot4.set_xlim([0.0,2000.0])
190    xlabel('Position')
191    ylabel('Ene') 
192   
193    print 'That took %.2f seconds'%(time.time()-t0)
194    show()
195       
196    filename = "%s%03i%s" %("2steadyflow_", i, ".png")
197    #filename = "%s" %("steady_flow.png")
198    #savefig(filename)
199    #finaltime = finaltime + 0.25
200    #print "finaltime=", finaltime
201    i = i + 1
202    #print "The domain.limiter is", domain.limiter
203    #print 'That took %.2f seconds'%(time.time()-t0)
Note: See TracBrowser for help on using the repository browser.