Changeset 4990

Timestamp:

Feb 5, 2008, 5:25:32 PM (17 years ago)

Author:

steve

Message:

Combining calculation of gradient, limiting edge by neighbour
and extrapolation in one procedure in quantity

Location:

anuga_core/source/anuga/abstract_2d_finite_volumes

Files:

• quantity.py (modified) (1 diff)
• quantity_ext.c (modified) (9 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -1404 +1404 @@
         extrapolate_from_gradient(self)
 
+    def extrapolate_second_order_and_limit(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        extrapolate_second_order_and_limit(self)        
+
     def backup_centroid_values(self):
         #Call correct module function

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

@@ -16 +16 @@
 //Shared code snippets
 #include "util_ext.h"
-//#include "quantity_ext.h" //Obsolete
-
-
+
+
+//-------------------------------------------
+// Low level routines (called from wrappers)
+//------------------------------------------
 
 int _compute_gradients(int N,
@@ -90 +92 @@
         return 0;
 }
-
 
 
@@ -135 +136 @@
         return 0;
 }
+
+
+int _extrapolate_and_limit_from_gradient(int N,double beta,
+                                         double* centroids,
+                                         long*   neighbours,
+                                         double* centroid_values,
+                                         double* vertex_coordinates,
+                                         double* vertex_values,
+                                         double* edge_values,
+                                         double* x_gradient,
+                                         double* y_gradient) {
+
+        int i, k, k2, k3, k6;
+        double x, y, x0, y0, x1, y1, x2, y2;
+        long n;
+        double qmin, qmax, qn, qc;
+        double dq, dqa[3], phi, r;
+
+
+        for (k=0; k<N; k++){
+                k6 = 6*k;
+                k3 = 3*k;
+                k2 = 2*k;
+
+                // Centroid coordinates
+                x = centroids[k2]; y = centroids[k2+1];
+
+                // vertex coordinates
+                // x0, y0, x1, y1, x2, y2 = X[k,:]
+                x0 = vertex_coordinates[k6 + 0];
+                y0 = vertex_coordinates[k6 + 1];
+                x1 = vertex_coordinates[k6 + 2];
+                y1 = vertex_coordinates[k6 + 3];
+                x2 = vertex_coordinates[k6 + 4];
+                y2 = vertex_coordinates[k6 + 5];
+
+                //Centroid Value
+                qc = centroid_values[k];
+
+                // Extrapolate to Vertices
+                vertex_values[k3+0] = centroid_values[k] + x_gradient[k]*(x0-x) + y_gradient[k]*(y0-y);
+                vertex_values[k3+1] = centroid_values[k] + x_gradient[k]*(x1-x) + y_gradient[k]*(y1-y);
+                vertex_values[k3+2] = centroid_values[k] + x_gradient[k]*(x2-x) + y_gradient[k]*(y2-y);
+
+                // Extrapolate to Edges (midpoints)
+                edge_values[k3+0] = 0.5*(vertex_values[k3 + 1]+vertex_values[k3 + 2]);
+                edge_values[k3+1] = 0.5*(vertex_values[k3 + 2]+vertex_values[k3 + 0]);
+                edge_values[k3+2] = 0.5*(vertex_values[k3 + 0]+vertex_values[k3 + 1]);
+
+
+                phi = 1.0;
+
+                for (i=0; i<3; i++) {
+                    dq = edge_values[k3+i] - qc;      //Delta between edge and centroid values
+                    dqa[i] = dq;                      //Save dq for use in updating vertex values
+
+                    n = neighbours[k3+i];
+                    if (n >= 0) {
+                        qn = centroid_values[n];     //Neighbour's centroid value
+
+                        qmin = min(qc, qn);
+                        qmax = max(qc, qn);
+
+
+                        r = 1.0;
+      
+                        if (dq > 0.0) r = (qmax - qc)/dq;
+                        if (dq < 0.0) r = (qmin - qc)/dq;      
+                    
+                        phi = min( min(r*beta, 1.0), phi);    
+                    }
+                }
+
+
+                //Update gradient, edge and vertex values using phi limiter
+                x_gradient[k] = x_gradient[k]*phi;
+                y_gradient[k] = y_gradient[k]*phi;
+
+                edge_values[k3+0] = qc + phi*dqa[0];
+                edge_values[k3+1] = qc + phi*dqa[1];
+                edge_values[k3+2] = qc + phi*dqa[2];
+
+
+                vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
+                vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
+                vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
+
+        }
+
+        return 0;
+
+}
+
+
 
 
@@ -542 +637 @@
 
 
-/////////////////////////////////////////////////
-// Gateways to Python
+//-----------------------------------------------------
+// Python method Wrappers 
+//-----------------------------------------------------
+
 PyObject *update(PyObject *self, PyObject *args) {
   // FIXME (Ole): It would be great to turn this text into a Python DOC string
@@ -937 +1034 @@
         Py_DECREF(domain);
 
-        //Allocate space for return vectors a and b (don't DECREF)
-        //dimensions[0] = N;
-        //a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
-        //b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
-
-        //FIXME: Odd that I couldn't use normal arrays
-        //Allocate space for return vectors a and b (don't DECREF)
-        //a = (double*) malloc(N * sizeof(double));
-        //if (!a) return NULL;
-        //b = (double*) malloc(N * sizeof(double));
-        //if (!b) return NULL;
-
-
-/*      err = _compute_gradients(N, */
-/*                      (double*) centroids -> data, */
-/*                      (double*) centroid_values -> data, */
-/*                      (long*) number_of_boundaries -> data, */
-/*                      (long*) surrogate_neighbours -> data, */
-/*                      (double*) x_gradient -> data, */
-/*                      (double*) y_gradient -> data); */
-
-/*      if (err != 0) { */
-/*        PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed"); */
-/*        return NULL; */
-/*      } */
-
         err = _extrapolate_from_gradient(N,
                         (double*) centroids -> data,
@@ -993 +1064 @@
 
         return Py_BuildValue("");
+}
+
+
+PyObject *extrapolate_second_order_and_limit(PyObject *self, PyObject *args) {
+    /* Compute edge values using second order approximation and limit values 
+       so that edge values are limited by the two corresponding centroid values
+       
+       Python Call:
+       extrapolate_second_order_and_limit(domain,quantity,beta)
+    */
+
+    PyObject *quantity, *domain;
+        
+    PyArrayObject
+        *domain_centroids,              //Coordinates at centroids
+        *domain_vertex_coordinates,     //Coordinates at vertices
+        *domain_number_of_boundaries,   //Number of boundaries for each triangle
+        *domain_surrogate_neighbours,   //True neighbours or - if one missing - self
+        *domain_neighbours,             //True neighbours, or if negative a link to boundary
+
+        *quantity_centroid_values,      //Values at centroids   
+        *quantity_vertex_values,        //Values at vertices
+        *quantity_edge_values,          //Values at edges
+        *quantity_x_gradient,           //x gradient
+        *quantity_y_gradient;           //y gradient
+    
+
+  // Local variables
+  int ntri;
+  double beta;
+  int err;
+    
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOd", 
+                        &domain,
+                        &quantity,
+                        &beta)) {
+      PyErr_SetString(PyExc_RuntimeError, 
+                      "quantity_ext.c: extrapolate_second_order_and_limit could not parse input");      
+      return NULL;
+  }
+
+        
+  // Get pertinent variables
+  domain_centroids              = get_consecutive_array(domain,   "centroid_coordinates");
+  domain_surrogate_neighbours   = get_consecutive_array(domain,   "surrogate_neighbours");
+  domain_number_of_boundaries   = get_consecutive_array(domain,   "number_of_boundaries");
+  domain_vertex_coordinates     = get_consecutive_array(domain,   "vertex_coordinates");
+  domain_neighbours             = get_consecutive_array(domain,   "neighbours");
+
+  quantity_centroid_values      = get_consecutive_array(quantity, "centroid_values");
+  quantity_vertex_values        = get_consecutive_array(quantity, "vertex_values");
+  quantity_edge_values          = get_consecutive_array(quantity, "edge_values");
+  quantity_x_gradient           = get_consecutive_array(quantity, "x_gradient");
+  quantity_y_gradient           = get_consecutive_array(quantity, "y_gradient");
+
+  ntri = quantity_centroid_values -> dimensions[0];
+
+
+
+  err = _compute_gradients(ntri,
+                           (double*) domain_centroids -> data,
+                           (double*) quantity_centroid_values -> data,
+                           (long*)   domain_number_of_boundaries -> data,
+                           (long*)   domain_surrogate_neighbours -> data,
+                           (double*) quantity_x_gradient -> data,
+                           (double*) quantity_y_gradient -> data);
+
+  if (err != 0) {
+      PyErr_SetString(PyExc_RuntimeError,
+                          "quantity_ext.c: Internal function _compute_gradient failed");
+      return NULL;
+  }
+
+
+  err = _extrapolate_and_limit_from_gradient(ntri, beta,
+                                             (double*) domain_centroids -> data,
+                                             (long*)   domain_neighbours -> data,
+                                             (double*) quantity_centroid_values -> data,
+                                             (double*) domain_vertex_coordinates -> data,
+                                             (double*) quantity_vertex_values -> data,
+                                             (double*) quantity_edge_values -> data,
+                                             (double*) quantity_x_gradient -> data,
+                                             (double*) quantity_y_gradient -> data);
+
+  if (err != 0) {
+      PyErr_SetString(PyExc_RuntimeError,
+                          "quantity_ext.c: Internal function _extrapolate_and_limit_from_gradient failed");
+      return NULL;
+  }
+
+
+  // Release
+  Py_DECREF(domain_centroids);
+  Py_DECREF(domain_surrogate_neighbours);
+  Py_DECREF(domain_number_of_boundaries);
+  Py_DECREF(domain_vertex_coordinates);
+
+  Py_DECREF(quantity_centroid_values);
+  Py_DECREF(quantity_vertex_values);
+  Py_DECREF(quantity_edge_values);
+  Py_DECREF(quantity_x_gradient);
+  Py_DECREF(quantity_y_gradient);
+
+  return Py_BuildValue("");
 }
 
