Changeset 4200

Timestamp:

Jan 31, 2007, 12:07:32 PM (18 years ago)

Author:

ole

Message:

Added diagnosticts about which triangle is responsible for the largest computed speed.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/domain.py (modified) (4 diffs)
• abstract_2d_finite_volumes/neighbour_mesh.py (modified) (1 diff)
• abstract_2d_finite_volumes/show_balanced_limiters.py (modified) (1 diff)
• shallow_water/shallow_water_domain.py (modified) (2 diffs)
• shallow_water/shallow_water_ext.c (modified) (6 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

@@ -4 +4 @@
 
    Copyright 2004
-   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
+   Ole Nielsen, Stephen Roberts, Duncan Gray
    Geoscience Australia
 """
 
-from Numeric import allclose
+from Numeric import allclose, argmax
 from anuga.config import epsilon
 
@@ -216 +216 @@
         self.already_computed_flux = zeros((N, 3), Int)
 
+        # Storage for maximal speeds computed for each triangle by compute_fluxes
+        # This is used for diagnostics only
+        self.max_speed = zeros(N, Float)
+
         if mesh_filename is not None:
             # If the mesh file passed any quantity values
@@ -546 +550 @@
         ##assert hasattr(self, 'boundary_objects')
 
-    def write_time(self):
-        print self.timestepping_statistics()
-
-
-    def timestepping_statistics(self):
+    def write_time(self, track_location=False):
+        print self.timestepping_statistics(track_location)
+
+
+    def timestepping_statistics(self, track_location=False):
         """Return string with time stepping statistics for printing or logging
+
+        Optional boolean keyword track_location decides whether to report location of smallest timestep. 
         """
 
@@ -568 +574 @@
                      self.max_timestep, self.number_of_steps,
                      self.number_of_first_order_steps)
+
+        if track_location is True:
+            msg += '\n'
+            # Find index of largest computed flux speed
+            i = argmax(self.max_speed)
+
+            x, y = self.get_centroid_coordinates()[i]
+
+            s = 'Triangle #%d with centroid (%.4f, %.4f) ' %(i, x, y)
+            s += 'had the largest computed speed: %.4f m/s' %(self.max_speed[i])
+
+            # FIXME (Ole): Maybe add all quantity values at vertices
+
+            msg += s
 
         return msg

anuga_core/source/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

@@ -616 +616 @@
             v1 = [x2-x1, y2-y1]
             v2 = [x0-x2, y0-y2]
-
             a0 = anglediff(v1, v0)
             a1 = anglediff(v2, v1)

anuga_core/source/anuga/abstract_2d_finite_volumes/show_balanced_limiters.py

@@ -82 +82 @@
 #Evolve
 for t in domain.evolve(yieldstep = 0.1, finaltime = 30):
-    domain.write_time()
+    domain.write_time(track_location=True)
     domain.write_boundary_statistics(['stage'],'left')
 

anuga_core/source/anuga/shallow_water/shallow_water_domain.py

@@ -780 +780 @@
         Xmom.explicit_update[k] = flux[1]
         Ymom.explicit_update[k] = flux[2]
+        domain.max_speed[k] = max_speed
 
 
@@ -832 +833 @@
     
     domain.timestep = compute_fluxes_ext(timestep, domain.epsilon, domain.g,
-                                     domain.neighbours,
-                                     domain.neighbour_edges,
-                                     domain.normals,
-                                     domain.edgelengths,
-                                     domain.radii,
-                                     domain.areas,
-                                     domain.tri_full_flag,
-                                     Stage.edge_values,
-                                     Xmom.edge_values,
-                                     Ymom.edge_values,
-                                     Bed.edge_values,
-                                     Stage.boundary_values,
-                                     Xmom.boundary_values,
-                                     Ymom.boundary_values,
-                                     Stage.explicit_update,
-                                     Xmom.explicit_update,
-                                     Ymom.explicit_update,
-                                     domain.already_computed_flux)
+                                         domain.neighbours,
+                                         domain.neighbour_edges,
+                                         domain.normals,
+                                         domain.edgelengths,
+                                         domain.radii,
+                                         domain.areas,
+                                         domain.tri_full_flag,
+                                         Stage.edge_values,
+                                         Xmom.edge_values,
+                                         Ymom.edge_values,
+                                         Bed.edge_values,
+                                         Stage.boundary_values,
+                                         Xmom.boundary_values,
+                                         Ymom.boundary_values,
+                                         Stage.explicit_update,
+                                         Xmom.explicit_update,
+                                         Ymom.explicit_update,
+                                         domain.already_computed_flux,
+                                         domain.max_speed)
 
 

anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -1217 +1217 @@
     *xmom_explicit_update,
     *ymom_explicit_update,
-    *already_computed_flux;//tracks whether the flux across an edge has already been computed
-
+    *already_computed_flux,//tracks whether the flux across an edge has already been computed
+    *max_speed_array; //Keeps track of max speeds for each triangle
 
   //Local variables
@@ -1231 +1231 @@
 
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOOO",
+  if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOOOO",
                         &timestep,
                         &epsilon,
@@ -1250 +1250 @@
                         &xmom_explicit_update,
                         &ymom_explicit_update,
-                        &already_computed_flux)) {
+                        &already_computed_flux,
+                        &max_speed_array)) {
     PyErr_SetString(PyExc_RuntimeError, "Input arguments failed");
     return NULL;
@@ -1261 +1262 @@
     ((double *) stage_explicit_update -> data)[k]=0.0;
     ((double *) xmom_explicit_update -> data)[k]=0.0;
-    ((double *) ymom_explicit_update -> data)[k]=0.0;
-  }
+    ((double *) ymom_explicit_update -> data)[k]=0.0;  
+  }
+  
   //Loop through neighbours and compute edge flux for each
   for (k=0; k<number_of_elements; k++) {
@@ -1327 +1329 @@
     }
     } // end for i
+    
     //Normalise by area and store for when all conserved
     //quantities get updated
@@ -1332 +1335 @@
     ((double *) xmom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
     ((double *) ymom_explicit_update -> data)[k] /= ((double *) areas -> data)[k];
+    
+   
+    //Keep track of maximal speeds
+    ((double *) max_speed_array -> data)[k] = max_speed;    
+    
   } //end for k
   return Py_BuildValue("d", timestep);