Changeset 5742


Ignore:
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:
6 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/anuga_1d/dam_h_elevation.py

    r5741 r5742  
    11import os
    22from math import sqrt, sin, cos, pi, exp
    3 from shallow_water_domain_new import *
     3from shallow_water_domain import *
    44from Numeric import zeros, Float
    55#from analytic_dam_sudi import AnalyticDam
     
    3131
    3232L=2000.0
    33 N=400
     33N=200
    3434
    3535cell_len=L/N
     
    6161    y=zeros(len(x), Float)
    6262    for i in range(len(x)):
    63         if 1100.0 <= x[i] <=1200.0:
     63        if 800.0 <= x[i] <=1200.0:
    6464            y[i]=10.0
    6565        else:
    66             y[i]=10.0
     66            y[i]=5.0
    6767    return y
    6868
     
    108108
    109109domain.set_quantity('stage',stage_perturb)
    110 domain.set_quantity('elevation', elevation_parabol)
     110domain.set_quantity('elevation', elevation_box)
    111111#domain.set_quantity('xmomentum', xmom_sincos)
    112112domain.order=domain.default_order
     
    152152
    153153        plot(X,ElevationQ,X,HeightQ)
    154         plot1.set_ylim([9.99,10.01])
     154        #plot1.set_ylim([9.99,10.01])
    155155        xlabel('Position')
    156156        ylabel('Stage')
     
    180180
    181181raw_input("Press any key")
     182
  • anuga_work/development/anuga_1d/domain.py

    r5741 r5742  
    5757
    5858        self.max_speed_array = zeros(N, Float)
     59     
    5960
    6061        self.normals = zeros((N, 2), Float)
     
    875876                self.number_of_steps = 0
    876877                self.number_of_first_order_steps = 0
    877                 self.max_speed = zeros(N, Float)
     878                self.max_speed_array = 0.0
    878879
    879880
  • anuga_work/development/anuga_1d/quantity.py

    r5741 r5742  
    860860        # (either from this module or C-extension)
    861861        #backup_centroid_values(self)
    862         N = self.domain.number_of_elements
    863         for i in range(N):
    864             self.centroid_backup_values[i] = self.centroid_values[i]
     862
     863        self.centroid_backup_values[:] = (self.centroid_values).astype('f')
    865864
    866865    def saxpy_centroid_values(self,a,b):
     
    868867        # (either from this module or C-extension)
    869868        #saxpy_centroid_values(self,a,b)
    870         N = self.domain.number_of_elements
    871         for i in range(N):
    872             self.centroid_backup_values[i] = a*self.centroid_values[i] + b*self.centroid_backup_values[i]       
    873 
     869        self.centroid_backup_values[:] = (a*self.centroid_values + b*self.centroid_backup_values).astype('f')
     870       
    874871class Conserved_quantity(Quantity):
    875872    """Class conserved quantity adds to Quantity:
  • anuga_work/development/anuga_1d/shallow_water_domain_new.py

    r5741 r5742  
    4949#Shallow water domain
    5050class Domain(Generic_Domain):
    51 
     51   
    5252    def __init__(self, coordinates, boundary = None, tagged_elements = None):
    5353
     
    6060        from config import minimum_allowed_height, g
    6161        self.minimum_allowed_height = minimum_allowed_height
     62        self.h0 = minimum_allowed_height
    6263        self.g = g
    6364
     
    8889        self.quantities_to_be_stored = ['stage', 'xmomentum', 'elevation']
    8990
    90         self.__doc__ = 'shallow_water_domain_new'
    91 
     91        self.__doc__ = 'shallow_water_domain_new'
    9292
    9393    def set_quantities_to_be_stored(self, q):
  • anuga_work/development/anuga_1d/test_shallow_water.py

    r5741 r5742  
    4141        #print domain.quantities['stage'].explicit_update
    4242        #print domain.quantities['xmomentum'].explicit_update
    43         print stage_ud
     43        #print stage_ud
    4444
    4545        assert allclose( domain.quantities['stage'].explicit_update, stage_ud )
     
    101101        Compare shallow_water_domain gravity calculation
    102102        """
    103                
    104         def slope_one(x):
     103
     104        def slope_one(x):
    105105            return x
    106                        
     106
    107107        domain = Domain(self.points)
    108108        domain.set_quantity('stage',4.0)
    109         domain.set_quantity('elevation',slope_one)
     109        domain.set_quantity('elevation',slope_one)
    110110        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
    111111
    112112        gravity(domain)
    113113       
    114         #print domain.quantities['stage'].vertex_values
    115         #print domain.quantities['elevation'].vertex_values
     114        #print domain.quantities['stage'].vertex_values
     115        #print domain.quantities['elevation'].vertex_values
    116116        #print domain.quantities['xmomentum'].explicit_update       
    117117
     
    123123        Compare still lake solution for various versions of shallow_water_domain
    124124        """
    125                
    126         def slope_square(x):
     125       
     126        def slope_square(x):
    127127            return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
    128                        
     128
    129129        domain = Domain(self.points2)
    130130        domain.set_quantity('stage',10.0)
    131         domain.set_quantity('elevation',slope_square)
     131        domain.set_quantity('elevation',slope_square)
    132132        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
    133133
     
    141141            pass
    142142       
    143 ##      print domain.quantities['stage'].vertex_values
    144 ##      print domain.quantities['elevation'].vertex_values
    145 ##        print domain.quantities['xmomentum'].vertex_values
    146 ##
    147 ##
    148 ##        print domain.quantities['stage'].centroid_values
    149 ##      print domain.quantities['elevation'].centroid_values
    150 ##        print domain.quantities['xmomentum'].centroid_values   
     143        ##      print domain.quantities['stage'].vertex_values
     144        ##      print domain.quantities['elevation'].vertex_values
     145        ##        print domain.quantities['xmomentum'].vertex_values
     146        ## 
     147        ##
     148        ##        print domain.quantities['stage'].centroid_values
     149        ##      print domain.quantities['elevation'].centroid_values   
     150        ##        print domain.quantities['xmomentum'].centroid_values   
    151151
    152152        assert allclose( 10.0*ones(10), domain.quantities['stage'].centroid_values )
     
    158158        Compare still lake solution for various versions of shallow_water_domain
    159159        """
    160                
    161         def slope_square(x):
     160
     161        def slope_square(x):
    162162            return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
    163                        
     163
    164164        domain = Domain(self.points2)
    165165        domain.set_quantity('stage',10.0)
    166         domain.set_quantity('elevation',slope_square)
     166        domain.set_quantity('elevation',slope_square)
    167167        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
    168168
     
    182182        Compare still lake solution for various versions of shallow_water_domain
    183183        """
    184                
    185         def slope_square(x):
     184
     185        def slope_square(x):
    186186            return maximum(4.0-(x-5.0)*(x-5.0), 0.0)
    187                        
     187
    188188        domain = Domain(self.points2)
    189189        domain.set_quantity('stage',10.0)
    190         domain.set_quantity('elevation',slope_square)
     190        domain.set_quantity('elevation',slope_square)
    191191        domain.set_boundary({'exterior' : Reflective_boundary(domain)})
    192192
     
    238238    Stage = domain.quantities['stage']
    239239    Xmom = domain.quantities['xmomentum']
    240 #    Ymom = domain.quantities['ymomentum']
     240    #    Ymom = domain.quantities['ymomentum']
    241241    Bed = domain.quantities['elevation']
    242242
     
    244244    #stage = Stage.edge_values
    245245    #xmom =  Xmom.edge_values
    246  #   ymom =  Ymom.edge_values
     246    #   ymom =  Ymom.edge_values
    247247    #bed =   Bed.edge_values
    248248   
     
    259259    #print 'stage_bdry',stage_bdry
    260260    #print 'xmom_bdry', xmom_bdry
    261 #    ymom_bdry =  Ymom.boundary_values
    262 
    263 #    flux = zeros(3, Float) #Work array for summing up fluxes
     261    #    ymom_bdry =  Ymom.boundary_values
     262
     263    #    flux = zeros(3, Float) #Work array for summing up fluxes
    264264    flux = zeros(2, Float) #Work array for summing up fluxes
    265265    ql = zeros(2, Float)
  • anuga_work/development/anuga_1d/util_ext.h

    r5563 r5742  
    279279}
    280280
     281
     282
    281283double get_python_double(PyObject *O, char *name) {
    282284  PyObject *TObject;
     285  #define BUFFER_SIZE 80
     286  char buf[BUFFER_SIZE];
    283287  double tmp;
     288  int n;
    284289 
    285290
     
    287292  TObject = PyObject_GetAttrString(O, name);
    288293  if (!TObject) {
    289     PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_double could not obtain double from object");
     294        n =  snprintf(buf, BUFFER_SIZE, "util_ext.h: get_python_double could not obtain double %s.\n", name);
     295        //printf("name = %s",name);
     296    PyErr_SetString(PyExc_RuntimeError, buf);
     297
    290298    return 0.0;
    291299  } 
     
    298306}
    299307
     308
     309
     310
    300311int get_python_integer(PyObject *O, char *name) {
    301312  PyObject *TObject;
    302   int tmp;
     313  #define BUFFER_SIZE 80
     314  char buf[BUFFER_SIZE];
     315  long tmp;
     316  int n;
    303317 
    304318
     
    306320  TObject = PyObject_GetAttrString(O, name);
    307321  if (!TObject) {
    308     PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_integer could not obtain double from object");
     322        n =  snprintf(buf, BUFFER_SIZE, "util_ext.h: get_python_integer could not obtain double %s.\n", name);
     323        //printf("name = %s",name);
     324    PyErr_SetString(PyExc_RuntimeError, buf);
    309325    return 0;
    310326  } 
    311327 
    312   tmp = PyFloat_AsDouble(TObject);
     328  tmp = PyInt_AsLong(TObject);
    313329 
    314330  Py_DECREF(TObject);
Note: See TracChangeset for help on using the changeset viewer.