Changeset 205


Ignore:
Timestamp:
Aug 23, 2004, 2:45:42 PM (21 years ago)
Author:
ole
Message:

Finished and tested gradient and limiters
Added forcing terms (not tested yet)

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

Legend:

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

    r198 r205  
    5757
    5858use_extensions = True   #Try to use C-extensions
    59 use_extensions = False   #Do not use C-extensions
     59#use_extensions = False   #Do not use C-extensions
    6060
    6161use_psyco = True  #Use psyco optimisations
  • inundation/ga/storm_surge/pyvolution/domain.py

    r196 r205  
    3636        for name in self.conserved_quantities + other_quantities:
    3737            self.quantities[name] = Quantity(self)
     38
     39
     40        #Create an empty list for explicit forcing terms
     41        #
     42        # Explicit terms must have the form
     43        #
     44        #    G(q, t)
     45        #
     46        # and explicit scheme is
     47        #
     48        #   q^{(n+1}) = q^{(n)} + delta_t G(q^{n}, n delta_t)
     49        #
     50        #
     51        # FIXME: How to call and how function should look
     52
     53        self.explicit_forcing_terms = []
     54
     55
     56        #Create an empty list for semi implicit forcing terms if any
     57        #
     58        # Semi implicit forcing terms are assumed to have the form
     59        #
     60        #    G(q, t) = H(q, t) q
     61        #
     62        # and the semi implicit scheme will then be
     63        #
     64        #   q^{(n+1}) = q^{(n)} + delta_t H(q^{n}, n delta_t) q^{(n+1})
     65       
     66        self.semi_implicit_forcing_terms = []
     67
    3868
    3969        #Defaults
     
    6292        self.store=False
    6393        self.format = 'dat'   
    64         #self.smooth = False
    6594        self.smooth = True
    6695
     
    386415        they should be defined in Domain subclass
    387416        """
    388         pass
     417
     418        for f in self.explicit_forcing_terms:
     419           
     420            pass
     421
    389422
    390423
     
    395428
    396429        from Numeric import ones, sum, equal, Float
    397    
    398430           
    399431        N = self.number_of_elements
     
    402434        timestep = self.timestep
    403435
     436        #Reset semi_implicit_updates
     437        #(explicit updates are already set by compute_fluyzes)
     438        for i, name in enumerate(self.conserved_quantities):
     439            Q = self.quantities[name]
     440            Q.semi_implicit_update[:] = 0.
     441
     442        #Compute forcing terms   
    404443        self.compute_forcing_terms()   
    405444
     
    421460        for name in self.conserved_quantities:
    422461            Q = self.quantities[name]
    423             Q.distribute_to_vertices_and_edges()
    424 
     462            if self.order == 1:
     463                Q.first_order_extrapolator()
     464            elif self.order == 2:
     465                Q.second_order_extrapolator()
     466                Q.limiter()                               
     467            else:
     468                raise 'Unknown order'
     469            Q.interpolate_from_vertices_to_edges()                           
    425470
    426471
  • inundation/ga/storm_surge/pyvolution/quantity.py

    r198 r205  
    106106            self.centroid_values /= denominator
    107107
     108
    108109    def compute_gradients(self):
    109110        """Compute gradients of triangle surfaces defined by centroids of
     
    138139
    139140                k0 = k  #Self
    140                
    141141                #Find index of the one neighbour
    142142                for k1 in self.domain.neighbours[k,:]:
     
    144144                        break
    145145                assert k1 != k0
    146                 assert k1 >= 0                     
    147                    
     146                assert k1 >= 0
     147
    148148                #Get data
    149149                q0 = self.centroid_values[k0]
     
    159159                    b[k] = (x0*q1 - x1*q0)/det
    160160
     161
    161162            else:
    162163                #One or zero missing neighbours
     
    180181                #Gradient
    181182                a[k], b[k] = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
    182                                       #array([q0]), array([q1]), array([q2]))
    183183       
    184184        return a, b
    185185       
    186186
    187     def second_order_limiter(self):
    188         """Extrapolate conserved quantities from centroid to
    189         vertices for each volume using second order scheme.
    190 
     187    def limiter(self):
     188        """Limit slopes for each volume to eliminate artificial variance
     189        introduced by e.g. second order extrapolator
     190       
     191        This is an unsophisticated limiter as it does not take into
     192        account dependencies among quantities.
     193       
    191194        precondition:
    192195           vertex values are estimated from gradient
     
    242245
    243246       
    244     def first_order_limiter(self):
     247    def first_order_extrapolator(self):
    245248        """Extrapolate conserved quantities from centroid to
    246249        vertices for each volume using
     
    254257            qv[:,i] = qc
    255258           
    256 
    257     def distribute_to_vertices_and_edges(self, order=1):
     259           
     260    def second_order_extrapolator(self):
    258261        """Extrapolate conserved quantities from centroid to
    259         vertices and edge-midpoints for each volume using
    260         specified order.
    261 
    262         This is an unsophisticated limiter as it does not take into
    263         account dependencies among quantities.
    264         """
    265 
    266         if order == 1:
    267             self.first_order_limiter()
    268         elif order == 2:
    269             a, b = self.compute_gradients()
    270            
    271             self.second_order_limiter()
    272         else:
    273             msg = 'ERROR: unknown order %d' %order
    274             raise msg
    275 
    276         self.interpolate_from_vertices_to_edges()
    277 
     262        vertices for each volume using
     263        second order scheme.
     264        """
     265       
     266        a, b = self.compute_gradients()
     267
     268        V = self.domain.get_vertex_coordinates()
     269        qc = self.centroid_values
     270        qv = self.vertex_values
     271       
     272        #Check each triangle
     273        for k in range(self.domain.number_of_elements):
     274            #Centroid coordinates           
     275            x, y = self.domain.centroids[k]
     276
     277            #vertex coordinates
     278            x0, y0, x1, y1, x2, y2 = V[k,:]
     279
     280            #Extrapolate
     281            qv[k,0] = qc[k] + a[k]*(x0-x) + b[k]*(y0-y)
     282            qv[k,1] = qc[k] + a[k]*(x1-x) + b[k]*(y1-y)
     283            qv[k,2] = qc[k] + a[k]*(x2-x) + b[k]*(y2-y)           
     284
     285           
    278286
    279287    def interpolate_from_vertices_to_edges(self):
  • inundation/ga/storm_surge/pyvolution/shallow_water.py

    r198 r205  
    7171        compute_fluxes(self)
    7272
    73     def first_order_limiter(self):
     73    def distribute_to_vertices_and_edges(self):
    7474        #Call correct module function
    75         first_order_limiter(self)       
     75        #(either from this module or C-extension)       
     76        distribute_to_vertices_and_edges(self)
    7677       
    7778####################################
     
    259260        Xmom.explicit_update[k] = flux[1]
    260261        Ymom.explicit_update[k] = flux[2]
    261        
     262
    262263        timestep = min(timestep, optimal_timestep)
    263264
    264265    domain.timestep = timestep   
     266
    265267       
    266268
    267269####################################
    268 # Functions for gradient limiting (distribute_to_vertices_and_edges)
     270# Module functions for gradient limiting (distribute_to_vertices_and_edges)
     271
     272def distribute_to_vertices_and_edges(domain):
     273   
     274    protect_against_infinitesimal_heights_centroid(domain)
     275    if domain.order == 1:
     276        first_order_extrapolator(domain)
     277    elif domain.order == 2:
     278        second_order_extrapolator(domain)       
     279    else:
     280        raise 'Unknown order'
    269281
    270282def protect_against_infinitesimal_heights_centroid(domain):
     
    297309               
    298310
    299    
    300 
    301 def first_order_limiter(domain):
    302     """First order limiter function, specific to the shallow water wave
    303     equation.
     311
     312
     313
     314def first_order_extrapolator(domain):
     315    """First order extrapolator function, specific
     316    to the shallow water wave equation.
    304317   
    305318    It will ensure that h (w-z) is always non-negative even in the
     
    315328    """
    316329
    317     #FIXME: Might go elsewhere
    318     protect_against_infinitesimal_heights_centroid(domain)
    319330   
    320331    #Update conserved quantities using straight first order
    321332    for name in domain.conserved_quantities:
    322333        Q = domain.quantities[name]
    323         Q.first_order_limiter()
    324 
    325        
     334        Q.first_order_extrapolator()
     335        Q.interpolate_from_vertices_to_edges()       
     336     
     337 
    326338    #Water levels at centroids
    327339    wc = domain.quantities['level'].centroid_values
     
    377389            alpha = 1.0
    378390           
    379         #Update stage
     391        #Update water levels
    380392        for i in range(3):
    381393            #wv[k,i] = zv[k,i] + hc[k] + alpha*(zc[k]-zv[k,i])
     
    384396
    385397       
    386 def second_order_limiter(domain):
     398def second_order_extrapolator(domain):
    387399    """Second order limiter function, specific to the shallow water wave
    388400    equation.
     
    406418   
    407419    """
    408 
    409     #FIXME: Might go elsewhere
    410     protect_against_infinitesimal_heights_centroid(domain)
     420   
     421    #FIXME: first and second order migh merge
     422
     423    from Numeric import minimum, maximum
    411424   
    412425    #Update conserved quantities using straight second order
    413426    for name in domain.conserved_quantities:
    414427        Q = domain.quantities[name]
    415         Q.second_order_limiter()
    416 
    417 
     428
     429        Q.second_order_extrapolator()
     430        Q.limiter()
     431        Q.interpolate_from_vertices_to_edges()               
     432
     433 
    418434    #Water levels at centroids
    419435    wc = domain.quantities['level'].centroid_values
     
    425441    hc = wc - zc
    426442
    427     #Momentums at centroids
    428     xmomc = domain.quantities['xmomentum'].centroid_values
    429     ymomc = domain.quantities['ymomentum'].centroid_values       
    430 
    431443    #Water levels at vertices
    432     #wv = domain.quantities['level'].vertex_values
    433     #
     444    wv = domain.quantities['level'].vertex_values
     445
    434446    #Bed elevations at vertices
    435     #zv = domain.quantities['elevation'].vertex_values
    436     #
    437     #Momentums at vertices
    438     #xmomv = domain.quantities['xmomentum'].vertex_values
    439     #ymomv = domain.quantities['ymomentum'].vertex_values
    440    
     447    zv = domain.quantities['elevation'].vertex_values
     448               
     449    #Water depths at vertices
     450    hv = wv-zv
     451
     452
     453   
     454    #Computed linear combination between constant levels and and
     455    #levels parallel to the bed elevation.     
     456    for k in range(domain.number_of_elements):               
     457               
     458        #Compute maximal variation in bed elevation
     459        #  This quantitiy is
     460        #    dz = max_i abs(z_i - z_c)
     461        #  and it is independent of dimension
     462        #  In the 1d case zc = (z0+z1)/2
     463        #  In the 2d case zc = (z0+z1+z2)/3
     464        # 
     465       
     466        z_range = max(abs(zv[k,0]-zc[k]),
     467                      abs(zv[k,1]-zc[k]),
     468                      abs(zv[k,2]-zc[k]))
     469
     470
     471        hmin = minimum( hv[k, :] )
     472
     473        #Create alpha in [0,1], where alpha==0 means using shallow
     474        #first order scheme and alpha==1 means using the limited
     475        #second order scheme for the stage w
     476        #If hmin < 0 then alpha=0 reverrting to first order.
     477     
     478        if z_range > 0.0:
     479            alpha = max( min( 2*hmin/z_range, 1.0), 0.0)
     480        else:
     481            alpha = 1.0
     482
     483        #Update water levels
     484        #  (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) =
     485        #   zvi + hc + alpha*(hvi - hc)         
     486        for i in range(3):
     487            wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
     488           
     489
     490        #Momentums at centroids
     491        xmomc = domain.quantities['xmomentum'].centroid_values
     492        ymomc = domain.quantities['ymomentum'].centroid_values       
     493
     494        #Momentums at vertices
     495        xmomv = domain.quantities['xmomentum'].vertex_values
     496        ymomv = domain.quantities['ymomentum'].vertex_values               
     497
     498        # Update momentum as a linear combination of
     499        # xmomc and ymomc (shallow) and momentum
     500        # from extrapolator xmomv and ymomv (deep).
     501       
     502
     503        ##f##or i in range(3):
     504        #####    xmomv[k,i] = (1-alpha)*xmomc[k] + alpha*xmomv[k,i];
     505           
     506        xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:];           
     507           
     508
     509    #Finally, protect against infinitesimal heights and high speeds
     510    #Water depths at vertices
     511    hv = wv-zv
     512    hc = wc-zc   
    441513    for k in range(domain.number_of_elements):
    442         #Protect against infinitesimal heights and high velocities
    443         if hc[k] < domain.minimum_allowed_height:
    444             #Control level and height
     514        hmax = maximum(hv[k,:])
     515       
     516        if hmax < minimum_allowed_height:       
     517            #Reset negative heights to bed elevation       
    445518            if hc[k] < 0.0:
    446                 wc[k] = zc[k]; hc[k] = 0.0
    447 
    448             #Control momentum as well!
    449             xmomc[k] = ymomc[k] = 0.0
    450                
    451 
    452     #Compute second order limiter for each conserved quantity       
    453     for name in domain.conserved_quantities:
    454         domain.quantities[name].second_order_limiter()
    455 
    456  
    457 #   //Adjusted water heights at vertices   
    458 #   hv0 = qv0[0]-zv0;
    459 #   hv1 = qv1[0]-zv1;
    460 #   hv2 = qv2[0]-zv2;
    461  
    462 #   hmin = min( min(hv0, hv1), hv2 );
    463 
    464 #   // -----------------------------------
    465 #   //  Second run - Linear combination with
    466 #   //  first order scheme for shallow depths
    467 #   // -----------------------------------
    468  
    469  
    470 #   //Compute maximal variation in bed elevation
    471 #   //
    472 #   //This quantitiy is
    473 #   // dz = max_i abs(z_i - z_c)
    474 #   //and it is independent of dimension
    475 #   //In the 1d case zc = (z0+z1)/2
    476 #   //In the 2d case zc = (z0+z1+z2)/3
    477 #   //
    478 #   z_range = max( max(fabs(zv0-zc), fabs(zv1-zc)), fabs(zv2-zc) );
    479  
    480 #   //Create alpha in [0,1], where alpha==0 means using shallow
    481 #   //first order scheme and alpha==1 means using the limited
    482 #   //second order scheme for the stage w
    483 #   //If hmin < 0 then alpha=0 and first order method is
    484 #   if (z_range>0.0) {
    485 #     alpha = max( min(2*hmin/z_range, 1.0), 0.0);
    486 #   } else {
    487 #     alpha = 1.0;
    488 #   }
    489  
    490 #   if (alpha < 1) {
    491 #     //Update stage
    492 #     //
    493 #     //(1-alpha)*(zvi+hc) + alpha*(zvi+hvi) =
    494 #     //zvi + hc + alpha*(hvi - hc)
    495 #     //
    496 #     qv0[0] = zv0 + hc + alpha*(hv0-hc);
    497 #     qv1[0] = zv1 + hc + alpha*(hv1-hc);
    498 #     qv2[0] = zv2 + hc + alpha*(hv2-hc);
    499  
    500        
    501 #     //Update momentum as a linear combination of
    502 #     //qc[j] (shallow) and momentum from gradient limiter (deep).
    503 #     for (j=1; j<number_of_conserved_quantities; j++) {
    504 #       qv0[j] = (1-alpha)*qc[j] + alpha*qv0[j];
    505 #       qv1[j] = (1-alpha)*qc[j] + alpha*qv1[j];
    506 #       qv2[j] = (1-alpha)*qc[j] + alpha*qv2[j];             
    507 #     }
    508 #   }
    509 
    510 
    511 #   //Finally, protect against infinitesimal heights and high speeds
    512 #   hv0 = qv0[0]-zv0; hv1 = qv1[0]-zv1; hv2 = qv2[0]-zv2;
    513 #   hmax = max(hv0, max(hv1, hv2));
    514 #   if (hmax < minimum_allowed_height) {
    515  
    516 #     //Set negative heights to bed elevation       
    517 #     if (hc < 0.0)  {qc[0] = zc; hc = 0.0;} 
    518 #     if (hv0 < 0.0) {qv0[0] = zv0; hv0 = 0.0;}         
    519 #     if (hv1 < 0.0) {qv1[0] = zv1; hv1 = 0.0;}
    520 #     if (hv2 < 0.0) {qv2[0] = zv2; hv2 = 0.0;}
    521 
    522    
    523 #     //FIXME: Do we need this desperately??
    524 #     //FIXME: We definitely do encounter negative heights,
    525 #     //but what about momentum?
    526    
    527 #     //Control momentum
    528 #     //qc[1] = 0.0; qc[2] = 0.0;
    529 #     //qv0[1] = 0.0; qv0[2] = 0.0;   
    530 #     //qv1[1] = 0.0; qv1[2] = 0.0;   
    531 #     //qv2[1] = 0.0; qv2[2] = 0.0;           
    532 #   }
    533 # }
    534 
    535 
    536 
    537 
    538 def limiter_org(k):
    539     """limiter function, specific to the shallow water wave equation.
    540        It will ensure that h is always non-negative.
    541        This version works on the consecutive data structure
    542     """
    543 
    544     #This version should be used to reproduce the bad limiting when levels
    545     #are close to bed elevation
    546 
    547     import sys
    548     from config import beta
    549     from Numeric import ones, Float
    550 
    551     #conserved quantities at centroid
    552     qc = Volume.conserved_quantities_centroid[k, :]
    553     N = len(qc)
    554    
    555     zc = Volume.field_values_centroid[k,0]         #bed elevation at centroid
    556     wc = qc[0]                                     #Stage at centroid
    557     hc = wc - zc                                   #Water depth at centroid
    558 
    559     #conserved quant. at vertices
    560     qv0 = Volume.conserved_quantities_vertex0[k,:]
    561     qv1 = Volume.conserved_quantities_vertex1[k,:]
    562     qv2 = Volume.conserved_quantities_vertex2[k,:]
    563 
    564 
    565     #Bed elevation at vertices
    566     zv0 = Volume.field_values_vertex0[k,0]
    567     zv1 = Volume.field_values_vertex1[k,0]
    568     zv2 = Volume.field_values_vertex2[k,0]
    569 
    570     #deltas
    571     dq0 = qv0-qc
    572     dq1 = qv1-qc
    573     dq2 = qv2-qc   
    574 
    575            
    576     #-------------------------------------
    577     hmin = hc
    578     qmin = qmax = qc
    579     for i in range(3):
    580         n = Volume.neighbours[k,i]
    581         if n >= 0:
    582             qn = Volume.conserved_quantities_centroid[n,:]     
    583             zn = Volume.field_values_centroid[n,0]       
    584             hn = qn[0]-zn
    585 
    586             hmin = min(hmin, hn)
    587             qmin = minimum(qmin, qn)
    588             qmax = maximum(qmax, qn)
    589 
    590     hmin = max(hmin,0.0)
    591 
    592    
    593     #-----------------------------------
    594     # First run - limit on conserved quantities
    595     # based on neighbouring heights
    596     #-----------------------------------   
    597 
    598     dqmax = qmax - qc
    599     dqmin = qmin - qc
    600     phi = ones(len(qc), Float)
    601     for j in range(N):
    602         #First find the gradient limiter (phi)
    603        
    604         #Vertex0:
    605         r = 1.0               
    606         if dq0[j] > 0: r = dqmax[j]/dq0[j]
    607         if dq0[j] < 0: r = dqmin[j]/dq0[j]
    608         phi[j] = min( min(r*beta, 1), phi[j] )
    609 
    610         #Vertex1:
    611         r = 1.0               
    612         if dq1[j] > 0: r = dqmax[j]/dq1[j]
    613         if dq1[j] < 0: r = dqmin[j]/dq1[j]
    614         phi[j] = min( min(r*beta, 1), phi[j] )
    615 
    616         #Vertex2:
    617         r = 1.0               
    618         if dq2[j] > 0: r = dqmax[j]/dq2[j]
    619         if dq2[j] < 0: r = dqmin[j]/dq2[j]
    620         phi[j] = min( min(r*beta, 1), phi[j] )
    621 
    622         #The update using phi limiter
    623         qv0[j] = qc[j] + phi[j]*dq0[j] #Vertex0
    624         qv1[j] = qc[j] + phi[j]*dq1[j] #Vertex1
    625         qv2[j] = qc[j] + phi[j]*dq2[j] #Vertex2
    626        
    627 
    628        
    629     #-----------------------------------
    630     # Second run - lower limit on height
    631     # based on neighbouring heights for
    632     # surfaces constructed below bed
    633     #-----------------------------------
    634 
    635     #Adjusted water heights at vertices   
    636     hv0 = qv0[0]-zv0
    637     hv1 = qv1[0]-zv1
    638     hv2 = qv2[0]-zv2
    639 
    640     #Deltas
    641     dhmin = hmin - hc
    642     dh0 = hv0 - hc
    643     dh1 = hv1 - hc
    644     dh2 = hv2 - hc
    645 
    646 
    647     #Work out if update is needed and calculate phi limiter
    648     phi = 1.0
    649     update_needed = False   
    650     if hv0 < 0.0:
    651         update_needed = True
    652         r = 1.0
    653         if dh0 < 0: r = dhmin/dh0
    654         phi = min( min(beta*r, 1), phi )
    655        
    656     if hv1 < 0.0:
    657         update_needed = True
    658         r = 1.0
    659         if dh1 < 0: r = dhmin/dh1
    660         phi = min( min(beta*r, 1), phi )
    661 
    662     if hv2 < 0.0:
    663         update_needed = True
    664         r = 1.0
    665         if dh2 < 0: r = dhmin/dh2
    666         phi = min( min(beta*r, 1), phi )
    667 
    668 
    669     #Update stage at vertices   
    670     if update_needed:
    671         hv0 = hc + phi*dh0
    672         if hv0 < 0.0: raise 'hv0 = %12f hc = %12f ' % (hv0 , hc)
    673         qv0[0] = zv0 + hv0
    674 
    675         hv1 = hc + phi*dh1
    676         if hv1 < 0.0: raise 'hv1 = %12f hc = %12f ' % (hv1 , hc)
    677         qv1[0] = zv1 + hv1       
    678 
    679         hv2 = hc + phi*dh2
    680         if hv2 < 0.0: raise 'hv2 = %12f hc = %12f ' % (hv2 , hc)
    681         qv2[0] = zv2 + hv2
    682        
    683        
    684     #In either case protect against infinitesimal heights
    685     dry(k)
    686 
    687 def dry(k):
    688     """Protect agains infinitesimal heights (and high speeds)
    689     There are no states for any conserved quantity when h is
    690     between minimum allowed height and 0.
    691     """
    692 
    693     #volume = Volume.instances[k]
    694    
    695     maxhv = 0.0
    696 
    697     wv0 = Volume.conserved_quantities_vertex0[k,0]       
    698     zv0 = Volume.field_values_vertex0[k,0]
    699     hv0 = wv0-zv0
    700 
    701     wv1 = Volume.conserved_quantities_vertex1[k,0]       
    702     zv1 = Volume.field_values_vertex1[k,0]
    703     hv1 = wv1-zv1
    704 
    705     wv2 = Volume.conserved_quantities_vertex2[k,0]       
    706     zv2 = Volume.field_values_vertex2[k,0]
    707     hv2 = wv2-zv2   
    708        
    709     maxhv = max(hv0, hv1, hv2)
    710 
    711     if maxhv < minimum_allowed_height:
    712         #Set water level to bed elevation and zero momentums.
    713         Volume.conserved_quantities_vertex0[k,:] =\
    714                                                  array([Volume.field_values_vertex0[k,0], 0, 0])
    715 
    716         Volume.conserved_quantities_vertex1[k,:] =\
    717                                                  array([Volume.field_values_vertex1[k,0], 0, 0])
    718 
    719         Volume.conserved_quantities_vertex2[k,:] =\
    720                                                  array([Volume.field_values_vertex2[k,0], 0, 0])
    721 
    722         #Set only momentums to zero
    723         #Volume.conserved_quantities_vertex0[k,1:] = array([0., 0.])
    724         #Volume.conserved_quantities_vertex1[k,1:] = array([0., 0.])
    725         #Volume.conserved_quantities_vertex2[k,1:] = array([0., 0.])
    726 
    727 
    728 
    729 
    730 def distribute_to_vertices_and_edges(domain):
    731     """Extrapolate conserved quantities from centroid to
    732     vertices and edge-midpoints for each volume
    733     """
    734 
    735     from config import minimum_allowed_height
    736 
    737     from mesh import Point       
    738     order = domain.order
    739        
    740     limiter = domain.limiter
    741     if limiter is None:
    742         from shallow_water import limit_on_w_then_h as limiter
    743 
    744 
    745     for k in range(Volume.number_of_instances):
    746        
    747         #conserved quantities at centroid       
    748         q = Volume.conserved_quantities_centroid[k,:]
    749 
    750         if order == 1:
    751             ########################
    752             #NOTE: This is now specific to the shallow water wave equation
    753 
    754             zc = Volume.field_values_centroid[k,0] #bed elevation at centroid
    755             wc = q[0]                              #Stage at centroid
    756             hc = wc - zc                           #Water depth at centroid
    757 
    758             #Protect against infinitesimal heights and high velocities
    759             if hc < minimum_allowed_height:
    760                 #Set momentums to zero
    761                 q[0] = zc
    762                 q[1] = 0.0
    763                 q[2] = 0.0     
    764            
    765 
    766             #Bed elevation at vertices
    767             zv0 = Volume.field_values_vertex0[k,0]
    768             zv1 = Volume.field_values_vertex1[k,0]
    769             zv2 = Volume.field_values_vertex2[k,0]
    770            
    771             #Update conserved quant. at vertices
    772             Volume.conserved_quantities_vertex0[k,1:] = q[1:] #Momentum
    773             Volume.conserved_quantities_vertex1[k,1:] = q[1:] #Momentum
    774             Volume.conserved_quantities_vertex2[k,1:] = q[1:] #Momentum
    775 
    776             #Stage
    777             z_range = max( max(abs(zv0-zv1), abs(zv0-zv2)), abs(zv1-zv2) )
    778             if z_range>0:
    779                 alpha = min(hc/z_range, 1.0)
    780             else:
    781                 alpha = 1
    782                
    783             zz = hc+alpha*zc
    784             Volume.conserved_quantities_vertex0[k,0] = zz + (1-alpha)*zv0
    785             Volume.conserved_quantities_vertex1[k,0] = zz + (1-alpha)*zv1
    786             Volume.conserved_quantities_vertex2[k,0] = zz + (1-alpha)*zv2
    787 
    788         elif order == 2:
    789             a, b = compute_gradient(k)           
    790 
    791             #Compute extrapolated values at vertices
    792             x,y = Point.coordinates[Volume.points[k,3]] #Centroid
    793 
    794             #Get vertex coordinates
    795             i = Volume.points[k,0]
    796             x0 = Point.coordinates[i,0]
    797             y0 = Point.coordinates[i,1]
    798 
    799             i = Volume.points[k,1]
    800             x1 = Point.coordinates[i,0]
    801             y1 = Point.coordinates[i,1]
    802 
    803             i = Volume.points[k,2]             
    804             x2 = Point.coordinates[i,0]
    805             y2 = Point.coordinates[i,1]                       
    806 
    807 
    808             #Extrapolate
    809             Volume.conserved_quantities_vertex0[k, :] = q + a*(x0-x) + b*(y0-y)
    810             Volume.conserved_quantities_vertex1[k, :] = q + a*(x1-x) + b*(y1-y)
    811             Volume.conserved_quantities_vertex2[k, :] = q + a*(x2-x) + b*(y2-y)
    812            
    813 
    814             if limiter is not None:
    815                 limiter(k)
    816         else:
    817             raise 'Unknown order'
    818 
    819 
    820         #Compute conserved quantities as interpolated
    821         #quantities at each edge midpoint
    822         q0 = Volume.conserved_quantities_vertex0[k,:]
    823         q1 = Volume.conserved_quantities_vertex1[k,:]
    824         q2 = Volume.conserved_quantities_vertex2[k,:]
    825 
    826         Volume.conserved_quantities_face0[k,:] = 0.5*(q1 + q2)
    827         Volume.conserved_quantities_face1[k,:] = 0.5*(q0 + q2)
    828         Volume.conserved_quantities_face2[k,:] = 0.5*(q0 + q1)       
     519                wc[k] = zc[k]
     520                ###hc[k] = 0.0
     521            for i in range(3):   
     522                if hv[k,i] < 0.0:
     523                    wv[k,i] = zv[k,i]
     524                    ##hv0 = 0.0;}         
     525
     526   
    829527
    830528
     
    869567
    870568
     569#########################
     570#Standard forcing terms:
     571
     572def gravity(domain):
     573    """Implement forcing function for bed slope working with
     574    consecutive data structures of class Volume
     575    """
     576
     577    from config import g
     578    from util import gradient
     579    from Numeric import zeros, Float, array, sum
     580
     581    Level = domain.quantities['level']
     582    Xmom = domain.quantities['xmomentum']
     583    Ymom = domain.quantities['ymomentum']
     584
     585    Elevation = domain.quantities['elevation']   
     586    h = Level.edge_values - Elevation.edge_values
     587    V = domain.get_vertex_coordinates()
     588   
     589    for k in range(domain.number_of_elements):   
     590        avg_h = sum( h[k,:] )/3
     591
     592        #Compute bed slope
     593        x0, y0, x1, y1, x2, y2 = V[k,:]   
     594        z0, z1, z2 = Elevation.vertex_values[k,:]
     595        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
     596
     597        #Update momentum
     598        Xmom.explicit_update[k] += -g*zx*avg_h
     599        Ymom.explicit_update[k] += -g*zy*avg_h       
     600   
     601
     602def manning_friction(domain):
     603    """Apply (Manning) friction to water momentum
     604    """
     605
     606    import Numeric
     607    from math import sqrt
     608    from config import g, minimum_allowed_height
     609
     610    Level = domain.quantities['level']
     611    Xmom = domain.quantities['xmomentum']
     612    Ymom = domain.quantities['ymomentum']
     613
     614    Friction = domain.quantities['friction']   
     615   
     616    for k in range(domain.number_of_elements):   
     617        w = Level.centroid_values[k]
     618        uh = Xmom.centroid_values[k]
     619        vh = Ymom.centroid_values[k]
     620        manning = Friction.centroid_values[k]
     621       
     622        if w >= minimum_allowed_height:
     623            S = -g * manning**2 * sqrt((uh**2 + vh**2))
     624            S /= w**(7.0/3)
     625
     626            #Update momentum
     627            Xmom.explicit_update[k] += S*uh
     628            Ymom.explicit_update[k] += S*vh
     629           
    871630
    872631##############################################
  • inundation/ga/storm_surge/pyvolution/test_quantity.py

    r198 r205  
    148148        quantity = Quantity(self.mesh4)
    149149
    150         #Set up for a gradient of (1,0)
    151         quantity.set_values([2.0/3, 4.0/3, 8.0/3, 2.0/3],
     150        #Set up for a gradient of (3,0) at mid triangle
     151        quantity.set_values([2.0, 4.0, 8.0, 2.0],
    152152                            location = 'centroids')
    153153       
    154154        #Gradients
    155         quantity.compute_gradients()
    156 
    157         #FIXME:
    158 
    159         #HERTIL
    160        
    161         #Gradient of fitted pwl surface   
    162         #a, b = compute_gradient(v2.id)
    163 
    164         #assert a == 1.0
    165         #assert b == 0.0
    166 
    167 
    168         #And now for the second order stuff       
    169         # - the full second order extrapolation
    170         #domain.order = 2
    171         #domain.limiter = dummy_limiter
    172         #distribute_to_vertices_and_edges(domain)       
    173    
    174         #assert v2.conserved_quantities_vertex0 == 0.0
    175         #assert v2.conserved_quantities_vertex1 == 2.0       
    176         #assert v2.conserved_quantities_vertex2 == 2.0       
     155        a, b = quantity.compute_gradients()
     156
     157
     158        #gradient bewteen t0 and t1 is undefined as det==0
     159        assert a[0] == 0.0
     160        assert b[0] == 0.0
     161        #The others are OK
     162        for i in range(1,4):
     163            assert a[i] == 3.0
     164            assert b[i] == 0.0
     165
     166
     167        quantity.second_order_extrapolator()
     168       
     169        assert allclose(quantity.vertex_values, [[2., 2.,  2.],
     170                                                 [0., 6.,  6.],
     171                                                 [6., 6., 12.],
     172                                                 [0., 0.,  6.]])
    177173
    178174
     
    355351
    356352
    357     def test_first_order_limiter(self):
     353    def test_first_order_extrapolator(self):
    358354        quantity = Quantity(self.mesh4)
    359355
     
    363359
    364360        #Extrapolate
    365         quantity.first_order_limiter()
     361        quantity.first_order_extrapolator()
    366362
    367363        #Check vertices but not edge values
     
    370366
    371367
    372     def test_second_order_limiter(self):
    373         quantity = Quantity(self.mesh4)
    374 
    375         #Create a deliberate overshoot (e.g. from gradient computation)
    376         quantity.set_values([[0,3,3], [6,2,2], [3,8,5], [3,5,8]])
    377 
    378         #Limit and extrapolate
    379         quantity.second_order_limiter()
     368    def test_second_order_extrapolator(self):
     369        quantity = Quantity(self.mesh4)
     370
     371        #Set up for a gradient of (3,0) at mid triangle
     372        quantity.set_values([2.0, 4.0, 8.0, 2.0],
     373                            location = 'centroids')
     374       
     375
     376
     377        quantity.second_order_extrapolator()
     378        quantity.limiter()
    380379
    381380
    382381        #Assert that central triangle is limited by neighbours
    383         assert quantity.vertex_values[1,0] <= quantity.vertex_values[2,2]
    384         assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
    385        
    386         assert quantity.vertex_values[1,1] <= quantity.vertex_values[3,0]
     382        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
     383        assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
     384       
     385        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
    387386        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
    388387       
    389388        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
    390         assert quantity.vertex_values[1,2] >= quantity.vertex_values[0,1]
     389        assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,1]
    391390
    392391       
     
    398397       
    399398       
    400         #Check vertices but not edge values
    401         #assert allclose(quantity.vertex_values,
    402         #                [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
     399
     400       
     401
     402    def test_limiter(self):       
     403        quantity = Quantity(self.mesh4)
     404
     405        #Create a deliberate overshoot (e.g. from gradient computation)
     406        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
     407
     408
     409        #Limit
     410        quantity.limiter()
     411
     412        #Assert that central triangle is limited by neighbours
     413        assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
     414        assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
     415       
     416        assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
     417        assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
     418       
     419        assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
     420        assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,1]
     421
     422
     423       
     424        #Assert that quantities are conserved
     425        from Numeric import sum
     426        for k in range(quantity.centroid_values.shape[0]):
     427            assert allclose (quantity.centroid_values[k],
     428                             sum(quantity.vertex_values[k,:])/3)
     429       
    403430
    404431
     
    412439
    413440        #Extrapolate
    414         quantity.distribute_to_vertices_and_edges(order=1)
     441        quantity.first_order_extrapolator()
     442
     443        #Interpolate
     444        quantity.interpolate_from_vertices_to_edges()       
    415445
    416446        assert allclose(quantity.vertex_values,
  • inundation/ga/storm_surge/pyvolution/test_shallow_water.py

    r198 r205  
    536536
    537537
    538     def test_first_order_limiter_const_z(self):
     538    def test_first_order_extrapolator_const_z(self):
    539539       
    540540        a = [0.0, 0.0]
     
    564564
    565565
    566         domain.first_order_limiter()
     566        domain.order = 1
     567        domain.distribute_to_vertices_and_edges()
    567568
    568569        #Check that centroid values were distributed to vertices
     
    608609        #- so that the limiter has something to work with
    609610        assert not alltrue(alltrue(greater_equal(L,E-epsilon)))               
    610                            
    611         domain.first_order_limiter()
     611
     612        domain.order = 1
     613        domain.distribute_to_vertices_and_edges()                           
    612614
    613615        #Check that all levels are above elevation (within eps)
  • inundation/ga/storm_surge/pyvolution/test_util.py

    r195 r205  
    2222        x2 = 0.0; y2 = 1.0; z2 = 0.0
    2323       
    24         zx, zy = gradient(x0, y0, x1, y1, x2, y2, array([z0]), array([z1]), array([z2]))
     24        zx, zy = gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2)
    2525
    2626        assert zx == -1.0
     
    5050        x2 = 2.0/3; y2 = 8.0/3
    5151
    52         q0 = array([2.0+2.0/3])
    53         q1 = array([8.0+2.0/3])
    54         q2 = array([2.0+8.0/3])
     52        q0 = 2.0+2.0/3
     53        q1 = 8.0+2.0/3
     54        q2 = 2.0+8.0/3
    5555       
    5656        #Gradient of fitted pwl surface   
    5757        a, b = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
    5858       
    59         assert abs(a[0] - 3.0) < epsilon
    60         assert abs(b[0] - 1.0) < epsilon
     59        assert abs(a - 3.0) < epsilon
     60        assert abs(b - 1.0) < epsilon
    6161   
    6262
     
    8484        x2 = 2.0/3; y2 = 8.0/3
    8585
    86         q0 = array([2.0+2.0/3])
    87         q1 = array([8.0+2.0/3])
    88         q2 = array([2.0+8.0/3])
     86        q0 = 2.0+2.0/3
     87        q1 = 8.0+2.0/3
     88        q2 = 2.0+8.0/3
    8989       
    9090        #Gradient of fitted pwl surface   
    9191        a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2)
    9292       
    93         assert abs(a[0] - 3.0) < epsilon
    94         assert abs(b[0] - 1.0) < epsilon
    95 
    96     def test_gradient_C_extension_scalars(self):
    97         from util_ext import gradient as gradient_c
    98        
    99         from RandomArray import uniform, seed
    100         seed(17, 53)
    101 
    102         x0, x1, x2, y0, y1, y2 = uniform(0.0,3.0,6)
    103 
    104         q0 = 10.0
    105         q1 = 3.0
    106         q2 = 7.0
    107        
    108        
    109         #Gradient of fitted pwl surface
    110         from util import gradient_python       
    111         a_ref, b_ref = gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2)
    112 
    113         #print a_ref, b_ref
    114         try:
    115             a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2)
    116         except TypeError:
    117             pass
    118         else:
    119             raise 'gradient c-extension should throw exception when q''s are scalars'
     93        assert abs(a - 3.0) < epsilon
     94        assert abs(b - 1.0) < epsilon
    12095
    12196       
     
    132107        q2 = uniform(7.0, 20.0, 4)       
    133108       
    134        
    135         #Gradient of fitted pwl surface
    136         from util import gradient_python               
    137         a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2)
    138109
    139         #print a_ref, b_ref
    140         a, b = gradient_c(x0, y0, x1, y1, x2, y2, q0, q1, q2)
     110        for i in range(4):
     111            #Gradient of fitted pwl surface
     112            from util import gradient_python           
     113            a_ref, b_ref = gradient(x0, y0, x1, y1, x2, y2,
     114                                    q0[i], q1[i], q2[i])
    141115
    142         #print a, a_ref, b, b_ref
    143         for i in range(4):
    144             assert abs(a[i] - a_ref[i]) < epsilon
    145             assert abs(b[i] - b_ref[i]) < epsilon
     116            #print a_ref, b_ref
     117            a, b = gradient_c(x0, y0, x1, y1, x2, y2,
     118                              q0[i], q1[i], q2[i])
     119
     120            #print a, a_ref, b, b_ref
     121            assert abs(a - a_ref) < epsilon
     122            assert abs(b - b_ref) < epsilon
    146123       
    147124
  • inundation/ga/storm_surge/pyvolution/util_ext.c

    r198 r205  
    1818              double x1, double y1,
    1919              double x2, double y2,
    20               double *q0, double *q1, double *q2,
    21               double *a, double *b,
    22               int N) {
     20              double q0, double q1, double q2,
     21              double *a, double *b) {
    2322
    2423  double det;
    25   int i;
    2624 
    2725  det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0);
    28  
    29   for (i=0; i<N; i++) {
    30     a[i] = (y2-y0)*(q1[i]-q0[i]) - (y1-y0)*(q2[i]-q0[i]);
    31     a[i] /= det;
    3226
    33     b[i] = (x1-x0)*(q2[i]-q0[i]) - (x2-x0)*(q1[i]-q0[i]);
    34     b[i] /= det;
    35   } 
     27  *a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0);
     28  *a /= det;
     29
     30  *b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0);
     31  *b /= det;
    3632
    3733  return 0;
     
    8177  //
    8278 
    83   double x0, y0, x1, y1, x2, y2;
    84   PyObject *Q0, *Q1, *Q2, *result;
    85   PyArrayObject *q0, *q1, *q2, *a, *b;
    86   int dimensions[1];
    87   int N;
    88 
     79  double x0, y0, x1, y1, x2, y2, q0, q1, q2, a, b;
     80  PyObject *result;
    8981 
    9082  // Convert Python arguments to C 
    91   if (!PyArg_ParseTuple(args, "ddddddOOO", &x0, &y0, &x1, &y1, &x2, &y2,
    92                         &Q0, &Q1, &Q2)) 
     83  if (!PyArg_ParseTuple(args, "ddddddddd", &x0, &y0, &x1, &y1, &x2, &y2,
     84                        &q0, &q1, &q2)) 
    9385    return NULL;
    94 
    95   //Input check
    96   if PyArray_Check(Q0) {
    97     q0 = (PyArrayObject *)
    98       PyArray_ContiguousFromObject(Q0, PyArray_DOUBLE, 0, 0);
    99   } else {
    100     PyErr_SetString(PyExc_TypeError, "q0 must be a numeric array");
    101     return NULL;
    102   }
    103    
    104   if PyArray_Check(Q1) {
    105     q1 = (PyArrayObject *)
    106       PyArray_ContiguousFromObject(Q1, PyArray_DOUBLE, 0, 0);
    107   } else {
    108     PyErr_SetString(PyExc_TypeError, "q1 must be a numeric array");
    109     return NULL;
    110   }
    111  
    112   if PyArray_Check(Q2) {
    113     q2 = (PyArrayObject *)
    114       PyArray_ContiguousFromObject(Q2, PyArray_DOUBLE, 0, 0);
    115   } else {
    116     PyErr_SetString(PyExc_TypeError, "q2 must be a numeric array");
    117     return NULL;
    118   } 
    119  
    120   //Allocate space for return vectors a and b 
    121  
    122   N = q0 -> dimensions[0];
    123   dimensions[0] = N;
    124   a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
    125   b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); 
    12686
    12787 
    12888  // Call underlying routine
    12989  _gradient(x0, y0, x1, y1, x2, y2,
    130             (double *) q0 -> data,
    131             (double *) q1 -> data,
    132             (double *) q2 -> data,
    133             (double *) a -> data,
    134             (double *) b-> data, N);
     90            q0, q1, q2, &a, &b);
    13591
    136   Py_DECREF(q0);
    137   Py_DECREF(q1);
    138   Py_DECREF(q2);   
    139              
    140   // Return arrays a and b
    141   result = Py_BuildValue("OO", a, b);
    142   Py_DECREF(a);
    143   Py_DECREF(b); 
     92  // Return values a and b
     93  result = Py_BuildValue("dd", a, b);
    14494  return result;
    14595}
Note: See TracChangeset for help on using the changeset viewer.