Changeset 9232

Timestamp:

Jul 2, 2014, 10:16:31 AM (11 years ago)

Author:

davies

Message:

Adjustment to riverwall implementation which seems to give better stability

Location:

trunk

Files:

• anuga_core/source/anuga/shallow_water/swDE1_domain_ext.c (modified) (5 diffs)
• anuga_core/validation_tests/behaviour_only/weir_1/plot_results.py (modified) (1 diff)
• anuga_core/validation_tests/behaviour_only/weir_1/results.tex (modified) (1 diff)
• anuga_core/validation_tests/behaviour_only/weir_1/runup.py (modified) (2 diffs)
• anuga_work/development/gareth/tests/levee2/runup.py (modified) (2 diffs)
• anuga_work/development/gareth/tests/levee2/runuplot.py (modified) (1 diff)

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

@@ -505 +505 @@
         }else{
             // Can't divide by edgeflux
+            // FIXME: This is fluxing mass but not momentum
+            //        It might be ok, but,.........
+            //        it has potential be a problem in the
+            //        case that there was zero edgeflux but non-zero weir flux
+            //        [which can occur because the edgeflux is computed from
+            //        edge quantities, whereas weir-flux is computed from
+            //        centroid info]. Then if we remove mass but not momentum
+            //        from the headwater, velocity gets faster.
+            //        
+            //        However, that situation would usually occur for a short time with a small flux,
+            //        so it might never actually create a problem.        
             edgeflux[0] = newFlux;
             edgeflux[1]*=0.;
@@ -576 +587 @@
     double hle, hre, zc, zc_n, Qfactor, s1, s2, h1, h2; 
     double stage_edge_lim, outgoing_mass_edges, bedtop, bedbot, pressure_flux, hc, hc_n, tmp, tmp2;
+    double h_left_tmp, h_right_tmp;
     static long call = 1; // Static local variable flagging already computed flux
     double speed_max_last, vol, smooth, local_speed, weir_height;
@@ -657 +669 @@
             z_half=max(zl,zr);
 
-            // Account for riverwalls
+            //// Account for riverwalls
             if(edge_flux_type[ki]==1){
                 // Update counter of riverwall edges == index of
@@ -664 +676 @@
                 
                 // Set central bed to riverwall elevation
-                z_half=max(riverwall_elevation[RiverWall_count-1], max(zl, zr)) ;
-
-                if(min(ql[0], qr[0]) < z_half){
-                    // Since there is a wall blocking the flow connection, use first order extrapolation for this edge
-                    ql[0]=stage_centroid_values[k];
-                    ql[1]=xmom_centroid_values[k];
-                    ql[2]=ymom_centroid_values[k];
-                    hle=hc;
-                    zl=zc;
-
-                    if(n>=0){
-                      qr[0]=stage_centroid_values[n];
-                      qr[1]=xmom_centroid_values[n];
-                      qr[2]=ymom_centroid_values[n];
-                      hre=hc_n;
-                      zr = zc_n;
-                    }else{
-                      hre=hc;
-                      zr = zc;
-                    }
-                    // Re-set central bed to riverwall elevation
-                    z_half=max(riverwall_elevation[RiverWall_count-1], max(zl, zr)) ;
-                }
-                
+                z_half=max(riverwall_elevation[RiverWall_count-1], z_half) ;
+
+                //if(min(ql[0], qr[0]) < z_half){
+                //    // Since there is a wall blocking the flow connection, use first order extrapolation for this edge
+                //    ql[0]=stage_centroid_values[k];
+                //    ql[1]=xmom_centroid_values[k];
+                //    ql[2]=ymom_centroid_values[k];
+                //    hle=hc;
+                //    zl=zc;
+
+                //    if(n>=0){
+                //      qr[0]=stage_centroid_values[n];
+                //      qr[1]=xmom_centroid_values[n];
+                //      qr[2]=ymom_centroid_values[n];
+                //      hre=hc_n;
+                //      zr = zc_n;
+                //    }else{
+                //      hre=hc;
+                //      zr = zc;
+                //    }
+                //    // Re-set central bed to riverwall elevation
+                //    z_half=max(riverwall_elevation[RiverWall_count-1], max(zl, zr)) ;
+                //}
 
             }
@@ -728 +739 @@
                     ii+=1;
                     h2=riverwall_hydraulic_properties[ii];
-
-                    //printf("%e, %e, %e, %e, %e \n", Qfactor, s1, s2, h1, h2);
-
-                    adjust_edgeflux_with_weir(edgeflux, h_left, h_right, g, 
+                    
+                    //adjust_edgeflux_with_weir(edgeflux, h_left, h_right, g, 
+                    //                          weir_height, Qfactor, 
+                    //                          s1, s2, h1, h2);
+
+                    // Use first-order h's for weir -- as the 'upstream/downstream' heads are
+                    //  measured away from the weir itself
+                    h_left_tmp= max(stage_centroid_values[k]-z_half,0.);
+                    if(n>=0){
+                        h_right_tmp= max(stage_centroid_values[n]-z_half,0.);
+                    }else{
+                        h_right_tmp= max(hc_n+zr-z_half,0.);
+                    }
+
+                    adjust_edgeflux_with_weir(edgeflux, h_left_tmp, h_right_tmp, g, 
                                               weir_height, Qfactor, 
                                               s1, s2, h1, h2);

trunk/anuga_core/validation_tests/behaviour_only/weir_1/plot_results.py

@@ -76 +76 @@
 #C=C*(ft_to_m)**0.5 # Convert to m^0.5/s
 C=2./3.*(9.81*2./3.)**0.5 # Another standard coefficient.
-L1=9.0 # In the model, the L1 region is 8m with two 1m rises at each edge, leading to the equivalent of 9m width on average
+L1=9.0 # Lengths of the w1/w2 regions
 L2=91.
 def simple_weir(H, C, L):

trunk/anuga_core/validation_tests/behaviour_only/weir_1/results.tex

@@ -4 +4 @@
 In this problem water runs up a broad sloping floodplain, which is split into 2 by a thin weir (riverwall). The landward side of the floodplain is initially dry but becomes wet as the water overtops the weir. We use the change in mass on the landward side of the weir to compute the flux over the weir according to \anuga{}, and compare with direct computation of the weir equation. 
 
-\anuga{} allows computation of flow over a weir using the riverwall structure (supported for DE1 only as of 19/05/2014). The default method adjusts the edge flux over the weir to satisfy a basic weir relation with Villemonte's submergence correction. However at high submergence ratios, or when the depth of flow over the weir is large compared with the weir height, \anuga{} smoothly reverts to the shallow water solution (because the weir relations are not sensible in these situations). This is required so that e.g. a weir of 1cm height covered by flow of 1m is basically no different from the shallow water solution - weir relations by themselves will not achieve this. See the documentation of riverwalls for more information.
+\anuga{} allows computation of flow over a weir using the riverwall structure (supported for Discontinuous Elevation Algorithms only as of 01/07/2014). The default method adjusts the edge flux over the weir to satisfy a basic weir relation with Villemonte's submergence correction. However at high submergence ratios, or when the depth of flow over the weir is large compared with the weir height, \anuga{} smoothly reverts to the shallow water solution (because the weir relations are not sensible in these situations). This is required so that e.g. a weir of 1cm height covered by flow of 1m is basically no different from the shallow water solution - weir relations by themselves will not achieve this. See the documentation of riverwalls for more information.
 
 \subsection{Results}

trunk/anuga_core/validation_tests/behaviour_only/weir_1/runup.py

@@ -22 +22 @@
 higherResPolygon=[ [40., 40.], [40., 60.], [60., 60.], [60., 40.]]
 # Riverwall = list of lists, each with a set of x,y,z (and optional QFactor) values
-riverWall={ 'centralWall':
+#riverWall={ 'centralWall':
+#                [ [50., 0.0, -0.0],
+#                  [50., 45., -0.0],
+#                  [50., 46., -0.2],
+#                  [50., 54., -0.2],
+#                  [50., 55., -0.0], 
+#                  [50., 100.0, -0.0]] 
+#          }
+riverWall={ 'leftWall':
                 [ [50., 0.0, -0.0],
-                  [50., 45., -0.0],
-                  [50., 46., -0.2],
-                  [50., 54., -0.2],
-                  [50., 55., -0.0], 
+                  [50., 45.5, -0.0]],
+             'centralWall': 
+                [[50., 45.5, -0.2],
+                 [50., 54.5, -0.2]],
+             'rightWall':
+                [ [50., 54.5, -0.0], 
                   [50., 100.0, -0.0]] 
           }
 
-riverWall_Par={'centralWall':{'Qfactor':1.0}}
+#riverWall_Par={'centralWall':{'Qfactor':1.0}}
 # Try to avoid any shallow-water type solution -- becomes unstable
 #riverWall_Par={'centralWall':{'Qfactor':1.0, 's1': 0.999, 's2':0.9999, 'h1':100, 'h2':150}}
@@ -98 +108 @@
 domain = distribute(domain)
 
-domain.riverwallData.create_riverwalls(riverWall, riverWall_Par)
+domain.riverwallData.create_riverwalls(riverWall)
 
 #--------------------------

trunk/anuga_work/development/gareth/tests/levee2/runup.py

@@ -17 +17 @@
 higherResPolygon=[ [40., 40.], [40., 60.], [60., 60.], [60., 40.]]
 # Riverwall = list of lists, each with a set of x,y,z (and optional QFactor) values
-riverWall={ 'centralWall':
+riverWall={ 'leftWall':
                 [ [50., 0.0, -0.0],
-                  [50., 45., -0.0],
-                  [50., 46., -0.2],
-                  [50., 54., -0.2],
-                  [50., 55., -0.0], 
+                  [50., 45.5, -0.0]],
+             'centralWall': 
+                [[50., 45.5, -0.2],
+                 [50., 54.5, -0.2]],
+             'rightWall':
+                [ [50., 54.5, -0.0], 
                   [50., 100.0, -0.0]] 
           }
@@ -59 +61 @@
 domain.set_datadir('.')                         
 domain.set_flow_algorithm('DE1')
-domain.set_store_vertices_uniquely()
+domain.set_store_vertices_uniquely(True)
 domain.riverwallData.create_riverwalls(riverWall, riverWall_Par)
 

trunk/anuga_work/development/gareth/tests/levee2/runuplot.py

@@ -46 +46 @@
 
 ## Get reference stage points
-seaLev=(abs(p.x-51.)+abs(p.y-50.)).argmin() # Index near levee on sea side
-landLev=(abs(p.x-49.)+abs(p.y-50.)).argmin() # Index near levee on land side
+seaLev=(abs(p.x-50.7)+abs(p.y-50.)).argmin() # Index near levee on sea side
+landLev=(abs(p.x-49.3)+abs(p.y-50.)).argmin() # Index near levee on land side
 heightDiff=p.stage[:,seaLev]-p.stage[:,landLev]
 ## Get indices on 'landward' side of riverwall