Changeset 245

Timestamp:

Aug 30, 2004, 6:44:45 PM (21 years ago)

Author:

ole

Message:

Started to work on C implementation of phi limiter

Location:

inundation/ga/storm_surge/pyvolution

Files:

• quantity.py (modified) (2 diffs)
• shallow_water_ext.c (modified) (2 diffs)

inundation/ga/storm_surge/pyvolution/quantity.py

@@ -333 +333 @@
 
     def limit(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
-        """
-
-        from Numeric import zeros, Float
-
-        beta = self.domain.beta
-        
-        qc = self.centroid_values
-        qv = self.vertex_values
-        
-        #Deltas between vertex and centroid values
-        dq = zeros( qv.shape, Float)
-        for i in range(3):
-            dq[:,i] = qv[:,i] - qc
-        
-        #Find min and max of this and neighbour's centroid values
-        qmax = zeros(qc.shape, Float)
-        qmin = zeros(qc.shape, Float)        
-
-        for k in range(self.domain.number_of_elements):
-            qmax[k] = qmin[k] = qc[k]
-            for i in range(3):
-                n = self.domain.neighbours[k,i]
-                if n >= 0:
-                    qn = qc[n] #Neighbour's centroid value
-      
-                    qmin[k] = min(qmin[k], qn)
-                    qmax[k] = max(qmax[k], qn)
-      
-    
-        #Diffences between centroids and maxima/minima
-        dqmax = qmax - qc 
-        dqmin = qmin - qc
-  
-        for k in range(self.domain.number_of_elements):
-
-            #Find the gradient limiter (phi) across vertices
-            phi = 1.0
-            for i in range(3):
-                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(3):            
-                qv[k,i] = qc[k] + phi*dq[k,i]
-
+        #Call correct module function
+        #(either from this module or C-extension)
+        limit(self)
+        
         
     def extrapolate_first_order(self):
@@ -431 +377 @@
 
             
+
+def limit(quantity):        
+    """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
+    """
+
+    from Numeric import zeros, Float
+
+    N = quantity.domain.number_of_elements
+    
+    beta = quantity.domain.beta
+        
+    qc = quantity.centroid_values
+    qv = quantity.vertex_values
+        
+    #Find min and max of this and neighbour's centroid values
+    qmax = zeros(qc.shape, Float)
+    qmin = zeros(qc.shape, Float)        
+    
+    for k in range(N):
+        qmax[k] = qmin[k] = qc[k]
+        for i in range(3):
+            n = quantity.domain.neighbours[k,i]
+            if n >= 0:
+                qn = qc[n] #Neighbour's centroid value
+                
+                qmin[k] = min(qmin[k], qn)
+                qmax[k] = max(qmax[k], qn)
+      
+    
+    #Diffences between centroids and maxima/minima
+    dqmax = qmax - qc 
+    dqmin = qmin - qc
+        
+    #Deltas between vertex and centroid values
+    dq = zeros(qv.shape, Float)
+    for i in range(3):
+        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(3):
+            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(3):            
+            qv[k,i] = qc[k] + phi*dq[k,i]
+
+
+
+import compile
+if compile.can_use_C_extension('quantity_ext.c'):
+    #Replace python version with c implementations
+
+    pass
+    #from quantity_ext import limit
+

inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

@@ -17 +17 @@
 #include "math.h"
 
+//FIXME: Should live in util_ext.c
 double max(double x, double y) {  
   //Return maximum of two doubles
@@ -33 +34 @@
 
 // Computational function for rotation
-// FIXME: Should probably go into shallow_water_ext.c
 int _rotate(double *q, double n1, double n2) {
   /*Rotate the momentum component q (q[1], q[2])