Changeset 8017

Timestamp:

Sep 16, 2010, 10:38:49 AM (14 years ago)

Author:

steve

Message:

Added code to pickup if elevation is discontinuous. compute_fluxes now produces an error in this case

Location:

trunk/anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/quantity.py (modified) (1 diff)
• shallow_water/shallow_water_domain.py (modified) (1 diff)
• shallow_water/shallow_water_ext.c (modified) (22 diffs)
• utilities/polygon_ext.c (modified) (2 diffs)
• utilities/util_ext.h (modified) (1 diff)

trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -1309 +1309 @@
         self.set_vertex_values(A, use_cache=use_cache, verbose=verbose)
 
+
+
     ############################################################################
     # Methods for outputting model results

trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -243 +243 @@
 
 
-    def set_all_betas(self, beta):
-        """Shorthand to assign one constant value [0,1] to all limiters.
+    def set_beta(self, beta):
+        """Shorthand to assign one constant value [0,2] to all limiters.
         0 Corresponds to first order, where as larger values make use of
         the second order scheme.

trunk/anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -22 +22 @@
 // Shared code snippets
 #include "util_ext.h"
+
+
+
 
 
@@ -1045 +1048 @@
             &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;
+      report_python_error(AT, "could not parse input arguments");
+      return NULL;
   }
 
@@ -1090 +1093 @@
             &xmom, &ymom)) {
     //&epsilon)) {
-    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments");
-    return NULL;
+      report_python_error(AT, "could not parse input arguments");
+      return NULL;
   }
 
@@ -1154 +1157 @@
   if (!PyArg_ParseTuple(args, "ddOOOOOOOO",
                         &g, &eps, &x, &w, &uh, &vh, &z,  &eta, &xmom, &ymom)) {
-    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction_new could not parse input arguments");
-    return NULL;
+      report_python_error(AT, "could not parse input arguments");
+      return NULL;
   }
 
@@ -1196 +1199 @@
   if (!PyArg_ParseTuple(args, "ddOOOOOOOO",
                         &g, &eps, &x, &w, &uh, &vh, &z,  &chezy, &xmom, &ymom)) {
-    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: chezy_friction could not parse input arguments");
-    return NULL;
+      report_python_error(AT, "could not parse input arguments");
+      return NULL;
   }
 
@@ -1238 +1241 @@
   if (!PyArg_ParseTuple(args, "ddOOOOOOO",
                         &g, &eps, &w, &uh, &vh, &z,  &eta, &xmom, &ymom)) {
-    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction_old could not parse input arguments");
-    return NULL;
+      report_python_error(AT, "could not parse input arguments");
+      return NULL;
   }
 
@@ -1438 +1441 @@
       if (area2 <= 0) 
       {
-        PyErr_SetString(PyExc_RuntimeError, 
-                        "shallow_water_ext.c: Negative triangle area encountered. This might be due to a triangle being flat or it's vertices not order anti-clockwise.");
-    
-        return -1;
+          report_python_error(AT, "Negative triangle area");
+          return -1;
       }  
       
@@ -1630 +1631 @@
       {
         // If we didn't find an internal neighbour
-        PyErr_SetString(PyExc_RuntimeError, 
-                        "shallow_water_ext.c: Internal neighbour not found");      
+        report_python_error(AT, "Internal neighbour not found");
         return -1;
       }
@@ -1866 +1866 @@
             &elevation_vertex_values,
             &optimise_dry_cells)) {         
-            
-    PyErr_SetString(PyExc_RuntimeError, 
-            "Input arguments to extrapolate_second_order_sw failed");
+
+    report_python_error(AT, "could not parse input arguments");
     return NULL;
   }
@@ -1957 +1956 @@
   if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
                    &Q, &Normal, &direction)) {
-    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: rotate could not parse input arguments");
+    report_python_error(AT, "could not parse input arguments");
     return NULL;
   }  
