Changeset 7852

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:

• base/generic_mesh.py (modified) (1 diff)
• base/limiters_python.py (added)
• base/quantity.py (modified) (6 diffs)
• base/test_generic_domain.py (modified) (1 diff)
• base/util.py (modified) (1 diff)
• channel/channel_domain.py (modified) (2 diffs)
• channel/channel_flow_1_padarn.py (modified) (1 diff)
• channel/profile_channel.py (modified) (3 diffs)
• compile_all.py (modified) (1 diff)
• sww/sww_boundary_conditions.py (modified) (1 diff)
• sww/sww_domain.py (modified) (2 diffs)
• sww/sww_forcing_terms.py (modified) (1 diff)
• sww/sww_vel_domain.py (modified) (2 diffs)

anuga_work/development/2010-projects/anuga_1d/base/generic_mesh.py

@@ -10 +10 @@
 import numpy
 
-def interval_mesh(n, x_0=0.0, x_1=1.0):
+def uniform_mesh(n, x_0=0.0, x_1=1.0):
     """Create points, and boundary for a uniform mesh with n sub-interval
     ranging from x_0 to x_1

anuga_work/development/2010-projects/anuga_1d/base/quantity.py

@@ -838 +838 @@
         second order scheme.
         """
+
         if self.domain.limiter == "pyvolution":
-            #Z = self.gradients
-            #print "gradients 1",Z
-            self.compute_gradients()
-            #print "gradients 2",Z
-
-            #Z = self.gradients
-            #print "gradients 1",Z
-            #self.compute_minmod_gradients()
-            #print "gradients 2", Z
-
-            G = self.gradients
-            V = self.domain.vertices
-            qc = self.centroid_values
-            qv = self.vertex_values        
-
-            #Check each triangle
-            for k in range(self.domain.number_of_elements):
-                #Centroid coordinates
-                x = self.domain.centroids[k]
-
-                #vertex coordinates
-                x0, x1 = V[k,:]
-
-                #Extrapolate
-                qv[k,0] = qc[k] + G[k]*(x0-x)
-                qv[k,1] = qc[k] + G[k]*(x1-x)
             self.limit_pyvolution()
 
         elif self.domain.limiter == "minmod_steve":
             self.limit_minmod()
+
         else:
             self.limit_range()
         
-        
-
-    def limit_minmod(self):
-        #Z = self.gradients
-        #print "gradients 1",Z
-        self.compute_minmod_gradients()
-        #print "gradients 2", Z
-
-        G = self.gradients
-        V = self.domain.vertices
-        qc = self.centroid_values
-        qv = self.vertex_values        
-        
-        #Check each triangle
-        for k in range(self.domain.number_of_elements):
-            #Centroid coordinates
-            x = self.domain.centroids[k]
-
-            #vertex coordinates
-            x0, x1 = V[k,:]
-
-            #Extrapolate
-            qv[k,0] = qc[k] + G[k]*(x0-x)
-            qv[k,1] = qc[k] + G[k]*(x1-x)
-
- 
-    def limit_pyvolution(self):
-        """
-        Limit slopes for each volume to eliminate artificial variance
-        introduced by e.g. second order extrapolator
-
-        This is an unsophisticated limiter as it does not take into
-        account dependencies among quantities.
-
-        precondition:
-        vertex values are estimated from gradient
-        postcondition:
-        vertex values are updated
-        """
-
-
-        N = self.domain.number_of_elements
-        beta = self.domain.beta
-        #beta = 0.8
-
-        qc = self.centroid_values
-        qv = self.vertex_values
-
-        #Find min and max of this and neighbour's centroid values
+
+
+
+
+    def find_qmax_qmin(self):
+        """ Find min and max of this and neighbour's centroid values"""
+
+
         qmax = self.qmax
         qmin = self.qmin
+
+        qc = self.centroid_values
 
         for k in range(N):
