Changeset 9328

Timestamp:

Sep 9, 2014, 12:38:56 PM (11 years ago)

Author:

davies

Message:

Tweaks to riverwalls + minor changes to template script

Location:

trunk

Files:

• anuga_core/source/anuga/shallow_water/swDE1_domain_ext.c (modified) (3 diffs)
• anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/project.py (modified) (1 diff)
• anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/run_model.py (modified) (2 diffs)
• anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/setup_initial_conditions.py (modified) (1 diff)

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

@@ -678 +678 @@
         if(fabs(edgeflux[0])>1.0e-100){ 
             scaleFlux=newFlux/edgeflux[0];
-
+        }else{
+            scaleFlux = 0.;
+        }
+
+        if(scaleFlux>0.){
+            //if(fabs(scaleFlux) > 10.){
+            //    printf("Very large scaleFlux %e, %e \n", scaleFlux, edgeflux[0]);
+            //}
             // FINAL ADJUSTED FLUX
             edgeflux[0]*=scaleFlux;
-            edgeflux[1]*=scaleFlux;
-            edgeflux[2]*=scaleFlux;
+
+            // FIXME: Do this in a cleaner way
+            // IDEA: Compute momentum flux implied by weir relations, and use
+            //       those in a weighted average (rather than the rescaling trick here)
+            // If we allow the scaling to momentum to be unbounded,
+            // velocity spikes can arise for very-shallow-flooded walls
+            edgeflux[1]*=min(scaleFlux,10.);
+            edgeflux[2]*=min(scaleFlux,10.);
         }else{
             // Can't divide by edgeflux, so enforce 'newFlux' directly
@@ -704 +717 @@
 
     // Adjust the max speed
-    *max_speed_local = sqrt(g*maxhd) + abs(edgeflux[0])/(maxhd+1.0e-100);
+    //*max_speed_local = sqrt(g*(maxhd)) + abs(edgeflux[0])/(maxhd+1.0e-100);
+    *max_speed_local = sqrt(g*(maxhd+weir_height)) + abs(edgeflux[0])/(maxhd+1.0e-100);
 
     return 0;
@@ -856 +870 @@
                     }
 
-                    //////////////////////////////////////////////////////////////////////////////////
-                    // Get Qfactor index - multiply the idealised weir discharge by this constant factor
-                    ii = D->riverwall_rowIndex[RiverWall_count-1] * D->ncol_riverwall_hydraulic_properties;
-                    Qfactor = D->riverwall_hydraulic_properties[ii];
-                    // Get s1, submergence ratio at which we start blending with the shallow water solution 
-                    ii+=1;
-                    s1= D->riverwall_hydraulic_properties[ii];
-                    // Get s2, submergence ratio at which we entirely use the shallow water solution 
-                    ii+=1;
-                    s2= D->riverwall_hydraulic_properties[ii];
-                    // Get h1, tailwater head / weir height at which we start blending with the shallow water solution 
-                    ii+=1;
-                    h1= D->riverwall_hydraulic_properties[ii];
-                    // Get h2, tailwater head / weir height at which we entirely use the shallow water solution 
-                    ii+=1;
-                    h2= D->riverwall_hydraulic_properties[ii];
-                    
-                    // Weir flux adjustment 
-                    adjust_edgeflux_with_weir(edgeflux, h_left_tmp, h_right_tmp, D->g, 
-                                              weir_height, Qfactor, 
-                                              s1, s2, h1, h2, &max_speed_local);
+                    if(h_left_tmp > 0. || h_right_tmp > 0.){
+
+                        //////////////////////////////////////////////////////////////////////////////////
+                        // Get Qfactor index - multiply the idealised weir discharge by this constant factor
+                        ii = D->riverwall_rowIndex[RiverWall_count-1] * D->ncol_riverwall_hydraulic_properties;
+                        Qfactor = D->riverwall_hydraulic_properties[ii];
+                        // Get s1, submergence ratio at which we start blending with the shallow water solution 
+                        ii+=1;
+                        s1= D->riverwall_hydraulic_properties[ii];
+                        // Get s2, submergence ratio at which we entirely use the shallow water solution 
+                        ii+=1;
+                        s2= D->riverwall_hydraulic_properties[ii];
+                        // Get h1, tailwater head / weir height at which we start blending with the shallow water solution 
+                        ii+=1;
+                        h1= D->riverwall_hydraulic_properties[ii];
+                        // Get h2, tailwater head / weir height at which we entirely use the shallow water solution 
+                        ii+=1;
+                        h2= D->riverwall_hydraulic_properties[ii];
+                        
+                        // Weir flux adjustment 
+                        adjust_edgeflux_with_weir(edgeflux, h_left_tmp, h_right_tmp, D->g, 
+                                                  weir_height, Qfactor, 
+                                                  s1, s2, h1, h2, &max_speed_local);
+                    }
                 }
             }

trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/project.py

@@ -34 +34 @@
 ## Total model runtime in seconds (starts at time=0.)
 #
-finalTime=0.01*3600.0
+finaltime=0.01*3600.0
 
 ## Model time between output file timesteps in seconds
 #
-yieldStep=20.0
+yieldstep=20.0
 
 ## Spatial Information

trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/run_model.py

@@ -55 +55 @@
 ##################################################################################
 print 'Making initial conditions'
-domain=setup_initial_conditions.setup_initial_conditions(domain,project)
+setup_initial_conditions.setup_initial_conditions(domain,project)
 
 ####################################################################################
@@ -96 +96 @@
 
 barrier()
-for t in domain.evolve(yieldstep=project.yieldStep, finaltime=project.finalTime):
+for t in domain.evolve(yieldstep=project.yieldstep, finaltime=project.finaltime):
     if(myid==0): print domain.timestepping_statistics()
    

trunk/anuga_work/development/gareth/template_scenarios/tsunami/tsunami1/scripts/setup_initial_conditions.py

@@ -230 +230 @@
         domain.riverwallData.export_riverwalls_to_text(output_dir=project.output_dir+'/'+project.spatial_txt_outputDir)
 
-    return(domain)
-
-
+    return
+
+