Changeset 7855


Ignore:
Timestamp:
Jun 17, 2010, 5:34:13 PM (10 years ago)
Author:
steve
Message:

Changing for loop to numpy.where

Location:
anuga_work/development/2010-projects/anuga_1d
Files:
6 edited
1 moved

Legend:

Unmodified
Added
Removed
  • anuga_work/development/2010-projects/anuga_1d/base/limiters_python.py

    r7852 r7855  
    2020    from numpy import sign, abs, minimum, where
    2121
    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
    2425
    2526
  • anuga_work/development/2010-projects/anuga_1d/base/quantity.py

    r7852 r7855  
    2626        #Initialise Quantity using optional vertex values.
    2727       
    28         from generic_domain import Generic_domain
     28        from anuga_1d.base.generic_domain import Generic_domain
    2929
    3030        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)'
    3232        assert isinstance(domain, Generic_domain), msg
    3333
     
    7878        self.qmin = numpy.zeros(self.centroid_values.shape, numpy.float)
    7979
    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)
    8190
    8291
     
    602611
    603612
    604 
    605         N = self.centroid_values.shape[0]
    606 
    607613        #Explicit updates
    608614        self.centroid_values += timestep*self.explicit_update
    609615       
    610616        #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
    619629
    620630
     
    625635
    626636
    627 
    628 
    629637        N = self.centroid_values.shape[0]
    630638
     
    634642        X = self.domain.centroids
    635643
    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
    709705
    710706
     
    719715       
    720716        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))
    722719
    723720        def xmic(t,a,b):
     
    733730        X = self.domain.centroids
    734731
    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
    820799       
    821800
     
    855834        """ Find min and max of this and neighbour's centroid values"""
    856835
     836        from numpy import maximum, minimum
    857837
    858838        qmax = self.qmax
     
    860840
    861841        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)
    872865
    873866
     
    990983        import sys
    991984
    992         from limiters_python import minmod, minmod_kurganov, maxmod, vanleer, vanalbada
     985        from limiters_python import minmod, minmod_kurganov, minmod_kurganov_old, maxmod, vanleer, vanalbada
    993986
    994987        limiter = self.domain.limiter
     
    1003996        x1 = self.domain.vertices[:,1]
    1004997
    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
    10241003
    10251004
     
    10271006        beta_m[:-1] = beta_p[1:]
    10281007        beta_x[1:-1] = (qc[2:]-qc[:-2])/(xc[2:]-xc[:-2])
    1029                
    1030         dx = numpy.zeros(qv.shape, numpy.float)
    10311008
    10321009        dx[:,0] = x0 - xc
    10331010        dx[:,1] = x1 - xc
    1034 
    1035 
    1036         phi = numpy.zeros(qc.shape, numpy.float)
    10371011
    10381012        phi[0] = (qc[1] - qc[0])/(xc[1] - xc[0])
     
    10531027
    10541028        elif limiter == "superbee":
    1055 
    10561029            slope1 = minmod(beta_m[1:-1],2.0*beta_p[1:-1])
    10571030            slope2 = minmod(2.0*beta_m[1:-1],beta_p[1:-1])
    10581031            phi[1:-1] = maxmod(slope1,slope2)
    1059 
    10601032
    10611033        else:
  • anuga_work/development/2010-projects/anuga_1d/base/test_generic_domain.py

    r7852 r7855  
    33import unittest
    44from math import sqrt, pi
    5 
     5import numpy
    66
    77from anuga_1d.base.generic_domain import *
    8 import numpy
     8
    99
    1010
     
    4747
    4848
     49
    4950#-------------------------------------------------------------
    5051if __name__ == "__main__":
  • anuga_work/development/2010-projects/anuga_1d/base/test_quantity_1d.py

    r7851 r7855  
    66
    77
    8 from generic_domain import Generic_domain as Domain
     8from anuga_1d.base.generic_domain import Generic_domain as Domain
    99#from shallow_water_domain import flux_function as domain_flux_function
    1010
    11 from quantity import *
     11from anuga_1d.base.quantity import *
    1212
    1313
     
    3737        self.domain6 = Domain(self.points6)
    3838
    39 
    4039    def tearDown(self):
    4140        pass
     
    4645
    4746        assert self.domain4.boundary  == {(0, 0): 'left', (3, 1): 'right'}
     47
    4848       
    4949    def test_creation(self):
     
    497497
    498498
     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
    499515    def test_distribute_first_order(self):
    500516        quantity = Quantity(self.domain4)
  • anuga_work/development/2010-projects/anuga_1d/channel/channel_domain.py

    r7852 r7855  
    202202
    203203    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
    208207    h0 = 1.0e-12
    209208
     
    217216    w_C[:] = h_C + z_C
    218217   
    219 #    for i in range(N):
    220 #
    221 #        if a_C[i] <= h0:
    222 #           a_C[i] = 0.0
    223 #           h_C[i] = 0.0
    224 #            d_C[i] = 0.0
    225 #            u_C[i] = 0.0
    226 #            w_C[i] = z_C[i]
    227 #        else:
    228 #
    229 #            if b_C[i]<=h0:
    230 #                a_C[i] = 0.0
    231 #                h_C[i] = 0.0
    232 #                d_C[i] = 0.0
    233 #                u_C[i] = 0.0
    234 #                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])
    240218               
    241219
     
    260238    b_V  = width.vertex_values
    261239
    262    
    263     #FIXME (SGR): replace with numpy.where
    264 
    265240
    266241    h_V[:] = w_V-z_V
     
    271246    d_V[:] = u_V*a_V
    272247
    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]=0
    278 #                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]
    281248   
    282249    return
  • anuga_work/development/2010-projects/anuga_1d/channel/profile_channel.py

    r7852 r7855  
    9292#p.strip_dirs().sort_stats(-1).print_stats(20)
    9393
    94 p.sort_stats('cumulative').print_stats(10)
     94p.sort_stats('cumulative').print_stats(30)
     95
     96print p
    9597
    9698
  • anuga_work/development/2010-projects/anuga_1d/test_all.py

    r7840 r7855  
    2323
    2424# 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']
     25exclude_dirs = ['channel',  '.svn']
     26               
    2827
    2928
     
    178177# @brief Check that the environment is sane.
    179178# @note Stops here if there is an error.
    180 def check_anuga_import():
     179def check_anuga_1d_import():
    181180    try:
    182181        # importing something that loads quickly
    183         import anuga.anuga_exceptions
     182        import anuga_1d.config
    184183    except ImportError:
    185         print "Python cannot import ANUGA module."
     184        print "Python cannot import ANUGA_1D module."
    186185        print "Check you have followed all steps of its installation."
    187186        import sys
     
    190189
    191190if __name__ == '__main__':
    192     check_anuga_import()
     191    check_anuga_1d_import()
     192
     193    print 'test-all'
    193194
    194195    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
Note: See TracChangeset for help on using the changeset viewer.