@@ -1969 +1968 @@
 
   if (normal -> dimensions[0] != 2) {
-    PyErr_SetString(PyExc_RuntimeError, "Normal vector must have 2 components");
+    report_python_error(AT, "Normal vector must have 2 components");
     return NULL;
   }
@@ -2000 +1999 @@
 
 // Computational function for flux computation
+
 double _compute_fluxes_central(int number_of_elements,
-                               double timestep,
-                               double epsilon,
-                               double H0,
-                               double g,
-                               long* neighbours,
-                               long* neighbour_edges,
-                               double* normals,
-                               double* edgelengths, 
-                               double* radii, 
-                               double* areas,
-                               long* tri_full_flag,
-                               double* stage_edge_values,
-                               double* xmom_edge_values,
-                               double* ymom_edge_values,
-                               double* bed_edge_values,
-                               double* stage_boundary_values,
-                               double* xmom_boundary_values,
-                               double* ymom_boundary_values,
-                               double* stage_explicit_update,
-                               double* xmom_explicit_update,
-                               double* ymom_explicit_update,
-                               long* already_computed_flux,
-                               double* max_speed_array,
-                               int optimise_dry_cells) 
-{                   
-  // Local variables
-  double max_speed, length, inv_area, zl, zr;
-  double h0 = H0*H0; // This ensures a good balance when h approaches H0.
-
-  double limiting_threshold = 10*H0; // Avoid applying limiter below this
-                                     // threshold for performance reasons.
-                                     // See ANUGA manual under flux limiting 
-  int k, i, m, n;
-  int ki, nm=0, ki2; // Index shorthands
- 
-  
-  // Workspace (making them static actually made function slightly slower (Ole))  
-  double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
-
-  static long call = 1; // Static local variable flagging already computed flux
- 
-  // Start computation
-  call++; // Flag 'id' of flux calculation for this timestep 
-
-  // Set explicit_update to zero for all conserved_quantities.
-  // This assumes compute_fluxes called before forcing terms
-  memset((char*) stage_explicit_update, 0, number_of_elements*sizeof(double));
-  memset((char*) xmom_explicit_update, 0, number_of_elements*sizeof(double));
-  memset((char*) ymom_explicit_update, 0, number_of_elements*sizeof(double));    
-  
-  // For all triangles
-  for (k = 0; k < number_of_elements; k++) 
-  {  
-    // Loop through neighbours and compute edge flux for each  
-    for (i = 0; i < 3; i++) 
-    {
-      ki = k*3 + i; // Linear index to edge i of triangle k
-      
-      if (already_computed_flux[ki] == call)
-      {
-        // We've already computed the flux across this edge
-        continue;
-      }
-    
-      // Get left hand side values from triangle k, edge i 
-      ql[0] = stage_edge_values[ki];
-      ql[1] = xmom_edge_values[ki];
-      ql[2] = ymom_edge_values[ki];
-      zl = bed_edge_values[ki];
-
-      // Get right hand side values either from neighbouring triangle
-      // or from boundary array (Quantities at neighbour on nearest face).
-      n = neighbours[ki];
-      if (n < 0) 
-      {
-        // Neighbour is a boundary condition
-        m = -n - 1; // Convert negative flag to boundary index
-    
-        qr[0] = stage_boundary_values[m];
-        qr[1] = xmom_boundary_values[m];
-        qr[2] = ymom_boundary_values[m];
-        zr = zl; // Extend bed elevation to boundary
-      } 
-      else 
-      {
-        // Neighbour is a real triangle
-        m = neighbour_edges[ki];
-        nm = n*3 + m; // Linear index (triangle n, edge m)
-        
-        qr[0] = stage_edge_values[nm];
-        qr[1] = xmom_edge_values[nm];
-        qr[2] = ymom_edge_values[nm];
-        zr = bed_edge_values[nm];
-      }
-      
-      // Now we have values for this edge - both from left and right side.
-      
-      if (optimise_dry_cells) 
-      {     
-        // Check if flux calculation is necessary across this edge
-        // This check will exclude dry cells.
-        // This will also optimise cases where zl != zr as 
-        // long as both are dry
-
-        if (fabs(ql[0] - zl) < epsilon && 
-            fabs(qr[0] - zr) < epsilon) 
-        {
-          // Cell boundary is dry
-          
-          already_computed_flux[ki] = call; // #k Done  
-          if (n >= 0)
-          {
-            already_computed_flux[nm] = call; // #n Done
-          }
-        
-          max_speed = 0.0;
-          continue;
-        }
-      }
-      
-      // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
-      ki2 = 2*ki; //k*6 + i*2
-
-      // Edge flux computation (triangle k, edge i)
-      _flux_function_central(ql, qr, zl, zr,
-                             normals[ki2], normals[ki2+1],
-                             epsilon, h0, limiting_threshold, g,
-                             edgeflux, &max_speed);
-      
-      
-      // Multiply edgeflux by edgelength
-      length = edgelengths[ki];
-      edgeflux[0] *= length;            
-      edgeflux[1] *= length;            
-      edgeflux[2] *= length;                        
-      
-      
-      // Update triangle k with flux from edge i
-      stage_explicit_update[k] -= edgeflux[0];
-      xmom_explicit_update[k] -= edgeflux[1];
-      ymom_explicit_update[k] -= edgeflux[2];
-      
-      already_computed_flux[ki] = call; // #k Done
-      
-      
-      // Update neighbour n with same flux but reversed sign
-      if (n >= 0) 
-      {
-        stage_explicit_update[n] += edgeflux[0];
-        xmom_explicit_update[n] += edgeflux[1];
-        ymom_explicit_update[n] += edgeflux[2];
-        
-        already_computed_flux[nm] = call; // #n Done
-      }
-
-      // Update timestep based on edge i and possibly neighbour n
-      if (tri_full_flag[k] == 1) 
-      {
-        if (max_speed > epsilon) 
-        {
-          // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
-          
-          // CFL for triangle k
-          timestep = min(timestep, radii[k]/max_speed);
-          
-          if (n >= 0)
-          {
-            // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
-            timestep = min(timestep, radii[n]/max_speed);
-          }
-
-          // Ted Rigby's suggested less conservative version
-          //if (n>=0) {         
-          //  timestep = min(timestep, (radii[k]+radii[n])/max_speed);
-          //} else {
-          //  timestep = min(timestep, radii[k]/max_speed);
-          // }
-        }
-      }
-      
-    } // End edge i (and neighbour n)
-    
-    
-    // Normalise triangle k by area and store for when all conserved
-    // quantities get updated
-    inv_area = 1.0/areas[k];
-    stage_explicit_update[k] *= inv_area;
-    xmom_explicit_update[k] *= inv_area;
-    ymom_explicit_update[k] *= inv_area;
-    
-   
-    // Keep track of maximal speeds
-    max_speed_array[k] = max_speed;
-    
-  } // End triangle k
-             
-  return timestep;
+        double timestep,
+        double epsilon,
+        double H0,
+        double g,
+        long* neighbours,
+        long* neighbour_edges,
+        double* normals,
+        double* edgelengths,
+        double* radii,
+        double* areas,
+        long* tri_full_flag,
+        double* stage_edge_values,
+        double* xmom_edge_values,
+        double* ymom_edge_values,
+        double* bed_edge_values,
+        double* stage_boundary_values,
+        double* xmom_boundary_values,
+        double* ymom_boundary_values,
+        double* stage_explicit_update,
+        double* xmom_explicit_update,
+        double* ymom_explicit_update,
+        long* already_computed_flux,
+        double* max_speed_array,
+        int optimise_dry_cells) {
+    // Local variables
+    double max_speed, length, inv_area, zl, zr;
+    double h0 = H0*H0; // This ensures a good balance when h approaches H0.
+
+    double limiting_threshold = 10 * H0; // Avoid applying limiter below this
+    // threshold for performance reasons.
+    // See ANUGA manual under flux limiting
+    int k, i, m, n;
+    int ki, nm = 0, ki2; // Index shorthands
+
+
+    // Workspace (making them static actually made function slightly slower (Ole))
+    double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
+
+    static long call = 1; // Static local variable flagging already computed flux
+
+    // Start computation
+    call++; // Flag 'id' of flux calculation for this timestep
+
+    // Set explicit_update to zero for all conserved_quantities.
+    // This assumes compute_fluxes called before forcing terms
+    memset((char*) stage_explicit_update, 0, number_of_elements * sizeof (double));
+    memset((char*) xmom_explicit_update, 0, number_of_elements * sizeof (double));
+    memset((char*) ymom_explicit_update, 0, number_of_elements * sizeof (double));
+
+    // For all triangles
+    for (k = 0; k < number_of_elements; k++) {
+        // Loop through neighbours and compute edge flux for each
+        for (i = 0; i < 3; i++) {
+            ki = k * 3 + i; // Linear index to edge i of triangle k
+
+            if (already_computed_flux[ki] == call) {
+                // We've already computed the flux across this edge
+                continue;
+            }
+
+            // Get left hand side values from triangle k, edge i
+            ql[0] = stage_edge_values[ki];
+            ql[1] = xmom_edge_values[ki];
+            ql[2] = ymom_edge_values[ki];
+            zl = bed_edge_values[ki];
+
+            // Get right hand side values either from neighbouring triangle
+            // or from boundary array (Quantities at neighbour on nearest face).
+            n = neighbours[ki];
+            if (n < 0) {
+                // Neighbour is a boundary condition
+                m = -n - 1; // Convert negative flag to boundary index
+
+                qr[0] = stage_boundary_values[m];
+                qr[1] = xmom_boundary_values[m];
+                qr[2] = ymom_boundary_values[m];
+                zr = zl; // Extend bed elevation to boundary
+            }
+            else {
+                // Neighbour is a real triangle
+                m = neighbour_edges[ki];
+                nm = n * 3 + m; // Linear index (triangle n, edge m)
+
+                qr[0] = stage_edge_values[nm];
+                qr[1] = xmom_edge_values[nm];
+                qr[2] = ymom_edge_values[nm];
+                zr = bed_edge_values[nm];
+            }
+
+            // Now we have values for this edge - both from left and right side.
+
+            if (optimise_dry_cells) {
+                // Check if flux calculation is necessary across this edge
+                // This check will exclude dry cells.
+                // This will also optimise cases where zl != zr as
+                // long as both are dry
+
+                if (fabs(ql[0] - zl) < epsilon &&
+                        fabs(qr[0] - zr) < epsilon) {
+                    // Cell boundary is dry
+
+                    already_computed_flux[ki] = call; // #k Done
+                    if (n >= 0) {
+                        already_computed_flux[nm] = call; // #n Done
+                    }
+
+                    max_speed = 0.0;
+                    continue;
+                }
+            }
+
+
+            if (fabs(zl-zr)>1.0e-10) {
+                report_python_error(AT,"Discontinuous Elevation");
+                return 0.0;
+            }
+            
+            // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
+            ki2 = 2 * ki; //k*6 + i*2
+
+            // Edge flux computation (triangle k, edge i)
+            _flux_function_central(ql, qr, zl, zr,
+                    normals[ki2], normals[ki2 + 1],
+                    epsilon, h0, limiting_threshold, g,
+                    edgeflux, &max_speed);
+
+
+            // Multiply edgeflux by edgelength
+            length = edgelengths[ki];
+            edgeflux[0] *= length;
+            edgeflux[1] *= length;
+            edgeflux[2] *= length;
+
+
+            // Update triangle k with flux from edge i
+            stage_explicit_update[k] -= edgeflux[0];
+            xmom_explicit_update[k] -= edgeflux[1];
+            ymom_explicit_update[k] -= edgeflux[2];
+
+            already_computed_flux[ki] = call; // #k Done
+
+
+            // Update neighbour n with same flux but reversed sign
+            if (n >= 0) {
+                stage_explicit_update[n] += edgeflux[0];
+                xmom_explicit_update[n] += edgeflux[1];
+                ymom_explicit_update[n] += edgeflux[2];
+
+                already_computed_flux[nm] = call; // #n Done
+            }
+
+            // Update timestep based on edge i and possibly neighbour n
+            if (tri_full_flag[k] == 1) {
+                if (max_speed > epsilon) {
+                    // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
+
+                    // CFL for triangle k
+                    timestep = min(timestep, radii[k] / max_speed);
+
+                    if (n >= 0) {
+                        // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
+                        timestep = min(timestep, radii[n] / max_speed);
+                    }
+
+                    // Ted Rigby's suggested less conservative version
+                    //if (n>=0) {
+                    //  timestep = min(timestep, (radii[k]+radii[n])/max_speed);
+                    //} else {
+                    //  timestep = min(timestep, radii[k]/max_speed);
+                    // }
+                }
+            }
+
+        } // End edge i (and neighbour n)
+
+
+        // Normalise triangle k by area and store for when all conserved
+        // quantities get updated
+        inv_area = 1.0 / areas[k];
+        stage_explicit_update[k] *= inv_area;
+        xmom_explicit_update[k] *= inv_area;
+        ymom_explicit_update[k] *= inv_area;
+
+
+        // Keep track of maximal speeds
+        max_speed_array[k] = max_speed;
+
+    } // End triangle k
+
+    return timestep;
 }
 
