source: trunk/anuga_work/development/sudi/sw_1d/shock_detector/shm/sf_RUN.py @ 7909

Last change on this file since 7909 was 7909, checked in by mungkasi, 14 years ago

Add 1d codes.

File size: 3.0 KB
Line 
1import os
2from math import sqrt
3from sww_domain_shm 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.28 #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=1.0
48t0=time.time()
49i=1
50
51while 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
Note: See TracBrowser for help on using the repository browser.