Changeset 7624
 Timestamp:
 Feb 16, 2010, 9:33:23 AM (13 years ago)
 Location:
 anuga_validation/analytical solutions
 Files:

 1 edited
 2 moved
Legend:
 Unmodified
 Added
 Removed

anuga_validation/analytical solutions/Analytical_solution_Yoon_cross_mesh.py
r3846 r7624 1 1 """Example of shallow water wave equation analytical solution 2 2 consists of a parabolic profile in a parabolic basin. Analytical 3 solutiuon to this problem was derived by Carrier and Greenspan3 solutiuon to this problem was derived by Thacker 4 4 and used by Yoon and Chou. 5 5 … … 14 14 #from os import sep 15 15 #sys.path.append('..'+sep+'pyvolution') 16 from anuga. pyvolution.shallow_waterimport Domain, Transmissive_boundary, Reflective_boundary,\16 from anuga.shallow_water_balanced.swb_domain import Domain, Transmissive_boundary, Reflective_boundary,\ 17 17 Dirichlet_boundary 18 19 #from anuga.interface import Domain, Transmissive_boundary, Reflective_boundary,\ 20 # Dirichlet_boundary 18 21 from math import sqrt, cos, sin, pi 19 from anuga. pyvolution.mesh_factoryimport rectangular_cross20 from anuga.utilities.polygon import inside_polygon 21 from Numericimport asarray22 from anuga.pyvolution.least_squares import Interpolation 22 from anuga.interface import rectangular_cross 23 from anuga.utilities.polygon import inside_polygon, is_inside_triangle 24 from numpy import asarray 25 23 26 24 27 # 25 28 # Domain 26 n = 100 27 m = 100 28 delta_x = 80.0 29 delta_y = 80.0 30 lenx = delta_x*n 31 leny = delta_y*m 29 n = 200 30 m = 200 31 lenx = 8000.0 32 leny = 8000.0 32 33 origin = (4000.0, 4000.0) 33 34 … … 37 38 # 38 39 # Order of scheme 39 domain.default_order = 1 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) 40 47 41 48 domain.smooth = True … … 54 61 # 55 62 # Reduction operation for get_vertex_values 56 from anuga. pyvolution.utilimport mean63 from anuga.utilities.numerical_tools import mean 57 64 domain.reduction = mean #domain.reduction = min #Looks better near steep slopes 58 65 … … 107 114 # Find triangle that contains the point points 108 115 # and print to file 109 points = [0.,0.] 116 117 118 119 points = (0.0, 0.0) 110 120 for n in range(len(domain.triangles)): 111 t = domain.triangles[n] 112 tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]] 121 tri = domain.get_vertex_coordinates(n) 113 122 114 if i nside_polygon(points,tri):115 print 'Point is within triangle with vertices '+'%s'%tri123 if is_inside_triangle(points,tri): 124 #print 'Point is within triangle with vertices '+'%s'%tri 116 125 n_point = n 117 126 118 127 print 'n_point = ',n_point 119 128 t = domain.triangles[n_point] 120 tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]] 129 print 't = ', t 130 tri = domain.get_vertex_coordinates(n) 121 131 122 132 filename=domain.get_name() … … 137 147 domain.write_time() 138 148 139 tri_array = asarray(tri)140 t_array = asarray([[0,1,2]])141 interp = Interpolation(tri_array,t_array,[points])149 #tri_array = asarray(tri) 150 #t_array = asarray([[0,1,2]]) 151 #interp = Interpolation(tri_array,t_array,[points]) 142 152 143 153 144 stage = Stage.get_values(location='centroids',indices=[n_point])[0] 145 xmomentum = Xmomentum.get_values(location='centroids',indices=[n_point])[0] 146 ymomentum = Ymomentum.get_values(location='centroids',indices=[n_point])[0] 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) 147 158 file.write( '%10.6f %10.6f %10.6f %10.6f\n'%(t,stage,xmomentum,ymomentum) ) 148 159 … … 161 172 xlabel('time(s)') 162 173 ylabel('stage (m)') 163 legend(('Observed', 'Modelled'), shadow=True, loc='upper left')174 #legend(('Observed', 'Modelled'), shadow=True, loc='upper left') 164 175 #savefig(name, dpi = 300) 165 176
Note: See TracChangeset
for help on using the changeset viewer.