Changeset 7818 for anuga_work/development/anuga_1d
- Timestamp:
- Jun 11, 2010, 10:27:11 AM (14 years ago)
- Location:
- anuga_work/development/anuga_1d
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/development/anuga_1d/quantity.py
r7815 r7818 232 232 except ValueError: 233 233 msg = ("Illegal type for variable 'numeric': %s" % type(numeric)) 234 raise Exception , msg234 raise Exception(msg) 235 235 self.set_values_from_constant(numeric) 236 236 … … 299 299 if location not in ['vertices', 'centroids', 'unique_vertices']: 300 300 msg = 'Invalid location: %s, (possible choices vertices, centroids, unique_vertices)' %location 301 raise msg301 raise Exception(msg) 302 302 303 303 if X is None: 304 304 msg = 'Given values are None' 305 raise msg305 raise Exception(msg) 306 306 307 307 import types … … 441 441 if location not in ['vertices', 'centroids', 'unique vertices']: 442 442 msg = 'Invalid location: %s' %location 443 raise msg443 raise Exception(msg) 444 444 445 445 import types, Numeric … … 463 463 if cells is None: 464 464 msg = 'Unique vertex not associated with cells' 465 raise msg465 raise Exception(msg) 466 466 467 467 # Go through all cells, vertex pairs … … 614 614 if sum(equal(denominator, 0.0)) > 0.0: 615 615 msg = 'Zero division in semi implicit update. Call Stephen :-)' 616 raise msg616 raise Exception(msg) 617 617 else: 618 618 #Update conserved_quantities from semi implicit updates -
anuga_work/development/anuga_1d/test_quantity.py
r7815 r7818 233 233 quantity2.set_values(2*quantity1) 234 234 235 print quantity2.vertex_values236 235 assert allclose(quantity2.vertex_values, 237 236 [[0,2], [2,4], [4,6], [6,8]]) … … 470 469 471 470 #Assert that quantities are conserved 472 from numpyimport sum471 from Numeric import sum 473 472 for k in range(quantity.centroid_values.shape[0]): 474 473 assert allclose (quantity.centroid_values[k], … … 562 561 quantity.update(timestep) 563 562 564 sem = array([1.,1.,1.,1.]) /array([1, 2, 3, 4])563 sem = array([1.,1.,1.,1.]) 565 564 denom = ones(4, Float)-timestep*sem 566 565 … … 594 593 assert allclose( quantity.centroid_values, x) 595 594 596 #Test smoothing597 def test_smoothing(self):598 599 600 from shallow_water import Domain, Transmissive_boundary601 from Numeric import zeros, Float602 from anuga.utilities.numerical_tools import mean603 604 605 #Create shallow water domain606 domain = Domain(points10)607 domain.default_order=2608 domain.reduction = mean609 610 611 #Set some field values612 domain.set_quantity('elevation', lambda x: x)613 domain.set_quantity('friction', 0.03)614 615 616 ######################617 # Boundary conditions618 B = Transmissive_boundary(domain)619 domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})620 621 622 ######################623 #Initial condition - with jumps624 625 bed = domain.quantities['elevation'].vertex_values626 stage = zeros(bed.shape, Float)627 628 h = 0.03629 for i in range(stage.shape[0]):630 if i % 2 == 0:631 stage[i,:] = bed[i,:] + h632 else:633 stage[i,:] = bed[i,:]634 635 domain.set_quantity('stage', stage)636 637 stage = domain.quantities['stage']638 639 #Get smoothed stage640 A, V = stage.get_vertex_values(xy=False, smooth=True)641 Q = stage.vertex_values642 643 644 assert A.shape[0] == 9645 assert V.shape[0] == 8646 assert V.shape[1] == 3647 648 #First four points649 assert allclose(A[0], (Q[0,2] + Q[1,1])/2)650 assert allclose(A[1], (Q[1,0] + Q[3,1] + Q[2,2])/3)651 assert allclose(A[2], Q[3,0])652 assert allclose(A[3], (Q[0,0] + Q[5,1] + Q[4,2])/3)653 654 #Center point655 assert allclose(A[4], (Q[0,1] + Q[1,2] + Q[2,0] +\656 Q[5,0] + Q[6,2] + Q[7,1])/6)657 658 659 #Check V660 assert allclose(V[0,:], [3,4,0])661 assert allclose(V[1,:], [1,0,4])662 assert allclose(V[2,:], [4,5,1])663 assert allclose(V[3,:], [2,1,5])664 assert allclose(V[4,:], [6,7,3])665 assert allclose(V[5,:], [4,3,7])666 assert allclose(V[6,:], [7,8,4])667 assert allclose(V[7,:], [5,4,8])668 669 #Get smoothed stage with XY670 X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=True)671 672 assert allclose(A, A1)673 assert allclose(V, V1)674 675 #Check XY676 assert allclose(X[4], 0.5)677 assert allclose(Y[4], 0.5)678 679 assert allclose(X[7], 1.0)680 assert allclose(Y[7], 0.5)681 682 683 684 685 def test_vertex_values_no_smoothing(self):686 687 from mesh_factory import rectangular688 from shallow_water import Domain, Transmissive_boundary689 from Numeric import zeros, Float690 from anuga.utilities.numerical_tools import mean691 692 693 #Create basic mesh694 points, vertices, boundary = rectangular(2, 2)695 696 #Create shallow water domain697 domain = Domain(points, vertices, boundary)698 domain.default_order=2699 domain.reduction = mean700 701 702 #Set some field values703 domain.set_quantity('elevation', lambda x,y: x)704 domain.set_quantity('friction', 0.03)705 706 707 ######################708 #Initial condition - with jumps709 710 bed = domain.quantities['elevation'].vertex_values711 stage = zeros(bed.shape, Float)712 713 h = 0.03714 for i in range(stage.shape[0]):715 if i % 2 == 0:716 stage[i,:] = bed[i,:] + h717 else:718 stage[i,:] = bed[i,:]719 720 domain.set_quantity('stage', stage)721 722 #Get stage723 stage = domain.quantities['stage']724 A, V = stage.get_vertex_values(xy=False, smooth=False)725 Q = stage.vertex_values.flat726 727 for k in range(8):728 assert allclose(A[k], Q[k])729 730 731 for k in range(8):732 assert V[k, 0] == 3*k733 assert V[k, 1] == 3*k+1734 assert V[k, 2] == 3*k+2735 736 737 738 X, Y, A1, V1 = stage.get_vertex_values(xy=True, smooth=False)739 740 741 assert allclose(A, A1)742 assert allclose(V, V1)743 744 #Check XY745 assert allclose(X[1], 0.5)746 assert allclose(Y[1], 0.5)747 assert allclose(X[4], 0.0)748 assert allclose(Y[4], 0.0)749 assert allclose(X[12], 1.0)750 assert allclose(Y[12], 0.0)751 752 753 595 754 596 755 597 #------------------------------------------------------------- 756 598 if __name__ == "__main__": 757 suite = unittest.makeSuite(Test_Quantity, 'test _set')599 suite = unittest.makeSuite(Test_Quantity, 'test') 758 600 runner = unittest.TextTestRunner() 759 601 runner.run(suite)
Note: See TracChangeset
for help on using the changeset viewer.