Changeset 2566


Ignore:
Timestamp:
Mar 21, 2006, 2:25:40 PM (18 years ago)
Author:
ole
Message:

Created maximum_allowed_speed in config and passed it through.
Also ditched epsilon and used h<=0 instead.

Location:
inundation
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • inundation/pyvolution/config.py

    r2528 r2566  
    3333
    3434
    35 v_max = 100 #For use in domain_ext.c
     35#v_max = 100 #For use in domain_ext.c
    3636sound_speed = 500
    3737
     
    106106#Specific to shallow water W.E.
    107107minimum_allowed_height = 1.0e-3 #Water depth below which it is considered to be 0
     108maximum_allowed_speed = 100.0 #Maximal particle speed of water
  • inundation/pyvolution/shallow_water.py

    r2516 r2566  
    7979                                tagged_elements, geo_reference, use_inscribed_circle)
    8080
    81         from config import minimum_allowed_height, g
     81        from config import minimum_allowed_height, maximum_allowed_speed, g
    8282        self.minimum_allowed_height = minimum_allowed_height
     83        self.maximum_allowed_speed = maximum_allowed_speed     
    8384        self.g = g
    8485
     
    667668    hc = wc - zc  #Water depths at centroids
    668669
    669     #Update
     670    #Update
     671    #FIXME: Modify accroditg to c-version - or discard altogether.
    670672    for k in range(domain.number_of_elements):
    671673
     
    691693    from shallow_water_ext import protect
    692694
    693     protect(domain.minimum_allowed_height, domain.epsilon,
    694             wc, zc, xmomc, ymomc)
     695    protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
     696            domain.epsilon, wc, zc, xmomc, ymomc)
    695697
    696698
  • inundation/pyvolution/shallow_water_ext.c

    r2561 r2566  
    356356int _protect(int N,
    357357             double minimum_allowed_height,
     358             double maximum_allowed_speed,           
    358359             double epsilon,
    359360             double* wc,
     
    364365  int k;
    365366  double hc;
    366   double vx, vy; 
    367   double maximum_allowed_speed = 100.0;  //FIXME (Ole): Pass in
     367  double u, v, tmp_speed; 
    368368
    369369  //Protect against initesimal and negative heights
     
    373373
    374374    if (hc < minimum_allowed_height) {
    375       if (hc < epsilon) {
     375    //if (1) {
     376      //if (hc < epsilon) {
     377      if (hc <= 0.0) {
    376378        wc[k] = zc[k]; //Contain 'lost mass' error
    377379        xmomc[k] = 0.0;
     
    380382        //Try to see if calculated speeds would blow up
    381383        //FIXME (Ole): This has not been written in Python
    382         vx = xmomc[k]/hc;
    383         vy = ymomc[k]/hc;       
    384         if (fabs(vx) > maximum_allowed_speed) {
    385           printf("WARNING: Speed has %.3f has been reduced to %.3f\n",
    386                  vx, maximum_allowed_speed);   
    387           xmomc[k] = maximum_allowed_speed * hc * vx/fabs(vx);
     384        u = xmomc[k]/hc;
     385        if (fabs(u) > maximum_allowed_speed) {
     386          tmp_speed = maximum_allowed_speed * u/fabs(u);
     387          //printf("WARNING: Speed (u) has been reduced from %.3f to %.3f\n",
     388          //     u, tmp_speed);
     389          xmomc[k] = tmp_speed * hc;
    388390        }
    389391       
    390         if (fabs(vy) > maximum_allowed_speed) {
    391           printf("WARNING: Speed has %.3f has been reduced to %.3f\n",
    392                  vy, maximum_allowed_speed);   
    393           ymomc[k] = maximum_allowed_speed * hc * vy/fabs(vy);
     392        v = ymomc[k]/hc;       
     393        if (fabs(v) > maximum_allowed_speed) {
     394          tmp_speed = maximum_allowed_speed * v/fabs(v);       
     395          //printf("WARNING: Speed (v) has been reduced from %.3f to %.3f\n",   
     396          //     v, maximum_allowed_speed);     
     397          ymomc[k] = tmp_speed * hc;     
    394398        }       
    395399       
     
    10431047PyObject *protect(PyObject *self, PyObject *args) {
    10441048  //
    1045   //    protect(minimum_allowed_height, wc, zc, xmomc, ymomc)
     1049  //    protect(minimum_allowed_height, maximum_allowed_speed, wc, zc, xmomc, ymomc)
    10461050
    10471051
     
    10541058
    10551059  int N;
    1056   double minimum_allowed_height, epsilon;
     1060  double minimum_allowed_height, maximum_allowed_speed, epsilon;
    10571061
    10581062  // Convert Python arguments to C
    1059   if (!PyArg_ParseTuple(args, "ddOOOO",
     1063  if (!PyArg_ParseTuple(args, "dddOOOO",
    10601064                        &minimum_allowed_height,
     1065                        &maximum_allowed_speed,                 
    10611066                        &epsilon,
    10621067                        &wc, &zc, &xmomc, &ymomc))
     
    10671072  _protect(N,
    10681073           minimum_allowed_height,
     1074           maximum_allowed_speed,
    10691075           epsilon,
    10701076           (double*) wc -> data,
  • inundation/pyvolution/test_shallow_water.py

    r2565 r2566  
    21742174        os.remove(domain.filename + '.sww')
    21752175
     2176       
     2177    def test_flatbed_second_order_vmax_0(self):
     2178        from mesh_factory import rectangular
     2179        from shallow_water import Domain,\
     2180             Reflective_boundary, Dirichlet_boundary
     2181
     2182        from Numeric import array
     2183
     2184        #Create basic mesh
     2185        N = 8
     2186        points, vertices, boundary = rectangular(N, N)
     2187
     2188        #Create shallow water domain
     2189        domain = Domain(points, vertices, boundary)
     2190        domain.smooth = False
     2191        domain.visualise = False
     2192        domain.default_order=2
     2193        domain.maximum_allowed_speed = 0.0 #Makes it like the 'oldstyle'
     2194
     2195        # Boundary conditions
     2196        Br = Reflective_boundary(domain)
     2197        Bd = Dirichlet_boundary([0.2,0.,0.])
     2198
     2199        domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br})
     2200        domain.check_integrity()
     2201
     2202        #Evolution
     2203        for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03):
     2204            pass
     2205
     2206
     2207        assert allclose(domain.min_timestep, 0.0210448446782)
     2208        assert allclose(domain.max_timestep, 0.0210448446782)
     2209
     2210        #FIXME: These numbers were from version before 21/3/6 -
     2211        #could be recreated by setting maximum_allowed_speed to 0 maybe
     2212        assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
     2213                        [ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
     2214       
     2215
     2216        assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
     2217                        [-0.00139463,0.0006156,-0.00060364,0.00061827])
     2218
     2219
     2220        os.remove(domain.filename + '.sww')
     2221
     2222       
    21762223
    21772224    def test_flatbed_second_order_distribute(self):
  • inundation/test_all.py

    r2563 r2566  
    2323exclude_dirs = ['pypar_dist', #Special requirements
    2424                'props', 'wcprops', 'prop-base', 'text-base', '.svn', #Svn
    25                 'tmp']
    26                 #'pmesh']     #Name conflict on my home machine (Ole)       
     25                'tmp',
     26                'pmesh']     #Name conflict on my home machine (Ole)       
    2727
    2828
Note: See TracChangeset for help on using the changeset viewer.