Changeset 5742

Timestamp:

Sep 6, 2008, 4:30:16 PM (16 years ago)

Author:

steve

Message:

Found error in get_python_double

Location:

anuga_work/development/anuga_1d

Files:

• dam_h_elevation.py (modified) (6 diffs)
• domain.py (modified) (2 diffs)
• quantity.py (modified) (2 diffs)
• shallow_water_domain_new.py (modified) (3 diffs)
• test_shallow_water.py (modified) (9 diffs)
• util_ext.h (modified) (4 diffs)

anuga_work/development/anuga_1d/dam_h_elevation.py

@@ -1 +1 @@
 import os
 from math import sqrt, sin, cos, pi, exp
-from shallow_water_domain_new import *
+from shallow_water_domain import *
 from Numeric import zeros, Float
 #from analytic_dam_sudi import AnalyticDam
@@ -31 +31 @@
 
 L=2000.0
-N=400
+N=200
 
 cell_len=L/N
@@ -61 +61 @@
     y=zeros(len(x), Float)
     for i in range(len(x)):
-        if 1100.0 <= x[i] <=1200.0:
+        if 800.0 <= x[i] <=1200.0:
             y[i]=10.0
         else:
-            y[i]=10.0
+            y[i]=5.0
     return y
 
@@ -108 +108 @@
 
 domain.set_quantity('stage',stage_perturb)
-domain.set_quantity('elevation', elevation_parabol)
+domain.set_quantity('elevation', elevation_box)
 #domain.set_quantity('xmomentum', xmom_sincos)
 domain.order=domain.default_order
@@ -152 +152 @@
 
         plot(X,ElevationQ,X,HeightQ)
-        plot1.set_ylim([9.99,10.01])
+        #plot1.set_ylim([9.99,10.01])
         xlabel('Position')
         ylabel('Stage')
@@ -180 +180 @@
 
 raw_input("Press any key")
+

anuga_work/development/anuga_1d/domain.py

@@ -57 +57 @@
 
         self.max_speed_array = zeros(N, Float)
+      
 
         self.normals = zeros((N, 2), Float)
@@ -875 +876 @@
                 self.number_of_steps = 0
                 self.number_of_first_order_steps = 0
-                self.max_speed = zeros(N, Float)
+                self.max_speed_array = 0.0
 
 

anuga_work/development/anuga_1d/quantity.py

@@ -860 +860 @@
         # (either from this module or C-extension)
         #backup_centroid_values(self)
-        N = self.domain.number_of_elements
-        for i in range(N):
-            self.centroid_backup_values[i] = self.centroid_values[i]
+
+        self.centroid_backup_values[:] = (self.centroid_values).astype('f')
 
     def saxpy_centroid_values(self,a,b):
@@ -868 +867 @@
         # (either from this module or C-extension)
         #saxpy_centroid_values(self,a,b)
-        N = self.domain.number_of_elements
-        for i in range(N):
-            self.centroid_backup_values[i] = a*self.centroid_values[i] + b*self.centroid_backup_values[i]        
-
+        self.centroid_backup_values[:] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f')
+        
 class Conserved_quantity(Quantity):
     """Class conserved quantity adds to Quantity:

anuga_work/development/anuga_1d/shallow_water_domain_new.py

@@ -49 +49 @@
 #Shallow water domain
 class Domain(Generic_Domain):
-
+    
     def __init__(self, coordinates, boundary = None, tagged_elements = None):
 
@@ -60 +60 @@
         from config import minimum_allowed_height, g
         self.minimum_allowed_height = minimum_allowed_height
+        self.h0 = minimum_allowed_height
         self.g = g
 
@@ -88 +89 @@
         self.quantities_to_be_stored = ['stage', 'xmomentum', 'elevation']
 
-        self.__doc__ = 'shallow_water_domain_new'
-
+        self.__doc__ = 'shallow_water_domain_new'
 
     def set_quantities_to_be_stored(self, q):

anuga_work/development/anuga_1d/test_shallow_water.py

@@ -41 +41 @@
         #print domain.quantities['stage'].explicit_update
         #print domain.quantities['xmomentum'].explicit_update
-        print stage_ud
+        #print stage_ud
 
         assert allclose( domain.quantities['stage'].explicit_update, stage_ud )
@@ -101 +101 @@
         Compare shallow_water_domain gravity calculation
         """
-                
-        def slope_one(x):
+
+        def slope_one(x):
             return x
-                        
+
         domain = Domain(self.points)
         domain.set_quantity('stage',4.0)
-        domain.set_quantity('elevation',slope_one)
+        domain.set_quantity('elevation',slope_one)
         domain.set_boundary({'exterior' : Reflective_boundary(domain)})
 
         gravity(domain)
         
-        #print domain.quantities['stage'].vertex_values 
-        #print domain.quantities['elevation'].vertex_values 
+        #print domain.quantities['stage'].vertex_values 
+        #print domain.quantities['elevation'].vertex_values 
         #print domain.quantities['xmomentum'].explicit_update       
 
