Changeset 234


Ignore:
Timestamp:
Aug 27, 2004, 4:33:24 PM (21 years ago)
Author:
ole
Message:

Played with protection against mnimal heights and
balanced limiting. Looks good.
The demo show_balanced_limiters now runs!

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

Legend:

Unmodified
Added
Removed
  • inundation/ga/storm_surge/pyvolution/HUSK.txt

    r232 r234  
     1
     2Perhaps control old/new style with minimum_allowed height
     3       
    14Test friction term!
    25
  • inundation/ga/storm_surge/pyvolution/config.py

    r221 r234  
    55
    66default_boundary_tag = 'exterior'
     7
     8newstyle = False #Flase meqans that we are equivalent to previous version pyvolution2
    79
    810
  • inundation/ga/storm_surge/pyvolution/domain.py

    r232 r234  
    7070
    7171        #Defaults
    72         from config import max_smallsteps, beta
     72        from config import max_smallsteps, beta, newstyle
    7373        self.beta = beta
     74        self.newstyle = newstyle
    7475        self.default_order = 1
    7576        self.order = self.default_order
  • inundation/ga/storm_surge/pyvolution/quantity.py

    r229 r234  
    183183                x2, y2 = self.domain.centroids[k2] #V2 centroid
    184184
    185                 #if k < 4:
    186                 #    print 'k, k0, k1, k2 = %d, %d, %d, %d' %(k, k0, k1, k2)
    187                 #    #print 'k = %d: (%f, %f), (%f, %f), (%f, %f)'\
    188                 #    #      %(k, x0, y0, x1, y1, x2, y2)
    189                 #    print 'k = %d: (%f, %f, %f)' %(k, q0, q1, q2)   
    190                
    191185                #Gradient
    192186                a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
     
    220214            dq[:,i] = qv[:,i] - qc
    221215       
    222         #Find min and max at this and neighbour's centroids
     216        #Find min and max of this and neighbour's centroid values
    223217        qmax = zeros(qc.shape, Float)
    224218        qmin = zeros(qc.shape, Float)       
     
    275269       
    276270        a, b = self.compute_gradients()
    277         ##for k in range(4):
    278         ##    print 'Gradients %d: %16.12f, %16.12f' %(k, a[k], b[k])
    279271
    280272        V = self.domain.get_vertex_coordinates()
  • inundation/ga/storm_surge/pyvolution/shallow_water.py

    r232 r234  
    272272def distribute_to_vertices_and_edges(domain):
    273273
    274     ##print 'Calling distrib within SW'
     274    protect_against_negative_heights_centroid(domain)
    275275    if domain.order == 1:
    276         #FIXME: This should be cleaned up, but we try to follow
    277         #pyvolution 2 strictly for now
    278         protect_against_negative_heights_centroid(domain)       
    279         protect_against_infinitesimal_heights_centroid(domain)       
    280276        extrapolate_first_order(domain)
    281277    elif domain.order == 2:
    282         #protect_against_negative_heights_centroid(domain)       
    283         extrapolate_second_order(domain)       
     278        extrapolate_second_order(domain)
    284279    else:
    285280        raise 'Unknown order'
    286 
     281   
    287282    #Compute edge values
    288283    for name in domain.conserved_quantities:
     
    290285        Q.interpolate_from_vertices_to_edges()           
    291286
    292 
    293 def protect_against_infinitesimal_heights_centroid(domain):
    294     """Adjust height and momentum at centroid if height is less than
     287def protect_against_infinitesimal_heights(domain):
     288    """Protect against infinitesimal heights and high speeds   
     289    Adjust height and momentum at centroid if height is less than
    295290    minimum allowed height
    296291    """
    297     #FIXME: Used only in first order
    298    
    299     #Water levels at centroids
     292
     293    #FIXME: May be unnecessary
     294    #Shortcuts
     295    wv = domain.quantities['level'].vertex_values
     296    zv = domain.quantities['elevation'].vertex_values
    300297    wc = domain.quantities['level'].centroid_values
    301 
    302     #Bed elevations at centroids
    303298    zc = domain.quantities['elevation'].centroid_values
    304 
    305     #Water depths at centroids   
     299    hv = wv-zv
    306300    hc = wc - zc
    307301
    308     #Momentums at centroids
     302    for k in range(domain.number_of_elements):
     303        hmax = max(hv[k,:])
     304       
     305        if hmax < domain.minimum_allowed_height:
     306           
     307            #Reset negative heights to bed elevation       
     308            if hc[k] < 0.0:
     309                print 'C'
     310                wc[k] = zc[k]
     311            for i in range(3):   
     312                if hv[k,i] < 0.0:
     313                    print 'V'                   
     314                    wv[k,i] = zv[k,i]
     315
     316
     317def protect_against_negative_heights_centroid(domain):
     318    """Protect against infinitesimal heights and associated high velocities
     319    """
     320
     321    #FIXME: Change name appropriately
     322   
     323    #Shortcuts
     324    wc = domain.quantities['level'].centroid_values
     325    zc = domain.quantities['elevation'].centroid_values
    309326    xmomc = domain.quantities['xmomentum'].centroid_values
    310     ymomc = domain.quantities['ymomentum'].centroid_values       
    311 
     327    ymomc = domain.quantities['ymomentum'].centroid_values           
     328    hc = wc - zc  #Water depths at centroids   
     329
     330    #Update
    312331    for k in range(domain.number_of_elements):
    313         #Protect against infinitesimal heights and high velocities
    314         if hc[k] < domain.minimum_allowed_height:
    315             #Control level and height
     332
     333
     334        if domain.newstyle:
     335            if hc[k] < domain.minimum_allowed_height:       
     336                if hc[k] < 0.0:
     337                    #Control level and height
     338                    wc[k] = zc[k]; hc[k] = 0.0
     339
     340                #Control momentum
     341                xmomc[k] = ymomc[k] = 0.0
     342        else:       
    316343            if hc[k] < 0.0:
     344                #Control level and height
    317345                wc[k] = zc[k]; hc[k] = 0.0
    318346
    319             #Control momentum
    320             xmomc[k] = ymomc[k] = 0.0
    321                
    322 
    323 def protect_against_negative_heights_centroid(domain):
    324     """Adjust height and momentum at centroid if height is less than zero
    325     """
    326 
    327    
    328     #Water levels at centroids
    329     wc = domain.quantities['level'].centroid_values
    330 
    331     #Bed elevations at centroids
    332     zc = domain.quantities['elevation'].centroid_values
    333 
    334     #Water depths at centroids   
    335     hc = wc - zc
    336 
    337     #Momentums at centroids
    338     xmomc = domain.quantities['xmomentum'].centroid_values
    339     ymomc = domain.quantities['ymomentum'].centroid_values       
    340 
    341     for k in range(domain.number_of_elements):
    342         #Protect against infinitesimal heights and high velocities
    343         if hc[k] < 0.0:
    344             #Control level and height
    345             wc[k] = zc[k]; hc[k] = 0.0
    346 
    347             #Control momentum
    348             xmomc[k] = ymomc[k] = 0.0
    349                
    350 
     347                #Control momentum
     348                xmomc[k] = ymomc[k] = 0.0
    351349
    352350
     
    367365    """
    368366
     367    #Shortcuts
     368    wc = domain.quantities['level'].centroid_values
     369    zc = domain.quantities['elevation'].centroid_values
     370    xmomc = domain.quantities['xmomentum'].centroid_values
     371    ymomc = domain.quantities['ymomentum'].centroid_values           
     372    hc = wc - zc   #Water depths at centroids
     373
     374
     375    if not domain.newstyle:
     376        #Protect against infinitesimal heights and high speeds at_centroid
     377        #FIXME: Try without when using the second balance function
     378       
     379        for k in range(domain.number_of_elements):
     380            #Protect against infinitesimal heights and high velocities
     381            if hc[k] < domain.minimum_allowed_height:
     382                #Control level and height
     383                if hc[k] < 0.0:
     384                    wc[k] = zc[k]; hc[k] = 0.0#
     385         
     386                #Control momentum
     387                xmomc[k] = ymomc[k] = 0.0
     388               
    369389   
    370390    #Update conserved quantities using straight first order
     
    373393        Q.extrapolate_first_order()
    374394
    375      
    376  
    377     #Water levels at centroids
    378     wc = domain.quantities['level'].centroid_values
    379 
    380     #Bed elevations at centroids
    381     zc = domain.quantities['elevation'].centroid_values
    382 
    383     #Water depths at centroids   
    384     hc = wc - zc
    385 
    386     #Water levels at vertices
    387     wv = domain.quantities['level'].vertex_values
    388 
    389     #Bed elevations at vertices
    390     zv = domain.quantities['elevation'].vertex_values
    391                
    392     #Water depths at vertices
    393     hv = wv-zv
    394    
    395    
    396     #Computed weighted balance between constant levels and and
    397     #levels parallel to the bed elevation.     
    398     for k in range(domain.number_of_elements):               
    399                
    400         #Compute maximal variation in bed elevation
    401         z_range = max(abs(zv[k,0]-zc[k]),
    402                       abs(zv[k,1]-zc[k]),
    403                       abs(zv[k,2]-zc[k]))
    404 
    405 
    406         #Weighted balance between stage parallel to bed elevation
    407         #(wvi = zvi + hc) and constant stage (wvi = wc = zc+hc)
    408         #where i=0,1,2 denotes the vertex ids
    409         #
    410         #It follows that
    411         #  wvi = (1-alpha)*(zvi+hc) + alpha*(zc+hc) =
    412         #  (1-alpha)*zvi + alpha*zc + hc =
    413         #  zvi + hc + alpha*(zc-zvi)
    414         #
    415         #where alpha in [0,1] and defined as the ratio between hc and
    416         #the maximal difference from zc to zv0, zv1 and zv2
    417         #
    418         #Mathematically the following can be continued on using hc as
    419         #  wvi =
    420         #  zvi + hc + alpha*(zc+hc-zvi-hc) =
    421         #  zvi + hc + alpha*(hvi-hc)
    422         #since hvi = zc+hc-zvi in the constant case
    423 
    424        
    425         if z_range > 0.0:
    426             alpha = min(hc[k]/z_range, 1.0)
    427         else:
    428             alpha = 1.0
    429            
    430         #Update water levels
    431         for i in range(3):
    432             #FIXME: Use the original first-order one first, then switch
    433             wv[k,i] = zv[k,i] + hc[k] + alpha*(zc[k]-zv[k,i])
    434             #wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
    435 
    436         #FIXME: What about alpha weighting of momentum??
     395
     396
     397    if domain.newstyle:
     398        balance_deep_and_shallow(domain)                      #This will be better
     399    else:   
     400        old_first_order_balancing_of_deep_and_shallow(domain) #Tests will pass
     401
     402   
    437403           
    438404       
     
    460426    """
    461427   
    462     #FIXME: first and second order might merge
    463 
    464     #Update conserved quantities using straight second order
     428
     429    #Update conserved quantities using straight second order with limiter
    465430    for name in domain.conserved_quantities:
    466431        Q = domain.quantities[name]
    467432        Q.extrapolate_second_order()
    468        
    469     #print 'y1', Q.vertex_values[1,:]   #OK
    470    
    471     #FIXME - like pyvolution 2 .......................
    472     protect_against_negative_heights_centroid(domain)   
    473 
    474     #print 'y1', Q.vertex_values[1,:]   #OK   
    475    
    476     for name in domain.conserved_quantities:
    477         Q = domain.quantities[name]     
    478         Q.limit()
    479 
    480  
    481     #print 'y1', Q.vertex_values[1,:]   #OK   
    482      
    483     #Water levels at centroids
     433        Q.limit()       
     434
     435
     436    balance_deep_and_shallow(domain)
     437    ##protect_against_infinitesimal_heights(domain)  #FIXME: Probably not necessary
     438
     439
     440
     441def balance_deep_and_shallow(domain):
     442    #FIXME: Grap math comments from old_first_order_balance....(below)
     443   
     444    #Shortcuts
    484445    wc = domain.quantities['level'].centroid_values
    485 
    486     #Bed elevations at centroids
    487446    zc = domain.quantities['elevation'].centroid_values
    488 
    489     #Water depths at centroids   
    490447    hc = wc - zc
    491 
    492     #Water levels at vertices
     448   
    493449    wv = domain.quantities['level'].vertex_values
    494 
    495     #Bed elevations at vertices
    496450    zv = domain.quantities['elevation'].vertex_values
    497                
    498     #Water depths at vertices
    499451    hv = wv-zv
    500452
    501 
    502    
    503453    #Computed linear combination between constant levels and and
    504454    #levels parallel to the bed elevation.     
     
    530480            alpha = 1.0
    531481   
    532         ##if k==1: print 'alpha', alpha #OK
    533        
    534482        #Update water levels
    535483        #  (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) =
     
    551499            # xmomc and ymomc (shallow) and momentum
    552500            # from extrapolator xmomv and ymomv (deep).
    553        
    554501            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:];           
    555502            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:];                     
    556503           
    557     #print 'y1', Q.vertex_values[1,:]   #OK   
    558            
    559     #Finally, protect against infinitesimal heights and high speeds
    560     #Water depths at vertices
    561     hv = wv-zv
    562     hc = wc-zc   
    563     for k in range(domain.number_of_elements):
    564         hmax = max(hv[k,:])
    565        
    566         if hmax < domain.minimum_allowed_height:       
    567             #Reset negative heights to bed elevation       
    568             if hc[k] < 0.0:
    569                 wc[k] = zc[k]
    570             for i in range(3):   
    571                 if hv[k,i] < 0.0:
    572                     wv[k,i] = zv[k,i]
    573 
    574504   
    575505
     
    674604
    675605            #Update momentum
    676             Xmom.explicit_update[k] += S*uh
    677             Ymom.explicit_update[k] += S*vh
     606            Xmom.semi_implicit_update[k] += S*uh
     607            Ymom.semi_implicit_update[k] += S*vh
    678608           
    679609
     
    847777            return self.W(x,y) + self.h
    848778
     779#STUFF
     780
     781
     782def old_first_order_balancing_of_deep_and_shallow(domain):
     783    #FIXME: Will be obsolete, but keep comments from this one.
     784   
     785    #Computed weighted balance between constant levels and and
     786    #levels parallel to the bed elevation.
     787    wc = domain.quantities['level'].centroid_values
     788    zc = domain.quantities['elevation'].centroid_values
     789    xmomc = domain.quantities['xmomentum'].centroid_values
     790    ymomc = domain.quantities['ymomentum'].centroid_values           
     791    hc = wc - zc   #Water depths at centroids
     792   
     793    wv = domain.quantities['level'].vertex_values
     794    zv = domain.quantities['elevation'].vertex_values
     795   
     796    for k in range(domain.number_of_elements):               
     797               
     798        #Compute maximal variation in bed elevation
     799        z_range = max(abs(zv[k,0]-zc[k]),
     800                      abs(zv[k,1]-zc[k]),
     801                      abs(zv[k,2]-zc[k]))
     802
     803
     804        #Weighted balance between stage parallel to bed elevation
     805        #(wvi = zvi + hc) and constant stage (wvi = wc = zc+hc)
     806        #where i=0,1,2 denotes the vertex ids
     807        #
     808        #It follows that
     809        #  wvi = (1-alpha)*(zvi+hc) + alpha*(zc+hc) =
     810        #  (1-alpha)*zvi + alpha*zc + hc =
     811        #  zvi + hc + alpha*(zc-zvi)
     812        #
     813        #where alpha in [0,1] and defined as the ratio between hc and
     814        #the maximal difference from zc to zv0, zv1 and zv2
     815        #
     816        #Mathematically the following could be continued on using hc as
     817        #  wvi =
     818        #  zvi + hc + alpha*(zc+hc-zvi-hc) =
     819        #  zvi + hc + alpha*(hvi-hc)
     820        #since hvi = zc+hc-zvi in the constant case
     821
     822       
     823        if z_range > 0.0:
     824            alpha = min(hc[k]/z_range, 1.0)
     825        else:
     826            alpha = 1.0
     827           
     828        #Update water levels
     829        for i in range(3):
     830            #FIXME: Use the original first-order one first, then switch
     831            wv[k,i] = zv[k,i] + hc[k] + alpha*(zc[k]-zv[k,i])
     832            #wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k]) #Requires updated hv!!
     833
     834        #FIXME: What about alpha weighting of momentum??
    849835
    850836
  • inundation/ga/storm_surge/pyvolution/show_balanced_limiters.py

    r233 r234  
    3636domain.visualise = False
    3737domain.visualise = True
    38 domain.default_order=2
     38domain.default_order = 2
     39domain.newstyle = True  #True is better
    3940
    4041
     
    4647print 'Field values'
    4748domain.set_quantity('elevation', Z)
    48 #domain.set_quantity('friction', manning)
     49domain.set_quantity('friction', manning)
    4950
    5051
     
    7778
    7879from Numeric import allclose
     80
    7981#Evolve
    80 for t in domain.evolve(yieldstep = 0.02, finaltime = 0.7):
    81 #for t in domain.evolve(yieldstep = 0.1, finaltime = 20):   
     82for t in domain.evolve(yieldstep = 0.1, finaltime = 20):
    8283    domain.write_time()
    83     print domain.quantities['level'].centroid_values[:6]
    84 
    85 if domain.time == 0.2:
    86     print 'Testing'
    87     assert allclose(domain.quantities['level'].centroid_values,
    88                     [ 3.96042007e-002,  5.61081828e-002,  4.66578380e-002,  5.73165329e-002,
    89                       4.72542001e-002,  5.74770060e-002,  4.74459150e-002,  5.77550805e-002,
    90                       4.80791695e-002,  5.85754074e-002,  4.90681598e-002,  6.02682664e-002,
    91                       1.16686827e-002,  1.75422685e-002,  1.17014731e-002,  2.15810992e-002,
    92                       1.30549421e-002,  2.14416681e-002,  1.31212934e-002,  2.15486277e-002,
    93                       1.34996488e-002,  2.24053139e-002,  1.50195194e-002,  2.22306851e-002,
    94                       -7.26762170e-003, -1.35071582e-003, -7.88143638e-003, -2.18165245e-003,
    95                       -7.81749271e-003, -1.06732807e-003, -7.76399231e-003, -1.00580353e-003,
    96                       -7.81877765e-003, -9.81086203e-004, -7.42500800e-003, -2.41412070e-004,
    97                       1.86244453e-001,  8.79324341e-002,  1.86232625e-001,  8.78313615e-002,
    98                       6.12537452e-002, -3.73125664e-002, -6.37550753e-002, -3.73269705e-002,
    99                       6.12145063e-002,  8.77700677e-002,  1.86257693e-001,  8.79121535e-002,
    100                       -4.83328632e-002,  1.18336527e-001, -4.83328400e-002,  1.18337005e-001,
    101                       -4.83328850e-002, -6.65776472e-003, -1.73331646e-001, -1.31654218e-001,
    102                       -1.73332232e-001, -6.66097985e-003, -4.83323869e-002,  1.18339536e-001,
    103                       -2.48333331e-001, -2.31666623e-001, -2.48333332e-001, -2.31666628e-001,
    104                       -2.48333332e-001, -2.31666627e-001, -2.48333330e-001, -2.31666575e-001,
    105                       -2.48333330e-001, -2.31666597e-001, -2.48333329e-001, -2.31666584e-001,
    106                       -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
    107                       -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
    108                       -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
    109                       -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
    110                       -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
    111                       -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
    112                       -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
    113                       -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
    114                       -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
    115                       -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
    116                       -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
    117                       -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
    118                       -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
    119                       -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
    120                       -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
    121                       -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
    122                       -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
    123                       -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
    124 
     84    #print domain.quantities['level'].centroid_values[:6]
    12585print 'Done'   
    12686   
Note: See TracChangeset for help on using the changeset viewer.