Changeset 7168


Ignore:
Timestamp:
Jun 10, 2009, 12:07:25 PM (13 years ago)
Author:
rwilson
Message:

Back-merge from Numeric.

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

Legend:

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

    r7035 r7168  
    233233    def test_sww_range(self):
    234234        """Test that constant sww information can be written correctly
    235         (non smooth)
     235        Use non-smooth to be able to compare to quantity values.
    236236        """
     237
    237238        self.domain.set_name('datatest' + str(id(self)))
    238239        self.domain.format = 'sww'
    239         self.domain.smooth = True
    240 
    241         self.domain.tight_slope_limiters = 0 # Backwards compatibility
    242         self.domain.use_centroid_velocities = 0 # Backwards compatibility (7/5/8)
    243        
     240        self.domain.set_store_vertices_uniquely(True)
    244241       
    245242        sww = get_dataobject(self.domain)       
    246243
     244        dqs = self.domain.get_quantity('stage')
     245        dqx = self.domain.get_quantity('xmomentum')
     246        dqy = self.domain.get_quantity('ymomentum')       
     247        xmom_min = ymom_min = stage_min = sys.maxint
     248        xmom_max = ymom_max = stage_max = -stage_min       
    247249        for t in self.domain.evolve(yieldstep = 1, finaltime = 1):
    248             pass
     250            wmax = max(dqs.get_values().flatten())
     251            if wmax > stage_max: stage_max = wmax
     252            wmin = min(dqs.get_values().flatten())
     253            if wmin < stage_min: stage_min = wmin           
     254           
     255            uhmax = max(dqx.get_values().flatten())
     256            if uhmax > xmom_max: xmom_max = uhmax
     257            uhmin = min(dqx.get_values().flatten())
     258            if uhmin < xmom_min: xmom_min = uhmin                       
     259           
     260            vhmax = max(dqy.get_values().flatten())
     261            if vhmax > ymom_max: ymom_max = vhmax
     262            vhmin = min(dqy.get_values().flatten())
     263            if vhmin < ymom_min: ymom_min = vhmin                                   
     264           
     265           
    249266           
    250267        # Get NetCDF
     
    253270        # Get the variables
    254271        range = fid.variables['stage_range'][:]
     272        assert num.allclose(range,[stage_min, stage_max])
     273
     274        range = fid.variables['xmomentum_range'][:]
    255275        #print range
    256         assert num.allclose(range,[-0.93519, 0.15]) or\
    257                num.allclose(range,[-0.9352743, 0.15]) or\
    258                num.allclose(range,[-0.93522203, 0.15000001]) # Old slope limiters
    259        
    260         range = fid.variables['xmomentum_range'][:]
    261         ##print range
    262         assert num.allclose(range,[0,0.4695096]) or \
    263                num.allclose(range,[0,0.47790655]) or\
    264                num.allclose(range,[0,0.46941957]) or\
    265                num.allclose(range,[0,0.47769409])
    266 
     276        assert num.allclose(range, [xmom_min, xmom_max])
    267277       
    268278        range = fid.variables['ymomentum_range'][:]
    269         ##print range
    270         assert num.allclose(range,[0,0.02174380]) or\
    271                num.allclose(range,[0,0.02174439]) or\
    272                num.allclose(range,[0,0.02283983]) or\
    273                num.allclose(range,[0,0.0217342]) or\
    274                num.allclose(range,[0,0.0227564]) # Old slope limiters
     279        #print range
     280        assert num.allclose(range, [ymom_min, ymom_max])       
     281
     282
    275283       
    276284        fid.close()
     
    344352
    345353        extrema = fid.variables['xmomentum.extrema'][:]
    346         assert num.allclose(extrema,[-0.06062178, 0.47873023]) or num.allclose(extrema, [-0.06062178, 0.47847986])
     354        assert num.allclose(extrema,[-0.06062178, 0.47873023]) or\
     355            num.allclose(extrema, [-0.06062178, 0.47847986]) or\
     356            num.allclose(extrema, [-0.06062178, 0.47848481]) # 27/5/9           
    347357       
    348358        extrema = fid.variables['ymomentum.extrema'][:]
     
    1108211092        for t in range(t_end+1):
    1108311093            for i in range(3):
    11084                 assert num.allclose(f(t, i), [1, 2, 0])
     11094                assert num.allclose(f(t, i), [1, 2, 0], atol=1.0e-6)
    1108511095           
    1108611096
     
    1111911129
    1112011130        #print i, Q
    11121         assert num.allclose(Q, 0)       
     11131        assert num.allclose(Q, 0, atol=1.0e-5)       
    1112211132
    1112311133
     
    1121611226            for i in range(3):
    1121711227                #print i, t, f(t, i)           
    11218                 assert num.allclose(f(t, i), [w, uh, 0])
     11228                assert num.allclose(f(t, i), [w, uh, 0], atol=1.0e-6)
    1121911229           
    1122011230
     
    1131611326            for i in range(3):
    1131711327                #print i, t, f(t, i)
    11318                 assert num.allclose(f(t, i), [w, uh, 0])
     11328                assert num.allclose(f(t, i), [w, uh, 0], atol=1.0e-6)
    1131911329           
    1132011330
  • branches/numpy/anuga/shallow_water/test_shallow_water_domain.py

    r7035 r7168  
    260260        edgeflux = num.zeros(3, num.float)
    261261        zl = zr = 0.
    262         H0 = 0.0
    263 
     262        H0 = 1.0e-3 # As suggested in the manual
    264263       
    265264        max_speed = flux_function(normal, ql, qr, zl, zr, edgeflux, epsilon, g, H0)
     
    15431542            assert num.allclose(w_t, w)
    15441543            assert num.allclose(uh_t, uh)
    1545             assert num.allclose(vh_t, 0.0)
     1544            assert num.allclose(vh_t, 0.0, atol=1.0e-6)
    15461545
    15471546            # Check flows through the middle
     
    16271626            assert num.allclose(w_t, w)
    16281627            assert num.allclose(uh_t, uh)
    1629             assert num.allclose(vh_t, 0.0)
     1628            assert num.allclose(vh_t, 0.0, atol=1.0e-6)
    16301629
    16311630            # Check energies through the middle
     
    16651664        points, vertices, boundary = rectangular_cross(10, 10) # Basic mesh
    16661665        domain = Domain(points, vertices, boundary) # Create domain
    1667         domain.set_default_order(1)       
     1666        domain.set_default_order(2)
    16681667        domain.set_quantities_to_be_stored(None)
    1669         domain.set_maximum_allowed_speed(100) #FIXME (Ole): try to remove this
    1670 
    1671         # FIXME (Ole): Need tests where this is commented out
    1672         domain.tight_slope_limiters = 0    # Backwards compatibility (14/4/7)
    1673         domain.H0 = 0    # Backwards compatibility (6/2/7)
    1674         domain.use_centroid_velocities = 0    # Backwards compatibility (7/5/8)
     1668        domain.H0 = 1.0e-3
    16751669
    16761670        #-----------------------------------------------------------------
     
    17111705
    17121706        #Reference (nautilus 26/6/2008)
    1713         G0 = [-0.20000000000000001, -0.20000000000000001, -0.19920600846161715,
    1714               -0.19153647344085376, -0.19127622768281194, -0.1770671909675095,
    1715               -0.16739412133181927, -0.16196038919122191, -0.15621633053131384,
    1716               -0.15130021599977705, -0.13930978857215484, -0.19349274358263582,
    1717               -0.19975307598803765, -0.19999897143103357, -0.1999999995532111,
    1718               -0.19999999999949952, -0.19999999999949952, -0.19999999999949952,
    1719               -0.19997270012494556, -0.19925805948554556, -0.19934513778450533,
    1720               -0.19966484196394893, -0.1997352860102084,  -0.19968260481750394,
    1721               -0.19980280797303882, -0.19998804881822749, -0.19999999778075916,
    1722               -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
    1723               -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
    1724               -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
    1725               -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
    1726               -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
    1727               -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
    1728               -0.19999999999966167, -0.19999999999966167, -0.19999999999966167,
    1729               -0.19999999999966167, -0.19999999999966167, -0.19999999999966167]
    1730 
    1731         G1 = [-0.29999999999999993, -0.29999588068034899, -0.29250047332330331,
    1732               -0.28335081844518584, -0.26142206997410805, -0.22656028856329835,
    1733               -0.21224087216745585, -0.19934324109114465, -0.1889857939783175,
    1734               -0.18146311603911383, -0.17401078727434263, -0.15419361061257214,
    1735               -0.16225060576782063, -0.19010941396999181, -0.20901161407004412,
    1736               -0.21670683975774699, -0.21771386270738891, -0.21481284465869752,
    1737               -0.21063120869004387, -0.20669243364582401, -0.20320707386714859,
    1738               -0.19984087691926442, -0.19725417448019505, -0.19633783049069981,
    1739               -0.19650494599999785, -0.19708316838336942, -0.19779309449413818,
    1740               -0.19853070294429562, -0.19917342167307153, -0.19964814677795845,
    1741               -0.19991627610824922, -0.20013162970144974, -0.20029864969405509,
    1742               -0.20036259676501131, -0.20030682824965193, -0.20016105135750167,
    1743               -0.19997664501985918, -0.19980185871568762, -0.19966836175417696,
    1744               -0.19958856744312226, -0.19955954696194517, -0.19956950051110917,
    1745               -0.19960377086336181, -0.19964885299433241, -0.19969427478531132,
    1746               -0.19973301547655564, -0.19976121574277764, -0.19977765285688653,
    1747               -0.19978315117522441, -0.19977994634841739, -0.19977101394878494]
    1748 
    1749         G2 = [-0.40000000000000002, -0.39077401254732241, -0.33350466136630474,
    1750               -0.29771023004255281, -0.27605439066140897, -0.25986156218997497,
    1751               -0.24502185018573647, -0.231792624329521,   -0.21981564668803993,
    1752               -0.20870707082936543, -0.19877739883776599, -0.18980922837977957,
    1753               -0.17308011674005838, -0.16306400164013773, -0.17798470933304333,
    1754               -0.1929554075869116,  -0.20236705191987037, -0.20695767560655007,
    1755               -0.20841025876092567, -0.20792102174869989, -0.20655350005579293,
    1756               -0.20492002526259828, -0.20310627026780645, -0.20105983335287836,
    1757               -0.19937394565794653, -0.19853917506699659, -0.19836389977624452,
    1758               -0.19850305023602796, -0.19877764028836831, -0.19910928131034669,
    1759               -0.19943705712418805, -0.19970344172958865, -0.19991076989870474,
    1760               -0.20010020127747646, -0.20025937787100062, -0.20035087292905965,
    1761               -0.20035829921463297, -0.20029606557316171, -0.20019606915365515,
    1762               -0.20009096093399206, -0.20000371608204368, -0.19994495432920584,
    1763               -0.19991535665176338, -0.19990981826533513, -0.19992106419898723,
    1764               -0.19994189853516578, -0.19996624091229293, -0.19998946016985167,
    1765               -0.20000842303470234, -0.20002144460718174, -0.20002815561337187]
    1766 
    1767         G3 = [-0.45000000000000001, -0.37631169657400332, -0.33000044342859486,
    1768               -0.30586045469008522, -0.28843572253009941, -0.27215308978603808,
    1769               -0.25712951540331219, -0.2431608296216613,  -0.23032023651386374,
    1770               -0.2184546873456619,  -0.20735123704254332, -0.19740397194806389,
    1771               -0.1859829564064375,  -0.16675980728362105, -0.16951575032846536,
    1772               -0.1832860872609344,  -0.19485758939241243, -0.20231368291811427,
    1773               -0.20625610376074754, -0.20758116241495619, -0.20721445402086161,
    1774               -0.20603406830353785, -0.20450262808396991, -0.2026769581185151,
    1775               -0.2007401212066364,  -0.19931160535777592, -0.19863606301128725,
    1776               -0.19848511940572691, -0.19860091042948352, -0.19885490669377764,
    1777               -0.19916542732701112, -0.19946678238611959, -0.19971209594104697,
    1778               -0.19991912886512292, -0.2001058430788881,  -0.20024959409472989,
    1779               -0.20032160254609382, -0.20031583165752354, -0.20025051539293123,
    1780               -0.2001556115816068,  -0.20005952955420872, -0.1999814429561611,
    1781               -0.19992977821558131, -0.19990457708664208, -0.19990104785490476,
    1782               -0.19991257153954825, -0.19993258231880562, -0.19995548502882532,
    1783               -0.19997700760919687, -0.19999429663503748, -0.20000588800248761]
    1784 
    1785         #FIXME (DSG):This is a hack so the anuga install, not precompiled
    1786         # works on DSG's win2000, python 2.3
    1787         #The problem is the gauge_values[X] are 52 long, not 51.
    1788         #
    1789         # This was probably fixed by Stephen in changeset:3804
    1790         #if len(gauge_values[0]) == 52: gauge_values[0].pop()
    1791         #if len(gauge_values[1]) == 52: gauge_values[1].pop()
    1792         #if len(gauge_values[2]) == 52: gauge_values[2].pop()
    1793         #if len(gauge_values[3]) == 52: gauge_values[3].pop()
    1794 
     1707       
     1708        G0 = [-0.20000000000000001, -0.20000000000000001, -0.20000000000000001, -0.1958465301767274, -0.19059602372752493, -0.18448466250700923, -0.16979321333876071, -0.15976372740651074, -0.1575649333345176, -0.15710373731900584, -0.1530922283220747, -0.18836084336565725, -0.19921529311644628, -0.19923451799698919, -0.19923795414410964, -0.19923178806924047, -0.19925157557666154, -0.19930747801697429, -0.1993266428576112, -0.19932004930281799, -0.19929691326931867, -0.19926285267313795, -0.19916645449780995, -0.1988942593296438, -0.19900620256621993, -0.19914327423060865, -0.19918708440899577, -0.19921557252449132, -0.1992404368022069, -0.19925070370697717, -0.19925975477892274, -0.1992671090445659, -0.19927254203777162, -0.19927631910959256, -0.19927843552031504, -0.19927880339239365, -0.19927763204453783, -0.19927545249577633, -0.19927289590622824, -0.19927076261495152, -0.19926974334736983, -0.19927002562760332, -0.19927138236272213, -0.1992734501064522, -0.19927573251318192, -0.19927778936001547, -0.1992793776883893, -0.19928040577720926, -0.19928092586206753, -0.19928110982948721, -0.19928118887248453]
     1709
     1710        G1 = [-0.29999999999999993, -0.29999999999999993, -0.29139135018319512, -0.27257394456094503, -0.24141437432643265, -0.22089173942479151, -0.20796171092975532, -0.19874580192293825, -0.19014580508752857, -0.18421165368665365, -0.18020808282748838, -0.17518824759550247, -0.16436633464497749, -0.18714479115225544, -0.2045242886738807, -0.21011244240826329, -0.21151316017424124, -0.21048112933621732, -0.20772920477355789, -0.20489184334204144, -0.20286043930678221, -0.20094305756540246, -0.19948172752345467, -0.19886725178309209, -0.1986680808256765, -0.19860718133373548, -0.19862076543539733, -0.19888949069732539, -0.19932190310819023, -0.19982482967777454, -0.20036045468470615, -0.20064263130562704, -0.2007255389410077, -0.20068815669152493, -0.20057471332984647, -0.20042203940851802, -0.20026779013499779, -0.20015995671464712, -0.2000684005446525, -0.20001486753189174, -0.20000743467898013, -0.20003739771775905, -0.20008784600912621, -0.20013758305093884, -0.20017277554845025, -0.20018629233766885, -0.20018106288462198, -0.20016648079299326, -0.20015155958426531, -0.20014259747382254, -0.20014096648362509]
     1711       
     1712       
     1713        G2 = [-0.40000000000000002, -0.38885199453206343, -0.33425057028323962, -0.30154253721772117, -0.27624597383474103, -0.26037811196890087, -0.24415404585285019, -0.22941383121091052, -0.21613496492144549, -0.20418199946908885, -0.19506212965221825, -0.18851924999737435, -0.18271143344718843, -0.16910750701722474, -0.17963775224176018, -0.19442870510406052, -0.20164216917300118, -0.20467219452758528, -0.20608246104917802, -0.20640259931640445, -0.2054749739152594, -0.20380549124050265, -0.20227296931678532, -0.20095834856297176, -0.20000430919304379, -0.19946673053844086, -0.1990733499211611, -0.19882136174363013, -0.19877442300323914, -0.19905182154377868, -0.19943266521643804, -0.19988524183849191, -0.20018905307631765, -0.20031895675727809, -0.20033991149804931, -0.20031574232920274, -0.20027004750680638, -0.20020472427796293, -0.20013382447039607, -0.2000635018536408, -0.20001515436367642, -0.19998427691514989, -0.19997263083178107, -0.19998545383896535, -0.20000134502238734, -0.2000127243362736, -0.20001564474711939, -0.20001267360809977, -0.20002707579781318, -0.20004087961702843, -0.20004212947389177]
     1714       
     1715        G3 = [-0.45000000000000001, -0.38058172993544187, -0.33756059941741273, -0.31015371357441396, -0.29214769368562965, -0.27545447937118606, -0.25871585649808154, -0.24254276680581988, -0.22758633129006092, -0.21417276895743134, -0.20237184768790789, -0.19369491041576814, -0.18721625333717057, -0.1794243868465818, -0.17052113574042196, -0.18534300640363346, -0.19601184621026671, -0.20185028431829469, -0.20476187496918136, -0.20602933256960082, -0.20598569228739247, -0.20492643756666848, -0.20339695828762758, -0.20196440373022595, -0.20070304052919338, -0.19986227854052355, -0.19933020476408528, -0.19898034831018035, -0.19878317651286193, -0.19886879323961787, -0.19915860801206181, -0.19953675278099042, -0.19992828019602107, -0.20015957043092364, -0.20025268671087426, -0.20028559516444974, -0.20027084877341045, -0.20022991487243985, -0.20016234295579871, -0.20009131445092507, -0.20003149397006523, -0.19998473356473795, -0.19996011913447218, -0.19995647168667186, -0.19996526209120422, -0.19996600297827097, -0.19997268800221216, -0.19998682275066659, -0.20000372259781876, -0.20001628681983963, -0.2000173314086407]
     1716       
    17951717        assert num.allclose(gauge_values[0], G0)
    1796         # FIXME(Ole): Disabled when ticket:314 was resolved.
    1797         # The slight change might in fact be for the better.
    1798         #assert num.allclose(gauge_values[1], G1)
     1718        assert num.allclose(gauge_values[1], G1)
    17991719        assert num.allclose(gauge_values[2], G2)
    18001720        assert num.allclose(gauge_values[3], G3)
     
    40553975        domain.beta_vh = 0.9
    40563976        domain.beta_vh_dry = 0.9
    4057         domain.H0 = 0    # Backwards compatibility (6/2/7)
     3977        domain.H0 = 1.0e-3
    40583978
    40593979        # Boundary conditions
     
    40723992        assert num.allclose(domain.max_timestep, 0.0396825396825)
    40733993
     3994        msg = ("domain.quantities['stage'].centroid_values[:12]=%s"
     3995               % str(domain.quantities['stage'].centroid_values[:12]))
    40743996        assert num.allclose(domain.quantities['stage'].centroid_values[:12],
    40753997                            [0.00171396, 0.02656103, 0.00241523, 0.02656103,
    40763998                             0.00241523, 0.02656103, 0.00241523, 0.02656103,
    4077                              0.00241523, 0.02656103, 0.00241523, 0.0272623])
     3999                             0.00241523, 0.02656103, 0.00241523, 0.0272623]), msg
    40784000
    40794001        domain.distribute_to_vertices_and_edges()
    40804002
    40814003        assert num.allclose(domain.quantities['stage'].vertex_values[:12,0],
    4082                             [0.0001714, 0.02656103, 0.00024152,
    4083                              0.02656103, 0.00024152, 0.02656103,
    4084                              0.00024152, 0.02656103, 0.00024152,
    4085                              0.02656103, 0.00024152, 0.0272623])
    4086 
     4004                            [0.001, 0.02656103, 0.001, 0.02656103, 0.001, 0.02656103,
     4005                             0.001, 0.02656103, 0.001, 0.02656103, 0.001, 0.0272623])     
     4006                             
    40874007        assert num.allclose(domain.quantities['stage'].vertex_values[:12,1],
    4088                             [0.00315012, 0.02656103, 0.00024152, 0.02656103,
    4089                              0.00024152, 0.02656103, 0.00024152, 0.02656103,
    4090                              0.00024152, 0.02656103, 0.00040506, 0.0272623])
     4008                            [0.00237867, 0.02656103, 0.001, 0.02656103, 0.001,
     4009                             0.02656103, 0.001, 0.02656103, 0.001, 0.02656103,
     4010                             0.00110647, 0.0272623])
    40914011
    40924012        assert num.allclose(domain.quantities['stage'].vertex_values[:12,2],
    4093                             [0.00182037, 0.02656103, 0.00676264,
    4094                              0.02656103, 0.00676264, 0.02656103,
    4095                              0.00676264, 0.02656103, 0.00676264,
    4096                              0.02656103, 0.0065991, 0.0272623])
     4013                            [0.00176321, 0.02656103, 0.00524568,
     4014                             0.02656103, 0.00524568, 0.02656103,
     4015                             0.00524568, 0.02656103, 0.00524568,
     4016                             0.02656103, 0.00513921, 0.0272623])
    40974017
    40984018        assert num.allclose(domain.quantities['xmomentum'].centroid_values[:12],
     
    41524072        domain.smooth = False
    41534073        domain.default_order = 1
    4154         domain.H0 = 0    # Backwards compatibility (6/2/7)
     4074        domain.H0 = 1.0e-3 # As suggested in the manual
    41554075
    41564076        # Boundary conditions
     
    41974117        domain.beta_vh = 0.9
    41984118        domain.beta_vh_dry = 0.9
    4199         #Makes it like the 'oldstyle' balance
    4200         #domain.minimum_allowed_height = 0.0
    4201         domain.H0 = 0    # Backwards compatibility (6/2/7)
     4119
     4120        domain.H0 = 1.0e-3 # As suggested in the manual
    42024121        domain.use_centroid_velocities = False # Backwards compatibility (8/5/8)
    42034122        domain.set_maximum_allowed_speed(1.0)
     
    42244143        #                [0.00101913,0.05352143,0.00104852,0.05354394])
    42254144
    4226         #FIXME: These numbers were from version before 21/3/6 -
    4227         #could be recreated by setting maximum_allowed_speed to 0 maybe
    4228         #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
    4229         #                [ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
    4230 
     4145        # Slight change due to flux limiter optimisation 28/5/9
     4146        assert num.allclose(domain.quantities['stage'].vertex_values[:4,0],
     4147                            [0.001, 0.05350407, 0.00106768, 0.05352525])
    42314148        assert num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
    4232                             [ 0.00090581, 0.03685719, 0.00088303, 0.03687313])
    4233 
    4234         #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
    4235         #                [0.00090581,0.03685719,0.00088303,0.03687313])
     4149                            [0.0008628, 0.03684647, 0.00087764, 0.03686007])
    42364150
    42374151        assert num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
    4238                             [-0.00139463,0.0006156,-0.00060364,0.00061827])
     4152                            [-0.00142114, 0.00061557, -0.00062362, 0.00061896])
    42394153
    42404154        os.remove(domain.get_name() + '.sww')
     4155
    42414156
    42424157    def test_flatbed_second_order_vmax_0(self):
     
    42584173        domain.beta_vh_dry = 0.9
    42594174        domain.maximum_allowed_speed = 0.0    # Makes it like the 'oldstyle'
    4260         domain.H0 = 0    # Backwards compatibility (6/2/7)
     4175        domain.H0 = 1.0e-3 # As suggested in the manual
    42614176        domain.use_centroid_velocities = False # Backwards compatibility (8/5/8)
    42624177
     
    42784193        #could be recreated by setting maximum_allowed_speed to 0 maybe
    42794194        assert num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
    4280                             [ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
    4281 
     4195                            [0.00066963, 0.03684647, 0.00085288, 0.03686007])
     4196       
    42824197
    42834198        assert num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
    4284                             [-0.00139463,0.0006156,-0.00060364,0.00061827])
     4199                            [-0.00142114, 0.00061557, -0.00062362, 0.00061896])
    42854200
    42864201        os.remove(domain.get_name() + '.sww')
     
    43064221        domain.beta_vh = 0.9
    43074222        domain.beta_vh_dry = 0.9
    4308         domain.H0 = 0    # Backwards compatibility (6/2/7)
    4309         domain.use_centroid_velocities = False # Backwards compatibility (8/5/8)
    4310         domain.set_maximum_allowed_speed(1.0)
     4223
     4224        domain.H0 = 1.0e-3 # ANUGA manual 28/5/9
    43114225
    43124226        # Boundary conditions
     
    43714285                assert num.allclose(domain.max_timestep, 0.0210448446782)
    43724286
    4373             # Centroids were correct but not vertices.
    4374             # Hence the check of distribute below.
    4375             assert num.allclose(domain.quantities['stage'].centroid_values[:4],
    4376                                 [0.00721206, 0.05352143,
    4377                                  0.01009108, 0.05354394])
    4378 
    4379             assert num.allclose(domain.quantities['xmomentum'].\
    4380                                     centroid_values[:4],
    4381                                 [0.00648352, 0.03685719,
    4382                                  0.00850733, 0.03687313])
    4383 
    4384             assert num.allclose(domain.quantities['ymomentum'].\
    4385                                     centroid_values[:4],
    4386                                 [-0.00139463, 0.0006156,
    4387                                  -0.00060364, 0.00061827])
    4388 
    4389             #assert allclose(domain.quantities['xmomentum'].\
    4390             #                    centroid_values[17],0.00028606084)
     4287            #print domain.quantities['stage'].centroid_values[:4]
     4288            #print domain.quantities['xmomentum'].centroid_values[:4]
     4289            #print domain.quantities['ymomentum'].centroid_values[:4]                       
     4290               
     4291            #Centroids were correct but not vertices.
     4292            #Hence the check of distribute below.
     4293
    43914294            if not V:
    4392                 #FIXME: These numbers were from version before 21/3/6 -
    4393                 #could be recreated by setting maximum_allowed_speed to 0 maybe
    4394                 #assert allclose(domain.quantities['xmomentum'].\
    4395                 #                    centroid_values[17], 0.0)
    4396                 assert num.allclose(domain.quantities['xmomentum'].\
    4397                                         centroid_values[17],
    4398                                     0.000286060839592)
    4399 
     4295                assert num.allclose(domain.quantities['stage'].centroid_values[:4],
     4296                                    [0.00725574, 0.05350737, 0.01008413, 0.0535293])           
     4297                assert num.allclose(domain.quantities['xmomentum'].centroid_values[:4],
     4298                                    [0.00654964, 0.03684904, 0.00852561, 0.03686323])
     4299
     4300                assert num.allclose(domain.quantities['ymomentum'].centroid_values[:4],
     4301                                    [-0.00143169, 0.00061027, -0.00062325, 0.00061408])
     4302
     4303                                   
     4304           
     4305                assert num.allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)               
    44004306            else:
    44014307                assert num.allclose(domain.quantities['xmomentum'].\
    44024308                                        centroid_values[17],
    44034309                                    0.00028606084)
     4310                return #FIXME - Bailout for V True
    44044311
    44054312            import copy
     
    44114318            domain.distribute_to_vertices_and_edges()
    44124319
    4413             #assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
    4414             #assert allclose(domain.quantities['xmomentum'].centroid_values[17],
    4415             #                0.0)
    4416 
    4417             assert num.allclose(domain.quantities['xmomentum'].\
    4418                                     centroid_values[17],
    4419                                 0.000286060839592)
    4420 
    4421             #FIXME: These numbers were from version before 25/10
    4422             #assert allclose(domain.quantities['stage'].vertex_values[:4,0],
    4423             #                [0.00101913,0.05352143,0.00104852,0.05354394])
    4424 
    4425             assert num.allclose(domain.quantities['ymomentum'].\
    4426                                     vertex_values[:4,0],
    4427                                 [-0.00139463, 0.0006156,
    4428                                  -0.00060364, 0.00061827])
    4429 
    4430             assert num.allclose(domain.quantities['xmomentum'].\
    4431                                     vertex_values[:4,0],
    4432                                 [0.00090581, 0.03685719,
    4433                                  0.00088303, 0.03687313])
     4320            assert num.allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)
     4321
     4322
     4323
     4324            assert num.allclose(domain.quantities['ymomentum'].vertex_values[:4,0], [-0.00019732, 0.00061027, -0.00062325, 0.00061408])
     4325
     4326            assert num.allclose(domain.quantities['xmomentum'].vertex_values[:4,0], [0.00090268, 0.03684904, 0.00085256, 0.03686323])
     4327
     4328
     4329
    44344330
    44354331        os.remove(domain.get_name() + '.sww')
     
    46024498
    46034499        assert num.allclose(domain.quantities['stage'].centroid_values,
    4604                             [ 0.01290985,  0.02356019,  0.01619096,
    4605                               0.02356019,  0.01619096,  0.02356019,
    4606                               0.01619096,  0.02356019,  0.01619096,
    4607                               0.02356019,  0.01619096,  0.0268413,
    4608                              -0.04411074, -0.0248011,  -0.04186556,
    4609                              -0.0248011,  -0.04186556, -0.0248011,
    4610                              -0.04186556, -0.0248011,  -0.04186556,
    4611                              -0.0248011,  -0.04186556, -0.02255593,
    4612                              -0.09966629, -0.08035666, -0.09742112,
    4613                              -0.08035666, -0.09742112, -0.08035666,
    4614                              -0.09742112, -0.08035666, -0.09742112,
    4615                              -0.08035666, -0.09742112, -0.07811149,
    4616                              -0.15522185, -0.13591222, -0.15297667,
    4617                              -0.13591222, -0.15297667, -0.13591222,
    4618                              -0.15297667, -0.13591222, -0.15297667,
    4619                              -0.13591222, -0.15297667, -0.13366704,
    4620                              -0.2107774,  -0.19146777, -0.20853223,
    4621                              -0.19146777, -0.20853223, -0.19146777,
    4622                              -0.20853223, -0.19146777, -0.20853223,
    4623                              -0.19146777, -0.20853223, -0.1892226,
    4624                              -0.26120669, -0.24776246, -0.25865535,
    4625                              -0.24776246, -0.25865535, -0.24776246,
    4626                              -0.25865535, -0.24776246, -0.25865535,
    4627                              -0.24776246, -0.25865535, -0.24521113])
     4500                            [ 0.01290985,  0.02356019,  0.01619096,  0.02356019,  0.01619096,
     4501                              0.02356019,  0.01619096,  0.02356019,  0.01619096,  0.02356019,
     4502                              0.01619096,  0.0268413,  -0.04411074, -0.0248011,  -0.04186556,
     4503                             -0.0248011,  -0.04186556, -0.0248011,  -0.04186556, -0.0248011,
     4504                             -0.04186556, -0.0248011,  -0.04186556, -0.02255593,
     4505                             -0.09966629, -0.08035666, -0.09742112, -0.08035666,
     4506                             -0.09742112, -0.08035666, -0.09742112, -0.08035666,
     4507                             -0.09742112, -0.08035666, -0.09742112, -0.07811149,
     4508                             -0.15522185, -0.13591222, -0.15297667, -0.13591222,
     4509                             -0.15297667, -0.13591222, -0.15297667, -0.13591222,
     4510                             -0.15297667, -0.13591222, -0.15297667, -0.13366704,
     4511                             -0.2107774,  -0.19146777, -0.20853223, -0.19146777,
     4512                             -0.20853223, -0.19146777, -0.20853223, -0.19146777,
     4513                             -0.20853223, -0.19146777, -0.20853223, -0.1892226,
     4514                             -0.26120669, -0.24776246, -0.25865535, -0.24776246,
     4515                             -0.25865535, -0.24776246, -0.25865535, -0.24776246,
     4516                             -0.25865535, -0.24776246, -0.25865535, -0.24521113])
    46284517
    46294518        os.remove(domain.get_name() + '.sww')
     
    47004589
    47014590        assert num.allclose(domain.quantities['stage'].centroid_values,
    4702                             [ 0.00855788,  0.01575204,  0.00994606,
    4703                               0.01717072,  0.01005985,  0.01716362,
    4704                               0.01005985,  0.01716299,  0.01007098,
    4705                               0.01736248,  0.01216452,  0.02026776,
    4706                              -0.04462374, -0.02479045, -0.04199789,
    4707                              -0.0229465,  -0.04184033, -0.02295693,
    4708                              -0.04184013, -0.02295675, -0.04184486,
    4709                              -0.0228168,  -0.04028876, -0.02036486,
    4710                              -0.10029444, -0.08170809, -0.09772846,
    4711                              -0.08021704, -0.09760006, -0.08022143,
    4712                              -0.09759984, -0.08022124, -0.09760261,
    4713                              -0.08008893, -0.09603914, -0.07758209,
    4714                              -0.15584152, -0.13723138, -0.15327266,
    4715                              -0.13572906, -0.15314427, -0.13573349,
    4716                              -0.15314405, -0.13573331, -0.15314679,
    4717                              -0.13560104, -0.15158523, -0.13310701,
    4718                              -0.21208605, -0.19283913, -0.20955631,
    4719                              -0.19134189, -0.20942821, -0.19134598,
    4720                              -0.20942799, -0.1913458,  -0.20943005,
    4721                              -0.19120952, -0.20781177, -0.18869401,
    4722                              -0.25384082, -0.2463294,  -0.25047649,
    4723                              -0.24464654, -0.25031159, -0.24464253,
    4724                              -0.25031112, -0.24464253, -0.25031463,
    4725                              -0.24454764, -0.24885323, -0.24286438])
     4591                            [ 0.00855788,  0.01575204,  0.00994606,  0.01717072,
     4592                              0.01005985,  0.01716362,  0.01005985,  0.01716299,
     4593                              0.01007098,  0.01736248,  0.01216452,  0.02026776,
     4594                             -0.04462374, -0.02479045, -0.04199789, -0.0229465,
     4595                             -0.04184033, -0.02295693, -0.04184013, -0.02295675,
     4596                             -0.04184486, -0.0228168,  -0.04028876, -0.02036486,
     4597                             -0.10029444, -0.08170809, -0.09772846, -0.08021704,
     4598                             -0.09760006, -0.08022143, -0.09759984, -0.08022124,
     4599                             -0.09760261, -0.08008893, -0.09603914, -0.07758209,
     4600                             -0.15584152, -0.13723138, -0.15327266, -0.13572906,
     4601                             -0.15314427, -0.13573349, -0.15314405, -0.13573331,
     4602                             -0.15314679, -0.13560104, -0.15158523, -0.13310701,
     4603                             -0.21208605, -0.19283913, -0.20955631, -0.19134189,
     4604                             -0.20942821, -0.19134598, -0.20942799, -0.1913458,
     4605                             -0.20943005, -0.19120952, -0.20781177, -0.18869401,
     4606                             -0.25384082, -0.2463294,  -0.25047649, -0.24464654,
     4607                             -0.25031159, -0.24464253, -0.25031112, -0.24464253,
     4608                             -0.25031463, -0.24454764, -0.24885323, -0.24286438])
     4609
    47264610
    47274611        os.remove(domain.get_name() + '.sww')
     
    47964680
    47974681        assert num.allclose(domain.quantities['stage'].centroid_values,
    4798                             [ 0.00855788,  0.01575204,  0.00994606,
    4799                               0.01717072,  0.01005985,  0.01716362,
    4800                               0.01005985,  0.01716299,  0.01007098,
    4801                               0.01736248,  0.01216452,  0.02026776,
    4802                              -0.04462374, -0.02479045, -0.04199789,
    4803                              -0.0229465,  -0.04184033, -0.02295693,
    4804                              -0.04184013, -0.02295675, -0.04184486,
    4805                              -0.0228168,  -0.04028876, -0.02036486,
    4806                              -0.10029444, -0.08170809, -0.09772846,
    4807                              -0.08021704, -0.09760006, -0.08022143,
    4808                              -0.09759984, -0.08022124, -0.09760261,
    4809                              -0.08008893, -0.09603914, -0.07758209,
    4810                              -0.15584152, -0.13723138, -0.15327266,
    4811                              -0.13572906, -0.15314427, -0.13573349,
    4812                              -0.15314405, -0.13573331, -0.15314679,
    4813                              -0.13560104, -0.15158523, -0.13310701,
    4814                              -0.21208605, -0.19283913, -0.20955631,
    4815                              -0.19134189, -0.20942821, -0.19134598,
    4816                              -0.20942799, -0.1913458,  -0.20943005,
    4817                              -0.19120952, -0.20781177, -0.18869401,
    4818                              -0.25384082, -0.2463294,  -0.25047649,
    4819                              -0.24464654, -0.25031159, -0.24464253,
    4820                              -0.25031112, -0.24464253, -0.25031463,
    4821                              -0.24454764, -0.24885323, -0.24286438])
     4682                            [ 0.00855788,  0.01575204,  0.00994606,  0.01717072,  0.01005985,
     4683                              0.01716362,  0.01005985,  0.01716299,  0.01007098,  0.01736248,
     4684                              0.01216452,  0.02026776, -0.04462374, -0.02479045, -0.04199789,
     4685                             -0.0229465,  -0.04184033, -0.02295693, -0.04184013,
     4686                             -0.02295675, -0.04184486, -0.0228168,  -0.04028876,
     4687                             -0.02036486, -0.10029444, -0.08170809, -0.09772846,
     4688                             -0.08021704, -0.09760006, -0.08022143, -0.09759984,
     4689                             -0.08022124, -0.09760261, -0.08008893, -0.09603914,
     4690                             -0.07758209, -0.15584152, -0.13723138, -0.15327266,
     4691                             -0.13572906, -0.15314427, -0.13573349, -0.15314405,
     4692                             -0.13573331, -0.15314679, -0.13560104, -0.15158523,
     4693                             -0.13310701, -0.21208605, -0.19283913, -0.20955631,
     4694                             -0.19134189, -0.20942821, -0.19134598, -0.20942799,
     4695                             -0.1913458,  -0.20943005, -0.19120952, -0.20781177,
     4696                             -0.18869401, -0.25384082, -0.2463294,  -0.25047649,
     4697                             -0.24464654, -0.25031159, -0.24464253, -0.25031112,
     4698                             -0.24464253, -0.25031463, -0.24454764, -0.24885323,
     4699                             -0.24286438])
    48224700
    48234701        os.remove(domain.get_name() + '.sww')
     
    48924770
    48934771        assert num.allclose(domain.quantities['stage'].centroid_values,
    4894                             [-0.02907028, -0.01475478, -0.02973417,
    4895                              -0.01447186, -0.02932665, -0.01428285,
    4896                              -0.02901975, -0.0141361,  -0.02898816,
    4897                              -0.01418135, -0.02961409, -0.01403487,
    4898                              -0.07597998, -0.06252591, -0.07664854,
    4899                              -0.06312532, -0.07638287, -0.06265139,
    4900                              -0.07571145, -0.06235231, -0.0756817,
    4901                              -0.06245309, -0.07652292, -0.06289946,
    4902                              -0.12367464, -0.11088981, -0.12237277,
    4903                              -0.11115338, -0.1218934,  -0.1107174,
    4904                              -0.12081485, -0.11000491, -0.12038451,
    4905                              -0.11010335, -0.12102113, -0.11012105,
    4906                              -0.16909116, -0.15831543, -0.16730214,
    4907                              -0.15786249, -0.1665493,  -0.15697919,
    4908                              -0.16496618, -0.15559852, -0.16338679,
    4909                              -0.15509088, -0.16364092, -0.15424423,
    4910                              -0.18771107, -0.19903904, -0.18903759,
    4911                              -0.19858437, -0.18701552, -0.19697797,
    4912                              -0.1833593,  -0.19505871, -0.1818806,
    4913                              -0.19418042, -0.18586159, -0.19576946,
    4914                              -0.13986873, -0.14170053, -0.14132188,
    4915                              -0.14560674, -0.14095617, -0.14373292,
    4916                              -0.13785933, -0.14033364, -0.13592955,
    4917                              -0.13936356, -0.13596008, -0.14216296])
     4772     [-0.02907028, -0.01475478, -0.02973417, -0.01447186, -0.02932665, -0.01428285,
     4773      -0.02901975, -0.0141361,  -0.02898816, -0.01418135, -0.02961409, -0.01403487,
     4774      -0.07597998, -0.06252591, -0.07664854, -0.06312532, -0.07638287, -0.06265139,
     4775      -0.07571145, -0.06235231, -0.0756817,  -0.06245309, -0.07652292, -0.06289946,
     4776      -0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934,  -0.1107174,
     4777      -0.12081485, -0.11000491, -0.12038451, -0.11010335, -0.12102113, -0.11012105,
     4778      -0.16909116, -0.15831543, -0.16730214, -0.15786249, -0.1665493,  -0.15697919,
     4779      -0.16496618, -0.15559852, -0.16338679, -0.15509088, -0.16364092, -0.15424423,
     4780      -0.18771107, -0.19903904, -0.18903759, -0.19858437, -0.18701552, -0.19697797,
     4781      -0.1833593,  -0.19505871, -0.1818806,  -0.19418042, -0.18586159, -0.19576946,
     4782      -0.13986873, -0.14170053, -0.14132188, -0.14560674, -0.14095617, -0.14373292,
     4783      -0.13785933, -0.14033364, -0.13592955, -0.13936356, -0.13596008, -0.14216296])
    49184784
    49194785        assert num.allclose(domain.quantities['xmomentum'].centroid_values,
    4920                             [ 0.00831121,  0.00317948,  0.00731797,
    4921                               0.00334939,  0.00764717,  0.00348053,
    4922                               0.00788729,  0.00356522,  0.00780649,
    4923                               0.00341919,  0.00693525,  0.00310375,
    4924                               0.02166196,  0.01421475,  0.02017737,
    4925                               0.01316839,  0.02037015,  0.01368659,
    4926                               0.02106,     0.01399161,  0.02074514,
    4927                               0.01354935,  0.01887407,  0.0123113,
    4928                               0.03775083,  0.02855197,  0.03689337,
    4929                               0.02759782,  0.03732848,  0.02812072,
    4930                               0.03872545,  0.02913348,  0.03880939,
    4931                               0.02803804,  0.03546499,  0.0260039,
    4932                               0.0632131,   0.04730634,  0.0576324,
    4933                               0.04592336,  0.05790921,  0.04690514,
    4934                               0.05986467,  0.04871165,  0.06170068,
    4935                               0.04811572,  0.05657041,  0.04416292,
    4936                               0.08489642,  0.07188097,  0.07835261,
    4937                               0.06843406,  0.07986412,  0.0698247,
    4938                               0.08201071,  0.07216756,  0.08378418,
    4939                               0.07273624,  0.080399,    0.06645841,
    4940                               0.01631548,  0.04691608,  0.0206632,
    4941                               0.044409,    0.02115518,  0.04560305,
    4942                               0.02160608,  0.04663725,  0.02174734,
    4943                               0.04795559,  0.02281427,  0.05667111])
     4786     [ 0.00831121,  0.00317948,  0.00731797,  0.00334939,  0.00764717,  0.00348053,
     4787       0.00788729,  0.00356522,  0.00780649,  0.00341919,  0.00693525,  0.00310375,
     4788       0.02166196,  0.01421475,  0.02017737,  0.01316839,  0.02037015,  0.01368659,
     4789       0.02106,     0.01399161,  0.02074514,  0.01354935,  0.01887407,  0.0123113,
     4790       0.03775083,  0.02855197,  0.03689337,  0.02759782,  0.03732848,  0.02812072,
     4791       0.03872545,  0.02913348,  0.03880939,  0.02803804,  0.03546499,  0.0260039,
     4792       0.0632131,   0.04730634,  0.0576324,   0.04592336,  0.05790921,  0.04690514,
     4793       0.05986467,  0.04871165,  0.06170068,  0.04811572,  0.05657041,  0.04416292,
     4794       0.08489642,  0.07188097,  0.07835261,  0.06843406,  0.07986412,  0.0698247,
     4795       0.08201071,  0.07216756,  0.08378418,  0.07273624,  0.080399,    0.06645841,
     4796       0.01631548,  0.04691608,  0.0206632,   0.044409,    0.02115518,  0.04560305,
     4797       0.02160608,  0.04663725,  0.02174734,  0.04795559,  0.02281427,  0.05667111])
    49444798
    49454799
    49464800        assert num.allclose(domain.quantities['ymomentum'].centroid_values,
    4947                             [ 1.45876601e-004, -3.24627393e-004,
    4948                              -1.57572719e-004, -2.92790187e-004,
    4949                              -9.90988382e-005, -3.06677335e-004,
    4950                              -1.62493106e-004, -3.71310004e-004,
    4951                              -1.99445058e-004, -3.28493467e-004,
    4952                               6.68217349e-005, -8.42042805e-006,
    4953                               5.05093371e-004, -1.42842214e-004,
    4954                              -6.81454718e-005, -5.02084057e-004,
    4955                              -8.50583861e-005, -4.65443981e-004,
    4956                              -1.96406564e-004, -5.88889562e-004,
    4957                              -2.70160173e-004, -5.35485454e-004,
    4958                               2.60780997e-004,  3.12145471e-005,
    4959                               5.16189608e-004,  1.07069062e-004,
    4960                               9.29989252e-005, -3.71211119e-004,
    4961                               1.16350246e-004, -3.82407830e-004,
    4962                              -1.62077969e-004, -6.30906636e-004,
    4963                              -4.74025708e-004, -6.94463009e-004,
    4964                               6.15092843e-005,  2.22106820e-004,
    4965                              -6.29589294e-004,  2.43611937e-004,
    4966                              -5.88125094e-004, -6.94293192e-005,
    4967                              -4.17914641e-004,  6.64609019e-005,
    4968                              -7.68334577e-004, -3.40232101e-004,
    4969                              -1.67424308e-003, -7.39485066e-004,
    4970                              -1.59966988e-003,  5.68262838e-005,
    4971                              -1.48470633e-003, -1.84554882e-003,
    4972                              -2.27200099e-003, -1.67506848e-003,
    4973                              -1.95610258e-003, -1.47638801e-003,
    4974                              -1.73779477e-003, -1.85498791e-003,
    4975                              -2.01357843e-003, -2.17675471e-003,
    4976                              -1.65783870e-003, -1.15818681e-003,
    4977                              -1.18663036e-003, -2.94229849e-003,
    4978                              -3.59309018e-003, -5.13496584e-003,
    4979                              -6.17359400e-003, -5.98761937e-003,
    4980                              -6.00540116e-003, -5.01121966e-003,
    4981                              -4.50964850e-003, -3.06319963e-003,
    4982                               6.08950810e-004, -4.79537921e-004])
     4801     [ 1.45876601e-004, -3.24627393e-004, -1.57572719e-004, -2.92790187e-004,
     4802      -9.90988382e-005, -3.06677335e-004, -1.62493106e-004, -3.71310004e-004,
     4803      -1.99445058e-004, -3.28493467e-004,  6.68217349e-005, -8.42042805e-006,
     4804       5.05093371e-004, -1.42842214e-004, -6.81454718e-005, -5.02084057e-004,
     4805      -8.50583861e-005, -4.65443981e-004, -1.96406564e-004, -5.88889562e-004,
     4806      -2.70160173e-004, -5.35485454e-004,  2.60780997e-004,  3.12145471e-005,
     4807       5.16189608e-004,  1.07069062e-004,  9.29989252e-005, -3.71211119e-004,
     4808       1.16350246e-004, -3.82407830e-004, -1.62077969e-004, -6.30906636e-004,
     4809      -4.74025708e-004, -6.94463009e-004,  6.15092843e-005,  2.22106820e-004,
     4810      -6.29589294e-004,  2.43611937e-004, -5.88125094e-004, -6.94293192e-005,
     4811      -4.17914641e-004,  6.64609019e-005, -7.68334577e-004, -3.40232101e-004,
     4812      -1.67424308e-003, -7.39485066e-004, -1.59966988e-003,  5.68262838e-005,
     4813      -1.48470633e-003, -1.84554882e-003, -2.27200099e-003, -1.67506848e-003,
     4814      -1.95610258e-003, -1.47638801e-003, -1.73779477e-003, -1.85498791e-003,
     4815      -2.01357843e-003, -2.17675471e-003, -1.65783870e-003, -1.15818681e-003,
     4816      -1.18663036e-003, -2.94229849e-003, -3.59309018e-003, -5.13496584e-003,
     4817      -6.17359400e-003, -5.98761937e-003, -6.00540116e-003, -5.01121966e-003,
     4818      -4.50964850e-003, -3.06319963e-003,  6.08950810e-004, -4.79537921e-004])
    49834819
    49844820        os.remove(domain.get_name() + '.sww')
     
    50544890
    50554891        assert num.allclose(domain.quantities['stage'].centroid_values,
    5056                             [-0.03348416, -0.01749303, -0.03299091,
    5057                              -0.01739241, -0.03246447, -0.01732016,
    5058                              -0.03205390, -0.01717833, -0.03146383,
    5059                              -0.01699831, -0.03076577, -0.01671795,
    5060                              -0.07952656, -0.06684763, -0.07721455,
    5061                              -0.06668388, -0.07632976, -0.06600113,
    5062                              -0.07523678, -0.06546373, -0.07447040,
    5063                              -0.06508861, -0.07438723, -0.06359288,
    5064                              -0.12526729, -0.11205668, -0.12179433,
    5065                              -0.11068104, -0.12048395, -0.10968948,
    5066                              -0.11912023, -0.10862628, -0.11784090,
    5067                              -0.10803744, -0.11790629, -0.10742354,
    5068                              -0.16859613, -0.15427413, -0.16664444,
    5069                              -0.15464452, -0.16570816, -0.15327556,
    5070                              -0.16409162, -0.15204092, -0.16264608,
    5071                              -0.15102139, -0.16162736, -0.14969205,
    5072                              -0.18736511, -0.19874036, -0.18811230,
    5073                              -0.19758289, -0.18590182, -0.19580301,
    5074                              -0.18234588, -0.19423215, -0.18100376,
    5075                              -0.19380116, -0.18509710, -0.19501636,
    5076                              -0.13982382, -0.14166819, -0.14132775,
    5077                              -0.14528694, -0.14096905, -0.14351126,
    5078                              -0.13800356, -0.14027920, -0.13613538,
    5079                              -0.13936795, -0.13621902, -0.14204982])
    5080 
     4892         [-0.03348416, -0.01749303, -0.03299091, -0.01739241, -0.03246447, -0.01732016,
     4893          -0.03205390, -0.01717833, -0.03146383, -0.01699831, -0.03076577, -0.01671795,
     4894          -0.07952656, -0.06684763, -0.07721455, -0.06668388, -0.07632976, -0.06600113,
     4895          -0.07523678, -0.06546373, -0.07447040, -0.06508861, -0.07438723, -0.06359288,
     4896          -0.12526729, -0.11205668, -0.12179433, -0.11068104, -0.12048395, -0.10968948,
     4897          -0.11912023, -0.10862628, -0.11784090, -0.10803744, -0.11790629, -0.10742354,
     4898          -0.16859613, -0.15427413, -0.16664444, -0.15464452, -0.16570816, -0.15327556,
     4899          -0.16409162, -0.15204092, -0.16264608, -0.15102139, -0.16162736, -0.14969205,
     4900          -0.18736511, -0.19874036, -0.18811230, -0.19758289, -0.18590182, -0.19580301,
     4901          -0.18234588, -0.19423215, -0.18100376, -0.19380116, -0.18509710, -0.19501636,
     4902          -0.13982382, -0.14166819, -0.14132775, -0.14528694, -0.14096905, -0.14351126,
     4903          -0.13800356, -0.14027920, -0.13613538, -0.13936795, -0.13621902, -0.14204982])
     4904
     4905                     
    50814906        assert num.allclose(domain.quantities['xmomentum'].centroid_values,
    5082                             [0.00600290,  0.00175780,  0.00591905,
    5083                              0.00190903,  0.00644462,  0.00203095,
    5084                              0.00684561,  0.00225089,  0.00708208,
    5085                              0.00236235,  0.00649095,  0.00222343,
    5086                              0.02068693,  0.01164034,  0.01983343,
    5087                              0.01159526,  0.02044611,  0.01233252,
    5088                              0.02135685,  0.01301289,  0.02161290,
    5089                              0.01260280,  0.01867612,  0.01133078,
    5090                              0.04091313,  0.02668283,  0.03634781,
    5091                              0.02733469,  0.03767692,  0.02836840,
    5092                              0.03906338,  0.02958073,  0.04025669,
    5093                              0.02953292,  0.03665616,  0.02583565,
    5094                              0.06314558,  0.04830935,  0.05663609,
    5095                              0.04564362,  0.05756200,  0.04739673,
    5096                              0.05967379,  0.04919083,  0.06124330,
    5097                              0.04965808,  0.05879240,  0.04629319,
    5098                              0.08220739,  0.06924725,  0.07713556,
    5099                              0.06782640,  0.07909499,  0.06992544,
    5100                              0.08116621,  0.07210181,  0.08281548,
    5101                              0.07222669,  0.07941059,  0.06755612,
    5102                              0.01581588,  0.04533609,  0.02017939,
    5103                              0.04342565,  0.02073232,  0.04476108,
    5104                              0.02117439,  0.04573358,  0.02129473,
    5105                              0.04694267,  0.02220398,  0.05533458])
     4907      [0.00600290,  0.00175780,  0.00591905,  0.00190903,  0.00644462,  0.00203095,
     4908       0.00684561,  0.00225089,  0.00708208,  0.00236235,  0.00649095,  0.00222343,
     4909       0.02068693,  0.01164034,  0.01983343,  0.01159526,  0.02044611,  0.01233252,
     4910       0.02135685,  0.01301289,  0.02161290,  0.01260280,  0.01867612,  0.01133078,
     4911       0.04091313,  0.02668283,  0.03634781,  0.02733469,  0.03767692,  0.02836840,
     4912       0.03906338,  0.02958073,  0.04025669,  0.02953292,  0.03665616,  0.02583565,
     4913       0.06314558,  0.04830935,  0.05663609,  0.04564362,  0.05756200,  0.04739673,
     4914       0.05967379,  0.04919083,  0.06124330,  0.04965808,  0.05879240,  0.04629319,
     4915       0.08220739,  0.06924725,  0.07713556,  0.06782640,  0.07909499,  0.06992544,
     4916       0.08116621,  0.07210181,  0.08281548,  0.07222669,  0.07941059,  0.06755612,
     4917       0.01581588,  0.04533609,  0.02017939,  0.04342565,  0.02073232,  0.04476108,
     4918       0.02117439,  0.04573358,  0.02129473,  0.04694267,  0.02220398,  0.05533458])
    51064919
    51074920
    51084921        assert num.allclose(domain.quantities['ymomentum'].centroid_values,
    5109                             [-7.65882069e-005, -1.46087080e-004,
    5110                              -1.09630102e-004, -7.80950424e-005,
    5111                              -1.15922807e-005, -9.09134899e-005,
    5112                              -1.35994542e-004, -1.95673476e-004,
    5113                              -4.25779199e-004, -2.95890312e-004,
    5114                              -4.00060341e-004, -9.42021290e-005,
    5115                              -3.41372596e-004, -1.54560195e-004,
    5116                              -2.94810038e-004, -1.08844546e-004,
    5117                              -6.97240892e-005,  3.50299623e-005,
    5118                              -2.40159184e-004, -2.01805883e-004,
    5119                              -7.60732405e-004, -5.10897642e-004,
    5120                              -1.00940001e-003, -1.38037759e-004,
    5121                              -1.06169131e-003, -3.12307760e-004,
    5122                              -9.90602307e-004, -4.21634250e-005,
    5123                              -6.02424239e-004,  1.52230578e-004,
    5124                              -7.63833035e-004, -1.10273481e-004,
    5125                              -1.40187071e-003, -5.57831837e-004,
    5126                              -1.63988285e-003, -2.48018092e-004,
    5127                              -1.83309840e-003, -6.19360836e-004,
    5128                              -1.29955242e-003, -3.76237145e-004,
    5129                              -1.00613007e-003, -8.63641918e-005,
    5130                              -1.13604124e-003, -3.90589728e-004,
    5131                              -1.91457355e-003, -9.43783961e-004,
    5132                              -2.28090840e-003, -5.79107025e-004,
    5133                              -1.54091533e-003, -2.39785792e-003,
    5134                              -2.47947427e-003, -2.02694009e-003,
    5135                              -2.10441194e-003, -1.82082650e-003,
    5136                              -1.80229336e-003, -2.10418336e-003,
    5137                              -1.93104408e-003, -2.23200334e-003,
    5138                              -1.57239706e-003, -1.31486358e-003,
    5139                              -1.17564993e-003, -2.85846494e-003,
    5140                              -3.52956754e-003, -5.12658193e-003,
    5141                              -6.24238960e-003, -6.01820113e-003,
    5142                              -6.09602201e-003, -5.04787190e-003,
    5143                              -4.59373845e-003, -3.01393146e-003,
    5144                              5.08550095e-004, -4.35896549e-004])
     4922     [-7.65882069e-005, -1.46087080e-004, -1.09630102e-004, -7.80950424e-005,
     4923      -1.15922807e-005, -9.09134899e-005, -1.35994542e-004, -1.95673476e-004,
     4924      -4.25779199e-004, -2.95890312e-004, -4.00060341e-004, -9.42021290e-005,
     4925      -3.41372596e-004, -1.54560195e-004, -2.94810038e-004, -1.08844546e-004,
     4926      -6.97240892e-005,  3.50299623e-005, -2.40159184e-004, -2.01805883e-004,
     4927      -7.60732405e-004, -5.10897642e-004, -1.00940001e-003, -1.38037759e-004,
     4928      -1.06169131e-003, -3.12307760e-004, -9.90602307e-004, -4.21634250e-005,
     4929      -6.02424239e-004,  1.52230578e-004, -7.63833035e-004, -1.10273481e-004,
     4930      -1.40187071e-003, -5.57831837e-004, -1.63988285e-003, -2.48018092e-004,
     4931      -1.83309840e-003, -6.19360836e-004, -1.29955242e-003, -3.76237145e-004,
     4932      -1.00613007e-003, -8.63641918e-005, -1.13604124e-003, -3.90589728e-004,
     4933      -1.91457355e-003, -9.43783961e-004, -2.28090840e-003, -5.79107025e-004,
     4934      -1.54091533e-003, -2.39785792e-003, -2.47947427e-003, -2.02694009e-003,
     4935      -2.10441194e-003, -1.82082650e-003, -1.80229336e-003, -2.10418336e-003,
     4936      -1.93104408e-003, -2.23200334e-003, -1.57239706e-003, -1.31486358e-003,
     4937      -1.17564993e-003, -2.85846494e-003, -3.52956754e-003, -5.12658193e-003,
     4938      -6.24238960e-003, -6.01820113e-003, -6.09602201e-003, -5.04787190e-003,
     4939      -4.59373845e-003, -3.01393146e-003,  5.08550095e-004, -4.35896549e-004])
    51454940
    51464941        os.remove(domain.get_name() + '.sww')
     
    65486343
    65496344
     6345    def test_volume_conservation_inflow(self):
     6346        """test_volume_conservation
     6347       
     6348        Test that total volume in domain is as expected, based on questions
     6349        raised by Petar Milevski in May 2009.
     6350       
     6351        This test adds inflow at a known rate and verifies that the total
     6352        terminal volume is as expected.
     6353       
     6354        """
     6355
     6356        verbose = False
     6357       
     6358
     6359        #---------------------------------------------------------------------
     6360        # Import necessary modules
     6361        #---------------------------------------------------------------------
     6362        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
     6363        from anuga.shallow_water import Domain
     6364        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
     6365        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
     6366        from anuga.shallow_water.shallow_water_domain import Inflow
     6367        from anuga.shallow_water.data_manager import get_flow_through_cross_section
     6368
     6369        #----------------------------------------------------------------------
     6370        # Setup computational domain
     6371        #----------------------------------------------------------------------
     6372        finaltime = 200.0
     6373
     6374        length = 300.
     6375        width  = 20.
     6376        dx = dy = 5       # Resolution: of grid on both axes
     6377       
     6378
     6379        points, vertices, boundary = rectangular_cross(int(length/dx),
     6380                                                       int(width/dy),
     6381                                                       len1=length, len2=width)
     6382
     6383
     6384        domain = Domain(points, vertices, boundary)   
     6385        domain.set_name('Inflow_volume_test')              # Output name
     6386               
     6387
     6388        #----------------------------------------------------------------------
     6389        # Setup initial conditions
     6390        #----------------------------------------------------------------------
     6391        slope = 0.0
     6392        def topography(x, y):
     6393            z=-x * slope
     6394            return z
     6395
     6396        domain.set_quantity('elevation', topography) # Use function for elevation
     6397        domain.set_quantity('friction', 0.0)         # Constant friction
     6398               
     6399        domain.set_quantity('stage',
     6400                            expression='elevation')  # Dry initially
     6401                           
     6402
     6403        #--------------------------------------------------------------
     6404        # Setup Inflow
     6405        #--------------------------------------------------------------
     6406
     6407        # Fixed Flowrate onto Area
     6408        fixed_inflow = Inflow(domain,
     6409                              center=(10.0, 10.0),
     6410                              radius=5.00,
     6411                              rate=10.00)                               
     6412                           
     6413        domain.forcing_terms.append(fixed_inflow)                           
     6414       
     6415        #----------------------------------------------------------------------
     6416        # Setup boundary conditions
     6417        #----------------------------------------------------------------------
     6418
     6419        Br = Reflective_boundary(domain) # Solid reflective wall
     6420        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
     6421
     6422       
     6423        #----------------------------------------------------------------------
     6424        # Evolve system through time
     6425        #----------------------------------------------------------------------
     6426        ref_volume = 0.0
     6427        ys = 10.0  # Yieldstep
     6428        for t in domain.evolve(yieldstep=ys, finaltime=finaltime):
     6429       
     6430            # Check volume
     6431            assert num.allclose(domain.compute_total_volume(), ref_volume)
     6432       
     6433            if verbose :
     6434                print domain.timestepping_statistics()
     6435                print domain.volumetric_balance_statistics()
     6436                print 'reference volume', ref_volume
     6437           
     6438           
     6439            # Update reference volume
     6440            ref_volume += ys * fixed_inflow.rate
     6441
     6442
     6443    def test_volume_conservation_rain(self):
     6444        """test_volume_conservation
     6445       
     6446        Test that total volume in domain is as expected, based on questions
     6447        raised by Petar Milevski in May 2009.
     6448       
     6449        This test adds rain at a known rate and verifies that the total
     6450        terminal volume is as expected.
     6451       
     6452        """
     6453
     6454        verbose = False
     6455       
     6456
     6457        #---------------------------------------------------------------------
     6458        # Import necessary modules
     6459        #---------------------------------------------------------------------
     6460        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
     6461        from anuga.shallow_water import Domain
     6462        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
     6463        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
     6464        from anuga.shallow_water.shallow_water_domain import Rainfall
     6465        from anuga.shallow_water.data_manager import get_flow_through_cross_section
     6466
     6467        #----------------------------------------------------------------------
     6468        # Setup computational domain
     6469        #----------------------------------------------------------------------
     6470        finaltime = 200.0
     6471
     6472        length = 300.
     6473        width  = 20.
     6474        dx = dy = 5       # Resolution: of grid on both axes
     6475       
     6476
     6477        points, vertices, boundary = rectangular_cross(int(length/dx),
     6478                                                       int(width/dy),
     6479                                                       len1=length, len2=width)
     6480
     6481
     6482        domain = Domain(points, vertices, boundary)   
     6483        domain.set_name('Rain_volume_test')              # Output name
     6484               
     6485
     6486        #----------------------------------------------------------------------
     6487        # Setup initial conditions
     6488        #----------------------------------------------------------------------
     6489        slope = 0.0
     6490        def topography(x, y):
     6491            z=-x * slope
     6492            return z
     6493
     6494        domain.set_quantity('elevation', topography) # Use function for elevation
     6495        domain.set_quantity('friction', 0.0)         # Constant friction
     6496               
     6497        domain.set_quantity('stage',
     6498                            expression='elevation')  # Dry initially
     6499                           
     6500
     6501        #--------------------------------------------------------------
     6502        # Setup rain
     6503        #--------------------------------------------------------------
     6504
     6505        # Fixed rain onto small circular area
     6506        fixed_rain = Rainfall(domain,
     6507                              center=(10.0, 10.0),
     6508                              radius=5.00,
     6509                              rate=10.00)   # 10 mm/s                           
     6510                           
     6511        domain.forcing_terms.append(fixed_rain)                           
     6512       
     6513        #----------------------------------------------------------------------
     6514        # Setup boundary conditions
     6515        #----------------------------------------------------------------------
     6516
     6517        Br = Reflective_boundary(domain) # Solid reflective wall
     6518        domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
     6519
     6520       
     6521        #----------------------------------------------------------------------
     6522        # Evolve system through time
     6523        #----------------------------------------------------------------------
     6524        ref_volume = 0.0
     6525        ys = 10.0  # Yieldstep
     6526        for t in domain.evolve(yieldstep=ys, finaltime=finaltime):
     6527       
     6528            # Check volume
     6529            V = domain.compute_total_volume()
     6530            msg = 'V = %e, Ref = %e' % (V, ref_volume)
     6531            assert num.allclose(V, ref_volume), msg
     6532       
     6533            if verbose :
     6534                print domain.timestepping_statistics()
     6535                print domain.volumetric_balance_statistics()
     6536                print 'reference volume', ref_volume
     6537                print V
     6538           
     6539           
     6540            # Update reference volume.
     6541            # FIXME: Note that rate has now been redefined
     6542            # as m/s internally. This is a little confusing
     6543            # when it was specfied as mm/s.
     6544           
     6545            delta_V = fixed_rain.rate*fixed_rain.exchange_area
     6546            ref_volume += ys * delta_V
     6547
     6548
     6549    def Xtest_rain_conservation_and_runoff(self):
     6550        """test_rain_conservation_and_runoff
     6551       
     6552        Test that total volume in domain is as expected, based on questions
     6553        raised by Petar Milevski in May 2009.
     6554       
     6555        This test adds rain at a known rate and verifies that the total
     6556        volume and outflows are as expected.
     6557       
     6558        """
     6559
     6560        # FIXME (Ole): Does not work yet. Investigate boundary flows
     6561       
     6562        verbose = True #False
     6563       
     6564
     6565        #---------------------------------------------------------------------
     6566        # Import necessary modules
     6567        #---------------------------------------------------------------------
     6568        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
     6569        from anuga.shallow_water import Domain
     6570        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
     6571        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
     6572        from anuga.shallow_water.shallow_water_domain import Rainfall
     6573        from anuga.shallow_water.data_manager import get_flow_through_cross_section
     6574
     6575        #----------------------------------------------------------------------
     6576        # Setup computational domain
     6577        #----------------------------------------------------------------------
     6578        finaltime = 500.0
     6579
     6580        length = 300.
     6581        width  = 20.
     6582        dx = dy = 5       # Resolution: of grid on both axes
     6583       
     6584
     6585        points, vertices, boundary = rectangular_cross(int(length/dx),
     6586                                                       int(width/dy),
     6587                                                       len1=length, len2=width)
     6588
     6589
     6590        domain = Domain(points, vertices, boundary)   
     6591        domain.set_name('Rain_volume_runoff_test')         # Output name
     6592               
     6593
     6594        #----------------------------------------------------------------------
     6595        # Setup initial conditions
     6596        #----------------------------------------------------------------------
     6597        slope = 0.0
     6598        def topography(x, y):
     6599            z=-x * slope
     6600            return z
     6601
     6602        domain.set_quantity('elevation', topography) # Use function for elevation
     6603        domain.set_quantity('friction', 0.0)         # Constant friction
     6604               
     6605        domain.set_quantity('stage',
     6606                            expression='elevation')  # Dry initially
     6607                           
     6608
     6609        #--------------------------------------------------------------
     6610        # Setup rain
     6611        #--------------------------------------------------------------
     6612
     6613        # Fixed rain onto small circular area
     6614        fixed_rain = Rainfall(domain,
     6615                              center=(10.0, 10.0),
     6616                              radius=5.00,
     6617                              rate=10.00)   # 10 mm/s                           
     6618                           
     6619        domain.forcing_terms.append(fixed_rain)                           
     6620       
     6621        #----------------------------------------------------------------------
     6622        # Setup boundary conditions
     6623        #----------------------------------------------------------------------
     6624
     6625        Br = Reflective_boundary(domain) # Solid reflective wall
     6626        Bt = Transmissive_stage_zero_momentum_boundary(domain)
     6627        Bd = Dirichlet_boundary([-10, 0, 0])
     6628        domain.set_boundary({'left': Bt, 'right': Bd, 'top': Bt, 'bottom': Bt})
     6629
     6630       
     6631        #----------------------------------------------------------------------
     6632        # Evolve system through time
     6633        #----------------------------------------------------------------------
     6634        ref_volume = 0.0
     6635        ys = 10.0  # Yieldstep
     6636        for t in domain.evolve(yieldstep=ys, finaltime=finaltime):
     6637       
     6638            # Check volume
     6639            V = domain.compute_total_volume()
     6640            msg = 'V = %e, Ref = %e' % (V, ref_volume)
     6641            #assert num.allclose(V, ref_volume) or V < ref_volume, msg
     6642       
     6643            if verbose:
     6644                print domain.timestepping_statistics()
     6645                print domain.volumetric_balance_statistics()
     6646                print 'reference volume', ref_volume
     6647                print V
     6648           
     6649           
     6650            # Update reference volume.
     6651            # FIXME: Note that rate has now been redefined
     6652            # as m/s internally. This is a little confusing
     6653            # when it was specfied as mm/s.
     6654           
     6655            delta_V = fixed_rain.rate*fixed_rain.exchange_area
     6656            ref_volume += ys * delta_V
     6657       
     6658            # Compute outflow at right hand downstream boundary
     6659            boundary_flows, inflow , outflow = domain.compute_boundary_flows()
     6660            net_outflow = outflow - inflow
     6661       
     6662            outflow = boundary_flows['right']
     6663            if verbose:
     6664                print 'Outflow', outflow
     6665                print 'Net outflow', net_outflow
     6666       
     6667            # Update reference volume
     6668            ref_volume += ys * outflow           
     6669
     6670       
     6671               
     6672       
    65506673    def test_inflow_using_flowline(self):
    65516674        """test_inflow_using_flowline
     
    67076830                msg = ('Predicted flow was %f, should have been %f'
    67086831                       % (q, ref_flow))
     6832                assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
     6833
     6834       
     6835    def Xtest_inflow_boundary_using_flowline(self):
     6836        """test_inflow_boundary_using_flowline
     6837        Test the ability of a flowline to match inflow above the flowline by
     6838        creating constant inflow into the boundary at the head of a 20m
     6839        wide by 300m long plane dipping at various slopes with a
     6840        perpendicular flowline and gauge downstream of the inflow and
     6841        a 45 degree flowlines at 200m downstream
     6842       
     6843       
     6844        """
     6845
     6846        # FIXME (Ole): Work in progress
     6847       
     6848        verbose = False
     6849       
     6850
     6851        #----------------------------------------------------------------------
     6852        # Import necessary modules
     6853        #----------------------------------------------------------------------
     6854        from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
     6855        from anuga.shallow_water import Domain
     6856        from anuga.shallow_water.shallow_water_domain import Reflective_boundary
     6857        from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
     6858        from anuga.shallow_water.shallow_water_domain import Inflow_boundary
     6859        from anuga.shallow_water.data_manager import get_flow_through_cross_section
     6860        from anuga.abstract_2d_finite_volumes.util import sww2csv_gauges, csv2timeseries_graphs
     6861
     6862
     6863        #----------------------------------------------------------------------
     6864        # Setup computational domain
     6865        #----------------------------------------------------------------------
     6866
     6867        finaltime = 500 #700.0 # If this is too short, steady state will not be achieved
     6868
     6869        length = 250.
     6870        width  = 20.
     6871        dx = dy = 5          # Resolution: of grid on both axes
     6872
     6873        points, vertices, boundary = rectangular_cross(int(length/dx), int(width/dy),
     6874                                                       len1=length, len2=width)
     6875
     6876        for mannings_n in [0.1, 0.01]:
     6877            # Loop over a range of roughnesses             
     6878           
     6879            for slope in [1.0/300, 1.0/100]:
     6880                # Loop over a range of bedslopes representing sub to super critical flows
     6881
     6882
     6883                domain = Domain(points, vertices, boundary)   
     6884                domain.set_name('inflow_boundary_flowline_test')
     6885               
     6886
     6887                #-------------------------------------------------------------
     6888                # Setup initial conditions
     6889                #-------------------------------------------------------------
     6890
     6891                def topography(x, y):
     6892                    z=-x * slope
     6893                    return z
     6894
     6895                domain.set_quantity('elevation', topography)
     6896                domain.set_quantity('friction', mannings_n)
     6897                domain.set_quantity('stage',
     6898                                    expression='elevation')
     6899
     6900
     6901                   
     6902                #--------------------------------------------------------------
     6903                # Setup boundary conditions
     6904                #--------------------------------------------------------------
     6905               
     6906
     6907
     6908                ref_flow = 10.00
     6909
     6910                # Compute normal depth on plane using Mannings equation
     6911                # v=1/n*(r^2/3)*(s^0.5) or r=(Q*n/(s^0.5*W))^0.6
     6912                normal_depth=(ref_flow*mannings_n/(slope**0.5*width))**0.6
     6913                if verbose:
     6914                    print
     6915                    print 'Slope:', slope, 'Mannings n:', mannings_n
     6916                   
     6917               
     6918               
     6919                Bi = Inflow_boundary(domain, rate=ref_flow)
     6920
     6921                Br = Reflective_boundary(domain)
     6922               
     6923                # Define downstream boundary based on predicted depth
     6924                def normal_depth_stage_downstream(t):
     6925                    return (-slope*length) + normal_depth
     6926               
     6927                Bt = Transmissive_momentum_set_stage_boundary(domain=domain,
     6928                                                              function=normal_depth_stage_downstream)
     6929               
     6930
     6931               
     6932
     6933                domain.set_boundary({'left': Bi,
     6934                                     'right': Bt,
     6935                                     'top': Br,
     6936                                     'bottom': Br})
     6937
     6938
     6939
     6940                #--------------------------------------------------------------
     6941                # Evolve system through time
     6942                #--------------------------------------------------------------
     6943
     6944
     6945                for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
     6946                    pass
     6947                    #if verbose :
     6948                    #    print domain.timestepping_statistics()
     6949                    #    print domain.volumetric_balance_statistics()
     6950                                                       
     6951
     6952
     6953                #--------------------------------------------------------------
     6954                # Compute flow thru flowlines ds of inflow
     6955                #--------------------------------------------------------------
     6956                   
     6957                # Square on flowline at 200m
     6958                q=domain.get_flow_through_cross_section([[200.0,0.0],[200.0,20.0]])
     6959                msg = 'Predicted flow was %f, should have been %f' % (q, ref_flow)
     6960                if verbose:
     6961                    print '90 degree flowline: ANUGA = %f, Ref = %f' % (q, ref_flow)
     6962                assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
     6963
     6964                           
     6965                # 45 degree flowline at 200m
     6966                q=domain.get_flow_through_cross_section([[200.0,0.0],[220.0,20.0]])
     6967                msg = 'Predicted flow was %f, should have been %f' % (q, ref_flow)
     6968                if verbose:
     6969                    print '45 degree flowline: ANUGA = %f, Ref = %f' % (q, ref_flow)
     6970                   
    67096971                assert num.allclose(q, ref_flow, rtol=1.0e-2), msg         
    67106972
Note: See TracChangeset for help on using the changeset viewer.