@@ -1046 +1222 @@
         Py_DECREF(domain);
 
-        //Allocate space for return vectors a and b (don't DECREF)
-        //dimensions[0] = N;
-        //a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
-        //b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
-
-        //FIXME: Odd that I couldn't use normal arrays
-        //Allocate space for return vectors a and b (don't DECREF)
-        //a = (double*) malloc(N * sizeof(double));
-        //if (!a) return NULL;
-        //b = (double*) malloc(N * sizeof(double));
-        //if (!b) return NULL;
-
 
         err = _compute_gradients(N,
@@ -1071 +1235 @@
           return NULL;
         }
-
-/*      err = _extrapolate(N, */
-/*                      (double*) centroids -> data, */
-/*                      (double*) centroid_values -> data, */
-/*                      (double*) vertex_coordinates -> data, */
-/*                      (double*) vertex_values -> data, */
-/*                      (double*) edge_values -> data, */
-/*                      (double*) x_gradient -> data, */
-/*                      (double*) y_gradient -> data); */
-
-
-/*      if (err != 0) { */
-/*        PyErr_SetString(PyExc_RuntimeError, */
-/*                        "Internal function _extrapolate failed"); */
-/*        return NULL; */
-/*      } */
 
 
@@ -1558 +1706 @@
         {"extrapolate_from_gradient", extrapolate_from_gradient,
                 METH_VARARGS, "Print out"},
+        {"extrapolate_second_order_and_limit", extrapolate_second_order_and_limit,
+                METH_VARARGS, "Print out"},
         {"interpolate_from_vertices_to_edges",
                 interpolate_from_vertices_to_edges,