Changeset 8254


Ignore:
Timestamp:
Nov 22, 2011, 1:59:47 PM (12 years ago)
Author:
paul
Message:

Cleaned up sqpipe code

Location:
trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe
Files:
2 added
22 deleted
3 edited

Legend:

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

    r8238 r8254  
    107107    ztplot, = pylab.plot(x, z+t)
    108108
    109     plot1.set_ylim([-1,40])
     109    plot1.set_ylim([-1,50])
    110110    pylab.xlabel('Position')
    111111    pylab.ylabel('Stage')
  • trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_domain.py

    r8253 r8254  
    4646
    4747from anuga_1d.base.generic_domain import Generic_domain
     48from anuga_1d.sqpipe.sqpipe_forcing_terms import *
     49from anuga_1d.sqpipe.sqpipe_boundary_conditions import *
    4850
    4951#Shallow water domain
    50 class Sqpipe_domain(Generic_domain):
    51 
    52     def __init__(self, coordinates, conserved_quantities, forcing_terms = [], boundary = None, state = None, update_state_flag = True, bulk_modulus = 1.0, tagged_elements = None):
    53 
     52class Domain(Generic_domain):
     53
     54    def __init__(self, coordinates, boundary = None, forcing_terms = [], state = None, update_state_flag = True, bulk_modulus = 1.0, tagged_elements = None):
     55        conserved_quantities = ['area', 'discharge']
    5456        evolved_quantities = ['area', 'discharge', 'elevation', 'height', 'velocity','width','top','stage']
    5557        other_quantities = ['friction']
     
    149151
    150152        # Import flux method and call it
    151         from anuga_1d.sqpipe.sqpipe_dual_fixed_width_a_d_comp_flux import compute_fluxes_ext
     153        from anuga_1d.sqpipe.sqpipe_comp_flux import compute_fluxes_ext
    152154        self.flux_timestep = compute_fluxes_ext(timestep,self)       
    153155
     
    171173
    172174    def distribute_to_vertices_and_edges(self):
    173         # Should be implemented by the derived class
    174         raise Exception("Not implemented")
     175        # Shortcuts
     176        h0 = self.h0
     177        epsilon = self.epsilon
     178
     179        area      = self.quantities['area']
     180        discharge = self.quantities['discharge']
     181        bed       = self.quantities['elevation']
     182        height    = self.quantities['height']
     183        velocity  = self.quantities['velocity']
     184        width     = self.quantities['width']
     185        top       = self.quantities['top']
     186        stage     = self.quantities['stage']
     187
     188        #Arrays   
     189        a_C   = area.centroid_values
     190        d_C   = discharge.centroid_values
     191        z_C   = bed.centroid_values
     192        h_C   = height.centroid_values
     193        u_C   = velocity.centroid_values
     194        b_C   = width.centroid_values
     195        t_C   = top.centroid_values
     196        w_C   = stage.centroid_values
     197
     198        #Calculate height (and fix negatives)better be non-negative!
     199        h_C[:] = a_C/(b_C + h0/(b_C + h0)) # Do we need to protect against small b? Make a b0 instead of h0 or use epsilon?
     200        w_C[:] = z_C + h_C
     201        u_C[:]  = d_C/(h_C + h0/(h_C + h0))
     202
     203        for name in ['velocity', 'stage' ]:
     204            Q = self.quantities[name]
     205            if self.order == 1:
     206                Q.extrapolate_first_order()
     207            elif self.order == 2:
     208                Q.extrapolate_second_order()
     209            else:
     210                raise 'Unknown order'
     211
     212        # These have been extrapolated
     213        w_V = self.quantities['stage'].vertex_values
     214        u_V = self.quantities['velocity'].vertex_values
     215
     216        # These are the given geometry and remain fixed
     217        z_V  = bed.vertex_values
     218        b_V  = width.vertex_values
     219        t_V  = top.vertex_values
     220
     221        # These need to be updated
     222        a_V  = area.vertex_values
     223        d_V  = discharge.vertex_values
     224        h_V  = height.vertex_values
     225
     226        h_V[:,:]    = w_V - z_V
     227
     228        # Fix up small heights
     229        h_0 = numpy.where(h_V[:,0] < 0.0, 0.0, h_V[:,0])
     230        h_1 = numpy.where(h_V[:,0] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,1])
     231
     232        h_V[:,0] = h_0
     233        h_V[:,1] = h_1
     234
     235        h_0 = numpy.where(h_V[:,1] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,0])
     236        h_1 = numpy.where(h_V[:,1] < 0.0, 0.0, h_V[:,1])
     237
     238        h_V[:,0] = h_0
     239        h_V[:,1] = h_1
     240
     241        # It may still be possible that h is small
     242        # If we set h to zero, we should also set u to 0
     243        h_V[:,:] = numpy.where(h_V < epsilon, 0.0, h_V)
     244        u_V[:,:] = numpy.where(h_V < epsilon, 0.0, u_V)
     245
     246        # Reconstruct conserved quantities and make everything consistent
     247        # with new h values
     248        w_V[:,:] = z_V + h_V
     249        a_V[:,:] = h_V * b_V
     250        d_V[:,:] = u_V * a_V
     251
     252
     253        return
    175254       
    176255    def evolve(self, yieldstep = None, finaltime = None, duration = None,
  • trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/test_sqpipe_domain.py

    r8236 r8254  
    11import anuga_1d.sqpipe.parabolic_canal
    22import time
    3 import anuga_1d.sqpipe.sqpipe_test_domain as dom
     3import anuga_1d.sqpipe.sqpipe_domain as dom
    44
    55
Note: See TracChangeset for help on using the changeset viewer.