Changeset 8866

Timestamp:

May 13, 2013, 5:09:25 PM (12 years ago)

Author:

davies

Message:

Updates to bal_dev

Location:

trunk/anuga_work/development/gareth

Files:

• experimental/balanced_dev/swb2_domain.py (modified) (8 diffs)
• experimental/balanced_dev/swb2_domain_ext.c (modified) (20 diffs)
• tests/runup_sinusoid/runup_sinusoid.py (modified) (1 diff)

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

@@ -66 +66 @@
         self.set_CFL(1.0)
         self.set_use_kinematic_viscosity(False)
-        self.timestepping_method='rk2'#'rk2'#'euler'#'rk2' 
+        self.timestepping_method='rk2'#'rk2'#'rk2'#'euler'#'rk2' 
         # The following allows storage of the negative depths associated with this method
         self.minimum_storable_height=-99999999999.0 
@@ -74 +74 @@
         self.extrapolate_velocity_second_order=True
 
-        # Note that the extrapolation method used in quantity_ext.c (e.g.
-        # extrapolate_second_order_and_limit_by_edge) uses a constant value for
-        # all the betas.  
         self.beta_w=0.0
         self.beta_w_dry=0.0
@@ -83 +80 @@
         self.beta_vh=0.0
         self.beta_vh_dry=0.0
+
+        #self.epsilon=1.0e-100
 
         #self.optimise_dry_cells=True
@@ -111 +110 @@
         self.boundary_flux_sum=numpy.ndarray(1)
         self.boundary_flux_sum[0]=0.
+
+        self.call=1 # Integer counting how many times we call compute_fluxes_central
+
+        # Integer recording the order of the time-stepping scheme
+        if(self.timestepping_method=='rk2'):
+          self.timestep_order=2
+        elif(self.timestepping_method=='euler'):
+          self.timestep_order=1
+        elif(self.timestepping_method=='rk3'):
+          self.timestep_order=3
+        else:
+          err_mess='ERROR: timestepping_method= ' + self.timestepping_method +' not supported in this solver'
+          raise Exception, err_mess
 
         print '##########################################################################'
@@ -252 +264 @@
         #    raise Exception, err_mess 
 
-        if(domain.timestepping_method=='rk2'):
-          timestep_order=2
-        elif(domain.timestepping_method=='euler'):
-          timestep_order=1
+        #if(domain.timestepping_method=='rk2'):
+        #  timestep_order=2
+        #elif(domain.timestepping_method=='euler'):
+        #  timestep_order=1
         #elif(domain.timestepping_method=='rk3'):
         #  timestep_order=3
-        else:
-          err_mess='ERROR: domain.timestepping_method= ' + domain.timestepping_method +' not supported in this solver'
-          raise Exception, err_mess
-
-        #print 'timestep_order=', timestep_order
+        #else:
+        #  err_mess='ERROR: domain.timestepping_method= ' + domain.timestepping_method +' not supported in this solver'
+        #  raise Exception, err_mess
+
+        ##print 'timestep_order=', timestep_order
 
         Stage = domain.quantities['stage']
@@ -271 +283 @@
         timestep = float(sys.maxint)
 
+        domain.call+=1
         flux_timestep = compute_fluxes_ext(timestep,
                                            domain.epsilon,
@@ -297 +310 @@
                                            domain.max_speed,
                                            int(domain.optimise_dry_cells),
-                                           timestep_order,
+                                           domain.timestep_order,
                                            Stage.centroid_values,
                                            Xmom.centroid_values,
@@ -304 +317 @@
                                            Bed.vertex_values)
 
-        #import pdb
-        #pdb.set_trace()
-        #print 'flux timestep: ', flux_timestep, domain.timestep
+
+        # Update the boundary flux integral
         if(domain.timestepping_method=='rk2'):
-          if(flux_timestep == float(sys.maxint)):
+            if(domain.call%2==1):
               domain.boundary_flux_integral[0]= domain.boundary_flux_integral[0] +\
                                                 domain.boundary_flux_sum[0]*domain.timestep*0.5
               #print 'dbfi ', domain.boundary_flux_integral, domain.boundary_flux_sum
               domain.boundary_flux_sum[0]=0.
