Changeset 7855 for anuga_work/development/2010-projects/anuga_1d
- Timestamp:
- Jun 17, 2010, 5:34:13 PM (14 years ago)
- Location:
- anuga_work/development/2010-projects/anuga_1d
- Files:
-
- 6 edited
- 1 moved
-
base/limiters_python.py (modified) (1 diff)
-
base/quantity.py (modified) (13 diffs)
-
base/test_generic_domain.py (modified) (2 diffs)
-
base/test_quantity_1d.py (moved) (moved from anuga_work/development/2010-projects/anuga_1d/base/test_quantity.py) (4 diffs)
-
channel/channel_domain.py (modified) (4 diffs)
-
channel/profile_channel.py (modified) (1 diff)
-
test_all.py (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/development/2010-projects/anuga_1d/base/limiters_python.py
r7852 r7855 20 20 from numpy import sign, abs, minimum, where 21 21 22 return where( (sign(a)*sign(b) > 0.0) & (sign(a)*sign(c)>0.0), 23 sign(a)*minimum(abs(a),abs(b),abs(c)), 0.0 ) 22 return where( (sign(a)*sign(b) > 0.0) & (sign(a)*sign(c)>0.0), 23 sign(a)*minimum(minimum(abs(a),abs(b)),abs(c)), 0.0 ) 24 24 25 25 26 -
anuga_work/development/2010-projects/anuga_1d/base/quantity.py
r7852 r7855 26 26 #Initialise Quantity using optional vertex values. 27 27 28 from generic_domain import Generic_domain28 from anuga_1d.base.generic_domain import Generic_domain 29 29 30 30 msg = 'First argument in Quantity.__init__ ' 31 msg += 'must be of class Domain (or a subclass thereof)'31 msg += 'must be of class Generic_domain (or a subclass thereof)' 32 32 assert isinstance(domain, Generic_domain), msg 33 33 … … 78 78 self.qmin = numpy.zeros(self.centroid_values.shape, numpy.float) 79 79 80 self.beta = domain.beta 80 self.beta = domain.beta 81 82 83 self.beta_p = numpy.zeros(N,numpy.float) 84 self.beta_m = numpy.zeros(N,numpy.float) 85 self.beta_x = numpy.zeros(N,numpy.float) 86 87 88 self.dx = numpy.zeros((N,2), numpy.float) 89 self.phi = numpy.zeros(N, numpy.float) 81 90 82 91 … … 602 611 603 612 604 605 N = self.centroid_values.shape[0]606 607 613 #Explicit updates 608 614 self.centroid_values += timestep*self.explicit_update 609 615 610 616 #Semi implicit updates 611 denominator = numpy.ones(N, numpy.float)-timestep*self.semi_implicit_update 612 613 if sum(numpy.equal(denominator, 0.0)) > 0.0: 614 msg = 'Zero division in semi implicit update. Call Stephen :-)' 615 raise Exception(msg) 616 else: 617 #Update conserved_quantities from semi implicit updates 618 self.centroid_values /= denominator 617 denominator = 1.0-timestep*self.semi_implicit_update 618 619 # if sum(numpy.equal(denominator, 0.0)) > 0.0: 620 # msg = 'Zero division in semi implicit update. Call Stephen :-)' 621 # raise Exception(msg) 622 # else: 623 # #Update conserved_quantities from semi implicit updates 624 # self.centroid_values /= denominator 625 # 626 627 #Update conserved_quantities from semi implicit updates 628 self.centroid_values /= denominator 619 629 620 630 … … 625 635 626 636 627 628 629 637 N = self.centroid_values.shape[0] 630 638 … … 634 642 X = self.domain.centroids 635 643 636 for k in range(N): 637 638 # first and last elements have boundaries 639 640 if k == 0: 641 642 #Get data 643 k0 = k 644 k1 = k+1 645 k2 = k+2 646 647 q0 = Q[k0] 648 q1 = Q[k1] 649 q2 = Q[k2] 650 651 x0 = X[k0] #V0 centroid 652 x1 = X[k1] #V1 centroid 653 x2 = X[k2] 654 655 #Gradient 656 #G[k] = (q1 - q0)/(x1 - x0) 657 658 G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0) 659 G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1) 660 661 elif k == N-1: 662 663 #Get data 664 k0 = k 665 k1 = k-1 666 k2 = k-2 667 668 q0 = Q[k0] 669 q1 = Q[k1] 670 q2 = Q[k2] 671 672 x0 = X[k0] #V0 centroid 673 x1 = X[k1] #V1 centroid 674 x2 = X[k2] 675 676 #Gradient 677 #G[k] = (q1 - q0)/(x1 - x0) 678 679 G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0) 680 G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1) 681 682 ## q0 = Q[k0] 683 ## q1 = Q[k1] 684 ## 685 ## x0 = X[k0] #V0 centroid 686 ## x1 = X[k1] #V1 centroid 687 ## 688 ## #Gradient 689 ## G[k] = (q1 - q0)/(x1 - x0) 690 691 else: 692 #Interior Volume (2 neighbours) 693 694 #Get data 695 k0 = k-1 696 k2 = k+1 697 698 q0 = Q[k0] 699 q1 = Q[k] 700 q2 = Q[k2] 701 702 x0 = X[k0] #V0 centroid 703 x1 = X[k] #V1 centroid (Self) 704 x2 = X[k2] #V2 centroid 705 706 #Gradient 707 #G[k] = (q2-q0)/(x2-x0) 708 G[k] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0) 644 # first element 645 646 k = 0 647 648 #Get data 649 k0 = k 650 k1 = k+1 651 k2 = k+2 652 653 q0 = Q[k0] 654 q1 = Q[k1] 655 q2 = Q[k2] 656 657 x0 = X[k0] #V0 centroid 658 x1 = X[k1] #V1 centroid 659 x2 = X[k2] 660 661 #Gradient 662 #G[k] = (q1 - q0)/(x1 - x0) 663 664 G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0) 665 G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1) 666 667 #last element 668 k = N-1 669 670 671 k0 = k 672 k1 = k-1 673 k2 = k-2 674 675 q0 = Q[k0] 676 q1 = Q[k1] 677 q2 = Q[k2] 678 679 x0 = X[k0] #V0 centroid 680 x1 = X[k1] #V1 centroid 681 x2 = X[k2] 682 683 #Gradient 684 #G[k] = (q1 - q0)/(x1 - x0) 685 686 G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0) 687 G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1) 688 689 690 #Interior Volume (2 neighbours) 691 692 693 q0 = Q[0:-2] 694 q1 = Q[1:-1] 695 q2 = Q[2:] 696 697 x0 = X[0:-2] #V0 centroid 698 x1 = X[1:-1] #V1 centroid (Self) 699 x2 = X[2:] #V2 centroid 700 701 #Gradient 702 #G[k] = (q2-q0)/(x2-x0) 703 G[1:-1] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0) 704 709 705 710 706 … … 719 715 720 716 def xmin(a,b): 721 return 0.5*(sign(a)+sign(b))*min(abs(a),abs(b)) 717 from numpy import sign, minimum 718 return 0.5*(sign(a)+sign(b))*minimum(abs(a),abs(b)) 722 719 723 720 def xmic(t,a,b): … … 733 730 X = self.domain.centroids 734 731 735 for k in range(N): 736 737 # first and last elements have boundaries 738 739 if k == 0: 740 741 #Get data 742 k0 = k 743 k1 = k+1 744 k2 = k+2 745 746 q0 = Q[k0] 747 q1 = Q[k1] 748 q2 = Q[k2] 749 750 x0 = X[k0] #V0 centroid 751 x1 = X[k1] #V1 centroid 752 x2 = X[k2] 753 754 #Gradient 755 #G[k] = (q1 - q0)/(x1 - x0) 756 757 G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0) 758 G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1) 759 760 elif k == N-1: 761 762 #Get data 763 k0 = k 764 k1 = k-1 765 k2 = k-2 766 767 q0 = Q[k0] 768 q1 = Q[k1] 769 q2 = Q[k2] 770 771 x0 = X[k0] #V0 centroid 772 x1 = X[k1] #V1 centroid 773 x2 = X[k2] 774 775 #Gradient 776 #G[k] = (q1 - q0)/(x1 - x0) 777 778 G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0) 779 G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1) 780 781 ## #Get data 782 ## k0 = k 783 ## k1 = k-1 784 ## 785 ## q0 = Q[k0] 786 ## q1 = Q[k1] 787 ## 788 ## x0 = X[k0] #V0 centroid 789 ## x1 = X[k1] #V1 centroid 790 ## 791 ## #Gradient 792 ## G[k] = (q1 - q0)/(x1 - x0) 793 794 elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2): 795 G[k] = 0.0 796 797 else: 798 #Interior Volume (2 neighbours) 799 800 #Get data 801 k0 = k-1 802 k2 = k+1 803 804 q0 = Q[k0] 805 q1 = Q[k] 806 q2 = Q[k2] 807 808 x0 = X[k0] #V0 centroid 809 x1 = X[k] #V1 centroid (Self) 810 x2 = X[k2] #V2 centroid 811 812 # assuming uniform grid 813 d1 = (q1 - q0)/(x1-x0) 814 d2 = (q2 - q1)/(x2-x1) 815 816 #Gradient 817 #G[k] = (d1+d2)*0.5 818 #G[k] = (d1*(x2-x1) - d2*(x0-x1))/(x2-x0) 819 G[k] = xmic( self.domain.beta, d1, d2 ) 732 #----------------- 733 #first element 734 #----------------- 735 k = 0 736 737 k0 = k 738 k1 = k+1 739 k2 = k+2 740 741 q0 = Q[k0] 742 q1 = Q[k1] 743 q2 = Q[k2] 744 745 x0 = X[k0] #V0 centroid 746 x1 = X[k1] #V1 centroid 747 x2 = X[k2] 748 749 #Gradient 750 #G[k] = (q1 - q0)/(x1 - x0) 751 752 G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0) 753 G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1) 754 755 #------------------- 756 # Last element 757 #------------------- 758 k = N-1 759 760 k0 = k 761 k1 = k-1 762 k2 = k-2 763 764 q0 = Q[k0] 765 q1 = Q[k1] 766 q2 = Q[k2] 767 768 x0 = X[k0] #V0 centroid 769 x1 = X[k1] #V1 centroid 770 x2 = X[k2] 771 772 #Gradient 773 #G[k] = (q1 - q0)/(x1 - x0) 774 775 G[k] = (q1 - q0)*(x2 - x0)*(x2 - x0) - (q2 - q0)*(x1 - x0)*(x1 - x0) 776 G[k] /= (x1 - x0)*(x2 - x0)*(x2 - x1) 777 778 779 780 #------------------------------ 781 #Interior Volume (2 neighbours) 782 #------------------------------ 783 784 q0 = Q[0:-2] 785 q1 = Q[1:-1] 786 q2 = Q[2:] 787 788 x0 = X[0:-2] #V0 centroid 789 x1 = X[1:-1] #V1 centroid (Self) 790 x2 = X[2:] #V2 centroid 791 792 # assuming uniform grid 793 d1 = (q1 - q0)/(x1-x0) 794 d2 = (q2 - q1)/(x2-x1) 795 796 #Gradient 797 G[1:-1] = xmic( self.domain.beta, d1, d2 ) 798 820 799 821 800 … … 855 834 """ Find min and max of this and neighbour's centroid values""" 856 835 836 from numpy import maximum, minimum 857 837 858 838 qmax = self.qmax … … 860 840 861 841 qc = self.centroid_values 862 863 for k in range(N): 864 qmax[k] = qmin[k] = qc[k] 865 for i in range(2): 866 n = self.domain.neighbours[k,i] 867 if n >= 0: 868 qn = qc[n] #Neighbour's centroid value 869 870 qmin[k] = min(qmin[k], qn) 871 qmax[k] = max(qmax[k], qn) 842 843 qmax[:] = qc 844 qmin[:] = qc 845 846 # check left 847 qmax[1:] = maximum(qmax[1:], qc[0:-1]) 848 qmin[1:] = minimum(qmin[1:], qc[0:-1]) 849 850 # check right 851 qmax[0:-1] = maximum(qmax[0:-1], qc[1:]) 852 qmin[0:-1] = minimum(qmin[0:-1], qc[1:]) 853 854 855 856 # for k in range(N): 857 # qmax[k] = qmin[k] = qc[k] 858 # for i in range(2): 859 # n = self.domain.neighbours[k,i] 860 # if n >= 0: 861 # qn = qc[n] #Neighbour's centroid value 862 # 863 # qmin[k] = min(qmin[k], qn) 864 # qmax[k] = max(qmax[k], qn) 872 865 873 866 … … 990 983 import sys 991 984 992 from limiters_python import minmod, minmod_kurganov, m axmod, vanleer, vanalbada985 from limiters_python import minmod, minmod_kurganov, minmod_kurganov_old, maxmod, vanleer, vanalbada 993 986 994 987 limiter = self.domain.limiter … … 1003 996 x1 = self.domain.vertices[:,1] 1004 997 1005 beta_p = numpy.zeros(N,numpy.float) 1006 beta_m = numpy.zeros(N,numpy.float) 1007 beta_x = numpy.zeros(N,numpy.float) 1008 1009 # for k in range(N): 1010 # 1011 # n0 = self.domain.neighbours[k,0] 1012 # n1 = self.domain.neighbours[k,1] 1013 # 1014 # if ( n0 >= 0) & (n1 >= 0): 1015 # #SLOPE DERIVATIVE LIMIT 1016 # beta_p[k] = (qc[k]-qc[k-1])/(C[k]-C[k-1]) 1017 # beta_m[k] = (qc[k+1]-qc[k])/(C[k+1]-C[k]) 1018 # beta_x[k] = (qc[k+1]-qc[k-1])/(C[k+1]-C[k-1]) 1019 1020 1021 1022 n0 = self.domain.neighbours[:,0] 1023 n1 = self.domain.neighbours[:,1] 998 beta_p = self.beta_p 999 beta_m = self.beta_m 1000 beta_x = self.beta_x 1001 phi = self.phi 1002 dx = self.dx 1024 1003 1025 1004 … … 1027 1006 beta_m[:-1] = beta_p[1:] 1028 1007 beta_x[1:-1] = (qc[2:]-qc[:-2])/(xc[2:]-xc[:-2]) 1029 1030 dx = numpy.zeros(qv.shape, numpy.float)1031 1008 1032 1009 dx[:,0] = x0 - xc 1033 1010 dx[:,1] = x1 - xc 1034 1035 1036 phi = numpy.zeros(qc.shape, numpy.float)1037 1011 1038 1012 phi[0] = (qc[1] - qc[0])/(xc[1] - xc[0]) … … 1053 1027 1054 1028 elif limiter == "superbee": 1055 1056 1029 slope1 = minmod(beta_m[1:-1],2.0*beta_p[1:-1]) 1057 1030 slope2 = minmod(2.0*beta_m[1:-1],beta_p[1:-1]) 1058 1031 phi[1:-1] = maxmod(slope1,slope2) 1059 1060 1032 1061 1033 else: -
anuga_work/development/2010-projects/anuga_1d/base/test_generic_domain.py
r7852 r7855 3 3 import unittest 4 4 from math import sqrt, pi 5 5 import numpy 6 6 7 7 from anuga_1d.base.generic_domain import * 8 import numpy 8 9 9 10 10 … … 47 47 48 48 49 49 50 #------------------------------------------------------------- 50 51 if __name__ == "__main__": -
anuga_work/development/2010-projects/anuga_1d/base/test_quantity_1d.py
r7851 r7855 6 6 7 7 8 from generic_domain import Generic_domain as Domain8 from anuga_1d.base.generic_domain import Generic_domain as Domain 9 9 #from shallow_water_domain import flux_function as domain_flux_function 10 10 11 from quantity import *11 from anuga_1d.base.quantity import * 12 12 13 13 … … 37 37 self.domain6 = Domain(self.points6) 38 38 39 40 39 def tearDown(self): 41 40 pass … … 46 45 47 46 assert self.domain4.boundary == {(0, 0): 'left', (3, 1): 'right'} 47 48 48 49 49 def test_creation(self): … … 497 497 498 498 499 def test_find_qmax_qmin(self): 500 quantity = Quantity(self.domain4) 501 502 503 quantity.set_values([1.0, 4.0, 8.0, 2.0], 504 location = 'centroids') 505 506 507 508 quantity.find_qmax_qmin() 509 510 511 assert allclose(quantity.qmax, [4.0, 8.0, 8.0, 8.0]) 512 assert allclose(quantity.qmin, [1.0, 1.0, 2.0, 2.0]) 513 514 499 515 def test_distribute_first_order(self): 500 516 quantity = Quantity(self.domain4) -
anuga_work/development/2010-projects/anuga_1d/channel/channel_domain.py
r7852 r7855 202 202 203 203 if domain.discontinousb: 204 domain.quantities['width'].extrapolate_second_order() 205 206 207 # FIXME (SGR): Replace with numpy.where to speed up code 204 width.extrapolate_second_order() 205 206 208 207 h0 = 1.0e-12 209 208 … … 217 216 w_C[:] = h_C + z_C 218 217 219 # for i in range(N):220 #221 # if a_C[i] <= h0:222 # a_C[i] = 0.0223 # h_C[i] = 0.0224 # d_C[i] = 0.0225 # u_C[i] = 0.0226 # w_C[i] = z_C[i]227 # else:228 #229 # if b_C[i]<=h0:230 # a_C[i] = 0.0231 # h_C[i] = 0.0232 # d_C[i] = 0.0233 # u_C[i] = 0.0234 # w_C[i] = z_C[i]235 #236 # else:237 # h_C[i] = a_C[i]/(b_C[i]+h0/b_C[i])238 # w_C[i] = h_C[i]+z_C[i]239 # u_C[i] = d_C[i]/(a_C[i]+h0/a_C[i])240 218 241 219 … … 260 238 b_V = width.vertex_values 261 239 262 263 #FIXME (SGR): replace with numpy.where264 265 240 266 241 h_V[:] = w_V-z_V … … 271 246 d_V[:] = u_V*a_V 272 247 273 # for i in range(len(h_C)):274 # for j in range(2):275 # h_V[i,j] = w_V[i,j]-z_V[i,j]276 # if h_V[i,j]<h0:277 # h_V[i,j]=0278 # w_V[i,j]=z_V[i,j]279 # a_V[i,j] = b_V[i,j]*h_V[i,j]280 # d_V[i,j]=u_V[i,j]*a_V[i,j]281 248 282 249 return -
anuga_work/development/2010-projects/anuga_1d/channel/profile_channel.py
r7852 r7855 92 92 #p.strip_dirs().sort_stats(-1).print_stats(20) 93 93 94 p.sort_stats('cumulative').print_stats(10) 94 p.sort_stats('cumulative').print_stats(30) 95 96 print p 95 97 96 98 -
anuga_work/development/2010-projects/anuga_1d/test_all.py
r7840 r7855 23 23 24 24 # Directories that should not be searched for test files. 25 exclude_dirs = ['channel', # Special requirements 26 '.svn', # subversion 27 'props', 'wcprops', 'prop-base', 'text-base', 'tmp'] 25 exclude_dirs = ['channel', '.svn'] 26 28 27 29 28 … … 178 177 # @brief Check that the environment is sane. 179 178 # @note Stops here if there is an error. 180 def check_anuga_ import():179 def check_anuga_1d_import(): 181 180 try: 182 181 # importing something that loads quickly 183 import anuga .anuga_exceptions182 import anuga_1d.config 184 183 except ImportError: 185 print "Python cannot import ANUGA module."184 print "Python cannot import ANUGA_1D module." 186 185 print "Check you have followed all steps of its installation." 187 186 import sys … … 190 189 191 190 if __name__ == '__main__': 192 check_anuga_import() 191 check_anuga_1d_import() 192 193 print 'test-all' 193 194 194 195 if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
Note: See TracChangeset
for help on using the changeset viewer.