@@ -935 +872 @@
 
 
+
+    def limit_minmod(self):
+        #Z = self.gradients
+        #print "gradients 1",Z
+        self.compute_minmod_gradients()
+        #print "gradients 2", Z
+
+        G = self.gradients
+        V = self.domain.vertices
+        qc = self.centroid_values
+        qv = self.vertex_values        
+
+        x = self.domain.centroids
+
+        x0 = V[:,0]
+        x1 = V[:,1]
+
+        #Extrapolate
+        qv[:,0] = qc + G*(x0-x)
+        qv[:,1] = qc + G*(x1-x)
+
+#        #Check each triangle
+#        for k in range(self.domain.number_of_elements):
+#            #Centroid coordinates
+#            x = self.domain.centroids[k]
+#
+#            #vertex coordinates
+#            x0, x1 = V[k,:]
+#
+#            #Extrapolate
+#            qv[k,0] = qc[k] + G[k]*(x0-x)
+#            qv[k,1] = qc[k] + G[k]*(x1-x)
+
+ 
+    def limit_pyvolution(self):
+        """
+        Limit slopes for each volume to eliminate artificial variance
+        introduced by e.g. second order extrapolator
+
+        This is an unsophisticated limiter as it does not take into
+        account dependencies among quantities.
+
+        precondition:
+        vertex values are estimated from gradient
+        postcondition:
+        vertex values are updated
+        """
+
+
+        N = self.domain.number_of_elements
+        beta = self.domain.beta
+        #beta = 0.8
+
+        self.compute_gradients()
+
+
+        G = self.gradients
+        V = self.domain.vertices
+        C = self.domain.centroids
+        qc = self.centroid_values
+        qv = self.vertex_values
+
+        V0 = V[:,0]
+        V1 = V[:,1]
+
+        # Extrapolate to vertices
+        qv[:,0] = qc + G*(V0-C)
+        qv[:,1] = qc + G*(V1-C)
+
+
+        # Find max and min values
+        self.find_qmax_qmin()
+
+        qmax = self.qmax
+        qmin = self.qmin
+
         #Diffences between centroids and maxima/minima
         dqmax = qmax - qc
@@ -941 +954 @@
         #Deltas between vertex and centroid values
         dq = numpy.zeros(qv.shape, numpy.float)
-        for i in range(2):
-            dq[:,i] = qv[:,i] - qc
-
-        #Phi limiter
-        for k in range(N):
-
-            #Find the gradient limiter (phi) across vertices
-            phi = 1.0
-            for i in range(2):
-                r = 1.0
-                if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
-                if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
-
-                phi = min( min(r*beta, 1), phi )
-
-            #Then update using phi limiter
-            for i in range(2):
-                qv[k,i] = qc[k] + phi*dq[k,i]
+
+        dq[:,0] = qv[:,0] - qc
+        dq[:,1] = qv[:,1] - qc
+
+        phi = numpy.ones(qc.shape, numpy.float)
+
+        r0 = numpy.where(dq[:,0]>0.0,dqmax/dq[:,0],1.0)
+        r0 = numpy.where(dq[:,0]<0.0,dqmin/dq[:,0],r0)
+
+        r1 = numpy.where(dq[:,1]>0.0,dqmax/dq[:,1],1.0)
+        r1 = numpy.where(dq[:,1]<0.0,dqmin/dq[:,1],r1)
+
+        phi = numpy.min(r0*beta,numpy.min(r1*beta,1.0))
+
+        qv[:,0] = qc + phi*dq[:,0]
+        qv[:,1] = qc + phi*dq[:,1]
+
+#        #Phi limiter
+#        for k in range(N):
+#
+#            #Find the gradient limiter (phi) across vertices
+#            phi = 1.0
+#            for i in range(2):
+#                r = 1.0
+#                if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
+#                if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
+#
+#                phi = min( min(r*beta, 1), phi )
+#
+#            #Then update using phi limiter
+#            for i in range(2):
+#                qv[k,i] = qc[k] + phi*dq[k,i]
 
     def limit_range(self):
         import sys
 
-        from util import minmod, minmod_kurganov, maxmod, vanleer, vanalbada
+        from limiters_python import minmod, minmod_kurganov, maxmod, vanleer, vanalbada
 
         limiter = self.domain.limiter
@@ -972 +999 @@
         qc = self.centroid_values
         qv = self.vertex_values
-        C = self.domain.centroids
-        X = self.domain.vertices
+        xc = self.domain.centroids
+        x0 = self.domain.vertices[:,0]
+        x1 = self.domain.vertices[:,1]
+
         beta_p = numpy.zeros(N,numpy.float)
         beta_m = numpy.zeros(N,numpy.float)
         beta_x = numpy.zeros(N,numpy.float)
         