@@ -2264 +2257 @@
     // Convert Python arguments to C
     if (!PyArg_ParseTuple(args, "dOOOO", &timestep, &domain, &stage, &xmom, &ymom, &bed )) {
-    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
-    return NULL;
+      report_python_error(AT, "could not parse input arguments");
+      return NULL;
     }
 
@@ -2450 +2443 @@
             &max_speed_array,
             &optimise_dry_cells)) {
-    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
+    report_python_error(AT, "could not parse input arguments");
     return NULL;
   }
@@ -2610 +2603 @@
             &ymom_explicit_update,
             &already_computed_flux)) {
-    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
+    report_python_error(AT, "could not parse input arguments");
     return NULL;
   }
@@ -2716 +2709 @@
             &epsilon,
             &wc, &zc, &xmomc, &ymomc)) {
-    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: protect could not parse input arguments");
+    report_python_error(AT, "could not parse input arguments");
     return NULL;
   }  
@@ -2775 +2768 @@
             &wv, &zv, &hvbar,
             &xmomc, &ymomc, &xmomv, &ymomv)) {
-    PyErr_SetString(PyExc_RuntimeError, 
-            "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
+    report_python_error(AT, "could not parse input arguments");
     return NULL;
   }  
@@ -2788 +2780 @@
   Tmp = PyObject_GetAttrString(domain, "alpha_balance");
   if (!Tmp) {
-    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object alpha_balance from domain");
+    report_python_error(AT, "could not obtain object alpha_balance from domain");
     return NULL;
   }  
@@ -2797 +2789 @@
   Tmp = PyObject_GetAttrString(domain, "H0");
   if (!Tmp) {
-    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object H0 from domain");
+    report_python_error(AT, "could not obtain object H0 from domain");
     return NULL;
   }  
@@ -2806 +2798 @@
   Tmp = PyObject_GetAttrString(domain, "tight_slope_limiters");
   if (!Tmp) {
-    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object tight_slope_limiters from domain");
+    report_python_error(AT, "could not obtain object tight_slope_limiters from domain");
     return NULL;
   }  
@@ -2815 +2807 @@
   Tmp = PyObject_GetAttrString(domain, "use_centroid_velocities");
   if (!Tmp) {
-    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: balance_deep_and_shallow could not obtain object use_centroid_velocities from domain");
+    report_python_error(AT, "could not obtain object use_centroid_velocities from domain");
     return NULL;
   }  
@@ -2869 +2861 @@
             &s_vec, &phi_vec,
             &cw)) {
-    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: assign_windfield_values could not parse input arguments");
+    report_python_error(AT, "could not parse input arguments");
     return NULL;
   }  

