Changeset 1913
 Timestamp:
 Oct 14, 2005, 7:57:07 AM (18 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

inundation/analytical solutions/Analytical solution_oblique_shock_01.py
r1909 r1913 31 31 from shallow_water import Domain, Constant_height 32 32 from shallow_water import Transmissive_boundary, Reflective_boundary,\ 33 Dirichlet_boundary 33 Dirichlet_boundary , Transmissive_Momentum_Set_Stage_boundary 34 34 35 35 from math import sqrt, cos, sin, pi … … 42 42 leny = 30. 43 43 lenx = 40. 44 n = 50 45 m = 60 46 n=n*4 47 m=m*4 44 f = 4 45 n = 60*f 46 m = 80*f 48 47 49 48 points, elements, boundary = oblique_cross(m, n, lenx, leny) … … 52 51 # Order of solver 53 52 domain.default_order=2 53 domain.beta_w = 0.8 54 54 55 55 # Store output … … 61 61 62 62 # Provide file name for storing output 63 domain.filename="oblique "63 domain.filename="oblique_cross" 64 64 65 65 # Visualization smoothing 66 66 domain.smooth=True 67 domain.visualise=True67 #domain.visualise=True 68 68 69 69 ####################### … … 71 71 def x_slope(x, y): 72 72 return 0*x 73 74 def shock_hb(t): 75 return 2.5 76 77 def shock_hh(t): 78 return 2.5 79 80 def shock_pb(t): 81 return 10 82 83 def shock_pp(t): 84 return 10 85 86 87 def shock_h(x,y): 88 n = x.shape[0] 89 w = 0*x 90 for i in range(n): 91 if x[i]<0: 92 w[i] = shock_hh(0.0) 93 else: 94 w[i] = 1.0 95 return w 96 97 98 def shock_p(x,y): 99 n = x.shape[0] 100 w = 0*x 101 for i in range(n): 102 if x[i]<0: 103 w[i] = shock_pp(0.0) 104 else: 105 w[i] = 0.0 106 return w 107 73 108 74 109 domain.set_quantity('elevation', x_slope) … … 80 115 R = Reflective_boundary(domain) 81 116 T = Transmissive_boundary(domain) 82 D = Dirichlet_boundary([ 1.82, 8.49, 0.0])83 117 D = Dirichlet_boundary([shock_hb(0.0) , shock_pb(0.0), 0.0]) 118 Bts = Transmissive_Momentum_Set_Stage_boundary(domain, shock_hb) 84 119 85 120 domain.set_boundary({'left': D, 'right': T, 'top': R, 'bottom': R}) … … 87 122 ###################### 88 123 #Initial condition 89 h = 0.5 90 domain.set_quantity('stage', Constant_height(x_slope, h) ) 124 125 domain.set_quantity('stage', shock_h ) 126 domain.set_quantity('xmomentum',shock_p ) 91 127 92 128 … … 99 135 import time 100 136 t0 = time.time() 101 for t in domain.evolve(yieldstep = 0. 1, finaltime = 5):137 for t in domain.evolve(yieldstep = 0.5, finaltime = 5.0): 102 138 domain.write_time() 103 print domain.quantities['stage'].get_values(location='centroids',indices=[3399]) 104 print domain.quantities['xmomentum'].get_values(location='centroids',indices=[3399]) 105 print domain.quantities['ymomentum'].get_values(location='centroids',indices=[3399]) 139 id = 3399 140 print domain.quantities['stage'].get_values(location='centroids',indices=[id]) 141 print domain.quantities['xmomentum'].get_values(location='centroids',indices=[id]) 142 print domain.quantities['ymomentum'].get_values(location='centroids',indices=[id]) 143 id = 12719 144 print domain.quantities['stage'].get_values(location='centroids',indices=[id]) 145 print domain.quantities['xmomentum'].get_values(location='centroids',indices=[id]) 146 print domain.quantities['ymomentum'].get_values(location='centroids',indices=[id]) 106 147 107 148 print 'That took %.2f seconds' %(time.time()t0) … … 114 155 #im = ImageGrab.grab() 115 156 #im.save("ccube.eps") 116
Note: See TracChangeset
for help on using the changeset viewer.