1 | import os |
---|
2 | from math import sqrt, sin, cos, pi, exp |
---|
3 | from shallow_water_domain_suggestion3 import * |
---|
4 | from Numeric import zeros, Float |
---|
5 | |
---|
6 | print "TEST 1D-SOLUTION I" |
---|
7 | |
---|
8 | L=2.0 |
---|
9 | N=200 |
---|
10 | cell_len=L/N |
---|
11 | points=zeros(N+1, Float) |
---|
12 | for i in range(N+1): |
---|
13 | points[i]=i*cell_len |
---|
14 | |
---|
15 | domain=Domain(points) |
---|
16 | domain.order = 2 |
---|
17 | domain.set_timestepping_method('rk2') |
---|
18 | domain.cfl = 1.0 |
---|
19 | domain.limiter = "minmod" |
---|
20 | |
---|
21 | def 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 | |
---|
30 | def 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 | |
---|
39 | domain.set_quantity('stage',stage_perturb) |
---|
40 | domain.set_quantity('elevation', elevation_parabol) |
---|
41 | domain.set_boundary({'exterior':Reflective_boundary(domain)}) |
---|
42 | #X=domain.vertices |
---|
43 | #C=domain.centroids |
---|
44 | |
---|
45 | import time |
---|
46 | yieldstep=finaltime=0.25/100 |
---|
47 | t0=time.time() |
---|
48 | i=0 |
---|
49 | print "integral", domain.quantities['stage'].get_integral() |
---|
50 | while 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 | |
---|