Changeset 4978 for anuga_core/source/anuga/advection

Timestamp:

Jan 28, 2008, 6:14:47 PM (17 years ago)

Author:

steve

Message:

Added a C extension for the advection directory (gave up on using f2py!)

Location:

anuga_core/source/anuga/advection

Files:

• advection.py (modified) (3 diffs)
• advection_ext.c (modified) (6 diffs)
• advection_ext.pyf (deleted)
• compile_all.py (deleted)
• setup.py (deleted)
• test_advection.py (modified) (1 diff)
• test_all.py (added)

anuga_core/source/anuga/advection/advection.py

@@ -66 +66 @@
                                 numproc=numproc)
 
-
+        import Numeric
         if velocity is None:
-            self.velocity = [1,0]
+            self.velocity = Numeric.array([1,0],'d')
         else:
-            self.velocity = velocity
+            self.velocity = Numeric.array(velocity,'d')
 
         #Only first is implemented for advection
@@ -108 +108 @@
         return flux, max_speed
 
+
+
     def compute_fluxes(self):
-
-        
-        self.compute_fluxes_ext()
-##         try:
-##             self.compute_fluxes_ext()
-##         except:
-##             self.compute_fluxes_python()
- 
+        """Compute all fluxes and the timestep suitable for all volumes
+        in domain.
+
+        Compute total flux for each conserved quantity using "flux_function"
+
+        Fluxes across each edge are scaled by edgelengths and summed up
+        Resulting flux is then scaled by area and stored in
+        domain.explicit_update
+
+        The maximal allowable speed computed by the flux_function
+        for each volume
+        is converted to a timestep that must not be exceeded. The minimum of
+        those is computed as the next overall timestep.
+
+        Post conditions:
+        domain.explicit_update is reset to computed flux values
+        domain.timestep is set to the largest step satisfying all volumes.
+        """
+
+        import sys
+        from Numeric import zeros, Float
+        from anuga.config import max_timestep
+
+
+        huge_timestep = float(sys.maxint)
+        Stage = self.quantities['stage']
+
+        """
+        print "======================================"
+        print "BEFORE compute_fluxes"
+        print "stage_update",Stage.explicit_update
+        print "stage_edge",Stage.edge_values
+        print "stage_bdry",Stage.boundary_values
+        print "neighbours",self.neighbours
+        print "neighbour_edges",self.neighbour_edges
+        print "normals",self.normals
+        print "areas",self.areas
+        print "radii",self.radii
+        print "edgelengths",self.edgelengths
+        print "tri_full_flag",self.tri_full_flag
+        print "huge_timestep",huge_timestep
+        print "max_timestep",max_timestep
+        print "velocity",self.velocity
+        """
+
+        import advection_ext            
+        self.timestep = advection_ext.compute_fluxes(self, Stage, huge_timestep, max_timestep)
+
+
+
+    def evolve(self,
+               yieldstep = None,
+               finaltime = None,
+               duration = None,
+               skip_initial_step = False):
+
+        """Specialisation of basic evolve method from parent class
+        """
+
+        #Call basic machinery from parent class
+        for t in Generic_domain.evolve(self,
+                                       yieldstep=yieldstep,
+                                       finaltime=finaltime,
+                                       duration=duration,
+                                       skip_initial_step=skip_initial_step):
+
+            #Pass control on to outer loop for more specific actions
+            yield(t)
+
+
 
 
@@ -207 +271 @@
         self.timestep = timestep
 
