Changeset 8375

Timestamp:

Mar 31, 2012, 2:47:12 PM (13 years ago)

Author:

davies

Message:

balanced_dev: fixed stack overflow problem, & new 'transect plot'
routine in util.py

Location:

trunk/anuga_work/development/gareth

Files:

• experimental/balanced_dev/swb2_domain.py (modified) (5 diffs)
• experimental/balanced_dev/swb2_domain_ext.c (modified) (6 diffs)
• tests/channel_floodplain/plotme.py (modified) (1 diff)
• tests/util.py (modified) (2 diffs)

trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain.py

@@ -144 +144 @@
         # self.check_integrity()
 
+        #print 'Into Evolve'
+
         msg = 'Attribute self.beta_w must be in the interval [0, 2]'
         assert 0 <= self.beta_w <= 2.0, msg
@@ -155 +157 @@
         if self.store is True and self.get_time() == self.get_starttime():
             self.initialise_storage()
-
+        #print 'Into Generic_Domain Evolve'
         # Call basic machinery from parent class
         for t in Generic_Domain.evolve(self, yieldstep=yieldstep,
                                        finaltime=finaltime, duration=duration,
                                        skip_initial_step=skip_initial_step):
+            #print 'Out of Generic_Domain Evolve'
             # Store model data, e.g. for subsequent visualisation
             if self.store is True:
@@ -208 +211 @@
         from swb2_domain_ext import compute_fluxes_ext_central \
                                       as compute_fluxes_ext
+
+        #print "."
 
         # Shortcuts
@@ -364 +369 @@
 
         # Remove very thin layers of water
+        #print 'Before protect'
         domain.protect_against_infinitesimal_and_negative_heights()
+        #print 'After protect'
 
         #for name in domain.conserved_quantities:
@@ -378 +385 @@
         #        raise Exception('Unknown order')
         
+        #print 'Before extrapolate'
         domain.extrapolate_second_order_edge_sw()
+        #print 'After extrapolate'
 
         #balance_deep_and_shallow(domain)

trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c

@@ -787 +787 @@
   double dqv[3], qmin, qmax, hmin, hmax, bedmax, stagemin;
   double hc, h0, h1, h2, beta_tmp, hfactor;
-  double xmom_centroid_store[number_of_elements], ymom_centroid_store[number_of_elements], dk, dv0, dv1, dv2, de[3];
-  double stage_centroid_store[number_of_elements];
-    
+  double dk, dv0, dv1, dv2, de[3];
+  //double xmom_centroid_store[number_of_elements], ymom_centroid_store[number_of_elements]
+  //double stage_centroid_store[number_of_elements];
+  
+  double *xmom_centroid_store, *ymom_centroid_store, *stage_centroid_store;
+
+  // Use malloc to avoid putting these variables on the stack, which can cause
+  // segfaults in large model runs
+  xmom_centroid_store = malloc(number_of_elements*sizeof(double));
+  ymom_centroid_store = malloc(number_of_elements*sizeof(double));
+  stage_centroid_store = malloc(number_of_elements*sizeof(double));
+ 
   if(extrapolate_velocity_second_order==1){
       // Replace momentum centroid with velocity centroid to allow velocity
@@ -819 +828 @@
 
     if (number_of_boundaries[k]==3)
+    //if (0==0)
     {
       // No neighbours, set gradient on the triangle to zero
@@ -1366 +1376 @@
       } 
   }
+
+  free(xmom_centroid_store);
+  free(ymom_centroid_store);
+  free(stage_centroid_store);
 
   return 0;
@@ -1859 +1873 @@
   number_of_elements = stage_centroid_values -> dimensions[0];  
 
+  //printf("In C before Extrapolate");
+  //e=1;
   // Call underlying computational routine
   e = _extrapolate_second_order_edge_sw(number_of_elements,
@@ -1884 +1900 @@
                    extrapolate_velocity_second_order);
 
+  //printf("In C before edges-to-vertices");
 
   //Extrapolate from edges to vertices
@@ -1896 +1913 @@
                                       (double *) elevation_vertex_values -> data, 
                                       extrapolate_velocity_second_order);
