Changeset 7519
- Timestamp:
- Sep 20, 2009, 6:35:14 PM (14 years ago)
- Location:
- anuga_core/source
- Files:
-
- 16 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py
r7498 r7519 65 65 boundary=None, 66 66 conserved_quantities=None, 67 evolved_quantities=None, 67 68 other_quantities=None, 68 69 tagged_elements=None, … … 90 91 conserved_quantities: List of quantity names entering the 91 92 conservation equations 93 evolved_quantities: List of all quantities that evolve 92 94 other_quantities: List of other quantity names 93 95 … … 154 156 self.conserved_quantities = conserved_quantities 155 157 158 if evolved_quantities is None: 159 self.evolved_quantities = self.conserved_quantities 160 else: 161 self.evolved_quantities = evolved_quantities 162 156 163 # List of other quantity names 157 164 if other_quantities is None: … … 160 167 self.other_quantities = other_quantities 161 168 169 # Test that conserved_quantities are stored in the first entries of evolved_quantities 170 for i, quantity in enumerate(self.conserved_quantities): 171 msg = 'The conserved quantities must be the first entries of evolved_quantities' 172 assert quantity == self.evolved_quantities[i], msg 173 174 162 175 # Build dictionary of Quantity instances keyed by quantity names 163 176 self.quantities = {} 164 177 165 for name in self. conserved_quantities:178 for name in self.evolved_quantities: 166 179 self.quantities[name] = Quantity(self) 167 180 for name in self.other_quantities: … … 325 338 return self.mesh.get_number_of_nodes(*args, **kwargs) 326 339 340 def get_number_of_triangles(self, *args, **kwargs): 341 return self.mesh.get_number_of_triangles(*args, **kwargs) 342 327 343 def get_normal(self, *args, **kwargs): 328 344 return self.mesh.get_normal(*args, **kwargs) … … 415 431 416 432 return q 433 434 ## 435 # @brief Get evolved quantities for a volume. 436 # @param vol_id ID of the volume we want the conserved quantities for. 437 # @param vertex If specified, use as index for edge values. 438 # @param edge If specified, use as index for edge values. 439 # @return Vector of conserved quantities. 440 # @note If neither 'vertex' or 'edge' specified, use centroid values. 441 # @note If both 'vertex' and 'edge' specified, raise exception. 442 def get_evolved_quantities(self, vol_id, 443 vertex=None, 444 edge=None): 445 """Get evolved quantities at volume vol_id. 446 447 If vertex is specified use it as index for vertex values 448 If edge is specified use it as index for edge values 449 If neither are specified use centroid values 450 If both are specified an exeception is raised 451 452 Return value: Vector of length == number_of_conserved quantities 453 """ 454 455 if not (vertex is None or edge is None): 456 msg = 'Values for both vertex and edge was specified.' 457 msg += 'Only one (or none) is allowed.' 458 raise Exception, msg 459 460 q = num.zeros(len(self.evolved_quantities), num.float) 461 462 for i, name in enumerate(self.evolved_quantities): 463 Q = self.quantities[name] 464 if vertex is not None: 465 q[i] = Q.vertex_values[vol_id, vertex] 466 elif edge is not None: 467 q[i] = Q.edge_values[vol_id, edge] 468 else: 469 q[i] = Q.centroid_values[vol_id] 470 471 return q 472 473 ## 474 # @brief 475 # @param flag 476 def set_CFL(self, cfl=1.0): 477 """Set CFL parameter, warn if greater than 1.0 478 """ 479 if cfl > 1.0: 480 self.CFL = cfl 481 log.warn('Setting CFL > 1.0') 482 483 assert cfl > 0.0 484 self.CFL = cfl 485 486 417 487 418 488 ## … … 876 946 assert quantity in self.quantities, msg 877 947 878 ##assert hasattr(self, 'boundary_objects') 948 949 for i, quantity in enumerate(self.conserved_quantities): 950 msg = 'Conserved quantities must be the first entries of evolved_quantities' 951 assert quantity == self.evolved_quantities[i], msg 952 879 953 880 954 ## … … 1654 1728 1655 1729 ## 1656 # @brief Backup conserved quantities .1730 # @brief Backup conserved quantities 1657 1731 def backup_conserved_quantities(self): 1658 N = len(self) # Number_of_triangles1659 1732 1660 1733 # Backup conserved_quantities centroid values … … 1664 1737 1665 1738 ## 1666 # @brief ??1667 # @param a ??1668 # @param b ??1739 # @brief Combines current C and saved centroid values S as C = aC + bS 1740 # @param a factor in combination 1741 # @param b factor in combination 1669 1742 def saxpy_conserved_quantities(self, a, b): 1670 N = len(self) #number_of_triangles1671 1743 1672 1744 # Backup conserved_quantities centroid values … … 1675 1747 Q.saxpy_centroid_values(a, b) 1676 1748 1749 1750 1751 1752 ## 1753 # @brief Mapping between conserved quantites and evolved quantities 1754 # @param q_cons array of conserved quantity values 1755 # @param q_evolved array of current evolved quantity values 1756 def conserved_to_evolved(self, q_cons, q_evolved): 1757 """Needs to be overridden by Domain subclass 1758 """ 1759 1760 if len(q_cons) == len(q_evolved): 1761 q_evolved[:] = q_cons 1762 else: 1763 msg = 'Method conserved_to_evolved must be overridden by Domain subclass' 1764 raise Exception, msg 1765 1677 1766 ## 1678 1767 # @brief Update boundary values for all conserved quantities. … … 1692 1781 log.critical('WARNING: Ignored boundary segment (None)') 1693 1782 else: 1694 q = B.evaluate(vol_id, edge_id) 1695 1696 for j, name in enumerate(self.conserved_quantities): 1783 q_cons = B.evaluate(vol_id, edge_id) 1784 1785 if len(q_cons) == len(self.evolved_quantities): 1786 # conserved and evolved quantities are the same 1787 q_evol = q_cons 1788 elif len(q_cons) == len(self.conserved_quantities): 1789 # boundary just returns conserved quantities 1790 # Need to calculate all the evolved quantities 1791 # Use default conversion 1792 1793 q_evol = self.get_evolved_quantities(vol_id, edge = edge_id) 1794 1795 self.conserved_to_evolved(q_cons, q_evol) 1796 else: 1797 msg = 'Boundary must return array of either conserved or evolved quantities' 1798 raise Exception, msg 1799 1800 for j, name in enumerate(self.evolved_quantities): 1697 1801 Q = self.quantities[name] 1698 Q.boundary_values[i] = q [j]1802 Q.boundary_values[i] = q_evol[j] 1699 1803 1700 1804 ## -
anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py
r7498 r7519 235 235 236 236 237 def get_number_of_triangles(self): 238 return self.number_of_triangles 239 240 237 241 def get_number_of_nodes(self): 238 242 return self.number_of_nodes 243 239 244 240 245 def get_nodes(self, absolute=False): -
anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c
r7276 r7519 712 712 713 713 714 // Explicit updates 715 for (k=0; k<N; k++) { 716 centroid_values[k] += timestep*explicit_update[k]; 717 } 718 719 720 714 721 // Semi implicit updates 715 722 for (k=0; k<N; k++) { 716 723 denominator = 1.0 - timestep*semi_implicit_update[k]; 717 if (denominator == 0.0) {724 if (denominator <= 0.0) { 718 725 return -1; 719 726 } else { … … 724 731 725 732 726 // Explicit updates 727 for (k=0; k<N; k++) { 728 centroid_values[k] += timestep*explicit_update[k]; 729 } 733 730 734 731 735 -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py
r7317 r7519 53 53 54 54 conserved_quantities = ['stage', 'xmomentum', 'ymomentum'] 55 evolved_quantities = ['stage', 'xmomentum', 'ymomentum', 'xvelocity'] 56 55 57 other_quantities = ['elevation', 'friction'] 56 58 57 59 domain = Domain(points, vertices, None, 58 conserved_quantities, other_quantities)60 conserved_quantities, evolved_quantities, other_quantities) 59 61 domain.check_integrity() 60 62 … … 65 67 assert num.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.) 66 68 69 70 71 def test_CFL(self): 72 a = [0.0, 0.0] 73 b = [0.0, 2.0] 74 c = [2.0,0.0] 75 d = [0.0, 4.0] 76 e = [2.0, 2.0] 77 f = [4.0,0.0] 78 79 points = [a, b, c, d, e, f] 80 #bac, bce, ecf, dbe, daf, dae 81 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] 82 83 conserved_quantities = ['stage', 'xmomentum', 'ymomentum'] 84 evolved_quantities = ['stage', 'xmomentum', 'ymomentum', 'xvelocity'] 85 86 other_quantities = ['elevation', 'friction'] 87 88 domain = Domain(points, vertices, None, 89 conserved_quantities, evolved_quantities, other_quantities) 90 91 try: 92 domain.set_CFL(-0.1) 93 except: 94 pass 95 else: 96 msg = 'Should have caught a negative cfl' 97 raise Exception, msg 98 99 100 101 try: 102 domain.set_CFL(2.0) 103 except: 104 pass 105 else: 106 msg = 'Should have warned of cfl>1.0' 107 raise Exception, msg 108 109 assert domain.CFL == 2.0 110 67 111 68 112 def test_conserved_quantities(self): … … 511 555 512 556 557 def test_conserved_evolved_boundary_conditions(self): 558 559 a = [0.0, 0.0] 560 b = [0.0, 2.0] 561 c = [2.0,0.0] 562 d = [0.0, 4.0] 563 e = [2.0, 2.0] 564 f = [4.0,0.0] 565 566 points = [a, b, c, d, e, f] 567 #bac, bce, ecf, dbe 568 vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 569 boundary = { (0, 0): 'First', 570 (0, 2): 'First', 571 (2, 0): 'Second', 572 (2, 1): 'Second', 573 (3, 1): 'Second', 574 (3, 2): 'Second'} 575 576 577 578 try: 579 domain = Domain(points, vertices, boundary, 580 conserved_quantities = ['stage', 'xmomentum', 'ymomentum'], 581 evolved_quantities =\ 582 ['stage', 'xmomentum', 'xvelocity', 'ymomentum', 'yvelocity']) 583 except: 584 pass 585 else: 586 msg = 'Should have caught the evolved quantities not being in order' 587 raise Exception, msg 588 589 590 domain = Domain(points, vertices, boundary, 591 conserved_quantities = ['stage', 'xmomentum', 'ymomentum'], 592 evolved_quantities =\ 593 ['stage', 'xmomentum', 'ymomentum', 'xvelocity', 'yvelocity']) 594 595 596 domain.set_quantity('stage', [[1,2,3], [5,5,5], 597 [0,0,9], [6, -3, 3]]) 598 599 600 domain.set_boundary( {'First': Dirichlet_boundary([5,2,1,4,6]), 601 'Second': Transmissive_boundary(domain)} ) 602 603 try: 604 domain.update_boundary() 605 except: 606 pass 607 else: 608 msg = 'Should have caught the lack of conserved_to_evolved member function' 609 raise Exception, msg 610 611 612 def conserved_to_evolved(q_cons, q_evol): 613 614 q_evol[0:3] = q_cons 615 q_evol[3] = q_cons[1]/q_cons[0] 616 q_evol[4] = q_cons[2]/q_cons[0] 617 618 domain.conserved_to_evolved = conserved_to_evolved 619 620 domain.update_boundary() 621 622 623 assert domain.quantities['stage'].boundary_values[0] == 5. #Dirichlet 624 assert domain.quantities['stage'].boundary_values[1] == 5. #Dirichlet 625 assert domain.quantities['xvelocity'].boundary_values[0] == 4. #Dirichlet 626 assert domain.quantities['yvelocity'].boundary_values[1] == 6. #Dirichlet 627 628 q_cons = domain.get_conserved_quantities(2, edge=0) #Transmissive 629 assert domain.quantities['stage' ].boundary_values[2] == q_cons[0] 630 assert domain.quantities['xmomentum'].boundary_values[2] == q_cons[1] 631 assert domain.quantities['ymomentum'].boundary_values[2] == q_cons[2] 632 assert domain.quantities['xvelocity'].boundary_values[2] == q_cons[1]/q_cons[0] 633 assert domain.quantities['yvelocity'].boundary_values[2] == q_cons[2]/q_cons[0] 634 635 q_cons = domain.get_conserved_quantities(2, edge=1) #Transmissive 636 assert domain.quantities['stage' ].boundary_values[3] == q_cons[0] 637 assert domain.quantities['xmomentum'].boundary_values[3] == q_cons[1] 638 assert domain.quantities['ymomentum'].boundary_values[3] == q_cons[2] 639 assert domain.quantities['xvelocity'].boundary_values[3] == q_cons[1]/q_cons[0] 640 assert domain.quantities['yvelocity'].boundary_values[3] == q_cons[2]/q_cons[0] 641 642 643 q_cons = domain.get_conserved_quantities(3, edge=1) #Transmissive 644 assert domain.quantities['stage' ].boundary_values[4] == q_cons[0] 645 assert domain.quantities['xmomentum'].boundary_values[4] == q_cons[1] 646 assert domain.quantities['ymomentum'].boundary_values[4] == q_cons[2] 647 assert domain.quantities['xvelocity'].boundary_values[4] == q_cons[1]/q_cons[0] 648 assert domain.quantities['yvelocity'].boundary_values[4] == q_cons[2]/q_cons[0] 649 650 651 q_cons = domain.get_conserved_quantities(3, edge=2) #Transmissive 652 assert domain.quantities['stage' ].boundary_values[5] == q_cons[0] 653 assert domain.quantities['xmomentum'].boundary_values[5] == q_cons[1] 654 assert domain.quantities['ymomentum'].boundary_values[5] == q_cons[2] 655 assert domain.quantities['xvelocity'].boundary_values[5] == q_cons[1]/q_cons[0] 656 assert domain.quantities['yvelocity'].boundary_values[5] == q_cons[2]/q_cons[0] 657 658 513 659 def test_distribute_first_order(self): 514 660 """Domain implements a default first order gradient limiter … … 598 744 #Set explicit_update 599 745 746 600 747 for name in domain.conserved_quantities: 601 748 domain.quantities[name].explicit_update = num.array([4.,3.,2.,1.]) … … 614 761 615 762 x = num.array([1., 2., 3., 4.]) 763 x += domain.timestep*num.array( [4,3,2,1] ) 616 764 x /= denom 617 x += domain.timestep*num.array( [4,3,2,1] ) 765 618 766 619 767 for name in domain.conserved_quantities: … … 906 1054 907 1055 if __name__ == "__main__": 908 suite = unittest.makeSuite(Test_Domain,'test ')1056 suite = unittest.makeSuite(Test_Domain,'test_conserved_evolved') 909 1057 runner = unittest.TextTestRunner() 910 1058 runner.run(suite) -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py
r7505 r7519 68 68 69 69 domain.conserved_quantities = ['stage', 'ymomentum'] 70 domain.evolved_quantities = ['stage', 'ymomentum'] 70 71 domain.quantities['stage'] =\ 71 72 Quantity(domain, [[1,2,3], [5,5,5], … … 133 134 134 135 domain.conserved_quantities = ['stage', 'ymomentum'] 136 domain.evolved_quantities = ['stage', 'ymomentum'] 135 137 domain.quantities['stage'] =\ 136 138 Quantity(domain, [[1,2,3], [5,5,5], … … 178 180 179 181 180 q = T.evaluate(1, 1) #Vol= 0, edge=2182 q = T.evaluate(1, 1) #Vol=1, edge=1 181 183 assert num.allclose(q, domain.get_edge_midpoint_coordinate(1,1)) 182 184 … … 204 206 205 207 domain.conserved_quantities = ['stage', 'ymomentum'] 208 domain.evolved_quantities = ['stage', 'ymomentum'] 206 209 domain.quantities['stage'] =\ 207 210 Quantity(domain, [[1,2,3], [5,5,5], … … 362 365 domain = Domain(points, elements) 363 366 domain.conserved_quantities = ['stage', 'xmomentum', 'ymomentum'] 367 domain.evolved_quantities = ['stage', 'xmomentum', 'ymomentum'] 364 368 domain.quantities['stage'] =\ 365 369 Quantity(domain, [[1,2,3], [5,5,5], -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_ghost.py
r7276 r7519 35 35 36 36 domain = Domain(points, vertices, None, 37 conserved_quantities, other_quantities)37 conserved_quantities, None, other_quantities) 38 38 domain.check_integrity() 39 39 -
anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py
r7276 r7519 1867 1867 1868 1868 x = num.array([1., 2., 3., 4.]) 1869 x += timestep*num.array( [4.0, 3.0, 2.0, 1.0] ) 1869 1870 x /= denom 1870 x += timestep*num.array( [4.0, 3.0, 2.0, 1.0] ) 1871 1871 1872 1872 1873 assert num.allclose( quantity.centroid_values, x) -
anuga_core/source/anuga/interface.py
r7452 r7519 18 18 19 19 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Time_boundary 20 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Time_space_boundary 20 21 from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Transmissive_boundary 21 22 … … 39 40 points, vertices, boundary = rectangular_cross(*args, **kwargs) 40 41 return Domain(points, vertices, boundary) 41 42 42 43 #---------------------------- 43 44 # Create domain from file -
anuga_core/source/anuga/shallow_water/__init__.py
r7326 r7519 16 16 Transmissive_n_momentum_zero_t_momentum_set_stage_boundary,\ 17 17 AWI_boundary 18 18 19 #from shallow_water_balanced_domain import Swb_domain 19 20 20 21 -
anuga_core/source/anuga/shallow_water/data_manager.py
r7509 r7519 2484 2484 #P = interp.mesh.get_boundary_polygon() 2485 2485 #inside_indices = inside_polygon(grid_points, P) 2486 2487 # change printoptions so that a long string of zeros in not 2488 # summarized as [0.0, 0.0, 0.0, ... 0.0, 0.0, 0.0] 2489 printoptions = num.get_printoptions() 2490 num.set_printoptions(threshold=sys.maxint) 2491 2486 2492 for i in range(nrows): 2487 2493 if verbose and i % ((nrows+10)/10) == 0: … … 2492 2498 slice = grid_values[base_index:base_index+ncols] 2493 2499 #s = array2string(slice, max_line_width=sys.maxint) 2500 2494 2501 s = num.array2string(slice, max_line_width=sys.maxint, 2495 2502 precision=number_of_decimal_places) 2496 2503 ascid.write(s[1:-1] + '\n') 2497 2504 2505 num.set_printoptions(threshold=printoptions['threshold']) 2506 2498 2507 #Close 2499 2508 ascid.close() … … 5882 5891 # brief Instantiate the SWW writer class. 5883 5892 def __init__(self, static_quantities, dynamic_quantities): 5884 """Initialise SWW writerwith two list af quantity names:5893 """Initialise Write_sww with two list af quantity names: 5885 5894 5886 5895 static_quantities (e.g. elevation or friction): -
anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r7452 r7519 147 147 # @param use_cache 148 148 # @param verbose 149 # @param evolved_quantities 149 150 # @param full_send_dict 150 151 # @param ghost_recv_dict … … 163 164 use_cache=False, 164 165 verbose=False, 166 evolved_quantities = None, 167 other_quantities = None, 165 168 full_send_dict=None, 166 169 ghost_recv_dict=None, … … 171 174 172 175 # Define quantities for the shallow_water domain 173 conserved_quantities = ['stage', 'xmomentum', 'ymomentum'] 174 other_quantities = ['elevation', 'friction'] 176 conserved_quantities = ['stage', 'xmomentum', 'ymomentum'] 177 178 if other_quantities == None: 179 other_quantities = ['elevation', 'friction'] 175 180 176 181 Generic_Domain.__init__(self, … … 179 184 boundary, 180 185 conserved_quantities, 186 evolved_quantities, 181 187 other_quantities, 182 188 tagged_elements, … … 208 214 self.optimise_dry_cells = optimise_dry_cells 209 215 216 self.use_new_mannings = False 210 217 self.forcing_terms.append(manning_friction_implicit) 211 218 self.forcing_terms.append(gravity) … … 225 232 self.use_centroid_velocities = use_centroid_velocities 226 233 234 235 ## 236 # @brief 237 # @param flag 238 def set_new_mannings_function(self, flag=True): 239 """Cludge to allow unit test to pass, but to 240 also introduce new mannings friction function 241 which takes into account the slope of the bed. 242 The flag is tested in the python wrapper 243 mannings_friction_implicit 244 """ 245 if flag: 246 self.use_new_mannings = True 247 else: 248 self.use_new_mannings = False 249 250 251 ## 252 # @brief 253 # @param flag 254 def set_use_edge_limiter(self, flag=True): 255 """Cludge to allow unit test to pass, but to 256 also introduce new edge limiting. The flag is 257 tested in distribute_to_vertices_and_edges 258 """ 259 if flag: 260 self.use_edge_limiter = True 261 else: 262 self.use_edge_limiter = False 263 264 265 227 266 ## 228 267 # @brief … … 488 527 489 528 490 491 492 ##493 # @brief Get the total flow through an arbitrary poly line.494 # @param polyline Representation of desired cross section.495 # @param verbose True if this method is to be verbose.496 # @note 'polyline' may contain multiple sections allowing complex shapes.497 # @note Assume absolute UTM coordinates.498 def old_get_flow_through_cross_section(self, polyline, verbose=False):499 """Get the total flow through an arbitrary poly line.500 501 This is a run-time equivalent of the function with same name502 in data_manager.py503 504 Input:505 polyline: Representation of desired cross section - it may contain506 multiple sections allowing for complex shapes. Assume507 absolute UTM coordinates.508 Format [[x0, y0], [x1, y1], ...]509 510 Output:511 Q: Total flow [m^3/s] across given segments.512 """513 514 # Find all intersections and associated triangles.515 segments = self.get_intersecting_segments(polyline, use_cache=True,516 verbose=verbose)517 518 # Get midpoints519 midpoints = segment_midpoints(segments)520 521 # Make midpoints Geospatial instances522 midpoints = ensure_geospatial(midpoints, self.geo_reference)523 524 # Compute flow525 if verbose:526 log.critical('Computing flow through specified cross section')527 528 # Get interpolated values529 xmomentum = self.get_quantity('xmomentum')530 ymomentum = self.get_quantity('ymomentum')531 532 uh = xmomentum.get_values(interpolation_points=midpoints,533 use_cache=True)534 vh = ymomentum.get_values(interpolation_points=midpoints,535 use_cache=True)536 537 # Compute and sum flows across each segment538 total_flow = 0539 for i in range(len(uh)):540 # Inner product of momentum vector with segment normal [m^2/s]541 normal = segments[i].normal542 normal_momentum = uh[i]*normal[0] + vh[i]*normal[1]543 544 # Flow across this segment [m^3/s]545 segment_flow = normal_momentum*segments[i].length546 547 # Accumulate548 total_flow += segment_flow549 550 return total_flow551 529 552 530 ## … … 589 567 590 568 return cross_section.get_energy_through_cross_section(kind) 591 592 593 ##594 # @brief595 # @param polyline Representation of desired cross section.596 # @param kind Select energy type to compute ('specific' or 'total').597 # @param verbose True if this method is to be verbose.598 # @note 'polyline' may contain multiple sections allowing complex shapes.599 # @note Assume absolute UTM coordinates.600 def old_get_energy_through_cross_section(self, polyline,601 kind='total',602 verbose=False):603 """Obtain average energy head [m] across specified cross section.604 605 Inputs:606 polyline: Representation of desired cross section - it may contain607 multiple sections allowing for complex shapes. Assume608 absolute UTM coordinates.609 Format [[x0, y0], [x1, y1], ...]610 kind: Select which energy to compute.611 Options are 'specific' and 'total' (default)612 613 Output:614 E: Average energy [m] across given segments for all stored times.615 616 The average velocity is computed for each triangle intersected by617 the polyline and averaged weighted by segment lengths.618 619 The typical usage of this function would be to get average energy of620 flow in a channel, and the polyline would then be a cross section621 perpendicular to the flow.622 623 #FIXME (Ole) - need name for this energy reflecting that its dimension624 is [m].625 """626 627 from anuga.config import g, epsilon, velocity_protection as h0628 629 # Find all intersections and associated triangles.630 segments = self.get_intersecting_segments(polyline, use_cache=True,631 verbose=verbose)632 633 # Get midpoints634 midpoints = segment_midpoints(segments)635 636 # Make midpoints Geospatial instances637 midpoints = ensure_geospatial(midpoints, self.geo_reference)638 639 # Compute energy640 if verbose: log.critical('Computing %s energy' % kind)641 642 # Get interpolated values643 stage = self.get_quantity('stage')644 elevation = self.get_quantity('elevation')645 xmomentum = self.get_quantity('xmomentum')646 ymomentum = self.get_quantity('ymomentum')647 648 w = stage.get_values(interpolation_points=midpoints, use_cache=True)649 z = elevation.get_values(interpolation_points=midpoints, use_cache=True)650 uh = xmomentum.get_values(interpolation_points=midpoints,651 use_cache=True)652 vh = ymomentum.get_values(interpolation_points=midpoints,653 use_cache=True)654 h = w-z # Depth655 656 # Compute total length of polyline for use with weighted averages657 total_line_length = 0.0658 for segment in segments:659 total_line_length += segment.length660 661 # Compute and sum flows across each segment662 average_energy = 0.0663 for i in range(len(w)):664 # Average velocity across this segment665 if h[i] > epsilon:666 # Use protection against degenerate velocities667 u = uh[i]/(h[i] + h0/h[i])668 v = vh[i]/(h[i] + h0/h[i])669 else:670 u = v = 0.0671 672 speed_squared = u*u + v*v673 kinetic_energy = 0.5*speed_squared/g674 675 if kind == 'specific':676 segment_energy = h[i] + kinetic_energy677 elif kind == 'total':678 segment_energy = w[i] + kinetic_energy679 else:680 msg = 'Energy kind must be either "specific" or "total".'681 msg += ' I got %s' %kind682 683 # Add to weighted average684 weigth = segments[i].length/total_line_length685 average_energy += segment_energy*weigth686 687 return average_energy688 569 689 570 … … 1866 1747 from shallow_water_ext import gravity as gravity_c 1867 1748 1868 xmom = domain.quantities['xmomentum'].explicit_update1869 ymom = domain.quantities['ymomentum'].explicit_update1749 xmom_update = domain.quantities['xmomentum'].explicit_update 1750 ymom_update = domain.quantities['ymomentum'].explicit_update 1870 1751 1871 1752 stage = domain.quantities['stage'] … … 1878 1759 g = domain.g 1879 1760 1880 gravity_c(g, h, z, x, xmom , ymom) #, 1.0e-6)1761 gravity_c(g, h, z, x, xmom_update, ymom_update) #, 1.0e-6) 1881 1762 1882 1763 ## … … 1889 1770 """ 1890 1771 1891 from shallow_water_ext import manning_friction as manning_friction_c 1772 from shallow_water_ext import manning_friction_old 1773 from shallow_water_ext import manning_friction_new 1892 1774 1893 1775 xmom = domain.quantities['xmomentum'] 1894 1776 ymom = domain.quantities['ymomentum'] 1895 1777 1778 x = domain.get_vertex_coordinates() 1779 1896 1780 w = domain.quantities['stage'].centroid_values 1897 z = domain.quantities['elevation']. centroid_values1781 z = domain.quantities['elevation'].vertex_values 1898 1782 1899 1783 uh = xmom.centroid_values … … 1908 1792 g = domain.g 1909 1793 1910 manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update) 1794 if domain.use_new_mannings: 1795 manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, ymom_update) 1796 else: 1797 manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, ymom_update) 1798 1911 1799 1912 1800 ## … … 1919 1807 """ 1920 1808 1921 from shallow_water_ext import manning_friction as manning_friction_c 1809 from shallow_water_ext import manning_friction_old 1810 from shallow_water_ext import manning_friction_new 1922 1811 1923 1812 xmom = domain.quantities['xmomentum'] 1924 1813 ymom = domain.quantities['ymomentum'] 1925 1814 1815 x = domain.get_vertex_coordinates() 1816 1926 1817 w = domain.quantities['stage'].centroid_values 1927 z = domain.quantities['elevation']. centroid_values1818 z = domain.quantities['elevation'].vertex_values 1928 1819 1929 1820 uh = xmom.centroid_values … … 1938 1829 g = domain.g 1939 1830 1940 manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update) 1831 1832 if domain.use_new_mannings: 1833 manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, ymom_update) 1834 else: 1835 manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, ymom_update) 1836 1941 1837 1942 1838 -
anuga_core/source/anuga/shallow_water/shallow_water_ext.c
r7276 r7519 516 516 517 517 518 void _manning_friction (double g, double eps, int N,519 double* w, double* z ,518 void _manning_friction_old(double g, double eps, int N, 519 double* w, double* zv, 520 520 double* uh, double* vh, 521 521 double* eta, double* xmom, double* ymom) { 522 522 523 int k ;524 double S, h ;523 int k, k3; 524 double S, h, z, z0, z1, z2; 525 525 526 526 for (k=0; k<N; k++) { 527 527 if (eta[k] > eps) { 528 h = w[k]-z[k]; 528 k3 = 3*k; 529 // Get bathymetry 530 z0 = zv[k3 + 0]; 531 z1 = zv[k3 + 1]; 532 z2 = zv[k3 + 2]; 533 z = (z0+z1+z2)/3.0; 534 h = w[k]-z; 529 535 if (h >= eps) { 530 536 S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k])); … … 542 548 } 543 549 550 551 void _manning_friction_new(double g, double eps, int N, 552 double* x, double* w, double* zv, 553 double* uh, double* vh, 554 double* eta, double* xmom_update, double* ymom_update) { 555 556 int k, k3, k6; 557 double S, h, z, z0, z1, z2, zs, zx, zy; 558 double x0,y0,x1,y1,x2,y2; 559 560 for (k=0; k<N; k++) { 561 if (eta[k] > eps) { 562 k3 = 3*k; 563 // Get bathymetry 564 z0 = zv[k3 + 0]; 565 z1 = zv[k3 + 1]; 566 z2 = zv[k3 + 2]; 567 568 // Compute bed slope 569 k6 = 6*k; // base index 570 571 x0 = x[k6 + 0]; 572 y0 = x[k6 + 1]; 573 x1 = x[k6 + 2]; 574 y1 = x[k6 + 3]; 575 x2 = x[k6 + 4]; 576 y2 = x[k6 + 5]; 577 578 _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy); 579 580 zs = sqrt(1.0 + zx*zx + zy*zy); 581 z = (z0+z1+z2)/3.0; 582 h = w[k]-z; 583 if (h >= eps) { 584 S = -g * eta[k]*eta[k] * zs * sqrt((uh[k]*uh[k] + vh[k]*vh[k])); 585 S /= pow(h, 7.0/3); //Expensive (on Ole's home computer) 586 //S /= exp((7.0/3.0)*log(h)); //seems to save about 15% over manning_friction 587 //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion 588 589 590 //Update momentum 591 xmom_update[k] += S*uh[k]; 592 ymom_update[k] += S*vh[k]; 593 } 594 } 595 } 596 } 544 597 545 598 /* … … 1043 1096 1044 1097 1045 PyObject *manning_friction (PyObject *self, PyObject *args) {1098 PyObject *manning_friction_new(PyObject *self, PyObject *args) { 1046 1099 // 1047 // manning_friction(g, eps, h, uh, vh, eta, xmom_update, ymom_update) 1100 // manning_friction_new(g, eps, x, h, uh, vh, z, eta, xmom_update, ymom_update) 1101 // 1102 1103 1104 PyArrayObject *x, *w, *z, *uh, *vh, *eta, *xmom, *ymom; 1105 int N; 1106 double g, eps; 1107 1108 if (!PyArg_ParseTuple(args, "ddOOOOOOOO", 1109 &g, &eps, &x, &w, &uh, &vh, &z, &eta, &xmom, &ymom)) { 1110 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction_new could not parse input arguments"); 1111 return NULL; 1112 } 1113 1114 // check that numpy array objects arrays are C contiguous memory 1115 CHECK_C_CONTIG(x); 1116 CHECK_C_CONTIG(w); 1117 CHECK_C_CONTIG(z); 1118 CHECK_C_CONTIG(uh); 1119 CHECK_C_CONTIG(vh); 1120 CHECK_C_CONTIG(eta); 1121 CHECK_C_CONTIG(xmom); 1122 CHECK_C_CONTIG(ymom); 1123 1124 N = w -> dimensions[0]; 1125 1126 _manning_friction_new(g, eps, N, 1127 (double*) x -> data, 1128 (double*) w -> data, 1129 (double*) z -> data, 1130 (double*) uh -> data, 1131 (double*) vh -> data, 1132 (double*) eta -> data, 1133 (double*) xmom -> data, 1134 (double*) ymom -> data); 1135 1136 return Py_BuildValue(""); 1137 } 1138 1139 1140 PyObject *manning_friction_old(PyObject *self, PyObject *args) { 1141 // 1142 // manning_friction_old(g, eps, h, uh, vh, z, eta, xmom_update, ymom_update) 1048 1143 // 1049 1144 … … 1054 1149 1055 1150 if (!PyArg_ParseTuple(args, "ddOOOOOOO", 1056 &g, &eps, &w, &z, &uh, &vh, &eta, 1057 &xmom, &ymom)) { 1058 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction could not parse input arguments"); 1151 &g, &eps, &w, &uh, &vh, &z, &eta, &xmom, &ymom)) { 1152 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction_old could not parse input arguments"); 1059 1153 return NULL; 1060 1154 } … … 1070 1164 1071 1165 N = w -> dimensions[0]; 1072 _manning_friction(g, eps, N, 1073 (double*) w -> data, 1074 (double*) z -> data, 1075 (double*) uh -> data, 1076 (double*) vh -> data, 1077 (double*) eta -> data, 1078 (double*) xmom -> data, 1079 (double*) ymom -> data); 1166 1167 _manning_friction_old(g, eps, N, 1168 (double*) w -> data, 1169 (double*) z -> data, 1170 (double*) uh -> data, 1171 (double*) vh -> data, 1172 (double*) eta -> data, 1173 (double*) xmom -> data, 1174 (double*) ymom -> data); 1080 1175 1081 1176 return Py_BuildValue(""); … … 2721 2816 {"compute_fluxes_ext_kinetic", compute_fluxes_ext_kinetic, METH_VARARGS, "Print out"}, 2722 2817 {"gravity", gravity, METH_VARARGS, "Print out"}, 2723 {"manning_friction", manning_friction, METH_VARARGS, "Print out"}, 2818 {"manning_friction_old", manning_friction_old, METH_VARARGS, "Print out"}, 2819 {"manning_friction_new", manning_friction_new, METH_VARARGS, "Print out"}, 2724 2820 {"flux_function_central", flux_function_central, METH_VARARGS, "Print out"}, 2725 2821 {"balance_deep_and_shallow", balance_deep_and_shallow, -
anuga_core/source/anuga/shallow_water/test_data_manager.py
r7509 r7519 352 352 353 353 extrema = fid.variables['xmomentum.extrema'][:] 354 assert num.allclose(extrema,[-0.06062178, 0.47873023]) or\ 355 num.allclose(extrema, [-0.06062178, 0.47847986]) or\ 356 num.allclose(extrema, [-0.06062178, 0.47848481]) # 27/5/9 354 assert num.allclose(extrema,[-0.06062178, 0.47873023]) or \ 355 num.allclose(extrema, [-0.06062178, 0.47847986]) or \ 356 num.allclose(extrema, [-0.06062178, 0.47848481]) or \ 357 num.allclose(extrema, [-0.06062178, 0.47763887]) # 18/09/09 357 358 358 359 extrema = fid.variables['ymomentum.extrema'][:] … … 2096 2097 fid.close() 2097 2098 2099 #Cleanup 2100 os.remove(prjfile) 2101 os.remove(ascfile) 2102 os.remove(swwfile) 2103 2104 2105 2106 def test_sww2dem_larger_zero(self): 2107 """Test that sww information can be converted correctly to asc/prj 2108 format readable by e.g. ArcView. Here: 2109 2110 ncols 11 2111 nrows 11 2112 xllcorner 308500 2113 yllcorner 6189000 2114 cellsize 10.000000 2115 NODATA_value -9999 2116 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200 2117 -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 2118 -80 -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 2119 -70 -80 -90 -100 -110 -120 -130 -140 -150 -160 -170 2120 -60 -70 -80 -90 -100 -110 -120 -130 -140 -150 -160 2121 -50 -60 -70 -80 -90 -100 -110 -120 -130 -140 -150 2122 -40 -50 -60 -70 -80 -90 -100 -110 -120 -130 -140 2123 -30 -40 -50 -60 -70 -80 -90 -100 -110 -120 -130 2124 -20 -30 -40 -50 -60 -70 -80 -90 -100 -110 -120 2125 -10 -20 -30 -40 -50 -60 -70 -80 -90 -100 -110 2126 0 -10 -20 -30 -40 -50 -60 -70 -80 -90 -100 2127 2128 """ 2129 2130 import time, os 2131 from Scientific.IO.NetCDF import NetCDFFile 2132 2133 #Setup 2134 2135 from mesh_factory import rectangular 2136 2137 #Create basic mesh (100m x 100m) 2138 points, vertices, boundary = rectangular(2, 2, 100, 100) 2139 2140 #Create shallow water domain 2141 domain = Domain(points, vertices, boundary) 2142 domain.default_order = 1 2143 2144 domain.set_name('datatest') 2145 2146 prjfile = domain.get_name() + '_elevation.prj' 2147 ascfile = domain.get_name() + '_elevation.asc' 2148 swwfile = domain.get_name() + '.sww' 2149 2150 domain.set_datadir('.') 2151 domain.format = 'sww' 2152 domain.smooth = True 2153 domain.geo_reference = Geo_reference(56, 308500, 6189000) 2154 2155 # 2156 domain.set_quantity('elevation', 0) 2157 domain.set_quantity('stage', 0) 2158 2159 B = Transmissive_boundary(domain) 2160 domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B}) 2161 2162 2163 # 2164 sww = SWW_file(domain) 2165 sww.store_connectivity() 2166 sww.store_timestep() 2167 2168 domain.tight_slope_limiters = 1 2169 domain.evolve_to_end(finaltime = 0.01) 2170 sww.store_timestep() 2171 2172 # Set theshold for printoptions to somrthing small 2173 # to pickup Rudy's error caused a long sequence of zeros printing 2174 # as [ 0.0, 0.0, 0.0, ... 0.0, 0.0, 0.0] by num.array2string 2175 printoptions = num.get_printoptions() 2176 num.set_printoptions(threshold=9) 2177 2178 cellsize = 10 #10m grid 2179 2180 2181 #Check contents 2182 #Get NetCDF 2183 2184 fid = NetCDFFile(sww.filename, netcdf_mode_r) 2185 2186 # Get the variables 2187 x = fid.variables['x'][:] 2188 y = fid.variables['y'][:] 2189 z = fid.variables['elevation'][:] 2190 time = fid.variables['time'][:] 2191 stage = fid.variables['stage'][:] 2192 2193 2194 #Export to ascii/prj files 2195 sww2dem(domain.get_name(), 2196 quantity = 'elevation', 2197 cellsize = cellsize, 2198 number_of_decimal_places = 9, 2199 verbose = self.verbose, 2200 format = 'asc', 2201 block_size=2) 2202 2203 2204 #Check prj (meta data) 2205 prjid = open(prjfile) 2206 lines = prjid.readlines() 2207 prjid.close() 2208 2209 L = lines[0].strip().split() 2210 assert L[0].strip().lower() == 'projection' 2211 assert L[1].strip().lower() == 'utm' 2212 2213 L = lines[1].strip().split() 2214 assert L[0].strip().lower() == 'zone' 2215 assert L[1].strip().lower() == '56' 2216 2217 L = lines[2].strip().split() 2218 assert L[0].strip().lower() == 'datum' 2219 assert L[1].strip().lower() == 'wgs84' 2220 2221 L = lines[3].strip().split() 2222 assert L[0].strip().lower() == 'zunits' 2223 assert L[1].strip().lower() == 'no' 2224 2225 L = lines[4].strip().split() 2226 assert L[0].strip().lower() == 'units' 2227 assert L[1].strip().lower() == 'meters' 2228 2229 L = lines[5].strip().split() 2230 assert L[0].strip().lower() == 'spheroid' 2231 assert L[1].strip().lower() == 'wgs84' 2232 2233 L = lines[6].strip().split() 2234 assert L[0].strip().lower() == 'xshift' 2235 assert L[1].strip().lower() == '500000' 2236 2237 L = lines[7].strip().split() 2238 assert L[0].strip().lower() == 'yshift' 2239 assert L[1].strip().lower() == '10000000' 2240 2241 L = lines[8].strip().split() 2242 assert L[0].strip().lower() == 'parameters' 2243 2244 2245 #Check asc file 2246 ascid = open(ascfile) 2247 lines = ascid.readlines() 2248 ascid.close() 2249 2250 L = lines[0].strip().split() 2251 assert L[0].strip().lower() == 'ncols' 2252 assert L[1].strip().lower() == '11' 2253 2254 L = lines[1].strip().split() 2255 assert L[0].strip().lower() == 'nrows' 2256 assert L[1].strip().lower() == '11' 2257 2258 L = lines[2].strip().split() 2259 assert L[0].strip().lower() == 'xllcorner' 2260 assert num.allclose(float(L[1].strip().lower()), 308500) 2261 2262 L = lines[3].strip().split() 2263 assert L[0].strip().lower() == 'yllcorner' 2264 assert num.allclose(float(L[1].strip().lower()), 6189000) 2265 2266 L = lines[4].strip().split() 2267 assert L[0].strip().lower() == 'cellsize' 2268 assert num.allclose(float(L[1].strip().lower()), cellsize) 2269 2270 L = lines[5].strip().split() 2271 assert L[0].strip() == 'NODATA_value' 2272 assert L[1].strip().lower() == '-9999' 2273 2274 #Check grid values (FIXME: Use same strategy for other sww2dem tests) 2275 for i, line in enumerate(lines[6:]): 2276 for j, value in enumerate( line.split() ): 2277 assert num.allclose(float(value), 0.0, 2278 atol=1.0e-12, rtol=1.0e-12) 2279 2280 # Note: Equality can be obtained in this case, 2281 # but it is better to use allclose. 2282 #assert float(value) == -(10-i+j)*cellsize 2283 2284 2285 fid.close() 2286 2287 2288 num.set_printoptions(threshold=printoptions['threshold']) 2098 2289 #Cleanup 2099 2290 os.remove(prjfile) … … 11712 11903 if __name__ == "__main__": 11713 11904 #suite = unittest.makeSuite(Test_Data_Manager, 'test_sww2domain2') 11714 suite = unittest.makeSuite(Test_Data_Manager, 'test ')11905 suite = unittest.makeSuite(Test_Data_Manager, 'test_sww') 11715 11906 11716 11907 -
anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py
r7492 r7519 2067 2067 4*S) 2068 2068 2069 2070 2071 def test_manning_friction_old(self): 2072 from anuga.config import g 2073 2074 a = [0.0, 0.0] 2075 b = [0.0, 2.0] 2076 c = [2.0, 0.0] 2077 d = [0.0, 4.0] 2078 e = [2.0, 2.0] 2079 f = [4.0, 0.0] 2080 2081 points = [a, b, c, d, e, f] 2082 # bac, bce, ecf, dbe 2083 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2084 2085 domain = Domain(points, vertices) 2086 2087 #Set up for a gradient of (3,0) at mid triangle (bce) 2088 def slope(x, y): 2089 return 3*x 2090 2091 h = 0.1 2092 def stage(x, y): 2093 return slope(x, y) + h 2094 2095 eta = 0.07 2096 domain.set_quantity('elevation', slope) 2097 domain.set_quantity('stage', stage) 2098 domain.set_quantity('friction', eta) 2099 2100 for name in domain.conserved_quantities: 2101 assert num.allclose(domain.quantities[name].explicit_update, 0) 2102 assert num.allclose(domain.quantities[name].semi_implicit_update, 0) 2103 2104 domain.compute_forcing_terms() 2105 2106 assert num.allclose(domain.quantities['stage'].explicit_update, 0) 2107 assert num.allclose(domain.quantities['xmomentum'].explicit_update, 2108 -g*h*3) 2109 assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0) 2110 2111 assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0) 2112 assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, 2113 0) 2114 assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 2115 0) 2116 2117 #Create some momentum for friction to work with 2118 domain.set_quantity('xmomentum', 1) 2119 S = -g*eta**2 / h**(7.0/3) 2120 2121 domain.compute_forcing_terms() 2122 assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0) 2123 assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, 2124 S) 2125 assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 2126 0) 2127 2128 #A more complex example 2129 domain.quantities['stage'].semi_implicit_update[:] = 0.0 2130 domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0 2131 domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0 2132 2133 domain.set_quantity('xmomentum', 3) 2134 domain.set_quantity('ymomentum', 4) 2135 2136 S = -g*eta**2*5 / h**(7.0/3) 2137 2138 domain.compute_forcing_terms() 2139 2140 assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0) 2141 assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, 2142 3*S) 2143 assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 2144 4*S) 2145 2146 2147 def test_manning_friction_new(self): 2148 from anuga.config import g 2149 2150 a = [0.0, 0.0] 2151 b = [0.0, 2.0] 2152 c = [2.0, 0.0] 2153 d = [0.0, 4.0] 2154 e = [2.0, 2.0] 2155 f = [4.0, 0.0] 2156 2157 points = [a, b, c, d, e, f] 2158 # bac, bce, ecf, dbe 2159 vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]] 2160 2161 domain = Domain(points, vertices) 2162 2163 # Use the new function which takes into account the extra 2164 # wetted area due to slope of bed 2165 domain.set_new_mannings_function(True) 2166 2167 #Set up for a gradient of (3,0) at mid triangle (bce) 2168 def slope(x, y): 2169 return 3*x 2170 2171 h = 0.1 2172 def stage(x, y): 2173 return slope(x, y) + h 2174 2175 eta = 0.07 2176 domain.set_quantity('elevation', slope) 2177 domain.set_quantity('stage', stage) 2178 domain.set_quantity('friction', eta) 2179 2180 for name in domain.conserved_quantities: 2181 assert num.allclose(domain.quantities[name].explicit_update, 0) 2182 assert num.allclose(domain.quantities[name].semi_implicit_update, 0) 2183 2184 domain.compute_forcing_terms() 2185 2186 assert num.allclose(domain.quantities['stage'].explicit_update, 0) 2187 assert num.allclose(domain.quantities['xmomentum'].explicit_update, 2188 -g*h*3) 2189 assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0) 2190 2191 assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0) 2192 assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, 2193 0) 2194 assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 2195 0) 2196 2197 #Create some momentum for friction to work with 2198 domain.set_quantity('xmomentum', 1) 2199 S = -g*eta**2 / h**(7.0/3) * sqrt(10) 2200 2201 domain.compute_forcing_terms() 2202 assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0) 2203 assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, 2204 S) 2205 assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 2206 0) 2207 2208 #A more complex example 2209 domain.quantities['stage'].semi_implicit_update[:] = 0.0 2210 domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0 2211 domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0 2212 2213 domain.set_quantity('xmomentum', 3) 2214 domain.set_quantity('ymomentum', 4) 2215 2216 S = -g*eta**2*5 / h**(7.0/3) * sqrt(10.0) 2217 2218 domain.compute_forcing_terms() 2219 2220 assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0) 2221 assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, 2222 3*S) 2223 assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 2224 4*S) 2225 2069 2226 def test_constant_wind_stress(self): 2070 2227 from anuga.config import rho_a, rho_w, eta_w … … 7334 7491 7335 7492 if __name__ == "__main__": 7336 suite = unittest.makeSuite(Test_Shallow_Water, 'test_volumetric_balance')7337 #suite = unittest.makeSuite(Test_Shallow_Water, 'test')7493 #suite = unittest.makeSuite(Test_Shallow_Water, 'test_volumetric_balance') 7494 suite = unittest.makeSuite(Test_Shallow_Water, 'test') 7338 7495 runner = unittest.TextTestRunner(verbosity=1) 7339 7496 runner.run(suite) -
anuga_core/source/anuga/utilities/log.py
r7276 r7519 11 11 log.console_logging_level = log.INFO 12 12 log.log_logging_level = log.DEBUG 13 log.log_filename ('./my.log')13 log.log_filename = './my.log' 14 14 15 15 # log away! -
anuga_core/source/anuga_parallel/test_parallel.py
r7507 r7519 169 169 assert(num.allclose(submesh['ghost_triangles'][1],true_ghost_triangles[1])) 170 170 assert num.allclose(submesh['ghost_triangles'][2],true_ghost_triangles[2]), ParallelException('X') 171 #if not num.allclose(submesh['ghost_triangles'][2],true_ghost_triangles[2]):172 # import pypar173 # raise Exception, 'Here it is'174 # pypar.abort()175 176 171 177 172 … … 295 290 296 291 #------------------------------------------------------------- 292 297 293 if __name__=="__main__": 298 294 runner = unittest.TextTestRunner()
Note: See TracChangeset
for help on using the changeset viewer.