-        for k in range(N):
-        
-            n0 = self.domain.neighbours[k,0]
-            n1 = self.domain.neighbours[k,1]
+#        for k in range(N):
+#
+#            n0 = self.domain.neighbours[k,0]
+#            n1 = self.domain.neighbours[k,1]
+#
+#            if ( n0 >= 0) & (n1 >= 0):
+#                #SLOPE DERIVATIVE LIMIT
+#                beta_p[k] = (qc[k]-qc[k-1])/(C[k]-C[k-1])
+#                beta_m[k] = (qc[k+1]-qc[k])/(C[k+1]-C[k])
+#                beta_x[k] = (qc[k+1]-qc[k-1])/(C[k+1]-C[k-1])
+
+
+
+        n0 = self.domain.neighbours[:,0]
+        n1 = self.domain.neighbours[:,1]
+
+
+        beta_p[1:]  = (qc[1:]-qc[:-1])/(xc[1:]-xc[:-1])
+        beta_m[:-1] = beta_p[1:]
+        beta_x[1:-1] = (qc[2:]-qc[:-2])/(xc[2:]-xc[:-2])
+                
+        dx = numpy.zeros(qv.shape, numpy.float)
+
+        dx[:,0] = x0 - xc
+        dx[:,1] = x1 - xc
+
+
+        phi = numpy.zeros(qc.shape, numpy.float)
+
+        phi[0] = (qc[1] - qc[0])/(xc[1] - xc[0])
+        phi[N-1] = (qc[N-1] - qc[N-2])/(xc[N-1] - xc[N-2])
+
+        if limiter == "minmod":
+            phi[1:-1] = minmod(beta_p[1:-1],beta_m[1:-1])
+
+        elif limiter == "vanleer":
+            phi[1:-1] = vanleer(beta_p[1:-1],beta_m[1:-1])
             
-            if ( n0 >= 0) & (n1 >= 0):
-                #SLOPE DERIVATIVE LIMIT
-                beta_p[k] = (qc[k]-qc[k-1])/(C[k]-C[k-1])
-                beta_m[k] = (qc[k+1]-qc[k])/(C[k+1]-C[k])
-                beta_x[k] = (qc[k+1]-qc[k-1])/(C[k+1]-C[k-1])
-                
-        dq = numpy.zeros(qv.shape, numpy.float)
-        for i in range(2):
-            dq[:,i] =self.domain.vertices[:,i]-self.domain.centroids
-            
+        elif limiter == "vanalbada":
+            phi[1:-1] = vanalbada(beta_p[1:-1],beta_m[1:-1])
+
+        elif limiter == "minmod_kurganov":
+            theta = 2.0
+            phi[1:-1] = minmod_kurganov(theta*beta_p[1:-1],theta*beta_m[1:-1], beta_x[1:-1])
+
+        elif limiter == "superbee":
+
+            slope1 = minmod(beta_m[1:-1],2.0*beta_p[1:-1])
+            slope2 = minmod(2.0*beta_m[1:-1],beta_p[1:-1])
+            phi[1:-1] = maxmod(slope1,slope2)
+
+
+        else:
+            msg = 'Unknown limiter'
+            raise Exception, msg
+    
+        qv[:,0] = qc + phi*dx[:,0]
+        qv[:,1] = qc + phi*dx[:,1]
+
+
         #Phi limiter
