Changeset 8251


Ignore:
Timestamp:
Nov 18, 2011, 7:55:23 PM (6 years ago)
Author:
steve
Message:

Tracking down NaN in 1d anuga code

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_core/source/anuga/version.py

    r6878 r8251  
    1 ## Subversion keywords:
    2        
    3 author_string = '$Author$'
    4 date_string = '$Date$'
    5 revision_string = '$LastChangedRevision$'
     1#! /usr/bin/python
    62
    7 #$LastChangedDate$
    8 #$LastChangedRevision$
    9 #$LastChangedBy$
     3# To change this template, choose Tools | Templates
     4# and open the template in the editor.
    105
    11 def get_version():
    12     """Extract version info from keyword strings
    13     """
    14    
    15     pass
     6__author__="steve"
     7__date__ ="$18/11/2011 11:35:22 AM$"
     8
     9if __name__ == "__main__":
     10    import anuga
     11    print 'Anuga version:',anuga.__version__
  • trunk/anuga_work/development/2010-projects/anuga_1d/base/generic_domain.py

    r7884 r8251  
    10591059        """
    10601060
     1061        import numpy as np
     1062        Stage      = self.quantities['stage']
     1063        print 'w_ex_update', np.any(np.isnan(Stage.explicit_update))
    10611064
    10621065        # Compute fluxes across each element edge
    10631066        self.compute_fluxes()
    10641067
     1068        print 'w_ex_update', np.any(np.isnan(Stage.explicit_update))
     1069
    10651070        # Update timestep to fit yieldstep and finaltime
    1066         self.update_timestep(yieldstep, finaltime)
     1071        self.update_timestep(yieldstep, finaltime)     
    10671072
    10681073        # Update conserved quantities
     
    13901395
    13911396
     1397        Stage      = self.quantities['stage']
     1398        Xmom       = self.quantities['xmomentum']
     1399        Bed        = self.quantities['elevation']
     1400        Height     = self.quantities['height']
     1401        Velocity   = self.quantities['velocity']
     1402
     1403        #Arrays
     1404        w_C   = Stage.centroid_values
     1405        uh_C  = Xmom.centroid_values
     1406        z_C   = Bed.centroid_values
     1407        h_C   = Height.centroid_values
     1408        u_C   = Velocity.centroid_values
     1409
     1410        import numpy as np
     1411        print '== before forcing ======='
     1412        print 'w_C', np.any(np.isnan(w_C))
     1413        print 'uh_C', np.any(np.isnan(uh_C))
     1414        print 'z_C', np.any(np.isnan(z_C))
     1415        print 'h_C', np.any(np.isnan(h_C))
     1416        print 'u_C', np.any(np.isnan(u_C))
     1417        print 'w_ex_update', np.any(np.isnan(Stage.explicit_update))
     1418
    13921419        timestep = self.timestep
     1420
     1421
    13931422
    13941423        #Compute forcing terms
    13951424        self.compute_forcing_terms()
     1425
     1426        print '==after forcing ======='
     1427        print 'w_C', np.any(np.isnan(w_C))
     1428        print 'uh_C', np.any(np.isnan(uh_C))
     1429        print 'z_C', np.any(np.isnan(z_C))
     1430        print 'h_C', np.any(np.isnan(h_C))
     1431        print 'u_C', np.any(np.isnan(u_C))
     1432        print 'w_ex_update', np.any(np.isnan(Stage.explicit_update))
    13961433
    13971434        #Update conserved_quantities
     
    14001437            Q.update(timestep)
    14011438
    1402 
     1439        print '==after quantity update  ======='
     1440        print 'w_C', np.any(np.isnan(w_C))
     1441        print 'uh_C', np.any(np.isnan(uh_C))
     1442        print 'z_C', np.any(np.isnan(z_C))
     1443        print 'h_C', np.any(np.isnan(h_C))
     1444        print 'u_C', np.any(np.isnan(u_C))
     1445        print 'w_ex_update', np.any(np.isnan(Stage.explicit_update))
    14031446
    14041447if __name__ == "__main__":
  • trunk/anuga_work/development/2010-projects/anuga_1d/sww/sww_domain.py

    r8225 r8251  
    9090        self.quantities_to_be_stored = ['stage','xmomentum']
    9191
    92         self.__doc__ = 'sww_domain'
     92        self.__doc__ = 'sww_domain'
    9393
    9494        self.set_spatial_order(2)
     
    166166
    167167
    168 
    169 
    170168    def distribute_to_vertices_and_edges(self):
    171169        #Call correct module function
    172170        #(either from this module or C-extension)
    173171        distribute_to_vertices_and_edges(self)
     172
    174173       
    175174    def evolve(self, yieldstep = None, finaltime = None, duration = None,
     
    223222
    224223    import sys
     224    import numpy as np
    225225    from anuga_1d.config import epsilon, h0
    226226       
     
    240240    h_C   = Height.centroid_values
    241241    u_C   = Velocity.centroid_values
     242
     243    w_V = domain.quantities['stage'].vertex_values
     244    z_V   = domain.quantities['elevation'].vertex_values
     245    h_V     = domain.quantities['height'].vertex_values
     246    u_V     = domain.quantities['velocity'].vertex_values
     247    uh_V  = domain.quantities['xmomentum'].vertex_values
     248
     249
     250#    print 'w_C', np.any(np.isnan(w_C))
     251#    print 'uh_C', np.any(np.isnan(uh_C))
     252#    print 'z_C', np.any(np.isnan(z_C))
     253#    print 'h_C', np.any(np.isnan(h_C))
     254#    print 'u_C', np.any(np.isnan(u_C))
     255#
     256#    print 'w_V', np.any(np.isnan(w_V))
     257#    print 'uh_V', np.any(np.isnan(uh_V))
     258#    print 'z_V', np.any(np.isnan(z_V))
     259#    print 'h_V', np.any(np.isnan(h_V))
     260#    print 'u_V', np.any(np.isnan(u_V))
    242261       
    243262    #Calculate height (and fix negatives)better be non-negative!
     
    245264    h_C[:] = w_C - z_C
    246265
    247     #h0 = 1.0e-12
    248     u_C[:]  = uh_C/(h_C + h0/(h_C + 1.0e-16))
     266    u_C[:] = numpy.where(h_C <= 1.0e-14, 0.0, uh_C/h_C)
     267   
    249268
    250269    #print 'domain.order', domain.order
     
    263282
    264283
    265     stage_V = domain.quantities['stage'].vertex_values                 
    266     bed_V   = domain.quantities['elevation'].vertex_values     
    267     h_V     = domain.quantities['height'].vertex_values
    268     u_V     = domain.quantities['velocity'].vertex_values
    269     xmom_V  = domain.quantities['xmomentum'].vertex_values
    270 
    271     h_V[:,:]    = stage_V - bed_V
    272 
     284#    print 'w_C', np.any(np.isnan(w_C))
     285#    print 'uh_C', np.any(np.isnan(uh_C))
     286#    print 'z_C', np.any(np.isnan(z_C))
     287#    print 'h_C', np.any(np.isnan(h_C))
     288#    print 'u_C', np.any(np.isnan(u_C))
     289
     290
     291
     292    h_V[:,:]    = w_V - z_V
    273293
    274294    #print 'any' ,numpy.any( h_V[:,0] < 0.0)
     
    306326#
    307327#
    308     stage_V[:,:] = bed_V + h_V
     328    w_V[:,:] = z_V + h_V
    309329    #bed_V[:,:] = stage_V - h_V
    310     xmom_V[:,:] = u_V * h_V
     330    uh_V[:,:] = u_V * h_V
     331
     332
     333#    print 'w_V', np.any(np.isnan(w_V))
     334#    print 'uh_V', np.any(np.isnan(uh_V))
     335#    print 'z_V', np.any(np.isnan(z_V))
     336#    print 'h_V', np.any(np.isnan(h_V))
     337#    print 'u_V', np.any(np.isnan(u_V))
    311338   
    312339    return
Note: See TracChangeset for help on using the changeset viewer.