Changeset 7852


Ignore:
Timestamp:
Jun 16, 2010, 5:30:17 PM (14 years ago)
Author:
steve
Message:

Moving calculation of limiters to numpy calculations

Location:
anuga_work/development/2010-projects/anuga_1d
Files:
1 added
12 edited

Legend:

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

    r7842 r7852  
    1010import numpy
    1111
    12 def interval_mesh(n, x_0=0.0, x_1=1.0):
     12def uniform_mesh(n, x_0=0.0, x_1=1.0):
    1313    """Create points, and boundary for a uniform mesh with n sub-interval
    1414    ranging from x_0 to x_1
  • anuga_work/development/2010-projects/anuga_1d/base/quantity.py

    r7850 r7852  
    838838        second order scheme.
    839839        """
     840
    840841        if self.domain.limiter == "pyvolution":
    841             #Z = self.gradients
    842             #print "gradients 1",Z
    843             self.compute_gradients()
    844             #print "gradients 2",Z
    845 
    846             #Z = self.gradients
    847             #print "gradients 1",Z
    848             #self.compute_minmod_gradients()
    849             #print "gradients 2", Z
    850 
    851             G = self.gradients
    852             V = self.domain.vertices
    853             qc = self.centroid_values
    854             qv = self.vertex_values       
    855 
    856             #Check each triangle
    857             for k in range(self.domain.number_of_elements):
    858                 #Centroid coordinates
    859                 x = self.domain.centroids[k]
    860 
    861                 #vertex coordinates
    862                 x0, x1 = V[k,:]
    863 
    864                 #Extrapolate
    865                 qv[k,0] = qc[k] + G[k]*(x0-x)
    866                 qv[k,1] = qc[k] + G[k]*(x1-x)
    867842            self.limit_pyvolution()
    868843
    869844        elif self.domain.limiter == "minmod_steve":
    870845            self.limit_minmod()
     846
    871847        else:
    872848            self.limit_range()
    873849       
    874        
    875 
    876     def limit_minmod(self):
    877         #Z = self.gradients
    878         #print "gradients 1",Z
    879         self.compute_minmod_gradients()
    880         #print "gradients 2", Z
    881 
    882         G = self.gradients
    883         V = self.domain.vertices
    884         qc = self.centroid_values
    885         qv = self.vertex_values       
    886        
    887         #Check each triangle
    888         for k in range(self.domain.number_of_elements):
    889             #Centroid coordinates
    890             x = self.domain.centroids[k]
    891 
    892             #vertex coordinates
    893             x0, x1 = V[k,:]
    894 
    895             #Extrapolate
    896             qv[k,0] = qc[k] + G[k]*(x0-x)
    897             qv[k,1] = qc[k] + G[k]*(x1-x)
    898 
    899  
    900     def limit_pyvolution(self):
    901         """
    902         Limit slopes for each volume to eliminate artificial variance
    903         introduced by e.g. second order extrapolator
    904 
    905         This is an unsophisticated limiter as it does not take into
    906         account dependencies among quantities.
    907 
    908         precondition:
    909         vertex values are estimated from gradient
    910         postcondition:
    911         vertex values are updated
    912         """
    913 
    914 
    915         N = self.domain.number_of_elements
    916         beta = self.domain.beta
    917         #beta = 0.8
    918 
    919         qc = self.centroid_values
    920         qv = self.vertex_values
    921 
    922         #Find min and max of this and neighbour's centroid values
     850
     851
     852
     853
     854    def find_qmax_qmin(self):
     855        """ Find min and max of this and neighbour's centroid values"""
     856
     857
    923858        qmax = self.qmax
    924859        qmin = self.qmin
     860
     861        qc = self.centroid_values
    925862
    926863        for k in range(N):
     
    935872
    936873
     874
     875    def limit_minmod(self):
     876        #Z = self.gradients
     877        #print "gradients 1",Z
     878        self.compute_minmod_gradients()
     879        #print "gradients 2", Z
     880
     881        G = self.gradients
     882        V = self.domain.vertices
     883        qc = self.centroid_values
     884        qv = self.vertex_values       
     885
     886        x = self.domain.centroids
     887
     888        x0 = V[:,0]
     889        x1 = V[:,1]
     890
     891        #Extrapolate
     892        qv[:,0] = qc + G*(x0-x)
     893        qv[:,1] = qc + G*(x1-x)
     894
     895#        #Check each triangle
     896#        for k in range(self.domain.number_of_elements):
     897#            #Centroid coordinates
     898#            x = self.domain.centroids[k]
     899#
     900#            #vertex coordinates
     901#            x0, x1 = V[k,:]
     902#
     903#            #Extrapolate
     904#            qv[k,0] = qc[k] + G[k]*(x0-x)
     905#            qv[k,1] = qc[k] + G[k]*(x1-x)
     906
     907 
     908    def limit_pyvolution(self):
     909        """
     910        Limit slopes for each volume to eliminate artificial variance
     911        introduced by e.g. second order extrapolator
     912
     913        This is an unsophisticated limiter as it does not take into
     914        account dependencies among quantities.
     915
     916        precondition:
     917        vertex values are estimated from gradient
     918        postcondition:
     919        vertex values are updated
     920        """
     921
     922
     923        N = self.domain.number_of_elements
     924        beta = self.domain.beta
     925        #beta = 0.8
     926
     927        self.compute_gradients()
     928
     929
     930        G = self.gradients
     931        V = self.domain.vertices
     932        C = self.domain.centroids
     933        qc = self.centroid_values
     934        qv = self.vertex_values
     935
     936        V0 = V[:,0]
     937        V1 = V[:,1]
     938
     939        # Extrapolate to vertices
     940        qv[:,0] = qc + G*(V0-C)
     941        qv[:,1] = qc + G*(V1-C)
     942
     943
     944        # Find max and min values
     945        self.find_qmax_qmin()
     946
     947        qmax = self.qmax
     948        qmin = self.qmin
     949
    937950        #Diffences between centroids and maxima/minima
    938951        dqmax = qmax - qc
     
    941954        #Deltas between vertex and centroid values
    942955        dq = numpy.zeros(qv.shape, numpy.float)
    943         for i in range(2):
    944             dq[:,i] = qv[:,i] - qc
    945 
    946         #Phi limiter
    947         for k in range(N):
    948 
    949             #Find the gradient limiter (phi) across vertices
    950             phi = 1.0
    951             for i in range(2):
    952                 r = 1.0
    953                 if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
    954                 if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
    955 
    956                 phi = min( min(r*beta, 1), phi )
    957 
    958             #Then update using phi limiter
    959             for i in range(2):
    960                 qv[k,i] = qc[k] + phi*dq[k,i]
     956
     957        dq[:,0] = qv[:,0] - qc
     958        dq[:,1] = qv[:,1] - qc
     959
     960        phi = numpy.ones(qc.shape, numpy.float)
     961
     962        r0 = numpy.where(dq[:,0]>0.0,dqmax/dq[:,0],1.0)
     963        r0 = numpy.where(dq[:,0]<0.0,dqmin/dq[:,0],r0)
     964
     965        r1 = numpy.where(dq[:,1]>0.0,dqmax/dq[:,1],1.0)
     966        r1 = numpy.where(dq[:,1]<0.0,dqmin/dq[:,1],r1)
     967
     968        phi = numpy.min(r0*beta,numpy.min(r1*beta,1.0))
     969
     970        qv[:,0] = qc + phi*dq[:,0]
     971        qv[:,1] = qc + phi*dq[:,1]
     972
     973#        #Phi limiter
     974#        for k in range(N):
     975#
     976#            #Find the gradient limiter (phi) across vertices
     977#            phi = 1.0
     978#            for i in range(2):
     979#                r = 1.0
     980#                if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
     981#                if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
     982#
     983#                phi = min( min(r*beta, 1), phi )
     984#
     985#            #Then update using phi limiter
     986#            for i in range(2):
     987#                qv[k,i] = qc[k] + phi*dq[k,i]
    961988
    962989    def limit_range(self):
    963990        import sys
    964991
    965         from util import minmod, minmod_kurganov, maxmod, vanleer, vanalbada
     992        from limiters_python import minmod, minmod_kurganov, maxmod, vanleer, vanalbada
    966993
    967994        limiter = self.domain.limiter
     
    972999        qc = self.centroid_values
    9731000        qv = self.vertex_values
    974         C = self.domain.centroids
    975         X = self.domain.vertices
     1001        xc = self.domain.centroids
     1002        x0 = self.domain.vertices[:,0]
     1003        x1 = self.domain.vertices[:,1]
     1004
    9761005        beta_p = numpy.zeros(N,numpy.float)
    9771006        beta_m = numpy.zeros(N,numpy.float)
    9781007        beta_x = numpy.zeros(N,numpy.float)
    9791008       
    980         for k in range(N):
    981        
    982             n0 = self.domain.neighbours[k,0]
    983             n1 = self.domain.neighbours[k,1]
     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]
     1024
     1025
     1026        beta_p[1:]  = (qc[1:]-qc[:-1])/(xc[1:]-xc[:-1])
     1027        beta_m[:-1] = beta_p[1:]
     1028        beta_x[1:-1] = (qc[2:]-qc[:-2])/(xc[2:]-xc[:-2])
     1029               
     1030        dx = numpy.zeros(qv.shape, numpy.float)
     1031
     1032        dx[:,0] = x0 - xc
     1033        dx[:,1] = x1 - xc
     1034
     1035
     1036        phi = numpy.zeros(qc.shape, numpy.float)
     1037
     1038        phi[0] = (qc[1] - qc[0])/(xc[1] - xc[0])
     1039        phi[N-1] = (qc[N-1] - qc[N-2])/(xc[N-1] - xc[N-2])
     1040
     1041        if limiter == "minmod":
     1042            phi[1:-1] = minmod(beta_p[1:-1],beta_m[1:-1])
     1043
     1044        elif limiter == "vanleer":
     1045            phi[1:-1] = vanleer(beta_p[1:-1],beta_m[1:-1])
    9841046           
    985             if ( n0 >= 0) & (n1 >= 0):
    986                 #SLOPE DERIVATIVE LIMIT
    987                 beta_p[k] = (qc[k]-qc[k-1])/(C[k]-C[k-1])
    988                 beta_m[k] = (qc[k+1]-qc[k])/(C[k+1]-C[k])
    989                 beta_x[k] = (qc[k+1]-qc[k-1])/(C[k+1]-C[k-1])
    990                
    991         dq = numpy.zeros(qv.shape, numpy.float)
    992         for i in range(2):
    993             dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
    994            
     1047        elif limiter == "vanalbada":
     1048            phi[1:-1] = vanalbada(beta_p[1:-1],beta_m[1:-1])
     1049
     1050        elif limiter == "minmod_kurganov":
     1051            theta = 2.0
     1052            phi[1:-1] = minmod_kurganov(theta*beta_p[1:-1],theta*beta_m[1:-1], beta_x[1:-1])
     1053
     1054        elif limiter == "superbee":
     1055
     1056            slope1 = minmod(beta_m[1:-1],2.0*beta_p[1:-1])
     1057            slope2 = minmod(2.0*beta_m[1:-1],beta_p[1:-1])
     1058            phi[1:-1] = maxmod(slope1,slope2)
     1059
     1060
     1061        else:
     1062            msg = 'Unknown limiter'
     1063            raise Exception, msg
     1064   
     1065        qv[:,0] = qc + phi*dx[:,0]
     1066        qv[:,1] = qc + phi*dx[:,1]
     1067
     1068
    9951069        #Phi limiter
    996         for k in range(N):
    997             n0 = self.domain.neighbours[k,0]
    998             n1 = self.domain.neighbours[k,1]
    999             if n0 < 0:
    1000                 phi = (qc[k+1] - qc[k])/(C[k+1] - C[k])
    1001             elif n1 < 0:
    1002                 phi = (qc[k] - qc[k-1])/(C[k] - C[k-1])
    1003             #elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
    1004             #    phi = 0.0
    1005             else:
    1006                 if limiter == "minmod":
    1007                     phi = minmod(beta_p[k],beta_m[k])
    1008 
    1009                 elif limiter == "minmod_kurganov":#Change this
    1010                     # Also known as monotonized central difference limiter
    1011                     # if theta = 2.0
    1012                     theta = 2.0
    1013                     phi = minmod_kurganov(theta*beta_p[k],theta*beta_m[k],beta_x[k])
    1014                 elif limiter == "superbee":
    1015                     slope1 = minmod(beta_m[k],2.0*beta_p[k])
    1016                     slope2 = minmod(2.0*beta_m[k],beta_p[k])
    1017                     phi = maxmod(slope1,slope2)
    1018 
    1019                 elif limiter == "vanleer":
    1020                     phi = vanleer(beta_p[k],beta_m[k])
    1021 
    1022                 elif limiter == "vanalbada":
    1023                     phi = vanalbada(beta_m[k],beta_p[k])
    1024            
    1025             for i in range(2):
    1026                 qv[k,i] = qc[k] + phi*dq[k,i]
     1070#        for k in range(N):
     1071#            n0 = self.domain.neighbours[k,0]
     1072#            n1 = self.domain.neighbours[k,1]
     1073#            if n0 < 0:
     1074#                phi = (qc[k+1] - qc[k])/(xc[k+1] - xc[k])
     1075#            elif n1 < 0:
     1076#                phi = (qc[k] - qc[k-1])/(xc[k] - xc[k-1])
     1077#            #elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
     1078#            #    phi = 0.0
     1079#            else:
     1080#                if limiter == "minmod":
     1081#                    phi = minmod(beta_p[k],beta_m[k])
     1082#
     1083#                elif limiter == "minmod_kurganov":#Change this
     1084#                    # Also known as monotonized central difference limiter
     1085#                    # if theta = 2.0
     1086#                    theta = 2.0
     1087#                    phi = minmod_kurganov(theta*beta_p[k],theta*beta_m[k],beta_x[k])
     1088#
     1089#                elif limiter == "superbee":
     1090#                    slope1 = minmod(beta_m[k],2.0*beta_p[k])
     1091#                    slope2 = minmod(2.0*beta_m[k],beta_p[k])
     1092#                    phi = maxmod(slope1,slope2)
     1093#
     1094#                elif limiter == "vanleer":
     1095#                    phi = vanleer(beta_p[k],beta_m[k])
     1096#
     1097#                elif limiter == "vanalbada":
     1098#                    phi = vanalbada(beta_m[k],beta_p[k])
     1099#
     1100#            for i in range(2):
     1101#                qv[k,i] = qc[k] + phi*dx[k,i]
    10271102
    10281103    def limit_steve_slope(self):   
     
    11071182
    11081183
    1109 def newLinePlot(title='Simple Plot'):
    1110     import pylab as g
    1111     g.ion()
    1112     g.hold(False)
    1113     g.title(title)
    1114     g.xlabel('x')
    1115     g.ylabel('y')
    1116    
    1117 
    1118 def linePlot(x,y):
    1119     import pylab as g
    1120     g.plot(x.flat,y.flat)
    1121 
    1122 
    1123 def closePlots():
    1124     import pylab as g
    1125     g.close('all')
     1184
    11261185   
    11271186if __name__ == "__main__":
    11281187    #from domain import Domain
    11291188    from generic_domain import Generic_domain as Domain
     1189
     1190
     1191    def newLinePlot(title='Simple Plot'):
     1192        import pylab as g
     1193        g.ion()
     1194        g.hold(False)
     1195        g.title(title)
     1196        g.xlabel('x')
     1197        g.ylabel('y')
     1198
     1199
     1200    def linePlot(x,y):
     1201        import pylab as g
     1202        g.plot(x.flat,y.flat)
     1203
     1204
     1205    def closePlots():
     1206        import pylab as g
     1207        g.close('all')
     1208
     1209       
    11301210   
    11311211    points1 = [0.0, 1.0, 2.0, 3.0]
     
    12131293
    12141294    raw_input('press return to quit')
    1215 closePlots()
     1295
     1296    closePlots()
  • anuga_work/development/2010-projects/anuga_1d/base/test_generic_domain.py

    r7840 r7852  
    55
    66
    7 from anuga_1d.generic.generic_domain import *
     7from anuga_1d.base.generic_domain import *
    88import numpy
    99
  • anuga_work/development/2010-projects/anuga_1d/base/util.py

    r7839 r7852  
    1919    return a
    2020
    21 def  minmod(beta_p,beta_m):
    22     if (abs(beta_p) < abs(beta_m)) & (beta_p*beta_m > 0.0):
    23         phi = beta_p
    24     elif (abs(beta_m) < abs(beta_p)) & (beta_p*beta_m > 0.0):
    25         phi = beta_m
    26     else:
    27         phi = 0.0
    28     return phi
    2921
    30 def  minmod_kurganov(a,b,c):
    31     from numpy import sign
    32     if sign(a)==sign(b)==sign(c):
    33         return sign(a)*min(abs(a),abs(b),abs(c))
    34     else:
    35         return 0.0
    36 
    37 def  maxmod(a,b):
    38     if (abs(a) > abs(b)) & (a*b > 0.0):
    39         phi = a
    40     elif (abs(b) > abs(a)) & (a*b > 0.0):
    41         phi = b
    42     else:
    43         phi = 0.0
    44     return phi
    45 
    46 def vanleer(a,b):
    47     if abs(a)+abs(b) > 1e-12:
    48         return (a*abs(b)+abs(a)*b)/(abs(a)+abs(b))
    49     else:
    50         return 0.0
    51 
    52 def vanalbada(a,b):
    53     if a*a+b*b > 1e-12:
    54         return (a*a*b+a*b*b)/(a*a+b*b)
    55     else:
    56         return 0.0
    57 
  • anuga_work/development/2010-projects/anuga_1d/channel/channel_domain.py

    r7842 r7852  
    5353
    5454
    55 from anuga_1d.generic.generic_domain import *
     55from anuga_1d.base.generic_domain import *
    5656import numpy
    5757
     
    9595       
    9696        #Reduction operation for get_vertex_values
    97         from anuga_1d.generic.util import mean
     97        from anuga_1d.base.util import mean
    9898        self.reduction = mean
    9999        #self.reduction = min  #Looks better near steep slopes
  • anuga_work/development/2010-projects/anuga_1d/channel/channel_flow_1_padarn.py

    r7842 r7852  
    66from anuga_1d.channel.channel_domain import *
    77from anuga_1d.config import g, epsilon
    8 from anuga_1d.generic.generic_mesh import interval_mesh
     8from anuga_1d.base.generic_mesh import interval_mesh
    99
    1010
  • anuga_work/development/2010-projects/anuga_1d/channel/profile_channel.py

    r7842 r7852  
    1515from anuga_1d.channel.channel_domain import *
    1616from anuga_1d.config import g, epsilon
    17 from anuga_1d.generic.generic_mesh import interval_mesh
     17from anuga_1d.base.generic_mesh import uniform_mesh
    1818
    1919
     
    5454
    5555    # Create domain with centroid points as defined above
    56     domain = Domain(*interval_mesh(N))
     56    domain = Domain(*uniform_mesh(N))
    5757
    5858
     
    7272    domain.set_CFL(1.0)
    7373    domain.set_limiter("vanleer")
     74    #domain.set_limiter("minmod")
    7475    #domain.h0=0.0001
    7576
  • anuga_work/development/2010-projects/anuga_1d/compile_all.py

    r7840 r7852  
    99
    1010os.chdir('..')
    11 os.chdir('generic')
     11os.chdir('base')
    1212execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
    1313
  • anuga_work/development/2010-projects/anuga_1d/sww/sww_boundary_conditions.py

    r7840 r7852  
    77__date__ ="$05/06/2010 5:44:05 PM$"
    88
    9 from anuga_1d.generic.generic_domain import *
     9from anuga_1d.base.generic_domain import *
    1010
    1111class Dirichlet_boundary(Boundary):
  • anuga_work/development/2010-projects/anuga_1d/sww/sww_domain.py

    r7840 r7852  
    4545import numpy
    4646
    47 from anuga_1d.generic.generic_domain import Generic_domain
     47from anuga_1d.base.generic_domain import Generic_domain
    4848from sww_boundary_conditions import *
    4949from sww_forcing_terms import *
     
    8484       
    8585        #Reduction operation for get_vertex_values
    86         from anuga_1d.generic.util import mean
     86        from anuga_1d.base.util import mean
    8787        self.reduction = mean
    8888        #self.reduction = min  #Looks better near steep slopes
  • anuga_work/development/2010-projects/anuga_1d/sww/sww_forcing_terms.py

    r7840 r7852  
    1414    """
    1515
    16     from anuga_1d.generic.util import gradient
     16    from anuga_1d.base.util import gradient
    1717
    1818    xmom  = domain.quantities['xmomentum'].explicit_update
  • anuga_work/development/2010-projects/anuga_1d/sww/sww_vel_domain.py

    r7840 r7852  
    4444
    4545
    46 from anuga_1d.generic.generic_domain import *
     46from anuga_1d.base.generic_domain import *
    4747from sww_boundary_conditions import *
    4848from sww_forcing_terms import *
     
    8282       
    8383        #Reduction operation for get_vertex_values
    84         from anuga_1d.generic.util import mean
     84        from anuga_1d.base.util import mean
    8585        self.reduction = mean
    8686        #self.reduction = min  #Looks better near steep slopes
Note: See TracChangeset for help on using the changeset viewer.