source: anuga_work/development/MISG_2008/sensitivity_study.py @ 5245

Last change on this file since 5245 was 5245, checked in by sexton, 16 years ago

(1) update of production processes document (2) test for generating points for URS output (3) update of graduate proposal

File size: 7.5 KB
Line 
1"""Simple water flow example using ANUGA
2
3Water driven up a linear slope and time varying boundary,
4similar to a beach environment
5
6Test case for MISG 2008, University of Wollongong
7"""
8
9
10#------------------------------------------------------------------------------
11# Import necessary modules
12#------------------------------------------------------------------------------
13
14from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
15from anuga.shallow_water import Domain
16from anuga.shallow_water import Reflective_boundary
17from anuga.shallow_water import Dirichlet_boundary
18from anuga.shallow_water import Time_boundary
19from anuga.shallow_water import Transmissive_boundary
20from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
21
22from anuga.pmesh.mesh_interface import create_mesh_from_regions
23
24from Numeric import zeros, Float
25from RandomArray import normal, seed
26
27seed(13, 17) # Ensure random number sequence is reproducible
28
29#------------------------------------------------------------------------------
30# Setup computational domain
31#------------------------------------------------------------------------------
32y_extent=2000.0
33bounding_polygon = [[0,0],[2000,0],[2000,2000],[0,2000]]
34beach_polygon = [[5,900],[225,900],[225,1100],[5,1100]]
35interior_regions = [[beach_polygon, 50]]
36create_mesh_from_regions(bounding_polygon,
37                         boundary_tags={'bottom':[0],
38                                        'right':[1],
39                                        'top':[2],
40                                        'left':[3]},
41                         maximum_triangle_area=2000,
42                         interior_regions= interior_regions,
43                         filename='misg.msh',
44                         use_cache=False,
45                         verbose=True)
46domain = Domain('misg.msh',use_cache=False, verbose=True)
47
48
49domain.set_name('sensitivity')              # Output to file runup.sww
50domain.set_datadir('.')                     # Use current directory for output
51domain.tight_slope_limiters = 1
52domain.beta_h = 0
53domain.set_minimum_storable_height(0.01)
54
55print 'number of elements: ', len(domain)
56
57#------------------------------------------------------------------------------
58# Setup initial conditions
59#------------------------------------------------------------------------------
60
61def topography(x,y):
62
63    N = len(x)
64    z = zeros(N, Float)
65   
66    for i in range(N):
67        if x[i] <= (depth-c)/m:
68            z[i] = m*x[i]+c
69        else:
70        # linear bed slope           
71            #z[i] = slope*x[i] + depth-(slope)*(depth-c)/m
72            z[i] = slope*x[i] + depth-(slope)*(depth-c)/m + 1.0*sin((x[i]-(depth-c)/m)/T) + 1.0*sin((y[i]-(depth-c)/m)/T)
73
74            # IID noise
75            #z[i] += normal(-100.0, 5.0)
76       
77    return z
78
79#------------------------------------------------------------------------------
80# Setup boundary conditions
81#------------------------------------------------------------------------------
82
83from math import sin, pi, exp
84Br = Reflective_boundary(domain)      # Solid reflective wall
85Bt = Transmissive_boundary(domain)    # Continue all values on boundary
86Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values
87Bw = Time_boundary(domain=domain,     # Time dependent boundary 
88                   f=lambda t: [0.0, 0.0, 0.0])
89
90from Numeric import zeros, Float, arange
91start=-1.0/40.0
92finish=-1.8/40.0
93n_vec=arange(start,finish,-2.0/400.0)
94N=len(n_vec)
95#N=1 # use for 2.5% and 5% change tests
96max0=zeros(N, Float)
97max1=zeros(N, Float)
98max2=zeros(N, Float)
99max3=zeros(N, Float)
100mmax0=zeros(N, Float)
101mmax1=zeros(N, Float)
102mmax2=zeros(N, Float)
103mmax3=zeros(N, Float)
104output_file ='misg_influence_data.csv'
105domain.store=True
106fid_out = open(output_file, 'w')
107s = 'variable, h0, h1, h2, h3, v0, v1, v2, v3\n'
108fid_out.write(s)
109
110for i in range(N):
111    manning = 0.01*1.0
112    A=1.0*1.0
113    w=1.0/100.0
114    m=-1.0/20.0
115    depth=-10.0
116    c=10.0
117    T=100.0*1.0
118    #slope=-1.0/40.0*1.0
119    slope = n_vec[i]
120
121    var = slope
122    domain.set_quantity('elevation', topography) # Use function for elevation
123    domain.set_quantity('friction', manning)
124    domain.set_quantity('stage', 0.0)            # Constant negative initial stage
125    domain.set_quantity('xmomentum', 0.0) # Use function for elevation
126    domain.set_quantity('ymomentum', 0.0) # Use function for elevation
127    swwname = 'sensitivity_2d_wavey_bottom'#+str(i)
128    domain.set_name(swwname)
129    Btd = Transmissive_Momentum_Set_Stage_boundary(domain=domain,
130                                                   function=lambda t: [(0<t<pi/w)*A*sin(w*t), 0.0, 0.0])
131    domain.set_boundary({'left': Br, 'right': Btd, 'top': Br, 'bottom': Br})
132
133    # locations to record wave height and velocity
134    tri_id0 = domain.get_triangle_containing_point([-c/m,y_extent*0.5])
135    tri_id1 = domain.get_triangle_containing_point([-c/m-10,y_extent*0.5])
136    tri_id2 = domain.get_triangle_containing_point([-c/m-20,y_extent*0.5])
137    tri_id3 = domain.get_triangle_containing_point([-c/m+50,y_extent*0.5])
138
139    print 'hello', -c/m, -c/m-10, -c/m-20, -c/m+50
140   
141    stage_centroid     = domain.quantities['stage'].centroid_values
142    xmom_centroid      = domain.quantities['xmomentum'].centroid_values
143    elevation_centroid = domain.quantities['elevation'].centroid_values
144   
145    elev0 = elevation_centroid[tri_id0]
146    elev1 = elevation_centroid[tri_id1]
147    elev2 = elevation_centroid[tri_id2]
148    elev3 = elevation_centroid[tri_id3]
149
150    print 'hello again', elev0, elev1, elev2, elev3
151    import time
152    t0 = time.time()
153    #------------------------------------------------------------------------------
154    # Evolve system through time
155    #------------------------------------------------------------------------------
156    domain.time=0.0
157    finaltime=4.0 #pi/w+150.0
158    max=0
159    for t in domain.evolve(yieldstep = 1.0, finaltime = finaltime):
160
161        if  stage_centroid[tri_id0]-elev0>max0[i]: 
162            max0[i]=stage_centroid[tri_id0]-elev0
163        if  stage_centroid[tri_id1]-elev1>max1[i]:
164            max1[i]=stage_centroid[tri_id1]-elev1
165        if  stage_centroid[tri_id2]-elev2>max2[i]:
166            max2[i]=stage_centroid[tri_id2]-elev2
167        if  stage_centroid[tri_id3]-elev3>max3[i]:
168            max3[i]=stage_centroid[tri_id3]-elev3
169           
170        if  abs(xmom_centroid[tri_id0])>mmax0[i]: 
171            mmax0[i]=abs(xmom_centroid[tri_id0])
172        if  abs(xmom_centroid[tri_id1])>mmax1[i]:
173            mmax1[i]=abs(xmom_centroid[tri_id1])
174        if  abs(xmom_centroid[tri_id2])>mmax2[i]:
175            mmax2[i]=abs(xmom_centroid[tri_id2])
176        if  abs(xmom_centroid[tri_id3])>mmax3[i]:
177            mmax3[i]=abs(xmom_centroid[tri_id3])
178
179    print 'finished cycle: ', i
180    print 'That took %.2f seconds' %(time.time()-t0)
181
182
183    s = '%.6f, %.6f, %.6f, %.6f, %.6f, %.6f, %.6f, %.6f, %.6f\n' %(var, max0[i], max1[i], max2[i], max3[i], mmax0[i], mmax1[i], mmax2[i], mmax3[i])
184    fid_out.write(s)   
185
186from pylab import plot, ion, savefig, xlabel, ylabel, title, close, legend,hist,subplot
187ion()
188plot(n_vec,max2,'bo-',n_vec,max1,'go-',n_vec,max0,'ro-',n_vec,max3,'co-')
189title('MISG testing: varying slope')
190xlabel('slope')
191ylabel('depth (m)')
192#legend(['loc2','loc1','loc0','loc3'])
193name='bed_slope'
194savefig(name)
195'''
196#title('MISG testing: varying slope PDF')
197subplot(2,2,1)
198[h0,bins,patches]=hist(max0,10)
199subplot(2,2,2)
200[h1,bins,pa3ches]=hist(max1,10)
201subplot(2,2,3)
202[h2,bins,patches]=hist(max2,10)
203subplot(2,2,4)
204[h3,bins,patches]=hist(max3,10)
205#plot(bins,h0,'bo-',bins,h1,'ro-',bins,h2,'go-',bins,h3,'co-')
206xlabel('Runup')
207ylabel('Frequency')
208name='bed_slope_PDF'
209savefig(name)
210'''
211
212close('all')
Note: See TracBrowser for help on using the repository browser.