Changeset 1509

Timestamp:

Jun 9, 2005, 4:09:13 PM (19 years ago)

Author:

matthew

Message:

The new method h_limiter_sw is an alternative to h_limiter that does not use malloc and free. It is defined as a PyObject?* method in shallow_water_ext.c. In shallow_water.py, we import h_limiter_sw as h_limiter so that the rest of the python code effectively uses h_limiter_sw instead of h_limiter.

For run_profile.py, with N=128, the previous version of h_limiter consumed 5.48 seconds, averaged over three runs. Importing h_limiter_sw as h_limiter, h_limiter now consumes 4.42 seconds.

Location:

inundation/ga/storm_surge/pyvolution

Files:

• shallow_water.py (modified) (1 diff)
• shallow_water_ext.c (modified) (2 diffs)

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -775 +775 @@
 
     #Call C-extension
-    from shallow_water_ext import h_limiter
+    from shallow_water_ext import h_limiter_sw as h_limiter
     hvbar = h_limiter(domain, hc, hv)
 

inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

@@ -1178 +1178 @@
 }
 
-
-
+PyObject *h_limiter_sw(PyObject *self, PyObject *args) {
+  //a faster version of h_limiter above
+  PyObject *domain, *Tmp;
+  PyArrayObject
+    *hv, *hc, //Depth at vertices and centroids
+    *hvbar,   //Limited depth at vertices (return values)
+    *neighbours;
+
+  int k, i, N, k3,k0,k1,k2;
+  int dimensions[2];
+  double beta_h; //Safety factor (see config.py)
+  double hmin, hmax, dh[3];
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv))
+    return NULL;
+  neighbours = get_consecutive_array(domain, "neighbours");
+  
+  //Get safety factor beta_h
+  Tmp = PyObject_GetAttrString(domain, "beta_h");
+  if (!Tmp)
+    return NULL;
+  beta_h = PyFloat_AsDouble(Tmp);
+
+  Py_DECREF(Tmp);
+
+  N = hc -> dimensions[0];
+
+  //Create hvbar
+  dimensions[0] = N;
+  dimensions[1] = 3;
+  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
+  for (k=0;k<N;k++){
+    k3=k*3;
+    //get the ids of the neighbours
+    k0 = ((long*) neighbours -> data)[k3];
+    k1 = ((long*) neighbours -> data)[k3+1];
+    k2 = ((long*) neighbours -> data)[k3+2];
+    //set hvbar provisionally
+    for (i=0;i<3;i++){
+      ((double*) hvbar -> data)[k3+i] = ((double*) hv -> data)[k3+i];
+      dh[i]=((double*) hvbar -> data)[k3+i]-((double*) hc -> data)[k];
+    }
+    hmin=((double*) hc -> data)[k];
+    hmax=hmin;
+    if (k0>=0){
+      hmin=min(hmin,((double*) hc -> data)[k0]);
+      hmax=max(hmax,((double*) hc -> data)[k0]);
+    }
+    if (k1>=0){
+      hmin=min(hmin,((double*) hc -> data)[k1]);
+      hmax=max(hmax,((double*) hc -> data)[k1]);
+    }
+    if (k2>=0){
+      hmin=min(hmin,((double*) hc -> data)[k2]);
+      hmax=max(hmax,((double*) hc -> data)[k2]);
+    }
+    hmin-=((double*) hc -> data)[k];
+    hmax-=((double*) hc -> data)[k];
+    limit_gradient(dh,hmin,hmax,beta_h);
+    for (i=0;i<3;i++)
+      ((double*) hvbar -> data)[k3+i] = ((double*) hc -> data)[k]+dh[i];
+  }
+  return PyArray_Return(hvbar);
+}
 
 PyObject *assign_windfield_values(PyObject *self, PyObject *args) {
@@ -1237 +1299 @@
    METH_VARARGS, "Print out"},
   {"h_limiter", h_limiter,
+   METH_VARARGS, "Print out"},
+  {"h_limiter_sw", h_limiter_sw,
    METH_VARARGS, "Print out"},
   {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},