Changeset 9148

Timestamp:

Jun 12, 2014, 7:55:33 AM (11 years ago)

Author:

steve

Message:

Working on validation tests

Location:

trunk/anuga_core/source/anuga

Files:

• config.py (modified) (1 diff)
• file/sww.py (modified) (2 diffs)
• shallow_water/swDE1_domain_ext.c (modified) (2 diffs)
• structures/boyd_box_operator.py (modified) (3 diffs)
• utilities/plot_utils.py (modified) (5 diffs)
• utilities/sww_merge.py (modified) (5 diffs)
• validation_utilities/parameters.py (modified) (1 diff)
• validation_utilities/run_validation.py (modified) (1 diff)

trunk/anuga_core/source/anuga/config.py

@@ -72 +72 @@
 tight_slope_limiters = True
 
-use_edge_limiter = False  # The edge limiter is better, but most runs have been
-                          # using vertex limiting. Validations passed with this
-                          # one True 9th May 2008, but many unit tests need
-                          # backward compatibility flag set FIXME(Ole).
+use_edge_limiter = False    # The edge limiter is better, but most runs have been
+                            # using vertex limiting. Validations passed with this
+                            # one True 9th May 2008, but many unit tests need
+                            # backward compatibility flag set FIXME(Ole).
 
 # Use centroid velocities to reconstruct momentum at vertices in

trunk/anuga_core/source/anuga/file/sww.py

@@ -7 +7 @@
 class NewQuantity(exceptions.Exception): pass
 class DataDomainError(exceptions.Exception): pass
-class DataMissingValuesError(exceptions.Exception): pass
 class DataTimeError(exceptions.Exception): pass
 
@@ -811 +810 @@
 
 
-         # dimension definitions
+        # dimension definitions
         #outfile.createDimension('number_of_volumes', number_of_volumes)
         #outfile.createDimension('number_of_vertices', 3)

trunk/anuga_core/source/anuga/shallow_water/swDE1_domain_ext.c

@@ -96 +96 @@
 //
 //}
