- Timestamp:
- Jun 17, 2010, 5:34:13 PM (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
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:
Note: See TracChangeset
for help on using the changeset viewer.
