Changeset 8867

Timestamp:

May 16, 2013, 8:57:00 PM (11 years ago)

Author:

davies

Message:

Updated to bal_dev

Location:

trunk/anuga_work/development/gareth

Files:

• experimental/balanced_dev/swb2_domain_ext.c (modified) (18 diffs)
• tests/channel_floodplain/plotme.py (modified) (1 diff)
• tests/dam_break/plotme.py (modified) (1 diff)
• tests/more_plot_util.py (added)
• tests/runup_sinusoid/runup_sinusoid.py (modified) (2 diffs)
• tests/wave_runup/plotme.py (modified) (4 diffs)

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

@@ -108 +108 @@
                            double *edgeflux, double *max_speed, 
                            double *pressure_flux, double hc, 
-                           double hc_n, 
-                           double speed_max_last) 
+                           double hc_n, double speed_max_last, 
+                           double smooth) 
 {
 
@@ -129 +129 @@
   double s_min, s_max, soundspeed_left, soundspeed_right;
   double denom, inverse_denominator, z;
-  double uint, t1, t2, t3, min_speed;
+  double uint, t1, t2, t3, min_speed, tmp;
   // Workspace (allocate once, use many)
   static double q_left_rotated[3], q_right_rotated[3], flux_right[3], flux_left[3];
@@ -277 +277 @@
       //edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]);
 
+      // Standard smoothing term, scaled
+      //edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i])*smooth;
+
       // Adjustment to the scheme by Kurganov and Lin 2007 Communications in Computational
       // Physics 2:141-163
-      //uint = (s_max*q_right_rotated[i] - s_min*q_left_rotated[i] - (flux_right[i] - flux_left[i]))*inverse_denominator;
-      //t1 = (q_right_rotated[i] - uint);
-      //t2 = (-q_left_rotated[i] + uint);
-      //t3 = _minmod(t1,t2);
-      //edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]-t3);
+      uint = (s_max*q_right_rotated[i] - s_min*q_left_rotated[i] - (flux_right[i] - flux_left[i]))*inverse_denominator;
+      t1 = (q_right_rotated[i] - uint);
+      t2 = (-q_left_rotated[i] + uint);
+      t3 = _minmod(t1,t2);
+      edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]-t3);
+
+      // test this
+      //edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i])*pow(min(hc/(hc_n+epsilon), hc_n/(hc+epsilon)),1.0);
+    
+      // test this 
+      //tmp=min(fabs(q_right_rotated[i])/(fabs(q_left_rotated[i])+epsilon) , fabs(q_left_rotated[i])/(fabs(q_right_rotated[i])+epsilon));
+      //edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]-t3)*pow(tmp,1.0);
 
       //GD HACK: Here, we supress the 'smoothing' part of the flux -- like scaling diffusion by local wave speed
+      // This is a bad idea! E.g. oscillations in landslide wave runup case
+      // However, it does supress the growth of some wet-dry instabilities (e.g. PNG).
       //if(i!=2){
-      edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i])*
-                     (min(*max_speed/(max(speed_max_last,1.0e-100)),1.0));//(min(hc/h_left,min(hc_n/h_right,1.0)));
+      //edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i])*
+      //               (min(*max_speed/(max(speed_max_last,1.0e-100)),1.0));//(min(hc/h_left,min(hc_n/h_right,1.0)));
       //}
 
@@ -357 +369 @@
     double stage_edge_lim, outgoing_mass_edges, bedtop, bedbot, pressure_flux, hc, hc_n, tmp, tmp2;
     static long call = 1; // Static local variable flagging already computed flux
-    double speed_max_last, vol;
-
-    double *max_bed_edgevalue, *min_bed_edgevalue, *edgeflux_store, *pressuregrad_store;
-
-    max_bed_edgevalue = malloc(number_of_elements*sizeof(double));
+    double speed_max_last, vol, smooth;
+
+    double *count_wet_edges, *count_wet_neighbours, *edgeflux_store, *pressuregrad_store;
+
+    count_wet_neighbours = malloc(number_of_elements*sizeof(double));
+    count_wet_edges = malloc(number_of_elements*sizeof(double));
     edgeflux_store = malloc(number_of_elements*9*sizeof(double));
     pressuregrad_store = malloc(number_of_elements*3*sizeof(double));