-        for k in range(N):
-            n0 = self.domain.neighbours[k,0]
-            n1 = self.domain.neighbours[k,1]
-            if n0 < 0:
-                phi = (qc[k+1] - qc[k])/(C[k+1] - C[k])
-            elif n1 < 0:
-                phi = (qc[k] - qc[k-1])/(C[k] - C[k-1])
-            #elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
-            #    phi = 0.0
-            else:
-                if limiter == "minmod":
-                    phi = minmod(beta_p[k],beta_m[k])
-
-                elif limiter == "minmod_kurganov":#Change this
-                    # Also known as monotonized central difference limiter
-                    # if theta = 2.0
-                    theta = 2.0 
-                    phi = minmod_kurganov(theta*beta_p[k],theta*beta_m[k],beta_x[k])
-                elif limiter == "superbee":
-                    slope1 = minmod(beta_m[k],2.0*beta_p[k])
-                    slope2 = minmod(2.0*beta_m[k],beta_p[k])
-                    phi = maxmod(slope1,slope2)
-
-                elif limiter == "vanleer":
-                    phi = vanleer(beta_p[k],beta_m[k])
-
-                elif limiter == "vanalbada":
-                    phi = vanalbada(beta_m[k],beta_p[k])
-            
-            for i in range(2):
-                qv[k,i] = qc[k] + phi*dq[k,i]
+#        for k in range(N):
+#            n0 = self.domain.neighbours[k,0]
+#            n1 = self.domain.neighbours[k,1]
+#            if n0 < 0:
+#                phi = (qc[k+1] - qc[k])/(xc[k+1] - xc[k])
+#            elif n1 < 0:
+#                phi = (qc[k] - qc[k-1])/(xc[k] - xc[k-1])
+#            #elif (self.domain.wet_nodes[k,0] == 2) & (self.domain.wet_nodes[k,1] == 2):
+#            #    phi = 0.0
+#            else:
+#                if limiter == "minmod":
+#                    phi = minmod(beta_p[k],beta_m[k])
+#
+#                elif limiter == "minmod_kurganov":#Change this
+#                    # Also known as monotonized central difference limiter
+#                    # if theta = 2.0
+#                    theta = 2.0
+#                    phi = minmod_kurganov(theta*beta_p[k],theta*beta_m[k],beta_x[k])
+#
+#                elif limiter == "superbee":
+#                    slope1 = minmod(beta_m[k],2.0*beta_p[k])
+#                    slope2 = minmod(2.0*beta_m[k],beta_p[k])
+#                    phi = maxmod(slope1,slope2)
+#
+#                elif limiter == "vanleer":
+#                    phi = vanleer(beta_p[k],beta_m[k])
+#
+#                elif limiter == "vanalbada":
+#                    phi = vanalbada(beta_m[k],beta_p[k])
+#
+#            for i in range(2):
+#                qv[k,i] = qc[k] + phi*dx[k,i]
 
     def limit_steve_slope(self):    
@@ -1107 +1182 @@
 
 
-def newLinePlot(title='Simple Plot'):
-    import pylab as g
-    g.ion()
-    g.hold(False)
-    g.title(title)
-    g.xlabel('x')
-    g.ylabel('y')
-    
-
-def linePlot(x,y):
-    import pylab as g
-    g.plot(x.flat,y.flat)
-
-
-def closePlots():
-    import pylab as g
-    g.close('all')
+
     
 if __name__ == "__main__":
     #from domain import Domain
     from generic_domain import Generic_domain as Domain
+
+
+    def newLinePlot(title='Simple Plot'):
+        import pylab as g
+        g.ion()
+        g.hold(False)
+        g.title(title)
+        g.xlabel('x')
+        g.ylabel('y')
+
+
+    def linePlot(x,y):
+        import pylab as g
+        g.plot(x.flat,y.flat)
+
+
+    def closePlots():
+        import pylab as g
+        g.close('all')
+
+        
     
     points1 = [0.0, 1.0, 2.0, 3.0]
@@ -1213 +1293 @@
 
     raw_input('press return to quit')
-closePlots()
+
+    closePlots()

anuga_work/development/2010-projects/anuga_1d/base/test_generic_domain.py

@@ -5 +5 @@
 
 
-from anuga_1d.generic.generic_domain import *
+from anuga_1d.base.generic_domain import *
 import numpy
 

anuga_work/development/2010-projects/anuga_1d/base/util.py

@@ -19 +19 @@
     return a
 
-def  minmod(beta_p,beta_m):
-    if (abs(beta_p) < abs(beta_m)) & (beta_p*beta_m > 0.0):
-        phi = beta_p
-    elif (abs(beta_m) < abs(beta_p)) & (beta_p*beta_m > 0.0):
-        phi = beta_m
-    else:
-        phi = 0.0
-    return phi
 
