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