[7837] | 1 | """Example of shallow water wave equation analytical solution |
---|
| 2 | consists of a parabolic profile in a parabolic basin. Analytical |
---|
| 3 | solutiuon to this problem was derived by Thacker |
---|
| 4 | and used by Yoon and Chou. |
---|
| 5 | |
---|
| 6 | Copyright 2005 |
---|
| 7 | Christopher Zoppou, Stephen Roberts, ANU, Geoscience Australia |
---|
| 8 | |
---|
| 9 | """ |
---|
| 10 | |
---|
| 11 | #--------------- |
---|
| 12 | # Module imports |
---|
[7886] | 13 | from anuga import Domain, Transmissive_boundary, Reflective_boundary,\ |
---|
| 14 | Dirichlet_boundary, rectangular_cross |
---|
| 15 | |
---|
| 16 | #from anuga.interface import Domain, Transmissive_boundary, Reflective_boundary,\ |
---|
[7837] | 17 | # Dirichlet_boundary |
---|
| 18 | from math import sqrt, cos, sin, pi |
---|
[7886] | 19 | #from anuga.interface import rectangular_cross |
---|
| 20 | #from anuga.utilities.polygon import inside_polygon, is_inside_triangle |
---|
[7837] | 21 | from numpy import asarray |
---|
| 22 | |
---|
| 23 | |
---|
| 24 | #------------------------------- |
---|
| 25 | # Domain |
---|
[7886] | 26 | n = 4#200 |
---|
| 27 | m = 4#200 |
---|
[7837] | 28 | lenx = 8000.0 |
---|
| 29 | leny = 8000.0 |
---|
| 30 | origin = (-4000.0, -4000.0) |
---|
| 31 | |
---|
| 32 | points, elements, boundary = rectangular_cross(m, n, lenx, leny, origin) |
---|
[7886] | 33 | |
---|
| 34 | print points |
---|
| 35 | print elements |
---|
| 36 | print boundary |
---|
| 37 | |
---|
[7837] | 38 | domain = Domain(points, elements, boundary) |
---|
| 39 | |
---|
| 40 | #---------------- |
---|
| 41 | # Order of scheme |
---|
| 42 | # Good compromise between |
---|
| 43 | # limiting and CFL |
---|
| 44 | #--------------- |
---|
| 45 | domain.set_default_order(2) |
---|
| 46 | domain.set_timestepping_method(2) |
---|
| 47 | domain.set_beta(0.7) |
---|
| 48 | domain.set_CFL(0.6) |
---|
| 49 | |
---|
| 50 | domain.smooth = True |
---|
| 51 | |
---|
| 52 | #------------------------------------- |
---|
| 53 | # Provide file name for storing output |
---|
| 54 | domain.store = False |
---|
| 55 | domain.format = 'sww' |
---|
| 56 | domain.set_name('yoon_mesh_second_order_cross') |
---|
| 57 | print 'Number of triangles = ', len(domain) |
---|
| 58 | |
---|
| 59 | #---------------------------------------------------------- |
---|
| 60 | # Decide which quantities are to be stored at each timestep |
---|
| 61 | domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum'] |
---|
| 62 | |
---|
| 63 | #------------------------------------------ |
---|
| 64 | # Reduction operation for get_vertex_values |
---|
| 65 | from anuga.utilities.numerical_tools import mean |
---|
| 66 | domain.reduction = mean #domain.reduction = min #Looks better near steep slopes |
---|
| 67 | |
---|
| 68 | #------------------ |
---|
| 69 | # Initial condition |
---|
| 70 | print 'Initial condition' |
---|
| 71 | t = 0.0 |
---|
| 72 | D0 = 1. |
---|
| 73 | L = 2500. |
---|
| 74 | R0 = 2000. |
---|
| 75 | g = 9.81 |
---|
| 76 | |
---|
| 77 | A = (L**4 - R0**4)/(L**4 + R0**4) |
---|
| 78 | omega = 2./L*sqrt(2.*g*D0) |
---|
| 79 | T = pi/omega |
---|
| 80 | |
---|
| 81 | #------------------ |
---|
| 82 | # Set bed elevation |
---|
| 83 | def x_slope(x,y): |
---|
| 84 | n = x.shape[0] |
---|
| 85 | z = 0*x |
---|
| 86 | for i in range(n): |
---|
| 87 | r = sqrt(x[i]*x[i] + y[i]*y[i]) |
---|
| 88 | z[i] = -D0*(1.-r*r/L/L) |
---|
| 89 | return z |
---|
| 90 | domain.set_quantity('elevation', x_slope) |
---|
| 91 | |
---|
| 92 | #---------------------------- |
---|
| 93 | # Set the initial water level |
---|
| 94 | def level(x,y): |
---|
| 95 | z = x_slope(x,y) |
---|
| 96 | n = x.shape[0] |
---|
| 97 | h = 0*x |
---|
| 98 | for i in range(n): |
---|
| 99 | r = sqrt(x[i]*x[i] + y[i]*y[i]) |
---|
| 100 | h[i] = D0*((sqrt(1-A*A))/(1.-A*cos(omega*t)) |
---|
| 101 | -1.-r*r/L/L*((1.-A*A)/((1.-A*cos(omega*t))**2)-1.)) |
---|
| 102 | if h[i] < z[i]: |
---|
| 103 | h[i] = z[i] |
---|
| 104 | return h |
---|
| 105 | domain.set_quantity('stage', level) |
---|
| 106 | |
---|
| 107 | #--------- |
---|
| 108 | # Boundary |
---|
| 109 | print 'Boundary conditions' |
---|
| 110 | R = Reflective_boundary(domain) |
---|
| 111 | T = Transmissive_boundary(domain) |
---|
| 112 | D = Dirichlet_boundary([0.0, 0.0, 0.0]) |
---|
| 113 | domain.set_boundary({'left': D, 'right': D, 'top': D, 'bottom': D}) |
---|
| 114 | |
---|
| 115 | #--------------------------------------------- |
---|
| 116 | # Find triangle that contains the point points |
---|
| 117 | # and print to file |
---|
| 118 | |
---|
| 119 | |
---|
| 120 | |
---|
| 121 | points = (0.0, 0.0) |
---|
| 122 | for n in range(len(domain.triangles)): |
---|
| 123 | tri = domain.get_vertex_coordinates(n) |
---|
| 124 | |
---|
| 125 | if is_inside_triangle(points,tri): |
---|
| 126 | #print 'Point is within triangle with vertices '+'%s'%tri |
---|
| 127 | n_point = n |
---|
| 128 | |
---|
| 129 | print 'n_point = ',n_point |
---|
| 130 | t = domain.triangles[n_point] |
---|
| 131 | print 't = ', t |
---|
| 132 | tri = domain.get_vertex_coordinates(n) |
---|
| 133 | |
---|
| 134 | filename=domain.get_name() |
---|
| 135 | file = open(filename,'w') |
---|
| 136 | |
---|
| 137 | #---------- |
---|
| 138 | # Evolution |
---|
| 139 | import time |
---|
| 140 | t0 = time.time() |
---|
| 141 | |
---|
| 142 | time_array = [] |
---|
| 143 | stage_array = [] |
---|
| 144 | Stage = domain.quantities['stage'] |
---|
| 145 | Xmomentum = domain.quantities['xmomentum'] |
---|
| 146 | Ymomentum = domain.quantities['ymomentum'] |
---|
| 147 | |
---|
| 148 | for t in domain.evolve(yieldstep = 20.0, finaltime = 17700.0 ): |
---|
| 149 | domain.write_time() |
---|
| 150 | |
---|
| 151 | #tri_array = asarray(tri) |
---|
| 152 | #t_array = asarray([[0,1,2]]) |
---|
| 153 | #interp = Interpolation(tri_array,t_array,[points]) |
---|
| 154 | |
---|
| 155 | |
---|
| 156 | stage = Stage.get_values(location='centroids',indices=[n_point]) |
---|
| 157 | xmomentum = Xmomentum.get_values(location='centroids',indices=[n_point]) |
---|
| 158 | ymomentum = Ymomentum.get_values(location='centroids',indices=[n_point]) |
---|
| 159 | #print '%10.6f %10.6f %10.6f %10.6f\n'%(t,stage,xmomentum,ymomentum) |
---|
| 160 | file.write( '%10.6f %10.6f %10.6f %10.6f\n'%(t,stage,xmomentum,ymomentum) ) |
---|
| 161 | |
---|
| 162 | time_array.append(t) |
---|
| 163 | stage_array.append(stage) |
---|
| 164 | |
---|
| 165 | file.close() |
---|
| 166 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
| 167 | |
---|
| 168 | |
---|
| 169 | from pylab import * |
---|
| 170 | ion() |
---|
| 171 | hold(False) |
---|
| 172 | plot(time_array, stage_array, 'r.-') |
---|
| 173 | #title('Gauge %s' %name) |
---|
| 174 | xlabel('time(s)') |
---|
| 175 | ylabel('stage (m)') |
---|
| 176 | #legend(('Observed', 'Modelled'), shadow=True, loc='upper left') |
---|
| 177 | #savefig(name, dpi = 300) |
---|
| 178 | |
---|
| 179 | raw_input('Next') |
---|
| 180 | show() |
---|
| 181 | |
---|
| 182 | |
---|