-def  minmod_kurganov(a,b,c):
-    from numpy import sign
-    if sign(a)==sign(b)==sign(c):
-        return sign(a)*min(abs(a),abs(b),abs(c))
-    else:
-        return 0.0
-
-def  maxmod(a,b):
-    if (abs(a) > abs(b)) & (a*b > 0.0):
-        phi = a
-    elif (abs(b) > abs(a)) & (a*b > 0.0):
-        phi = b
-    else:
-        phi = 0.0
-    return phi
-
-def vanleer(a,b):
-    if abs(a)+abs(b) > 1e-12:
-        return (a*abs(b)+abs(a)*b)/(abs(a)+abs(b))
-    else:
-        return 0.0
-
-def vanalbada(a,b):
-    if a*a+b*b > 1e-12:
-        return (a*a*b+a*b*b)/(a*a+b*b)
-    else:
-        return 0.0
-

anuga_work/development/2010-projects/anuga_1d/channel/channel_domain.py

@@ -53 +53 @@
 
 
-from anuga_1d.generic.generic_domain import *
+from anuga_1d.base.generic_domain import *
 import numpy
 
@@ -95 +95 @@
         
         #Reduction operation for get_vertex_values
-        from anuga_1d.generic.util import mean
+        from anuga_1d.base.util import mean
         self.reduction = mean
         #self.reduction = min  #Looks better near steep slopes

anuga_work/development/2010-projects/anuga_1d/channel/channel_flow_1_padarn.py

@@ -6 +6 @@
 from anuga_1d.channel.channel_domain import *
 from anuga_1d.config import g, epsilon
-from anuga_1d.generic.generic_mesh import interval_mesh
+from anuga_1d.base.generic_mesh import interval_mesh
 
 

anuga_work/development/2010-projects/anuga_1d/channel/profile_channel.py

@@ -15 +15 @@
 from anuga_1d.channel.channel_domain import *
 from anuga_1d.config import g, epsilon
-from anuga_1d.generic.generic_mesh import interval_mesh
+from anuga_1d.base.generic_mesh import uniform_mesh
 
 
@@ -54 +54 @@
 
     # Create domain with centroid points as defined above
-    domain = Domain(*interval_mesh(N))
+    domain = Domain(*uniform_mesh(N))
 
 
@@ -72 +72 @@
     domain.set_CFL(1.0)
     domain.set_limiter("vanleer")
+    #domain.set_limiter("minmod")
     #domain.h0=0.0001
 

anuga_work/development/2010-projects/anuga_1d/compile_all.py

@@ -9 +9 @@
 
 os.chdir('..')
-os.chdir('generic')
+os.chdir('base')
 execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
 

anuga_work/development/2010-projects/anuga_1d/sww/sww_boundary_conditions.py

@@ -7 +7 @@
 __date__ ="$05/06/2010 5:44:05 PM$"
 
-from anuga_1d.generic.generic_domain import *
+from anuga_1d.base.generic_domain import *
 
 class Dirichlet_boundary(Boundary):

anuga_work/development/2010-projects/anuga_1d/sww/sww_domain.py

@@ -45 +45 @@
 import numpy
 
-from anuga_1d.generic.generic_domain import Generic_domain
+from anuga_1d.base.generic_domain import Generic_domain
 from sww_boundary_conditions import *
 from sww_forcing_terms import *
@@ -84 +84 @@
         
         #Reduction operation for get_vertex_values
-        from anuga_1d.generic.util import mean
+        from anuga_1d.base.util import mean
         self.reduction = mean
         #self.reduction = min  #Looks better near steep slopes

anuga_work/development/2010-projects/anuga_1d/sww/sww_forcing_terms.py

@@ -14 +14 @@
     """
 
-    from anuga_1d.generic.util import gradient
+    from anuga_1d.base.util import gradient
 
     xmom  = domain.quantities['xmomentum'].explicit_update

anuga_work/development/2010-projects/anuga_1d/sww/sww_vel_domain.py

@@ -44 +44 @@
 
 
-from anuga_1d.generic.generic_domain import *
+from anuga_1d.base.generic_domain import *
 from sww_boundary_conditions import *
 from sww_forcing_terms import *
@@ -82 +82 @@
         
         #Reduction operation for get_vertex_values
-        from anuga_1d.generic.util import mean
+        from anuga_1d.base.util import mean
         self.reduction = mean
         #self.reduction = min  #Looks better near steep slopes