Changeset 8253


Ignore:
Timestamp:
Nov 21, 2011, 1:36:41 PM (12 years ago)
Author:
paul
Message:

sqpipe update_state now uses neighbouring states

Location:
trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_domain.py

    r8238 r8253  
    5050class Sqpipe_domain(Generic_domain):
    5151
    52     def __init__(self, coordinates, conserved_quantities, forcing_terms = [], boundary = None, state = None, update_state = True, bulk_modulus = 1.0, tagged_elements = None):
     52    def __init__(self, coordinates, conserved_quantities, forcing_terms = [], boundary = None, state = None, update_state_flag = True, bulk_modulus = 1.0, tagged_elements = None):
    5353
    5454        evolved_quantities = ['area', 'discharge', 'elevation', 'height', 'velocity','width','top','stage']
     
    7070            self.state = state           
    7171
    72         self.update_state = update_state
     72        self.update_state_flag = update_state_flag
    7373
    7474        # Get computational parameters
     
    126126            q = [q] # Turn argument into a list
    127127
    128         #Check correcness   
     128        #Check correctness   
    129129        for quantity_name in q:
    130130            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
     
    145145
    146146        # Update state before computing flux (if necessary)
    147         if self.update_state :
    148             h = self.quantities['height'].centroid_values
    149             t = self.quantities['top'].centroid_values
    150             self.state = numpy.where(h>t, 1, 0)
     147        if self.update_state_flag:
     148            self.update_state()
    151149
    152150        # Import flux method and call it
     
    154152        self.flux_timestep = compute_fluxes_ext(timestep,self)       
    155153
     154    def update_state(self):
     155        h = self.quantities['height'].centroid_values
     156        t = self.quantities['top'].centroid_values
     157
     158        new_state = numpy.zeros_like(self.state)
     159
     160        # Update boundary state
     161        new_state[0] = update_cell_state(0, [0, 1], h, t, self.state)
     162        N = self.number_of_elements-1
     163        new_state[N] = update_cell_state(N, [N-1, N], h, t, self.state)
     164
     165        # Update interior states
     166        for i in range(1, N-1):
     167            new_state[i] = update_cell_state(i, [i-1, i+1], h, t, self.state)
     168
     169        self.state = new_state
     170        #self.state = numpy.where(h>=t, 1, 0)
    156171
    157172    def distribute_to_vertices_and_edges(self):
     
    171186            yield(t)
    172187
     188# Auxillary methods
     189
     190# Method to update the state of a cell
     191# cell is the index of the cell to update
     192# neighbours is a list of the indices of it's neighbours
     193# (for boundary cells, this could include the cell itself)
     194def update_cell_state(cell, neighbours, h, t, s):
     195    # If height is bigger than top, then pressurise
     196    if h[cell] >= t[cell]:
     197        return 1
     198    else:
     199        # If the cell was pressurised, and all it's neighbours
     200        # are pressurised, it remains pressurised even though
     201        # height is less than top
     202        # Otherwise, it's free surface
     203        if all([s[i] for i in neighbours]):
     204            return s[cell]
     205        else:
     206            return 0
     207
    173208#=============== End of Shallow Water Domain ===============================
  • trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_test_domain.py

    r8238 r8253  
    1717                               state                = state,
    1818                               bulk_modulus         = bulk_modulus,
    19         #                       update_state         = False,
     19        #                       update_state_flag         = False,
    2020                               tagged_elements      = tagged_elements)
    2121
Note: See TracChangeset for help on using the changeset viewer.