Changeset 255

Timestamp:

Aug 31, 2004, 11:47:24 AM (21 years ago)

Author:

ole

Message:

 

Location:

inundation/ga/storm_surge/pyvolution

Files:

• quantity.py (modified) (2 diffs)
• quantity_ext.c (modified) (4 diffs)

inundation/ga/storm_surge/pyvolution/quantity.py

@@ -352 +352 @@
             
     def extrapolate_second_order(self):
-        """Extrapolate conserved quantities from centroid to
-        vertices for each volume using
-        second order scheme.
-        """
-        
-        a, b = self.compute_gradients()
-
-        V = self.domain.get_vertex_coordinates()
-        qc = self.centroid_values
-        qv = self.vertex_values
-        
-        #Check each triangle
-        for k in range(self.domain.number_of_elements):
-            #Centroid coordinates            
-            x, y = self.domain.centroids[k]
-
-            #vertex coordinates
-            x0, y0, x1, y1, x2, y2 = V[k,:]
-
-            #Extrapolate
-            qv[k,0] = qc[k] + a[k]*(x0-x) + b[k]*(y0-y)
-            qv[k,1] = qc[k] + a[k]*(x1-x) + b[k]*(y1-y)
-            qv[k,2] = qc[k] + a[k]*(x2-x) + b[k]*(y2-y)            
+        #Call correct module function
+        #(either from this module or C-extension)
+        extrapolate_second_order(self)
+
+
+def extrapolate_second_order(self):        
+    """Extrapolate conserved quantities from centroid to
+    vertices for each volume using
+    second order scheme.
+    """
+        
+    a, b = self.compute_gradients()
+
+    V = self.domain.get_vertex_coordinates()
+    qc = self.centroid_values
+    qv = self.vertex_values
+    
+    #Check each triangle
+    for k in range(self.domain.number_of_elements):
+        #Centroid coordinates            
+        x, y = self.domain.centroids[k]
+        
+        #vertex coordinates
+        x0, y0, x1, y1, x2, y2 = V[k,:]
+        
+        #Extrapolate
+        qv[k,0] = qc[k] + a[k]*(x0-x) + b[k]*(y0-y)
+        qv[k,1] = qc[k] + a[k]*(x1-x) + b[k]*(y1-y)
+        qv[k,2] = qc[k] + a[k]*(x2-x) + b[k]*(y2-y)            
 
             
@@ -446 +452 @@
     #Replace python version with c implementations
 
-    pass
-    #from quantity_ext import limit
-
+    from quantity_ext import limit #, extrapolate_second_order        

inundation/ga/storm_surge/pyvolution/quantity_ext.c

