Changeset 9016 for trunk


Ignore:
Timestamp:
Nov 7, 2013, 10:05:45 AM (11 years ago)
Author:
davies
Message:

Cleaning dev_audusse -- some problems remain

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

Legend:

Unmodified
Added
Removed

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

    r9014 r9016  
    6666        self.set_CFL(1.00)
    6767        self.set_use_kinematic_viscosity(False)
    68         self.timestepping_method='rk2'#'rk3'#'euler'#'rk2'
     68        self.timestepping_method='rk2'#'rk2'#'rk3'#'euler'#'rk2'
    6969        # The following allows storage of the negative depths associated with this method
    7070        self.minimum_storable_height=-99999999999.0

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

    r9015 r9016  
    5151}
    5252
    53 // Function to obtain speed from momentum and depth.
    54 // This is used by flux functions
    55 // Input parameters uh and h may be modified by this function.
    56 // Tried to inline, but no speedup was achieved 27th May 2009 (Ole)
    57 //static inline double _compute_speed(double *uh,
    58 double _compute_speed(double *uh,
    59                       double *h,
    60                       double epsilon,
    61                       double h0,
    62                       double limiting_threshold) {
    63  
    64   double u;
    65 
    66   if (*h < limiting_threshold) {   
    67     // Apply limiting of speeds according to the ANUGA manual
    68     if (*h < epsilon) {
    69       //*h = max(0.0,*h);  // Could have been negative
    70       u = 0.0;
    71     } else {
    72       u = *uh/(*h + h0/ *h);   
    73     }
    74  
    75 
    76     // Adjust momentum to be consistent with speed
    77     *uh = u * *h;
    78   } else {
    79     // We are in deep water - no need for limiting
    80     u = *uh/ *h;
    81   }
    82  
    83   return u;
    84 }
    85 
    86 // minmod limiter
    87 int _minmod(double a, double b){
    88     // Compute minmod
    89 
    90     if(sign(a)!=sign(b)){
    91         return 0.0;
    92     }else{
    93         return min(fabs(a), fabs(b))*sign(a);
    94     }
    95 
    96 
    97 }
     53//// Function to obtain speed from momentum and depth.
     54//// This is used by flux functions
     55//// Input parameters uh and h may be modified by this function.
     56//// Tried to inline, but no speedup was achieved 27th May 2009 (Ole)
     57////static inline double _compute_speed(double *uh,
     58//double _compute_speed(double *uh,
     59//                    double *h,
     60//                    double epsilon,
     61//                    double h0,
     62//                    double limiting_threshold) {
     63// 
     64//  double u;
     65//
     66//  if (*h < limiting_threshold) {   
     67//    // Apply limiting of speeds according to the ANUGA manual
     68//    if (*h < epsilon) {
     69//      //*h = max(0.0,*h);  // Could have been negative
     70//      u = 0.0;
     71//    } else {
     72//      u = *uh/(*h + h0/ *h);   
     73//    }
     74// 
     75//
     76//    // Adjust momentum to be consistent with speed
     77//    *uh = u * *h;
     78//  } else {
     79//    // We are in deep water - no need for limiting
     80//    u = *uh/ *h;
     81//  }
     82// 
     83//  return u;
     84//}
     85//
     86//// minmod limiter
     87//int _minmod(double a, double b){
     88//    // Compute minmod
     89//
     90//    if(sign(a)!=sign(b)){
     91//        return 0.0;
     92//    }else{
     93//        return min(fabs(a), fabs(b))*sign(a);
     94//    }
     95//
     96//
     97//}
    9898
    9999// Innermost flux function (using stage w=z+h)
    … …  
    149149  uh_left=q_left_rotated[1];
    150150  vh_left=q_left_rotated[2];
    151   if(hle>1.0e-03){
     151  if(hle>0.0e-06){
    152152    u_left = uh_left/(hle+0.0e-03) ; //max(h_left, 1.0e-06);
    153153    uh_left=h_left*u_left;
    … …  
    165165  uh_right = q_right_rotated[1];
    166166  vh_right = q_right_rotated[2];
    167   if(hre>1.0e-03){
     167  if(hre>0.0e-06){
    168168    u_right = uh_right/(hre+0.0e-03);//max(h_right, 1.0e-06);
    169169    uh_right=h_right*u_right;
    … …  
    241241  // Flux computation
    242242  denom = s_max - s_min;
    243   //if (denom < epsilon)
    244   if (0==1)
     243  if (denom < epsilon)
     244  //if (0==1)
    245245  {
    246246    // Both wave speeds are very small
    … …  
    263263
    264264      //if(i==0){
    265       //    // Standard smoothing term
    266       //    edgeflux[i] += (s_max*s_min)*(h_right - h_left);
     265      //    edgeflux[i] += (s_max*s_min)*(hre - hle);
    267266      //}else{
    268267          // Standard smoothing term
    … …  
    361360    // Compute previous call max_speed
    362361    speed_max_last=0.0;
    363     for (k = 0; k < number_of_elements; k++){
    364        
    365         speed_max_last=max(speed_max_last, max_speed_array[k]);
    366     }
     362    //for (k = 0; k < number_of_elements; k++){
     363    //   
     364    //    speed_max_last=max(speed_max_last, max_speed_array[k]);
     365    //}
    367366
    368367    // For all triangles
    … …  
    406405            } else {
    407406                // Neighbour is a real triangle
    408                 hc_n = max(stage_centroid_values[n] - bed_centroid_values[n],0.0);
     407                hc_n = height_centroid_values[n];//max(stage_centroid_values[n] - bed_centroid_values[n],0.0);
    409408                zc_n = bed_centroid_values[n];
    410409                m = neighbour_edges[ki];
    … …  
    426425            h_right= max(hre+zr-z_half,0.);
    427426
    428             //if (fabs(zl-zr)>1.0e-10) {
    429             //    report_python_error(AT,"Discontinuous Elevation");
    430             //    return 0.0;
    431             //}
    432            
    433 
    434            
    435427            // Edge flux computation (triangle k, edge i)
    436428            _flux_function_central(ql, qr,
    … …  
    448440            edgeflux[2] *= length;
    449441
    450             // Don't allow an outward advective flux if the cell centroid stage
    451             // is < the edge value
    452             //if(((stage_centroid_values[k] < z_half)| hle<H0) && edgeflux[0] > 0.){
     442            //// Don't allow an outward advective flux if the cell centroid stage
     443            //// is < the edge value
     444            //if((hc<H0) && edgeflux[0] > 0.){
    453445            //    edgeflux[0] = 0.;
    454446            //    edgeflux[1] = 0.;
    … …  
    458450            //}
    459451            ////
    460             //if(((stage_centroid_values[n] < z_half)| hre<H0) && edgeflux[0] < 0.){
     452            //if((hc_n<H0) && edgeflux[0] < 0.){
    461453            //    edgeflux[0] = 0.;
    462454            //    edgeflux[1] = 0.;
    … …  
    472464
    473465            // bedslope_work contains all gravity related terms -- weighting of
    474             bedslope_work=length*(-g*0.5*(h_left*h_left - hle*hle -(hle+hc)*(zl-zc))+pressure_flux);
     466            //if(h_left>0. || h_right>0.){
     467                bedslope_work=length*(-g*0.5*(h_left*h_left - hle*hle -(hle+hc)*(zl-zc))+pressure_flux);
     468            //}else{
     469            //    bedslope_work=0.;
     470            //}
    475471
    476472            pressuregrad_store[ki]=bedslope_work;
    … …  
    499495                edgeflux_store[nm3 + 1 ] = edgeflux[1];
    500496                edgeflux_store[nm3 + 2 ] = edgeflux[2];
    501 
    502                 bedslope_work=length*(-g*0.5*(h_right*h_right-hre*hre-(hre+hc_n)*(zr-zc_n))+pressure_flux);
     497                //if(h_right>0. || h_left >0.){
     498                    bedslope_work=length*(-g*0.5*(h_right*h_right-hre*hre-(hre+hc_n)*(zr-zc_n))+pressure_flux);
     499                //}else{
     500                //    bedslope_work=0.;
     501                //}
    503502                pressuregrad_store[nm]=bedslope_work;
    504503
    … …  
    529528                }
    530529            }
     530       
     531        // Keep track of maximal speeds
     532        //max_speed_array[k] = max(max_speed_array[k],max_speed);
    531533
    532534        } // End edge i (and neighbour n)
    533 
    534535        // Keep track of maximal speeds
     536        // FIXME: Doesn't this only store edge 2?? Not the max-speed of the
     537        // triangle as a whole. Not sure if it is used anyway
    535538        max_speed_array[k] = max_speed;
     539
    536540
    537541    } // End triangle k
    … …  
    541545    for(k=0; k< number_of_elements; k++){
    542546            //continue;
    543             hc = max(stage_centroid_values[k] - bed_centroid_values[k],0.);
     547            hc = height_centroid_values[k]; //max(stage_centroid_values[k] - bed_centroid_values[k],0.);
    544548            // Loop over every edge
    545549            for(i = 0; i<3; i++){
    546550                if(i==0){
    547                     // Add up the outgoing flux through the cell
     551                    // Add up the outgoing flux through the cell -- only do this once (i==0)
    548552                    outgoing_mass_edges=0.0;
    549553                    for(useint=0; useint<3; useint++){
    … …  
    624628            // GD HACK
    625629            // Compute bed slope term
    626             //if(hc > H0*0.5){
     630            //if(hc > H0){
    627631                xmom_explicit_update[k] -= normals[ki2]*pressuregrad_store[ki];
    628632                ymom_explicit_update[k] -= normals[ki2+1]*pressuregrad_store[ki];
    … …  
    640644        xmom_explicit_update[k] *= inv_area;
    641645        ymom_explicit_update[k] *= inv_area;
     646   
     647        //GD HACK
     648        //if(k==30) printf("%e \n", edgeflux_store[ki3]);
    642649    }  // end cell k
     650
    643651
    644652    if((call-2)%timestep_order==0) timestep=local_timestep; // Hack to ensure we only update the timestep on the first call within each rk2/rk3 step
    … …  
    728736  // distance between the bed_centroid_value and the max bed_edge_value of
    729737  // every triangle.
    730   double minimum_relative_height=0.05;
     738  //double minimum_relative_height=0.05;
    731739  int mass_added=0;
    732740
    … …  
    735743    for (k=0; k<N; k++) {
    736744      hc = wc[k] - zc[k];
    737       // Definine the maximum bed edge value on triangle k.
    738       //bmax = 0.5*max((zv[3*k]+zv[3*k+1]),max((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
    739       //bmax = max(zv[3*k],max(zv[3*k+2],zv[3*k+1]));
    740 
    741       //if (hc < max(minimum_relative_height*(bmax-zc[k]), minimum_allowed_height)) {
    742       //  xmomc[k]*=0.1;
    743       //  ymomc[k]*=0.1;
    744       //}
    745 
    746 
    747       //if (hc < minimum_allowed_height*2.0 | hc < minimum_relative_height*(bmax-zc[k]) ){
    748       if (hc < minimum_allowed_height*2.0 ){
    749       //if (hc < minimum_allowed_height*2.0 ){
     745      if (hc < minimum_allowed_height*1.0 ){
    750746            // Set momentum to zero and ensure h is non negative
    751             //xmomc[k] = 0.;//xmomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
    752             //ymomc[k] = 0.;//ymomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
    753             //xmomc[k] = xmomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
    754             //ymomc[k] = ymomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
    755             //xmomc[k] = xmomc[k]*max(hc-minimum_allowed_height,0.)/(minimum_allowed_height);// 0.0;
    756             //ymomc[k] = ymomc[k]*max(hc-minimum_allowed_height,0.)/(minimum_allowed_height);// 0.0;
    757 
     747            xmomc[k] = 0.;//xmomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
     748            ymomc[k] = 0.;//ymomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
    758749        if (hc <= 0.0){
    759              // Definine the minimum bed edge value on triangle k.
    760              //bmin = 0.5*min((zv[3*k]+zv[3*k+1]),min((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
    761              //bmin =0.5*bmin + 0.5*min(zv[3*k],min(zv[3*k+1],zv[3*k+2]));
    762              //bmin =min(zv[3*k],min(zv[3*k+1],zv[3*k+2])) -minimum_allowed_height;
    763              bmin=zc[k]-minimum_allowed_height;//-1.0e-03;
     750             bmin=zc[k];//minimum_allowed_height;//-1.0e-03;
    764751             // Minimum allowed stage = bmin
    765752
    … …  
    769756                 mass_added=1; //Flag to warn of added mass               
    770757
    771                  wc[k] = max(wc[k], bmin);
     758                 wc[k] = bmin;
    772759                 //printf(" mass_add, %f, %f, %f,%f,%f,%d\n", xc[k], yc[k], wc[k],zc[k],wc[k]-zc[k], k) ;
    773760             
    … …  
    777764                 // time step, for 'dry' areas where the designated stage is
    778765                 // less than the bed centroid value
    779                  //wv[3*k] = max(wv[3*k], min(bmin, zv[3*k]-minimum_allowed_height));
    780                  //wv[3*k+1] = max(wv[3*k+1], min(bmin, zv[3*k+1]-minimum_allowed_height));
    781                  //wv[3*k+2] = max(wv[3*k+2], min(bmin, zv[3*k+2]-minimum_allowed_height));
    782                  wv[3*k] = min(bmin, wc[3*k]); //zv[3*k]-minimum_allowed_height);
    783                  wv[3*k+1] = min(bmin, wc[3*k+1]); //zv[3*k+1]-minimum_allowed_height);
    784                  wv[3*k+2] = min(bmin, wc[3*k+2]); //zv[3*k+2]-minimum_allowed_height);
     766                 wv[3*k] = min(bmin, wc[k]); //zv[3*k]-minimum_allowed_height);
     767                 wv[3*k+1] = min(bmin, wc[k]); //zv[3*k+1]-minimum_allowed_height);
     768                 wv[3*k+2] = min(bmin, wc[k]); //zv[3*k+2]-minimum_allowed_height);
    785769            }
    786770        }
    … …  
    928912  }
    929913
    930   // Alternative 'PROTECT' step
    931   for(k=0; k<number_of_elements;k++){
    932     if((cell_wetness_scale[k]==0. ) ){
    933         xmom_centroid_values[k]=0.;
    934         ymom_centroid_values[k]=0.;
    935         //xmom_centroid_values[k]*=0.9;
    936         //ymom_centroid_values[k]*=0.9;
    937 
    938     }
    939 
    940    
    941     ////// Try some Froude-number limiting in shallow depths
    942     //dpth=max(stage_centroid_values[k] - elevation_centroid_values[k], 0.0);
    943     ////
    944     //if(dpth< max(max_elevation_edgevalue[k]-min_elevation_edgevalue[k],10*minimum_allowed_height)){
    945     //    // momnorm = momentum^2
    946     //    momnorm=(xmom_centroid_values[k]*xmom_centroid_values[k]+
    947     //             ymom_centroid_values[k]*ymom_centroid_values[k]);
    948     //   
    949     //    // vel^2 < constant*g*dpth [-- then multiply both sides by dpth^2]
    950     //    if(momnorm > 1*9.81*dpth*dpth*dpth){
    951     //        // Down-scale momentum so that Froude number < constant
    952     //        tmp=sqrt((1*9.81*dpth*dpth*dpth)/momnorm);
    953     //        xmom_centroid_values[k] *=tmp;
    954     //        ymom_centroid_values[k] *=tmp;
    955     //    }
    956     //}
    957   }
     914  //// Alternative 'PROTECT' step
     915  //for(k=0; k<number_of_elements;k++){
     916  //  //if((cell_wetness_scale[k]==0. ) ){
     917  //  //    xmom_centroid_values[k]=0.;
     918  //  //    ymom_centroid_values[k]=0.;
     919  //  //    //xmom_centroid_values[k]*=0.9;
     920  //  //    //ymom_centroid_values[k]*=0.9;
     921
     922  //  //}
     923
     924  // 
     925  //  ////// Try some Froude-number limiting in shallow depths
     926  //  //dpth=max(stage_centroid_values[k] - elevation_centroid_values[k], 0.0);
     927  //  ////
     928  //  //if(dpth< max(max_elevation_edgevalue[k]-min_elevation_edgevalue[k],10*minimum_allowed_height)){
     929  //  //    // momnorm = momentum^2
     930  //  //    momnorm=(xmom_centroid_values[k]*xmom_centroid_values[k]+
     931  //  //             ymom_centroid_values[k]*ymom_centroid_values[k]);
     932  //  //   
     933  //  //    // vel^2 < constant*g*dpth [-- then multiply both sides by dpth^2]
     934  //  //    if(momnorm > 4*9.81*dpth*dpth*dpth){
     935  //  //        // Down-scale momentum so that Froude number < constant
     936  //  //        tmp=sqrt((4*9.81*dpth*dpth*dpth)/momnorm);
     937  //  //        xmom_centroid_values[k] *=tmp;
     938  //  //        ymom_centroid_values[k] *=tmp;
     939  //  //    }
     940  //  //}
     941  //}
    958942
    959943 
    … …  
    963947      // extrapolation This will be changed back at the end of the routine
    964948      for (k=0; k<number_of_elements; k++){
    965 
    966           dk = max(stage_centroid_values[k]-elevation_centroid_values[k],minimum_allowed_height);
    967           xmom_centroid_store[k] = xmom_centroid_values[k];
    968           xmom_centroid_values[k] = xmom_centroid_values[k]/dk;
    969 
    970           ymom_centroid_store[k] = ymom_centroid_values[k];
    971           ymom_centroid_values[k] = ymom_centroid_values[k]/dk;
    972 
    973 
    974           // SET HEIGHT IN PLACE
     949         
    975950          height_centroid_values[k] = max(stage_centroid_values[k] - elevation_centroid_values[k], 0.);
    976            
    977           }
    978 
     951
     952          dk = height_centroid_values[k]; //max(stage_centroid_values[k]-elevation_centroid_values[k],minimum_allowed_height);
     953          if(dk>minimum_allowed_height){
     954              xmom_centroid_store[k] = xmom_centroid_values[k];
     955              xmom_centroid_values[k] = xmom_centroid_values[k]/dk;
     956
     957              ymom_centroid_store[k] = ymom_centroid_values[k];
     958              ymom_centroid_values[k] = ymom_centroid_values[k]/dk;
     959          }else{
     960              xmom_centroid_store[k] = 0.;
     961              xmom_centroid_values[k] = 0.;
     962              ymom_centroid_store[k] = 0.;
     963              ymom_centroid_values[k] = 0.;
     964
     965         }
    979966      }
     967  }
    980968
    981969  // Begin extrapolation routine
    … …  
    1004992      ymom_edge_values[k3+1]  = ymom_centroid_values[k];
    1005993      ymom_edge_values[k3+2]  = ymom_centroid_values[k];
    1006       dk = max(stage_centroid_values[k] - elevation_centroid_values[k], 0.);
     994
     995      dk = height_centroid_values[k];//max(stage_centroid_values[k] - elevation_centroid_values[k], 0.);
    1007996      height_edge_values[k3] = dk;
    1008997      height_edge_values[k3+1] = dk;
    … …  
    10581047      k1 = surrogate_neighbours[k3 + 1];
    10591048      k2 = surrogate_neighbours[k3 + 2];
     1049
     1050      //if(cell_wetness_scale[k0]==0 && cell_wetness_scale[k1]==0 && cell_wetness_scale[k2]==0){
     1051      //    xmom_centroid_store[k]=0.;
     1052      //    ymom_centroid_store[k]=0.;
     1053      //    xmom_centroid_values[k] = 0.;
     1054      //    ymom_centroid_values[k] = 0.;
     1055      //}
    10601056     
    10611057      // Test to see whether we accept the surrogate neighbours
    … …  
    10641060      // used
    10651061      // FIXME: Remove cell_wetness_scale if you don't need it
    1066       if(k2<0 || (cell_wetness_scale[k2]==-10.0 && stage_centroid_values[k]<max_elevation_edgevalue[k])){
    1067           k2 = k ;
    1068       }
    1069       if(k0 < 0 || (cell_wetness_scale[k0]==-10.0 && stage_centroid_values[k]<max_elevation_edgevalue[k])){
    1070           k0 = k ;
    1071       }
    1072       if(k1 < 0 || (cell_wetness_scale[k1]==-10.0 && stage_centroid_values[k]<max_elevation_edgevalue[k])){
    1073           k1 = k ;
    1074       }
     1062      //if( (cell_wetness_scale[k2]==0.0 && stage_centroid_values[k]<max_elevation_edgevalue[k]+100.)){
     1063      //    k2 = k ;
     1064      //}
     1065      //if((cell_wetness_scale[k0]==0.0 && stage_centroid_values[k]<max_elevation_edgevalue[k])+100.){
     1066      //    k0 = k ;
     1067      //}
     1068      //if((cell_wetness_scale[k1]==0.0 && stage_centroid_values[k]<max_elevation_edgevalue[k])+100.){
     1069      //    k1 = k ;
     1070      //}
    10751071
    10761072      // Take note if the max neighbour bed elevation is greater than the min
    … …  
    11131109      // The triangle is guaranteed to be counter-clockwise     
    11141110      area2 = dy2*dx1 - dy1*dx2;
     1111      //if(k==54) printf("K=54\n");
    11151112       
    1116       //// Treat triangles with zero or 1 wet neighbours (area2 <=0.)
    1117       if ((area2 <= 0.) )//|| (cell_wetness_scale[k]==0.)) //|(count_wet_neighbours[k]==0))
     1113      //// Treat triangles with no neighbours (area2 <=0.)
     1114      if ((area2 <= 0.))//|( cell_wetness_scale[k2]==0. && cell_wetness_scale[k1]==0. && cell_wetness_scale[k0]==0.))//|| (cell_wetness_scale[k]==0.)) //|(count_wet_neighbours[k]==0))
    11181115      {
    11191116
     1117
    11201118          // Isolated wet cell -- constant stage/depth extrapolation
    1121           //stage_edge_values[k3]   = stage_centroid_values[k];
    1122           //stage_edge_values[k3+1] = stage_centroid_values[k];
    1123           //stage_edge_values[k3+2] = stage_centroid_values[k];
     1119          stage_edge_values[k3]   = stage_centroid_values[k];
     1120          stage_edge_values[k3+1] = stage_centroid_values[k];
     1121          stage_edge_values[k3+2] = stage_centroid_values[k];
    11241122
    11251123          dk=height_centroid_values[k]; //max(stage_centroid_values[k]-elevation_centroid_values[k],0.);
    … …  
    11281126          height_edge_values[k3+2] = dk;
    11291127         
    1130           stage_edge_values[k3]   = elevation_centroid_values[k]+dk;
    1131           stage_edge_values[k3+1] = elevation_centroid_values[k]+dk;
    1132           stage_edge_values[k3+2] = elevation_centroid_values[k]+dk;
     1128          //stage_edge_values[k3]   = elevation_centroid_values[k]+dk;
     1129          //stage_edge_values[k3+1] = elevation_centroid_values[k]+dk;
     1130          //stage_edge_values[k3+2] = elevation_centroid_values[k]+dk;
    11331131          //stage_edge_values[k3] = elevation_edge_values[k3]+dk;
    11341132          //stage_edge_values[k3+1] = elevation_edge_values[k3+1]+dk;
    11351133          //stage_edge_values[k3+2] = elevation_edge_values[k3+2]+dk;
     1134
     1135          //xmom_centroid_values[k]=0.;
     1136          //ymom_centroid_values[k]=0.;
     1137          //xmom_centroid_store[k] = 0.;
     1138          //ymom_centroid_store[k] = 0.;
    11361139
    11371140          xmom_edge_values[k3]    = xmom_centroid_values[k];
    … …  
    11461149     
    11471150      // Calculate heights of neighbouring cells
    1148       hc = stage_centroid_values[k]  - elevation_centroid_values[k];
    1149       h0 = stage_centroid_values[k0] - elevation_centroid_values[k0];
    1150       h1 = stage_centroid_values[k1] - elevation_centroid_values[k1];
    1151       h2 = stage_centroid_values[k2] - elevation_centroid_values[k2];
     1151      hc = height_centroid_values[k];//stage_centroid_values[k]  - elevation_centroid_values[k];
     1152      h0 = height_centroid_values[k0];// - elevation_centroid_values[k0];
     1153      h1 = height_centroid_values[k1];// - elevation_centroid_values[k1];
     1154      h2 = height_centroid_values[k2];// - elevation_centroid_values[k2];
    11521155     
    11531156      hmin = min(min(h0, min(h1, h2)), hc);
    … …  
    11551158
    11561159      // Look for strong changes in cell depth, or shallow cell depths, as an indicator of near-wet-dry
    1157       hfactor= max(0., min(2.0*max(hmin,0.0)/max(hc,1.0e-06)-0.1,
    1158                            min(2.0*max(hc,0.)/max(hmax,1.0e-06)-0.1, 1.0))
     1160      // Reduce beta linearly from 1-0 between depth ratio of 0.3-0.1
     1161      hfactor= max(0., min(5.0*max(hmin,0.0)/max(hc,1.0e-06)-0.5,
     1162                           min(5.0*max(hc,0.)/max(hmax,1.0e-06)-0.5, 1.0))
    11591163                  );
     1164      //hfactor= max(0., min(5.0*max(hmin,0.0)/max(hmax,1.0e-06)-0.5, 1.0)
     1165      //            );
    11601166      //hfactor=1.0;
    11611167      hfactor=min( max(hmin,0.)/(max(hmin,0.)+1.*minimum_allowed_height), hfactor);
    … …  
    12331239      // from the centroid to the min and max
    12341240      find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax);
    1235    
    1236       beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
    1237       //beta_tmp=beta_w;
    1238      
     1241   
    12391242      // Limit the gradient
    12401243      limit_gradient(dqv, qmin, qmax, beta_tmp);
    … …  
    12461249
    12471250      // REDEFINE hfactor for momentum terms -- make MORE first order
    1248       hfactor= max(0., min(1.6*max(hmin,0.0)/max(hc,1.0e-06)-0.5,
    1249                            min(1.6*max(hc,0.)/max(hmax,1.0e-06)-0.5, 1.0))
     1251      // Reduce beta linearly from 1-0 between depth ratio of 0.6-0.4
     1252      hfactor= max(0., min(5*max(hmin,0.0)/max(hc,1.0e-06)-2.0,
     1253                           min(5*max(hc,0.)/max(hmax,1.0e-06)-2.0, 1.0))
    12501254                  );
     1255      //hfactor= max(0., min(5*max(hmin,0.0)/max(hmax,1.0e-06)-2.0,
     1256      //                      1.0));
    12511257      hfactor=min( max(hmin,0.)/(max(hmin,0.)+10.*minimum_allowed_height), hfactor);
    12521258     
    … …  
    14911497      limit_gradient(dqv, qmin, qmax, beta_w);
    14921498     
    1493       //for (i=0;i<3;i++)
    1494       //xmom_edge_values[k3] = xmom_centroid_values[k] + dqv[0];
    1495       //xmom_edge_values[k3 + 1] = xmom_centroid_values[k] + dqv[1];
    1496       //xmom_edge_values[k3 + 2] = xmom_centroid_values[k] + dqv[2];
    1497      
    14981499      for (i = 0; i < 3;i++)
    14991500      {
    … …  
    15551556      height_vertex_values[k3+2] =  height_edge_values[k3] + height_edge_values[k3+1]-height_edge_values[k3+2];
    15561557
    1557      // Compute xmom vertex values
    1558       xmom_vertex_values[k3] = xmom_edge_values[k3+1] + xmom_edge_values[k3+2] -xmom_edge_values[k3] ;
    1559       xmom_vertex_values[k3+1] =  xmom_edge_values[k3] + xmom_edge_values[k3+2]-xmom_edge_values[k3+1];
    1560       xmom_vertex_values[k3+2] =  xmom_edge_values[k3] + xmom_edge_values[k3+1]-xmom_edge_values[k3+2];
    1561 
    1562       // Compute ymom vertex values
    1563       ymom_vertex_values[k3] = ymom_edge_values[k3+1] + ymom_edge_values[k3+2] -ymom_edge_values[k3] ;
    1564       ymom_vertex_values[k3+1] =  ymom_edge_values[k3] + ymom_edge_values[k3+2]-ymom_edge_values[k3+1];
    1565       ymom_vertex_values[k3+2] =  ymom_edge_values[k3] + ymom_edge_values[k3+1]-ymom_edge_values[k3+2];
    1566 
    15671558      // If needed, convert from velocity to momenta
    15681559      if(extrapolate_velocity_second_order==1){
    … …  
    15701561          xmom_centroid_values[k] = xmom_centroid_store[k];
    15711562          ymom_centroid_values[k] = ymom_centroid_store[k];
    1572 
     1563     
    15731564          // Re-compute momenta at edges
    15741565          for (i=0; i<3; i++){
    1575               de[i] = height_edge_values[k3+i]; //max(stage_vertex_values[k3+i]-elevation_vertex_values[k3+i],0.0);
    1576               // Method to perhaps slightly reduce momenta
    1577               //de[i] = min(height_edge_values[k3+i], max(stage_edge_values[k3+i]-elevation_edge_values[k3+i],0.)); //max(stage_edge_values[k3+i]-elevation_edge_values[k3+i],0.0);
    1578               //if(de[i]<minimum_allowed_height){
    1579               //  de[i]=0.;
     1566              //if(hfactor>=0.8){
     1567              dk= height_edge_values[k3+i];
     1568              //}else{
     1569              //   de[i] = height_centroid_values[k];
    15801570              //}
    1581               xmom_edge_values[k3+i]=xmom_edge_values[k3+i]*de[i];
    1582               ymom_edge_values[k3+i]=ymom_edge_values[k3+i]*de[i];
     1571              xmom_edge_values[k3+i]=xmom_edge_values[k3+i]*dk;
     1572              ymom_edge_values[k3+i]=ymom_edge_values[k3+i]*dk;
    15831573          }
    1584 
    1585           // Re-compute momenta at vertices
    1586           for (i=0; i<3; i++){
    1587               de[i] = height_vertex_values[k3+i]; //max(stage_vertex_values[k3+i]-elevation_vertex_values[k3+i],0.0);
    1588               // Method to perhaps slightly reduce momenta
    1589               //de[i] = min(height_vertex_values[k3+i], max(stage_vertex_values[k3+i]-elevation_vertex_values[k3+i],0.)); //max(stage_vertex_values[k3+i]-elevation_vertex_values[k3+i],0.0);
    1590               xmom_vertex_values[k3+i]=xmom_vertex_values[k3+i]*de[i];
    1591               ymom_vertex_values[k3+i]=ymom_vertex_values[k3+i]*de[i];
    1592           }
    1593 
    15941574      }
     1575      // Compute momenta at vertices
     1576      xmom_vertex_values[k3] = xmom_edge_values[k3+1] + xmom_edge_values[k3+2] -xmom_edge_values[k3] ;
     1577      xmom_vertex_values[k3+1] =  xmom_edge_values[k3] + xmom_edge_values[k3+2]-xmom_edge_values[k3+1];
     1578      xmom_vertex_values[k3+2] =  xmom_edge_values[k3] + xmom_edge_values[k3+1]-xmom_edge_values[k3+2];
     1579      ymom_vertex_values[k3] = ymom_edge_values[k3+1] + ymom_edge_values[k3+2] -ymom_edge_values[k3] ;
     1580      ymom_vertex_values[k3+1] =  ymom_edge_values[k3] + ymom_edge_values[k3+2]-ymom_edge_values[k3+1];
     1581      ymom_vertex_values[k3+2] =  ymom_edge_values[k3] + ymom_edge_values[k3+1]-ymom_edge_values[k3+2];
     1582
     1583     
    15951584
    15961585      // Compute new bed elevation
    … …  
    15981587      elevation_edge_values[k3+1]=stage_edge_values[k3+1]-height_edge_values[k3+1];
    15991588      elevation_edge_values[k3+2]=stage_edge_values[k3+2]-height_edge_values[k3+2];
    1600       elevation_vertex_values[k3]=stage_vertex_values[k3]-height_vertex_values[k3];
    1601       elevation_vertex_values[k3+1]=stage_vertex_values[k3+1]-height_vertex_values[k3+1];
    1602       elevation_vertex_values[k3+2]=stage_vertex_values[k3+2]-height_vertex_values[k3+2];
     1589      elevation_vertex_values[k3] = elevation_edge_values[k3+1] + elevation_edge_values[k3+2] -elevation_edge_values[k3] ;
     1590      elevation_vertex_values[k3+1] =  elevation_edge_values[k3] + elevation_edge_values[k3+2]-elevation_edge_values[k3+1];
     1591      elevation_vertex_values[k3+2] =  elevation_edge_values[k3] + elevation_edge_values[k3+1]-elevation_edge_values[k3+2];
    16031592
    16041593   

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

    r8751 r9016  
    2424import anuga
    2525import project
    26 from balanced_dev import *
     26#from balanced_dev import *
     27from bal_and import *
    2728#from anuga_tsunami import *
    2829

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

    r9014 r9016  
    8686    return boundary_ws*scale_me
    8787    #return -0.2*scale_me #-0.1 #(0.0*sin(t*2*pi)-0.1)*exp(-t)-0.1
     88def waveform2(t):
     89    return [ (0.1*sin(2*pi*t)-0.3)*exp(-t), 0., 0.]
     90
    8891Bt2=anuga.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain,waveform)
    89 #Bw=anuga.Time_boundary(domain=domain,
    90 #       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.
     92Bw=anuga.Time_boundary(domain=domain, waveform2) # Time varying boundary -- get rid of the 0.0 to do a runup.
    9193
    9294#----------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.