Changeset 7519


Ignore:
Timestamp:
Sep 20, 2009, 6:35:14 PM (14 years ago)
Author:
steve
Message:

Have been playing with using a slope limited velocity to calculate
fluxes (hence the addition of evolved_quantities as well as conserved
quantities.

But the commit is to fix a problem Rudy found with sww2dem. Seems
numpy.array2string is a little too clever, in that it summarizes
output if there is a long sequence of zeros to
[0.0, 0.0, 0.0, ... 0.0, 0.0 ,0.0] To get around this I have added
a call to numpy.set_options(threshold=sys.max_int) to turn this
behaviour off!

Location:
anuga_core/source
Files:
16 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

    r7498 r7519  
    6565                       boundary=None,
    6666                       conserved_quantities=None,
     67                       evolved_quantities=None,
    6768                       other_quantities=None,
    6869                       tagged_elements=None,
    … …  
    9091          conserved_quantities: List of quantity names entering the
    9192                                conservation equations
     93          evolved_quantities:   List of all quantities that evolve
    9294          other_quantities:     List of other quantity names
    9395
    … …  
    154156            self.conserved_quantities = conserved_quantities
    155157
     158        if evolved_quantities is None:
     159            self.evolved_quantities = self.conserved_quantities
     160        else:
     161            self.evolved_quantities = evolved_quantities
     162           
    156163        # List of other quantity names
    157164        if other_quantities is None:
    … …  
    160167            self.other_quantities = other_quantities
    161168
     169        # Test that conserved_quantities are stored in the first entries of evolved_quantities
     170        for i, quantity in enumerate(self.conserved_quantities):
     171            msg = 'The conserved quantities must be the first entries of evolved_quantities'
     172            assert quantity == self.evolved_quantities[i], msg
     173           
     174
    162175        # Build dictionary of Quantity instances keyed by quantity names
    163176        self.quantities = {}
    164177
    165         for name in self.conserved_quantities:
     178        for name in self.evolved_quantities:
    166179            self.quantities[name] = Quantity(self)
    167180        for name in self.other_quantities:
    … …  
    325338        return self.mesh.get_number_of_nodes(*args, **kwargs)
    326339
     340    def get_number_of_triangles(self, *args, **kwargs):
     341        return self.mesh.get_number_of_triangles(*args, **kwargs)   
     342
    327343    def get_normal(self, *args, **kwargs):
    328344        return self.mesh.get_normal(*args, **kwargs)
    … …  
    415431
    416432        return q
     433
     434    ##
     435    # @brief Get evolved quantities for a volume.
     436    # @param vol_id ID of the volume we want the conserved quantities for.
     437    # @param vertex If specified, use as index for edge values.
     438    # @param edge If specified, use as index for edge values.
     439    # @return Vector of conserved quantities.
     440    # @note If neither 'vertex' or 'edge' specified, use centroid values.
     441    # @note If both 'vertex' and 'edge' specified, raise exception.
     442    def get_evolved_quantities(self, vol_id,
     443                               vertex=None,
     444                               edge=None):
     445        """Get evolved quantities at volume vol_id.
     446
     447        If vertex is specified use it as index for vertex values
     448        If edge is specified use it as index for edge values
     449        If neither are specified use centroid values
     450        If both are specified an exeception is raised
     451
     452        Return value: Vector of length == number_of_conserved quantities
     453        """
     454
     455        if not (vertex is None or edge is None):
     456            msg = 'Values for both vertex and edge was specified.'
     457            msg += 'Only one (or none) is allowed.'
     458            raise Exception, msg
     459
     460        q = num.zeros(len(self.evolved_quantities), num.float)
     461
     462        for i, name in enumerate(self.evolved_quantities):
     463            Q = self.quantities[name]
     464            if vertex is not None:
     465                q[i] = Q.vertex_values[vol_id, vertex]
     466            elif edge is not None:
     467                q[i] = Q.edge_values[vol_id, edge]
     468            else:
     469                q[i] = Q.centroid_values[vol_id]
     470
     471        return q
     472
     473  ##
     474    # @brief
     475    # @param flag
     476    def set_CFL(self, cfl=1.0):
     477        """Set CFL parameter, warn if greater than 1.0
     478        """
     479        if cfl > 1.0:
     480            self.CFL = cfl
     481            log.warn('Setting CFL > 1.0')
     482
     483        assert cfl > 0.0
     484        self.CFL = cfl
     485
     486
    417487
    418488    ##
    … …  
    876946            assert quantity in self.quantities, msg
    877947
    878         ##assert hasattr(self, 'boundary_objects')
     948
     949        for i, quantity in enumerate(self.conserved_quantities):
     950            msg = 'Conserved quantities must be the first entries of evolved_quantities'
     951            assert quantity == self.evolved_quantities[i], msg
     952 
    879953
    880954    ##
    … …  
    16541728
    16551729    ##
    1656     # @brief Backup conserved quantities.
     1730    # @brief Backup conserved quantities 
    16571731    def backup_conserved_quantities(self):
    1658         N = len(self) # Number_of_triangles
    16591732
    16601733        # Backup conserved_quantities centroid values
    … …  
    16641737
    16651738    ##
    1666     # @brief ??
    1667     # @param a ??
    1668     # @param b ??
     1739    # @brief Combines current C and saved centroid values S as C = aC + bS
     1740    # @param a factor in combination
     1741    # @param b factor in combination
    16691742    def saxpy_conserved_quantities(self, a, b):
    1670         N = len(self) #number_of_triangles
    16711743
    16721744        # Backup conserved_quantities centroid values
    … …  
    16751747            Q.saxpy_centroid_values(a, b)
    16761748
     1749           
     1750
     1751
     1752    ##
     1753    # @brief Mapping between conserved quantites and evolved quantities
     1754    # @param q_cons array of conserved quantity values
     1755    # @param q_evolved array of current evolved quantity values
     1756    def  conserved_to_evolved(self, q_cons, q_evolved):
     1757        """Needs to be overridden by Domain subclass
     1758        """
     1759
     1760        if len(q_cons) == len(q_evolved):
     1761            q_evolved[:] = q_cons
     1762        else:
     1763            msg = 'Method conserved_to_evolved must be overridden by Domain subclass'
     1764            raise Exception, msg
     1765   
    16771766    ##
    16781767    # @brief Update boundary values for all conserved quantities.
    … …  
    16921781                log.critical('WARNING: Ignored boundary segment (None)')
    16931782            else:
    1694                 q = B.evaluate(vol_id, edge_id)
    1695 
    1696                 for j, name in enumerate(self.conserved_quantities):
     1783                q_cons = B.evaluate(vol_id, edge_id)
     1784
     1785                if len(q_cons) == len(self.evolved_quantities):
     1786                    # conserved and evolved quantities are the same
     1787                    q_evol = q_cons
     1788                elif len(q_cons) == len(self.conserved_quantities):
     1789                    # boundary just returns conserved quantities
     1790                    # Need to calculate all the evolved quantities
     1791                    # Use default conversion
     1792
     1793                    q_evol = self.get_evolved_quantities(vol_id, edge = edge_id)
     1794
     1795                    self.conserved_to_evolved(q_cons, q_evol)
     1796                else:
     1797                    msg = 'Boundary must return array of either conserved or evolved quantities'
     1798                    raise Exception, msg
     1799               
     1800                for j, name in enumerate(self.evolved_quantities):
    16971801                    Q = self.quantities[name]
    1698                     Q.boundary_values[i] = q[j]
     1802                    Q.boundary_values[i] = q_evol[j]
    16991803
    17001804    ##

  • anuga_core/source/anuga/abstract_2d_finite_volumes/general_mesh.py

    r7498 r7519  
    235235               
    236236
     237    def get_number_of_triangles(self):
     238        return self.number_of_triangles
     239
     240   
    237241    def get_number_of_nodes(self):
    238242        return self.number_of_nodes
     243
    239244
    240245    def get_nodes(self, absolute=False):

  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

    r7276 r7519  
    712712
    713713
     714        // Explicit updates
     715        for (k=0; k<N; k++) {
     716                centroid_values[k] += timestep*explicit_update[k];
     717        }
     718
     719
     720
    714721        // Semi implicit updates
    715722        for (k=0; k<N; k++) {
    716723                denominator = 1.0 - timestep*semi_implicit_update[k];
    717                 if (denominator == 0.0) {
     724                if (denominator <= 0.0) {
    718725                        return -1;
    719726                } else {
    … …  
    724731
    725732
    726         // Explicit updates
    727         for (k=0; k<N; k++) {
    728                 centroid_values[k] += timestep*explicit_update[k];
    729         }
     733
    730734
    731735

  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_domain.py

    r7317 r7519  
    5353
    5454        conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
     55        evolved_quantities = ['stage', 'xmomentum', 'ymomentum', 'xvelocity']
     56       
    5557        other_quantities = ['elevation', 'friction']
    5658
    5759        domain = Domain(points, vertices, None,
    58                         conserved_quantities, other_quantities)
     60                        conserved_quantities, evolved_quantities, other_quantities)
    5961        domain.check_integrity()
    6062
    … …  
    6567        assert num.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.)
    6668
     69
     70
     71    def test_CFL(self):
     72        a = [0.0, 0.0]
     73        b = [0.0, 2.0]
     74        c = [2.0,0.0]
     75        d = [0.0, 4.0]
     76        e = [2.0, 2.0]
     77        f = [4.0,0.0]
     78
     79        points = [a, b, c, d, e, f]
     80        #bac, bce, ecf, dbe, daf, dae
     81        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     82
     83        conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
     84        evolved_quantities = ['stage', 'xmomentum', 'ymomentum', 'xvelocity']
     85       
     86        other_quantities = ['elevation', 'friction']
     87
     88        domain = Domain(points, vertices, None,
     89                        conserved_quantities, evolved_quantities, other_quantities)
     90
     91        try:
     92            domain.set_CFL(-0.1)
     93        except:
     94            pass
     95        else:
     96            msg = 'Should have caught a negative cfl'
     97            raise Exception, msg
     98
     99
     100       
     101        try:
     102            domain.set_CFL(2.0)
     103        except:
     104            pass
     105        else:
     106            msg = 'Should have warned of cfl>1.0'
     107            raise Exception, msg
     108
     109        assert domain.CFL == 2.0
     110       
    67111
    68112    def test_conserved_quantities(self):
    … …  
    511555
    512556
     557    def test_conserved_evolved_boundary_conditions(self):
     558
     559        a = [0.0, 0.0]
     560        b = [0.0, 2.0]
     561        c = [2.0,0.0]
     562        d = [0.0, 4.0]
     563        e = [2.0, 2.0]
     564        f = [4.0,0.0]
     565
     566        points = [a, b, c, d, e, f]
     567        #bac, bce, ecf, dbe
     568        vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
     569        boundary = { (0, 0): 'First',
     570                     (0, 2): 'First',
     571                     (2, 0): 'Second',
     572                     (2, 1): 'Second',
     573                     (3, 1): 'Second',
     574                     (3, 2): 'Second'}
     575
     576
     577 
     578        try:
     579            domain = Domain(points, vertices, boundary,
     580                            conserved_quantities = ['stage', 'xmomentum', 'ymomentum'],
     581                            evolved_quantities =\
     582                                   ['stage', 'xmomentum', 'xvelocity', 'ymomentum', 'yvelocity'])
     583        except:
     584            pass
     585        else:
     586            msg = 'Should have caught the evolved quantities not being in order'
     587            raise Exception, msg           
     588
     589
     590        domain = Domain(points, vertices, boundary,
     591                        conserved_quantities = ['stage', 'xmomentum', 'ymomentum'],
     592                        evolved_quantities =\
     593                        ['stage', 'xmomentum', 'ymomentum', 'xvelocity', 'yvelocity'])
     594
     595
     596        domain.set_quantity('stage', [[1,2,3], [5,5,5],
     597                                      [0,0,9], [6, -3, 3]])
     598
     599
     600        domain.set_boundary( {'First': Dirichlet_boundary([5,2,1,4,6]),
     601                              'Second': Transmissive_boundary(domain)} )
     602
     603        try:
     604            domain.update_boundary()
     605        except:
     606            pass
     607        else:
     608            msg = 'Should have caught the lack of conserved_to_evolved member function'
     609            raise Exception, msg
     610
     611
     612        def  conserved_to_evolved(q_cons, q_evol):
     613
     614            q_evol[0:3] = q_cons
     615            q_evol[3] = q_cons[1]/q_cons[0]
     616            q_evol[4] = q_cons[2]/q_cons[0]
     617
     618        domain.conserved_to_evolved = conserved_to_evolved
     619
     620        domain.update_boundary()
     621
     622
     623        assert domain.quantities['stage'].boundary_values[0] == 5. #Dirichlet
     624        assert domain.quantities['stage'].boundary_values[1] == 5. #Dirichlet
     625        assert domain.quantities['xvelocity'].boundary_values[0] == 4. #Dirichlet
     626        assert domain.quantities['yvelocity'].boundary_values[1] == 6. #Dirichlet
     627
     628        q_cons = domain.get_conserved_quantities(2, edge=0) #Transmissive
     629        assert domain.quantities['stage'    ].boundary_values[2] == q_cons[0]
     630        assert domain.quantities['xmomentum'].boundary_values[2] == q_cons[1]
     631        assert domain.quantities['ymomentum'].boundary_values[2] == q_cons[2]
     632        assert domain.quantities['xvelocity'].boundary_values[2] == q_cons[1]/q_cons[0]
     633        assert domain.quantities['yvelocity'].boundary_values[2] == q_cons[2]/q_cons[0]
     634
     635        q_cons = domain.get_conserved_quantities(2, edge=1) #Transmissive
     636        assert domain.quantities['stage'    ].boundary_values[3] == q_cons[0]
     637        assert domain.quantities['xmomentum'].boundary_values[3] == q_cons[1]
     638        assert domain.quantities['ymomentum'].boundary_values[3] == q_cons[2]
     639        assert domain.quantities['xvelocity'].boundary_values[3] == q_cons[1]/q_cons[0]
     640        assert domain.quantities['yvelocity'].boundary_values[3] == q_cons[2]/q_cons[0]       
     641
     642
     643        q_cons = domain.get_conserved_quantities(3, edge=1) #Transmissive
     644        assert domain.quantities['stage'    ].boundary_values[4] == q_cons[0]
     645        assert domain.quantities['xmomentum'].boundary_values[4] == q_cons[1]
     646        assert domain.quantities['ymomentum'].boundary_values[4] == q_cons[2]
     647        assert domain.quantities['xvelocity'].boundary_values[4] == q_cons[1]/q_cons[0]
     648        assert domain.quantities['yvelocity'].boundary_values[4] == q_cons[2]/q_cons[0]               
     649
     650
     651        q_cons = domain.get_conserved_quantities(3, edge=2) #Transmissive
     652        assert domain.quantities['stage'    ].boundary_values[5] == q_cons[0]
     653        assert domain.quantities['xmomentum'].boundary_values[5] == q_cons[1]
     654        assert domain.quantities['ymomentum'].boundary_values[5] == q_cons[2]
     655        assert domain.quantities['xvelocity'].boundary_values[5] == q_cons[1]/q_cons[0]
     656        assert domain.quantities['yvelocity'].boundary_values[5] == q_cons[2]/q_cons[0]
     657 
     658
    513659    def test_distribute_first_order(self):
    514660        """Domain implements a default first order gradient limiter
    … …  
    598744        #Set explicit_update
    599745
     746
    600747        for name in domain.conserved_quantities:
    601748            domain.quantities[name].explicit_update = num.array([4.,3.,2.,1.])
    … …  
    614761
    615762        x = num.array([1., 2., 3., 4.])
     763        x += domain.timestep*num.array( [4,3,2,1] )
    616764        x /= denom
    617         x += domain.timestep*num.array( [4,3,2,1] )
     765
    618766
    619767        for name in domain.conserved_quantities:
    … …  
    9061054
    9071055if __name__ == "__main__":
    908     suite = unittest.makeSuite(Test_Domain,'test')
     1056    suite = unittest.makeSuite(Test_Domain,'test_conserved_evolved')
    9091057    runner = unittest.TextTestRunner()
    9101058    runner.run(suite)

  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py

    r7505 r7519  
    6868
    6969        domain.conserved_quantities = ['stage', 'ymomentum']
     70        domain.evolved_quantities = ['stage', 'ymomentum']       
    7071        domain.quantities['stage'] =\
    7172                                   Quantity(domain, [[1,2,3], [5,5,5],
    … …  
    133134
    134135        domain.conserved_quantities = ['stage', 'ymomentum']
     136        domain.evolved_quantities = ['stage', 'ymomentum']       
    135137        domain.quantities['stage'] =\
    136138                                   Quantity(domain, [[1,2,3], [5,5,5],
    … …  
    178180
    179181
    180         q = T.evaluate(1, 1)  #Vol=0, edge=2
     182        q = T.evaluate(1, 1)  #Vol=1, edge=1
    181183        assert num.allclose(q, domain.get_edge_midpoint_coordinate(1,1))       
    182184
    … …  
    204206
    205207        domain.conserved_quantities = ['stage', 'ymomentum']
     208        domain.evolved_quantities = ['stage', 'ymomentum']       
    206209        domain.quantities['stage'] =\
    207210                                   Quantity(domain, [[1,2,3], [5,5,5],
    … …  
    362365        domain = Domain(points, elements)
    363366        domain.conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
     367        domain.evolved_quantities = ['stage', 'xmomentum', 'ymomentum']
    364368        domain.quantities['stage'] =\
    365369                                   Quantity(domain, [[1,2,3], [5,5,5],

  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_ghost.py

    r7276 r7519  
    3535
    3636        domain = Domain(points, vertices, None,
    37                         conserved_quantities, other_quantities)
     37                        conserved_quantities, None, other_quantities)
    3838        domain.check_integrity()
    3939

  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

    r7276 r7519  
    18671867
    18681868        x = num.array([1., 2., 3., 4.])
     1869        x += timestep*num.array( [4.0, 3.0, 2.0, 1.0] )       
    18691870        x /= denom
    1870         x += timestep*num.array( [4.0, 3.0, 2.0, 1.0] )
     1871
    18711872
    18721873        assert num.allclose( quantity.centroid_values, x)

  • anuga_core/source/anuga/interface.py

    r7452 r7519  
    1818
    1919from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Time_boundary
     20from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Time_space_boundary
    2021from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import Transmissive_boundary
    2122
    … …  
    3940    points, vertices, boundary = rectangular_cross(*args, **kwargs)
    4041    return Domain(points, vertices, boundary)
    41    
     42
    4243#----------------------------
    4344# Create domain from file

  • anuga_core/source/anuga/shallow_water/__init__.py

    r7326 r7519  
    1616     Transmissive_n_momentum_zero_t_momentum_set_stage_boundary,\
    1717     AWI_boundary
    18      
     18
     19#from shallow_water_balanced_domain import Swb_domain
    1920
    2021

  • anuga_core/source/anuga/shallow_water/data_manager.py

    r7509 r7519  
    24842484        #P = interp.mesh.get_boundary_polygon()
    24852485        #inside_indices = inside_polygon(grid_points, P)
     2486
     2487        # change printoptions so that a long string of zeros in not
     2488        # summarized as [0.0, 0.0, 0.0, ... 0.0, 0.0, 0.0]
     2489        printoptions = num.get_printoptions()
     2490        num.set_printoptions(threshold=sys.maxint)
     2491
    24862492        for i in range(nrows):
    24872493            if verbose and i % ((nrows+10)/10) == 0:
    … …  
    24922498            slice = grid_values[base_index:base_index+ncols]
    24932499            #s = array2string(slice, max_line_width=sys.maxint)
     2500
    24942501            s = num.array2string(slice, max_line_width=sys.maxint,
    24952502                                 precision=number_of_decimal_places)
    24962503            ascid.write(s[1:-1] + '\n')
    24972504
     2505        num.set_printoptions(threshold=printoptions['threshold'])
     2506       
    24982507        #Close
    24992508        ascid.close()
    … …  
    58825891    # brief Instantiate the SWW writer class.
    58835892    def __init__(self, static_quantities, dynamic_quantities):
    5884         """Initialise SWW writer with two list af quantity names:
     5893        """Initialise Write_sww with two list af quantity names:
    58855894       
    58865895        static_quantities (e.g. elevation or friction):

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

    r7452 r7519  
    147147    # @param use_cache
    148148    # @param verbose
     149    # @param evolved_quantities
    149150    # @param full_send_dict
    150151    # @param ghost_recv_dict
    … …  
    163164                 use_cache=False,
    164165                 verbose=False,
     166                 evolved_quantities = None,
     167                 other_quantities = None,
    165168                 full_send_dict=None,
    166169                 ghost_recv_dict=None,
    … …  
    171174
    172175        # Define quantities for the shallow_water domain         
    173         conserved_quantities = ['stage', 'xmomentum', 'ymomentum']         
    174         other_quantities = ['elevation', 'friction']
     176        conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
     177
     178        if other_quantities == None:
     179            other_quantities = ['elevation', 'friction']
    175180       
    176181        Generic_Domain.__init__(self,
    … …  
    179184                                boundary,
    180185                                conserved_quantities,
     186                                evolved_quantities,                               
    181187                                other_quantities,
    182188                                tagged_elements,
    … …  
    208214        self.optimise_dry_cells = optimise_dry_cells
    209215
     216        self.use_new_mannings = False
    210217        self.forcing_terms.append(manning_friction_implicit)
    211218        self.forcing_terms.append(gravity)
    … …  
    225232        self.use_centroid_velocities = use_centroid_velocities
    226233
     234
     235    ##
     236    # @brief
     237    # @param flag
     238    def set_new_mannings_function(self, flag=True):
     239        """Cludge to allow unit test to pass, but to
     240        also introduce new mannings friction function
     241        which takes into account the slope of the bed.
     242        The flag is tested in the python wrapper
     243        mannings_friction_implicit
     244        """
     245        if flag:
     246            self.use_new_mannings = True
     247        else:
     248            self.use_new_mannings = False
     249
     250
     251    ##
     252    # @brief
     253    # @param flag
     254    def set_use_edge_limiter(self, flag=True):
     255        """Cludge to allow unit test to pass, but to
     256        also introduce new edge limiting. The flag is
     257        tested in distribute_to_vertices_and_edges
     258        """
     259        if flag:
     260            self.use_edge_limiter = True
     261        else:
     262            self.use_edge_limiter = False
     263
     264
     265         
    227266    ##
    228267    # @brief
    … …  
    488527
    489528
    490 
    491 
    492     ##
    493     # @brief Get the total flow through an arbitrary poly line.
    494     # @param polyline Representation of desired cross section.
    495     # @param verbose True if this method is to be verbose.
    496     # @note 'polyline' may contain multiple sections allowing complex shapes.
    497     # @note Assume absolute UTM coordinates.
    498     def old_get_flow_through_cross_section(self, polyline, verbose=False):
    499         """Get the total flow through an arbitrary poly line.
    500 
    501         This is a run-time equivalent of the function with same name
    502         in data_manager.py
    503 
    504         Input:
    505             polyline: Representation of desired cross section - it may contain
    506                       multiple sections allowing for complex shapes. Assume
    507                       absolute UTM coordinates.
    508                       Format [[x0, y0], [x1, y1], ...]
    509 
    510         Output:
    511             Q: Total flow [m^3/s] across given segments.
    512         """
    513 
    514         # Find all intersections and associated triangles.
    515         segments = self.get_intersecting_segments(polyline, use_cache=True,
    516                                                   verbose=verbose)
    517 
    518         # Get midpoints
    519         midpoints = segment_midpoints(segments)
    520 
    521         # Make midpoints Geospatial instances
    522         midpoints = ensure_geospatial(midpoints, self.geo_reference)
    523 
    524         # Compute flow
    525         if verbose:
    526             log.critical('Computing flow through specified cross section')
    527 
    528         # Get interpolated values
    529         xmomentum = self.get_quantity('xmomentum')
    530         ymomentum = self.get_quantity('ymomentum')
    531 
    532         uh = xmomentum.get_values(interpolation_points=midpoints,
    533                                   use_cache=True)
    534         vh = ymomentum.get_values(interpolation_points=midpoints,
    535                                   use_cache=True)
    536 
    537         # Compute and sum flows across each segment
    538         total_flow = 0
    539         for i in range(len(uh)):
    540             # Inner product of momentum vector with segment normal [m^2/s]
    541             normal = segments[i].normal
    542             normal_momentum = uh[i]*normal[0] + vh[i]*normal[1]
    543 
    544             # Flow across this segment [m^3/s]
    545             segment_flow = normal_momentum*segments[i].length
    546 
    547             # Accumulate
    548             total_flow += segment_flow
    549 
    550         return total_flow
    551529
    552530    ##
    … …  
    589567
    590568        return cross_section.get_energy_through_cross_section(kind)
    591 
    592 
    593     ##
    594     # @brief
    595     # @param polyline Representation of desired cross section.
    596     # @param kind Select energy type to compute ('specific' or 'total').
    597     # @param verbose True if this method is to be verbose.
    598     # @note 'polyline' may contain multiple sections allowing complex shapes.
    599     # @note Assume absolute UTM coordinates.
    600     def old_get_energy_through_cross_section(self, polyline,
    601                                          kind='total',
    602                                          verbose=False):
    603         """Obtain average energy head [m] across specified cross section.
    604 
    605         Inputs:
    606             polyline: Representation of desired cross section - it may contain
    607                       multiple sections allowing for complex shapes. Assume
    608                       absolute UTM coordinates.
    609                       Format [[x0, y0], [x1, y1], ...]
    610             kind:     Select which energy to compute.
    611                       Options are 'specific' and 'total' (default)
    612 
    613         Output:
    614             E: Average energy [m] across given segments for all stored times.
    615 
    616         The average velocity is computed for each triangle intersected by
    617         the polyline and averaged weighted by segment lengths.
    618 
    619         The typical usage of this function would be to get average energy of
    620         flow in a channel, and the polyline would then be a cross section
    621         perpendicular to the flow.
    622 
    623         #FIXME (Ole) - need name for this energy reflecting that its dimension
    624         is [m].
    625         """
    626 
    627         from anuga.config import g, epsilon, velocity_protection as h0
    628 
    629         # Find all intersections and associated triangles.
    630         segments = self.get_intersecting_segments(polyline, use_cache=True,
    631                                                   verbose=verbose)
    632 
    633         # Get midpoints
    634         midpoints = segment_midpoints(segments)
    635 
    636         # Make midpoints Geospatial instances
    637         midpoints = ensure_geospatial(midpoints, self.geo_reference)
    638 
    639         # Compute energy
    640         if verbose: log.critical('Computing %s energy' % kind)
    641 
    642         # Get interpolated values
    643         stage = self.get_quantity('stage')
    644         elevation = self.get_quantity('elevation')
    645         xmomentum = self.get_quantity('xmomentum')
    646         ymomentum = self.get_quantity('ymomentum')
    647 
    648         w = stage.get_values(interpolation_points=midpoints, use_cache=True)
    649         z = elevation.get_values(interpolation_points=midpoints, use_cache=True)
    650         uh = xmomentum.get_values(interpolation_points=midpoints,
    651                                   use_cache=True)
    652         vh = ymomentum.get_values(interpolation_points=midpoints,
    653                                   use_cache=True)
    654         h = w-z                # Depth
    655 
    656         # Compute total length of polyline for use with weighted averages
    657         total_line_length = 0.0
    658         for segment in segments:
    659             total_line_length += segment.length
    660 
    661         # Compute and sum flows across each segment
    662         average_energy = 0.0
    663         for i in range(len(w)):
    664             # Average velocity across this segment
    665             if h[i] > epsilon:
    666                 # Use protection against degenerate velocities
    667                 u = uh[i]/(h[i] + h0/h[i])
    668                 v = vh[i]/(h[i] + h0/h[i])
    669             else:
    670                 u = v = 0.0
    671 
    672             speed_squared = u*u + v*v
    673             kinetic_energy = 0.5*speed_squared/g
    674 
    675             if kind == 'specific':
    676                 segment_energy = h[i] + kinetic_energy
    677             elif kind == 'total':
    678                 segment_energy = w[i] + kinetic_energy
    679             else:
    680                 msg = 'Energy kind must be either "specific" or "total".'
    681                 msg += ' I got %s' %kind
    682 
    683             # Add to weighted average
    684             weigth = segments[i].length/total_line_length
    685             average_energy += segment_energy*weigth
    686 
    687         return average_energy
    688569
    689570
    … …  
    18661747    from shallow_water_ext import gravity as gravity_c
    18671748
    1868     xmom = domain.quantities['xmomentum'].explicit_update
    1869     ymom = domain.quantities['ymomentum'].explicit_update
     1749    xmom_update = domain.quantities['xmomentum'].explicit_update
     1750    ymom_update = domain.quantities['ymomentum'].explicit_update
    18701751
    18711752    stage = domain.quantities['stage']
    … …  
    18781759    g = domain.g
    18791760
    1880     gravity_c(g, h, z, x, xmom, ymom)    #, 1.0e-6)
     1761    gravity_c(g, h, z, x, xmom_update, ymom_update)    #, 1.0e-6)
    18811762
    18821763##
    … …  
    18891770    """
    18901771
    1891     from shallow_water_ext import manning_friction as manning_friction_c
     1772    from shallow_water_ext import manning_friction_old
     1773    from shallow_water_ext import manning_friction_new
    18921774
    18931775    xmom = domain.quantities['xmomentum']
    18941776    ymom = domain.quantities['ymomentum']
    18951777
     1778    x = domain.get_vertex_coordinates()
     1779   
    18961780    w = domain.quantities['stage'].centroid_values
    1897     z = domain.quantities['elevation'].centroid_values
     1781    z = domain.quantities['elevation'].vertex_values
    18981782
    18991783    uh = xmom.centroid_values
    … …  
    19081792    g = domain.g
    19091793
    1910     manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
     1794    if domain.use_new_mannings:
     1795        manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, ymom_update)
     1796    else:
     1797        manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, ymom_update)
     1798   
    19111799
    19121800##
    … …  
    19191807    """
    19201808
    1921     from shallow_water_ext import manning_friction as manning_friction_c
     1809    from shallow_water_ext import manning_friction_old
     1810    from shallow_water_ext import manning_friction_new
    19221811
    19231812    xmom = domain.quantities['xmomentum']
    19241813    ymom = domain.quantities['ymomentum']
    19251814
     1815    x = domain.get_vertex_coordinates()
     1816   
    19261817    w = domain.quantities['stage'].centroid_values
    1927     z = domain.quantities['elevation'].centroid_values
     1818    z = domain.quantities['elevation'].vertex_values
    19281819
    19291820    uh = xmom.centroid_values
    … …  
    19381829    g = domain.g
    19391830
    1940     manning_friction_c(g, eps, w, z, uh, vh, eta, xmom_update, ymom_update)
     1831
     1832    if domain.use_new_mannings:
     1833        manning_friction_new(g, eps, x, w, uh, vh, z, eta, xmom_update, ymom_update)
     1834    else:
     1835        manning_friction_old(g, eps, w, uh, vh, z, eta, xmom_update, ymom_update)
     1836
    19411837
    19421838

  • anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r7276 r7519  
    516516
    517517
    518 void _manning_friction(double g, double eps, int N,
    519                double* w, double* z,
     518void _manning_friction_old(double g, double eps, int N,
     519               double* w, double* zv,
    520520               double* uh, double* vh,
    521521               double* eta, double* xmom, double* ymom) {
    522522
    523   int k;
    524   double S, h;
     523  int k, k3;
     524  double S, h, z, z0, z1, z2;
    525525
    526526  for (k=0; k<N; k++) {
    527527    if (eta[k] > eps) {
    528       h = w[k]-z[k];
     528      k3 = 3*k;
     529      // Get bathymetry
     530      z0 = zv[k3 + 0];
     531      z1 = zv[k3 + 1];
     532      z2 = zv[k3 + 2];
     533      z = (z0+z1+z2)/3.0;
     534      h = w[k]-z;
    529535      if (h >= eps) {
    530536        S = -g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
    … …  
    542548}
    543549
     550
     551void _manning_friction_new(double g, double eps, int N,
     552                           double* x, double* w, double* zv,
     553                           double* uh, double* vh,
     554                           double* eta, double* xmom_update, double* ymom_update) {
     555
     556  int k, k3, k6;
     557  double S, h, z, z0, z1, z2, zs, zx, zy;
     558  double x0,y0,x1,y1,x2,y2;
     559
     560  for (k=0; k<N; k++) {
     561    if (eta[k] > eps) {
     562      k3 = 3*k;
     563      // Get bathymetry
     564      z0 = zv[k3 + 0];
     565      z1 = zv[k3 + 1];
     566      z2 = zv[k3 + 2];
     567
     568      // Compute bed slope
     569      k6 = 6*k;  // base index
     570 
     571      x0 = x[k6 + 0];
     572      y0 = x[k6 + 1];
     573      x1 = x[k6 + 2];
     574      y1 = x[k6 + 3];
     575      x2 = x[k6 + 4];
     576      y2 = x[k6 + 5];
     577
     578      _gradient(x0, y0, x1, y1, x2, y2, z0, z1, z2, &zx, &zy);
     579
     580      zs = sqrt(1.0 + zx*zx + zy*zy);
     581      z = (z0+z1+z2)/3.0;
     582      h = w[k]-z;
     583      if (h >= eps) {
     584        S = -g * eta[k]*eta[k] * zs * sqrt((uh[k]*uh[k] + vh[k]*vh[k]));
     585        S /= pow(h, 7.0/3);      //Expensive (on Ole's home computer)
     586        //S /= exp((7.0/3.0)*log(h));      //seems to save about 15% over manning_friction
     587        //S /= h*h*(1 + h/3.0 - h*h/9.0); //FIXME: Could use a Taylor expansion
     588
     589
     590        //Update momentum
     591        xmom_update[k] += S*uh[k];
     592        ymom_update[k] += S*vh[k];
     593      }
     594    }
     595  }
     596}
    544597
    545598/*
    … …  
    10431096
    10441097
    1045 PyObject *manning_friction(PyObject *self, PyObject *args) {
     1098PyObject *manning_friction_new(PyObject *self, PyObject *args) {
    10461099  //
    1047   // manning_friction(g, eps, h, uh, vh, eta, xmom_update, ymom_update)
     1100  // manning_friction_new(g, eps, x, h, uh, vh, z, eta, xmom_update, ymom_update)
     1101  //
     1102
     1103
     1104  PyArrayObject *x, *w, *z, *uh, *vh, *eta, *xmom, *ymom;
     1105  int N;
     1106  double g, eps;
     1107
     1108  if (!PyArg_ParseTuple(args, "ddOOOOOOOO",
     1109                        &g, &eps, &x, &w, &uh, &vh, &z,  &eta, &xmom, &ymom)) {
     1110    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction_new could not parse input arguments");
     1111    return NULL;
     1112  }
     1113
     1114  // check that numpy array objects arrays are C contiguous memory
     1115  CHECK_C_CONTIG(x);
     1116  CHECK_C_CONTIG(w);
     1117  CHECK_C_CONTIG(z);
     1118  CHECK_C_CONTIG(uh);
     1119  CHECK_C_CONTIG(vh);
     1120  CHECK_C_CONTIG(eta);
     1121  CHECK_C_CONTIG(xmom);
     1122  CHECK_C_CONTIG(ymom);
     1123
     1124  N = w -> dimensions[0];
     1125
     1126  _manning_friction_new(g, eps, N,
     1127                        (double*) x -> data,
     1128                        (double*) w -> data,
     1129                        (double*) z -> data,
     1130                        (double*) uh -> data,
     1131                        (double*) vh -> data,
     1132                        (double*) eta -> data,
     1133                        (double*) xmom -> data,
     1134                        (double*) ymom -> data);
     1135
     1136  return Py_BuildValue("");
     1137}
     1138
     1139
     1140PyObject *manning_friction_old(PyObject *self, PyObject *args) {
     1141  //
     1142  // manning_friction_old(g, eps, h, uh, vh, z, eta, xmom_update, ymom_update)
    10481143  //
    10491144
    … …  
    10541149
    10551150  if (!PyArg_ParseTuple(args, "ddOOOOOOO",
    1056             &g, &eps, &w, &z, &uh, &vh, &eta,
    1057             &xmom, &ymom)) {
    1058     PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction could not parse input arguments");
     1151                        &g, &eps, &w, &uh, &vh, &z,  &eta, &xmom, &ymom)) {
     1152    PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction_old could not parse input arguments");
    10591153    return NULL;
    10601154  }
    … …  
    10701164
    10711165  N = w -> dimensions[0];
    1072   _manning_friction(g, eps, N,
    1073             (double*) w -> data,
    1074             (double*) z -> data,
    1075             (double*) uh -> data,
    1076             (double*) vh -> data,
    1077             (double*) eta -> data,
    1078             (double*) xmom -> data,
    1079             (double*) ymom -> data);
     1166
     1167  _manning_friction_old(g, eps, N,
     1168                        (double*) w -> data,
     1169                        (double*) z -> data,
     1170                        (double*) uh -> data,
     1171                        (double*) vh -> data,
     1172                        (double*) eta -> data,
     1173                        (double*) xmom -> data,
     1174                        (double*) ymom -> data);
    10801175
    10811176  return Py_BuildValue("");
    … …  
    27212816  {"compute_fluxes_ext_kinetic", compute_fluxes_ext_kinetic, METH_VARARGS, "Print out"},
    27222817  {"gravity", gravity, METH_VARARGS, "Print out"},
    2723   {"manning_friction", manning_friction, METH_VARARGS, "Print out"},
     2818  {"manning_friction_old", manning_friction_old, METH_VARARGS, "Print out"},
     2819  {"manning_friction_new", manning_friction_new, METH_VARARGS, "Print out"},
    27242820  {"flux_function_central", flux_function_central, METH_VARARGS, "Print out"}, 
    27252821  {"balance_deep_and_shallow", balance_deep_and_shallow,

  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r7509 r7519  
    352352
    353353        extrema = fid.variables['xmomentum.extrema'][:]
    354         assert num.allclose(extrema,[-0.06062178, 0.47873023]) or\
    355             num.allclose(extrema, [-0.06062178, 0.47847986]) or\
    356             num.allclose(extrema, [-0.06062178, 0.47848481]) # 27/5/9           
     354        assert num.allclose(extrema,[-0.06062178, 0.47873023]) or \
     355            num.allclose(extrema, [-0.06062178, 0.47847986]) or \
     356            num.allclose(extrema, [-0.06062178, 0.47848481]) or \
     357            num.allclose(extrema, [-0.06062178, 0.47763887]) # 18/09/09
    357358       
    358359        extrema = fid.variables['ymomentum.extrema'][:]
    … …  
    20962097        fid.close()
    20972098
     2099        #Cleanup
     2100        os.remove(prjfile)
     2101        os.remove(ascfile)
     2102        os.remove(swwfile)
     2103
     2104
     2105
     2106    def test_sww2dem_larger_zero(self):
     2107        """Test that sww information can be converted correctly to asc/prj
     2108        format readable by e.g. ArcView. Here:
     2109
     2110        ncols         11
     2111        nrows         11
     2112        xllcorner     308500
     2113        yllcorner     6189000
     2114        cellsize      10.000000
     2115        NODATA_value  -9999
     2116        -100 -110 -120 -130 -140 -150 -160 -170 -180 -190 -200
     2117         -90 -100 -110 -120 -130 -140 -150 -160 -170 -180 -190
     2118         -80  -90 -100 -110 -120 -130 -140 -150 -160 -170 -180
     2119         -70  -80  -90 -100 -110 -120 -130 -140 -150 -160 -170
     2120         -60  -70  -80  -90 -100 -110 -120 -130 -140 -150 -160
     2121         -50  -60  -70  -80  -90 -100 -110 -120 -130 -140 -150
     2122         -40  -50  -60  -70  -80  -90 -100 -110 -120 -130 -140
     2123         -30  -40  -50  -60  -70  -80  -90 -100 -110 -120 -130
     2124         -20  -30  -40  -50  -60  -70  -80  -90 -100 -110 -120
     2125         -10  -20  -30  -40  -50  -60  -70  -80  -90 -100 -110
     2126           0  -10  -20  -30  -40  -50  -60  -70  -80  -90 -100
     2127
     2128        """
     2129
     2130        import time, os
     2131        from Scientific.IO.NetCDF import NetCDFFile
     2132
     2133        #Setup
     2134
     2135        from mesh_factory import rectangular
     2136
     2137        #Create basic mesh (100m x 100m)
     2138        points, vertices, boundary = rectangular(2, 2, 100, 100)
     2139
     2140        #Create shallow water domain
     2141        domain = Domain(points, vertices, boundary)
     2142        domain.default_order = 1
     2143
     2144        domain.set_name('datatest')
     2145
     2146        prjfile = domain.get_name() + '_elevation.prj'
     2147        ascfile = domain.get_name() + '_elevation.asc'
     2148        swwfile = domain.get_name() + '.sww'
     2149
     2150        domain.set_datadir('.')
     2151        domain.format = 'sww'
     2152        domain.smooth = True
     2153        domain.geo_reference = Geo_reference(56, 308500, 6189000)
     2154
     2155        #
     2156        domain.set_quantity('elevation', 0)
     2157        domain.set_quantity('stage', 0)
     2158
     2159        B = Transmissive_boundary(domain)
     2160        domain.set_boundary( {'left': B, 'right': B, 'top': B, 'bottom': B})
     2161
     2162
     2163        #
     2164        sww = SWW_file(domain)
     2165        sww.store_connectivity()
     2166        sww.store_timestep()
     2167       
     2168        domain.tight_slope_limiters = 1
     2169        domain.evolve_to_end(finaltime = 0.01)
     2170        sww.store_timestep()
     2171
     2172        # Set theshold for printoptions to somrthing small
     2173        # to pickup Rudy's error caused a long sequence of zeros printing
     2174        # as [ 0.0, 0.0, 0.0, ... 0.0, 0.0, 0.0] by num.array2string
     2175        printoptions = num.get_printoptions()
     2176        num.set_printoptions(threshold=9)
     2177       
     2178        cellsize = 10  #10m grid
     2179
     2180
     2181        #Check contents
     2182        #Get NetCDF
     2183
     2184        fid = NetCDFFile(sww.filename, netcdf_mode_r)
     2185
     2186        # Get the variables
     2187        x = fid.variables['x'][:]
     2188        y = fid.variables['y'][:]
     2189        z = fid.variables['elevation'][:]
     2190        time = fid.variables['time'][:]
     2191        stage = fid.variables['stage'][:]
     2192
     2193
     2194        #Export to ascii/prj files
     2195        sww2dem(domain.get_name(),
     2196                quantity = 'elevation',
     2197                cellsize = cellsize,
     2198                number_of_decimal_places = 9,
     2199                verbose = self.verbose,
     2200                format = 'asc',
     2201                block_size=2)
     2202
     2203
     2204        #Check prj (meta data)
     2205        prjid = open(prjfile)
     2206        lines = prjid.readlines()
     2207        prjid.close()
     2208
     2209        L = lines[0].strip().split()
     2210        assert L[0].strip().lower() == 'projection'
     2211        assert L[1].strip().lower() == 'utm'
     2212
     2213        L = lines[1].strip().split()
     2214        assert L[0].strip().lower() == 'zone'
     2215        assert L[1].strip().lower() == '56'
     2216
     2217        L = lines[2].strip().split()
     2218        assert L[0].strip().lower() == 'datum'
     2219        assert L[1].strip().lower() == 'wgs84'
     2220
     2221        L = lines[3].strip().split()
     2222        assert L[0].strip().lower() == 'zunits'
     2223        assert L[1].strip().lower() == 'no'
     2224
     2225        L = lines[4].strip().split()
     2226        assert L[0].strip().lower() == 'units'
     2227        assert L[1].strip().lower() == 'meters'
     2228
     2229        L = lines[5].strip().split()
     2230        assert L[0].strip().lower() == 'spheroid'
     2231        assert L[1].strip().lower() == 'wgs84'
     2232
     2233        L = lines[6].strip().split()
     2234        assert L[0].strip().lower() == 'xshift'
     2235        assert L[1].strip().lower() == '500000'
     2236
     2237        L = lines[7].strip().split()
     2238        assert L[0].strip().lower() == 'yshift'
     2239        assert L[1].strip().lower() == '10000000'
     2240
     2241        L = lines[8].strip().split()
     2242        assert L[0].strip().lower() == 'parameters'
     2243
     2244
     2245        #Check asc file
     2246        ascid = open(ascfile)
     2247        lines = ascid.readlines()
     2248        ascid.close()
     2249
     2250        L = lines[0].strip().split()
     2251        assert L[0].strip().lower() == 'ncols'
     2252        assert L[1].strip().lower() == '11'
     2253
     2254        L = lines[1].strip().split()
     2255        assert L[0].strip().lower() == 'nrows'
     2256        assert L[1].strip().lower() == '11'
     2257
     2258        L = lines[2].strip().split()
     2259        assert L[0].strip().lower() == 'xllcorner'
     2260        assert num.allclose(float(L[1].strip().lower()), 308500)
     2261
     2262        L = lines[3].strip().split()
     2263        assert L[0].strip().lower() == 'yllcorner'
     2264        assert num.allclose(float(L[1].strip().lower()), 6189000)
     2265
     2266        L = lines[4].strip().split()
     2267        assert L[0].strip().lower() == 'cellsize'
     2268        assert num.allclose(float(L[1].strip().lower()), cellsize)
     2269
     2270        L = lines[5].strip().split()
     2271        assert L[0].strip() == 'NODATA_value'
     2272        assert L[1].strip().lower() == '-9999'
     2273
     2274        #Check grid values (FIXME: Use same strategy for other sww2dem tests)
     2275        for i, line in enumerate(lines[6:]):
     2276            for j, value in enumerate( line.split() ):
     2277                assert num.allclose(float(value), 0.0,
     2278                                    atol=1.0e-12, rtol=1.0e-12)
     2279
     2280                # Note: Equality can be obtained in this case,
     2281                # but it is better to use allclose.
     2282                #assert float(value) == -(10-i+j)*cellsize
     2283
     2284
     2285        fid.close()
     2286
     2287
     2288        num.set_printoptions(threshold=printoptions['threshold'])
    20982289        #Cleanup
    20992290        os.remove(prjfile)
    … …  
    1171211903if __name__ == "__main__":
    1171311904    #suite = unittest.makeSuite(Test_Data_Manager, 'test_sww2domain2')
    11714     suite = unittest.makeSuite(Test_Data_Manager, 'test')
     11905    suite = unittest.makeSuite(Test_Data_Manager, 'test_sww')
    1171511906   
    1171611907   

  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r7492 r7519  
    20672067                            4*S)
    20682068
     2069
     2070
     2071    def test_manning_friction_old(self):
     2072        from anuga.config import g
     2073
     2074        a = [0.0, 0.0]
     2075        b = [0.0, 2.0]
     2076        c = [2.0, 0.0]
     2077        d = [0.0, 4.0]
     2078        e = [2.0, 2.0]
     2079        f = [4.0, 0.0]
     2080
     2081        points = [a, b, c, d, e, f]
     2082        #             bac,     bce,     ecf,     dbe
     2083        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     2084
     2085        domain = Domain(points, vertices)
     2086
     2087        #Set up for a gradient of (3,0) at mid triangle (bce)
     2088        def slope(x, y):
     2089            return 3*x
     2090
     2091        h = 0.1
     2092        def stage(x, y):
     2093            return slope(x, y) + h
     2094
     2095        eta = 0.07
     2096        domain.set_quantity('elevation', slope)
     2097        domain.set_quantity('stage', stage)
     2098        domain.set_quantity('friction', eta)
     2099
     2100        for name in domain.conserved_quantities:
     2101            assert num.allclose(domain.quantities[name].explicit_update, 0)
     2102            assert num.allclose(domain.quantities[name].semi_implicit_update, 0)
     2103
     2104        domain.compute_forcing_terms()
     2105
     2106        assert num.allclose(domain.quantities['stage'].explicit_update, 0)
     2107        assert num.allclose(domain.quantities['xmomentum'].explicit_update,
     2108                            -g*h*3)
     2109        assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
     2110
     2111        assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
     2112        assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update,
     2113                            0)
     2114        assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update,
     2115                            0)
     2116
     2117        #Create some momentum for friction to work with
     2118        domain.set_quantity('xmomentum', 1)
     2119        S = -g*eta**2 / h**(7.0/3)
     2120
     2121        domain.compute_forcing_terms()
     2122        assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
     2123        assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update,
     2124                            S)
     2125        assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update,
     2126                            0)
     2127
     2128        #A more complex example
     2129        domain.quantities['stage'].semi_implicit_update[:] = 0.0
     2130        domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0
     2131        domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0
     2132
     2133        domain.set_quantity('xmomentum', 3)
     2134        domain.set_quantity('ymomentum', 4)
     2135
     2136        S = -g*eta**2*5 / h**(7.0/3)
     2137
     2138        domain.compute_forcing_terms()
     2139
     2140        assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
     2141        assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update,
     2142                            3*S)
     2143        assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update,
     2144                            4*S)
     2145
     2146
     2147    def test_manning_friction_new(self):
     2148        from anuga.config import g
     2149
     2150        a = [0.0, 0.0]
     2151        b = [0.0, 2.0]
     2152        c = [2.0, 0.0]
     2153        d = [0.0, 4.0]
     2154        e = [2.0, 2.0]
     2155        f = [4.0, 0.0]
     2156
     2157        points = [a, b, c, d, e, f]
     2158        #             bac,     bce,     ecf,     dbe
     2159        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
     2160
     2161        domain = Domain(points, vertices)
     2162
     2163        # Use the new function which takes into account the extra
     2164        # wetted area due to slope of bed
     2165        domain.set_new_mannings_function(True)
     2166       
     2167        #Set up for a gradient of (3,0) at mid triangle (bce)
     2168        def slope(x, y):
     2169            return 3*x
     2170
     2171        h = 0.1
     2172        def stage(x, y):
     2173            return slope(x, y) + h
     2174
     2175        eta = 0.07
     2176        domain.set_quantity('elevation', slope)
     2177        domain.set_quantity('stage', stage)
     2178        domain.set_quantity('friction', eta)
     2179
     2180        for name in domain.conserved_quantities:
     2181            assert num.allclose(domain.quantities[name].explicit_update, 0)
     2182            assert num.allclose(domain.quantities[name].semi_implicit_update, 0)
     2183
     2184        domain.compute_forcing_terms()
     2185
     2186        assert num.allclose(domain.quantities['stage'].explicit_update, 0)
     2187        assert num.allclose(domain.quantities['xmomentum'].explicit_update,
     2188                            -g*h*3)
     2189        assert num.allclose(domain.quantities['ymomentum'].explicit_update, 0)
     2190
     2191        assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
     2192        assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update,
     2193                            0)
     2194        assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update,
     2195                            0)
     2196
     2197        #Create some momentum for friction to work with
     2198        domain.set_quantity('xmomentum', 1)
     2199        S = -g*eta**2 / h**(7.0/3) * sqrt(10)
     2200
     2201        domain.compute_forcing_terms()
     2202        assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
     2203        assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update,
     2204                            S)
     2205        assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update,
     2206                            0)
     2207
     2208        #A more complex example
     2209        domain.quantities['stage'].semi_implicit_update[:] = 0.0
     2210        domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0
     2211        domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0
     2212
     2213        domain.set_quantity('xmomentum', 3)
     2214        domain.set_quantity('ymomentum', 4)
     2215
     2216        S = -g*eta**2*5 / h**(7.0/3) * sqrt(10.0)
     2217
     2218        domain.compute_forcing_terms()
     2219
     2220        assert num.allclose(domain.quantities['stage'].semi_implicit_update, 0)
     2221        assert num.allclose(domain.quantities['xmomentum'].semi_implicit_update,
     2222                            3*S)
     2223        assert num.allclose(domain.quantities['ymomentum'].semi_implicit_update,
     2224                            4*S)
     2225       
    20692226    def test_constant_wind_stress(self):
    20702227        from anuga.config import rho_a, rho_w, eta_w
    … …  
    73347491
    73357492if __name__ == "__main__":
    7336     suite = unittest.makeSuite(Test_Shallow_Water, 'test_volumetric_balance')
    7337     #suite = unittest.makeSuite(Test_Shallow_Water, 'test')
     7493    #suite = unittest.makeSuite(Test_Shallow_Water, 'test_volumetric_balance')
     7494    suite = unittest.makeSuite(Test_Shallow_Water, 'test')
    73387495    runner = unittest.TextTestRunner(verbosity=1)
    73397496    runner.run(suite)

  • anuga_core/source/anuga/utilities/log.py

    r7276 r7519  
    1111    log.console_logging_level = log.INFO
    1212    log.log_logging_level = log.DEBUG
    13     log.log_filename('./my.log')
     13    log.log_filename = './my.log'
    1414
    1515    # log away!

  • anuga_core/source/anuga_parallel/test_parallel.py

    r7507 r7519  
    169169            assert(num.allclose(submesh['ghost_triangles'][1],true_ghost_triangles[1]))
    170170            assert num.allclose(submesh['ghost_triangles'][2],true_ghost_triangles[2]), ParallelException('X')
    171             #if not num.allclose(submesh['ghost_triangles'][2],true_ghost_triangles[2]):
    172             #    import pypar
    173             #    raise Exception, 'Here it is'
    174             #    pypar.abort()
    175 
    176171
    177172           
    … …  
    295290
    296291#-------------------------------------------------------------
     292
    297293if __name__=="__main__":
    298294    runner = unittest.TextTestRunner()
Note: See TracChangeset for help on using the changeset viewer.