@@ -123 +123 @@
         Compare still lake solution for various versions of shallow_water_domain
         """
-                
-        def slope_square(x):
+        
+        def slope_square(x):
             return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
-                        
+
         domain = Domain(self.points2)
         domain.set_quantity('stage',10.0)
-        domain.set_quantity('elevation',slope_square)
+        domain.set_quantity('elevation',slope_square)
         domain.set_boundary({'exterior' : Reflective_boundary(domain)})
 
@@ -141 +141 @@
             pass
         
-##      print domain.quantities['stage'].vertex_values 
-##      print domain.quantities['elevation'].vertex_values 
-##        print domain.quantities['xmomentum'].vertex_values
-##
-##
-##        print domain.quantities['stage'].centroid_values 
-##      print domain.quantities['elevation'].centroid_values 
-##        print domain.quantities['xmomentum'].centroid_values    
+        ##      print domain.quantities['stage'].vertex_values 
+        ##      print domain.quantities['elevation'].vertex_values 
+        ##        print domain.quantities['xmomentum'].vertex_values
+        ##  
+        ##
+        ##        print domain.quantities['stage'].centroid_values  
+        ##      print domain.quantities['elevation'].centroid_values    
+        ##        print domain.quantities['xmomentum'].centroid_values    
 
         assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
@@ -158 +158 @@
         Compare still lake solution for various versions of shallow_water_domain
         """
-                
-        def slope_square(x):
+
+        def slope_square(x):
             return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
-                        
+
         domain = Domain(self.points2)
         domain.set_quantity('stage',10.0)
-        domain.set_quantity('elevation',slope_square)
+        domain.set_quantity('elevation',slope_square)
         domain.set_boundary({'exterior' : Reflective_boundary(domain)})
 
@@ -182 +182 @@
         Compare still lake solution for various versions of shallow_water_domain
         """
-                
-        def slope_square(x):
+
+        def slope_square(x):
             return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
-                        
+
         domain = Domain(self.points2)
         domain.set_quantity('stage',10.0)
-        domain.set_quantity('elevation',slope_square)
+        domain.set_quantity('elevation',slope_square)
         domain.set_boundary({'exterior' : Reflective_boundary(domain)})
 
@@ -238 +238 @@
     Stage = domain.quantities['stage']
     Xmom = domain.quantities['xmomentum']
-#    Ymom = domain.quantities['ymomentum']
+    #    Ymom = domain.quantities['ymomentum']
     Bed = domain.quantities['elevation']
 
@@ -244 +244 @@
     #stage = Stage.edge_values
     #xmom =  Xmom.edge_values
- #   ymom =  Ymom.edge_values
+    #   ymom =  Ymom.edge_values
     #bed =   Bed.edge_values
     
@@ -259 +259 @@
     #print 'stage_bdry',stage_bdry
     #print 'xmom_bdry', xmom_bdry
-#    ymom_bdry =  Ymom.boundary_values
-
-#    flux = zeros(3, Float) #Work array for summing up fluxes
+    #    ymom_bdry =  Ymom.boundary_values
+
+    #    flux = zeros(3, Float) #Work array for summing up fluxes
     flux = zeros(2, Float) #Work array for summing up fluxes
     ql = zeros(2, Float)

anuga_work/development/anuga_1d/util_ext.h

@@ -279 +279 @@
 }
 
+
+
 double get_python_double(PyObject *O, char *name) {
   PyObject *TObject;
+  #define BUFFER_SIZE 80
+  char buf[BUFFER_SIZE];
   double tmp;
+  int n;
   
 
@@ -287 +292 @@
   TObject = PyObject_GetAttrString(O, name);
   if (!TObject) {
-    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_double could not obtain double from object");
+        n =  snprintf(buf, BUFFER_SIZE, "util_ext.h: get_python_double could not obtain double %s.\n", name);
+        //printf("name = %s",name);
+    PyErr_SetString(PyExc_RuntimeError, buf);
+
     return 0.0;
   }  
@@ -298 +306 @@
 }
 
+
+
+
 int get_python_integer(PyObject *O, char *name) {
   PyObject *TObject;
-  int tmp;
+  #define BUFFER_SIZE 80
+  char buf[BUFFER_SIZE];
+  long tmp;
+  int n;
   
 
@@ -306 +320 @@
   TObject = PyObject_GetAttrString(O, name);
   if (!TObject) {
-    PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_integer could not obtain double from object");
+        n =  snprintf(buf, BUFFER_SIZE, "util_ext.h: get_python_integer could not obtain double %s.\n", name);
+        //printf("name = %s",name);
+    PyErr_SetString(PyExc_RuntimeError, buf);
     return 0;
   }  
   
-  tmp = PyFloat_AsDouble(TObject);
+  tmp = PyInt_AsLong(TObject);
   
   Py_DECREF(TObject);