Changeset 9257


Ignore:
Timestamp:
Jul 10, 2014, 5:42:24 PM (11 years ago)
Author:
davies
Message:

Cutting a boundary_flux_integral operator into the code

Location:
trunk
Files:
2 added
3 edited

Legend:

Unmodified
Added
Removed

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

    r9256 r9257  
    250250        self.riverwallData=anuga.structures.riverwall.RiverWall(self)
    251251
     252        # Keep track of the fluxes through the boundaries
     253        # Only works for DE algorithms at present
     254        max_time_substeps=3 # Maximum number of substeps supported by any timestepping method
     255        # boundary_flux_sum holds boundary fluxes on each sub-step
     256        self.boundary_flux_sum=num.array([0.]*max_time_substeps)
     257        from anuga.operators.boundary_flux_integral_operator import boundary_flux_integral_operator
     258        self.boundary_flux_integral=boundary_flux_integral_operator(self)
     259        # Make an integer counting how many times we call compute_fluxes_central -- so we know which substep we are on
     260        self.call=1
    252261
    253262
    … …  
    480489        self.maximum_allowed_speed=0.0
    481490
    482         ## FIXME: Should implement tracking of boundary fluxes
    483         ## Keep track of the fluxes through the boundaries
    484         self.boundary_flux_integral=num.ndarray(1)
    485         self.boundary_flux_integral[0]=0.
    486         self.boundary_flux_sum=num.ndarray(1)
    487         self.boundary_flux_sum[0]=0.
    488 
    489         self.call=1 # Integer counting how many times we call compute_fluxes_central
    490491
    491492        if self.processor == 0 and self.verbose:
    … …  
    547548        self.maximum_allowed_speed=0.0
    548549
    549         ## FIXME: Should implement tracking of boundary fluxes
    550         ## Keep track of the fluxes through the boundaries
    551         self.boundary_flux_integral=num.ndarray(1)
    552         self.boundary_flux_integral[0]=0.
    553         self.boundary_flux_sum=num.ndarray(1)
    554         self.boundary_flux_sum[0]=0.
    555 
    556         self.call=1 # Integer counting how many times we call compute_fluxes_central
    557 
    558550        if self.processor == 0 and self.verbose:
    559551            print '##########################################################################'
    … …  
    615607        self.maximum_allowed_speed=0.0
    616608
    617         ## FIXME: Should implement tracking of boundary fluxes
    618         ## Keep track of the fluxes through the boundaries
    619         self.boundary_flux_integral=num.ndarray(1)
    620         self.boundary_flux_integral[0]=0.
    621         self.boundary_flux_sum=num.ndarray(1)
    622         self.boundary_flux_sum[0]=0.
    623 
    624         self.call=1 # Integer counting how many times we call compute_fluxes_central
    625 
    626609        if self.processor == 0 and self.verbose:
    627610            print '##########################################################################'
    … …  
    682665
    683666        self.maximum_allowed_speed=0.0
    684 
    685         ## FIXME: Should implement tracking of boundary fluxes
    686         ## Keep track of the fluxes through the boundaries
    687         self.boundary_flux_integral=num.ndarray(1)
    688         self.boundary_flux_integral[0]=0.
    689         self.boundary_flux_sum=num.ndarray(1)
    690         self.boundary_flux_sum[0]=0.
    691 
    692         self.call=1 # Integer counting how many times we call compute_fluxes_central
    693667
    694668        if self.processor == 0 and self.verbose:
    … …  
    13371311        return water_volume
    13381312
     1313    def get_boundary_flux_integral(self):
     1314        """
     1315            Compute the boundary flux integral.
     1316            Should work in parallel
     1317        """
     1318   
     1319        from anuga import myid, numprocs, send, receive, barrier
     1320
     1321        if not self.compute_fluxes_method=='DE':
     1322            msg='Boundary flux integral only supported for DE fluxes '+\
     1323                '(because computation of boundary_flux_sum is only implemented there)'
     1324            raise Exception, msg
     1325           
     1326        flux_integral = self.boundary_flux_integral.boundary_flux_integral
     1327           
     1328        if myid == 0:
     1329            for i in range(1,numprocs):
     1330                remote_flux_integral = receive(i)
     1331                flux_integral = flux_integral + remote_flux_integral
     1332        else:
     1333            send(flux_integral,0)
     1334       
     1335        barrier()
     1336   
     1337        if myid == 0:
     1338            for i in range(1,numprocs):
     1339                send(flux_integral,i)
     1340        else:
     1341            flux_integral = receive(0)
     1342       
     1343        return flux_integral
    13391344
    13401345    def get_flow_through_cross_section(self, polyline, verbose=False):
    … …  
    15871592                                               self.riverwallData.hydraulic_properties)
    15881593
    1589             ## FIXME: This won't work in parallel since the timestep has not been updated to include other processors.
    1590             ## Update the boundary flux integral
    1591             ## FIXME GD: This should be moved to its own routine -- something like this
    1592             ##        could be included in all solvers.
    1593             #if(self.timestepping_method=='rk2'):
    1594             #    if(self.call%2==1):
    1595             #      self.boundary_flux_integral[0]= self.boundary_flux_integral[0] +\
    1596             #                                        self.boundary_flux_sum[0]*self.timestep*0.5
    1597             #      #print 'dbfi ', self.boundary_flux_integral, self.boundary_flux_sum
    1598             #      self.boundary_flux_sum[0]=0.
    1599             #elif(self.timestepping_method=='euler'):
    1600             #    # FIXME GD: Unsupported at present.
    1601             #    self.boundary_flux_integral=0.
    1602             #    # This presently doesn't work -- this section of code may need to go elsewhere
    1603             #    #self.boundary_flux_integral[0]= self.boundary_flux_integral[0] +\
    1604             #    #                                  self.boundary_flux_sum[0]*self.timestep
    1605             #    #self.boundary_flux_sum[0]=0.
    1606             #elif(self.timestepping_method=='rk3'):
    1607             #    # FIXME GD: Unsupported at present.
    1608             #    self.boundary_flux_integral=0.
    1609             #    # FIXME GD: This needs to be implemented
    1610             #else:
    1611             #    mess='ERROR: domain.timestepping_method', self.timestepping_method,' method not recognised'
    1612             #    raise Exception, mess
    1613 
    16141594            self.flux_timestep = flux_timestep
    16151595
    … …  
    23382318        These calculations are only approximate since they don't use the
    23392319        flux calculation used in evolve
     2320
     2321        See get_boundary_flux_integral for an exact computation
    23402322        """
    23412323               

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

    r9232 r9257  
    415415      // Smoothing by stage alone can cause high velocities / slow draining for nearly dry cells
    416416      if(i==0) edgeflux[i] += (s_max*s_min)*(max(q_right_rotated[i],ze) - max(q_left_rotated[i],ze));
     417      //if(i==0) edgeflux[i] += (s_max*s_min)*(h_right - h_left);
    417418      if(i==1) edgeflux[i] += (s_max*s_min)*(uh_right - uh_left);
    418419      if(i==2) edgeflux[i] += (s_max*s_min)*(vh_right - vh_left);
    … …  
    584585    double bedslope_work, local_timestep;
    585586    int neighbours_wet[3];//Work array
    586     int RiverWall_count;
     587    int RiverWall_count, substep_count;
    587588    double hle, hre, zc, zc_n, Qfactor, s1, s2, h1, h2;
    588589    double stage_edge_lim, outgoing_mass_edges, bedtop, bedbot, pressure_flux, hc, hc_n, tmp, tmp2;
    … …  
    606607    local_timestep=timestep;
    607608    RiverWall_count=0;
     609    substep_count=(call-2)%timestep_fluxcalls;
    608610
    609611
    … …  
    790792            edgeflux_store[ki3 + 2 ] = -edgeflux[2];
    791793
    792             // bedslope_work contains all gravity related terms -- weighting of
     794            // bedslope_work contains all gravity related terms
    793795            bedslope_work=length*(-g*0.5*(h_left*h_left - hle*hle -(hle+hc)*(zl-zc))+pressure_flux);
    794796
    … …  
    836838        } // End edge i (and neighbour n)
    837839        // Keep track of maximal speeds
    838         if((call-2)%timestep_fluxcalls==0) max_speed_array[k] = speed_max_last; //max_speed;
     840        if(substep_count==0) max_speed_array[k] = speed_max_last; //max_speed;
    839841
    840842
    … …  
    843845    //// GD HACK
    844846    //// Limit edgefluxes, for mass conservation near wet/dry cells
     847    //// This doesn't seem to be needed anymore
    845848    //for(k=0; k< number_of_elements; k++){
    846849    //    //continue;
    … …  
    925928            // boundary_flux_integral
    926929            if( (n<0 & tri_full_flag[k]==1) | ( n>=0 && (tri_full_flag[k]==1 & tri_full_flag[n]==0)) ){
    927                 boundary_flux_sum[0] += edgeflux_store[ki3];
     930                // boundary_flux_sum is an array with length = timestep_fluxcalls
     931                // For each sub-step, we put the boundary flux sum in.
     932                boundary_flux_sum[substep_count] += edgeflux_store[ki3];
    928933            }
    929 
    930934            // GD HACK
    931935            // Compute bed slope term
    … …  
    950954
    951955    // Hack to ensure we only update the timestep on the first call within each rk2/rk3 step
    952     if((call-2)%timestep_fluxcalls==0) timestep=local_timestep;
    953 
     956    if(substep_count==0) timestep=local_timestep;
     957           
    954958    free(edgeflux_store);
    955959    free(pressuregrad_store);

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

    r9200 r9257  
    99
    1010from math import sin, pi, exp
    11 #from anuga.shallow_water.shallow_water_domain import Domain as Domain
    12 #from anuga.shallow_water_balanced2.swb2_domain import Domain as Domain
    13 #path.append('/home/gareth/storage/anuga_clean/anuga_jan12/trunk/anuga_work/development/gareth/balanced_basic')
    14 #from swb2_domain import *
    15 #from balanced_basic import *
    16 #from balanced_dev import *
    1711from anuga import Domain as Domain
    18 #from bal_and import *
    19 #from anuga_tsunami import *
     12from anuga.operators.boundary_flux_integral_operator import boundary_flux_integral_operator
    2013
    2114#---------
    … …  
    2720domain.set_name('runup_sinusoid_v2')                         # Output to file runup.sww
    2821#domain.set_store_centroids(True)
    29 #domain.set_timestepping_method('euler')
    30 #domain.set_flow_algorithm('DE1')
    31 #domain.set_flow_algorithm('tsunami')
    3222domain.set_flow_algorithm('DE0')
    3323domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2,
    3424         'ymomentum': 2, 'elevation': 2})
    35 #domain.minimum_storable_height=1.0e-03
    3625domain.set_store_vertices_uniquely()
    3726#------------------
    … …  
    7059
    7160print '#-#', domain.minimum_allowed_height
    72 #domain.minimum_allowed_height=0.001
    73 #def frict_change(x,y):
    74 #       return 0.2*(x>0.5)+0.1*(x<=0.5)
    75 #
    76 #domain.set_quantity('friction',frict_change)
    7761
    7862domain.set_quantity('stage', stagefun)              # Constant negative initial stage
    … …  
    10892#Evolve the system through time
    10993#------------------------------
    110 
    11194for t in domain.evolve(yieldstep=tstep,finaltime=lasttime):
    11295#for t in domain.evolve(yieldstep=0.2,finaltime=60.):
    … …  
    121104    vv = vv*(dd>1.0e-03)
    122105    print 'Peak velocity is: ', vv.max(), vv.argmax(), dd[vv.argmax()]
    123     print 'Volume is', sum(dd_raw*domain.areas)   
     106    vol=sum(dd_raw*domain.areas)
     107    print 'Volume is', vol   
     108    flux_integral=domain.get_boundary_flux_integral()
     109    print 'Flux integral is: ', flux_integral
     110    print 'Their difference is:', vol-flux_integral
    124111    #print 'Volume less flux int', sum(dd_raw*domain.areas) - domain.boundary_flux_integral
    125 
     112   
    126113
    127114vel_final_inds=(vv>1.0e-01).nonzero()
Note: See TracChangeset for help on using the changeset viewer.