Changeset 266


Ignore:
Timestamp:
Sep 1, 2004, 6:31:42 PM (21 years ago)
Author:
ole
Message:
 
Location:
inundation/ga/storm_surge/pyvolution
Files:
2 edited

Legend:

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

    r263 r266  
    474474
    475475def balance_deep_and_shallow(domain):
    476     """Compute linear combination between stage as computed by gradient-limiters and
    477     stage computed as constant height above bed.
    478     The former takes precedence when heights are large compared to the bed slope while the latter
    479     takes precedence when heights are relatively small.
    480     Anything in between is computed as a balanced linear combination in order to avoid numerical
    481     disturbances which would otherwise appear as a result of hard switching between modes.
     476    """Compute linear combination between stage as computed by
     477    gradient-limiters and stage computed as constant height above bed.
     478    The former takes precedence when heights are large compared to the
     479    bed slope while the latter takes precedence when heights are
     480    relatively small.  Anything in between is computed as a balanced
     481    linear combination in order to avoid numerical disturbances which
     482    would otherwise appear as a result of hard switching between
     483    modes.
    482484    """
    483485   
     
    553555            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:];                     
    554556           
     557
     558def balance_deep_and_shallow_c(domain):
     559    """Wrapper for C implementation
     560    """
     561   
     562    #Shortcuts
     563    wc = domain.quantities['level'].centroid_values
     564    zc = domain.quantities['elevation'].centroid_values
     565    hc = wc - zc
     566   
     567    wv = domain.quantities['level'].vertex_values
     568    zv = domain.quantities['elevation'].vertex_values
     569    hv = wv-zv
     570
     571    #Momentums at centroids
     572    xmomc = domain.quantities['xmomentum'].centroid_values
     573    ymomc = domain.quantities['ymomentum'].centroid_values       
     574
     575    #Momentums at vertices
     576    xmomv = domain.quantities['xmomentum'].vertex_values
     577    ymomv = domain.quantities['ymomentum'].vertex_values         
     578
     579   
     580
     581    from shallow_water_ext import balance_deep_and_shallow
     582    balance_deep_and_shallow(domain.number_of_elements,
     583                             wc, zc, wv, zv, hv,
     584                             xmomc, ymomc, xmomv, ymomv)
     585
    555586   
    556587
     
    954985    gravity = gravity_c
    955986    manning_friction = manning_friction_c
    956    
     987    #balance_deep_and_shallow = balance_deep_and_shallow_c
    957988   
    958989    #distribute_to_vertices_and_edges = distribute_to_vertices_and_edges_c
  • inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

    r260 r266  
    176176           
    177177
     178/*
     179int _balance_deep_and_shallow()     
     180
     181    #Computed linear combination between constant levels and and
     182    #levels parallel to the bed elevation.     
     183    for k in range(domain.number_of_elements):
     184        #Compute maximal variation in bed elevation
     185        #  This quantitiy is
     186        #    dz = max_i abs(z_i - z_c)
     187        #  and it is independent of dimension
     188        #  In the 1d case zc = (z0+z1)/2
     189        #  In the 2d case zc = (z0+z1+z2)/3
     190
     191        dz = max(abs(zv[k,0]-zc[k]),
     192                 abs(zv[k,1]-zc[k]),
     193                 abs(zv[k,2]-zc[k]))
     194
     195       
     196        hmin = min( hv[k,:] )
     197
     198        #Create alpha in [0,1], where alpha==0 means using shallow
     199        #first order scheme and alpha==1 means using the stage w as
     200        #computed by the gradient limiter (1st or 2nd order)
     201        #
     202        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
     203        #If hmin < 0 then alpha = 0 reverting to constant height above bed.
     204     
     205        if dz > 0.0:
     206            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
     207        else:
     208            #Flat bed
     209            alpha = 1.0
     210
     211
     212        #Weighted balance between stage parallel to bed elevation
     213        #(wvi = zvi + hc) and stage as computed by 1st or 2nd
     214        #order gradient limiter
     215        #(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids
     216        #
     217        #It follows that the updated wvi is
     218        #  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) =
     219        #  zvi + hc + alpha*(hvi - hc)
     220        #
     221        #Note that hvi = zc+hc-zvi in the first order case (constant).
     222
     223        if alpha < 1:         
     224            for i in range(3):
     225                wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k])
     226           
     227
     228            # Update momentum as a linear combination of
     229            # xmomc and ymomc (shallow) and momentum
     230            # from extrapolator xmomv and ymomv (deep).
     231            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:];           
     232            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:];                     
     233*/           
     234           
    178235       
    179236///////////////////////////////////////////////////////////////////
     
    482539
    483540
     541PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) {
     542  //
     543  //    balance_deep_and_shallow(domain.number_of_elements,
     544  //                             wc, zc, wv, zv, hv,
     545  //                             xmomc, ymomc, xmomv, ymomv)
     546 
     547
     548  PyArrayObject
     549    *wc,            //Level at centroids
     550    *zc,            //Elevation at centroids   
     551    *wv,            //Level at vertices
     552    *zv,            //Elevation at vertices
     553    *hv,            //Heights at vertices   
     554    *xmomc,         //Momentums at centroids and vertices
     555    *ymomc,
     556    *xmomv,
     557    *ymomv;   
     558   
     559  int N; //, err;
     560 
     561  // Convert Python arguments to C 
     562  if (!PyArg_ParseTuple(args, "dOOOOOOOOO", &N, &wc, &zc, &wv, &zv, &hv,
     563                        &xmomc, &ymomc, &xmomv, &ymomv))
     564    return NULL;
     565
     566 
     567  /*
     568  _balance_deep_and_shallow(N,
     569                            (double*) wc -> data,
     570                            (double*) zc -> data,
     571                            (double*) wv -> data,
     572                            (double*) zv -> data,
     573                            (double*) hv -> data,
     574                            (double*) xmomc -> data,
     575                            (double*) ymomc -> data,
     576                            (double*) xmomv -> data,
     577                            (double*) ymomv -> data); 
     578  */
     579 
     580  return Py_BuildValue(""); 
     581}
    484582
    485583
     
    499597  {"gravity", gravity, METH_VARARGS, "Print out"},     
    500598  {"manning_friction", manning_friction, METH_VARARGS, "Print out"},       
     599  {"balance_deep_and_shallow", balance_deep_and_shallow,
     600   METH_VARARGS, "Print out"},         
    501601  //{"distribute_to_vertices_and_edges",
    502602  // distribute_to_vertices_and_edges, METH_VARARGS},   
Note: See TracChangeset for help on using the changeset viewer.