+
+        elif(domain.timestepping_method=='euler'):
+            domain.boundary_flux_integral=0.
+            # This presently doesn't work -- this section of code may need to go elsewhere
+            #domain.boundary_flux_integral[0]= domain.boundary_flux_integral[0] +\
+            #                                  domain.boundary_flux_sum[0]*domain.timestep
+            #domain.boundary_flux_sum[0]=0.
+
+        elif(domain.timestepping_method=='rk3'):
+            domain.boundary_flux_integral=0.
+            # This needs to be designed
+        else:
+            mess='ERROR: domain.timestepping_method', domain.timestepping_method,' method not recognised'
+            raise Exception, mess
 
         domain.flux_timestep = flux_timestep

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

@@ -85 +85 @@
 }
 
+// minmod limiter
+int _minmod(double a, double b){
+    // Compute minmod
+
+    if(sign(a)!=sign(b)){
+        return 0.0;
+    }else{
+        return min(fabs(a), fabs(b))*sign(a);
+    }
+
+
+}
+
 // Innermost flux function (using stage w=z+h)
 int _flux_function_central(double *q_left, double *q_right,
@@ -116 +129 @@
   double s_min, s_max, soundspeed_left, soundspeed_right;
   double denom, inverse_denominator, z;
-  double uint, t1, t2, t3;
+  double uint, t1, t2, t3, min_speed;
   // Workspace (allocate once, use many)
   static double q_left_rotated[3], q_right_rotated[3], flux_right[3], flux_left[3];
@@ -145 +158 @@
   if(h_left>h0){
     u_left = uh_left/max(h_left, 1.0e-06);
+    uh_left=h_left*u_left;
   }else{
     u_left=0.;
+    uh_left=h_left*u_left;
   }
   //u_left = _compute_speed(&uh_left, &h_left, 
@@ -158 +173 @@
   if(h_right>h0){
     u_right = uh_right/max(h_right, 1.0e-06);
+    uh_right=h_right*u_right;
   }else{
     u_right=0.;
+    uh_right=h_right*u_right;
   }
 
@@ -236 +253 @@
   // Flux computation
   denom = s_max - s_min;
-  if (denom < epsilon) 
+  //if (denom < -epsilon) 
+  if (0==1) 
   { 
     // Both wave speeds are very small
@@ -242 +260 @@
 
     *max_speed = 0.0;
-    *pressure_flux = 0.0;
-    //*pressure_flux = 0.5*g*0.5*(h_left*h_left+h_right*h_right);
+    //*pressure_flux = 0.0;
+    *pressure_flux = 0.5*g*0.5*(h_left*h_left+h_right*h_right);
   } 
   else 
@@ -249 +267 @@
     // Maximal wavespeed
     *max_speed = max(s_max, -s_min);
-    
-    inverse_denominator = 1.0/denom;
+    //min_speed=min(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] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]);
+
+      // 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);
+
       //GD HACK: Here, we supress the 'smoothing' part of the flux -- like scaling diffusion by local wave speed
+      //if(i!=2){
       edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i])*
-                      (min(*max_speed/(max(speed_max_last,1e-06)),1.0));//(min(hc/h_left,min(hc_n/h_right,1.0)));
+                     (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] *= inverse_denominator;
     }
@@ -449 +481 @@
             // bedslope_work contains all gravity related terms -- weighting of
             // conservative and non-conservative versions.