@@ -379 +392 @@
     speed_max_last=0.0;
     for (k = 0; k < number_of_elements; k++){
-        max_bed_edgevalue[k] = max(bed_edge_values[3*k], 
-                                   max(bed_edge_values[3*k+1], bed_edge_values[3*k+2]));
+        
+        count_wet_edges[k] = 0.;
+        count_wet_neighbours[k] = 0.;
+        for (i =0; i<3; i++){
+            ki=3*k+i;
+            if(stage_edge_values[ki] > bed_edge_values[ki]+epsilon) count_wet_edges[k]+=1.0;
+            
+            n=neighbours[ki];
+            if(n<0 || stage_centroid_values[n] > bed_centroid_values[n]+h0+epsilon) count_wet_neighbours[k]+=1.0;
+        }
+
         speed_max_last=max(speed_max_last, max_speed_array[k]);
-
-        //if(stage_centroid_values[k]<bed_centroid_values[k]){
-        //    printf("Stage < bed");
-        //}
     }
-
-    //printf("SML: %.12lf, %.12lf, %.12lf \n", speed_max_last, fabs(-speed_max_last), max(speed_max_last*2.0, 0.5*speed_max_last)/2.0);
-
-    //printf("SPM: %lf \n", speed_max_last);
 
     // For all triangles
@@ -459 +473 @@
             //}
 
+            smooth=1.0;
+            //if((count_wet_edges[k]<=1.0 | count_wet_neighbours[k]<=1.0) & 
+            //    (n>=0 && (count_wet_edges[n]<=1.0 | count_wet_neighbours[n]<=1.0))){
+            //    smooth=0.0;
+            //}
+            
             // Edge flux computation (triangle k, edge i)
             _flux_function_central(ql, qr, zl, zr,
@@ -464 +484 @@
                     epsilon, h0, limiting_threshold, g,
                     edgeflux, &max_speed, &pressure_flux, hc, hc_n, 
-                    speed_max_last);
+                    speed_max_last, smooth);
 
            
@@ -478 +498 @@
 
             // Update triangle k with flux from edge i
-
             // bedslope_work contains all gravity related terms -- weighting of
             // conservative and non-conservative versions.
             if(hc> h0){
+            //if(stage_centroid_values[k] > count_wet_edges[k]){
               bedslope_work = length*(g*(hc*ql[0]-0.5*max(ql[0]-zl,0.)*(ql[0]-zl)) + pressure_flux);
-              //bedslope_work = length*(g*(hc*stage_centroid_values[k]-0.5*max(ql[0]-zl,0.)*(ql[0]-zl)) + pressure_flux);
+              //bedslope_work = length*(g*(hc*stage_centroid_values[k]*0.-
+              //                           0.5*max(stage_centroid_values[k]-zl,0.)*(stage_centroid_values[k]-zl)) 
+              //                        + pressure_flux);
               //bedslope_work+= (1.0-tmp)*length*g*hc*ql[0]; // Non-conservative water surface slope 
               pressuregrad_store[ki] = bedslope_work;
             }else{
+              // NOTE: When hc<h0, pressure_flux is entirely computed from the
+              // other side of the edge. Hence, we can't necessarily satisfy
+              // well-balancing by just assuming that the 2nd and 3rd terms in
+              // bedslope_work cancel. So, we force it here -- equivalent to
+              // using the water surface gravity term
               bedslope_work = length*(g*(hc*ql[0]-0.0*max(ql[0]-zl,0.)*(ql[0]-zl)) + 0.0*pressure_flux);
               //bedslope_work = length*(g*(hc*stage_centroid_values[k]-0.0*max(ql[0]-zl,0.)*(ql[0]-zl)) + 0.0*pressure_flux);
@@ -492 +519 @@
               pressuregrad_store[ki] = bedslope_work;
             }
-            
+
+            //if(count_wet_edges[k]<=1.0) pressuregrad_store[ki]=0.0;
+
             already_computed_flux[ki] = call; // #k Done
 
@@ -504 +533 @@
 
                 if(hc_n> h0){
+                //if(stage_centroid_values[n] > count_wet_edges[n]){
                   bedslope_work = length*(g*(hc_n*qr[0]-0.5*max(qr[0]-zr,0.)*(qr[0]-zr)) + pressure_flux);
-                  //bedslope_work = length*(g*(hc_n*stage_centroid_values[n]-0.5*max(qr[0]-zr,0.)*(qr[0]-zr)) + pressure_flux);
+                  //bedslope_work = length*(g*(hc_n*stage_centroid_values[n]*0.-0.5*max(stage_centroid_values[n]-zr,0.)*
+                  //                                (stage_centroid_values[n]-zr)) + pressure_flux);
                   //bedslope_work+= (1.-tmp)*length*g*hc_n*qr[0];  
                   pressuregrad_store[nm] = bedslope_work;
                 }else{
+                  // NOTE: When hc<h0, pressure_flux is entirely computed from the
+                  // other side of the edge. Hence, we can't necessarily satisfy
+                  // well-balancing by just assuming that the 2nd and 3rd terms in
+                  // bedslope_work cancel. So, we force it here -- equivalent to
+                  // using the water surface gravity term
                   bedslope_work = length*(g*(hc_n*qr[0]-0.0*max(qr[0]-zr,0.)*(qr[0]-zr)) + 0.0*pressure_flux);
                   //bedslope_work = length*(g*(hc_n*stage_centroid_values[n]-0.0*max(qr[0]-zr,0.)*(qr[0]-zr)) + 0.0*pressure_flux);
@@ -514 +550 @@
                   pressuregrad_store[nm] = bedslope_work;
                 }