@@ -31 +31 @@
 
 
-void _limit(int N, double beta, double* qc, double* qv, double* qmin, double* qmax) { 
-
-  int k, i, ki;
-  double dq, phi, r;
+void _limit(int N, double beta, double* qc, double* qv, 
+            double* qmin, double* qmax) { 
+
+  int k, i, k3;
+  double dq, dqa[3], phi, r;
   
   for (k=0; k<N; k++) {
-  
+    k3 = k*3;
+    
     //Find the gradient limiter (phi) across vertices  
     phi = 1.0;
     for (i=0; i<3; i++) {    
-      ki = k*3+i;
       r = 1.0;
       
-      dq = qv[ki] - qc[k];    //Delta between vertex and centroid values
-    
+      dq = qv[k3+i] - qc[k];    //Delta between vertex and centroid values
+      dqa[i] = dq;              //Save dq for use in the next loop
+      
       if (dq > 0.0) r = (qmax[k] - qc[k])/dq;
       if (dq < 0.0) r = (qmin[k] - qc[k])/dq;      
   
+  
       phi = min( min(r*beta, 1.0), phi);    
     }
+    
     //Then update using phi limiter
     for (i=0; i<3; i++) {    
-      qv[ki] = qc[k] + phi*(qv[ki] - qc[k]); //FIXME: Could precompute dq (3 elements)
+      qv[k3+i] = qc[k] + phi*dqa[i];
     }
   }
-
-     
-  /*            
-    #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]
-            
-  */
-  
-}
-
-// Gateway to Python
-//FIXME: DOES NOT WORK YET
-PyObject *limit(PyObject *self, PyObject *args) {
-  
-  PyObject *quantity;
-  
-  PyObject *Tmp;
-  
+}
+
+/////////////////////////////////////////////////
+// Gateways to Python
+
+/*
+PyObject *extrapolate_second_order(PyObject *self, PyObject *args) {
+
+  PyObject *quantity, *domain, *Tmp;
   PyArrayObject 
     *qv, //Conserved quantities at vertices 
@@ -101 +74 @@
     *neighbours;
     
-  int k, i, n, N;
+  int k, i, n, N, k3;
   double beta; //Safety factor
   double *qmin, *qmax, qn;
@@ -108 +81 @@
   if (!PyArg_ParseTuple(args, "O", &quantity)) 
     return NULL;
-  
-  //Get safety factor beta
-  Tmp = PyObject_GetAttrString(quantity, "beta"); 
-  beta = PyFloat_AsDouble(Tmp);  
-  Py_DECREF(Tmp);         
-  
-  
+
+  domain = PyObject_GetAttrString(quantity, "domain");     
+  if (!domain) 
+    return NULL;  
+    
+  neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours");
+
   qc = (PyArrayObject*) PyObject_GetAttrString(quantity, "centroid_values"); 
   qv = (PyArrayObject*) PyObject_GetAttrString(quantity, "vertex_values");   
-  
-  Tmp = PyObject_GetAttrString(quantity, "domain");     
-  neighbours = (PyArrayObject*) PyObject_GetAttrString(Tmp, "neighbours");     
-  Py_DECREF(Tmp);           
   N = qc -> dimensions[0];
-
-  
+*/
+  /*
+def extrapolate_second_order(self):        
+    """Extrapolate conserved quantities from centroid to
+    vertices for each volume using
+    second order scheme.
+    """
+        
+    a, b = self.compute_gradients()
+
+    V = self.domain.get_vertex_coordinates()
+    qc = self.centroid_values
+    qv = self.vertex_values
+    
+    #Check each triangle
+    for k in range(self.domain.number_of_elements):
+        #Centroid coordinates            
+        x, y = self.domain.centroids[k]
+        
+        #vertex coordinates
+        x0, y0, x1, y1, x2, y2 = V[k,:]
+        
+        #Extrapolate
+        qv[k,0] = qc[k] + a[k]*(x0-x) + b[k]*(y0-y)
+        qv[k,1] = qc[k] + a[k]*(x1-x) + b[k]*(y1-y)
+        qv[k,2] = qc[k] + a[k]*(x2-x) + b[k]*(y2-y)            
+  */
+//}    
+
+
+
+PyObject *limit(PyObject *self, PyObject *args) {
+  
+  PyObject *quantity, *domain, *Tmp;
+  PyArrayObject 
+    *qv, //Conserved quantities at vertices 
+    *qc, //Conserved quantities at centroids
+    *neighbours;
+    
+  int k, i, n, N, k3;
+  double beta; //Safety factor
+  double *qmin, *qmax, qn;
+  
+  // Convert Python arguments to C  
+  if (!PyArg_ParseTuple(args, "O", &quantity)) 
+    return NULL;
+
+  domain = PyObject_GetAttrString(quantity, "domain");     
+  if (!domain) 
+    return NULL;  
+    
+  neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours");
+
+  //Get safety factor beta
+  Tmp = PyObject_GetAttrString(domain, "beta"); 
+  if (!Tmp) 
+    return NULL;
+      
+  beta = PyFloat_AsDouble(Tmp);  
+  
+  Py_DECREF(Tmp);    
+  Py_DECREF(domain);      
+
+  qc = (PyArrayObject*) PyObject_GetAttrString(quantity, "centroid_values"); 
+  qv = (PyArrayObject*) PyObject_GetAttrString(quantity, "vertex_values");   
+  N = qc -> dimensions[0];
+    
   //Find min and max of this and neighbour's centroid values
   qmin = malloc(N * sizeof(double));
   qmax = malloc(N * sizeof(double));  
   for (k=0; k<N; k++) {
-    qmax[k] = qmin[k] = ((double*) qc -> data)[k];
+    k3=k*3;
+    
+    qmin[k] = ((double*) qc -> data)[k];
+    qmax[k] = qmin[k];
+    
     for (i=0; i<3; i++) {
-      n = ((int*) neighbours -> data)[k*3+i];
+      n = ((int*) neighbours -> data)[k3+i];
       if (n >= 0) {
         qn = ((double*) qc -> data)[n]; //Neighbour's centroid value
@@ -141 +179 @@
   
   // Call underlying routine
-  //_limit(N, beta, (double*) qc -> data, (double*) qv -> data, qmin, qmax);
+  _limit(N, beta, (double*) qc -> data, (double*) qv -> data, qmin, qmax);
          
   free(qmin);