+
+
+// Innermost flux function (using stage w=z+h)
+int _flux_function_toro(double *q_left, double *q_right,
+                           double h_left, double h_right,
+                           double hle, double hre,
+                           double n1, double n2,
+                           double epsilon,
+                           double ze,
+                           double limiting_threshold,
+                           double g,
+                           double *edgeflux, double *max_speed,
+                           double *pressure_flux, double hc,
+                           double hc_n, double speed_max_last)
+{
+
+  /*Compute fluxes between volumes for the shallow water wave equation
+    cast in terms of the 'stage', w = h+z using
+
+        HLL scheme of Fraccarollo and Toro Experimental and numerical assessment of the shallow
+water model for two-dimensional dam-break type. Journal of Computational Physics,
+33(6):843–864, 1995.
+
+    FIXME: Several variables in this interface are no longer used, clean up
+  */
+
+  int i;
+
+  double w_left,  uh_left, vh_left, u_left;
+  double w_right, uh_right, vh_right, u_right;
+  double s_min, s_max, soundspeed_left, soundspeed_right;
+  double u_m, h_m, soundspeed_m, s_m;
+  double denom, inverse_denominator;
+  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];
+
+  // Copy conserved quantities to protect from modification
+  q_left_rotated[0] = q_left[0];
+  q_right_rotated[0] = q_right[0];
+  q_left_rotated[1] = q_left[1];
+  q_right_rotated[1] = q_right[1];
+  q_left_rotated[2] = q_left[2];
+  q_right_rotated[2] = q_right[2];
+
+  // Align x- and y-momentum with x-axis
+  _rotate(q_left_rotated, n1, n2);
+  _rotate(q_right_rotated, n1, n2);
+
+
+  // Compute speeds in x-direction
+  //w_left = q_left_rotated[0];
+  uh_left=q_left_rotated[1];
+  vh_left=q_left_rotated[2];
+  if(hle>0.0){
+    u_left = uh_left/hle ; //max(h_left, 1.0e-06);
+    uh_left=h_left*u_left;
+    vh_left=h_left/(hle)*vh_left;
+  }else{
+    u_left=0.;
+    uh_left=0.;
+    vh_left=0.;
+  }
+
+  //u_left = _compute_speed(&uh_left, &hle,
+  //            epsilon, h0, limiting_threshold);
+
+  //w_right = q_right_rotated[0];
+  uh_right = q_right_rotated[1];
+  vh_right = q_right_rotated[2];
+  if(hre>0.0){
+    u_right = uh_right/hre;//max(h_right, 1.0e-06);
+    uh_right=h_right*u_right;
+    vh_right=h_right/hre*vh_right;
+  }else{
+    u_right=0.;
+    uh_right=0.;
+    vh_right=0.;
+  }
+  //u_right = _compute_speed(&uh_right, &hre,
+  //              epsilon, h0, limiting_threshold);
+
+
+
+  // Maximal and minimal wave speeds
+  soundspeed_left  = sqrt(g*h_left);
+  soundspeed_right = sqrt(g*h_right);
+
+  // Toro for shallow water
+  u_m = 0.5*(u_left + u_right) + soundspeed_left - soundspeed_right;
+  h_m = (u_left + 2.0*soundspeed_left - u_right + 2.0*soundspeed_right);
+  h_m = h_m*h_m/(16.0*g);
+  soundspeed_m = sqrt(g*h_m);
+
+
+  if(h_left < 1.0e-10){
+          s_min = u_right - 2.0*soundspeed_right;
+          s_max = u_right + soundspeed_right;
+          s_m   = s_min;
+  }
+  else if (h_right < 1.0e-10){
+          s_min = u_left - soundspeed_left;
+          s_max = u_left + 2.0*soundspeed_left;
+          s_m = s_max;
+  }
+  else {
+          s_max = max(u_right + soundspeed_right, u_m + soundspeed_right);
+          s_min = min(u_left - soundspeed_left, u_m - soundspeed_m);
+  }
+
+  if (s_max < 0.0)
+  {
+    s_max = 0.0;
+  }
+
+  if (s_min > 0.0)
+  {
+    s_min = 0.0;
+  }
+
+
+  // Flux formulas
+  flux_left[0] = u_left*h_left;
+  flux_left[1] = u_left*uh_left; //+ 0.5*g*h_left*h_left;
+  flux_left[2] = u_left*vh_left;
+
+  flux_right[0] = u_right*h_right;
+  flux_right[1] = u_right*uh_right ; //+ 0.5*g*h_right*h_right;
+  flux_right[2] = u_right*vh_right;
+
+  // Flux computation
+  denom = s_max - s_min;
+  if (denom < epsilon)
+  {
+    // Both wave speeds are very small
+    memset(edgeflux, 0, 3*sizeof(double));
+
+    *max_speed = 0.0;
+    //*pressure_flux = 0.0;
+    *pressure_flux = 0.5*g*0.5*(h_left*h_left+h_right*h_right);
+  }
+  else
+  {
+    // Maximal wavespeed
+    *max_speed = max(s_max, -s_min);
+
+    inverse_denominator = 1.0/max(denom,1.0e-100);
+    for (i = 0; i < 3; i++)
+    {
+      edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i];
+
+      // Standard smoothing term
+      // edgeflux[i] += 1.0*(s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]);
+      // Smoothing by stage alone can cause high velocities / slow draining for nearly dry cells
+      if(i==0) edgeflux[i] += (s_max*s_min)*(max(q_right_rotated[i],ze) - max(q_left_rotated[i],ze));
+      if(i==1) edgeflux[i] += (s_max*s_min)*(uh_right - uh_left);
+      if(i==2) edgeflux[i] += (s_max*s_min)*(vh_right - vh_left);
+
+      edgeflux[i] *= inverse_denominator;
+    }
+    // Separate pressure flux, so we can apply different wet-dry hacks to it
+    *pressure_flux = 0.5*g*( s_max*h_left*h_left -s_min*h_right*h_right)*inverse_denominator;
+
+
+    // Rotate back
+    _rotate(edgeflux, n1, -n2);
+  }
+
+  return 0;
+}
+
+
+
+
 
 // Innermost flux function (using stage w=z+h)
