Changeset 6451


Ignore:
Timestamp:
Mar 4, 2009, 3:04:42 PM (16 years ago)
Author:
rwilson
Message:

Fixed bug introduced in numpy changeover.

Location:
branches/numpy/anuga/shallow_water
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/numpy/anuga/shallow_water/shallow_water_domain.py

    r6441 r6451  
    16021602            msg = 'Function %s could not be executed:\n%s' %(f, e)
    16031603            # FIXME: Reconsider this semantics
    1604             raise msg
     1604            raise Exception, msg
    16051605
    16061606        try:
    16071607            q = num.array(q, num.float)
    16081608        except:
    1609             msg = 'Return value from vector function %s could ' %f
    1610             msg += 'not be converted into a numeric array of floats.\n'
    1611             msg += 'Specified function should return either list or array.'
     1609            msg = ('Return value from vector function %s could not '
     1610                   'be converted into a numeric array of floats.\nSpecified '
     1611                   'function should return either list or array.' % f)
    16121612            raise Exception, msg
    16131613
     
    16161616        func_info = (f.func_name, f.func_code.co_filename,
    16171617                     f.func_code.co_firstlineno)
    1618         msg = ('Function %s() must return vector (defined in %s, line %d)'
    1619                % func_info)
    1620         assert hasattr(q, 'len'), msg
    1621 
    1622         msg = ('Return vector from function %s() must have same '
    1623                'length as input vectors\nq=%s' % (f.func_name, str(q)))
    1624         assert len(q) == N, msg
     1618        func_msg = 'Function %s (defined in %s, line %d)' % func_info
     1619        try:
     1620            result_len = len(q)
     1621        except:
     1622            msg = '%s must return vector' % func_msg
     1623            self.fail(msg)
     1624        msg = '%s must return vector of length %d' % (func_msg, N)
     1625        assert result_len == N, msg
    16251626    else:
    16261627        try:
  • branches/numpy/anuga/shallow_water/test_shallow_water_domain.py

    <
    r6441 r6451  
    77from anuga.config import g, epsilon
    88from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    9 import numpy as num
    109from anuga.utilities.numerical_tools import mean
    1110from anuga.utilities.polygon import is_inside_polygon
     
    1413from anuga.geospatial_data.geospatial_data import Geospatial_data
    1514from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
     15
    1616from shallow_water_domain import *
     17
     18import numpy as num
    1719
    1820# Get gateway to C implementation of flux function for direct testing
    1921from shallow_water_ext import flux_function_central as flux_function
     22
    2023
    2124# For test_fitting_using_shallow_water_domain example
    2225def linear_function(point):
    2326    point = num.array(point)
    24     return point[:,0]+point[:,1]
     27    return point[:,0] + point[:,1]
     28
    2529
    2630class Weir:
    2731    """Set a bathymetry for weir with a hole and a downstream gutter
     32
    2833    x,y are assumed to be in the unit square
    2934    """
     
    3338
    3439    def __call__(self, x, y):
    35 
    3640        N = len(x)
    3741        assert N == len(y)
     
    3943        z = num.zeros(N, num.float)
    4044        for i in range(N):
    41             z[i] = -x[i]/2  #General slope
    42 
    43             #Flattish bit to the left
     45            z[i] = -x[i] / 2    # General slope
     46
     47            # Flattish bit to the left
    4448            if x[i] < 0.3:
    4549                z[i] = -x[i]/10
    4650
    47             #Weir
     51            # Weir
    4852            if x[i] >= 0.3 and x[i] < 0.4:
    4953                z[i] = -x[i]+0.9
    5054
    51             #Dip
     55            # Dip
    5256            x0 = 0.6
    53             #depth = -1.3
    5457            depth = -1.0
    55             #plateaux = -0.9
    5658            plateaux = -0.6
    5759            if y[i] < 0.7:
    5860                if x[i] > x0 and x[i] < 0.9:
    5961                    z[i] = depth
    60 
    61                 #RHS plateaux
     62                # RHS plateaux
    6263                if x[i] >= 0.9:
    6364                    z[i] = plateaux
    64 
    65 
    6665            elif y[i] >= 0.7 and y[i] < 1.5:
    67                 #Restrict and deepen
     66                # Restrict and deepen
    6867                if x[i] >= x0 and x[i] < 0.8:
    69                     z[i] = depth-(y[i]/3-0.3)
    70                     #z[i] = depth-y[i]/5
    71                     #z[i] = depth
     68                    z[i] = depth - (y[i]/3 - 0.3)
    7269                elif x[i] >= 0.8:
    73                     #RHS plateaux
     70                    # RHS plateaux
    7471                    z[i] = plateaux
    75 
    7672            elif y[i] >= 1.5:
    7773                if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2:
    78                     #Widen up and stay at constant depth
    79                     z[i] = depth-1.5/5
    80                 elif x[i] >= 0.8 + (y[i]-1.5)/1.2:
    81                     #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
    8278                    z[i] = plateaux
    8379
    84 
    85             #Hole in weir (slightly higher than inflow condition)
     80            # Hole in weir (slightly higher than inflow condition)
    8681            if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
    87                 z[i] = -x[i]+self.inflow_stage + 0.02
    88 
    89             #Channel behind weir
     82                z[i] = -x[i] + self.inflow_stage + 0.02
     83
     84            # Channel behind weir
    9085            x0 = 0.5
    9186            if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4:
    92                 z[i] = -x[i]+self.inflow_stage + 0.02
     87                z[i] = -x[i] + self.inflow_stage + 0.02
    9388
    9489            if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4:
    95                 #Flatten it out towards the end
    96                 z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5
    97 
    98             #Hole to the east
    99             x0 = 1.1; y0 = 0.35
    100             #if x[i] < -0.2 and y < 0.5:
     90                # Flatten it out towards the end
     91                z[i] = -x0 + self.inflow_stage + 0.02 + (x0 - x[i])/5
     92
     93            # Hole to the east
     94            x0 = 1.1
     95            y0 = 0.35
    10196            if num.sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
    102                 z[i] = num.sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-1.0
    103 
    104             #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
    105100            if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6:
    106                 z[i] = -0.9 #North south
     101                z[i] = -0.9    # North south
    107102
    108103            if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65:
    109                 z[i] = -1.0 + (x[i]-0.9)/3 #East west
    110 
    111 
    112 
    113             #Stuff not in use
    114 
    115             #Upward slope at inlet to the north west
    116             #if x[i] < 0.0: # and y[i] > 0.5:
     104                z[i] = -1.0 + (x[i]-0.9)/3    # East west
     105
     106            # Stuff not in use
     107
     108            # Upward slope at inlet to the north west
     109            # if x[i] < 0.0: # and y[i] > 0.5:
    117110            #    #z[i] = -y[i]+0.5  #-x[i]/2
    118111            #    z[i] = x[i]/4 - y[i]**2 + 0.5
    119112
    120             #Hole to the west
    121             #x0 = -0.4; y0 = 0.35 # center
    122             #if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
     113            # Hole to the west
     114            # x0 = -0.4; y0 = 0.35 # center
     115            # if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
    123116            #    z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-0.2
    124117
    125 
    126 
    127 
    128 
    129118        return z/2
     119
    130120
    131121class Weir_simple:
    132122    """Set a bathymetry for weir with a hole and a downstream gutter
     123
    133124    x,y are assumed to be in the unit square
    134125    """
     
    138129
    139130    def __call__(self, x, y):
    140 
    141131        N = len(x)
    142132        assert N == len(y)
     
    144134        z = num.zeros(N, num.float)
    145135        for i in range(N):
    146             z[i] = -x[i]  #General slope
    147 
    148             #Flat bit to the left
     136            z[i] = -x[i]    # General slope
     137
     138            # Flat bit to the left
    149139            if x[i] < 0.3:
    150                 z[i] = -x[i]/10  #General slope
    151 
    152             #Weir
     140                z[i] = -x[i]/10    # General slope
     141
     142            # Weir
    153143            if x[i] > 0.3 and x[i] < 0.4:
    154                 z[i] = -x[i]+0.9
    155 
    156             #Dip
     144                z[i] = -x[i] + 0.9
     145
     146            # Dip
    157147            if x[i] > 0.6 and x[i] < 0.9:
    158                 z[i] = -x[i]-0.5  #-y[i]/5
    159 
    160             #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)
    161151            if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
    162                 z[i] = -x[i]+self.inflow_stage + 0.05
    163 
     152                z[i] = -x[i] + self.inflow_stage + 0.05
    164153
    165154        return z/2
    166155
    167156
    168 
    169 
    170 #Variable windfield implemented using functions
    171 def speed(t,x,y):
     157# Variable windfield implemented using functions
     158def speed(t, x, y):
    172159    """Large speeds halfway between center and edges
     160
    173161    Low speeds at center and edges
    174162    """
     
    180168
    181169    N = len(x)
    182     s = 0*x  #New array
     170    s = 0 * x    # New array
    183171
    184172    for k in range(N):
    185 
    186173        r = num.sqrt(x[k]**2 + y[k]**2)
    187 
    188         factor = exp( -(r-0.15)**2 )
    189 
     174        factor = exp(-(r-0.15)**2)
    190175        s[k] = 4000 * factor * (cos(t*2*pi/150) + 2)
    191176
    192177    return s
    193178
    194 
    195 def scalar_func(t,x,y):
     179def scalar_func(t, x, y):
    196180    """Function that returns a scalar.
     181
    197182    Used to test error message when numeric array is expected
    198183    """
     
    200185    return 17.7
    201186
    202 
    203 def scalar_func_list(t,x,y):
     187def scalar_func_list(t, x, y):
    204188    """Function that returns a scalar.
     189
    205190    Used to test error message when numeric array is expected
    206191    """
     
    208193    return [17.7]
    209194
    210 
    211 def angle(t,x,y):
     195def angle(t, x, y):
    212196    """Rotating field
    213197    """
     
    218202
    219203    N = len(x)
    220     a = 0*x  #New array
     204    a = 0 * x    # New array
    221205
    222206    for k in range(N):
     
    226210
    227211        if x[k] < 0:
    228             angle+=pi #atan in ]-pi/2; pi/2[
    229 
    230         #Take normal direction
     212            angle += pi
     213
     214        # Take normal direction
    231215        angle -= pi/2
    232216
    233         #Ensure positive radians
     217        # Ensure positive radians
    234218        if angle < 0:
    235219            angle += 2*pi
     
    248232
    249233    def test_rotate(self):
    250         normal = num.array([0.0,-1.0])
    251 
    252         q = num.array([1.0,2.0,3.0])
     234        normal = num.array([0.0, -1.0])
     235
     236        q = num.array([1.0, 2.0, 3.0])
    253237
    254238        r = rotate(q, normal, direction = 1)
     
    260244        assert num.allclose(w, q)
    261245
    262         #Check error check
     246        # Check error check
    263247        try:
    264             rotate(r, num.array([1,1,1]))
     248            rotate(r, num.array([1, 1, 1]))
    265249        except:
    266250            pass
    267251        else:
    268             raise 'Should have raised an exception'
    269 
     252            raise Exception, 'Should have raised an exception'
    270253
    271254    # Individual flux tests
    272255    def test_flux_zero_case(self):
    273         ql = num.zeros( 3, num.float )
    274         qr = num.zeros( 3, num.float )
    275         normal = num.zeros( 2, num.float )
    276         edgeflux = num.zeros( 3, num.float )
     256        ql = num.zeros(3, num.float)
     257        qr = num.zeros(3, num.float)
     258        normal = num.zeros(2, num.float)
     259        edgeflux = num.zeros(3, num.float)
    277260        zl = zr = 0.
    278261        H0 = 0.0
    279        
    280         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)
    281 
    282         assert num.allclose(edgeflux, [0,0,0])
     262
     263        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
     264                                  epsilon, g, H0)
     265
     266        assert num.allclose(edgeflux, [0, 0, 0])
    283267        assert max_speed == 0.
    284268
     
    286270        w = 2.0
    287271
    288         normal = num.array([1.,0])
     272        normal = num.array([1., 0])
    289273        ql = num.array([w, 0, 0])
    290274        qr = num.array([w, 0, 0])
    291         edgeflux = num.zeros(3, num.float)       
     275        edgeflux = num.zeros(3, num.float)
    292276        zl = zr = 0.
    293277        h = w - (zl+zr)/2
    294278        H0 = 0.0
    295279
    296         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)       
     280        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
     281                                  epsilon, g, H0)
     282
    297283        assert num.allclose(edgeflux, [0., 0.5*g*h**2, 0.])
    298284        assert max_speed == num.sqrt(g*h)
     
    302288    #    w = 2.0
    303289    #
    304     #    normal = array([1.,0])
     290    #    normal = array([1., 0])
    305291    #    ql = array([w, 0, 0])
    306292    #    qr = array([w, 0, 0])
     
    313299    #    assert max_speed == sqrt(g*h)
    314300
    315 
    316301    def test_flux1(self):
    317         #Use data from previous version of abstract_2d_finite_volumes
    318         normal = num.array([1.,0])
     302        # Use data from previous version of abstract_2d_finite_volumes
     303        normal = num.array([1., 0])
    319304        ql = num.array([-0.2, 2, 3])
    320305        qr = num.array([-0.2, 2, 3])
    321306        zl = zr = -0.5
    322         edgeflux = num.zeros(3, num.float)               
     307        edgeflux = num.zeros(3, num.float)
    323308
    324309        H0 = 0.0
    325310
    326         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)       
    327 
    328         assert num.allclose(edgeflux, [2.,13.77433333, 20.])
     311        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
     312                                  epsilon, g, H0)
     313
     314        assert num.allclose(edgeflux, [2., 13.77433333, 20.])
    329315        assert num.allclose(max_speed, 8.38130948661)
    330316
    331 
    332317    def test_flux2(self):
    333         #Use data from previous version of abstract_2d_finite_volumes
     318        # Use data from previous version of abstract_2d_finite_volumes
    334319        normal = num.array([0., -1.])
    335320        ql = num.array([-0.075, 2, 3])
     
    337322        zl = zr = -0.375
    338323
    339         edgeflux = num.zeros(3, num.float)               
     324        edgeflux = num.zeros(3, num.float)
    340325        H0 = 0.0
    341         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)       
    342 
    343         assert num.allclose(edgeflux, [-3.,-20.0, -30.441])
     326        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
     327                                  epsilon, g, H0)
     328
     329        assert num.allclose(edgeflux, [-3., -20.0, -30.441])
    344330        assert num.allclose(max_speed, 11.7146428199)
    345331
    346332    def test_flux3(self):
    347         #Use data from previous version of abstract_2d_finite_volumes
     333        # Use data from previous version of abstract_2d_finite_volumes
    348334        normal = num.array([-sqrt(2)/2, sqrt(2)/2])
    349335        ql = num.array([-0.075, 2, 3])
     
    351337        zl = zr = -0.375
    352338
    353         edgeflux = num.zeros(3, num.float)               
     339        edgeflux = num.zeros(3, num.float)
    354340        H0 = 0.0
    355         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)       
     341        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
     342                                  epsilon, g, H0)
    356343
    357344        assert num.allclose(edgeflux, [sqrt(2)/2, 4.40221112, 7.3829019])
     
    359346
    360347    def test_flux4(self):
    361         #Use data from previous version of abstract_2d_finite_volumes
     348        # Use data from previous version of abstract_2d_finite_volumes
    362349        normal = num.array([-sqrt(2)/2, sqrt(2)/2])
    363350        ql = num.array([-0.34319278, 0.10254161, 0.07273855])
     
    365352        zl = zr = -0.375
    366353
    367         edgeflux = num.zeros(3, num.float)               
     354        edgeflux = num.zeros(3, num.float)
    368355        H0 = 0.0
    369         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)               
     356        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
     357                                  epsilon, g, H0)
    370358
    371359        assert num.allclose(edgeflux, [-0.04072676, -0.07096636, -0.01604364])
    372360        assert num.allclose(max_speed, 1.31414103233)
    373361
    374     def test_flux_computation(self):   
    375         """test_flux_computation - test flux calculation (actual C implementation)
    376         This one tests the constant case where only the pressure term contributes to each edge and cancels out
    377         once the total flux has been summed up.
     362    def test_flux_computation(self):
     363        """test flux calculation (actual C implementation)
     364
     365        This one tests the constant case where only the pressure term
     366        contributes to each edge and cancels out once the total flux has
     367        been summed up.
    378368        """
    379                
     369
    380370        a = [0.0, 0.0]
    381371        b = [0.0, 2.0]
    382         c = [2.0,0.0]
     372        c = [2.0, 0.0]
    383373        d = [0.0, 4.0]
    384374        e = [2.0, 2.0]
    385         f = [4.0,0.0]
     375        f = [4.0, 0.0]
    386376
    387377        points = [a, b, c, d, e, f]
    388         #bac, bce, ecf, dbe, daf, dae
     378        #              bac,     bce,     ecf,     dbe
    389379        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    390380
     
    392382        domain.check_integrity()
    393383
    394         # The constant case             
     384        # The constant case
    395385        domain.set_quantity('elevation', -1)
    396         domain.set_quantity('stage', 1) 
    397        
     386        domain.set_quantity('stage', 1)
     387
    398388        domain.compute_fluxes()
    399         assert num.allclose(domain.get_quantity('stage').explicit_update[1], 0) # Central triangle
    400        
    401 
    402         # The more general case                 
    403         def surface(x,y):
    404             return -x/2                   
    405        
     389        # Central triangle
     390        assert num.allclose(domain.get_quantity('stage').explicit_update[1], 0)
     391
     392        # The more general case
     393        def surface(x, y):
     394            return -x/2
     395
    406396        domain.set_quantity('elevation', -10)
    407         domain.set_quantity('stage', surface)   
    408         domain.set_quantity('xmomentum', 1)             
    409        
     397        domain.set_quantity('stage', surface)
     398        domain.set_quantity('xmomentum', 1)
     399
    410400        domain.compute_fluxes()
    411        
     401
    412402        #print domain.get_quantity('stage').explicit_update
    413403        # FIXME (Ole): TODO the general case
    414         #assert allclose(domain.get_quantity('stage').explicit_update[1], ........??)
    415                
    416        
    417                
     404        #assert allclose(domain.get_quantity('stage').explicit_update[1], ...??)
     405
    418406    def test_sw_domain_simple(self):
    419407        a = [0.0, 0.0]
    420408        b = [0.0, 2.0]
    421         c = [2.0,0.0]
     409        c = [2.0, 0.0]
    422410        d = [0.0, 4.0]
    423411        e = [2.0, 2.0]
    424         f = [4.0,0.0]
     412        f = [4.0, 0.0]
    425413
    426414        points = [a, b, c, d, e, f]
    427         #bac, bce, ecf, dbe, daf, dae
    428         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    429 
     415        #             bac,     bce,     ecf,     dbe
     416        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    430417
    431418        #from anuga.abstract_2d_finite_volumes.domain import Domain as Generic_domain
     
    443430        assert num.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.)
    444431
    445 
    446432    def test_boundary_conditions(self):
    447 
    448433        a = [0.0, 0.0]
    449434        b = [0.0, 2.0]
    450         c = [2.0,0.0]
     435        c = [2.0, 0.0]
    451436        d = [0.0, 4.0]
    452437        e = [2.0, 2.0]
    453         f = [4.0,0.0]
     438        f = [4.0, 0.0]
    454439
    455440        points = [a, b, c, d, e, f]
    456         #bac, bce, ecf, dbe
    457         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
    458         boundary = { (0, 0): 'Third',
    459                      (0, 2): 'First',
    460                      (2, 0): 'Second',
    461                      (2, 1): 'Second',
    462                      (3, 1): 'Second',
    463                      (3, 2): 'Third'}
    464 
     441        #             bac,     bce,     ecf,     dbe
     442        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     443        boundary = {(0, 0): 'Third',
     444                    (0, 2): 'First',
     445                    (2, 0): 'Second',
     446                    (2, 1): 'Second',
     447                    (3, 1): 'Second',
     448                    (3, 2): 'Third'}
    465449
    466450        domain = Domain(points, vertices, boundary)
    467451        domain.check_integrity()
    468452
    469 
    470         domain.set_quantity('stage', [[1,2,3], [5,5,5],
    471                                       [0,0,9], [-6, 3, 3]])
    472 
    473         domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
    474                                           [3,3,3], [4, 4, 4]])
     453        domain.set_quantity('stage', [[1,2,3], [5,5,5], [0,0,9], [-6,3,3]])
     454
     455        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], [3,3,3], [4,4,4]])
    475456
    476457        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
    477                                           [30,30,30], [40, 40, 40]])
    478 
    479 
    480         D = Dirichlet_boundary([5,2,1])
     458                                          [30,30,30], [40,40,40]])
     459
     460        D = Dirichlet_boundary([5, 2, 1])
    481461        T = Transmissive_boundary(domain)
    482462        R = Reflective_boundary(domain)
    483         domain.set_boundary( {'First': D, 'Second': T, 'Third': R})
     463        domain.set_boundary({'First': D, 'Second': T, 'Third': R})
    484464
    485465        domain.update_boundary()
    486466
    487         #Stage
     467        # Stage
    488468        assert domain.quantities['stage'].boundary_values[0] == 2.5
    489         assert domain.quantities['stage'].boundary_values[0] ==\
    490                domain.get_conserved_quantities(0, edge=0)[0] #Reflective (2.5)
    491         assert domain.quantities['stage'].boundary_values[1] == 5. #Dirichlet
    492         assert domain.quantities['stage'].boundary_values[2] ==\
    493                domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5)
    494         assert domain.quantities['stage'].boundary_values[3] ==\
    495                domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5)
    496         assert domain.quantities['stage'].boundary_values[4] ==\
    497                domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5)
    498         assert domain.quantities['stage'].boundary_values[5] ==\
    499                domain.get_conserved_quantities(3, edge=2)[0] #Reflective (-1.5)
    500 
    501         #Xmomentum
    502         assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 #Reflective
    503         assert domain.quantities['xmomentum'].boundary_values[1] == 2. #Dirichlet
    504         assert domain.quantities['xmomentum'].boundary_values[2] ==\
    505                domain.get_conserved_quantities(2, edge=0)[1] #Transmissive
    506         assert domain.quantities['xmomentum'].boundary_values[3] ==\
    507                domain.get_conserved_quantities(2, edge=1)[1] #Transmissive
    508         assert domain.quantities['xmomentum'].boundary_values[4] ==\
    509                domain.get_conserved_quantities(3, edge=1)[1] #Transmissive
    510         assert domain.quantities['xmomentum'].boundary_values[5] == -4.0  #Reflective
    511 
    512         #Ymomentum
    513         assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 #Reflective
    514         assert domain.quantities['ymomentum'].boundary_values[1] == 1.  #Dirichlet
    515         assert domain.quantities['ymomentum'].boundary_values[2] == 30. #Transmissive
    516         assert domain.quantities['ymomentum'].boundary_values[3] == 30. #Transmissive
    517         assert domain.quantities['ymomentum'].boundary_values[4] == 40. #Transmissive
    518         assert domain.quantities['ymomentum'].boundary_values[5] == 40. #Reflective
    519 
     469        # Reflective (2.5)
     470        assert (domain.quantities['stage'].boundary_values[0] ==
     471                domain.get_conserved_quantities(0, edge=0)[0])
     472        # Dirichlet
     473        assert domain.quantities['stage'].boundary_values[1] == 5.
     474        # Transmissive (4.5)
     475        assert (domain.quantities['stage'].boundary_values[2] ==
     476                domain.get_conserved_quantities(2, edge=0)[0])
     477        # Transmissive (4.5)
     478        assert (domain.quantities['stage'].boundary_values[3] ==
     479                domain.get_conserved_quantities(2, edge=1)[0])
     480        # Transmissive (-1.5)
     481        assert (domain.quantities['stage'].boundary_values[4] ==
     482                domain.get_conserved_quantities(3, edge=1)[0])
     483        # Reflective (-1.5)
     484        assert (domain.quantities['stage'].boundary_values[5] ==
     485                domain.get_conserved_quantities(3, edge=2)[0])
     486
     487        # Xmomentum
     488        # Reflective
     489        assert domain.quantities['xmomentum'].boundary_values[0] == 1.0
     490        # Dirichlet
     491        assert domain.quantities['xmomentum'].boundary_values[1] == 2.
     492        # Transmissive
     493        assert (domain.quantities['xmomentum'].boundary_values[2] ==
     494                domain.get_conserved_quantities(2, edge=0)[1])
     495        # Transmissive
     496        assert (domain.quantities['xmomentum'].boundary_values[3] ==
     497                domain.get_conserved_quantities(2, edge=1)[1])
     498        # Transmissive
     499        assert (domain.quantities['xmomentum'].boundary_values[4] ==
     500                domain.get_conserved_quantities(3, edge=1)[1])
     501        # Reflective
     502        assert domain.quantities['xmomentum'].boundary_values[5] == -4.0
     503
     504        # Ymomentum
     505        # Reflective
     506        assert domain.quantities['ymomentum'].boundary_values[0] == -10.0
     507        # Dirichlet
     508        assert domain.quantities['ymomentum'].boundary_values[1] == 1.
     509        # Transmissive
     510        assert domain.quantities['ymomentum'].boundary_values[2] == 30.
     511        # Transmissive
     512        assert domain.quantities['ymomentum'].boundary_values[3] == 30.
     513        # Transmissive
     514        assert domain.quantities['ymomentum'].boundary_values[4] == 40.
     515        # Reflective
     516        assert domain.quantities['ymomentum'].boundary_values[5] == 40.
    520517
    521518    def test_boundary_conditionsII(self):
    522 
    523519        a = [0.0, 0.0]
    524520        b = [0.0, 2.0]
    525         c = [2.0,0.0]
     521        c = [2.0, 0.0]
    526522        d = [0.0, 4.0]
    527523        e = [2.0, 2.0]
    528         f = [4.0,0.0]
     524        f = [4.0, 0.0]
    529525
    530526        points = [a, b, c, d, e, f]
    531         #bac, bce, ecf, dbe
    532         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
    533         boundary = { (0, 0): 'Third',
    534                      (0, 2): 'First',
    535                      (2, 0): 'Second',
    536                      (2, 1): 'Second',
    537                      (3, 1): 'Second',
    538                      (3, 2): 'Third',
    539                      (0, 1): 'Internal'}
    540 
     527        #             bac,     bce,     ecf,     dbe
     528        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     529        boundary = {(0, 0): 'Third',
     530                    (0, 2): 'First',
     531                    (2, 0): 'Second',
     532                    (2, 1): 'Second',
     533                    (3, 1): 'Second',
     534                    (3, 2): 'Third',
     535                    (0, 1): 'Internal'}
    541536
    542537        domain = Domain(points, vertices, boundary)
    543538        domain.check_integrity()
    544539
    545 
    546         domain.set_quantity('stage', [[1,2,3], [5,5,5],
    547                                       [0,0,9], [-6, 3, 3]])
    548 
    549         domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
    550                                           [3,3,3], [4, 4, 4]])
     540        domain.set_quantity('stage', [[1,2,3], [5,5,5], [0,0,9], [-6,3,3]])
     541
     542        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], [3,3,3], [4,4,4]])
    551543
    552544        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
    553                                           [30,30,30], [40, 40, 40]])
    554 
    555 
    556         D = Dirichlet_boundary([5,2,1])
     545                                          [30,30,30], [40,40,40]])
     546
     547        D = Dirichlet_boundary([5, 2, 1])
    557548        T = Transmissive_boundary(domain)
    558549        R = Reflective_boundary(domain)
    559         domain.set_boundary( {'First': D, 'Second': T,
    560                               'Third': R, 'Internal': None})
     550        domain.set_boundary({'First': D, 'Second': T,
     551                             'Third': R, 'Internal': None})
    561552
    562553        domain.update_boundary()
    563554        domain.check_integrity()
    564555
    565        
    566 
    567        
    568556    def test_boundary_conditionsIII(self):
    569557        """test_boundary_conditionsIII
    570        
     558
    571559        Test Transmissive_stage_zero_momentum_boundary
    572560        """
    573        
     561
    574562        a = [0.0, 0.0]
    575563        b = [0.0, 2.0]
    576         c = [2.0,0.0]
     564        c = [2.0, 0.0]
    577565        d = [0.0, 4.0]
    578566        e = [2.0, 2.0]
    579         f = [4.0,0.0]
     567        f = [4.0, 0.0]
    580568
    581569        points = [a, b, c, d, e, f]
    582         #bac, bce, ecf, dbe
    583         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
    584         boundary = { (0, 0): 'Third',
    585                      (0, 2): 'First',
    586                      (2, 0): 'Second',
    587                      (2, 1): 'Second',
    588                      (3, 1): 'Second',
    589                      (3, 2): 'Third'}
    590 
     570        #             bac,     bce,     ecf,     dbe
     571        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     572        boundary = {(0, 0): 'Third',
     573                    (0, 2): 'First',
     574                    (2, 0): 'Second',
     575                    (2, 1): 'Second',
     576                    (3, 1): 'Second',
     577                    (3, 2): 'Third'}
    591578
    592579        domain = Domain(points, vertices, boundary)
    593580        domain.check_integrity()
    594581
    595 
    596         domain.set_quantity('stage', [[1,2,3], [5,5,5],
    597                                       [0,0,9], [-6, 3, 3]])
    598 
    599         domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
    600                                           [3,3,3], [4, 4, 4]])
     582        domain.set_quantity('stage', [[1,2,3], [5,5,5], [0,0,9], [-6,3,3]])
     583
     584        domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], [3,3,3], [4,4,4]])
    601585
    602586        domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
    603                                           [30,30,30], [40, 40, 40]])
    604 
    605 
    606         D = Dirichlet_boundary([5,2,1])
     587                                          [30,30,30], [40,40,40]])
     588
     589        D = Dirichlet_boundary([5, 2, 1])
    607590        T = Transmissive_stage_zero_momentum_boundary(domain)
    608591        R = Reflective_boundary(domain)
    609         domain.set_boundary( {'First': D, 'Second': T, 'Third': R})
     592        domain.set_boundary({'First': D, 'Second': T, 'Third': R})
    610593
    611594        domain.update_boundary()
    612595
    613596        # Stage
     597        # Reflective (2.5)
    614598        assert domain.quantities['stage'].boundary_values[0] == 2.5
    615         assert domain.quantities['stage'].boundary_values[0] ==\
    616                domain.get_conserved_quantities(0, edge=0)[0] #Reflective (2.5)
    617         assert domain.quantities['stage'].boundary_values[1] == 5. #Dirichlet
    618         assert domain.quantities['stage'].boundary_values[2] ==\
    619                domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5)
    620         assert domain.quantities['stage'].boundary_values[3] ==\
    621                domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5)
    622         assert domain.quantities['stage'].boundary_values[4] ==\
    623                domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5)
    624         assert domain.quantities['stage'].boundary_values[5] ==\
    625                domain.get_conserved_quantities(3, edge=2)[0] #Reflective (-1.5)
     599        assert (domain.quantities['stage'].boundary_values[0] ==
     600                domain.get_conserved_quantities(0, edge=0)[0])
     601        # Dirichlet
     602        assert domain.quantities['stage'].boundary_values[1] == 5.
     603        # Transmissive (4.5)
     604        assert (domain.quantities['stage'].boundary_values[2] ==
     605                domain.get_conserved_quantities(2, edge=0)[0])
     606        # Transmissive (4.5)
     607        assert (domain.quantities['stage'].boundary_values[3] ==
     608                domain.get_conserved_quantities(2, edge=1)[0])
     609        # Transmissive (-1.5)
     610        assert (domain.quantities['stage'].boundary_values[4] ==
     611                domain.get_conserved_quantities(3, edge=1)[0])
     612        # Reflective (-1.5)
     613        assert (domain.quantities['stage'].boundary_values[5] ==
     614                domain.get_conserved_quantities(3, edge=2)[0])
    626615
    627616        # Xmomentum
    628         assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 #Reflective
    629         assert domain.quantities['xmomentum'].boundary_values[1] == 2. #Dirichlet
    630         assert domain.quantities['xmomentum'].boundary_values[2] == 0.0
    631         assert domain.quantities['xmomentum'].boundary_values[3] == 0.0
    632         assert domain.quantities['xmomentum'].boundary_values[4] == 0.0
    633         assert domain.quantities['xmomentum'].boundary_values[5] == -4.0  #Reflective
     617        # Reflective
     618        assert domain.quantities['xmomentum'].boundary_values[0] == 1.0
     619        # Dirichlet
     620        assert domain.quantities['xmomentum'].boundary_values[1] == 2.
     621        assert domain.quantities['xmomentum'].boundary_values[2] == 0.0
     622        assert domain.quantities['xmomentum'].boundary_values[3] == 0.0
     623        assert domain.quantities['xmomentum'].boundary_values[4] == 0.0
     624        # Reflective
     625        assert domain.quantities['xmomentum'].boundary_values[5] == -4.0
    634626
    635627        # Ymomentum
    636         assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 #Reflective
    637         assert domain.quantities['ymomentum'].boundary_values[1] == 1.  #Dirichlet
    638         assert domain.quantities['ymomentum'].boundary_values[2] == 0.0
    639         assert domain.quantities['ymomentum'].boundary_values[3] == 0.0
    640         assert domain.quantities['ymomentum'].boundary_values[4] == 0.0
    641         assert domain.quantities['ymomentum'].boundary_values[5] == 40. #Reflective
    642 
    643                
    644                
    645        
     628        # Reflective
     629        assert domain.quantities['ymomentum'].boundary_values[0] == -10.0
     630        # Dirichlet
     631        assert domain.quantities['ymomentum'].boundary_values[1] == 1.
     632        assert domain.quantities['ymomentum'].boundary_values[2] == 0.0
     633        assert domain.quantities['ymomentum'].boundary_values[3] == 0.0
     634        assert domain.quantities['ymomentum'].boundary_values[4] == 0.0
     635        # Reflective
     636        assert domain.quantities['ymomentum'].boundary_values[5] == 40.
     637
    646638    def test_boundary_condition_time(self):
    647639        """test_boundary_condition_time
    648        
     640
    649641        This tests that boundary conditions are evaluated
    650642        using the right time from domain.
    651 
    652643        """
    653        
     644
    654645        # Setup computational domain
    655         from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
    656 
     646        from anuga.abstract_2d_finite_volumes.mesh_factory \
     647                import rectangular_cross
    657648
    658649        #--------------------------------------------------------------
     
    660651        #--------------------------------------------------------------
    661652        N = 5
    662         points, vertices, boundary = rectangular_cross(N, N) 
     653        points, vertices, boundary = rectangular_cross(N, N)
    663654        domain = Domain(points, vertices, boundary)
    664655
     
    670661        domain.set_quantity('stage', 0.0)
    671662
    672 
    673663        #--------------------------------------------------------------
    674664        # Setup boundary conditions
    675665        #--------------------------------------------------------------
    676         Bt = Time_boundary(domain=domain,             # Time dependent boundary 
    677                           f=lambda t: [t, 0.0, 0.0])
    678        
     666        # Time dependent boundary
     667        Bt = Time_boundary(domain=domain, f=lambda t: [t, 0.0, 0.0])
     668
    679669        Br = Reflective_boundary(domain)              # Reflective wall
    680670
    681671        domain.set_boundary({'left': Bt, 'right': Br, 'top': Br, 'bottom': Br})
    682        
     672
    683673        for t in domain.evolve(yieldstep = 10, finaltime = 20.0):
    684674            q = Bt.evaluate()
    685    
    686             # FIXME (Ole): This test would not have passed in 
     675
     676            # FIXME (Ole): This test would not have passed in
    687677            # changeset:5846.
    688678            msg = 'Time boundary not evaluated correctly'
    689679            assert num.allclose(t, q[0]), msg
    690            
    691        
    692680
    693681    def test_compute_fluxes0(self):
     
    697685        a = [0.0, 0.0]
    698686        b = [0.0, 2.0]
    699         c = [2.0,0.0]
     687        c = [2.0, 0.0]
    700688        d = [0.0, 4.0]
    701689        e = [2.0, 2.0]
    702         f = [4.0,0.0]
     690        f = [4.0, 0.0]
    703691
    704692        points = [a, b, c, d, e, f]
    705         #bac, bce, ecf, dbe
    706         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     693        #             bac,     bce,     ecf,    dbe
     694        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    707695
    708696        domain = Domain(points, vertices)
    709         domain.set_quantity('stage', [[2,2,2], [2,2,2],
    710                                       [2,2,2], [2,2,2]])
     697        domain.set_quantity('stage', [[2,2,2], [2,2,2], [2,2,2], [2,2,2]])
    711698        domain.check_integrity()
    712699
    713         assert num.allclose(domain.neighbours, [[-1,1,-1], [2,3,0], [-1,-1,1],[1,-1,-1]])
    714         assert num.allclose(domain.neighbour_edges, [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]])
    715 
    716         zl=zr=0. # Assume flat bed
    717 
    718         edgeflux = num.zeros(3, num.float)       
     700        assert num.allclose(domain.neighbours,
     701                            [[-1,1,-1], [2,3,0], [-1,-1,1],[1,-1,-1]])
     702        assert num.allclose(domain.neighbour_edges,
     703                            [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]])
     704
     705        zl = zr = 0.     # Assume flat bed
     706
     707        edgeflux = num.zeros(3, num.float)
    719708        edgeflux0 = num.zeros(3, num.float)
    720709        edgeflux1 = num.zeros(3, num.float)
    721         edgeflux2 = num.zeros(3, num.float)                               
    722         H0 = 0.0       
     710        edgeflux2 = num.zeros(3, num.float)
     711        H0 = 0.0
    723712
    724713        # Flux across right edge of volume 1
    725         normal = domain.get_normal(1,0)
     714        normal = domain.get_normal(1, 0)
    726715        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
    727716        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
    728         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, epsilon, g, H0)                       
     717        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0,
     718                                  epsilon, g, H0)
    729719
    730720        # Check that flux seen from other triangles is inverse
    731         tmp = qr; qr=ql; ql=tmp
    732         normal = domain.get_normal(2,2)
    733         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)                               
     721        (ql, qr) = (qr, ql)
     722        #tmp = qr
     723        #qr = ql
     724        #ql = tmp
     725        normal = domain.get_normal(2, 2)
     726        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
     727                                  epsilon, g, H0)
    734728
    735729        assert num.allclose(edgeflux0 + edgeflux, 0.)
    736730
    737731        # Flux across upper edge of volume 1
    738         normal = domain.get_normal(1,1)
     732        normal = domain.get_normal(1, 1)
    739733        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
    740734        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
    741         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, epsilon, g, H0)                                       
     735        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1,
     736                                  epsilon, g, H0)
    742737
    743738        # Check that flux seen from other triangles is inverse
    744         tmp = qr; qr=ql; ql=tmp
    745         normal = domain.get_normal(3,0)
    746         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)                                               
    747 
    748         assert num.allclose(edgeflux1 + edgeflux, 0.)       
    749        
     739        (ql, qr) = (qr, ql)
     740        #tmp = qr
     741        #qr = ql
     742        #ql = tmp
     743        normal = domain.get_normal(3, 0)
     744        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
     745                                  epsilon, g, H0)
     746
     747        assert num.allclose(edgeflux1 + edgeflux, 0.)
    750748
    751749        # Flux across lower left hypotenuse of volume 1
    752         normal = domain.get_normal(1,2)
     750        normal = domain.get_normal(1, 2)
    753751        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
    754752        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
    755         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, epsilon, g, H0)                                                               
     753        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2,
     754                                  epsilon, g, H0)
    756755
    757756        # Check that flux seen from other triangles is inverse
    758         tmp = qr; qr=ql; ql=tmp
    759         normal = domain.get_normal(0,1)
    760         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)                                                       
     757        (ql, qr) = (qr, ql)
     758        #tmp = qr
     759        #qr=ql
     760        #ql=tmp
     761        normal = domain.get_normal(0, 1)
     762        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux,
     763                                  epsilon, g, H0)
    761764        assert num.allclose(edgeflux2 + edgeflux, 0.)
    762 
    763765
    764766        # Scale by edgelengths, add up anc check that total flux is zero
     
    767769        e2 = domain.edgelengths[1, 2]
    768770
    769         assert num.allclose(e0*edgeflux0+e1*edgeflux1+e2*edgeflux2, 0.)
     771        assert num.allclose(e0*edgeflux0 + e1*edgeflux1 + e2*edgeflux2, 0.)
    770772
    771773        # Now check that compute_flux yields zeros as well
     
    773775
    774776        for name in ['stage', 'xmomentum', 'ymomentum']:
    775             #print name, domain.quantities[name].explicit_update
    776777            assert num.allclose(domain.quantities[name].explicit_update[1], 0)
    777 
    778 
    779778
    780779    def test_compute_fluxes1(self):
    781780        #Use values from previous version
    782 
    783781        a = [0.0, 0.0]
    784782        b = [0.0, 2.0]
    785         c = [2.0,0.0]
     783        c = [2.0, 0.0]
    786784        d = [0.0, 4.0]
    787785        e = [2.0, 2.0]
    788         f = [4.0,0.0]
     786        f = [4.0, 0.0]
    789787
    790788        points = [a, b, c, d, e, f]
    791         #bac, bce, ecf, dbe
    792         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     789        #             bac,     bce,     ecf,    dbe
     790        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    793791
    794792        domain = Domain(points, vertices)
    795         val0 = 2.+2.0/3
    796         val1 = 4.+4.0/3
    797         val2 = 8.+2.0/3
    798         val3 = 2.+8.0/3
     793        val0 = 2. + 2.0/3
     794        val1 = 4. + 4.0/3
     795        val2 = 8. + 2.0/3
     796        val3 = 2. + 8.0/3
    799797
    800798        domain.set_quantity('stage', [[val0, val0, val0], [val1, val1, val1],
     
    802800        domain.check_integrity()
    803801
    804         zl=zr=0. #Assume flat bed
    805 
    806         edgeflux = num.zeros(3, num.float)       
     802        zl = zr = 0.    # Assume flat bed
     803
     804        edgeflux = num.zeros(3, num.float)
    807805        edgeflux0 = num.zeros(3, num.float)
    808806        edgeflux1 = num.zeros(3, num.float)
    809         edgeflux2 = num.zeros(3, num.float)                               
    810         H0 = 0.0       
    811        
     807        edgeflux2 = num.zeros(3, num.float)
     808        H0 = 0.0
    812809
    813810        # Flux across right edge of volume 1
    814         normal = domain.get_normal(1,0) #Get normal 0 of triangle 1
     811        normal = domain.get_normal(1, 0)    # Get normal 0 of triangle 1
    815812        assert num.allclose(normal, [1, 0])
    816        
     813
    817814        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
    818815        assert num.allclose(ql, [val1, 0, 0])
    819        
     816
    820817        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
    821818        assert num.allclose(qr, [val2, 0, 0])
    822819
    823         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, epsilon, g, H0)                                                       
     820        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0,
     821                                  epsilon, g, H0)
    824822
    825823        # Flux across edge in the east direction (as per normal vector)
     
    827825        assert num.allclose(max_speed, 9.21592824046)
    828826
    829 
    830827        #Flux across edge in the west direction (opposite sign for xmomentum)
    831         normal_opposite = domain.get_normal(2,2) #Get normal 2 of triangle 2
     828        normal_opposite = domain.get_normal(2, 2)   # Get normal 2 of triangle 2
    832829        assert num.allclose(normal_opposite, [-1, 0])
    833830
    834         max_speed = flux_function(normal_opposite, ql, qr, zl, zr, edgeflux, epsilon, g, H0)                                             
    835         #flux_opposite, max_speed = flux_function([-1, 0], ql, qr, zl, zr)
     831        max_speed = flux_function(normal_opposite, ql, qr, zl, zr, edgeflux,
     832                                  epsilon, g, H0)
    836833        assert num.allclose(edgeflux, [-15.3598804, -253.71111111, 0.])
    837        
    838834
    839835        #Flux across upper edge of volume 1
    840         normal = domain.get_normal(1,1)
     836        normal = domain.get_normal(1, 1)
    841837        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
    842838        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
    843         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, epsilon, g, H0)                                                               
     839        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1,
     840                                  epsilon, g, H0)
    844841
    845842        assert num.allclose(edgeflux1, [2.4098563, 0., 123.04444444])
     
    847844
    848845        #Flux across lower left hypotenuse of volume 1
    849         normal = domain.get_normal(1,2)
     846        normal = domain.get_normal(1, 2)
    850847        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
    851848        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
    852         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, epsilon, g, H0)       
     849        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2,
     850                                  epsilon, g, H0)
    853851
    854852        assert num.allclose(edgeflux2, [9.63942522, -61.59685738, -61.59685738])
     
    860858        e2 = domain.edgelengths[1, 2]
    861859
    862         total_flux = -(e0*edgeflux0+e1*edgeflux1+e2*edgeflux2)/domain.areas[1]
     860        total_flux = -(e0*edgeflux0 +
     861                       e1*edgeflux1 +
     862                       e2*edgeflux2) / domain.areas[1]
     863
    863864        assert num.allclose(total_flux, [-0.68218178, -166.6, -35.93333333])
    864865
    865 
    866866        domain.compute_fluxes()
    867867
    868         #assert num.allclose(total_flux, domain.explicit_update[1,:])
    869868        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
    870869            assert num.allclose(total_flux[i],
    871870                                domain.quantities[name].explicit_update[1])
    872 
    873         #assert allclose(domain.explicit_update, [
    874         #    [0., -69.68888889, -69.68888889],
    875         #    [-0.68218178, -166.6, -35.93333333],
    876         #    [-111.77316251, 69.68888889, 0.],
    877         #    [-35.68522449, 0., 69.68888889]])
    878871
    879872        assert num.allclose(domain.quantities['stage'].explicit_update,
     
    884877                            [-69.68888889, -35.93333333, 0., 69.68888889])
    885878
    886 
    887         #assert allclose(domain.quantities[name].explicit_update
    888 
    889 
    890 
    891 
    892 
    893879    def test_compute_fluxes2(self):
    894880        #Random values, incl momentum
    895 
    896881        a = [0.0, 0.0]
    897882        b = [0.0, 2.0]
    898         c = [2.0,0.0]
     883        c = [2.0, 0.0]
    899884        d = [0.0, 4.0]
    900885        e = [2.0, 2.0]
    901         f = [4.0,0.0]
     886        f = [4.0, 0.0]
    902887
    903888        points = [a, b, c, d, e, f]
    904         #bac, bce, ecf, dbe
    905         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     889        #             bac,     bce,     ecf,    dbe
     890        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    906891
    907892        domain = Domain(points, vertices)
     893        val0 = 2. + 2.0/3
     894        val1 = 4. + 4.0/3
     895        val2 = 8. + 2.0/3
     896        val3 = 2. + 8.0/3
     897
     898        zl = zr = 0    # Assume flat zero bed
     899        edgeflux = num.zeros(3, num.float)
     900        edgeflux0 = num.zeros(3, num.float)
     901        edgeflux1 = num.zeros(3, num.float)
     902        edgeflux2 = num.zeros(3, num.float)
     903        H0 = 0.0
     904
     905        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
     906
     907        domain.set_quantity('stage', [[val0, val0-1, val0-2],
     908                                      [val1, val1+1, val1],
     909                                      [val2, val2-2, val2],
     910                                      [val3-0.5, val3, val3]])
     911
     912        domain.set_quantity('xmomentum',
     913                            [[1,2,3], [3,4,5], [1,-1,0], [0,-2,2]])
     914
     915        domain.set_quantity('ymomentum',
     916                            [[1,-1,0], [0,-3,2], [0,1,0], [-1,2,2]])
     917
     918        domain.check_integrity()
     919
     920        # Flux across right edge of volume 1
     921        normal = domain.get_normal(1, 0)
     922        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
     923        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
     924        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0,
     925                                  epsilon, g, H0)
     926
     927        # Flux across upper edge of volume 1
     928        normal = domain.get_normal(1, 1)
     929        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
     930        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
     931        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1,
     932                                  epsilon, g, H0)
     933
     934        # Flux across lower left hypotenuse of volume 1
     935        normal = domain.get_normal(1, 2)
     936        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
     937        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
     938        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2,
     939                                  epsilon, g, H0)
     940
     941        # Scale, add up and check that compute_fluxes is correct for vol 1
     942        e0 = domain.edgelengths[1, 0]
     943        e1 = domain.edgelengths[1, 1]
     944        e2 = domain.edgelengths[1, 2]
     945
     946        total_flux = -(e0*edgeflux0 +
     947                       e1*edgeflux1 +
     948                       e2*edgeflux2) / domain.areas[1]
     949
     950        domain.compute_fluxes()
     951
     952        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
     953            assert num.allclose(total_flux[i],
     954                                domain.quantities[name].explicit_update[1])
     955
     956    # FIXME (Ole): Need test like this for fluxes in very shallow water.
     957    def test_compute_fluxes3(self):
     958        #Random values, incl momentum
     959        a = [0.0, 0.0]
     960        b = [0.0, 2.0]
     961        c = [2.0, 0.0]
     962        d = [0.0, 4.0]
     963        e = [2.0, 2.0]
     964        f = [4.0, 0.0]
     965
     966        points = [a, b, c, d, e, f]
     967        #             bac,     bce,     ecf,     dbe
     968        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     969
     970        domain = Domain(points, vertices)
     971
    908972        val0 = 2.+2.0/3
    909973        val1 = 4.+4.0/3
     
    911975        val3 = 2.+8.0/3
    912976
    913         zl=zr=0 #Assume flat zero bed
    914         edgeflux = num.zeros(3, num.float)       
     977        zl = zr = -3.75    # Assume constant bed (must be less than stage)
     978        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
     979
     980        edgeflux = num.zeros(3, num.float)
    915981        edgeflux0 = num.zeros(3, num.float)
    916982        edgeflux1 = num.zeros(3, num.float)
    917         edgeflux2 = num.zeros(3, num.float)                               
    918         H0 = 0.0       
    919        
    920 
    921         domain.set_quantity('elevation', zl*num.ones( (4,3), num.int )) #array default#
    922 
     983        edgeflux2 = num.zeros(3, num.float)
     984        H0 = 0.0
    923985
    924986        domain.set_quantity('stage', [[val0, val0-1, val0-2],
     
    928990
    929991        domain.set_quantity('xmomentum',
    930                             [[1, 2, 3], [3, 4, 5],
    931                              [1, -1, 0], [0, -2, 2]])
     992                            [[1,2,3], [3,4,5], [1,-1,0], [0,-2,2]])
    932993
    933994        domain.set_quantity('ymomentum',
    934                             [[1, -1, 0], [0, -3, 2],
    935                              [0, 1, 0], [-1, 2, 2]])
    936 
     995                            [[1,-1,0], [0,-3,2], [0,1,0], [-1,2,2]])
    937996
    938997        domain.check_integrity()
    939998
    940 
    941 
    942         #Flux across right edge of volume 1
    943         normal = domain.get_normal(1,0)
     999        # Flux across right edge of volume 1
     1000        normal = domain.get_normal(1, 0)
    9441001        ql = domain.get_conserved_quantities(vol_id=1, edge=0)
    9451002        qr = domain.get_conserved_quantities(vol_id=2, edge=2)
    946         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, epsilon, g, H0)               
    947 
    948         #Flux across upper edge of volume 1
    949         normal = domain.get_normal(1,1)
     1003        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0,
     1004                                  epsilon, g, H0)
     1005
     1006        # Flux across upper edge of volume 1
     1007        normal = domain.get_normal(1, 1)
    9501008        ql = domain.get_conserved_quantities(vol_id=1, edge=1)
    9511009        qr = domain.get_conserved_quantities(vol_id=3, edge=0)
    952         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, epsilon, g, H0)                       
    953 
    954         #Flux across lower left hypotenuse of volume 1
    955         normal = domain.get_normal(1,2)
     1010        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1,
     1011                                  epsilon, g, H0)
     1012
     1013        # Flux across lower left hypotenuse of volume 1
     1014        normal = domain.get_normal(1, 2)
    9561015        ql = domain.get_conserved_quantities(vol_id=1, edge=2)
    9571016        qr = domain.get_conserved_quantities(vol_id=0, edge=1)
    958         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, epsilon, g, H0)               
    959 
    960         #Scale, add up and check that compute_fluxes is correct for vol 1
     1017        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2,
     1018                                  epsilon, g, H0)
     1019
     1020        # Scale, add up and check that compute_fluxes is correct for vol 1
    9611021        e0 = domain.edgelengths[1, 0]
    9621022        e1 = domain.edgelengths[1, 1]
    9631023        e2 = domain.edgelengths[1, 2]
    9641024
    965         total_flux = -(e0*edgeflux0+e1*edgeflux1+e2*edgeflux2)/domain.areas[1]
    966 
     1025        total_flux = -(e0*edgeflux0 +
     1026                       e1*edgeflux1 +
     1027                       e2*edgeflux2) / domain.areas[1]
    9671028
    9681029        domain.compute_fluxes()
     1030
    9691031        for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
    9701032            assert num.allclose(total_flux[i],
    9711033                                domain.quantities[name].explicit_update[1])
    972         #assert allclose(total_flux, domain.explicit_update[1,:])
    973 
    974 
    975     # FIXME (Ole): Need test like this for fluxes in very shallow water.   
    976     def test_compute_fluxes3(self):
    977         #Random values, incl momentum
     1034
     1035    def xtest_catching_negative_heights(self):
     1036        #OBSOLETE
    9781037
    9791038        a = [0.0, 0.0]
    9801039        b = [0.0, 2.0]
    981         c = [2.0,0.0]
     1040        c = [2.0, 0.0]
    9821041        d = [0.0, 4.0]
    9831042        e = [2.0, 2.0]
    984         f = [4.0,0.0]
     1043        f = [4.0, 0.0]
    9851044
    9861045        points = [a, b, c, d, e, f]
    987         #bac, bce, ecf, dbe
    988         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     1046        #             bac,     bce,     ecf,    dbe
     1047        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    9891048
    9901049        domain = Domain(points, vertices)
    991         val0 = 2.+2.0/3
    992         val1 = 4.+4.0/3
    993         val2 = 8.+2.0/3
    994         val3 = 2.+8.0/3
    995 
    996         zl=zr=-3.75 #Assume constant bed (must be less than stage)
    997         domain.set_quantity('elevation', zl*num.ones( (4,3), num.int )) #array default#
    998 
    999 
    1000         edgeflux = num.zeros(3, num.float)       
    1001         edgeflux0 = num.zeros(3, num.float)
    1002         edgeflux1 = num.zeros(3, num.float)
    1003         edgeflux2 = num.zeros(3, num.float)                               
    1004         H0 = 0.0       
    1005        
    1006 
    1007 
    1008         domain.set_quantity('stage', [[val0, val0-1, val0-2],
    1009                                       [val1, val1+1, val1],
    1010                                       [val2, val2-2, val2],
    1011                                       [val3-0.5, val3, val3]])
    1012 
    1013         domain.set_quantity('xmomentum',
    1014                             [[1, 2, 3], [3, 4, 5],
    1015                              [1, -1, 0], [0, -2, 2]])
    1016 
    1017         domain.set_quantity('ymomentum',
    1018                             [[1, -1, 0], [0, -3, 2],
    1019                              [0, 1, 0], [-1, 2, 2]])
    1020 
    1021 
    1022         domain.check_integrity()
    1023 
    1024 
    1025 
    1026         #Flux across right edge of volume 1
    1027         normal = domain.get_normal(1,0)
    1028         ql = domain.get_conserved_quantities(vol_id=1, edge=0)
    1029         qr = domain.get_conserved_quantities(vol_id=2, edge=2)
    1030         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux0, epsilon, g, H0)
    1031 
    1032         #Flux across upper edge of volume 1
    1033         normal = domain.get_normal(1,1)
    1034         ql = domain.get_conserved_quantities(vol_id=1, edge=1)
    1035         qr = domain.get_conserved_quantities(vol_id=3, edge=0)
    1036         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux1, epsilon, g, H0)       
    1037 
    1038         #Flux across lower left hypotenuse of volume 1
    1039         normal = domain.get_normal(1,2)
    1040         ql = domain.get_conserved_quantities(vol_id=1, edge=2)
    1041         qr = domain.get_conserved_quantities(vol_id=0, edge=1)
    1042         max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux2, epsilon, g, H0)       
    1043 
    1044         #Scale, add up and check that compute_fluxes is correct for vol 1
    1045         e0 = domain.edgelengths[1, 0]
    1046         e1 = domain.edgelengths[1, 1]
    1047         e2 = domain.edgelengths[1, 2]
    1048 
    1049         total_flux = -(e0*edgeflux0+e1*edgeflux1+e2*edgeflux2)/domain.areas[1]
    1050 
    1051         domain.compute_fluxes()
    1052         for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
    1053             assert num.allclose(total_flux[i],
    1054                                 domain.quantities[name].explicit_update[1])
    1055 
    1056 
    1057 
    1058     def xtest_catching_negative_heights(self):
    1059 
    1060         #OBSOLETE
    1061        
    1062         a = [0.0, 0.0]
    1063         b = [0.0, 2.0]
    1064         c = [2.0,0.0]
    1065         d = [0.0, 4.0]
    1066         e = [2.0, 2.0]
    1067         f = [4.0,0.0]
    1068 
    1069         points = [a, b, c, d, e, f]
    1070         #bac, bce, ecf, dbe
    1071         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    1072 
    1073         domain = Domain(points, vertices)
    1074         val0 = 2.+2.0/3
    1075         val1 = 4.+4.0/3
    1076         val2 = 8.+2.0/3
    1077         val3 = 2.+8.0/3
    1078 
    1079         zl=zr=4  #Too large
    1080         domain.set_quantity('elevation', zl*num.ones( (4,3), num.int )) #array default#
     1050        val0 = 2. + 2.0/3
     1051        val1 = 4. + 4.0/3
     1052        val2 = 8. + 2.0/3
     1053        val3 = 2. + 8.0/3
     1054
     1055        zl = zr = 4    # Too large
     1056        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
    10811057        domain.set_quantity('stage', [[val0, val0-1, val0-2],
    10821058                                      [val1, val1+1, val1],
     
    10901066            pass
    10911067
    1092 
    1093 
    10941068    def test_get_wet_elements(self):
    1095 
    10961069        a = [0.0, 0.0]
    10971070        b = [0.0, 2.0]
    1098         c = [2.0,0.0]
     1071        c = [2.0, 0.0]
    10991072        d = [0.0, 4.0]
    11001073        e = [2.0, 2.0]
    1101         f = [4.0,0.0]
     1074        f = [4.0, 0.0]
    11021075
    11031076        points = [a, b, c, d, e, f]
    1104         #bac, bce, ecf, dbe
    1105         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     1077        #             bac,     bce,     ecf,    dbe
     1078        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    11061079
    11071080        domain = Domain(points, vertices)
    1108         val0 = 2.+2.0/3
    1109         val1 = 4.+4.0/3
    1110         val2 = 8.+2.0/3
    1111         val3 = 2.+8.0/3
    1112 
    1113         zl=zr=5
    1114         domain.set_quantity('elevation', zl*num.ones( (4,3), num.int )) #array default#
     1081
     1082        val0 = 2. + 2.0/3
     1083        val1 = 4. + 4.0/3
     1084        val2 = 8. + 2.0/3
     1085        val3 = 2. + 8.0/3
     1086
     1087        zl = zr = 5
     1088        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
    11151089        domain.set_quantity('stage', [[val0, val0-1, val0-2],
    11161090                                      [val1, val1+1, val1],
     
    11181092                                      [val3-0.5, val3, val3]])
    11191093
    1120 
    1121 
    1122         #print domain.get_quantity('elevation').get_values(location='centroids')
    1123         #print domain.get_quantity('stage').get_values(location='centroids')       
    11241094        domain.check_integrity()
    11251095
    11261096        indices = domain.get_wet_elements()
    1127         assert num.allclose(indices, [1,2])
    1128 
    1129         indices = domain.get_wet_elements(indices=[0,1,3])
     1097        assert num.allclose(indices, [1, 2])
     1098
     1099        indices = domain.get_wet_elements(indices=[0, 1, 3])
    11301100        assert num.allclose(indices, [1])
    1131        
    1132 
    11331101
    11341102    def test_get_maximum_inundation_1(self):
    1135 
    11361103        a = [0.0, 0.0]
    11371104        b = [0.0, 2.0]
    1138         c = [2.0,0.0]
     1105        c = [2.0, 0.0]
    11391106        d = [0.0, 4.0]
    11401107        e = [2.0, 2.0]
    1141         f = [4.0,0.0]
     1108        f = [4.0, 0.0]
    11421109
    11431110        points = [a, b, c, d, e, f]
    1144         #bac, bce, ecf, dbe
    1145         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     1111        #             bac,     bce,     ecf,    dbe
     1112        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    11461113
    11471114        domain = Domain(points, vertices)
    11481115
    1149         domain.set_quantity('elevation', lambda x, y: x+2*y) #2 4 4 6
     1116        domain.set_quantity('elevation', lambda x, y: x+2*y)    # 2 4 4 6
    11501117        domain.set_quantity('stage', 3)
    11511118
     
    11561123
    11571124        q = domain.get_maximum_inundation_elevation()
    1158         assert num.allclose(q, domain.get_quantity('elevation').get_values(location='centroids')[0])
     1125        assert num.allclose(q, domain.get_quantity('elevation').\
     1126                                   get_values(location='centroids')[0])
    11591127
    11601128        x, y = domain.get_maximum_inundation_location()
    11611129        assert num.allclose([x, y], domain.get_centroid_coordinates()[0])
    11621130
    1163 
    11641131    def test_get_maximum_inundation_2(self):
    11651132        """test_get_maximum_inundation_2(self)
     
    11671134        Test multiple wet cells with same elevation
    11681135        """
    1169        
     1136
    11701137        a = [0.0, 0.0]
    11711138        b = [0.0, 2.0]
    1172         c = [2.0,0.0]
     1139        c = [2.0, 0.0]
    11731140        d = [0.0, 4.0]
    11741141        e = [2.0, 2.0]
    1175         f = [4.0,0.0]
     1142        f = [4.0, 0.0]
    11761143
    11771144        points = [a, b, c, d, e, f]
    1178         #bac, bce, ecf, dbe
    1179         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     1145        #             bac,     bce,     ecf,    dbe
     1146        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    11801147
    11811148        domain = Domain(points, vertices)
    11821149
    1183         domain.set_quantity('elevation', lambda x, y: x+2*y) #2 4 4 6
     1150        domain.set_quantity('elevation', lambda x, y: x+2*y)    # 2 4 4 6
    11841151        domain.set_quantity('stage', 4.1)
    11851152
     
    11871154
    11881155        indices = domain.get_wet_elements()
    1189         assert num.allclose(indices, [0,1,2])
     1156        assert num.allclose(indices, [0, 1, 2])
    11901157
    11911158        q = domain.get_maximum_inundation_elevation()
    1192         assert num.allclose(q, 4)       
     1159        assert num.allclose(q, 4)
    11931160
    11941161        x, y = domain.get_maximum_inundation_location()
    1195         assert num.allclose([x, y], domain.get_centroid_coordinates()[1])       
    1196        
     1162        assert num.allclose([x, y], domain.get_centroid_coordinates()[1])
    11971163
    11981164    def test_get_maximum_inundation_3(self):
     
    12021168        """
    12031169
    1204         from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
     1170        from anuga.abstract_2d_finite_volumes.mesh_factory \
     1171                import rectangular_cross
    12051172
    12061173        initial_runup_height = -0.4
    12071174        final_runup_height = -0.3
    1208 
    12091175
    12101176        #--------------------------------------------------------------
     
    12121178        #--------------------------------------------------------------
    12131179        N = 5
    1214         points, vertices, boundary = rectangular_cross(N, N) 
     1180        points, vertices, boundary = rectangular_cross(N, N)
    12151181        domain = Domain(points, vertices, boundary)
    1216         domain.set_maximum_allowed_speed(1.0)       
     1182        domain.set_maximum_allowed_speed(1.0)
    12171183
    12181184        #--------------------------------------------------------------
    12191185        # Setup initial conditions
    12201186        #--------------------------------------------------------------
    1221         def topography(x,y):
     1187        def topography(x, y):
    12221188            return -x/2                             # linear bed slope
    1223            
    1224 
    1225         domain.set_quantity('elevation', topography)       # Use function for elevation
    1226         domain.set_quantity('friction', 0.)                # Zero friction 
    1227         domain.set_quantity('stage', initial_runup_height) # Constant negative initial stage
    1228 
     1189
     1190        # Use function for elevation
     1191        domain.set_quantity('elevation', topography)
     1192        domain.set_quantity('friction', 0.)                # Zero friction
     1193        # Constant negative initial stage
     1194        domain.set_quantity('stage', initial_runup_height)
    12291195
    12301196        #--------------------------------------------------------------
    12311197        # Setup boundary conditions
    12321198        #--------------------------------------------------------------
    1233         Br = Reflective_boundary(domain)              # Reflective wall
    1234         Bd = Dirichlet_boundary([final_runup_height,  # Constant inflow
    1235                                  0,
    1236                                  0])
    1237 
    1238         # All reflective to begin with (still water)
     1199        Br = Reflective_boundary(domain)                       # Reflective wall
     1200        Bd = Dirichlet_boundary([final_runup_height, 0, 0])    # Constant inflow
     1201
     1202        # All reflective to begin with (still water)
    12391203        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
    1240 
    12411204
    12421205        #--------------------------------------------------------------
     
    12451208
    12461209        indices = domain.get_wet_elements()
    1247         z = domain.get_quantity('elevation').\
    1248             get_values(location='centroids', indices=indices)
    1249         assert num.alltrue(z<initial_runup_height)
     1210        z = domain.get_quantity('elevation').get_values(location='centroids',
     1211                                                        indices=indices)
     1212        assert num.alltrue(z < initial_runup_height)
    12501213
    12511214        q = domain.get_maximum_inundation_elevation()
    1252         assert num.allclose(q, initial_runup_height, rtol = 1.0/N) # First order accuracy
     1215        # First order accuracy
     1216        assert num.allclose(q, initial_runup_height, rtol=1.0/N)
    12531217
    12541218        x, y = domain.get_maximum_inundation_location()
    12551219
    1256         qref = domain.get_quantity('elevation').get_values(interpolation_points = [[x, y]])
     1220        qref = domain.get_quantity('elevation').\
     1221                     get_values(interpolation_points=[[x, y]])
    12571222        assert num.allclose(q, qref)
    12581223
    1259 
    12601224        wet_elements = domain.get_wet_elements()
    1261         wet_elevations = domain.get_quantity('elevation').get_values(location='centroids',
    1262                                                                      indices=wet_elements)
    1263         assert num.alltrue(wet_elevations<initial_runup_height)
    1264         assert num.allclose(wet_elevations, z)       
    1265 
    1266 
    1267         #print domain.get_quantity('elevation').get_maximum_value(indices=wet_elements)
    1268         #print domain.get_quantity('elevation').get_maximum_location(indices=wet_elements)
    1269         #print domain.get_quantity('elevation').get_maximum_index(indices=wet_elements)
    1270 
    1271        
     1225        wet_elevations = domain.get_quantity('elevation').\
     1226                                    get_values(location='centroids',
     1227                                               indices=wet_elements)
     1228        assert num.alltrue(wet_elevations < initial_runup_height)
     1229        assert num.allclose(wet_elevations, z)
     1230
    12721231        #--------------------------------------------------------------
    12731232        # Let triangles adjust
     
    12761235            pass
    12771236
    1278 
    12791237        #--------------------------------------------------------------
    12801238        # Test inundation height again
    12811239        #--------------------------------------------------------------
    1282 
    12831240        indices = domain.get_wet_elements()
    1284         z = domain.get_quantity('elevation').\
    1285             get_values(location='centroids', indices=indices)
    1286 
    1287         assert num.alltrue(z<initial_runup_height)
     1241        z = domain.get_quantity('elevation').get_values(location='centroids',
     1242                                                        indices=indices)
     1243
     1244        assert num.alltrue(z < initial_runup_height)
    12881245
    12891246        q = domain.get_maximum_inundation_elevation()
    1290         assert num.allclose(q, initial_runup_height, rtol = 1.0/N) # First order accuracy
    1291        
     1247        # First order accuracy
     1248        assert num.allclose(q, initial_runup_height, rtol=1.0/N)
     1249
    12921250        x, y = domain.get_maximum_inundation_location()
    1293         qref = domain.get_quantity('elevation').get_values(interpolation_points = [[x, y]])
    1294         assert num.allclose(q, qref)       
    1295 
     1251        qref = domain.get_quantity('elevation').\
     1252                        get_values(interpolation_points=[[x, y]])
     1253        assert num.allclose(q, qref)
    12961254
    12971255        #--------------------------------------------------------------
     
    13001258        domain.set_boundary({'right': Bd})
    13011259
    1302        
    13031260        #--------------------------------------------------------------
    13041261        # Evolve system through time
    13051262        #--------------------------------------------------------------
    13061263        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0):
    1307             #print domain.timestepping_statistics(track_speeds=True)
    1308             #domain.write_time()
    13091264            pass
    1310    
     1265
    13111266        #--------------------------------------------------------------
    13121267        # Test inundation height again
    13131268        #--------------------------------------------------------------
    1314 
    13151269        indices = domain.get_wet_elements()
    13161270        z = domain.get_quantity('elevation').\
    1317             get_values(location='centroids', indices=indices)
    1318 
    1319         assert num.alltrue(z<final_runup_height)
     1271                    get_values(location='centroids', indices=indices)
     1272
     1273        assert num.alltrue(z < final_runup_height)
    13201274
    13211275        q = domain.get_maximum_inundation_elevation()
    1322         assert num.allclose(q, final_runup_height, rtol = 1.0/N) # First order accuracy
     1276        # First order accuracy
     1277        assert num.allclose(q, final_runup_height, rtol=1.0/N)
    13231278
    13241279        x, y = domain.get_maximum_inundation_location()
    1325         qref = domain.get_quantity('elevation').get_values(interpolation_points = [[x, y]])
    1326         assert num.allclose(q, qref)       
    1327 
     1280        qref = domain.get_quantity('elevation').\
     1281                        get_values(interpolation_points=[[x, y]])
     1282        assert num.allclose(q, qref)
    13281283
    13291284        wet_elements = domain.get_wet_elements()
    1330         wet_elevations = domain.get_quantity('elevation').get_values(location='centroids',
    1331                                                                      indices=wet_elements)
    1332         assert num.alltrue(wet_elevations<final_runup_height)
    1333         assert num.allclose(wet_elevations, z)       
    1334        
    1335 
     1285        wet_elevations = domain.get_quantity('elevation').\
     1286                             get_values(location='centroids',
     1287                                        indices=wet_elements)
     1288        assert num.alltrue(wet_elevations < final_runup_height)
     1289        assert num.allclose(wet_elevations, z)
    13361290
    13371291    def test_get_maximum_inundation_from_sww(self):
     
    13401294        Test of get_maximum_inundation_elevation()
    13411295        and get_maximum_inundation_location() from data_manager.py
    1342        
     1296
    13431297        This is based on test_get_maximum_inundation_3(self) but works with the
    13441298        stored results instead of with the internal data structure.
     
    13471301        """
    13481302
    1349         from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
     1303        from anuga.abstract_2d_finite_volumes.mesh_factory \
     1304                import rectangular_cross
    13501305        from data_manager import get_maximum_inundation_elevation
    13511306        from data_manager import get_maximum_inundation_location
    13521307        from data_manager import get_maximum_inundation_data
    1353        
    13541308
    13551309        initial_runup_height = -0.4
    13561310        final_runup_height = -0.3
    1357 
    13581311
    13591312        #--------------------------------------------------------------
     
    13611314        #--------------------------------------------------------------
    13621315        N = 10
    1363         points, vertices, boundary = rectangular_cross(N, N) 
     1316        points, vertices, boundary = rectangular_cross(N, N)
    13641317        domain = Domain(points, vertices, boundary)
    13651318        domain.set_name('runup_test')
    13661319        domain.set_maximum_allowed_speed(1.0)
    13671320
    1368         domain.tight_slope_limiters = 0 # FIXME: This works better with old limiters so far
     1321        # FIXME: This works better with old limiters so far
     1322        domain.tight_slope_limiters = 0
    13691323
    13701324        #--------------------------------------------------------------
    13711325        # Setup initial conditions
    13721326        #--------------------------------------------------------------
    1373         def topography(x,y):
     1327        def topography(x, y):
    13741328            return -x/2                             # linear bed slope
    1375            
    1376 
    1377         domain.set_quantity('elevation', topography)       # Use function for elevation
    1378         domain.set_quantity('friction', 0.)                # Zero friction 
    1379         domain.set_quantity('stage', initial_runup_height) # Constant negative initial stage
    1380 
     1329
     1330        # Use function for elevation
     1331        domain.set_quantity('elevation', topography)
     1332        domain.set_quantity('friction', 0.)                # Zero friction
     1333        # Constant negative initial stage
     1334        domain.set_quantity('stage', initial_runup_height)
    13811335
    13821336        #--------------------------------------------------------------
    13831337        # Setup boundary conditions
    13841338        #--------------------------------------------------------------
    1385         Br = Reflective_boundary(domain)              # Reflective wall
    1386         Bd = Dirichlet_boundary([final_runup_height,  # Constant inflow
    1387                                  0,
    1388                                  0])
    1389 
    1390         # All reflective to begin with (still water)
     1339        Br = Reflective_boundary(domain)                       # Reflective wall
     1340        Bd = Dirichlet_boundary([final_runup_height, 0, 0])    # Constant inflow
     1341
     1342        # All reflective to begin with (still water)
    13911343        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
    1392 
    13931344
    13941345        #--------------------------------------------------------------
    13951346        # Test initial inundation height
    13961347        #--------------------------------------------------------------
    1397 
    13981348        indices = domain.get_wet_elements()
    13991349        z = domain.get_quantity('elevation').\
    1400             get_values(location='centroids', indices=indices)
    1401         assert num.alltrue(z<initial_runup_height)
     1350                get_values(location='centroids', indices=indices)
     1351        assert num.alltrue(z < initial_runup_height)
    14021352
    14031353        q_ref = domain.get_maximum_inundation_elevation()
    1404         assert num.allclose(q_ref, initial_runup_height, rtol = 1.0/N) # First order accuracy
    1405 
    1406        
     1354        # First order accuracy
     1355        assert num.allclose(q_ref, initial_runup_height, rtol=1.0/N)
     1356
    14071357        #--------------------------------------------------------------
    14081358        # Let triangles adjust
     
    14111361            pass
    14121362
    1413 
    14141363        #--------------------------------------------------------------
    14151364        # Test inundation height again
    14161365        #--------------------------------------------------------------
    1417        
    14181366        q_ref = domain.get_maximum_inundation_elevation()
    14191367        q = get_maximum_inundation_elevation('runup_test.sww')
    1420         msg = 'We got %f, should have been %f' %(q, q_ref)
     1368        msg = 'We got %f, should have been %f' % (q, q_ref)
    14211369        assert num.allclose(q, q_ref, rtol=1.0/N), msg
    14221370
    1423 
    14241371        q = get_maximum_inundation_elevation('runup_test.sww')
    1425         msg = 'We got %f, should have been %f' %(q, initial_runup_height)
    1426         assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
    1427 
     1372        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
     1373        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
    14281374
    14291375        # Test error condition if time interval is out
     
    14391385        # Check correct time interval
    14401386        q, loc = get_maximum_inundation_data('runup_test.sww',
    1441                                              time_interval=[0.0, 3.0])       
    1442         msg = 'We got %f, should have been %f' %(q, initial_runup_height)
     1387                                             time_interval=[0.0, 3.0])
     1388        msg = 'We got %f, should have been %f' % (q, initial_runup_height)
    14431389        assert num.allclose(q, initial_runup_height, rtol = 1.0/N), msg
    1444         assert num.allclose(-loc[0]/2, q) # From topography formula         
    1445        
     1390        assert num.allclose(-loc[0]/2, q)    # From topography formula
    14461391
    14471392        #--------------------------------------------------------------
     
    14501395        domain.set_boundary({'right': Bd})
    14511396
    1452        
    14531397        #--------------------------------------------------------------
    14541398        # Evolve system through time
     
    14561400        q_max = None
    14571401        for t in domain.evolve(yieldstep = 0.1, finaltime = 3.0,
    1458                                skip_initial_step = True): 
     1402                               skip_initial_step = True):
    14591403            q = domain.get_maximum_inundation_elevation()
    1460             if q > q_max: q_max = q
    1461 
    1462    
     1404            if q > q_max:
     1405                q_max = q
     1406
    14631407        #--------------------------------------------------------------
    14641408        # Test inundation height again
    14651409        #--------------------------------------------------------------
    1466 
    14671410        indices = domain.get_wet_elements()
    14681411        z = domain.get_quantity('elevation').\
    1469             get_values(location='centroids', indices=indices)
    1470 
    1471         assert num.alltrue(z<final_runup_height)
     1412                get_values(location='centroids', indices=indices)
     1413
     1414        assert num.alltrue(z < final_runup_height)
    14721415
    14731416        q = domain.get_maximum_inundation_elevation()
    1474         assert num.allclose(q, final_runup_height, rtol = 1.0/N) # First order accuracy
    1475 
    1476         q, loc = get_maximum_inundation_data('runup_test.sww', time_interval=[3.0, 3.0])
    1477         msg = 'We got %f, should have been %f' %(q, final_runup_height)
     1417        # First order accuracy
     1418        assert num.allclose(q, final_runup_height, rtol=1.0/N)
     1419
     1420        q, loc = get_maximum_inundation_data('runup_test.sww',
     1421                                             time_interval=[3.0, 3.0])
     1422        msg = 'We got %f, should have been %f' % (q, final_runup_height)
    14781423        assert num.allclose(q, final_runup_height, rtol=1.0/N), msg
    1479         #print 'loc', loc, q       
    1480         assert num.allclose(-loc[0]/2, q) # From topography formula         
     1424        assert num.allclose(-loc[0]/2, q)    # From topography formula
    14811425
    14821426        q = get_maximum_inundation_elevation('runup_test.sww')
    1483         loc = get_maximum_inundation_location('runup_test.sww')       
    1484         msg = 'We got %f, should have been %f' %(q, q_max)
     1427        loc = get_maximum_inundation_location('runup_test.sww')
     1428        msg = 'We got %f, should have been %f' % (q, q_max)
    14851429        assert num.allclose(q, q_max, rtol=1.0/N), msg
    1486         #print 'loc', loc, q
    1487         assert num.allclose(-loc[0]/2, q) # From topography formula
    1488 
    1489        
    1490 
    1491         q = get_maximum_inundation_elevation('runup_test.sww', time_interval=[0, 3])
    1492         msg = 'We got %f, should have been %f' %(q, q_max)
     1430        assert num.allclose(-loc[0]/2, q)    # From topography formula
     1431
     1432        q = get_maximum_inundation_elevation('runup_test.sww',
     1433                                             time_interval=[0, 3])
     1434        msg = 'We got %f, should have been %f' % (q, q_max)
    14931435        assert num.allclose(q, q_max, rtol=1.0/N), msg
    14941436
    1495 
    14961437        # Check polygon mode
    1497         polygon = [[0.3, 0.0], [0.9, 0.0], [0.9, 1.0], [0.3, 1.0]] # Runup region
     1438        # Runup region
     1439        polygon = [[0.3, 0.0], [0.9, 0.0], [0.9, 1.0], [0.3, 1.0]]
    14981440        q = get_maximum_inundation_elevation('runup_test.sww',
    14991441                                             polygon = polygon,
    15001442                                             time_interval=[0, 3])
    1501         msg = 'We got %f, should have been %f' %(q, q_max)
     1443        msg = 'We got %f, should have been %f' % (q, q_max)
    15021444        assert num.allclose(q, q_max, rtol=1.0/N), msg
    15031445
    1504        
    1505         polygon = [[0.9, 0.0], [1.0, 0.0], [1.0, 1.0], [0.9, 1.0]] # Offshore region
     1446        # Offshore region
     1447        polygon = [[0.9, 0.0], [1.0, 0.0], [1.0, 1.0], [0.9, 1.0]]
    15061448        q, loc = get_maximum_inundation_data('runup_test.sww',
    15071449                                             polygon = polygon,
    15081450                                             time_interval=[0, 3])
    1509         msg = 'We got %f, should have been %f' %(q, -0.475)
     1451        msg = 'We got %f, should have been %f' % (q, -0.475)
    15101452        assert num.allclose(q, -0.475, rtol=1.0/N), msg
    15111453        assert is_inside_polygon(loc, polygon)
    1512         assert num.allclose(-loc[0]/2, q) # From topography formula         
    1513 
    1514 
    1515         polygon = [[0.0, 0.0], [0.4, 0.0], [0.4, 1.0], [0.0, 1.0]] # Dry region
     1454        assert num.allclose(-loc[0]/2, q)    # From topography formula
     1455
     1456        # Dry region
     1457        polygon = [[0.0, 0.0], [0.4, 0.0], [0.4, 1.0], [0.0, 1.0]]
    15161458        q, loc = get_maximum_inundation_data('runup_test.sww',
    15171459                                             polygon = polygon,
    15181460                                             time_interval=[0, 3])
    1519         msg = 'We got %s, should have been None' %(q)
     1461        msg = 'We got %s, should have been None' % (q)
    15201462        assert q is None, msg
    1521         msg = 'We got %s, should have been None' %(loc)
    1522         assert loc is None, msg       
     1463        msg = 'We got %s, should have been None' % (loc)
     1464        assert loc is None, msg
    15231465
    15241466        # Check what happens if no time point is within interval
    15251467        try:
    1526             q = get_maximum_inundation_elevation('runup_test.sww', time_interval=[2.75, 2.75])
     1468            q = get_maximum_inundation_elevation('runup_test.sww',
     1469                                                 time_interval=[2.75, 2.75])
    15271470        except AssertionError:
    15281471            pass
    15291472        else:
    15301473            msg = 'Time interval should have raised an exception'
    1531             raise msg
     1474            raise Exception, msg
    15321475
    15331476        # Cleanup
     
    15371480            pass
    15381481            #FIXME(Ole): Windows won't allow removal of this
    1539            
    1540        
    1541        
     1482
    15421483    def test_get_flow_through_cross_section_with_geo(self):
    15431484        """test_get_flow_through_cross_section(self):
     
    15451486        Test that the total flow through a cross section can be
    15461487        correctly obtained at run-time from the ANUGA domain.
    1547        
     1488
    15481489        This test creates a flat bed with a known flow through it and tests
    15491490        that the function correctly returns the expected flow.
     
    15571498        q = u*h*w = 12 m^3/s
    15581499
    1559 
    15601500        This run tries it with georeferencing and with elevation = -1
    1561        
    15621501        """
    15631502
    15641503        import time, os
    15651504        from Scientific.IO.NetCDF import NetCDFFile
    1566 
    1567         # Setup
    15681505        from mesh_factory import rectangular
    15691506
     
    15721509        length = 20
    15731510        t_end = 1
    1574         points, vertices, boundary = rectangular(length, width,
    1575                                                  length, width)
     1511        points, vertices, boundary = rectangular(length, width, length, width)
    15761512
    15771513        # Create shallow water domain
    15781514        domain = Domain(points, vertices, boundary,
    1579                         geo_reference=Geo_reference(56,308500,6189000))
     1515                        geo_reference=Geo_reference(56, 308500, 6189000))
    15801516
    15811517        domain.default_order = 2
    15821518        domain.set_quantities_to_be_stored(None)
    1583 
    15841519
    15851520        e = -1.0
     
    15891524        uh = u*h
    15901525
    1591         Br = Reflective_boundary(domain)     # Side walls
    1592         Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet:
    1593 
     1526        Br = Reflective_boundary(domain)       # Side walls
     1527        Bd = Dirichlet_boundary([w, uh, 0])    # 2 m/s across the 3 m inlet:
    15941528
    15951529        # Initial conditions
     
    15971531        domain.set_quantity('stage', w)
    15981532        domain.set_quantity('xmomentum', uh)
    1599         domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
    1600        
    1601        
     1533        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
     1534
    16021535        # Interpolation points down the middle
    16031536        I = [[0, width/2.],
    16041537             [length/2., width/2.],
    16051538             [length, width/2.]]
    1606         interpolation_points = domain.geo_reference.get_absolute(I)       
    1607        
     1539        interpolation_points = domain.geo_reference.get_absolute(I)
     1540
    16081541        # Shortcuts to quantites
    1609         stage = domain.get_quantity('stage')       
    1610         xmomentum = domain.get_quantity('xmomentum')       
    1611         ymomentum = domain.get_quantity('ymomentum')               
    1612 
    1613         for t in domain.evolve(yieldstep=0.1, finaltime = t_end):
     1542        stage = domain.get_quantity('stage')
     1543        xmomentum = domain.get_quantity('xmomentum')
     1544        ymomentum = domain.get_quantity('ymomentum')
     1545
     1546        for t in domain.evolve(yieldstep=0.1, finaltime=t_end):
    16141547            # Check that quantities are they should be in the interior
    1615 
    1616             w_t = stage.get_values(interpolation_points)           
     1548            w_t = stage.get_values(interpolation_points)
    16171549            uh_t = xmomentum.get_values(interpolation_points)
    1618             vh_t = ymomentum.get_values(interpolation_points)           
    1619            
     1550            vh_t = ymomentum.get_values(interpolation_points)
     1551
    16201552            assert num.allclose(w_t, w)
    1621             assert num.allclose(uh_t, uh)           
    1622             assert num.allclose(vh_t, 0.0)                       
    1623            
    1624            
     1553            assert num.allclose(uh_t, uh)
     1554            assert num.allclose(vh_t, 0.0)
     1555
    16251556            # Check flows through the middle
    16261557            for i in range(5):
    1627                 x = length/2. + i*0.23674563 # Arbitrary
     1558                x = length/2. + i*0.23674563    # Arbitrary
    16281559                cross_section = [[x, 0], [x, width]]
    16291560
    1630                 cross_section = domain.geo_reference.get_absolute(cross_section)           
     1561                cross_section = domain.geo_reference.get_absolute(cross_section)
    16311562                Q = domain.get_flow_through_cross_section(cross_section,
    16321563                                                          verbose=False)
     
    16341565                assert num.allclose(Q, uh*width)
    16351566
    1636 
    1637        
    16381567    def test_get_energy_through_cross_section_with_geo(self):
    16391568        """test_get_energy_through_cross_section(self):
     
    16411570        Test that the total and specific energy through a cross section can be
    16421571        correctly obtained at run-time from the ANUGA domain.
    1643        
     1572
    16441573        This test creates a flat bed with a known flow through it and tests
    16451574        that the function correctly returns the expected energies.
     
    16531582        q = u*h*w = 12 m^3/s
    16541583
    1655 
    16561584        This run tries it with georeferencing and with elevation = -1
    1657        
    16581585        """
    16591586
    16601587        import time, os
    16611588        from Scientific.IO.NetCDF import NetCDFFile
    1662 
    1663         # Setup
    16641589        from mesh_factory import rectangular
    16651590
     
    16681593        length = 20
    16691594        t_end = 1
    1670         points, vertices, boundary = rectangular(length, width,
    1671                                                  length, width)
     1595        points, vertices, boundary = rectangular(length, width, length, width)
    16721596
    16731597        # Create shallow water domain
    16741598        domain = Domain(points, vertices, boundary,
    1675                         geo_reference=Geo_reference(56,308500,6189000))
     1599                        geo_reference=Geo_reference(56, 308500, 6189000))
    16761600
    16771601        domain.default_order = 2
    16781602        domain.set_quantities_to_be_stored(None)
    1679 
    16801603
    16811604        e = -1.0
     
    16851608        uh = u*h
    16861609
    1687         Br = Reflective_boundary(domain)     # Side walls
    1688         Bd = Dirichlet_boundary([w, uh, 0])  # 2 m/s across the 3 m inlet:
    1689 
     1610        Br = Reflective_boundary(domain)       # Side walls
     1611        Bd = Dirichlet_boundary([w, uh, 0])    # 2 m/s across the 3 m inlet:
    16901612
    16911613        # Initial conditions
     
    16931615        domain.set_quantity('stage', w)
    16941616        domain.set_quantity('xmomentum', uh)
    1695         domain.set_boundary( {'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
    1696        
    1697        
     1617        domain.set_boundary({'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
     1618
    16981619        # Interpolation points down the middle
    16991620        I = [[0, width/2.],
    17001621             [length/2., width/2.],
    17011622             [length, width/2.]]
    1702         interpolation_points = domain.geo_reference.get_absolute(I)       
    1703        
     1623        interpolation_points = domain.geo_reference.get_absolute(I)
     1624
    17041625        # Shortcuts to quantites
    1705         stage = domain.get_quantity('stage')       
    1706         xmomentum = domain.get_quantity('xmomentum')       
    1707         ymomentum = domain.get_quantity('ymomentum')               
    1708 
    1709         for t in domain.evolve(yieldstep=0.1, finaltime = t_end):
     1626        stage = domain.get_quantity('stage')
     1627        xmomentum = domain.get_quantity('xmomentum')
     1628        ymomentum = domain.get_quantity('ymomentum')
     1629
     1630        for t in domain.evolve(yieldstep=0.1, finaltime=t_end):
    17101631            # Check that quantities are they should be in the interior
    1711 
    1712             w_t = stage.get_values(interpolation_points)           
     1632            w_t = stage.get_values(interpolation_points)
    17131633            uh_t = xmomentum.get_values(interpolation_points)
    1714             vh_t = ymomentum.get_values(interpolation_points)           
    1715            
     1634            vh_t = ymomentum.get_values(interpolation_points)
     1635
    17161636            assert num.allclose(w_t, w)
    1717             assert num.allclose(uh_t, uh)           
    1718             assert num.allclose(vh_t, 0.0)                       
    1719            
    1720            
     1637            assert num.allclose(uh_t, uh)
     1638            assert num.allclose(vh_t, 0.0)
     1639
    17211640            # Check energies through the middle
    17221641            for i in range(5):
    1723                 x = length/2. + i*0.23674563 # Arbitrary
     1642                x = length/2. + i*0.23674563    # Arbitrary
    17241643                cross_section = [[x, 0], [x, width]]
    17251644
    1726                 cross_section = domain.geo_reference.get_absolute(cross_section)   
     1645                cross_section = domain.geo_reference.get_absolute(cross_section)
    17271646                Es = domain.get_energy_through_cross_section(cross_section,
    17281647                                                             kind='specific',
    17291648                                                             verbose=False)
    1730                                                      
     1649
    17311650                assert num.allclose(Es, h + 0.5*u*u/g)
    1732            
     1651
    17331652                Et = domain.get_energy_through_cross_section(cross_section,
    17341653                                                             kind='total',
    17351654                                                             verbose=False)
    1736                 assert num.allclose(Et, w + 0.5*u*u/g)           
    1737 
    1738            
    1739        
    1740        
     1655                assert num.allclose(Et, w + 0.5*u*u/g)
    17411656
    17421657    def test_another_runup_example(self):
     
    17471662        """
    17481663
    1749         #-----------------------------------------------------------------
    1750         # Import necessary modules
    1751         #-----------------------------------------------------------------
    1752 
    17531664        from anuga.pmesh.mesh_interface import create_mesh_from_regions
    1754         from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
     1665        from anuga.abstract_2d_finite_volumes.mesh_factory \
     1666                import rectangular_cross
    17551667        from anuga.shallow_water import Domain
    17561668        from anuga.shallow_water import Reflective_boundary
    17571669        from anuga.shallow_water import Dirichlet_boundary
    17581670
    1759 
    17601671        #-----------------------------------------------------------------
    17611672        # Setup computational domain
    17621673        #-----------------------------------------------------------------
    1763         points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh
    1764         domain = Domain(points, vertices, boundary) # Create domain
     1674        points, vertices, boundary = rectangular_cross(10, 10)    # Basic mesh
     1675        domain = Domain(points, vertices, boundary)    # Create domain
    17651676        domain.set_quantities_to_be_stored(None)
    17661677        domain.set_maximum_allowed_speed(100) #FIXME (Ole): try to remove this
    1767        
     1678
    17681679        # FIXME (Ole): Need tests where this is commented out
    1769         domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)
    1770         domain.H0 = 0 # Backwards compatibility (6/2/7)
    1771         domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
    1772        
     1680        domain.tight_slope_limiters = 0    # Backwards compatibility (14/4/7)
     1681        domain.H0 = 0    # Backwards compatibility (6/2/7)
     1682        domain.use_centroid_velocities = 0    # Backwards compatibility (7/5/8)
    17731683
    17741684        #-----------------------------------------------------------------
    17751685        # Setup initial conditions
    17761686        #-----------------------------------------------------------------
    1777 
    1778         def topography(x,y):
     1687        def topography(x, y):
    17791688            return -x/2                              # linear bed slope
    17801689
    1781         domain.set_quantity('elevation', topography)
    1782         domain.set_quantity('friction', 0.0)         
    1783         domain.set_quantity('stage', expression='elevation')           
    1784 
     1690        domain.set_quantity('elevation', topography)
     1691        domain.set_quantity('friction', 0.0)
     1692        domain.set_quantity('stage', expression='elevation')
    17851693
    17861694        #----------------------------------------------------------------
    17871695        # Setup boundary conditions
    17881696        #----------------------------------------------------------------
    1789 
    1790         Br = Reflective_boundary(domain)      # Solid reflective wall
    1791         Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values
     1697        Br = Reflective_boundary(domain)           # Solid reflective wall
     1698        Bd = Dirichlet_boundary([-0.2, 0., 0.])    # Constant boundary values
    17921699        domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom': Br})
    1793 
    17941700
    17951701        #----------------------------------------------------------------
    17961702        # Evolve system through time
    17971703        #----------------------------------------------------------------
    1798 
    17991704        interpolation_points = [[0.4,0.5], [0.6,0.5], [0.8,0.5], [0.9,0.5]]
    18001705        gauge_values = []
     
    18031708
    18041709        time = []
    1805         for t in domain.evolve(yieldstep = 0.1, finaltime = 5.0):
     1710        for t in domain.evolve(yieldstep=0.1, finaltime=5.0):
    18061711            # Record time series at known points
    18071712            time.append(domain.get_time())
    1808            
     1713
    18091714            stage = domain.get_quantity('stage')
    18101715            w = stage.get_values(interpolation_points=interpolation_points)
    1811            
     1716
    18121717            for i, _ in enumerate(interpolation_points):
    18131718                gauge_values[i].append(w[i])
    18141719
    1815 
    1816         #print
    1817         #print time
    1818         #print
    1819         #for i, (x,y) in enumerate(interpolation_points):
    1820         #    print i, gauge_values[i]
    1821         #    print
    1822 
    18231720        #Reference (nautilus 26/6/2008)
    1824        
    1825         G0 = [-0.20000000000000001, -0.20000000000000001, -0.19920600846161715, -0.19153647344085376, -0.19127622768281194, -0.1770671909675095, -0.16739412133181927, -0.16196038919122191, -0.15621633053131384, -0.15130021599977705, -0.13930978857215484, -0.19349274358263582, -0.19975307598803765, -0.19999897143103357, -0.1999999995532111, -0.19999999999949952, -0.19999999999949952, -0.19999999999949952, -0.19997270012494556, -0.19925805948554556, -0.19934513778450533, -0.19966484196394893, -0.1997352860102084, -0.19968260481750394, -0.19980280797303882, -0.19998804881822749, -0.19999999778075916, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167, -0.19999999999966167]
    1826 
    1827         G1 = [-0.29999999999999993, -0.29999588068034899, -0.29250047332330331, -0.28335081844518584, -0.26142206997410805, -0.22656028856329835, -0.21224087216745585, -0.19934324109114465, -0.1889857939783175, -0.18146311603911383, -0.17401078727434263, -0.15419361061257214, -0.16225060576782063, -0.19010941396999181, -0.20901161407004412, -0.21670683975774699, -0.21771386270738891, -0.21481284465869752, -0.21063120869004387, -0.20669243364582401, -0.20320707386714859, -0.19984087691926442, -0.19725417448019505, -0.19633783049069981, -0.19650494599999785, -0.19708316838336942, -0.19779309449413818, -0.19853070294429562, -0.19917342167307153, -0.19964814677795845, -0.19991627610824922, -0.20013162970144974, -0.20029864969405509, -0.20036259676501131, -0.20030682824965193, -0.20016105135750167, -0.19997664501985918, -0.19980185871568762, -0.19966836175417696, -0.19958856744312226, -0.19955954696194517, -0.19956950051110917, -0.19960377086336181, -0.19964885299433241, -0.19969427478531132, -0.19973301547655564, -0.19976121574277764, -0.19977765285688653, -0.19978315117522441, -0.19977994634841739, -0.19977101394878494]
    1828        
    1829         G2 = [-0.40000000000000002, -0.39077401254732241, -0.33350466136630474, -0.29771023004255281, -0.27605439066140897, -0.25986156218997497, -0.24502185018573647, -0.231792624329521, -0.21981564668803993, -0.20870707082936543, -0.19877739883776599, -0.18980922837977957, -0.17308011674005838, -0.16306400164013773, -0.17798470933304333, -0.1929554075869116, -0.20236705191987037, -0.20695767560655007, -0.20841025876092567, -0.20792102174869989, -0.20655350005579293, -0.20492002526259828, -0.20310627026780645, -0.20105983335287836, -0.19937394565794653, -0.19853917506699659, -0.19836389977624452, -0.19850305023602796, -0.19877764028836831, -0.19910928131034669, -0.19943705712418805, -0.19970344172958865, -0.19991076989870474, -0.20010020127747646, -0.20025937787100062, -0.20035087292905965, -0.20035829921463297, -0.20029606557316171, -0.20019606915365515, -0.20009096093399206, -0.20000371608204368, -0.19994495432920584, -0.19991535665176338, -0.19990981826533513, -0.19992106419898723, -0.19994189853516578, -0.19996624091229293, -0.19998946016985167, -0.20000842303470234, -0.20002144460718174, -0.20002815561337187]
    1830        
    1831         G3 = [-0.45000000000000001, -0.37631169657400332, -0.33000044342859486, -0.30586045469008522, -0.28843572253009941, -0.27215308978603808, -0.25712951540331219, -0.2431608296216613, -0.23032023651386374, -0.2184546873456619, -0.20735123704254332, -0.19740397194806389, -0.1859829564064375, -0.16675980728362105, -0.16951575032846536, -0.1832860872609344, -0.19485758939241243, -0.20231368291811427, -0.20625610376074754, -0.20758116241495619, -0.20721445402086161, -0.20603406830353785, -0.20450262808396991, -0.2026769581185151, -0.2007401212066364, -0.19931160535777592, -0.19863606301128725, -0.19848511940572691, -0.19860091042948352, -0.19885490669377764, -0.19916542732701112, -0.19946678238611959, -0.19971209594104697, -0.19991912886512292, -0.2001058430788881, -0.20024959409472989, -0.20032160254609382, -0.20031583165752354, -0.20025051539293123, -0.2001556115816068, -0.20005952955420872, -0.1999814429561611, -0.19992977821558131, -0.19990457708664208, -0.19990104785490476, -0.19991257153954825, -0.19993258231880562, -0.19995548502882532, -0.19997700760919687, -0.19999429663503748, -0.20000588800248761]
    1832        
     1721        G0 = [-0.20000000000000001, -0.20000000000000001, -0.19920600846161715,
     1722              -0.19153647344085376, -0.19127622768281194, -0.1770671909675095,
     1723              -0.16739412133181927, -0.16196038919122191, -0.15621633053131384,
     1724              -0.15130021599977705, -0.13930978857215484, -0.19349274358263582,
     1725              -0.19975307598803765, -0.19999897143103357, -0.1999999995532111,
     1726              -0.19999999999949952, -0.19999999999949952, -0.19999999999949952,
     1727              -0.19997270012494556, -0.19925805948554556, -0.19934513778450533,
     1728              -0.19966484196394893, -0.1997352860102084,  -0.19968260481750394,
     1729              -0.19980280797303882, -0.19998804881822749, -0.19999999778075916,
     1730              -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
     1731              -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
     1732              -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
     1733              -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
     1734              -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
     1735              -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
     1736              -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
     1737              -0.19999999999966167, -0.19999999999966167, -0.19999999999966167]
     1738
     1739        G1 = [-0.29999999999999993, -0.29999588068034899, -0.29250047332330331,
     1740              -0.28335081844518584, -0.26142206997410805, -0.22656028856329835,
     1741              -0.21224087216745585, -0.19934324109114465, -0.1889857939783175,
     1742              -0.18146311603911383, -0.17401078727434263, -0.15419361061257214,
     1743              -0.16225060576782063, -0.19010941396999181, -0.20901161407004412,
     1744              -0.21670683975774699, -0.21771386270738891, -0.21481284465869752,
     1745              -0.21063120869004387, -0.20669243364582401, -0.20320707386714859,
     1746              -0.19984087691926442, -0.19725417448019505, -0.19633783049069981,
     1747              -0.19650494599999785, -0.19708316838336942, -0.19779309449413818,
     1748              -0.19853070294429562, -0.19917342167307153, -0.19964814677795845,
     1749              -0.19991627610824922, -0.20013162970144974, -0.20029864969405509,
     1750              -0.20036259676501131, -0.20030682824965193, -0.20016105135750167,
     1751              -0.19997664501985918, -0.19980185871568762, -0.19966836175417696,
     1752              -0.19958856744312226, -0.19955954696194517, -0.19956950051110917,
     1753              -0.19960377086336181, -0.19964885299433241, -0.19969427478531132,
     1754              -0.19973301547655564, -0.19976121574277764, -0.19977765285688653,
     1755              -0.19978315117522441, -0.19977994634841739, -0.19977101394878494]
     1756
     1757        G2 = [-0.40000000000000002, -0.39077401254732241, -0.33350466136630474,
     1758              -0.29771023004255281, -0.27605439066140897, -0.25986156218997497,
     1759              -0.24502185018573647, -0.231792624329521,   -0.21981564668803993,
     1760              -0.20870707082936543, -0.19877739883776599, -0.18980922837977957,
     1761              -0.17308011674005838, -0.16306400164013773, -0.17798470933304333,
     1762              -0.1929554075869116,  -0.20236705191987037, -0.20695767560655007,
     1763              -0.20841025876092567, -0.20792102174869989, -0.20655350005579293,
     1764              -0.20492002526259828, -0.20310627026780645, -0.20105983335287836,
     1765              -0.19937394565794653, -0.19853917506699659, -0.19836389977624452,
     1766              -0.19850305023602796, -0.19877764028836831, -0.19910928131034669,
     1767              -0.19943705712418805, -0.19970344172958865, -0.19991076989870474,
     1768              -0.20010020127747646, -0.20025937787100062, -0.20035087292905965,
     1769              -0.20035829921463297, -0.20029606557316171, -0.20019606915365515,
     1770              -0.20009096093399206, -0.20000371608204368, -0.19994495432920584,
     1771              -0.19991535665176338, -0.19990981826533513, -0.19992106419898723,
     1772              -0.19994189853516578, -0.19996624091229293, -0.19998946016985167,
     1773              -0.20000842303470234, -0.20002144460718174, -0.20002815561337187]
     1774
     1775        G3 = [-0.45000000000000001, -0.37631169657400332, -0.33000044342859486,
     1776              -0.30586045469008522, -0.28843572253009941, -0.27215308978603808,
     1777              -0.25712951540331219, -0.2431608296216613,  -0.23032023651386374,
     1778              -0.2184546873456619,  -0.20735123704254332, -0.19740397194806389,
     1779              -0.1859829564064375,  -0.16675980728362105, -0.16951575032846536,
     1780              -0.1832860872609344,  -0.19485758939241243, -0.20231368291811427,
     1781              -0.20625610376074754, -0.20758116241495619, -0.20721445402086161,
     1782              -0.20603406830353785, -0.20450262808396991, -0.2026769581185151,
     1783              -0.2007401212066364,  -0.19931160535777592, -0.19863606301128725,
     1784              -0.19848511940572691, -0.19860091042948352, -0.19885490669377764,
     1785              -0.19916542732701112, -0.19946678238611959, -0.19971209594104697,
     1786              -0.19991912886512292, -0.2001058430788881,  -0.20024959409472989,
     1787              -0.20032160254609382, -0.20031583165752354, -0.20025051539293123,
     1788              -0.2001556115816068,  -0.20005952955420872, -0.1999814429561611,
     1789              -0.19992977821558131, -0.19990457708664208, -0.19990104785490476,
     1790              -0.19991257153954825, -0.19993258231880562, -0.19995548502882532,
     1791              -0.19997700760919687, -0.19999429663503748, -0.20000588800248761]
     1792
    18331793        #FIXME (DSG):This is a hack so the anuga install, not precompiled
    18341794        # works on DSG's win2000, python 2.3
     
    18411801        #if len(gauge_values[3]) == 52: gauge_values[3].pop()
    18421802
    1843 ##         print len(G0), len(gauge_values[0])
    1844 ##         print len(G1), len(gauge_values[1])
    1845        
    1846         #print gauge_values[3]
    1847         #print G0[:4]
    1848         #print array(gauge_values[0])-array(G0)
    1849        
    1850        
    18511803        assert num.allclose(gauge_values[0], G0)
    18521804        assert num.allclose(gauge_values[1], G1)
    18531805        assert num.allclose(gauge_values[2], G2)
    1854         assert num.allclose(gauge_values[3], G3)       
    1855 
    1856 
    1857 
    1858 
    1859 
    1860 
     1806        assert num.allclose(gauge_values[3], G3)
    18611807
    18621808    #####################################################
     
    18641810    def test_flux_optimisation(self):
    18651811        """test_flux_optimisation
     1812
    18661813        Test that fluxes are correctly computed using
    18671814        dry and still cell exclusions
     
    18791826
    18801827        points = [a, b, c, d, e, f]
    1881         #bac, bce, ecf, dbe
    1882         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     1828        #             bac,     bce,     ecf,    dbe
     1829        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    18831830
    18841831        domain = Domain(points, vertices)
     
    18891836
    18901837        h = 0.1
    1891         def stage(x,y):
    1892             return slope(x,y)+h
     1838        def stage(x, y):
     1839            return slope(x, y) + h
    18931840
    18941841        domain.set_quantity('elevation', slope)
    18951842        domain.set_quantity('stage', stage)
    18961843
    1897         # Allow slope limiters to work (FIXME (Ole): Shouldn't this be automatic in ANUGA?)     
    1898         domain.distribute_to_vertices_and_edges()       
     1844        # Allow slope limiters to work
     1845        # (FIXME (Ole): Shouldn't this be automatic in ANUGA?)
     1846        domain.distribute_to_vertices_and_edges()
    18991847
    19001848        initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
     
    19021850        domain.set_boundary({'exterior': Reflective_boundary(domain)})
    19031851
    1904 
    1905         #  Check that update arrays are initialised to zero
     1852        #  Check that update arrays are initialised to zero
    19061853        assert num.allclose(domain.get_quantity('stage').explicit_update, 0)
    19071854        assert num.allclose(domain.get_quantity('xmomentum').explicit_update, 0)
    1908         assert num.allclose(domain.get_quantity('ymomentum').explicit_update, 0)               
    1909 
     1855        assert num.allclose(domain.get_quantity('ymomentum').explicit_update, 0)
    19101856
    19111857        # Get true values
     
    19141860        stage_ref = copy.copy(domain.get_quantity('stage').explicit_update)
    19151861        xmom_ref = copy.copy(domain.get_quantity('xmomentum').explicit_update)
    1916         ymom_ref = copy.copy(domain.get_quantity('ymomentum').explicit_update)       
     1862        ymom_ref = copy.copy(domain.get_quantity('ymomentum').explicit_update)
    19171863
    19181864        # Try with flux optimisation
     
    19201866        domain.compute_fluxes()
    19211867
    1922         assert num.allclose(stage_ref, domain.get_quantity('stage').explicit_update)
    1923         assert num.allclose(xmom_ref, domain.get_quantity('xmomentum').explicit_update)
    1924         assert num.allclose(ymom_ref, domain.get_quantity('ymomentum').explicit_update)
    1925        
    1926    
    1927        
     1868        assert num.allclose(stage_ref,
     1869                            domain.get_quantity('stage').explicit_update)
     1870        assert num.allclose(xmom_ref,
     1871                            domain.get_quantity('xmomentum').explicit_update)
     1872        assert num.allclose(ymom_ref,
     1873                            domain.get_quantity('ymomentum').explicit_update)
     1874
    19281875    def test_initial_condition(self):
    19291876        """test_initial_condition
     1877
    19301878        Test that initial condition is output at time == 0 and that
    19311879        computed values change as system evolves
     
    19431891
    19441892        points = [a, b, c, d, e, f]
    1945         #bac, bce, ecf, dbe
    1946         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     1893        #             bac,     bce,     ecf,    dbe
     1894        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    19471895
    19481896        domain = Domain(points, vertices)
     
    19531901
    19541902        h = 0.1
    1955         def stage(x,y):
    1956             return slope(x,y)+h
     1903        def stage(x, y):
     1904            return slope(x, y) + h
    19571905
    19581906        domain.set_quantity('elevation', slope)
    19591907        domain.set_quantity('stage', stage)
    19601908
    1961         # Allow slope limiters to work (FIXME (Ole): Shouldn't this be automatic in ANUGA?)     
    1962         domain.distribute_to_vertices_and_edges()       
     1909        # Allow slope limiters to work
     1910        # (FIXME (Ole): Shouldn't this be automatic in ANUGA?)
     1911        domain.distribute_to_vertices_and_edges()
    19631912
    19641913        initial_stage = copy.copy(domain.quantities['stage'].vertex_values)
     
    19671916
    19681917        domain.optimise_dry_cells = True
     1918
    19691919        #Evolution
    1970         for t in domain.evolve(yieldstep = 0.5, finaltime = 2.0):
     1920        for t in domain.evolve(yieldstep=0.5, finaltime=2.0):
    19711921            stage = domain.quantities['stage'].vertex_values
    19721922
     
    19761926                assert not num.allclose(stage, initial_stage)
    19771927
    1978 
    19791928        os.remove(domain.get_name() + '.sww')
    19801929
    1981 
    1982 
    19831930    #####################################################
     1931
    19841932    def test_gravity(self):
    19851933        #Assuming no friction
     
    19951943
    19961944        points = [a, b, c, d, e, f]
    1997         #bac, bce, ecf, dbe
    1998         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     1945        #             bac,     bce,     ecf,    dbe
     1946        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    19991947
    20001948        domain = Domain(points, vertices)
     
    20051953
    20061954        h = 0.1
    2007         def stage(x,y):
    2008             return slope(x,y)+h
     1955        def stage(x, y):
     1956            return slope(x, y) + h
    20091957
    20101958        domain.set_quantity('elevation', slope)
     
    20181966
    20191967        assert num.allclose(domain.quantities['stage'].explicit_update, 0)
    2020         assert num.allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3)
     1968        assert num.allclose(domain.quantities['xmomentum'].explicit_update,
     1969                            -g*h*3)
    20211970        assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
    2022 
    20231971
    20241972    def test_manning_friction(self):
     
    20331981
    20341982        points = [a, b, c, d, e, f]
    2035         #bac, bce, ecf, dbe
    2036         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     1983        #             bac,     bce,     ecf,    dbe
     1984        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    20371985
    20381986        domain = Domain(points, vertices)
     
    20431991
    20441992        h = 0.1
    2045         def stage(x,y):
    2046             return slope(x,y)+h
     1993        def stage(x, y):
     1994            return slope(x, y) + h
    20471995
    20481996        eta = 0.07
     
    20582006
    20592007        assert num.allclose(domain.quantities['stage'].explicit_update, 0)
    2060         assert num.allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3)
     2008        assert num.allclose(domain.quantities['xmomentum'].explicit_update,
     2009                            -g*h*3)
    20612010        assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
    20622011
    20632012        assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
    2064         assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, 0)
    2065         assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)
     2013        assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update,
     2014                            0)
     2015        assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update,
     2016                            0)
    20662017
    20672018        #Create some momentum for friction to work with
    20682019        domain.set_quantity('xmomentum', 1)
    2069         S = -g * eta**2 / h**(7.0/3)
     2020        S = -g*eta**2 / h**(7.0/3)
    20702021
    20712022        domain.compute_forcing_terms()
    20722023        assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
    2073         assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, S)
    2074         assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 0)
     2024        assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update,
     2025                            S)
     2026        assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update,
     2027                            0)
    20752028
    20762029        #A more complex example
     
    20822035        domain.set_quantity('ymomentum', 4)
    20832036
    2084         S = -g * eta**2 * 5 / h**(7.0/3)
    2085 
     2037        S = -g*eta**2*5 / h**(7.0/3)
    20862038
    20872039        domain.compute_forcing_terms()
    20882040
    20892041        assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
    2090         assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update, 3*S)
    2091         assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update, 4*S)
     2042        assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update,
     2043                            3*S)
     2044        assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update,
     2045                            4*S)
    20922046
    20932047    def test_constant_wind_stress(self):
     
    21032057
    21042058        points = [a, b, c, d, e, f]
    2105         #bac, bce, ecf, dbe
    2106         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    2107 
     2059        #             bac,     bce,     ecf,     dbe
     2060        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    21082061
    21092062        domain = Domain(points, vertices)
     
    21212074        phi = 135
    21222075        domain.forcing_terms = []
    2123         domain.forcing_terms.append( Wind_stress(s, phi) )
     2076        domain.forcing_terms.append(Wind_stress(s, phi))
    21242077
    21252078        domain.compute_forcing_terms()
    21262079
    2127 
    2128         const = eta_w*rho_a/rho_w
     2080        const = eta_w*rho_a / rho_w
    21292081
    21302082        #Convert to radians
    2131         phi = phi*pi/180
     2083        phi = phi*pi / 180
    21322084
    21332085        #Compute velocity vector (u, v)
     
    21412093        assert num.allclose(domain.quantities['xmomentum'].explicit_update, S*u)
    21422094        assert num.allclose(domain.quantities['ymomentum'].explicit_update, S*v)
    2143 
    21442095
    21452096    def test_variable_wind_stress(self):
     
    21552106
    21562107        points = [a, b, c, d, e, f]
    2157         #bac, bce, ecf, dbe
    2158         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     2108        #             bac,     bce,     ecf,    dbe
     2109        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    21592110
    21602111        domain = Domain(points, vertices)
     
    21682119        domain.set_boundary({'exterior': Br})
    21692120
    2170 
    2171         domain.time = 5.54 #Take a random time (not zero)
     2121        domain.time = 5.54    # Take a random time (not zero)
    21722122
    21732123        #Setup only one forcing term, constant wind stress
     
    21752125        phi = 135
    21762126        domain.forcing_terms = []
    2177         domain.forcing_terms.append( Wind_stress(s = speed, phi = angle) )
     2127        domain.forcing_terms.append(Wind_stress(s=speed, phi=angle))
    21782128
    21792129        domain.compute_forcing_terms()
    21802130
    21812131        #Compute reference solution
    2182         const = eta_w*rho_a/rho_w
    2183 
    2184         N = len(domain) # number_of_triangles
     2132        const = eta_w*rho_a / rho_w
     2133
     2134        N = len(domain)    # number_of_triangles
    21852135
    21862136        xc = domain.get_centroid_coordinates()
     
    21922142        phi_vec = angle(t,x,y)
    21932143
    2194 
    21952144        for k in range(N):
    2196             #Convert to radians
    2197             phi = phi_vec[k]*pi/180
     2145            # Convert to radians
     2146            phi = phi_vec[k]*pi / 180
    21982147            s = s_vec[k]
    21992148
    2200             #Compute velocity vector (u, v)
     2149            # Compute velocity vector (u, v)
    22012150            u = s*cos(phi)
    22022151            v = s*sin(phi)
    22032152
    2204             #Compute wind stress
     2153            # Compute wind stress
    22052154            S = const * num.sqrt(u**2 + v**2)
    22062155
    2207             assert num.allclose(domain.quantities['stage'].explicit_update[k], 0)
    2208             assert num.allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
    2209             assert num.allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
    2210 
    2211 
    2212 
    2213 
    2214 
     2156            assert num.allclose(domain.quantities['stage'].explicit_update[k],
     2157                                0)
     2158            assert num.allclose(domain.quantities['xmomentum'].\
     2159                                     explicit_update[k],
     2160                                S*u)
     2161            assert num.allclose(domain.quantities['ymomentum'].\
     2162                                     explicit_update[k],
     2163                                S*v)
    22152164
    22162165    def test_windfield_from_file(self):
     2166        import time
    22172167        from anuga.config import rho_a, rho_w, eta_w
    22182168        from math import pi, cos, sin
    22192169        from anuga.config import time_format
    22202170        from anuga.abstract_2d_finite_volumes.util import file_function
    2221         import time
    2222 
    22232171
    22242172        a = [0.0, 0.0]
     
    22302178
    22312179        points = [a, b, c, d, e, f]
    2232         #bac, bce, ecf, dbe
    2233         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     2180        #             bac,     bce,     ecf,    dbe
     2181        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    22342182
    22352183        domain = Domain(points, vertices)
    22362184
    2237         #Flat surface with 1m of water
     2185        # Flat surface with 1m of water
    22382186        domain.set_quantity('elevation', 0)
    22392187        domain.set_quantity('stage', 1.0)
     
    22432191        domain.set_boundary({'exterior': Br})
    22442192
    2245 
    2246         domain.time = 7 #Take a time that is represented in file (not zero)
    2247 
    2248         #Write wind stress file (ensure that domain.time is covered)
    2249         #Take x=1 and y=0
     2193        domain.time = 7    # Take a time that is represented in file (not zero)
     2194
     2195        # Write wind stress file (ensure that domain.time is covered)
     2196        # Take x=1 and y=0
    22502197        filename = 'test_windstress_from_file'
    22512198        start = time.mktime(time.strptime('2000', '%Y'))
    22522199        fid = open(filename + '.txt', 'w')
    2253         dt = 1  #One second interval
     2200        dt = 1    # One second interval
    22542201        t = 0.0
    22552202        while t <= 10.0:
    22562203            t_string = time.strftime(time_format, time.gmtime(t+start))
    22572204
    2258             fid.write('%s, %f %f\n' %(t_string,
    2259                                       speed(t,[1],[0])[0],
    2260                                       angle(t,[1],[0])[0]))
     2205            fid.write('%s, %f %f\n' %
     2206                      (t_string, speed(t,[1],[0])[0], angle(t,[1],[0])[0]))
    22612207            t += dt
    22622208
    22632209        fid.close()
    22642210
    2265 
    2266         #Convert ASCII file to NetCDF (Which is what we really like!)
    2267         from data_manager import timefile2netcdf       
     2211        # Convert ASCII file to NetCDF (Which is what we really like!)
     2212        from data_manager import timefile2netcdf
     2213
    22682214        timefile2netcdf(filename)
    22692215        os.remove(filename + '.txt')
    22702216
    2271        
    2272         #Setup wind stress
    2273         F = file_function(filename + '.tms', quantities = ['Attribute0',
    2274                                                            'Attribute1'])
     2217        # Setup wind stress
     2218        F = file_function(filename + '.tms',
     2219                          quantities=['Attribute0', 'Attribute1'])
    22752220        os.remove(filename + '.tms')
    2276        
    2277 
    2278         #print 'F(5)', F(5)
    2279 
    2280         #print 'F(5,x,y)', F(5,x=zeros(3),y=zeros(3))
    2281        
    2282         #print dir(F)
    2283         #print F.T
    2284         #print F.precomputed_values
    2285         #
    2286         #F = file_function(filename + '.txt')       
    2287         #
    2288         #print dir(F)
    2289         #print F.T       
    2290         #print F.Q
    2291        
     2221
    22922222        W = Wind_stress(F)
    22932223
     
    22972227        domain.compute_forcing_terms()
    22982228
    2299         #Compute reference solution
    2300         const = eta_w*rho_a/rho_w
    2301 
    2302         N = len(domain) # number_of_triangles
     2229        # Compute reference solution
     2230        const = eta_w*rho_a / rho_w
     2231
     2232        N = len(domain)    # number_of_triangles
    23032233
    23042234        t = domain.time
    23052235
    2306         s = speed(t,[1],[0])[0]
    2307         phi = angle(t,[1],[0])[0]
    2308 
    2309         #Convert to radians
    2310         phi = phi*pi/180
    2311 
    2312 
    2313         #Compute velocity vector (u, v)
     2236        s = speed(t, [1], [0])[0]
     2237        phi = angle(t, [1], [0])[0]
     2238
     2239        # Convert to radians
     2240        phi = phi*pi / 180
     2241
     2242        # Compute velocity vector (u, v)
    23142243        u = s*cos(phi)
    23152244        v = s*sin(phi)
    23162245
    2317         #Compute wind stress
     2246        # Compute wind stress
    23182247        S = const * num.sqrt(u**2 + v**2)
    23192248
    23202249        for k in range(N):
    2321             assert num.allclose(domain.quantities['stage'].explicit_update[k], 0)
    2322             assert num.allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
    2323             assert num.allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
    2324 
     2250            assert num.allclose(domain.quantities['stage'].explicit_update[k],
     2251                                0)
     2252            assert num.allclose(domain.quantities['xmomentum'].\
     2253                                    explicit_update[k],
     2254                                S*u)
     2255            assert num.allclose(domain.quantities['ymomentum'].\
     2256                                    explicit_update[k],
     2257                                S*v)
    23252258
    23262259    def test_windfield_from_file_seconds(self):
     2260        import time
    23272261        from anuga.config import rho_a, rho_w, eta_w
    23282262        from math import pi, cos, sin
    23292263        from anuga.config import time_format
    23302264        from anuga.abstract_2d_finite_volumes.util import file_function
    2331         import time
    2332 
    23332265
    23342266        a = [0.0, 0.0]
     
    23402272
    23412273        points = [a, b, c, d, e, f]
    2342         #bac, bce, ecf, dbe
    2343         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     2274        #             bac,     bce,     ecf,    dbe
     2275        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    23442276
    23452277        domain = Domain(points, vertices)
    23462278
    2347         #Flat surface with 1m of water
     2279        # Flat surface with 1m of water
    23482280        domain.set_quantity('elevation', 0)
    23492281        domain.set_quantity('stage', 1.0)
     
    23532285        domain.set_boundary({'exterior': Br})
    23542286
    2355 
    2356         domain.time = 7 #Take a time that is represented in file (not zero)
    2357 
    2358         #Write wind stress file (ensure that domain.time is covered)
    2359         #Take x=1 and y=0
     2287        domain.time = 7    # Take a time that is represented in file (not zero)
     2288
     2289        # Write wind stress file (ensure that domain.time is covered)
     2290        # Take x=1 and y=0
    23602291        filename = 'test_windstress_from_file'
    23612292        start = time.mktime(time.strptime('2000', '%Y'))
    23622293        fid = open(filename + '.txt', 'w')
    2363         dt = 0.5 #1  #One second interval
     2294        dt = 0.5    # Half second interval
    23642295        t = 0.0
    23652296        while t <= 10.0:
    2366             fid.write('%s, %f %f\n' %(str(t),
    2367                                       speed(t,[1],[0])[0],
    2368                                       angle(t,[1],[0])[0]))
     2297            fid.write('%s, %f %f\n'
     2298                      % (str(t), speed(t, [1], [0])[0], angle(t, [1], [0])[0]))
    23692299            t += dt
    23702300
    23712301        fid.close()
    23722302
    2373 
    2374         #Convert ASCII file to NetCDF (Which is what we really like!)
    2375         from data_manager import timefile2netcdf       
     2303        # Convert ASCII file to NetCDF (Which is what we really like!)
     2304        from data_manager import timefile2netcdf
     2305
    23762306        timefile2netcdf(filename, time_as_seconds=True)
    23772307        os.remove(filename + '.txt')
    23782308
    2379        
    2380         #Setup wind stress
    2381         F = file_function(filename + '.tms', quantities = ['Attribute0',
    2382                                                            'Attribute1'])
     2309        # Setup wind stress
     2310        F = file_function(filename + '.tms',
     2311                          quantities=['Attribute0', 'Attribute1'])
    23832312        os.remove(filename + '.tms')
    2384        
    2385 
    2386         #print 'F(5)', F(5)
    2387 
    2388         #print 'F(5,x,y)', F(5,x=zeros(3),y=zeros(3))
    2389        
    2390         #print dir(F)
    2391         #print F.T
    2392         #print F.precomputed_values
    2393         #
    2394         #F = file_function(filename + '.txt')       
    2395         #
    2396         #print dir(F)
    2397         #print F.T       
    2398         #print F.Q
    2399        
     2313
    24002314        W = Wind_stress(F)
    24012315
     
    24052319        domain.compute_forcing_terms()
    24062320
    2407         #Compute reference solution
    2408         const = eta_w*rho_a/rho_w
    2409 
    2410         N = len(domain) # number_of_triangles
     2321        # Compute reference solution
     2322        const = eta_w*rho_a / rho_w
     2323
     2324        N = len(domain)    # number_of_triangles
    24112325
    24122326        t = domain.time
    24132327
    2414         s = speed(t,[1],[0])[0]
    2415         phi = angle(t,[1],[0])[0]
    2416 
    2417         #Convert to radians
    2418         phi = phi*pi/180
    2419 
    2420 
    2421         #Compute velocity vector (u, v)
     2328        s = speed(t, [1], [0])[0]
     2329        phi = angle(t, [1], [0])[0]
     2330
     2331        # Convert to radians
     2332        phi = phi*pi / 180
     2333
     2334        # Compute velocity vector (u, v)
    24222335        u = s*cos(phi)
    24232336        v = s*sin(phi)
    24242337
    2425         #Compute wind stress
     2338        # Compute wind stress
    24262339        S = const * num.sqrt(u**2 + v**2)
    24272340
    24282341        for k in range(N):
    2429             assert num.allclose(domain.quantities['stage'].explicit_update[k], 0)
    2430             assert num.allclose(domain.quantities['xmomentum'].explicit_update[k], S*u)
    2431             assert num.allclose(domain.quantities['ymomentum'].explicit_update[k], S*v)
    2432 
    2433 
    2434        
     2342            assert num.allclose(domain.quantities['stage'].explicit_update[k],
     2343                                0)
     2344            assert num.allclose(domain.quantities['xmomentum'].\
     2345                                    explicit_update[k],
     2346                                S*u)
     2347            assert num.allclose(domain.quantities['ymomentum'].\
     2348                                    explicit_update[k],
     2349                                S*v)
    24352350
    24362351    def test_wind_stress_error_condition(self):
     
    24392354        """
    24402355
     2356        from math import pi, cos, sin
    24412357        from anuga.config import rho_a, rho_w, eta_w
    2442         from math import pi, cos, sin
    24432358
    24442359        a = [0.0, 0.0]
     
    24502365
    24512366        points = [a, b, c, d, e, f]
    2452         #bac, bce, ecf, dbe
    2453         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     2367        #             bac,     bce,     ecf,    dbe
     2368        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    24542369
    24552370        domain = Domain(points, vertices)
    24562371
    2457         #Flat surface with 1m of water
     2372        # Flat surface with 1m of water
    24582373        domain.set_quantity('elevation', 0)
    24592374        domain.set_quantity('stage', 1.0)
     
    24632378        domain.set_boundary({'exterior': Br})
    24642379
    2465 
    2466         domain.time = 5.54 #Take a random time (not zero)
    2467 
    2468         #Setup only one forcing term, bad func
     2380        domain.time = 5.54    # Take a random time (not zero)
     2381
     2382        # Setup only one forcing term, bad func
    24692383        domain.forcing_terms = []
    24702384
    24712385        try:
    2472             domain.forcing_terms.append(Wind_stress(s = scalar_func_list,
    2473                                                     phi = angle))
     2386            domain.forcing_terms.append(Wind_stress(s=scalar_func_list,
     2387                                                    phi=angle))
    24742388        except AssertionError:
    24752389            pass
    24762390        else:
    24772391            msg = 'Should have raised exception'
    2478             raise msg
    2479 
     2392            raise Exception, msg
    24802393
    24812394        try:
    2482             domain.forcing_terms.append(Wind_stress(s = speed,
    2483                                                     phi = scalar_func))
     2395            domain.forcing_terms.append(Wind_stress(s=speed, phi=scalar_func))
    24842396        except Exception:
    24852397            pass
    24862398        else:
    24872399            msg = 'Should have raised exception'
    2488             raise msg
     2400            raise Exception, msg
    24892401
    24902402        try:
    2491             domain.forcing_terms.append(Wind_stress(s = speed,
    2492                                                     phi = 'xx'))
     2403            domain.forcing_terms.append(Wind_stress(s=speed, phi='xx'))
    24932404        except:
    24942405            pass
    24952406        else:
    24962407            msg = 'Should have raised exception'
    2497             raise msg
    2498 
    2499 
     2408            raise Exception, msg
    25002409
    25012410    def test_rainfall(self):
     
    25102419
    25112420        points = [a, b, c, d, e, f]
    2512         #bac, bce, ecf, dbe
    2513         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    2514 
     2421        #             bac,     bce,     ecf,     dbe
     2422        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    25152423
    25162424        domain = Domain(points, vertices)
    25172425
    2518         #Flat surface with 1m of water
     2426        # Flat surface with 1m of water
    25192427        domain.set_quantity('elevation', 0)
    25202428        domain.set_quantity('stage', 1.0)
     
    25262434        # Setup only one forcing term, constant rainfall
    25272435        domain.forcing_terms = []
    2528         domain.forcing_terms.append( Rainfall(domain, rate=2.0) )
     2436        domain.forcing_terms.append(Rainfall(domain, rate=2.0))
    25292437
    25302438        domain.compute_forcing_terms()
    2531         assert num.allclose(domain.quantities['stage'].explicit_update, 2.0/1000)
    2532 
    2533 
     2439        assert num.allclose(domain.quantities['stage'].explicit_update,
     2440                            2.0/1000)
    25342441
    25352442    def test_rainfall_restricted_by_polygon(self):
     
    25442451
    25452452        points = [a, b, c, d, e, f]
    2546         #bac, bce, ecf, dbe
    2547         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    2548 
     2453        #             bac,     bce,     ecf,     dbe
     2454        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    25492455
    25502456        domain = Domain(points, vertices)
    25512457
    2552         #Flat surface with 1m of water
     2458        # Flat surface with 1m of water
    25532459        domain.set_quantity('elevation', 0)
    25542460        domain.set_quantity('stage', 1.0)
     
    25582464        domain.set_boundary({'exterior': Br})
    25592465
    2560         # Setup only one forcing term, constant rainfall restricted to a polygon enclosing triangle #1 (bce)
     2466        # Setup only one forcing term, constant rainfall
     2467        # restricted to a polygon enclosing triangle #1 (bce)
    25612468        domain.forcing_terms = []
    2562         R = Rainfall(domain,
    2563                      rate=2.0,
    2564                      polygon = [[1,1], [2,1], [2,2], [1,2]])
     2469        R = Rainfall(domain, rate=2.0, polygon=[[1,1], [2,1], [2,2], [1,2]])
    25652470
    25662471        assert num.allclose(R.exchange_area, 1)
    2567        
     2472
    25682473        domain.forcing_terms.append(R)
    25692474
    25702475        domain.compute_forcing_terms()
    2571         #print domain.quantities['stage'].explicit_update
    2572        
     2476
    25732477        assert num.allclose(domain.quantities['stage'].explicit_update[1],
    25742478                            2.0/1000)
    25752479        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
    2576         assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)       
    2577        
    2578 
     2480        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
    25792481
    25802482    def test_time_dependent_rainfall_restricted_by_polygon(self):
    2581 
    25822483        a = [0.0, 0.0]
    25832484        b = [0.0, 2.0]
     
    25882489
    25892490        points = [a, b, c, d, e, f]
    2590         #bac, bce, ecf, dbe
    2591         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    2592 
     2491        #             bac,     bce,     ecf,     dbe
     2492        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    25932493
    25942494        domain = Domain(points, vertices)
    25952495
    2596         #Flat surface with 1m of water
     2496        # Flat surface with 1m of water
    25972497        domain.set_quantity('elevation', 0)
    25982498        domain.set_quantity('stage', 1.0)
     
    26022502        domain.set_boundary({'exterior': Br})
    26032503
    2604         # Setup only one forcing term, time dependent rainfall restricted to a polygon enclosing triangle #1 (bce)
     2504        # Setup only one forcing term, time dependent rainfall
     2505        # restricted to a polygon enclosing triangle #1 (bce)
    26052506        domain.forcing_terms = []
    2606         R = Rainfall(domain,
    2607                      rate=lambda t: 3*t + 7,
    2608                      polygon = [[1,1], [2,1], [2,2], [1,2]])
     2507        R = Rainfall(domain, rate=lambda t: 3*t + 7,
     2508                     polygon=[[1,1], [2,1], [2,2], [1,2]])
    26092509
    26102510        assert num.allclose(R.exchange_area, 1)
    2611        
     2511
    26122512        domain.forcing_terms.append(R)
    26132513
    2614 
    26152514        domain.time = 10.
    26162515
    26172516        domain.compute_forcing_terms()
    2618         #print domain.quantities['stage'].explicit_update
    2619        
     2517
    26202518        assert num.allclose(domain.quantities['stage'].explicit_update[1],
    2621                             (3*domain.time+7)/1000)
     2519                            (3*domain.time + 7)/1000)
    26222520        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
    2623         assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)       
    2624        
    2625 
    2626 
     2521        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
    26272522
    26282523    def test_time_dependent_rainfall_using_starttime(self):
    2629    
    2630         rainfall_poly = ensure_numeric([[1,1], [2,1], [2,2], [1,2]], num.float)   
     2524        rainfall_poly = ensure_numeric([[1,1], [2,1], [2,2], [1,2]], num.float)
    26312525
    26322526        a = [0.0, 0.0]
     
    26382532
    26392533        points = [a, b, c, d, e, f]
    2640         #bac, bce, ecf, dbe
    2641         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    2642 
     2534        #             bac,     bce,     ecf,     dbe
     2535        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    26432536
    26442537        domain = Domain(points, vertices)
    26452538
    2646         #Flat surface with 1m of water
     2539        # Flat surface with 1m of water
    26472540        domain.set_quantity('elevation', 0)
    26482541        domain.set_quantity('stage', 1.0)
     
    26522545        domain.set_boundary({'exterior': Br})
    26532546
    2654         # Setup only one forcing term, time dependent rainfall restricted to a polygon enclosing triangle #1 (bce)
     2547        # Setup only one forcing term, time dependent rainfall
     2548        # restricted to a polygon enclosing triangle #1 (bce)
    26552549        domain.forcing_terms = []
    2656         R = Rainfall(domain,
    2657                      rate=lambda t: 3*t + 7,
    2658                      polygon=rainfall_poly)                     
     2550        R = Rainfall(domain, rate=lambda t: 3*t + 7,
     2551                     polygon=rainfall_poly)
    26592552
    26602553        assert num.allclose(R.exchange_area, 1)
    2661        
     2554
    26622555        domain.forcing_terms.append(R)
    26632556
     
    26692562
    26702563        domain.compute_forcing_terms()
    2671         #print domain.quantities['stage'].explicit_update
    2672 
    2673         #print domain.get_time()
     2564
    26742565        assert num.allclose(domain.quantities['stage'].explicit_update[1],
    2675                             (3*domain.get_time()+7)/1000)
     2566                            (3*domain.get_time() + 7)/1000)
    26762567        assert num.allclose(domain.quantities['stage'].explicit_update[1],
    2677                             (3*(domain.time + domain.starttime)+7)/1000)
     2568                            (3*(domain.time + domain.starttime) + 7)/1000)
    26782569
    26792570        # Using internal time her should fail
    26802571        assert not num.allclose(domain.quantities['stage'].explicit_update[1],
    2681                                 (3*domain.time+7)/1000)               
     2572                                (3*domain.time + 7)/1000)
    26822573
    26832574        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
    2684         assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)       
    2685        
    2686 
    2687 
    2688        
     2575        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
     2576
    26892577    def test_time_dependent_rainfall_using_georef(self):
    26902578        """test_time_dependent_rainfall_using_georef
    2691        
     2579
    26922580        This will also test the General forcing term using georef
    26932581        """
    2694        
    2695         #Mesh in zone 56 (absolute coords)
    2696 
     2582
     2583        # Mesh in zone 56 (absolute coords)
    26972584        x0 = 314036.58727982
    26982585        y0 = 6224951.2960092
    26992586
    2700        
    27012587        rainfall_poly = ensure_numeric([[1,1], [2,1], [2,2], [1,2]], num.float)
    27022588        rainfall_poly += [x0, y0]
     
    27102596
    27112597        points = [a, b, c, d, e, f]
    2712         #bac, bce, ecf, dbe
    2713         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    2714 
     2598        #             bac,     bce,     ecf,     dbe
     2599        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    27152600
    27162601        domain = Domain(points, vertices,
    2717                         geo_reference = Geo_reference(56, x0, y0))
    2718 
    2719         #Flat surface with 1m of water
     2602                        geo_reference=Geo_reference(56, x0, y0))
     2603
     2604        # Flat surface with 1m of water
    27202605        domain.set_quantity('elevation', 0)
    27212606        domain.set_quantity('stage', 1.0)
     
    27252610        domain.set_boundary({'exterior': Br})
    27262611
    2727         # Setup only one forcing term, time dependent rainfall restricted to a polygon enclosing triangle #1 (bce)
     2612        # Setup only one forcing term, time dependent rainfall
     2613        # restricted to a polygon enclosing triangle #1 (bce)
    27282614        domain.forcing_terms = []
    2729         R = Rainfall(domain,
    2730                      rate=lambda t: 3*t + 7,
    2731                      polygon=rainfall_poly)
     2615        R = Rainfall(domain, rate=lambda t: 3*t + 7, polygon=rainfall_poly)
    27322616
    27332617        assert num.allclose(R.exchange_area, 1)
    2734        
     2618
    27352619        domain.forcing_terms.append(R)
    27362620
     
    27422626
    27432627        domain.compute_forcing_terms()
    2744         #print domain.quantities['stage'].explicit_update
    2745 
    2746         #print domain.get_time()
     2628
    27472629        assert num.allclose(domain.quantities['stage'].explicit_update[1],
    2748                             (3*domain.get_time()+7)/1000)
     2630                            (3*domain.get_time() + 7)/1000)
    27492631        assert num.allclose(domain.quantities['stage'].explicit_update[1],
    2750                             (3*(domain.time + domain.starttime)+7)/1000)
     2632                            (3*(domain.time + domain.starttime) + 7)/1000)
    27512633
    27522634        # Using internal time her should fail
    27532635        assert not num.allclose(domain.quantities['stage'].explicit_update[1],
    2754                                 (3*domain.time+7)/1000)               
     2636                                (3*domain.time + 7)/1000)
    27552637
    27562638        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
    2757         assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)       
    2758        
    2759 
    2760        
    2761        
    2762 
     2639        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
    27632640
    27642641    def test_time_dependent_rainfall_restricted_by_polygon_with_default(self):
    2765         """test_time_dependent_rainfall_restricted_by_polygon_with_default
    2766 
     2642        """
    27672643        Test that default rainfall can be used when given rate runs out of data.
    27682644        """
     2645
    27692646        a = [0.0, 0.0]
    27702647        b = [0.0, 2.0]
     
    27752652
    27762653        points = [a, b, c, d, e, f]
    2777         #bac, bce, ecf, dbe
    2778         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    2779 
     2654        #             bac,     bce,     ecf,     dbe
     2655        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    27802656
    27812657        domain = Domain(points, vertices)
    27822658
    2783         #Flat surface with 1m of water
     2659        # Flat surface with 1m of water
    27842660        domain.set_quantity('elevation', 0)
    27852661        domain.set_quantity('stage', 1.0)
     
    27892665        domain.set_boundary({'exterior': Br})
    27902666
    2791         # Setup only one forcing term, time dependent rainfall that expires at t==20
     2667        # Setup only one forcing term, time dependent rainfall
     2668        # that expires at t==20
    27922669        from anuga.fit_interpolate.interpolate import Modeltime_too_late
     2670
    27932671        def main_rate(t):
    27942672            if t > 20:
     
    27972675            else:
    27982676                return 3*t + 7
    2799        
     2677
    28002678        domain.forcing_terms = []
    2801         R = Rainfall(domain,
    2802                      rate=main_rate,
    2803                      polygon = [[1,1], [2,1], [2,2], [1,2]],
    2804                      default_rate=5.0)
     2679        R = Rainfall(domain, rate=main_rate,
     2680                     polygon=[[1,1], [2,1], [2,2], [1,2]], default_rate=5.0)
    28052681
    28062682        assert num.allclose(R.exchange_area, 1)
    2807        
     2683
    28082684        domain.forcing_terms.append(R)
    28092685
    2810 
    28112686        domain.time = 10.
    28122687
    28132688        domain.compute_forcing_terms()
    2814         #print domain.quantities['stage'].explicit_update
    2815        
    2816         assert num.allclose(domain.quantities['stage'].explicit_update[1], (3*domain.time+7)/1000)
     2689
     2690        assert num.allclose(domain.quantities['stage'].explicit_update[1],
     2691                            (3*domain.time+7)/1000)
    28172692        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
    2818         assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)       
    2819 
     2693        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
    28202694
    28212695        domain.time = 100.
    2822         domain.quantities['stage'].explicit_update[:] = 0.0  # Reset
     2696        domain.quantities['stage'].explicit_update[:] = 0.0     # Reset
    28232697        domain.compute_forcing_terms()
    2824         #print domain.quantities['stage'].explicit_update
    2825        
    2826         assert num.allclose(domain.quantities['stage'].explicit_update[1], 5.0/1000) # Default value
     2698
     2699        assert num.allclose(domain.quantities['stage'].explicit_update[1],
     2700                            5.0/1000) # Default value
    28272701        assert num.allclose(domain.quantities['stage'].explicit_update[0], 0)
    2828         assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)       
    2829        
    2830 
    2831 
    2832 
    2833        
    2834 
     2702        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
    28352703
    28362704    def test_rainfall_forcing_with_evolve(self):
     
    28392707        Test how forcing terms are called within evolve
    28402708        """
    2841        
     2709
    28422710        # FIXME(Ole): This test is just to experiment
    2843        
     2711
    28442712        a = [0.0, 0.0]
    28452713        b = [0.0, 2.0]
     
    28502718
    28512719        points = [a, b, c, d, e, f]
    2852         #bac, bce, ecf, dbe
    2853         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    2854 
     2720        #             bac,     bce,     ecf,     dbe
     2721        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    28552722
    28562723        domain = Domain(points, vertices)
    28572724
    2858         #Flat surface with 1m of water
     2725        # Flat surface with 1m of water
    28592726        domain.set_quantity('elevation', 0)
    28602727        domain.set_quantity('stage', 1.0)
     
    28642731        domain.set_boundary({'exterior': Br})
    28652732
    2866         # Setup only one forcing term, time dependent rainfall that expires at t==20
     2733        # Setup only one forcing term, time dependent rainfall
     2734        # that expires at t==20
    28672735        from anuga.fit_interpolate.interpolate import Modeltime_too_late
     2736
    28682737        def main_rate(t):
    28692738            if t > 20:
     
    28722741            else:
    28732742                return 3*t + 7
    2874        
     2743
    28752744        domain.forcing_terms = []
    2876         R = Rainfall(domain,
    2877                      rate=main_rate,
    2878                      polygon=[[1,1], [2,1], [2,2], [1,2]],
    2879                      default_rate=5.0)
     2745        R = Rainfall(domain, rate=main_rate,
     2746                     polygon=[[1,1], [2,1], [2,2], [1,2]], default_rate=5.0)
    28802747
    28812748        assert num.allclose(R.exchange_area, 1)
    2882        
     2749
    28832750        domain.forcing_terms.append(R)
    28842751
    28852752        for t in domain.evolve(yieldstep=1, finaltime=25):
    28862753            pass
    2887            
    2888             #print t, domain.quantities['stage'].explicit_update, (3*t+7)/1000
    2889            
    28902754            #FIXME(Ole):  A test here is hard because explicit_update also
    28912755            # receives updates from the flux calculation.
    2892 
    2893        
    2894        
    28952756
    28962757    def test_inflow_using_circle(self):
     
    29052766
    29062767        points = [a, b, c, d, e, f]
    2907         #bac, bce, ecf, dbe
    2908         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    2909 
     2768        #             bac,     bce,     ecf,     dbe
     2769        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    29102770
    29112771        domain = Domain(points, vertices)
     
    29192779        domain.set_boundary({'exterior': Br})
    29202780
    2921         # Setup only one forcing term, constant inflow of 2 m^3/s on a circle affecting triangles #0 and #1 (bac and bce)
     2781        # Setup only one forcing term, constant inflow of 2 m^3/s
     2782        # on a circle affecting triangles #0 and #1 (bac and bce)
    29222783        domain.forcing_terms = []
    2923         domain.forcing_terms.append( Inflow(domain, rate=2.0, center=(1,1), radius=1) )
     2784        domain.forcing_terms.append(Inflow(domain, rate=2.0,
     2785                                           center=(1,1), radius=1))
    29242786
    29252787        domain.compute_forcing_terms()
    2926         #print domain.quantities['stage'].explicit_update
    2927        
    2928         assert num.allclose(domain.quantities['stage'].explicit_update[1], 2.0/pi)
    2929         assert num.allclose(domain.quantities['stage'].explicit_update[0], 2.0/pi)
    2930         assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)       
    2931 
     2788
     2789        assert num.allclose(domain.quantities['stage'].explicit_update[1],
     2790                            2.0/pi)
     2791        assert num.allclose(domain.quantities['stage'].explicit_update[0],
     2792                            2.0/pi)
     2793        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
    29322794
    29332795    def test_inflow_using_circle_function(self):
     
    29422804
    29432805        points = [a, b, c, d, e, f]
    2944         #bac, bce, ecf, dbe
    2945         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    2946 
     2806        #             bac,     bce,     ecf,     dbe
     2807        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    29472808
    29482809        domain = Domain(points, vertices)
     
    29562817        domain.set_boundary({'exterior': Br})
    29572818
    2958         # Setup only one forcing term, time dependent inflow of 2 m^3/s on a circle affecting triangles #0 and #1 (bac and bce)
     2819        # Setup only one forcing term, time dependent inflow of 2 m^3/s
     2820        # on a circle affecting triangles #0 and #1 (bac and bce)
    29592821        domain.forcing_terms = []
    2960         domain.forcing_terms.append( Inflow(domain, rate=lambda t: 2., center=(1,1), radius=1) )
     2822        domain.forcing_terms.append(Inflow(domain, rate=lambda t: 2.,
     2823                                           center=(1,1), radius=1))
    29612824
    29622825        domain.compute_forcing_terms()
    2963        
    2964         assert num.allclose(domain.quantities['stage'].explicit_update[1], 2.0/pi)
    2965         assert num.allclose(domain.quantities['stage'].explicit_update[0], 2.0/pi)
    2966         assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)       
    2967        
    2968 
    2969 
     2826
     2827        assert num.allclose(domain.quantities['stage'].explicit_update[1],
     2828                            2.0/pi)
     2829        assert num.allclose(domain.quantities['stage'].explicit_update[0],
     2830                            2.0/pi)
     2831        assert num.allclose(domain.quantities['stage'].explicit_update[2:], 0)
    29702832
    29712833    def test_inflow_catch_too_few_triangles(self):
    2972         """test_inflow_catch_too_few_triangles
    2973        
    2974         Test that exception is thrown if no triangles are covered by the inflow area
    29752834        """
     2835        Test that exception is thrown if no triangles are covered
     2836        by the inflow area
     2837        """
     2838
    29762839        from math import pi, cos, sin
    29772840
     
    29842847
    29852848        points = [a, b, c, d, e, f]
    2986         #bac, bce, ecf, dbe
    2987         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    2988 
     2849        #             bac,     bce,     ecf,     dbe
     2850        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    29892851
    29902852        domain = Domain(points, vertices)
     
    29982860        domain.set_boundary({'exterior': Br})
    29992861
    3000         # Setup only one forcing term, constant inflow of 2 m^3/s on a circle affecting triangles #0 and #1 (bac and bce)
    3001 
     2862        # Setup only one forcing term, constant inflow of 2 m^3/s
     2863        # on a circle affecting triangles #0 and #1 (bac and bce)
    30022864        try:
    30032865            Inflow(domain, rate=2.0, center=(1,1.1), radius=0.01)
     
    30082870            raise Exception, msg
    30092871
    3010 
    3011 
    3012 
    30132872    def Xtest_inflow_outflow_conservation(self):
    3014         """test_inflow_outflow_conservation
    3015        
    3016         Test what happens if water is abstracted from one area and
     2873        """
     2874        Test what happens if water is abstracted from one area and
    30172875        injected into another - especially if there is not enough
    3018         water to match the abstraction. 
     2876        water to match the abstraction.
    30192877        This tests that the total volume is kept constant under a range of
    30202878        scenarios.
    3021        
     2879
    30222880        This test will fail as the problem was only fixed for culverts.
    30232881        """
    3024        
     2882
    30252883        from math import pi, cos, sin
    3026        
     2884
    30272885        length = 20.
    30282886        width = 10.
    30292887
    3030         dx = dy = 2  # 1 or 2 OK
     2888        dx = dy = 2    # 1 or 2 OK
    30312889        points, vertices, boundary = rectangular_cross(int(length/dx),
    30322890                                                       int(width/dy),
    3033                                                        len1=length, 
     2891                                                       len1=length,
    30342892                                                       len2=width)
    3035         domain = Domain(points, vertices, boundary)   
    3036         domain.set_name('test_inflow_conservation')  # Output name
     2893        domain = Domain(points, vertices, boundary)
     2894        domain.set_name('test_inflow_conservation')    # Output name
    30372895        domain.set_default_order(2)
    3038        
    30392896
    30402897        # Flat surface with 1m of water
     
    30472904        domain.set_boundary({'left': Br, 'right': Br, 'bottom': Br, 'top': Br})
    30482905
    3049         # Setup one forcing term, constant inflow of 2 m^3/s on a circle 
     2906        # Setup one forcing term, constant inflow of 2 m^3/s on a circle
    30502907        domain.forcing_terms = []
    3051         domain.forcing_terms.append(Inflow(domain, rate=2.0, center=(5,5), radius=1))
     2908        domain.forcing_terms.append(Inflow(domain, rate=2.0,
     2909                                           center=(5,5), radius=1))
    30522910
    30532911        domain.compute_forcing_terms()
    3054         #print domain.quantities['stage'].explicit_update
    3055        
     2912
    30562913        # Check that update values are correct
    30572914        for x in domain.quantities['stage'].explicit_update:
    30582915            assert num.allclose(x, 2.0/pi) or num.allclose(x, 0.0)
    30592916
    3060        
    30612917        # Check volumes without inflow
    3062         domain.forcing_terms = []       
     2918        domain.forcing_terms = []
    30632919        initial_volume = domain.quantities['stage'].get_integral()
    3064        
     2920
    30652921        assert num.allclose(initial_volume, width*length*stage)
    3066        
     2922
    30672923        for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
    3068             volume =  domain.quantities['stage'].get_integral()
    3069             assert num.allclose (volume, initial_volume)
    3070            
    3071            
     2924            volume = domain.quantities['stage'].get_integral()
     2925            assert num.allclose(volume, initial_volume)
     2926
    30722927        # Now apply the inflow and check volumes for a range of stage values
    30732928        for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]:
    30742929            domain.time = 0.0
    3075             domain.set_quantity('stage', stage)       
    3076                    
    3077             domain.forcing_terms = []       
    3078             domain.forcing_terms.append(Inflow(domain, rate=2.0, center=(5,5), radius=1))       
     2930            domain.set_quantity('stage', stage)
     2931            domain.forcing_terms = []
     2932            domain.forcing_terms.append(Inflow(domain, rate=2.0,
     2933                                               center=(5,5), radius=1))
    30792934            initial_volume = domain.quantities['stage'].get_integral()
    30802935            predicted_volume = initial_volume
    30812936            dt = 0.05
    3082             for t in domain.evolve(yieldstep = dt, finaltime = 5.0):
     2937            for t in domain.evolve(yieldstep=dt, finaltime=5.0):
    30832938                volume = domain.quantities['stage'].get_integral()
    3084                
    3085                 assert num.allclose (volume, predicted_volume)           
     2939                assert num.allclose (volume, predicted_volume)
    30862940                predicted_volume = predicted_volume + 2.0/pi/100/dt # Why 100?
    3087            
    3088            
    3089         # Apply equivalent outflow only and check volumes for a range of stage values
     2941
     2942        # Apply equivalent outflow only and check volumes
     2943        # for a range of stage values
    30902944        for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]:
    30912945            print stage
    3092            
     2946
    30932947            domain.time = 0.0
    3094             domain.set_quantity('stage', stage)       
    3095             domain.forcing_terms = []       
    3096             domain.forcing_terms.append(Inflow(domain, rate=-2.0, center=(15,5), radius=1))       
     2948            domain.set_quantity('stage', stage)
     2949            domain.forcing_terms = []
     2950            domain.forcing_terms.append(Inflow(domain, rate=-2.0,
     2951                                               center=(15,5), radius=1))
    30972952            initial_volume = domain.quantities['stage'].get_integral()
    30982953            predicted_volume = initial_volume
    30992954            dt = 0.05
    3100             for t in domain.evolve(yieldstep = dt, finaltime = 5.0):
     2955            for t in domain.evolve(yieldstep=dt, finaltime=5.0):
    31012956                volume = domain.quantities['stage'].get_integral()
    3102                
    31032957                print t, volume, predicted_volume
    3104                 assert num.allclose (volume, predicted_volume)           
    3105                 predicted_volume = predicted_volume - 2.0/pi/100/dt # Why 100?           
    3106            
    3107            
    3108         # Apply both inflow and outflow and check volumes being constant for a
    3109         # range of stage values
    3110         for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]:       
     2958                assert num.allclose (volume, predicted_volume)
     2959                predicted_volume = predicted_volume - 2.0/pi/100/dt # Why 100?
     2960
     2961        # Apply both inflow and outflow and check volumes being constant for a
     2962        # range of stage values
     2963        for stage in [2.0, 1.0, 0.5, 0.25, 0.1, 0.0]:
    31112964            print stage
    3112            
     2965
    31132966            domain.time = 0.0
    3114             domain.set_quantity('stage', stage)       
    3115             domain.forcing_terms = []       
    3116             domain.forcing_terms.append(Inflow(domain, rate=2.0, center=(5,5), radius=1))       
    3117             domain.forcing_terms.append(Inflow(domain, rate=-2.0, center=(15,5), radius=1))               
     2967            domain.set_quantity('stage', stage)
     2968            domain.forcing_terms = []
     2969            domain.forcing_terms.append(Inflow(domain, rate=2.0,
     2970                                               center=(5,5), radius=1))
     2971            domain.forcing_terms.append(Inflow(domain, rate=-2.0,
     2972                                               center=(15,5), radius=1))
    31182973            initial_volume = domain.quantities['stage'].get_integral()
    31192974
    31202975            dt = 0.05
    3121             for t in domain.evolve(yieldstep = dt, finaltime = 5.0):
     2976            for t in domain.evolve(yieldstep=dt, finaltime=5.0):
    31222977                volume = domain.quantities['stage'].get_integral()
    3123                
     2978
    31242979                print t, volume
    3125                 assert num.allclose (volume, initial_volume)           
    3126 
    3127            
    3128            
     2980                assert num.allclose(volume, initial_volume)
    31292981
    31302982    #####################################################
     2983
    31312984    def test_first_order_extrapolator_const_z(self):
    3132 
    31332985        a = [0.0, 0.0]
    31342986        b = [0.0, 2.0]
     
    31392991
    31402992        points = [a, b, c, d, e, f]
    3141         #bac, bce, ecf, dbe
    3142         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     2993        #             bac,     bce,     ecf,    dbe
     2994        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    31432995
    31442996        domain = Domain(points, vertices)
    3145         val0 = 2.+2.0/3
    3146         val1 = 4.+4.0/3
    3147         val2 = 8.+2.0/3
    3148         val3 = 2.+8.0/3
    3149 
    3150         zl=zr=-3.75 #Assume constant bed (must be less than stage)
    3151         domain.set_quantity('elevation', zl*num.ones( (4,3), num.int )) #array default#
     2997        val0 = 2. + 2.0/3
     2998        val1 = 4. + 4.0/3
     2999        val2 = 8. + 2.0/3
     3000        val3 = 2. + 8.0/3
     3001
     3002        zl = zr = -3.75    # Assume constant bed (must be less than stage)
     3003        domain.set_quantity('elevation', zl*num.ones((4, 3), num.int)) #array default#
    31523004        domain.set_quantity('stage', [[val0, val0-1, val0-2],
    31533005                                      [val1, val1+1, val1],
     
    31553007                                      [val3-0.5, val3, val3]])
    31563008
    3157 
    3158 
    31593009        domain._order_ = 1
    31603010        domain.distribute_to_vertices_and_edges()
     
    31633013        C = domain.quantities['stage'].centroid_values
    31643014        for i in range(3):
    3165             assert num.allclose( domain.quantities['stage'].vertex_values[:,i], C)
    3166 
     3015            assert num.allclose(domain.quantities['stage'].vertex_values[:,i],
     3016                                C)
    31673017
    31683018    def test_first_order_limiter_variable_z(self):
    3169         #Check that first order limiter follows bed_slope
     3019        '''Check that first order limiter follows bed_slope'''
     3020
    31703021        from anuga.config import epsilon
    31713022
    31723023        a = [0.0, 0.0]
    31733024        b = [0.0, 2.0]
    3174         c = [2.0,0.0]
     3025        c = [2.0, 0.0]
    31753026        d = [0.0, 4.0]
    31763027        e = [2.0, 2.0]
    3177         f = [4.0,0.0]
     3028        f = [4.0, 0.0]
    31783029
    31793030        points = [a, b, c, d, e, f]
    3180         #bac, bce, ecf, dbe
    3181         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     3031        #             bac,     bce,     ecf,    dbe
     3032        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    31823033
    31833034        domain = Domain(points, vertices)
    3184         val0 = 2.+2.0/3
    3185         val1 = 4.+4.0/3
    3186         val2 = 8.+2.0/3
    3187         val3 = 2.+8.0/3
     3035        val0 = 2. + 2.0/3
     3036        val1 = 4. + 4.0/3
     3037        val2 = 8. + 2.0/3
     3038        val3 = 2. + 8.0/3
    31883039
    31893040        domain.set_quantity('elevation', [[0,0,0], [6,0,0],
     
    31973048        L = domain.quantities['stage'].vertex_values
    31983049
    3199 
    3200         #Check that some stages are not above elevation (within eps)
    3201         #- so that the limiter has something to work with
    3202         assert not num.alltrue(num.alltrue(num.greater_equal(L,E-epsilon)))
     3050        # Check that some stages are not above elevation (within eps) -
     3051        # so that the limiter has something to work with
     3052        assert not num.alltrue(num.alltrue(num.greater_equal(L, E-epsilon)))
    32033053
    32043054        domain._order_ = 1
     
    32063056
    32073057        #Check that all stages are above elevation (within eps)
    3208         assert num.alltrue(num.alltrue(num.greater_equal(L,E-epsilon)))
    3209 
     3058        assert num.alltrue(num.alltrue(num.greater_equal(L, E-epsilon)))
    32103059
    32113060    #####################################################
     3061
    32123062    def test_distribute_basic(self):
    32133063        #Using test data generated by abstract_2d_finite_volumes-2
     
    32223072
    32233073        points = [a, b, c, d, e, f]
    3224         #bac, bce, ecf, dbe
    3225         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     3074        #             bac,     bce,     ecf,    dbe
     3075        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    32263076
    32273077        domain = Domain(points, vertices)
     
    32363086        L = domain.quantities['stage'].vertex_values
    32373087
    3238         #First order
     3088        # First order
    32393089        domain._order_ = 1
    32403090        domain.distribute_to_vertices_and_edges()
    32413091        assert num.allclose(L[1], val1)
    32423092
    3243         #Second order
     3093        # Second order
    32443094        domain._order_ = 2
    3245         domain.beta_w      = 0.9
    3246         domain.beta_w_dry  = 0.9
    3247         domain.beta_uh     = 0.9
     3095        domain.beta_w = 0.9
     3096        domain.beta_w_dry = 0.9
     3097        domain.beta_uh = 0.9
    32483098        domain.beta_uh_dry = 0.9
    3249         domain.beta_vh     = 0.9
     3099        domain.beta_vh = 0.9
    32503100        domain.beta_vh_dry = 0.9
    32513101        domain.distribute_to_vertices_and_edges()
    32523102        assert num.allclose(L[1], [2.2, 4.9, 4.9])
    3253 
    3254 
    32553103
    32563104    def test_distribute_away_from_bed(self):
     
    32663114
    32673115        points = [a, b, c, d, e, f]
    3268         #bac, bce, ecf, dbe
    3269         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     3116        #             bac,     bce,     ecf,    dbe
     3117        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    32703118
    32713119        domain = Domain(points, vertices)
    32723120        L = domain.quantities['stage'].vertex_values
    32733121
    3274         def stage(x,y):
     3122        def stage(x, y):
    32753123            return x**2
    32763124
     
    32803128
    32813129        a, b = domain.quantities['stage'].get_gradients()
    3282                
     3130
    32833131        assert num.allclose(a[1], 3.33333334)
    32843132        assert num.allclose(b[1], 0.0)
     
    32893137
    32903138        domain._order_ = 2
    3291         domain.beta_w      = 0.9
    3292         domain.beta_w_dry  = 0.9
    3293         domain.beta_uh     = 0.9
     3139        domain.beta_w = 0.9
     3140        domain.beta_w_dry = 0.9
     3141        domain.beta_uh = 0.9
    32943142        domain.beta_uh_dry = 0.9
    3295         domain.beta_vh     = 0.9
     3143        domain.beta_vh = 0.9
    32963144        domain.beta_vh_dry = 0.9
    32973145        domain.distribute_to_vertices_and_edges()
    32983146        assert num.allclose(L[1], [0.57777777, 2.37777778, 2.37777778])
    3299 
    3300 
    33013147
    33023148    def test_distribute_away_from_bed1(self):
     
    33123158
    33133159        points = [a, b, c, d, e, f]
    3314         #bac, bce, ecf, dbe
    3315         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     3160        #             bac,     bce,     ecf,    dbe
     3161        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    33163162
    33173163        domain = Domain(points, vertices)
    33183164        L = domain.quantities['stage'].vertex_values
    33193165
    3320         def stage(x,y):
    3321             return x**4+y**2
     3166        def stage(x, y):
     3167            return x**4 + y**2
    33223168
    33233169        domain.set_quantity('stage', stage, location='centroids')
    3324         #print domain.quantities['stage'].centroid_values
    33253170
    33263171        domain.quantities['stage'].compute_gradients()
     
    33343179
    33353180        domain._order_ = 2
    3336         domain.beta_w      = 0.9
    3337         domain.beta_w_dry  = 0.9
    3338         domain.beta_uh     = 0.9
     3181        domain.beta_w = 0.9
     3182        domain.beta_w_dry = 0.9
     3183        domain.beta_uh = 0.9
    33393184        domain.beta_uh_dry = 0.9
    3340         domain.beta_vh     = 0.9
     3185        domain.beta_vh = 0.9
    33413186        domain.beta_vh_dry = 0.9
    33423187        domain.distribute_to_vertices_and_edges()
    33433188        assert num.allclose(L[1], [1.07160494, 6.46058131, 7.28262855])
    33443189
    3345 
    3346 
    33473190    def test_distribute_near_bed(self):
    3348 
    33493191        a = [0.0, 0.0]
    33503192        b = [0.0, 2.0]
     
    33553197
    33563198        points = [a, b, c, d, e, f]
    3357         #bac, bce, ecf, dbe
    3358         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     3199        #             bac,     bce,     ecf,    dbe
     3200        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    33593201
    33603202        domain = Domain(points, vertices)
    33613203
    3362 
    3363         #Set up for a gradient of (10,0) at mid triangle (bce)
     3204        # Set up for a gradient of (10,0) at mid triangle (bce)
    33643205        def slope(x, y):
    33653206            return 10*x
     
    33713212        domain.set_quantity('elevation', slope)
    33723213        domain.set_quantity('stage', stage, location='centroids')
    3373 
    3374         #print domain.quantities['elevation'].centroid_values
    3375         #print domain.quantities['stage'].centroid_values
    33763214
    33773215        E = domain.quantities['elevation'].vertex_values
     
    33823220        for i in range(len(L)):
    33833221            volumes.append(num.sum(L[i])/3)
    3384             assert num.allclose(volumes[i], domain.quantities['stage'].centroid_values[i]) 
    3385        
    3386        
     3222            assert num.allclose(volumes[i],
     3223                                domain.quantities['stage'].centroid_values[i])
     3224
    33873225        domain._order_ = 1
    3388        
     3226
    33893227        domain.tight_slope_limiters = 0
    33903228        domain.distribute_to_vertices_and_edges()
    33913229        assert num.allclose(L[1], [0.1, 20.1, 20.1])
    33923230        for i in range(len(L)):
    3393             assert num.allclose(volumes[i], num.sum(L[i])/3)                   
    3394        
    3395         domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
     3231            assert num.allclose(volumes[i], num.sum(L[i])/3)
     3232
     3233        # Allow triangle to be flatter (closer to bed)
     3234        domain.tight_slope_limiters = 1
     3235
    33963236        domain.distribute_to_vertices_and_edges()
    33973237        assert num.allclose(L[1], [0.298, 20.001, 20.001])
    33983238        for i in range(len(L)):
    3399             assert num.allclose(volumes[i], num.sum(L[i])/3)   
     3239            assert num.allclose(volumes[i], num.sum(L[i])/3)
    34003240
    34013241        domain._order_ = 2
    3402        
     3242
    34033243        domain.tight_slope_limiters = 0
    34043244        domain.distribute_to_vertices_and_edges()
    3405         assert num.allclose(L[1], [0.1, 20.1, 20.1])       
     3245        assert num.allclose(L[1], [0.1, 20.1, 20.1])
    34063246        for i in range(len(L)):
    3407             assert num.allclose(volumes[i], num.sum(L[i])/3)           
    3408        
    3409         domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
     3247            assert num.allclose(volumes[i], num.sum(L[i])/3)
     3248
     3249        # Allow triangle to be flatter (closer to bed)
     3250        domain.tight_slope_limiters = 1
     3251
    34103252        domain.distribute_to_vertices_and_edges()
    34113253        assert num.allclose(L[1], [0.298, 20.001, 20.001])
    34123254        for i in range(len(L)):
    3413             assert num.allclose(volumes[i], num.sum(L[i])/3)   
    3414        
    3415 
     3255            assert num.allclose(volumes[i], num.sum(L[i])/3)
    34163256
    34173257    def test_distribute_near_bed1(self):
    3418 
    34193258        a = [0.0, 0.0]
    34203259        b = [0.0, 2.0]
     
    34253264
    34263265        points = [a, b, c, d, e, f]
    3427         #bac, bce, ecf, dbe
    3428         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     3266        #             bac,     bce,     ecf,    dbe
     3267        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    34293268
    34303269        domain = Domain(points, vertices)
    34313270
    3432 
    3433         #Set up for a gradient of (8,2) at mid triangle (bce)
     3271        # Set up for a gradient of (8,2) at mid triangle (bce)
    34343272        def slope(x, y):
    3435             return x**4+y**2
     3273            return x**4 + y**2
    34363274
    34373275        h = 0.1
    3438         def stage(x,y):
    3439             return slope(x,y)+h
     3276        def stage(x, y):
     3277            return slope(x, y) + h
    34403278
    34413279        domain.set_quantity('elevation', slope)
    34423280        domain.set_quantity('stage', stage)
    3443 
    3444         #print domain.quantities['elevation'].centroid_values
    3445         #print domain.quantities['stage'].centroid_values
    34463281
    34473282        E = domain.quantities['elevation'].vertex_values
     
    34523287        for i in range(len(L)):
    34533288            volumes.append(num.sum(L[i])/3)
    3454             assert num.allclose(volumes[i], domain.quantities['stage'].centroid_values[i]) 
    3455        
    3456         #print E
     3289            assert num.allclose(volumes[i],
     3290                                domain.quantities['stage'].centroid_values[i])
     3291
    34573292        domain._order_ = 1
    3458        
     3293
    34593294        domain.tight_slope_limiters = 0
    34603295        domain.distribute_to_vertices_and_edges()
    3461         assert num.allclose(L[1], [4.1, 16.1, 20.1])       
     3296        assert num.allclose(L[1], [4.1, 16.1, 20.1])
    34623297        for i in range(len(L)):
    34633298            assert num.allclose(volumes[i], num.sum(L[i])/3)
    3464        
    3465                
    3466         domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
     3299
     3300        # Allow triangle to be flatter (closer to bed)
     3301        domain.tight_slope_limiters = 1
     3302
    34673303        domain.distribute_to_vertices_and_edges()
    34683304        assert num.allclose(L[1], [4.2386, 16.0604, 20.001])
    34693305        for i in range(len(L)):
    3470             assert num.allclose(volumes[i], num.sum(L[i])/3)   
    3471        
     3306            assert num.allclose(volumes[i], num.sum(L[i])/3)
    34723307
    34733308        domain._order_ = 2
    3474        
    3475         domain.tight_slope_limiters = 0   
     3309
     3310        domain.tight_slope_limiters = 0
    34763311        domain.distribute_to_vertices_and_edges()
    34773312        assert num.allclose(L[1], [4.1, 16.1, 20.1])
    34783313        for i in range(len(L)):
    3479             assert num.allclose(volumes[i], num.sum(L[i])/3)   
    3480        
    3481         domain.tight_slope_limiters = 1 # Allow triangle to be flatter (closer to bed)
     3314            assert num.allclose(volumes[i], num.sum(L[i])/3)
     3315
     3316        # Allow triangle to be flatter (closer to bed)
     3317        domain.tight_slope_limiters = 1
     3318
    34823319        domain.distribute_to_vertices_and_edges()
    3483         #print L[1]
    3484         assert num.allclose(L[1], [4.23370103, 16.06529897, 20.001]) or\
    3485                num.allclose(L[1], [4.18944138, 16.10955862, 20.001]) or\
    3486                num.allclose(L[1], [4.19351461, 16.10548539, 20.001]) # old limiters
    3487        
     3320        assert (num.allclose(L[1], [4.23370103, 16.06529897, 20.001]) or
     3321                num.allclose(L[1], [4.18944138, 16.10955862, 20.001]) or
     3322                num.allclose(L[1], [4.19351461, 16.10548539, 20.001]))
     3323        # old limiters
     3324
    34883325        for i in range(len(L)):
    34893326            assert num.allclose(volumes[i], num.sum(L[i])/3)
    3490 
    34913327
    34923328    def test_second_order_distribute_real_data(self):
     
    35033339
    35043340        points = [a, b, c, d, e, f, g]
    3505         #bae, efb, cbf, feg
    3506         vertices = [ [1,0,4], [4,5,1], [2,1,5], [5,4,6]]
     3341        #             bae,     efb,     cbf,    feg
     3342        vertices = [[1,0,4], [4,5,1], [2,1,5], [5,4,6]]
    35073343
    35083344        domain = Domain(points, vertices)
     
    35303366        Y = domain.quantities['ymomentum'].vertex_values
    35313367
    3532         #print E
    35333368        domain._order_ = 2
    3534         domain.beta_w      = 0.9
    3535         domain.beta_w_dry  = 0.9
    3536         domain.beta_uh     = 0.9
     3369        domain.beta_w = 0.9
     3370        domain.beta_w_dry = 0.9
     3371        domain.beta_uh = 0.9
    35373372        domain.beta_uh_dry = 0.9
    3538         domain.beta_vh     = 0.9
     3373        domain.beta_vh = 0.9
    35393374        domain.beta_vh_dry = 0.9
    3540        
     3375
    35413376        # FIXME (Ole): Need tests where this is commented out
    3542         domain.tight_slope_limiters = 0 # Backwards compatibility (14/4/7)                 
    3543         domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
    3544        
    3545                
     3377        domain.tight_slope_limiters = 0       # Backwards compatibility (14/4/7)
     3378        domain.use_centroid_velocities = 0    # Backwards compatibility (7/5/8)
     3379
    35463380        domain.distribute_to_vertices_and_edges()
    3547 
    3548         #print L[1,:]
    3549         #print X[1,:]
    3550         #print Y[1,:]
    35513381
    35523382        assert num.allclose(L[1,:], [-0.00825735775384,
     
    35603390                                     -0.000151505429018])
    35613391
    3562 
    3563 
    35643392    def test_balance_deep_and_shallow(self):
    35653393        """Test that balanced limiters preserve conserved quantites.
    35663394        This test is using old depth based balanced limiters
    35673395        """
     3396
    35683397        import copy
    35693398
     
    35763405
    35773406        points = [a, b, c, d, e, f]
    3578 
    3579         #bac, bce, ecf, dbe
    3580         elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
     3407        #             bac,     bce,     ecf,     dbe
     3408        elements = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
    35813409
    35823410        domain = Domain(points, elements)
    35833411        domain.check_integrity()
    35843412
    3585         #Create a deliberate overshoot
     3413        # Create a deliberate overshoot
    35863414        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
    3587         domain.set_quantity('elevation', 0) #Flat bed
     3415        domain.set_quantity('elevation', 0)    # Flat bed
    35883416        stage = domain.quantities['stage']
    35893417
    3590         ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
    3591 
    3592         #Limit
    3593         domain.tight_slope_limiters = 0               
     3418        ref_centroid_values = copy.copy(stage.centroid_values[:])    # Copy
     3419
     3420        # Limit
     3421        domain.tight_slope_limiters = 0
    35943422        domain.distribute_to_vertices_and_edges()
    35953423
    3596         #Assert that quantities are conserved
     3424        # Assert that quantities are conserved
    35973425        for k in range(len(domain)):
    3598             assert num.allclose (ref_centroid_values[k],
    3599                                  num.sum(stage.vertex_values[k,:])/3)
    3600 
    3601 
    3602         #Now try with a non-flat bed - closely hugging initial stage in places
    3603         #This will create alphas in the range [0, 0.478260, 1]
     3426            assert num.allclose(ref_centroid_values[k],
     3427                                num.sum(stage.vertex_values[k,:])/3)
     3428
     3429        # Now try with a non-flat bed - closely hugging initial stage in places
     3430        # This will create alphas in the range [0, 0.478260, 1]
    36043431        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
    36053432        domain.set_quantity('elevation', [[0,0,0],
    3606                                         [1.8,1.9,5.9],
    3607                                         [4.6,0,0],
    3608                                         [0,2,4]])
     3433                                          [1.8,1.9,5.9],
     3434                                          [4.6,0,0],
     3435                                          [0,2,4]])
    36093436        stage = domain.quantities['stage']
    36103437
    3611         ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
    3612         ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
    3613 
    3614         #Limit
    3615         domain.tight_slope_limiters = 0       
     3438        ref_centroid_values = copy.copy(stage.centroid_values[:])    # Copy
     3439        ref_vertex_values = copy.copy(stage.vertex_values[:])        # Copy
     3440
     3441        # Limit
     3442        domain.tight_slope_limiters = 0
    36163443        domain.distribute_to_vertices_and_edges()
    36173444
    3618 
    3619         #Assert that all vertex quantities have changed
     3445        # Assert that all vertex quantities have changed
    36203446        for k in range(len(domain)):
    3621             #print ref_vertex_values[k,:], stage.vertex_values[k,:]
    3622             assert not num.allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])
    3623         #and assert that quantities are still conserved
     3447            assert not num.allclose(ref_vertex_values[k,:],
     3448                                    stage.vertex_values[k,:])
     3449        # and assert that quantities are still conserved
    36243450        for k in range(len(domain)):
    3625             assert num.allclose (ref_centroid_values[k],
    3626                                  num.sum(stage.vertex_values[k,:])/3)
    3627 
     3451            assert num.allclose(ref_centroid_values[k],
     3452                                num.sum(stage.vertex_values[k,:])/3)
    36283453
    36293454        # Check actual results
    3630         assert num.allclose (stage.vertex_values,
    3631                              [[2,2,2],
    3632                               [1.93333333, 2.03333333, 6.03333333],
    3633                               [6.93333333, 4.53333333, 4.53333333],
    3634                               [5.33333333, 5.33333333, 5.33333333]])
    3635 
     3455        assert num.allclose(stage.vertex_values,
     3456                            [[2,2,2],
     3457                             [1.93333333, 2.03333333, 6.03333333],
     3458                             [6.93333333, 4.53333333, 4.53333333],
     3459                             [5.33333333, 5.33333333, 5.33333333]])
    36363460
    36373461    def test_balance_deep_and_shallow_tight_SL(self):
     
    36393463        This test is using Tight Slope Limiters
    36403464        """
     3465
    36413466        import copy
    36423467
     
    36493474
    36503475        points = [a, b, c, d, e, f]
    3651 
    3652         #bac, bce, ecf, dbe
    3653         elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
     3476        #             bac,     bce,     ecf,     dbe
     3477        elements = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
    36543478
    36553479        domain = Domain(points, elements)
    36563480        domain.check_integrity()
    36573481
    3658         #Create a deliberate overshoot
     3482        # Create a deliberate overshoot
    36593483        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
    3660         domain.set_quantity('elevation', 0) #Flat bed
     3484        domain.set_quantity('elevation', 0)    # Flat bed
    36613485        stage = domain.quantities['stage']
    36623486
    3663         ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
    3664 
    3665         #Limit
    3666         domain.tight_slope_limiters = 1               
     3487        ref_centroid_values = copy.copy(stage.centroid_values[:])    # Copy
     3488
     3489        # Limit
     3490        domain.tight_slope_limiters = 1
    36673491        domain.distribute_to_vertices_and_edges()
    36683492
    3669         #Assert that quantities are conserved
     3493        # Assert that quantities are conserved
    36703494        for k in range(len(domain)):
    36713495            assert num.allclose (ref_centroid_values[k],
    36723496                                 num.sum(stage.vertex_values[k,:])/3)
    36733497
    3674 
    3675         #Now try with a non-flat bed - closely hugging initial stage in places
    3676         #This will create alphas in the range [0, 0.478260, 1]
     3498        # Now try with a non-flat bed - closely hugging initial stage in places
     3499        # This will create alphas in the range [0, 0.478260, 1]
    36773500        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
    36783501        domain.set_quantity('elevation', [[0,0,0],
    3679                                         [1.8,1.9,5.9],
    3680                                         [4.6,0,0],
    3681                                         [0,2,4]])
     3502                                          [1.8,1.9,5.9],
     3503                                          [4.6,0,0],
     3504                                          [0,2,4]])
    36823505        stage = domain.quantities['stage']
    36833506
    3684         ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
    3685         ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
    3686 
    3687         #Limit
    3688         domain.tight_slope_limiters = 1       
     3507        ref_centroid_values = copy.copy(stage.centroid_values[:])    # Copy
     3508        ref_vertex_values = copy.copy(stage.vertex_values[:])        # Copy
     3509
     3510        # Limit
     3511        domain.tight_slope_limiters = 1
    36893512        domain.distribute_to_vertices_and_edges()
    36903513
    3691 
    3692         #Assert that all vertex quantities have changed
     3514        # Assert that all vertex quantities have changed
    36933515        for k in range(len(domain)):
    3694             #print ref_vertex_values[k,:], stage.vertex_values[k,:]
    3695             assert not num.allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])
    3696         #and assert that quantities are still conserved
     3516            assert not num.allclose(ref_vertex_values[k,:],
     3517                                    stage.vertex_values[k,:])
     3518        # and assert that quantities are still conserved
    36973519        for k in range(len(domain)):
    3698             assert num.allclose (ref_centroid_values[k],
    3699                                  num.sum(stage.vertex_values[k,:])/3)
    3700 
    3701 
    3702         #Also check that Python and C version produce the same
    3703         # No longer applicable if tight_slope_limiters == 1
    3704         #print stage.vertex_values
    3705         #assert allclose (stage.vertex_values,
    3706         #                 [[2,2,2],
    3707         #                  [1.93333333, 2.03333333, 6.03333333],
    3708         #                  [6.93333333, 4.53333333, 4.53333333],
    3709         #                  [5.33333333, 5.33333333, 5.33333333]])
    3710 
    3711 
     3520            assert num.allclose(ref_centroid_values[k],
     3521                                num.sum(stage.vertex_values[k,:])/3)
    37123522
    37133523    def test_balance_deep_and_shallow_Froude(self):
     
    37163526        This test is using tight slope limiters.
    37173527        """
     3528
    37183529        import copy
    37193530
     
    37263537
    37273538        points = [a, b, c, d, e, f]
    3728 
    3729         # bac, bce, ecf, dbe
    3730         elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
     3539        #             bac,     bce,     ecf,     dbe
     3540        elements = [[1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
    37313541
    37323542        domain = Domain(points, elements)
    37333543        domain.check_integrity()
    37343544        domain.tight_slope_limiters = True
    3735         domain.use_centroid_velocities = True               
     3545        domain.use_centroid_velocities = True
    37363546
    37373547        # Create non-flat bed - closely hugging initial stage in places
     
    37393549        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
    37403550        domain.set_quantity('elevation', [[0,0,0],
    3741                                         [1.8,1.999,5.999],
    3742                                         [4.6,0,0],
    3743                                         [0,2,4]])
     3551                                          [1.8,1.999,5.999],
     3552                                          [4.6,0,0],
     3553                                          [0,2,4]])
    37443554
    37453555        # Create small momenta, that nonetheless will generate large speeds
     
    37483558        domain.set_quantity('ymomentum', 0.0890)
    37493559
    3750 
    3751 
    3752        
    37533560        stage = domain.quantities['stage']
    37543561        elevation = domain.quantities['elevation']
    37553562        xmomentum = domain.quantities['xmomentum']
    3756         ymomentum = domain.quantities['ymomentum']       
     3563        ymomentum = domain.quantities['ymomentum']
    37573564
    37583565        # Setup triangle #1 to mimick real Froude explosion observed
    37593566        # in the Onslow example 13 Nov 2007.
    3760 
    37613567        stage.vertex_values[1,:] = [1.6385, 1.6361, 1.2953]
    3762         elevation.vertex_values[1,:] = [1.6375, 1.6336, 0.4647]       
     3568        elevation.vertex_values[1,:] = [1.6375, 1.6336, 0.4647]
    37633569        xmomentum.vertex_values[1,:] = [-0.0058, -0.0050, -0.0066]
    37643570        ymomentum.vertex_values[1,:] = [0.0890, 0.0890, 0.0890]
    37653571
    37663572        xmomentum.interpolate()
    3767         ymomentum.interpolate()       
    3768         stage.interpolate()       
     3573        ymomentum.interpolate()
     3574        stage.interpolate()
    37693575        elevation.interpolate()
    37703576
     
    37803586        v = ymomentum/depth
    37813587
    3782         denom = (depth*g)**0.5 
     3588        denom = (depth*g)**0.5
    37833589        Fx = u/denom
    37843590        Fy = v/denom
    3785        
    3786    
     3591
    37873592        # Verify against Onslow example (14 Nov 2007)
    37883593        assert num.allclose(depth.centroid_values[1], 0.278033)
     
    38023607        assert num.allclose(Fy.centroid_values[1], 0.193924048435)
    38033608
    3804 
    38053609        # But Froude numbers are huge at some vertices and edges
    38063610        assert num.allclose(Fx.vertex_values[1,:], [-5.85888475e+01,
     
    38193623                                                  1.06035244e-01,
    38203624                                                  3.88346947e+02])
    3821        
    3822        
     3625
     3626
    38233627        # The task is now to arrange the limiters such that Froude numbers
    38243628        # remain under control whil at the same time obeying the conservation
    38253629        # laws.
    3826 
    3827        
    3828         ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
    3829         ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
     3630        ref_centroid_values = copy.copy(stage.centroid_values[:])    # Copy
     3631        ref_vertex_values = copy.copy(stage.vertex_values[:])        # Copy
    38303632
    38313633        # Limit (and invoke balance_deep_and_shallow)
     
    38343636
    38353637        # Redo derived quantities
    3836         depth = stage-elevation
     3638        depth = stage - elevation
    38373639        u = xmomentum/depth
    38383640        v = ymomentum/depth
     
    38403642        # Assert that all vertex velocities stay within one
    38413643        # order of magnitude of centroid velocities.
    3842         #print u.vertex_values[1,:]
    3843         #print u.centroid_values[1]
    3844        
    3845         assert num.alltrue(num.absolute(u.vertex_values[1,:]) <= num.absolute(u.centroid_values[1])*10)
    3846         assert num.alltrue(num.absolute(v.vertex_values[1,:]) <= num.absolute(v.centroid_values[1])*10)
    3847 
    3848         denom = (depth*g)**0.5
     3644        assert num.alltrue(num.absolute(u.vertex_values[1,:]) <=
     3645                           num.absolute(u.centroid_values[1])*10)
     3646        assert num.alltrue(num.absolute(v.vertex_values[1,:]) <=
     3647                           num.absolute(v.centroid_values[1])*10)
     3648
     3649        denom = (depth*g)**0.5
    38493650        Fx = u/denom
    38503651        Fy = v/denom
    3851 
    38523652
    38533653        # Assert that Froude numbers are less than max value (TBA)
    38543654        # at vertices, edges and centroids.
    38553655        from anuga.config import maximum_froude_number
    3856         assert num.alltrue(num.absolute(Fx.vertex_values[1,:]) < maximum_froude_number)
    3857         assert num.alltrue(num.absolute(Fy.vertex_values[1,:]) < maximum_froude_number)
    3858 
     3656
     3657        assert num.alltrue(num.absolute(Fx.vertex_values[1,:]) <
     3658                           maximum_froude_number)
     3659        assert num.alltrue(num.absolute(Fy.vertex_values[1,:]) <
     3660                           maximum_froude_number)
    38593661
    38603662        # Assert that all vertex quantities have changed
    38613663        for k in range(len(domain)):
    3862             #print ref_vertex_values[k,:], stage.vertex_values[k,:]
    3863             assert not num.allclose (ref_vertex_values[k,:],
    3864                                      stage.vertex_values[k,:])
    3865            
     3664            assert not num.allclose(ref_vertex_values[k,:],
     3665                                    stage.vertex_values[k,:])
     3666
    38663667        # Assert that quantities are still conserved
    38673668        for k in range(len(domain)):
    3868             assert num.allclose (ref_centroid_values[k],
    3869                                  num.sum(stage.vertex_values[k,:])/3)
    3870 
    3871 
    3872        
     3669            assert num.allclose(ref_centroid_values[k],
     3670                                num.sum(stage.vertex_values[k,:])/3)
     3671
    38733672        return
    3874    
     3673
    38753674        qwidth = 12
    3876         for k in [1]: #range(len(domain)):
    3877             print 'Triangle %d (C, V, E)' %k
    3878            
    3879             print 'stage'.ljust(qwidth), stage.centroid_values[k],\
    3880                   stage.vertex_values[k,:], stage.edge_values[k,:]
    3881             print 'elevation'.ljust(qwidth), elevation.centroid_values[k],\
    3882                   elevation.vertex_values[k,:], elevation.edge_values[k,:]
    3883             print 'depth'.ljust(qwidth), depth.centroid_values[k],\
    3884                   depth.vertex_values[k,:], depth.edge_values[k,:]
    3885             print 'xmomentum'.ljust(qwidth), xmomentum.centroid_values[k],\
    3886                   xmomentum.vertex_values[k,:], xmomentum.edge_values[k,:]
    3887             print 'ymomentum'.ljust(qwidth), ymomentum.centroid_values[k],\
    3888                   ymomentum.vertex_values[k,:], ymomentum.edge_values[k,:]
    3889             print 'u'.ljust(qwidth),u.centroid_values[k],\
    3890                   u.vertex_values[k,:], u.edge_values[k,:]
    3891             print 'v'.ljust(qwidth), v.centroid_values[k],\
    3892                   v.vertex_values[k,:], v.edge_values[k,:]
    3893             print 'Fx'.ljust(qwidth), Fx.centroid_values[k],\
    3894                   Fx.vertex_values[k,:], Fx.edge_values[k,:]
    3895             print 'Fy'.ljust(qwidth), Fy.centroid_values[k],\
    3896                   Fy.vertex_values[k,:], Fy.edge_values[k,:]
    3897            
    3898            
    3899        
    3900 
    3901 
     3675        for k in [1]:    # range(len(domain)):
     3676            print 'Triangle %d (C, V, E)' % k
     3677
     3678            print ('stage'.ljust(qwidth), stage.centroid_values[k],
     3679                   stage.vertex_values[k,:], stage.edge_values[k,:])
     3680            print ('elevation'.ljust(qwidth), elevation.centroid_values[k],
     3681                   elevation.vertex_values[k,:], elevation.edge_values[k,:])
     3682            print ('depth'.ljust(qwidth), depth.centroid_values[k],
     3683                   depth.vertex_values[k,:], depth.edge_values[k,:])
     3684            print ('xmomentum'.ljust(qwidth), xmomentum.centroid_values[k],
     3685                   xmomentum.vertex_values[k,:], xmomentum.edge_values[k,:])
     3686            print ('ymomentum'.ljust(qwidth), ymomentum.centroid_values[k],
     3687                   ymomentum.vertex_values[k,:], ymomentum.edge_values[k,:])
     3688            print ('u'.ljust(qwidth),u.centroid_values[k],
     3689                   u.vertex_values[k,:], u.edge_values[k,:])
     3690            print ('v'.ljust(qwidth), v.centroid_values[k],
     3691                   v.vertex_values[k,:], v.edge_values[k,:])
     3692            print ('Fx'.ljust(qwidth), Fx.centroid_values[k],
     3693                   Fx.vertex_values[k,:], Fx.edge_values[k,:])
     3694            print ('Fy'.ljust(qwidth), Fy.centroid_values[k],
     3695                   Fy.vertex_values[k,:], Fy.edge_values[k,:])
    39023696
    39033697    def test_conservation_1(self):
     
    39073701        initial condition
    39083702        """
     3703
    39093704        from mesh_factory import rectangular
    39103705
    3911         #Create basic mesh
     3706        # Create basic mesh
    39123707        points, vertices, boundary = rectangular(6, 6)
    39133708
    3914         #Create shallow water domain
     3709        # Create shallow water domain
    39153710        domain = Domain(points, vertices, boundary)
    39163711        domain.smooth = False
    3917         domain.default_order=2
    3918 
    3919         #IC
     3712        domain.default_order = 2
     3713
     3714        # IC
    39203715        def x_slope(x, y):
    39213716            return x/3
     
    39343729        initial_xmom = domain.quantities['xmomentum'].get_integral()
    39353730
    3936         #print initial_xmom
    3937 
    3938         #Evolution
    3939         for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
    3940             volume =  domain.quantities['stage'].get_integral()
    3941             assert num.allclose (volume, initial_volume)
     3731        # Evolution
     3732        for t in domain.evolve(yieldstep=0.05, finaltime=5.0):
     3733            volume = domain.quantities['stage'].get_integral()
     3734            assert num.allclose(volume, initial_volume)
    39423735
    39433736            #I don't believe that the total momentum should be the same
     
    39493742        os.remove(domain.get_name() + '.sww')
    39503743
    3951 
    39523744    def test_conservation_2(self):
    39533745        """Test that stage is conserved globally
     
    39563748        initial condition
    39573749        """
     3750
    39583751        from mesh_factory import rectangular
    39593752
    3960         #Create basic mesh
     3753        # Create basic mesh
    39613754        points, vertices, boundary = rectangular(6, 6)
    39623755
    3963         #Create shallow water domain
     3756        # Create shallow water domain
    39643757        domain = Domain(points, vertices, boundary)
    39653758        domain.smooth = False
    3966         domain.default_order=2
    3967 
    3968         #IC
     3759        domain.default_order = 2
     3760
     3761        # IC
    39693762        def x_slope(x, y):
    39703763            return x/3
     
    39723765        domain.set_quantity('elevation', x_slope)
    39733766        domain.set_quantity('friction', 0)
    3974         domain.set_quantity('stage', 0.4) #Steady
     3767        domain.set_quantity('stage', 0.4)    # Steady
    39753768
    39763769        # Boundary conditions (reflective everywhere)
     
    39833776        initial_xmom = domain.quantities['xmomentum'].get_integral()
    39843777
    3985         #print initial_xmom
    3986 
    3987         #Evolution
    3988         for t in domain.evolve(yieldstep = 0.05, finaltime = 5.0):
    3989             volume =  domain.quantities['stage'].get_integral()
    3990             assert num.allclose (volume, initial_volume)
     3778        # Evolution
     3779        for t in domain.evolve(yieldstep=0.05, finaltime=5.0):
     3780            volume = domain.quantities['stage'].get_integral()
     3781            assert num.allclose(volume, initial_volume)
    39913782
    39923783            #FIXME: What would we expect from momentum
     
    40003791        """Test that stage is conserved globally
    40013792
    4002         This one uses a larger grid, convoluted bed, reflective bdries and a suitable
    4003         initial condition
     3793        This one uses a larger grid, convoluted bed, reflective boundaries
     3794        and a suitable initial condition
    40043795        """
     3796
    40053797        from mesh_factory import rectangular
    40063798
    4007         #Create basic mesh
     3799        # Create basic mesh
    40083800        points, vertices, boundary = rectangular(2, 1)
    40093801
    4010         #Create shallow water domain
     3802        # Create shallow water domain
    40113803        domain = Domain(points, vertices, boundary)
    40123804        domain.smooth = False
     
    40143806        domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
    40153807
    4016         #IC
     3808        # IC
    40173809        def x_slope(x, y):
    40183810            z = 0*x
     
    40283820            return z
    40293821
    4030 
    4031 
    40323822        domain.set_quantity('elevation', x_slope)
    40333823        domain.set_quantity('friction', 0)
     
    40443834
    40453835        import copy
    4046         ref_centroid_values =\
    4047                  copy.copy(domain.quantities['stage'].centroid_values)
    4048 
    4049         #print 'ORG', domain.quantities['stage'].centroid_values
     3836
     3837        ref_centroid_values = copy.copy(domain.quantities['stage'].\
     3838                                            centroid_values)
     3839
    40503840        domain.distribute_to_vertices_and_edges()
    40513841
    4052 
    4053         #print domain.quantities['stage'].centroid_values
    40543842        assert num.allclose(domain.quantities['stage'].centroid_values,
    40553843                            ref_centroid_values)
    40563844
    4057 
    4058         #Check that initial limiter doesn't violate cons quan
     3845        # Check that initial limiter doesn't violate cons quan
    40593846        assert num.allclose(domain.quantities['stage'].get_integral(),
    40603847                            initial_volume)
    40613848
    4062         #Evolution
    4063         for t in domain.evolve(yieldstep = 0.05, finaltime = 10):
     3849        # Evolution
     3850        for t in domain.evolve(yieldstep=0.05, finaltime=10):
    40643851            volume =  domain.quantities['stage'].get_integral()
    4065             #print t, volume, initial_volume
    40663852            assert num.allclose (volume, initial_volume)
    40673853
     
    40713857        """Test that stage is conserved globally
    40723858
    4073         This one uses a larger grid, convoluted bed, reflective bdries and a suitable
    4074         initial condition
     3859        This one uses a larger grid, convoluted bed, reflective boundaries
     3860        and a suitable initial condition
    40753861        """
     3862
    40763863        from mesh_factory import rectangular
    40773864
    4078         #Create basic mesh
     3865        # Create basic mesh
    40793866        points, vertices, boundary = rectangular(6, 6)
    40803867
    4081         #Create shallow water domain
     3868        # Create shallow water domain
    40823869        domain = Domain(points, vertices, boundary)
    40833870        domain.smooth = False
    4084         domain.default_order=2
    4085 
    4086         #IC
     3871        domain.default_order = 2
     3872
     3873        # IC
    40873874        def x_slope(x, y):
    40883875            z = 0*x
     
    40933880                    z[i] = -0.5
    40943881                if 0.5 <= x[i] < 0.7:
    4095                     #z[i] = 0.3 #OK with beta == 0.2
    4096                     z[i] = 0.34 #OK with beta == 0.0
    4097                     #z[i] = 0.35#Fails after 80 timesteps with an error
    4098                                 #of the order 1.0e-5
     3882                    #z[i] = 0.3     # OK with beta == 0.2
     3883                    z[i] = 0.34     # OK with beta == 0.0
     3884                    #z[i] = 0.35    # Fails after 80 timesteps with an error
     3885                                    # of the order 1.0e-5
    40993886                if 0.7 <= x[i]:
    41003887                    z[i] = x[i]/3
    41013888            return z
    41023889
    4103 
    4104 
    41053890        domain.set_quantity('elevation', x_slope)
    41063891        domain.set_quantity('friction', 0)
     
    41173902
    41183903        import copy
    4119         ref_centroid_values =\
    4120                  copy.copy(domain.quantities['stage'].centroid_values)
    4121 
    4122         #Test limiter by itself