-
-    def compute_fluxes_ext(self):
-        """Compute all fluxes and the timestep suitable for all volumes
-        in domain.
-
-        Compute total flux for each conserved quantity using "flux_function"
-
-        Fluxes across each edge are scaled by edgelengths and summed up
-        Resulting flux is then scaled by area and stored in
-        domain.explicit_update
-
-        The maximal allowable speed computed by the flux_function
-        for each volume
-        is converted to a timestep that must not be exceeded. The minimum of
-        those is computed as the next overall timestep.
-
-        Post conditions:
-        domain.explicit_update is reset to computed flux values
-        domain.timestep is set to the largest step satisfying all volumes.
-        """
-
-        import sys
-        from Numeric import zeros, Float
-        from anuga.config import max_timestep
-
-        ntri = len(self)
-
-        timestep = max_timestep
-
-        #Shortcuts
-        Stage = self.quantities['stage']
-
-        #Arrays
-        neighbours      = self.neighbours
-        neighbour_edges = self.neighbour_edges
-        normals         = self.normals
-        areas           = self.areas
-        radii           = self.radii
-        edgelengths     = self.edgelengths
-        tri_full_flag   = self.tri_full_flag
-
-        stage_edge      = Stage.edge_values
-        stage_bdry      = Stage.boundary_values
-        stage_update    = Stage.explicit_update
-
-        huge_timestep = float(sys.maxint)
-
-        v = self.velocity
-
-        nbdry = len(stage_bdry)
-
-        from advection_ext import compute_fluxes
-
-        print "stage_update",stage_update
-        print "stage_edge",stage_edge
-        print "stage_bdry",stage_bdry
-        print "neighbours",neighbours
-        print "neighbour_edges",neighbour_edges
-        print "normals",normals
-        print "areas",areas
-        print "radii",radii
-        print "edgelengths",edgelengths
-        print "tri_full_flag",tri_full_flag
-        print "huge_timestep",huge_timestep
-        print "max_timestep",max_timestep
-        print "v",v
-        print "ntri",ntri
-        print "nbdry",nbdry
-
-
-                
-        timestep = compute_fluxes(stage_update,stage_edge,stage_bdry,
-                                  neighbours,neighbour_edges,normals,
-                                  areas,radii,edgelengths,
-                                  tri_full_flag,
-                                  huge_timestep,max_timestep,v,ntri,nbdry)
-
-        print "stage_update",stage_update        
-        
-        #print 'timestep out2 =',timestep
-
-        self.timestep = timestep
-
-
-    def evolve(self,
-               yieldstep = None,
-               finaltime = None,
-               duration = None,
-               skip_initial_step = False):
-
-        """Specialisation of basic evolve method from parent class
-        """
-
-        #Call basic machinery from parent class
-        for t in Generic_domain.evolve(self,
-                                       yieldstep=yieldstep,
-                                       finaltime=finaltime,
-                                       duration=duration,
-                                       skip_initial_step=skip_initial_step):
-
-            #Pass control on to outer loop for more specific actions
-            yield(t)

anuga_core/source/anuga/advection/advection_ext.c

@@ +1 @@
+// Python - C extension module for advection.py
+//
+// To compile (Python2.X):
+// use python ../utilities/compile.py
+//
+// See the module advection.py
+//
+//
+// Steve Roberts, ANU 2008
+
+
+#include "Python.h"
+#include "Numeric/arrayobject.h"
 #include "math.h"
 #include "stdio.h"
 
