Changeset 8073


Ignore:
Timestamp:
Nov 19, 2010, 2:53:46 PM (13 years ago)
Author:
steve
Message:

Added in a unit test for inlet operator

Location:
trunk
Files:
2 added
21 edited
1 moved

Legend:

Unmodified
Added
Removed

  • trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_domain.py

    r8068 r8073  
    13211321    # @note Raises exception of method not known.
    13221322    def set_timestepping_method(self, timestepping_method):
    1323         methods = ['euler', 'rk2', 'rk3']   
     1323        methods = ['euler', 'rk2', 'rk3']
    13241324        if timestepping_method in methods:
    13251325            self.timestepping_method = timestepping_method
    13261326            return
    13271327        if timestepping_method in [1,2,3]:
    1328             self.timetepping_method = methods[timestepping_method-1]
     1328            self.timestepping_method = methods[timestepping_method-1]
    13291329            return
    13301330

  • trunk/anuga_core/source/anuga/geometry/polygon.py

    r8062 r8073  
    44
    55import numpy as num
     6import math
    67
    78from anuga.utilities.numerical_tools import ensure_numeric
    … …  
    274275
    275276    return indices[:count]
     277
     278
     279def line_length(line):
     280    """Determine the length of the line
     281    """
     282   
     283    l12 = line[1]-line[0]
     284
     285    return math.sqrt(num.dot(l12,l12))
    276286   
    277287def not_line_intersect(triangles, line, verbose=False):

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

    r7967 r8073  
    278278        from math import pi, cos, sin
    279279
     280        if domain.numproc > 1:
     281            msg = 'Not implemented to run in parallel'
     282            assert self.__parallel_safe(), msg
     283
    280284        if center is None:
    281285            msg = 'I got radius but no center.'
    … …  
    510514
    511515
     516    def __parallel_safe(self):
     517
     518        return false
    512519##
    513520# @brief A class for rainfall forcing function.

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

    r8068 r8073  
    173173
    174174        # Stored output
    175         self.store = True
     175        self.set_store(True)
    176176        self.set_store_vertices_uniquely(False)
    177177        self.quantities_to_be_stored = {'elevation': 1,
    … …  
    222222
    223223
     224    def set_store(self, flag=True):
     225        """Set whether data saved to sww file.
     226        """
     227
     228        self.store = flag
     229
     230       
    224231    def set_sloped_mannings_function(self, flag=True):
    225232        """Set mannings friction function to use the sloped

  • trunk/anuga_core/source/anuga/shallow_water_balanced/test_swb_forcing_terms.py

    r8072 r8073  
    99from anuga.config import g, epsilon
    1010from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    11 from anuga.utilities.numerical_tools import mean
     11from anuga.utilities.numerical_tools import mean, ensure_numeric
    1212from anuga.geometry.polygon import is_inside_polygon
    1313from anuga.coordinate_transforms.geo_reference import Geo_reference
    … …  
    13231323
    13241324       
    1325     def test_inflow_using_flowline(self):
     1325    def xtest_inflow_using_flowline(self):
    13261326        """test_inflow_using_flowline
    13271327
    … …  
    14301430
    14311431                for t in domain.evolve(yieldstep=100.0, finaltime=finaltime):
    1432                     pass
     1432                    print domain.timestepping_statistics()
    14331433                    #if verbose :
    14341434                    #    print domain.timestepping_statistics()
    … …  
    16101610if __name__ == "__main__":
    16111611    suite = unittest.makeSuite(Test_swb_forcing_terms, 'test')
    1612     runner = unittest.TextTestRunner(verbosity=1)
     1612    runner = unittest.TextTestRunner(verbosity=2)
    16131613    runner.run(suite)

  • trunk/anuga_core/source/anuga/structures/boyd_box_operator.py

    r8056 r8073  
    1010    compute_discharge method for specific subclasses)
    1111   
    12     Input: Two points, pipe_size (either diameter or width, height), 
     12    Input: Two points, pipe_size (either diameter or width, height),
    1313    mannings_rougness,
    1414    """
    … …  
    7777        self.case = 'N/A'
    7878
    79    
     79
     80
    8081    def discharge_routine(self):
    8182
    8283        local_debug ='false'
    83        
    84         if self.inflow.get_enquiry_height() > 0.01: #this value was 0.01:
     84
     85        if self.use_velocity_head:
     86            self.delta_total_energy = self.inlets[0].get_enquiry_total_energy() - self.inlets[1].get_enquiry_total_energy()
     87        else:
     88            self.delta_total_energy = self.inlets[0].get_enquiry_stage() - self.inlets[1].get_enquiry_stage()
     89
     90        self.inflow  = self.inlets[0]
     91        self.outflow = self.inlets[1]
     92
     93        if self.delta_total_energy < 0:
     94            self.inflow  = self.inlets[1]
     95            self.outflow = self.inlets[0]
     96            self.delta_total_energy = -self.delta_total_energy
     97
     98
     99
     100        if self.inflow.get_enquiry_depth() > 0.01: #this value was 0.01:
    85101            if local_debug =='true':
    86102                anuga.log.critical('Specific E & Deltat Tot E = %s, %s'
    … …  
    97113                self.driving_energy = self.inflow.get_enquiry_specific_energy()
    98114            else:
    99                 self.driving_energy = self.inflow.get_enquiry_height()
    100 
    101             height = self.culvert_height
     115                self.driving_energy = self.inflow.get_enquiry_depth()
     116
     117            depth = self.culvert_height
    102118            width = self.culvert_width
    103119            flow_width = self.culvert_width
    … …  
    106122            # but ensure the correct flow area and wetted perimeter are used
    107123            Q_inlet_unsubmerged = 0.544*anuga.g**0.5*width*self.driving_energy**1.50 # Flow based on Inlet Ctrl Inlet Unsubmerged
    108             Q_inlet_submerged = 0.702*anuga.g**0.5*width*height**0.89*self.driving_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
     124            Q_inlet_submerged = 0.702*anuga.g**0.5*width*depth**0.89*self.driving_energy**0.61  # Flow based on Inlet Ctrl Inlet Submerged
    109125
    110126            # FIXME(Ole): Are these functions really for inlet control?
    … …  
    112128                Q = Q_inlet_unsubmerged
    113129                dcrit = (Q**2/anuga.g/width**2)**0.333333
    114                 if dcrit > height:
    115                     dcrit = height
     130                if dcrit > depth:
     131                    dcrit = depth
    116132                    flow_area = width*dcrit
    117133                    perimeter= 2.0*(width+dcrit)
    118                 else: # dcrit < height
     134                else: # dcrit < depth
    119135                    flow_area = width*dcrit
    120136                    perimeter= 2.0*dcrit+width
    … …  
    124140                Q = Q_inlet_submerged
    125141                dcrit = (Q**2/anuga.g/width**2)**0.333333
    126                 if dcrit > height:
    127                     dcrit = height
     142                if dcrit > depth:
     143                    dcrit = depth
    128144                    flow_area = width*dcrit
    129145                    perimeter= 2.0*(width+dcrit)
    130                 else: # dcrit < height
     146                else: # dcrit < depth
    131147                    flow_area = width*dcrit
    132148                    perimeter= 2.0*dcrit+width
    … …  
    137153            # May not need this .... check if same is done above
    138154            outlet_culvert_depth = dcrit
    139             if outlet_culvert_depth > height:
    140                 outlet_culvert_depth = height  # Once again the pipe is flowing full not partfull
    141                 flow_area = width*height  # Cross sectional area of flow in the culvert
    142                 perimeter = 2*(width+height)
     155            if outlet_culvert_depth > depth:
     156                outlet_culvert_depth = depth  # Once again the pipe is flowing full not partfull
     157                flow_area = width*depth  # Cross sectional area of flow in the culvert
     158                perimeter = 2*(width+depth)
    143159                self.case = 'Inlet CTRL Outlet unsubmerged PIPE PART FULL'
    144160            else:
    … …  
    157173
    158174                # Determine the depth at the outlet relative to the depth of flow in the Culvert
    159                 if self.outflow.get_enquiry_height() > height:        # The Outlet is Submerged
    160                     outlet_culvert_depth=height
    161                     flow_area=width*height       # Cross sectional area of flow in the culvert
    162                     perimeter=2.0*(width+height)
     175                if self.outflow.get_enquiry_depth() > depth:        # The Outlet is Submerged
     176                    outlet_culvert_depth=depth
     177                    flow_area=width*depth       # Cross sectional area of flow in the culvert
     178                    perimeter=2.0*(width+depth)
    163179                    self.case = 'Outlet submerged'
    164180                else:   # Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    165181                    dcrit = (Q**2/anuga.g/width**2)**0.333333
    166182                    outlet_culvert_depth=dcrit   # For purpose of calculation assume the outlet depth = Critical Depth
    167                     if outlet_culvert_depth > height:
    168                         outlet_culvert_depth=height
    169                         flow_area=width*height
    170                         perimeter=2.0*(width+height)
     183                    if outlet_culvert_depth > depth:
     184                        outlet_culvert_depth=depth
     185                        flow_area=width*depth
     186                        perimeter=2.0*(width+depth)
    171187                        self.case = 'Outlet is Flowing Full'
    172188                    else:
    … …  
    201217        # END CODE BLOCK for DEPTH  > Required depth for CULVERT Flow
    202218
    203         else: # self.inflow.get_enquiry_height() < 0.01:
     219        else: # self.inflow.get_enquiry_depth() < 0.01:
    204220            Q = barrel_velocity = outlet_culvert_depth = 0.0
    205221

  • trunk/anuga_core/source/anuga/structures/boyd_pipe_operator.py

    r8056 r8073  
    3838                                          enquiry_points,
    3939                                          width=diameter,
    40                                           height=None,
     40                                          depth=None,
    4141                                          apron=apron,
    4242                                          manning=manning,
    … …  
    7373        self.case = 'N/A'
    7474       
    75    
     75
     76    def __determine_inflow_outflow(self):
     77        # Determine flow direction based on total energy difference
     78
     79        if self.use_velocity_head:
     80            self.delta_total_energy = self.inlets[0].get_enquiry_total_energy() - self.inlets[1].get_enquiry_total_energy()
     81        else:
     82            self.delta_total_energy = self.inlets[0].get_enquiry_stage() - self.inlets[1].get_enquiry_stage()
     83
     84        self.inflow  = self.inlets[0]
     85        self.outflow = self.inlets[1]
     86
     87        if self.delta_total_energy < 0:
     88            self.inflow  = self.inlets[1]
     89            self.outflow = self.inlets[0]
     90            self.delta_total_energy = -self.delta_total_energy
     91
    7692    def discharge_routine(self):
     93
     94        self.__determine_inflow_outflow()
    7795
    7896        local_debug ='false'
    … …  
    8199        #pdb.set_trace()
    82100       
    83         if self.inflow.get_enquiry_height() > 0.01: #this value was 0.01: Remember this needs to be compared to the Invert Lvl
     101        if self.inflow.get_enquiry_depth() > 0.01: #this value was 0.01: Remember this needs to be compared to the Invert Lvl
    84102            if local_debug =='true':
    85103                anuga.log.critical('Specific E & Deltat Tot E = %s, %s'
    … …  
    96114                self.driving_energy = self.inflow.get_enquiry_specific_energy()
    97115            else:
    98                 self.driving_energy = self.inflow.get_enquiry_height()
     116                self.driving_energy = self.inflow.get_enquiry_depth()
    99117                """
    100118        For a circular pipe the Boyd method reviews 3 conditions
    … …  
    109127
    110128        local_debug ='false'
    111         if self.inflow.get_average_height() > 0.01: #this should test against invert
     129        if self.inflow.get_average_depth() > 0.01: #this should test against invert
    112130            if local_debug =='true':
    113131                anuga.log.critical('Specific E & Deltat Tot E = %s, %s'
    … …  
    167185
    168186                # Determine the depth at the outlet relative to the depth of flow in the Culvert
    169                 if self.outflow.get_average_height() > diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
     187                if self.outflow.get_average_depth() > diameter:       # Outlet is submerged Assume the end of the Pipe is flowing FULL
    170188                    outlet_culvert_depth=diameter
    171189                    flow_area = (diameter/2)**2 * math.pi  # Cross sectional area of flow in the culvert
    … …  
    176194                        anuga.log.critical('Outlet submerged')
    177195                else:   # Culvert running PART FULL for PART OF ITS LENGTH   Here really should use the Culvert Slope to calculate Actual Culvert Depth & Velocity
    178                     # IF  self.outflow.get_average_height() < diameter
     196                    # IF  self.outflow.get_average_depth() < diameter
    179197                    dcrit1 = diameter/1.26*(Q/anuga.g**0.5*diameter**2.5)**(1/3.75)
    180198                    dcrit2 = diameter/0.95*(Q/anuga.g**0.5*diameter**2.5)**(1/1.95)
    … …  
    243261            barrel_velocity = Q/(flow_area + anuga.velocity_protection/flow_area)
    244262
    245         else: # self.inflow.get_average_height() < 0.01:
     263        else: # self.inflow.get_average_depth() < 0.01:
    246264            Q = barrel_velocity = outlet_culvert_depth = 0.0
    247265

  • trunk/anuga_core/source/anuga/structures/inlet.py

    r8056 r8073  
     1import anuga.geometry.polygon
    12from anuga.geometry.polygon import inside_polygon, is_inside_polygon, line_intersect
    23from anuga.config import velocity_protection, g
    … …  
    1819        self.compute_triangle_indices()
    1920        self.compute_area()
     21        self.compute_inlet_length()
     22
    2023
    2124
    … …  
    6164
    6265
     66    def compute_inlet_length(self):
     67        """ Compute the length of the inlet (as
     68        defined by the input line
     69        """
     70
     71        point0 = self.line[0]
     72        point1 = self.line[1]
     73
     74        self.inlet_length = anuga.geometry.polygon.line_length(self.line)
     75
     76
     77    def get_inlet_length(self):
     78
     79        return self.inlet_length
     80
     81    def get_line(self):
     82
     83        return self.line
     84       
    6385    def get_area(self):
    6486
    … …  
    110132   
    111133
    112     def get_heights(self):
     134    def get_depths(self):
    113135   
    114136        return self.get_stages() - self.get_elevations()
    … …  
    117139    def get_total_water_volume(self):
    118140       
    119        return num.sum(self.get_heights()*self.get_areas())
     141       return num.sum(self.get_depths()*self.get_areas())
    120142 
    121143
    122     def get_average_height(self):
     144    def get_average_depth(self):
    123145   
    124146        return self.get_total_water_volume()/self.area
    … …  
    127149    def get_velocities(self):
    128150       
    129             heights = self.get_heights()
    130             u = self.get_xmoms()/(heights + velocity_protection/heights)
    131             v = self.get_ymoms()/(heights + velocity_protection/heights)
     151            depths = self.get_depths()
     152            u = self.get_xmoms()/(depths + velocity_protection/depths)
     153            v = self.get_ymoms()/(depths + velocity_protection/depths)
    132154           
    133155            return u, v
    … …  
    136158    def get_xvelocities(self):
    137159
    138             heights = self.get_heights()
    139             return self.get_xmoms()/(heights + velocity_protection/heights)
     160            depths = self.get_depths()
     161            return self.get_xmoms()/(depths + velocity_protection/depths)
    140162
    141163    def get_yvelocities(self):
    142164
    143             heights = self.get_heights()
    144             return self.get_ymoms()/(heights + velocity_protection/heights)
     165            depths = self.get_depths()
     166            return self.get_ymoms()/(depths + velocity_protection/depths)
    145167           
    146168           
    … …  
    167189    def get_average_specific_energy(self):
    168190       
    169         return self.get_average_velocity_head() + self.get_average_height()
    170 
    171 
    172 
    173     def set_heights(self,height):
    174 
    175         self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, self.get_elevations() + height)
     191        return self.get_average_velocity_head() + self.get_average_depth()
     192
     193
     194
     195    def set_depths(self,depth):
     196
     197        self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, self.get_elevations() + depth)
    176198
    177199

  • trunk/anuga_core/source/anuga/structures/inlet_enquiry.py

    r8056 r8073  
    4545            raise Exception, msg
    4646
     47    def get_enquiry_position(self):
    4748
     49        return self.domain.get_centroid_coordinates(absolute=True)[self.enquiry_index]
    4850
    4951    def get_enquiry_stage(self):
    … …  
    6567        return self.domain.quantities['elevation'].centroid_values[self.enquiry_index]
    6668
    67     def get_enquiry_height(self):
     69    def get_enquiry_depth(self):
    6870
    6971        return self.get_enquiry_stage() - self.get_enquiry_elevation()
    … …  
    7274    def get_enquiry_velocity(self):
    7375
    74             height = self.get_enquiry_height()
    75             u = self.get_enquiry_xmom()/(height + velocity_protection/height)
    76             v = self.get_enquiry_ymom()/(height + velocity_protection/height)
     76            depth = self.get_enquiry_depth()
     77            u = self.get_enquiry_xmom()/(depth + velocity_protection/depth)
     78            v = self.get_enquiry_ymom()/(depth + velocity_protection/depth)
    7779
    7880            return u, v
    … …  
    8183    def get_enquiry_xvelocity(self):
    8284
    83             height = self.get_enquiry_height()
    84             return self.get_enquiry_xmom()/(height + velocity_protection/height)
     85            depth = self.get_enquiry_depth()
     86            return self.get_enquiry_xmom()/(depth + velocity_protection/depth)
    8587
    8688    def get_enquiry_yvelocity(self):
    8789
    88             height = self.get_enquiry_height()
    89             return self.get_enquiry_ymom()/(height + velocity_protection/height)
     90            depth = self.get_enquiry_depth()
     91            return self.get_enquiry_ymom()/(depth + velocity_protection/depth)
    9092
    9193
    … …  
    109111    def get_enquiry_specific_energy(self):
    110112
    111         return self.get_enquiry_velocity_head() + self.get_enquiry_height()
     113        return self.get_enquiry_velocity_head() + self.get_enquiry_depth()
    112114
    113115

  • trunk/anuga_core/source/anuga/structures/inlet_operator.py

    r8050 r8073  
    88
    99class Inlet_operator:
    10     """Inlet Operator - add water to one inlet.
     10    """Inlet Operator - add water to an inlet.
    1111    Sets up the geometry of problem
    1212   
    … …  
    1414    discharge_routine method for specific subclasses)
    1515   
    16     Input: Two points, pipe_size (either diameter or width, height),
    17     mannings_rougness,
     16    Input: domain, Two points
    1817    """
    1918
    … …  
    7473
    7574        # FIXME (SR): Might be nice to spread the over the inlet so that it is flat
    76         new_inlet_height = self.inlet.get_average_height() + (Q*timestep/self.inlet.get_area())
    77         self.inlet.set_heights(new_inlet_height)
     75        new_inlet_depth = self.inlet.get_average_depth() + (Q*timestep/self.inlet.get_area())
     76        self.inlet.set_depths(new_inlet_depth)
    7877
    7978           

  • trunk/anuga_core/source/anuga/structures/structure_operator.py

    r8056 r8073  
    55
    66from anuga.utilities.system_tools import log_to_file
     7from anuga.utilities.numerical_tools import ensure_numeric
     8
    79
    810
    … …  
    1416    discharge_routine method for specific subclasses)
    1517   
    16     Input: Two points, pipe_size (either diameter or width, height),
     18    Input: Two points, pipe_size (either diameter or width, depth),
    1719    mannings_rougness,
    1820    """
    … …  
    3840        self.domain = domain
    3941        self.domain.set_fractional_step_operator(self)
    40         self.end_points = end_points
    41         self.exchange_lines = exchange_lines
    42         self.enquiry_points = enquiry_points
     42        self.end_points = ensure_numeric(end_points)
     43        self.exchange_lines = ensure_numeric(exchange_lines)
     44        self.enquiry_points = ensure_numeric(enquiry_points)
     45
     46
     47        if domain.numproc > 1:
     48            msg = 'Not implemented to run in parallel'
     49            assert self.__parallel_safe(), msg
     50
    4351       
    4452        if height is None:
    … …  
    106114        timestep = self.domain.get_timestep()
    107115       
    108         self.determine_inflow_outflow()
    109        
    110116        Q, barrel_speed, outlet_depth = self.discharge_routine()
    111117
    112         old_inflow_height = self.inflow.get_average_height()
     118        old_inflow_depth = self.inflow.get_average_depth()
     119        old_inflow_stage = self.inflow.get_average_stage()
    113120        old_inflow_xmom = self.inflow.get_average_xmom()
    114121        old_inflow_ymom = self.inflow.get_average_ymom()
    115122
     123
    116124        # Implement the update of flow over a timestep by
    117125        # using a semi-implict update. This ensures that
    118         # the update does not create a negative height
    119         if old_inflow_height > 0.0 :
    120                 Q_star = Q/old_inflow_height
     126        # the update does not create a negative depth
     127        if old_inflow_depth > 0.0 :
     128                Q_star = Q/old_inflow_depth
    121129        else:
    122130                Q_star = 0.0
    … …  
    124132        factor = 1.0/(1.0 + Q_star*timestep/self.inflow.get_area())
    125133
    126         new_inflow_height = old_inflow_height*factor
     134        new_inflow_depth = old_inflow_depth*factor
    127135        new_inflow_xmom = old_inflow_xmom*factor
    128136        new_inflow_ymom = old_inflow_ymom*factor
    129137           
    130         self.inflow.set_heights(new_inflow_height)
     138        self.inflow.set_depths(new_inflow_depth)
    131139
    132140        #inflow.set_xmoms(Q/inflow.get_area())
    … …  
    136144        self.inflow.set_ymoms(new_inflow_ymom)
    137145
    138         loss = (old_inflow_height - new_inflow_height)*self.inflow.get_area()
     146        loss = (old_inflow_depth - new_inflow_depth)*self.inflow.get_area()
    139147
    140148        # set outflow
    141         if old_inflow_height > 0.0 :
    142                 timestep_star = timestep*new_inflow_height/old_inflow_height
     149        if old_inflow_depth > 0.0 :
     150                timestep_star = timestep*new_inflow_depth/old_inflow_depth
    143151        else:
    144152            timestep_star = 0.0
    145153
    146         outflow_extra_height = Q*timestep_star/self.outflow.get_area()
     154
     155
     156
     157        outflow_extra_depth = Q*timestep_star/self.outflow.get_area()
    147158        outflow_direction = - self.outflow.outward_culvert_vector
    148         outflow_extra_momentum = outflow_extra_height*barrel_speed*outflow_direction
    149            
    150         gain = outflow_extra_height*self.outflow.get_area()
     159        outflow_extra_momentum = outflow_extra_depth*barrel_speed*outflow_direction
     160           
     161        gain = outflow_extra_depth*self.outflow.get_area()
    151162           
    152163        #print Q, Q*timestep, barrel_speed, outlet_depth, Qstar, factor, timestep_star
    … …  
    154165
    155166        # Stats
    156         self.discharge  = Q#outflow_extra_height*self.outflow.get_area()/timestep
     167        self.discharge  = Q#outflow_extra_depth*self.outflow.get_area()/timestep
    157168        self.velocity = barrel_speed#self.discharge/outlet_depth/self.width
    158169
    159         new_outflow_height = self.outflow.get_average_height() + outflow_extra_height
     170        new_outflow_depth = self.outflow.get_average_depth() + outflow_extra_depth
    160171
    161172        if self.use_momentum_jet :
    … …  
    164175            #new_outflow_ymom = outflow.get_average_ymom() + outflow_extra_momentum[1]
    165176
    166             new_outflow_xmom = barrel_speed*new_outflow_height*outflow_direction[0]
    167             new_outflow_ymom = barrel_speed*new_outflow_height*outflow_direction[1]
     177            new_outflow_xmom = barrel_speed*new_outflow_depth*outflow_direction[0]
     178            new_outflow_ymom = barrel_speed*new_outflow_depth*outflow_direction[1]
    168179
    169180        else:
    … …  
    174185            new_outflow_ymom = 0.0
    175186
    176         self.outflow.set_heights(new_outflow_height)
     187        self.outflow.set_depths(new_outflow_depth)
    177188        self.outflow.set_xmoms(new_outflow_xmom)
    178189        self.outflow.set_ymoms(new_outflow_ymom)
    179190
    180191
    181     def determine_inflow_outflow(self):
    182         # Determine flow direction based on total energy difference
    183 
    184         if self.use_velocity_head:
    185             self.delta_total_energy = self.inlets[0].get_enquiry_total_energy() - self.inlets[1].get_enquiry_total_energy()
    186         else:
    187             self.delta_total_energy = self.inlets[0].get_enquiry_stage() - self.inlets[1].get_enquiry_stage()
    188 
    189         self.inflow  = self.inlets[0]
    190         self.outflow = self.inlets[1]
    191 
    192         if self.delta_total_energy < 0:
    193             self.inflow  = self.inlets[1]
    194             self.outflow = self.inlets[0]
    195             self.delta_total_energy = -self.delta_total_energy
     192
    196193
    197194
    … …  
    262259            self.enquiry_points.append(centre_point1 + gap)
    263260           
    264        
     261
     262    def __parallel_safe(self):
     263
     264        return False
     265
    265266    def discharge_routine(self):
    266        
    267         pass
     267
     268        msg = 'Need to impelement '
     269        raise
    268270           
    269271

  • trunk/anuga_core/source/anuga/structures/test_boyd_box_operator.py

    r8056 r8073  
    2424    def tearDown(self):
    2525        pass
    26 
     26   
     27   
     28    def _create_domain(self,d_length,
     29                            d_width,
     30                            dx,
     31                            dy,
     32                            elevation_0,
     33                            elevation_1,
     34                            stage_0,
     35                            stage_1):
     36       
     37        points, vertices, boundary = rectangular_cross(int(d_length/dx), int(d_width/dy),
     38                                                        len1=d_length, len2=d_width)
     39        domain = Domain(points, vertices, boundary)   
     40        domain.set_name('Test_Outlet_Inlet')                 # Output name
     41        domain.set_store()
     42        domain.set_default_order(2)
     43        domain.H0 = 0.01
     44        domain.tight_slope_limiters = 1
     45
     46        #print 'Size', len(domain)
     47
     48        #------------------------------------------------------------------------------
     49        # Setup initial conditions
     50        #------------------------------------------------------------------------------
     51
     52        def elevation(x, y):
     53            """Set up a elevation
     54            """
     55           
     56            z = numpy.zeros(x.shape,dtype='d')
     57            z[:] = elevation_0
     58           
     59            numpy.putmask(z, x > d_length/2, elevation_1)
     60   
     61            return z
     62           
     63        def stage(x,y):
     64            """Set up stage
     65            """
     66            z = numpy.zeros(x.shape,dtype='d')
     67            z[:] = stage_0
     68           
     69            numpy.putmask(z, x > d_length/2, stage_1)
     70
     71            return z
     72           
     73        #print 'Setting Quantities....'
     74        domain.set_quantity('elevation', elevation)  # Use function for elevation
     75        domain.set_quantity('stage',  stage)   # Use function for elevation
     76       
     77        return domain
    2778
    2879    def test_boyd_non_skew(self):
    … …  
    3788        elevation_0 = 10.0
    3889        elevation_1 = 10.0
     90
     91        domain_length = 200.0
     92        domain_width = 200.0
    3993
    4094        culvert_length = 20.0
    … …  
    46100        culvert_apron = 0.0
    47101        enquiry_gap = 10.0
    48        
    49         expected_Q = 4.55
    50         expected_v = 2.3
    51         expected_d = 0.54
    52        
    53 
    54         # Probably no need to change below here
    55        
    56         domain_length = 200.  #x-Dir
    57         domain_width  = 200.  #y-dir
    58         dx = dy = 5.0          # Resolution: Length of subdivisions on both axes
    59 
    60 
    61         points, vertices, boundary = rectangular_cross(int(domain_length/dx), int(domain_width/dy),
    62                                                         len1=domain_length, len2=domain_width)
    63         domain = Domain(points, vertices, boundary)   
    64         domain.set_name('Test_Outlet_Inlet')                 # Output name
    65         domain.set_default_order(2)
    66         domain.H0 = 0.01
    67         domain.tight_slope_limiters = 1
    68 
    69         print 'Size', len(domain)
    70 
    71         #------------------------------------------------------------------------------
    72         # Setup initial conditions
    73         #------------------------------------------------------------------------------
    74 
    75         def elevation(x, y):
    76             """Set up a elevation
    77             """
    78            
    79             z = numpy.zeros(x.shape,dtype='d')
    80             z[:] = elevation_0
    81            
    82             numpy.putmask(z, x > domain_length/2, elevation_1)
    83    
    84             return z
    85            
    86         def stage(x,y):
    87             """Set up stage
    88             """
    89             z = numpy.zeros(x.shape,dtype='d')
    90             z[:] = stage_0
    91            
    92             numpy.putmask(z, x > domain_length/2, stage_1)
    93 
    94             return z
    95            
    96         print 'Setting Quantities....'
    97         domain.set_quantity('elevation', elevation)  # Use function for elevation
    98         domain.set_quantity('stage',  stage)   # Use function for elevation
    99 
    100 
    101         print 'Defining Structures'
     102
     103       
     104        expected_Q = 6.23
     105        expected_v = 2.55
     106        expected_d = 0.66
     107       
     108
     109        domain = self._create_domain(d_length=domain_length,
     110                                     d_width=domain_width,
     111                                     dx = 5.0,
     112                                     dy = 5.0,
     113                                     elevation_0 = elevation_0,
     114                                     elevation_1 = elevation_1,
     115                                     stage_0 = stage_0,
     116                                     stage_1 = stage_1)
     117 
     118
     119        #print 'Defining Structures'
    102120       
    103121        ep0 = numpy.array([domain_length/2-culvert_length/2, 100.0])
    … …  
    118136                                    verbose=False)
    119137
    120         culvert.determine_inflow_outflow()
     138        #culvert.determine_inflow_outflow()
    121139       
    122140        ( Q, v, d ) = culvert.discharge_routine()
    123141       
    124         print 'test_boyd_non_skew'
    125         print 'Q: ', Q, 'expected_Q: ', expected_Q
     142        #print 'test_boyd_non_skew'
     143        #print 'Q: ', Q, 'expected_Q: ', expected_Q
    126144       
    127145
    … …  
    142160        elevation_0 = 10.0
    143161        elevation_1 = 10.0
     162
     163        domain_length = 200.0
     164        domain_width = 200.0
    144165
    145166        culvert_length = 20.0
    … …  
    152173        enquiry_gap = 10.0
    153174       
    154         expected_Q = 4.55
    155         expected_v = 2.3
    156         expected_d = 0.54
    157        
    158 
    159         # Probably no need to change below here
    160        
    161         domain_length = 200.  #x-Dir
    162         domain_width  = 200.  #y-dir
    163         dx = dy = 5.0          # Resolution: Length of subdivisions on both axes
    164 
    165 
    166         points, vertices, boundary = rectangular_cross(int(domain_length/dx), int(domain_width/dy),
    167                                                         len1=domain_length, len2=domain_width)
    168         domain = Domain(points, vertices, boundary)   
    169         domain.set_name('Test_Outlet_Inlet')                 # Output name
    170         domain.set_default_order(2)
    171         domain.H0 = 0.01
    172         domain.tight_slope_limiters = 1
    173 
    174         print 'Size', len(domain)
    175 
    176         #------------------------------------------------------------------------------
    177         # Setup initial conditions
    178         #------------------------------------------------------------------------------
    179 
    180         def elevation(x, y):
    181             """Set up a elevation
    182             """
    183            
    184             z = numpy.zeros(x.shape,dtype='d')
    185             z[:] = elevation_0
    186            
    187             numpy.putmask(z, x > domain_length/2, elevation_1)
    188    
    189             return z
    190            
    191         def stage(x,y):
    192             """Set up stage
    193             """
    194             z = numpy.zeros(x.shape,dtype='d')
    195             z[:] = stage_0
    196            
    197             numpy.putmask(z, x > domain_length/2, stage_1)
    198 
    199             return z
    200            
    201         print 'Setting Quantities....'
    202         domain.set_quantity('elevation', elevation)  # Use function for elevation
    203         domain.set_quantity('stage',  stage)   # Use function for elevation
    204 
    205 
    206         print 'Defining Structures'
     175        expected_Q = 6.23
     176        expected_v = 2.55
     177        expected_d = 0.66
     178       
     179
     180        domain = self._create_domain(d_length=domain_length,
     181                                     d_width=domain_width,
     182                                     dx = 5.0,
     183                                     dy = 5.0,
     184                                     elevation_0 = elevation_0,
     185                                     elevation_1 = elevation_1,
     186                                     stage_0 = stage_0,
     187                                     stage_1 = stage_1)
     188
     189        #print 'Defining Structures'
    207190       
    208191        a = domain_length/2 - culvert_length/2
    … …  
    225208                                    verbose=False)
    226209
    227         culvert.determine_inflow_outflow()
     210        #culvert.determine_inflow_outflow()
    228211       
    229212        ( Q, v, d ) = culvert.discharge_routine()
    230213       
    231         print 'test_boyd_skew'
    232         print 'Q: ', Q, 'expected_Q: ', expected_Q
     214        #print 'test_boyd_skew'
     215        #print 'Q: ', Q, 'expected_Q: ', expected_Q
    233216
    234217        assert numpy.allclose(Q, expected_Q, rtol=1.0e-2) #inflow
    … …  
    248231        elevation_0 = 10.0
    249232        elevation_1 = 10.0
     233
     234        domain_length = 200.0
     235        domain_width = 200.0
    250236
    251237        culvert_length = 20.0
    … …  
    258244        enquiry_gap = 10.0
    259245       
    260         expected_Q = 4.55
    261         expected_v = 2.3
    262         expected_d = 0.54
     246        expected_Q = 6.23
     247        expected_v = 2.55
     248        expected_d = 0.66
    263249       
    264250
    265251        # Probably no need to change below here
    266252       
    267         domain_length = 200.  #x-Dir
    268         domain_width  = 200.  #y-dir
    269         dx = dy = 5.0          # Resolution: Length of subdivisions on both axes
    270 
    271 
    272         points, vertices, boundary = rectangular_cross(int(domain_length/dx), int(domain_width/dy),
    273                                                         len1=domain_length, len2=domain_width)
    274         domain = Domain(points, vertices, boundary)   
    275         domain.set_name('Test_Outlet_Inlet')                 # Output name
    276         domain.set_default_order(2)
    277         domain.H0 = 0.01
    278         domain.tight_slope_limiters = 1
    279 
    280         print 'Size', len(domain)
    281 
    282         #------------------------------------------------------------------------------
    283         # Setup initial conditions
    284         #------------------------------------------------------------------------------
    285 
    286         def elevation(x, y):
    287             """Set up a elevation
    288             """
    289            
    290             z = numpy.zeros(x.shape,dtype='d')
    291             z[:] = elevation_0
    292            
    293             numpy.putmask(z, x > domain_length/2, elevation_1)
    294    
    295             return z
    296            
    297         def stage(x,y):
    298             """Set up stage
    299             """
    300             z = numpy.zeros(x.shape,dtype='d')
    301             z[:] = stage_0
    302            
    303             numpy.putmask(z, x > domain_length/2, stage_1)
    304 
    305             return z
    306            
    307         print 'Setting Quantities....'
    308         domain.set_quantity('elevation', elevation)  # Use function for elevation
    309         domain.set_quantity('stage',  stage)   # Use function for elevation
    310 
    311 
    312         print 'Defining Structures'
     253        domain = self._create_domain(d_length=domain_length,
     254                                     d_width=domain_width,
     255                                     dx = 5.0,
     256                                     dy = 5.0,
     257                                     elevation_0 = elevation_0,
     258                                     elevation_1 = elevation_1,
     259                                     stage_0 = stage_0,
     260                                     stage_1 = stage_1)
     261
     262
     263        #print 'Defining Structures'
    313264       
    314265        a = domain_length/2 - culvert_length/2
    … …  
    334285                                    verbose=False)
    335286
    336         culvert.determine_inflow_outflow()
     287        #culvert.determine_inflow_outflow()
    337288       
    338289        ( Q, v, d ) = culvert.discharge_routine()
    339290       
    340         print test_boyd_non_skew_enquiry_points
    341         print 'Q: ', Q, 'expected_Q: ', expected_Q
     291        #print 'test_boyd_non_skew_enquiry_points'
     292        #print 'Q: ', Q, 'expected_Q: ', expected_Q
     293        #print 'v: ', v, 'expected_v: ', expected_v
     294        #print 'd: ', d, 'expected_d: ', expected_d
    342295
    343296        assert numpy.allclose(Q, expected_Q, rtol=1.0e-2) #inflow

  • trunk/anuga_core/source/anuga/structures/testing_open_slot_wide_bridge.py

    r7980 r8073  
    154154                       end_point0=[40.0, 75.0],
    155155                       end_point1=[50.0, 75.0],
    156                        width=50.0,height=5.0,
     156                       width=50.0,depth=5.0,
    157157                       culvert_routine=boyd_generalised_culvert_model,       
    158158                       number_of_barrels=1,

  • trunk/anuga_core/source/anuga/structures/testing_outlet_control.py

    r8067 r8073  
    150150#                            end_point1=[50.0, 75.0],
    151151#                            width=50.0,
    152 #                            height=10.0,
     152#                            depth=10.0,
    153153#                            apron=5.0,
    154154#                            verbose=False)
    … …  
    164164    culvert_width = 50.0/number_of_culverts
    165165    y = 100-i*culvert_width - culvert_width/2.0
    166     ep0 = [40.0, y]
    167     ep1 = [50.0, y]
     166    ep0 = num.array([40.0, y])
     167    ep1 = num.array([50.0, y])
     168   
    168169    losses = {'inlet':0.5, 'outlet':1, 'bend':0, 'grate':0, 'pier': 0, 'other': 0}
    169170#    culverts.append(Boyd_pipe_operator(domain,
    … …  
    181182
    182183
    183     Boyd_pipe_operator(domain,
    184                             end_point0=ep0,
    185                             end_point1=ep1,
    186                             losses=losses,
    187                             diameter=1.5, #culvert_width, #3.658,
    188                             apron=6.0,
    189                             use_momentum_jet=True,
    190                             use_velocity_head=True,
    191                             manning=0.013,
    192                             logging=True,
    193                             label='pipe_culvert',
    194                             verbose=False)
     184
     185
    195186
    196187    Boyd_box_operator(domain,
    197                             end_point0=ep0,
    198                             end_point1=ep1,
     188                            end_points=[ep0,ep1],
    199189                            losses=losses,
    200190                            width=culvert_width,
    201                             height=10.0,
     191                            depth=10.0,
    202192                            apron=6.0,
    203193                            use_momentum_jet=True,
    … …  
    216206                            #losses,
    217207                            #width=25.0,
    218                             #height=10.0,
     208                            #depth=10.0,
    219209                            #apron=5.0,
    220210                            #manning=0.013,

  • trunk/anuga_core/source/anuga/structures/testing_wide_bridge.py

    r8008 r8073  
    171171#                            end_point1=[50.0, 75.0],
    172172#                            width=50.0,
    173 #                            height=10.0,
     173#                            depth=10.0,
    174174#                            apron=5.0,
    175175#                            verbose=False)
    … …  
    191191                            losses=1.5,
    192192                            width=3.658,
    193                             height=3.658,
     193                            depth=3.658,
    194194                            apron=5.0,
    195195                            use_momentum_jet=True,
    … …  
    206206#                            end_point1=[50.0, 62.5],
    207207#                            width=25.0,
    208 #                            height=10.0,
     208#                            depth=10.0,
    209209#                            apron=5.0,
    210210#                            manning=0.013,
    … …  
    219219    end_point0=[40.0, 75.0],
    220220    end_point1=[50.0, 75.0],
    221     width=50.0,height=5.0,
     221    width=50.0,depth=5.0,
    222222    culvert_routine=boyd_generalised_culvert_model,  #culvert_routine=weir_orifice_channel_culvert_model,     
    223223    number_of_barrels=1,

  • trunk/anuga_core/source/anuga/structures/testing_wide_bridge_old.py

    r7980 r8073  
    168168    end_point0=[40.0, 75.0],
    169169    end_point1=[50.0, 75.0],
    170     width=50.0,height=5.0,
     170    width=50.0,depth=5.0,
    171171    culvert_routine=boyd_generalised_culvert_model,  #culvert_routine=weir_orifice_channel_culvert_model,     
    172172    number_of_barrels=1,

  • trunk/anuga_core/source/anuga/utilities/numerical_tools.py

    r7690 r8073  
    251251        A: String.   Array of ASCII values (numpy can't handle this)
    252252
     253        A:None.      Return None
     254
    253255        typecode:    numeric type. If specified, use this in the conversion.
    254256                     If not, let numeric package decide.
    … …  
    264266#        msg = 'Sorry, cannot handle strings in ensure_numeric()'
    265267#        raise Exception, msg
     268
     269    if A is None:
     270        return None
    266271
    267272    if typecode is None:

  • trunk/anuga_core/source/anuga_parallel/run_parallel_sw_merimbula.py

    r8014 r8073  
    4747
    4848#mesh_filename = "merimbula_10785_1.tsh" ; x0 = 756000.0 ; x1 = 756500.0
    49 mesh_filename = "merimbula_43200.tsh"   ; x0 = 756000.0 ; x1 = 756500.0
    50 #mesh_filename = "test-100.tsh" ; x0 = 0.25 ; x1 = 0.5
     49#mesh_filename = "merimbula_43200.tsh"   ; x0 = 756000.0 ; x1 = 756500.0
     50mesh_filename = "test-100.tsh" ; x0 = 0.25 ; x1 = 0.5
    5151yieldstep = 5
    5252finaltime = 200

  • trunk/anuga_work/development/2010-projects/anuga_1d/sww/run_dry_dam.py

    r7884 r8073  
    7575
    7676
    77 finaltime = 4.0
    78 yieldstep = 0.1
     77finaltime = 20.0
     78yieldstep = 1.0
    7979L = 2000.0     # Length of channel (m)
    8080
    8181k = 0
    8282
    83 N = 3200
     83N = 10000
    8484print "Evaluating domain with %d cells" %N
    8585domain = Domain(*uniform_mesh(N,x_1=L))

  • trunk/anuga_work/development/2010-projects/anuga_1d/sww/run_parabolic_canal.py

    r7884 r8073  
    2828
    2929L_x = 2500.0     # Length of channel (m)
    30 N = 400         # Number of compuational cells
     30N = 1000         # Number of compuational cells
    3131cell_len = 4*L_x/N   # Origin = 0.0
    3232

  • trunk/anuga_work/development/2010-projects/kv/profiling.py

    r8051 r8073  
    1 from anuga.abstract_2d_finite_volumes.domain import Domain
     1from anuga.abstract_2d_finite_volumes.generic_domain import Generic_Domain
    22import anuga.abstract_2d_finite_volumes.mesh_factory as mesh_factory
    33import anuga.abstract_2d_finite_volumes.neighbour_mesh as neighbour_mesh
    … …  
    77import anuga.utilities.log as log
    88
    9 def evolve():
    10     print "Setting up"
    11     grid_rows = 40
    12     grid_cols = 40
    13     #Set up a rectangular grid
    14     points,elements,boundary = mesh_factory.rectangular_cross(grid_rows,grid_cols)
    15     mesh = neighbour_mesh.Mesh(points,elements)
    16     n = len(mesh)
    17     print "n =",n
    18     boundary_map = {}
    19     for (vol_id,edge) in boundary.keys():
    20         x = mesh.get_edge_midpoint_coordinate(vol_id,edge)[0]
    21         y = mesh.get_edge_midpoint_coordinate(vol_id,edge)[1]
    22         #Define boundary conditions: u=x^2-y^2, v=x-y+2
    23         boundary_map[(vol_id,edge)] = Dirichlet_boundary([1,x*x-y*y,x-y+2])
    24     #Now define an operator
    25     domain = Domain(source=points,triangles=elements,boundary=boundary_map)
    26     KV_Operator = Kinematic_Viscosity_Operator(domain)
    27     print "Applying stage heights"
    28     #Set up the operator and RHS
    29     h = num.ones((n,),num.float)
     9
     10print "Setting up"
     11grid_rows = 40
     12grid_cols = 40
     13#Set up a rectangular grid
     14points,elements,boundary = mesh_factory.rectangular_cross(grid_rows,grid_cols)
     15mesh = neighbour_mesh.Mesh(points,elements)
     16n = len(mesh)
     17print "n =",n
     18boundary_map = {}
     19for (vol_id,edge) in boundary.keys():
     20    x = mesh.get_edge_midpoint_coordinate(vol_id,edge)[0]
     21    y = mesh.get_edge_midpoint_coordinate(vol_id,edge)[1]
     22    #Define boundary conditions: u=x^2-y^2, v=x-y+2
     23    boundary_map[(vol_id,edge)] = Dirichlet_boundary([1,x*x-y*y,x-y+2])
     24#Now define an operator
     25domain = Generic_Domain(source=points,triangles=elements,boundary=boundary_map)
     26KV_Operator = Kinematic_Viscosity_Operator(domain)
     27print "Applying stage heights"
     28#Set up the operator and RHS
     29h = num.ones((n,),num.float)
     30   
     31def solve():   
    3032    KV_Operator.apply_stage_heights(h)
    3133    rhs = num.zeros((n,2),num.float) #we will solve for u and v
    … …  
    3840import profile, pstats
    3941FN = 'kv_profile.dat'
    40 print "Starting..."
    41 profile.run('evolve()', FN)
     42print "Starting solve..."
     43profile.run('solve()', FN)
    4244S = pstats.Stats(FN)
    4345s = S.sort_stats('cumulative').print_stats(30)
Note: See TracChangeset for help on using the changeset viewer.