Changeset 8251 for trunk/anuga_work

Timestamp:

Nov 18, 2011, 7:55:23 PM (13 years ago)

Author:

steve

Message:

Tracking down NaN in 1d anuga code

Location:

trunk/anuga_work/development/2010-projects/anuga_1d

Files:

• base/generic_domain.py (modified) (3 diffs)
• sww/sww_domain.py (modified) (7 diffs)

trunk/anuga_work/development/2010-projects/anuga_1d/base/generic_domain.py

@@ -1059 +1059 @@
         """
 
+        import numpy as np
+        Stage      = self.quantities['stage']
+        print 'w_ex_update', np.any(np.isnan(Stage.explicit_update))
 
         # Compute fluxes across each element edge
         self.compute_fluxes()
 
+        print 'w_ex_update', np.any(np.isnan(Stage.explicit_update))
+
         # Update timestep to fit yieldstep and finaltime
-        self.update_timestep(yieldstep, finaltime)
+        self.update_timestep(yieldstep, finaltime)     
 
         # Update conserved quantities
@@ -1390 +1395 @@
 
 
+        Stage      = self.quantities['stage']
+        Xmom       = self.quantities['xmomentum']
+        Bed        = self.quantities['elevation']
+        Height     = self.quantities['height']
+        Velocity   = self.quantities['velocity']
+
+        #Arrays
+        w_C   = Stage.centroid_values
+        uh_C  = Xmom.centroid_values
+        z_C   = Bed.centroid_values
+        h_C   = Height.centroid_values
+        u_C   = Velocity.centroid_values
+
+        import numpy as np
+        print '== before forcing ======='
+        print 'w_C', np.any(np.isnan(w_C))
+        print 'uh_C', np.any(np.isnan(uh_C))
+        print 'z_C', np.any(np.isnan(z_C))
+        print 'h_C', np.any(np.isnan(h_C))
+        print 'u_C', np.any(np.isnan(u_C))
+        print 'w_ex_update', np.any(np.isnan(Stage.explicit_update))
+
         timestep = self.timestep
+
+
 
         #Compute forcing terms
         self.compute_forcing_terms()
+
+        print '==after forcing ======='
+        print 'w_C', np.any(np.isnan(w_C))
+        print 'uh_C', np.any(np.isnan(uh_C))
+        print 'z_C', np.any(np.isnan(z_C))
+        print 'h_C', np.any(np.isnan(h_C))
+        print 'u_C', np.any(np.isnan(u_C))
+        print 'w_ex_update', np.any(np.isnan(Stage.explicit_update))
 
         #Update conserved_quantities
@@ -1400 +1437 @@
             Q.update(timestep)
 
-
+        print '==after quantity update  ======='
+        print 'w_C', np.any(np.isnan(w_C))
+        print 'uh_C', np.any(np.isnan(uh_C))
+        print 'z_C', np.any(np.isnan(z_C))
+        print 'h_C', np.any(np.isnan(h_C))
+        print 'u_C', np.any(np.isnan(u_C))
+        print 'w_ex_update', np.any(np.isnan(Stage.explicit_update))
 
 if __name__ == "__main__":

trunk/anuga_work/development/2010-projects/anuga_1d/sww/sww_domain.py

@@ -90 +90 @@
         self.quantities_to_be_stored = ['stage','xmomentum']
 
-        self.__doc__ = 'sww_domain'
+        self.__doc__ = 'sww_domain'
 
         self.set_spatial_order(2)
@@ -166 +166 @@
 
 
-
-
     def distribute_to_vertices_and_edges(self):
         #Call correct module function
         #(either from this module or C-extension)
         distribute_to_vertices_and_edges(self)
+
         
     def evolve(self, yieldstep = None, finaltime = None, duration = None,
@@ -223 +222 @@
 
     import sys
+    import numpy as np
     from anuga_1d.config import epsilon, h0
         
@@ -240 +240 @@
     h_C   = Height.centroid_values
     u_C   = Velocity.centroid_values
+
+    w_V = domain.quantities['stage'].vertex_values
+    z_V   = domain.quantities['elevation'].vertex_values
+    h_V     = domain.quantities['height'].vertex_values
+    u_V     = domain.quantities['velocity'].vertex_values
+    uh_V  = domain.quantities['xmomentum'].vertex_values
+
+
+#    print 'w_C', np.any(np.isnan(w_C))
+#    print 'uh_C', np.any(np.isnan(uh_C))
+#    print 'z_C', np.any(np.isnan(z_C))
+#    print 'h_C', np.any(np.isnan(h_C))
+#    print 'u_C', np.any(np.isnan(u_C))
+#
+#    print 'w_V', np.any(np.isnan(w_V))
+#    print 'uh_V', np.any(np.isnan(uh_V))
+#    print 'z_V', np.any(np.isnan(z_V))
+#    print 'h_V', np.any(np.isnan(h_V))
+#    print 'u_V', np.any(np.isnan(u_V))
         
     #Calculate height (and fix negatives)better be non-negative!
@@ -245 +264 @@
     h_C[:] = w_C - z_C
 
-    #h0 = 1.0e-12
-    u_C[:]  = uh_C/(h_C + h0/(h_C + 1.0e-16))
+    u_C[:] = numpy.where(h_C <= 1.0e-14, 0.0, uh_C/h_C)
+    
 
     #print 'domain.order', domain.order
@@ -263 +282 @@
 
 
-    stage_V = domain.quantities['stage'].vertex_values                  
-    bed_V   = domain.quantities['elevation'].vertex_values      
-    h_V     = domain.quantities['height'].vertex_values
-    u_V     = domain.quantities['velocity'].vertex_values
-    xmom_V  = domain.quantities['xmomentum'].vertex_values
-
-    h_V[:,:]    = stage_V - bed_V
-
+#    print 'w_C', np.any(np.isnan(w_C))
+#    print 'uh_C', np.any(np.isnan(uh_C))
+#    print 'z_C', np.any(np.isnan(z_C))
+#    print 'h_C', np.any(np.isnan(h_C))
+#    print 'u_C', np.any(np.isnan(u_C))
+
+
+
+    h_V[:,:]    = w_V - z_V
 
     #print 'any' ,numpy.any( h_V[:,0] < 0.0)
@@ -306 +326 @@
 #
 #
-    stage_V[:,:] = bed_V + h_V
+    w_V[:,:] = z_V + h_V
     #bed_V[:,:] = stage_V - h_V
-    xmom_V[:,:] = u_V * h_V
+    uh_V[:,:] = u_V * h_V
+
+
+#    print 'w_V', np.any(np.isnan(w_V))
+#    print 'uh_V', np.any(np.isnan(uh_V))
+#    print 'z_V', np.any(np.isnan(z_V))
+#    print 'h_V', np.any(np.isnan(h_V))
+#    print 'u_V', np.any(np.isnan(u_V))
     
     return