Ignore:
Timestamp:
Feb 5, 2010, 5:24:39 PM (13 years ago)
Author:
steve
Message:

Updating to relocated repository

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water_balanced/swb_domain.py

    r7575 r7616  
    8080        # set some defaults
    8181        #---------------------
    82         self.set_timestepping_method(1)
     82        self.set_timestepping_method(2)
    8383        self.set_default_order(2)
    8484        self.set_new_mannings_function(True)
    8585        self.set_centroid_transmissive_bc(True)
     86        self.set_CFL(1.0)
     87        self.set_beta(1.0)
    8688
    8789    ##
     
    114116        #Call correct module function (either from this module or C-extension)
    115117
     118        #compute_fluxes(self)
     119
     120        #return
    116121        from swb_domain_ext import compute_fluxes_c
    117122
     
    129134        self.flux_timestep = compute_fluxes_c(timestep, self, W, UH, VH, H, Z, U, V)
    130135
     136        #print self.flux_timestep
     137
    131138    ##
    132139    # @brief
     
    150157        """
    151158
    152 
     159        #Sww_domain.distribute_to_vertices_and_edges(self)
     160        #return
     161   
    153162        #Shortcuts
    154163        W  = self.quantities['stage']
     
    178187        num.putmask(uh_C, h_C < 1.0e-15, 0.0)
    179188        num.putmask(vh_C, h_C < 1.0e-15, 0.0)
    180         num.putmask(h_C, h_C < 1.0e-15, 1.0e-15)       
    181        
    182         u_C[:]  = uh_C/h_C
    183         v_C[:]  = vh_C/h_C
    184        
    185         for name in [ 'height', 'elevation', 'xvelocity', 'yvelocity' ]:
     189        num.putmask(h_C, h_C < 1.0e-15, 1.0e-16)       
     190
     191        H0 = 1.0e-8
     192        u_C[:]  = uh_C/(h_C + H0/h_C)
     193        v_C[:]  = vh_C/(h_C + H0/h_C)
     194
     195        num.putmask(h_C, h_C < 1.0e-15, 0.0)
     196       
     197        for name in [ 'stage', 'height', 'xvelocity', 'yvelocity' ]:
    186198            Q = self.quantities[name]
    187199            if self._order_ == 1:
     
    189201            elif self._order_ == 2:
    190202                Q.extrapolate_second_order_and_limit_by_edge()
     203                #Q.extrapolate_second_order_and_limit_by_vertex()
    191204            else:
    192205                raise 'Unknown order'
     
    202215
    203216
    204         w_E[:]   = z_E + h_E
    205 
    206         #num.putmask(u_E, h_temp <= 0.0, 0.0)
    207         #num.putmask(v_E, h_temp <= 0.0, 0.0)
    208         #num.putmask(w_E, h_temp <= 0.0, z_E+h_E)
     217        #minh_E = num.min(h_E)
     218        #msg = 'min h_E = %g ' % minh_E
     219        #assert minh_E >= -1.0e-15, msg
     220
     221        z_E[:]   = w_E - h_E
     222
     223        num.putmask(h_E, h_E <= 1.0e-15, 0.0)
     224        num.putmask(u_E, h_E <= 1.0e-15, 0.0)
     225        num.putmask(v_E, h_E <= 1.0e-15, 0.0)
     226        num.putmask(w_E, h_E <= 1.0e-15, z_E)
    209227        #num.putmask(h_E, h_E <= 0.0, 0.0)
    210228       
     
    246264
    247265
     266
     267
     268    ##
     269    # @brief
     270    def distribute_to_vertices_and_edges_h(self):
     271        """Distribution from centroids to edges specific to the SWW eqn.
     272
     273        It will ensure that h (w-z) is always non-negative even in the
     274        presence of steep bed-slopes by taking a weighted average between shallow
     275        and deep cases.
     276
     277        In addition, all conserved quantities get distributed as per either a
     278        constant (order==1) or a piecewise linear function (order==2).
     279       
     280
     281        Precondition:
     282        All conserved quantities defined at centroids and bed elevation defined at
     283        edges.
     284       
     285        Postcondition
     286        Evolved quantities defined at vertices and edges
     287        """
     288
     289
     290        #Shortcuts
     291        W  = self.quantities['stage']
     292        UH = self.quantities['xmomentum']
     293        VH = self.quantities['ymomentum']
     294        H  = self.quantities['height']
     295        Z  = self.quantities['elevation']
     296        U  = self.quantities['xvelocity']
     297        V  = self.quantities['yvelocity']
     298
     299        #Arrays   
     300        w_C   = W.centroid_values   
     301        uh_C  = UH.centroid_values
     302        vh_C  = VH.centroid_values   
     303        z_C   = Z.centroid_values
     304        h_C   = H.centroid_values
     305        u_C   = U.centroid_values
     306        v_C   = V.centroid_values
     307
     308        w_C[:] = num.maximum(w_C, z_C)
     309       
     310        h_C[:] = w_C - z_C
     311
     312
     313        assert num.min(h_C) >= 0
     314               
     315        num.putmask(uh_C, h_C < 1.0e-15, 0.0)
     316        num.putmask(vh_C, h_C < 1.0e-15, 0.0)
     317        num.putmask(h_C, h_C < 1.0e-15, 1.0e-16)       
     318       
     319        u_C[:]  = uh_C/h_C
     320        v_C[:]  = vh_C/h_C
     321
     322        num.putmask(h_C, h_C < 1.0e-15, 0.0)
     323       
     324        for name in [ 'stage', 'height', 'xvelocity', 'yvelocity' ]:
     325            Q = self.quantities[name]
     326            if self._order_ == 1:
     327                Q.extrapolate_first_order()
     328            elif self._order_ == 2:
     329                Q.extrapolate_second_order_and_limit_by_edge()
     330                #Q.extrapolate_second_order_and_limit_by_vertex()
     331            else:
     332                raise 'Unknown order'
     333
     334
     335        w_E     = W.edge_values
     336        uh_E    = UH.edge_values
     337        vh_E    = VH.edge_values       
     338        h_E     = H.edge_values
     339        z_E     = Z.edge_values
     340        u_E     = U.edge_values
     341        v_E     = V.edge_values         
     342
     343
     344        #minh_E = num.min(h_E)
     345        #msg = 'min h_E = %g ' % minh_E
     346        #assert minh_E >= -1.0e-15, msg
     347
     348        z_E[:]   = w_E - h_E
     349
     350        num.putmask(h_E, h_E <= 1.0e-8, 0.0)
     351        num.putmask(u_E, h_E <= 1.0e-8, 0.0)
     352        num.putmask(v_E, h_E <= 1.0e-8, 0.0)
     353        num.putmask(w_E, h_E <= 1.0e-8, z_E)
     354        #num.putmask(h_E, h_E <= 0.0, 0.0)
     355       
     356        uh_E[:] = u_E * h_E
     357        vh_E[:] = v_E * h_E
     358
     359        """
     360        print '=========================================================='
     361        print 'Time ', self.get_time()
     362        print h_E
     363        print uh_E
     364        print vh_E
     365        """
     366       
     367        # Compute vertex values by interpolation
     368        for name in self.evolved_quantities:
     369            Q = self.quantities[name]
     370            Q.interpolate_from_edges_to_vertices()
     371
     372
     373        w_V     = W.vertex_values
     374        uh_V    = UH.vertex_values
     375        vh_V    = VH.vertex_values     
     376        z_V     = Z.vertex_values       
     377        h_V     = H.vertex_values
     378        u_V     = U.vertex_values
     379        v_V     = V.vertex_values               
     380
     381
     382        #w_V[:]    = z_V + h_V
     383
     384        #num.putmask(u_V, h_V <= 0.0, 0.0)
     385        #num.putmask(v_V, h_V <= 0.0, 0.0)
     386        #num.putmask(w_V, h_V <= 0.0, z_V)       
     387        #num.putmask(h_V, h_V <= 0.0, 0.0)
     388       
     389        uh_V[:] = u_V * h_V
     390        vh_V[:] = v_V * h_V
     391
     392
     393
     394
    248395    ##
    249396    # @brief Code to let us use old shallow water domain BCs
Note: See TracChangeset for help on using the changeset viewer.