Changeset 7616
- Timestamp:
- Feb 5, 2010, 5:24:39 PM (14 years ago)
- Location:
- anuga_core/source/anuga
- Files:
-
- 5 edited
-
abstract_2d_finite_volumes/quantity_ext.c (modified) (5 diffs)
-
compile_all.py (modified) (1 diff)
-
shallow_water_balanced/swb_domain.py (modified) (8 diffs)
-
shallow_water_balanced/test_swb_balance.py (modified) (5 diffs)
-
shallow_water_balanced/test_swb_boundary_condition.py (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c
r7573 r7616 362 362 363 363 sign = 0.0; 364 if (qmin > 0 ) {364 if (qmin > 0.0) { 365 365 sign = 1.0; 366 366 } else if (qmax < 0) { … … 373 373 dqa[i] = dq; //Save dq for use in updating vertex values 374 374 375 // FIXME SR 20091125 This caused problems in shallow_water_balanced 376 // commenting out problem 375 377 376 // Just limit non boundary edges so that we can reconstruct a linear function 378 if (neighbours[k3+i] >= 0) { 377 // FIXME Problem with stability on edges 378 //if (neighbours[k3+i] >= 0) { 379 379 r = 1.0; 380 380 … … 383 383 384 384 phi = min( min(r*beta, 1.0), phi); 385 } 386 387 if (neighbours[k3+i] < 0) { 388 r = 1.0; 389 390 if (dq > 0.0 && sign == -1.0 ) r = (0.0 - qc)/dq; 391 if (dq < 0.0 && sign == 1.0 ) r = (0.0 - qc)/dq; 385 // } 386 387 // 388 /* if (neighbours[k3+i] < 0) { */ 389 /* r = 1.0; */ 390 391 /* if (dq > 0.0 && (sign == -1.0 || sign == 0.0 )) r = (0.0 - qc)/dq; */ 392 /* if (dq < 0.0 && (sign == 1.0 || sign == 0.0 )) r = (0.0 - qc)/dq; */ 392 393 393 phi = min( min(r*beta, 1.0), phi); 394 } 394 /* phi = min( min(r*beta, 1.0), phi); */ 395 /* } */ 395 396 396 397 } … … 1194 1195 if (!PyArg_ParseTuple(args, "O",&quantity)) { 1195 1196 PyErr_SetString(PyExc_RuntimeError, 1196 "quantity_ext.c: extrapolate_second_order_and_limit could not parse input");1197 "quantity_ext.c: extrapolate_second_order_and_limit_by_edge could not parse input"); 1197 1198 return NULL; 1198 1199 } … … 1201 1202 if (!domain) { 1202 1203 PyErr_SetString(PyExc_RuntimeError, 1203 "quantity_ext.c: extrapolate_second_order_and_limit could not obtain domain object from quantity");1204 "quantity_ext.c: extrapolate_second_order_and_limit_by_edge could not obtain domain object from quantity"); 1204 1205 return NULL; 1205 1206 } -
anuga_core/source/anuga/compile_all.py
r6718 r7616 21 21 22 22 os.chdir('..') 23 os.chdir('shallow_water_balanced') 24 execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py') 25 26 os.chdir('..') 23 27 os.chdir('mesh_engine') 24 28 execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py') -
anuga_core/source/anuga/shallow_water_balanced/swb_domain.py
r7575 r7616 80 80 # set some defaults 81 81 #--------------------- 82 self.set_timestepping_method( 1)82 self.set_timestepping_method(2) 83 83 self.set_default_order(2) 84 84 self.set_new_mannings_function(True) 85 85 self.set_centroid_transmissive_bc(True) 86 self.set_CFL(1.0) 87 self.set_beta(1.0) 86 88 87 89 ## … … 114 116 #Call correct module function (either from this module or C-extension) 115 117 118 #compute_fluxes(self) 119 120 #return 116 121 from swb_domain_ext import compute_fluxes_c 117 122 … … 129 134 self.flux_timestep = compute_fluxes_c(timestep, self, W, UH, VH, H, Z, U, V) 130 135 136 #print self.flux_timestep 137 131 138 ## 132 139 # @brief … … 150 157 """ 151 158 152 159 #Sww_domain.distribute_to_vertices_and_edges(self) 160 #return 161 153 162 #Shortcuts 154 163 W = self.quantities['stage'] … … 178 187 num.putmask(uh_C, h_C < 1.0e-15, 0.0) 179 188 num.putmask(vh_C, h_C < 1.0e-15, 0.0) 180 num.putmask(h_C, h_C < 1.0e-15, 1.0e-15) 181 182 u_C[:] = uh_C/h_C 183 v_C[:] = vh_C/h_C 184 185 for name in [ 'height', 'elevation', 'xvelocity', 'yvelocity' ]: 189 num.putmask(h_C, h_C < 1.0e-15, 1.0e-16) 190 191 H0 = 1.0e-8 192 u_C[:] = uh_C/(h_C + H0/h_C) 193 v_C[:] = vh_C/(h_C + H0/h_C) 194 195 num.putmask(h_C, h_C < 1.0e-15, 0.0) 196 197 for name in [ 'stage', 'height', 'xvelocity', 'yvelocity' ]: 186 198 Q = self.quantities[name] 187 199 if self._order_ == 1: … … 189 201 elif self._order_ == 2: 190 202 Q.extrapolate_second_order_and_limit_by_edge() 203 #Q.extrapolate_second_order_and_limit_by_vertex() 191 204 else: 192 205 raise 'Unknown order' … … 202 215 203 216 204 w_E[:] = z_E + h_E 205 206 #num.putmask(u_E, h_temp <= 0.0, 0.0) 207 #num.putmask(v_E, h_temp <= 0.0, 0.0) 208 #num.putmask(w_E, h_temp <= 0.0, z_E+h_E) 217 #minh_E = num.min(h_E) 218 #msg = 'min h_E = %g ' % minh_E 219 #assert minh_E >= -1.0e-15, msg 220 221 z_E[:] = w_E - h_E 222 223 num.putmask(h_E, h_E <= 1.0e-15, 0.0) 224 num.putmask(u_E, h_E <= 1.0e-15, 0.0) 225 num.putmask(v_E, h_E <= 1.0e-15, 0.0) 226 num.putmask(w_E, h_E <= 1.0e-15, z_E) 209 227 #num.putmask(h_E, h_E <= 0.0, 0.0) 210 228 … … 246 264 247 265 266 267 268 ## 269 # @brief 270 def distribute_to_vertices_and_edges_h(self): 271 """Distribution from centroids to edges specific to the SWW eqn. 272 273 It will ensure that h (w-z) is always non-negative even in the 274 presence of steep bed-slopes by taking a weighted average between shallow 275 and deep cases. 276 277 In addition, all conserved quantities get distributed as per either a 278 constant (order==1) or a piecewise linear function (order==2). 279 280 281 Precondition: 282 All conserved quantities defined at centroids and bed elevation defined at 283 edges. 284 285 Postcondition 286 Evolved quantities defined at vertices and edges 287 """ 288 289 290 #Shortcuts 291 W = self.quantities['stage'] 292 UH = self.quantities['xmomentum'] 293 VH = self.quantities['ymomentum'] 294 H = self.quantities['height'] 295 Z = self.quantities['elevation'] 296 U = self.quantities['xvelocity'] 297 V = self.quantities['yvelocity'] 298 299 #Arrays 300 w_C = W.centroid_values 301 uh_C = UH.centroid_values 302 vh_C = VH.centroid_values 303 z_C = Z.centroid_values 304 h_C = H.centroid_values 305 u_C = U.centroid_values 306 v_C = V.centroid_values 307 308 w_C[:] = num.maximum(w_C, z_C) 309 310 h_C[:] = w_C - z_C 311 312 313 assert num.min(h_C) >= 0 314 315 num.putmask(uh_C, h_C < 1.0e-15, 0.0) 316 num.putmask(vh_C, h_C < 1.0e-15, 0.0) 317 num.putmask(h_C, h_C < 1.0e-15, 1.0e-16) 318 319 u_C[:] = uh_C/h_C 320 v_C[:] = vh_C/h_C 321 322 num.putmask(h_C, h_C < 1.0e-15, 0.0) 323 324 for name in [ 'stage', 'height', 'xvelocity', 'yvelocity' ]: 325 Q = self.quantities[name] 326 if self._order_ == 1: 327 Q.extrapolate_first_order() 328 elif self._order_ == 2: 329 Q.extrapolate_second_order_and_limit_by_edge() 330 #Q.extrapolate_second_order_and_limit_by_vertex() 331 else: 332 raise 'Unknown order' 333 334 335 w_E = W.edge_values 336 uh_E = UH.edge_values 337 vh_E = VH.edge_values 338 h_E = H.edge_values 339 z_E = Z.edge_values 340 u_E = U.edge_values 341 v_E = V.edge_values 342 343 344 #minh_E = num.min(h_E) 345 #msg = 'min h_E = %g ' % minh_E 346 #assert minh_E >= -1.0e-15, msg 347 348 z_E[:] = w_E - h_E 349 350 num.putmask(h_E, h_E <= 1.0e-8, 0.0) 351 num.putmask(u_E, h_E <= 1.0e-8, 0.0) 352 num.putmask(v_E, h_E <= 1.0e-8, 0.0) 353 num.putmask(w_E, h_E <= 1.0e-8, z_E) 354 #num.putmask(h_E, h_E <= 0.0, 0.0) 355 356 uh_E[:] = u_E * h_E 357 vh_E[:] = v_E * h_E 358 359 """ 360 print '==========================================================' 361 print 'Time ', self.get_time() 362 print h_E 363 print uh_E 364 print vh_E 365 """ 366 367 # Compute vertex values by interpolation 368 for name in self.evolved_quantities: 369 Q = self.quantities[name] 370 Q.interpolate_from_edges_to_vertices() 371 372 373 w_V = W.vertex_values 374 uh_V = UH.vertex_values 375 vh_V = VH.vertex_values 376 z_V = Z.vertex_values 377 h_V = H.vertex_values 378 u_V = U.vertex_values 379 v_V = V.vertex_values 380 381 382 #w_V[:] = z_V + h_V 383 384 #num.putmask(u_V, h_V <= 0.0, 0.0) 385 #num.putmask(v_V, h_V <= 0.0, 0.0) 386 #num.putmask(w_V, h_V <= 0.0, z_V) 387 #num.putmask(h_V, h_V <= 0.0, 0.0) 388 389 uh_V[:] = u_V * h_V 390 vh_V[:] = v_V * h_V 391 392 393 394 248 395 ## 249 396 # @brief Code to let us use old shallow water domain BCs -
anuga_core/source/anuga/shallow_water_balanced/test_swb_balance.py
r7573 r7616 63 63 64 64 # Limit 65 domain.tight_slope_limiters = 066 65 domain.distribute_to_vertices_and_edges() 67 66 … … 72 71 73 72 # Now try with a non-flat bed - closely hugging initial stage in places 74 # This will create alphas in the range [0, 0.478260, 1]75 73 domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) 76 74 domain.set_quantity('elevation', [[0,0,0], … … 79 77 [0,2,4]]) 80 78 stage = domain.quantities['stage'] 79 elevation = domain.quantities['elevation'] 80 height = domain.quantities['height'] 81 81 82 82 ref_centroid_values = copy.copy(stage.centroid_values[:]) # Copy … … 84 84 85 85 # Limit 86 domain.tight_slope_limiters = 087 86 domain.distribute_to_vertices_and_edges() 88 87 … … 98 97 # Check actual results 99 98 assert num.allclose(stage.vertex_values, 100 [[ 2., 2., 2. ], 101 [ 1.93333333, 2.03333333, 6.03333333], 102 [ 8.4, 3.8, 3.8 ], 103 [ 3.33333333, 5.33333333, 7.33333333]]) or \ 104 num.allclose(stage.vertex_values, 105 [[ 1.06666667, 3.86666667, 1.06666667], 106 [ 1.93333333, 2.03333333, 6.03333333], 107 [ 5.46666667, 3.06666667, 7.46666667], 108 [ 6.53333333, 4.69333333, 4.77333333]]) 99 [[ 2.66666667, 0.66666667, 2.66666667], 100 [ 3.33333333, 3.33333333, 3.33333333], 101 [ 3.73333333, 4.93333333, 7.33333333], 102 [ 7.33333333, 4.93333333, 3.73333333]]) 109 103 110 104 -
anuga_core/source/anuga/shallow_water_balanced/test_swb_boundary_condition.py
r7573 r7616 871 871 872 872 if __name__ == "__main__": 873 suite = unittest.makeSuite(Test_swb_boundary_condition, 'test ')873 suite = unittest.makeSuite(Test_swb_boundary_condition, 'test_boundary_condition_time') 874 874 runner = unittest.TextTestRunner(verbosity=1) 875 875 runner.run(suite)
Note: See TracChangeset
for help on using the changeset viewer.
