Changeset 2632


Ignore:
Timestamp:
Mar 29, 2006, 4:33:30 PM (19 years ago)
Author:
ole
Message:

Work on creep demo

Location:
inundation
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • inundation/examples/island.py

    r2621 r2632  
    77
    88
    9 #-------------------------------------------------------------------------------
     9#------------------------------------------------------------------------------
    1010# Import necessary modules
    11 #-------------------------------------------------------------------------------
     11#------------------------------------------------------------------------------
    1212
    1313# Standard modules
     
    2121
    2222
    23 #-------------------------------------------------------------------------------   
     23#------------------------------------------------------------------------------
    2424# Setup computational domain
    25 #-------------------------------------------------------------------------------
     25#------------------------------------------------------------------------------
    2626
    2727#Create basic mesh
     
    3737                          maximum_triangle_area = 25,
    3838                          filename = 'island.msh',
    39                           interior_regions=[ ([[55,25], [75,25],
    40                                                [75,75], [55,75]], 3)]
     39                          interior_regions=[ ([[50,25], [70,25],
     40                                               [70,75], [50,75]], 3)]
    4141                          )
    4242                         
     
    6363
    6464
    65 #-------------------------------------------------------------------------------
     65#------------------------------------------------------------------------------
    6666# Setup initial conditions
    67 #-------------------------------------------------------------------------------
     67#------------------------------------------------------------------------------
    6868
    6969def island(x, y):
     
    7373    return z
    7474
    75 domain.set_quantity('friction', 0.0)
     75domain.set_quantity('friction', 2.0)
    7676domain.set_quantity('elevation', island)
    7777domain.set_quantity('stage', 1)
    7878
    7979
    80 #-------------------------------------------------------------------------------
     80#------------------------------------------------------------------------------
    8181# Setup boundary conditions (all reflective)
    82 #-------------------------------------------------------------------------------
     82#------------------------------------------------------------------------------
    8383
    8484Br = Reflective_boundary(domain)
     
    8888
    8989
    90 #-------------------------------------------------------------------------------
     90#------------------------------------------------------------------------------
    9191# Evolve system through time
    92 #-------------------------------------------------------------------------------
     92#------------------------------------------------------------------------------
    9393
    9494import time
    95 for t in domain.evolve(yieldstep = 0.5, finaltime = 100):
     95for t in domain.evolve(yieldstep = 1, finaltime = 100):
    9696    domain.write_time()
     97#for t in domain.evolve(yieldstep = 0.05, finaltime = 5):
     98#    domain.write_time()
     99#    #print '    Volume: ', domain.get_quantity('stage').get_integral()
    97100
    98101print 'Done'
  • inundation/examples/island_timeseries.py

    r2620 r2632  
    1010
    1111swwfile = 'island.sww'
    12 #gauges = [[0,0.5], [0.2,0.5], [0.4,0.5], [0.6,0.5], [0.8,0.5], [1,0.5]]
    13 #gauges = [[20, 50], [30, 50], [40, 50], [50, 50]]
    14 gauges = [[42, 48], [58, 48]]     
     12#gauges = [[56, 48], [58, 48], [60, 48], [62, 48], [64, 48]]
     13gauges = [[60, 48], [62, 48], [64, 48]]   
    1514
     15
     16store = False
    1617
    1718#Read model output
     
    103104    hold(False)
    104105
    105     if elevations[0] < -10:
     106    if elevations[0] < 0:
    106107        plot(model_time, stages, '-b')
    107108    else:   
    108109        plot(model_time, stages, '-b',
    109110             model_time, elevations, '-k')
    110     axis([0, 100, z-0.01, z+0.005])
     111    #axis([0, 100, z-0.01, z+0.005])
    111112       
    112113    #name = 'Gauge_%d: (%.1f, %.1f)' %(k, g[0], g[1])
     
    120121           shadow=True,
    121122           loc='upper right')
    122     #savefig('Gauge_%d_stage' %k)
    123     savefig('Gauge_%s_stage' %myloc)
     123
     124    if store is True: savefig('Gauge_%s_stage' %myloc)
    124125
    125126    raw_input('Next')
     
    136137    ylabel('sqrt( uh^2 + vh^2 ) [m^2/s]')   
    137138    #savefig('Gauge_%d_momentum' %k)
    138     savefig('Gauge_%s_momentum' %myloc)
     139    if store is True: savefig('Gauge_%s_momentum' %myloc)
    139140   
    140141    raw_input('Next')
     
    160161    ylabel(' atan(vh/uh) [degrees from North]')   
    161162    #savefig('Gauge_%d_bearing' %k)
    162     savefig('Gauge_%s_bearing' %myloc)
     163    if store is True: savefig('Gauge_%s_bearing' %myloc)
    163164   
    164165    raw_input('Next')
     
    176177    ylabel('sqrt( uh^2 + vh^2 ) / depth [m/s]')   
    177178    #savefig('Gauge_%d_speed' %k)
    178     savefig('Gauge_%s_speed' %myloc)
     179    if store is True: savefig('Gauge_%s_speed' %myloc)
    179180   
    180181    raw_input('Next')
  • inundation/pyvolution/shallow_water_ext.c

    r2621 r2632  
    246246      if (h >= eps) {
    247247        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
    248         //S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
    249         S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
     248        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
     249        //S /= exp(7.0/3.0*log(h));      //seems to save about 15% over manning_friction
    250250        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
    251251
     
    372372
    373373    if (hc < minimum_allowed_height) {
     374   
     375      //Old code: Set momentum to zero and ensure h is non negative
     376      /*
     377      xmomc[k] = 0.0;
     378      ymomc[k] = 0.0;     
     379      //wc[k] = zc[k];
     380      if (hc <= 0.0) wc[k] = zc[k];
     381      */
     382
     383      //New code: Adjust momentum to guarantee speeds are physical
     384      //          ensure h is non negative     
     385   
    374386      if (hc <= 0.0) {
    375         wc[k] = zc[k]; //Contain 'lost mass' error
     387        wc[k] = zc[k];
    376388        xmomc[k] = 0.0;
    377389        ymomc[k] = 0.0;
     
    383395          reduced_speed = maximum_allowed_speed * u/fabs(u);
    384396          //printf("Speed (u) has been reduced from %.3f to %.3f\n",
    385           //        u, reduced_speed); 
     397          //     u, reduced_speed);     
    386398          xmomc[k] = reduced_speed * hc;
    387399        }
     
    391403          reduced_speed = maximum_allowed_speed * v/fabs(v);   
    392404          //printf("Speed (v) has been reduced from %.3f to %.3f\n",   
    393           //        v, reduced_speed); 
     405          //     v, reduced_speed);     
    394406          ymomc[k] = reduced_speed * hc;         
    395407        }       
    396408      }
    397      
    398       //Old code
    399       //wc[k] = zc[k]; //Contain 'lost mass' error
    400       //printf("Speed has been reduced to %.3f\n", 0.0);             
    401       //xmomc[k] = 0.0;
    402       //ymomc[k] = 0.0;     
    403409    }
    404410  }
Note: See TracChangeset for help on using the changeset viewer.