source: anuga_work/development/anuga_1d/anuga_1d_suggestion1_2_3/perturbation1.py @ 7576

Last change on this file since 7576 was 7576, checked in by sudi, 14 years ago
File size: 2.6 KB
Line 
1import os
2from math import sqrt, sin, cos, pi, exp
3from shallow_water_domain_suggestion3 import *
4from Numeric import zeros, Float
5
6print "TEST 1D-SOLUTION I"
7
8L=2.0
9N=200
10cell_len=L/N
11points=zeros(N+1, Float)
12for i in range(N+1):
13    points[i]=i*cell_len
14
15domain=Domain(points)
16domain.order = 2
17domain.set_timestepping_method('rk2')
18domain.cfl = 1.0
19domain.limiter = "minmod"
20
21def stage_perturb(x):
22    y=zeros(len(x), Float)
23    for i in range(len(x)):
24        if 1.1 <= x[i] <=1.2:
25            y[i]=1.001
26        else:
27            y[i]=1.0
28    return y
29
30def elevation_parabol(x):
31    y=zeros(len(x), Float)
32    for i in range(len(x)):
33        if 1.4 <= x[i] <= 1.6:
34            y[i] = 0.25*(1.0+cos(10*pi*(x[i]-0.5)))
35        else:
36            y[i] = 0.0
37    return y
38
39domain.set_quantity('stage',stage_perturb)
40domain.set_quantity('elevation', elevation_parabol)
41domain.set_boundary({'exterior':Reflective_boundary(domain)})
42#X=domain.vertices
43#C=domain.centroids
44
45import time
46yieldstep=finaltime=0.25/100
47t0=time.time()
48i=0
49print "integral", domain.quantities['stage'].get_integral()
50while finaltime < 0.25001/100:
51    for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
52        domain.write_time() 
53        X = domain.vertices
54        StageQ = domain.quantities['stage'].vertex_values
55        XmomQ = domain.quantities['xmomentum'].vertex_values
56        VelQ = domain.quantities['velocity'].vertex_values
57        ElevQ = domain.quantities['elevation'].vertex_values
58
59        from pylab import plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
60        hold(False)
61           
62        plot1 = subplot(311)
63        plot(X,StageQ,X,ElevQ) #(X,ElevQ,X,StageQ)
64        plot1.set_ylim([-0.1,1.1])
65        #plot1.set_ylim([0.9999,1.001])
66        legend(('Stage','Bed Elevation'),'center left', shadow=False)
67        ylabel('Stage')
68           
69        plot2 = subplot(312)
70        plot(X,XmomQ)
71        ylabel('Momentum')
72   
73        plot3 = subplot(313)
74        plot(X,VelQ)
75        #legend(('Numerical Solution',' '),'lower right', shadow=False)
76        xlabel('Position')
77        ylabel('Velocity')
78   
79        #print "The domain.limiter is", domain.limiter
80        #print "domain order", domain.order
81        #print 'That took %.2f seconds' %(time.time()-t0)
82        #print '============================================================================='
83 
84        #file = "perturb"
85        #file += str(finaltime)
86        #file += ".png"
87        #show()
88        #savefig(file)
89       
90        filename = "%s%04i%s" %("perturbation", i, ".eps")
91        savefig(filename)
92
93        finaltime = finaltime + 0.25/100
94        i=i+1
95
Note: See TracBrowser for help on using the repository browser.