Changeset 6902


Ignore:
Timestamp:
Apr 24, 2009, 5:22:14 PM (15 years ago)
Author:
rwilson
Message:

Back-merge from Numeric trunk to numpy branch.

Files:
20 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/utilities/log_test.py

    r6587 r6902  
    3131        if os.path.exists(STDOUT_LOG_NAME):
    3232            os.remove(STDOUT_LOG_NAME)
    33         pass
    3433
    3534    ##
  • branches/numpy/anuga/abstract_2d_finite_volumes/domain.py

    r6553 r6902  
    9292        """
    9393
     94        number_of_full_nodes=None
     95        number_of_full_triangles=None
     96       
    9497        # Determine whether source is a mesh filename or coordinates
    9598        if type(source) == types.StringType:
     
    117120
    118121        # Expose Mesh attributes (FIXME: Maybe turn into methods)
     122        self.triangles = self.mesh.triangles       
    119123        self.centroid_coordinates = self.mesh.centroid_coordinates
    120124        self.vertex_coordinates = self.mesh.vertex_coordinates
     
    205209        # the ghost triangles.
    206210        if not num.allclose(self.tri_full_flag[:self.number_of_full_nodes], 1):
    207             print ('WARNING:  '
    208                    'Not all full triangles are store before ghost triangles')
     211            if self.numproc>1:
     212                print ('WARNING:  '
     213                       'Not all full triangles are store before ghost triangles')
    209214
    210215        # Defaults
     
    306311    def get_normal(self, *args, **kwargs):
    307312        return self.mesh.get_normal(*args, **kwargs)
     313    def get_triangle_containing_point(self, *args, **kwargs):
     314        return self.mesh.get_triangle_containing_point(*args, **kwargs)
    308315
    309316    def get_intersecting_segments(self, *args, **kwargs):
     
    17451752
    17461753    ##
    1747     # @brief ??
     1754    # @brief Sequential update of ghost cells
    17481755    def update_ghosts(self):
    1749         pass
    1750 
     1756        # We must send the information from the full cells and
     1757        # receive the information for the ghost cells
     1758        # We have a list with ghosts expecting updates
     1759
     1760        from Numeric import take,put
     1761
     1762
     1763        #Update of ghost cells
     1764        iproc = self.processor
     1765        if self.full_send_dict.has_key(iproc):
     1766
     1767            # now store full as local id, global id, value
     1768            Idf  = self.full_send_dict[iproc][0]
     1769
     1770            # now store ghost as local id, global id, value
     1771            Idg = self.ghost_recv_dict[iproc][0]
     1772
     1773            for i, q in enumerate(self.conserved_quantities):
     1774                Q_cv =  self.quantities[q].centroid_values
     1775                put(Q_cv,     Idg, take(Q_cv,     Idf))
     1776
     1777 
    17511778    ##
    17521779    # @brief Extrapolate conserved quantities from centroid to vertices
  • branches/numpy/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

    r6553 r6902  
    8282    """Time dependent boundary returns values for the
    8383    conserved quantities as a function of time.
    84     Must specify domain to get access to model time and a function f(t)
    85     which must return conserved quantities as a function time
     84    Must specify domain to get access to model time and a function of t
     85    which must return conserved quantities as a function time.
     86   
     87    Example:
     88      B = Time_boundary(domain,
     89                        function=lambda t: [(60<t<3660)*2, 0, 0])
     90     
     91      This will produce a boundary condition with is a 2m high square wave
     92      starting 60 seconds into the simulation and lasting one hour.
     93      Momentum applied will be 0 at all times.
     94                       
    8695    """
    8796
     
    98107        self.verbose = verbose
    99108
     109       
     110        # FIXME: Temporary code to deal with both f and function
     111        if function is not None and f is not None:
     112            raise Exception, 'Specify either function or f to Time_boundary'
     113           
     114        if function is None and f is None:
     115            raise Exception, 'You must specify a function to Time_boundary'
     116           
     117        if f is None:
     118            f = function
     119        #####
     120       
    100121        try:
    101122            q = f(0.0)
  • branches/numpy/anuga/abstract_2d_finite_volumes/mesh_factory.py

    r6304 r6902  
    201201
    202202    return points, elements, boundary
     203
     204
     205def rectangular_periodic(m_g, n_g, len1_g=1.0, len2_g=1.0, origin_g = (0.0, 0.0)):
     206
     207
     208    """Setup a rectangular grid of triangles
     209    with m+1 by n+1 grid points
     210    and side lengths len1, len2. If side lengths are omitted
     211    the mesh defaults to the unit square, divided between all the
     212    processors
     213
     214    len1: x direction (left to right)
     215    len2: y direction (bottom to top)
     216
     217    """
     218
     219    processor = 0
     220    numproc   = 1
     221
     222   
     223    n = n_g
     224    m_low  = -1
     225    m_high = m_g +1
     226
     227    m = m_high - m_low
     228
     229    delta1 = float(len1_g)/m_g
     230    delta2 = float(len2_g)/n_g
     231
     232    len1 = len1_g*float(m)/float(m_g)
     233    len2 = len2_g
     234    origin = ( origin_g[0]+float(m_low)/float(m_g)*len1_g, origin_g[1] )
     235
     236    #Calculate number of points
     237    Np = (m+1)*(n+1)
     238
     239    class VIndex:
     240
     241        def __init__(self, n,m):
     242            self.n = n
     243            self.m = m
     244
     245        def __call__(self, i,j):
     246            return j+i*(self.n+1)
     247
     248    class EIndex:
     249
     250        def __init__(self, n,m):
     251            self.n = n
     252            self.m = m
     253
     254        def __call__(self, i,j):
     255            return 2*(j+i*self.n)
     256
     257
     258    I = VIndex(n,m)
     259    E = EIndex(n,m)
     260
     261    points = num.zeros( (Np,2), num.Float)
     262
     263    for i in range(m+1):
     264        for j in range(n+1):
     265
     266            points[I(i,j),:] = [i*delta1 + origin[0], j*delta2 + origin[1]]
     267
     268    #Construct 2 triangles per rectangular element and assign tags to boundary
     269    #Calculate number of triangles
     270    Nt = 2*m*n
     271
     272
     273    elements = num.zeros( (Nt,3), num.Int)
     274    boundary = {}
     275    Idgl = []
     276    Idfl = []
     277    Idgr = []
     278    Idfr = []
     279
     280    full_send_dict = {}
     281    ghost_recv_dict = {}
     282    nt = -1
     283    for i in range(m):
     284        for j in range(n):
     285
     286            i1 = I(i,j+1)
     287            i2 = I(i,j)
     288            i3 = I(i+1,j+1)
     289            i4 = I(i+1,j)
     290
     291            #Lower Element
     292            nt = E(i,j)
     293            if i == 0:
     294                Idgl.append(nt)
     295
     296            if i == 1:
     297                Idfl.append(nt)
     298
     299            if i == m-2:
     300                Idfr.append(nt)
     301
     302            if i == m-1:
     303                Idgr.append(nt)
     304
     305            if i == m-1:
     306                if processor == numproc-1:
     307                    boundary[nt, 2] = 'right'
     308                else:
     309                    boundary[nt, 2] = 'ghost'
     310       
     311            if j == 0:
     312                boundary[nt, 1] = 'bottom'
     313            elements[nt,:] = [i4,i3,i2]
     314
     315            #Upper Element
     316            nt = E(i,j)+1
     317            if i == 0:
     318                Idgl.append(nt)
     319
     320            if i == 1:
     321                Idfl.append(nt)
     322
     323            if i == m-2:
     324                Idfr.append(nt)
     325
     326            if i == m-1:
     327                Idgr.append(nt)
     328
     329            if i == 0:
     330                if processor == 0:
     331                    boundary[nt, 2] = 'left'
     332                else:
     333                    boundary[nt, 2] = 'ghost'
     334            if j == n-1:
     335                boundary[nt, 1] = 'top'
     336            elements[nt,:] = [i1,i2,i3]
     337
     338    Idfl.extend(Idfr)
     339    Idgr.extend(Idgl)
     340
     341    Idfl = num.array(Idfl, num.Int)
     342    Idgr = num.array(Idgr, num.Int)
     343   
     344    full_send_dict[processor]  = [Idfl, Idfl]
     345    ghost_recv_dict[processor] = [Idgr, Idgr]
     346
     347
     348    return  points, elements, boundary, full_send_dict, ghost_recv_dict
     349
    203350
    204351def oblique(m, n, lenx = 1.0, leny = 1.0, theta = 8.95, origin = (0.0, 0.0)):
  • branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py

    r6689 r6902  
    11781178        # Edges have already been deprecated in set_values, see changeset:5521,
    11791179        # but *might* be useful in get_values. Any thoughts anyone?
     1180        # YES (Ole): Edge values are necessary for volumetric balance check
    11801181
    11811182        if location not in ['vertices', 'centroids',
  • branches/numpy/anuga/abstract_2d_finite_volumes/test_domain.py

    r6553 r6902  
    397397                                     
    398398    def test_setting_timestepping_method(self):
    399         """test_set_quanitities_to_be_monitored
     399        """test_setting_timestepping_method
    400400        """
    401401
     
    803803                             [ 11.0,  11.0,  11.0]])
    804804                             
     805    def test_rectangular_periodic_and_ghosts(self):
     806
     807        from mesh_factory import rectangular_periodic
     808       
     809
     810        M=5
     811        N=2
     812        points, vertices, boundary, full_send_dict, ghost_recv_dict = rectangular_periodic(M, N)
     813
     814        assert num.allclose(ghost_recv_dict[0][0], [24, 25, 26, 27,  0,  1,  2,  3])
     815        assert num.allclose(full_send_dict[0][0] , [ 4,  5,  6,  7, 20, 21, 22, 23])
     816
     817        #print 'HERE'
     818       
     819        conserved_quantities = ['quant1', 'quant2']
     820        domain = Domain(points, vertices, boundary, conserved_quantities,
     821                        full_send_dict=full_send_dict,
     822                        ghost_recv_dict=ghost_recv_dict)
     823
     824
     825
     826
     827        assert num.allclose(ghost_recv_dict[0][0], [24, 25, 26, 27,  0,  1,  2,  3])
     828        assert num.allclose(full_send_dict[0][0] , [ 4,  5,  6,  7, 20, 21, 22, 23])
     829
     830        def xylocation(x,y):
     831            return 15*x + 9*y
     832
     833       
     834        domain.set_quantity('quant1',xylocation,location='centroids')
     835        domain.set_quantity('quant2',xylocation,location='centroids')
     836
     837
     838        assert num.allclose(domain.quantities['quant1'].centroid_values,
     839                            [  0.5,   1.,   5.,    5.5,   3.5,   4.,    8.,    8.5,   6.5,  7.,   11.,   11.5,   9.5,
     840                               10.,   14.,   14.5,  12.5,  13.,   17.,   17.5,  15.5,  16.,   20.,   20.5,
     841                               18.5,  19.,   23.,   23.5])
     842
     843
     844
     845        assert num.allclose(domain.quantities['quant2'].centroid_values,
     846                            [  0.5,   1.,   5.,    5.5,   3.5,   4.,    8.,    8.5,   6.5,  7.,   11.,   11.5,   9.5,
     847                               10.,   14.,   14.5,  12.5,  13.,   17.,   17.5,  15.5,  16.,   20.,   20.5,
     848                               18.5,  19.,   23.,   23.5])
     849
     850        domain.update_ghosts()
     851
     852
     853        assert num.allclose(domain.quantities['quant1'].centroid_values,
     854                            [  15.5,  16.,   20.,   20.5,   3.5,   4.,    8.,    8.5,   6.5,  7.,   11.,   11.5,   9.5,
     855                               10.,   14.,   14.5,  12.5,  13.,   17.,   17.5,  15.5,  16.,   20.,   20.5,
     856                                3.5,   4.,    8.,    8.5])
     857
     858
     859
     860        assert num.allclose(domain.quantities['quant2'].centroid_values,
     861                            [  15.5,  16.,   20.,   20.5,   3.5,   4.,    8.,    8.5,   6.5,  7.,   11.,   11.5,   9.5,
     862                               10.,   14.,   14.5,  12.5,  13.,   17.,   17.5,  15.5,  16.,   20.,   20.5,
     863                                3.5,   4.,    8.,    8.5])
     864
     865       
     866        assert num.allclose(domain.tri_full_flag, [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     867                                                   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0])
     868
     869        #Test that points are arranged in a counter clock wise order
     870        domain.check_integrity()
     871
     872
    805873    def test_that_mesh_methods_exist(self):
    806874        """test_that_mesh_methods_exist
     
    828896        domain.get_number_of_nodes()
    829897        domain.get_normal(0,0)
     898        domain.get_triangle_containing_point([0.4,0.5])
    830899        domain.get_intersecting_segments([[0.0, 0.0], [0.0, 1.0]])
    831900        domain.get_disconnected_triangles()
  • branches/numpy/anuga/coordinate_transforms/redfearn.py

    r6553 r6902  
    4848    If zone is specified reproject lat and long to specified zone instead of
    4949    standard zone
     50
    5051    If meridian is specified, reproject lat and lon to that instead of zone. In this case
    5152    zone will be set to -1 to indicate non-UTM projection
  • branches/numpy/anuga/culvert_flows/culvert_class.py

    r6553 r6902  
    325325            if time - self.last_update > self.update_interval or time == 0.0:
    326326                update = True
    327 
    328327           
    329328        if self.log_filename is not None:       
  • branches/numpy/anuga/culvert_flows/culvert_routines.py

    r6689 r6902  
    217217    else: # inlet_depth < 0.01:
    218218        Q = barrel_velocity = outlet_culvert_depth = 0.0
     219
    219220    # Temporary flow limit
    220221    if barrel_velocity > max_velocity:
  • branches/numpy/anuga/culvert_flows/test_culvert_routines.py

    r6790 r6902  
    1919
    2020
    21     def test_boyd_1(self):
    22         """test_boyd_1
    23        
    24         This tests the Boyd routine with data obtained from ??? by Petar Milevski   
    25         """
    26         # FIXME(Ole): This test fails (20 Feb 2009)
    27 
     21    def test_boyd_0(self):
     22        """test_boyd_0
     23       
     24        This tests the Boyd routine with data obtained from ??? by Petar Milevski
     25        This test is the only one that passed in late February 2009
     26        """
     27     
    2828        g=9.81
    2929        culvert_slope=0.1  # Downward
     
    6767       
    6868
    69     def test_boyd_2(self):
    70         """test_boyd_2
     69    def Xtest_boyd_00(self):
     70        """test_boyd_00
    7171       
    7272        This tests the Boyd routine with data obtained from ??? by Petar Milevski   
     
    108108       
    109109        #print Q, v, d
    110         #assert num.allclose(Q, 0.185, rtol=1.0e-3)
     110        assert num.allclose(Q, 0.185, rtol=1.0e-3)
    111111        #assert num.allclose(v, 0.93)
    112112        #assert num.allclose(d, 0.0)
    113113       
     114    def Xtest_boyd_1(self):
     115        """test_boyd_1
     116       
     117        This tests the Boyd routine with data obtained from ??? by Petar Milevski   
     118        """
     119        # FIXME(Ole): This test fails (20 Feb 2009)
     120
     121        g=9.81
     122        culvert_slope=0.01  # Downward
     123
     124        inlet_depth=0.263
     125        outlet_depth=0.0
     126
     127        culvert_length=4.0
     128        culvert_width=0.75
     129        culvert_height=0.75
     130       
     131        culvert_type='pipe'
     132        manning=0.013
     133        sum_loss=1.5
     134
     135        inlet_specific_energy=inlet_depth #+0.5*v**2/g
     136        z_in = 0.0
     137        z_out = -culvert_length*culvert_slope/100
     138        E_in = z_in+inlet_depth  #+ 0.5*v**2/g
     139        E_out = z_out+outlet_depth  #+ 0.5*v**2/g
     140        delta_total_energy = E_in-E_out
     141
     142        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
     143                                                 outlet_depth,
     144                                                 inlet_specific_energy,
     145                                                 delta_total_energy,
     146                                                 g,
     147                                                 culvert_length,
     148                                                 culvert_width,
     149                                                 culvert_height,
     150                                                 culvert_type,
     151                                                 manning,
     152                                                 sum_loss)
     153       
     154        print Q, v, d
     155        assert num.allclose(Q, 0.10, rtol=1.0e-2) #inflow
     156        assert num.allclose(v, 1.13, rtol=1.0e-2) #outflow velocity
     157        assert num.allclose(d, 0.15, rtol=1.0e-2) #depth at outlet used to calc v
     158       
     159    def Xtest_boyd_2(self):
     160        """test_boyd_2
     161       
     162        This tests the Boyd routine with data obtained from ??? by Petar Milevski   
     163        """
     164        # FIXME(Ole): This test fails (20 Feb 2009)
     165
     166        g=9.81
     167        culvert_slope=0.01  # Downward
     168
     169        inlet_depth=1.135
     170        outlet_depth=0.0
     171
     172        culvert_length=4.0
     173        culvert_width=0.75
     174        culvert_height=0.75
     175       
     176        culvert_type='pipe'
     177        manning=0.013
     178        sum_loss=1.5
     179
     180        inlet_specific_energy=inlet_depth #+0.5*v**2/g
     181        z_in = 0.0
     182        z_out = -culvert_length*culvert_slope/100
     183        E_in = z_in+inlet_depth  #+ 0.5*v**2/g
     184        E_out = z_out+outlet_depth  #+ 0.5*v**2/g
     185        delta_total_energy = E_in-E_out
     186
     187        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
     188                                                 outlet_depth,
     189                                                 inlet_specific_energy,
     190                                                 delta_total_energy,
     191                                                 g,
     192                                                 culvert_length,
     193                                                 culvert_width,
     194                                                 culvert_height,
     195                                                 culvert_type,
     196                                                 manning,
     197                                                 sum_loss)
     198       
     199        print Q, v, d
     200        assert num.allclose(Q, 1.00, rtol=1.0e-2) #inflow
     201        assert num.allclose(v, 2.59, rtol=1.0e-2) #outflow velocity
     202        assert num.allclose(d, 0.563, rtol=1.0e-2) #depth at outlet used to calc v 
     203
     204    def Xtest_boyd_3(self):
     205        """test_boyd_3
     206       
     207        This tests the Boyd routine with data obtained from ??? by Petar Milevski   
     208        """
     209        # FIXME(Ole): This test fails (20 Feb 2009)
     210
     211        g=9.81
     212        culvert_slope=0.01  # Downward
     213
     214        inlet_depth=12.747
     215        outlet_depth=0.0
     216
     217        culvert_length=4.0
     218        culvert_width=0.75
     219        culvert_height=0.75
     220       
     221        culvert_type='pipe'
     222        manning=0.013
     223        sum_loss=1.5
     224
     225        inlet_specific_energy=inlet_depth #+0.5*v**2/g
     226        z_in = 0.0
     227        z_out = -culvert_length*culvert_slope/100
     228        E_in = z_in+inlet_depth  #+ 0.5*v**2/g
     229        E_out = z_out+outlet_depth  #+ 0.5*v**2/g
     230        delta_total_energy = E_in-E_out
     231
     232        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
     233                                                 outlet_depth,
     234                                                 inlet_specific_energy,
     235                                                 delta_total_energy,
     236                                                 g,
     237                                                 culvert_length,
     238                                                 culvert_width,
     239                                                 culvert_height,
     240                                                 culvert_type,
     241                                                 manning,
     242                                                 sum_loss)
     243       
     244        print Q, v, d
     245        assert num.allclose(Q, 5.00, rtol=1.0e-2) #inflow
     246        assert num.allclose(v, 11.022, rtol=1.0e-2) #outflow velocity
     247        assert num.allclose(d, 0.72, rtol=1.0e-2) #depth at outlet used to calc v
     248
     249    def Xtest_boyd_4(self):
     250        """test_boyd_4
     251       
     252        This tests the Boyd routine with data obtained from ??? by Petar Milevski   
     253        """
     254        # FIXME(Ole): This test fails (20 Feb 2009)
     255
     256        g=9.81
     257        culvert_slope=0.01  # Downward
     258
     259        inlet_depth=1.004
     260        outlet_depth=1.00
     261
     262        culvert_length=4.0
     263        culvert_width=0.75
     264        culvert_height=0.75
     265       
     266        culvert_type='pipe'
     267        manning=0.013
     268        sum_loss=1.5
     269
     270        inlet_specific_energy=inlet_depth #+0.5*v**2/g
     271        z_in = 0.0
     272        z_out = -culvert_length*culvert_slope/100
     273        E_in = z_in+inlet_depth  #+ 0.5*v**2/g
     274        E_out = z_out+outlet_depth  #+ 0.5*v**2/g
     275        delta_total_energy = E_in-E_out
     276
     277        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
     278                                                 outlet_depth,
     279                                                 inlet_specific_energy,
     280                                                 delta_total_energy,
     281                                                 g,
     282                                                 culvert_length,
     283                                                 culvert_width,
     284                                                 culvert_height,
     285                                                 culvert_type,
     286                                                 manning,
     287                                                 sum_loss)
     288       
     289        print Q, v, d
     290        assert num.allclose(Q, 0.10, rtol=1.0e-2) #inflow
     291        assert num.allclose(v, 0.22, rtol=1.0e-2) #outflow velocity
     292        assert num.allclose(d, 0.76, rtol=1.0e-2) #depth at outlet used to calc v
     293
     294    def Xtest_boyd_5(self):
     295        """test_boyd_5
     296       
     297        This tests the Boyd routine with data obtained from ??? by Petar Milevski   
     298        """
     299        # FIXME(Ole): This test fails (20 Feb 2009)
     300
     301        g=9.81
     302        culvert_slope=0.01  # Downward
     303
     304        inlet_depth=1.401
     305        outlet_depth=1.00
     306
     307        culvert_length=4.0
     308        culvert_width=0.75
     309        culvert_height=0.75
     310       
     311        culvert_type='pipe'
     312        manning=0.013
     313        sum_loss=1.5
     314
     315        inlet_specific_energy=inlet_depth #+0.5*v**2/g
     316        z_in = 0.0
     317        z_out = -culvert_length*culvert_slope/100
     318        E_in = z_in+inlet_depth  #+ 0.5*v**2/g
     319        E_out = z_out+outlet_depth  #+ 0.5*v**2/g
     320        delta_total_energy = E_in-E_out
     321
     322        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
     323                                                 outlet_depth,
     324                                                 inlet_specific_energy,
     325                                                 delta_total_energy,
     326                                                 g,
     327                                                 culvert_length,
     328                                                 culvert_width,
     329                                                 culvert_height,
     330                                                 culvert_type,
     331                                                 manning,
     332                                                 sum_loss)
     333       
     334        print Q, v, d
     335        assert num.allclose(Q, 1.00, rtol=1.0e-2) #inflow
     336        assert num.allclose(v, 2.204, rtol=1.0e-2) #outflow velocity
     337        assert num.allclose(d, 0.76, rtol=1.0e-2) #depth at outlet used to calc v
     338
     339
     340    def Xtest_boyd_6(self):
     341        """test_boyd_5
     342       
     343        This tests the Boyd routine with data obtained from ??? by Petar Milevski   
     344        """
     345        # FIXME(Ole): This test fails (20 Feb 2009)
     346
     347        g=9.81
     348        culvert_slope=0.01  # Downward
     349
     350        inlet_depth=12.747
     351        outlet_depth=1.00
     352
     353        culvert_length=4.0
     354        culvert_width=0.75
     355        culvert_height=0.75
     356       
     357        culvert_type='pipe'
     358        manning=0.013
     359        sum_loss=1.5
     360
     361        inlet_specific_energy=inlet_depth #+0.5*v**2/g
     362        z_in = 0.0
     363        z_out = -culvert_length*culvert_slope/100
     364        E_in = z_in+inlet_depth  #+ 0.5*v**2/g
     365        E_out = z_out+outlet_depth  #+ 0.5*v**2/g
     366        delta_total_energy = E_in-E_out
     367
     368        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
     369                                                 outlet_depth,
     370                                                 inlet_specific_energy,
     371                                                 delta_total_energy,
     372                                                 g,
     373                                                 culvert_length,
     374                                                 culvert_width,
     375                                                 culvert_height,
     376                                                 culvert_type,
     377                                                 manning,
     378                                                 sum_loss)
     379       
     380        print Q, v, d
     381        assert num.allclose(Q, 5.00, rtol=1.0e-2) #inflow
     382        assert num.allclose(v, 11.022, rtol=1.0e-2) #outflow velocity
     383        assert num.allclose(d, 0.76, rtol=1.0e-2) #depth at outlet used to calc v
     384
     385
     386    def Xtest_boyd_7(self):
     387        """test_boyd_7
     388       
     389        This tests the Boyd routine with data obtained from ??? by Petar Milevski   
     390        """
     391        # FIXME(Ole): This test fails (20 Feb 2009)
     392
     393        g=9.81
     394        culvert_slope=0.1  # Downward
     395
     396        inlet_depth=0.303
     397        outlet_depth=0.00
     398
     399        culvert_length=4.0
     400        culvert_width=0.75
     401        culvert_height=0.75
     402       
     403        culvert_type='pipe'
     404        manning=0.013
     405        sum_loss=1.5
     406
     407        inlet_specific_energy=inlet_depth #+0.5*v**2/g
     408        z_in = 0.0
     409        z_out = -culvert_length*culvert_slope/100
     410        E_in = z_in+inlet_depth  #+ 0.5*v**2/g
     411        E_out = z_out+outlet_depth  #+ 0.5*v**2/g
     412        delta_total_energy = E_in-E_out
     413
     414        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
     415                                                 outlet_depth,
     416                                                 inlet_specific_energy,
     417                                                 delta_total_energy,
     418                                                 g,
     419                                                 culvert_length,
     420                                                 culvert_width,
     421                                                 culvert_height,
     422                                                 culvert_type,
     423                                                 manning,
     424                                                 sum_loss)
     425       
     426        print Q, v, d
     427        assert num.allclose(Q, 0.10, rtol=1.0e-2) #inflow
     428        assert num.allclose(v, 1.13, rtol=1.0e-2) #outflow velocity
     429        assert num.allclose(d, 0.19, rtol=1.0e-2) #depth at outlet used to calc v
     430
     431
     432    def Xtest_boyd_8(self):
     433        """test_boyd_8
     434       
     435        This tests the Boyd routine with data obtained from ??? by Petar Milevski   
     436        """
     437        # FIXME(Ole): This test fails (20 Feb 2009)
     438
     439        g=9.81
     440        culvert_slope=0.1  # Downward
     441
     442        inlet_depth=1.135
     443        outlet_depth=0.00
     444
     445        culvert_length=4.0
     446        culvert_width=0.75
     447        culvert_height=0.75
     448       
     449        culvert_type='pipe'
     450        manning=0.013
     451        sum_loss=1.5
     452
     453        inlet_specific_energy=inlet_depth #+0.5*v**2/g
     454        z_in = 0.0
     455        z_out = -culvert_length*culvert_slope/100
     456        E_in = z_in+inlet_depth  #+ 0.5*v**2/g
     457        E_out = z_out+outlet_depth  #+ 0.5*v**2/g
     458        delta_total_energy = E_in-E_out
     459
     460        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
     461                                                 outlet_depth,
     462                                                 inlet_specific_energy,
     463                                                 delta_total_energy,
     464                                                 g,
     465                                                 culvert_length,
     466                                                 culvert_width,
     467                                                 culvert_height,
     468                                                 culvert_type,
     469                                                 manning,
     470                                                 sum_loss)
     471       
     472        print Q, v, d
     473        assert num.allclose(Q, 1.00, rtol=1.0e-2) #inflow
     474        assert num.allclose(v, 2.204, rtol=1.0e-2) #outflow velocity
     475        assert num.allclose(d, 0.76, rtol=1.0e-2) #depth at outlet used to calc v
     476
     477    def Xtest_boyd_9(self):
     478        """test_boyd_9
     479       
     480        This tests the Boyd routine with data obtained from ??? by Petar Milevski   
     481        """
     482        # FIXME(Ole): This test fails (20 Feb 2009)
     483
     484        g=9.81
     485        culvert_slope=0.1  # Downward
     486
     487        inlet_depth=1.1504
     488        outlet_depth=1.5
     489
     490        culvert_length=4.0
     491        culvert_width=0.75
     492        culvert_height=0.75
     493       
     494        culvert_type='pipe'
     495        manning=0.013
     496        sum_loss=1.5
     497
     498        inlet_specific_energy=inlet_depth #+0.5*v**2/g
     499        z_in = 0.0
     500        z_out = -culvert_length*culvert_slope/100
     501        E_in = z_in+inlet_depth  #+ 0.5*v**2/g
     502        E_out = z_out+outlet_depth  #+ 0.5*v**2/g
     503        delta_total_energy = E_in-E_out
     504
     505        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
     506                                                 outlet_depth,
     507                                                 inlet_specific_energy,
     508                                                 delta_total_energy,
     509                                                 g,
     510                                                 culvert_length,
     511                                                 culvert_width,
     512                                                 culvert_height,
     513                                                 culvert_type,
     514                                                 manning,
     515                                                 sum_loss)
     516       
     517        print Q, v, d
     518        assert num.allclose(Q, 0.10, rtol=1.0e-2) #inflow
     519        assert num.allclose(v, 0.22, rtol=1.0e-2) #outflow velocity
     520        assert num.allclose(d, 0.76, rtol=1.0e-2) #depth at outlet used to calc v
     521
     522
     523    def Xtest_boyd_10(self):
     524        """test_boyd_9
     525       
     526        This tests the Boyd routine with data obtained from ??? by Petar Milevski   
     527        """
     528        # FIXME(Ole): This test fails (20 Feb 2009)
     529
     530        g=9.81
     531        culvert_slope=0.1  # Downward
     532
     533        inlet_depth=1.901
     534        outlet_depth=1.5
     535
     536        culvert_length=4.0
     537        culvert_width=0.75
     538        culvert_height=0.75
     539       
     540        culvert_type='pipe'
     541        manning=0.013
     542        sum_loss=1.5
     543
     544        inlet_specific_energy=inlet_depth #+0.5*v**2/g
     545        z_in = 0.0
     546        z_out = -culvert_length*culvert_slope/100
     547        E_in = z_in+inlet_depth  #+ 0.5*v**2/g
     548        E_out = z_out+outlet_depth  #+ 0.5*v**2/g
     549        delta_total_energy = E_in-E_out
     550
     551        Q, v, d = boyd_generalised_culvert_model(inlet_depth,
     552                                                 outlet_depth,
     553                                                 inlet_specific_energy,
     554                                                 delta_total_energy,
     555                                                 g,
     556                                                 culvert_length,
     557                                                 culvert_width,
     558                                                 culvert_height,
     559                                                 culvert_type,
     560                                                 manning,
     561                                                 sum_loss)
     562       
     563        print Q, v, d
     564        assert num.allclose(Q, 1.00, rtol=1.0e-2) #inflow
     565        assert num.allclose(v, 2.204, rtol=1.0e-2) #outflow velocity
     566        assert num.allclose(d, 0.76, rtol=1.0e-2) #depth at outlet used to calc v
    114567   
    115568               
  • branches/numpy/anuga/geospatial_data/geospatial_data.py

    r6768 r6902  
    862862                 self.file_pointer) = _read_csv_file_blocking(self.file_pointer,
    863863                                                              self.header[:],
    864                                                               max_read_lines= \
     864                                                              max_read_lines=
    865865                                                           self.max_read_lines,
    866                                                               verbose= \
     866                                                              verbose=
    867867                                                                  self.verbose)
    868868
  • branches/numpy/anuga/load_mesh/loadASCII.py

    r6553 r6902  
    817817        mesh['vertex_attributes'] = None
    818818
    819     mesh['vertex_attribute_titles'] = []
     819#    mesh['vertex_attribute_titles'] = []
    820820    try:
    821821        titles = fid.variables['vertex_attribute_titles'][:]
     
    829829        mesh['segments'] = num.array([], num.int)      #array default#
    830830
    831     mesh['segment_tags'] = []
     831#    mesh['segment_tags'] = []
    832832    try:
    833833        tags = fid.variables['segment_tags'][:]
     
    844844        mesh['triangle_neighbors'] = num.array([], num.int)     #array default#
    845845
    846     mesh['triangle_tags'] = []
     846#    mesh['triangle_tags'] = []
    847847    try:
    848848        tags = fid.variables['triangle_tags'][:]
  • branches/numpy/anuga/load_mesh/test_loadASCII.py

    r6553 r6902  
    410410        self.failUnless(num.alltrue(loaded_dict['vertex_attributes'] ==
    411411                                    dict['vertex_attributes']) or
    412                         (loaded_dict['vertex_attributes'] == None and
    413                          dict['vertex_attributes'] == []),
     412                                    (loaded_dict['vertex_attributes'] ==
     413                                    None and
     414                                    dict['vertex_attributes'] == []),
    414415                        fail_string + ' failed!! Test vertex_attributes')
    415416
     
    423424            self.failUnless(seg == ldseg, fail_string + ' failed\n' + msg)
    424425
    425         dict_array = num.array(dict['vertex_attribute_titles'])
    426         loaded_dict_array = num.array(loaded_dict['vertex_attribute_titles'])
     426#        dict_array = num.array(dict['vertex_attribute_titles'])
     427#        loaded_dict_array = num.array(loaded_dict['vertex_attribute_titles'])
    427428        try:
    428429            assert num.allclose(dict_array, loaded_dict_array)
  • branches/numpy/anuga/shallow_water/data_manager.py

    r6689 r6902  
    35753575    other_quantities.remove('z')
    35763576    other_quantities.remove('volumes')
    3577     #try:
    3578     #    other_quantities.remove('stage_range')
    3579     #    other_quantities.remove('xmomentum_range')
    3580     #    other_quantities.remove('ymomentum_range')
    3581     #    other_quantities.remove('elevation_range')
    3582     #except:
    3583     #    pass
     3577    try:
     3578        other_quantities.remove('stage_range')
     3579        other_quantities.remove('xmomentum_range')
     3580        other_quantities.remove('ymomentum_range')
     3581        other_quantities.remove('elevation_range')
     3582    except:
     3583        pass
    35843584
    35853585    conserved_quantities.remove('time')
     
    53385338        j += 1
    53395339
    5340     #if verbose: sww.verbose_quantities(outfile)
     5340    if verbose: sww.verbose_quantities(outfile)
    53415341
    53425342    outfile.close()
     
    58555855    # @param smoothing True if smoothing is to be used.
    58565856    # @param order
    5857     # @param sww_precision Data type of the quantitiy written (netcdf constant)
     5857    # @param sww_precision Data type of the quantity written (netcdf constant)
    58585858    # @param verbose True if this function is to be verbose.
    58595859    # @note If 'times' is a list, the info will be made relative.
     
    59365936                               ('number_of_points',))
    59375937        q = 'elevation'
    5938         #outfile.createVariable(q + Write_sww.RANGE, sww_precision,
    5939         #                       ('numbers_in_range',))
    5940 
    5941 
    5942         ## Initialise ranges with small and large sentinels.
    5943         ## If this was in pure Python we could have used None sensibly
    5944         #outfile.variables[q+Write_sww.RANGE][0] = max_float  # Min
    5945         #outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
     5938        outfile.createVariable(q + Write_sww.RANGE, sww_precision,
     5939                               ('numbers_in_range',))
     5940
     5941        # Initialise ranges with small and large sentinels.
     5942        # If this was in pure Python we could have used None sensibly
     5943        outfile.variables[q+Write_sww.RANGE][0] = max_float  # Min
     5944        outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
    59465945
    59475946        # FIXME: Backwards compatibility
     
    59585957            outfile.createVariable(q, sww_precision, ('number_of_timesteps',
    59595958                                                      'number_of_points'))
    5960             #outfile.createVariable(q + Write_sww.RANGE, sww_precision,
    5961             #                       ('numbers_in_range',))
     5959            outfile.createVariable(q + Write_sww.RANGE, sww_precision,
     5960                                   ('numbers_in_range',))
    59625961
    59635962            # Initialise ranges with small and large sentinels.
    59645963            # If this was in pure Python we could have used None sensibly
    5965             #outfile.variables[q+Write_sww.RANGE][0] = max_float  # Min
    5966             #outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
     5964            outfile.variables[q+Write_sww.RANGE][0] = max_float  # Min
     5965            outfile.variables[q+Write_sww.RANGE][1] = -max_float # Max
    59675966
    59685967        if isinstance(times, (list, num.ndarray)):
     
    60736072        q = 'elevation'
    60746073        # This updates the _range values
    6075         #outfile.variables[q + Write_sww.RANGE][0] = min(elevation)
    6076         #outfile.variables[q + Write_sww.RANGE][1] = max(elevation)
     6074        outfile.variables[q + Write_sww.RANGE][0] = num.min(elevation)
     6075        outfile.variables[q + Write_sww.RANGE][1] = num.max(elevation)
    60776076
    60786077
     
    61266125                                q_values.astype(sww_precision)
    61276126       
    6128         #        # This updates the _range values
    6129         #        q_range = outfile.variables[q + Write_sww.RANGE][:]
    6130         #        q_values_min = min(q_values)
    6131         #        if q_values_min < q_range[0]:
    6132         #            outfile.variables[q + Write_sww.RANGE][0] = q_values_min
    6133         #        q_values_max = max(q_values)
    6134         #        if q_values_max > q_range[1]:
    6135         #            outfile.variables[q + Write_sww.RANGE][1] = q_values_max
     6127                # This updates the _range values
     6128                q_range = outfile.variables[q + Write_sww.RANGE][:]
     6129                q_values_min = num.min(q_values)
     6130                if q_values_min < q_range[0]:
     6131                    outfile.variables[q + Write_sww.RANGE][0] = q_values_min
     6132                q_values_max = num.max(q_values)
     6133                if q_values_max > q_range[1]:
     6134                    outfile.variables[q + Write_sww.RANGE][1] = q_values_max
    61366135
    61376136    ##
    61386137    # @brief Print the quantities in the underlying file.
    61396138    # @param outfile UNUSED.
    6140     #def verbose_quantities(self, outfile):
    6141     #    print '------------------------------------------------'
    6142     #    print 'More Statistics:'
    6143     #    for q in Write_sww.sww_quantities:
    6144     #        print '  %s in [%f, %f]' % (q,
    6145     #                                    outfile.variables[q+Write_sww.RANGE][0],
    6146     #                                    outfile.variables[q+Write_sww.RANGE][1])
    6147     #    print '------------------------------------------------'
     6139    def verbose_quantities(self, outfile):
     6140        print '------------------------------------------------'
     6141        print 'More Statistics:'
     6142        for q in Write_sww.sww_quantities:
     6143            print '  %s in [%f, %f]' % (q,
     6144                                        outfile.variables[q+Write_sww.RANGE][0],
     6145                                        outfile.variables[q+Write_sww.RANGE][1])
     6146        print '------------------------------------------------'
    61486147
    61496148
     
    64966495
    64976496                # This updates the _range values
    6498                 #q_range = outfile.variables[q + Write_sts.RANGE][:]
    6499                 #q_values_min = min(q_values)
    6500                 #if q_values_min < q_range[0]:
    6501                 #    outfile.variables[q + Write_sts.RANGE][0] = q_values_min
    6502                 #q_values_max = max(q_values)
    6503                 #if q_values_max > q_range[1]:
    6504                 #    outfile.variables[q + Write_sts.RANGE][1] = q_values_max
     6497                q_range = outfile.variables[q + Write_sts.RANGE][:]
     6498                q_values_min = num.min(q_values)
     6499                if q_values_min < q_range[0]:
     6500                    outfile.variables[q + Write_sts.RANGE][0] = q_values_min
     6501                q_values_max = num.max(q_values)
     6502                if q_values_max > q_range[1]:
     6503                    outfile.variables[q + Write_sts.RANGE][1] = q_values_max
    65056504
    65066505
     
    66736672
    66746673
    6675 # FIXME (DSG): Add unit test, make general, not just 2 files,
    6676 # but any number of files.
    66776674##
    66786675# @brief Copy a file to a directory, and optionally append another file to it.
     
    66806677# @param filename Path to file to copy to directory 'dir_name'.
    66816678# @param filename2 Optional path to file to append to copied file.
    6682 def copy_code_files(dir_name, filename1, filename2=None):
     6679# @param verbose True if this function is to be verbose.
     6680# @note Allow filenames to be either a string or sequence of strings.
     6681def copy_code_files(dir_name, filename1, filename2=None, verbose=False):
    66836682    """Copies "filename1" and "filename2" to "dir_name".
    6684     Very useful for information management
    6685     filename1 and filename2 are both absolute pathnames
    6686     """
    6687 
     6683
     6684    Each 'filename' may be a string or list of filename strings.
     6685
     6686    Filenames must be absolute pathnames
     6687    """
     6688
     6689    ##
     6690    # @brief copies a file or sequence to destination directory.
     6691    # @param dest The destination directory to copy to.
     6692    # @param file A filename string or sequence of filename strings.
     6693    def copy_file_or_sequence(dest, file):
     6694        if hasattr(file, '__iter__'):
     6695            for f in file:
     6696                shutil.copy(f, dir_name)
     6697                if verbose:
     6698                    print 'File %s copied' % (f)
     6699        else:
     6700            shutil.copy(file, dir_name)
     6701            if verbose:
     6702                print 'File %s copied' % (file)
     6703
     6704    # check we have a destination directory, create if necessary
    66886705    if access(dir_name, F_OK) == 0:
    6689         print 'Make directory %s' % dir_name
    6690         mkdir (dir_name,0777)
    6691 
    6692     shutil.copy(filename1, dir_name + sep + basename(filename1))
    6693 
    6694     if filename2 != None:
    6695         shutil.copy(filename2, dir_name + sep + basename(filename2))
    6696         print 'Files %s and %s copied' %(filename1, filename2)
    6697     else:
    6698         print 'File %s copied' %(filename1)
     6706        if verbose:
     6707            print 'Make directory %s' % dir_name
     6708        mkdir(dir_name, 0777)
     6709
     6710    copy_file_or_sequence(dir_name, filename1)
     6711
     6712    if not filename2 is None:
     6713        copy_file_or_sequence(dir_name, filename2)
    66996714
    67006715
  • branches/numpy/anuga/shallow_water/shallow_water_ext.c

    r6780 r6902  
    2828
    2929// Computational function for rotation
     30// FIXME: Perhaps inline this and profile
    3031int _rotate(double *q, double n1, double n2) {
    3132  /*Rotate the momentum component q (q[1], q[2])
     
    5253
    5354int find_qmin_and_qmax(double dq0, double dq1, double dq2,
    54                        double *qmin, double *qmax){
     55               double *qmin, double *qmax){
    5556  // Considering the centroid of an FV triangle and the vertices of its
    5657  // auxiliary triangle, find
     
    6768    if (dq1>=dq2){
    6869      if (dq1>=0.0)
    69         *qmax=dq0+dq1;
     70    *qmax=dq0+dq1;
    7071      else
    71         *qmax=dq0;
    72        
     72    *qmax=dq0;
     73   
    7374      *qmin=dq0+dq2;
    7475      if (*qmin>=0.0) *qmin = 0.0;
     
    7677    else{// dq1<dq2
    7778      if (dq2>0)
    78         *qmax=dq0+dq2;
     79    *qmax=dq0+dq2;
    7980      else
    80         *qmax=dq0;
    81        
    82       *qmin=dq0+dq1;   
     81    *qmax=dq0;
     82   
     83      *qmin=dq0+dq1;   
    8384      if (*qmin>=0.0) *qmin=0.0;
    8485    }
     
    8788    if (dq1<=dq2){
    8889      if (dq1<0.0)
    89         *qmin=dq0+dq1;
     90    *qmin=dq0+dq1;
    9091      else
    91         *qmin=dq0;
    92        
    93       *qmax=dq0+dq2;   
     92    *qmin=dq0;
     93   
     94      *qmax=dq0+dq2;   
    9495      if (*qmax<=0.0) *qmax=0.0;
    9596    }
    9697    else{// dq1>dq2
    9798      if (dq2<0.0)
    98         *qmin=dq0+dq2;
     99    *qmin=dq0+dq2;
    99100      else
    100         *qmin=dq0;
    101        
     101    *qmin=dq0;
     102   
    102103      *qmax=dq0+dq1;
    103104      if (*qmax<=0.0) *qmax=0.0;
     
    133134 
    134135  phi=min(r*beta_w,1.0);
    135   for (i=0;i<3;i++)
    136     dqv[i]=dqv[i]*phi;
     136  //for (i=0;i<3;i++)
     137  dqv[0]=dqv[0]*phi;
     138  dqv[1]=dqv[1]*phi;
     139  dqv[2]=dqv[2]*phi;
    137140   
    138141  return 0;
     
    141144
    142145double compute_froude_number(double uh,
    143                              double h,
    144                              double g,
    145                              double epsilon) {
    146                          
     146                 double h,
     147                 double g,
     148                 double epsilon) {
     149             
    147150  // Compute Froude number; v/sqrt(gh)
    148151  // FIXME (Ole): Not currently in use
     
    166169// This is used by flux functions
    167170// Input parameters uh and h may be modified by this function.
     171
     172// FIXME: Perhaps inline this and profile
    168173double _compute_speed(double *uh,
    169                       double *h,
    170                       double epsilon,
    171                       double h0) {
     174              double *h,
     175              double epsilon,
     176              double h0) {
    172177 
    173178  double u;
     
    188193}
    189194
     195
     196// Optimised squareroot computation
     197//
     198//See
     199//http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf
     200//and http://mail.opensolaris.org/pipermail/tools-gcc/2005-August/000047.html
     201float fast_squareroot_approximation(float number) {
     202  float x;
     203  const float f = 1.5;
     204
     205  // Allow i and y to occupy the same memory
     206  union
     207  {
     208    long i;
     209    float y;
     210  } u;
     211 
     212  // Good initial guess 
     213  x = number * 0.5; 
     214  u.y  = number;
     215  u.i  = 0x5f3759df - ( u.i >> 1 );
     216 
     217  // Take a few iterations
     218  u.y  = u.y * ( f - ( x * u.y * u.y ) );
     219  u.y  = u.y * ( f - ( x * u.y * u.y ) );
     220   
     221  return number * u.y;
     222}
     223
     224
     225
     226// Optimised squareroot computation (double version, slower)
     227double Xfast_squareroot_approximation(double number) {
     228  double x;
     229  const double f = 1.5;
     230   
     231  // Allow i and y to occupy the same memory   
     232  union
     233  {
     234    long long i;
     235    double y;
     236  } u;
     237
     238  // Good initial guess
     239  x = number * 0.5;
     240  u.y  = number;
     241  u.i  = 0x5fe6ec85e7de30daLL - ( u.i >> 1 );
     242 
     243  // Take a few iterations 
     244  u.y  = u.y * ( f - ( x * u.y * u.y ) );
     245  u.y  = u.y * ( f - ( x * u.y * u.y ) );
     246
     247  return number * u.y;
     248}
     249
     250
    190251// Innermost flux function (using stage w=z+h)
    191252int _flux_function_central(double *q_left, double *q_right,
    192                            double z_left, double z_right,
    193                            double n1, double n2,
    194                            double epsilon, double H0, double g,
    195                            double *edgeflux, double *max_speed) {
     253                           double z_left, double z_right,
     254                           double n1, double n2,
     255                           double epsilon, double H0, double g,
     256                           double *edgeflux, double *max_speed)
     257{
    196258
    197259  /*Compute fluxes between volumes for the shallow water wave equation
     
    212274  double v_left, v_right; 
    213275  double s_min, s_max, soundspeed_left, soundspeed_right;
    214   double denom, z;
     276  double denom, inverse_denominator, z;
    215277
    216278  // Workspace (allocate once, use many)
     
    219281  double h0 = H0*H0; // This ensures a good balance when h approaches H0.
    220282                     // But evidence suggests that h0 can be as little as
    221                      // epsilon!
     283             // epsilon!
    222284 
    223285  // Copy conserved quantities to protect from modification
    224   for (i=0; i<3; i++) {
    225     q_left_rotated[i] = q_left[i];
    226     q_right_rotated[i] = q_right[i];
    227   }
     286  q_left_rotated[0] = q_left[0];
     287  q_right_rotated[0] = q_right[0];
     288  q_left_rotated[1] = q_left[1];
     289  q_right_rotated[1] = q_right[1];
     290  q_left_rotated[2] = q_left[2];
     291  q_right_rotated[2] = q_right[2];
    228292
    229293  // Align x- and y-momentum with x-axis
     
    231295  _rotate(q_right_rotated, n1, n2);
    232296
    233   z = (z_left+z_right)/2; // Average elevation values.
    234                           // Even though this will nominally allow for discontinuities
    235                           // in the elevation data, there is currently no numerical
    236                           // support for this so results may be strange near jumps in the bed.
     297  z = 0.5*(z_left + z_right); // Average elevation values.
     298                            // Even though this will nominally allow for discontinuities
     299                            // in the elevation data, there is currently no numerical
     300                            // support for this so results may be strange near jumps in the bed.
    237301
    238302  // Compute speeds in x-direction
    239303  w_left = q_left_rotated[0];         
    240   h_left = w_left-z;
     304  h_left = w_left - z;
    241305  uh_left = q_left_rotated[1];
    242306  u_left = _compute_speed(&uh_left, &h_left, epsilon, h0);
    243307
    244308  w_right = q_right_rotated[0];
    245   h_right = w_right-z;
     309  h_right = w_right - z;
    246310  uh_right = q_right_rotated[1];
    247311  u_right = _compute_speed(&uh_right, &h_right, epsilon, h0); 
     
    257321  // Maximal and minimal wave speeds
    258322  soundspeed_left  = sqrt(g*h_left);
    259   soundspeed_right = sqrt(g*h_right);
    260 
    261   s_max = max(u_left+soundspeed_left, u_right+soundspeed_right);
    262   if (s_max < 0.0) s_max = 0.0;
    263 
    264   s_min = min(u_left-soundspeed_left, u_right-soundspeed_right);
    265   if (s_min > 0.0) s_min = 0.0;
    266 
     323  soundspeed_right = sqrt(g*h_right); 
     324 
     325  // Code to use fast square root optimisation if desired.
     326  //soundspeed_left  = fast_squareroot_approximation(g*h_left);
     327  //soundspeed_right = fast_squareroot_approximation(g*h_right);
     328
     329  s_max = max(u_left + soundspeed_left, u_right + soundspeed_right);
     330  if (s_max < 0.0)
     331  {
     332    s_max = 0.0;
     333  }
     334
     335  s_min = min(u_left - soundspeed_left, u_right - soundspeed_right);
     336  if (s_min > 0.0)
     337  {
     338    s_min = 0.0;
     339  }
    267340 
    268341  // Flux formulas
     
    275348  flux_right[2] = u_right*vh_right;
    276349
    277  
    278350  // Flux computation
    279   denom = s_max-s_min;
    280   if (denom < epsilon) { // FIXME (Ole): Try using H0 here
    281     for (i=0; i<3; i++) edgeflux[i] = 0.0;
     351  denom = s_max - s_min;
     352  if (denom < epsilon)
     353  { // FIXME (Ole): Try using H0 here
     354    memset(edgeflux, 0, 3*sizeof(double));
    282355    *max_speed = 0.0;
    283   } else {
    284     for (i=0; i<3; i++) {
     356  }
     357  else
     358  {
     359    inverse_denominator = 1.0/denom;
     360    for (i = 0; i < 3; i++)
     361    {
    285362      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
    286       edgeflux[i] += s_max*s_min*(q_right_rotated[i]-q_left_rotated[i]);
    287       edgeflux[i] /= denom;
     363      edgeflux[i] += s_max*s_min*(q_right_rotated[i] - q_left_rotated[i]);
     364      edgeflux[i] *= inverse_denominator;
    288365    }
    289366
     
    315392// Computational function for flux computation (using stage w=z+h)
    316393int flux_function_kinetic(double *q_left, double *q_right,
    317                   double z_left, double z_right,
    318                   double n1, double n2,
    319                   double epsilon, double H0, double g,
    320                   double *edgeflux, double *max_speed) {
     394          double z_left, double z_right,
     395          double n1, double n2,
     396          double epsilon, double H0, double g,
     397          double *edgeflux, double *max_speed) {
    321398
    322399  /*Compute fluxes between volumes for the shallow water wave equation
     
    413490
    414491void _manning_friction(double g, double eps, int N,
    415                        double* w, double* z,
    416                        double* uh, double* vh,
    417                        double* eta, double* xmom, double* ymom) {
     492               double* w, double* z,
     493               double* uh, double* vh,
     494               double* eta, double* xmom, double* ymom) {
    418495
    419496  int k;
     
    426503        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
    427504        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
    428         //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
     505        //S /= exp((7.0/3.0)*log(h));      //seems to save about 15% over manning_friction
    429506        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
    430507
     
    441518/*
    442519void _manning_friction_explicit(double g, double eps, int N,
    443                        double* w, double* z,
    444                        double* uh, double* vh,
    445                        double* eta, double* xmom, double* ymom) {
     520               double* w, double* z,
     521               double* uh, double* vh,
     522               double* eta, double* xmom, double* ymom) {
    446523
    447524  int k;
     
    452529      h = w[k]-z[k];
    453530      if (h >= eps) {
    454         S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
    455         S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
    456         //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
    457         //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
    458 
    459 
    460         //Update momentum
    461         xmom[k] += S*uh[k];
    462         ymom[k] += S*vh[k];
     531    S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
     532    S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
     533    //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
     534    //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
     535
     536
     537    //Update momentum
     538    xmom[k] += S*uh[k];
     539    ymom[k] += S*vh[k];
    463540      }
    464541    }
     
    471548
    472549double velocity_balance(double uh_i,
    473                         double uh,
    474                         double h_i,
    475                         double h,
    476                         double alpha,
    477                         double epsilon) {
     550            double uh,
     551            double h_i,
     552            double h,
     553            double alpha,
     554            double epsilon) {
    478555  // Find alpha such that speed at the vertex is within one
    479   // order of magnitude of the centroid speed   
     556  // order of magnitude of the centroid speed   
    480557
    481558  // FIXME(Ole): Work in progress
     
    486563 
    487564  printf("alpha = %f, uh_i=%f, uh=%f, h_i=%f, h=%f\n",
    488         alpha, uh_i, uh, h_i, h);
     565    alpha, uh_i, uh, h_i, h);
    489566     
    490567 
     
    500577  if (h < epsilon) {
    501578    b = 1.0e10; // Limit
    502   } else {     
     579  } else { 
    503580    b = m*fabs(h_i - h)/h;
    504581  }
     
    507584
    508585  if (a > b) {
    509     estimate = (m-1)/(a-b);             
     586    estimate = (m-1)/(a-b);     
    510587   
    511588    printf("Alpha %f, estimate %f\n",
    512            alpha, estimate);   
    513            
     589       alpha, estimate);   
     590       
    514591    if (alpha < estimate) {
    515592      printf("Adjusting alpha from %f to %f\n",
    516              alpha, estimate);
     593         alpha, estimate);
    517594      alpha = estimate;
    518595    }
     
    520597 
    521598    if (h < h_i) {
    522       estimate = (m-1)/(a-b);                
     599      estimate = (m-1)/(a-b);            
    523600   
    524601      printf("Alpha %f, estimate %f\n",
    525              alpha, estimate);   
    526            
     602         alpha, estimate);   
     603       
    527604      if (alpha < estimate) {
    528         printf("Adjusting alpha from %f to %f\n",
    529                alpha, estimate);
    530         alpha = estimate;
     605    printf("Adjusting alpha from %f to %f\n",
     606           alpha, estimate);
     607    alpha = estimate;
    531608      }   
    532609    }
     
    540617
    541618int _balance_deep_and_shallow(int N,
    542                               double* wc,
    543                               double* zc,
    544                               double* wv,
    545                               double* zv,
    546                               double* hvbar, // Retire this
    547                               double* xmomc,
    548                               double* ymomc,
    549                               double* xmomv,
    550                               double* ymomv,
    551                               double H0,
    552                               int tight_slope_limiters,
    553                               int use_centroid_velocities,
    554                               double alpha_balance) {
     619                  double* wc,
     620                  double* zc,
     621                  double* wv,
     622                  double* zv,
     623                  double* hvbar, // Retire this
     624                  double* xmomc,
     625                  double* ymomc,
     626                  double* xmomv,
     627                  double* ymomv,
     628                  double H0,
     629                  int tight_slope_limiters,
     630                  int use_centroid_velocities,
     631                  double alpha_balance) {
    555632
    556633  int k, k3, i;
     
    579656      // FIXME: Try with this one precomputed
    580657      for (i=0; i<3; i++) {
    581         dz = max(dz, fabs(zv[k3+i]-zc[k]));
     658    dz = max(dz, fabs(zv[k3+i]-zc[k]));
    582659      }
    583660    }
     
    592669
    593670    //if (hmin < 0.0 ) {
    594     //  printf("hmin = %f\n",hmin);
     671    //  printf("hmin = %f\n",hmin);
    595672    //}
    596673   
     
    607684     
    608685      if (dz > 0.0) {
    609         alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );     
     686    alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 );     
    610687      } else {
    611         alpha = 1.0;  // Flat bed
     688    alpha = 1.0;  // Flat bed
    612689      }
    613690      //printf("Using old style limiter\n");
     
    622699   
    623700      if (hmin < H0) {
    624         alpha = 1.0;
    625         for (i=0; i<3; i++) {
    626          
    627           h_diff = hc_k - hv[i];         
    628           if (h_diff <= 0) {
    629             // Deep water triangle is further away from bed than
    630             // shallow water (hbar < h). Any alpha will do
    631          
    632           } else { 
    633             // Denominator is positive which means that we need some of the
    634             // h-limited stage.
    635            
    636             alpha = min(alpha, (hc_k - H0)/h_diff);        
    637           }
    638         }
    639 
    640         // Ensure alpha in [0,1]
    641         if (alpha>1.0) alpha=1.0;
    642         if (alpha<0.0) alpha=0.0;
    643        
     701    alpha = 1.0;
     702    for (i=0; i<3; i++) {
     703     
     704      h_diff = hc_k - hv[i];     
     705      if (h_diff <= 0) {
     706        // Deep water triangle is further away from bed than
     707        // shallow water (hbar < h). Any alpha will do
     708     
     709      } else { 
     710        // Denominator is positive which means that we need some of the
     711        // h-limited stage.
     712       
     713        alpha = min(alpha, (hc_k - H0)/h_diff);    
     714      }
     715    }
     716
     717    // Ensure alpha in [0,1]
     718    if (alpha>1.0) alpha=1.0;
     719    if (alpha<0.0) alpha=0.0;
     720   
    644721      } else {
    645         // Use w-limited stage exclusively in deeper water.
    646         alpha = 1.0;       
     722    // Use w-limited stage exclusively in deeper water.
     723    alpha = 1.0;       
    647724      }
    648725    }
    649726
    650727
    651     //  Let
     728    //  Let
    652729    //
    653     //    wvi be the w-limited stage (wvi = zvi + hvi)
    654     //    wvi- be the h-limited state (wvi- = zvi + hvi-)
     730    //    wvi be the w-limited stage (wvi = zvi + hvi)
     731    //    wvi- be the h-limited state (wvi- = zvi + hvi-)
    655732    //
    656733    //
    657     //  where i=0,1,2 denotes the vertex ids
     734    //  where i=0,1,2 denotes the vertex ids
    658735    //
    659736    //  Weighted balance between w-limited and h-limited stage is
    660737    //
    661     //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
     738    //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
    662739    //
    663740    //  It follows that the updated wvi is
    664741    //    wvi := zvi + (1-alpha)*hvi- + alpha*hvi
    665742    //
    666     //  Momentum is balanced between constant and limited
     743    //  Momentum is balanced between constant and limited
    667744
    668745
     
    670747      for (i=0; i<3; i++) {
    671748     
    672         wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i];     
    673 
    674         // Update momentum at vertices
    675         if (use_centroid_velocities == 1) {
    676           // This is a simple, efficient and robust option
    677           // It uses first order approximation of velocities, but retains
    678           // the order used by stage.
    679        
    680           // Speeds at centroids
    681           if (hc_k > epsilon) {
    682             uc = xmomc[k]/hc_k;
    683             vc = ymomc[k]/hc_k;
    684           } else {
    685             uc = 0.0;
    686             vc = 0.0;
    687           }
    688          
    689           // Vertex momenta guaranteed to be consistent with depth guaranteeing
    690           // controlled speed
    691           hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
    692           xmomv[k3+i] = uc*hv[i];
    693           ymomv[k3+i] = vc*hv[i];
    694          
    695         } else {
    696           // Update momentum as a linear combination of
    697           // xmomc and ymomc (shallow) and momentum
    698           // from extrapolator xmomv and ymomv (deep).
    699           // This assumes that values from xmomv and ymomv have
    700           // been established e.g. by the gradient limiter.
    701 
    702           // FIXME (Ole): I think this should be used with vertex momenta
    703           // computed above using centroid_velocities instead of xmomc
    704           // and ymomc as they'll be more representative first order
    705           // values.
    706          
    707           xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
    708           ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
    709        
    710         }
     749    wv[k3+i] = zv[k3+i] + (1-alpha)*hc_k + alpha*hv[i];
     750
     751    // Update momentum at vertices
     752    if (use_centroid_velocities == 1) {
     753      // This is a simple, efficient and robust option
     754      // It uses first order approximation of velocities, but retains
     755      // the order used by stage.
     756   
     757      // Speeds at centroids
     758      if (hc_k > epsilon) {
     759        uc = xmomc[k]/hc_k;
     760        vc = ymomc[k]/hc_k;
     761      } else {
     762        uc = 0.0;
     763        vc = 0.0;
     764      }
     765     
     766      // Vertex momenta guaranteed to be consistent with depth guaranteeing
     767      // controlled speed
     768      hv[i] = wv[k3+i] - zv[k3+i]; // Recompute (balanced) vertex depth
     769      xmomv[k3+i] = uc*hv[i];
     770      ymomv[k3+i] = vc*hv[i];
     771     
     772    } else {
     773      // Update momentum as a linear combination of
     774      // xmomc and ymomc (shallow) and momentum
     775      // from extrapolator xmomv and ymomv (deep).
     776      // This assumes that values from xmomv and ymomv have
     777      // been established e.g. by the gradient limiter.
     778
     779      // FIXME (Ole): I think this should be used with vertex momenta
     780      // computed above using centroid_velocities instead of xmomc
     781      // and ymomc as they'll be more representative first order
     782      // values.
     783     
     784      xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
     785      ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
     786   
     787    }
    711788      }
    712789    }
     
    719796
    720797int _protect(int N,
    721              double minimum_allowed_height,
    722              double maximum_allowed_speed,
    723              double epsilon,
    724              double* wc,
    725              double* zc,
    726              double* xmomc,
    727              double* ymomc) {
     798         double minimum_allowed_height,
     799         double maximum_allowed_speed,
     800         double epsilon,
     801         double* wc,
     802         double* zc,
     803         double* xmomc,
     804         double* ymomc) {
    728805
    729806  int k;
     
    738815
    739816      if (hc < minimum_allowed_height) {
    740        
    741         // Set momentum to zero and ensure h is non negative
    742         xmomc[k] = 0.0;
    743         ymomc[k] = 0.0;
    744         if (hc <= 0.0) wc[k] = zc[k];
     817       
     818    // Set momentum to zero and ensure h is non negative
     819    xmomc[k] = 0.0;
     820    ymomc[k] = 0.0;
     821    if (hc <= 0.0) wc[k] = zc[k];
    745822      }
    746823    }
     
    757834             
    758835        if (hc <= 0.0) {
    759                 wc[k] = zc[k];
    760         xmomc[k] = 0.0;
    761         ymomc[k] = 0.0;
     836            wc[k] = zc[k];
     837        xmomc[k] = 0.0;
     838        ymomc[k] = 0.0;
    762839        } else {
    763840          //Reduce excessive speeds derived from division by small hc
    764         //FIXME (Ole): This may be unnecessary with new slope limiters
    765         //in effect.
     841        //FIXME (Ole): This may be unnecessary with new slope limiters
     842        //in effect.
    766843         
    767844          u = xmomc[k]/hc;
    768           if (fabs(u) > maximum_allowed_speed) {
    769             reduced_speed = maximum_allowed_speed * u/fabs(u);
    770             //printf("Speed (u) has been reduced from %.3f to %.3f\n",
    771             //  u, reduced_speed);
    772             xmomc[k] = reduced_speed * hc;
    773           }
     845      if (fabs(u) > maximum_allowed_speed) {
     846        reduced_speed = maximum_allowed_speed * u/fabs(u);
     847        //printf("Speed (u) has been reduced from %.3f to %.3f\n",
     848        //  u, reduced_speed);
     849        xmomc[k] = reduced_speed * hc;
     850      }
    774851
    775852          v = ymomc[k]/hc;
    776           if (fabs(v) > maximum_allowed_speed) {
    777             reduced_speed = maximum_allowed_speed * v/fabs(v);
    778             //printf("Speed (v) has been reduced from %.3f to %.3f\n",
    779             //  v, reduced_speed);
    780             ymomc[k] = reduced_speed * hc;
    781           }
     853      if (fabs(v) > maximum_allowed_speed) {
     854        reduced_speed = maximum_allowed_speed * v/fabs(v);
     855        //printf("Speed (v) has been reduced from %.3f to %.3f\n",
     856        //  v, reduced_speed);
     857        ymomc[k] = reduced_speed * hc;
     858      }
    782859        }
    783860      }
     
    790867
    791868int _assign_wind_field_values(int N,
    792                               double* xmom_update,
    793                               double* ymom_update,
    794                               double* s_vec,
    795                               double* phi_vec,
    796                               double cw) {
     869                  double* xmom_update,
     870                  double* ymom_update,
     871                  double* s_vec,
     872                  double* phi_vec,
     873                  double cw) {
    797874
    798875  // Assign windfield values to momentum updates
     
    839916
    840917  if (!PyArg_ParseTuple(args, "OOOddOddd",
    841                         &normal, &ql, &qr, &zl, &zr, &edgeflux,
    842                         &epsilon, &g, &H0)) {
     918            &normal, &ql, &qr, &zl, &zr, &edgeflux,
     919            &epsilon, &g, &H0)) {
    843920    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: flux_function_central could not parse input arguments");
    844921    return NULL;
     
    847924 
    848925  _flux_function_central((double*) ql -> data,
    849                         (double*) qr -> data,
    850                         zl,
    851                          zr,                                           
    852                         ((double*) normal -> data)[0],
    853                          ((double*) normal -> data)[1],                 
    854                         epsilon, H0, g,
    855                         (double*) edgeflux -> data,
    856                         &max_speed);
     926            (double*) qr -> data,
     927            zl,
     928             zr,                       
     929            ((double*) normal -> data)[0],
     930             ((double*) normal -> data)[1],         
     931            epsilon, H0, g,
     932            (double*) edgeflux -> data,
     933            &max_speed);
    857934 
    858935  return Py_BuildValue("d", max_speed); 
     
    874951
    875952  if (!PyArg_ParseTuple(args, "dOOOOO",
    876                         &g, &h, &v, &x,
    877                         &xmom, &ymom)) {
     953            &g, &h, &v, &x,
     954            &xmom, &ymom)) {
    878955    //&epsilon)) {
    879956    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments");
     
    9401017
    9411018  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
    942                         &g, &eps, &w, &z, &uh, &vh, &eta,
    943                         &xmom, &ymom)) {
     1019            &g, &eps, &w, &z, &uh, &vh, &eta,
     1020            &xmom, &ymom)) {
    9441021    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction could not parse input arguments");
    9451022    return NULL;
     
    9571034  N = w -> dimensions[0];
    9581035  _manning_friction(g, eps, N,
    959                     (double*) w -> data,
    960                     (double*) z -> data,
    961                     (double*) uh -> data,
    962                     (double*) vh -> data,
    963                     (double*) eta -> data,
    964                     (double*) xmom -> data,
    965                     (double*) ymom -> data);
     1036            (double*) w -> data,
     1037            (double*) z -> data,
     1038            (double*) uh -> data,
     1039            (double*) vh -> data,
     1040            (double*) eta -> data,
     1041            (double*) xmom -> data,
     1042            (double*) ymom -> data);
    9661043
    9671044  return Py_BuildValue("");
     
    9811058
    9821059  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
    983                         &g, &eps, &w, &z, &uh, &vh, &eta,
    984                         &xmom, &ymom))
     1060            &g, &eps, &w, &z, &uh, &vh, &eta,
     1061            &xmom, &ymom))
    9851062    return NULL;
    9861063
     
    9961073  N = w -> dimensions[0];
    9971074  _manning_friction_explicit(g, eps, N,
    998                     (double*) w -> data,
    999                     (double*) z -> data,
    1000                     (double*) uh -> data,
    1001                     (double*) vh -> data,
    1002                     (double*) eta -> data,
    1003                     (double*) xmom -> data,
    1004                     (double*) ymom -> data);
     1075            (double*) w -> data,
     1076            (double*) z -> data,
     1077            (double*) uh -> data,
     1078            (double*) vh -> data,
     1079            (double*) eta -> data,
     1080            (double*) xmom -> data,
     1081            (double*) ymom -> data);
    10051082
    10061083  return Py_BuildValue("");
     
    10121089// Computational routine
    10131090int _extrapolate_second_order_sw(int number_of_elements,
    1014                                   double epsilon,
    1015                                   double minimum_allowed_height,
    1016                                   double beta_w,
    1017                                   double beta_w_dry,
    1018                                   double beta_uh,
    1019                                   double beta_uh_dry,
    1020                                   double beta_vh,
    1021                                   double beta_vh_dry,
    1022                                   long* surrogate_neighbours,
    1023                                   long* number_of_boundaries,
    1024                                   double* centroid_coordinates,
    1025                                   double* stage_centroid_values,
    1026                                   double* xmom_centroid_values,
    1027                                   double* ymom_centroid_values,
    1028                                   double* elevation_centroid_values,
    1029                                   double* vertex_coordinates,
    1030                                   double* stage_vertex_values,
    1031                                   double* xmom_vertex_values,
    1032                                   double* ymom_vertex_values,
    1033                                   double* elevation_vertex_values,
    1034                                   int optimise_dry_cells) {
    1035                                  
    1036                                  
     1091                                 double epsilon,
     1092                                 double minimum_allowed_height,
     1093                                 double beta_w,
     1094                                 double beta_w_dry,
     1095                                 double beta_uh,
     1096                                 double beta_uh_dry,
     1097                                 double beta_vh,
     1098                                 double beta_vh_dry,
     1099                                 long* surrogate_neighbours,
     1100                                 long* number_of_boundaries,
     1101                                 double* centroid_coordinates,
     1102                                 double* stage_centroid_values,
     1103                                 double* xmom_centroid_values,
     1104                                 double* ymom_centroid_values,
     1105                                 double* elevation_centroid_values,
     1106                                 double* vertex_coordinates,
     1107                                 double* stage_vertex_values,
     1108                                 double* xmom_vertex_values,
     1109                                 double* ymom_vertex_values,
     1110                                 double* elevation_vertex_values,
     1111                                 int optimise_dry_cells) {
     1112                 
     1113                 
    10371114
    10381115  // Local variables
    10391116  double a, b; // Gradient vector used to calculate vertex values from centroids
    1040   int k,k0,k1,k2,k3,k6,coord_index,i;
    1041   double x,y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2; // Vertices of the auxiliary triangle
    1042   double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;
     1117  int k, k0, k1, k2, k3, k6, coord_index, i;
     1118  double x, y, x0, y0, x1, y1, x2, y2, xv0, yv0, xv1, yv1, xv2, yv2; // Vertices of the auxiliary triangle
     1119  double dx1, dx2, dy1, dy2, dxv0, dxv1, dxv2, dyv0, dyv1, dyv2, dq0, dq1, dq2, area2, inv_area2;
    10431120  double dqv[3], qmin, qmax, hmin, hmax;
    10441121  double hc, h0, h1, h2, beta_tmp, hfactor;
    10451122
    1046        
    1047   for (k=0; k<number_of_elements; k++) {
     1123   
     1124  for (k = 0; k < number_of_elements; k++)
     1125  {
    10481126    k3=k*3;
    10491127    k6=k*6;
    10501128
    1051    
    1052     if (number_of_boundaries[k]==3){
     1129    if (number_of_boundaries[k]==3)
     1130    {
    10531131      // No neighbours, set gradient on the triangle to zero
    10541132     
     
    10651143      continue;
    10661144    }
    1067     else {
     1145    else
     1146    {
    10681147      // Triangle k has one or more neighbours.
    10691148      // Get centroid and vertex coordinates of the triangle
    10701149     
    10711150      // Get the vertex coordinates
    1072       xv0 = vertex_coordinates[k6];   yv0=vertex_coordinates[k6+1];
    1073       xv1 = vertex_coordinates[k6+2]; yv1=vertex_coordinates[k6+3];
    1074       xv2 = vertex_coordinates[k6+4]; yv2=vertex_coordinates[k6+5];
     1151      xv0 = vertex_coordinates[k6];   
     1152      yv0 = vertex_coordinates[k6+1];
     1153      xv1 = vertex_coordinates[k6+2];
     1154      yv1 = vertex_coordinates[k6+3];
     1155      xv2 = vertex_coordinates[k6+4];
     1156      yv2 = vertex_coordinates[k6+5];
    10751157     
    10761158      // Get the centroid coordinates
    1077       coord_index=2*k;
    1078       x=centroid_coordinates[coord_index];
    1079       y=centroid_coordinates[coord_index+1];
     1159      coord_index = 2*k;
     1160      x = centroid_coordinates[coord_index];
     1161      y = centroid_coordinates[coord_index+1];
    10801162     
    10811163      // Store x- and y- differentials for the vertices of
    10821164      // triangle k relative to the centroid
    1083       dxv0=xv0-x; dxv1=xv1-x; dxv2=xv2-x;
    1084       dyv0=yv0-y; dyv1=yv1-y; dyv2=yv2-y;
     1165      dxv0 = xv0 - x;
     1166      dxv1 = xv1 - x;
     1167      dxv2 = xv2 - x;
     1168      dyv0 = yv0 - y;
     1169      dyv1 = yv1 - y;
     1170      dyv2 = yv2 - y;
    10851171    }
    10861172
    10871173
    10881174           
    1089     if (number_of_boundaries[k]<=1){
    1090    
     1175    if (number_of_boundaries[k]<=1)
     1176    {
    10911177      //==============================================
    10921178      // Number of boundaries <= 1
     
    10991185      // from this centroid and its two neighbours
    11001186     
    1101       k0=surrogate_neighbours[k3];
    1102       k1=surrogate_neighbours[k3+1];
    1103       k2=surrogate_neighbours[k3+2];
     1187      k0 = surrogate_neighbours[k3];
     1188      k1 = surrogate_neighbours[k3 + 1];
     1189      k2 = surrogate_neighbours[k3 + 2];
    11041190     
    11051191      // Get the auxiliary triangle's vertex coordinates
    11061192      // (really the centroids of neighbouring triangles)
    1107       coord_index=2*k0;
    1108       x0=centroid_coordinates[coord_index];
    1109       y0=centroid_coordinates[coord_index+1];
    1110      
    1111       coord_index=2*k1;
    1112       x1=centroid_coordinates[coord_index];
    1113       y1=centroid_coordinates[coord_index+1];
    1114      
    1115       coord_index=2*k2;
    1116       x2=centroid_coordinates[coord_index];
    1117       y2=centroid_coordinates[coord_index+1];
     1193      coord_index = 2*k0;
     1194      x0 = centroid_coordinates[coord_index];
     1195      y0 = centroid_coordinates[coord_index+1];
     1196     
     1197      coord_index = 2*k1;
     1198      x1 = centroid_coordinates[coord_index];
     1199      y1 = centroid_coordinates[coord_index+1];
     1200     
     1201      coord_index = 2*k2;
     1202      x2 = centroid_coordinates[coord_index];
     1203      y2 = centroid_coordinates[coord_index+1];
    11181204     
    11191205      // Store x- and y- differentials for the vertices
    11201206      // of the auxiliary triangle
    1121       dx1=x1-x0; dx2=x2-x0;
    1122       dy1=y1-y0; dy2=y2-y0;
     1207      dx1 = x1 - x0;
     1208      dx2 = x2 - x0;
     1209      dy1 = y1 - y0;
     1210      dy2 = y2 - y0;
    11231211     
    11241212      // Calculate 2*area of the auxiliary triangle
     
    11281216      // If the mesh is 'weird' near the boundary,
    11291217      // the triangle might be flat or clockwise:
    1130       if (area2<=0) {
    1131         PyErr_SetString(PyExc_RuntimeError,
    1132                         "shallow_water_ext.c: negative triangle area encountered");
    1133         return -1;
     1218      if (area2 <= 0)
     1219      {
     1220        PyErr_SetString(PyExc_RuntimeError,
     1221                        "shallow_water_ext.c: negative triangle area encountered");
     1222   
     1223        return -1;
    11341224      } 
    11351225     
     
    11391229      h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
    11401230      h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
    1141       hmin = min(min(h0,min(h1,h2)),hc);
     1231      hmin = min(min(h0, min(h1, h2)), hc);
    11421232      //hfactor = hc/(hc + 1.0);
    11431233
    11441234      hfactor = 0.0;
    1145       if (hmin > 0.001 ) {
    1146           hfactor = (hmin-0.001)/(hmin+0.004);
    1147       }
    1148      
    1149       if (optimise_dry_cells) {     
    1150         // Check if linear reconstruction is necessary for triangle k
    1151         // This check will exclude dry cells.
    1152 
    1153         hmax = max(h0,max(h1,h2));     
    1154         if (hmax < epsilon) {
    1155           continue;
    1156         }
    1157       }
    1158 
    1159            
     1235      if (hmin > 0.001)
     1236      {
     1237        hfactor = (hmin - 0.001)/(hmin + 0.004);
     1238      }
     1239     
     1240      if (optimise_dry_cells)
     1241      {     
     1242        // Check if linear reconstruction is necessary for triangle k
     1243        // This check will exclude dry cells.
     1244
     1245        hmax = max(h0, max(h1, h2));     
     1246        if (hmax < epsilon)
     1247        {
     1248            continue;
     1249        }
     1250      }
     1251   
    11601252      //-----------------------------------
    11611253      // stage
     
    11641256      // Calculate the difference between vertex 0 of the auxiliary
    11651257      // triangle and the centroid of triangle k
    1166       dq0=stage_centroid_values[k0]-stage_centroid_values[k];
     1258      dq0 = stage_centroid_values[k0] - stage_centroid_values[k];
    11671259     
    11681260      // Calculate differentials between the vertices
    11691261      // of the auxiliary triangle (centroids of neighbouring triangles)
    1170       dq1=stage_centroid_values[k1]-stage_centroid_values[k0];
    1171       dq2=stage_centroid_values[k2]-stage_centroid_values[k0];
    1172      
     1262      dq1 = stage_centroid_values[k1] - stage_centroid_values[k0];
     1263      dq2 = stage_centroid_values[k2] - stage_centroid_values[k0];
     1264     
     1265          inv_area2 = 1.0/area2;
    11731266      // Calculate the gradient of stage on the auxiliary triangle
    11741267      a = dy2*dq1 - dy1*dq2;
    1175       a /= area2;
     1268      a *= inv_area2;
    11761269      b = dx1*dq2 - dx2*dq1;
    1177       b /= area2;
     1270      b *= inv_area2;
    11781271     
    11791272      // Calculate provisional jumps in stage from the centroid
    11801273      // of triangle k to its vertices, to be limited
    1181       dqv[0]=a*dxv0+b*dyv0;
    1182       dqv[1]=a*dxv1+b*dyv1;
    1183       dqv[2]=a*dxv2+b*dyv2;
     1274      dqv[0] = a*dxv0 + b*dyv0;
     1275      dqv[1] = a*dxv1 + b*dyv1;
     1276      dqv[2] = a*dxv2 + b*dyv2;
    11841277     
    11851278      // Now we want to find min and max of the centroid and the
    11861279      // vertices of the auxiliary triangle and compute jumps
    11871280      // from the centroid to the min and max
    1188       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
     1281      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
    11891282     
    11901283      // Playing with dry wet interface
     
    11931286      //if (hmin>minimum_allowed_height)
    11941287      beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
    1195        
     1288   
    11961289      //printf("min_alled_height = %f\n",minimum_allowed_height);
    11971290      //printf("hmin = %f\n",hmin);
     
    11991292      //printf("beta_tmp = %f\n",beta_tmp);
    12001293      // Limit the gradient
    1201       limit_gradient(dqv,qmin,qmax,beta_tmp);
    1202      
    1203       for (i=0;i<3;i++)
    1204         stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
     1294      limit_gradient(dqv, qmin, qmax, beta_tmp);
     1295     
     1296      //for (i=0;i<3;i++)
     1297      stage_vertex_values[k3+0] = stage_centroid_values[k] + dqv[0];
     1298          stage_vertex_values[k3+1] = stage_centroid_values[k] + dqv[1];
     1299          stage_vertex_values[k3+2] = stage_centroid_values[k] + dqv[2];
    12051300     
    12061301     
     
    12111306      // Calculate the difference between vertex 0 of the auxiliary
    12121307      // triangle and the centroid of triangle k     
    1213       dq0=xmom_centroid_values[k0]-xmom_centroid_values[k];
     1308      dq0 = xmom_centroid_values[k0] - xmom_centroid_values[k];
    12141309     
    12151310      // Calculate differentials between the vertices
    12161311      // of the auxiliary triangle
    1217       dq1=xmom_centroid_values[k1]-xmom_centroid_values[k0];
    1218       dq2=xmom_centroid_values[k2]-xmom_centroid_values[k0];
     1312      dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k0];
     1313      dq2 = xmom_centroid_values[k2] - xmom_centroid_values[k0];
    12191314     
    12201315      // Calculate the gradient of xmom on the auxiliary triangle
    12211316      a = dy2*dq1 - dy1*dq2;
    1222       a /= area2;
     1317      a *= inv_area2;
    12231318      b = dx1*dq2 - dx2*dq1;
    1224       b /= area2;
     1319      b *= inv_area2;
    12251320     
    12261321      // Calculate provisional jumps in stage from the centroid
    12271322      // of triangle k to its vertices, to be limited     
    1228       dqv[0]=a*dxv0+b*dyv0;
    1229       dqv[1]=a*dxv1+b*dyv1;
    1230       dqv[2]=a*dxv2+b*dyv2;
     1323      dqv[0] = a*dxv0+b*dyv0;
     1324      dqv[1] = a*dxv1+b*dyv1;
     1325      dqv[2] = a*dxv2+b*dyv2;
    12311326     
    12321327      // Now we want to find min and max of the centroid and the
    12331328      // vertices of the auxiliary triangle and compute jumps
    12341329      // from the centroid to the min and max
    1235       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
     1330      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
    12361331      //beta_tmp = beta_uh;
    12371332      //if (hmin<minimum_allowed_height)
     
    12401335
    12411336      // Limit the gradient
    1242       limit_gradient(dqv,qmin,qmax,beta_tmp);
    1243 
    1244       for (i=0;i<3;i++)
    1245         xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i];
    1246      
     1337      limit_gradient(dqv, qmin, qmax, beta_tmp);
     1338
     1339      for (i=0; i < 3; i++)
     1340      {
     1341        xmom_vertex_values[k3+i] = xmom_centroid_values[k] + dqv[i];
     1342      }
    12471343     
    12481344      //-----------------------------------
     
    12521348      // Calculate the difference between vertex 0 of the auxiliary
    12531349      // triangle and the centroid of triangle k
    1254       dq0=ymom_centroid_values[k0]-ymom_centroid_values[k];
     1350      dq0 = ymom_centroid_values[k0] - ymom_centroid_values[k];
    12551351     
    12561352      // Calculate differentials between the vertices
    12571353      // of the auxiliary triangle
    1258       dq1=ymom_centroid_values[k1]-ymom_centroid_values[k0];
    1259       dq2=ymom_centroid_values[k2]-ymom_centroid_values[k0];
     1354      dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k0];
     1355      dq2 = ymom_centroid_values[k2] - ymom_centroid_values[k0];
    12601356     
    12611357      // Calculate the gradient of xmom on the auxiliary triangle
    12621358      a = dy2*dq1 - dy1*dq2;
    1263       a /= area2;
     1359      a *= inv_area2;
    12641360      b = dx1*dq2 - dx2*dq1;
    1265       b /= area2;
     1361      b *= inv_area2;
    12661362     
    12671363      // Calculate provisional jumps in stage from the centroid
    12681364      // of triangle k to its vertices, to be limited
    1269       dqv[0]=a*dxv0+b*dyv0;
    1270       dqv[1]=a*dxv1+b*dyv1;
    1271       dqv[2]=a*dxv2+b*dyv2;
     1365      dqv[0] = a*dxv0 + b*dyv0;
     1366      dqv[1] = a*dxv1 + b*dyv1;
     1367      dqv[2] = a*dxv2 + b*dyv2;
    12721368     
    12731369      // Now we want to find min and max of the centroid and the
    12741370      // vertices of the auxiliary triangle and compute jumps
    12751371      // from the centroid to the min and max
    1276       find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax);
     1372      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
    12771373     
    12781374      //beta_tmp = beta_vh;
     
    12801376      //if (hmin<minimum_allowed_height)
    12811377      //beta_tmp = beta_vh_dry;
    1282       beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;       
     1378      beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;   
    12831379
    12841380      // Limit the gradient
    1285       limit_gradient(dqv,qmin,qmax,beta_tmp);
     1381      limit_gradient(dqv, qmin, qmax, beta_tmp);
    12861382     
    12871383      for (i=0;i<3;i++)
    1288         ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i];
    1289        
     1384      {
     1385        ymom_vertex_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
     1386      }
    12901387    } // End number_of_boundaries <=1
    1291     else{
     1388    else
     1389    {
    12921390
    12931391      //==============================================
     
    12981396     
    12991397      // Find the only internal neighbour (k1?)
    1300       for (k2=k3;k2<k3+3;k2++){
    1301         // Find internal neighbour of triangle k     
    1302         // k2 indexes the edges of triangle k   
    1303      
    1304         if (surrogate_neighbours[k2]!=k)
    1305           break;
    1306       }
    1307      
    1308       if ((k2==k3+3)) {
    1309         // If we didn't find an internal neighbour
    1310         PyErr_SetString(PyExc_RuntimeError,
    1311                         "shallow_water_ext.c: Internal neighbour not found");     
    1312         return -1;
    1313       }
    1314      
    1315       k1=surrogate_neighbours[k2];
     1398      for (k2 = k3; k2 < k3 + 3; k2++)
     1399      {
     1400      // Find internal neighbour of triangle k     
     1401      // k2 indexes the edges of triangle k   
     1402     
     1403          if (surrogate_neighbours[k2] != k)
     1404          {
     1405             break;
     1406          }
     1407      }
     1408     
     1409      if ((k2 == k3 + 3))
     1410      {
     1411        // If we didn't find an internal neighbour
     1412        PyErr_SetString(PyExc_RuntimeError,
     1413                        "shallow_water_ext.c: Internal neighbour not found");     
     1414        return -1;
     1415      }
     1416     
     1417      k1 = surrogate_neighbours[k2];
    13161418     
    13171419      // The coordinates of the triangle are already (x,y).
    13181420      // Get centroid of the neighbour (x1,y1)
    1319       coord_index=2*k1;
    1320       x1=centroid_coordinates[coord_index];
    1321       y1=centroid_coordinates[coord_index+1];
     1421      coord_index = 2*k1;
     1422      x1 = centroid_coordinates[coord_index];
     1423      y1 = centroid_coordinates[coord_index + 1];
    13221424     
    13231425      // Compute x- and y- distances between the centroid of
    13241426      // triangle k and that of its neighbour
    1325       dx1=x1-x; dy1=y1-y;
     1427      dx1 = x1 - x;
     1428      dy1 = y1 - y;
    13261429     
    13271430      // Set area2 as the square of the distance
    1328       area2=dx1*dx1+dy1*dy1;
     1431      area2 = dx1*dx1 + dy1*dy1;
    13291432     
    13301433      // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2)
     
    13321435      // respectively correspond to the x- and y- gradients
    13331436      // of the conserved quantities
    1334       dx2=1.0/area2;
    1335       dy2=dx2*dy1;
    1336       dx2*=dx1;
     1437      dx2 = 1.0/area2;
     1438      dy2 = dx2*dy1;
     1439      dx2 *= dx1;
    13371440     
    13381441     
     
    13421445
    13431446      // Compute differentials
    1344       dq1=stage_centroid_values[k1]-stage_centroid_values[k];
     1447      dq1 = stage_centroid_values[k1] - stage_centroid_values[k];
    13451448     
    13461449      // Calculate the gradient between the centroid of triangle k
    13471450      // and that of its neighbour
    1348       a=dq1*dx2;
    1349       b=dq1*dy2;
     1451      a = dq1*dx2;
     1452      b = dq1*dy2;
    13501453     
    13511454      // Calculate provisional vertex jumps, to be limited
    1352       dqv[0]=a*dxv0+b*dyv0;
    1353       dqv[1]=a*dxv1+b*dyv1;
    1354       dqv[2]=a*dxv2+b*dyv2;
     1455      dqv[0] = a*dxv0 + b*dyv0;
     1456      dqv[1] = a*dxv1 + b*dyv1;
     1457      dqv[2] = a*dxv2 + b*dyv2;
    13551458     
    13561459      // Now limit the jumps
    1357       if (dq1>=0.0){
    1358         qmin=0.0;
    1359         qmax=dq1;
    1360       }
    1361       else{
    1362         qmin=dq1;
    1363         qmax=0.0;
     1460      if (dq1>=0.0)
     1461      {
     1462        qmin=0.0;
     1463        qmax=dq1;
     1464      }
     1465      else
     1466      {
     1467        qmin = dq1;
     1468        qmax = 0.0;
    13641469      }
    13651470     
    13661471      // Limit the gradient
    1367       limit_gradient(dqv,qmin,qmax,beta_w);
    1368      
    1369       for (i=0;i<3;i++)
    1370         stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i];
    1371      
     1472      limit_gradient(dqv, qmin, qmax, beta_w);
     1473     
     1474      //for (i=0; i < 3; i++)
     1475      //{
     1476      stage_vertex_values[k3] = stage_centroid_values[k] + dqv[0];
     1477      stage_vertex_values[k3 + 1] = stage_centroid_values[k] + dqv[1];
     1478      stage_vertex_values[k3 + 2] = stage_centroid_values[k] + dqv[2];
     1479      //}
    13721480     
    13731481      //-----------------------------------
     
    13761484     
    13771485      // Compute differentials
    1378       dq1=xmom_centroid_values[k1]-xmom_centroid_values[k];
     1486      dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k];
    13791487     
    13801488      // Calculate the gradient between the centroid of triangle k
    13811489      // and that of its neighbour
    1382       a=dq1*dx2;
    1383       b=dq1*dy2;
     1490      a = dq1*dx2;
     1491      b = dq1*dy2;
    13841492     
    13851493      // Calculate provisional vertex jumps, to be limited
    1386       dqv[0]=a*dxv0+b*dyv0;
    1387       dqv[1]=a*dxv1+b*dyv1;
    1388       dqv[2]=a*dxv2+b*dyv2;
     1494      dqv[0] = a*dxv0+b*dyv0;
     1495      dqv[1] = a*dxv1+b*dyv1;
     1496      dqv[2] = a*dxv2+b*dyv2;
    13891497     
    13901498      // Now limit the jumps
    1391       if (dq1>=0.0){
    1392         qmin=0.0;
    1393         qmax=dq1;
    1394       }
    1395       else{
    1396         qmin=dq1;
    1397         qmax=0.0;
     1499      if (dq1 >= 0.0)
     1500      {
     1501        qmin = 0.0;
     1502        qmax = dq1;
     1503      }
     1504      else
     1505      {
     1506        qmin = dq1;
     1507        qmax = 0.0;
    13981508      }
    13991509     
    14001510      // Limit the gradient     
    1401       limit_gradient(dqv,qmin,qmax,beta_w);
    1402      
    1403       for (i=0;i<3;i++)
    1404         xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i];
    1405      
     1511      limit_gradient(dqv, qmin, qmax, beta_w);
     1512     
     1513      //for (i=0;i<3;i++)
     1514      //xmom_vertex_values[k3] = xmom_centroid_values[k] + dqv[0];
     1515      //xmom_vertex_values[k3 + 1] = xmom_centroid_values[k] + dqv[1];
     1516      //xmom_vertex_values[k3 + 2] = xmom_centroid_values[k] + dqv[2];
     1517     
     1518      for (i = 0; i < 3;i++)
     1519      {
     1520          xmom_vertex_values[k3 + i] = xmom_centroid_values[k] + dqv[i];
     1521      }
    14061522     
    14071523      //-----------------------------------
     
    14101526
    14111527      // Compute differentials
    1412       dq1=ymom_centroid_values[k1]-ymom_centroid_values[k];
     1528      dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k];
    14131529     
    14141530      // Calculate the gradient between the centroid of triangle k
    14151531      // and that of its neighbour
    1416       a=dq1*dx2;
    1417       b=dq1*dy2;
     1532      a = dq1*dx2;
     1533      b = dq1*dy2;
    14181534     
    14191535      // Calculate provisional vertex jumps, to be limited
    1420       dqv[0]=a*dxv0+b*dyv0;
    1421       dqv[1]=a*dxv1+b*dyv1;
    1422       dqv[2]=a*dxv2+b*dyv2;
     1536      dqv[0] = a*dxv0 + b*dyv0;
     1537      dqv[1] = a*dxv1 + b*dyv1;
     1538      dqv[2] = a*dxv2 + b*dyv2;
    14231539     
    14241540      // Now limit the jumps
    1425       if (dq1>=0.0){
    1426         qmin=0.0;
    1427         qmax=dq1;
    1428       }
    1429       else{
    1430         qmin=dq1;
    1431         qmax=0.0;
     1541      if (dq1>=0.0)
     1542      {
     1543        qmin = 0.0;
     1544        qmax = dq1;
     1545      }
     1546      else
     1547      {
     1548        qmin = dq1;
     1549        qmax = 0.0;
    14321550      }
    14331551     
    14341552      // Limit the gradient
    1435       limit_gradient(dqv,qmin,qmax,beta_w);
    1436      
    1437       for (i=0;i<3;i++)
    1438         ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i];
    1439        
     1553      limit_gradient(dqv, qmin, qmax, beta_w);
     1554     
     1555      //for (i=0;i<3;i++)
     1556      //ymom_vertex_values[k3] = ymom_centroid_values[k] + dqv[0];
     1557      //ymom_vertex_values[k3 + 1] = ymom_centroid_values[k] + dqv[1];
     1558      //ymom_vertex_values[k3 + 2] = ymom_centroid_values[k] + dqv[2];
     1559
     1560for (i=0;i<3;i++)
     1561        {
     1562ymom_vertex_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
     1563        }
     1564//ymom_vertex_values[k3] = ymom_centroid_values[k] + dqv[0];
     1565//ymom_vertex_values[k3 + 1] = ymom_centroid_values[k] + dqv[1];
     1566//ymom_vertex_values[k3 + 2] = ymom_centroid_values[k] + dqv[2];
    14401567    } // else [number_of_boundaries==2]
    14411568  } // for k=0 to number_of_elements-1
    14421569 
    14431570  return 0;
    1444 }                       
     1571}           
    14451572
    14461573
     
    14771604    Post conditions:
    14781605            The vertices of each triangle have values from a
    1479             limited linear reconstruction
    1480             based on centroid values
     1606        limited linear reconstruction
     1607        based on centroid values
    14811608
    14821609  */
     
    15051632  // Convert Python arguments to C
    15061633  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi",
    1507                         &domain,
    1508                         &surrogate_neighbours,
    1509                         &number_of_boundaries,
    1510                         &centroid_coordinates,
    1511                         &stage_centroid_values,
    1512                         &xmom_centroid_values,
    1513                         &ymom_centroid_values,
    1514                         &elevation_centroid_values,
    1515                         &vertex_coordinates,
    1516                         &stage_vertex_values,
    1517                         &xmom_vertex_values,
    1518                         &ymom_vertex_values,
    1519                         &elevation_vertex_values,
    1520                         &optimise_dry_cells)) {                 
    1521                        
     1634            &domain,
     1635            &surrogate_neighbours,
     1636            &number_of_boundaries,
     1637            &centroid_coordinates,
     1638            &stage_centroid_values,
     1639            &xmom_centroid_values,
     1640            &ymom_centroid_values,
     1641            &elevation_centroid_values,
     1642            &vertex_coordinates,
     1643            &stage_vertex_values,
     1644            &xmom_vertex_values,
     1645            &ymom_vertex_values,
     1646            &elevation_vertex_values,
     1647            &optimise_dry_cells)) {         
     1648           
    15221649    PyErr_SetString(PyExc_RuntimeError,
    1523                     "Input arguments to extrapolate_second_order_sw failed");
     1650            "Input arguments to extrapolate_second_order_sw failed");
    15241651    return NULL;
    15251652  }
     
    15571684  // Call underlying computational routine
    15581685  e = _extrapolate_second_order_sw(number_of_elements,
    1559                                    epsilon,
    1560                                    minimum_allowed_height,
    1561                                    beta_w,
    1562                                    beta_w_dry,
    1563                                    beta_uh,
    1564                                    beta_uh_dry,
    1565                                    beta_vh,
    1566                                    beta_vh_dry,
    1567                                    (long*) surrogate_neighbours -> data,
    1568                                    (long*) number_of_boundaries -> data,
    1569                                    (double*) centroid_coordinates -> data,
    1570                                    (double*) stage_centroid_values -> data,
    1571                                    (double*) xmom_centroid_values -> data,
    1572                                    (double*) ymom_centroid_values -> data,
    1573                                    (double*) elevation_centroid_values -> data,
    1574                                    (double*) vertex_coordinates -> data,
    1575                                    (double*) stage_vertex_values -> data,
    1576                                    (double*) xmom_vertex_values -> data,
    1577                                    (double*) ymom_vertex_values -> data,
    1578                                    (double*) elevation_vertex_values -> data,
    1579                                    optimise_dry_cells);
     1686                   epsilon,
     1687                   minimum_allowed_height,
     1688                   beta_w,
     1689                   beta_w_dry,
     1690                   beta_uh,
     1691                   beta_uh_dry,
     1692                   beta_vh,
     1693                   beta_vh_dry,
     1694                   (long*) surrogate_neighbours -> data,
     1695                   (long*) number_of_boundaries -> data,
     1696                   (double*) centroid_coordinates -> data,
     1697                   (double*) stage_centroid_values -> data,
     1698                   (double*) xmom_centroid_values -> data,
     1699                   (double*) ymom_centroid_values -> data,
     1700                   (double*) elevation_centroid_values -> data,
     1701                   (double*) vertex_coordinates -> data,
     1702                   (double*) stage_vertex_values -> data,
     1703                   (double*) xmom_vertex_values -> data,
     1704                   (double*) ymom_vertex_values -> data,
     1705                   (double*) elevation_vertex_values -> data,
     1706                   optimise_dry_cells);
    15801707  if (e == -1) {
    15811708    // Use error string set inside computational routine
    15821709    return NULL;
    1583   }                              
     1710  }                  
    15841711 
    15851712 
     
    16091736  // Convert Python arguments to C
    16101737  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
    1611                                    &Q, &Normal, &direction)) {
     1738                   &Q, &Normal, &direction)) {
    16121739    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: rotate could not parse input arguments");
    16131740    return NULL;
     
    16541781// Computational function for flux computation
    16551782double _compute_fluxes_central(int number_of_elements,
    1656                                double timestep,
    1657                                double epsilon,
    1658                                double H0,
    1659                                double g,
    1660                                long* neighbours,
    1661                                long* neighbour_edges,
    1662                                double* normals,
    1663                                double* edgelengths,
    1664                                double* radii,
    1665                                double* areas,
    1666                                long* tri_full_flag,
    1667                                double* stage_edge_values,
    1668                                double* xmom_edge_values,
    1669                                double* ymom_edge_values,
    1670                                double* bed_edge_values,
    1671                                double* stage_boundary_values,
    1672                                double* xmom_boundary_values,
    1673                                double* ymom_boundary_values,
    1674                                double* stage_explicit_update,
    1675                                double* xmom_explicit_update,
    1676                                double* ymom_explicit_update,
    1677                                long* already_computed_flux,
    1678                                double* max_speed_array,
    1679                                int optimise_dry_cells) {
    1680                                
     1783                               double timestep,
     1784                               double epsilon,
     1785                               double H0,
     1786                               double g,
     1787                               long* neighbours,
     1788                               long* neighbour_edges,
     1789                               double* normals,
     1790                               double* edgelengths,
     1791                               double* radii,
     1792                               double* areas,
     1793                               long* tri_full_flag,
     1794                               double* stage_edge_values,
     1795                               double* xmom_edge_values,
     1796                               double* ymom_edge_values,
     1797                               double* bed_edge_values,
     1798                               double* stage_boundary_values,
     1799                               double* xmom_boundary_values,
     1800                               double* ymom_boundary_values,
     1801                               double* stage_explicit_update,
     1802                               double* xmom_explicit_update,
     1803                               double* ymom_explicit_update,
     1804                               long* already_computed_flux,
     1805                               double* max_speed_array,
     1806                               int optimise_dry_cells)
     1807{                   
    16811808  // Local variables
    1682   double max_speed, length, area, zl, zr;
     1809  double max_speed, length, inv_area, zl, zr;
    16831810  int k, i, m, n;
    16841811  int ki, nm=0, ki2; // Index shorthands
     
    16871814  double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes
    16881815
    1689   static long call=1; // Static local variable flagging already computed flux
    1690                                
    1691        
     1816  static long call = 1; // Static local variable flagging already computed flux
     1817                   
    16921818  // Start computation
    16931819  call++; // Flag 'id' of flux calculation for this timestep
    1694  
     1820
    16951821  // Set explicit_update to zero for all conserved_quantities.
    16961822  // This assumes compute_fluxes called before forcing terms
    1697   for (k=0; k<number_of_elements; k++) {
    1698     stage_explicit_update[k]=0.0;
    1699     xmom_explicit_update[k]=0.0;
    1700     ymom_explicit_update[k]=0.0; 
    1701   }
    1702 
     1823  memset((char*) stage_explicit_update, 0, number_of_elements*sizeof(double));
     1824  memset((char*) xmom_explicit_update, 0, number_of_elements*sizeof(double));
     1825  memset((char*) ymom_explicit_update, 0, number_of_elements*sizeof(double));   
     1826 
    17031827  // For all triangles
    1704   for (k=0; k<number_of_elements; k++) {
    1705    
     1828  for (k = 0; k < number_of_elements; k++)
     1829  { 
    17061830    // Loop through neighbours and compute edge flux for each 
    1707     for (i=0; i<3; i++) {
    1708       ki = k*3+i; // Linear index (triangle k, edge i)
     1831    for (i = 0; i < 3; i++)
     1832    {
     1833      ki = k*3 + i; // Linear index to edge i of triangle k
    17091834     
    17101835      if (already_computed_flux[ki] == call)
     1836      {
    17111837        // We've already computed the flux across this edge
    1712         continue;
    1713        
    1714        
     1838        continue;
     1839      }
     1840   
     1841      // Get left hand side values from triangle k, edge i
    17151842      ql[0] = stage_edge_values[ki];
    17161843      ql[1] = xmom_edge_values[ki];
     
    17181845      zl = bed_edge_values[ki];
    17191846
    1720       // Quantities at neighbour on nearest face
     1847      // Get right hand side values either from neighbouring triangle
     1848      // or from boundary array (Quantities at neighbour on nearest face).
    17211849      n = neighbours[ki];
    1722       if (n < 0) {
     1850      if (n < 0)
     1851      {
    17231852        // Neighbour is a boundary condition
    1724         m = -n-1; // Convert negative flag to boundary index
    1725        
    1726         qr[0] = stage_boundary_values[m];
    1727         qr[1] = xmom_boundary_values[m];
    1728         qr[2] = ymom_boundary_values[m];
    1729         zr = zl; // Extend bed elevation to boundary
    1730       } else {
    1731         // Neighbour is a real element
    1732         m = neighbour_edges[ki];
    1733         nm = n*3+m; // Linear index (triangle n, edge m)
    1734        
    1735         qr[0] = stage_edge_values[nm];
    1736         qr[1] = xmom_edge_values[nm];
    1737         qr[2] = ymom_edge_values[nm];
    1738         zr = bed_edge_values[nm];
    1739       }
    1740      
    1741      
    1742       if (optimise_dry_cells) {     
    1743         // Check if flux calculation is necessary across this edge
    1744         // This check will exclude dry cells.
    1745         // This will also optimise cases where zl != zr as
    1746         // long as both are dry
    1747 
    1748         if ( fabs(ql[0] - zl) < epsilon &&
    1749              fabs(qr[0] - zr) < epsilon ) {
    1750           // Cell boundary is dry
    1751          
    1752           already_computed_flux[ki] = call; // #k Done 
    1753           if (n>=0)
    1754             already_computed_flux[nm] = call; // #n Done
    1755        
    1756           max_speed = 0.0;
    1757           continue;
    1758         }
    1759       }
    1760    
     1853        m = -n - 1; // Convert negative flag to boundary index
     1854   
     1855        qr[0] = stage_boundary_values[m];
     1856        qr[1] = xmom_boundary_values[m];
     1857        qr[2] = ymom_boundary_values[m];
     1858        zr = zl; // Extend bed elevation to boundary
     1859      }
     1860      else
     1861      {
     1862        // Neighbour is a real triangle
     1863        m = neighbour_edges[ki];
     1864        nm = n*3 + m; // Linear index (triangle n, edge m)
     1865       
     1866        qr[0] = stage_edge_values[nm];
     1867        qr[1] = xmom_edge_values[nm];
     1868        qr[2] = ymom_edge_values[nm];
     1869        zr = bed_edge_values[nm];
     1870      }
     1871     
     1872      // Now we have values for this edge - both from left and right side.
     1873     
     1874      if (optimise_dry_cells)
     1875      {     
     1876        // Check if flux calculation is necessary across this edge
     1877        // This check will exclude dry cells.
     1878        // This will also optimise cases where zl != zr as
     1879        // long as both are dry
     1880
     1881        if (fabs(ql[0] - zl) < epsilon &&
     1882            fabs(qr[0] - zr) < epsilon)
     1883        {
     1884          // Cell boundary is dry
     1885         
     1886          already_computed_flux[ki] = call; // #k Done 
     1887          if (n >= 0)
     1888          {
     1889            already_computed_flux[nm] = call; // #n Done
     1890          }
     1891       
     1892          max_speed = 0.0;
     1893          continue;
     1894        }
     1895      }
    17611896     
    17621897      // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
     
    17651900      // Edge flux computation (triangle k, edge i)
    17661901      _flux_function_central(ql, qr, zl, zr,
    1767                              normals[ki2], normals[ki2+1],
    1768                              epsilon, H0, g,
    1769                              edgeflux, &max_speed);
     1902                             normals[ki2], normals[ki2+1],
     1903                             epsilon, H0, g,
     1904                             edgeflux, &max_speed);
    17701905     
    17711906     
     
    17861921     
    17871922      // Update neighbour n with same flux but reversed sign
    1788       if (n>=0) {
    1789         stage_explicit_update[n] += edgeflux[0];
    1790         xmom_explicit_update[n] += edgeflux[1];
    1791         ymom_explicit_update[n] += edgeflux[2];
    1792        
    1793         already_computed_flux[nm] = call; // #n Done
    1794       }
    1795 
     1923      if (n >= 0)
     1924      {
     1925        stage_explicit_update[n] += edgeflux[0];
     1926        xmom_explicit_update[n] += edgeflux[1];
     1927        ymom_explicit_update[n] += edgeflux[2];
     1928       
     1929        already_computed_flux[nm] = call; // #n Done
     1930      }
    17961931
    17971932      // Update timestep based on edge i and possibly neighbour n
    1798       if (tri_full_flag[k] == 1) {
    1799         if (max_speed > epsilon) {
    1800        
    1801           // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
    1802          
    1803           // CFL for triangle k
    1804           timestep = min(timestep, radii[k]/max_speed);
    1805          
    1806           if (n>=0)
    1807             // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
    1808             timestep = min(timestep, radii[n]/max_speed);
    1809          
    1810           // Ted Rigby's suggested less conservative version
    1811           //if (n>=0) {             
    1812           //  timestep = min(timestep, (radii[k]+radii[n])/max_speed);
    1813           //} else {
    1814           //  timestep = min(timestep, radii[k]/max_speed);
    1815           // }
    1816         }
     1933      if (tri_full_flag[k] == 1)
     1934      {
     1935        if (max_speed > epsilon)
     1936        {
     1937          // Apply CFL condition for triangles joining this edge (triangle k and triangle n)
     1938         
     1939          // CFL for triangle k
     1940          timestep = min(timestep, radii[k]/max_speed);
     1941         
     1942          if (n >= 0)
     1943          {
     1944            // Apply CFL condition for neigbour n (which is on the ith edge of triangle k)
     1945            timestep = min(timestep, radii[n]/max_speed);
     1946          }
     1947
     1948          // Ted Rigby's suggested less conservative version
     1949          //if (n>=0) {         
     1950          //  timestep = min(timestep, (radii[k]+radii[n])/max_speed);
     1951          //} else {
     1952          //  timestep = min(timestep, radii[k]/max_speed);
     1953          // }
     1954        }
    18171955      }
    18181956     
     
    18221960    // Normalise triangle k by area and store for when all conserved
    18231961    // quantities get updated
    1824     area = areas[k];
    1825     stage_explicit_update[k] /= area;
    1826     xmom_explicit_update[k] /= area;
    1827     ymom_explicit_update[k] /= area;
     1962    inv_area = 1.0/areas[k];
     1963    stage_explicit_update[k] *= inv_area;
     1964    xmom_explicit_update[k] *= inv_area;
     1965    ymom_explicit_update[k] *= inv_area;
    18281966   
    18291967   
     
    18321970   
    18331971  } // End triangle k
    1834        
    1835                                
    1836                                
     1972             
    18371973  return timestep;
    18381974}
     
    18692005
    18702006    PyObject
    1871         *domain,
    1872         *stage,
    1873         *xmom,
    1874         *ymom,
    1875         *bed;
     2007    *domain,
     2008    *stage,
     2009    *xmom,
     2010    *ymom,
     2011    *bed;
    18762012
    18772013    PyArrayObject
    1878         *neighbours,
    1879         *neighbour_edges,
    1880         *normals,
    1881         *edgelengths,
    1882         *radii,
    1883         *areas,
    1884         *tri_full_flag,
    1885         *stage_edge_values,
    1886         *xmom_edge_values,
    1887         *ymom_edge_values,
    1888         *bed_edge_values,
    1889         *stage_boundary_values,
    1890         *xmom_boundary_values,
    1891         *ymom_boundary_values,
    1892         *stage_explicit_update,
    1893         *xmom_explicit_update,
    1894         *ymom_explicit_update,
    1895         *already_computed_flux, //Tracks whether the flux across an edge has already been computed
    1896         *max_speed_array; //Keeps track of max speeds for each triangle
     2014    *neighbours,
     2015    *neighbour_edges,
     2016    *normals,
     2017    *edgelengths,
     2018    *radii,
     2019    *areas,
     2020    *tri_full_flag,
     2021    *stage_edge_values,
     2022    *xmom_edge_values,
     2023    *ymom_edge_values,
     2024    *bed_edge_values,
     2025    *stage_boundary_values,
     2026    *xmom_boundary_values,
     2027    *ymom_boundary_values,
     2028    *stage_explicit_update,
     2029    *xmom_explicit_update,
     2030    *ymom_explicit_update,
     2031    *already_computed_flux, //Tracks whether the flux across an edge has already been computed
     2032    *max_speed_array; //Keeps track of max speeds for each triangle
    18972033
    18982034   
     
    19022038    // Convert Python arguments to C
    19032039    if (!PyArg_ParseTuple(args, "dOOOO", &timestep, &domain, &stage, &xmom, &ymom, &bed )) {
    1904         PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
    1905         return NULL;
     2040    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
     2041    return NULL;
    19062042    }
    19072043
     
    19402076  // the explicit update arrays
    19412077  timestep = _compute_fluxes_central(number_of_elements,
    1942                                      timestep,
    1943                                      epsilon,
    1944                                      H0,
    1945                                      g,
    1946                                      (long*) neighbours -> data,
    1947                                      (long*) neighbour_edges -> data,
    1948                                      (double*) normals -> data,
    1949                                      (double*) edgelengths -> data,
    1950                                      (double*) radii -> data,
    1951                                      (double*) areas -> data,
    1952                                      (long*) tri_full_flag -> data,
    1953                                      (double*) stage_edge_values -> data,
    1954                                      (double*) xmom_edge_values -> data,
    1955                                      (double*) ymom_edge_values -> data,
    1956                                      (double*) bed_edge_values -> data,
    1957                                      (double*) stage_boundary_values -> data,
    1958                                      (double*) xmom_boundary_values -> data,
    1959                                      (double*) ymom_boundary_values -> data,
    1960                                      (double*) stage_explicit_update -> data,
    1961                                      (double*) xmom_explicit_update -> data,
    1962                                      (double*) ymom_explicit_update -> data,
    1963                                      (long*) already_computed_flux -> data,
    1964                                      (double*) max_speed_array -> data,
    1965                                      optimise_dry_cells);
     2078                     timestep,
     2079                     epsilon,
     2080                     H0,
     2081                     g,
     2082                     (long*) neighbours -> data,
     2083                     (long*) neighbour_edges -> data,
     2084                     (double*) normals -> data,
     2085                     (double*) edgelengths -> data,
     2086                     (double*) radii -> data,
     2087                     (double*) areas -> data,
     2088                     (long*) tri_full_flag -> data,
     2089                     (double*) stage_edge_values -> data,
     2090                     (double*) xmom_edge_values -> data,
     2091                     (double*) ymom_edge_values -> data,
     2092                     (double*) bed_edge_values -> data,
     2093                     (double*) stage_boundary_values -> data,
     2094                     (double*) xmom_boundary_values -> data,
     2095                     (double*) ymom_boundary_values -> data,
     2096                     (double*) stage_explicit_update -> data,
     2097                     (double*) xmom_explicit_update -> data,
     2098                     (double*) ymom_explicit_update -> data,
     2099                     (long*) already_computed_flux -> data,
     2100                     (double*) max_speed_array -> data,
     2101                     optimise_dry_cells);
    19662102
    19672103  Py_DECREF(neighbours);
     
    20132149    domain.timestep = compute_fluxes(timestep,
    20142150                                     domain.epsilon,
    2015                                      domain.H0,
     2151                     domain.H0,
    20162152                                     domain.g,
    20172153                                     domain.neighbours,
     
    20332169                                     Ymom.explicit_update,
    20342170                                     already_computed_flux,
    2035                                      optimise_dry_cells)                                     
     2171                     optimise_dry_cells)                     
    20362172
    20372173
     
    20662202  // Convert Python arguments to C
    20672203  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi",
    2068                         &timestep,
    2069                         &epsilon,
    2070                         &H0,
    2071                         &g,
    2072                         &neighbours,
    2073                         &neighbour_edges,
    2074                         &normals,
    2075                         &edgelengths, &radii, &areas,
    2076                         &tri_full_flag,
    2077                         &stage_edge_values,
    2078                         &xmom_edge_values,
    2079                         &ymom_edge_values,
    2080                         &bed_edge_values,
    2081                         &stage_boundary_values,
    2082                         &xmom_boundary_values,
    2083                         &ymom_boundary_values,
    2084                         &stage_explicit_update,
    2085                         &xmom_explicit_update,
    2086                         &ymom_explicit_update,
    2087                         &already_computed_flux,
    2088                         &max_speed_array,
    2089                         &optimise_dry_cells)) {
     2204            &timestep,
     2205            &epsilon,
     2206            &H0,
     2207            &g,
     2208            &neighbours,
     2209            &neighbour_edges,
     2210            &normals,
     2211            &edgelengths, &radii, &areas,
     2212            &tri_full_flag,
     2213            &stage_edge_values,
     2214            &xmom_edge_values,
     2215            &ymom_edge_values,
     2216            &bed_edge_values,
     2217            &stage_boundary_values,
     2218            &xmom_boundary_values,
     2219            &ymom_boundary_values,
     2220            &stage_explicit_update,
     2221            &xmom_explicit_update,
     2222            &ymom_explicit_update,
     2223            &already_computed_flux,
     2224            &max_speed_array,
     2225            &optimise_dry_cells)) {
    20902226    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
    20912227    return NULL;
     
    21182254  // the explicit update arrays
    21192255  timestep = _compute_fluxes_central(number_of_elements,
    2120                                      timestep,
    2121                                      epsilon,
    2122                                      H0,
    2123                                      g,
    2124                                      (long*) neighbours -> data,
    2125                                      (long*) neighbour_edges -> data,
    2126                                      (double*) normals -> data,
    2127                                      (double*) edgelengths -> data,
    2128                                      (double*) radii -> data,
    2129                                      (double*) areas -> data,
    2130                                      (long*) tri_full_flag -> data,
    2131                                      (double*) stage_edge_values -> data,
    2132                                      (double*) xmom_edge_values -> data,
    2133                                      (double*) ymom_edge_values -> data,
    2134                                      (double*) bed_edge_values -> data,
    2135                                      (double*) stage_boundary_values -> data,
    2136                                      (double*) xmom_boundary_values -> data,
    2137                                      (double*) ymom_boundary_values -> data,
    2138                                      (double*) stage_explicit_update -> data,
    2139                                      (double*) xmom_explicit_update -> data,
    2140                                      (double*) ymom_explicit_update -> data,
    2141                                      (long*) already_computed_flux -> data,
    2142                                      (double*) max_speed_array -> data,
    2143                                      optimise_dry_cells);
     2256                     timestep,
     2257                     epsilon,
     2258                     H0,
     2259                     g,
     2260                     (long*) neighbours -> data,
     2261                     (long*) neighbour_edges -> data,
     2262                     (double*) normals -> data,
     2263                     (double*) edgelengths -> data,
     2264                     (double*) radii -> data,
     2265                     (double*) areas -> data,
     2266                     (long*) tri_full_flag -> data,
     2267                     (double*) stage_edge_values -> data,
     2268                     (double*) xmom_edge_values -> data,
     2269                     (double*) ymom_edge_values -> data,
     2270                     (double*) bed_edge_values -> data,
     2271                     (double*) stage_boundary_values -> data,
     2272                     (double*) xmom_boundary_values -> data,
     2273                     (double*) ymom_boundary_values -> data,
     2274                     (double*) stage_explicit_update -> data,
     2275                     (double*) xmom_explicit_update -> data,
     2276                     (double*) ymom_explicit_update -> data,
     2277                     (long*) already_computed_flux -> data,
     2278                     (double*) max_speed_array -> data,
     2279                     optimise_dry_cells);
    21442280 
    21452281  // Return updated flux timestep
     
    22282364  // Convert Python arguments to C
    22292365  if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOO",
    2230                         &timestep,
    2231                         &epsilon,
    2232                         &H0,
    2233                         &g,
    2234                         &neighbours,
    2235                         &neighbour_edges,
    2236                         &normals,
    2237                         &edgelengths, &radii, &areas,
    2238                         &tri_full_flag,
    2239                         &stage_edge_values,
    2240                         &xmom_edge_values,
    2241                         &ymom_edge_values,
    2242                         &bed_edge_values,
    2243                         &stage_boundary_values,
    2244                         &xmom_boundary_values,
    2245                         &ymom_boundary_values,
    2246                         &stage_explicit_update,
    2247                         &xmom_explicit_update,
    2248                         &ymom_explicit_update,
    2249                         &already_computed_flux)) {
     2366            &timestep,
     2367            &epsilon,
     2368            &H0,
     2369            &g,
     2370            &neighbours,
     2371            &neighbour_edges,
     2372            &normals,
     2373            &edgelengths, &radii, &areas,
     2374            &tri_full_flag,
     2375            &stage_edge_values,
     2376            &xmom_edge_values,
     2377            &ymom_edge_values,
     2378            &bed_edge_values,
     2379            &stage_boundary_values,
     2380            &xmom_boundary_values,
     2381            &ymom_boundary_values,
     2382            &stage_explicit_update,
     2383            &xmom_explicit_update,
     2384            &ymom_explicit_update,
     2385            &already_computed_flux)) {
    22502386    PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
    22512387    return NULL;
     
    22652401      ki = k*3+i;
    22662402      if (((long *) already_computed_flux->data)[ki]==call)//we've already computed the flux across this edge
    2267         continue;
     2403    continue;
    22682404      ql[0] = ((double *) stage_edge_values -> data)[ki];
    22692405      ql[1] = ((double *) xmom_edge_values -> data)[ki];
     
    22742410      n = ((long *) neighbours -> data)[ki];
    22752411      if (n < 0) {
    2276         m = -n-1; //Convert negative flag to index
    2277         qr[0] = ((double *) stage_boundary_values -> data)[m];
    2278         qr[1] = ((double *) xmom_boundary_values -> data)[m];
    2279         qr[2] = ((double *) ymom_boundary_values -> data)[m];
    2280         zr = zl; //Extend bed elevation to boundary
     2412    m = -n-1; //Convert negative flag to index
     2413    qr[0] = ((double *) stage_boundary_values -> data)[m];
     2414    qr[1] = ((double *) xmom_boundary_values -> data)[m];
     2415    qr[2] = ((double *) ymom_boundary_values -> data)[m];
     2416    zr = zl; //Extend bed elevation to boundary
    22812417      } else {
    2282         m = ((long *) neighbour_edges -> data)[ki];
    2283         nm = n*3+m;
    2284         qr[0] = ((double *) stage_edge_values -> data)[nm];
    2285         qr[1] = ((double *) xmom_edge_values -> data)[nm];
    2286         qr[2] = ((double *) ymom_edge_values -> data)[nm];
    2287         zr =    ((double *) bed_edge_values -> data)[nm];
     2418    m = ((long *) neighbour_edges -> data)[ki];
     2419    nm = n*3+m;
     2420    qr[0] = ((double *) stage_edge_values -> data)[nm];
     2421    qr[1] = ((double *) xmom_edge_values -> data)[nm];
     2422    qr[2] = ((double *) ymom_edge_values -> data)[nm];
     2423    zr =    ((double *) bed_edge_values -> data)[nm];
    22882424      }
    22892425      // Outward pointing normal vector
     
    22942430      //Edge flux computation
    22952431      flux_function_kinetic(ql, qr, zl, zr,
    2296                     normal[0], normal[1],
    2297                     epsilon, H0, g,
    2298                     edgeflux, &max_speed);
     2432            normal[0], normal[1],
     2433            epsilon, H0, g,
     2434            edgeflux, &max_speed);
    22992435      //update triangle k
    23002436      ((long *) already_computed_flux->data)[ki]=call;
     
    23042440      //update the neighbour n
    23052441      if (n>=0){
    2306         ((long *) already_computed_flux->data)[nm]=call;
    2307         ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
    2308         ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
    2309         ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
     2442    ((long *) already_computed_flux->data)[nm]=call;
     2443    ((double *) stage_explicit_update -> data)[n] += edgeflux[0]*((double *) edgelengths -> data)[nm];
     2444    ((double *) xmom_explicit_update -> data)[n] += edgeflux[1]*((double *) edgelengths -> data)[nm];
     2445    ((double *) ymom_explicit_update -> data)[n] += edgeflux[2]*((double *) edgelengths -> data)[nm];
    23102446      }
    23112447      ///for (j=0; j<3; j++) {
    2312         ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
    2313         ///}
    2314         //Update timestep
    2315         //timestep = min(timestep, domain.radii[k]/max_speed)
    2316         //FIXME: SR Add parameter for CFL condition
     2448    ///flux[j] -= edgeflux[j]*((double *) edgelengths -> data)[ki];
     2449    ///}
     2450    //Update timestep
     2451    //timestep = min(timestep, domain.radii[k]/max_speed)
     2452    //FIXME: SR Add parameter for CFL condition
    23172453    if ( ((long *) tri_full_flag -> data)[k] == 1) {
    2318             if (max_speed > epsilon) {
    2319                 timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
    2320                 //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
    2321                 if (n>=0)
    2322                     timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
    2323             }
     2454        if (max_speed > epsilon) {
     2455            timestep = min(timestep, ((double *) radii -> data)[k]/max_speed);
     2456            //maxspeed in flux_function is calculated as max(|u+a|,|u-a|)
     2457            if (n>=0)
     2458                timestep = min(timestep, ((double *) radii -> data)[n]/max_speed);
     2459        }
    23242460    }
    23252461    } // end for i
     
    23502486  // Convert Python arguments to C
    23512487  if (!PyArg_ParseTuple(args, "dddOOOO",
    2352                         &minimum_allowed_height,
    2353                         &maximum_allowed_speed,
    2354                         &epsilon,
    2355                         &wc, &zc, &xmomc, &ymomc)) {
     2488            &minimum_allowed_height,
     2489            &maximum_allowed_speed,
     2490            &epsilon,
     2491            &wc, &zc, &xmomc, &ymomc)) {
    23562492    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: protect could not parse input arguments");
    23572493    return NULL;
     
    23612497
    23622498  _protect(N,
    2363            minimum_allowed_height,
    2364            maximum_allowed_speed,
    2365            epsilon,
    2366            (double*) wc -> data,
    2367            (double*) zc -> data,
    2368            (double*) xmomc -> data,
    2369            (double*) ymomc -> data);
     2499       minimum_allowed_height,
     2500       maximum_allowed_speed,
     2501       epsilon,
     2502       (double*) wc -> data,
     2503       (double*) zc -> data,
     2504       (double*) xmomc -> data,
     2505       (double*) ymomc -> data);
    23702506
    23712507  return Py_BuildValue("");
     
    24092545  // Convert Python arguments to C
    24102546  if (!PyArg_ParseTuple(args, "OOOOOOOOOO",
    2411                         &domain,
    2412                         &wc, &zc,
    2413                         &wv, &zv, &hvbar,
    2414                         &xmomc, &ymomc, &xmomv, &ymomv)) {
     2547            &domain,
     2548            &wc, &zc,
     2549            &wv, &zv, &hvbar,
     2550            &xmomc, &ymomc, &xmomv, &ymomv)) {
    24152551    PyErr_SetString(PyExc_RuntimeError,
    2416                     "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
     2552            "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments");
    24172553    return NULL;
    24182554  } 
    2419          
    2420          
     2555     
     2556     
    24212557  // FIXME (Ole): I tested this without GetAttrString and got time down
    24222558  // marginally from 4.0s to 3.8s. Consider passing everything in
     
    24632599  N = wc -> dimensions[0];
    24642600  _balance_deep_and_shallow(N,
    2465                             (double*) wc -> data,
    2466                             (double*) zc -> data,
    2467                             (double*) wv -> data,
    2468                             (double*) zv -> data,
    2469                             (double*) hvbar -> data,
    2470                             (double*) xmomc -> data,
    2471                             (double*) ymomc -> data,
    2472                             (double*) xmomv -> data,
    2473                             (double*) ymomv -> data,
    2474                             H0,
    2475                             (int) tight_slope_limiters,
    2476                             (int) use_centroid_velocities,                         
    2477                             alpha_balance);
     2601                (double*) wc -> data,
     2602                (double*) zc -> data,
     2603                (double*) wv -> data,
     2604                (double*) zv -> data,
     2605                (double*) hvbar -> data,
     2606                (double*) xmomc -> data,
     2607                (double*) ymomc -> data,
     2608                (double*) xmomv -> data,
     2609                (double*) ymomv -> data,
     2610                H0,
     2611                (int) tight_slope_limiters,
     2612                (int) use_centroid_velocities,             
     2613                alpha_balance);
    24782614
    24792615
     
    25032639  // Convert Python arguments to C
    25042640  if (!PyArg_ParseTuple(args, "OOOOd",
    2505                         &xmom_update,
    2506                         &ymom_update,
    2507                         &s_vec, &phi_vec,
    2508                         &cw)) {
     2641            &xmom_update,
     2642            &ymom_update,
     2643            &s_vec, &phi_vec,
     2644            &cw)) {
    25092645    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: assign_windfield_values could not parse input arguments");
    25102646    return NULL;
    25112647  } 
    2512                        
     2648           
    25132649
    25142650  N = xmom_update -> dimensions[0];
    25152651
    25162652  _assign_wind_field_values(N,
    2517            (double*) xmom_update -> data,
    2518            (double*) ymom_update -> data,
    2519            (double*) s_vec -> data,
    2520            (double*) phi_vec -> data,
    2521            cw);
     2653       (double*) xmom_update -> data,
     2654       (double*) ymom_update -> data,
     2655       (double*) s_vec -> data,
     2656       (double*) phi_vec -> data,
     2657       cw);
    25222658
    25232659  return Py_BuildValue("");
  • branches/numpy/anuga/shallow_water/test_data_manager.py

    r6689 r6902  
    1515from struct import pack, unpack
    1616from sets import ImmutableSet
     17import shutil
    1718
    1819from anuga.shallow_water import *
     
    4748       
    4849        self.verbose = Test_Data_Manager.verbose
    49         #Create basic mesh
     50        # Create basic mesh
    5051        points, vertices, boundary = rectangular(2, 2)
    5152
    52         #Create shallow water domain
     53        # Create shallow water domain
    5354        domain = Domain(points, vertices, boundary)
    5455        domain.default_order = 2
    5556
    56         #Set some field values
     57        # Set some field values
    5758        domain.set_quantity('elevation', lambda x,y: -x)
    5859        domain.set_quantity('friction', 0.03)
     
    252253
    253254        # Get the variables
    254         #range = fid.variables['stage_range'][:]
     255        range = fid.variables['stage_range'][:]
    255256        #print range
    256         #assert num.allclose(range,[-0.93519, 0.15]) or\
    257         #       num.allclose(range,[-0.9352743, 0.15]) or\
    258         #       num.allclose(range,[-0.93522203, 0.15000001]) # Old slope limiters
    259        
    260         #range = fid.variables['xmomentum_range'][:]
     257        assert num.allclose(range,[-0.93519, 0.15]) or\
     258               num.allclose(range,[-0.9352743, 0.15]) or\
     259               num.allclose(range,[-0.93522203, 0.15000001]) # Old slope limiters
     260       
     261        range = fid.variables['xmomentum_range'][:]
    261262        ##print range
    262         #assert num.allclose(range,[0,0.4695096]) or \
    263         #       num.allclose(range,[0,0.47790655]) or\
    264         #       num.allclose(range,[0,0.46941957]) or\
    265         #       num.allclose(range,[0,0.47769409])
    266 
    267        
    268         #range = fid.variables['ymomentum_range'][:]
     263        assert num.allclose(range,[0,0.4695096]) or \
     264               num.allclose(range,[0,0.47790655]) or\
     265               num.allclose(range,[0,0.46941957]) or\
     266               num.allclose(range,[0,0.47769409])
     267
     268       
     269        range = fid.variables['ymomentum_range'][:]
    269270        ##print range
    270         #assert num.allclose(range,[0,0.02174380]) or\
    271         #       num.allclose(range,[0,0.02174439]) or\
    272         #       num.allclose(range,[0,0.02283983]) or\
    273         #       num.allclose(range,[0,0.0217342]) or\
    274         #       num.allclose(range,[0,0.0227564]) # Old slope limiters
     271        assert num.allclose(range,[0,0.02174380]) or\
     272               num.allclose(range,[0,0.02174439]) or\
     273               num.allclose(range,[0,0.02283983]) or\
     274               num.allclose(range,[0,0.0217342]) or\
     275               num.allclose(range,[0,0.0227564]) # Old slope limiters
    275276       
    276277        fid.close()
     
    12321233        from Scientific.IO.NetCDF import NetCDFFile
    12331234
    1234         #Setup
     1235        # Setup
    12351236        self.domain.set_name('datatest')
    12361237
     
    12451246        self.domain.set_quantity('stage', 1.0)
    12461247
    1247         self.domain.geo_reference = Geo_reference(56,308500,6189000)
     1248        self.domain.geo_reference = Geo_reference(56, 308500, 6189000)
    12481249
    12491250        sww = get_dataobject(self.domain)
     
    55635564            quantities_init[0].append(this_ha) # HA
    55645565            quantities_init[1].append(this_ua) # UA
    5565             quantities_init[2].append(this_va) # VA
     5566            quantities_init[2].append(this_va) # VA 
    55665567               
    55675568        file_handle, base_name = tempfile.mkstemp("")
     
    61046105                        #print 'writing', time, point_i, q_time[time, point_i]
    61056106                        f.write(pack('f', q_time[time, point_i]))
    6106 
    61076107            f.close()
    61086108
     
    62156215
    62166216        msg='Incorrect gauge depths returned'
    6217         assert num.allclose(elevation,-depth),msg
     6217        assert num.allclose(elevation, -depth),msg
    62186218        msg='incorrect gauge height time series returned'
    6219         assert num.allclose(stage,ha)
     6219        assert num.allclose(stage, ha)
    62206220        msg='incorrect gauge ua time series returned'
    6221         assert num.allclose(xvelocity,ua)
     6221        assert num.allclose(xvelocity, ua)
    62226222        msg='incorrect gauge va time series returned'
    62236223        assert num.allclose(yvelocity, -va) # South is positive in MUX
     
    63976397                if j == 2: assert num.allclose(data[i][:parameters_index], -va0[permutation[i], :])
    63986398       
     6399        self.delete_mux(filesI)
    63996400
    64006401
     
    64076408        wrong values Win32
    64086409
    6409         This test does not pass on Windows but test_read_mux_platform_problem1 does
     6410        This test does not pass on Windows but test_read_mux_platform_problem1
     6411        does
    64106412        """
    64116413       
     
    64226424        times_ref = num.arange(0, time_step_count*time_step, time_step)
    64236425       
    6424         lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115), (-21.,115.), (-22., 117.)]
     6426        lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115),
     6427                           (-21.,115.), (-22., 117.)]
    64256428        n = len(lat_long_points)
    64266429       
     
    64666469             for j in range(0, first_tstep[i]-1) + range(last_tstep[i], time_step_count):
    64676470                 # For timesteps before and after recording range
    6468                  ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0                                  
     6471                 ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0
    64696472
    64706473
     
    64746477        #print 'va', va1[0,:]
    64756478       
    6476         # Write second mux file to be combined by urs2sts                                            
     6479        # Write second mux file to be combined by urs2sts
    64776480        base_nameII, filesII = self.write_mux2(lat_long_points,
    64786481                                               time_step_count, time_step,
     
    66056608
    66066609            f.close()
    6607 
    6608 
    6609 
    6610 
    6611 
    6612 
    6613        
    66146610                                               
    66156611        # Create ordering file
    66166612        permutation = ensure_numeric([4,0,2])
    66176613
    6618         _, ordering_filename = tempfile.mkstemp('')
    6619         order_fid = open(ordering_filename, 'w') 
    6620         order_fid.write('index, longitude, latitude\n')
    6621         for index in permutation:
    6622             order_fid.write('%d, %f, %f\n' %(index,
    6623                                              lat_long_points[index][1],
    6624                                              lat_long_points[index][0]))
    6625         order_fid.close()
    6626        
    6627        
    6628 
     6614       #  _, ordering_filename = tempfile.mkstemp('')
     6615#         order_fid = open(ordering_filename, 'w') 
     6616#         order_fid.write('index, longitude, latitude\n')
     6617#         for index in permutation:
     6618#             order_fid.write('%d, %f, %f\n' %(index,
     6619#                                              lat_long_points[index][1],
     6620#                                              lat_long_points[index][0]))
     6621#         order_fid.close()
     6622       
    66296623        # -------------------------------------
    66306624        # Now read files back and check values
     
    66386632        for j, file in enumerate(filesII):
    66396633            # Read stage, u, v enumerated as j
    6640 
    6641            
    66426634            #print 'Reading', j, file
    66436635            data = read_mux2(1, [file], weights, file_params, permutation, verbose)
     
    66456637            #print 'Data received by Python'
    66466638            #print data[1][8]
    6647 
    6648            
    66496639            number_of_selected_stations = data.shape[0]
    66506640
     
    66596649                #print i, parameters_index
    66606650                #print quantity[i][:]
    6661 
    6662                
    66636651                if j == 0: assert num.allclose(data[i][:parameters_index], ha1[permutation[i], :])
    66646652                if j == 1: assert num.allclose(data[i][:parameters_index], ua1[permutation[i], :])
     
    66696657                    #print j, i
    66706658                    #print 'Input'
    6671                     #print 'u', ua1[permutation[i], 8]                                        
     6659                    #print 'u', ua1[permutation[i], 8]       
    66726660                    #print 'v', va1[permutation[i], 8]
    66736661               
    66746662                    #print 'Output'
    6675                     #print 'v ', data[i][:parameters_index][8]                   
     6663                    #print 'v ', data[i][:parameters_index][8] 
     6664
     6665                    # South is positive in MUX
     6666                    #print "data[i][:parameters_index]", data[i][:parameters_index]
     6667                    #print "-va1[permutation[i], :]", -va1[permutation[i], :]
     6668                    assert num.allclose(data[i][:parameters_index], -va1[permutation[i], :])
     6669       
     6670        self.delete_mux(filesII)
     6671           
     6672    def test_read_mux_platform_problem3(self):
     6673       
     6674        # This is to test a situation where read_mux returned
     6675        # wrong values Win32
     6676
     6677       
     6678        from urs_ext import read_mux2
     6679       
     6680        from anuga.config import single_precision as epsilon       
     6681       
     6682        verbose = False
    66766683               
    6677                     # South is positive in MUX
     6684        tide = 1.5
     6685        time_step_count = 10
     6686        time_step = 0.02
     6687
     6688        '''
     6689        Win results
     6690        time_step = 0.2000001
     6691        This is OK       
     6692        '''
     6693       
     6694        '''
     6695        Win results
     6696        time_step = 0.20000001
     6697
     6698        ======================================================================
     6699ERROR: test_read_mux_platform_problem3 (__main__.Test_Data_Manager)
     6700----------------------------------------------------------------------
     6701Traceback (most recent call last):
     6702  File "test_data_manager.py", line 6718, in test_read_mux_platform_problem3
     6703    ha1[0]=num.sin(times_ref)
     6704ValueError: matrices are not aligned for copy
     6705
     6706        '''
     6707
     6708        '''
     6709        Win results
     6710        time_step = 0.200000001
     6711        FAIL
     6712         assert num.allclose(data[i][:parameters_index],
     6713         -va1[permutation[i], :])
     6714        '''
     6715        times_ref = num.arange(0, time_step_count*time_step, time_step)
     6716        #print "times_ref", times_ref
     6717       
     6718        lat_long_points = [(-21.5,114.5), (-21,114.5), (-21.5,115),
     6719                           (-21.,115.), (-22., 117.)]
     6720        stations = len(lat_long_points)
     6721       
     6722        # Create different timeseries starting and ending at different times
     6723        first_tstep=num.ones(stations, num.int)
     6724        first_tstep[0]+=2   # Point 0 starts at 2
     6725        first_tstep[1]+=4   # Point 1 starts at 4       
     6726        first_tstep[2]+=3   # Point 2 starts at 3
     6727       
     6728        last_tstep=(time_step_count)*num.ones(stations, num.int)
     6729        last_tstep[0]-=1    # Point 0 ends 1 step early
     6730        last_tstep[1]-=2    # Point 1 ends 2 steps early               
     6731        last_tstep[4]-=3    # Point 4 ends 3 steps early       
     6732       
     6733        # Create varying elevation data (positive values for seafloor)
     6734        gauge_depth=20*num.ones(stations, num.float)
     6735        for i in range(stations):
     6736            gauge_depth[i] += i**2
     6737           
     6738        # Create data to be written to second mux file       
     6739        ha1=num.ones((stations,time_step_count), num.float)
     6740        ha1[0]=num.sin(times_ref)
     6741        ha1[1]=2*num.sin(times_ref - 3)
     6742        ha1[2]=5*num.sin(4*times_ref)
     6743        ha1[3]=num.sin(times_ref)
     6744        ha1[4]=num.sin(2*times_ref-0.7)
     6745               
     6746        ua1=num.zeros((stations,time_step_count),num.float)
     6747        ua1[0]=3*num.cos(times_ref)       
     6748        ua1[1]=2*num.sin(times_ref-0.7)   
     6749        ua1[2]=num.arange(3*time_step_count,4*time_step_count)
     6750        ua1[4]=2*num.ones(time_step_count)
     6751       
     6752        va1=num.zeros((stations,time_step_count),num.float)
     6753        va1[0]=2*num.cos(times_ref-0.87)       
     6754        va1[1]=3*num.ones(time_step_count)
     6755        va1[3]=2*num.sin(times_ref-0.71)       
     6756        #print "va1[0]", va1[0]  # The 8th element is what will go bad.
     6757        # Ensure data used to write mux file to be zero when gauges are
     6758        # not recording
     6759        for i in range(stations):
     6760             # For each point
     6761             for j in range(0, first_tstep[i]-1) + range(last_tstep[i],
     6762                                                         time_step_count):
     6763                 # For timesteps before and after recording range
     6764                 ha1[i][j] = ua1[i][j] = va1[i][j] = 0.0
     6765
     6766
     6767        #print 'Second station to be written to MUX'
     6768        #print 'ha', ha1[0,:]
     6769        #print 'ua', ua1[0,:]
     6770        #print 'va', va1[0,:]
     6771       
     6772        # Write second mux file to be combined by urs2sts
     6773        base_nameII, filesII = self.write_mux2(lat_long_points,
     6774                                               time_step_count, time_step,
     6775                                               first_tstep, last_tstep,
     6776                                               depth=gauge_depth,
     6777                                               ha=ha1,
     6778                                               ua=ua1,
     6779                                               va=va1)
     6780        #print "filesII", filesII
     6781
     6782
     6783
     6784
     6785        # Read mux file back and verify it's correcness
     6786
     6787        ####################################################
     6788        # FIXME (Ole): This is where the test should
     6789        # verify that the MUX files are correct.
     6790
     6791        #JJ: It appears as though
     6792        #that certain quantities are not being stored with enough precision
     6793        #inn muxfile or more likely that they are being cast into a
     6794        #lower precision when read in using read_mux2 Time step and q_time
     6795        # are equal but only to approx 1e-7
     6796        ####################################################
     6797
     6798        #define information as it should be stored in mus2 files
     6799        points_num=len(lat_long_points)
     6800        depth=gauge_depth
     6801        ha=ha1
     6802        ua=ua1
     6803        va=va1
     6804       
     6805        quantities = ['HA','UA','VA']
     6806        mux_names = [WAVEHEIGHT_MUX2_LABEL,
     6807                     EAST_VELOCITY_MUX2_LABEL,
     6808                     NORTH_VELOCITY_MUX2_LABEL]
     6809        quantities_init = [[],[],[]]
     6810        latlondeps = []
     6811        #irrelevant header information
     6812        ig=ilon=ilat=0
     6813        mcolat=mcolon=centerlat=centerlon=offset=az=baz=id=0.0
     6814        # urs binary is latitude fastest
     6815        for i,point in enumerate(lat_long_points):
     6816            lat = point[0]
     6817            lon = point[1]
     6818            _ , e, n = redfearn(lat, lon)
     6819            if depth is None:
     6820                this_depth = n
     6821            else:
     6822                this_depth = depth[i]
     6823            latlondeps.append([lat, lon, this_depth])
     6824
     6825            if ha is None:
     6826                this_ha = e
     6827                quantities_init[0].append(num.ones(time_step_count,
     6828                                                   num.float)*this_ha) # HA
     6829            else:
     6830                quantities_init[0].append(ha[i])
     6831            if ua is None:
     6832                this_ua = n
     6833                quantities_init[1].append(num.ones(time_step_count,
     6834                                                   num.float)*this_ua) # UA
     6835            else:
     6836                quantities_init[1].append(ua[i])
     6837            if va is None:
     6838                this_va = e
     6839                quantities_init[2].append(num.ones(time_step_count,
     6840                                                   num.float)*this_va) #
     6841            else:
     6842                quantities_init[2].append(va[i])
     6843
     6844        for i, q in enumerate(quantities):
     6845            #print
     6846            #print i, q
     6847           
     6848            q_time = num.zeros((time_step_count, points_num), num.float)
     6849            quantities_init[i] = ensure_numeric(quantities_init[i])
     6850            for time in range(time_step_count):
     6851                #print i, q, time, quantities_init[i][:,time]
     6852                q_time[time,:] = quantities_init[i][:,time]
     6853                #print i, q, time, q_time[time, :]
     6854
     6855           
     6856            filename = base_nameII + mux_names[i]
     6857            f = open(filename, 'rb')
     6858            if self.verbose: print 'Reading' + filename
     6859            assert abs(points_num-unpack('i',f.read(4))[0])<epsilon
     6860            #write mux 2 header
     6861            for latlondep in latlondeps:
     6862                assert abs(latlondep[0]-unpack('f',f.read(4))[0])<epsilon
     6863                assert abs(latlondep[1]-unpack('f',f.read(4))[0])<epsilon
     6864                assert abs(mcolat-unpack('f',f.read(4))[0])<epsilon
     6865                assert abs(mcolon-unpack('f',f.read(4))[0])<epsilon
     6866                assert abs(ig-unpack('i',f.read(4))[0])<epsilon
     6867                assert abs(ilon-unpack('i',f.read(4))[0])<epsilon
     6868                assert abs(ilat-unpack('i',f.read(4))[0])<epsilon
     6869                assert abs(latlondep[2]-unpack('f',f.read(4))[0])<epsilon
     6870                assert abs(centerlat-unpack('f',f.read(4))[0])<epsilon
     6871                assert abs(centerlon-unpack('f',f.read(4))[0])<epsilon
     6872                assert abs(offset-unpack('f',f.read(4))[0])<epsilon
     6873                assert abs(az-unpack('f',f.read(4))[0])<epsilon
     6874                assert abs(baz-unpack('f',f.read(4))[0])<epsilon
     6875               
     6876                x = unpack('f', f.read(4))[0]
     6877                #print time_step
     6878                #print x
     6879                assert abs(time_step-x)<epsilon
     6880                assert abs(time_step_count-unpack('i',f.read(4))[0])<epsilon
     6881                for j in range(4): # identifier
     6882                    assert abs(id-unpack('i',f.read(4))[0])<epsilon
     6883
     6884            #first_tstep=1
     6885            #last_tstep=time_step_count
     6886            for i,latlondep in enumerate(latlondeps):
     6887                assert abs(first_tstep[i]-unpack('i',f.read(4))[0])<epsilon
     6888            for i,latlondep in enumerate(latlondeps):
     6889                assert abs(last_tstep[i]-unpack('i',f.read(4))[0])<epsilon
     6890
     6891            # Find when first station starts recording
     6892            min_tstep = min(first_tstep)
     6893            # Find when all stations have stopped recording
     6894            max_tstep = max(last_tstep)
     6895
     6896            #for time in  range(time_step_count):
     6897            for time in range(min_tstep-1,max_tstep):
     6898                assert abs(time*time_step-unpack('f',f.read(4))[0])<epsilon
     6899                for point_i in range(points_num):
     6900                    if time+1>=first_tstep[point_i] and time+1<=last_tstep[point_i]:
     6901                        x = unpack('f',f.read(4))[0]
     6902                        #print time, x, q_time[time, point_i]
     6903                        if q == 'VA': x = -x # South is positive in MUX
     6904                        #print q+" q_time[%d, %d] = %f" %(time, point_i,
     6905                                                      #q_time[time, point_i])
     6906                        assert abs(q_time[time, point_i]-x)<epsilon
     6907
     6908            f.close()
     6909                           
     6910        permutation = ensure_numeric([4,0,2])
     6911                   
     6912        # Create ordering file
     6913#         _, ordering_filename = tempfile.mkstemp('')
     6914#         order_fid = open(ordering_filename, 'w') 
     6915#         order_fid.write('index, longitude, latitude\n')
     6916#         for index in permutation:
     6917#             order_fid.write('%d, %f, %f\n' %(index,
     6918#                                              lat_long_points[index][1],
     6919#                                              lat_long_points[index][0]))
     6920#         order_fid.close()
     6921       
     6922        # -------------------------------------
     6923        # Now read files back and check values
     6924        weights = ensure_numeric([1.0])
     6925
     6926        # For each quantity read the associated list of source mux2 file with
     6927        # extention associated with that quantity
     6928        file_params=-1*num.ones(3,num.float) # [nsta,dt,nt]
     6929        OFFSET = 5
     6930
     6931        for j, file in enumerate(filesII):
     6932            # Read stage, u, v enumerated as j
     6933            #print 'Reading', j, file
     6934            #print "file", file
     6935            data = read_mux2(1, [file], weights, file_params,
     6936                             permutation, verbose)
     6937            #print str(j) + "data", data
     6938
     6939            #print 'Data received by Python'
     6940            #print data[1][8]
     6941            number_of_selected_stations = data.shape[0]
     6942            #print "number_of_selected_stations", number_of_selected_stations
     6943            #print "stations", stations
     6944
     6945            # Index where data ends and parameters begin
     6946            parameters_index = data.shape[1]-OFFSET         
     6947                 
     6948            for i in range(number_of_selected_stations):
     6949       
     6950                #print i, parameters_index
     6951                if j == 0:
     6952                    assert num.allclose(data[i][:parameters_index],
     6953                                        ha1[permutation[i], :])
     6954                   
     6955                if j == 1: assert num.allclose(data[i][:parameters_index], ua1[permutation[i], :])
     6956                if j == 2:
    66786957                    assert num.allclose(data[i][:parameters_index], -va1[permutation[i], :])
    6679                    
     6958       
     6959        self.delete_mux(filesII)           
    66806960       
    66816961    def test_urs2sts0(self):
     
    68197099
    68207100        base_name, files = self.write_mux2(lat_long_points,
    6821                                       time_step_count, time_step,
    6822                                       first_tstep, last_tstep,
    6823                                       depth=gauge_depth,
    6824                                       ha=ha,
    6825                                       ua=ua,
    6826                                       va=va)
     7101                                           time_step_count, time_step,
     7102                                           first_tstep, last_tstep,
     7103                                           depth=gauge_depth,
     7104                                           ha=ha,
     7105                                           ua=ua,
     7106                                           va=va)
    68277107
    68287108        urs2sts(base_name,
     
    68507130        y = points[:,1]
    68517131
    6852         #Check that all coordinate are correctly represented       
    6853         #Using the non standard projection (50)
     7132        # Check that all coordinate are correctly represented       
     7133        # Using the non standard projection (50)
    68547134        for i in range(4):
    68557135            zone, e, n = redfearn(lat_long_points[i][0],
     
    68587138            assert num.allclose([x[i],y[i]], [e,n])
    68597139            assert zone==-1
     7140       
     7141        self.delete_mux(files)
    68607142
    68617143           
     
    69157197        y = points[:,1]
    69167198
    6917         #Check that all coordinate are correctly represented       
    6918         #Using the non standard projection (50)
     7199        # Check that all coordinate are correctly represented       
     7200        # Using the non standard projection (50)
    69197201        for i in range(4):
    69207202            zone, e, n = redfearn(lat_long_points[i][0], lat_long_points[i][1],
     
    69227204            assert num.allclose([x[i],y[i]], [e,n])
    69237205            assert zone==geo_reference.zone
     7206       
     7207        self.delete_mux(files)
    69247208
    69257209           
     
    76527936            raise Exception, msg
    76537937
     7938       
     7939        self.delete_mux(filesI)
     7940        self.delete_mux(filesII)
    76547941
    76557942       
     
    80348321        os.remove(sts_file)
    80358322       
    8036        
    8037        
    80388323        #----------------------
    80398324        # Then read the mux files together and test
     
    81488433        os.remove(sts_file)
    81498434       
    8150 
    8151        
    81528435        #---------------
    81538436        # "Manually" add the timeseries up with weights and test
     
    81578440        stage_man = weights[0]*(stage0-tide) + weights[1]*(stage1-tide) + tide
    81588441        assert num.allclose(stage_man, stage)
    8159        
    8160        
    8161        
    8162        
    8163        
    81648442               
    81658443       
     
    82888566                       
    82898567        assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values,
    8290                             domain_drchlt.quantities['xmomentum'].vertex_values)                        
     8568                            domain_drchlt.quantities['xmomentum'].vertex_values)
    82918569                       
    82928570        assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values,
    8293                             domain_drchlt.quantities['ymomentum'].vertex_values)                                               
     8571                            domain_drchlt.quantities['ymomentum'].vertex_values)
    82948572       
    82958573       
    82968574        os.remove(sts_file+'.sts')
    82978575        os.remove(meshname)
    8298        
    8299        
    8300        
    83018576               
    83028577       
     
    84198694
    84208695
    8421 
    8422        
    8423        
    8424        
    8425            
    8426 
    8427        
    84288696        domain_drchlt = Domain(meshname)
    84298697        domain_drchlt.set_quantity('stage', tide)
     
    84478715                       
    84488716        assert num.allclose(domain_fbound.quantities['xmomentum'].vertex_values,
    8449                             domain_drchlt.quantities['xmomentum'].vertex_values)                        
     8717                            domain_drchlt.quantities['xmomentum'].vertex_values)
    84508718                       
    84518719        assert num.allclose(domain_fbound.quantities['ymomentum'].vertex_values,
    8452                             domain_drchlt.quantities['ymomentum'].vertex_values)                                               
    8453        
     8720                            domain_drchlt.quantities['ymomentum'].vertex_values)
    84548721       
    84558722        os.remove(sts_file+'.sts')
    84568723        os.remove(meshname)
    8457 
    8458 
    8459        
    84608724               
    84618725       
     
    1093411198        domain.set_quantity('xmomentum', uh)
    1093511199        domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
     11200
    1093611201        for t in domain.evolve(yieldstep=1, finaltime = t_end):
    1093711202            pass
     
    1135711622
    1135811623           
     11624    def test_copy_code_files(self):
     11625        '''test that the copy_code_files() function is sane.'''
     11626
     11627        def create_file(f):
     11628            fd = open(f, 'w')
     11629            fd.write('%s\n' % f)
     11630            fd.close()
     11631
     11632        # create working directories and test files
     11633        work_dir = tempfile.mkdtemp()
     11634        dst_dir = tempfile.mkdtemp(dir=work_dir)
     11635        src_dir = tempfile.mkdtemp(dir=work_dir)
     11636
     11637        f1 = 'file1'       
     11638        filename1 = os.path.join(src_dir, f1)
     11639        create_file(filename1)
     11640        f2 = 'file2'       
     11641        filename2 = os.path.join(src_dir, f2)
     11642        create_file(filename2)
     11643        f3 = 'file3'       
     11644        filename3 = os.path.join(src_dir, f3)
     11645        create_file(filename3)
     11646        f4 = 'file4'       
     11647        filename4 = os.path.join(src_dir, f4)
     11648        create_file(filename4)
     11649        f5 = 'file5'       
     11650        filename5 = os.path.join(src_dir, f5)
     11651        create_file(filename5)
     11652
     11653        # exercise the copy function
     11654        copy_code_files(dst_dir, filename1)
     11655        copy_code_files(dst_dir, filename1, filename2)
     11656        copy_code_files(dst_dir, (filename4, filename5, filename3))
     11657
     11658        # test that files were actually copied
     11659        self.failUnless(access(os.path.join(dst_dir, f1), F_OK))
     11660        self.failUnless(access(os.path.join(dst_dir, f2), F_OK))
     11661        self.failUnless(access(os.path.join(dst_dir, f3), F_OK))
     11662        self.failUnless(access(os.path.join(dst_dir, f4), F_OK))
     11663        self.failUnless(access(os.path.join(dst_dir, f5), F_OK))
     11664
     11665        # clean up
     11666        shutil.rmtree(work_dir)
    1135911667           
    1136011668#-------------------------------------------------------------
     
    1136211670if __name__ == "__main__":
    1136311671    suite = unittest.makeSuite(Test_Data_Manager,'test')
    11364 
    11365     # FIXME (Ole): This is the test that fails under Windows
    11366     #suite = unittest.makeSuite(Test_Data_Manager,'test_read_mux_platform_problem2')
    11367     #suite = unittest.makeSuite(Test_Data_Manager,'test_file_boundary_stsIV')
    1136811672   
    1136911673    if len(sys.argv) > 1 and sys.argv[1][0].upper() == 'V':
  • branches/numpy/anuga/shallow_water/test_shallow_water_domain.py

    r6689 r6902  
    2727def linear_function(point):
    2828    point = num.array(point)
    29     return point[:,0] + point[:,1]
    30 
     29    return point[:,0]+point[:,1]
    3130
    3231class Weir:
    3332    """Set a bathymetry for weir with a hole and a downstream gutter
    34 
    3533    x,y are assumed to be in the unit square
    3634    """
     
    4543        z = num.zeros(N, num.float)
    4644        for i in range(N):
    47             z[i] = -x[i] / 2    # General slope
    48 
    49             # Flattish bit to the left
     45            z[i] = -x[i]/2  #General slope
     46
     47            #Flattish bit to the left
    5048            if x[i] < 0.3:
    5149                z[i] = -x[i]/10
    5250
    53             # Weir
     51            #Weir
    5452            if x[i] >= 0.3 and x[i] < 0.4:
    5553                z[i] = -x[i]+0.9
    5654
    57             # Dip
     55            #Dip
    5856            x0 = 0.6
    5957            depth = -1.0
     
    6260                if x[i] > x0 and x[i] < 0.9:
    6361                    z[i] = depth
    64                 # RHS plateaux
     62                #RHS plateaux
    6563                if x[i] >= 0.9:
    6664                    z[i] = plateaux
    6765            elif y[i] >= 0.7 and y[i] < 1.5:
    68                 # Restrict and deepen
     66                #Restrict and deepen
    6967                if x[i] >= x0 and x[i] < 0.8:
    7068                    z[i] = depth - (y[i]/3 - 0.3)
    7169                elif x[i] >= 0.8:
    72                     # RHS plateaux
     70                    #RHS plateaux
    7371                    z[i] = plateaux
    7472            elif y[i] >= 1.5:
    7573                if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2:
    76                     # Widen up and stay at constant depth
    77                     z[i] = depth - 1.5/5
    78                 elif x[i] >= 0.8 + (y[i] - 1.5)/1.2:
    79                     # RHS plateaux
     74                    #Widen up and stay at constant depth
     75                    z[i] = depth-1.5/5
     76                elif x[i] >= 0.8 + (y[i]-1.5)/1.2:
     77                    #RHS plateaux
    8078                    z[i] = plateaux
    8179
    82             # Hole in weir (slightly higher than inflow condition)
     80            #Hole in weir (slightly higher than inflow condition)
    8381            if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
    84                 z[i] = -x[i] + self.inflow_stage + 0.02
    85 
    86             # Channel behind weir
     82                z[i] = -x[i]+self.inflow_stage + 0.02
     83
     84            #Channel behind weir
    8785            x0 = 0.5
    8886            if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4:
    89                 z[i] = -x[i] + self.inflow_stage + 0.02
     87                z[i] = -x[i]+self.inflow_stage + 0.02
    9088
    9189            if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4:
    92                 # Flatten it out towards the end
    93                 z[i] = -x0 + self.inflow_stage + 0.02 + (x0 - x[i])/5
     90                #Flatten it out towards the end
     91                z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5
    9492
    9593            # Hole to the east
     
    9795            y0 = 0.35
    9896            if num.sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
    99                 z[i] = num.sqrt((x[i]-x0)**2 + (y[i]-y0)**2) - 1.0
    100 
    101             # Tiny channel draining hole
     97                z[i] = num.sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-1.0
     98
     99            #Tiny channel draining hole
    102100            if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6:
    103                 z[i] = -0.9    # North south
     101                z[i] = -0.9 #North south
    104102
    105103            if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65:
    106                 z[i] = -1.0 + (x[i]-0.9)/3    # East west
     104                z[i] = -1.0 + (x[i]-0.9)/3 #East west
    107105
    108106            # Stuff not in use
     
    136134        z = num.zeros(N, num.float)
    137135        for i in range(N):
    138             z[i] = -x[i]    # General slope
    139 
    140             # Flat bit to the left
     136            z[i] = -x[i]  #General slope
     137
     138            #Flat bit to the left
    141139            if x[i] < 0.3:
    142                 z[i] = -x[i]/10    # General slope
    143 
    144             # Weir
     140                z[i] = -x[i]/10  #General slope
     141
     142            #Weir
    145143            if x[i] > 0.3 and x[i] < 0.4:
    146                 z[i] = -x[i] + 0.9
    147 
    148             # Dip
     144                z[i] = -x[i]+0.9
     145
     146            #Dip
    149147            if x[i] > 0.6 and x[i] < 0.9:
    150                 z[i] = -x[i] - 0.5    # -y[i]/5
    151 
    152             # Hole in weir (slightly higher than inflow condition)
     148                z[i] = -x[i]-0.5  #-y[i]/5
     149
     150            #Hole in weir (slightly higher than inflow condition)
    153151            if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
    154                 z[i] = -x[i] + self.inflow_stage + 0.05
     152                z[i] = -x[i]+self.inflow_stage + 0.05
     153
    155154
    156155        return z/2
     
    170169
    171170    N = len(x)
    172     s = 0 * x    # New array
     171    s = 0*x  #New array
    173172
    174173    for k in range(N):
     
    263262        H0 = 0.0
    264263
    265         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
    266                                   epsilon, g, H0)
    267 
    268         assert num.allclose(edgeflux, [0, 0, 0])
     264       
     265        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)
     266
     267        assert num.allclose(edgeflux, [0,0,0])
    269268        assert max_speed == 0.
    270269
     
    272271        w = 2.0
    273272
    274         normal = num.array([1., 0])
     273        normal = num.array([1.,0])
    275274        ql = num.array([w, 0, 0])
    276275        qr = num.array([w, 0, 0])
     
    280279        H0 = 0.0
    281280
    282         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
    283                                   epsilon, g, H0)
    284 
     281        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)       
    285282        assert num.allclose(edgeflux, [0., 0.5*g*h**2, 0.])
    286283        assert max_speed == num.sqrt(g*h)
     
    290287    #    w = 2.0
    291288    #
    292     #    normal = array([1., 0])
     289    #    normal = array([1.,0])
    293290    #    ql = array([w, 0, 0])
    294291    #    qr = array([w, 0, 0])
     
    15171514        uh = u*h
    15181515
    1519         Br = Reflective_boundary(domain)       # Side walls
    1520         Bd = Dirichlet_boundary([w, uh, 0])    # 2 m/s across the 3 m inlet:
     1516        Br = Reflective_boundary(domain)     # Side walls
     1517        Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet:
     1518
    15211519
    15221520        # Initial conditions
  • branches/numpy/anuga/shallow_water/urs_ext.c

    r6780 r6902  
    104104                memcpy(data + it, muxData + offset, sizeof(float));
    105105               
    106                 //printf("%d: muxdata=%f\n", it, muxData[offset]);                     
    107                 //printf("data[%d]=%f, offset=%d\n", it, data[it], offset);       
     106                //printf("%d: muxdata=%f\n", it, muxData[offset]); 
     107                //printf("data[%d]=%f, offset=%d\n", it, data[it], offset);
    108108                offset++;
    109109            }
     
    221221
    222222        // Open the mux file
    223         if((fp = fopen(muxFileName, "r")) == NULL)
     223        if((fp = fopen(muxFileName, "rb")) == NULL)
    224224        {
    225225            char *err_msg = strerror(errno);
     
    237237            }
    238238       
    239             fros = (int*) malloc(*total_number_of_stations*numSrc*sizeof(int)); 
     239            fros = (int*) malloc(*total_number_of_stations*numSrc*sizeof(int));
    240240            lros = (int*) malloc(*total_number_of_stations*numSrc*sizeof(int));
    241241     
     
    447447    temp_sts_data = (float*) calloc(len_sts_data, sizeof(float));
    448448
    449     muxData = (float*) malloc(numDataMax*sizeof(float));
     449    muxData = (float*) calloc(numDataMax, sizeof(float));
    450450   
    451451    // Loop over all sources
     
    460460        // Read in data block from mux2 file
    461461        muxFileName = muxFileNameArray[isrc];
    462         if((fp = fopen(muxFileName, "r")) == NULL)
     462        if((fp = fopen(muxFileName, "rb")) == NULL)
    463463        {
    464464            fprintf(stderr, "cannot open file %s\n", muxFileName);
     
    470470        }
    471471
    472         offset = sizeof(int) + total_number_of_stations*(sizeof(struct tgsrwg) + 2*sizeof(int));
    473         fseek(fp, offset, 0);
     472        offset = (long int)sizeof(int) + total_number_of_stations*(sizeof(struct tgsrwg) + 2*sizeof(int));
     473        //printf("\n offset %i ", (long int)offset);
     474        fseek(fp, offset, 0);
    474475
    475476        numData = getNumData(fros_per_source,
    476477                             lros_per_source,
    477478                             total_number_of_stations);
    478                              
    479                              
    480                                      
    481         elements_read = fread(muxData, ((int) numData)*sizeof(float), 1, fp);
     479        // Note numData is larger than what it has to be.                   
     480        //elements_read = fread(muxData, ((int) numData)*sizeof(float), 1, fp);
     481        elements_read = fread(muxData, (size_t) sizeof(float), (size_t) numData, fp);
     482        //printf("\n elements_read  %d, ", (int)elements_read);
     483        //printf("\n ferror(fp)  %d, ", (int)ferror(fp));
    482484        if ((int) elements_read == 0 && ferror(fp)) {
    483485          fprintf(stderr, "Error reading mux data\n");
    484486          return NULL;
    485         }
    486                
    487         fclose(fp);
     487        }       
    488488       
    489 
    490         // FIXME (Ole): This is where Nariman and Ole traced the platform dependent
    491         // difference on 11 November 2008. We don't think the problem lies in the
    492         // C code. Maybe it is a problem with the MUX files written by the unit test
    493         // that fails on Windows but works OK on Linux. JJ's test on 17th November shows
    494         // that as far as Python is concerned this file should be OK on both platforms.
    495        
    496        
    497        
    498         // These printouts are enough to show the problem when compared
    499         // on the two platforms
    500         //printf("\nRead %d elements, ", (int) numData);
    501         //printf("muxdata[%d]=%f\n", 39, muxData[39]);         
    502        
    503        
    504         /*
    505         In Windows we get
    506        
    507         Read 85 elements, muxdata[39]=0.999574
    508         Read 85 elements, muxdata[39]=-0.087599
    509         Read 85 elements, muxdata[39]=-0.087599
    510        
    511        
    512         In Linux we get (the correct)
    513        
    514         Read 85 elements, muxdata[39]=0.999574
    515         Read 85 elements, muxdata[39]=-0.087599
    516         Read 85 elements, muxdata[39]=1.490349
    517         */
    518        
    519        
    520        
    521        
     489        fclose(fp); 
    522490       
    523491        // loop over stations present in the permutation array
  • branches/numpy/anuga/utilities/log.py

    r6888 r6902  
    6060NOTSET = logging.NOTSET
    6161
     62# get True if python version 2.5 or later
     63(version_major, version_minor, _, _, _) = sys.version_info
     64new_python = ((version_major == 2 and version_minor >= 5) or version_major > 2)
     65
    6266
    6367################################################################################
     
    7983    global _setup, log_logging_level
    8084
    81     # get running python version for later
    82     (version_major, version_minor, _, _, _) = sys.version_info
    83 
    8485    # have we been setup?
    8586    if not _setup:
     
    8990
    9091        # setup the file logging system
    91         if version_major >= 2 and version_minor >= 5:
     92        if new_python:
    9293            fmt = '%(asctime)s %(levelname)-8s %(mname)25s:%(lnum)-4d|%(message)s'
    9394        else:
     
    109110                        logging.getLevelName(log_logging_level),
    110111                        logging.getLevelName(console_logging_level)))
    111         if version_major >= 2 and version_minor >= 5:
     112        if new_python:
    112113            logging.log(logging.CRITICAL, start_msg,
    113114                        extra={'mname': __name__, 'lnum': 0})
     
    127128            break
    128129
    129     if version_major >= 2 and version_minor >= 5:
     130    if new_python:
    130131        logging.log(level, msg, extra={'mname': fname, 'lnum': lnum})
    131132    else:
     
    139140# @brief Shortcut for log(DEBUG, msg).
    140141# @param msg Message string to log at logging.DEBUG level.
    141 def debug(msg):
     142def debug(msg=''):
    142143    log(logging.DEBUG, msg)
    143144
     
    145146# @brief Shortcut for log(INFO, msg).
    146147# @param msg Message string to log at logging.INFO level.
    147 def info(msg):
     148def info(msg=''):
    148149    log(logging.INFO, msg)
    149150
     
    151152# @brief Shortcut for log(WARNING, msg).
    152153# @param msg Message string to log at logging.WARNING level.
    153 def warning(msg):
     154def warning(msg=''):
    154155    log(logging.WARNING, msg)
    155156
     
    157158# @brief Shortcut for log(ERROR, msg).
    158159# @param msg Message string to log at logging.ERROR level.
    159 def error(msg):
     160def error(msg=''):
    160161    log(logging.ERROR, msg)
    161162
     
    163164# @brief Shortcut for log(CRITICAL, msg).
    164165# @param msg Message string to log at logging.CRITICAL level.
    165 def critical(msg):
     166def critical(msg=''):
    166167    log(logging.CRITICAL, msg)
    167168
     
    219220    else:
    220221        pass
    221 
    222 if __name__ == "__main__":
    223     debug('DEBUG message')
    224     info('INFO message')
    225     warning('WARNING message')
  • branches/numpy/anuga/utilities/log_test.py

    r6702 r6902  
    2727
    2828    def tearDown(self):
    29         pass
    30 #        if os.path.exists(LOGFILE_NAME):
    31 #            os.remove(LOGFILE_NAME)
    32 #        if os.path.exists(STDOUT_LOG_NAME):
    33 #            os.remove(STDOUT_LOG_NAME)
     29        if os.path.exists(LOGFILE_NAME):
     30            os.remove(LOGFILE_NAME)
     31        if os.path.exists(STDOUT_LOG_NAME):
     32            os.remove(STDOUT_LOG_NAME)
    3433
    3534    ##
Note: See TracChangeset for help on using the changeset viewer.