Changeset 8994


Ignore:
Timestamp:
Sep 27, 2013, 10:36:23 AM (11 years ago)
Author:
steve
Message:

Changing inlets to take into account the apron.

Location:
trunk/anuga_core/source
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r8813 r8994  
    27572757    double dk, dv0, dv1, dv2;
    27582758
    2759     double *xmom_centroid_store, *ymom_centroid_store, *stage_centroid_store;
    2760 
    2761 
    2762     // Associate memory location of Domain varialbe with local aliases
     2759    double *xmom_centroid_store;
     2760    double *ymom_centroid_store;
     2761    //double *stage_centroid_store;
     2762
     2763
     2764    // Associate memory location of Domain varibles with local aliases
    27632765    number_of_elements     = D->number_of_elements;
    27642766    epsilon                = D->epsilon;
     
    28362838    xmom_centroid_store = malloc(number_of_elements*sizeof(double));
    28372839    ymom_centroid_store = malloc(number_of_elements*sizeof(double));
    2838     stage_centroid_store = malloc(number_of_elements*sizeof(double));
     2840    // stage_centroid_store = malloc(number_of_elements*sizeof(double));
    28392841
    28402842    if (extrapolate_velocity_second_order == 1) {
     
    33183320    free(xmom_centroid_store);
    33193321    free(ymom_centroid_store);
    3320     free(stage_centroid_store);
     3322    //free(stage_centroid_store);
    33213323
    33223324
  • trunk/anuga_core/source/anuga/structures/boyd_box_operator.py

    r8985 r8994  
    155155                print 'delta_total_energy ',self.delta_total_energy
    156156                print 'outlet_enquiry_depth ',self.outflow.get_enquiry_depth()
     157                print 'inflow_enquiry_depth ',self.inflow.get_enquiry_depth()
     158                print 'outlet_enquiry_speed ',self.outflow.get_enquiry_speed()
     159                print 'inflow_enquiry_speed ',self.inflow.get_enquiry_speed()
    157160                print 'sum_loss ',self.sum_loss
    158161                print 'manning ',self.manning
  • trunk/anuga_core/source/anuga/structures/inlet.py

    r8873 r8994  
    4545        else: # poly is a polygon
    4646
    47             self.triangle_indices = inside_polygon(domain_centroids, self.poly)
     47            tris_0 = line_intersect(vertex_coordinates, [self.poly[0],self.poly[1]])
     48            tris_1 = inside_polygon(domain_centroids, self.poly)
     49            #print 40*"="
     50            #print tris_0
     51            #print tris_1
     52            self.triangle_indices = num.union1d(tris_0, tris_1)
     53            #print self.triangle_indices
    4854
    4955        if len(self.triangle_indices) == 0:
    5056            msg = 'Inlet poly=%s ' % (self.poly)
    51             msg += 'No triangles intersecting line '
     57            msg += 'No triangle centroids intersecting poly '
    5258            raise Exception, msg
    5359
  • trunk/anuga_core/source/anuga/structures/inlet_enquiry.py

    r8876 r8994  
    5454        if self.enquiry_index in self.triangle_indices:
    5555            msg = 'Enquiry point %s' % (self.enquiry_pt)
    56             msg += 'is in an inlet triangle'
     56            msg += ' is in an inlet triangle'
    5757            import warnings
    5858            warnings.warn(msg)
  • trunk/anuga_core/source/anuga/structures/run_culvert.py

    r8361 r8994  
    4242domain = anuga.Domain(points, vertices, boundary)   
    4343domain.set_starttime(10)
    44 domain.set_name('Test_culvert')                 # Output name
     44domain.set_name()                 # Output name
    4545domain.set_default_order(2)
    4646#domain.set_beta(1.5)
  • trunk/anuga_core/source/anuga/structures/structure_operator.py

    r8950 r8994  
    101101
    102102        # Slots for recording current statistics
     103        self.accumulated_flow = 0.0
    103104        self.discharge = 0.0
    104105        self.velocity = 0.0
     106        self.outlet_depth = 0.0
    105107        self.delta_total_energy = 0.0
    106108        self.driving_energy = 0.0
     
    116118        self.inlets = []
    117119        line0 = self.exchange_lines[0] #self.inlet_lines[0]
     120        if self.apron is None:
     121            poly0 = line0
     122        else:
     123            offset = -self.apron*self.outward_vector_0
     124            #print line0
     125            #print offset
     126            poly0 = num.array([ line0[0], line0[1], line0[1]+offset, line0[0]+offset])
     127            #print poly0
    118128        if self.invert_elevations is None:
    119129            invert_elevation0 = None
     
    126136        self.inlets.append(inlet_enquiry.Inlet_enquiry(
    127137                           self.domain,
    128                            line0,
     138                           poly0,
    129139                           enquiry_point0,
    130140                           invert_elevation = invert_elevation0,
     
    132142                           verbose = self.verbose))
    133143
    134 
     144        tris_0 = self.inlets[0].triangle_indices
     145        #print tris_0
     146        #print self.domain.centroid_coordinates[tris_0]
    135147
    136148        line1 = self.exchange_lines[1]
     149        if self.apron is None:
     150            poly1 = line1
     151        else:
     152            offset = -self.apron*self.outward_vector_1
     153            #print line1
     154            #print offset
     155            poly1 = num.array([ line1[0], line1[1], line1[1]+offset, line1[0]+offset])
     156            #print poly1
     157           
    137158        if self.invert_elevations is None:
    138159            invert_elevation1 = None
     
    144165        self.inlets.append(inlet_enquiry.Inlet_enquiry(
    145166                           self.domain,
    146                            line1,
     167                           poly1,
    147168                           enquiry_point1,
    148169                           invert_elevation = invert_elevation1,
     
    150171                           verbose = self.verbose))
    151172
     173
     174        tris_1 = self.inlets[1].triangle_indices
     175        #print tris_1
     176        #print self.domain.centroid_coordinates[tris_1]       
     177       
    152178        self.set_logging(logging)
    153179
     
    166192        old_inflow_ymom = self.inflow.get_average_ymom()
    167193
    168 
    169         # Implement the update of flow over a timestep by
    170         # using a semi-implict update. This ensures that
    171         # the update does not create a negative depth
    172         if old_inflow_depth > 0.0 :
    173                 Q_star = Q/old_inflow_depth
    174         else:
    175                 Q_star = 0.0
    176 
    177         factor = 1.0/(1.0 + Q_star*timestep/self.inflow.get_area())
    178 
    179         new_inflow_depth = old_inflow_depth*factor
    180         new_inflow_xmom = old_inflow_xmom*factor
    181         new_inflow_ymom = old_inflow_ymom*factor
    182            
    183         self.inflow.set_depths(new_inflow_depth)
    184 
    185         #inflow.set_xmoms(Q/inflow.get_area())
    186         #inflow.set_ymoms(0.0)
    187 
    188         self.inflow.set_xmoms(new_inflow_xmom)
    189         self.inflow.set_ymoms(new_inflow_ymom)
    190 
    191         loss = (old_inflow_depth - new_inflow_depth)*self.inflow.get_area()
    192 
    193         # set outflow
    194         if old_inflow_depth > 0.0 :
    195             timestep_star = timestep*new_inflow_depth/old_inflow_depth
    196         else:
    197             timestep_star = 0.0
    198 
    199 
    200 
    201 
    202         outflow_extra_depth = Q*timestep_star/self.outflow.get_area()
    203         outflow_direction = - self.outflow.outward_culvert_vector
    204         outflow_extra_momentum = outflow_extra_depth*barrel_speed*outflow_direction
    205            
    206         gain = outflow_extra_depth*self.outflow.get_area()
    207            
    208         #print Q, Q*timestep, barrel_speed, outlet_depth, Qstar, factor, timestep_star
    209         #print '  ', loss, gain
    210 
     194        semi_implicit = True
     195        if semi_implicit:   
     196            # Implement the update of flow over a timestep by
     197            # using a semi-implict update. This ensures that
     198            # the update does not create a negative depth
     199            if old_inflow_depth > 0.0 :
     200                    Q_star = Q/old_inflow_depth
     201            else:
     202                    Q_star = 0.0
     203   
     204            factor = 1.0/(1.0 + Q_star*timestep/self.inflow.get_area())
     205   
     206            new_inflow_depth = old_inflow_depth*factor
     207            new_inflow_xmom = old_inflow_xmom*factor
     208            new_inflow_ymom = old_inflow_ymom*factor
     209               
     210            self.inflow.set_depths(new_inflow_depth)
     211   
     212            #inflow.set_xmoms(Q/inflow.get_area())
     213            #inflow.set_ymoms(0.0)
     214   
     215            self.inflow.set_xmoms(new_inflow_xmom)
     216            self.inflow.set_ymoms(new_inflow_ymom)
     217   
     218            loss = (old_inflow_depth - new_inflow_depth)*self.inflow.get_area()
     219   
     220            # set outflow
     221            if old_inflow_depth > 0.0 :
     222                timestep_star = timestep*new_inflow_depth/old_inflow_depth
     223            else:
     224                timestep_star = 0.0
     225   
     226   
     227   
     228   
     229            outflow_extra_depth = Q*timestep_star/self.outflow.get_area()
     230            outflow_direction = - self.outflow.outward_culvert_vector
     231            outflow_extra_momentum = outflow_extra_depth*barrel_speed*outflow_direction
     232               
     233            gain = outflow_extra_depth*self.outflow.get_area()
     234           
     235            #print gain, loss
     236            assert num.allclose(gain-loss, 0.0)
     237               
     238            #print Q, Q*timestep, barrel_speed, outlet_depth, Qstar, factor, timestep_star
     239            #print '  ', loss, gain
     240        else:
     241
     242   
     243            factor = Q*timestep/self.inflow.get_area()
     244   
     245            new_inflow_depth = old_inflow_depth - factor
     246            new_inflow_xmom = old_inflow_xmom/old_inflow_depth*new_inflow_depth
     247            new_inflow_ymom = old_inflow_ymom/old_inflow_depth*new_inflow_depth
     248               
     249            self.inflow.set_depths(new_inflow_depth)   
     250            self.inflow.set_xmoms(new_inflow_xmom)
     251            self.inflow.set_ymoms(new_inflow_ymom)
     252   
     253   
     254            loss = (old_inflow_depth - new_inflow_depth)*self.inflow.get_area()
     255   
     256
     257   
     258            outflow_extra_depth = Q*timestep/self.outflow.get_area()
     259            outflow_direction = - self.outflow.outward_culvert_vector
     260            outflow_extra_momentum = outflow_extra_depth*barrel_speed*outflow_direction
     261               
     262            gain = outflow_extra_depth*self.outflow.get_area()
     263           
     264            assert num.allclose(gain-loss, 0.0)
     265               
     266            #print Q, Q*timestep, barrel_speed, outlet_depth, Qstar, factor, timestep_star
     267            #print '  ', loss, gain
     268
     269   
    211270        # Stats
    212         self.discharge  = Q#outflow_extra_depth*self.outflow.get_area()/timestep
    213         self.velocity = barrel_speed#self.discharge/outlet_depth/self.width
     271       
     272        self.accumulated_flow += gain
     273        self.discharge  = outflow_extra_depth*self.outflow.get_area()/timestep #Q
     274        self.velocity =   barrel_speed # self.discharge/outlet_depth/self.width
     275        self.outlet_depth = outlet_depth
    214276
    215277        new_outflow_depth = self.outflow.get_average_depth() + outflow_extra_depth
     
    217279        if self.use_momentum_jet :
    218280            # FIXME (SR) Review momentum to account for possible hydraulic jumps at outlet
    219             #new_outflow_xmom = outflow.get_average_xmom() + outflow_extra_momentum[0]
    220             #new_outflow_ymom = outflow.get_average_ymom() + outflow_extra_momentum[1]
     281            #new_outflow_xmom = self.outflow.get_average_xmom() + outflow_extra_momentum[0]
     282            #new_outflow_ymom = self.outflow.get_average_ymom() + outflow_extra_momentum[1]
    221283            new_outflow_xmom = barrel_speed*new_outflow_depth*outflow_direction[0]
    222284            new_outflow_ymom = barrel_speed*new_outflow_depth*outflow_direction[1]
     
    396458        message += '--------------------------\n'
    397459        message += 'Type: %s\n' % self.structure_type
     460       
     461        message += 'inlets[0]_enquiry_depth [m]:  %.2f\n' %self.inlets[0].get_enquiry_depth()
     462        message += 'inlets[0]_enquiry_speed [m/s]:  %.2f\n' %self.inlets[0].get_enquiry_speed()
     463        message += 'inlets[0]_enquiry_stage [m]:  %.2f\n' %self.inlets[0].get_enquiry_stage()
     464        message += 'inlets[0]_enquiry_elevation [m]:  %.2f\n' %self.inlets[0].get_enquiry_elevation()
     465        message += 'inlets[0]_average_depth [m]:  %.2f\n' %self.inlets[0].get_average_depth()
     466        message += 'inlets[0]_average_speed [m/s]:  %.2f\n' %self.inlets[0].get_average_speed()
     467        message += 'inlets[0]_average_stage [m]:  %.2f\n' %self.inlets[0].get_average_stage()
     468        message += 'inlets[0]_average_elevation [m]:  %.2f\n' %self.inlets[0].get_average_elevation()
     469
     470        message += '\n'
     471       
     472        message += 'inlets[1]_enquiry_depth [m]:  %.2f\n' %self.inlets[1].get_enquiry_depth()
     473        message += 'inlets[1]_enquiry_speed [m/s]:  %.2f\n' %self.inlets[1].get_enquiry_speed()
     474        message += 'inlets[1]_enquiry_stage [m]:  %.2f\n' %self.inlets[1].get_enquiry_stage()
     475        message += 'inlets[1]_enquiry_elevation [m]:  %.2f\n' %self.inlets[1].get_enquiry_elevation()
     476
     477        message += 'inlets[1]_average_depth [m]:  %.2f\n' %self.inlets[1].get_average_depth()
     478        message += 'inlets[1]_average_speed [m/s]:  %.2f\n' %self.inlets[1].get_average_speed()
     479        message += 'inlets[1]_average_stage [m]:  %.2f\n' %self.inlets[1].get_average_stage()
     480        message += 'inlets[1]_average_elevation [m]:  %.2f\n' %self.inlets[1].get_average_elevation()
     481
     482       
    398483        message += 'Discharge [m^3/s]: %.2f\n' % self.discharge
    399484        message += 'Velocity  [m/s]: %.2f\n' % self.velocity
     485        message += 'Outlet Depth  [m]: %.2f\n' % self.outlet_depth
     486        message += 'Accumulated Flow [m^3]: %.2f\n' % self.accumulated_flow
    400487        message += 'Inlet Driving Energy %.2f\n' % self.driving_energy
    401488        message += 'Delta Total Energy %.2f\n' % self.delta_total_energy
    402489        message += 'Control at this instant: %s\n' % self.case
     490       
     491
     492
    403493
    404494        print message
     
    423513        message += '%.5f, ' % self.discharge
    424514        message += '%.5f, ' % self.velocity
     515        message += '%.5f, ' % self.accumulated_flow
    425516        message += '%.5f, ' % self.driving_energy
    426517        message += '%.5f' % self.delta_total_energy
     518       
     519
    427520
    428521        return message
  • trunk/anuga_core/source/anuga_parallel/run_parallel_gate_operator.py

    r8877 r8994  
    22import sys
    33
    4 from anuga.utilities.system_tools import get_pathname_from_package
    5 from anuga.geometry.polygon_function import Polygon_function
     4#from anuga.utilities.system_tools import get_pathname_from_package
     5#from anuga.geometry.polygon_function import Polygon_function
    66       
    7 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
    8 from anuga.abstract_2d_finite_volumes.quantity import Quantity
     7#from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
     8#from anuga.abstract_2d_finite_volumes.quantity import Quantity
    99
    1010import anuga
     
    1616import numpy as num
    1717
    18 from anuga_parallel import distribute, myid, numprocs, finalize, barrier
     18from anuga import distribute, myid, numprocs, finalize, barrier
    1919
    20 from parallel_operator_factory import Inlet_operator, Boyd_box_operator
     20from anuga import Inlet_operator, Boyd_box_operator
    2121
    2222
     
    5656    return z
    5757
    58 
    59 """test_that_culvert_runs_rating
    60 
    61 This test exercises the culvert and checks values outside rating curve
    62 are dealt with       
    63 """
    64 
    65 path = get_pathname_from_package('anuga.culvert_flows')   
     58 
    6659
    6760length = 40.
     
    7366if myid == 0:
    7467
    75     points, vertices, boundary = rectangular_cross(int(length/dx),
    76                                                int(width/dy),
    77                                                len1=length,
    78                                                len2=width)
     68    points, vertices, boundary = anuga.rectangular_cross(int(length/dx),
     69                                                         int(width/dy),
     70                                                         len1=length,
     71                                                         len2=width)
    7972    domain = anuga.Domain(points, vertices, boundary)
    80     domain.set_name('run_gate_operator')                 # Output name
     73    domain.set_name()                 # Output name
    8174    domain.set_flow_algorithm('2_0')
    8275
Note: See TracChangeset for help on using the changeset viewer.