Changeset 4769 for anuga_core/source/anuga/shallow_water/shallow_water_ext.c

Timestamp:

Oct 30, 2007, 5:13:17 PM (16 years ago)

Author:

ole

Message:

Retired all Python code which is superseded by C implementations.
This has made the code base much smaller and there is less confusion when people try to understand the algorithms.

File:

• anuga_core/source/anuga/shallow_water/shallow_water_ext.c (modified) (13 diffs)

anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -194 +194 @@
 
 // Innermost flux function (using stage w=z+h)
-int flux_function_central(double *q_left, double *q_right,
-                          double z_left, double z_right,
-                          double n1, double n2,
-                          double epsilon, double H0, double g,
-                          double *edgeflux, double *max_speed) {
+int _flux_function_central(double *q_left, double *q_right,
+                           double z_left, double z_right,
+                           double n1, double n2,
+                           double epsilon, double H0, double g,
+                           double *edgeflux, double *max_speed) {
 
   /*Compute fluxes between volumes for the shallow water wave equation
@@ -737 +737 @@
 ///////////////////////////////////////////////////////////////////
 // Gateways to Python
+
+
+
+PyObject *flux_function_central(PyObject *self, PyObject *args) {
+  //
+  // Gateway to innermost flux function.
+  // This is only used by the unit tests as the c implementation is
+  // normally called by compute_fluxes in this module.
+
+
+  PyArrayObject *normal, *ql, *qr,  *edgeflux;
+  double g, epsilon, max_speed, H0, zl, zr;
+
+  if (!PyArg_ParseTuple(args, "OOOddOddd",
+                        &normal, &ql, &qr, &zl, &zr, &edgeflux,
+                        &epsilon, &g, &H0)) {
+    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: flux_function_central could not parse input arguments");
+    return NULL;
+  }
+
+  
+  _flux_function_central((double*) ql -> data, 
+                         (double*) qr -> data, 
+                         zl, 
+                         zr,                                             
+                         ((double*) normal -> data)[0],                                                                          
+                         ((double*) normal -> data)[1],                  
+                         epsilon, H0, g,
+                         (double*) edgeflux -> data, 
+                         &max_speed);
+  
+  return Py_BuildValue("d", max_speed);  
+}
+
+
 
 PyObject *gravity(PyObject *self, PyObject *args) {
@@ -1494 +1529 @@
   }  
 
-  //Input checks (convert sequences into numeric arrays)
+  // Input checks (convert sequences into numeric arrays)
   q = (PyArrayObject *)
     PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
@@ -1506 +1541 @@
   }
 
-  //Allocate space for return vector r (don't DECREF)
+  // Allocate space for return vector r (don't DECREF)
   dimensions[0] = 3;
   r = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
 
-  //Copy
+  // Copy
   for (i=0; i<3; i++) {
     ((double *) (r -> data))[i] = ((double *) (q -> data))[i];
   }
 
-  //Get normal and direction
+  // Get normal and direction
   n1 = ((double *) normal -> data)[0];
   n2 = ((double *) normal -> data)[1];
   if (direction == -1) n2 = -n2;
 
-  //Rotate
+  // Rotate
   _rotate((double *) r -> data, n1, n2);
 
-  //Release numeric arrays
+  // Release numeric arrays
   Py_DECREF(q);
   Py_DECREF(normal);
 
-  //return result using PyArray to avoid memory leak
+  // Return result using PyArray to avoid memory leak
   return PyArray_Return(r);
 }
@@ -1642 +1677 @@
 
       // Edge flux computation (triangle k, edge i)
-      flux_function_central(ql, qr, zl, zr,
-                            normals[ki2], normals[ki2+1],
-                            epsilon, H0, g,
-                            edgeflux, &max_speed);
+      _flux_function_central(ql, qr, zl, zr,
+                             normals[ki2], normals[ki2+1],
+                             epsilon, H0, g,
+                             edgeflux, &max_speed);
       
       
@@ -1898 +1933 @@
     *xmom_explicit_update,
     *ymom_explicit_update,
-    *already_computed_flux;//tracks whether the flux across an edge has already been computed
-
-
-  //Local variables
+    *already_computed_flux; // Tracks whether the flux across an edge has already been computed
+
+
+  // Local variables
   double timestep, max_speed, epsilon, g, H0;
   double normal[2], ql[3], qr[3], zl, zr;
@@ -2060 +2095 @@
 
 PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) {
+  // Compute linear combination between stage as computed by
+  // gradient-limiters limiting using w, and stage computed by
+  // gradient-limiters limiting using h (h-limiter).
+  // The former takes precedence when heights are large compared to the
+  // bed slope while the latter takes precedence when heights are
+  // relatively small.  Anything in between is computed as a balanced
+  // linear combination in order to avoid numerical disturbances which
+  // would otherwise appear as a result of hard switching between
+  // modes.
   //
   //    balance_deep_and_shallow(beta_h, wc, zc, wv, zv,
@@ -2160 +2204 @@
   int k, i, n, N, k3;
   int dimensions[2];
-  double beta_h; //Safety factor (see config.py)
+  double beta_h; // Safety factor (see config.py)
   double *hmin, *hmax, hn;
 
@@ -2201 +2245 @@
       n = ((long*) neighbours -> data)[k3+i];
 
-      //Initialise hvbar with values from hv
+      // Initialise hvbar with values from hv
       ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
 
       if (n >= 0) {
-        hn = ((double*) hc -> data)[n]; //Neighbour's centroid value
+        hn = ((double*) hc -> data)[n]; // Neighbour's centroid value
 
         hmin[k] = min(hmin[k], hn);
@@ -2216 +2260 @@
   _limit(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax);
 
-  // // //Py_DECREF(domain); //FIXME: NEcessary?
+  // // //Py_DECREF(domain); //FIXME: Necessary?
   free(hmin);
   free(hmax);
 
-  //return result using PyArray to avoid memory leak
+  // Return result using PyArray to avoid memory leak
   return PyArray_Return(hvbar);
-  //return Py_BuildValue("");
 }
 
@@ -2337 +2380 @@
 
 
-//////////////////////////////////////////
+//-------------------------------
 // Method table for python module
+//-------------------------------
 static struct PyMethodDef MethodTable[] = {
   /* The cast of the function is necessary since PyCFunction values
@@ -2351 +2395 @@
   {"gravity", gravity, METH_VARARGS, "Print out"},
   {"manning_friction", manning_friction, METH_VARARGS, "Print out"},
+  {"flux_function_central", flux_function_central, METH_VARARGS, "Print out"},  
   {"balance_deep_and_shallow", balance_deep_and_shallow,
    METH_VARARGS, "Print out"},
@@ -2360 +2405 @@
   {"assign_windfield_values", assign_windfield_values,
    METH_VARARGS | METH_KEYWORDS, "Print out"},
-  //{"distribute_to_vertices_and_edges",
-  // distribute_to_vertices_and_edges, METH_VARARGS},
-  //{"update_conserved_quantities",
-  // update_conserved_quantities, METH_VARARGS},
-  //{"set_initialcondition",
-  // set_initialcondition, METH_VARARGS},
   {NULL, NULL}
 };