[4937] | 1 | """Simple water flow example using ANUGA |
---|
| 2 | |
---|
| 3 | Water driven up a linear slope and time varying boundary, |
---|
| 4 | similar to a beach environment |
---|
[5066] | 5 | |
---|
| 6 | Test case for MISG 2008, University of Wollongong |
---|
[4937] | 7 | """ |
---|
| 8 | |
---|
| 9 | |
---|
| 10 | #------------------------------------------------------------------------------ |
---|
| 11 | # Import necessary modules |
---|
| 12 | #------------------------------------------------------------------------------ |
---|
| 13 | |
---|
| 14 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross |
---|
| 15 | from anuga.shallow_water import Domain |
---|
| 16 | from anuga.shallow_water import Reflective_boundary |
---|
| 17 | from anuga.shallow_water import Dirichlet_boundary |
---|
| 18 | from anuga.shallow_water import Time_boundary |
---|
| 19 | from anuga.shallow_water import Transmissive_boundary |
---|
[5066] | 20 | from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary |
---|
[4937] | 21 | |
---|
[5066] | 22 | from anuga.pmesh.mesh_interface import create_mesh_from_regions |
---|
| 23 | |
---|
[4937] | 24 | from Numeric import zeros, Float |
---|
| 25 | from RandomArray import normal, seed |
---|
| 26 | |
---|
| 27 | seed(13, 17) # Ensure random number sequence is reproducible |
---|
| 28 | |
---|
| 29 | #------------------------------------------------------------------------------ |
---|
| 30 | # Setup computational domain |
---|
| 31 | #------------------------------------------------------------------------------ |
---|
[5066] | 32 | y_extent=2000.0 |
---|
| 33 | bounding_polygon = [[0,0],[2000,0],[2000,2000],[0,2000]] |
---|
| 34 | beach_polygon = [[5,900],[225,900],[225,1100],[5,1100]] |
---|
| 35 | interior_regions = [[beach_polygon, 50]] |
---|
| 36 | create_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) |
---|
| 46 | domain = Domain('misg.msh',use_cache=False, verbose=True) |
---|
[4937] | 47 | |
---|
| 48 | |
---|
| 49 | domain.set_name('sensitivity') # Output to file runup.sww |
---|
| 50 | domain.set_datadir('.') # Use current directory for output |
---|
[5066] | 51 | domain.set_minimum_storable_height(0.01) |
---|
[4937] | 52 | |
---|
[5066] | 53 | print 'number of elements: ', len(domain) |
---|
| 54 | |
---|
[4937] | 55 | #------------------------------------------------------------------------------ |
---|
| 56 | # Setup initial conditions |
---|
| 57 | #------------------------------------------------------------------------------ |
---|
| 58 | |
---|
| 59 | def topography(x,y): |
---|
| 60 | |
---|
| 61 | N = len(x) |
---|
| 62 | z = zeros(N, Float) |
---|
[5066] | 63 | |
---|
[4937] | 64 | for i in range(N): |
---|
[5066] | 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) |
---|
[4937] | 71 | |
---|
[5066] | 72 | # IID noise |
---|
| 73 | #z[i] += normal(-100.0, 5.0) |
---|
[4937] | 74 | |
---|
| 75 | return z |
---|
| 76 | |
---|
| 77 | #------------------------------------------------------------------------------ |
---|
| 78 | # Setup boundary conditions |
---|
| 79 | #------------------------------------------------------------------------------ |
---|
| 80 | |
---|
| 81 | from math import sin, pi, exp |
---|
| 82 | Br = Reflective_boundary(domain) # Solid reflective wall |
---|
| 83 | Bt = Transmissive_boundary(domain) # Continue all values on boundary |
---|
| 84 | Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values |
---|
| 85 | Bw = Time_boundary(domain=domain, # Time dependent boundary |
---|
[5066] | 86 | f=lambda t: [0.0, 0.0, 0.0]) |
---|
[4937] | 87 | |
---|
[5066] | 88 | from Numeric import zeros, Float, arange |
---|
| 89 | start=-1.0/40.0 |
---|
| 90 | finish=-1.8/40.0 |
---|
| 91 | n_vec=arange(start,finish,-2.0/400.0) |
---|
| 92 | N=len(n_vec) |
---|
| 93 | #N=1 # use for 2.5% and 5% change tests |
---|
| 94 | max0=zeros(N, Float) |
---|
| 95 | max1=zeros(N, Float) |
---|
| 96 | max2=zeros(N, Float) |
---|
| 97 | max3=zeros(N, Float) |
---|
| 98 | mmax0=zeros(N, Float) |
---|
| 99 | mmax1=zeros(N, Float) |
---|
| 100 | mmax2=zeros(N, Float) |
---|
| 101 | mmax3=zeros(N, Float) |
---|
| 102 | output_file ='misg_influence_data.csv' |
---|
| 103 | domain.store=True |
---|
| 104 | fid_out = open(output_file, 'w') |
---|
| 105 | s = 'variable, h0, h1, h2, h3, v0, v1, v2, v3\n' |
---|
| 106 | fid_out.write(s) |
---|
[4937] | 107 | |
---|
[5066] | 108 | for 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] |
---|
[4937] | 118 | |
---|
[5066] | 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}) |
---|
[4937] | 130 | |
---|
[5066] | 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]) |
---|
[5245] | 136 | |
---|
| 137 | print 'hello', -c/m, -c/m-10, -c/m-20, -c/m+50 |
---|
| 138 | |
---|
[5066] | 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] |
---|
[4937] | 147 | |
---|
[5245] | 148 | print 'hello again', elev0, elev1, elev2, elev3 |
---|
[5066] | 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): |
---|
[4937] | 158 | |
---|
[5066] | 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 | |
---|
| 184 | from pylab import plot, ion, savefig, xlabel, ylabel, title, close, legend,hist,subplot |
---|
| 185 | ion() |
---|
| 186 | plot(n_vec,max2,'bo-',n_vec,max1,'go-',n_vec,max0,'ro-',n_vec,max3,'co-') |
---|
| 187 | title('MISG testing: varying slope') |
---|
| 188 | xlabel('slope') |
---|
| 189 | ylabel('depth (m)') |
---|
| 190 | #legend(['loc2','loc1','loc0','loc3']) |
---|
| 191 | name='bed_slope' |
---|
| 192 | savefig(name) |
---|
| 193 | ''' |
---|
| 194 | #title('MISG testing: varying slope PDF') |
---|
| 195 | subplot(2,2,1) |
---|
| 196 | [h0,bins,patches]=hist(max0,10) |
---|
| 197 | subplot(2,2,2) |
---|
| 198 | [h1,bins,pa3ches]=hist(max1,10) |
---|
| 199 | subplot(2,2,3) |
---|
| 200 | [h2,bins,patches]=hist(max2,10) |
---|
| 201 | subplot(2,2,4) |
---|
| 202 | [h3,bins,patches]=hist(max3,10) |
---|
| 203 | #plot(bins,h0,'bo-',bins,h1,'ro-',bins,h2,'go-',bins,h3,'co-') |
---|
| 204 | xlabel('Runup') |
---|
| 205 | ylabel('Frequency') |
---|
| 206 | name='bed_slope_PDF' |
---|
| 207 | savefig(name) |
---|
| 208 | ''' |
---|
| 209 | |
---|
| 210 | close('all') |
---|