Changeset 8446


Ignore:
Timestamp:
Jun 14, 2012, 12:31:32 AM (12 years ago)
Author:
davies
Message:

bal_dev: Updates

Location:
trunk/anuga_work/development/gareth
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_work/development/gareth/experimental/balanced_dev/__init__.py

    r8396 r8446  
    99from balanced_dev.swb2_domain import Domain as Domain
    1010from balanced_dev.swb2_domain import Domain
    11 from balanced_dev.swb2_boundary_conditions import zero_mass_flux_transmissive_momentum_flux_boundary as zero_mass_flux_transmissive_momentum_flux_boundary
     11#from balanced_dev.swb2_boundary_conditions import zero_mass_flux_transmissive_momentum_flux_boundary as zero_mass_flux_transmissive_momentum_flux_boundary
    1212#, \
    1313#     Transmissive_boundary,\
  • trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_boundary_conditions.py

    r8396 r8446  
    66#####
    77
    8 class  zero_mass_flux_transmissive_momentum_flux_boundary(Boundary):
    9     """ Boundary which operates directly on the fluxes
    10         CAN BE UNSTABLE -- CAREFUL
    11         Sets boundary values of h, uh, vh as in transmissive boundary, but
    12         also ensures:
    13         flux[0] = 0.0,
     8#class  zero_mass_flux_transmissive_momentum_flux_boundary(Boundary):
     9#    """ Boundary which operates directly on the fluxes
     10#        CAN BE UNSTABLE -- CAREFUL
     11#        Sets boundary values of h, uh, vh as in transmissive boundary, but
     12#        also ensures:
     13#        flux[0] = 0.0,
     14#    """
     15#    def __init__(self, domain = None):
     16#        Boundary.__init__(self)
     17#
     18#        if domain is None:
     19#            msg = 'Domain must be specified for zero_mass_flux_transmissive_momentum_flux boundary'
     20#            raise Exception, msg
     21#
     22#        self.domain = domain
     23#
     24#    def __repr__(self):
     25#        return 'zero_mass_flux_transmissive_momentum_flux_boundary(%s)' %self.domain
     26#
     27#    def evaluate(self, vol_id, edge_id):
     28#        """Transmissive boundaries return the edge values
     29#        of the volume they serve.
     30#        """
     31#
     32#        if self.domain.get_centroid_transmissive_bc() :
     33#            q = self.domain.get_evolved_quantities(vol_id)
     34#        else:
     35#            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
     36#        return q
     37#
     38######
     39#
     40#class  zero_mass_flux_zero_momentum_boundary(Boundary):
     41#    """ Boundary which operates directly on the fluxes
     42#        Sets boundary values of h=neighbour_h, uh=0, vh=0.
     43#        and ensure that:
     44#        flux[0] = 0.0
     45#    """
     46#    def __init__(self, domain = None):
     47#        Boundary.__init__(self)
     48#
     49#        if domain is None:
     50#            msg = 'Domain must be specified for zero_mass_flux_zero_momentum_boundary'
     51#            raise Exception, msg
     52#
     53#        self.domain = domain
     54#
     55#    def __repr__(self):
     56#        return 'zero_mass_flux_zero_momentum_boundary(%s)' %self.domain
     57#
     58#    def evaluate(self, vol_id, edge_id):
     59#        """ Apply chosen values to q[1], q[2]
     60#        """
     61#
     62#        if self.domain.get_centroid_transmissive_bc() :
     63#            q = self.domain.get_evolved_quantities(vol_id)
     64#        else:
     65#            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
     66#
     67#        normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
     68#
     69#        r = rotate(q, normal, direction = 1)
     70#        #r[1] = -r[1]
     71#        tdamp=0.0
     72#        ndamp=0.0
     73#        r[1]=r[1]*ndamp
     74#        r[2]=r[2]*tdamp
     75#        q = rotate(r, normal, direction = -1)
     76#
     77#        # q[1]=0.
     78#        # q[2]=0.
     79#
     80#        return q
     81#
     82#
     83#class  zero_mass_flux_zero_n_transmissive_t_momentum_boundary(Boundary):
     84#    """ Boundary which operates directly on the fluxes
     85#        Sets boundary values of h=neighbour_h, uh=0, vh=vh
     86#        and ensure that:
     87#        flux[0] = 0.0
     88#    """
     89#    def __init__(self, domain = None):
     90#        Boundary.__init__(self)
     91#
     92#        if domain is None:
     93#            msg = 'Domain must be specified for zero_mass_flux_zero_n_transmissive_t_momentum_boundary'
     94#            raise Exception, msg
     95#
     96#        self.domain = domain
     97#
     98#    def __repr__(self):
     99#        return 'zero_mass_flux_zero_n_transmissive_t_momentum_boundary(%s)' %self.domain
     100#
     101#    def evaluate(self, vol_id, edge_id):
     102#        """ Apply chosen values to q[1], q[2]
     103#        """
     104#
     105#        if self.domain.get_centroid_transmissive_bc() :
     106#            q = self.domain.get_evolved_quantities(vol_id)
     107#        else:
     108#            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
     109#
     110#        normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
     111#
     112#        r = rotate(q, normal, direction = 1)
     113#        #r[1] = -r[1]
     114#        tdamp=1.0
     115#        ndamp=0.0
     116#        r[1]=r[1]*ndamp
     117#        r[2]=r[2]*tdamp
     118#        q = rotate(r, normal, direction = -1)
     119#
     120#        # q[1]=0.
     121#        # q[2]=0.
     122#
     123#        return q
     124#
     125#
     126#class  zero_mass_flux_zero_t_transmissive_n_momentum_boundary(Boundary):
     127#    """ Boundary which operates directly on the fluxes
     128#        Sets boundary values of h=neighbour_h, uh=neighbour_uh, vh=0.
     129#        and ensure that:
     130#        flux[0] = 0.0
     131#    """
     132#    def __init__(self, domain = None):
     133#        Boundary.__init__(self)
     134#
     135#        if domain is None:
     136#            msg = 'Domain must be specified for zero_mass_flux_zero_t_transmissive_n_momentum_boundary'
     137#            raise Exception, msg
     138#
     139#        self.domain = domain
     140#
     141#    def __repr__(self):
     142#        return 'zero_mass_flux_zero_t_transmissive_n_momentum_boundary(%s)' %self.domain
     143#
     144#    def evaluate(self, vol_id, edge_id):
     145#        """ Apply chosen values to q[1], q[2]
     146#        """
     147#
     148#        if self.domain.get_centroid_transmissive_bc() :
     149#            q = self.domain.get_evolved_quantities(vol_id)
     150#        else:
     151#            q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
     152#
     153#        normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
     154#
     155#        r = rotate(q, normal, direction = 1)
     156#        #r[1] = -r[1]
     157#        tdamp=0.0
     158#        ndamp=1.0
     159#        r[1]=r[1]*ndamp
     160#        r[2]=r[2]*tdamp
     161#        q = rotate(r, normal, direction = -1)
     162#
     163#        # q[1]=0.
     164#        # q[2]=0.
     165#
     166#        return q
     167#
     168#
     169class Characteristic_boundary(Boundary):
    14170    """
    15     def __init__(self, domain = None):
     171
     172    Example:
     173
     174    def waveform(t):
     175        return sea_level + normalized_amplitude/cosh(t-25)**2
     176
     177    Bf = Nudge_boundary(domain, waveform)
     178
     179    Underlying domain must be specified when boundary is instantiated
     180    """
     181
     182    def __init__(self, domain=None, function=None):
     183        """ Instantiate a
     184            Nudge_boundary.
     185            domain is the domain containing the boundary
     186            function is the function to apply
     187        """
     188
    16189        Boundary.__init__(self)
    17190
    18191        if domain is None:
    19             msg = 'Domain must be specified for zero_mass_flux_transmissive_momentum_flux boundary'
     192            msg = 'Domain must be specified for this type boundary'
    20193            raise Exception, msg
    21194
     195        if function is None:
     196            msg = 'Function must be specified for this type boundary'
     197            raise Exception, msg
     198
    22199        self.domain = domain
     200        self.function = function
     201
    23202
    24203    def __repr__(self):
    25         return 'zero_mass_flux_transmissive_momentum_flux_boundary(%s)' %self.domain
     204        """ Return a representation of this instance. """
     205        msg = 'Nudge_boundary'
     206        msg += '(%s)' % self.domain
     207        return msg
     208
    26209
    27210    def evaluate(self, vol_id, edge_id):
    28         """Transmissive boundaries return the edge values
    29         of the volume they serve.
    30211        """
    31 
    32         if self.domain.get_centroid_transmissive_bc() :
    33             q = self.domain.get_evolved_quantities(vol_id)
     212        """
     213
     214        q = self.domain.get_conserved_quantities(vol_id, edge = edge_id)
     215        bed = self.domain.quantities['elevation'].centroid_values[vol_id]
     216        depth=max(q[0]-bed,0.0)
     217        dt=self.domain.timestep
     218
     219        normal = self.domain.get_normal(vol_id, edge_id)
     220
     221
     222        t = self.domain.get_time()
     223
     224        if hasattr(self.function, 'time'):
     225            # Roll boundary over if time exceeds
     226            while t > self.function.time[-1]:
     227                msg = 'WARNING: domain time %.2f has exceeded' % t
     228                msg += 'time provided in '
     229                msg += 'Flather_boundary object.\n'
     230                msg += 'I will continue, reusing the object from t==0'
     231                log.critical(msg)
     232                t -= self.function.time[-1]
     233
     234        value = self.function(t)
     235        try:
     236            x = float(value)
     237        except:
     238            x = float(value[0])
     239
     240        ##
     241        ndotq = (normal[0]*q[1] + normal[1]*q[2]) # momentum perpendicular to the boundary
     242
     243        if(depth==0.):
     244            q[0] = x
     245            q[1] = 0.
     246            q[2] = 0.
     247
    34248        else:
    35             q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
     249
     250            # Asssume sub-critical flow. Set the values of the characteristics as
     251            # appropriate, depending on whether we have inflow or outflow
     252            tmp = (9.8/depth)**0.5
     253            if(ndotq>0.):
     254                # Only need to impose one characteristic
     255                # Assume zero external velocity (dangerous?)
     256                # FIXME: Allow external velocities to be specified in the function
     257                # This would be needed to use the boundary for real inflows (e.g. run_wave.py)
     258                # At the moment, it seems to work okay as a 'passive' outflow boundary (e.g. run_wave.py)
     259                # Here it attenuates the wave at the outflow a bit, but doesn't do too much damage compared
     260                # with the other boundary conditions
     261                w1 = 0. - tmp*x
     262                w2 = (+normal[1]*q[1] -normal[0]*q[2])/depth
     263                w3 = ndotq/depth + tmp*q[0]
     264               
     265            else:
     266                # Need to set 2 characteristics
     267                # Assume zero external velocity (dangerous?)
     268                # FIXME: Allow this to be specified with a function
     269                w1 = 0. - tmp*x
     270                w2 = 0.
     271                #w2 = (+normal[1]*q[1] -normal[0]*q[2])/depth
     272                w3 = ndotq/depth + tmp*q[0]
     273
     274
     275            q[0] = (w3-w1)/(2*tmp)
     276            qperp= (w3+w1)/2. *depth
     277            qpar=  w2*depth
     278
     279            # So q[1], q[2] = qperp*(normal[0], normal[1]) + qpar*(-normal[1], normal[0])
     280
     281            q[1] = qperp*normal[0] + qpar*normal[1]
     282            q[2] = qperp*normal[1] -qpar*normal[0]
     283
    36284        return q
    37 
    38 #####
    39 
    40 class  zero_mass_flux_zero_momentum_boundary(Boundary):
    41     """ Boundary which operates directly on the fluxes
    42         Sets boundary values of h=neighbour_h, uh=0, vh=0.
    43         and ensure that:
    44         flux[0] = 0.0
    45     """
    46     def __init__(self, domain = None):
    47         Boundary.__init__(self)
    48 
    49         if domain is None:
    50             msg = 'Domain must be specified for zero_mass_flux_zero_momentum_boundary'
    51             raise Exception, msg
    52 
    53         self.domain = domain
    54 
    55     def __repr__(self):
    56         return 'zero_mass_flux_zero_momentum_boundary(%s)' %self.domain
    57 
    58     def evaluate(self, vol_id, edge_id):
    59         """ Apply chosen values to q[1], q[2]
    60         """
    61 
    62         if self.domain.get_centroid_transmissive_bc() :
    63             q = self.domain.get_evolved_quantities(vol_id)
    64         else:
    65             q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
    66 
    67         normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
    68 
    69         r = rotate(q, normal, direction = 1)
    70         #r[1] = -r[1]
    71         tdamp=0.0
    72         ndamp=0.0
    73         r[1]=r[1]*ndamp
    74         r[2]=r[2]*tdamp
    75         q = rotate(r, normal, direction = -1)
    76 
    77         # q[1]=0.
    78         # q[2]=0.
    79 
    80         return q
    81 
    82 
    83 class  zero_mass_flux_zero_n_transmissive_t_momentum_boundary(Boundary):
    84     """ Boundary which operates directly on the fluxes
    85         Sets boundary values of h=neighbour_h, uh=0, vh=vh
    86         and ensure that:
    87         flux[0] = 0.0
    88     """
    89     def __init__(self, domain = None):
    90         Boundary.__init__(self)
    91 
    92         if domain is None:
    93             msg = 'Domain must be specified for zero_mass_flux_zero_n_transmissive_t_momentum_boundary'
    94             raise Exception, msg
    95 
    96         self.domain = domain
    97 
    98     def __repr__(self):
    99         return 'zero_mass_flux_zero_n_transmissive_t_momentum_boundary(%s)' %self.domain
    100 
    101     def evaluate(self, vol_id, edge_id):
    102         """ Apply chosen values to q[1], q[2]
    103         """
    104 
    105         if self.domain.get_centroid_transmissive_bc() :
    106             q = self.domain.get_evolved_quantities(vol_id)
    107         else:
    108             q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
    109 
    110         normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
    111 
    112         r = rotate(q, normal, direction = 1)
    113         #r[1] = -r[1]
    114         tdamp=1.0
    115         ndamp=0.0
    116         r[1]=r[1]*ndamp
    117         r[2]=r[2]*tdamp
    118         q = rotate(r, normal, direction = -1)
    119 
    120         # q[1]=0.
    121         # q[2]=0.
    122 
    123         return q
    124 
    125 
    126 class  zero_mass_flux_zero_t_transmissive_n_momentum_boundary(Boundary):
    127     """ Boundary which operates directly on the fluxes
    128         Sets boundary values of h=neighbour_h, uh=neighbour_uh, vh=0.
    129         and ensure that:
    130         flux[0] = 0.0
    131     """
    132     def __init__(self, domain = None):
    133         Boundary.__init__(self)
    134 
    135         if domain is None:
    136             msg = 'Domain must be specified for zero_mass_flux_zero_t_transmissive_n_momentum_boundary'
    137             raise Exception, msg
    138 
    139         self.domain = domain
    140 
    141     def __repr__(self):
    142         return 'zero_mass_flux_zero_t_transmissive_n_momentum_boundary(%s)' %self.domain
    143 
    144     def evaluate(self, vol_id, edge_id):
    145         """ Apply chosen values to q[1], q[2]
    146         """
    147 
    148         if self.domain.get_centroid_transmissive_bc() :
    149             q = self.domain.get_evolved_quantities(vol_id)
    150         else:
    151             q = self.domain.get_evolved_quantities(vol_id, edge = edge_id)
    152 
    153         normal = self.domain.normals[vol_id, 2*edge_id:2*edge_id+2]
    154 
    155         r = rotate(q, normal, direction = 1)
    156         #r[1] = -r[1]
    157         tdamp=0.0
    158         ndamp=1.0
    159         r[1]=r[1]*ndamp
    160         r[2]=r[2]*tdamp
    161         q = rotate(r, normal, direction = -1)
    162 
    163         # q[1]=0.
    164         # q[2]=0.
    165 
    166         return q
    167 
    168 
    169285###
    170286def rotate(q,normal,direction=1):
  • trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain.py

    r8397 r8446  
    7575        # Note that the extrapolation method used in quantity_ext.c (e.g.
    7676        # extrapolate_second_order_and_limit_by_edge) uses a constant value for
    77         # all the betas == beta_rk2 in our case
    78         self.beta_w=2.0
    79         self.beta_w_dry=2.0
    80         self.beta_uh=2.0
    81         self.beta_uh_dry=2.0
    82         self.beta_vh=2.0
    83         self.beta_vh_dry=2.0
     77        # all the betas
     78        self.beta_w=1.0
     79        self.beta_w_dry=1.0
     80        self.beta_uh=1.0
     81        self.beta_uh_dry=1.0
     82        self.beta_vh=1.0
     83        self.beta_vh_dry=1.0
    8484
    8585        #self.optimise_dry_cells=True
  • trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c

    r8397 r8446  
    227227      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
    228228      //if(i<2) {
    229         edgeflux[i] += s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3);
     229      edgeflux[i] += s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3);
    230230      //}else{
    231231      //  edgeflux[i] += 0.5*s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3);
     
    307307    int useint;
    308308    double stage_edge_lim, scale_factor_shallow, bedtop, bedbot, pressure_flux, hc, hc_n;
    309     double max_bed_edgevalue[number_of_elements];
    310     //double depthtemp, max_bed_edge_value;
     309    //double min_bed_edgevalue[number_of_elements]; // FIXME: Set memory explicitly
    311310    static long call = 1; // Static local variable flagging already computed flux
    312311
     312    double *min_bed_edgevalue;
     313
     314    min_bed_edgevalue = malloc(number_of_elements*sizeof(double));
    313315    // Start computation
    314316    call++; // Flag 'id' of flux calculation for this timestep
     
    321323
    322324
    323     // Compute maximum bed edge value on each triangle
    324     //for (k = 0; k < number_of_elements; k++){
    325     //    max_bed_edgevalue[k] = max(bed_edge_values[3*k],
    326     //                               max(bed_edge_values[3*k+1], bed_edge_values[3*k+2]));
    327     //
    328     //}
     325    // Compute minimum bed edge value on each triangle
     326    for (k = 0; k < number_of_elements; k++){
     327        min_bed_edgevalue[k] = min(bed_edge_values[3*k],
     328                                   min(bed_edge_values[3*k+1], bed_edge_values[3*k+2]));
     329        //min_bed_edgevalue[k] = bed_centroid_values[k];
     330   
     331    }
    329332
    330333
     
    402405            // NOTE: This also means no bedslope term -- which is arguably
    403406            // appropriate, since the depth is zero at the cell centroid
    404             if(( hc == 0.0)&((hc_n == 0.0)|((n<0)&(qr[0]<zl)))){
     407            if(( hc == 0.0)&(hc_n == 0.0)){
     408                already_computed_flux[ki] = call; // #k Done
     409                already_computed_flux[nm] = call; // #n Done
     410                max_speed = 0.0;
    405411                continue;
    406412            }
     
    410416            if(n>=0){
    411417                if(hc==0.0){
    412                     ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
    413                     //ql[0]=qr[0]; // max(min(qr[0],stage_centroid_values[k]),zl);
     418                    //ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
     419                    //ql[0]=max(min(qr[0],0.5*(stage_centroid_values[k]+stage_centroid_values[n])),zl);
     420                    ql[0]=qr[0]; // max(min(qr[0],stage_centroid_values[k]),zl);
    414421                }
    415422                if(hc_n==0.0){
    416                     qr[0]=max(min(ql[0],stage_centroid_values[n]),zr);
    417                     //qr[0]=ql[0];
     423                    //qr[0]=max(min(ql[0],stage_centroid_values[n]),zr);
     424                    //qr[0]=max(min(ql[0],0.5*(stage_centroid_values[n]+stage_centroid_values[k])),zr);
     425                    qr[0]=ql[0];
    418426                }
    419427            }else{
     
    421429                if((hc==0.0)){
    422430                    ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
     431                    //ql[0]=max(min(qr[0],ql[0]),zl);
    423432                }
    424433            }
     
    432441                    epsilon, h0, limiting_threshold, g,
    433442                    edgeflux, &max_speed, &pressure_flux);
    434             //_flux_function_fraccaro(ql, qr, zl, zr,
    435             //        normals[ki2], normals[ki2 + 1],
    436             //        epsilon, h0, limiting_threshold, g,
    437             //        edgeflux, &max_speed);
    438    
     443
     444            // Prevent outflow from 'seriously' dry cells
     445            if(stage_centroid_values[k]<min_bed_edgevalue[k]){
     446                if(edgeflux[0]>0.0){
     447                    edgeflux[0]=0.0;
     448                    edgeflux[1]=0.0;
     449                    edgeflux[2]=0.0;
     450                }
     451            }
     452            if(n>=0){
     453                if(stage_centroid_values[n]<min_bed_edgevalue[n]){
     454                    if(edgeflux[0]<0.0){
     455                        edgeflux[0]=0.0;
     456                        edgeflux[1]=0.0;
     457                        edgeflux[2]=0.0;
     458                    }
     459                }
     460            }
     461
    439462            // Multiply edgeflux by edgelength
    440463            length = edgelengths[ki];
     
    443466            edgeflux[2] *= length;
    444467
    445             if(n<0){
    446                 if(boundary_flux_type[m]>0){
    447                     // IMPOSE BOUNDARY CONDITIONS DIRECTLY ON THE FLUX
    448                     if(boundary_flux_type[m]==1){
    449                         // Zero_mass_flux_transmissive_momentum_flux boundary, or
    450                         // zero_mass_flux_zero_momentum_boundary
    451                         // Just need to set the mass flux to zero, as the
    452                         // transmissive momentum is ensured by the values of
    453                         // qr[0], qr[1], qr[2]
    454                         edgeflux[0] = 0.0;
    455                     }
    456                 }
    457             }
     468            // GET RID OF THIS: EXPERIMENTAL
     469            //if(n<0){
     470            //    if(boundary_flux_type[m]>0){
     471            //        // IMPOSE BOUNDARY CONDITIONS DIRECTLY ON THE FLUX
     472            //        if(boundary_flux_type[m]==1){
     473            //            // Zero_mass_flux_transmissive_momentum_flux boundary, or
     474            //            // zero_mass_flux_zero_momentum_boundary
     475            //            // Just need to set the mass flux to zero, as the
     476            //            // transmissive momentum is ensured by the values of
     477            //            // qr[0], qr[1], qr[2]
     478            //            printf("Warning: WEIRD EDGEFLUX THING IS ACTING");
     479            //            edgeflux[0] = 0.0;
     480            //        }
     481            //    }
     482            //}
    458483
    459484            // Update triangle k with flux from edge i
     
    462487            ymom_explicit_update[k] -= edgeflux[2];
    463488
    464             // Calculate the bed slope term, using the 'divergence form' from
    465             // Alessandro Valiani and Lorenzo Begnudelli, Divergence Form for Bed Slope Source Term in Shallow Water Equations
    466             // Journal of Hydraulic Engineering, 2006, 132:652-665
    467             if(hc>0.0e-03){
     489            // Compute bed slope term
     490            if(hc>-9999.0){
    468491                //Bedslope approx 1:
    469492                //bedslope_work = g*hc*(ql[0])*length-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length;
     
    486509            }else{
    487510                // Treat nearly dry cells
    488                 //bedslope_work =-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length; //
     511                bedslope_work =-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length; //
    489512                //
    490                 bedslope_work = -pressure_flux*length;
     513                //bedslope_work = -pressure_flux*length;
    491514                xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
    492515                ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
     
    506529                ymom_explicit_update[n] += edgeflux[2];
    507530                //Add bed slope term here
    508                 if(hc_n>0.0e-03){
     531                if(hc_n>-9999.0){
    509532                //if(stage_centroid_values[n] > max_bed_edgevalue[n]){
    510533                    // Bedslope approx 1:
     
    528551                }else{
    529552                    // Treat nearly dry cells
    530                     //bedslope_work = -0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length; //-pressure_flux*length; //-0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length;
    531                     bedslope_work = -pressure_flux*length;
     553                    bedslope_work = -0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length; //-pressure_flux*length; //-0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length;
     554                    //bedslope_work = -pressure_flux*length;
    532555                    xmom_explicit_update[n] += normals[ki2]*bedslope_work;
    533556                    ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
     
    578601    } // End triangle k
    579602
     603    free(min_bed_edgevalue);
    580604
    581605    return timestep;
     
    599623
    600624  // This acts like minimum_allowed height, but scales with the vertical
    601   // distance between the bed_centroid_value and the max bed_edge_value of every triangle.
    602   double minimum_relative_height=0.05;
     625  // distance between the bed_centroid_value and the max bed_edge_value of
     626  // every triangle.
     627  double minimum_relative_height=0.1;
    603628
    604629
     
    619644        if (hc <= 0.0){
    620645             // Definine the minimum bed edge value on triangle k.
    621              bmin = 0.5*min((zv[3*k]+zv[3*k+1]),min((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
     646             // Setting = minimum edge value can lead to mass conservation problems
     647             //bmin = 0.5*min((zv[3*k]+zv[3*k+1]),min((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
     648             //bmin =0.5*bmin + 0.5*min(zv[3*k],min(zv[3*k+1],zv[3*k+2]));
     649             // Setting = minimum vertex value seems better, but might tend to be less smooth
     650             bmin =min(zv[3*k],min(zv[3*k+1],zv[3*k+2]));
    622651             // Minimum allowed stage = bmin
     652             // WARNING: ADDING MASS if wc[k]<bmin
     653             if(wc[k]<bmin){
     654                 printf("Adding mass to dry cell\n");
     655             }
    623656             wc[k] = max(wc[k], bmin) ;
     657
     658             // Set vertex values as well. Seems that this shouldn't be
     659             // needed. However from memory this is important at the first
     660             // time step, for 'dry' areas where the designated stage is
     661             // less than the bed centroid value
    624662             wv[3*k] = max(wv[3*k], bmin);
    625663             wv[3*k+1] = max(wv[3*k+1], bmin);
     
    683721  // dq2=q(vertex2)-q(vertex0)
    684722
    685   // This is a simple implementation -- the one below is the
    686   // original version -- I wonder if it is any faster?
     723  // This is a simple implementation
    687724  *qmax = max(max(dq0, max(dq0+dq1, dq0+dq2)), 0.0) ;
    688725  *qmin = min(min(dq0, min(dq0+dq1, dq0+dq2)), 0.0) ;
    689726 
    690   //if (dq0>=0.0){
    691   //  if (dq1>=dq2){
    692   //    if (dq1>=0.0)
    693   //  *qmax=dq0+dq1;
    694   //    else
    695   //  *qmax=dq0;
    696   //   
    697   //    *qmin=dq0+dq2;
    698   //    if (*qmin>=0.0) *qmin = 0.0;
    699   //  }
    700   //  else{// dq1<dq2
    701   //    if (dq2>0)
    702   //  *qmax=dq0+dq2;
    703   //    else
    704   //  *qmax=dq0;
    705   //   
    706   //    *qmin=dq0+dq1;   
    707   //    if (*qmin>=0.0) *qmin=0.0;
    708   //  }
    709   //}
    710   //else{//dq0<0
    711   //  if (dq1<=dq2){
    712   //    if (dq1<0.0)
    713   //  *qmin=dq0+dq1;
    714   //    else
    715   //  *qmin=dq0;
    716   //   
    717   //    *qmax=dq0+dq2;   
    718   //    if (*qmax<=0.0) *qmax=0.0;
    719   //  }
    720   //  else{// dq1>dq2
    721   //    if (dq2<0.0)
    722   //  *qmin=dq0+dq2;
    723   //    else
    724   //  *qmin=dq0;
    725   //   
    726   //    *qmax=dq0+dq1;
    727   //    if (*qmax<=0.0) *qmax=0.0;
    728   //  }
    729   //}
    730727  return 0;
    731728}
     
    759756  phi=min(r*beta_w,1.0);
    760757  //phi=1.;
    761   //for (i=0;i<3;i++)
    762758  dqv[0]=dqv[0]*phi;
    763759  dqv[1]=dqv[1]*phi;
     
    766762  return 0;
    767763}
     764
     765//int limit_gradient2(double *dqv, double qmin, double qmax, double beta_w, double r0scale){
     766//  // EXPERIMENTAL LIMITER, DOESN'T WORK
     767//  // Given provisional jumps dqv from the FV triangle centroid to its
     768//  // vertices and jumps qmin (qmax) between the centroid of the FV
     769//  // triangle and the minimum (maximum) of the values at the centroid of
     770//  // the FV triangle and the auxiliary triangle vertices,
     771//  // calculate a multiplicative factor phi by which the provisional
     772//  // vertex jumps are to be limited
     773// 
     774//  int i;
     775//  double r=1000.0, r0=1.0, phi=1.0;
     776//  static double TINY = 1.0e-100; // to avoid machine accuracy problems.
     777//  // FIXME: Perhaps use the epsilon used elsewhere.
     778// 
     779//  for (i=0;i<3;i++){
     780//    if (dqv[i]<-TINY)
     781//      r0=qmin/dqv[i]*r0scale*2.0;
     782//     
     783//    if (dqv[i]>TINY)
     784//      r0=qmax/dqv[i]*r0scale*2.0;
     785//     
     786//    r=min(r0,r);
     787//  }
     788// 
     789//  phi=min(r*beta_w,1.0);
     790//  //phi=1.;
     791//  //for (i=0;i<3;i++)
     792//  dqv[0]=dqv[0]*phi;
     793//  dqv[1]=dqv[1]*phi;
     794//  dqv[2]=dqv[2]*phi;
     795//   
     796//  return 0;
     797//}
    768798
    769799// Computational routine
     
    789819                                 double* ymom_edge_values,
    790820                                 double* elevation_edge_values,
     821                                 double* stage_vertex_values,
     822                                 double* xmom_vertex_values,
     823                                 double* ymom_vertex_values,
     824                                 double* elevation_vertex_values,
    791825                                 int optimise_dry_cells,
    792826                                 int extrapolate_velocity_second_order) {
    793827                 
    794   // Note that although this routine refers to EDGE values, it can equally well
    795   // be applied to vertex values -- it is a modified version of such a routine.
    796   // Also note that momentum can be extrapolated using velocity, via
    797   // modifications at the start and end of the routine.                 
    798 
    799828  // Local variables
    800829  double a, b; // Gradient vector used to calculate edge values from centroids
     
    804833  double dqv[3], qmin, qmax, hmin, hmax, bedmax, stagemin;
    805834  double hc, h0, h1, h2, beta_tmp, hfactor;
    806   double dk, dv0, dv1, dv2, de[3];
    807   //double xmom_centroid_store[number_of_elements], ymom_centroid_store[number_of_elements]
    808   //double stage_centroid_store[number_of_elements];
     835  double dk, dv0, dv1, dv2, de[3], demin, dcmax, r0scale;
    809836 
    810837  double *xmom_centroid_store, *ymom_centroid_store, *stage_centroid_store;
     
    828855          ymom_centroid_values[k] = ymom_centroid_values[k]/dk;
    829856
    830           // Experimental -- replace stage with a weighting of 'stage' and
    831           // 'depth', where the weights depend on the froude number **2
    832           // stage_centroid_store[k]=stage_centroid_values[k];
    833           // hfactor =(xmom_centroid_values[k]*xmom_centroid_values[k] +
    834           //          ymom_centroid_values[k]*ymom_centroid_values[k])/20.0;
    835           //          //(9.8*max(stage_centroid_values[k] - elevation_centroid_values[k],minimum_allowed_height));
    836           //stage_centroid_values[k] -= min(1.0, hfactor)*elevation_centroid_values[k];
    837857      }
    838858  }
     
    887907      dyv1 = yv1 - y;
    888908      dyv2 = yv2 - y;
     909      // Compute the minimum distance from the centroid to an edge
     910      //demin=min(dxv0*dxv0 +dyv0*dyv0, min(dxv1*dxv1+dyv1*dyv1, dxv2*dxv2+dyv2*dyv2));
     911      //demin=sqrt(demin);
    889912    }
    890913
     
    911934      // triangle area will be zero, and a first order extrapolation will be
    912935      // used
    913       //if(stage_centroid_values[k2]<=elevation_centroid_values[k2]){
    914       //    k2 = k ;
    915       //}
    916       //if(stage_centroid_values[k0]<=elevation_centroid_values[k0]){
    917       //    k0 = k ;
    918       //}
    919       //if(stage_centroid_values[k1]<=elevation_centroid_values[k1]){
    920       //    k1 = k ;
    921       //}
     936      if(stage_centroid_values[k2]<=elevation_centroid_values[k2]){
     937          k2 = k ;
     938      }
     939      if(stage_centroid_values[k0]<=elevation_centroid_values[k0]){
     940          k0 = k ;
     941      }
     942      if(stage_centroid_values[k1]<=elevation_centroid_values[k1]){
     943          k1 = k ;
     944      }
    922945      // Alternative approach (BRUTAL) -- if the max neighbour bed elevation is greater
    923946      // than the min neighbour stage, then we use first order extrapolation
    924       bedmax = max(elevation_centroid_values[k],
    925                    max(elevation_centroid_values[k0],
    926                        max(elevation_centroid_values[k1], elevation_centroid_values[k2])));
    927       stagemin =min(stage_centroid_values[k],
    928                    min(stage_centroid_values[k0],
    929                        min(stage_centroid_values[k1], stage_centroid_values[k2])));
    930  
    931       if(stagemin < bedmax){
    932          // This will cause first order extrapolation
    933          k2 = k;
    934          k0 = k;
    935          k1 = k;
    936       }
     947      //bedmax = max(elevation_centroid_values[k],
     948      //             max(elevation_centroid_values[k0],
     949      //                 max(elevation_centroid_values[k1], elevation_centroid_values[k2])));
     950      //stagemin = min(stage_centroid_values[k],
     951      //             min(stage_centroid_values[k0],
     952      //                 min(stage_centroid_values[k1], stage_centroid_values[k2])));
     953      //
     954      //if(stagemin < bedmax){
     955      //   // This will cause first order extrapolation
     956      //   k2 = k;
     957      //   k0 = k;
     958      //   k1 = k;
     959      //}
    937960     
    938961      // Get the auxiliary triangle's vertex coordinates
     
    949972      x2 = centroid_coordinates[coord_index];
    950973      y2 = centroid_coordinates[coord_index+1];
    951      
     974
     975      // compute the maximum distance from the centroid to a neighbouring
     976      // centroid
     977      //dcmax=max( (x0-x)*(x0-x) + (y0-y)*(y0-y),     
     978      //           max((x1-x)*(x1-x) + (y1-y)*(y1-y),
     979      //               (x2-x)*(x2-x) + (y2-y)*(y2-y)));
     980      //dcmax=sqrt(dcmax);
     981      //// Ratio of centroid to edge distance -- useful in attempting to adapt limiter
     982      //if(dcmax>0.){
     983      //    r0scale=demin/dcmax;
     984      //    //printf("%f \n", r0scale);
     985      //}else{
     986      //    r0scale=0.5;
     987      //}
     988
    952989      // Store x- and y- differentials for the vertices
    953990      // of the auxiliary triangle
     
    9701007          //return -1;
    9711008
    972           //stage_edge_values[k3]   = 0.5*(stage_centroid_values[k]+stage_centroid_values[k0]);
    973           //stage_edge_values[k3+1] = 0.5*(stage_centroid_values[k]+stage_centroid_values[k1]);
    974           //stage_edge_values[k3+2] = 0.5*(stage_centroid_values[k]+stage_centroid_values[k2]);
    9751009          stage_edge_values[k3]   = stage_centroid_values[k];
    9761010          stage_edge_values[k3+1] = stage_centroid_values[k];
     
    10021036      {
    10031037        hfactor = 1.0 ;//hmin/(hmin + 0.004);
     1038        //hfactor=hmin/(hmin + 0.004);
    10041039      }
    10051040     
     
    10481083      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
    10491084     
    1050       // Playing with dry wet interface
    1051       //hmin = qmin;
    1052       //beta_tmp = beta_w_dry;
    1053       //if (hmin>minimum_allowed_height)
    10541085      beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
    10551086   
    1056       //printf("min_alled_height = %f\n",minimum_allowed_height);
    1057       //printf("hmin = %f\n",hmin);
    1058       //printf("beta_w = %f\n",beta_w);
    1059       //printf("beta_tmp = %f\n",beta_tmp);
    10601087      // Limit the gradient
    10611088      limit_gradient(dqv, qmin, qmax, beta_tmp);
     1089      //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale);
    10621090     
    10631091      //for (i=0;i<3;i++)
     
    10661094      stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2];
    10671095     
    1068       //if(k==300){
    1069       //  printf("Stage centroid values: %f, v1, %f, v2, %f, v3, %f \n", stage_centroid_values[k],
    1070       //          stage_edge_values[k3+0], stage_edge_values[k3+1], stage_edge_values[k3+2]);
    1071       //}     
    1072  
    10731096      //-----------------------------------
    10741097      // xmomentum
     
    11021125      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
    11031126
    1104       //beta_tmp = beta_uh;
    1105       //if (hmin<minimum_allowed_height)
    1106       //beta_tmp = beta_uh_dry;
    11071127      beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor;
    1108 
    1109       //if(k==300){
    1110       //  printf("dqv0: %f, dqv1, %f, dqv2, %f \n", dqv[0],dqv[1], dqv[2]) ;
    1111 
    1112       //}
    11131128
    11141129      // Limit the gradient
    11151130      limit_gradient(dqv, qmin, qmax, beta_tmp);
    1116      
    1117       //if(k==300){
    1118       //  printf("dqv0: %f, dqv1, %f, dqv2, %f \n", dqv[0],dqv[1], dqv[2]) ;
    1119 
    1120       //}
     1131      //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale);
     1132     
    11211133
    11221134      for (i=0; i < 3; i++)
     
    11241136        xmom_edge_values[k3+i] = xmom_centroid_values[k] + dqv[i];
    11251137      }
    1126      
    1127       //if(k==300){
    1128       //  printf("xmom centroid values: %f, v1, %f, v2, %f, v3, %f \n", xmom_centroid_values[k],
    1129       //          xmom_edge_values[k3+0], xmom_edge_values[k3+1], xmom_edge_values[k3+2]);
    1130       //}     
    11311138     
    11321139      //-----------------------------------
     
    11611168      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
    11621169     
    1163       //beta_tmp = beta_vh;
    1164       //
    1165       //if (hmin<minimum_allowed_height)
    1166       //beta_tmp = beta_vh_dry;
    11671170      beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;   
    11681171
    11691172      // Limit the gradient
    11701173      limit_gradient(dqv, qmin, qmax, beta_tmp);
     1174      //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale);
    11711175     
    11721176      for (i=0;i<3;i++)
     
    11751179      }
    11761180     
    1177       //if(k==300){
    1178       //  printf("ymom centroid values: %f, v1, %f, v2, %f, v3, %f \n", ymom_centroid_values[k],
    1179       //          ymom_edge_values[k3+0], ymom_edge_values[k3+1], ymom_edge_values[k3+2]);
    1180       //}     
    11811181    } // End number_of_boundaries <=1
    11821182    else
     
    13461346      limit_gradient(dqv, qmin, qmax, beta_w);
    13471347     
    1348       //for (i=0;i<3;i++)
    1349       //ymom_edge_values[k3] = ymom_centroid_values[k] + dqv[0];
    1350       //ymom_edge_values[k3 + 1] = ymom_centroid_values[k] + dqv[1];
    1351       //ymom_edge_values[k3 + 2] = ymom_centroid_values[k] + dqv[2];
    1352 
    1353 for (i=0;i<3;i++)
    1354         {
    1355         ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
    1356         }
    1357 //ymom_edge_values[k3] = ymom_centroid_values[k] + dqv[0];
    1358 //ymom_edge_values[k3 + 1] = ymom_centroid_values[k] + dqv[1];
    1359 //ymom_edge_values[k3 + 2] = ymom_centroid_values[k] + dqv[2];
     1348      for (i=0;i<3;i++)
     1349              {
     1350              ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
     1351              }
    13601352    } // else [number_of_boundaries==2]
    13611353  } // for k=0 to number_of_elements-1
    13621354 
    1363   if(extrapolate_velocity_second_order==1){
    1364       // Convert back from velocity to momentum
    1365      
    1366       for (k=0; k<number_of_elements; k++){
    1367           k3=3*k;
    1368 
    1369           // Convert energy back to stage
    1370           //stage_centroid_values[k] = stage_centroid_store[k];
    1371 
    1372           //Convert velocity back to momenta
     1355  // Compute vertex values of quantities
     1356  for (k=0; k<number_of_elements; k++){
     1357      k3=3*k;
     1358
     1359      // Compute stage vertex values
     1360      stage_vertex_values[k3] = stage_edge_values[k3+1] + stage_edge_values[k3+2] -stage_edge_values[k3] ;
     1361      stage_vertex_values[k3+1] =  stage_edge_values[k3] + stage_edge_values[k3+2]-stage_edge_values[k3+1];
     1362      stage_vertex_values[k3+2] =  stage_edge_values[k3] + stage_edge_values[k3+1]-stage_edge_values[k3+2];
     1363     
     1364     // Compute xmom vertex values
     1365      xmom_vertex_values[k3] = xmom_edge_values[k3+1] + xmom_edge_values[k3+2] -xmom_edge_values[k3] ;
     1366      xmom_vertex_values[k3+1] =  xmom_edge_values[k3] + xmom_edge_values[k3+2]-xmom_edge_values[k3+1];
     1367      xmom_vertex_values[k3+2] =  xmom_edge_values[k3] + xmom_edge_values[k3+1]-xmom_edge_values[k3+2];
     1368
     1369      // Compute ymom vertex values
     1370      ymom_vertex_values[k3] = ymom_edge_values[k3+1] + ymom_edge_values[k3+2] -ymom_edge_values[k3] ;
     1371      ymom_vertex_values[k3+1] =  ymom_edge_values[k3] + ymom_edge_values[k3+2]-ymom_edge_values[k3+1];
     1372      ymom_vertex_values[k3+2] =  ymom_edge_values[k3] + ymom_edge_values[k3+1]-ymom_edge_values[k3+2];
     1373
     1374      // If needed, convert from velocity to momenta
     1375      if(extrapolate_velocity_second_order==1){
     1376          //Convert velocity back to momenta at centroids
    13731377          xmom_centroid_values[k] = xmom_centroid_store[k];
    13741378          ymom_centroid_values[k] = ymom_centroid_store[k];
    13751379
     1380          // Re-compute momenta at edges
    13761381          for (i=0; i<3; i++){
    1377               //stage_edge_values[k3+i] -=(xmom_edge_values[k3+i]*xmom_edge_values[k3+i]
    1378               //                           + ymom_edge_values[k3+i]*ymom_edge_values[k3+i])*0.0;
    1379               //stage_edge_values[k3+i] += elevation_edge_values[k3+i];
    1380               //de[i] = max(stage_centroid_values[k]-elevation_centroid_values[k],minimum_allowed_height);
    1381               //for(ii=0; ii<1; ii++){
    1382               //hfactor =(xmom_edge_values[k3+i]*xmom_edge_values[k3+i] +
    1383               //          ymom_edge_values[k3+i]*ymom_edge_values[k3+i])/20.0; // 20 ~= 2*g
    1384                         //(9.8*de[i]);
    1385               //hfactor = 0.0; //max(hfactor/20.0 - 0.2, 0.0);
    1386               //stage_edge_values[k3+i] = stage_edge_values[k3+i];//+ min(hfactor, 1.0)*elevation_edge_values[k3+i];
    1387               //de[i] = max(stage_edge_values[k3+i]-elevation_edge_values[k3+i],1.0e-03);
    1388               //}
    1389 
    13901382              de[i] = max(stage_edge_values[k3+i]-elevation_edge_values[k3+i],0.0);
    13911383              xmom_edge_values[k3+i]=xmom_edge_values[k3+i]*de[i];
     
    13931385          }
    13941386
    1395          
    1396       }
    1397   }
     1387          // Re-compute momenta at vertices
     1388          for (i=0; i<3; i++){
     1389              de[i] = max(stage_vertex_values[k3+i]-elevation_vertex_values[k3+i],0.0);
     1390              xmom_vertex_values[k3+i]=xmom_vertex_values[k3+i]*de[i];
     1391              ymom_vertex_values[k3+i]=ymom_vertex_values[k3+i]*de[i];
     1392          }
     1393
     1394      }
     1395
     1396   
     1397  }
    13981398
    13991399  free(xmom_centroid_store);
     
    14041404}           
    14051405
    1406 int extrapolate_from_edges_to_vertices(int number_of_elements,
    1407                                         double* stage_edge_values,
    1408                                         double* xmom_edge_values,
    1409                                         double* ymom_edge_values,
    1410                                         double* elevation_edge_values,
    1411                                         double* stage_vertex_values,
    1412                                         double* xmom_vertex_values,
    1413                                         double* ymom_vertex_values,
    1414                                         double* elevation_vertex_values,
    1415                                         int extrapolate_velocity_second_order) {
    1416 
    1417   int i, i3;
    1418   double v0, v1, v2;
    1419 
    1420   for(i=0; i<number_of_elements; i++){
    1421         i3 = 3*i;
    1422         stage_vertex_values[i3] = stage_edge_values[i3+1] + stage_edge_values[i3+2] -stage_edge_values[i3] ;
    1423         stage_vertex_values[i3+1] =  stage_edge_values[i3] + stage_edge_values[i3+2]-stage_edge_values[i3+1];
    1424         stage_vertex_values[i3+2] =  stage_edge_values[i3] + stage_edge_values[i3+1]-stage_edge_values[i3+2];
    1425 
    1426         if(extrapolate_velocity_second_order==1){
    1427             // Compute x-velocity, carefully to avoid division by zero
    1428             //v0 = xmom_edge_values[i3+0]/(stage_edge_values[i3+0]-elevation_edge_values[i3+0]);
    1429             //v1 = xmom_edge_values[i3+1]/(stage_edge_values[i3+1]-elevation_edge_values[i3+1]);
    1430             //v2 = xmom_edge_values[i3+2]/(stage_edge_values[i3+2]-elevation_edge_values[i3+2]);
    1431             if(stage_edge_values[i3]>elevation_edge_values[i3]){
    1432                 v0 = xmom_edge_values[i3]/(stage_edge_values[i3]-elevation_edge_values[i3]);
    1433             }else{
    1434                 v0 = 0.0;
    1435             }
    1436             if(stage_edge_values[i3+1]>elevation_edge_values[i3+1]){
    1437                 v1 = xmom_edge_values[i3+1]/(stage_edge_values[i3+1]-elevation_edge_values[i3+1]);
    1438             }else{
    1439                 v1 = 0.0;
    1440             }
    1441 
    1442             if(stage_edge_values[i3+2]>elevation_edge_values[i3+2]){
    1443                 v2 = xmom_edge_values[i3+2]/(stage_edge_values[i3+2]-elevation_edge_values[i3+2]);
    1444             }else{
    1445                 v2 = 0.0;
    1446             }
    1447             //v1 = xmom_edge_values[i3+1]/(stage_edge_values[i3+1]-elevation_edge_values[i3+1]);
    1448             //v2 = xmom_edge_values[i3+2]/(stage_edge_values[i3+2]-elevation_edge_values[i3+2]);
    1449 
    1450             xmom_vertex_values[i3] = (-v0 + v1 +v2)*(stage_vertex_values[i3]-elevation_vertex_values[i3]);
    1451             xmom_vertex_values[i3+1] = (-v1 + v0 +v2)*(stage_vertex_values[i3+1]-elevation_vertex_values[i3+1]);
    1452             xmom_vertex_values[i3+2] = (-v2 + v0 +v1)*(stage_vertex_values[i3+2]-elevation_vertex_values[i3+2]);
    1453            
    1454             // Compute y-velocity, carefully to avoid division by zero
    1455             //v0 = ymom_edge_values[i3+0]/(stage_edge_values[i3+0]-elevation_edge_values[i3+0]);
    1456             //v1 = ymom_edge_values[i3+1]/(stage_edge_values[i3+1]-elevation_edge_values[i3+1]);
    1457             //v2 = ymom_edge_values[i3+2]/(stage_edge_values[i3+2]-elevation_edge_values[i3+2]);
    1458             if(stage_edge_values[i3]>elevation_edge_values[i3]){
    1459                 v0 = ymom_edge_values[i3]/(stage_edge_values[i3]-elevation_edge_values[i3]);
    1460             }else{
    1461                 v0 = 0.0;
    1462             }
    1463             if(stage_edge_values[i3+1]>elevation_edge_values[i3+1]){
    1464                 v1 = ymom_edge_values[i3+1]/(stage_edge_values[i3+1]-elevation_edge_values[i3+1]);
    1465             }else{
    1466                 v1 = 0.0;
    1467             }
    1468 
    1469             if(stage_edge_values[i3+2]>elevation_edge_values[i3+2]){
    1470                 v2 = ymom_edge_values[i3+2]/(stage_edge_values[i3+2]-elevation_edge_values[i3+2]);
    1471             }else{
    1472                 v2 = 0.0;
    1473             }
    1474             //v0 = ymom_edge_values[i3]/(stage_edge_values[i3]-elevation_edge_values[i3]);
    1475             //v1 = ymom_edge_values[i3+1]/(stage_edge_values[i3+1]-elevation_edge_values[i3+1]);
    1476             //v2 = ymom_edge_values[i3+2]/(stage_edge_values[i3+2]-elevation_edge_values[i3+2]);
    1477 
    1478             ymom_vertex_values[i3] = (-v0 + v1 +v2)*(stage_vertex_values[i3]-elevation_vertex_values[i3]);
    1479             ymom_vertex_values[i3+1] = (-v1 + v0 +v2)*(stage_vertex_values[i3+1]-elevation_vertex_values[i3+1]);
    1480             ymom_vertex_values[i3+2] = (-v2 + v0 +v1)*(stage_vertex_values[i3+2]-elevation_vertex_values[i3+2]);
    1481 
    1482         }else{
    1483            
    1484             xmom_vertex_values[i3] = xmom_edge_values[i3+1] + xmom_edge_values[i3+2]-xmom_edge_values[i3] ;
    1485             xmom_vertex_values[i3+1] =  xmom_edge_values[i3] + xmom_edge_values[i3+2]-xmom_edge_values[i3+1];
    1486             xmom_vertex_values[i3+2] =  xmom_edge_values[i3] + xmom_edge_values[i3+1]-xmom_edge_values[i3+2];
    1487 
    1488             ymom_vertex_values[i3] =  ymom_edge_values[i3+1] + ymom_edge_values[i3+2]-ymom_edge_values[i3];
    1489             ymom_vertex_values[i3+1] =  ymom_edge_values[i3] + ymom_edge_values[i3+2] -ymom_edge_values[i3+1];
    1490             ymom_vertex_values[i3+2] =  ymom_edge_values[i3] + ymom_edge_values[i3+1]-ymom_edge_values[i3+2] ;
    1491 
    1492         }
    1493 
    1494     }
    1495     return 0 ;
    1496 }
    14971406//=========================================================================
    14981407// Python Glue
    14991408//=========================================================================
    15001409
    1501 
    1502 //PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {
    1503 //  //
    1504 //  // r = rotate(q, normal, direction=1)
    1505 //  //
    1506 //  // Where q is assumed to be a Float numeric array of length 3 and
    1507 //  // normal a Float numeric array of length 2.
    1508 //
    1509 //  // FIXME(Ole): I don't think this is used anymore
    1510 //
    1511 //  PyObject *Q, *Normal;
    1512 //  PyArrayObject *q, *r, *normal;
    1513 //
    1514 //  static char *argnames[] = {"q", "normal", "direction", NULL};
    1515 //  int dimensions[1], i, direction=1;
    1516 //  double n1, n2;
    1517 //
    1518 //  // Convert Python arguments to C
    1519 //  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,
    1520 //                   &Q, &Normal, &direction)) {
    1521 //    report_python_error(AT, "could not parse input arguments");
    1522 //    return NULL;
    1523 //  } 
    1524 //
    1525 //  // Input checks (convert sequences into numeric arrays)
    1526 //  q = (PyArrayObject *)
    1527 //    PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);
    1528 //  normal = (PyArrayObject *)
    1529 //    PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);
    1530 //
    1531 //
    1532 //  if (normal -> dimensions[0] != 2) {
    1533 //    report_python_error(AT, "Normal vector must have 2 components");
    1534 //    return NULL;
    1535 //  }
    1536 //
    1537 //  // Allocate space for return vector r (don't DECREF)
    1538 //  dimensions[0] = 3;
    1539 //  r = (PyArrayObject *) anuga_FromDims(1, dimensions, PyArray_DOUBLE);
    1540 //
    1541 //  // Copy
    1542 //  for (i=0; i<3; i++) {
    1543 //    ((double *) (r -> data))[i] = ((double *) (q -> data))[i];
    1544 //  }
    1545 //
    1546 //  // Get normal and direction
    1547 //  n1 = ((double *) normal -> data)[0];
    1548 //  n2 = ((double *) normal -> data)[1];
    1549 //  if (direction == -1) n2 = -n2;
    1550 //
    1551 //  // Rotate
    1552 //  _rotate((double *) r -> data, n1, n2);
    1553 //
    1554 //  // Release numeric arrays
    1555 //  Py_DECREF(q);
    1556 //  Py_DECREF(normal);
    1557 //
    1558 //  // Return result using PyArray to avoid memory leak
    1559 //  return PyArray_Return(r);
    1560 //}
    15611410
    15621411//========================================================================
     
    19811830                   (double*) ymom_edge_values -> data,
    19821831                   (double*) elevation_edge_values -> data,
     1832                   (double*) stage_vertex_values -> data,
     1833                   (double*) xmom_vertex_values -> data,
     1834                   (double*) ymom_vertex_values -> data,
     1835                   (double*) elevation_vertex_values -> data,
    19831836                   optimise_dry_cells,
    19841837                   extrapolate_velocity_second_order);
    1985 
    1986   //printf("In C before edges-to-vertices");
    1987 
    1988   //Extrapolate from edges to vertices
    1989   e2 = extrapolate_from_edges_to_vertices(number_of_elements,
    1990                                       (double *) stage_edge_values -> data,
    1991                                       (double *) xmom_edge_values -> data,
    1992                                       (double *) ymom_edge_values -> data,
    1993                                       (double *) elevation_edge_values -> data,
    1994                                       (double *) stage_vertex_values -> data,
    1995                                       (double *) xmom_vertex_values -> data,
    1996                                       (double *) ymom_vertex_values -> data,
    1997                                       (double *) elevation_vertex_values -> data,
    1998                                       extrapolate_velocity_second_order);
    1999   //printf("In C after edges-to-vertices");
    20001838
    20011839  if (e == -1) {
  • trunk/anuga_work/development/gareth/tests/channel_floodplain/channel_floodplain1.py

    r8353 r8446  
    1010import numpy
    1111from anuga.structures.inlet_operator import Inlet_operator
    12 from anuga import *
     12#from anuga import *
    1313#from swb_domain import domain
    1414#from anuga import *
     
    8282
    8383
    84 domain.set_name('channel_floodplain1_bal_dev_lowvisc') # Output name
     84domain.set_name('channel_floodplain1') # Output name
    8585#domain.set_store_vertices_uniquely(True)
    8686#domain.use_edge_limiter=False
     
    142142
    143143# Define inlet operator
    144 flow_in_yval=100.0
     144flow_in_yval=10.0
    145145if True:
    146146    line1 = [ [floodplain_width/2. - chan_width/2., flow_in_yval],\
     
    245245    if( numpy.floor(t/100.) == t/100. ):
    246246        print '#### COMPUTING FLOW THROUGH CROSS-SECTIONS########'
     247        s0 = domain.get_flow_through_cross_section([[0., 10.0], [floodplain_width, 10.]])
    247248        s1 = domain.get_flow_through_cross_section([[0., floodplain_length-300.0], [floodplain_width, floodplain_length-300.0]])
    248249        s2 = domain.get_flow_through_cross_section([[0., floodplain_length-1.0], [floodplain_width, floodplain_length-1.0]])
    249250       
    250         print 'Cross sectional flows: ', s1, s2
     251        print 'Cross sectional flows: ', s0, s1, s2
    251252        print '$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
    252253
  • trunk/anuga_work/development/gareth/tests/channel_floodplain/plotme.py

    r8384 r8446  
    44
    55# Time-index to plot outputs from
    6 index=53
     6index=60
    77
    88#p2 = util.get_output('channel_floodplain1_bal_dev.sww', minimum_allowed_height=0.01)
     
    1010#p=util.get_centroids(p2, velocity_extrapolation=True)
    1111
    12 p2 = util.get_output('channel_floodplain1_bal_dev_lowvisc.sww', 0.01)
     12#p2 = util.get_output('channel_floodplain1_bal_dev_lowvisc.sww', 0.01)
     13p2 = util.get_output('channel_floodplain1.sww', 0.01)
    1314p=util.get_centroids(p2, velocity_extrapolation=True)
    1415
  • trunk/anuga_work/development/gareth/tests/dam_break/plotme.py

    r8384 r8446  
    99p2_st=util.get_centroids(p_st)
    1010
    11 
    12 p_dev = util.get_output('dam_break_20120404_230353/dam_break.sww', 0.001)
     11p_dev = util.get_output('dam_break_20120517_003225/dam_break.sww', 0.001)
    1312p2_dev=util.get_centroids(p_dev, velocity_extrapolation=True)
    1413
  • trunk/anuga_work/development/gareth/tests/runup/runup.py

    r8353 r8446  
    6767#------------------------------
    6868
    69 #xwrite=open("xvel.out","wb")
    70 #ywrite=open("yvel.out","wb")
    71 ## Set print options to be compatible with file writing via the 'print' statement
    72 #numpy.set_printoptions(threshold=numpy.nan, linewidth=numpy.nan)
    73 
    74 for t in domain.evolve(yieldstep=0.2,finaltime=30.0):
     69for t in domain.evolve(yieldstep=0.2,finaltime=60.0):
    7570    print domain.timestepping_statistics()
    76 
     71    uh=domain.quantities['xmomentum'].centroid_values
     72    vh=domain.quantities['ymomentum'].centroid_values
     73    depth=domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values
     74    depth=depth*(depth>1.0e-06) + 1.0e-06
     75    vel=((uh/depth)**2 + (vh/depth)**2)**0.5
     76    print 'peak speed is', vel.max()
    7777
    7878print 'Finished'
  • trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py

    r8359 r8446  
    6565Bt=anuga.Transmissive_boundary(domain)          # Continue all values of boundary -- not used in this example
    6666Bd=anuga.Dirichlet_boundary([-0.1*scale_me,0.,0.])       # Constant boundary values -- not used in this example
     67def waveform(t):
     68    return (0.0*sin(t*2*pi)-0.1)*exp(-t)-0.1
     69Bt2=anuga.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain,waveform)
    6770#Bw=anuga.Time_boundary(domain=domain,
    6871#       f=lambda t: [(0.0*sin(t*2*pi)-0.1)*exp(-t)-0.1,0.0,0.0]) # Time varying boundary -- get rid of the 0.0 to do a runup.
     
    8285    yy = domain.quantities['ymomentum'].centroid_values
    8386    dd = domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values
    84     dd = (dd)*(dd>1.0e-03)+1.0e-06
     87    dd = (dd)*(dd>1.0e-06)+1.0e-03
    8588    vv = ( (xx/dd)**2 + (yy/dd)**2 )**0.5
    8689    vv = vv*(dd>1.0e-03)
  • trunk/anuga_work/development/gareth/tests/wave/run_wave.py

    r8396 r8446  
    8181Br = anuga.Reflective_boundary(domain)      # Solid reflective wall
    8282Bt = anuga.Transmissive_boundary(domain)    # Continue all values on boundary
    83 Bz = swb2_boundary_conditions.zero_mass_flux_zero_momentum_boundary(domain) # Strong reflections
    84 Bz1 = swb2_boundary_conditions.zero_mass_flux_zero_t_transmissive_n_momentum_boundary(domain) # Strong reflections
    85 Bz2 = swb2_boundary_conditions.zero_mass_flux_zero_n_transmissive_t_momentum_boundary(domain) # Strong reflections
     83#Bz = swb2_boundary_conditions.zero_mass_flux_zero_momentum_boundary(domain) # Strong reflections
     84#Bz1 = swb2_boundary_conditions.zero_mass_flux_zero_t_transmissive_n_momentum_boundary(domain) # Strong reflections
     85#Bz2 = swb2_boundary_conditions.zero_mass_flux_zero_n_transmissive_t_momentum_boundary(domain) # Strong reflections
    8686Bs = anuga.Transmissive_stage_zero_momentum_boundary(domain) # Strong reflections
    8787Bd = anuga.Dirichlet_boundary([1,0.,0.]) # Constant boundary values
     
    100100    return -amplitude*sin((1./wave_length)*t*2*pi)
    101101
     102def waveform2(t):
     103    return 0.
     104
    102105Bw2 = anuga.shallow_water.boundaries.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain, waveform)
     106Bf_in = swb2_boundary_conditions.Characteristic_boundary(domain, waveform)
    103107#Bw3 = swb2_boundary_conditions.Transmissive_momentum_nudge_stage_boundary(domain, waveform)
     108Bf=swb2_boundary_conditions.Characteristic_boundary(domain,waveform2)
    104109                   
    105110
    106111
    107112# Associate boundary tags with boundary objects
    108 domain.set_boundary({'left': Bw2, 'right': Bt, 'top': Br, 'bottom': Br})
     113domain.set_boundary({'left': Bw2, 'right': Bf, 'top': Br, 'bottom': Br})
    109114
    110115
Note: See TracChangeset for help on using the changeset viewer.