Changeset 7868 for anuga_work/development/2010-projects/anuga_1d/base/quantity.py

Timestamp:

Jun 22, 2010, 5:30:32 PM (13 years ago)

Author:

steve

Message:

Added in some more c based limiters

File:

• anuga_work/development/2010-projects/anuga_1d/base/quantity.py (modified) (9 diffs)

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

@@ -70 +70 @@
         #flux calculations and forcing functions
         
-        N = domain.number_of_elements
+        self.N = domain.number_of_elements
+        N = self.N
         self.explicit_update = numpy.zeros(N, numpy.float )
         self.semi_implicit_update = numpy.zeros(N, numpy.float )
@@ -78 +79 @@
         self.qmin = numpy.zeros(self.centroid_values.shape, numpy.float)
 
-        self.beta = domain.beta
+        #These are taken from domain but can be set for each quantity
+        # if so desired
+        self.beta = 0.0
+        self.limiter = 'vanleer'
 
 
@@ -604 +608 @@
         return integral
 
+    def get_beta(self,beta):
+        """Get limiting parameter
+        """
+        
+        return self.beta
+
+    def set_beta(self,beta):
+        """Set limiting parameter
+        """
+
+        #Probably should test that it is not too large
+        self.beta = beta
+
+
+    def get_limiter(self):
+        return self.limiter
+
+    def set_limiter(self,limiter):
+
+        possible_limiters = \
+        ['pyvolution', 'minmod_steve', 'minmod', 'minmod_kurganov', 'superbee', 'vanleer', 'vanalbada']
+
+        if limiter in possible_limiters:
+            self.limiter = limiter
+            return
+
+        msg = '%s is an incorrect limiter type.\n'% limiter
+        msg += 'Possible types are: '+ ", ".join(["%s" % el for el in possible_limiters])
+        raise Exception, msg
 
     def update(self, timestep):
@@ -795 +828 @@
 
         #Gradient
-        G[1:-1] = xmic( self.domain.beta, d1, d2 )      
+        G[1:-1] = xmic( self.beta, d1, d2 )      
 
         
@@ -818 +851 @@
         """
 
-        if self.domain.limiter == "pyvolution":
+        if self.limiter == "pyvolution":
             self.limit_pyvolution()
 
-        elif self.domain.limiter == "minmod_steve":
+        elif self.limiter == "minmod_steve":
             self.limit_minmod()
 
@@ -985 +1018 @@
         from limiters_python import minmod, minmod_kurganov, minmod_kurganov_old, maxmod, vanleer, vanalbada
 
-        limiter = self.domain.limiter
+        limiter = self.get_limiter()
         #print limiter
         
         #print 'limit_range'
-        N = self.domain.number_of_elements
+        N = self.N
         qc = self.centroid_values
         qv = self.vertex_values
@@ -1013 +1046 @@
         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])
@@ -1023 +1057 @@
 
         elif limiter == "minmod_kurganov":
-            theta = 2.0
+            theta = self.beta
             phi[1:-1] = minmod_kurganov(theta*beta_p[1:-1],theta*beta_m[1:-1], beta_x[1:-1])
 
@@ -1034 +1068 @@
             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])/(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):