Changeset 907


Ignore:
Timestamp:
Feb 17, 2005, 5:36:29 PM (19 years ago)
Author:
ole
Message:

Addressed problem with artificial momentum generated by discontinuous water depths in the presence of steep slopes.
Now very shallow water is limited with a separate h-limiter controlled by beta_h
(see config.py) for details.

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

Legend:

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

    r844 r907  
    4545
    4646
     47#Betas [0;1] control the allowed steepness of gradient for second order
     48#extrapolations. Values of 1 allow the steepes gradients while
     49#lower values are more conservative. Values of 0 correspond to
     50#1'st order extrapolations.
     51#
     52# Large values of beta_h may cause simulations to require more timesteps
     53# as surface will 'hug' closer to the bed.
     54# Small values of beta_h will make code faster, but one may experience
     55# artificial momenta caused by discontinuities in water depths in
     56# the presence of steep slopes. One example of this would be
     57# stationary water 'lapping' upwards to a higher point on the coast.
     58#
     59#
     60#
     61#There are separate betas for the w-limiter and the h-limiter
     62#
     63#
     64#
     65#
     66#Good values are:
     67#beta_w = 0.9
     68#beta_h = 0.2
    4769
    48 beta = 0.9  #]0;1] (e.g. 0.9). 1 allows the steepest gradients,
    49             #lower values are the more conservative
    50             #the limiter being a first order method.
     70
     71
     72beta_w = 0.9
     73beta_h = 0.2
     74
    5175
    5276pmesh_filename = '.\\pmesh'
  • inundation/ga/storm_surge/pyvolution/domain.py

    r851 r907  
    5252
    5353        #Defaults
    54         from config import max_smallsteps, beta, epsilon
    55         self.beta = beta
     54        from config import max_smallsteps, beta_w, beta_h, epsilon
     55        self.beta_w = beta_w
     56        self.beta_h = beta_h       
    5657        self.epsilon = epsilon
     58
     59        #FIXME: Maybe have separate orders for h-limiter and w-limiter?
     60        #Or maybe get rid of order altogether and use beta_w and beta_h
    5761        self.default_order = 1
    5862        self.order = self.default_order
  • inundation/ga/storm_surge/pyvolution/netherlands.py

    r901 r907  
    9999#N = 140
    100100
    101 N = 15
     101#N = 15
    102102
    103103print 'Creating domain'
     
    110110domain.check_integrity()
    111111domain.default_order = 2
     112#domain.beta_h=0
    112113
    113114#Output params
  • inundation/ga/storm_surge/pyvolution/quantity.py

    r826 r907  
    813813    N = quantity.domain.number_of_elements
    814814   
    815     beta = quantity.domain.beta
     815    beta_w = quantity.domain.beta_w
    816816       
    817817    qc = quantity.centroid_values
     
    852852            if (dq[k,i] < 0): r = dqmin[k]/dq[k,i]
    853853           
    854             phi = min( min(r*beta, 1), phi )   
     854            phi = min( min(r*beta_w, 1), phi )   
    855855
    856856        #Then update using phi limiter
  • inundation/ga/storm_surge/pyvolution/quantity_ext.c

    r766 r907  
    1515
    1616//Shared code snippets
    17 #include "util_ext.h"
    18 #include "quantity_ext.h"
     17#include "util_ext.h" 
     18//#include "quantity_ext.h" //Obsolete
    1919
    2020
     
    450450   
    451451  int k, i, n, N, k3;
    452   double beta; //Safety factor
     452  double beta_w; //Safety factor
    453453  double *qmin, *qmax, qn;
    454454 
     
    464464  neighbours = get_consecutive_array(domain, "neighbours");
    465465
    466   //Get safety factor beta
    467   Tmp = PyObject_GetAttrString(domain, "beta");
     466  //Get safety factor beta_w
     467  Tmp = PyObject_GetAttrString(domain, "beta_w");
    468468  if (!Tmp)
    469469    return NULL;
    470470     
    471   beta = PyFloat_AsDouble(Tmp); 
     471  beta_w = PyFloat_AsDouble(Tmp); 
    472472 
    473473  Py_DECREF(Tmp);   
     
    502502 
    503503  // Call underlying routine
    504   _limit(N, beta, (double*) qc -> data, (double*) qv -> data, qmin, qmax);
     504  _limit(N, beta_w, (double*) qc -> data, (double*) qv -> data, qmin, qmax);
    505505         
    506506  free(qmin);
  • inundation/ga/storm_surge/pyvolution/quantity_ext.h

    r258 r907  
    1 
    2 void _limit(int N, double beta, double* qc, double* qv,
    3             double* qmin, double* qmax) {
    4 
    5   int k, i, k3;
    6   double dq, dqa[3], phi, r;
    7  
    8   for (k=0; k<N; k++) {
    9     k3 = k*3;
    10    
    11     //Find the gradient limiter (phi) across vertices 
    12     phi = 1.0;
    13     for (i=0; i<3; i++) {   
    14       r = 1.0;
    15      
    16       dq = qv[k3+i] - qc[k];    //Delta between vertex and centroid values
    17       dqa[i] = dq;              //Save dq for use in the next loop
    18      
    19       if (dq > 0.0) r = (qmax[k] - qc[k])/dq;
    20       if (dq < 0.0) r = (qmin[k] - qc[k])/dq;     
    21  
    22  
    23       phi = min( min(r*beta, 1.0), phi);   
    24     }
    25    
    26     //Then update using phi limiter
    27     for (i=0; i<3; i++) {   
    28       qv[k3+i] = qc[k] + phi*dqa[i];
    29     }
    30   }
    31 }
  • inundation/ga/storm_surge/pyvolution/shallow_water.py

    r901 r907  
    521521            raise 'Unknown order'
    522522
    523     #Take bed elevation into account when water heights are small   
     523    #Take bed elevation into account when water heights are small
    524524    balance_deep_and_shallow(domain)
    525525
     
    594594    protect(domain.minimum_allowed_height, wc, zc, xmomc, ymomc)
    595595   
    596    
    597 def balance_deep_and_shallow_old(domain):
    598     """Compute linear combination between stage as computed by
    599     gradient-limiters and stage computed as constant height above bed.
    600     The former takes precedence when heights are large compared to the
    601     bed slope while the latter takes precedence when heights are
    602     relatively small.  Anything in between is computed as a balanced
    603     linear combination in order to avoid numerical disturbances which
    604     would otherwise appear as a result of hard switching between
    605     modes.
    606     """
    607    
    608     #Shortcuts
    609     wc = domain.quantities['stage'].centroid_values
    610     zc = domain.quantities['elevation'].centroid_values
    611     hc = wc - zc
    612    
    613     wv = domain.quantities['stage'].vertex_values
    614     zv = domain.quantities['elevation'].vertex_values
    615     hv = wv-zv
    616 
    617 
    618     #Computed linear combination between constant stages and and
    619     #stages parallel to the bed elevation.     
    620     for k in range(domain.number_of_elements):
    621         #Compute maximal variation in bed elevation
    622         #  This quantitiy is
    623         #    dz = max_i abs(z_i - z_c)
    624         #  and it is independent of dimension
    625         #  In the 1d case zc = (z0+z1)/2
    626         #  In the 2d case zc = (z0+z1+z2)/3
    627 
    628         dz = max(abs(zv[k,0]-zc[k]),
    629                  abs(zv[k,1]-zc[k]),
    630                  abs(zv[k,2]-zc[k]))
    631 
    632        
    633         hmin = min( hv[k,:] )
    634 
    635         #Create alpha in [0,1], where alpha==0 means using shallow
    636         #first order scheme and alpha==1 means using the stage w as
    637         #computed by the gradient limiter (1st or 2nd order)
    638         #
    639         #If hmin > dz/2 then alpha = 1 and the bed will have no effect
    640         #If hmin < 0 then alpha = 0 reverting to constant height above bed.
    641      
    642         if dz > 0.0:
    643             alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
    644         else:
    645             #Flat bed
    646             alpha = 1.0
    647 
    648         #Weighted balance between stage parallel to bed elevation
    649         #(wvi = zvi + hc) and stage as computed by 1st or 2nd
    650         #order gradient limiter
    651         #(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids
    652         #
    653         #It follows that the updated wvi is
    654         #  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) =
    655         #  zvi + hc + alpha*(hvi - hc)
    656         #
    657         #Note that hvi = zc+hc-zvi in the first order case (constant).
    658 
    659         if alpha < 1:         
    660             for i in range(3):
    661                 wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
    662            
    663 
    664             #Momentums at centroids
    665             xmomc = domain.quantities['xmomentum'].centroid_values
    666             ymomc = domain.quantities['ymomentum'].centroid_values       
    667 
    668             #Momentums at vertices
    669             xmomv = domain.quantities['xmomentum'].vertex_values
    670             ymomv = domain.quantities['ymomentum'].vertex_values         
    671 
    672             # Update momentum as a linear combination of
    673             # xmomc and ymomc (shallow) and momentum
    674             # from extrapolator xmomv and ymomv (deep).
    675             xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]           
    676             ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
    677            
    678 
    679 def limit_on_h(domain):       
     596
     597
     598def h_limiter(domain):       
    680599    """Limit slopes for each volume to eliminate artificial variance
    681600    introduced by e.g. second order extrapolator
    682601   
    683602    limit on h = w-z
     603
     604    This limiter depends on two quantities (w,z) so it resides within
     605    this module rather than within quantity.py
    684606    """
    685607
     
    687609
    688610    N = domain.number_of_elements
    689     beta = domain.beta   
     611    beta_h = domain.beta_h
    690612   
    691613    #Shortcuts
     
    734656            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]
    735657           
    736             phi = min( min(r*beta, 1), phi )   
     658            phi = min( min(r*beta_h, 1), phi )   
    737659
    738660        #Then update using phi limiter
     
    741663
    742664    return hvbar
    743    
     665
     666
     667
     668def h_limiter_c(domain):       
     669    """Limit slopes for each volume to eliminate artificial variance
     670    introduced by e.g. second order extrapolator
     671   
     672    limit on h = w-z
     673
     674    This limiter depends on two quantities (w,z) so it resides within
     675    this module rather than within quantity.py
     676
     677    Wrapper for c-extension
     678    """
     679
     680    from Numeric import zeros, Float
     681
     682    N = domain.number_of_elements
     683    beta_h = domain.beta_h
     684   
     685    #Shortcuts
     686    wc = domain.quantities['stage'].centroid_values
     687    zc = domain.quantities['elevation'].centroid_values
     688    hc = wc - zc
     689   
     690    wv = domain.quantities['stage'].vertex_values
     691    zv = domain.quantities['elevation'].vertex_values
     692    hv = wv-zv
     693
     694    #Call C-extension
     695    from shallow_water_ext import h_limiter
     696    hvbar = h_limiter(domain, hc, hv)
     697
     698    return hvbar
     699
    744700                                           
    745701def balance_deep_and_shallow(domain):
     
    776732    hv = wv-zv
    777733
     734    #Limit h
     735    hvbar = h_limiter(domain)
     736
    778737    for k in range(domain.number_of_elements):
    779738        #Compute maximal variation in bed elevation
     
    791750        hmin = min( hv[k,:] )
    792751
    793         #Create alpha in [0,1], where alpha==0 means using shallow
    794         #first order scheme and alpha==1 means using the stage w as
    795         #computed by the gradient limiter (1st or 2nd order)
    796         #
     752        #Create alpha in [0,1], where alpha==0 means using the h-limited
     753        #stage and alpha==1 means using the w-limited stage as
     754        #computed by the gradient limiter (both 1st or 2nd order)
     755   
    797756        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
    798757        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
    799      
     758
    800759        if dz > 0.0:
    801760            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
     
    805764
    806765        #Let
    807         #
     766        #
     767        #  wvi be the w-limited stage (wvi = zvi + hvi)       
    808768        #  wvi- be the h-limited state (wvi- = zvi + hvi-)
    809         #  wvi~ be the w-limited stage (wvi~ = zvi + hvi~)
    810769        # 
    811770        #   
     
    814773        #Weighted balance between w-limited and h-limited stage is 
    815774        #
    816         #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi~)   
     775        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)   
    817776        #
    818777        #It follows that the updated wvi is
    819         #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi~
     778        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
    820779        #
    821         # Momentum is balanced between constant and limited ???
    822 
    823         #print 'Alpha = ', alpha
    824         #FIXME: if alpha == 0 compute only h-limited slope,
    825         #if alpha == 1 compute only w-limited slope
    826         #
     780        # Momentum is balanced between constant and limited
     781
    827782        if alpha < 1:         
    828783       
    829             #Limit h
    830             hvbar = limit_on_h(domain)
    831    
    832784            for i in range(3):
    833785                wv[k,i] = zv[k,i] + (1-alpha)*hvbar[k,i] + alpha*hv[k,i]
     
    848800                                   
    849801           
    850 
    851802def balance_deep_and_shallow_c(domain):
    852803    """Wrapper for C implementation
     
    870821    ymomv = domain.quantities['ymomentum'].vertex_values         
    871822
    872    
     823    #Limit h
     824    hvbar = h_limiter(domain)   
    873825
    874826    from shallow_water_ext import balance_deep_and_shallow
    875     balance_deep_and_shallow(wc, zc, hc, wv, zv, hv,
    876                              xmomc, ymomc, xmomv, ymomv)
     827    #balance_deep_and_shallow(wc, zc, hc, wv, zv, hv,
     828    #                         xmomc, ymomc, xmomv, ymomv)
     829
     830    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv, hvbar,
     831                             xmomc, ymomc, xmomv, ymomv)   
    877832
    878833   
     
    12721227        ymom_update[k] += S*v       
    12731228           
     1229
     1230
     1231##############################
     1232#OBSOLETE STUFF
     1233   
     1234def balance_deep_and_shallow_old(domain):
     1235    """Compute linear combination between stage as computed by
     1236    gradient-limiters and stage computed as constant height above bed.
     1237    The former takes precedence when heights are large compared to the
     1238    bed slope while the latter takes precedence when heights are
     1239    relatively small.  Anything in between is computed as a balanced
     1240    linear combination in order to avoid numerical disturbances which
     1241    would otherwise appear as a result of hard switching between
     1242    modes.
     1243    """
     1244
     1245    #OBSOLETE
     1246   
     1247    #Shortcuts
     1248    wc = domain.quantities['stage'].centroid_values
     1249    zc = domain.quantities['elevation'].centroid_values
     1250    hc = wc - zc
     1251   
     1252    wv = domain.quantities['stage'].vertex_values
     1253    zv = domain.quantities['elevation'].vertex_values
     1254    hv = wv-zv
     1255
     1256
     1257    #Computed linear combination between constant stages and and
     1258    #stages parallel to the bed elevation.     
     1259    for k in range(domain.number_of_elements):
     1260        #Compute maximal variation in bed elevation
     1261        #  This quantitiy is
     1262        #    dz = max_i abs(z_i - z_c)
     1263        #  and it is independent of dimension
     1264        #  In the 1d case zc = (z0+z1)/2
     1265        #  In the 2d case zc = (z0+z1+z2)/3
     1266
     1267        dz = max(abs(zv[k,0]-zc[k]),
     1268                 abs(zv[k,1]-zc[k]),
     1269                 abs(zv[k,2]-zc[k]))
     1270
     1271       
     1272        hmin = min( hv[k,:] )
     1273
     1274        #Create alpha in [0,1], where alpha==0 means using shallow
     1275        #first order scheme and alpha==1 means using the stage w as
     1276        #computed by the gradient limiter (1st or 2nd order)
     1277        #
     1278        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
     1279        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
     1280     
     1281        if dz > 0.0:
     1282            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
     1283        else:
     1284            #Flat bed
     1285            alpha = 1.0
     1286
     1287        #Weighted balance between stage parallel to bed elevation
     1288        #(wvi = zvi + hc) and stage as computed by 1st or 2nd
     1289        #order gradient limiter
     1290        #(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids
     1291        #
     1292        #It follows that the updated wvi is
     1293        #  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) =
     1294        #  zvi + hc + alpha*(hvi - hc)
     1295        #
     1296        #Note that hvi = zc+hc-zvi in the first order case (constant).
     1297
     1298        if alpha < 1:         
     1299            for i in range(3):
     1300                wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
     1301           
     1302
     1303            #Momentums at centroids
     1304            xmomc = domain.quantities['xmomentum'].centroid_values
     1305            ymomc = domain.quantities['ymomentum'].centroid_values       
     1306
     1307            #Momentums at vertices
     1308            xmomv = domain.quantities['xmomentum'].vertex_values
     1309            ymomv = domain.quantities['ymomentum'].vertex_values         
     1310
     1311            # Update momentum as a linear combination of
     1312            # xmomc and ymomc (shallow) and momentum
     1313            # from extrapolator xmomv and ymomv (deep).
     1314            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]           
     1315            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
     1316           
    12741317
    12751318
     
    14571500    gravity = gravity_c
    14581501    manning_friction = manning_friction_c
     1502    h_limiter = h_limiter_c
    14591503    balance_deep_and_shallow = balance_deep_and_shallow_c
    1460     #balance_deep_and_shallow = balance_deep_and_shallow_old
    14611504    protect_against_infinitesimal_and_negative_heights = protect_against_infinitesimal_and_negative_heights_c
    14621505
  • inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

    r773 r907  
    190190                              double* zv,
    191191                              double* hv,
     192                              double* hvbar,                         
    192193                              double* xmomc,
    193194                              double* ymomc,
     
    198199  double dz, hmin, alpha;
    199200 
    200   //Compute linear combination between constant stages and and
    201   //stages parallel to the bed elevation.     
     201  //Compute linear combination between w-limited stages and
     202  //h-limited stages close to the bed elevation.     
    202203 
    203204  for (k=0; k<N; k++) {
     
    220221
    221222   
    222     //Create alpha in [0,1], where alpha==0 means using shallow
    223     //first order scheme and alpha==1 means using the stage w as
    224     //computed by the gradient limiter (1st or 2nd order)
     223    //Create alpha in [0,1], where alpha==0 means using the h-limited
     224    //stage and alpha==1 means using the w-limited stage as
     225    //computed by the gradient limiter (both 1st or 2nd order)
    225226    //
    226227    //If hmin > dz/2 then alpha = 1 and the bed will have no effect
     
    232233      alpha = 1.0;  //Flat bed
    233234
    234      
    235     //Weighted balance between stage parallel to bed elevation
    236     //(wvi = zvi + hc) and stage as computed by 1st or 2nd
    237     //order gradient limiter
    238     //(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids
    239     //
    240     //It follows that the updated wvi is
    241     //  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) =
    242     //  zvi + hc + alpha*(hvi - hc)
    243     //
    244     //Note that hvi = zc+hc-zvi in the first order case (constant).
     235
     236    //  Let
     237    // 
     238    //    wvi be the w-limited stage (wvi = zvi + hvi)       
     239    //    wvi- be the h-limited state (wvi- = zvi + hvi-)
     240    //     
     241    //     
     242    //  where i=0,1,2 denotes the vertex ids
     243    // 
     244    //  Weighted balance between w-limited and h-limited stage is 
     245    //   
     246    //    wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)   
     247    // 
     248    //  It follows that the updated wvi is
     249    //    wvi := zvi + (1-alpha)*hvi- + alpha*hvi
     250    //   
     251    //   Momentum is balanced between constant and limited
    245252
    246253    if (alpha < 1) {         
    247254      for (i=0; i<3; i++) {
    248         wv[k3+i] = zv[k3+i] + hc[k] + alpha*(hv[k3+i]-hc[k]);
     255        wv[k3+i] = zv[k3+i] + (1-alpha)*hvbar[k3+i] + alpha*hv[k3+i];
    249256           
    250257     
     
    286293  return 0;
    287294}
     295
     296
     297
     298
     299
    288300
    289301
     
    683695    *wv,            //Stage at vertices
    684696    *zv,            //Elevation at vertices
    685     *hv,            //Heights at vertices   
     697    *hv,            //Depths at vertices   
     698    *hvbar,         //h-Limited depths at vertices       
    686699    *xmomc,         //Momentums at centroids and vertices
    687700    *ymomc,
     
    692705 
    693706  // Convert Python arguments to C 
    694   if (!PyArg_ParseTuple(args, "OOOOOOOOOO", 
     707  if (!PyArg_ParseTuple(args, "OOOOOOOOOOO", 
    695708                        &wc, &zc, &hc,
    696                         &wv, &zv, &hv,
     709                        &wv, &zv, &hv, &hvbar,
    697710                        &xmomc, &ymomc, &xmomv, &ymomv))
    698711    return NULL;
     
    707720                            (double*) zv -> data,
    708721                            (double*) hv -> data,
     722                            (double*) hvbar -> data,                       
    709723                            (double*) xmomc -> data,
    710724                            (double*) ymomc -> data,
     
    718732
    719733
     734PyObject *h_limiter(PyObject *self, PyObject *args) {
     735 
     736  PyObject *domain, *Tmp;
     737  PyArrayObject
     738    *hv, *hc, //Depth at vertices and centroids       
     739    *hvbar,   //Limited depth at vertices (return values)
     740    *neighbours;
     741   
     742  int k, i, n, N, k3;
     743  int dimensions[2];
     744  double beta_h; //Safety factor (see config.py)
     745  double *hmin, *hmax, hn;
     746 
     747  // Convert Python arguments to C 
     748  if (!PyArg_ParseTuple(args, "OOO", &domain, &hc, &hv))
     749    return NULL;
     750
     751  neighbours = get_consecutive_array(domain, "neighbours");
     752
     753  //Get safety factor beta_h
     754  Tmp = PyObject_GetAttrString(domain, "beta_h");
     755  if (!Tmp)
     756    return NULL;
     757     
     758  beta_h = PyFloat_AsDouble(Tmp); 
     759 
     760  Py_DECREF(Tmp);   
     761
     762  N = hc -> dimensions[0];
     763 
     764  //Create hvbar
     765  dimensions[0] = N;
     766  dimensions[1] = 3; 
     767  hvbar = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
     768
     769     
     770  //Find min and max of this and neighbour's centroid values
     771  hmin = malloc(N * sizeof(double));
     772  hmax = malloc(N * sizeof(double)); 
     773  for (k=0; k<N; k++) {
     774    k3=k*3;
     775   
     776    hmin[k] = ((double*) hc -> data)[k];
     777    hmax[k] = hmin[k];
     778   
     779    for (i=0; i<3; i++) {
     780      n = ((long*) neighbours -> data)[k3+i];
     781      if (n >= 0) {
     782        hn = ((double*) hc -> data)[n]; //Neighbour's centroid value
     783         
     784        hmin[k] = min(hmin[k], hn);
     785        hmax[k] = max(hmax[k], hn);
     786      }
     787    }
     788  }
     789 
     790 
     791  // Call underlying routine
     792  _limit(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax);
     793         
     794  // // //Py_DECREF(domain); //FIXME: NEcessary?         
     795  free(hmin);
     796  free(hmax); 
     797 
     798  //return result using PyArray to avoid memory leak 
     799  return PyArray_Return(hvbar);
     800}
    720801
    721802
     
    776857  {"balance_deep_and_shallow", balance_deep_and_shallow,
    777858   METH_VARARGS, "Print out"},         
     859  {"h_limiter", h_limiter,
     860   METH_VARARGS, "Print out"},             
    778861  {"protect", protect, METH_VARARGS | METH_KEYWORDS, "Print out"},   
    779862  {"assign_windfield_values", assign_windfield_values,
  • inundation/ga/storm_surge/pyvolution/test_data_manager.py

    r849 r907  
    414414        self.domain.smooth = False
    415415        self.domain.store = True
     416        self.domain.beta_h = 0
    416417
    417418        #Evolution
  • inundation/ga/storm_surge/pyvolution/test_interpolate_sww.py

    r867 r907  
    2828        domain = Domain(points, vertices, boundary)
    2929        domain.default_order=2
     30        domain.beta_h = 0
    3031
    3132       
     
    4344        ######################
    4445        #Initial condition - with jumps
    45 
    4646
    4747        bed = domain.quantities['elevation'].vertex_values
  • inundation/ga/storm_surge/pyvolution/test_shallow_water.py

    r901 r907  
    14341434
    14351435        #print E
    1436         domain.order = 2       
     1436        domain.order = 2
     1437        domain.beta_h = 0.0 #Use first order in h-limiter       
    14371438        domain.distribute_to_vertices_and_edges()
    14381439
     
    18371838        domain.smooth = False
    18381839        domain.default_order=1
     1840        domain.beta_h = 0.0 #Use first order in h-limiter
    18391841
    18401842        #Bed-slope and friction
     
    19861988        domain.smooth = False
    19871989        domain.default_order=2
     1990        domain.beta_h = 0.0 #Use first order in h-limiter       
    19881991
    19891992        #Bed-slope and friction at vertices (and interpolated elsewhere)
     
    20782081        domain.smooth = False
    20792082        domain.default_order=2
     2083        domain.beta_h = 0.0 #Use first order in h-limiter       
    20802084
    20812085        #Bed-slope and friction at vertices (and interpolated elsewhere)
     
    21662170        domain.smooth = False
    21672171        domain.default_order=2
     2172        domain.beta_h = 0.0 #Use first order in h-limiter       
    21682173
    21692174        #Bed-slope and friction at vertices (and interpolated elsewhere)
     
    22822287        domain.smooth = False
    22832288        domain.default_order=2
     2289        domain.beta_h = 0.0 #Use first order in h-limiter       
    22842290
    22852291        #Bed-slope and friction at vertices (and interpolated elsewhere)
  • inundation/ga/storm_surge/pyvolution/util_ext.h

    r753 r907  
    5151
    5252
     53void _limit(int N, double beta, double* qc, double* qv,
     54            double* qmin, double* qmax) {
     55
     56  //N are the number of elements
     57  int k, i, k3;
     58  double dq, dqa[3], phi, r;
     59 
     60  //printf("INSIDE\n");
     61  for (k=0; k<N; k++) {
     62    k3 = k*3;
     63   
     64    //Find the gradient limiter (phi) across vertices 
     65    phi = 1.0;
     66    for (i=0; i<3; i++) {   
     67      r = 1.0;
     68     
     69      dq = qv[k3+i] - qc[k];    //Delta between vertex and centroid values
     70      dqa[i] = dq;              //Save dq for use in the next loop
     71     
     72      if (dq > 0.0) r = (qmax[k] - qc[k])/dq;
     73      if (dq < 0.0) r = (qmin[k] - qc[k])/dq;     
     74 
     75 
     76      phi = min( min(r*beta, 1.0), phi);   
     77    }
     78   
     79    //Then update using phi limiter
     80    for (i=0; i<3; i++) {   
     81      qv[k3+i] = qc[k] + phi*dqa[i];
     82    }
     83  }
     84}
     85
     86
    5387void print_numeric_array(PyArrayObject *x) { 
    5488  int i, j;
Note: See TracChangeset for help on using the changeset viewer.