Changeset 6694


Ignore:
Timestamp:
Apr 1, 2009, 8:41:37 PM (16 years ago)
Author:
sudi
Message:

Revised codes for discontinuous bed.

Location:
anuga_work/development/anuga_1d
Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c

    r5832 r6694  
    33#include "math.h"
    44#include <stdio.h>
     5#include "util_ext.h"
    56const double pi = 3.14159265358979;
    67
    7 
    8 // Shared code snippets
    9 #include "util_ext.h"
    10 
    11 
    12 //double max(double a, double b) {
    13 //      double z;
    14 //      z=(a>b)?a:b;
    15 //      return z;}
    16        
    17 //double min(double a, double b) {
    18 //      double z;
    19 //      z=(a<b)?a:b;
    20 //      return z;}
    21 
    22 
    23 
    24 
    25 
    268//Innermost flux function (using w=z+h)
    27 int _flux_function_wellbalanced(double *q_left, double *q_right,
    28                 double normals, double g, double epsilon,
    29                 double *edgeflux, double *max_speed) {
    30                
     9int _flux_function_wellbalanced(double *q_left,
     10                                                                double *q_right,
     11                                                                double normals,
     12                                                                double g,
     13                                                                double epsilon,
     14                                                                double *edgeflux,
     15                                                                double *max_speed) {           
    3116                int i;
    3217        double z_left, z_right;
    3318                double flux_left[2], flux_right[2];
    3419                double z_star, w_left, h_left, h_left_star, uh_left, soundspeed_left, u_left;
    35                 double w_right, h_right, h_right_star, uh_right, soundspeed_right, u_right;                          //h_right,
     20                double w_right, h_right, h_right_star, uh_right, soundspeed_right, u_right;
    3621                double s_max, s_min, denom;
    3722               
    38                 w_left = q_left[0];
     23                w_left  = q_left[0];
    3924                uh_left = q_left[1];
    4025                uh_left = uh_left*normals;
    41         z_left = q_left[2];
    42        
    43                 w_right = q_right[0];
     26        z_left  = q_left[2];
     27                h_left  = q_left[3];
     28                u_left  = q_left[4];
     29       
     30                w_right  = q_right[0];
    4431                uh_right = q_right[1];
    4532                uh_right = uh_right*normals;
    46         z_right = q_right[2];
     33        z_right  = q_right[2];
     34                h_right  = q_right[3];
     35                u_right  = q_right[4];
    4736         
    4837                if (z_left-z_right != 0.0 ) printf("z_l - z_r = %f \n",z_left-z_right);
    4938       
    50                 //z = (z_left+z_right)/2.0;
    5139                z_star = max(z_left, z_right);                                                                                                          //equation(7) 
    5240                                   
    53                 // Compute speeds in x-direction
    54                 h_left = w_left-z_left;                                                                         //This is used for computing the edgeflux[1].
    55                 h_left_star = max(0, w_left-z_star);                                                                                                            //(8)
    56    
     41                // Compute speeds in x-direction.
     42                h_left_star = max(0.0, w_left-z_star);                                                                                                          //equation(8)
     43               
    5744                if (h_left_star < epsilon) {
    5845                        u_left = 0.0; h_left_star = 0.0;
     
    6148                        u_left = uh_left/h_left_star;
    6249                }
    63        
    64                 h_right = w_right-z_right;                                                                                                                                                                                                     
    65                 h_right_star = max(0, w_right-z_star);                                                                                                          //(9)
     50                                                                                                                                                                                                   
     51                h_right_star = max(0.0, w_right-z_star);                                                                                                        //equation(9)
    6652
    6753                if (h_right_star < epsilon) {
     
    9278                // Flux computation
    9379                denom = s_max-s_min;
    94                 if (denom < epsilon) {
     80                if (denom < epsilon && z_right > z_left) {
    9581                        for (i=0; i<2; i++) edgeflux[i] = 0.0;
     82                        // Maximal wavespeed                   
     83                        edgeflux[1] = edgeflux[1] + (0.5*g*h_left*h_left - 0.5*g*h_left_star*h_left_star );
     84                        edgeflux[1] *= normals;
    9685                        *max_speed = 0.0;
    97                 } else {
     86                } else if (denom < epsilon){
     87                        for (i=0; i<2; i++) edgeflux[i] = 0.0;
     88                        // Maximal wavespeed
     89                        *max_speed = 0.0;               
     90                }else {
    9891                        edgeflux[0] = s_max*flux_left[0] - s_min*flux_right[0];
    9992                        edgeflux[0] += s_max*s_min*(h_right_star-h_left_star);                   //s_max*s_min*(qr[0]-ql[0]);
     
    10295                        edgeflux[1] += s_max*s_min*(flux_right[0]-flux_left[0]);                 //s_max*s_min*(qr[1]-ql[1]);
    10396                        edgeflux[1] /= denom;
     97                        edgeflux[1] = edgeflux[1] + (0.5*g*h_left*h_left - 0.5*g*h_left_star*h_left_star );
    10498                        edgeflux[1] *= normals;
    105    
    106                 //edgeflux[1] = edgeflux[1] + 0.5*g*0.5*(h_left+h_right)*0.5*(h_left+h_right) - 0.5*g*h_left_star*h_left_star;
    107                 //edgeflux[1] = edgeflux[1] + 0.5*g*h_right*h_right - 0.5*g*h_right_star*h_right_star;
    108                 edgeflux[1] = edgeflux[1] + 0.5*g*h_left*h_left - 0.5*g*h_left_star*h_left_star;           //what about the h_right and h_right_star??????????????????????
    109                 // Maximal wavespeed
    110         *max_speed = max(fabs(s_max), fabs(s_min));
     99                        // Maximal wavespeed
     100                        *max_speed = max(fabs(s_max), fabs(s_min));
    111101                } 
    112                 //printf("max_speed = %g \n",*max_speed);
     102               
    113103    return 0;           
    114104        }
     
    121111                double epsilon,
    122112                double g,
     113               
    123114                long* neighbours,
    124115                long* neighbour_vertices,
    125116                double* normals,
    126117                double* areas,
     118               
    127119                double* stage_edge_values,
    128120                double* xmom_edge_values,
    129121                double* bed_edge_values,
     122                double* height_edge_values,
     123                double* vel_edge_values,
     124               
    130125                double* stage_boundary_values,
    131126                double* xmom_boundary_values,
     127                double* bed_boundary_values,
     128                double* height_boundary_values,
     129                double* vel_boundary_values,
     130               
    132131                double* stage_explicit_update,
    133132                double* xmom_explicit_update,
     133                double* bed_explicit_update,
     134                double* height_explicit_update,
     135                double* vel_explicit_update,
     136               
    134137                int number_of_elements,
    135138                double* max_speed_array) {
    136139               
    137                 double flux[2], ql[3], qr[3], edgeflux[2];
     140                double flux[2], ql[5], qr[5], edgeflux[2];
    138141                double max_speed, normal;
    139142                int k, i, ki, n, m, nm=0;
     
    149152                                ql[1] = xmom_edge_values[ki];
    150153                                ql[2] = bed_edge_values[ki];
    151                                 //ql[3] = velocity_edge_values[ki];
    152                 //ql[4] = height_edge_values[ki];
     154                ql[3] = height_edge_values[ki];
     155                                ql[4] = vel_edge_values[ki];
    153156               
    154157                                n = neighbours[ki];
     158                                printf("neighbours[ki] = %li \n",neighbours[ki]);
    155159                                if (n<0) {
    156160                                        m = -n-1;
    157161                                        qr[0] = stage_boundary_values[m];
    158162                                        qr[1] = xmom_boundary_values[m];
    159                                         qr[2] = ql[2];
     163                                        qr[2] = ql[2];                                                                                                                                                 
     164                                        qr[3] = height_edge_values[m];                                                                                                         
     165                                        qr[4] = vel_edge_values[m];                                                                                                                             
    160166                   
    161167                                } else {
     
    164170                                        qr[0] = stage_edge_values[nm];
    165171                                        qr[1] = xmom_edge_values[nm];
    166                                         qr[2] = bed_edge_values[nm];                           
     172                                        qr[2] = bed_edge_values[nm];
     173                                        qr[3] = height_edge_values[nm];                                                                                                                 
     174                                        qr[4] = vel_edge_values[nm];                                                                                                                   
    167175                                }
    168176                               
     
    175183                                if (max_speed > epsilon) {
    176184                        // Original CFL calculation
    177                                        
    178                                         timestep = min(timestep, 0.5*areas[k]/max_speed); //Here, CFL=1.0 is assumed. ?????????????????????????????????????????????
     185                                        timestep = min(timestep, 0.5*areas[k]/max_speed); //Here, CFL=1.0 is assumed.
    179186                                        if (n>=0) {
    180                                                 timestep = min(timestep, 0.5*areas[n]/max_speed); //Here, CFL=1.0 is assumed. ?????????????????????????????????????????????
     187                                                timestep = min(timestep, 0.5*areas[n]/max_speed); //Here, CFL=1.0 is assumed.
    181188                                        }
    182189                                }
     
    190197                        max_speed_array[k]=max_speed;
    191198                }
    192                 //printf("timestep = %f \n",timestep);
     199               
    193200                return timestep;       
    194201        }
     
    209216        *normals,
    210217        *areas,
     218       
    211219        *stage_edge_values,
    212220        *xmom_edge_values,
    213221        *bed_edge_values,
     222        *height_edge_values,
     223        *vel_edge_values,
     224       
    214225        *stage_boundary_values,
    215226        *xmom_boundary_values,
     227        *bed_boundary_values,
     228        *height_boundary_values,
     229        *vel_boundary_values,
     230       
    216231        *stage_explicit_update,
    217232        *xmom_explicit_update,
     233        *bed_explicit_update,
     234        *height_explicit_update,
     235        *vel_explicit_update,
     236       
    218237        *max_speed_array;
    219238   
     
    222241   
    223242  // Convert Python arguments to C
    224   if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOiO",
     243  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOOOOiO",
    225244                        &timestep,
    226245                        &epsilon,
    227246                        &g,
     247                       
    228248                        &neighbours,
    229249                        &neighbour_vertices,
    230250                        &normals,
    231251                        &areas,
     252                       
    232253                        &stage_edge_values,
    233254                        &xmom_edge_values,
    234255                        &bed_edge_values,
     256                        &height_edge_values,
     257                        &vel_edge_values,
     258                       
    235259                        &stage_boundary_values,
    236260                        &xmom_boundary_values,
     261                        &bed_boundary_values,
     262                        &height_boundary_values,
     263                        &vel_boundary_values,
     264                       
    237265                        &stage_explicit_update,
    238266                        &xmom_explicit_update,
     267                        &bed_explicit_update,
     268                        &height_explicit_update,
     269                        &vel_explicit_update,
     270                       
    239271                        &number_of_elements,
    240272                        &max_speed_array)) {
    241     PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext could not parse input");
     273    PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext_wellbalanced.c: compute_fluxes_ext_wellbalanced could not parse input");
    242274    return NULL;
    243275  }
     
    249281                                 epsilon,
    250282                                 g,
     283                                 
    251284                                 (long*) neighbours -> data,
    252285                                 (long*) neighbour_vertices -> data,
    253286                                 (double*) normals -> data,
    254287                                 (double*) areas -> data,
     288                                 
    255289                                 (double*) stage_edge_values -> data,
    256290                                 (double*) xmom_edge_values -> data,
    257291                                 (double*) bed_edge_values -> data,
     292                                 (double*) height_edge_values -> data,
     293                                 (double*) vel_edge_values -> data,
     294                                 
    258295                                 (double*) stage_boundary_values -> data,
    259296                                 (double*) xmom_boundary_values -> data,
     297                                 (double*) bed_boundary_values -> data,
     298                                 (double*) height_boundary_values -> data,
     299                                 (double*) vel_boundary_values -> data,
     300                                 
    260301                                 (double*) stage_explicit_update -> data,
    261302                                 (double*) xmom_explicit_update -> data,
     303                                 (double*) bed_explicit_update -> data,
     304                                 (double*) height_explicit_update -> data,
     305                                 (double*) vel_explicit_update -> data,
     306                                 
    262307                                 number_of_elements,
    263308                                 (double*) max_speed_array -> data);
     
    266311  return Py_BuildValue("d", timestep);
    267312}
    268 
    269 
    270313
    271314
     
    285328  import_array();
    286329}
    287        
    288        
  • anuga_work/development/anuga_1d/config.py

    r5827 r6694  
    109109
    110110#Specific to shallow water W.E.
    111 minimum_allowed_height = 1.0e-6 #Water depth below which it is considered to be 0
     111h_min = minimum_allowed_height = 1.0e-6 #1.0e-6 #Water depth below which it is considered to be 0
    112112maximum_allowed_speed = 1000.0
  • anuga_work/development/anuga_1d/domain.py

    r6453 r6694  
    916916        # Update ghosts
    917917        self.update_ghosts()
    918        
     918
    919919        # Initial update of vertex and edge values
    920         self.distribute_to_vertices_and_edges() 
    921        
     920        self.distribute_to_vertices_and_edges()
     921
    922922        # Update extrema if necessary (for reporting)
    923923        self.update_extrema()
     
    15021502                Q = self.quantities[name]
    15031503                Q.boundary_values[i] = q[j]
     1504                #print 'Q=',Q
    15041505   
    15051506    def update_timestep(self, yieldstep, finaltime):
  • anuga_work/development/anuga_1d/shallow_water_domain_suggestion2.py

    r5827 r6694  
    5353
    5454        conserved_quantities = ['stage', 'xmomentum']
    55         other_quantities = ['elevation', 'friction', 'height', 'velocity']
     55        evolved_quantities = ['stage', 'xmomentum', 'elevation', 'height', 'velocity']
     56        other_quantities = ['friction']
    5657        Generic_Domain.__init__(self, coordinates, boundary,
    57                                 conserved_quantities, other_quantities,
     58                                conserved_quantities, evolved_quantities, other_quantities,
    5859                                tagged_elements)
    5960       
     
    9192        self.quantities_to_be_stored = ['stage','xmomentum']
    9293
    93         self.__doc__ = 'shallow_water_domain'
     94        self.__doc__ = 'shallow_water_domain_suggestion2'
     95        self.check_integrity()
    9496
    9597
     
    148150        msg = 'Second conserved quantity must be "xmomentum"'
    149151        assert self.conserved_quantities[1] == 'xmomentum', msg
     152               
     153        msg = 'First evolved quantity must be "stage"'
     154        assert self.evolved_quantities[0] == 'stage', msg
     155        msg = 'Second evolved quantity must be "xmomentum"'
     156        assert self.evolved_quantities[1] == 'xmomentum', msg
     157        msg = 'Third evolved quantity must be "elevation"'
     158        assert self.evolved_quantities[2] == 'elevation', msg
     159        msg = 'Fourth evolved quantity must be "height"'
     160        assert self.evolved_quantities[3] == 'height', msg
     161        msg = 'Fifth evolved quantity must be "velocity"'
     162        assert self.evolved_quantities[4] == 'velocity', msg
    150163
    151164    def extrapolate_second_order_sw(self):
     
    168181        distribute_to_vertices_and_edges(self)
    169182       
    170     def evolve(self, yieldstep = None, finaltime = None, duration = None,
    171                skip_initial_step = False):
    172         """Specialisation of basic evolve method from parent class
    173         """
    174 
    175         #Call check integrity here rather than from user scripts
    176         #self.check_integrity()
    177 
    178         #msg = 'Parameter beta_h must be in the interval [0, 1)'
    179         #assert 0 <= self.beta_h < 1.0, msg
    180         #msg = 'Parameter beta_w must be in the interval [0, 1)'
    181         #assert 0 <= self.beta_w < 1.0, msg
    182 
    183 
    184         #Initial update of vertex and edge values before any storage
    185         #and or visualisation
    186        
    187         #self.distribute_to_vertices_and_edges() ?????????????????????????????????
    188        
    189         #Initialise real time viz if requested
    190         #if self.visualise is True and self.time == 0.0:
    191         #    if self.visualiser is None:
    192         #        self.initialise_visualiser()
    193         #
    194         #    self.visualiser.update_timer()
    195         #    self.visualiser.setup_all()
    196 
    197         #Store model data, e.g. for visualisation
    198         #if self.store is True and self.time == 0.0:
    199         #    self.initialise_storage()
    200         #    #print 'Storing results in ' + self.writer.filename
    201         #else:
    202         #    pass
    203         #    #print 'Results will not be stored.'
    204         #    #print 'To store results set domain.store = True'
    205         #    #FIXME: Diagnostic output should be controlled by
    206         #    # a 'verbose' flag living in domain (or in a parent class)
    207 
     183    def evolve(self, yieldstep = None, finaltime = None, duration = None, skip_initial_step = False):
    208184        #Call basic machinery from parent class
    209         for t in Generic_Domain.evolve(self, yieldstep, finaltime,duration,
    210                                        skip_initial_step):
    211             #Real time viz
    212         #    if self.visualise is True:
    213         #        self.visualiser.update_all()
    214         #        self.visualiser.update_timer()
    215 
    216 
    217             #Store model data, e.g. for subsequent visualisation
    218         #    if self.store is True:
    219         #        self.store_timestep(self.quantities_to_be_stored)
    220 
    221             #FIXME: Could maybe be taken from specified list
    222             #of 'store every step' quantities
    223 
    224             #Pass control on to outer loop for more specific actions
     185        for t in Generic_Domain.evolve(self, yieldstep, finaltime, duration, skip_initial_step):
    225186            yield(t)
    226187
     
    352313       
    353314
    354     #We have got   h   and    u   at vertex, then  the following is the calculation of fluxes!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     315    #We have got h and u at vertex, then the following is the calculation of fluxes!
    355316    soundspeed_left  = sqrt(g*h_left)
    356317    soundspeed_right = sqrt(g*h_right)
     
    656617    import sys
    657618    from Numeric import zeros, Float
    658     from config import epsilon, h0
    659619   
    660620    N = domain.number_of_elements
    661621    timestep = float(sys.maxint)
    662     #epsilon = domain.epsilon
     622    epsilon = domain.epsilon
    663623    g = domain.g
    664624    neighbours = domain.neighbours
     
    667627    areas = domain.areas
    668628   
     629    Stage    = domain.quantities['stage']
     630    Xmom     = domain.quantities['xmomentum']
     631    Bed      = domain.quantities['elevation']
     632    Height   = domain.quantities['height']
     633    Velocity = domain.quantities['velocity']
     634
     635
     636    stage_boundary_values = Stage.boundary_values
     637    xmom_boundary_values  = Xmom.boundary_values
     638    bed_boundary_values   = Bed.boundary_values
     639    height_boundary_values= Height.boundary_values
     640    vel_boundary_values   = Velocity.boundary_values
     641   
     642    stage_explicit_update = Stage.explicit_update
     643    xmom_explicit_update  = Xmom.explicit_update
     644    bed_explicit_values   = Bed.explicit_update
     645    height_explicit_values= Height.explicit_update
     646    vel_explicit_values   = Velocity.explicit_update
     647   
     648    max_speed_array = domain.max_speed_array
     649   
     650    domain.distribute_to_vertices_and_edges()
     651    domain.update_boundary()
     652   
     653    stage_V = Stage.vertex_values
     654    xmom_V  = Xmom.vertex_values
     655    bed_V   = Bed.vertex_values
     656    height_V= Height.vertex_values
     657    vel_V   = Velocity.vertex_values
     658
     659    number_of_elements = len(stage_V)
     660
     661    #from config import epsilon, h0, h_min
     662    # ##Check this:
     663    #for i in range(N):
     664    #    height_V[i] = stage_V[i] - bed_V[i]                                           
     665    #    if height_V[i] <= h_min:  #epsilon :
     666    #       height_V[i] = 0.0
     667    #       stage_V[i]  = bed_V[i]
     668    #       xmom_V[i]   = 0.0
     669    #       vel_V[i]    = 0.0
     670    #   else:
     671    #        vel_V[i]  = xmom_V[i]/(height_V[i] +  h0/height_V[i])
     672    # ##End of part to be checked.
     673
     674   
     675    #print 'neighbours=',neighbours
     676    #print 'neighbour_vertices',neighbour_vertices
     677    #print 'normals=',normals
     678    #print 'areas=',areas
     679    from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced
     680    domain.flux_timestep = compute_fluxes_ext_wellbalanced(timestep,
     681                                  epsilon,
     682                                  g,
     683                                                           
     684                                  neighbours,
     685                                  neighbour_vertices,
     686                                  normals,
     687                                  areas,
     688                                                           
     689                                  stage_V,
     690                                  xmom_V,
     691                                  bed_V,
     692                                  height_V,
     693                                  vel_V,
     694                                                           
     695                                  stage_boundary_values,
     696                                  xmom_boundary_values,
     697                                  bed_boundary_values,
     698                                  height_boundary_values,
     699                                  vel_boundary_values,
     700                                                           
     701                                  stage_explicit_update,
     702                                  xmom_explicit_update,
     703                                  bed_explicit_values,
     704                                  height_explicit_values,
     705                                  vel_explicit_values,
     706                                                           
     707                                  number_of_elements,
     708                                  max_speed_array)
     709 
     710
     711# ###################################
     712
     713
     714
     715
     716
     717
     718# Module functions for gradient limiting (distribute_to_vertices_and_edges)
     719
     720def distribute_to_vertices_and_edges(domain):
     721    """Distribution from centroids to vertices specific to the
     722    shallow water wave
     723    equation.
     724
     725    It will ensure that h (w-z) is always non-negative even in the
     726    presence of steep bed-slopes by taking a weighted average between shallow
     727    and deep cases.
     728
     729    In addition, all conserved quantities get distributed as per either a
     730    constant (order==1) or a piecewise linear function (order==2).
     731
     732    FIXME: more explanation about removal of artificial variability etc
     733
     734    Precondition:
     735      All quantities defined at centroids and bed elevation defined at
     736      vertices.
     737
     738    Postcondition
     739      Conserved quantities defined at vertices
     740
     741    """
     742
     743    #from config import optimised_gradient_limiter
     744
     745    #Remove very thin layers of water
     746    #protect_against_infinitesimal_and_negative_heights(domain) 
     747
     748    import sys
     749    from Numeric import zeros, Float
     750    from config import epsilon, h0, h_min
     751
     752    N = domain.number_of_elements
     753
     754    #Shortcuts
    669755    Stage = domain.quantities['stage']
    670     Xmom  = domain.quantities['xmomentum']
    671     Bed   = domain.quantities['elevation']
     756    Xmom = domain.quantities['xmomentum']
     757    Bed = domain.quantities['elevation']
    672758    Height = domain.quantities['height']
    673759    Velocity = domain.quantities['velocity']
    674760
     761    #Arrays   
    675762    w_C   = Stage.centroid_values   
    676763    uh_C  = Xmom.centroid_values   
    677     z_C   = Bed.centroid_values
     764    z_C   = Bed.centroid_values 
    678765    h_C   = Height.centroid_values
    679766    u_C   = Velocity.centroid_values
     
    681768    for i in range(N):
    682769        h_C[i] = w_C[i] - z_C[i]                                               
    683         if h_C[i] < epsilon :
     770        if h_C[i] <= h_min:  #epsilon :
    684771            h_C[i] = 0.0
    685772            w_C[i] = z_C[i]
    686773            uh_C[i] = 0.0
    687 
    688     for i in range(len(h_C)):
    689         if h_C[i] < epsilon:
    690             u_C[i] = 0.0  #Could have been negative
    691             h_C[i] = 0.0
     774            u_C[i] = 0.0
    692775        else:
    693             u_C[i]  = uh_C[i]/(h_C[i] +  h0/h_C[i])
    694 
    695     for name in ['height', 'velocity']:
     776            u_C[i]  = uh_C[i]/(h_C[i] +  h0/h_C[i])
     777
     778    for name in ['stage', 'height', 'velocity']:
    696779        Q = domain.quantities[name]
    697780        if domain.order == 1:
     
    705788            raise 'Unknown order'
    706789
    707     w_V = domain.quantities['stage'].vertex_values                     
    708     h_V = domain.quantities['height'].vertex_values
    709     u_V = domain.quantities['velocity'].vertex_values           
    710     uh_V = domain.quantities['xmomentum'].vertex_values
    711     z_V   = w_V - h_V
    712 
     790    stage_V = domain.quantities['stage'].vertex_values
     791    bed_V   = domain.quantities['elevation'].vertex_values
     792    h_V     = domain.quantities['height'].vertex_values
     793    u_V     = domain.quantities['velocity'].vertex_values               
     794    xmom_V  = domain.quantities['xmomentum'].vertex_values
     795
     796    bed_V[:] = stage_V - h_V
     797    xmom_V[:] = u_V * h_V
    713798   
    714     stage_boundary_values = Stage.boundary_values
    715     xmom_boundary_values =  Xmom.boundary_values
    716     stage_explicit_update = Stage.explicit_update
    717     xmom_explicit_update = Xmom.explicit_update
    718     max_speed_array = domain.max_speed_array
    719     #number_of_elements = len(w_V)
     799    return
     800
    720801   
    721     #domain.distribute_to_vertices_and_edges()
    722     #domain.update_boundary()
    723     #stage_V = Stage.vertex_values
    724     #xmom_V  = Xmom.vertex_values
    725     #bed_V   = Bed.vertex_values
    726    
    727        
    728     from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced  #from comp_flux_ext import compute_fluxes_ext
    729        
    730     domain.flux_timestep = compute_fluxes_ext_wellbalanced(timestep,
    731                                   epsilon,
    732                                   g,
    733                                   neighbours,
    734                                   neighbour_vertices,
    735                                   normals,
    736                                   areas,
    737                                   w_V,
    738                                   uh_V,
    739                                   z_V,
    740                                   stage_boundary_values,
    741                                   xmom_boundary_values,
    742                                   stage_explicit_update,
    743                                   xmom_explicit_update,
    744                                   N,
    745                                   max_speed_array)
    746 
    747 # ###################################
    748 
    749 
    750 
    751 
    752 
    753 
    754 # Module functions for gradient limiting (distribute_to_vertices_and_edges)
    755 
    756 def distribute_to_vertices_and_edges(domain):
    757     """Distribution from centroids to vertices specific to the
    758     shallow water wave
    759     equation.
    760 
    761     It will ensure that h (w-z) is always non-negative even in the
    762     presence of steep bed-slopes by taking a weighted average between shallow
    763     and deep cases.
    764 
    765     In addition, all conserved quantities get distributed as per either a
    766     constant (order==1) or a piecewise linear function (order==2).
    767 
    768     FIXME: more explanation about removal of artificial variability etc
    769 
    770     Precondition:
    771       All quantities defined at centroids and bed elevation defined at
    772       vertices.
    773 
    774     Postcondition
    775       Conserved quantities defined at vertices
    776 
    777     """
    778 
    779     #from config import optimised_gradient_limiter
    780 
    781     #Remove very thin layers of water
    782     protect_against_infinitesimal_and_negative_heights(domain) 
    783 
    784     import sys
    785     from Numeric import zeros, Float
    786 
    787     N = domain.number_of_elements
    788 
    789     #Shortcuts
    790     Stage = domain.quantities['stage']
    791     Xmom = domain.quantities['xmomentum']
    792     Bed = domain.quantities['elevation']
    793     Height = domain.quantities['height']
    794     Velocity = domain.quantities['velocity']
    795 
    796     #Arrays   
    797     w_C   = Stage.centroid_values   
    798     uh_C  = Xmom.centroid_values   
    799     z_C   = Bed.centroid_values
    800     h_C   = Height.centroid_values
    801     u_C   = Velocity.centroid_values
    802        
    803     #print id(h_C)
    804     for i in range(N):
    805         h_C[i] = w_C[i] - z_C[i]                                                #This is h at centroids:  OK
    806    
    807     from config import epsilon, h0
    808        
    809     for i in range(len(h_C)):
    810         if h_C[i] < epsilon:
    811             u_C[i] = 0.0  #Could have been negative
    812             h_C[i] = 0.0
    813         else:
    814             u_C[i]  = uh_C[i]/(h_C[i] +  h0/h_C[i])
    815        
    816     for name in ['stage', 'velocity', 'height' ]:
    817         Q = domain.quantities[name]
    818         if domain.order == 1:
    819             Q.extrapolate_first_order()
    820         elif domain.order == 2:
    821             #print "add extrapolate second order to shallow water"
    822             #if name != 'height':
    823             Q.extrapolate_second_order()
    824             #Q.limit()
    825         else:
    826             raise 'Unknown order'
    827 
    828     stage  = domain.quantities['stage'].vertex_values                   #This is w at vertices:   OK
    829     h_V    = domain.quantities['height'].vertex_values
    830     bed    = stage - h_V
    831     u_V    = domain.quantities['velocity'].vertex_values               
    832     xmom_V = domain.quantities['xmomentum'].vertex_values       
    833        
    834     #h_V[:]    = stage - bed
    835     xmom_V[:] = u_V * h_V
    836     return stage, xmom_V, bed, h_V, u_V
    837802#
    838803def protect_against_infinitesimal_and_negative_heights(domain):
     
    10641029        #self.xmom    = domain.quantities['xmomentum'].edge_values
    10651030        #self.ymom    = domain.quantities['ymomentum'].edge_values
    1066         self.normals = domain.normals
    1067         self.stage   = domain.quantities['stage'].vertex_values
    1068         self.xmom    = domain.quantities['xmomentum'].vertex_values
     1031        self.normals  = domain.normals
     1032        self.stage    = domain.quantities['stage'].vertex_values
     1033        self.xmom     = domain.quantities['xmomentum'].vertex_values
     1034        self.bed      = domain.quantities['elevation'].vertex_values
     1035        self.height   = domain.quantities['height'].vertex_values
     1036        self.velocity = domain.quantities['velocity'].vertex_values
    10691037
    10701038        from Numeric import zeros, Float
    10711039        #self.conserved_quantities = zeros(3, Float)
    1072         self.conserved_quantities = zeros(2, Float)
     1040        #self.conserved_quantities = zeros(2, Float)
     1041        self.evolved_quantities = zeros(5, Float)
    10731042
    10741043    def __repr__(self):
    10751044        return 'Reflective_boundary'
    10761045
    1077 
     1046    """
    10781047    def evaluate(self, vol_id, edge_id):
    1079         """Reflective boundaries reverses the outward momentum
    1080         of the volume they serve.
    1081         """
    1082 
    1083         q = self.conserved_quantities
     1048        #Reflective boundaries reverses the outward momentum
     1049        #of the volume they serve.
     1050       
     1051        #q = self.conserved_quantities
     1052        q = self.evolved_quantities
    10841053        q[0] = self.stage[vol_id, edge_id]
    10851054        q[1] = self.xmom[vol_id, edge_id]
     
    11011070
    11021071        return q
     1072    """
     1073    def evaluate(self, vol_id, edge_id):
     1074        """Reflective boundaries reverses the outward momentum
     1075        of the volume they serve.
     1076        """
     1077        q = self.evolved_quantities
     1078        q[0] = self.stage[vol_id, edge_id]
     1079        q[1] = -self.xmom[vol_id, edge_id]
     1080        q[2] = self.bed[vol_id, edge_id]
     1081        q[3] = self.height[vol_id, edge_id]
     1082        q[4] = -self.velocity[vol_id, edge_id]
     1083        return q
     1084
     1085
     1086   
    11031087
    11041088class Dirichlet_boundary(Boundary):
Note: See TracChangeset for help on using the changeset viewer.