Changeset 215


Ignore:
Timestamp:
Aug 24, 2004, 6:11:26 PM (21 years ago)
Author:
ole
Message:

More testing of gravity, forcing terms, distributors

Location:
inundation/ga/storm_surge/pyvolution
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • inundation/ga/storm_surge/pyvolution/domain.py

    r212 r215  
    438438        timestep = self.timestep
    439439
    440         #Reset semi_implicit_updates
    441         #(explicit updates are already set by compute_fluyzes)
    442         for i, name in enumerate(self.conserved_quantities):
    443             Q = self.quantities[name]
    444             Q.semi_implicit_update[:] = 0.
    445 
    446440        #Compute forcing terms   
    447441        self.compute_forcing_terms()   
    448442
    449443        #Update conserved_quantities from explicit updates
    450         for i, name in enumerate(self.conserved_quantities):
     444        for name in self.conserved_quantities:
    451445            Q = self.quantities[name]
    452446            Q.update(timestep)
     447
     448            #Clean up
     449            Q.semi_implicit_update[:] = 0.0
     450            Q.explicit_update[:] = 0.0 #Not necessary as fluxes will set it
    453451           
    454452
  • inundation/ga/storm_surge/pyvolution/quantity.py

    r209 r215  
    328328            msg = 'Invalid location: %s' %location
    329329            raise msg
     330
     331        if X is None:
     332            msg = 'Given values are None'
     333            raise msg           
    330334       
    331335        import types
  • inundation/ga/storm_surge/pyvolution/shallow_water.py

    r212 r215  
    586586    for k in range(domain.number_of_elements):   
    587587        avg_h = sum( h[k,:] )/3
    588 
     588       
    589589        #Compute bed slope
    590590        x0, y0, x1, y1, x2, y2 = V[k,:]   
    591591        z0, z1, z2 = Elevation.vertex_values[k,:]
     592
    592593        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
    593594
     
    595596        Xmom.explicit_update[k] += -g*zx*avg_h
    596597        Ymom.explicit_update[k] += -g*zy*avg_h       
    597    
     598
    598599
    599600def manning_friction(domain):
     
    627628
    628629###########################
     630###########################
     631#Geometries
     632
     633
     634#FIXME: Rethink this way of creating values.
     635
     636
     637class Weir:
     638    """Set a bathymetry for weir with a hole and a downstream gutter
     639    x,y are assumed to be in the unit square
     640    """
     641   
     642    def __init__(self, stage):
     643        self.inflow_stage = stage
     644
     645    def __call__(self, x, y):   
     646        from Numeric import zeros, Float
     647        from math import sqrt
     648   
     649        N = len(x)
     650        assert N == len(y)   
     651
     652        z = zeros(N, Float)
     653        for i in range(N):
     654            z[i] = -x[i]/2  #General slope
     655
     656            #Flattish bit to the left
     657            if x[i] < 0.3:
     658                z[i] = -x[i]/10
     659           
     660            #Weir
     661            if x[i] >= 0.3 and x[i] < 0.4:
     662                z[i] = -x[i]+0.9
     663
     664            #Dip
     665            x0 = 0.6
     666            #depth = -1.3
     667            depth = -1.0
     668            #plateaux = -0.9
     669            plateaux = -0.6           
     670            if y[i] < 0.7:
     671                if x[i] > x0 and x[i] < 0.9:
     672                    z[i] = depth
     673
     674                #RHS plateaux
     675                if x[i] >= 0.9:
     676                    z[i] = plateaux
     677                   
     678                   
     679            elif y[i] >= 0.7 and y[i] < 1.5:
     680                #Restrict and deepen
     681                if x[i] >= x0 and x[i] < 0.8:
     682                    z[i] = depth-(y[i]/3-0.3)
     683                    #z[i] = depth-y[i]/5
     684                    #z[i] = depth
     685                elif x[i] >= 0.8:   
     686                    #RHS plateaux
     687                    z[i] = plateaux
     688                   
     689            elif y[i] >= 1.5:
     690                if x[i] >= x0 and x[i] < 0.8 + (y[i]-1.5)/1.2:
     691                    #Widen up and stay at constant depth
     692                    z[i] = depth-1.5/5
     693                elif x[i] >= 0.8 + (y[i]-1.5)/1.2:                       
     694                    #RHS plateaux
     695                    z[i] = plateaux                                   
     696                   
     697
     698            #Hole in weir (slightly higher than inflow condition)
     699            if x[i] >= 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
     700                z[i] = -x[i]+self.inflow_stage + 0.02
     701
     702            #Channel behind weir
     703            x0 = 0.5
     704            if x[i] >= 0.4 and x[i] < x0 and y[i] > 0.2 and y[i] < 0.4:
     705                z[i] = -x[i]+self.inflow_stage + 0.02
     706               
     707            if x[i] >= x0 and x[i] < 0.6 and y[i] > 0.2 and y[i] < 0.4:
     708                #Flatten it out towards the end
     709                z[i] = -x0+self.inflow_stage + 0.02 + (x0-x[i])/5
     710
     711            #Hole to the east
     712            x0 = 1.1; y0 = 0.35
     713            #if x[i] < -0.2 and y < 0.5:
     714            if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
     715                z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-1.0
     716
     717            #Tiny channel draining hole
     718            if x[i] >= 1.14 and x[i] < 1.2 and y[i] >= 0.4 and y[i] < 0.6:
     719                z[i] = -0.9 #North south
     720               
     721            if x[i] >= 0.9 and x[i] < 1.18 and y[i] >= 0.58 and y[i] < 0.65:
     722                z[i] = -1.0 + (x[i]-0.9)/3 #East west
     723           
     724           
     725
     726            #Stuff not in use
     727           
     728            #Upward slope at inlet to the north west
     729            #if x[i] < 0.0: # and y[i] > 0.5:
     730            #    #z[i] = -y[i]+0.5  #-x[i]/2
     731            #    z[i] = x[i]/4 - y[i]**2 + 0.5
     732
     733            #Hole to the west
     734            #x0 = -0.4; y0 = 0.35 # center
     735            #if sqrt((2*(x[i]-x0))**2 + (2*(y[i]-y0))**2) < 0.2:
     736            #    z[i] = sqrt(((x[i]-x0))**2 + ((y[i]-y0))**2)-0.2
     737
     738
     739
     740
     741
     742        return z/2
     743
     744class Weir_simple:
     745    """Set a bathymetry for weir with a hole and a downstream gutter
     746    x,y are assumed to be in the unit square
     747    """
     748   
     749    def __init__(self, stage):
     750        self.inflow_stage = stage
     751
     752    def __call__(self, x, y):   
     753        from Numeric import zeros, Float
     754   
     755        N = len(x)
     756        assert N == len(y)   
     757
     758        z = zeros(N, Float)
     759        for i in range(N):
     760            z[i] = -x[i]  #General slope
     761
     762            #Flat bit to the left
     763            if x[i] < 0.3:
     764                z[i] = -x[i]/10  #General slope
     765           
     766            #Weir
     767            if x[i] > 0.3 and x[i] < 0.4:
     768                z[i] = -x[i]+0.9
     769
     770            #Dip
     771            if x[i] > 0.6 and x[i] < 0.9:
     772                z[i] = -x[i]-0.5  #-y[i]/5
     773
     774            #Hole in weir (slightly higher than inflow condition)
     775            if x[i] > 0.3 and x[i] < 0.4 and y[i] > 0.2 and y[i] < 0.4:
     776                z[i] = -x[i]+self.inflow_stage + 0.05       
     777           
     778
     779        return z/2
     780       
     781
     782           
    629783class Constant_height:
    630784    """Set an initial condition with constant water height, e.g
    631785    stage s = z+h
    632786    """
    633 
    634     #FIXME: Rethink this way of creating values.
    635 
    636787    def __init__(self, W, h):
    637788        self.W = W
  • inundation/ga/storm_surge/pyvolution/test_domain.py

    r195 r215  
    165165
    166166
    167     def test_distribute(self):
     167    def test_distribute_first_order(self):
    168168        """Domain implements a default first order gradient limiter
    169169        """
     
    206206
    207207
    208         #print domain.quantities['level'].centroid_values
    209208        domain.distribute_to_vertices_and_edges()
    210209
     
    244243
    245244
    246         domain.set_quantity('level', [[1,2,3], [5,5,5],
    247                                       [0,0,9], [-6, 3, 3]])
    248 
    249         domain.set_quantity('xmomentum', [[1,1,1], [2,2,2],
    250                                           [3,3,3], [4, 4, 4]])
    251 
    252         domain.set_quantity('ymomentum', [[10,10,10], [20,20,20],
    253                                           [30,30,30], [40, 40, 40]])       
    254 
     245        domain.set_quantity('level', [1,2,3,4], 'centroids')
     246        domain.set_quantity('xmomentum', [1,2,3,4], 'centroids')
     247        domain.set_quantity('ymomentum', [1,2,3,4], 'centroids')
     248                                         
    255249
    256250        #Assign some values to update vectors
    257         #domain.explicit_updates[]
     251        #Set explicit_update
     252
     253        for name in domain.conserved_quantities:
     254            domain.quantities[name].explicit_update = array([4.,3.,2.,1.])
     255            domain.quantities[name].semi_implicit_update = array([1.,1.,1.,1.])
     256           
     257
     258        #Update with given timestep (assuming no other forcing terms)
     259        domain.timestep = 0.1
     260        domain.update_conserved_quantities()
     261
     262        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )
     263        x /= array( [.9,.9,.9,.9] )
     264
     265        for name in domain.conserved_quantities:
     266            assert allclose(domain.quantities[name].centroid_values, x)       
     267
    258268       
    259269       
  • inundation/ga/storm_surge/pyvolution/test_shallow_water.py

    r212 r215  
    534534            pass
    535535
    536 
    537 
     536    #####################################################   
     537    def test_gravity(self):
     538        #Assuming no friction
     539
     540        from config import g
     541       
     542        a = [0.0, 0.0]
     543        b = [0.0, 2.0]
     544        c = [2.0, 0.0]
     545        d = [0.0, 4.0]
     546        e = [2.0, 2.0]
     547        f = [4.0, 0.0]
     548
     549        points = [a, b, c, d, e, f]
     550        #bac, bce, ecf, dbe
     551        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     552       
     553        domain = Domain(points, vertices)
     554
     555        #Set up for a gradient of (3,0) at mid triangle         
     556        def slope(x, y):
     557            return 3*x
     558
     559        h = 0.1
     560        def level(x,y):
     561            return slope(x,y)+h
     562
     563        domain.set_quantity('elevation', slope)
     564        domain.set_quantity('level', level)
     565
     566        for name in domain.conserved_quantities:
     567            assert allclose(domain.quantities[name].explicit_update, 0)
     568            assert allclose(domain.quantities[name].semi_implicit_update, 0)
     569           
     570        domain.compute_forcing_terms()                                   
     571
     572        assert allclose(domain.quantities['level'].explicit_update, 0)
     573        assert allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3)
     574        assert allclose(domain.quantities['ymomentum'].explicit_update, 0)
     575       
     576       
     577
     578    #####################################################
    538579    def test_first_order_extrapolator_const_z(self):
     580       
     581        a = [0.0, 0.0]
     582        b = [0.0, 2.0]
     583        c = [2.0, 0.0]
     584        d = [0.0, 4.0]
     585        e = [2.0, 2.0]
     586        f = [4.0, 0.0]
     587
     588        points = [a, b, c, d, e, f]
     589        #bac, bce, ecf, dbe
     590        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     591       
     592        domain = Domain(points, vertices)
     593        val0 = 2.+2.0/3
     594        val1 = 4.+4.0/3
     595        val2 = 8.+2.0/3
     596        val3 = 2.+8.0/3
     597
     598        zl=zr=-3.75 #Assume constant bed (must be less than level)       
     599        domain.set_quantity('elevation', zl*ones( (4,3) ))
     600        domain.set_quantity('level', [[val0, val0-1, val0-2],
     601                                      [val1, val1+1, val1],
     602                                      [val2, val2-2, val2],
     603                                      [val3-0.5, val3, val3]])
     604
     605
     606
     607        domain.order = 1
     608        domain.distribute_to_vertices_and_edges()
     609
     610        #Check that centroid values were distributed to vertices
     611        C = domain.quantities['level'].centroid_values
     612        for i in range(3):
     613            assert allclose( domain.quantities['level'].vertex_values[:,i], C)
     614       
     615
     616    def test_first_order_limiter_variable_z(self):
     617        #Check that first order limiter follows bed_slope
     618        from Numeric import alltrue, greater_equal
     619        from config import epsilon
    539620       
    540621        a = [0.0, 0.0]
     
    555636        val3 = 2.+8.0/3
    556637
    557         zl=zr=-3.75 #Assume constant bed (must be less than level)       
    558         domain.set_quantity('elevation', zl*ones( (4,3) ))
    559         domain.set_quantity('level', [[val0, val0-1, val0-2],
    560                                       [val1, val1+1, val1],
    561                                       [val2, val2-2, val2],
    562                                       [val3-0.5, val3, val3]])
    563 
    564 
    565 
    566         domain.order = 1
    567         domain.distribute_to_vertices_and_edges()
    568 
    569         #Check that centroid values were distributed to vertices
    570         C = domain.quantities['level'].centroid_values
    571         for i in range(3):
    572             assert allclose( domain.quantities['level'].vertex_values[:,i], C)
    573        
    574 
    575     def test_first_order_limiter_variable_z(self):
    576         #Check that first order limiter follows bed_slope
    577         from Numeric import alltrue, greater_equal
    578         from config import epsilon
    579        
    580         a = [0.0, 0.0]
    581         b = [0.0, 2.0]
    582         c = [2.0,0.0]
    583         d = [0.0, 4.0]
    584         e = [2.0, 2.0]
    585         f = [4.0,0.0]
    586 
    587         points = [a, b, c, d, e, f]
    588         #bac, bce, ecf, dbe
    589         vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
    590        
    591         domain = Domain(points, vertices)
    592         val0 = 2.+2.0/3
    593         val1 = 4.+4.0/3
    594         val2 = 8.+2.0/3
    595         val3 = 2.+8.0/3
    596 
    597638        domain.set_quantity('elevation', [[0,0,0], [6,0,0],
    598639                                          [6,6,6], [6,6,6]])
     
    616657        assert alltrue(alltrue(greater_equal(L,E-epsilon)))
    617658
     659
     660    #####################################################   
     661    def test_distribute_near_bed(self):
     662        #Assuming no friction
     663
     664        from config import g
     665       
     666        a = [0.0, 0.0]
     667        b = [0.0, 2.0]
     668        c = [2.0, 0.0]
     669        d = [0.0, 4.0]
     670        e = [2.0, 2.0]
     671        f = [4.0, 0.0]
     672
     673        points = [a, b, c, d, e, f]
     674        #bac, bce, ecf, dbe
     675        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     676       
     677        domain = Domain(points, vertices)
     678        domain.order = 2
     679
     680        #Set up for a gradient of (10,0) at mid triangle         
     681        def slope(x, y):
     682            return 10*x
     683
     684        h = 0.1
     685        def level(x,y):
     686            return slope(x,y)+h
     687
     688        domain.set_quantity('elevation', slope)
     689        domain.set_quantity('level', level, 'centroids')
     690
     691        E = domain.quantities['elevation'].vertex_values
     692        L = domain.quantities['level'].vertex_values
     693
     694        #print E
     695        domain.distribute_to_vertices_and_edges()
     696        #print L
     697
     698        #FIXME: HERTIL 24/8/4. Write similar test using pyvolution2 and
     699        #get test data from that
    618700
    619701       
     
    850932        domain.default_order=2
    851933
    852         #Bed-slope and friction
     934        #Bed-slope and friction at vertices (and interpolated elsewhere)
    853935        def x_slope(x, y):
    854936            return -x/3
Note: See TracChangeset for help on using the changeset viewer.