-void  print_double_array(char* name, double* array, int n, int m){
-
-    int k,i,km;
-
-    printf("%s = [",name);
-    for (k=0; k<n; k++){
-        km = m*k;
-        printf("[");
-        for (i=0; i<m ; i++){
-            printf("%g ",array[km+i]);
-        }
-        if (k==(n-1))
-            printf("]");
-        else
-            printf("]\n");
-    }
-    printf("]\n");
-}
-
-void  print_int_array(char* name, int* array, int n, int m){
-
-    int k,i,km;
-
-    printf("%s = [",name);
-    for (k=0; k<n; k++){
-        km = m*k;
-        printf("[");
-        for (i=0; i<m ; i++){
-            printf("%i ",array[km+i]);
-        }
-        if (k==(n-1))
-            printf("]");
-        else
-            printf("]\n");
-    }
-    printf("]\n");
-}
-
-
-double compute_fluxes(
-                    double* stage_update,
-                    double* stage_edge,
-                    double* stage_bdry,
-                    int* neighbours,
-                    int* neighbour_edges,
-                    double* normals,
-                    double* areas,
-                    double* radii,
-                    double* edgelengths,
-                    int*   tri_full_flag,
+// Shared code snippets
+#include "util_ext.h"
+
+
+
+double _compute_fluxes(
+                    double* quantity_update,
+                    double* quantity_edge,
+                    double* quantity_bdry,
+                    long*   domain_neighbours,
+                    long*   domain_neighbour_edges,
+                    double* domain_normals,
+                    double* domain_areas,
+                    double* domain_radii,
+                    double* domain_edgelengths,
+                    long*   domain_tri_full_flag,
+                    double* domain_velocity,
                     double  huge_timestep,
                     double  max_timestep,
-                    double* v,
                     int ntri,
                     int nbdry){
+
+ 
         //Loop
 
@@ -72 +53 @@
         timestep = max_timestep;
 
-
-        printf("======================================================\n");
-        printf("INSIDE compute_fluxes\n");
-
-
-        print_double_array( "stage_update",stage_update, ntri, 1);
-        print_double_array( "stage_edge",stage_edge, ntri, 3);
-        print_double_array( "stage_bdry",stage_bdry, nbdry, 1);
-        print_int_array( "neighbours",neighbours, ntri, 3);
-        print_int_array( "neighbour_edges",neighbour_edges, ntri, 3);
-        print_double_array( "normals",normals, ntri, 6);
-        print_double_array( "areas",areas, ntri, 1);
-        print_double_array( "radii",radii, ntri, 1);
-        print_double_array( "edgelengths",edgelengths, ntri, 3);
-        print_int_array( "tri_full_flag",tri_full_flag, ntri, 1);
-        printf("huge_timestep = %g\n",huge_timestep);
-        printf("max_timestep = %g\n",max_timestep);
-        print_double_array( "v",v, 2, 1);
-        printf("ntri = %i\n",ntri);
-        printf("nbdry = %i\n",nbdry);
-
-
         for (k=0; k<ntri; k++){
             optimal_timestep = huge_timestep;
@@ -101 +60 @@
                 k_i = k3+i;
                 //Quantities inside volume facing neighbour i
-                ql = stage_edge[k_i];
+                ql = quantity_edge[k_i];
 
 
                 //Quantities at neighbour on nearest face
-                n = neighbours[k_i];
+                n = domain_neighbours[k_i];
                 if (n < 0) {
                     m = -n-1; //Convert neg flag to index
-                    qr = stage_bdry[m];
+                    qr = quantity_bdry[m];
                 } else {
-                    m = neighbour_edges[k_i];
+                    m = domain_neighbour_edges[k_i];
                     n_m = 3*n+m;
-                    qr = stage_edge[n_m];
+                    qr = quantity_edge[n_m];
                 }
 
@@ -119 +78 @@
                 for (j=0; j<2; j++){
                     k_i_j = 6*k+2*i+j;
-                    normal[j] = normals[k_i_j];
+                    normal[j] = domain_normals[k_i_j];
                 }
 
 
                 //Flux computation using provided function
-                normal_velocity = v[0]*normal[0] + v[1]*normal[1];
+                normal_velocity = domain_velocity[0]*normal[0] + domain_velocity[1]*normal[1];
 
                 if (normal_velocity < 0) {
@@ -133 +92 @@
 
                 max_speed = fabs(normal_velocity);
-                flux = flux - edgeflux * edgelengths[k_i];
+                flux = flux - edgeflux * domain_edgelengths[k_i];
 
                 //Update optimal_timestep
-                if (tri_full_flag[k] == 1) {
+                if (domain_tri_full_flag[k] == 1) {
                     if (max_speed != 0.0) {
-                        optimal_timestep = (optimal_timestep>radii[k]/max_speed) ? radii[k]/max_speed : optimal_timestep;
+                        optimal_timestep = (optimal_timestep>domain_radii[k]/max_speed) ? 
+                          domain_radii[k]/max_speed : optimal_timestep;
                     }
                 }
@@ -146 +106 @@
             //Normalise by area and store for when all conserved
             //quantities get updated
-            stage_update[k] = flux/areas[k];
+            quantity_update[k] = flux/domain_areas[k];
 
             timestep = (timestep>optimal_timestep) ? optimal_timestep : timestep;
 
         }
-
-        //printf("timestep out = %g\n",timestep);
-
-
-
-        printf("INSIDE compute_fluxes, end \n");
-
-        print_double_array( "stage_update",stage_update, ntri, 1);
-
-        printf("FINISHED compute_fluxes\n");
-        printf("======================================================\n");
-
 
         return timestep;
 }
 
+
+//-----------------------------------------------------
+// Python method Wrappers 
+//-----------------------------------------------------
+
+
+
+PyObject *compute_fluxes(PyObject *self, PyObject *args) {
+  /*Compute all fluxes and the timestep suitable for all volumes
+    in domain.
+
+    Fluxes across each edge are scaled by edgelengths and summed up
+    Resulting flux is then scaled by area and stored in
+    explicit_update for the conserved quantity stage.
+
+    The maximal allowable speed computed for each volume
+    is converted to a timestep that must not be exceeded. The minimum of
+    those is computed as the next overall timestep.
+
+    Python call:
+    timestep = advection_ext.compute_fluxes(domain,quantity,huge_timestep,max_timestep)
+
+    Post conditions:
+      domain.explicit_update is reset to computed flux values
+      domain.timestep is set to the largest step satisfying all volumes.
+
+  */
+
+  PyObject *domain, *quantity;
+  PyArrayObject 
+    * quantity_update,
+    * quantity_edge,
+    * quantity_bdry,
+    * domain_neighbours,
+    * domain_neighbour_edges,
+    * domain_normals,
+    * domain_areas,
+    * domain_radii,
+    * domain_edgelengths,
+    * domain_tri_full_flag,
+    * domain_velocity;
+                  
+
+  // Local variables
+  int ntri, nbdry;
+  double huge_timestep, max_timestep;
+  double timestep;
+
+
+  // Convert Python arguments to C
+  if (!PyArg_ParseTuple(args, "OOdd",
+                        &domain,
+                        &quantity,
+                        &huge_timestep,
+                        &max_timestep)) {
+    PyErr_SetString(PyExc_RuntimeError, "advection_ext.c: compute_fluxes could not parse input");
+    return NULL;
+  }
+
+  quantity_edge          = get_consecutive_array(quantity, "edge_values");
+  quantity_bdry          = get_consecutive_array(quantity, "boundary_values");
+  quantity_update        = get_consecutive_array(quantity, "explicit_update");
+  domain_neighbours      = get_consecutive_array(domain,   "neighbours");
+  domain_neighbour_edges = get_consecutive_array(domain,   "neighbour_edges");
+  domain_normals         = get_consecutive_array(domain,   "normals");
+  domain_areas           = get_consecutive_array(domain,   "areas");
+  domain_radii           = get_consecutive_array(domain,   "radii");
+  domain_edgelengths     = get_consecutive_array(domain,   "edgelengths");
+  domain_tri_full_flag   = get_consecutive_array(domain,   "tri_full_flag");
+  domain_velocity        = get_consecutive_array(domain,   "velocity");  
+    
+  ntri  = quantity_edge -> dimensions[0];
+  nbdry = quantity_bdry -> dimensions[0]; 
+
+  timestep = _compute_fluxes((double*) quantity_update -> data,
+                             (double*) quantity_edge -> data,
+                             (double*) quantity_bdry -> data,
+                             (long*)   domain_neighbours -> data,
+                             (long*)   domain_neighbour_edges -> data,
+                             (double*) domain_normals -> data,
+                             (double*) domain_areas ->data,
+                             (double*) domain_radii -> data,
+                             (double*) domain_edgelengths -> data,
+                             (long*)   domain_tri_full_flag -> data,
+                             (double*) domain_velocity -> data,
+                             huge_timestep,
+                             max_timestep,
+                             ntri,
+                             nbdry);
+
+  // Release and return
+  Py_DECREF(quantity_update);
+  Py_DECREF(quantity_edge);
+  Py_DECREF(quantity_bdry);
+  Py_DECREF(domain_neighbours);
+  Py_DECREF(domain_neighbour_edges);
+  Py_DECREF(domain_normals);
+  Py_DECREF(domain_areas);
+  Py_DECREF(domain_radii);
+  Py_DECREF(domain_edgelengths);
+  Py_DECREF(domain_tri_full_flag);
+  Py_DECREF(domain_velocity);
+
+
+  return Py_BuildValue("d", timestep);
+}
+
+
+
+//-------------------------------
+// Method table for python module
+//-------------------------------
+static struct PyMethodDef MethodTable[] = {
+  /* The cast of the function is necessary since PyCFunction values
+   * only take two PyObject* parameters, and rotate() takes
+   * three.
+   */
+
+  {"compute_fluxes", compute_fluxes, METH_VARARGS, "Print out"},
+  {NULL, NULL}
+};
+
+// Module initialisation
+void initadvection_ext(void){
+  Py_InitModule("advection_ext", MethodTable);
+  import_array(); // Necessary for handling of NumPY structures
+}

anuga_core/source/anuga/advection/test_advection.py

@@ -60 +60 @@
         U = -domain.quantities['stage'].explicit_update
         R = -0.5/domain.areas[0]
-
-        print 'U ', U
-        print 'R ', R
-        
 
         assert U==R, '%s %s' %(U, R)