Changeset 4886

Timestamp:

Dec 12, 2007, 5:38:37 PM (17 years ago)

Author:

steve

Message:

Working towards an edge limiter (at least in quantity)

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/quantity.py (modified) (3 diffs)
• abstract_2d_finite_volumes/quantity_ext.c (modified) (6 diffs)
• abstract_2d_finite_volumes/test_quantity.py (modified) (5 diffs)
• utilities/util_ext.h (modified) (1 diff)

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -229 +229 @@
         #(either from this module or C-extension)
         interpolate_from_vertices_to_edges(self)
+
+    def interpolate_from_edges_to_vertices(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        interpolate_from_edges_to_vertices(self)
 
 
@@ -1400 +1405 @@
 
 
-##     def limit_by_vertex(self):
-##         #Call correct module function
-##         #(either from this module or C-extension)
-##         limit_by_vertex()
-
-
-##     def limit_by_edge(self):
-##         #Call correct module function
-##         #(either from this module or C-extension)
-##         limit_by_edge()        
+    def limit_by_vertex(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        limit_by_vertex(self)
+
+
+    def limit_by_edge(self):
+        #Call correct module function
+        #(either from this module or C-extension)
+        limit_by_edge(self)        
 
     def extrapolate_second_order(self):
@@ -1439 +1444 @@
          compute_gradients,\
          limit_old,\
+         limit_by_vertex,\
+         limit_by_edge,\
          extrapolate_second_order,\
          interpolate_from_vertices_to_edges,\
+         interpolate_from_edges_to_vertices,\
          update    
 else:

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

@@ -122 +122 @@
                 y2 = vertex_coordinates[k6 + 5];
 
-                // Extrapolate
+                // Extrapolate to Vertices
                 vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y);
                 vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y);
                 vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[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]);
