Ignore:
Timestamp:
Jun 22, 2010, 5:30:32 PM (14 years ago)
Author:
steve
Message:

Added in some more c based limiters

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/2010-projects/anuga_1d/sww/sww_domain.py

    r7852 r7868  
    9090        self.quantities_to_be_stored = ['stage','xmomentum']
    9191
    92         self.__doc__ = 'shallow_water_domain'
     92        self.__doc__ = 'sww_domain'
    9393
    9494
     
    178178            yield(t)
    179179
    180     def initialise_storage(self):
    181         """Create and initialise self.writer object for storing data.
    182         Also, save x and bed elevation
    183         """
    184 
    185         import data_manager
    186 
    187         #Initialise writer
    188         self.writer = data_manager.get_dataobject(self, mode = 'w')
    189 
    190         #Store vertices and connectivity
    191         self.writer.store_connectivity()
    192 
    193 
    194     def store_timestep(self, name):
    195         """Store named quantity and time.
    196 
    197         Precondition:
    198            self.write has been initialised
    199         """
    200         self.writer.store_timestep(name)
     180
     181
    201182
    202183
     
    253234    u_C   = Velocity.centroid_values
    254235       
    255     #print id(h_C)
    256     #FIXME replace with numpy.where
    257     for i in range(N):
    258         h_C[i] = w_C[i] - z_C[i]                                               
    259         if h_C[i] <= 0.0:
    260             #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
    261             h_C[i] = 1.0e-15
    262             w_C[i] = z_C[i]
    263             uh_C[i] = 0.0
    264                
    265    
    266 ##    for i in range(len(h_C)):
    267 ##      if h_C[i] < epsilon:
    268 ##          u_C[i] = 0.0  #Could have been negative
    269 ##          h_C[i] = 0.0
    270 ##      else:
    271 
    272     u_C[:,]  = uh_C/(h_C + h0/h_C)
     236    #Calculate height (and fix negatives)better be non-negative!
     237    h_C[:] = w_C - z_C
     238    u_C[:]  = uh_C/(h_C + h0/h_C)
    273239       
    274240    for name in [ 'velocity', 'stage' ]:
     
    284250            raise 'Unknown order'
    285251
     252
    286253    stage_V = domain.quantities['stage'].vertex_values                 
    287254    bed_V   = domain.quantities['elevation'].vertex_values     
     
    291258
    292259    h_V[:,:]    = stage_V - bed_V
    293     for i in range(len(h_C)):
    294         for j in range(2):
    295             if h_V[i,j] < 0.0 :
    296                 #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
    297                 dh = h_V[i,j]
    298                 h_V[i,j] = 0.0
    299                 stage_V[i,j] = bed_V[i,j]
    300                 h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
    301                 stage_V[i,(j+1)%2] = stage_V[i,(j+1)%2] + dh
    302                
     260
     261#    # protect from edge values going negative
     262#    h_V[:,1] = numpy.where(h_V[:,0] < 0.0 , h_V[:,1]-h_V[:,0], h_V[:,1])
     263#    h_V[:,0] = numpy.where(h_V[:,0] < 0.0 , 0.0, h_V[:,0])
     264#
     265#    h_V[:,0] = numpy.where(h_V[:,1] < 0.0 , h_V[:,0]-h_V[:,1], h_V[:,0])
     266#    h_V[:,1] = numpy.where(h_V[:,1] < 0.0 , 0.0, h_V[:,1])
     267#
     268#
     269#    stage_V[:,:] = bed_V + h_V
    303270    xmom_V[:,:] = u_V * h_V
    304271   
Note: See TracChangeset for help on using the changeset viewer.