source: trunk/anuga_work/development/MISG_2008/sensitivity_study.py @ 7884

Last change on this file since 7884 was 5442, checked in by ole, 16 years ago

Retired h-limiter and beta_h as per ticket:194.
All unit tests and validation tests pass.

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.set_minimum_storable_height(0.01)
52
53print 'number of elements: ', len(domain)
54
55#------------------------------------------------------------------------------
56# Setup initial conditions
57#------------------------------------------------------------------------------
58
59def topography(x,y):
60
61    N = len(x)
62    z = zeros(N, Float)
63   
64    for i in range(N):
65        if x[i] <= (depth-c)/m:
66            z[i] = m*x[i]+c
67        else:
68        # linear bed slope           
69            #z[i] = slope*x[i] + depth-(slope)*(depth-c)/m
70            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)
71
72            # IID noise
73            #z[i] += normal(-100.0, 5.0)
74       
75    return z
76
77#------------------------------------------------------------------------------
78# Setup boundary conditions
79#------------------------------------------------------------------------------
80
81from math import sin, pi, exp
82Br = Reflective_boundary(domain)      # Solid reflective wall
83Bt = Transmissive_boundary(domain)    # Continue all values on boundary
84Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values
85Bw = Time_boundary(domain=domain,     # Time dependent boundary 
86                   f=lambda t: [0.0, 0.0, 0.0])
87
88from Numeric import zeros, Float, arange
89start=-1.0/40.0
90finish=-1.8/40.0
91n_vec=arange(start,finish,-2.0/400.0)
92N=len(n_vec)
93#N=1 # use for 2.5% and 5% change tests
94max0=zeros(N, Float)
95max1=zeros(N, Float)
96max2=zeros(N, Float)
97max3=zeros(N, Float)
98mmax0=zeros(N, Float)
99mmax1=zeros(N, Float)
100mmax2=zeros(N, Float)
101mmax3=zeros(N, Float)
102output_file ='misg_influence_data.csv'
103domain.store=True
104fid_out = open(output_file, 'w')
105s = 'variable, h0, h1, h2, h3, v0, v1, v2, v3\n'
106fid_out.write(s)
107
108for i in range(N):
109    manning = 0.01*1.0
110    A=1.0*1.0
111    w=1.0/100.0
112    m=-1.0/20.0
113    depth=-10.0
114    c=10.0
115    T=100.0*1.0
116    #slope=-1.0/40.0*1.0
117    slope = n_vec[i]
118
119    var = slope
120    domain.set_quantity('elevation', topography) # Use function for elevation
121    domain.set_quantity('friction', manning)
122    domain.set_quantity('stage', 0.0)            # Constant negative initial stage
123    domain.set_quantity('xmomentum', 0.0) # Use function for elevation
124    domain.set_quantity('ymomentum', 0.0) # Use function for elevation
125    swwname = 'sensitivity_2d_wavey_bottom'#+str(i)
126    domain.set_name(swwname)
127    Btd = Transmissive_Momentum_Set_Stage_boundary(domain=domain,
128                                                   function=lambda t: [(0<t<pi/w)*A*sin(w*t), 0.0, 0.0])
129    domain.set_boundary({'left': Br, 'right': Btd, 'top': Br, 'bottom': Br})
130
131    # locations to record wave height and velocity
132    tri_id0 = domain.get_triangle_containing_point([-c/m,y_extent*0.5])
133    tri_id1 = domain.get_triangle_containing_point([-c/m-10,y_extent*0.5])
134    tri_id2 = domain.get_triangle_containing_point([-c/m-20,y_extent*0.5])
135    tri_id3 = domain.get_triangle_containing_point([-c/m+50,y_extent*0.5])
136
137    print 'hello', -c/m, -c/m-10, -c/m-20, -c/m+50
138   
139    stage_centroid     = domain.quantities['stage'].centroid_values
140    xmom_centroid      = domain.quantities['xmomentum'].centroid_values
141    elevation_centroid = domain.quantities['elevation'].centroid_values
142   
143    elev0 = elevation_centroid[tri_id0]
144    elev1 = elevation_centroid[tri_id1]
145    elev2 = elevation_centroid[tri_id2]
146    elev3 = elevation_centroid[tri_id3]
147
148    print 'hello again', elev0, elev1, elev2, elev3
149    import time
150    t0 = time.time()
151    #------------------------------------------------------------------------------
152    # Evolve system through time
153    #------------------------------------------------------------------------------
154    domain.time=0.0
155    finaltime=4.0 #pi/w+150.0
156    max=0
157    for t in domain.evolve(yieldstep = 1.0, finaltime = finaltime):
158
159        if  stage_centroid[tri_id0]-elev0>max0[i]: 
160            max0[i]=stage_centroid[tri_id0]-elev0
161        if  stage_centroid[tri_id1]-elev1>max1[i]:
162            max1[i]=stage_centroid[tri_id1]-elev1
163        if  stage_centroid[tri_id2]-elev2>max2[i]:
164            max2[i]=stage_centroid[tri_id2]-elev2
165        if  stage_centroid[tri_id3]-elev3>max3[i]:
166            max3[i]=stage_centroid[tri_id3]-elev3
167           
168        if  abs(xmom_centroid[tri_id0])>mmax0[i]: 
169            mmax0[i]=abs(xmom_centroid[tri_id0])
170        if  abs(xmom_centroid[tri_id1])>mmax1[i]:
171            mmax1[i]=abs(xmom_centroid[tri_id1])
172        if  abs(xmom_centroid[tri_id2])>mmax2[i]:
173            mmax2[i]=abs(xmom_centroid[tri_id2])
174        if  abs(xmom_centroid[tri_id3])>mmax3[i]:
175            mmax3[i]=abs(xmom_centroid[tri_id3])
176
177    print 'finished cycle: ', i
178    print 'That took %.2f seconds' %(time.time()-t0)
179
180
181    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])
182    fid_out.write(s)   
183
184from pylab import plot, ion, savefig, xlabel, ylabel, title, close, legend,hist,subplot
185ion()
186plot(n_vec,max2,'bo-',n_vec,max1,'go-',n_vec,max0,'ro-',n_vec,max3,'co-')
187title('MISG testing: varying slope')
188xlabel('slope')
189ylabel('depth (m)')
190#legend(['loc2','loc1','loc0','loc3'])
191name='bed_slope'
192savefig(name)
193'''
194#title('MISG testing: varying slope PDF')
195subplot(2,2,1)
196[h0,bins,patches]=hist(max0,10)
197subplot(2,2,2)
198[h1,bins,pa3ches]=hist(max1,10)
199subplot(2,2,3)
200[h2,bins,patches]=hist(max2,10)
201subplot(2,2,4)
202[h3,bins,patches]=hist(max3,10)
203#plot(bins,h0,'bo-',bins,h1,'ro-',bins,h2,'go-',bins,h3,'co-')
204xlabel('Runup')
205ylabel('Frequency')
206name='bed_slope_PDF'
207savefig(name)
208'''
209
210close('all')
Note: See TracBrowser for help on using the repository browser.