Changeset 9039


Ignore:
Timestamp:
Dec 9, 2013, 9:18:55 AM (11 years ago)
Author:
davies
Message:

Cleaning up DE1 implementation

Location:
trunk
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r9038 r9039  
    452452        self.maximum_allowed_speed=0.0
    453453
    454         # Keep track of the fluxes through the boundaries
     454        ## FIXME: Should implement tracking of boundary fluxes
     455        ## Keep track of the fluxes through the boundaries
    455456        self.boundary_flux_integral=num.ndarray(1)
    456457        self.boundary_flux_integral[0]=0.
    … …  
    12441245            from swDE1_domain_ext import compute_fluxes_ext_central \
    12451246                                      as compute_fluxes_ext
    1246             if self.timestepping_method!='rk2':
    1247                 raise Exception, 'DE1 only works with rk2 at present, due to some hard-coded statements in set_DE1_defaults'
     1247            if self.timestepping_method=='euler':
     1248                raise Exception, "DE1 doesn't seems to work with euler at present \
     1249                                 (boundary conditions?), and is mostly designed for rk2"
    12481250            Stage = self.quantities['stage']
    12491251            Xmom = self.quantities['xmomentum']
    … …  
    12551257
    12561258            # Keep track of number of calls to compute_fluxes_ext
    1257             # This is a hacky method of knowing whether we are in the
    1258             # first or subsequent sub-steps of our 2nd/3rd order timestepping schemes,
    1259             # which is important for the mass conservation enforcement.
    12601259            # Note that in ANUGA, all sub-steps within an rk2/rk3 timestep have the same timestep
    12611260            self.call+=1
    … …  
    12951294                                               Bed.vertex_values)
    12961295
    1297             # Update the boundary flux integral
    1298             # FIXME GD: This should be moved to its own routine -- something like this
    1299             #        could be included in all solvers.
    1300             if(self.timestepping_method=='rk2'):
    1301                 if(self.call%2==1):
    1302                   self.boundary_flux_integral[0]= self.boundary_flux_integral[0] +\
    1303                                                     self.boundary_flux_sum[0]*self.timestep*0.5
    1304                   #print 'dbfi ', self.boundary_flux_integral, self.boundary_flux_sum
    1305                   self.boundary_flux_sum[0]=0.
    1306             elif(self.timestepping_method=='euler'):
    1307                 # FIXME GD: Unsupported at present.
    1308                 self.boundary_flux_integral=0.
    1309                 # This presently doesn't work -- this section of code may need to go elsewhere
    1310                 #self.boundary_flux_integral[0]= self.boundary_flux_integral[0] +\
    1311                 #                                  self.boundary_flux_sum[0]*self.timestep
    1312                 #self.boundary_flux_sum[0]=0.
    1313             elif(self.timestepping_method=='rk3'):
    1314                 # FIXME GD: Unsupported at present.
    1315                 self.boundary_flux_integral=0.
    1316                 # FIXME GD: This needs to be implemented
    1317             else:
    1318                 mess='ERROR: domain.timestepping_method', self.timestepping_method,' method not recognised'
    1319                 raise Exception, mess
     1296            ## FIXME: This won't work in parallel since the timestep has not been updated to include other processors.
     1297            ## Update the boundary flux integral
     1298            ## FIXME GD: This should be moved to its own routine -- something like this
     1299            ##        could be included in all solvers.
     1300            #if(self.timestepping_method=='rk2'):
     1301            #    if(self.call%2==1):
     1302            #      self.boundary_flux_integral[0]= self.boundary_flux_integral[0] +\
     1303            #                                        self.boundary_flux_sum[0]*self.timestep*0.5
     1304            #      #print 'dbfi ', self.boundary_flux_integral, self.boundary_flux_sum
     1305            #      self.boundary_flux_sum[0]=0.
     1306            #elif(self.timestepping_method=='euler'):
     1307            #    # FIXME GD: Unsupported at present.
     1308            #    self.boundary_flux_integral=0.
     1309            #    # This presently doesn't work -- this section of code may need to go elsewhere
     1310            #    #self.boundary_flux_integral[0]= self.boundary_flux_integral[0] +\
     1311            #    #                                  self.boundary_flux_sum[0]*self.timestep
     1312            #    #self.boundary_flux_sum[0]=0.
     1313            #elif(self.timestepping_method=='rk3'):
     1314            #    # FIXME GD: Unsupported at present.
     1315            #    self.boundary_flux_integral=0.
     1316            #    # FIXME GD: This needs to be implemented
     1317            #else:
     1318            #    mess='ERROR: domain.timestepping_method', self.timestepping_method,' method not recognised'
     1319            #    raise Exception, mess
    13201320
    13211321            self.flux_timestep = flux_timestep

  • trunk/anuga_core/source/anuga/shallow_water/swDE1_domain_ext.c

    r9038 r9039  
    287287        double* max_speed_array,
    288288        int optimise_dry_cells,
    289         int timestep_order,
     289        int timestep_fluxcalls,
    290290        double* stage_centroid_values,
    291291        double* xmom_centroid_values,
    … …  
    450450            // steps, NOT within them (since a constant timestep is used within
    451451            // each rk2/rk3 sub-step)
    452             if ((tri_full_flag[k] == 1) ) {
    453 
    454                 // On the 2nd/3rd timsteps of rk2 / rk3 sequence, don't
    455                 // recompute the max-speed. ANUGA uses a constant timestep and
    456                 // we need to respect that in the volume-protection
    457                 // computations below (achieved by fixing max-speed)
    458                 if((call-2)%timestep_order!=0) max_speed = max_speed_array[k]; // HACK to Ensure that local timestep is the same as the last timestep
    459                 // On the 1st timestep of an rk2/rk3 sequence, recompute the max-speed
    460                 if((call-2)%timestep_order==0) speed_max_last=max(speed_max_last, max_speed);
     452            if ((tri_full_flag[k] == 1) & ( (call-2)%timestep_fluxcalls==0)) {
     453
     454                speed_max_last=max(speed_max_last, max_speed);
    461455
    462456                if (max_speed > epsilon) {
    … …  
    479473        } // End edge i (and neighbour n)
    480474        // Keep track of maximal speeds
    481         if((call-2)%timestep_order==0) max_speed_array[k] = speed_max_last; //max_speed;
     475        if((call-2)%timestep_fluxcalls==0) max_speed_array[k] = speed_max_last; //max_speed;
    482476
    483477
    484478    } // End triangle k
    485479 
    486     // GD HACK
    487     // Limit edgefluxes, for mass conservation near wet/dry cells
    488     for(k=0; k< number_of_elements; k++){
    489         //continue;
    490         hc = height_centroid_values[k];
    491         // Loop over every edge
    492         for(i = 0; i<3; i++){
    493             if(i==0){
    494                 // Add up the outgoing flux through the cell -- only do this once (i==0)
    495                 outgoing_mass_edges=0.0;
    496                 for(useint=0; useint<3; useint++){
    497                     if(edgeflux_store[3*(3*k+useint)]< 0.){
    498                         //outgoing_mass_edges+=1.0;
    499                         outgoing_mass_edges+=(edgeflux_store[3*(3*k+useint)]);
    500                     }
    501                 }
    502                 outgoing_mass_edges*=local_timestep;
    503             }
    504 
    505             ki=3*k+i;   
    506             ki2=ki*2;
    507             ki3 = ki*3;
    508            
    509             // Prevent outflow from 'seriously' dry cells
    510             // Idea: The cell will not go dry if:
    511             // total_outgoing_flux <= cell volume = Area_triangle*hc
    512             vol=areas[k]*hc;
    513             if((edgeflux_store[ki3]< 0.0) && (-outgoing_mass_edges> vol)){
    514                
    515                 // This bound could be improved (e.g. we could actually sum the
    516                 // + and - fluxes and check if they are too large).  However,
    517                 // the advantage of this method is that we don't have to worry
    518                 // about subsequent changes to the + edgeflux caused by
    519                 // constraints associated with neighbouring triangles.
    520                 tmp = vol/(-(outgoing_mass_edges)) ;
    521                 if(tmp< 1.0){
    522                     edgeflux_store[ki3+0]*=tmp;
    523                     edgeflux_store[ki3+1]*=tmp;
    524                     edgeflux_store[ki3+2]*=tmp;
    525 
    526                     // Compute neighbour edge index
    527                     n = neighbours[ki];
    528                     if(n>=0){
    529                         nm = 3*n + neighbour_edges[ki];
    530                         nm3 = nm*3;
    531                         edgeflux_store[nm3+0]*=tmp;
    532                         edgeflux_store[nm3+1]*=tmp;
    533                         edgeflux_store[nm3+2]*=tmp;
    534                     }
    535                 }
    536             }
    537         }
    538     }
    539 
    540     //printf("%e \n", edgeflux_store[3*30*3]);
     480    //// GD HACK
     481    //// Limit edgefluxes, for mass conservation near wet/dry cells
     482    //for(k=0; k< number_of_elements; k++){
     483    //    //continue;
     484    //    hc = height_centroid_values[k];
     485    //    // Loop over every edge
     486    //    for(i = 0; i<3; i++){
     487    //        if(i==0){
     488    //            // Add up the outgoing flux through the cell -- only do this once (i==0)
     489    //            outgoing_mass_edges=0.0;
     490    //            for(useint=0; useint<3; useint++){
     491    //                if(edgeflux_store[3*(3*k+useint)]< 0.){
     492    //                    //outgoing_mass_edges+=1.0;
     493    //                    outgoing_mass_edges+=(edgeflux_store[3*(3*k+useint)]);
     494    //                }
     495    //            }
     496    //            outgoing_mass_edges*=local_timestep;
     497    //        }
     498
     499    //        ki=3*k+i;   
     500    //        ki2=ki*2;
     501    //        ki3 = ki*3;
     502    //       
     503    //        // Prevent outflow from 'seriously' dry cells
     504    //        // Idea: The cell will not go dry if:
     505    //        // total_outgoing_flux <= cell volume = Area_triangle*hc
     506    //        vol=areas[k]*hc;
     507    //        if((edgeflux_store[ki3]< 0.0) && (-outgoing_mass_edges> vol)){
     508    //           
     509    //            // This bound could be improved (e.g. we could actually sum the
     510    //            // + and - fluxes and check if they are too large).  However,
     511    //            // the advantage of this method is that we don't have to worry
     512    //            // about subsequent changes to the + edgeflux caused by
     513    //            // constraints associated with neighbouring triangles.
     514    //            tmp = vol/(-(outgoing_mass_edges)) ;
     515    //            if(tmp< 1.0){
     516    //                edgeflux_store[ki3+0]*=tmp;
     517    //                edgeflux_store[ki3+1]*=tmp;
     518    //                edgeflux_store[ki3+2]*=tmp;
     519
     520    //                // Compute neighbour edge index
     521    //                n = neighbours[ki];
     522    //                if(n>=0){
     523    //                    nm = 3*n + neighbour_edges[ki];
     524    //                    nm3 = nm*3;
     525    //                    edgeflux_store[nm3+0]*=tmp;
     526    //                    edgeflux_store[nm3+1]*=tmp;
     527    //                    edgeflux_store[nm3+2]*=tmp;
     528    //                }
     529    //            }
     530    //        }
     531    //    }
     532    // }
     533
     534    ////printf("%e \n", edgeflux_store[3*30*3]);
    541535
    542536    // Now add up stage, xmom, ymom explicit updates
    … …  
    591585
    592586    // Hack to ensure we only update the timestep on the first call within each rk2/rk3 step
    593     if((call-2)%timestep_order==0) timestep=local_timestep;
     587    if((call-2)%timestep_fluxcalls==0) timestep=local_timestep;
    594588
    595589    free(edgeflux_store);
    … …  
    16131607   
    16141608  double timestep, epsilon, H0, g;
    1615   int optimise_dry_cells, timestep_order;
     1609  int optimise_dry_cells, timestep_fluxcalls;
    16161610   
    16171611  // Convert Python arguments to C
    … …  
    16421636            &max_speed_array,
    16431637            &optimise_dry_cells,
    1644             &timestep_order,
     1638            &timestep_fluxcalls,
    16451639            &stage_centroid_values,
    16461640            &xmom_centroid_values,
    … …  
    17141708                     (double*) max_speed_array -> data,
    17151709                     optimise_dry_cells,
    1716                      timestep_order,
     1710                     timestep_fluxcalls,
    17171711                     (double*) stage_centroid_values -> data,
    17181712                     (double*) xmom_centroid_values -> data,

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

    r9038 r9039  
    6262        self.set_CFL(1.00)
    6363        self.set_use_kinematic_viscosity(False)
    64         self.timestepping_method='rk2'#'rk3'#'euler'#'rk2'
     64        self.set_timestepping_method('rk2')#'rk3'#'euler'#'rk2'
    6565
    6666        # Don't place any restriction on the minimum storable height
    … …  
    7373
    7474        self.beta_w=1.0
    75         self.beta_w_dry=0.0
     75        self.beta_w_dry=0.4
    7676        self.beta_uh=1.0
    7777        self.beta_uh_dry=0.0
    … …  
    9595
    9696        # Keep track of the fluxes through the boundaries
     97        # FIXME: Need to implement this
    9798        self.boundary_flux_integral=numpy.ndarray(1)
    9899        self.boundary_flux_integral[0]=0.
    … …  
    101102
    102103        self.call=1 # Integer counting how many times we call compute_fluxes_central
    103 
    104         # Integer recording the order of the time-stepping scheme
    105         if(self.timestepping_method=='rk2'):
    106           self.timestep_order=2
    107         elif(self.timestepping_method=='euler'):
    108           self.timestep_order=1
    109         elif(self.timestepping_method=='rk3'):
    110           self.timestep_order=3
    111         else:
    112           err_mess='ERROR: timestepping_method= ' + self.timestepping_method +' not supported in this solver'
    113           raise Exception, err_mess
    114104
    115105        print '##########################################################################'
    … …  
    254244                                           domain.max_speed,
    255245                                           int(domain.optimise_dry_cells),
    256                                            domain.timestep_order,
     246                                           domain.timestep_fluxcalls,
    257247                                           Stage.centroid_values,
    258248                                           Xmom.centroid_values,
    … …  
    263253
    264254
    265         # Update the boundary flux integral
    266         if(domain.timestepping_method=='rk2'):
    267             if(domain.call%2==1):
    268               domain.boundary_flux_integral[0]= domain.boundary_flux_integral[0] +\
    269                                                 domain.boundary_flux_sum[0]*domain.timestep*0.5
    270               #print 'dbfi ', domain.boundary_flux_integral, domain.boundary_flux_sum
    271               domain.boundary_flux_sum[0]=0.
    272         elif(domain.timestepping_method=='euler'):
    273             domain.boundary_flux_integral=0.
    274             # This presently doesn't work -- this section of code may need to go elsewhere
    275             #domain.boundary_flux_integral[0]= domain.boundary_flux_integral[0] +\
    276             #                                  domain.boundary_flux_sum[0]*domain.timestep
    277             #domain.boundary_flux_sum[0]=0.
    278         elif(domain.timestepping_method=='rk3'):
    279             domain.boundary_flux_integral=0.
    280             # FIXME: This needs to be implemented
    281         else:
    282             mess='ERROR: domain.timestepping_method', domain.timestepping_method,' method not recognised'
    283             raise Exception, mess
     255        ## FIXME: This won't work in parallel
     256        ## Update the boundary flux integral
     257        #if(domain.timestepping_method=='rk2'):
     258        #    if(domain.call%2==1):
     259        #      domain.boundary_flux_integral[0]= domain.boundary_flux_integral[0] +\
     260        #                                        domain.boundary_flux_sum[0]*domain.timestep*0.5
     261        #      #print 'dbfi ', domain.boundary_flux_integral, domain.boundary_flux_sum
     262        #      domain.boundary_flux_sum[0]=0.
     263        #elif(domain.timestepping_method=='euler'):
     264        #    domain.boundary_flux_integral=0.
     265        #    # This presently doesn't work -- this section of code may need to go elsewhere
     266        #    #domain.boundary_flux_integral[0]= domain.boundary_flux_integral[0] +\
     267        #    #                                  domain.boundary_flux_sum[0]*domain.timestep
     268        #    #domain.boundary_flux_sum[0]=0.
     269        #elif(domain.timestepping_method=='rk3'):
     270        #    domain.boundary_flux_integral=0.
     271        #    # FIXME: This needs to be implemented
     272        #else:
     273        #    mess='ERROR: domain.timestepping_method', domain.timestepping_method,' method not recognised'
     274        #    raise Exception, mess
    284275
    285276        domain.flux_timestep = flux_timestep

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

    r9038 r9039  
    288288        double* max_speed_array,
    289289        int optimise_dry_cells,
    290         int timestep_order,
     290        int timestep_fluxcalls,
    291291        double* stage_centroid_values,
    292292        double* xmom_centroid_values,
    … …  
    451451            // steps, NOT within them (since a constant timestep is used within
    452452            // each rk2/rk3 sub-step)
    453             if ((tri_full_flag[k] == 1) ) {
    454 
    455                 // On the 2nd/3rd timsteps of rk2 / rk3 sequence, don't
    456                 // recompute the max-speed. ANUGA uses a constant timestep and
    457                 // we need to respect that in the volume-protection
    458                 // computations below (achieved by fixing max-speed)
    459                 if((call-2)%timestep_order!=0) max_speed = max_speed_array[k]; // HACK to Ensure that local timestep is the same as the last timestep
    460                 // On the 1st timestep of an rk2/rk3 sequence, recompute the max-speed
    461                 if((call-2)%timestep_order==0) speed_max_last=max(speed_max_last, max_speed);
     453            if ((tri_full_flag[k] == 1) & ((call-2)%timestep_fluxcalls==0) ) {
     454
     455                speed_max_last=max(speed_max_last, max_speed);
    462456
    463457                if (max_speed > epsilon) {
    … …  
    480474        } // End edge i (and neighbour n)
    481475        // Keep track of maximal speeds
    482         if((call-2)%timestep_order==0) max_speed_array[k] = speed_max_last; //max_speed;
     476        if((call-2)%timestep_fluxcalls==0) max_speed_array[k] = speed_max_last; //max_speed;
    483477
    484478
    485479    } // End triangle k
    486480 
    487     // GD HACK
    488     // Limit edgefluxes, for mass conservation near wet/dry cells
    489     for(k=0; k< number_of_elements; k++){
    490         //continue;
    491         hc = height_centroid_values[k];
    492         // Loop over every edge
    493         for(i = 0; i<3; i++){
    494             if(i==0){
    495                 // Add up the outgoing flux through the cell -- only do this once (i==0)
    496                 outgoing_mass_edges=0.0;
    497                 for(useint=0; useint<3; useint++){
    498                     if(edgeflux_store[3*(3*k+useint)]< 0.){
    499                         //outgoing_mass_edges+=1.0;
    500                         outgoing_mass_edges+=(edgeflux_store[3*(3*k+useint)]);
    501                     }
    502                 }
    503                 outgoing_mass_edges*=local_timestep;
    504             }
    505 
    506             ki=3*k+i;   
    507             ki2=ki*2;
    508             ki3 = ki*3;
    509            
    510             // Prevent outflow from 'seriously' dry cells
    511             // Idea: The cell will not go dry if:
    512             // total_outgoing_flux <= cell volume = Area_triangle*hc
    513             vol=areas[k]*hc;
    514             if((edgeflux_store[ki3]< 0.0) && (-outgoing_mass_edges> vol)){
    515                
    516                 // This bound could be improved (e.g. we could actually sum the
    517                 // + and - fluxes and check if they are too large).  However,
    518                 // the advantage of this method is that we don't have to worry
    519                 // about subsequent changes to the + edgeflux caused by
    520                 // constraints associated with neighbouring triangles.
    521                 tmp = vol/(-(outgoing_mass_edges)) ;
    522                 if(tmp< 1.0){
    523                     edgeflux_store[ki3+0]*=tmp;
    524                     edgeflux_store[ki3+1]*=tmp;
    525                     edgeflux_store[ki3+2]*=tmp;
    526 
    527                     // Compute neighbour edge index
    528                     n = neighbours[ki];
    529                     if(n>=0){
    530                         nm = 3*n + neighbour_edges[ki];
    531                         nm3 = nm*3;
    532                         edgeflux_store[nm3+0]*=tmp;
    533                         edgeflux_store[nm3+1]*=tmp;
    534                         edgeflux_store[nm3+2]*=tmp;
    535                     }
    536                 }
    537             }
    538         }
    539     }
     481    //// GD HACK
     482    //// Limit edgefluxes, for mass conservation near wet/dry cells
     483    //for(k=0; k< number_of_elements; k++){
     484    //    //continue;
     485    //    hc = height_centroid_values[k];
     486    //    // Loop over every edge
     487    //    for(i = 0; i<3; i++){
     488    //        if(i==0){
     489    //            // Add up the outgoing flux through the cell -- only do this once (i==0)
     490    //            outgoing_mass_edges=0.0;
     491    //            for(useint=0; useint<3; useint++){
     492    //                if(edgeflux_store[3*(3*k+useint)]< 0.){
     493    //                    //outgoing_mass_edges+=1.0;
     494    //                    outgoing_mass_edges+=(edgeflux_store[3*(3*k+useint)]);
     495    //                }
     496    //            }
     497    //            outgoing_mass_edges*=local_timestep;
     498    //        }
     499
     500    //        ki=3*k+i;   
     501    //        ki2=ki*2;
     502    //        ki3 = ki*3;
     503    //       
     504    //        // Prevent outflow from 'seriously' dry cells
     505    //        // Idea: The cell will not go dry if:
     506    //        // total_outgoing_flux <= cell volume = Area_triangle*hc
     507    //        vol=areas[k]*hc;
     508    //        if((edgeflux_store[ki3]< 0.0) && (-outgoing_mass_edges> vol)){
     509    //           
     510    //            // This bound could be improved (e.g. we could actually sum the
     511    //            // + and - fluxes and check if they are too large).  However,
     512    //            // the advantage of this method is that we don't have to worry
     513    //            // about subsequent changes to the + edgeflux caused by
     514    //            // constraints associated with neighbouring triangles.
     515    //            tmp = vol/(-(outgoing_mass_edges)) ;
     516    //            if(tmp< 1.0){
     517    //                edgeflux_store[ki3+0]*=tmp;
     518    //                edgeflux_store[ki3+1]*=tmp;
     519    //                edgeflux_store[ki3+2]*=tmp;
     520
     521    //                // Compute neighbour edge index
     522    //                n = neighbours[ki];
     523    //                if(n>=0){
     524    //                    nm = 3*n + neighbour_edges[ki];
     525    //                    nm3 = nm*3;
     526    //                    edgeflux_store[nm3+0]*=tmp;
     527    //                    edgeflux_store[nm3+1]*=tmp;
     528    //                    edgeflux_store[nm3+2]*=tmp;
     529    //                }
     530    //            }
     531    //        }
     532    //    }
     533    // }
    540534
    541535    //printf("%e \n", edgeflux_store[3*30*3]);
    … …  
    592586
    593587    // Hack to ensure we only update the timestep on the first call within each rk2/rk3 step
    594     if((call-2)%timestep_order==0) timestep=local_timestep;
     588    if((call-2)%timestep_fluxcalls==0) timestep=local_timestep;
    595589
    596590    free(edgeflux_store);
    … …  
    851845  // condition) set its momentum to zero too. Aim is to prevent local 'pits'
    852846  // with a non-zero water surface gradient building up momentum over time
    853   //for (k=0; k<number_of_elements;k++){
    854 
    855   //    k3=k*3;
    856   //    k0 = surrogate_neighbours[k3];
    857   //    k1 = surrogate_neighbours[k3 + 1];
    858   //    k2 = surrogate_neighbours[k3 + 2];
    859 
    860   //    if((height_centroid_values[k0] < minimum_allowed_height | k0==k) &
    861   //       (height_centroid_values[k1] < minimum_allowed_height | k1==k) &
    862   //       (height_centroid_values[k2] < minimum_allowed_height | k2==k)){
    863   //            xmom_centroid_store[k] = 0.;
    864   //            xmom_centroid_values[k] = 0.;
    865   //            ymom_centroid_store[k] = 0.;
    866   //            ymom_centroid_values[k] = 0.;
    867 
    868   //    }
    869   //}
     847  for (k=0; k<number_of_elements;k++){
     848
     849      k3=k*3;
     850      k0 = surrogate_neighbours[k3];
     851      k1 = surrogate_neighbours[k3 + 1];
     852      k2 = surrogate_neighbours[k3 + 2];
     853
     854      if((height_centroid_values[k0] < minimum_allowed_height | k0==k) &
     855         (height_centroid_values[k1] < minimum_allowed_height | k1==k) &
     856         (height_centroid_values[k2] < minimum_allowed_height | k2==k)){
     857              xmom_centroid_store[k] = 0.;
     858              xmom_centroid_values[k] = 0.;
     859              ymom_centroid_store[k] = 0.;
     860              ymom_centroid_values[k] = 0.;
     861
     862      }
     863  }
    870864
    871865  // Begin extrapolation routine
    … …  
    10701064      c_tmp=1.0/(a_tmp-b_tmp);
    10711065      d_tmp= 1.0-(c_tmp*a_tmp);
     1066      hfactor=1.0;
    10721067      hfactor= max(0., min(c_tmp*max(hmin,0.0)/max(hc,1.0e-06)+d_tmp,
    10731068                           min(c_tmp*max(hc,0.)/max(hmax,1.0e-06)+d_tmp, 1.0))
    … …  
    10801075      // Set hfactor to 0 smoothly as hmin --> minumum_allowed_height
    10811076      hfactor=min( 1.2*max(hmin-minimum_allowed_height,0.)/(max(hmin,0.)+1.*minimum_allowed_height), hfactor);
     1077      //hfactor=min( 1.2*max(hc-minimum_allowed_height,0.)/(max(hc,0.)+1.*minimum_allowed_height), hfactor);
    10821078
    10831079      //-----------------------------------
    … …  
    11121108   
    11131109      beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor;
     1110
     1111      //if(hmin>minimum_allowed_height) beta_tmp=1.0;
    11141112      //beta_tmp = beta_w_dry*0. + (beta_w - beta_w_dry*0.) * hfactor;
    11151113      //beta_tmp=1.0;
    … …  
    16131611   
    16141612  double timestep, epsilon, H0, g;
    1615   int optimise_dry_cells, timestep_order;
     1613  int optimise_dry_cells, timestep_fluxcalls;
    16161614   
    16171615  // Convert Python arguments to C
    … …  
    16421640            &max_speed_array,
    16431641            &optimise_dry_cells,
    1644             &timestep_order,
     1642            &timestep_fluxcalls,
    16451643            &stage_centroid_values,
    16461644            &xmom_centroid_values,
    … …  
    17141712                     (double*) max_speed_array -> data,
    17151713                     optimise_dry_cells,
    1716                      timestep_order,
     1714                     timestep_fluxcalls,
    17171715                     (double*) stage_centroid_values -> data,
    17181716                     (double*) xmom_centroid_values -> data,

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

    r9033 r9039  
    2828domain.set_name('runup_v2')                         # Output to file runup.sww
    2929domain.set_datadir('.')                          # Use current folder
    30 domain.set_store_centroids(True)
     30#domain.set_store_centroids(True)
    3131#domain.set_flow_algorithm('DE1')
    3232#domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})

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

    r9038 r9039  
    2121from anuga.structures.inlet_operator import Inlet_operator
    2222from anuga.shallow_water.shallow_water_domain import Domain as Domain
    23 #from bal_and import *
     23from bal_and import *
    2424#from balanced_dev import Domain as Domain
    2525#from anuga_tsunami import *
    … …  
    3232domain = Domain(points, vertices, boundary) # Create domain
    3333domain.set_name('channel_SU_2_v2') # Output name
    34 domain.set_flow_algorithm('DE1')
     34#domain.set_flow_algorithm('DE1')
    3535#------------------------------------------------------------------------------
    3636# Setup initial conditions
Note: See TracChangeset for help on using the changeset viewer.