source: trunk/anuga_work/development/sudi/sw_1d/shock_detector_shv/sf_interactive_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: 2.8 KB
Line 
1import os
2from math import sqrt
3from sww_domain_simplified import *
4from Numeric import Float
5from numpy import zeros
6from sf_parameters import *
7
8N = int(N) # number of cells
9print "number of cells=",N
10boundary = {(0,0):'left', (N-1,1): 'right'}
11domain = Domain(points,boundary)
12domain.order = 2
13domain.set_timestepping_method('rk2')
14domain.cfl = 1.0
15domain.limiter = "minmod"
16
17def 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.33
24    return y
25domain.set_quantity('stage',stage)
26
27def 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
35domain.set_quantity('elevation',elevation)
36
37### ================ Define the boundary function =========================
38# ['stage', 'xmomentum', 'elevation', 'height', 'velocity']
39D_left = Dirichlet_boundary([0.4125,  0.18,  0.0,  0.4125,  0.18/0.4125])
40D_right = Dirichlet_boundary([0.33,  0.18,  0.0,  0.33,  0.18/0.33])
41domain.set_boundary({'left':D_left,'right':D_right})
42### ================  End of the definition of boundary function ===========
43
44X=domain.vertices
45C=domain.centroids
46import time
47yieldstep=finaltime=0.1
48t0=time.time()
49i=1
50
51while finaltime < 0.101:
52    for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
53        domain.write_time()
54    finaltime = finaltime + 10.0 
55if 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    hold(False)
71
72    plot1 = subplot(211)
73    plot(X,StageQ, X,BedQ)
74    #plot1.set_ylim([-1,11])
75    #plot1.set_xlim([0.0,2000.0])
76    legend(('Numerical solution', 'Bed elevation'),
77            'upper left', shadow=False)   
78    #xlabel('Position')
79    ylabel('Stage')
80   
81
82    plot2 = subplot(212)
83    plot(points,SD)
84    #plot2.set_xlim([0.0,2000.0])
85    xlabel('Position')
86    ylabel('Smoothness indicator')
87   
88    print 'That took %.2f seconds'%(time.time()-t0)
89    show()
90       
91    #filename = "%s%04i%s" %("dam_", i, ".eps")
92    #savefig(filename)
93    #finaltime = finaltime + 0.25
94    #print "finaltime=", finaltime
95    #i = i + 1
96    #print "The domain.limiter is", domain.limiter
97    #print 'That took %.2f seconds'%(time.time()-t0)
98    #print '============================================================================='
Note: See TracBrowser for help on using the repository browser.