trunk/anuga_core/source/anuga/utilities/polygon_ext.c

@@ -20 +20 @@
 #include "util_ext.h"
 
+#define YES 1
+#define NO 0
+
 
 double dist(double x,
@@ -89 +92 @@
 
 
+//  public domain function by Darel Rex Finley, 2006
+//  http://www.alienryderflex.com/intersect/
+//
+//  Determines the intersection point of the line segment defined by points A and B
+//  with the line segment defined by points C and D.
+//
+//  Returns YES if the intersection point was found, and stores that point in X,Y.
+//  Returns NO if there is no determinable intersection point, in which case X,Y will
+//  be unmodified.
+
+int __lineSegmentIntersection(
+        double Ax, double Ay,
+        double Bx, double By,
+        double Cx, double Cy,
+        double Dx, double Dy,
+        double *X, double *Y) {
+
+    double distAB, theCos, theSin, newX, ABpos;
+
+    //  Fail if either line segment is zero-length.
+    if ( (Ax == Bx && Ay == By) || (Cx == Dx && Cy == Dy) ) return NO ;
+
+    //  Fail if the segments share an end-point.
+    if ( (Ax == Cx && Ay == Cy) || (Bx == Cx && By == Cy)
+            || (Ax == Dx && Ay == Dy) || (Bx == Dx && By == Dy) ) {
+        return NO;
+    }
+
+    //  (1) Translate the system so that point A is on the origin.
+    Bx -= Ax;
+    By -= Ay;
+    Cx -= Ax;
+    Cy -= Ay;
+    Dx -= Ax;
+    Dy -= Ay;
+
+    //  Discover the length of segment A-B.
+    distAB = sqrt(Bx * Bx + By * By);
+
+    //  (2) Rotate the system so that point B is on the positive X axis.
+    theCos = Bx / distAB;
+    theSin = By / distAB;
+    newX = Cx * theCos + Cy*theSin;
+    Cy = Cy * theCos - Cx*theSin;
+    Cx = newX;
+    newX = Dx * theCos + Dy*theSin;
+    Dy = Dy * theCos - Dx*theSin;
+    Dx = newX;
+
+    //  Fail if segment C-D doesn't cross line A-B.
+    if ( (Cy < 0. && Dy < 0.) || (Cy >= 0. && Dy >= 0.) ) return NO;
+
+    //  (3) Discover the position of the intersection point along line A-B.
+    ABpos = Dx + (Cx - Dx) * Dy / (Dy - Cy);
+
+    //  Fail if segment C-D crosses line A-B outside of segment A-B.
+    if (ABpos < 0. || ABpos > distAB) return NO;
+
+    //  (4) Apply the discovered position to line A-B in the original coordinate system.
+    *X = Ax + ABpos*theCos;
+    *Y = Ay + ABpos*theSin;
+
+    //  Success.
+    return YES;
+}
 
 /*

trunk/anuga_core/source/anuga/utilities/util_ext.h

@@ -14 +14 @@
 #include "math.h"
 
+#include <stdio.h>
+#define STRINGIFY(x) #x
+#define TOSTRING(x) STRINGIFY(x)
+#define AT __FILE__ ":" TOSTRING(__LINE__)
+#define P_ERROR_BUFFER_SIZE 65
+void report_python_error(const char *location, const char *msg)
+{
+
+    char buf[P_ERROR_BUFFER_SIZE];
+    int n;
+    n = snprintf(buf, P_ERROR_BUFFER_SIZE, "Error at %s: %s\n", location, msg);
+
+    PyErr_SetString(PyExc_RuntimeError, buf);
+}
+
+
 
 double max(double x, double y) {