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