[3223] | 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 |
---|
| 5 | """ |
---|
| 6 | |
---|
| 7 | |
---|
| 8 | #------------------------------------------------------------------------------ |
---|
| 9 | # Import necessary modules |
---|
| 10 | #------------------------------------------------------------------------------ |
---|
| 11 | |
---|
[7808] | 12 | import anuga |
---|
[3223] | 13 | |
---|
[3247] | 14 | |
---|
[3223] | 15 | #------------------------------------------------------------------------------ |
---|
| 16 | # Setup computational domain |
---|
| 17 | #------------------------------------------------------------------------------ |
---|
| 18 | |
---|
[3257] | 19 | waveheight = 1800.0 |
---|
[3247] | 20 | depth_east_edge = -4000. |
---|
[3257] | 21 | timeend = 1000.0 |
---|
[3247] | 22 | west = 0. |
---|
| 23 | east = 80000. |
---|
| 24 | south = 0. |
---|
| 25 | north = 5000. |
---|
[3223] | 26 | polygon = [[east,north],[west,north],[west,south],[east,south]] |
---|
| 27 | meshname = 'test.msh' |
---|
| 28 | |
---|
[7808] | 29 | anuga.create_mesh_from_regions(polygon, |
---|
[3223] | 30 | boundary_tags={'top': [0], |
---|
| 31 | 'left': [1], |
---|
| 32 | 'bottom': [2], |
---|
| 33 | 'right': [3]}, |
---|
[3254] | 34 | maximum_triangle_area=200000, |
---|
[3223] | 35 | filename=meshname) |
---|
| 36 | #interior_regions=interior_regions) |
---|
| 37 | |
---|
[7808] | 38 | domain = anuga.Domain(meshname,use_cache=True,verbose = True) |
---|
[3223] | 39 | domain.set_name('test') # Output to test.sww |
---|
| 40 | domain.set_datadir('.') # Use current directory for output |
---|
| 41 | domain.set_quantities_to_be_stored(['stage',# Store all conserved quantities |
---|
| 42 | 'xmomentum', |
---|
| 43 | 'ymomentum']) |
---|
| 44 | |
---|
| 45 | |
---|
| 46 | #------------------------------------------------------------------------------ |
---|
| 47 | # Setup initial conditions |
---|
| 48 | #------------------------------------------------------------------------------ |
---|
| 49 | |
---|
| 50 | def topography(x,y): |
---|
[3247] | 51 | slope = depth_east_edge/((east-west)/2.) |
---|
| 52 | c = -(west+east)/2.*slope |
---|
| 53 | return slope*x+c # linear bed slope |
---|
[3223] | 54 | |
---|
| 55 | domain.set_quantity('elevation', topography) # Use function for elevation |
---|
| 56 | domain.set_quantity('friction', 0. ) # Constant friction |
---|
| 57 | domain.set_quantity('stage', 0.0) # Constant initial stage |
---|
| 58 | |
---|
| 59 | #------------------------------------------------------------------------------ |
---|
| 60 | # Setup boundary conditions |
---|
| 61 | #------------------------------------------------------------------------------ |
---|
| 62 | |
---|
| 63 | from math import sin, pi |
---|
[7808] | 64 | Br = anuga.Reflective_boundary(domain) # Solid reflective wall |
---|
| 65 | Bt = anuga.Transmissive_boundary(domain) # Continue all values on boundary |
---|
| 66 | Bd = anuga.Dirichlet_boundary([0.,0.,0.]) # Constant boundary values |
---|
| 67 | Bw = anuga.Time_boundary(domain=domain, # Time dependent boundary |
---|
[3223] | 68 | f=lambda t: [((10<t<20)*waveheight), 0.0, 0.0]) |
---|
| 69 | |
---|
| 70 | # Associate boundary tags with boundary objects |
---|
| 71 | domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom': Br}) |
---|
| 72 | |
---|
| 73 | |
---|
| 74 | #------------------------------------------------------------------------------ |
---|
| 75 | # Evolve system through time |
---|
| 76 | #------------------------------------------------------------------------------ |
---|
| 77 | |
---|
| 78 | stagestep = [] |
---|
[3257] | 79 | for t in domain.evolve(yieldstep = 10, finaltime = timeend): |
---|
[3223] | 80 | domain.write_time() |
---|
| 81 | |
---|
[3247] | 82 | #----------------------------------------------------------------------------- |
---|
| 83 | # Interrogate solution |
---|
| 84 | #----------------------------------------------------------------------------- |
---|
| 85 | |
---|
[3223] | 86 | # Gauge line |
---|
| 87 | def gauge_line(west,east,north,south): |
---|
[7808] | 88 | from numpy import arange |
---|
[3258] | 89 | gaugex = arange(west,(west+east)/2.,1000.) |
---|
[3223] | 90 | gauges = [] |
---|
[3247] | 91 | gaugey = [] |
---|
[3223] | 92 | for x in gaugex: |
---|
| 93 | y = (north+south)/2. |
---|
[3247] | 94 | gaugey.append(y) |
---|
[3223] | 95 | gauges.append([x,y]) |
---|
[3247] | 96 | return gauges, gaugex, gaugey |
---|
[3223] | 97 | |
---|
[3247] | 98 | gauges, gaugex, gaugey = gauge_line(west,east,north,south) |
---|
[3223] | 99 | |
---|
[7808] | 100 | f = anuga.file_function('test.sww', |
---|
[3247] | 101 | quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'], |
---|
| 102 | interpolation_points = gauges, |
---|
| 103 | verbose = True, |
---|
| 104 | use_cache = True) |
---|
| 105 | |
---|
[3256] | 106 | from pylab import plot, xlabel, ylabel, title, ion, close, savefig, figure, axis |
---|
[3247] | 107 | ion() |
---|
| 108 | |
---|
| 109 | maxw = [] |
---|
| 110 | minw = [] |
---|
| 111 | maxd = [] |
---|
| 112 | count = 0 |
---|
[3223] | 113 | for k, g in enumerate(gauges): |
---|
[3247] | 114 | stage = [] |
---|
| 115 | elevation = [] |
---|
| 116 | depths = [] |
---|
[3254] | 117 | bed = topography(g[0],g[1]) |
---|
[3223] | 118 | for i, t in enumerate(f.get_time()): |
---|
[3247] | 119 | w = f(t, point_id = k)[0] |
---|
| 120 | elev = f(t, point_id = k)[1] |
---|
| 121 | depth = w-elev |
---|
| 122 | stage.append(w) |
---|
| 123 | elevation.append(elev) |
---|
| 124 | depths.append(elev) |
---|
[3256] | 125 | |
---|
| 126 | if max(stage)-bed <= 0.2: |
---|
| 127 | count+=1 |
---|
| 128 | posx=g[0] |
---|
| 129 | loc = k |
---|
[3247] | 130 | maxw.append(max(stage)) |
---|
| 131 | minw.append(min(stage)) |
---|
| 132 | maxd.append(max(depths)) |
---|
[3256] | 133 | |
---|
[3223] | 134 | |
---|
[3247] | 135 | filename = 'maxrunup'+str(waveheight)+str(depth_east_edge)+'.csv' |
---|
[3223] | 136 | fid = open(filename,'w') |
---|
[3247] | 137 | s = 'Depth at eastern boundary,Waveheight,Runup distance,Runup height\n' |
---|
[3223] | 138 | fid.write(s) |
---|
[3247] | 139 | runup_distance = posx |
---|
| 140 | runup_height = topography(runup_distance,(north+south)/2.) |
---|
| 141 | print 'run up height: ', runup_height |
---|
| 142 | if runup_distance < (east+west)/2.: |
---|
| 143 | runup_distance_from_coast = (east+west)/2.-runup_distance |
---|
| 144 | else: |
---|
| 145 | runup_distance_from_coast = runup_distance-(east+west)/2. |
---|
| 146 | print 'run up distance from coastline: ', runup_distance_from_coast |
---|
| 147 | s = '%.2f,%.2f,%.2f,%.2f\n' %(depth_east_edge,waveheight,runup_distance_from_coast,runup_height) |
---|
[3223] | 148 | fid.write(s) |
---|
| 149 | |
---|
[3247] | 150 | figure(1) |
---|
[3254] | 151 | plot(gaugex,maxw,'g+',gaugex,topography(gaugex,(north+south)/2.),'r-') |
---|
[3223] | 152 | xlabel('x') |
---|
| 153 | ylabel('stage') |
---|
[3224] | 154 | title('Maximum stage for gauge line') |
---|
[3256] | 155 | #axis([33000, 47000, -1000, 3000]) |
---|
[3247] | 156 | savefig('max_stage') |
---|
| 157 | |
---|
| 158 | close('all') |
---|