Ignore:
Timestamp:
Apr 1, 2009, 8:41:37 PM (15 years ago)
Author:
sudi
Message:

Revised codes for discontinuous bed.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/anuga_1d/shallow_water_domain_suggestion2.py

    r5827 r6694  
    5353
    5454        conserved_quantities = ['stage', 'xmomentum']
    55         other_quantities = ['elevation', 'friction', 'height', 'velocity']
     55        evolved_quantities = ['stage', 'xmomentum', 'elevation', 'height', 'velocity']
     56        other_quantities = ['friction']
    5657        Generic_Domain.__init__(self, coordinates, boundary,
    57                                 conserved_quantities, other_quantities,
     58                                conserved_quantities, evolved_quantities, other_quantities,
    5859                                tagged_elements)
    5960       
     
    9192        self.quantities_to_be_stored = ['stage','xmomentum']
    9293
    93         self.__doc__ = 'shallow_water_domain'
     94        self.__doc__ = 'shallow_water_domain_suggestion2'
     95        self.check_integrity()
    9496
    9597
     
    148150        msg = 'Second conserved quantity must be "xmomentum"'
    149151        assert self.conserved_quantities[1] == 'xmomentum', msg
     152               
     153        msg = 'First evolved quantity must be "stage"'
     154        assert self.evolved_quantities[0] == 'stage', msg
     155        msg = 'Second evolved quantity must be "xmomentum"'
     156        assert self.evolved_quantities[1] == 'xmomentum', msg
     157        msg = 'Third evolved quantity must be "elevation"'
     158        assert self.evolved_quantities[2] == 'elevation', msg
     159        msg = 'Fourth evolved quantity must be "height"'
     160        assert self.evolved_quantities[3] == 'height', msg
     161        msg = 'Fifth evolved quantity must be "velocity"'
     162        assert self.evolved_quantities[4] == 'velocity', msg
    150163
    151164    def extrapolate_second_order_sw(self):
     
    168181        distribute_to_vertices_and_edges(self)
    169182       
    170     def evolve(self, yieldstep = None, finaltime = None, duration = None,
    171                skip_initial_step = False):
    172         """Specialisation of basic evolve method from parent class
    173         """
    174 
    175         #Call check integrity here rather than from user scripts
    176         #self.check_integrity()
    177 
    178         #msg = 'Parameter beta_h must be in the interval [0, 1)'
    179         #assert 0 <= self.beta_h < 1.0, msg
    180         #msg = 'Parameter beta_w must be in the interval [0, 1)'
    181         #assert 0 <= self.beta_w < 1.0, msg
    182 
    183 
    184         #Initial update of vertex and edge values before any storage
    185         #and or visualisation
    186        
    187         #self.distribute_to_vertices_and_edges() ?????????????????????????????????
    188        
    189         #Initialise real time viz if requested
    190         #if self.visualise is True and self.time == 0.0:
    191         #    if self.visualiser is None:
    192         #        self.initialise_visualiser()
    193         #
    194         #    self.visualiser.update_timer()
    195         #    self.visualiser.setup_all()
    196 
    197         #Store model data, e.g. for visualisation
    198         #if self.store is True and self.time == 0.0:
    199         #    self.initialise_storage()
    200         #    #print 'Storing results in ' + self.writer.filename
    201         #else:
    202         #    pass
    203         #    #print 'Results will not be stored.'
    204         #    #print 'To store results set domain.store = True'
    205         #    #FIXME: Diagnostic output should be controlled by
    206         #    # a 'verbose' flag living in domain (or in a parent class)
    207 
     183    def evolve(self, yieldstep = None, finaltime = None, duration = None, skip_initial_step = False):
    208184        #Call basic machinery from parent class
    209         for t in Generic_Domain.evolve(self, yieldstep, finaltime,duration,
    210                                        skip_initial_step):
    211             #Real time viz
    212         #    if self.visualise is True:
    213         #        self.visualiser.update_all()
    214         #        self.visualiser.update_timer()
    215 
    216 
    217             #Store model data, e.g. for subsequent visualisation
    218         #    if self.store is True:
    219         #        self.store_timestep(self.quantities_to_be_stored)
    220 
    221             #FIXME: Could maybe be taken from specified list
    222             #of 'store every step' quantities
    223 
    224             #Pass control on to outer loop for more specific actions
     185        for t in Generic_Domain.evolve(self, yieldstep, finaltime, duration, skip_initial_step):
    225186            yield(t)
    226187
     
    352313       
    353314
    354     #We have got   h   and    u   at vertex, then  the following is the calculation of fluxes!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     315    #We have got h and u at vertex, then the following is the calculation of fluxes!
    355316    soundspeed_left  = sqrt(g*h_left)
    356317    soundspeed_right = sqrt(g*h_right)
     
    656617    import sys
    657618    from Numeric import zeros, Float
    658     from config import epsilon, h0
    659619   
    660620    N = domain.number_of_elements
    661621    timestep = float(sys.maxint)
    662     #epsilon = domain.epsilon
     622    epsilon = domain.epsilon
    663623    g = domain.g
    664624    neighbours = domain.neighbours
     
    667627    areas = domain.areas
    668628   
     629    Stage    = domain.quantities['stage']
     630    Xmom     = domain.quantities['xmomentum']
     631    Bed      = domain.quantities['elevation']
     632    Height   = domain.quantities['height']
     633    Velocity = domain.quantities['velocity']
     634
     635
     636    stage_boundary_values = Stage.boundary_values
     637    xmom_boundary_values  = Xmom.boundary_values
     638    bed_boundary_values   = Bed.boundary_values
     639    height_boundary_values= Height.boundary_values
     640    vel_boundary_values   = Velocity.boundary_values
     641   
     642    stage_explicit_update = Stage.explicit_update
     643    xmom_explicit_update  = Xmom.explicit_update
     644    bed_explicit_values   = Bed.explicit_update
     645    height_explicit_values= Height.explicit_update
     646    vel_explicit_values   = Velocity.explicit_update
     647   
     648    max_speed_array = domain.max_speed_array
     649   
     650    domain.distribute_to_vertices_and_edges()
     651    domain.update_boundary()
     652   
     653    stage_V = Stage.vertex_values
     654    xmom_V  = Xmom.vertex_values
     655    bed_V   = Bed.vertex_values
     656    height_V= Height.vertex_values
     657    vel_V   = Velocity.vertex_values
     658
     659    number_of_elements = len(stage_V)
     660
     661    #from config import epsilon, h0, h_min
     662    # ##Check this:
     663    #for i in range(N):
     664    #    height_V[i] = stage_V[i] - bed_V[i]                                           
     665    #    if height_V[i] <= h_min:  #epsilon :
     666    #       height_V[i] = 0.0
     667    #       stage_V[i]  = bed_V[i]
     668    #       xmom_V[i]   = 0.0
     669    #       vel_V[i]    = 0.0
     670    #   else:
     671    #        vel_V[i]  = xmom_V[i]/(height_V[i] +  h0/height_V[i])
     672    # ##End of part to be checked.
     673
     674   
     675    #print 'neighbours=',neighbours
     676    #print 'neighbour_vertices',neighbour_vertices
     677    #print 'normals=',normals
     678    #print 'areas=',areas
     679    from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced
     680    domain.flux_timestep = compute_fluxes_ext_wellbalanced(timestep,
     681                                  epsilon,
     682                                  g,
     683                                                           
     684                                  neighbours,
     685                                  neighbour_vertices,
     686                                  normals,
     687                                  areas,
     688                                                           
     689                                  stage_V,
     690                                  xmom_V,
     691                                  bed_V,
     692                                  height_V,
     693                                  vel_V,
     694                                                           
     695                                  stage_boundary_values,
     696                                  xmom_boundary_values,
     697                                  bed_boundary_values,
     698                                  height_boundary_values,
     699                                  vel_boundary_values,
     700                                                           
     701                                  stage_explicit_update,
     702                                  xmom_explicit_update,
     703                                  bed_explicit_values,
     704                                  height_explicit_values,
     705                                  vel_explicit_values,
     706                                                           
     707                                  number_of_elements,
     708                                  max_speed_array)
     709 
     710
     711# ###################################
     712
     713
     714
     715
     716
     717
     718# Module functions for gradient limiting (distribute_to_vertices_and_edges)
     719
     720def distribute_to_vertices_and_edges(domain):
     721    """Distribution from centroids to vertices specific to the
     722    shallow water wave
     723    equation.
     724
     725    It will ensure that h (w-z) is always non-negative even in the
     726    presence of steep bed-slopes by taking a weighted average between shallow
     727    and deep cases.
     728
     729    In addition, all conserved quantities get distributed as per either a
     730    constant (order==1) or a piecewise linear function (order==2).
     731
     732    FIXME: more explanation about removal of artificial variability etc
     733
     734    Precondition:
     735      All quantities defined at centroids and bed elevation defined at
     736      vertices.
     737
     738    Postcondition
     739      Conserved quantities defined at vertices
     740
     741    """
     742
     743    #from config import optimised_gradient_limiter
     744
     745    #Remove very thin layers of water
     746    #protect_against_infinitesimal_and_negative_heights(domain) 
     747
     748    import sys
     749    from Numeric import zeros, Float
     750    from config import epsilon, h0, h_min
     751
     752    N = domain.number_of_elements
     753
     754    #Shortcuts
    669755    Stage = domain.quantities['stage']
    670     Xmom  = domain.quantities['xmomentum']
    671     Bed   = domain.quantities['elevation']
     756    Xmom = domain.quantities['xmomentum']
     757    Bed = domain.quantities['elevation']
    672758    Height = domain.quantities['height']
    673759    Velocity = domain.quantities['velocity']
    674760
     761    #Arrays   
    675762    w_C   = Stage.centroid_values   
    676763    uh_C  = Xmom.centroid_values   
    677     z_C   = Bed.centroid_values
     764    z_C   = Bed.centroid_values 
    678765    h_C   = Height.centroid_values
    679766    u_C   = Velocity.centroid_values
     
    681768    for i in range(N):
    682769        h_C[i] = w_C[i] - z_C[i]                                               
    683         if h_C[i] < epsilon :
     770        if h_C[i] <= h_min:  #epsilon :
    684771            h_C[i] = 0.0
    685772            w_C[i] = z_C[i]
    686773            uh_C[i] = 0.0
    687 
    688     for i in range(len(h_C)):
    689         if h_C[i] < epsilon:
    690             u_C[i] = 0.0  #Could have been negative
    691             h_C[i] = 0.0
     774            u_C[i] = 0.0
    692775        else:
    693             u_C[i]  = uh_C[i]/(h_C[i] +  h0/h_C[i])
    694 
    695     for name in ['height', 'velocity']:
     776            u_C[i]  = uh_C[i]/(h_C[i] +  h0/h_C[i])
     777
     778    for name in ['stage', 'height', 'velocity']:
    696779        Q = domain.quantities[name]
    697780        if domain.order == 1:
     
    705788            raise 'Unknown order'
    706789
    707     w_V = domain.quantities['stage'].vertex_values                     
    708     h_V = domain.quantities['height'].vertex_values
    709     u_V = domain.quantities['velocity'].vertex_values           
    710     uh_V = domain.quantities['xmomentum'].vertex_values
    711     z_V   = w_V - h_V
    712 
     790    stage_V = domain.quantities['stage'].vertex_values
     791    bed_V   = domain.quantities['elevation'].vertex_values
     792    h_V     = domain.quantities['height'].vertex_values
     793    u_V     = domain.quantities['velocity'].vertex_values               
     794    xmom_V  = domain.quantities['xmomentum'].vertex_values
     795
     796    bed_V[:] = stage_V - h_V
     797    xmom_V[:] = u_V * h_V
    713798   
    714     stage_boundary_values = Stage.boundary_values
    715     xmom_boundary_values =  Xmom.boundary_values
    716     stage_explicit_update = Stage.explicit_update
    717     xmom_explicit_update = Xmom.explicit_update
    718     max_speed_array = domain.max_speed_array
    719     #number_of_elements = len(w_V)
     799    return
     800
    720801   
    721     #domain.distribute_to_vertices_and_edges()
    722     #domain.update_boundary()
    723     #stage_V = Stage.vertex_values
    724     #xmom_V  = Xmom.vertex_values
    725     #bed_V   = Bed.vertex_values
    726    
    727        
    728     from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced  #from comp_flux_ext import compute_fluxes_ext
    729        
    730     domain.flux_timestep = compute_fluxes_ext_wellbalanced(timestep,
    731                                   epsilon,
    732                                   g,
    733                                   neighbours,
    734                                   neighbour_vertices,
    735                                   normals,
    736                                   areas,
    737                                   w_V,
    738                                   uh_V,
    739                                   z_V,
    740                                   stage_boundary_values,
    741                                   xmom_boundary_values,
    742                                   stage_explicit_update,
    743                                   xmom_explicit_update,
    744                                   N,
    745                                   max_speed_array)
    746 
    747 # ###################################
    748 
    749 
    750 
    751 
    752 
    753 
    754 # Module functions for gradient limiting (distribute_to_vertices_and_edges)
    755 
    756 def distribute_to_vertices_and_edges(domain):
    757     """Distribution from centroids to vertices specific to the
    758     shallow water wave
    759     equation.
    760 
    761     It will ensure that h (w-z) is always non-negative even in the
    762     presence of steep bed-slopes by taking a weighted average between shallow
    763     and deep cases.
    764 
    765     In addition, all conserved quantities get distributed as per either a
    766     constant (order==1) or a piecewise linear function (order==2).
    767 
    768     FIXME: more explanation about removal of artificial variability etc
    769 
    770     Precondition:
    771       All quantities defined at centroids and bed elevation defined at
    772       vertices.
    773 
    774     Postcondition
    775       Conserved quantities defined at vertices
    776 
    777     """
    778 
    779     #from config import optimised_gradient_limiter
    780 
    781     #Remove very thin layers of water
    782     protect_against_infinitesimal_and_negative_heights(domain) 
    783 
    784     import sys
    785     from Numeric import zeros, Float
    786 
    787     N = domain.number_of_elements
    788 
    789     #Shortcuts
    790     Stage = domain.quantities['stage']
    791     Xmom = domain.quantities['xmomentum']
    792     Bed = domain.quantities['elevation']
    793     Height = domain.quantities['height']
    794     Velocity = domain.quantities['velocity']
    795 
    796     #Arrays   
    797     w_C   = Stage.centroid_values   
    798     uh_C  = Xmom.centroid_values   
    799     z_C   = Bed.centroid_values
    800     h_C   = Height.centroid_values
    801     u_C   = Velocity.centroid_values
    802        
    803     #print id(h_C)
    804     for i in range(N):
    805         h_C[i] = w_C[i] - z_C[i]                                                #This is h at centroids:  OK
    806    
    807     from config import epsilon, h0
    808        
    809     for i in range(len(h_C)):
    810         if h_C[i] < epsilon:
    811             u_C[i] = 0.0  #Could have been negative
    812             h_C[i] = 0.0
    813         else:
    814             u_C[i]  = uh_C[i]/(h_C[i] +  h0/h_C[i])
    815        
    816     for name in ['stage', 'velocity', 'height' ]:
    817         Q = domain.quantities[name]
    818         if domain.order == 1:
    819             Q.extrapolate_first_order()
    820         elif domain.order == 2:
    821             #print "add extrapolate second order to shallow water"
    822             #if name != 'height':
    823             Q.extrapolate_second_order()
    824             #Q.limit()
    825         else:
    826             raise 'Unknown order'
    827 
    828     stage  = domain.quantities['stage'].vertex_values                   #This is w at vertices:   OK
    829     h_V    = domain.quantities['height'].vertex_values
    830     bed    = stage - h_V
    831     u_V    = domain.quantities['velocity'].vertex_values               
    832     xmom_V = domain.quantities['xmomentum'].vertex_values       
    833        
    834     #h_V[:]    = stage - bed
    835     xmom_V[:] = u_V * h_V
    836     return stage, xmom_V, bed, h_V, u_V
    837802#
    838803def protect_against_infinitesimal_and_negative_heights(domain):
     
    10641029        #self.xmom    = domain.quantities['xmomentum'].edge_values
    10651030        #self.ymom    = domain.quantities['ymomentum'].edge_values
    1066         self.normals = domain.normals
    1067         self.stage   = domain.quantities['stage'].vertex_values
    1068         self.xmom    = domain.quantities['xmomentum'].vertex_values
     1031        self.normals  = domain.normals
     1032        self.stage    = domain.quantities['stage'].vertex_values
     1033        self.xmom     = domain.quantities['xmomentum'].vertex_values
     1034        self.bed      = domain.quantities['elevation'].vertex_values
     1035        self.height   = domain.quantities['height'].vertex_values
     1036        self.velocity = domain.quantities['velocity'].vertex_values
    10691037
    10701038        from Numeric import zeros, Float
    10711039        #self.conserved_quantities = zeros(3, Float)
    1072         self.conserved_quantities = zeros(2, Float)
     1040        #self.conserved_quantities = zeros(2, Float)
     1041        self.evolved_quantities = zeros(5, Float)
    10731042
    10741043    def __repr__(self):
    10751044        return 'Reflective_boundary'
    10761045
    1077 
     1046    """
    10781047    def evaluate(self, vol_id, edge_id):
    1079         """Reflective boundaries reverses the outward momentum
    1080         of the volume they serve.
    1081         """
    1082 
    1083         q = self.conserved_quantities
     1048        #Reflective boundaries reverses the outward momentum
     1049        #of the volume they serve.
     1050       
     1051        #q = self.conserved_quantities
     1052        q = self.evolved_quantities
    10841053        q[0] = self.stage[vol_id, edge_id]
    10851054        q[1] = self.xmom[vol_id, edge_id]
     
    11011070
    11021071        return q
     1072    """
     1073    def evaluate(self, vol_id, edge_id):
     1074        """Reflective boundaries reverses the outward momentum
     1075        of the volume they serve.
     1076        """
     1077        q = self.evolved_quantities
     1078        q[0] = self.stage[vol_id, edge_id]
     1079        q[1] = -self.xmom[vol_id, edge_id]
     1080        q[2] = self.bed[vol_id, edge_id]
     1081        q[3] = self.height[vol_id, edge_id]
     1082        q[4] = -self.velocity[vol_id, edge_id]
     1083        return q
     1084
     1085
     1086   
    11031087
    11041088class Dirichlet_boundary(Boundary):
Note: See TracChangeset for help on using the changeset viewer.