+  //printf("In C after edges-to-vertices");
 
   if (e == -1) {

trunk/anuga_work/development/gareth/tests/channel_floodplain/plotme.py

@@ -4 +4 @@
 
 # Time-index to plot outputs from
-index=75
+index=900
 
 #p2 = util.get_output('channel_floodplain1_bal_dev.sww', minimum_allowed_height=0.01)

trunk/anuga_work/development/gareth/tests/util.py

@@ -73 +73 @@
 ##############
 
-def animate_1D(time, var, x, ylab=' '): #, x=range(var.shape[1]), vmin=var.min(), vmax=var.max()):
-    # Input: time = one-dimensional time vector;
-    #        var =  array with first dimension = len(time) ;
-    #        x = (optional) vector width dimension equal to var.shape[1];
-    
-    import pylab
-    import numpy
-   
-    
-
-    pylab.close()
-    pylab.ion()
-
-    # Initial plot
-    vmin=var.min()
-    vmax=var.max()
-    line, = pylab.plot( (x.min(), x.max()), (vmin, vmax), 'o')
-
-    # Lots of plots
-    for i in range(len(time)):
-        line.set_xdata(x)
-        line.set_ydata(var[i,:])
-        pylab.draw()
-        pylab.xlabel('x')
-        pylab.ylabel(ylab)
-        pylab.title('time = ' + str(time[i]))
-    
-    return
 
 
@@ -198 +170 @@
 #pylab.quiver(x,y,xvel[n,:],yvel[n,:])
 #pylab.savefig('Plot1.png')
+
+def animate_1D(time, var, x, ylab=' '): #, x=range(var.shape[1]), vmin=var.min(), vmax=var.max()):
+    # Input: time = one-dimensional time vector;
+    #        var =  array with first dimension = len(time) ;
+    #        x = (optional) vector width dimension equal to var.shape[1];
+    
+    import pylab
+    import numpy
+   
+    
+
+    pylab.close()
+    pylab.ion()
+
+    # Initial plot
+    vmin=var.min()
+    vmax=var.max()
+    line, = pylab.plot( (x.min(), x.max()), (vmin, vmax), 'o')
+
+    # Lots of plots
+    for i in range(len(time)):
+        line.set_xdata(x)
+        line.set_ydata(var[i,:])
+        pylab.draw()
+        pylab.xlabel('x')
+        pylab.ylabel(ylab)
+        pylab.title('time = ' + str(time[i]))
+    
+    return
+
+def near_transect(p, point1, point2, tol=1.):
+    # Function to get the indices of points in p less than 'tol' from the line
+    # joining (x1,y1), and (x2,y2)
+    # p comes from util.get_output('mysww.sww')
+    #
+    # e.g.
+    # import util
+    # #import transect_interpolate
+    # from matplotlib import pyplot
+    # p=util.get_output('merewether_1m.sww',0.01)
+    # p2=util.get_centroids(p,velocity_extrapolation=True)
+    # #xxx=transect_interpolate.near_transect(p,[95., 85.], [120.,68.],tol=2.)
+    # xxx=util.near_transect(p,[95., 85.], [120.,68.],tol=2.)
+    # pyplot.scatter(xxx[1],p.vel[140,xxx[0]],color='red')
+
+    x1=point1[0]
+    y1=point1[1]
+
+    x2=point2[0]
+    y2=point2[1]
+
+    # Find line equation a*x + b*y + c = 0
+    # based on y=gradient*x +intercept
+    if x1!=x2:
+        gradient= (y2-y1)/(x2-x1)
+        intercept = y1 - gradient*x1
+
+        a = -gradient
+        b = 1.
+        c = -intercept
+    else:
+        #print 'FIXME: Still need to treat 0 and infinite gradients'
+        #assert 0==1
+        a=1.
+        b=0.
+        c=-x2 
+
+    # Distance formula
+    inv_denom = 1./(a**2 + b**2)**0.5
+    distp = abs(p.x*a + p.y*b + c)*inv_denom
+
+    near_points = (distp<tol).nonzero()[0]
+    
+    # Now find a 'local' coordinate for the point, projected onto the line
+    # g1 = unit vector parallel to the line
+    # g2 = vector joining (x1,y1) and (p.x,p.y)
+    g1x = x2-x1 
+    g1y = y2-y1
+    g1_norm = (g1x**2 + g1y**2)**0.5
+    g1x=g1x/g1_norm
+    g1y=g1x/g1_norm
+
+    g2x = p.x[near_points] - x1
+    g2y = p.y[near_points] - y1
+    
+    # Dot product = projected distance == a local coordinate
+    local_coord = g1x*g2x + g1y*g2y
+
+    return near_points, local_coord