+            
+                //if(count_wet_edges[n]<=1.0) pressuregrad_store[nm]=0.0;
 
                 already_computed_flux[nm] = call; // #n Done
             }
 
+            //if(n==576 | n==579) pressuregrad_store[nm]=0.;
             // Update timestep based on edge i and possibly neighbour n
             // GD HACK:
@@ -709 +748 @@
     //return -1;
 //
-    free(max_bed_edgevalue);
-    //free(min_bed_edgevalue);
+    free(count_wet_edges);
+    free(count_wet_neighbours);
     free(edgeflux_store);
     free(pressuregrad_store);
@@ -739 +778 @@
   // distance between the bed_centroid_value and the max bed_edge_value of
   // every triangle.
-  double minimum_relative_height=0.0; 
+  double minimum_relative_height=0.1; 
   int mass_added=0;
 
@@ -747 +786 @@
       hc = wc[k] - zc[k];
       // Definine the maximum bed edge value on triangle k.
-      //bmax = 0.5*max((zv[3*k]+zv[3*k+1]),max((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
+      bmax = 0.5*max((zv[3*k]+zv[3*k+1]),max((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k])));
 
       //if (hc < max(minimum_relative_height*(bmax-zc[k]), minimum_allowed_height)) {
@@ -755 +794 @@
 
 
-      //if (hc < max(minimum_relative_height*(bmax-zc[k]), minimum_allowed_height)) {
-      if (hc < minimum_allowed_height*2.0) {
+      //if (hc < minimum_allowed_height*2.0 | hc < minimum_relative_height*(bmax-zc[k]) ){
+      if (hc < minimum_allowed_height*2.0 ){
             // Set momentum to zero and ensure h is non negative
-            //xmomc[k] = 0.;//xmomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
-            //ymomc[k] = 0.;//ymomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
+            xmomc[k] = 0.;//xmomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
+            ymomc[k] = 0.;//ymomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
             //xmomc[k] = xmomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
             //ymomc[k] = ymomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0;
-            xmomc[k] = xmomc[k]*max(hc-minimum_allowed_height,0.)/(minimum_allowed_height);// 0.0;
-            ymomc[k] = ymomc[k]*max(hc-minimum_allowed_height,0.)/(minimum_allowed_height);// 0.0;
+            //xmomc[k] = xmomc[k]*max(hc-minimum_allowed_height,0.)/(minimum_allowed_height);// 0.0;
+            //ymomc[k] = ymomc[k]*max(hc-minimum_allowed_height,0.)/(minimum_allowed_height);// 0.0;
 
         if (hc <= 0.0){
@@ -1965 +2004 @@
   PyArrayObject *normal, *ql, *qr,  *edgeflux;
   double g, epsilon, max_speed, H0, zl, zr;
-  double h0, limiting_threshold, pressure_flux;
+  double h0, limiting_threshold, pressure_flux, smooth;
 
   if (!PyArg_ParseTuple(args, "OOOddOddd",
@@ -1981 +2020 @@
  
   pressure_flux = 0.0; // Store the water pressure related component of the flux 
+  smooth = 1.0 ; // term to scale diffusion in riemann solver
+
   _flux_function_central((double*) ql -> data, 
                          (double*) qr -> data, 
@@ -1994 +2035 @@
                          ((double*) normal -> data)[0],
                          ((double*) normal -> data)[1], 
-                         ((double*) normal -> data)[1] 
+                         ((double*) normal -> data)[1],
+             smooth 
              );
   

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

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

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

@@ -9 +9 @@
 p2_st=util.get_centroids(p_st)
 
-p_dev = util.get_output('dam_break_20130504_145918/dam_break.sww', 0.001)
+p_dev = util.get_output('dam_break_20130514_201704/dam_break.sww', 0.001)
 p2_dev=util.get_centroids(p_dev, velocity_extrapolation=True)
 

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

@@ -28 +28 @@
 # Define topography
 #------------------
-scale_me=1.0
-boundary_ws=-0.1
+
+### Pathological
+scale_me=100.0
+boundary_ws=-0.1999
 init_ws=-0.2
-bumpiness=10. # Higher = shorter wavelength oscillations in topography
+bumpiness=50. # Higher = shorter wavelength oscillations in topography
+tstep=0.002
+lasttime=1.1
+
+### Sensible
+#scale_me=1.0
+#boundary_ws=-0.1
+#init_ws=-0.2
+#bumpiness=50. # Higher = shorter wavelength oscillations in topography
+#tstep=0.2
+#lasttime=40.
+
 #domain.minimum_allowed_height=domain.minimum_allowed_height*scale_me # Seems needed to make the algorithms behave
 
@@ -82 +95 @@
 #------------------------------
 
-for t in domain.evolve(yieldstep=0.2,finaltime=120.00):
+for t in domain.evolve(yieldstep=tstep,finaltime=lasttime):
+#for t in domain.evolve(yieldstep=0.2,finaltime=60.):
     print domain.timestepping_statistics()
     #print domain.boundary_flux_integral

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

@@ -2 +2 @@
 import util
 from matplotlib import pyplot as pyplot
-import numpy
+import scipy, numpy
 
 p= util.get_output('runup_v2.sww')
@@ -39 +39 @@
 
     pyplot.plot(p2.x[v], p2.elev[v], '--', color='green')
-    pyplot.title('Time ' + str(round(t,2)))
+    pyplot.title('Stage (red is analytical): Time ' + str(round(t,2)))
     pyplot.ylim((-25,20))
     pyplot.savefig('figure_stage_'+str(i)+'.png')
@@ -52 +52 @@
         pyplot.plot(t220[:,0], t220[:,2],'-',color='red')
     pyplot.ylim((-25,20))
-    pyplot.title('Time ' + str(round(t,2)))
+    pyplot.title('Velocity (red is analytical): Time ' + str(round(t,2)))
     pyplot.savefig('figure_vel_'+str(i)+'.png')
     pyplot.clf()
@@ -59 +59 @@
 model_shore_x=p2.time*0.
 model_shore_u=p2.time*0.
+
+# Hacks to identify the location of the shoreline
+shoreline_depth_thresh1=0.01
+shoreline_depth_thresh2=0.4
+
 for i in range(len(model_shore_x)):
     # Compute index of shoreline
-    vtmp = p2.stage[i,:]>p2.elev+0.002
-    # FIXME -- introduced a bug here
+    vtmp = p2.stage[i,:]>p2.elev+shoreline_depth_thresh1
+    
     model_shore_x[i]=p2.x[vtmp].min()
-    #v_ind=p2.x[vtmp].argmin()
+    #v_ind=(p2.stage[i,vtmp]-p2.elev[vtmp]).argmin()
     # Store x and xvel
     #model_shore_x[i]=p2.x[vtmp[v_ind]]
     #model_shore_u[i]=p2.xvel[i,vtmp[v_ind]]
 
+    # Extract shoreline velocity. It is tricky, to avoid getting
+    # zero velocity points -- might be a better way to do it?
+    vtmp2= (p2.stage[i,:]>p2.elev+shoreline_depth_thresh1)*\
+           (p2.stage[i,:]<p2.elev+shoreline_depth_thresh2)
+    mloc=abs(p2.xvel[i,vtmp2]).argmax()
+    #print abs(p2.xvel[i,vtmp2]).max(), mloc
+    model_shore_u[i]=p2.xvel[i,vtmp2][mloc]
+    #vtmp2 = (p2.stage[i,:]>p2.elev+0.01)
+    #mloc= (vtmp2*(p2.stage[i,:]-p2.elev)).argmin()
+    #model_shore_u[i] = p2.xvel[i,mloc]
+
 pyplot.plot(p2.time, model_shore_x-200.,'-o')
 pyplot.plot(shoreline[:,0], shoreline[:,1],'-', color='red')
+pyplot.title('Shoreline position (where depth >'+str(shoreline_depth_thresh1)+')')
 pyplot.savefig('Shoreline_position.png')
 
 # Can plot velocity as well -- but notice that velocity at the shoreline
 # keeps being set to zero, so the model result is constant zeros
-# pyplot.plot(p2.time, model_shore_u,'-o')
-# pyplot.plot(shoreline[:,0], shoreline[:,2],'-',color='red')
+pyplot.clf()
+pyplot.plot(p2.time, model_shore_u,'-o')
+pyplot.plot(shoreline[:,0], shoreline[:,2],'-',color='red')
+pyplot.title('Shoreline velocity: Peak speed where depth >'\
+             + str(shoreline_depth_thresh1) + ' and depth < '\
+             + str(shoreline_depth_thresh2))
+pyplot.savefig('Shoreline_velocity.png')