Changeset 8178
- Timestamp:
- May 15, 2011, 4:22:40 PM (13 years ago)
- Location:
- trunk/anuga_core/source/anuga
- Files:
-
- 6 edited
-
shallow_water/shallow_water_domain.py (modified) (3 diffs)
-
shallow_water_balanced/run_step.py (modified) (4 diffs)
-
shallow_water_balanced/run_wave.py (modified) (1 diff)
-
shallow_water_balanced/run_wave_shelf.py (modified) (5 diffs)
-
shallow_water_balanced/swb_domain.py (modified) (8 diffs)
-
shallow_water_balanced/swb_domain_ext.c (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r8175 r8178 574 574 def distribute_to_vertices_and_edges(self): 575 575 """ Call correct module function """ 576 577 #Shortcuts 578 #W = self.quantities['stage'] 579 #Z = self.quantities['elevation'] 580 581 #Arrays 582 #w_C = W.centroid_values 583 #z_C = Z.centroid_values 584 585 #num_min = num.min(w_C-z_C) 586 #if num_min < 0.0: 587 # print 'num.min(w_C-z_C)', num_min 588 576 589 if self.use_edge_limiter: 577 590 distribute_using_edge_limiter(self) … … 1052 1065 Bed = domain.quantities['elevation'] 1053 1066 1067 1068 1054 1069 timestep = float(sys.maxint) 1055 1070 … … 1240 1255 """ 1241 1256 1242 from shallow_water_ext import balance_deep_and_shallow \1257 from anuga.shallow_water.shallow_water_ext import balance_deep_and_shallow \ 1243 1258 as balance_deep_and_shallow_c 1244 1259 -
trunk/anuga_core/source/anuga/shallow_water_balanced/run_step.py
r8175 r8178 27 27 time = strftime('%Y%m%d_%H%M%S',localtime()) 28 28 29 output_dir = 'wave_'+time 29 #output_dir = 'wave_'+time 30 output_dir = '.' 30 31 output_file = 'wave' 31 32 … … 33 34 #start_screen_catcher(output_dir+sep) 34 35 35 interactive_visualisation = False36 interactive_visualisation = True 36 37 37 38 #------------------------------------------------------------------------------ … … 73 74 def stage(x,y): 74 75 z = num.zeros_like(x) 75 num.putmask(z, x < L/2, 1.0)76 num.putmask(z, (2*L/5 < x) * (x < 3*L/5) , 1.0) 76 77 return z 77 78 … … 97 98 98 99 # Associate boundary tags with boundary objects 99 domain.set_boundary({'left': B d, 'right': Bt, 'top': Br, 'bottom': Br})100 domain.set_boundary({'left': Bt, 'right': Bt, 'top': Br, 'bottom': Br}) 100 101 101 102 -
trunk/anuga_core/source/anuga/shallow_water_balanced/run_wave.py
r7734 r8178 27 27 time = strftime('%Y%m%d_%H%M%S',localtime()) 28 28 29 output_dir = 'wave_'+time 29 #output_dir = 'wave_'+time 30 output_dir = '.' 30 31 output_file = 'wave' 31 32 -
trunk/anuga_core/source/anuga/shallow_water_balanced/run_wave_shelf.py
r7734 r8178 1 1 """Simple water flow example using ANUGA 2 2 3 Will Powers example of a simple sinusoidal wave which showed diffusive effects of 4 thefirst order and standard second order method. Problem resolved if "rk2" timestepping 5 and higher beta = 2 limiter used. Also new edge limiter with rk2 resolves problem 3 Oscillatory wave from off shore coming into shore 6 4 """ 7 5 … … 11 9 12 10 import sys 13 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross 11 import anuga 14 12 from swb_domain import * 15 13 … … 27 25 time = strftime('%Y%m%d_%H%M%S',localtime()) 28 26 29 output_dir = 'wave_shelf_'+time 27 #output_dir = 'wave_shelf_'+time 28 output_dir = '.' 30 29 output_file = 'wave_shelf' 31 30 … … 44 43 45 44 # structured mesh 46 points, vertices, boundary = rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, -W/2))45 points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, -W/2)) 47 46 48 domain = Domain(points, vertices, boundary) 47 domain = Domain(points, vertices, boundary) 49 48 50 49 domain.set_name(output_file) … … 91 90 #------------------------------------------------------------------------------ 92 91 from math import sin, pi, exp 93 Br = Reflective_boundary(domain) # Solid reflective wall94 Bt = Transmissive_boundary(domain) # Continue all values on boundary95 Bd = Dirichlet_boundary([1,0.,0.]) # Constant boundary values92 Br = anuga.Reflective_boundary(domain) # Solid reflective wall 93 Bt = anuga.Transmissive_boundary(domain) # Continue all values on boundary 94 Bd = anuga.Dirichlet_boundary([1,0.,0.]) # Constant boundary values 96 95 amplitude = 1 97 96 wave_length = 2000.0 98 Bw = Time_boundary(domain=domain, # Time dependent boundary97 Bw = anuga.Time_boundary(domain=domain, # Time dependent boundary 99 98 ## Sine wave 100 99 f=lambda t: [(amplitude*sin((1./wave_length)*t*2*pi)), 0.0, 0.0]) -
trunk/anuga_core/source/anuga/shallow_water_balanced/swb_domain.py
r8175 r8178 69 69 self.set_CFL(1.0) 70 70 self.set_beta(1.0) 71 self.quantities['height'].set_beta(1.0) 72 73 print 'Swb_Domain' 71 74 72 75 def check_integrity(self): … … 91 94 msg = 'First other quantity must be "friction"' 92 95 assert self.other_quantities[0] == 'friction', msg 96 msg = 'Second other quantity must be "x"' 97 assert self.other_quantities[1] == 'x', msg 98 msg = 'Third other quantity must be "y"' 99 assert self.other_quantities[1] == 'y', msg 100 93 101 94 102 def compute_fluxes(self): 95 #Call correct module function (either from this module or C-extension) 96 97 #compute_fluxes(self) 98 99 #return 103 """ 104 Call correct module function (either from this module or C-extension) 105 """ 106 100 107 from swb_domain_ext import compute_fluxes_c 101 108 … … 111 118 timestep = self.get_evolve_max_timestep() 112 119 113 self.flux_timestep = compute_fluxes_c(timestep, self, W, UH, VH, H, Z, U, V)114 115 #print self.flux_timestep 120 self.flux_timestep = \ 121 compute_fluxes_c(timestep, self, W, UH, VH, H, Z, U, V) 122 116 123 117 124 def distribute_to_vertices_and_edges(self): … … 134 141 """ 135 142 136 #Sww_domain.distribute_to_vertices_and_edges(self)137 #return143 from anuga.shallow_water.shallow_water_domain import \ 144 protect_against_infinitesimal_and_negative_heights as protect 138 145 139 146 #Shortcuts … … 155 162 v_C = V.centroid_values 156 163 164 num_min = num.min(w_C-z_C) 165 if num_min < 0.0: 166 print '**** num.min(w_C-z_C)', num_min 167 168 169 w_C[:] = num.maximum(w_C, z_C) 170 h_C[:] = w_C - z_C 171 172 173 assert num.min(h_C) >= 0.0 174 175 num.putmask(uh_C, h_C < 1.0e-15, 0.0) 176 num.putmask(vh_C, h_C < 1.0e-15, 0.0) 177 #num.putmask(h_C, h_C < 1.0e-15, 1.0e-16) 178 179 # Noelle has an alternative method for calculating velcities 180 # Check it out in the GPU shallow water paper 181 H0 = 1.0e-16 182 u_C[:] = uh_C/(h_C + H0/h_C) 183 v_C[:] = vh_C/(h_C + H0/h_C) 184 185 #num.putmask(h_C, h_C < 1.0e-15, 0.0) 186 187 for name in [ 'stage', 'xvelocity', 'yvelocity' ]: 188 Q = self.quantities[name] 189 if self._order_ == 1: 190 Q.extrapolate_first_order() 191 elif self._order_ == 2: 192 Q.extrapolate_second_order_and_limit_by_edge() 193 #Q.extrapolate_second_order_and_limit_by_vertex() 194 else: 195 raise Exception('Unknown order: %s' % str(self._order_)) 196 197 for name in [ 'height' ]: 198 Q = self.quantities[name] 199 if self._order_ == 1: 200 Q.extrapolate_first_order() 201 elif self._order_ == 2: 202 #Q.extrapolate_second_order_and_limit_by_edge() 203 Q.extrapolate_second_order_and_limit_by_vertex() 204 else: 205 raise Exception('Unknown order: %s' % str(self._order_)) 206 207 208 w_V = W.vertex_values 209 u_V = U.vertex_values 210 v_V = V.vertex_values 211 z_V = Z.vertex_values 212 213 h_V = H.vertex_values 214 uh_V = UH.vertex_values 215 vh_V = VH.vertex_values 216 217 218 # Update other quantities 219 #protect(self) 220 221 z_V[:] = w_V - h_V 222 uh_V[:] = u_V * h_V 223 vh_V[:] = v_V * h_V 224 225 226 num_min = num.min(h_V) 227 if num_min < 0.0: 228 print 'num.min(h_V)', num_min 229 230 231 # Compute edge values by interpolation 232 for name in ['xmomentum', 'ymomentum', 'elevation']: 233 Q = self.quantities[name] 234 Q.interpolate_from_vertices_to_edges() 235 236 237 238 239 240 def distribute_to_vertices_and_edges_h(self): 241 """Distribution from centroids to edges specific to the SWW eqn. 242 243 It will ensure that h (w-z) is always non-negative even in the 244 presence of steep bed-slopes by taking a weighted average between shallow 245 and deep cases. 246 247 In addition, all conserved quantities get distributed as per either a 248 constant (order==1) or a piecewise linear function (order==2). 249 250 251 Precondition: 252 All conserved quantities defined at centroids and bed elevation defined at 253 edges. 254 255 Postcondition 256 Evolved quantities defined at vertices and edges 257 """ 258 259 260 #Shortcuts 261 W = self.quantities['stage'] 262 UH = self.quantities['xmomentum'] 263 VH = self.quantities['ymomentum'] 264 H = self.quantities['height'] 265 Z = self.quantities['elevation'] 266 U = self.quantities['xvelocity'] 267 V = self.quantities['yvelocity'] 268 269 #Arrays 270 w_C = W.centroid_values 271 uh_C = UH.centroid_values 272 vh_C = VH.centroid_values 273 z_C = Z.centroid_values 274 h_C = H.centroid_values 275 u_C = U.centroid_values 276 v_C = V.centroid_values 277 157 278 w_C[:] = num.maximum(w_C, z_C) 158 279 … … 165 286 num.putmask(vh_C, h_C < 1.0e-15, 0.0) 166 287 num.putmask(h_C, h_C < 1.0e-15, 1.0e-16) 167 168 H0 = 1.0e-8 169 u_C[:] = uh_C/(h_C + H0/h_C) 170 v_C[:] = vh_C/(h_C + H0/h_C) 288 289 u_C[:] = uh_C/h_C 290 v_C[:] = vh_C/h_C 171 291 172 292 num.putmask(h_C, h_C < 1.0e-15, 0.0) … … 198 318 z_E[:] = w_E - h_E 199 319 200 num.putmask(h_E, h_E <= 1.0e- 15, 0.0)201 num.putmask(u_E, h_E <= 1.0e- 15, 0.0)202 num.putmask(v_E, h_E <= 1.0e- 15, 0.0)203 num.putmask(w_E, h_E <= 1.0e- 15, z_E)320 num.putmask(h_E, h_E <= 1.0e-8, 0.0) 321 num.putmask(u_E, h_E <= 1.0e-8, 0.0) 322 num.putmask(v_E, h_E <= 1.0e-8, 0.0) 323 num.putmask(w_E, h_E <= 1.0e-8, z_E) 204 324 #num.putmask(h_E, h_E <= 0.0, 0.0) 205 325 … … 243 363 244 364 245 def distribute_to_vertices_and_edges_h(self):246 """Distribution from centroids to edges specific to the SWW eqn.247 248 It will ensure that h (w-z) is always non-negative even in the249 presence of steep bed-slopes by taking a weighted average between shallow250 and deep cases.251 252 In addition, all conserved quantities get distributed as per either a253 constant (order==1) or a piecewise linear function (order==2).254 255 256 Precondition:257 All conserved quantities defined at centroids and bed elevation defined at258 edges.259 260 Postcondition261 Evolved quantities defined at vertices and edges262 """263 264 265 #Shortcuts266 W = self.quantities['stage']267 UH = self.quantities['xmomentum']268 VH = self.quantities['ymomentum']269 H = self.quantities['height']270 Z = self.quantities['elevation']271 U = self.quantities['xvelocity']272 V = self.quantities['yvelocity']273 274 #Arrays275 w_C = W.centroid_values276 uh_C = UH.centroid_values277 vh_C = VH.centroid_values278 z_C = Z.centroid_values279 h_C = H.centroid_values280 u_C = U.centroid_values281 v_C = V.centroid_values282 283 w_C[:] = num.maximum(w_C, z_C)284 285 h_C[:] = w_C - z_C286 287 288 assert num.min(h_C) >= 0289 290 num.putmask(uh_C, h_C < 1.0e-15, 0.0)291 num.putmask(vh_C, h_C < 1.0e-15, 0.0)292 num.putmask(h_C, h_C < 1.0e-15, 1.0e-16)293 294 u_C[:] = uh_C/h_C295 v_C[:] = vh_C/h_C296 297 num.putmask(h_C, h_C < 1.0e-15, 0.0)298 299 for name in [ 'stage', 'height', 'xvelocity', 'yvelocity' ]:300 Q = self.quantities[name]301 if self._order_ == 1:302 Q.extrapolate_first_order()303 elif self._order_ == 2:304 Q.extrapolate_second_order_and_limit_by_edge()305 #Q.extrapolate_second_order_and_limit_by_vertex()306 else:307 raise Exception('Unknown order: %s' % str(self._order_))308 309 310 w_E = W.edge_values311 uh_E = UH.edge_values312 vh_E = VH.edge_values313 h_E = H.edge_values314 z_E = Z.edge_values315 u_E = U.edge_values316 v_E = V.edge_values317 318 319 #minh_E = num.min(h_E)320 #msg = 'min h_E = %g ' % minh_E321 #assert minh_E >= -1.0e-15, msg322 323 z_E[:] = w_E - h_E324 325 num.putmask(h_E, h_E <= 1.0e-8, 0.0)326 num.putmask(u_E, h_E <= 1.0e-8, 0.0)327 num.putmask(v_E, h_E <= 1.0e-8, 0.0)328 num.putmask(w_E, h_E <= 1.0e-8, z_E)329 #num.putmask(h_E, h_E <= 0.0, 0.0)330 331 uh_E[:] = u_E * h_E332 vh_E[:] = v_E * h_E333 334 """335 print '=========================================================='336 print 'Time ', self.get_time()337 print h_E338 print uh_E339 print vh_E340 """341 342 # Compute vertex values by interpolation343 for name in self.evolved_quantities:344 Q = self.quantities[name]345 Q.interpolate_from_edges_to_vertices()346 347 348 w_V = W.vertex_values349 uh_V = UH.vertex_values350 vh_V = VH.vertex_values351 z_V = Z.vertex_values352 h_V = H.vertex_values353 u_V = U.vertex_values354 v_V = V.vertex_values355 356 357 #w_V[:] = z_V + h_V358 359 #num.putmask(u_V, h_V <= 0.0, 0.0)360 #num.putmask(v_V, h_V <= 0.0, 0.0)361 #num.putmask(w_V, h_V <= 0.0, z_V)362 #num.putmask(h_V, h_V <= 0.0, 0.0)363 364 uh_V[:] = u_V * h_V365 vh_V[:] = v_V * h_V366 367 368 369 370 365 def conserved_values_to_evolved_values(self, q_cons, q_evol): 371 366 """Mapping between conserved quantities and the evolved quantities. -
trunk/anuga_core/source/anuga/shallow_water_balanced/swb_domain_ext.c
r7734 r8178 1 1 // Python - C extension module for shallow_water.py 2 2 // 3 // To compile (Python2. 3):4 // gcc -c swb_domain_ext.c -I/usr/include/python2. 3-o domain_ext.o -Wall -O3 // To compile (Python2.6): 4 // gcc -c swb_domain_ext.c -I/usr/include/python2.6 -o domain_ext.o -Wall -O 5 5 // gcc -shared swb_domain_ext.o -o swb_domain_ext.so 6 6 // … … 150 150 denom = s_max - s_min; 151 151 152 memset(edgeflux, 0, 3*sizeof(double)); 152 edgeflux[0] = 0.0; 153 edgeflux[1] = 0.0; 154 edgeflux[2] = 0.0; 153 155 154 156 if (denom < epsilon) … … 170 172 171 173 // Balance the pressure term 174 // FIXME SR: I think we need to add a term which uses simpson's rule. 172 175 edgeflux[1] += 0.5*g*h_inside*h_inside - 0.5*g*h_inside_star*h_inside_star; 173 176
Note: See TracChangeset
for help on using the changeset viewer.
