Changeset 3424


Ignore:
Timestamp:
Jul 26, 2006, 12:12:43 PM (18 years ago)
Author:
steve
Message:

minmod limiter for shallow water

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

Legend:

Unmodified
Added
Removed

  • development/pyvolution-1d/analytic_dam.py

    r3401 r3424  
    44
    55    def __init__(self, h0 = 5.0, h1 = 10.0, L = 2000.0):
    6 
     6        from math import sqrt
     7       
    78        self.h0 = h0 # depth upstream (m)
    89        self.h1 = h1 # depth downstream (m)
    910        self.L  = L  # length of domain
    10 
    11     def __call__(self, C,t):
    12        
    13         from Numeric import zeros,Float
    14         from math import sqrt, pi
    15    
    16         #t  = 0.0     # time (s)
    17         h0 = self.h0   
    18         h1 = self.h1   
    19         L = self.L
    20         n = len(C)    # number of cells
    2111
    2212        g  = 9.81    # gravity (m/s^2)
    … …  
    2414        c0 = sqrt(g*h0) #left celerity
    2515        c1 = sqrt(g*h1) #right celerity
    26 
    27 
    28         u = zeros(n,Float)
    29         h = zeros(n,Float)
    30         uh = zeros(n,Float)
    31         x = C-L/2
    32         #x = zeros(n,Float)
    33         #for i in range(n):
    34         #    x[i] = C[i]-1000.0
    35 
    36         # Upstream and downstream boundary conditions are set to the intial water
    37         # depth for all time.
    38 
    39         # Calculate Shock Speed
    40         #h2 = 7.2692044
    4116       
    42         #S2 = 2*h2/(h2-h0)*(sqrt(g*h1)-sqrt(g*h2))
    43         #u2 = S2 - g*h0/(4*S2)*(1+sqrt(1+8*S2*S2/(g*h0)))
    44 
    4517        zmin=-100.0
    4618        zmax=101.0
    … …  
    5729        if( abs(z) > 99.0):
    5830            print 'no convergence'
     31
     32        self.u2 = u2
     33        self.c0 = c0
     34        self.c1 = c1
     35        self.c2 = c2
     36        self.g = g
     37        self.z = z
     38
     39
     40
     41    def __call__(self, C,t):
     42       
     43        from Numeric import zeros,Float
     44        from math import sqrt
     45   
     46        #t  = 0.0     # time (s)
     47        h0 = self.h0   
     48        h1 = self.h1   
     49        L = self.L
     50        n = len(C)    # number of cells
     51
     52        u2 = self.u2
     53        c0 = self.c0
     54        c1 = self.c1
     55        c2 = self.c2
     56       
     57        g = self.g
     58        z = self.z
     59
     60        u = zeros(n,Float)
     61        h = zeros(n,Float)
     62        uh = zeros(n,Float)
     63        x = C-L/2.0
     64        #x = zeros(n,Float)
     65        #for i in range(n):
     66        #    x[i] = C[i]-1000.0
     67
     68        # Upstream and downstream boundary conditions are set to the intial water
     69        # depth for all time.
     70
     71        # Calculate Shock Speed
     72        #h2 = 7.2692044
     73       
     74        #S2 = 2*h2/(h2-h0)*(sqrt(g*h1)-sqrt(g*h2))
     75        #u2 = S2 - g*h0/(4*S2)*(1+sqrt(1+8*S2*S2/(g*h0)))
     76
     77
    5978
    6079        h2=h0/(1.0-u2/z)

  • development/pyvolution-1d/dam.py

    r3370 r3424  
    55from Numeric import allclose, array, zeros, ones, Float, take, sqrt
    66from config import g, epsilon
     7from analytic_dam import AnalyticDam
    78
    8 def analytical_sol(C,t):
     9h0 = 0.1
     10h1 = 10.0
     11
     12analytical_sol = AnalyticDam(h0, h1)
     13
     14## def analytical_sol(C,t):
    915   
    10     #t  = 0.0     # time (s)
    11     g  = 9.81    # gravity (m/s^2)
    12     h1 = 10.0    # depth upstream (m)
    13     h0 = 5.0     # depth downstream (m)
    14     L = 2000.0   # length of stream/domain (m)
    15     n = len(C)    # number of cells
     16##     #t  = 0.0     # time (s)
     17##     g  = 9.81    # gravity (m/s^2)
     18##     h1 = 10.0    # depth upstream (m)
     19##     h0 = 5.0     # depth downstream (m)
     20##     L = 2000.0   # length of stream/domain (m)
     21##     n = len(C)    # number of cells
    1622
    17     u = zeros(n,Float)
    18     h = zeros(n,Float)
    19     x = C-L/2
    20     #x = zeros(n,Float)
    21     #for i in range(n):
    22     #    x[i] = C[i]-1000.0
     23##     u = zeros(n,Float)
     24##     h = zeros(n,Float)
     25##     x = C-L/2
     26##     #x = zeros(n,Float)
     27##     #for i in range(n):
     28##     #    x[i] = C[i]-1000.0
    2329
    24     # Upstream and downstream boundary conditions are set to the intial water
    25     # depth for all time.
     30##     # Upstream and downstream boundary conditions are set to the intial water
     31##     # depth for all time.
    2632
    27     # Calculate Shock Speed
    28     h2 = 7.2692044
    29     S2 = 2*h2/(h2-h0)*(sqrt(g*h1)-sqrt(g*h2))
    30     u2 = S2 - g*h0/(4*S2)*(1+sqrt(1+8*S2*S2/(g*h0)))
     33##     # Calculate Shock Speed
     34##     h2 = 7.2692044
     35##     S2 = 2*h2/(h2-h0)*(sqrt(g*h1)-sqrt(g*h2))
     36##     u2 = S2 - g*h0/(4*S2)*(1+sqrt(1+8*S2*S2/(g*h0)))
    3137   
    32     #t=50
    33     #x = (-L/2:L/2)
    34     for i in range(n):
    35         # Calculate Analytical Solution at time t > 0
    36         u3 = 2/3*(sqrt(g*h1)+x[i]/t)
    37         h3 = 4/(9*g)*(sqrt(g*h1)-x[i]/(2*t))*(sqrt(g*h1)-x[i]/(2*t))
     38##     #t=50
     39##     #x = (-L/2:L/2)
     40##     for i in range(n):
     41##         # Calculate Analytical Solution at time t > 0
     42##         u3 = 2/3*(sqrt(g*h1)+x[i]/t)
     43##         h3 = 4/(9*g)*(sqrt(g*h1)-x[i]/(2*t))*(sqrt(g*h1)-x[i]/(2*t))
    3844
    39         if ( x[i] <= -t*sqrt(g*h1) ):
    40             u[i] = 0.0
    41             h[i] = h1
    42         elif ( x[i] <= t*(u2-sqrt(g*h2)) ):
    43             u[i] = u3
    44             h[i] = h3
    45         elif ( x[i] < t*S2 ):
    46             u[i] = u2
    47             h[i] = h2
    48         else:
    49             u[i] = 0.0
    50             h[i] = h0
     45##         if ( x[i] <= -t*sqrt(g*h1) ):
     46##             u[i] = 0.0
     47##             h[i] = h1
     48##         elif ( x[i] <= t*(u2-sqrt(g*h2)) ):
     49##             u[i] = u3
     50##             h[i] = h3
     51##         elif ( x[i] < t*S2 ):
     52##             u[i] = u2
     53##             h[i] = h2
     54##         else:
     55##             u[i] = 0.0
     56##             h[i] = h0
    5157
    52     return h , u*h
     58##     return h , u*h
    5359
    5460
    5561def newLinePlot(title='Simple Plot'):
    5662    import Gnuplot
    57     g = Gnuplot.Gnuplot(persist=1)
     63    g = Gnuplot.Gnuplot()
    5864    g.title(title)
    5965    g('set data style linespoints')
    … …  
    7581
    7682L = 2000.0     # Length of channel (m)
    77 N = 100        # Number of compuational cells
     83N = 400        # Number of compuational cells
    7884cell_len = L/N # Origin = 0.0
    7985
    … …  
    8894    for i in range(len(x)):
    8995        if x[i]<=1000.0:
    90             y[i] = 10.0
     96            y[i] = h1
    9197        else:
    92             y[i] = 5.0
     98            y[i] = h0
    9399    return y
    94100
    95101domain.set_quantity('stage', stage)
    96102
    97 domain.order = 2
     103
    98104domain.default_order = 2
    99105domain.cfl = 0.8
    100106#domain.beta = 0.0
    101 print "domain.order", domain.order
     107print "domain.order", domain.default_order
    102108
    103109if (debug == True):
    … …  
    117123plot1 = newLinePlot("Stage")
    118124plot2 = newLinePlot("Momentum")
     125plot3 = newLinePlot("Velocity")
    119126
    120127import time
    … …  
    122129yieldstep = 1.0
    123130finaltime = 50.0
     131
     132print domain.neighbour_vertices
     133
    124134for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
    125135    domain.write_time()
    … …  
    129139        linePlot(plot1,X,StageQ,X,y)
    130140        MomentumQ = domain.quantities['xmomentum'].vertex_values
     141
     142        Velocity = MomentumQ/StageQ
    131143        linePlot(plot2,X,MomentumQ,X,my)
     144        linePlot(plot3,X,Velocity,X,my/y)
    132145    #raw_input('press_return')
    133146    #pass
    … …  
    144157    print C
    145158
     159raw_input('press return')
    146160#f = file('test_solution_I.out', 'w')
    147161#for i in range(N):

  • development/pyvolution-1d/quantity.py

    r3370 r3424  
    414414        self.centroid_values += timestep*self.explicit_update
    415415
    416         #Semi implicit updates
    417         denominator = ones(N, Float)-timestep*self.semi_implicit_update
    418 
    419         if sum(equal(denominator, 0.0)) > 0.0:
    420             msg = 'Zero division in semi implicit update. Call Stephen :-)'
    421             raise msg
    422         else:
    423             #Update conserved_quantities from semi implicit updates
    424             self.centroid_values /= denominator
     416##         #Semi implicit updates
     417##         denominator = ones(N, Float)-timestep*self.semi_implicit_update
     418
     419##         if sum(equal(denominator, 0.0)) > 0.0:
     420##             msg = 'Zero division in semi implicit update. Call Stephen :-)'
     421##             raise msg
     422##         else:
     423##             #Update conserved_quantities from semi implicit updates
     424##             self.centroid_values /= denominator
    425425
    426426
    … …  
    493493                G[k] = ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)
    494494
    495         return G
     495
     496    def compute_minmod_gradients(self):
     497        """Compute gradients of piecewise linear function defined by centroids of
     498        neighbouring volumes.
     499        """
     500
     501        from Numeric import array, zeros, Float,sign
     502       
     503        def xmin(a,b):
     504            return 0.5*(sign(a)+sign(b))*min(abs(a),abs(b))
     505
     506        def xmic(t,a,b):
     507            return xmin(t*xmin(a,b), 0.50*(a+b) )
     508
     509
     510
     511        N = self.centroid_values.shape[0]
     512
     513
     514        G = self.gradients
     515        Q = self.centroid_values
     516        X = self.domain.centroids
     517
     518        for k in range(N):
     519
     520            # first and last elements have boundaries
     521
     522            if k == 0:
     523
     524                #Get data
     525                k0 = k
     526                k1 = k+1
     527
     528                q0 = Q[k0]
     529                q1 = Q[k1]
     530
     531                x0 = X[k0] #V0 centroid
     532                x1 = X[k1] #V1 centroid
     533
     534                #Gradient
     535                G[k] = (q1 - q0)/(x1 - x0)
     536
     537            elif k == N-1:
     538
     539                #Get data
     540                k0 = k
     541                k1 = k-1
     542
     543                q0 = Q[k0]
     544                q1 = Q[k1]
     545
     546                x0 = X[k0] #V0 centroid
     547                x1 = X[k1] #V1 centroid
     548
     549                #Gradient
     550                G[k] = (q1 - q0)/(x1 - x0)
     551
     552            else:
     553                #Interior Volume (2 neighbours)
     554
     555                #Get data
     556                k0 = k-1
     557                k2 = k+1
     558
     559                q0 = Q[k0]
     560                q1 = Q[k]
     561                q2 = Q[k2]
     562
     563                x0 = X[k0] #V0 centroid
     564                x1 = X[k]  #V1 centroid (Self)
     565                x2 = X[k2] #V2 centroid
     566
     567                # assuming uniform grid
     568                d1 = (q1 - q0)/(x1-x0)
     569                d2 = (q2 - q1)/(x2-x1)
     570
     571                #Gradient
     572                #G[k] = (d1+d2)*0.5
     573                #G[k] = (d1*(x2-x1) - d2*(x0-x1))/(x2-x0)               
     574                G[k] = xmic( 0.9, d1, d2 )
     575       
    496576
    497577    def extrapolate_first_order(self):
    … …  
    517597        #print "gradients 1",Z
    518598       
    519         self.compute_gradients()
     599        self.compute_minmod_gradients()
    520600        #print "gradients 2", Z
    521601
    522602        G = self.gradients
    523603        V = self.domain.vertices
    524         Qc = self.centroid_values
    525         Qv = self.vertex_values
     604        qc = self.centroid_values
     605        qv = self.vertex_values       
    526606       
    527607        #Check each triangle
    … …  
    534614
    535615            #Extrapolate
    536             Qv[k,0] = Qc[k] + G[k]*(x0-x)
    537             Qv[k,1] = Qc[k] + G[k]*(x1-x)
     616            qv[k,0] = qc[k] + G[k]*(x0-x)
     617            qv[k,1] = qc[k] + G[k]*(x1-x)
    538618
    539619    """    def limit(self):

  • development/pyvolution-1d/shallow_water_1d.py

    r3370 r3424  
    1 u"""Class Domain -
     1"""Class Domain -
    221D interval domains for finite-volume computations of
    33the shallow water wave equation.
    … …  
    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
    … …  
    512512    domain.timestep = timestep
    513513
    514 #see comments in the corresponding method in shallow_water_ext.c
    515 def extrapolate_second_order_sw_c(domain):
    516     """Wrapper calling C version of extrapolate_second_order_sw
    517     """
    518     import sys
    519     from Numeric import zeros, Float
    520 
    521     N = domain.number_of_elements
    522 
    523     #Shortcuts
    524     Stage = domain.quantities['stage']
    525     Xmom = domain.quantities['xmomentum']
    526 #    Ymom = domain.quantities['ymomentum']
    527     from shallow_water_ext import extrapolate_second_order_sw
    528     extrapolate_second_order_sw(domain,domain.surrogate_neighbours,
    529 # cant find this in domain                      domain.number_of_boundaries,
    530                                 domain.centroid_coordinates,
    531                                 Stage.centroid_values,
    532                                 Xmom.centroid_values,
    533 #                                Ymom.centroid_values,
    534                                 domain.vertex_coordinates,
    535                                 Stage.vertex_values,
    536                                 Xmom.vertex_values)#,
    537 #                                Ymom.vertex_values)
    538 
    539 def compute_fluxes_c(domain):
    540     """Wrapper calling C version of compute fluxes
    541     """
    542 
    543     import sys
    544     from Numeric import zeros, Float
    545 
    546     N = domain.number_of_elements
    547 
    548     #Shortcuts
    549     Stage = domain.quantities['stage']
    550     Xmom = domain.quantities['xmomentum']
    551     #Ymom = domain.quantities['ymomentum']
    552     Bed = domain.quantities['elevation']
    553 
    554     timestep = float(sys.maxint)
    555     from shallow_water_ext import compute_fluxes
    556     domain.timestep = compute_fluxes(timestep, domain.epsilon, domain.g,
    557                                      domain.neighbours,
    558                                      #domain.neighbour_edges,
    559                                      domain.neighbour_vertices,
    560                                      domain.normals,
    561                                      #domain.edgelengths,
    562                                      #domain.radii,
    563                                      domain.areas,
    564                                      #Stage.edge_values,
    565                                      #Xmom.edge_values,
    566                                      #Ymom.edge_values,
    567                                      #Bed.edge_values,
    568                                      Stage.vertex_values,
    569                                      Xmom.vertex_values,
    570                                      Bed.vertex_values,
    571                                      Stage.boundary_values,
    572                                      Xmom.boundary_values,
    573                                      #Ymom.boundary_values,
    574                                      Stage.explicit_update,
    575                                      Xmom.explicit_update,
    576                                      #Ymom.explicit_update,
    577                                      domain.already_computed_flux)
     514## #see comments in the corresponding method in shallow_water_ext.c
     515## def extrapolate_second_order_sw_c(domain):
     516##     """Wrapper calling C version of extrapolate_second_order_sw
     517##     """
     518##     import sys
     519##     from Numeric import zeros, Float
     520
     521##     N = domain.number_of_elements
     522
     523##     #Shortcuts
     524##     Stage = domain.quantities['stage']
     525##     Xmom = domain.quantities['xmomentum']
     526## #    Ymom = domain.quantities['ymomentum']
     527##     from shallow_water_ext import extrapolate_second_order_sw
     528##     extrapolate_second_order_sw(domain,domain.surrogate_neighbours,
     529## # cant find this in domain                      domain.number_of_boundaries,
     530##                                 domain.centroid_coordinates,
     531##                                 Stage.centroid_values,
     532##                                 Xmom.centroid_values,
     533## #                                Ymom.centroid_values,
     534##                                 domain.vertex_coordinates,
     535##                                 Stage.vertex_values,
     536##                                 Xmom.vertex_values)#,
     537## #                                Ymom.vertex_values)
     538
     539## def compute_fluxes_c(domain):
     540##     """Wrapper calling C version of compute fluxes
     541##     """
     542
     543##     import sys
     544##     from Numeric import zeros, Float
     545
     546##     N = domain.number_of_elements
     547
     548##     #Shortcuts
     549##     Stage = domain.quantities['stage']
     550##     Xmom = domain.quantities['xmomentum']
     551##     #Ymom = domain.quantities['ymomentum']
     552##     Bed = domain.quantities['elevation']
     553
     554##     timestep = float(sys.maxint)
     555##     from shallow_water_ext import compute_fluxes
     556##     domain.timestep = compute_fluxes(timestep, domain.epsilon, domain.g,
     557##                                      domain.neighbours,
     558##                                      #domain.neighbour_edges,
     559##                                      domain.neighbour_vertices,
     560##                                      domain.normals,
     561##                                      #domain.edgelengths,
     562##                                      #domain.radii,
     563##                                      domain.areas,
     564##                                      #Stage.edge_values,
     565##                                      #Xmom.edge_values,
     566##                                      #Ymom.edge_values,
     567##                                      #Bed.edge_values,
     568##                                      Stage.vertex_values,
     569##                                      Xmom.vertex_values,
     570##                                      Bed.vertex_values,
     571##                                      Stage.boundary_values,
     572##                                      Xmom.boundary_values,
     573##                                      #Ymom.boundary_values,
     574##                                      Stage.explicit_update,
     575##                                      Xmom.explicit_update,
     576##                                      #Ymom.explicit_update,
     577##                                      domain.already_computed_flux)
    578578
    579579
    … …  
    608608
    609609    #Remove very thin layers of water
    610     protect_against_infinitesimal_and_negative_heights(domain)
     610    ##protect_against_infinitesimal_and_negative_heights(domain)
    611611
    612612    #Extrapolate all conserved quantities
    … …  
    632632        elif domain.order == 2:
    633633            Q.extrapolate_second_order()
    634             Q.limit()
     634            #Q.limit()
    635635        else:
    636636            raise 'Unknown order'
    637637
    638638    #Take bed elevation into account when water heights are small
    639     balance_deep_and_shallow(domain)
     639    #balance_deep_and_shallow(domain)
    640640
    641641    #Compute edge values by interpolation
Note: See TracChangeset for help on using the changeset viewer.