Changeset 7616


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

Updating to relocated repository

Location:
anuga_core/source/anuga
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

    r7573 r7616  
    362362   
    363363    sign = 0.0;
    364     if (qmin > 0) {
     364    if (qmin > 0.0) {
    365365      sign = 1.0;
    366366    } else if (qmax < 0) {
     
    373373      dqa[i] = dq;                      //Save dq for use in updating vertex values 
    374374     
    375       // FIXME SR 20091125 This caused problems in shallow_water_balanced
    376       // commenting out problem
     375
    377376      // Just limit non boundary edges so that we can reconstruct a linear function
    378       if (neighbours[k3+i] >= 0) {
     377      // FIXME Problem with stability on edges
     378      //if (neighbours[k3+i] >= 0) {
    379379        r = 1.0;
    380380     
     
    383383           
    384384        phi = min( min(r*beta, 1.0), phi);
    385         }
    386 
    387       if (neighbours[k3+i] < 0) {
    388         r = 1.0;
    389      
    390         if (dq > 0.0 && sign == -1.0 ) r = (0.0 - qc)/dq;
    391         if (dq < 0.0 && sign ==  1.0 ) r = (0.0 - qc)/dq;     
     385        //      }
     386
     387      //
     388      /* if (neighbours[k3+i] < 0) { */
     389      /*        r = 1.0; */
     390     
     391      /*        if (dq > 0.0 && (sign == -1.0 || sign == 0.0 )) r = (0.0 - qc)/dq; */
     392      /*        if (dq < 0.0 && (sign ==  1.0 || sign == 0.0 )) r = (0.0 - qc)/dq; */
    392393           
    393         phi = min( min(r*beta, 1.0), phi);
    394         }
     394      /*        phi = min( min(r*beta, 1.0), phi); */
     395      /*        } */
    395396   
    396397    }
     
    11941195  if (!PyArg_ParseTuple(args, "O",&quantity)) {
    11951196      PyErr_SetString(PyExc_RuntimeError,
    1196                       "quantity_ext.c: extrapolate_second_order_and_limit could not parse input");     
     1197                      "quantity_ext.c: extrapolate_second_order_and_limit_by_edge could not parse input");     
    11971198      return NULL;
    11981199  }
     
    12011202  if (!domain) {
    12021203      PyErr_SetString(PyExc_RuntimeError,
    1203                       "quantity_ext.c: extrapolate_second_order_and_limit could not obtain domain object from quantity");       
     1204                      "quantity_ext.c: extrapolate_second_order_and_limit_by_edge could not obtain domain object from quantity");       
    12041205      return NULL;
    12051206  }
  • anuga_core/source/anuga/compile_all.py

    r6718 r7616  
    2121
    2222os.chdir('..')
     23os.chdir('shallow_water_balanced')
     24execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
     25
     26os.chdir('..')
    2327os.chdir('mesh_engine')
    2428execfile('..' + os.sep + 'utilities' + os.sep + 'compile.py')
  • 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
  • anuga_core/source/anuga/shallow_water_balanced/test_swb_balance.py

    r7573 r7616  
    6363
    6464        # Limit
    65         domain.tight_slope_limiters = 0
    6665        domain.distribute_to_vertices_and_edges()
    6766
     
    7271
    7372        # Now try with a non-flat bed - closely hugging initial stage in places
    74         # This will create alphas in the range [0, 0.478260, 1]
    7573        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
    7674        domain.set_quantity('elevation', [[0,0,0],
     
    7977                                          [0,2,4]])
    8078        stage = domain.quantities['stage']
     79        elevation = domain.quantities['elevation']
     80        height = domain.quantities['height']
    8181
    8282        ref_centroid_values = copy.copy(stage.centroid_values[:])    # Copy
     
    8484
    8585        # Limit
    86         domain.tight_slope_limiters = 0
    8786        domain.distribute_to_vertices_and_edges()
    8887
     
    9897        # Check actual results
    9998        assert num.allclose(stage.vertex_values,
    100                             [[ 2.,          2.,          2.        ],
    101                              [ 1.93333333,  2.03333333,  6.03333333],
    102                              [ 8.4,         3.8,         3.8       ],
    103                              [ 3.33333333,  5.33333333,  7.33333333]]) or \
    104                num.allclose(stage.vertex_values,
    105                             [[ 1.06666667,  3.86666667,  1.06666667],
    106                              [ 1.93333333,  2.03333333,  6.03333333],
    107                              [ 5.46666667,  3.06666667,  7.46666667],
    108                              [ 6.53333333,  4.69333333,  4.77333333]])
     99                            [[ 2.66666667,  0.66666667,  2.66666667],
     100                             [ 3.33333333,  3.33333333,  3.33333333],
     101                             [ 3.73333333,  4.93333333,  7.33333333],
     102                             [ 7.33333333,  4.93333333,  3.73333333]])
    109103
    110104
  • anuga_core/source/anuga/shallow_water_balanced/test_swb_boundary_condition.py

    r7573 r7616  
    871871
    872872if __name__ == "__main__":
    873     suite = unittest.makeSuite(Test_swb_boundary_condition, 'test')
     873    suite = unittest.makeSuite(Test_swb_boundary_condition, 'test_boundary_condition_time')
    874874    runner = unittest.TextTestRunner(verbosity=1)
    875875    runner.run(suite)
Note: See TracChangeset for help on using the changeset viewer.