-            if(hc> -h0){
+            if(hc> h0){
               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);
@@ -471 +503 @@
                 edgeflux_store[nm3 + 2 ] = edgeflux[2];
 
-                if(hc_n> -h0){
+                if(hc_n> h0){
                   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);
@@ -493 +525 @@
             if ((tri_full_flag[k] == 1) ) {
 
-                if(call%timestep_order!=0) max_speed = max_speed_array[k]; // HACK to Ensure that local timestep is the same as the last timestep 
+                if((call-2)%timestep_order!=0) max_speed = max_speed_array[k]; // HACK to Ensure that local timestep is the same as the last timestep 
 
                 if (max_speed > epsilon) {
@@ -585 +617 @@
             // GD HACK
             // Option to limit advective fluxes
-            if(hc > h0){
+            if(hc > -h0){
                 stage_explicit_update[k] += edgeflux_store[ki3+0];
                 // Cut these terms for shallow flows, and protect stationary states from round-off error
@@ -598 +630 @@
             if(n<0){ 
                 boundary_flux_sum[0] += edgeflux_store[ki3];
-                //printf("%f, %f, %d \n", boundary_flux_sum[0], timestep, ki3);
             }
 
             // GD HACK
             // Compute bed slope term
-            if(hc > h0){
+            if(hc > -h0){
                 xmom_explicit_update[k] -= normals[ki2]*pressuregrad_store[ki];
                 ymom_explicit_update[k] -= normals[ki2+1]*pressuregrad_store[ki];
@@ -621 +652 @@
     }  // end cell k
 
-    if(call%timestep_order==0) timestep=local_timestep;
+    if((call-2)%timestep_order==0) timestep=local_timestep; // Hack to ensure we only update the timestep on the first call within each rk2/rk3 step
 
     // Look for 'dry' edges, and set momentum (& component of momentum_update)
@@ -729 +760 @@
             //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]*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;
 
         if (hc <= 0.0){
@@ -805 +838 @@
   // Any provisional jump with magnitude < TINY does not contribute to 
   // the limiting process.
+  //return 0;
   
   for (i=0;i<3;i++){
@@ -821 +855 @@
   dqv[1]=dqv[1]*phi;
   dqv[2]=dqv[2]*phi;
+
+  
+  //printf("%lf \n ", beta_w);
     
   return 0;
@@ -1120 +1157 @@
               }else{
                 // Constant stage extrapolation
-                stage_edge_values[k3]   = stage_centroid_values[k];
-                stage_edge_values[k3+1] = stage_centroid_values[k];
-                stage_edge_values[k3+2] = stage_centroid_values[k];
+                //stage_edge_values[k3]   = stage_centroid_values[k];
+                //stage_edge_values[k3+1] = stage_centroid_values[k];
+                //stage_edge_values[k3+2] = stage_centroid_values[k];
               }
             //}else{
@@ -1223 +1260 @@
               //}
               // First order momentum / velocity extrapolation
-              //stage_edge_values[k3]    = stage_centroid_values[k];
-              //stage_edge_values[k3+1]  = stage_centroid_values[k];
-              //stage_edge_values[k3+2]  = stage_centroid_values[k];
+              stage_edge_values[k3]    = stage_centroid_values[k];
+              stage_edge_values[k3+1]  = stage_centroid_values[k];
+              stage_edge_values[k3+2]  = stage_centroid_values[k];
               xmom_edge_values[k3]    = xmom_centroid_values[k];
               xmom_edge_values[k3+1]  = xmom_centroid_values[k];
@@ -1251 +1288 @@
             //hfactor=hmin/(hmin + 0.004);
           //}
-          hfactor= min(2.0*max(hmin,0.0)/max(hc,max(0.5*(bedmax-bedmin),minimum_allowed_height*10.)), 1.0);
+          hfactor= min(2.0*max(hmin,0.0)/max(hc,max(0.5*(bedmax-bedmin),minimum_allowed_height*2.)), 1.0);
           //hfactor= min(2.0*max(hmin,0.0)/max(hc,max(0.5*(bedmax-bedmin),minimum_allowed_height)), 1.0);
           
@@ -1317 +1354 @@
           //if(count_wet_neighbours[k]>0){
           limit_gradient(dqv, qmin, qmax, beta_tmp);
+          //printf("%lf \n ", beta_tmp);
           //}
           //}

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

@@ -94 +94 @@
     print 'Peak velocity is: ', vv.max(), vv.argmax(), dd[vv.argmax()]
     print 'Volume is', sum(dd_raw*domain.areas)    
-    #print 'Volume less flux int', sum(dd_raw*domain.areas) - domain.boundary_flux_integral
+    print 'Volume less flux int', sum(dd_raw*domain.areas) - domain.boundary_flux_integral