Changeset 7930

Timestamp:

Aug 6, 2010, 3:17:50 PM (14 years ago)

Author:

steve

Message:

Moving directories from contributed to development

Location:

trunk

Files:

• anuga_core/benchmarks/create_test_sww.py (modified) (1 diff)
• anuga_core/documentation/user_manual/anuga_user_manual.pdf (modified) (previous)
• anuga_work/development/2010-projects/anuga_1d/base/limiters_python.py (modified) (1 diff)
• anuga_work/development/2010-projects/anuga_1d/base/quantity_ext.c (modified) (3 diffs)
• anuga_work/development/2010-projects/anuga_1d/base/test_quantity_1d.py (modified) (1 diff)
• anuga_work/development/GIS (moved) (moved from trunk/anuga_core/source/contributed_code/GIS)
• anuga_work/development/flood (moved) (moved from trunk/anuga_core/source/contributed_code/flood)
• anuga_work/development/graphics (moved) (moved from trunk/anuga_core/source/contributed_code/graphics)
• anuga_work/development/tools (moved) (moved from trunk/anuga_core/source/contributed_code/tools)
• anuga_work/development/tsumami (moved) (moved from trunk/anuga_core/source/contributed_code/tsunami)

trunk/anuga_core/benchmarks/create_test_sww.py

@@ -10 +10 @@
         #------------------------------------------------------------------------------
         points, vertices, boundary = anuga.rectangular_cross(5, 5,
-                                                                                                   len1=50.0, len2=50.0) # Mesh
+                                                             len1=50.0, len2=50.0) # Mesh
 
         domain = anuga.Domain(points, vertices, boundary)  # Create domain

trunk/anuga_work/development/2010-projects/anuga_1d/base/limiters_python.py

@@ -46 +46 @@
     from numpy import abs, where
 
-    return where((a*a + b*b >= 1.0e-16), (a*a*b+a*b*b)/(a*a+b*b), 0.0)
+    return where((a*a + b*b >= 1.0e-32), (a*a*b+a*b*b)/(a*a+b*b), 0.0)
 
 

trunk/anuga_work/development/2010-projects/anuga_1d/base/quantity_ext.c

@@ -142 +142 @@
   // Convert Python arguments to C
   if (!PyArg_ParseTuple(args, "O", &quantity)) {
-      PyErr_SetString(PyExc_RuntimeError, "quantity_ext.c: limit_ minmod_kurganov_ext could not parse input");
+      PyErr_SetString(PyExc_RuntimeError, "quantity_ext.c: limit_minmod_kurganov_ext could not parse input");
       return NULL;
   }
@@ -320 +320 @@
 
 
+//====================================================================
+PyObject *limit_vanalbada_ext(PyObject *self, PyObject *args) {
+
+   PyObject
+           *domain,
+           *quantity;
+
+    PyArrayObject
+        *qco,
+        *qvo,
+        *xco,
+        *xvo;
+
+    double *qc,
+            *qv,
+            *xc,
+            *xv;
+
+  double a, b;
+  //double c;
+  double phi, dx0, dx1;
+  double theta;
+  int N, k, k2;
+
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "O", &quantity)) {
+      PyErr_SetString(PyExc_RuntimeError, "quantity_ext.c: limit_vanalbada_ext could not parse input");
+      return NULL;
+  }
+
+
+  domain = PyObject_GetAttrString(quantity, "domain");
+
+  //printf("B = %p\n",(void*)domain);
+  if (!domain) {
+    printf("quantity_ext.c: Could not obtain python object");
+    fflush(stdout);
+    PyErr_SetString(PyExc_RuntimeError, "quantity_ext.c: Could not obtain python object domain");
+    return NULL;
+  }
+
+
+
+  N  = get_python_integer(quantity,"N");
+  theta = get_python_double(quantity,"beta");
+
+
+  //printf("beta = %f",theta);
+  //fflush(stdout);
+
+  qco = get_consecutive_array(quantity, "centroid_values");
+  qvo = get_consecutive_array(quantity, "vertex_values");
+  xco = get_consecutive_array(domain, "centroids");
+  xvo = get_consecutive_array(domain, "vertices");
+
+  qc = (double *) qco -> data;
+  qv = (double *) qvo -> data;
+  xc = (double *) xco -> data;
+  xv = (double *) xvo -> data;
+
+
+  for (k=0; k<N; k++) {
+      k2 = 2*k;
+      if (k == 0) {
+        phi = (qc[1]-qc[0])/(xc[1] - xc[0]);
+      }
+      else if (k==N-1) {
+          phi = (qc[N-1] - qc[N-2])/(xc[N-1] - xc[N-2]);
+      }
+      else {
+          a  = (qc[k]-qc[k-1])/(xc[k]-xc[k-1]);
+          b  = (qc[k+1]-qc[k])/(xc[k+1]-xc[k]);
+          //c  = (qc[k+1]-qc[k-1])/(xc[k+1]-xc[k-1]);
+
+
+          phi = 0.0;
+          if (a*a + b*b >= 1.0e-32) {
+              phi = (a*a*b+a*b*b)/(a*a+b*b);
+          }
+
+
+      }
+
+      dx0 = xv[k2] - xc[k];
+      dx1 = xv[k2+1] - xc[k];
+
+
+      qv[k2]   = qc[k] + phi*dx0;
+      qv[k2+1] = qc[k] + phi*dx1;
+  }
+
+
+  Py_DECREF(qco);
+  Py_DECREF(qvo);
+  Py_DECREF(xco);
+  Py_DECREF(xvo);
+
+  // Return updated flux timestep
+  return Py_BuildValue("");
+}
+
 
 
@@ -330 +432 @@
   {"limit_minmod_kurganov_ext", limit_minmod_kurganov_ext, METH_VARARGS, "Print out"},
   {"limit_vanleer_ext",         limit_vanleer_ext, METH_VARARGS, "Print out"},
+  {"limit_vanalbada_ext",       limit_vanalbada_ext, METH_VARARGS, "Print out"},
   {NULL, NULL}
 };

trunk/anuga_work/development/2010-projects/anuga_1d/base/test_quantity_1d.py

@@ -595 +595 @@
                                                  [ 10.39393939,  11.60606061],
                                                  [ 11.4,         14.6       ]] )
-#
+
+    def test_limit_vanalbada(self):
+        quantity = Quantity(self.domain4)
+
+        #Create a deliberate overshoot (e.g. from gradient computation)
+        quantity.set_values([[0,0], [2,10], [9,13], [12,14]])
+
+        #Limit
+        quantity.limiter = 'vanalbada'
+        quantity.limit_range()
+
+
+        print quantity.vertex_values
+
+        assert allclose(quantity.vertex_values, [[ -2.4,          2.4       ],
+                                                 [  2.32805995,   9.67194005],
+                                                 [ 10.52104499,  11.47895501],
+                                                 [ 11.4,         14.6       ]])
+
+        from anuga_1d.base.quantity_ext import limit_vanalbada_ext
+
+        quantity.set_values([[0,0], [2,10], [9,13], [12,14]])
+
+        limit_vanalbada_ext(quantity)
+
+        print quantity.vertex_values
+
+
+
+
+        assert allclose(quantity.vertex_values, [[ -2.4,          2.4       ],
+                                                 [  2.32805995,   9.67194005],
+                                                 [ 10.52104499,  11.47895501],
+                                                 [ 11.4,         14.6       ]] )
+
+
+
 
     def test_find_qmax_qmin(self):