Changeset 8866 for trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c

Timestamp:

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

Author:

davies

Message:

Updates to bal_dev

File:

• trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c (modified) (20 diffs)

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);
           //}
           //}