Changeset 3293


Ignore:
Timestamp:
Jul 10, 2006, 11:23:11 AM (18 years ago)
Author:
jakeman
Message:

Updating lots of changes

Location:
development/pyvolution-1d
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • development/pyvolution-1d/domain.py

    r2791 r3293  
    706706
    707707        #Initialise interval of timestep sizes (for reporting only)
     708        # SEEMS WIERD
    708709        self.min_timestep = max_timestep
    709710        self.max_timestep = min_timestep
     
    747748            #Update boundary values
    748749            self.update_boundary()
     750
     751            #print 'timestep', self.timestep
    749752
    750753            #Update time
     
    801804            else:
    802805                raise 'Unknown order'
    803             Q.interpolate_from_vertices_to_edges()
     806            #Q.interpolate_from_vertices_to_edges()
    804807
    805808
  • development/pyvolution-1d/quantity.py

    r2791 r3293  
    518518        Qc = self.centroid_values
    519519        Qv = self.vertex_values
    520 
     520       
    521521        #Check each triangle
    522522        for k in range(self.domain.number_of_elements):
     
    526526            #vertex coordinates
    527527            x0, x1 = V[k,:]
    528 
    529528            #Extrapolate
    530529            Qv[k,0] = Qc[k] + G[k]*(x0-x)
    531530            Qv[k,1] = Qc[k] + G[k]*(x1-x)
    532 
    533 
    534 
    535531
    536532    def limit(self):
  • development/pyvolution-1d/shallow_water_1d.py

    r2791 r3293  
    6363
    6464        #forcing terms not included in 1d domain ?WHy?
    65         self.forcing_terms.append(gravity)
    66         self.forcing_terms.append(manning_friction)
     65        #self.forcing_terms.append(gravity)
     66        #self.forcing_terms.append(manning_friction)
    6767        #print "\nI have Removed forcing terms line 64 1dsw"
    6868
     
    176176
    177177        #Initialise real time viz if requested
    178         if self.visualise is True and self.time == 0.0:
    179             if self.visualiser is None:
    180                 self.initialise_visualiser()
    181 
    182             self.visualiser.update_timer()
    183             self.visualiser.setup_all()
     178        #if self.visualise is True and self.time == 0.0:
     179        #    if self.visualiser is None:
     180        #        self.initialise_visualiser()
     181        #
     182        #    self.visualiser.update_timer()
     183        #    self.visualiser.setup_all()
    184184
    185185        #Store model data, e.g. for visualisation
    186         if self.store is True and self.time == 0.0:
    187             self.initialise_storage()
    188             #print 'Storing results in ' + self.writer.filename
    189         else:
    190             pass
    191             #print 'Results will not be stored.'
    192             #print 'To store results set domain.store = True'
    193             #FIXME: Diagnostic output should be controlled by
    194             # a 'verbose' flag living in domain (or in a parent class)
     186        #if self.store is True and self.time == 0.0:
     187        #    self.initialise_storage()
     188        #    #print 'Storing results in ' + self.writer.filename
     189        #else:
     190        #    pass
     191        #    #print 'Results will not be stored.'
     192        #    #print 'To store results set domain.store = True'
     193        #    #FIXME: Diagnostic output should be controlled by
     194        #    # a 'verbose' flag living in domain (or in a parent class)
    195195
    196196        #Call basic machinery from parent class
     
    198198                                       skip_initial_step):
    199199            #Real time viz
    200             if self.visualise is True:
    201                 self.visualiser.update_all()
    202                 self.visualiser.update_timer()
     200        #    if self.visualise is True:
     201        #        self.visualiser.update_all()
     202        #        self.visualiser.update_timer()
    203203
    204204
    205205            #Store model data, e.g. for subsequent visualisation
    206             if self.store is True:
    207                 self.store_timestep(self.quantities_to_be_stored)
     206        #    if self.store is True:
     207        #        self.store_timestep(self.quantities_to_be_stored)
    208208
    209209            #FIXME: Could maybe be taken from specified list
     
    312312    q_right[1] = q_right[1]*normal
    313313
    314     z = (zl+zr)/2 #Take average of field values
     314    #z = (zl+zr)/2 #Take average of field values
     315    z = 0.5*(zl+zr) #Take average of field values
    315316
    316317    w_left  = q_left[0]   #w=h+z
     
    350351    #Flux computation
    351352
    352     #FIXME(Ole): Why is it again that we don't
    353     #use uh_left and uh_right directly in the first entries?
    354353    #flux_left  = array([u_left*h_left,
    355     #                    u_left*uh_left + 0.5*g*h_left**2,
    356     #                    u_left*vh_left])
     354    #                    u_left*uh_left + 0.5*g*h_left**2])
    357355    #flux_right = array([u_right*h_right,
    358     #                    u_right*uh_right + 0.5*g*h_right**2,
    359     #                    u_right*vh_right])
    360 
     356    #                    u_right*uh_right + 0.5*g*h_right**2])
    361357    flux_left  = array([u_left*h_left,
    362                         u_left*uh_left + 0.5*g*h_left**2])
     358                        u_left*uh_left + 0.5*g*h_left*h_left])
    363359    flux_right = array([u_right*h_right,
    364                         u_right*uh_right + 0.5*g*h_right**2])
     360                        u_right*uh_right + 0.5*g*h_right*h_right])
    365361
    366362    denom = s_max-s_min
    367363    if denom == 0.0:
    368         #edgeflux = array([0.0, 0.0, 0.0])
    369364        edgeflux = array([0.0, 0.0])
    370365        max_speed = 0.0
     
    373368        edgeflux += s_max*s_min*(q_right-q_left)/denom
    374369       
    375         #edgeflux[1] = -edgeflux[1]*normal
    376370        edgeflux[1] = edgeflux[1]*normal
    377371
    378         #edgeflux = rotate(edgeflux, normal, direction=-1)
    379372        max_speed = max(abs(s_max), abs(s_min))
    380373
     
    445438        for i in range(2):
    446439            #Quantities inside volume facing neighbour i
    447             #ql = [stage[k, i], xmom[k, i], ymom[k, i]]
    448             ql[0] = stage[k, i]
    449             ql[1] = xmom[k, i]
     440            #ql[0] = stage[k, i]
     441            #ql[1] = xmom[k, i]
     442            ql = [stage[k, i], xmom[k, i]]
    450443            zl = bed[k, i]
    451444
     
    454447            if n < 0:
    455448                m = -n-1 #Convert negative flag to index
    456                 #qr = [stage_bdry[m], xmom_bdry[m], ymom_bdry[m]]
    457449                qr[0] = stage_bdry[m]
    458450                qr[1] = xmom_bdry[m]
     
    469461            #Outward pointing normal vector
    470462            normal = domain.normals[k, i]
    471             #CHECK THIS LINE [k, 2*i:2*i+1]
    472 
     463       
    473464            #Flux computation using provided function
    474465            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
     
    486477            try:
    487478                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
    488                 timestep = min(timestep, 0.5/max_speed)
     479                #timestep = 0.01
     480                timestep = min(timestep, 0.5*domain.areas[k]/max_speed)
     481                if timestep < 0.00001:
     482                    #print 'max_speed', max_speed
     483                    s = domain.quantities['stage']
     484                    s = s.centroid_values
     485                    xm = domain.quantities['xmomentum']
     486                    xm = xm.centroid_values
     487                    #print 'h', s[k]
     488                    #print 'xm', xm[k]
     489                    #print 'u', xm[k]/s[k]
     490                    #break
     491                #timestep = 0.01
     492                #print 'areas', domain.areas[k]
     493                #print "timestep", timestep
    489494            except ZeroDivisionError:
    490495                pass
     
    492497        #Normalise by area and store for when all conserved
    493498        #quantities get updated
    494         flux /= domain.areas[k]
     499        #flux /= domain.areas[k]
     500        # ADD ABOVE LINE AGAIN
    495501        Stage.explicit_update[k] = flux[0]
    496502        Xmom.explicit_update[k] = flux[1]
     
    596602
    597603    #Remove very thin layers of water
    598     protect_against_infinitesimal_and_negative_heights(domain)
     604    #protect_against_infinitesimal_and_negative_heights(domain)
    599605
    600606    #Extrapolate all conserved quantities
    601     if optimised_gradient_limiter:
    602         #MH090605 if second order,
    603         #perform the extrapolation and limiting on
    604         #all of the conserved quantities
    605 
    606         if (domain.order == 1):
    607             for name in domain.conserved_quantities:
    608                 Q = domain.quantities[name]
    609                 Q.extrapolate_first_order()
     607    #if optimised_gradient_limiter:
     608    #    #MH090605 if second order,
     609    #    #perform the extrapolation and limiting on
     610    #    #all of the conserved quantities
     611
     612    #    if (domain.order == 1):
     613    #        for name in domain.conserved_quantities:
     614    #            Q = domain.quantities[name]
     615    #            Q.extrapolate_first_order()
     616    #    elif domain.order == 2:
     617    #        domain.extrapolate_second_order_sw()
     618    #    else:
     619    #        raise 'Unknown order'
     620    #else:
     621        #old code:
     622    for name in domain.conserved_quantities:
     623        Q = domain.quantities[name]
     624        if domain.order == 1:
     625            Q.extrapolate_first_order()
    610626        elif domain.order == 2:
    611             domain.extrapolate_second_order_sw()
     627            Q.extrapolate_second_order()
     628            Q.limit()
    612629        else:
    613630            raise 'Unknown order'
    614     else:
    615         #old code:
    616         for name in domain.conserved_quantities:
    617             Q = domain.quantities[name]
    618             if domain.order == 1:
    619                 Q.extrapolate_first_order()
    620             elif domain.order == 2:
    621                 Q.extrapolate_second_order()
    622                 Q.limit()
    623             else:
    624                 raise 'Unknown order'
    625 
    626631
    627632    #Take bed elevation into account when water heights are small
    628     balance_deep_and_shallow(domain)
     633    #balance_deep_and_shallow(domain)
    629634
    630635    #Compute edge values by interpolation
     
    664669    wc = domain.quantities['stage'].centroid_values
    665670    zc = domain.quantities['elevation'].centroid_values
    666 #    xmomc = domain.quantities['xmomentum'].centroid_values
    667     ymomc = domain.quantities['ymomentum'].centroid_values
     671    xmomc = domain.quantities['xmomentum'].centroid_values
     672#    ymomc = domain.quantities['ymomentum'].centroid_values
    668673
    669674    from shallow_water_ext import protect
  • development/pyvolution-1d/test_shallow_water_1d.py

    r2791 r3293  
    466466    def test_compute_fluxes1(self):
    467467        #Use values from previous version
    468 
    469         #a = [0.0, 0.0]
    470         #b = [0.0, 2.0]
    471         #c = [2.0,0.0]
    472         #d = [0.0, 4.0]
    473         #e = [2.0, 2.0]
    474         #f = [4.0,0.0]
    475         a=0.0
    476         b=2.0
    477         c=4.0
    478         d=6.0
    479         e=8.0
    480         #f=10.0
    481        
    482         #points = [a, b, c, d, e, f]
     468        print "test compute fluxes 1"
     469
     470        a=0.0
     471        b=2.0
     472        c=4.0
     473        d=6.0
     474        e=8.0
     475               
    483476        points = [a, b, c, d, e]
    484         #bac, bce, ecf, dbe
    485         #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    486 
    487         #domain = Domain(points, vertices)
     477       
    488478        domain = Domain(points)
    489479       
     
    493483        val3 = 2.+8.0/3
    494484
    495         #domain.set_quantity('stage', [[val0, val0, val0], [val1, val1, val1],
    496         #                              [val2, val2, val2], [val3, val3, val3]])
    497         domain.set_quantity('stage', [[val0, val0], [val1, val1],
    498                                       [val2, val2], [val3, val3]])
     485        domain.set_quantity('stage', [[val3, val3], [val1, val1],
     486                                      [val2, val2], [val0, val0]])
    499487        stage = domain.get_quantity('stage')
    500488        print stage
     
    504492
    505493        #Flux across right edge of volume 1
    506         #normal = domain.get_normal(1,0)
    507         #ql = domain.get_conserved_quantities(vol_id=1, edge=0)
    508         #qr = domain.get_conserved_quantities(vol_id=2, edge=2)
    509         #flux0, max_speed = flux_function(normal, ql, qr, zl, zr)
    510 
    511         #Flux across right edge of volume 1
    512         #ql = domain.get_conserved_quantities(vol_id=1, edge=1)
    513         #qr = domain.get_conserved_quantities(vol_id=2, edge=0)
    514494        ql = domain.get_conserved_quantities(vol_id=1, vertex=1)
    515495        qr = domain.get_conserved_quantities(vol_id=2, vertex=0)
     
    521501       
    522502        #Check that flux seen from other triangles is inverse
    523         #tmp = qr; qr=ql; ql=tmp
    524         #normal = domain.get_normal(3,0)
    525         #flux, max_speed = flux_function(normal, ql, qr, zl, zr)
    526503        flux, max_speed = flux_function(-1.0, qr, ql, zl, zr)
    527504        print flux0
     
    529506        assert allclose(flux + flux0, 0.)
    530507
    531         #print flux0
    532         #print max_speed
    533         #assert allclose(flux0, [-15.3598804, 253.71111111, 0.])
    534508        assert allclose(flux0, [-15.3598804, 253.71111111])
    535509        assert allclose(max_speed, 9.21592824046)
    536510
    537511
    538         #Flux across upper edge of volume 1
    539         #normal = domain.get_normal(1,1)
    540         #ql = domain.get_conserved_quantities(vol_id=1, edge=1)
    541         #qr = domain.get_conserved_quantities(vol_id=2, edge=0)
     512        #Flux across left edge of volume 1
    542513        ql = domain.get_conserved_quantities(vol_id=1, vertex=0)
    543514        qr = domain.get_conserved_quantities(vol_id=0, vertex=1)
    544         #flux1, max_speed = flux_function(normal, ql, qr, zl, zr)
    545515        flux1, max_speed = flux_function(-1.0, ql, qr, zl, zr)
    546         #assert allclose(flux1, [2.4098563, 0., 123.04444444])
    547516        print 'flux1', flux1
    548         assert allclose(flux1, [2.4098563, 0.])
     517        assert allclose(flux1, [2.4098563, -123.04444444])
    549518        assert allclose(max_speed, 7.22956891292)
    550519
    551         #Flux across lower left hypotenuse of volume 1
    552         #normal = domain.get_normal(1,2)
    553         #ql = domain.get_conserved_quantities(vol_id=1, edge=2)
    554         #qr = domain.get_conserved_quantities(vol_id=0, edge=1)
    555         #ql = domain.get_conserved_quantities(vol_id=1, vertex=2)
    556         #qr = domain.get_conserved_quantities(vol_id=0, vertex=1)
    557         #flux2, max_speed = flux_function(normal, ql, qr, zl, zr)
    558 
    559         #assert allclose(flux2, [9.63942522, -61.59685738, -61.59685738])
    560         #assert allclose(max_speed, 7.22956891292)
    561 
    562         #Scale, add up and check that compute_fluxes is correct for vol 1
    563         #e0 = domain.edgelengths[1, 0]
    564         #e1 = domain.edgelengths[1, 1]
    565         #e2 = domain.edgelengths[1, 2]
    566 
    567         #total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1]
     520        #Add up and check that compute_fluxes is correct for vol 1
    568521        total_flux = -(flux0+flux1)/domain.areas[1]
    569         #assert allclose(total_flux, [-0.68218178, -166.6, -35.93333333])
    570         assert allclose(total_flux, [-0.68218178, -166.6])
    571        
    572 
    573 
     522        print 'total flux', total_flux
     523        print 'domain area',domain.areas[1]
     524        assert allclose(total_flux, [6.47501205, -65.333333333])
     525       
    574526        domain.compute_fluxes()
    575527
    576         #assert allclose(total_flux, domain.explicit_update[1,:])
    577         for i, name in enumerate(['stage', 'xmomentum', 'ymomentum']):
     528        for i, name in enumerate(['stage', 'xmomentum']):
    578529            assert allclose(total_flux[i],
    579530                            domain.quantities[name].explicit_update[1])
     
    585536        #    [-35.68522449, 0., 69.68888889]])
    586537
     538        print domain.quantities['stage'].explicit_update
     539        print domain.quantities['xmomentum'].explicit_update
    587540        assert allclose(domain.quantities['stage'].explicit_update,
    588                         [0., -0.68218178, -111.77316251, -35.68522449])
     541                        [0., 6.47501205, -111.77316251, -35.68522449])
    589542        assert allclose(domain.quantities['xmomentum'].explicit_update,
    590                         [-69.68888889, -166.6, 69.68888889, 0])
     543                        [-69.68888889, -65.33333333, 69.68888889, 0])
    591544        assert allclose(domain.quantities['ymomentum'].explicit_update,
    592545                        [-69.68888889, -35.93333333, 0., 69.68888889])
     
    597550    def test_catching_negative_heights(self):
    598551
    599         #a = [0.0, 0.0]
    600         #b = [0.0, 2.0]
    601         #c = [2.0,0.0]
    602         #d = [0.0, 4.0]
    603         #e = [2.0, 2.0]
    604         #f = [4.0,0.0]
    605         a=0.0
    606         b=2.0
    607         c=4.0
    608         d=6.0
    609         e=8.0
    610         #f=10.0
    611 
    612         #points = [a, b, c, d, e, f]
     552        a=0.0
     553        b=2.0
     554        c=4.0
     555        d=6.0
     556        e=8.0
     557   
    613558        points = [a, b, c, d, e]
    614         #bac, bce, ecf, dbe
    615         #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    616 
    617         #domain = Domain(points, vertices)
    618559        domain = Domain(points)
    619560        val0 = 2.+2.0/3
     
    623564
    624565        zl=zr=4  #Too large
    625         #domain.set_quantity('elevation', zl*ones( (4,3) ))
    626         #domain.set_quantity('stage', [[val0, val0-1, val0-2],
    627         #                              [val1, val1+1, val1],
    628         #                              [val2, val2-2, val2],
    629         #                              [val3-0.5, val3, val3]])
    630566        domain.set_quantity('elevation', zl*ones( (4,2) ))
    631567        domain.set_quantity('stage', [[val0, val0-1],
     
    703639
    704640        from config import g
    705         #a = [0.0, 0.0]
    706         #b = [0.0, 2.0]
    707         #c = [2.0, 0.0]
    708         #d = [0.0, 4.0]
    709         #e = [2.0, 2.0]
    710         #f = [4.0, 0.0]
    711         a=0.0
    712         b=2.0
    713         c=4.0
    714         d=6.0
    715         e=8.0
    716         #f=10.0
    717 
    718         #points = [a, b, c, d, e, f]
     641        a=0.0
     642        b=2.0
     643        c=4.0
     644        d=6.0
     645        e=8.0
     646       
    719647        points = [a, b, c, d, e]
    720         #bac, bce, ecf, dbe
    721         #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    722 
    723         #domain = Domain(points, vertices)
     648       
    724649        domain = Domain(points)
    725650
     
    751676        from config import g
    752677
    753         #a = [0.0, 0.0]
    754         #b = [0.0, 2.0]
    755         #c = [2.0, 0.0]
    756         #d = [0.0, 4.0]
    757         #e = [2.0, 2.0]
    758         #f = [4.0, 0.0]
    759         a=0.0
    760         b=2.0
    761         c=4.0
    762         d=6.0
    763         e=8.0
    764         #f=10.0
    765 
    766         #points = [a, b, c, d, e, f]
     678        a=0.0
     679        b=2.0
     680        c=4.0
     681        d=6.0
     682        e=8.0
     683
    767684        points = [a, b, c, d, e]
    768         #bac, bce, ecf, dbe
    769         #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    770 
    771         #domain = Domain(points, vertices)
     685               
    772686        domain = Domain(points)
    773687
     
    10861000        #Using test data generated by pyvolution-2
    10871001        #Assuming no friction and flat bed (0.0)
    1088 
    1089         #a = [0.0, 0.0]
    1090         #b = [0.0, 2.0]
    1091         #c = [2.0,0.0]
    1092         #d = [0.0, 4.0]
    1093         #e = [2.0, 2.0]
    1094         #f = [4.0,0.0]
    1095         a=0.0
    1096         b=2.0
    1097         c=4.0
    1098         d=6.0
    1099         e=8.0
    1100         #f=10.0
    1101 
    1102         #points = [a, b, c, d, e, f]
     1002        print "\ntest distriute basic"
     1003
     1004        a=0.0
     1005        b=2.0
     1006        c=4.0
     1007        d=6.0
     1008        e=8.0
     1009
    11031010        points = [a, b, c, d, e]
    1104         #bac, bce, ecf, dbe
    1105         #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    1106 
    1107         #domain = Domain(points, vertices)
    11081011        domain = Domain(points)
    11091012
     
    11201023        domain.order = 1
    11211024        domain.distribute_to_vertices_and_edges()
     1025        print "First order L", L
    11221026        assert allclose(L[1], val1)
    11231027
    1124         #Second order
     1028        #Second order       
    11251029        domain.order = 2
     1030        a = domain.quantities['stage'].compute_gradients()
     1031        print 'gradients', a
     1032               
    11261033        domain.distribute_to_vertices_and_edges()
    1127         assert allclose(L[1], [2.2, 4.9, 4.9])
    1128 
    1129 
     1034        print "Second order L", L
     1035        print "L[1]", L[1]
     1036        #assert allclose(L[1], [2.2, 4.9, 4.9])
     1037        assert allclose(L[1], [2.5, 5.5])
    11301038
    11311039    def test_distribute_away_from_bed(self):
     
    11331041        #Assuming no friction and flat bed (0.0)
    11341042
    1135         #a = [0.0, 0.0]
    1136         #b = [0.0, 2.0]
    1137         #c = [2.0,0.0]
    1138         #d = [0.0, 4.0]
    1139         #e = [2.0, 2.0]
    1140         #f = [4.0,0.0]
     1043        print "\ntest distribute awaway from bed"
     1044
    11411045        a=0.0
    11421046        b=2.0
     
    11481052        #points = [a, b, c, d, e, f]
    11491053        points = [a, b, c, d, e]
    1150         #bac, bce, ecf, dbe
    1151         #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    1152 
    1153         #domain = Domain(points, vertices)
     1054
    11541055        domain = Domain(points)
    11551056
     
    11781079        assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778])
    11791080
    1180 
    1181 
    1182 
    1183 
     1081    def test_1d_solution_I(self):
     1082        print "TEST 1D-SOLUTION I"
     1083
     1084        L = 2000.0     # Length of channel (m)
     1085        N = 100        # Number of compuational cells
     1086        cell_len = L/N # Origin = 0.0
     1087       
     1088        points = zeros(N+1,Float)
     1089        for i in range(N+1):
     1090            points[i] = i*cell_len
     1091           
     1092        domain = Domain(points)
     1093           
     1094        def stage(x):
     1095            for i in range(len(x)):
     1096                if x[i]<=1000.0:
     1097                    x[i] = 10.0
     1098                else:
     1099                    x[i] = 5.0
     1100            return x
     1101               
     1102        domain.set_quantity('stage', stage)
     1103        #L = domain.quantities['stage'].vertex_values
     1104        #print "Initial Stage"
     1105        #print L
     1106        domain.set_boundary({'exterior': Reflective_boundary(domain)})
     1107
     1108        import time
     1109        t0 = time.time()
     1110        yieldstep = 10
     1111        finaltime = 50.0
     1112        for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
     1113            pass
     1114
     1115        print 'That took %.2f seconds' %(time.time()-t0)
     1116
     1117        #L = domain.quantities['stage'].vertex_values
     1118        #print "Final Stage"
     1119        #print L
     1120
     1121        C = domain.quantities['stage'].vertex_values
     1122        #print C
     1123
     1124        f = file('test_solution_I.out', 'w')
     1125        for i in range(N):
     1126            f.write(str(C[i,1]))
     1127            f.write("\n")
     1128        f.close
     1129"""
     1130    def test_1d_solution_II(self):
     1131        print "TEST 1D-SOLUTION II"
     1132       
     1133        L = 2000.0     # Length of channel (m)
     1134        N = 100        # Number of compuational cells
     1135        cell_len = L/N # Origin = 0.0
     1136       
     1137        points = zeros(N+1,Float)
     1138        for i in range(N+1):
     1139            points[i] = i*cell_len
     1140           
     1141        domain = Domain(points)
     1142           
     1143        def stage(x):
     1144            for i in range(len(x)):
     1145                if x[i]<=1000.0:
     1146                    x[i] = 10.0
     1147                else:
     1148                    x[i] = 0.1
     1149            return x
     1150                   
     1151
     1152        domain.set_quantity('stage', stage)
     1153        #L = domain.quantities['stage'].vertex_values
     1154        #print "Initial Stage"
     1155        #print L
     1156        domain.set_boundary({'exterior': Reflective_boundary(domain)})
     1157       
     1158        import time
     1159        t0 = time.time()
     1160        yieldstep = 1.0
     1161        finaltime = 50.0
     1162       
     1163        for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
     1164            a=0
     1165
     1166        print 'That took %.2f seconds' %(time.time()-t0)
     1167        #L = domain.quantities['stage'].vertex_values
     1168        #print "Final Stage"
     1169        #print L
     1170
     1171        C = domain.quantities['stage'].vertex_values
     1172        #print C
     1173       
     1174        f = file('test_solution_II.out', 'w')
     1175        for i in range(N):
     1176            f.write(str(C[i,1]))
     1177            f.write("\n")
     1178        f.close
     1179
     1180    def test_1d_solution_III(self):
     1181        print "TEST 1D-SOLUTION III"
     1182        L = 2000.0     # Length of channel (m)
     1183        N = 100        # Number of compuational cells
     1184        cell_len = L/N # Origin = 0.0
     1185       
     1186        points = zeros(N+1,Float)
     1187        for i in range(N+1):
     1188            points[i] = i*cell_len
     1189           
     1190        domain = Domain(points)
     1191
     1192        def stage(x):
     1193            for i in range(len(x)):
     1194                if x[i]<=1000.0:
     1195                    x[i] = 10.0
     1196                else:
     1197                    x[i] = 0.0
     1198            return x
     1199
     1200        domain.set_quantity('stage', stage)
     1201        #L = domain.quantities['stage'].vertex_values
     1202        #print "Initial Stage"
     1203        #print L
     1204        domain.set_boundary({'exterior': Reflective_boundary(domain)})
     1205       
     1206        import time
     1207        t0 = time.time()
     1208        yieldstep = 1.0
     1209        finaltime = 30.0
     1210       
     1211        for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
     1212            a=0
     1213           
     1214        print 'That took %.2f seconds' %(time.time()-t0)
     1215           
     1216        #L = domain.quantities['stage'].vertex_values
     1217        #print "Final Stage"
     1218        #print L
     1219       
     1220        C = domain.quantities['stage'].vertex_values
     1221        #print C
     1222
     1223        f = file('test_solution_III.out', 'w')
     1224        for i in range(N):
     1225            f.write(str(C[i,1]))
     1226            f.write("\n")
     1227        f.close
     1228   """                                         
    11841229#-------------------------------------------------------------
    11851230if __name__ == "__main__":
Note: See TracChangeset for help on using the changeset viewer.