@@ -516 +690 @@
 
             // Edge flux computation (triangle k, edge i)
-            _flux_function_central(ql, qr, 
+            _flux_function_central(ql, qr,
                     h_left, h_right,
                     hle, hre,

trunk/anuga_core/source/anuga/structures/boyd_box_operator.py

@@ -36 +36 @@
                  manning=0.013,
                  enquiry_gap=0.0,
+                 smoothing_timescale=0.0,
                  use_momentum_jet=True,
                  use_velocity_head=True,
@@ -90 +91 @@
         self.case = 'N/A'
         
+        # May/June 2014 -- allow 'smoothing ' of driving_energy, delta total energy, and outflow_enq_depth
+        self.smoothing_timescale=0.
+        self.smooth_delta_total_energy=0.
+        self.smooth_Q=0.
+        # Set them based on a call to the discharge routine with smoothing_timescale=0.
+        # [values of self.smooth_* are required in discharge_routine, hence dummy values above]
+        Qvd=self.discharge_routine()
+        self.smooth_delta_total_energy=1.0*self.delta_total_energy
+        self.smooth_Q=Qvd[0]
+        # Finally, set the smoothing timescale we actually want
+        self.smoothing_timescale=smoothing_timescale
+
+        
         if verbose:
             print self.get_culvert_slope()
@@ -129 +143 @@
             self.outflow = self.inlets[0]
             self.delta_total_energy = -self.delta_total_energy
+            
+            
+            
 
 

trunk/anuga_core/source/anuga/utilities/plot_utils.py

@@ -70 +70 @@
     which would make an object p2 that is like p, but holds information at centroids
     """
-    def __init__(self, filename_list, minimum_allowed_height=1.0e-03):
+    def __init__(self, filename_list, minimum_allowed_height=1.0e-03, verbose=False):
         #
         # Go through the sww files in 'filename_list', and combine them into one object.
@@ -76 +76 @@
 
         for i, filename in enumerate(filename_list):
-            print i, filename
+            if verbose: print i, filename
             # Store output from filename
             p_tmp = get_output(filename, minimum_allowed_height)
@@ -139 +139 @@
        p then contains most relevant information as e.g., p.stage, p.elev, p.xmom, etc 
     """
-    def __init__(self, filename, minimum_allowed_height=1.0e-03):
+    def __init__(self, filename, minimum_allowed_height=1.0e-03, verbose=False):
         self.x, self.y, self.time, self.vols, self.stage, \
                 self.height, self.elev, self.friction, self.xmom, self.ymom, \
                 self.xvel, self.yvel, self.vel, self.minimum_allowed_height,\
                 self.xllcorner, self.yllcorner = \
-                read_output(filename, minimum_allowed_height)
+                read_output(filename, minimum_allowed_height,verbose=verbose)
         self.filename=filename
-
-
-def read_output(filename, minimum_allowed_height):
+        self.verbose = verbose
+
+
+def read_output(filename, minimum_allowed_height,verbose):
     # Input: The name of an .sww file to read data from,
     #                    e.g. read_sww('channel3.sww')
@@ -240 +241 @@
          This is done because presently centroid elevations do not change over time.
     """
-    def __init__(self,p, velocity_extrapolation=False):
+    def __init__(self,p, velocity_extrapolation=False, verbose=False):
         
         self.time, self.x, self.y, self.stage, self.xmom,\
              self.ymom, self.height, self.elev, self.friction, self.xvel, \
              self.yvel, self.vel= \
-             get_centroid_values(p, velocity_extrapolation)
+             get_centroid_values(p, velocity_extrapolation,verbose=verbose)
                                  
 
-def get_centroid_values(p, velocity_extrapolation):
+def get_centroid_values(p, velocity_extrapolation, verbose=False):
     # Input: p is the result of e.g. p=util.get_output('mysww.sww'). See the get_output class defined above
     # Output: Values of x, y, Stage, xmom, ymom, elev, xvel, yvel, vel at centroids
@@ -332 +333 @@
     else:
         # Get centroid values from file 
-        print 'Reading centroids from file'
+        if verbose: print 'Reading centroids from file'
         stage_cent=fid.variables['stage_c'][:]
         elev_cent=fid.variables['elevation_c'][:]

trunk/anuga_core/source/anuga/utilities/sww_merge.py

@@ -356 +356 @@
         ftri_ids = num.where(tri_full_flag>0)
         ftri_l2g = num.compress(tri_full_flag, tri_l2g)
+        
+        #f_ids = num.argwhere(tri_full_flag==1).reshape(-1,)
+        #f_gids = tri_l2g[f_ids]
 
         #print l_volumes
@@ -361 +364 @@
         #print tri_l2g
         #print ftri_l2g
-
-        fg_volumes = num.compress(tri_full_flag,l_volumes,axis=0)
-        g_volumes[ftri_l2g] = fg_volumes
+        
+        f_volumes0 = num.compress(tri_full_flag,volumes[:,0])
+        f_volumes1 = num.compress(tri_full_flag,volumes[:,1])
+        f_volumes2 = num.compress(tri_full_flag,volumes[:,2])
+        
+        g_volumes[ftri_l2g,0] = node_l2g[f_volumes0]
+        g_volumes[ftri_l2g,1] = node_l2g[f_volumes1]
+        g_volumes[ftri_l2g,2] = node_l2g[f_volumes2]
+
+        #fg_volumes = num.compress(tri_full_flag,l_volumes,axis=0)
+        #g_volumes[ftri_l2g] = fg_volumes
 
 
@@ -662 +673 @@
 
         # Just pick out the full triangles
-        ftri_ids = num.where(tri_full_flag>0)
-        ftri_l2g = num.compress(tri_full_flag, tri_l2g)
-
-        assert num.allclose(f_ids, ftri_ids)
-        assert num.allclose(ftri_l2g, f_gids)
+        #ftri_ids = num.where(tri_full_flag>0)
+        #ftri_l2g = num.compress(tri_full_flag, tri_l2g)
+
+        #assert num.allclose(f_ids, ftri_ids)
+        #assert num.allclose(ftri_l2g, f_gids)
 
         g_vids = (3*f_gids.reshape(-1,1) + num.array([0,1,2])).reshape(-1,)
@@ -708 +719 @@
             #             num.array(fid.variables[quantity],dtype=num.float32)
             q = fid.variables[quantity]
-            out_s_c_quantities[quantity][ftri_l2g] = \
-                         num.array(q,dtype=num.float32)[ftri_ids]
+            out_s_c_quantities[quantity][f_gids] = \
+                         num.array(q,dtype=num.float32)[f_ids]
 
         
@@ -717 +728 @@
             #print q.shape
             for i in range(n_steps):
-                out_d_c_quantities[quantity][i][ftri_l2g] = \
-                           num.array(q[i],dtype=num.float32)[ftri_ids]
+                out_d_c_quantities[quantity][i][f_gids] = \
+                           num.array(q[i],dtype=num.float32)[f_ids]
 
         fid.close()

trunk/anuga_core/source/anuga/validation_utilities/parameters.py

@@ -7 +7 @@
 __date__ ="$20/08/2012 11:20:00 PM$"
 
-alg = 'DE1'
+alg = 'DE0'
 cfl = 1.0
 
+# If a file local_parameters.py exits in current directory
+# use the parameters defined there to override the parameters set here.
 try:
-    from anuga_validation_tests.local_parameters import *
+    from local_parameters import *
 except:
     pass

trunk/anuga_core/source/anuga/validation_utilities/run_validation.py

@@ -7 +7 @@
 
 def run_validation_script(script):
-    from anuga.validation_utilities.fabricate import run
+    #from anuga.validation_utilities.fabricate import run
     from anuga.validation_utilities.parameters import alg,cfl
-    run('python', script,  '-alg', alg, '-cfl', cfl)
+    
+    import subprocess
+    cmd = 'python %s -cfl %s -alg %s' % (script, str(cfl), str(alg))
+    print 50*'='
+    print 'Run '+cmd
+    print 50*'='
+    subprocess.call([cmd], shell=True)
+    #run('python', script,  '-alg', alg, '-cfl', cfl)