@@ -137 +137 @@
 
 
-
-
-int _interpolate(int N,
+int _limit_by_vertex(int N, double beta,
+                     double* centroid_values,
+                     double* vertex_values,
+                     double* edge_values,
+                     long*   neighbours) {
+
+        int i, k, k2, k3, k6;
+        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;
+
+                qc = centroid_values[k];
+                qmin = qc;
+                qmax = qc;
+
+                for (i=0; i<3; i++) {
+                    n = neighbours[k3+i];
+                    if (n >= 0) {
+                        qn = centroid_values[n]; //Neighbour's centroid value
+
+                        qmin = min(qmin, qn);
+                        qmax = max(qmax, qn);
+                    }
+                }
+
+                phi = 1.0;
+                for (i=0; i<3; i++) {    
+                    r = 1.0;
+      
+                    dq = vertex_values[k3+i] - qc;    //Delta between vertex and centroid values
+                    dqa[i] = dq;                      //Save dq for use in updating vertex values
+      
+                    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 vertex and edge values using phi limiter
+                vertex_values[k3+0] = qc + phi*dqa[0];
+                vertex_values[k3+1] = qc + phi*dqa[1];
+                vertex_values[k3+2] = qc + phi*dqa[2];
+                
+                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]);
+
+        }
+
+        return 0;
+}
+
+int _limit_by_edge(int N, double beta,
+                     double* centroid_values,
+                     double* vertex_values,
+                     double* edge_values,
+                     long*   neighbours) {
+
+        int i, k, k2, k3, k6;
+        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;
+
+                qc = centroid_values[k];
+                qmin = qc;
+                qmax = qc;
+
+                for (i=0; i<3; i++) {
+                    n = neighbours[k3+i];
+                    if (n >= 0) {
+                        qn = centroid_values[n]; //Neighbour's centroid value
+
+                        qmin = min(qmin, qn);
+                        qmax = max(qmax, qn);
+                    }
+                }
+
+                phi = 1.0;
+                for (i=0; i<3; i++) {    
+                    r = 1.0;
+      
+                    dq = edge_values[k3+i] - qc;     //Delta between edge and centroid values
+                    dqa[i] = dq;                      //Save dq for use in updating vertex values
+      
+                    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 vertex and edge values using phi limiter
+                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;
+}
+
+
+
+int _interpolate_from_vertices_to_edges(int N,
                  double* vertex_values,
                  double* edge_values) {
@@ -157 +273 @@
                 edge_values[k3 + 1] = 0.5*(q0+q2);
                 edge_values[k3 + 2] = 0.5*(q0+q1);
+        }
+        return 0;
+}
+
+
+int _interpolate_from_edges_to_vertices(int N,
+                 double* vertex_values,
+                 double* edge_values) {
+
+        int k, k3;
+        double e0, e1, e2;
+
+
+        for (k=0; k<N; k++) {
+                k3 = 3*k;
+
+                e0 = edge_values[k3 + 0];
+                e1 = edge_values[k3 + 1];
+                e2 = edge_values[k3 + 2];
+
+                vertex_values[k3 + 0] = e1 + e2 - e0;
+                vertex_values[k3 + 1] = e2 + e0 - e1;
+                vertex_values[k3 + 2] = e0 + e1 - e2;
         }
         return 0;
@@ -458 +597 @@
         N = vertex_values -> dimensions[0];
 
-        err = _interpolate(N,
+        err = _interpolate_from_vertices_to_edges(N,
+                           (double*) vertex_values -> data,
+                           (double*) edge_values -> data);
+
+        if (err != 0) {
+          PyErr_SetString(PyExc_RuntimeError, 
+                          "Interpolate could not be computed");
+          return NULL;
+        }
+
+        // Release and return
+        Py_DECREF(vertex_values);
+        Py_DECREF(edge_values);
+
+        return Py_BuildValue("");
+}
+
+
+PyObject *interpolate_from_edges_to_vertices(PyObject *self, PyObject *args) {
+        //
+        //Compute vertex values from edge values using linear interpolation
+        // 
+        
+        PyObject *quantity;
+        PyArrayObject *vertex_values, *edge_values;
+
+        int N, err;
+
+        // Convert Python arguments to C
+        if (!PyArg_ParseTuple(args, "O", &quantity)) {
+          PyErr_SetString(PyExc_RuntimeError, 
+                          "quantity_ext.c: interpolate_from_edges_to_vertices could not parse input");
+          return NULL;
+        }
+        
+        vertex_values = get_consecutive_array(quantity, "vertex_values");
+        edge_values = get_consecutive_array(quantity, "edge_values");
+
+        N = vertex_values -> dimensions[0];
+
+        err = _interpolate_from_edges_to_vertices(N,
                            (double*) vertex_values -> data,
                            (double*) edge_values -> data);
@@ -795 +974 @@
 
 
+PyObject *limit_by_vertex(PyObject *self, PyObject *args) {
+  //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 and edge values are updated
+  //
+
+        PyObject *quantity, *domain, *Tmp;
+        PyArrayObject
+            *vertex_values,   //Conserved quantities at vertices
+            *centroid_values, //Conserved quantities at centroids
+            *edge_values,     //Conserved quantities at edges
+            *neighbours;
+
+        double beta_w; //Safety factor
+        int N, err;
+
+
+        // Convert Python arguments to C
+        if (!PyArg_ParseTuple(args, "O", &quantity)) {
+          PyErr_SetString(PyExc_RuntimeError, 
+                          "quantity_ext.c: limit_by_vertex could not parse input");
+          return NULL;
+        }
+
+        domain = PyObject_GetAttrString(quantity, "domain");
+        if (!domain) {
+          PyErr_SetString(PyExc_RuntimeError, 
+                          "quantity_ext.c: limit_by_vertex could not obtain domain object from quantity");                      
+          
+          return NULL;
+        }
+
+        // Get safety factor beta_w
+        Tmp = PyObject_GetAttrString(domain, "beta_w");
+        if (!Tmp) {
+          PyErr_SetString(PyExc_RuntimeError, 
+                          "quantity_ext.c: limit_by_vertex could not obtain beta_w object from domain");                        
+          
+          return NULL;
+        }       
+
+
+        // Get pertinent variables
+        neighbours       = get_consecutive_array(domain, "neighbours");
+        centroid_values  = get_consecutive_array(quantity, "centroid_values");
+        vertex_values    = get_consecutive_array(quantity, "vertex_values");
+        edge_values      = get_consecutive_array(quantity, "edge_values");
+        beta_w           = PyFloat_AsDouble(Tmp);
+
+
+        N = centroid_values -> dimensions[0];
+
+        err = _limit_by_vertex(N, beta_w,
+                               (double*) centroid_values -> data,
+                               (double*) vertex_values -> data,
+                               (double*) edge_values -> data,
+                               (long*)   neighbours -> data);
+
+        if (err != 0) {
+          PyErr_SetString(PyExc_RuntimeError,
+                          "Internal function _limit_by_vertex failed");
+          return NULL;
+        }       
+
+
+        // Release
+        Py_DECREF(neighbours);
+        Py_DECREF(centroid_values);
+        Py_DECREF(vertex_values);
+        Py_DECREF(edge_values);
+        Py_DECREF(Tmp);
+
+
+        return Py_BuildValue("");
+}
+
+
+
+PyObject *limit_by_edge(PyObject *self, PyObject *args) {
+  //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 and edge values are updated
+  //
+
+        PyObject *quantity, *domain, *Tmp;
+        PyArrayObject
+            *vertex_values,   //Conserved quantities at vertices
+            *centroid_values, //Conserved quantities at centroids
+            *edge_values,     //Conserved quantities at edges
+            *neighbours;
+
+        double beta_w; //Safety factor
+        int N, err;
+
+
+        // Convert Python arguments to C
+        if (!PyArg_ParseTuple(args, "O", &quantity)) {
+          PyErr_SetString(PyExc_RuntimeError, 
+                          "quantity_ext.c: limit_by_edge could not parse input");
+          return NULL;
+        }
+
+        domain = PyObject_GetAttrString(quantity, "domain");
+        if (!domain) {
+          PyErr_SetString(PyExc_RuntimeError, 
+                          "quantity_ext.c: limit_by_edge could not obtain domain object from quantity");                        
+          
+          return NULL;
+        }
+
+        // Get safety factor beta_w
+        Tmp = PyObject_GetAttrString(domain, "beta_w");
+        if (!Tmp) {
+          PyErr_SetString(PyExc_RuntimeError, 
+                          "quantity_ext.c: limit_by_edge could not obtain beta_w object from domain");                  
+          
+          return NULL;
+        }       
+
+
+        // Get pertinent variables
+        neighbours       = get_consecutive_array(domain, "neighbours");
+        centroid_values  = get_consecutive_array(quantity, "centroid_values");
+        vertex_values    = get_consecutive_array(quantity, "vertex_values");
+        edge_values      = get_consecutive_array(quantity, "edge_values");
+        beta_w           = PyFloat_AsDouble(Tmp);
+
+
+        N = centroid_values -> dimensions[0];
+
+        err = _limit_by_edge(N, beta_w,
+                               (double*) centroid_values -> data,
+                               (double*) vertex_values -> data,
+                               (double*) edge_values -> data,
+                               (long*)   neighbours -> data);
+
+        if (err != 0) {
+          PyErr_SetString(PyExc_RuntimeError,
+                          "Internal function _limit_by_vertex failed");
+          return NULL;
+        }       
+
+
+        // Release
+        Py_DECREF(neighbours);
+        Py_DECREF(centroid_values);
+        Py_DECREF(vertex_values);
+        Py_DECREF(edge_values);
+        Py_DECREF(Tmp);
+
+
+        return Py_BuildValue("");
+}
+
+
 
 // Method table for python module
 static struct PyMethodDef MethodTable[] = {
         {"limit_old", limit_old, METH_VARARGS, "Print out"},
+        {"limit_by_vertex", limit_by_vertex, METH_VARARGS, "Print out"},
+        {"limit_by_edge", limit_by_edge, METH_VARARGS, "Print out"},
         {"update", update, METH_VARARGS, "Print out"},
         {"backup_centroid_values", backup_centroid_values, METH_VARARGS, "Print out"},
@@ -808 +1158 @@
                 interpolate_from_vertices_to_edges,
                 METH_VARARGS, "Print out"},
+        {"interpolate_from_edges_to_vertices",
+                interpolate_from_edges_to_vertices,
+                METH_VARARGS, "Print out"},
         {"average_vertex_values", average_vertex_values, METH_VARARGS, "Print out"},            
         {NULL, NULL, 0, NULL}   // sentinel

anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

@@ -1445 +1445 @@
 
 
-    def test_limiter(self):
+    def test_vertex_limiter(self):
         quantity = Conserved_quantity(self.mesh4)
 
@@ -1453 +1453 @@
 
         #Limit
-        quantity.limit()
+        quantity.limit_by_vertex()
 
         #Assert that central triangle is limited by neighbours
@@ -1474 +1474 @@
 
 
+
+    def test_edge_limiter(self):
+        quantity = Conserved_quantity(self.mesh4)
+
+        #Create a deliberate overshoot (e.g. from gradient computation)
+        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
+
+
+        #Limit
+        quantity.limit_by_edge()
+
+        #Assert that central triangle is limited by neighbours
+        assert quantity.edge_values[1,0] <= quantity.centroid_values[2]
+        assert quantity.edge_values[1,0] >= quantity.centroid_values[0]
+
+        assert quantity.edge_values[1,1] <= quantity.centroid_values[2]
+        assert quantity.edge_values[1,1] >= quantity.centroid_values[0]
+
+        assert quantity.edge_values[1,2] <= quantity.centroid_values[2]
+        assert quantity.edge_values[1,2] <= quantity.centroid_values[0]
+
+
+
+        #Assert that quantities are conserved
+        from Numeric import sum
+        for k in range(quantity.centroid_values.shape[0]):
+            assert allclose (quantity.centroid_values[k],
+                             sum(quantity.vertex_values[k,:])/3)
+
+
     def test_limiter2(self):
         """Taken from test_shallow_water
@@ -1518 +1548 @@
 
 
-        #Extrapolate
+        #Extrapolate from centroid to vertices and edges
         quantity.extrapolate_first_order()
 
         #Interpolate
-        quantity.interpolate_from_vertices_to_edges()
+        #quantity.interpolate_from_vertices_to_edges()
 
         assert allclose(quantity.vertex_values,
@@ -1528 +1558 @@
         assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
                                                [3,3,3], [4, 4, 4]])
+
+
+    def test_interpolate_from_vertices_to_edges(self):
+        quantity = Quantity(self.mesh4)
+
+        quantity.vertex_values = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]],Float)
+
+        quantity.interpolate_from_vertices_to_edges()
+
+        assert allclose(quantity.edge_values, [[1., 1.5, 0.5],
+                                               [3., 2.5, 1.5],
+                                               [3.5, 4.5, 3.],
+                                               [2.5, 3.5, 2]])
+
+
+    def test_interpolate_from_edges_to_vertices(self):
+        quantity = Quantity(self.mesh4)
+
+        quantity.edge_values = array([[1., 1.5, 0.5],
+                                [3., 2.5, 1.5],
+                                [3.5, 4.5, 3.],
+                                [2.5, 3.5, 2]],Float)
+
+        quantity.interpolate_from_edges_to_vertices()
+
+        assert allclose(quantity.vertex_values,
+                        [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
+
 
 

anuga_core/source/anuga/utilities/util_ext.h

@@ -166 +166 @@
 
 
-void _limit_by_edge(int N, double beta, double* qc, double* qv, 
-            double* qmin, double* qmax) { 
-
-  //N are the number of elements
-  int k, i, k3;
-  double dq, dqa[3], phi, r;
-  
-  //printf("INSIDE\n");
-  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++) {    
-      r = 1.0;
-      
-      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[k3+i] = qc[k] + phi*dqa[i];
-    }
-  }
-}
-
-void _limit_by_vertex(int N, double beta, double* qc, double* qv, 
-            double* qmin, double* qmax) { 
-
-  //N are the number of elements
-  int k, i, k3;
-  double dq, dqa[3], phi, r;
-  
-  //printf("INSIDE\n");
-  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++) {    
-      r = 1.0;
-      
-      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[k3+i] = qc[k] + phi*dqa[i];
-    }
-  }
-}