Changeset 8547


Ignore:
Timestamp:
Aug 31, 2012, 8:38:25 PM (13 years ago)
Author:
davies
Message:

Major experimental changes to balanced_dev

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

Legend:

Unmodified
Added
Removed

  • trunk/anuga_work/development/gareth/experimental/balanced_dev/runup.py

    r8353 r8547  
    1515
    1616from math import sin, pi, exp
     17from balanced_dev import *
     18#from anuga_tsunami import *
    1719#from anuga.shallow_water_balanced2.swb2_domain import Domain as Domain
    1820#from anuga.shallow_water.shallow_water_domain import Domain as Domain
    … …  
    2325#from swb_domain import *
    2426#path.append('/home/gareth/storage/anuga_clean/anuga_jan12/trunk/anuga_work/development/gareth/balanced_basic')
    25 from balanced_basic.swb2_domain import *
     27#from balanced_basic.swb2_domain import *
    2628#---------
    2729#Setup computational domain
    … …  
    6062Bt=anuga.Transmissive_boundary(domain)          # Continue all values of boundary -- not used in this example
    6163Bd=anuga.Dirichlet_boundary([-0.2,0.,0.])       # Constant boundary values -- not used in this example
    62 Bw=anuga.Time_boundary(domain=domain,
    63         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.
     64Bw=anuga.Time_boundary(domain=domain, function=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.
    6465
    6566#----------------------------------------------

  • trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain.py

    r8446 r8547  
    102102
    103103        self.maximum_allowed_speed=0.0
    104 
     104        #self.forcing_terms.append(manning_friction_explicit)
     105        #self.forcing_terms.remove(manning_friction_implicit)
    105106
    106107        print '##########################################################################'
    … …  
    120121        print "#"
    121122        print "# 2) In field scale applications (where the depth is typically > 1m), I suggest you set"
    122         print "# domain.minimum_allowed_height=0.01 (the default is 1.0e-3). With this setting, velocity"
    123         print "# seems to be very well behaved"
     123        print "# domain.minimum_allowed_height=0.01 (the default is 1.0e-3). "
    124124        print "#"
    125125        print "# 3) This solver is not expected to perform well in problems with very"
    … …  
    290290        xmomc = domain.quantities['xmomentum'].centroid_values
    291291        ymomc = domain.quantities['ymomentum'].centroid_values
     292        areas = domain.areas
     293        xc = domain.centroid_coordinates[:,0]
     294        yc = domain.centroid_coordinates[:,1]
    292295
    293296        protect(domain.minimum_allowed_height, domain.maximum_allowed_speed,
    294                 domain.epsilon, wc, wv, zc,zv, xmomc, ymomc)
     297                domain.epsilon, wc, wv, zc,zv, xmomc, ymomc, areas, xc, yc)
    295298
    296299    def conserved_values_to_evolved_values(self, q_cons, q_evol):
    … …  
    416419        #    Q = domain.quantities[name]
    417420        #    Q.interpolate_from_edges_to_vertices()
    418         #    #Q.interpolate_from_vertices_to_edges()
     421        #    Q.interpolate_from_vertices_to_edges()
    419422
    420423
    … …  
    521524        extrapol2(self,
    522525                  self.surrogate_neighbours,
     526                  self.neighbour_edges,
    523527                  self.number_of_boundaries,
    524528                  self.centroid_coordinates,

  • trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c

    r8466 r8547  
    6767    // Apply limiting of speeds according to the ANUGA manual
    6868    if (*h < epsilon) {
    69       //*h = max(0.0,*h);  // Could have been negative
    70       *h = 0.0;  // Could have been negative
     69      *h = max(0.0,*h);  // Could have been negative
     70      //*h = 0.0;  // Could have been negative
    7171      u = 0.0;
    7272    } else {
    … …  
    196196  // Flux formulas
    197197  flux_left[0] = u_left*h_left;
    198   flux_left[1] = u_left*uh_left + 0.5*g*h_left*h_left;
     198  flux_left[1] = u_left*uh_left + 0.0*g*h_left*h_left;
    199199  flux_left[2] = u_left*vh_left;
    200200
    201201  flux_right[0] = u_right*h_right;
    202   flux_right[1] = u_right*uh_right + 0.5*g*h_right*h_right;
     202  flux_right[1] = u_right*uh_right + 0.0*g*h_right*h_right;
    203203  flux_right[2] = u_right*vh_right;
    204204
    … …  
    208208  { // FIXME (Ole): Try using h0 here
    209209    memset(edgeflux, 0, 3*sizeof(double));
    210     //for(i=0;i<3;i++){
    211     //    edgeflux[i] = _minmod(flux_left[i], flux_right[i]);
    212     //}
     210
    213211    *max_speed = 0.0;
     212    *pressure_flux = 0.0;
    214213  }
    215214  else
    … …  
    218217    for (i = 0; i < 3; i++)
    219218    {
    220       // Adjustment to the scheme by Kurganov and Lin 2007 Communications in Computational
    221       // Physics 2:141-163
    222       //uint = (s_max*q_right_rotated[i] - s_min*q_left_rotated[i] - (flux_right[i] - flux_left[i]))*inverse_denominator;
    223       //t1 = (q_right_rotated[i] - uint);
    224       //t2 = (-q_left_rotated[i] + uint);
    225       //t3 = _minmod(t1,t2);
    226       t3 = 0.0;
    227219      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
    228       //if(i<2) {
    229       edgeflux[i] += s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3);
    230       //}else{
    231       //  edgeflux[i] += 0.5*s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3);
    232       //}
    233       //t1 = s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3);
    234       //if(fabs(edgeflux[i])>fabs(t1)){
    235       //      edgeflux[i]+=t1;
    236       //}else{
    237       //      edgeflux[i]+=fabs(edgeflux[i])*copysign(1.0,t1) ;
    238       //}
    239 
     220      edgeflux[i] += s_max*s_min*(q_right_rotated[i] - q_left_rotated[i]);
    240221      edgeflux[i] *= inverse_denominator;
    241222    }
    242223
    243224    *pressure_flux = 0.5*g*( s_max*h_left*h_left -s_min*h_right*h_right)*inverse_denominator;
    244     //if(edgeflux[0]<0.0){
    245     //    //edgeflux[2] = edgeflux[0]*0.5*(vh_right/h_right + vh_left/h_left);
    246     //    edgeflux[2] = edgeflux[0]*vh_right/h_right;
    247     //}else{
    248     //    edgeflux[2] = edgeflux[0]*vh_left/h_left;
    249     //}
    250 
     225   
    251226    // Maximal wavespeed
    252227    *max_speed = max(fabs(s_max), fabs(s_min));
    … …  
    306281    int neighbours_wet[3];//Work array
    307282    int useint;
    308     double stage_edge_lim, scale_factor_shallow, bedtop, bedbot, pressure_flux, hc, hc_n;
    309     //double min_bed_edgevalue[number_of_elements]; // FIXME: Set memory explicitly
     283    double stage_edge_lim, scale_factor_shallow, bedtop, bedbot, pressure_flux, hc, hc_n, tmp;
     284    //double max_bed_edgevalue[number_of_elements]; // FIXME: Set memory explicitly
    310285    static long call = 1; // Static local variable flagging already computed flux
    311286
    312     double *min_bed_edgevalue;
    313 
     287    double *max_bed_edgevalue, *min_bed_edgevalue;
     288
     289    max_bed_edgevalue = malloc(number_of_elements*sizeof(double));
    314290    min_bed_edgevalue = malloc(number_of_elements*sizeof(double));
    315291    // Start computation
    … …  
    325301    // Compute minimum bed edge value on each triangle
    326302    for (k = 0; k < number_of_elements; k++){
     303        max_bed_edgevalue[k] = max(bed_edge_values[3*k],
     304                                   max(bed_edge_values[3*k+1], bed_edge_values[3*k+2]));
     305        //max_bed_edgevalue[k]= 0.5*(max_bed_edgevalue[j]+bed_centroid_values[k]);
    327306        min_bed_edgevalue[k] = min(bed_edge_values[3*k],
    328307                                   min(bed_edge_values[3*k+1], bed_edge_values[3*k+2]));
    329         //min_bed_edgevalue[k] = bed_centroid_values[k];
    330308   
    331309    }
    … …  
    349327            zl = bed_edge_values[ki];
    350328            hc = max(stage_centroid_values[k] - bed_centroid_values[k],0.0);
     329
     330            //if(ql[0]<=zl){
     331            //    ql[1]=0.;
     332            //    ql[2]=0.;
     333            //}
    351334            // Get right hand side values either from neighbouring triangle
    352335            // or from boundary array (Quantities at neighbour on nearest face).
    … …  
    373356                zr = bed_edge_values[nm];
    374357            }
    375 
    376             // Now we have values for this edge - both from left and right side.
    377 
    378             if (optimise_dry_cells) {
    379                 // Check if flux calculation is necessary across this edge
    380                 // This check will exclude dry cells.
    381                 // This will also optimise cases where zl != zr as
    382                 // long as both are dry
    383 
    384                 if ((ql[0] - zl) < epsilon &&
    385                         (qr[0] - zr) < epsilon) {
    386                     // Cell boundary is dry
    387 
    388                     already_computed_flux[ki] = call; // #k Done
    389                     if (n >= 0) {
    390                         already_computed_flux[nm] = call; // #n Done
    391                     }
    392 
    393                     max_speed = 0.0;
    394                     continue;
    395                 }
    396             }
    397 
     358           
    398359
    399360            if (fabs(zl-zr)>1.0e-10) {
    … …  
    402363            }
    403364           
    404             // If both centroids are dry, then the cells should not exchange flux
    405             // NOTE: This also means no bedslope term -- which is arguably
    406             // appropriate, since the depth is zero at the cell centroid
    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;
    411                 continue;
    412             }
    413            
    414365            //If one centroid is dry, then extrapolate its edge values from the neighbouring centroid,
    415366            // unless the local centroid value is smaller
    416367            if(n>=0){
    417                 if(hc==0.0){
     368                //if(hc==0.0){
     369                if((hc<=H0)&&(hc_n>H0)){
     370                    ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
     371                    //ql[0]=max(min(qr[0],0.5*(stage_centroid_values[k]+stage_centroid_values[n])),zl);
    418372                    //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);
    421373                }
    422                 if(hc_n==0.0){
    423                     //qr[0]=max(min(ql[0],stage_centroid_values[n]),zr);
     374                //if(hc_n==0.0){
     375                if((hc_n<=H0)&&(hc>H0)){
     376                    qr[0]=max(min(ql[0],stage_centroid_values[n]),zr);
    424377                    //qr[0]=max(min(ql[0],0.5*(stage_centroid_values[n]+stage_centroid_values[k])),zr);
    425                     qr[0]=ql[0];
     378                    //qr[0]=ql[0];
    426379                }
     380
     381                //if((hc_n<=H0)&&(hc<=H0)){
     382                //    ql[0] = zl;
     383                //    qr[0] = zr;
     384                //}
     385
    427386            }else{
    428387                // Treat the boundary case
    429                 if((hc==0.0)){
    430                     ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
    431                     //ql[0]=max(min(qr[0],ql[0]),zl);
    432                 }
     388                //if((hc==0.0)){
     389                //    ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);
     390                //    //ql[0]=max(min(qr[0],ql[0]),zl);
     391                //}
    433392            }
     393            //if(k==229 | n==229){
     394            //    printf(" ---------------\n ");
     395            //    printf("k, ql, qr, zl, zr==%d, %20.20e, %20.20e, %20.20e, %20.20e \n", k, ql[0], qr[0], zl, zr);
     396
     397            //}
    434398           
    435399            // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2])
    … …  
    443407
    444408            // Prevent outflow from 'seriously' dry cells
    445             if((stage_centroid_values[k]<=min_bed_edgevalue[k])|
     409            // Idea: The cell will not go dry if:
     410            // mass_flux = edgeflux[0]*(dt*edgelength) <= Area_triangle*hc
     411            if((stage_centroid_values[k]<=max_bed_edgevalue[k])|
    446412               (ql[0]<=zl)){
    447413                if(edgeflux[0]>0.0){
    448                     edgeflux[0]=0.0;
    449                     edgeflux[1]=0.0;
    450                     edgeflux[2]=0.0;
     414                    tmp=min((1.0/3.0)*areas[k]*hc/(edgelengths[ki]*max(timestep,epsilon)), 1.0); // 28/7 -- Right answer for channel flow problem.
     415                    tmp = min(min(edgeflux[0], tmp)/edgeflux[0], 1.0);
     416                    if((tmp<0.)|tmp>1.0){
     417                        printf("Invalid tmp value\n");
     418                    }
     419                    edgeflux[0]*=tmp;
     420                    edgeflux[1]*=tmp;
     421                    edgeflux[2]*=tmp;
    451422                }
    452423            }
    453424            if(n>=0){
    454                 if((stage_centroid_values[n]<=min_bed_edgevalue[n])|
     425                if((stage_centroid_values[n]<=max_bed_edgevalue[n])|
    455426                   (qr[0]<=zr)){
    456427                    if(edgeflux[0]<0.0){
    457                         edgeflux[0]=0.0;
    458                         edgeflux[1]=0.0;
    459                         edgeflux[2]=0.0;
     428                        tmp=min((1.0/3.0)*areas[n]*hc_n/(edgelengths[ki]*max(timestep,epsilon)), 1.0); // 28/7 -- Right answer for channel flow problem.
     429                        tmp = min( max(edgeflux[0], -tmp)/edgeflux[0], 1.0);
     430                        if((tmp<0.)|tmp>1.0){
     431                            printf("Invalid tmp value\n");
     432                        }
     433                        edgeflux[0]*=tmp;
     434                        edgeflux[1]*=tmp;
     435                        edgeflux[2]*=tmp;
    460436                    }
    461437                }
    462438            }
     439           
     440            //if(k==229|n==229){
     441            //    printf("++edgeflux=, %20.20e, %20.20e, %20.20e \n", edgeflux[0],edgeflux[1], edgeflux[2]);
     442            //}
    463443
    464444            // Multiply edgeflux by edgelength
    … …  
    468448            edgeflux[2] *= length;
    469449
    470             // GET RID OF THIS: EXPERIMENTAL
    471             //if(n<0){
    472             //    if(boundary_flux_type[m]>0){
    473             //        // IMPOSE BOUNDARY CONDITIONS DIRECTLY ON THE FLUX
    474             //        if(boundary_flux_type[m]==1){
    475             //            // Zero_mass_flux_transmissive_momentum_flux boundary, or
    476             //            // zero_mass_flux_zero_momentum_boundary
    477             //            // Just need to set the mass flux to zero, as the
    478             //            // transmissive momentum is ensured by the values of
    479             //            // qr[0], qr[1], qr[2]
    480             //            printf("Warning: WEIRD EDGEFLUX THING IS ACTING");
    481             //            edgeflux[0] = 0.0;
    482             //        }
    483             //    }
    484             //}
    485450
    486451            // Update triangle k with flux from edge i
    … …  
    490455
    491456            // Compute bed slope term
    492             if(hc>-9999.0){
     457            //if(hc>-9999.0){
    493458                //Bedslope approx 1:
    494                 //bedslope_work = g*hc*(ql[0])*length-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length;
    495                 bedslope_work = g*length*( hc*(ql[0])-0.5*max(ql[0]-zl,0.)*(ql[0]-zl) );
     459            if(hc>H0){
     460                // Add the pressure gradient flux term back in
     461                // Do this to avoid the above mass-conservative limiting
     462                // Because limiting the pressure term would destroy balancing with the bedslope term
     463                //edgeflux[1]+=normals[ki2]*pressure_flux;
     464                //edgeflux[2]+=normals[ki2+1]*pressure_flux;
     465
     466                //bedslope_work = -0.5*g*max(stage_centroid_values[k]-zl,0.)*(stage_centroid_values[k]-zl)*length;
     467                bedslope_work = length*g*(hc*ql[0]-0.5*max(ql[0]-zl,0.)*(ql[0]-zl)) + pressure_flux*length;
     468                xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
     469                ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
     470            }else{
     471                bedslope_work = 0.;//-pressure_flux*length ;//g*length*( 0.*(ql[0])-0.5*max(ql[0]-zl,0.)*(ql[0]-zl) );
     472                xmom_explicit_update[k] *= 0.;
     473                ymom_explicit_update[k] *= 0.;
     474            }
     475           
     476            //if(k==229|n==229){
     477            //    printf("--k-pres, xmom, bedwrk, %20.20e, %20.20e, %20.20e \n", pressure_flux*length, xmom_explicit_update[k], bedslope_work);
     478            ////    printf(" Edgesum: %20.20e\n", edgelengths[3*k]*normals[6*k]*g*hc*stage_edge_values[3*k]
     479            ////                   + edgelengths[3*k+1]*normals[6*k+2]*g*hc*stage_edge_values[3*k+1]
     480            ////                   + edgelengths[3*k+2]*normals[6*k+4]*g*hc*stage_edge_values[3*k+2]);
     481            //}
    496482                //
    497483                // Bedslope approx 2
    … …  
    508494                //bedslope_work = -0.5*g*max(stage_centroid_values[k]-zl,0.)*(stage_centroid_values[k]-zl)*length;
    509495                //
    510                 xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
    511                 ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
    512             }else{
    513                 // Treat nearly dry cells
    514                 bedslope_work =-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length; //
    515                 //
    516                 //bedslope_work = -pressure_flux*length;
    517                 xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
    518                 ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
    519 
    520                 //xmom_explicit_update[k] = 0.0;
    521                 //ymom_explicit_update[k] = 0.0;
    522 
    523             }
     496            //}else{
     497            //   // Treat nearly dry cells
     498            //   bedslope_work =-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length; //
     499            //   //
     500            //   //bedslope_work = -pressure_flux*length;
     501            //   xmom_explicit_update[k] -= normals[ki2]*bedslope_work;
     502            //   ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work;
     503            //   //xmom_explicit_update[k] = 0.0;
     504            //   //ymom_explicit_update[k] = 0.0;
     505
     506            ////}
    524507
    525508            already_computed_flux[ki] = call; // #k Done
    … …  
    532515                ymom_explicit_update[n] += edgeflux[2];
    533516                //Add bed slope term here
    534                 if(hc_n>-9999.0){
     517                //if(hc_n>-9999.0){
    535518                //if(stage_centroid_values[n] > max_bed_edgevalue[n]){
    536519                    // Bedslope approx 1:
    537                     //bedslope_work = g*hc_n*(qr[0])*length-0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length;
    538                     bedslope_work = g*length*(hc_n*(qr[0])-0.5*max(qr[0]-zr,0.)*(qr[0]-zr));
     520                if(hc_n>H0){
     521                    //bedslope_work = -0.5*g*max(stage_centroid_values[n]-zr,0.)*(stage_centroid_values[n]-zr)*length;
     522                    bedslope_work = length*g*(hc_n*qr[0]-0.5*max(qr[0]-zr,0.)*(qr[0]-zr)) + pressure_flux*length;
     523                    xmom_explicit_update[n] += normals[ki2]*bedslope_work;
     524                    ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
     525                }else{
     526                    bedslope_work = 0.;//-pressure_flux*length;//g*length*(0.*(qr[0])-0.5*max(qr[0]-zr,0.)*(qr[0]-zr));
     527                    xmom_explicit_update[n] *= 0.;
     528                    ymom_explicit_update[n] *= 0.;
     529                }
     530                //if(k==229| n==229){
     531                //    printf("-n-pres,xmom,bedwrk, %20.20e, %20.20e, %20.20e \n", pressure_flux*length, xmom_explicit_update[n], bedslope_work);
     532                //}
     533
    539534                    //
    540535                    // Bedslope approx 2:
    … …  
    551546                    //bedslope_work = -0.5*g*max(stage_centroid_values[n]-zr,0.)*(stage_centroid_values[n]-zr)*length;
    552547                    //
    553                     xmom_explicit_update[n] += normals[ki2]*bedslope_work;
    554                     ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
    555                 }else{
    556                     // Treat nearly dry cells
    557                     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;
    558                     //bedslope_work = -pressure_flux*length;
    559                     xmom_explicit_update[n] += normals[ki2]*bedslope_work;
    560                     ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
    561 
    562                     //xmom_explicit_update[n] = 0.0;
    563                     //ymom_explicit_update[n] = 0.0;
    564                 }
     548                //}else{
     549                //    // Treat nearly dry cells
     550                //    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;
     551                //    //bedslope_work = -pressure_flux*length;
     552                //    xmom_explicit_update[n] += normals[ki2]*bedslope_work;
     553                //    ymom_explicit_update[n] += normals[ki2+1]*bedslope_work;
     554
     555                //    //xmom_explicit_update[n] = 0.0;
     556                //    //ymom_explicit_update[n] = 0.0;
     557                //}
    565558
    566559                already_computed_flux[nm] = call; // #n Done
    … …  
    591584        } // End edge i (and neighbour n)
    592585
     586        //if(k==2082){
     587        //    printf("xmom_explicit: %20.20e \n", xmom_explicit_update[k]);
     588        //}
    593589
    594590        // Normalise triangle k by area and store for when all conserved
    … …  
    599595        ymom_explicit_update[k] *= inv_area;
    600596
    601 
     597        //if(k==2082){
     598        //    printf(" inv_area: %20.20e \n", inv_area);
     599        //    printf("xmom_explicit: %20.20e \n", xmom_explicit_update[k]);
     600        //}
    602601        // Keep track of maximal speeds
    603602        max_speed_array[k] = max_speed;
    604603
    605604    } // End triangle k
    606 
     605 
     606     
     607    //bedtop=0.;
     608    //for(k=0; k<number_of_elements; k++){
     609    //      bedtop = bedtop + stage_explicit_update[k]*areas[k];
     610    //}
     611    //printf("Bedtop = %20.20e \n ", bedtop);
     612    //    for(i=0;i<number_of_elements;i++){
     613    //        bedtop=max(bedtop, fabs(ymom_edge_values[3*k+i]));
     614    //    }
     615    //   
     616    //   
     617    //    if((fabs(xmom_explicit_update[k])>1.0e-06)|fabs(ymom_explicit_update[k])>1.0e-06){
     618    //        printf(" %e, %e, %e,%e,%e,%e, %e,%e,%e,%d,%e,%e,%e\n", stage_explicit_update[k], xmom_explicit_update[k], ymom_explicit_update[k],
     619    //                                    stage_centroid_values[k], bed_centroid_values[k],
     620    //                                    stage_centroid_values[k]- bed_centroid_values[k],
     621    //                                    stage_edge_values[3*neighbours[3*k+0]+neighbour_edges[3*k+0]],
     622    //                                    stage_edge_values[3*neighbours[3*k+1]+neighbour_edges[3*k+1]],
     623    //                                    stage_edge_values[3*neighbours[3*k+2]+neighbour_edges[3*k+2]],
     624    //                                    k,
     625    //                                    stage_edge_values[3*k+0],
     626    //                                    stage_edge_values[3*k+1],
     627    //                                    stage_edge_values[3*k+2]);
     628    //    }
     629    //}
     630    //printf("max ymom edge value is: %e \n", bedtop);
     631    //report_python_error(AT, "Written out stage, xmom and ymom updates");
     632    //return -1;
     633//
     634    free(max_bed_edgevalue);
    607635    free(min_bed_edgevalue);
    608636
    … …  
    620648         double* zv,
    621649         double* xmomc,
    622          double* ymomc) {
     650         double* ymomc,
     651         double* areas,
     652         double* xc,
     653         double* yc) {
    623654
    624655  int k;
    625656  double hc, bmin, bmax;
    626657  double u, v, reduced_speed;
    627 
     658  static double mass_error=0.;
    628659  // This acts like minimum_allowed height, but scales with the vertical
    629660  // distance between the bed_centroid_value and the max bed_edge_value of
    630661  // every triangle.
    631   double minimum_relative_height=0.1;
    632 
     662  double minimum_relative_height=0.0;
     663  int mass_added=0;
    633664
    634665  // Protect against inifintesimal and negative heights 
    635   if (maximum_allowed_speed < epsilon) {
     666  //if (maximum_allowed_speed < epsilon) {
    636667    for (k=0; k<N; k++) {
    637668      hc = wc[k] - zc[k];
    … …  
    643674        // Set momentum to zero and ensure h is non negative
    644675        // NOTE: THIS IS IMPORTANT -- WE ARE SETTING MOMENTUM TO ZERO
    645         xmomc[k] = 0.0;
    646         ymomc[k] = 0.0;
     676        //if(hc<=epsilon){
     677            xmomc[k] = 0.0;
     678            ymomc[k] = 0.0;
     679        //}
    647680
    648681        if (hc <= 0.0){
    … …  
    652685             //bmin =0.5*bmin + 0.5*min(zv[3*k],min(zv[3*k+1],zv[3*k+2]));
    653686             // Setting = minimum vertex value seems better, but might tend to be less smooth
    654              bmin =min(zv[3*k],min(zv[3*k+1],zv[3*k+2])) -minimum_allowed_height;
     687             //bmin =min(zv[3*k],min(zv[3*k+1],zv[3*k+2])) -minimum_allowed_height;
     688             bmin=zc[k]-minimum_allowed_height;
    655689             // Minimum allowed stage = bmin
    656690             // WARNING: ADDING MASS if wc[k]<bmin
    657691             if(wc[k]<bmin){
    658                  printf("Adding mass to dry cell %d %f %f %f %f \n", k, zv[3*k], zv[3*k+1], zv[3*k+2], wc[k]);
    659              }
    660              wc[k] = max(wc[k], bmin);
    661 
    662              // Set vertex values as well. Seems that this shouldn't be
    663              // needed. However from memory this is important at the first
    664              // time step, for 'dry' areas where the designated stage is
    665              // less than the bed centroid value
    666              wv[3*k] = max(wv[3*k], bmin);
    667              wv[3*k+1] = max(wv[3*k+1], bmin);
    668              wv[3*k+2] = max(wv[3*k+2], bmin);
     692                 mass_error+=(bmin-wc[k])*areas[k];
     693                 mass_added=1; //Flag to warn of added mass               
     694
     695                 wc[k] = max(wc[k], bmin);
     696                 printf(" mass_add, %f, %f, %f,%d\n", xc[k], yc[k], wc[k], k) ;
     697             
     698
     699                 // Set vertex values as well. Seems that this shouldn't be
     700                 // needed. However from memory this is important at the first
     701                 // time step, for 'dry' areas where the designated stage is
     702                 // less than the bed centroid value
     703                 //wv[3*k] = max(wv[3*k], min(bmin, zv[3*k]-minimum_allowed_height));
     704                 //wv[3*k+1] = max(wv[3*k+1], min(bmin, zv[3*k+1]-minimum_allowed_height));
     705                 //wv[3*k+2] = max(wv[3*k+2], min(bmin, zv[3*k+2]-minimum_allowed_height));
     706                 wv[3*k] = min(bmin, wc[3*k]); //zv[3*k]-minimum_allowed_height);
     707                 wv[3*k+1] = min(bmin, wc[3*k+1]); //zv[3*k+1]-minimum_allowed_height);
     708                 wv[3*k+2] = min(bmin, wc[3*k+2]); //zv[3*k+2]-minimum_allowed_height);
     709            }
    669710        }
    670711      }
    671712    }
    672   } else {
    673    
    674     // Protect against infinitesimal and negative heights
    675     for (k=0; k<N; k++) {
    676       hc = wc[k] - zc[k];
    677      
    678       if (hc < minimum_allowed_height) {
    679 
    680         //New code: Adjust momentum to guarantee speeds are physical
    681         //          ensure h is non negative
    682              
    683         if (hc <= 0.0) {
    684             wc[k] = zc[k];
    685         xmomc[k] = 0.0;
    686         ymomc[k] = 0.0;
    687         } else {
    688           //Reduce excessive speeds derived from division by small hc
    689           //FIXME (Ole): This may be unnecessary with new slope limiters
    690           //in effect.
    691          
    692           u = xmomc[k]/hc;
    693           if (fabs(u) > maximum_allowed_speed) {
    694             reduced_speed = maximum_allowed_speed * u/fabs(u);
    695             //printf("Speed (u) has been reduced from %.3f to %.3f\n",
    696             //   u, reduced_speed);
    697             xmomc[k] = reduced_speed * hc;
    698           }
    699          
    700           v = ymomc[k]/hc;
    701           if (fabs(v) > maximum_allowed_speed) {
    702             reduced_speed = maximum_allowed_speed * v/fabs(v);
    703             //printf("Speed (v) has been reduced from %.3f to %.3f\n",
    704             //   v, reduced_speed);
    705             ymomc[k] = reduced_speed * hc;
    706           }
    707         }
    708       }
    709     }
     713
     714  if(mass_added==1){
     715     printf("Cumulative mass protection: %f m^3 \n", mass_error);
    710716  }
    711717  return 0;
    … …  
    767773}
    768774
    769 //int limit_gradient2(double *dqv, double qmin, double qmax, double beta_w, double r0scale){
    770 //  // EXPERIMENTAL LIMITER, DOESN'T WORK
    771 //  // Given provisional jumps dqv from the FV triangle centroid to its
    772 //  // vertices and jumps qmin (qmax) between the centroid of the FV
    773 //  // triangle and the minimum (maximum) of the values at the centroid of
    774 //  // the FV triangle and the auxiliary triangle vertices,
    775 //  // calculate a multiplicative factor phi by which the provisional
    776 //  // vertex jumps are to be limited
    777 // 
    778 //  int i;
    779 //  double r=1000.0, r0=1.0, phi=1.0;
    780 //  static double TINY = 1.0e-100; // to avoid machine accuracy problems.
    781 //  // FIXME: Perhaps use the epsilon used elsewhere.
    782 // 
    783 //  for (i=0;i<3;i++){
    784 //    if (dqv[i]<-TINY)
    785 //      r0=qmin/dqv[i]*r0scale*2.0;
    786 //     
    787 //    if (dqv[i]>TINY)
    788 //      r0=qmax/dqv[i]*r0scale*2.0;
    789 //     
    790 //    r=min(r0,r);
    791 //  }
    792 // 
    793 //  phi=min(r*beta_w,1.0);
    794 //  //phi=1.;
    795 //  //for (i=0;i<3;i++)
    796 //  dqv[0]=dqv[0]*phi;
    797 //  dqv[1]=dqv[1]*phi;
    798 //  dqv[2]=dqv[2]*phi;
    799 //   
    800 //  return 0;
    801 //}
    802 
    803775// Computational routine
    804776int _extrapolate_second_order_edge_sw(int number_of_elements,
    … …  
    812784                                 double beta_vh_dry,
    813785                                 long* surrogate_neighbours,
     786                                 long* neighbour_edges,
    814787                                 long* number_of_boundaries,
    815788                                 double* centroid_coordinates,
    … …  
    832805  // Local variables
    833806  double a, b; // Gradient vector used to calculate edge values from centroids
    834   int k, k0, k1, k2, k3, k6, coord_index, i, ii;
     807  int k, k0, k1, k2, k3, k6, coord_index, i, ii, ktmp, k_wetdry;
    835808  double x, y, x0, y0, x1, y1, x2, y2, xv0, yv0, xv1, yv1, xv2, yv2; // Vertices of the auxiliary triangle
    836809  double dx1, dx2, dy1, dy2, dxv0, dxv1, dxv2, dyv0, dyv1, dyv2, dq0, dq1, dq2, area2, inv_area2;
    837810  double dqv[3], qmin, qmax, hmin, hmax, bedmax, stagemin;
    838   double hc, h0, h1, h2, beta_tmp, hfactor;
    839   double dk, dv0, dv1, dv2, de[3], demin, dcmax, r0scale;
    840  
    841   double *xmom_centroid_store, *ymom_centroid_store, *stage_centroid_store;
     811  double hc, h0, h1, h2, beta_tmp, hfactor, xtmp, ytmp, weight;
     812  double dk, dv0, dv1, dv2, de[3], demin, dcmax, r0scale, vect_norm, l1, l2;
     813 
     814  double *xmom_centroid_store, *ymom_centroid_store, *stage_centroid_store, *min_elevation_edgevalue, *max_elevation_edgevalue;
     815  double *neighbour_wetness_scale;
     816  int *count_wet_neighbours;
    842817
    843818  // Use malloc to avoid putting these variables on the stack, which can cause
    … …  
    846821  ymom_centroid_store = malloc(number_of_elements*sizeof(double));
    847822  stage_centroid_store = malloc(number_of_elements*sizeof(double));
     823  min_elevation_edgevalue = malloc(number_of_elements*sizeof(double));
     824  max_elevation_edgevalue = malloc(number_of_elements*sizeof(double));
     825  neighbour_wetness_scale = malloc(number_of_elements*sizeof(double));
     826  count_wet_neighbours = malloc(number_of_elements*sizeof(int));
    848827 
    849828  if(extrapolate_velocity_second_order==1){
    … …  
    859838          ymom_centroid_values[k] = ymom_centroid_values[k]/dk;
    860839
     840          min_elevation_edgevalue[k] = min(elevation_edge_values[3*k],
     841                                           min(elevation_edge_values[3*k+1],
     842                                               elevation_edge_values[3*k+2]));
     843          max_elevation_edgevalue[k] = max(elevation_edge_values[3*k],
     844                                           max(elevation_edge_values[3*k+1],
     845                                               elevation_edge_values[3*k+2]));
     846           
     847          }
     848
     849      }
     850
     851  // Compute some useful statistics on wetness/dryness
     852  for(k=0; k<number_of_elements;k++){
     853      count_wet_neighbours[k]=0;
     854      neighbour_wetness_scale[k] = 0.;
     855      // Cell k is wet
     856      if(stage_centroid_values[k] > max_elevation_edgevalue[k] + minimum_allowed_height+epsilon){
     857          neighbour_wetness_scale[k] = 1.; 
     858      }
     859      // Count the wet neighbours
     860      for (i=0; i<3; i++){
     861        ktmp = surrogate_neighbours[3*k+i];             
     862        if(stage_centroid_values[ktmp] > max_elevation_edgevalue[ktmp]){
     863            count_wet_neighbours[k]+=1;
     864        }
     865       
    861866      }
    862867  }
    863868
    864   // Begin extrapolation routine
    865   for (k = 0; k < number_of_elements; k++)
    866   {
    867     k3=k*3;
    868     k6=k*6;
    869 
    870     if (number_of_boundaries[k]==3)
    871     //if (0==0)
    872     {
    873       // No neighbours, set gradient on the triangle to zero
    874      
    875       stage_edge_values[k3]   = stage_centroid_values[k];
    876       stage_edge_values[k3+1] = stage_centroid_values[k];
    877       stage_edge_values[k3+2] = stage_centroid_values[k];
    878       xmom_edge_values[k3]    = xmom_centroid_values[k];
    879       xmom_edge_values[k3+1]  = xmom_centroid_values[k];
    880       xmom_edge_values[k3+2]  = xmom_centroid_values[k];
    881       ymom_edge_values[k3]    = ymom_centroid_values[k];
    882       ymom_edge_values[k3+1]  = ymom_centroid_values[k];
    883       ymom_edge_values[k3+2]  = ymom_centroid_values[k];
    884      
    885       continue;
    886     }
    887     else
    888     {
    889       // Triangle k has one or more neighbours.
    890       // Get centroid and edge coordinates of the triangle
    891      
    892       // Get the edge coordinates
    893       xv0 = edge_coordinates[k6];   
    894       yv0 = edge_coordinates[k6+1];
    895       xv1 = edge_coordinates[k6+2];
    896       yv1 = edge_coordinates[k6+3];
    897       xv2 = edge_coordinates[k6+4];
    898       yv2 = edge_coordinates[k6+5];
    899      
    900       // Get the centroid coordinates
    901       coord_index = 2*k;
    902       x = centroid_coordinates[coord_index];
    903       y = centroid_coordinates[coord_index+1];
    904      
    905       // Store x- and y- differentials for the edges of
    906       // triangle k relative to the centroid
    907       dxv0 = xv0 - x;
    908       dxv1 = xv1 - x;
    909       dxv2 = xv2 - x;
    910       dyv0 = yv0 - y;
    911       dyv1 = yv1 - y;
    912       dyv2 = yv2 - y;
    913       // Compute the minimum distance from the centroid to an edge
    914       //demin=min(dxv0*dxv0 +dyv0*dyv0, min(dxv1*dxv1+dyv1*dyv1, dxv2*dxv2+dyv2*dyv2));
    915       //demin=sqrt(demin);
    916     }
    917 
    918 
    919            
    920     if (number_of_boundaries[k]<=1)
    921     {
    922       //==============================================
    923       // Number of boundaries <= 1
    924       //==============================================   
    925    
    926    
    927       // If no boundaries, auxiliary triangle is formed
    928       // from the centroids of the three neighbours
    929       // If one boundary, auxiliary triangle is formed
    930       // from this centroid and its two neighbours
    931      
    932       k0 = surrogate_neighbours[k3];
    933       k1 = surrogate_neighbours[k3 + 1];
    934       k2 = surrogate_neighbours[k3 + 2];
    935      
    936       // Test to see whether we accept the surrogate neighbours
    937       // Note that if ki is replaced with k in more than 1 neighbour, then the
    938       // triangle area will be zero, and a first order extrapolation will be
    939       // used
    940       if(stage_centroid_values[k2]<=elevation_centroid_values[k2]){
    941           k2 = k ;
     869 
     870  // First treat all 'fully wet' cells, then treat 'partially dry' cells
     871  for(k_wetdry=0; k_wetdry<2; k_wetdry++){
     872      //if(((stage_centroid_values[k] > max_elevation_edgevalue[k]))==k_wetdry){
     873      if(((stage_centroid_values[k] < elevation_centroid_values[k]+minimum_allowed_height))==k_wetdry){
     874        continue;
    942875      }
    943       if(stage_centroid_values[k0]<=elevation_centroid_values[k0]){
    944           k0 = k ;
    945       }
    946       if(stage_centroid_values[k1]<=elevation_centroid_values[k1]){
    947           k1 = k ;
    948       }
    949       // Alternative approach (BRUTAL) -- if the max neighbour bed elevation is greater
    950       // than the min neighbour stage, then we use first order extrapolation
    951       //bedmax = max(elevation_centroid_values[k],
    952       //             max(elevation_centroid_values[k0],
    953       //                 max(elevation_centroid_values[k1], elevation_centroid_values[k2])));
    954       //stagemin = min(stage_centroid_values[k],
    955       //             min(stage_centroid_values[k0],
    956       //                 min(stage_centroid_values[k1], stage_centroid_values[k2])));
    957       //
    958       //if(stagemin < bedmax){
    959       //   // This will cause first order extrapolation
    960       //   k2 = k;
    961       //   k0 = k;
    962       //   k1 = k;
    963       //}
    964      
    965       // Get the auxiliary triangle's vertex coordinates
    966       // (really the centroids of neighbouring triangles)
    967       coord_index = 2*k0;
    968       x0 = centroid_coordinates[coord_index];
    969       y0 = centroid_coordinates[coord_index+1];
    970      
    971       coord_index = 2*k1;
    972       x1 = centroid_coordinates[coord_index];
    973       y1 = centroid_coordinates[coord_index+1];
    974      
    975       coord_index = 2*k2;
    976       x2 = centroid_coordinates[coord_index];
    977       y2 = centroid_coordinates[coord_index+1];
    978 
    979       // compute the maximum distance from the centroid to a neighbouring
    980       // centroid
    981       //dcmax=max( (x0-x)*(x0-x) + (y0-y)*(y0-y),     
    982       //           max((x1-x)*(x1-x) + (y1-y)*(y1-y),
    983       //               (x2-x)*(x2-x) + (y2-y)*(y2-y)));
    984       //dcmax=sqrt(dcmax);
    985       //// Ratio of centroid to edge distance -- useful in attempting to adapt limiter
    986       //if(dcmax>0.){
    987       //    r0scale=demin/dcmax;
    988       //    //printf("%f \n", r0scale);
    989       //}else{
    990       //    r0scale=0.5;
    991       //}
    992 
    993       // Store x- and y- differentials for the vertices
    994       // of the auxiliary triangle
    995       dx1 = x1 - x0;
    996       dx2 = x2 - x0;
    997       dy1 = y1 - y0;
    998       dy2 = y2 - y0;
    999      
    1000       // Calculate 2*area of the auxiliary triangle
    1001       // The triangle is guaranteed to be counter-clockwise     
    1002       area2 = dy2*dx1 - dy1*dx2;
    1003      
    1004       // If the mesh is 'weird' near the boundary,
    1005       // the triangle might be flat or clockwise
    1006       // Default to zero gradient
    1007       if (area2 <= 0)
     876      // For partially wet cells, we now know that the edges of neighbouring
     877      // fully wet cells have been defined
     878
     879      // Begin extrapolation routine
     880      for (k = 0; k < number_of_elements; k++)
    1008881      {
    1009           //printf("Error negative triangle area \n");
    1010           //report_python_error(AT, "Negative triangle area");
    1011           //return -1;
    1012 
     882        k3=k*3;
     883        k6=k*6;
     884
     885        if (number_of_boundaries[k]==3)
     886        {
     887          // No neighbours, set gradient on the triangle to zero
     888         
    1013889          stage_edge_values[k3]   = stage_centroid_values[k];
    1014890          stage_edge_values[k3+1] = stage_centroid_values[k];
    … …  
    1020896          ymom_edge_values[k3+1]  = ymom_centroid_values[k];
    1021897          ymom_edge_values[k3+2]  = ymom_centroid_values[k];
    1022 
     898         
    1023899          continue;
    1024       } 
    1025      
    1026       // Calculate heights of neighbouring cells
    1027       hc = stage_centroid_values[k]  - elevation_centroid_values[k];
    1028       h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
    1029       h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
    1030       h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
    1031       hmin = min(min(h0, min(h1, h2)), hc);
    1032       //hmin = min(h0, min(h1, h2));
    1033       //hmin = max(hmin, 0.0);
    1034       //hfactor = hc/(hc + 1.0);
    1035 
    1036       hfactor = 0.0;
    1037       //if (hmin > 0.001)
    1038       if (hmin > 0.)
    1039       //if (hc>0.0)
    1040       {
    1041         hfactor = 1.0 ;//hmin/(hmin + 0.004);
    1042         //hfactor=hmin/(hmin + 0.004);
    1043       }
    1044      
    1045       if (optimise_dry_cells)
    1046       {     
    1047         // Check if linear reconstruction is necessary for triangle k
    1048         // This check will exclude dry cells.
    1049 
    1050         //hmax = max(h0, max(h1, max(h2, hc)));     
    1051         hmax = max(h0, max(h1, h2));     
    1052         if (hmax < epsilon)
     900        }
     901        else
    1053902        {
    1054             continue;
     903          // Triangle k has one or more neighbours.
     904          // Get centroid and edge coordinates of the triangle
     905         
     906          // Get the edge coordinates
     907          xv0 = edge_coordinates[k6];   
     908          yv0 = edge_coordinates[k6+1];
     909          xv1 = edge_coordinates[k6+2];
     910          yv1 = edge_coordinates[k6+3];
     911          xv2 = edge_coordinates[k6+4];
     912          yv2 = edge_coordinates[k6+5];
     913         
     914          // Get the centroid coordinates
     915          coord_index = 2*k;
     916          x = centroid_coordinates[coord_index];
     917          y = centroid_coordinates[coord_index+1];
     918         
     919          // Store x- and y- differentials for the edges of
     920          // triangle k relative to the centroid
     921          dxv0 = xv0 - x;
     922          dxv1 = xv1 - x;
     923          dxv2 = xv2 - x;
     924          dyv0 = yv0 - y;
     925          dyv1 = yv1 - y;
     926          dyv2 = yv2 - y;
    1055927        }
    1056       }
    1057    
    1058       //-----------------------------------
    1059       // stage
    1060       //-----------------------------------     
    1061      
    1062       // Calculate the difference between vertex 0 of the auxiliary
    1063       // triangle and the centroid of triangle k
    1064       dq0 = stage_centroid_values[k0] - stage_centroid_values[k];
    1065      
    1066       // Calculate differentials between the vertices
    1067       // of the auxiliary triangle (centroids of neighbouring triangles)
    1068       dq1 = stage_centroid_values[k1] - stage_centroid_values[k0];
    1069       dq2 = stage_centroid_values[k2] - stage_centroid_values[k0];
    1070      
    1071       inv_area2 = 1.0/area2;
    1072       // Calculate the gradient of stage on the auxiliary triangle
    1073       a = dy2*dq1 - dy1*dq2;
    1074       a *= inv_area2;
    1075       b = dx1*dq2 - dx2*dq1;
    1076       b *= inv_area2;
    1077      
    1078       // Calculate provisional jumps in stage from the centroid
    1079       // of triangle k to its vertices, to be limited
    1080       dqv[0] = a*dxv0 + b*dyv0;
    1081       dqv[1] = a*dxv1 + b*dyv1;
    1082       dqv[2] = a*dxv2 + b*dyv2;
    1083    
    1084       // Idea: New limiter, computed using only neighbouring centroid values and local centroid value
    1085       // Step1: Compute the value of stage at the centroid of triangle k ,
    1086       //       using only the neighbouring centroid stage values
    1087       //tmp = a*(x-x0) + b*(y-y0) + stage_centroid_values[k0];
    1088       // Step2: Now, compute a limiting factor, so that if gradients are multiplied by this,
    1089       //        then for a triangle with these limited gradients, based at the local centroid value
    1090       //        ,the re-constructed neighbour centroid values will not over shoot (if larger) or undershoot (if smaller)
    1091       //limiting_factor=1 - max_all_k0((tmp - stage_centroid_values[k])/(tmp - stage_centroid_values[k0]))
    1092       //limiting_factor=min(max(limiting_factor,0),1)
    1093       // Advantages: No limiting for linear problem. Worth a look anyway
    1094  
    1095       // Now we want to find min and max of the centroid and the
    1096       // vertices of the auxiliary triangle and compute jumps
    1097       // from the centroid to the min and max
    1098       find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
    1099      
    1100       beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
    1101    
    1102       // Limit the gradient
    1103       //if(hmin/hc<0.5){
    1104       limit_gradient(dqv, qmin, qmax, beta_tmp);
    1105       //}
    1106       //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale);
    1107      
    1108       //for (i=0;i<3;i++)
    1109       stage_edge_values[k3+0] = stage_centroid_values[k] + dqv[0];
    1110       stage_edge_values[k3+1] = stage_centroid_values[k] + dqv[1];
    1111       stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2];
    1112      
    1113       //-----------------------------------
    1114       // xmomentum
    1115       //-----------------------------------           
    1116 
    1117       // Calculate the difference between vertex 0 of the auxiliary
    1118       // triangle and the centroid of triangle k     
    1119       dq0 = xmom_centroid_values[k0] - xmom_centroid_values[k];
    1120      
    1121       // Calculate differentials between the vertices
    1122       // of the auxiliary triangle
    1123       dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k0];
    1124       dq2 = xmom_centroid_values[k2] - xmom_centroid_values[k0];
    1125      
    1126       // Calculate the gradient of xmom on the auxiliary triangle
    1127       a = dy2*dq1 - dy1*dq2;
    1128       a *= inv_area2;
    1129       b = dx1*dq2 - dx2*dq1;
    1130       b *= inv_area2;
    1131      
    1132       // Calculate provisional jumps in stage from the centroid
    1133       // of triangle k to its vertices, to be limited     
    1134       dqv[0] = a*dxv0+b*dyv0;
    1135       dqv[1] = a*dxv1+b*dyv1;
    1136       dqv[2] = a*dxv2+b*dyv2;
    1137      
    1138       // Now we want to find min and max of the centroid and the
    1139       // vertices of the auxiliary triangle and compute jumps
    1140       // from the centroid to the min and max
    1141       //
    1142       find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
    1143 
    1144       beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor;
    1145 
    1146       // Limit the gradient
    1147       limit_gradient(dqv, qmin, qmax, beta_tmp);
    1148       //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale);
    1149      
    1150 
    1151       for (i=0; i < 3; i++)
    1152       {
    1153         xmom_edge_values[k3+i] = xmom_centroid_values[k] + dqv[i];
    1154       }
    1155      
    1156       //-----------------------------------
    1157       // ymomentum
    1158       //-----------------------------------                 
    1159 
    1160       // Calculate the difference between vertex 0 of the auxiliary
    1161       // triangle and the centroid of triangle k
    1162       dq0 = ymom_centroid_values[k0] - ymom_centroid_values[k];
    1163      
    1164       // Calculate differentials between the vertices
    1165       // of the auxiliary triangle
    1166       dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k0];
    1167       dq2 = ymom_centroid_values[k2] - ymom_centroid_values[k0];
    1168      
    1169       // Calculate the gradient of xmom on the auxiliary triangle
    1170       a = dy2*dq1 - dy1*dq2;
    1171       a *= inv_area2;
    1172       b = dx1*dq2 - dx2*dq1;
    1173       b *= inv_area2;
    1174      
    1175       // Calculate provisional jumps in stage from the centroid
    1176       // of triangle k to its vertices, to be limited
    1177       dqv[0] = a*dxv0 + b*dyv0;
    1178       dqv[1] = a*dxv1 + b*dyv1;
    1179       dqv[2] = a*dxv2 + b*dyv2;
    1180      
    1181       // Now we want to find min and max of the centroid and the
    1182       // vertices of the auxiliary triangle and compute jumps
    1183       // from the centroid to the min and max
    1184       //
    1185       find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
    1186      
    1187       beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;   
    1188 
    1189       // Limit the gradient
    1190       limit_gradient(dqv, qmin, qmax, beta_tmp);
    1191       //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale);
    1192      
    1193       for (i=0;i<3;i++)
    1194       {
    1195         ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
    1196       }
    1197      
    1198     } // End number_of_boundaries <=1
    1199     else
    1200     {
    1201 
    1202       //==============================================
    1203       // Number of boundaries == 2
    1204       //==============================================       
     928
     929
     930               
     931        if (number_of_boundaries[k]<=1)
     932        {
     933          //==============================================
     934          // Number of boundaries <= 1
     935          //==============================================   
    1205936       
    1206       // One internal neighbour and gradient is in direction of the neighbour's centroid
    1207      
    1208       // Find the only internal neighbour (k1?)
    1209       for (k2 = k3; k2 < k3 + 3; k2++)
    1210       {
    1211       // Find internal neighbour of triangle k     
    1212       // k2 indexes the edges of triangle k   
    1213      
    1214           if (surrogate_neighbours[k2] != k)
     937       
     938          // If no boundaries, auxiliary triangle is formed
     939          // from the centroids of the three neighbours
     940          // If one boundary, auxiliary triangle is formed
     941          // from this centroid and its two neighbours
     942         
     943          k0 = surrogate_neighbours[k3];
     944          k1 = surrogate_neighbours[k3 + 1];
     945          k2 = surrogate_neighbours[k3 + 2];
     946         
     947          // Test to see whether we accept the surrogate neighbours
     948          // Note that if ki is replaced with k in more than 1 neighbour, then the
     949          // triangle area will be zero, and a first order extrapolation will be
     950          // used
     951          //if(stage_centroid_values[k2]<=max_elevation_edgevalue[k2]+epsilon){
     952          if(stage_centroid_values[k2]<=elevation_centroid_values[k2]+minimum_allowed_height+epsilon){
     953              k2 = k ;
     954          }
     955          //if(stage_centroid_values[k0]<=max_elevation_edgevalue[k0]+epsilon){
     956          if(stage_centroid_values[k0]<=elevation_centroid_values[k0]+minimum_allowed_height+epsilon){
     957              k0 = k ;
     958          }
     959          //if(stage_centroid_values[k1]<=max_elevation_edgevalue[k1]+epsilon){
     960          if(stage_centroid_values[k1]<=elevation_centroid_values[k1]+minimum_allowed_height+epsilon){
     961              k1 = k ;
     962          }
     963
     964          // Take note if the max neighbour bed elevation is greater than the min
     965          // neighbour stage -- suggests a 'steep' bed relative to the flow
     966          //bedmax = max(elevation_centroid_values[k],
     967          //             max(elevation_centroid_values[k0],
     968          //                 max(elevation_centroid_values[k1],
     969          //                     elevation_centroid_values[k2])));
     970          //stagemin = min(max(stage_centroid_values[k], elevation_centroid_values[k]),
     971          //               min(max(stage_centroid_values[k0], elevation_centroid_values[k0]),
     972          //                   min(max(stage_centroid_values[k1], elevation_centroid_values[k1]),
     973          //                       max(stage_centroid_values[k2], elevation_centroid_values[k2]))));
     974
     975          // Get the auxiliary triangle's vertex coordinates
     976          // (really the centroids of neighbouring triangles)
     977          coord_index = 2*k0;
     978          x0 = centroid_coordinates[coord_index];
     979          y0 = centroid_coordinates[coord_index+1];
     980         
     981          coord_index = 2*k1;
     982          x1 = centroid_coordinates[coord_index];
     983          y1 = centroid_coordinates[coord_index+1];
     984         
     985          coord_index = 2*k2;
     986          x2 = centroid_coordinates[coord_index];
     987          y2 = centroid_coordinates[coord_index+1];
     988
     989          // Store x- and y- differentials for the vertices
     990          // of the auxiliary triangle
     991          dx1 = x1 - x0;
     992          dx2 = x2 - x0;
     993          dy1 = y1 - y0;
     994          dy2 = y2 - y0;
     995         
     996          // Calculate 2*area of the auxiliary triangle
     997          // The triangle is guaranteed to be counter-clockwise     
     998          area2 = dy2*dx1 - dy1*dx2;
     999           
     1000          // Treat triangles with zero or 1 wet neighbours.
     1001          if ((area2 <= 0)) //|(count_wet_neighbours[k]==0))
    12151002          {
    1216              break;
     1003              if( ((k==k0) & (k==k1) & (k==k2))){
     1004                  // Isolated wet cell  -- use constant depth extrapolation for stage
     1005                  //stage_edge_values[k3]   = stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3];
     1006                  //stage_edge_values[k3+1] = stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+1];
     1007                  //stage_edge_values[k3+2] = stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+2];
     1008
     1009                  // Isolated wet cell -- constant depth extrapolation, but don't let water flow uphill
     1010                  //stage_edge_values[k3]   = min(stage_centroid_values[k], stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3]);
     1011                  //stage_edge_values[k3+1] = min(stage_centroid_values[k], stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+1]);
     1012                  //stage_edge_values[k3+2] = min(stage_centroid_values[k], stage_centroid_values[k]- elevation_centroid_values[k] + elevation_edge_values[k3+2]);
     1013                  // Isolated wet cell -- constant stage extrapolation
     1014                  stage_edge_values[k3]   = stage_centroid_values[k];
     1015                  stage_edge_values[k3+1] = stage_centroid_values[k];
     1016                  stage_edge_values[k3+2] = stage_centroid_values[k];
     1017              }else{
     1018                  // Cell with one wet neighbour.
     1019                  // Find which neighbour is wet [if k!=k0, then k0 is wet]
     1020                  if(k!=k0){
     1021                    ktmp=k0;
     1022                    ii=0;
     1023                    xtmp = x0;
     1024                    ytmp = y0;
     1025                  }else if(k!=k1){
     1026                    ktmp=k1;
     1027                    ii=1;
     1028                    xtmp = x1;
     1029                    ytmp = y1;
     1030                  }else if(k!=k2){
     1031                    ktmp=k2;
     1032                    ii=2;
     1033                    xtmp = x2;
     1034                    ytmp = y2;
     1035                  }else{
     1036                      printf("ERROR in extrapolation routine: Should have had one ki == k \n");
     1037                      report_python_error(AT, "Negative triangle area");
     1038                      return -1;
     1039                  }
     1040
     1041                  //if(stage_centroid_values[k]>elevation_centroid_values[k]+minimum_allowed_height){
     1042                  if(k_wetdry==0){
     1043                  //    // Cell is wet
     1044                  //    //if(count_wet_neighbours[ktmp]>0){
     1045                  //    //    stage_edge_values[k3]= stage_centroid_values[k];
     1046                  //    //    stage_edge_values[k3+1]= stage_centroid_values[k];
     1047                  //    //    stage_edge_values[k3+2]= stage_centroid_values[k];
     1048                  //    //}else{
     1049                  //        // Constant depth extrapolation
     1050                      //if(stage_centroid_values[k]<stage_centroid_values[ktmp]){
     1051                          //stage_edge_values[k3]= stage_centroid_values[k]-elevation_centroid_values[k]+elevation_edge_values[k3];
     1052                          //stage_edge_values[k3+1]= stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+1];
     1053                          //stage_edge_values[k3+2]= stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+2];
     1054                      stage_edge_values[k3] = stage_centroid_values[k];
     1055                      stage_edge_values[k3+1] = stage_centroid_values[k];
     1056                      stage_edge_values[k3+2] = stage_centroid_values[k];
     1057                      //}else{
     1058                      //    stage_edge_values[k3] = stage_centroid_values[ktmp];
     1059                      //    stage_edge_values[k3+1] = stage_centroid_values[ktmp];
     1060                      //    stage_edge_values[k3+2] = stage_centroid_values[ktmp];
     1061                      //    //stage_edge_values[k3]= max(stage_centroid_values[k]-elevation_centroid_values[k]+elevation_edge_values[k3], stage_centroid_values[ktmp]);
     1062                      //    //stage_edge_values[k3+1]= max(stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+1],stage_centroid_values[ktmp]);
     1063                      //    //stage_edge_values[k3+2]= max(stage_centroid_values[k] -elevation_centroid_values[k]+elevation_edge_values[k3+2], stage_centroid_values[ktmp]);
     1064                      //}
     1065
     1066                  }else{
     1067                     
     1068                  //    // This cell is 'dry'.
     1069                  //    // Extrapolate using the neighbouring wet cell if appropriate
     1070                  //    // This needs to preserve a wet-dry interface
     1071                      stage_edge_values[k3]= stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+ii]];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+i]], elevation_edge_values[k3]) ;
     1072                      stage_edge_values[k3+1]= stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+ii]];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+i]], elevation_edge_values[k3+1]);
     1073                      stage_edge_values[k3+2]= stage_edge_values[3*ktmp+neighbour_edges[3*ktmp+ii]];//max(stage_edge_values[3*ktmp+neighbour_edges[3*ktmp +i]], elevation_edge_values[k3+2]);
     1074                  //    stage_edge_values[k3]= stage_centroid_values[ktmp] ;
     1075                  //    stage_edge_values[k3+1]= stage_centroid_values[ktmp];
     1076                  //    stage_edge_values[k3+2]= stage_centroid_values[ktmp];
     1077                  //    stage_edge_values[k3]= stage_centroid_values[k] ;
     1078                  //    stage_edge_values[k3+1]= stage_centroid_values[k];
     1079                  //    stage_edge_values[k3+2]= stage_centroid_values[k];
     1080                  }
     1081
     1082              }
     1083              // First order momentum / velocity extrapolation
     1084              //stage_edge_values[k3]    = stage_centroid_values[k];
     1085              //stage_edge_values[k3+1]  = stage_centroid_values[k];
     1086              //stage_edge_values[k3+2]  = stage_centroid_values[k];
     1087              xmom_edge_values[k3]    = xmom_centroid_values[k];
     1088              xmom_edge_values[k3+1]  = xmom_centroid_values[k];
     1089              xmom_edge_values[k3+2]  = xmom_centroid_values[k];
     1090              ymom_edge_values[k3]    = ymom_centroid_values[k];
     1091              ymom_edge_values[k3+1]  = ymom_centroid_values[k];
     1092              ymom_edge_values[k3+2]  = ymom_centroid_values[k];
     1093
     1094              continue;
     1095          } 
     1096         
     1097          // Calculate heights of neighbouring cells
     1098          hc = stage_centroid_values[k]  - elevation_centroid_values[k];
     1099          h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
     1100          h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
     1101          h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
     1102          hmin = min(min(h0, min(h1, h2)), hc);
     1103
     1104          hfactor = 0.0;
     1105          //if (hmin > 0.001)
     1106          if (hmin > 0.)
     1107          //if (hc>0.0)
     1108          {
     1109            hfactor = 1.0 ;//hmin/(hmin + 0.004);
     1110            //hfactor=hmin/(hmin + 0.004);
    12171111          }
    1218       }
    1219      
    1220       if ((k2 == k3 + 3))
    1221       {
    1222         // If we didn't find an internal neighbour
    1223         report_python_error(AT, "Internal neighbour not found");
    1224         return -1;
    1225       }
    1226      
    1227       k1 = surrogate_neighbours[k2];
    1228      
    1229       // The coordinates of the triangle are already (x,y).
    1230       // Get centroid of the neighbour (x1,y1)
    1231       coord_index = 2*k1;
    1232       x1 = centroid_coordinates[coord_index];
    1233       y1 = centroid_coordinates[coord_index + 1];
    1234      
    1235       // Compute x- and y- distances between the centroid of
    1236       // triangle k and that of its neighbour
    1237       dx1 = x1 - x;
    1238       dy1 = y1 - y;
    1239      
    1240       // Set area2 as the square of the distance
    1241       area2 = dx1*dx1 + dy1*dy1;
    1242      
    1243       // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2)
    1244       // and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
    1245       // respectively correspond to the x- and y- gradients
    1246       // of the conserved quantities
    1247       dx2 = 1.0/area2;
    1248       dy2 = dx2*dy1;
    1249       dx2 *= dx1;
    1250      
    1251      
    1252       //-----------------------------------
    1253       // stage
    1254       //-----------------------------------           
    1255 
    1256       // Compute differentials
    1257       dq1 = stage_centroid_values[k1] - stage_centroid_values[k];
    1258      
    1259       // Calculate the gradient between the centroid of triangle k
    1260       // and that of its neighbour
    1261       a = dq1*dx2;
    1262       b = dq1*dy2;
    1263      
    1264       // Calculate provisional edge jumps, to be limited
    1265       dqv[0] = a*dxv0 + b*dyv0;
    1266       dqv[1] = a*dxv1 + b*dyv1;
    1267       dqv[2] = a*dxv2 + b*dyv2;
    1268      
    1269       // Now limit the jumps
    1270       if (dq1>=0.0)
    1271       {
    1272         qmin=0.0;
    1273         qmax=dq1;
    1274       }
    1275       else
    1276       {
    1277         qmin = dq1;
    1278         qmax = 0.0;
    1279       }
    1280      
    1281       // Limit the gradient
    1282       limit_gradient(dqv, qmin, qmax, beta_w);
    1283      
    1284       //for (i=0; i < 3; i++)
    1285       //{
    1286       stage_edge_values[k3] = stage_centroid_values[k] + dqv[0];
    1287       stage_edge_values[k3 + 1] = stage_centroid_values[k] + dqv[1];
    1288       stage_edge_values[k3 + 2] = stage_centroid_values[k] + dqv[2];
    1289       //}
    1290      
    1291       //-----------------------------------
    1292       // xmomentum
    1293       //-----------------------------------                       
    1294      
    1295       // Compute differentials
    1296       dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k];
    1297      
    1298       // Calculate the gradient between the centroid of triangle k
    1299       // and that of its neighbour
    1300       a = dq1*dx2;
    1301       b = dq1*dy2;
    1302      
    1303       // Calculate provisional edge jumps, to be limited
    1304       dqv[0] = a*dxv0+b*dyv0;
    1305       dqv[1] = a*dxv1+b*dyv1;
    1306       dqv[2] = a*dxv2+b*dyv2;
    1307      
    1308       // Now limit the jumps
    1309       if (dq1 >= 0.0)
    1310       {
    1311         qmin = 0.0;
    1312         qmax = dq1;
    1313       }
    1314       else
    1315       {
    1316         qmin = dq1;
    1317         qmax = 0.0;
    1318       }
    1319      
    1320       // Limit the gradient     
    1321       limit_gradient(dqv, qmin, qmax, beta_w);
    1322      
    1323       //for (i=0;i<3;i++)
    1324       //xmom_edge_values[k3] = xmom_centroid_values[k] + dqv[0];
    1325       //xmom_edge_values[k3 + 1] = xmom_centroid_values[k] + dqv[1];
    1326       //xmom_edge_values[k3 + 2] = xmom_centroid_values[k] + dqv[2];
    1327      
    1328       for (i = 0; i < 3;i++)
    1329       {
    1330           xmom_edge_values[k3 + i] = xmom_centroid_values[k] + dqv[i];
    1331       }
    1332      
    1333       //-----------------------------------
    1334       // ymomentum
    1335       //-----------------------------------                       
    1336 
    1337       // Compute differentials
    1338       dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k];
    1339      
    1340       // Calculate the gradient between the centroid of triangle k
    1341       // and that of its neighbour
    1342       a = dq1*dx2;
    1343       b = dq1*dy2;
    1344      
    1345       // Calculate provisional edge jumps, to be limited
    1346       dqv[0] = a*dxv0 + b*dyv0;
    1347       dqv[1] = a*dxv1 + b*dyv1;
    1348       dqv[2] = a*dxv2 + b*dyv2;
    1349      
    1350       // Now limit the jumps
    1351       if (dq1>=0.0)
    1352       {
    1353         qmin = 0.0;
    1354         qmax = dq1;
    1355       }
    1356       else
    1357       {
    1358         qmin = dq1;
    1359         qmax = 0.0;
    1360       }
    1361      
    1362       // Limit the gradient
    1363       limit_gradient(dqv, qmin, qmax, beta_w);
    1364      
    1365       for (i=0;i<3;i++)
     1112         
     1113          //-----------------------------------
     1114          // stage
     1115          //-----------------------------------     
     1116         
     1117          // Calculate the difference between vertex 0 of the auxiliary
     1118          // triangle and the centroid of triangle k
     1119          dq0 = stage_centroid_values[k0] - stage_centroid_values[k];
     1120         
     1121          // Calculate differentials between the vertices
     1122          // of the auxiliary triangle (centroids of neighbouring triangles)
     1123          dq1 = stage_centroid_values[k1] - stage_centroid_values[k0];
     1124          dq2 = stage_centroid_values[k2] - stage_centroid_values[k0];
     1125         
     1126          inv_area2 = 1.0/area2;
     1127          // Calculate the gradient of stage on the auxiliary triangle
     1128          a = dy2*dq1 - dy1*dq2;
     1129          a *= inv_area2;
     1130          b = dx1*dq2 - dx2*dq1;
     1131          b *= inv_area2;
     1132         
     1133          // Calculate provisional jumps in stage from the centroid
     1134          // of triangle k to its vertices, to be limited
     1135          dqv[0] = a*dxv0 + b*dyv0;
     1136          dqv[1] = a*dxv1 + b*dyv1;
     1137          dqv[2] = a*dxv2 + b*dyv2;
     1138       
     1139          // Now we want to find min and max of the centroid and the
     1140          // vertices of the auxiliary triangle and compute jumps
     1141          // from the centroid to the min and max
     1142          find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
     1143         
     1144          beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
     1145       
     1146         
     1147          // Limit the gradient
     1148          // This is more complex for stage than other quantities
     1149          if((count_wet_neighbours[k]>0)){ //&&(stage_centroid_values[k] > elevation_centroid_values[k])){
     1150              limit_gradient(dqv, qmin, qmax, beta_tmp);
     1151              stage_edge_values[k3+0] = stage_centroid_values[k] + dqv[0];
     1152              stage_edge_values[k3+1] = stage_centroid_values[k] + dqv[1];
     1153              stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2];
     1154          }else{
     1155          //if(count_wet_neighbours[k]==0){
     1156              //weight=max(neighbour_wetness_scale[k] - (stage_centroid_values[k]-elevation_centroid_values[k]), 0.);
     1157              //weight/=(max_elevation_edgevalue[k] - elevation_centroid_values[k]);
     1158              //weight=min(weight, 1.0);
     1159              //weight==0;
     1160              //weight=0.;
     1161              // 'Shallow flow on steep slope'
     1162              stage_edge_values[k3+0] =stage_centroid_values[k] + dqv[0];
     1163              stage_edge_values[k3+1] =stage_centroid_values[k] + dqv[1];
     1164              stage_edge_values[k3+2] =stage_centroid_values[k] + dqv[2];
     1165
     1166              // There exists a neighbouring bed cell with elevation > minimum
     1167              // neighbouring stage value
     1168
     1169              // We have to be careful in this situation. Limiting the stage gradient can
     1170              // cause problems for shallow flows down steep slopes (rainfall-runoff type problems)
     1171              // Previously 'balance_deep_and_shallow' was used to deal with this issue
     1172              //for(i=0; i<3; i++){
     1173              //    stage_edge_values[k3+i] = weight*stage_edge_values[k3+i] +
     1174              //                             (1.0-weight)*(stage_centroid_values[k] -elevation_centroid_values[k]
     1175              //                                           +elevation_edge_values[k3+i]);
     1176              //}
     1177          }
     1178           
     1179          //-----------------------------------
     1180          // xmomentum
     1181          //-----------------------------------           
     1182
     1183          // Calculate the difference between vertex 0 of the auxiliary
     1184          // triangle and the centroid of triangle k     
     1185          dq0 = xmom_centroid_values[k0] - xmom_centroid_values[k];
     1186         
     1187          // Calculate differentials between the vertices
     1188          // of the auxiliary triangle
     1189          dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k0];
     1190          dq2 = xmom_centroid_values[k2] - xmom_centroid_values[k0];
     1191         
     1192          // Calculate the gradient of xmom on the auxiliary triangle
     1193          a = dy2*dq1 - dy1*dq2;
     1194          a *= inv_area2;
     1195          b = dx1*dq2 - dx2*dq1;
     1196          b *= inv_area2;
     1197         
     1198          // Calculate provisional jumps in stage from the centroid
     1199          // of triangle k to its vertices, to be limited     
     1200          dqv[0] = a*dxv0+b*dyv0;
     1201          dqv[1] = a*dxv1+b*dyv1;
     1202          dqv[2] = a*dxv2+b*dyv2;
     1203         
     1204          // Now we want to find min and max of the centroid and the
     1205          // vertices of the auxiliary triangle and compute jumps
     1206          // from the centroid to the min and max
     1207          //
     1208          find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
     1209
     1210          beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor;
     1211
     1212          // Limit the gradient
     1213          //if((k!=k0)&(k!=k1)&(k!=k2))
     1214          //if(stagemin>=bedmax){
     1215          if((count_wet_neighbours[k]>0) &&
     1216            (stage_centroid_values[k]>(minimum_allowed_height+elevation_centroid_values[k]))){
     1217            limit_gradient(dqv, qmin, qmax, beta_tmp);
     1218          }else{
     1219            dqv[0]=0.;
     1220            dqv[1]=0.;
     1221            dqv[2]=0.;
     1222          //  limit_gradient(dqv, qmin, qmax, 0.);
     1223          }
     1224          //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale);
     1225         
     1226
     1227          for (i=0; i < 3; i++)
     1228          {
     1229            xmom_edge_values[k3+i] = xmom_centroid_values[k] + dqv[i];
     1230          }
     1231         
     1232          //-----------------------------------
     1233          // ymomentum
     1234          //-----------------------------------                 
     1235
     1236          // Calculate the difference between vertex 0 of the auxiliary
     1237          // triangle and the centroid of triangle k
     1238          dq0 = ymom_centroid_values[k0] - ymom_centroid_values[k];
     1239         
     1240          // Calculate differentials between the vertices
     1241          // of the auxiliary triangle
     1242          dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k0];
     1243          dq2 = ymom_centroid_values[k2] - ymom_centroid_values[k0];
     1244         
     1245          // Calculate the gradient of xmom on the auxiliary triangle
     1246          a = dy2*dq1 - dy1*dq2;
     1247          a *= inv_area2;
     1248          b = dx1*dq2 - dx2*dq1;
     1249          b *= inv_area2;
     1250         
     1251          // Calculate provisional jumps in stage from the centroid
     1252          // of triangle k to its vertices, to be limited
     1253          dqv[0] = a*dxv0 + b*dyv0;
     1254          dqv[1] = a*dxv1 + b*dyv1;
     1255          dqv[2] = a*dxv2 + b*dyv2;
     1256         
     1257          // Now we want to find min and max of the centroid and the
     1258          // vertices of the auxiliary triangle and compute jumps
     1259          // from the centroid to the min and max
     1260          //
     1261          find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
     1262         
     1263          beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor;   
     1264
     1265          // Limit the gradient
     1266          //if(stagemin>=bedmax){
     1267          if((count_wet_neighbours[k]>0)&&
     1268            (stage_centroid_values[k]>(minimum_allowed_height+elevation_centroid_values[k]))){
     1269            limit_gradient(dqv, qmin, qmax, beta_tmp);
     1270          }else{
     1271            dqv[0]=0.;
     1272            dqv[1]=0.;
     1273            dqv[2]=0.;
     1274          }
     1275         
     1276          for (i=0;i<3;i++)
     1277          {
     1278            ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
     1279          }
     1280         
     1281        } // End number_of_boundaries <=1
     1282        else
     1283        {
     1284
     1285          //==============================================
     1286          // Number of boundaries == 2
     1287          //==============================================       
     1288           
     1289          // One internal neighbour and gradient is in direction of the neighbour's centroid
     1290         
     1291          // Find the only internal neighbour (k1?)
     1292          for (k2 = k3; k2 < k3 + 3; k2++)
     1293          {
     1294          // Find internal neighbour of triangle k     
     1295          // k2 indexes the edges of triangle k   
     1296         
     1297              if (surrogate_neighbours[k2] != k)
    13661298              {
    1367               ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
     1299                 break;
    13681300              }
    1369     } // else [number_of_boundaries==2]
    1370   } // for k=0 to number_of_elements-1
    1371  
     1301          }
     1302         
     1303          if ((k2 == k3 + 3))
     1304          {
     1305            // If we didn't find an internal neighbour
     1306            report_python_error(AT, "Internal neighbour not found");
     1307            return -1;
     1308          }
     1309         
     1310          k1 = surrogate_neighbours[k2];
     1311         
     1312          // The coordinates of the triangle are already (x,y).
     1313          // Get centroid of the neighbour (x1,y1)
     1314          coord_index = 2*k1;
     1315          x1 = centroid_coordinates[coord_index];
     1316          y1 = centroid_coordinates[coord_index + 1];
     1317         
     1318          // Compute x- and y- distances between the centroid of
     1319          // triangle k and that of its neighbour
     1320          dx1 = x1 - x;
     1321          dy1 = y1 - y;
     1322         
     1323          // Set area2 as the square of the distance
     1324          area2 = dx1*dx1 + dy1*dy1;
     1325         
     1326          // Set dx2=(x1-x0)/((x1-x0)^2+(y1-y0)^2)
     1327          // and dy2=(y1-y0)/((x1-x0)^2+(y1-y0)^2) which
     1328          // respectively correspond to the x- and y- gradients
     1329          // of the conserved quantities
     1330          dx2 = 1.0/area2;
     1331          dy2 = dx2*dy1;
     1332          dx2 *= dx1;
     1333         
     1334         
     1335          //-----------------------------------
     1336          // stage
     1337          //-----------------------------------           
     1338
     1339          // Compute differentials
     1340          dq1 = stage_centroid_values[k1] - stage_centroid_values[k];
     1341         
     1342          // Calculate the gradient between the centroid of triangle k
     1343          // and that of its neighbour
     1344          a = dq1*dx2;
     1345          b = dq1*dy2;
     1346         
     1347          // Calculate provisional edge jumps, to be limited
     1348          dqv[0] = a*dxv0 + b*dyv0;
     1349          dqv[1] = a*dxv1 + b*dyv1;
     1350          dqv[2] = a*dxv2 + b*dyv2;
     1351         
     1352          // Now limit the jumps
     1353          if (dq1>=0.0)
     1354          {
     1355            qmin=0.0;
     1356            qmax=dq1;
     1357          }
     1358          else
     1359          {
     1360            qmin = dq1;
     1361            qmax = 0.0;
     1362          }
     1363         
     1364          // Limit the gradient
     1365          limit_gradient(dqv, qmin, qmax, beta_w);
     1366         
     1367          //for (i=0; i < 3; i++)
     1368          //{
     1369          stage_edge_values[k3] = stage_centroid_values[k] + dqv[0];
     1370          stage_edge_values[k3 + 1] = stage_centroid_values[k] + dqv[1];
     1371          stage_edge_values[k3 + 2] = stage_centroid_values[k] + dqv[2];
     1372          //}
     1373         
     1374          //-----------------------------------
     1375          // xmomentum
     1376          //-----------------------------------                       
     1377         
     1378          // Compute differentials
     1379          dq1 = xmom_centroid_values[k1] - xmom_centroid_values[k];
     1380         
     1381          // Calculate the gradient between the centroid of triangle k
     1382          // and that of its neighbour
     1383          a = dq1*dx2;
     1384          b = dq1*dy2;
     1385         
     1386          // Calculate provisional edge jumps, to be limited
     1387          dqv[0] = a*dxv0+b*dyv0;
     1388          dqv[1] = a*dxv1+b*dyv1;
     1389          dqv[2] = a*dxv2+b*dyv2;
     1390         
     1391          // Now limit the jumps
     1392          if (dq1 >= 0.0)
     1393          {
     1394            qmin = 0.0;
     1395            qmax = dq1;
     1396          }
     1397          else
     1398          {
     1399            qmin = dq1;
     1400            qmax = 0.0;
     1401          }
     1402         
     1403          // Limit the gradient     
     1404          limit_gradient(dqv, qmin, qmax, beta_w);
     1405         
     1406          //for (i=0;i<3;i++)
     1407          //xmom_edge_values[k3] = xmom_centroid_values[k] + dqv[0];
     1408          //xmom_edge_values[k3 + 1] = xmom_centroid_values[k] + dqv[1];
     1409          //xmom_edge_values[k3 + 2] = xmom_centroid_values[k] + dqv[2];
     1410         
     1411          for (i = 0; i < 3;i++)
     1412          {
     1413              xmom_edge_values[k3 + i] = xmom_centroid_values[k] + dqv[i];
     1414          }
     1415         
     1416          //-----------------------------------
     1417          // ymomentum
     1418          //-----------------------------------                       
     1419
     1420          // Compute differentials
     1421          dq1 = ymom_centroid_values[k1] - ymom_centroid_values[k];
     1422         
     1423          // Calculate the gradient between the centroid of triangle k
     1424          // and that of its neighbour
     1425          a = dq1*dx2;
     1426          b = dq1*dy2;
     1427         
     1428          // Calculate provisional edge jumps, to be limited
     1429          dqv[0] = a*dxv0 + b*dyv0;
     1430          dqv[1] = a*dxv1 + b*dyv1;
     1431          dqv[2] = a*dxv2 + b*dyv2;
     1432         
     1433          // Now limit the jumps
     1434          if (dq1>=0.0)
     1435          {
     1436            qmin = 0.0;
     1437            qmax = dq1;
     1438          }
     1439          else
     1440          {
     1441            qmin = dq1;
     1442            qmax = 0.0;
     1443          }
     1444         
     1445          // Limit the gradient
     1446          limit_gradient(dqv, qmin, qmax, beta_w);
     1447         
     1448          for (i=0;i<3;i++)
     1449                  {
     1450                  ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i];
     1451                  }
     1452        } // else [number_of_boundaries==2]
     1453      } // for k=0 to number_of_elements-1
     1454  }
     1455
     1456  //for(k=0;k<number_of_elements; k++){
     1457  //    // Hack for 'dry' cells
     1458  //    // Make sure edge value is not greater than neighbour edge value
     1459  //    if(stage_centroid_values[k] - elevation_centroid_values[k] < minimum_allowed_height+epsilon){
     1460  //      for(i=0; i<3;i++){
     1461  //          if(surrogate_neighbours[3*k+i] >= 0){
     1462  //              ktmp=3*surrogate_neighbours[3*k+i] + neighbour_edges[3*k+i];
     1463  //          //if(neighbour_wetness_scale[surrogate_neighbours[3*k+i]]>0.){
     1464  //              stage_edge_values[3*k+i] = max(stage_edge_values[3*k+i], stage_edge_values[ktmp]);
     1465  //          }
     1466  //      }
     1467  //    }
     1468  //}
     1469
    13721470  // Compute vertex values of quantities
    13731471  for (k=0; k<number_of_elements; k++){
    … …  
    13981496          for (i=0; i<3; i++){
    13991497              de[i] = max(stage_edge_values[k3+i]-elevation_edge_values[k3+i],0.0);
     1498              //if(de[i]<minimum_allowed_height){
     1499              //  de[i]=0.;
     1500              //}
    14001501              xmom_edge_values[k3+i]=xmom_edge_values[k3+i]*de[i];
    14011502              ymom_edge_values[k3+i]=ymom_edge_values[k3+i]*de[i];
    … …  
    14171518  free(ymom_centroid_store);
    14181519  free(stage_centroid_store);
     1520  free(min_elevation_edgevalue);
     1521  free(max_elevation_edgevalue);
     1522  free(neighbour_wetness_scale);
     1523  free(count_wet_neighbours);
    14191524
    14201525  return 0;
    … …  
    17401845  */
    17411846  PyArrayObject *surrogate_neighbours,
     1847    *neighbour_edges,
    17421848    *number_of_boundaries,
    17431849    *centroid_coordinates,
    … …  
    17641870 
    17651871  // Convert Python arguments to C
    1766   if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOOOOii",
     1872  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOOOOOii",
    17671873            &domain,
    17681874            &surrogate_neighbours,
     1875            &neighbour_edges,
    17691876            &number_of_boundaries,
    17701877            &centroid_coordinates,
    … …  
    17911898  // check that numpy array objects arrays are C contiguous memory
    17921899  CHECK_C_CONTIG(surrogate_neighbours);
     1900  CHECK_C_CONTIG(neighbour_edges);
    17931901  CHECK_C_CONTIG(number_of_boundaries);
    17941902  CHECK_C_CONTIG(centroid_coordinates);
    … …  
    18361944                   beta_vh_dry,
    18371945                   (long*) surrogate_neighbours -> data,
     1946                   (long*) neighbour_edges -> data,
    18381947                   (long*) number_of_boundaries -> data,
    18391948                   (double*) centroid_coordinates -> data,
    … …  
    18791988  *zv,            //Elevation at vertices
    18801989  *xmomc,         //Momentums at centroids
    1881   *ymomc;
    1882 
     1990  *ymomc,
     1991  *areas,         // Triangle areas
     1992  *xc,
     1993  *yc;
    18831994
    18841995  int N;
    … …  
    18861997
    18871998  // Convert Python arguments to C
    1888   if (!PyArg_ParseTuple(args, "dddOOOOOO",
     1999  if (!PyArg_ParseTuple(args, "dddOOOOOOOOO",
    18892000            &minimum_allowed_height,
    18902001            &maximum_allowed_speed,
    18912002            &epsilon,
    1892             &wc, &wv, &zc, &zv, &xmomc, &ymomc)) {
     2003            &wc, &wv, &zc, &zv, &xmomc, &ymomc, &areas, &xc, &yc)) {
    18932004    report_python_error(AT, "could not parse input arguments");
    18942005    return NULL;
    … …  
    19062017       (double*) zv -> data,
    19072018       (double*) xmomc -> data,
    1908        (double*) ymomc -> data);
     2019       (double*) ymomc -> data,
     2020       (double*) areas -> data,
     2021       (double*) xc -> data,
     2022       (double*) yc -> data );
    19092023
    19102024  return Py_BuildValue("");

  • trunk/anuga_work/development/gareth/tests/channel_floodplain/channel_floodplain1.py

    r8446 r8547  
    1414#from anuga import *
    1515#from balanced_basic import *
    16 from balanced_dev import *
     16#from balanced_dev import *
     17from anuga_tsunami import *
    1718#from balanced_basic.swb2_domain import *
    1819#from balanced_basic.swb2_domain import Domain

  • trunk/anuga_work/development/gareth/tests/channel_floodplain/plotme.py

    r8446 r8547  
    44
    55# Time-index to plot outputs from
    6 index=60
     6index=1200
    77
    88#p2 = util.get_output('channel_floodplain1_bal_dev.sww', minimum_allowed_height=0.01)

  • trunk/anuga_work/development/gareth/tests/dam_break/dam_break.py

    r8353 r8547  
    1515from time import localtime, strftime, gmtime
    1616from balanced_dev import *
     17#from anuga_tsunami import *
    1718
    1819

  • trunk/anuga_work/development/gareth/tests/dam_break/plotme.py

    r8466 r8547  
    99p2_st=util.get_centroids(p_st)
    1010
    11 p_dev = util.get_output('dam_break_20120627_231618/dam_break.sww', 0.001)
     11p_dev = util.get_output('dam_break_20120831_172309/dam_break.sww', 0.001)
    1212p2_dev=util.get_centroids(p_dev, velocity_extrapolation=True)
    1313

  • trunk/anuga_work/development/gareth/tests/merimbula_steve/run_sw_merimbula.py

    r8396 r8547  
    2929#from balanced_basic import *
    3030from balanced_dev import *
     31#from anuga_tsunami import *
    3132
    3233from anuga import Reflective_boundary

  • trunk/anuga_work/development/gareth/tests/okushiri/run_okushiri.py

    r8354 r8547  
    2424import anuga
    2525import project
    26 from balanced_dev import *
     26#from balanced_dev import *
     27from anuga_tsunami import *
    2728
    2829def main(elevation_in_mesh=False):

  • trunk/anuga_work/development/gareth/tests/parabolic/parabolic.py

    r8353 r8547  
    99
    1010from balanced_dev import *
     11#from anuga_tsunami import *
    1112#from balanced_basic import *
    1213#from anuga.shallow_water_balanced2.swb2_domain import Domain as Domain
    … …  
    2122domain.set_datadir('.')                          # Use current folder
    2223domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
    23 domain.set_minimum_allowed_height(0.01)
     24domain.set_minimum_allowed_height(0.001)
    2425# Time stepping
    2526#domain.set_timestepping_method('euler') # Default

  • trunk/anuga_work/development/gareth/tests/runup/runup.py

    r8446 r8547  
    1919#from balanced_basic import *
    2020from balanced_dev import *
     21#from anuga_tsunami import *
    2122#---------
    2223#Setup computational domain
    … …  
    4546domain.get_quantity('elevation').smooth_vertex_values() # Steve's fix -- without this, substantial artificial velcities are generated everywhere in the domain. With this fix, there are artificial velocities near the coast, but not elsewhere.
    4647
    47 domain.set_quantity('friction',0.00)             # Constant friction
     48domain.set_quantity('friction',0.10)             # Constant friction
    4849
    4950domain.set_quantity('stage', stagefun)              # Constant negative initial stage
    … …  
    6162# Associate boundary tags with boundary objects
    6263#----------------------------------------------
    63 domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom':Br})
     64domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom':Br})
    6465
    6566#------------------------------

  • trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py

    r8466 r8547  
    1515#from balanced_basic import *
    1616from balanced_dev import *
     17#from anuga_tsunami import *
    1718#---------
    1819#Setup computational domain
    … …  
    2526domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
    2627#domain.set_store_vertices_uniquely(True)
    27 domain.minimum_allowed_height=0.001
    28 
    2928#------------------
    3029# Define topography
    … …  
    3837
    3938def stagefun(x,y):
    40     stge=-0.2*scale_me #-0.1*(x>0.5) -0.2*(x<=0.5)
     39    stge=-0.2*scale_me # +0.01*(x>0.9)
    4140    #topo=topography(x,y)
    4241    return stge#*(stge>topo) + (topo)*(stge<=topo)
    … …  
    4443domain.set_quantity('elevation',topography)     # Use function for elevation
    4544domain.get_quantity('elevation').smooth_vertex_values()
    46 
    47 domain.set_quantity('friction',0.00)             # Constant friction
     45domain.set_quantity('friction',0.03)             # Constant friction
    4846
    4947#def frict_change(x,y):
    … …  
    7472# Associate boundary tags with boundary objects
    7573#----------------------------------------------
    76 domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom':Br})
     74domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom':Br})
    7775
    7876#------------------------------
    … …  
    8078#------------------------------
    8179
    82 for t in domain.evolve(yieldstep=0.1,finaltime=20.0):
     80for t in domain.evolve(yieldstep=0.2,finaltime=20.0):
    8381    print domain.timestepping_statistics()
    8482    xx = domain.quantities['xmomentum'].centroid_values
    8583    yy = domain.quantities['ymomentum'].centroid_values
    8684    dd = domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values
    87     dd = (dd)*(dd>1.0e-06)+1.0e-03
    88     vv = ( (xx/dd)**2 + (yy/dd)**2 )**0.5
     85    dd_raw=1.0*dd
     86    dd = (dd)*(dd>1.0e-03)+1.0e-03
     87    vv = ( (xx/dd)**2 + (yy/dd)**2)**0.5
    8988    vv = vv*(dd>1.0e-03)
    9089    print 'Peak velocity is: ', vv.max(), vv.argmax()
     90    print 'Volume is', sum(dd_raw*domain.areas)   
     91
     92
     93vel_final_inds=(vv>1.0e-01).nonzero()
     94print 'vel_final_inds', vel_final_inds
    9195
    9296print 'Finished'

  • trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_2plot.py

    r8466 r8547  
    1717
    1818#------------------
    19 # Select line
     19# Select line along which we plot
    2020#------------------
    2121#v=(y==2.5)
    22 v=(p2.y==p2.y[3])
     22tmp = abs(p2.y-50.)
     23n = tmp.argmin()
     24v=(p2.y==p2.y[n])
    2325
    2426#-------------------------------

  • trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_sparse.py

    r8359 r8547  
    33"""
    44import sys
    5 
     5import numpy
    66#---------------
    77# I am having trouble with my environmental variables -- this is a
    … …  
    2424from balanced_dev import *
    2525from balanced_dev import Domain as Domain
     26#from anuga_tsunami import *
     27#from anuga_tsunami import Domain as Domain
    2628#------------------------------------------------------------------------------
    2729# Setup computational domain
    … …  
    3739#------------------------------------------------------------------------------
    3840def topography(x, y):
    39         return -x/10. #+abs(y-2.5)/10. # linear bed slope
     41        return -x/10. + numpy.sin(x/10.) +abs(y-50.)/100. # linear bed slope
    4042
    4143def stagetopo(x,y):
    … …  
    4446    return stg#*(stg>topo) + topo*(stg<=topo)
    4547
    46 line1=[ [0.,0.], [0., 100.] ]
     48line1=[ [10.,0.], [10., 100.] ]
    4749Qin=20.
    4850Inlet_operator(domain, line1,Qin)
    … …  
    6062
    6163def constant_water(t):
    62     return -9.936903
     64    return -8.936903
    6365
    64 #Bo= anuga.shallow_water.boundaries.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain, constant_water) # Outflow for steady uniform flow with a depth integrated velocity of 0.1
    65 Bo = anuga.Transmissive_boundary(domain)
     66Bo= anuga.shallow_water.boundaries.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain, constant_water) # Outflow for steady uniform flow with a depth integrated velocity of 0.1
     67#Bo = anuga.Transmissive_boundary(domain)
     68#Bd = anuga.Dirichlet_boundary([-20., 0., 0.])
    6669domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br})
     70
    6771#------------------------------------------------------------------------------
    6872# Evolve system through time

  • trunk/anuga_work/development/gareth/tests/urban_flow/ideal_urban.py

    r8450 r8547  
    1212#from anuga.structures.inlet_operator import Inlet_operator
    1313
    14 from balanced_dev import *
    15 from balanced_dev import create_domain_from_regions as create_domain_from_regions
     14#from balanced_dev import *
     15#from balanced_dev import create_domain_from_regions as create_domain_from_regions
     16from anuga_tsunami import *
     17from anuga_tsunami import create_domain_from_regions as create_domain_from_regions
    1618#------------------------------------------------------------------------------
    1719# Useful parameters for controlling this case

  • trunk/anuga_work/development/gareth/tests/wave/run_wave.py

    r8446 r8547  
    1414#from anuga import Domain
    1515
    16 from balanced_dev import *
     16#from balanced_dev import *
     17from anuga_tsunami import *
    1718#from balanced_dev import Domain as Domain
    1819
    … …  
    104105
    105106Bw2 = anuga.shallow_water.boundaries.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain, waveform)
    106 Bf_in = swb2_boundary_conditions.Characteristic_boundary(domain, waveform)
     107#Bf_in = swb2_boundary_conditions.Characteristic_boundary(domain, waveform)
    107108#Bw3 = swb2_boundary_conditions.Transmissive_momentum_nudge_stage_boundary(domain, waveform)
    108 Bf=swb2_boundary_conditions.Characteristic_boundary(domain,waveform2)
     109#Bf=swb2_boundary_conditions.Characteristic_boundary(domain,waveform2)
    109110                   
    110111
    111112
    112113# Associate boundary tags with boundary objects
    113 domain.set_boundary({'left': Bw2, 'right': Bf, 'top': Br, 'bottom': Br})
     114domain.set_boundary({'left': Bw2, 'right': Bt, 'top': Br, 'bottom': Br})
    114115
    115116
Note: See TracChangeset for help on using the changeset viewer.