Changeset 245


Ignore:
Timestamp:
Aug 30, 2004, 6:44:45 PM (21 years ago)
Author:
ole
Message:

Started to work on C implementation of phi limiter

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

Legend:

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

    r242 r245  
    333333
    334334    def limit(self):
    335         """Limit slopes for each volume to eliminate artificial variance
    336         introduced by e.g. second order extrapolator
    337        
    338         This is an unsophisticated limiter as it does not take into
    339         account dependencies among quantities.
    340        
    341         precondition:
    342            vertex values are estimated from gradient
    343         postcondition:
    344            vertex values are updated
    345         """
    346 
    347         from Numeric import zeros, Float
    348 
    349         beta = self.domain.beta
    350        
    351         qc = self.centroid_values
    352         qv = self.vertex_values
    353        
    354         #Deltas between vertex and centroid values
    355         dq = zeros( qv.shape, Float)
    356         for i in range(3):
    357             dq[:,i] = qv[:,i] - qc
    358        
    359         #Find min and max of this and neighbour's centroid values
    360         qmax = zeros(qc.shape, Float)
    361         qmin = zeros(qc.shape, Float)       
    362 
    363         for k in range(self.domain.number_of_elements):
    364             qmax[k] = qmin[k] = qc[k]
    365             for i in range(3):
    366                 n = self.domain.neighbours[k,i]
    367                 if n >= 0:
    368                     qn = qc[n] #Neighbour's centroid value
    369      
    370                     qmin[k] = min(qmin[k], qn)
    371                     qmax[k] = max(qmax[k], qn)
    372      
    373    
    374         #Diffences between centroids and maxima/minima
    375         dqmax = qmax - qc
    376         dqmin = qmin - qc
    377  
    378         for k in range(self.domain.number_of_elements):
    379 
    380             #Find the gradient limiter (phi) across vertices
    381             phi = 1.0
    382             for i in range(3):
    383                 r = 1.0
    384                 if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
    385                 if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
    386                
    387                 phi = min( min(r*beta, 1), phi )   
    388 
    389             #Then update using phi limiter
    390             for i in range(3):           
    391                 qv[k,i] = qc[k] + phi*dq[k,i]
    392 
     335        #Call correct module function
     336        #(either from this module or C-extension)
     337        limit(self)
     338       
    393339       
    394340    def extrapolate_first_order(self):
     
    431377
    432378           
     379
     380def limit(quantity):       
     381    """Limit slopes for each volume to eliminate artificial variance
     382    introduced by e.g. second order extrapolator
     383   
     384    This is an unsophisticated limiter as it does not take into
     385    account dependencies among quantities.
     386   
     387    precondition:
     388    vertex values are estimated from gradient
     389    postcondition:
     390    vertex values are updated
     391    """
     392
     393    from Numeric import zeros, Float
     394
     395    N = quantity.domain.number_of_elements
     396   
     397    beta = quantity.domain.beta
     398       
     399    qc = quantity.centroid_values
     400    qv = quantity.vertex_values
     401       
     402    #Find min and max of this and neighbour's centroid values
     403    qmax = zeros(qc.shape, Float)
     404    qmin = zeros(qc.shape, Float)       
     405   
     406    for k in range(N):
     407        qmax[k] = qmin[k] = qc[k]
     408        for i in range(3):
     409            n = quantity.domain.neighbours[k,i]
     410            if n >= 0:
     411                qn = qc[n] #Neighbour's centroid value
     412               
     413                qmin[k] = min(qmin[k], qn)
     414                qmax[k] = max(qmax[k], qn)
     415     
     416   
     417    #Diffences between centroids and maxima/minima
     418    dqmax = qmax - qc
     419    dqmin = qmin - qc
     420       
     421    #Deltas between vertex and centroid values
     422    dq = zeros(qv.shape, Float)
     423    for i in range(3):
     424        dq[:,i] = qv[:,i] - qc
     425
     426    #Phi limiter   
     427    for k in range(N):
     428       
     429        #Find the gradient limiter (phi) across vertices
     430        phi = 1.0
     431        for i in range(3):
     432            r = 1.0
     433            if (dq[k,i] > 0): r = dqmax[k]/dq[k,i]
     434            if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
     435           
     436            phi = min( min(r*beta, 1), phi )   
     437
     438        #Then update using phi limiter
     439        for i in range(3):           
     440            qv[k,i] = qc[k] + phi*dq[k,i]
     441
     442
     443
     444import compile
     445if compile.can_use_C_extension('quantity_ext.c'):
     446    #Replace python version with c implementations
     447
     448    pass
     449    #from quantity_ext import limit
     450
  • inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

    r240 r245  
    1717#include "math.h"
    1818
     19//FIXME: Should live in util_ext.c
    1920double max(double x, double y) { 
    2021  //Return maximum of two doubles
     
    3334
    3435// Computational function for rotation
    35 // FIXME: Should probably go into shallow_water_ext.c
    3636int _rotate(double *q, double n1, double n2) {
    3737  /*Rotate the momentum component q (q[1], q[2])
Note: See TracChangeset for help on using the changeset viewer.