Changeset 267


Ignore:
Timestamp:
Sep 2, 2004, 3:00:07 AM (21 years ago)
Author:
ole
Message:

C implementation of balanced_deep_and_shallow.

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

Legend:

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

    r240 r267  
    11import os
     2
     3#Clean up
     4
     5for file in os.listdir('.'):
     6    if file[-1] == '~':
     7        os.remove(file)
     8
    29
    310os.system('python compile.py')
  • inundation/ga/storm_surge/pyvolution/shallow_water.py

    r266 r267  
    511511        hmin = min( hv[k,:] )
    512512
     513        print dz, hmin
     514       
    513515        #Create alpha in [0,1], where alpha==0 means using shallow
    514516        #first order scheme and alpha==1 means using the stage w as
     
    552554            # xmomc and ymomc (shallow) and momentum
    553555            # from extrapolator xmomv and ymomv (deep).
    554             xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:];           
    555             ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:];                     
     556            xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]           
     557            ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]
     558                                   
    556559           
    557560
     
    580583
    581584    from shallow_water_ext import balance_deep_and_shallow
    582     balance_deep_and_shallow(domain.number_of_elements,
    583                              wc, zc, wv, zv, hv,
     585    balance_deep_and_shallow(wc, zc, hc, wv, zv, hv,
    584586                             xmomc, ymomc, xmomv, ymomv)
    585587
     
    985987    gravity = gravity_c
    986988    manning_friction = manning_friction_c
    987     #balance_deep_and_shallow = balance_deep_and_shallow_c
     989    balance_deep_and_shallow = balance_deep_and_shallow_c
    988990   
    989991    #distribute_to_vertices_and_edges = distribute_to_vertices_and_edges_c
  • inundation/ga/storm_surge/pyvolution/shallow_water_ext.c

    r266 r267  
    176176           
    177177
    178 /*
    179 int _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.
     178
     179int _balance_deep_and_shallow(int N,
     180                              double* wc,
     181                              double* zc,
     182                              double* hc,                             
     183                              double* wv,
     184                              double* zv,
     185                              double* hv,
     186                              double* xmomc,
     187                              double* ymomc,
     188                              double* xmomv,
     189                              double* ymomv) {
     190 
     191  int k, k3, i;
     192  double dz, hmin, alpha;
     193 
     194  //Compute linear combination between constant levels and and
     195  //levels parallel to the bed elevation.     
     196 
     197  for (k=0; k<N; k++) {
     198    // Compute maximal variation in bed elevation
     199    // This quantitiy is
     200    //     dz = max_i abs(z_i - z_c)
     201    // and it is independent of dimension
     202    // In the 1d case zc = (z0+z1)/2
     203    // In the 2d case zc = (z0+z1+z2)/3
     204
     205    k3 = 3*k;
     206   
     207    //FIXME: Try with this one precomputed
     208    dz = 0.0;
     209    hmin = hv[k3];
     210    for (i=0; i<3; i++) {
     211      dz = max(dz, fabs(zv[k3+i]-zc[k]));
     212      hmin = min(hmin, hv[k3+i]);
     213    }
     214
     215   
     216    //Create alpha in [0,1], where alpha==0 means using shallow
     217    //first order scheme and alpha==1 means using the stage w as
     218    //computed by the gradient limiter (1st or 2nd order)
     219    //
     220    //If hmin > dz/2 then alpha = 1 and the bed will have no effect
     221    //If hmin < 0 then alpha = 0 reverting to constant height above bed.
     222   
     223    if (dz > 0.0)
     224      alpha = max( min( 2*hmin/dz, 1.0), 0.0 );
     225    else
     226      alpha = 1.0;  //Flat bed
     227
     228
     229    //Weighted balance between stage parallel to bed elevation
     230    //(wvi = zvi + hc) and stage as computed by 1st or 2nd
     231    //order gradient limiter
     232    //(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids
     233    //
     234    //It follows that the updated wvi is
     235    //  wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) =
     236    //  zvi + hc + alpha*(hvi - hc)
     237    //
     238    //Note that hvi = zc+hc-zvi in the first order case (constant).
     239
     240    if (alpha < 1) {         
     241      for (i=0; i<3; i++) {
     242        wv[k3+i] = zv[k3+i] + hc[k] + alpha*(hv[k3+i]-hc[k]);
     243           
    204244     
    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            
    235        
     245        //Update momentum as a linear combination of
     246        //xmomc and ymomc (shallow) and momentum
     247        //from extrapolator xmomv and ymomv (deep).
     248        xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];           
     249        ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
     250      }
     251    }
     252  }         
     253  return 0;
     254}
     255
    236256///////////////////////////////////////////////////////////////////
    237257// Gateways to Python
     
    542562  //
    543563  //    balance_deep_and_shallow(domain.number_of_elements,
    544   //                             wc, zc, wv, zv, hv,
     564  //                             wc, zc, hc, wv, zv, hv,
    545565  //                             xmomc, ymomc, xmomv, ymomv)
    546566 
     
    549569    *wc,            //Level at centroids
    550570    *zc,            //Elevation at centroids   
     571    *hc,            //Height at centroids       
    551572    *wv,            //Level at vertices
    552573    *zv,            //Elevation at vertices
     
    560581 
    561582  // Convert Python arguments to C 
    562   if (!PyArg_ParseTuple(args, "dOOOOOOOOO", &N, &wc, &zc, &wv, &zv, &hv,
     583  if (!PyArg_ParseTuple(args, "OOOOOOOOOO", 
     584                        &wc, &zc, &hc,
     585                        &wv, &zv, &hv,
    563586                        &xmomc, &ymomc, &xmomv, &ymomv))
    564587    return NULL;
    565588
    566  
    567   /*
     589  N = wc -> dimensions[0];
     590   
    568591  _balance_deep_and_shallow(N,
    569592                            (double*) wc -> data,
    570593                            (double*) zc -> data,
     594                            (double*) hc -> data,                           
    571595                            (double*) wv -> data,
    572596                            (double*) zv -> data,
     
    576600                            (double*) xmomv -> data,
    577601                            (double*) ymomv -> data); 
    578   */
     602 
    579603 
    580604  return Py_BuildValue(""); 
Note: See TracChangeset for help on using the changeset viewer.