Changeset 9287

Timestamp:

Aug 6, 2014, 9:14:12 AM (11 years ago)

Author:

steve

Message:

Fixed bug in volume reporting with negative rate operators

Location:

trunk/anuga_core

Files:

• source/anuga/operators/rate_operators.py (modified) (2 diffs)
• source/anuga/operators/test_rate_operators.py (modified) (3 diffs)
• validation_tests/case_studies/towradgi/simulation.py (modified) (1 diff)
• validation_tests/experimental_data/okushiri (moved) (moved from trunk/anuga_core/validation_tests/case_studies/okushiri)
• validation_tests/reports/all_tests_report.tex (modified) (3 diffs)

trunk/anuga_core/source/anuga/operators/rate_operators.py

@@ -117 +117 @@
         if self.verbose is True:
             log.critical('Rate of %s at time = %.2f = %f'
-                         % (self.quantity_name, domain.get_time(), rate))
+                         % (self.quantity_name, self.domain.get_time(), rate))
 
         fid = self.full_indices
@@ -123 +123 @@
             # Record the local flux for mass conservation tracking
             if indices is None:
-                self.local_influx=factor*timestep*(rate*self.areas)[fid].sum()
-                self.stage_c[:] = self.stage_c[:]  \
-                       + factor*rate*timestep
-            else:
-                self.local_influx=factor*timestep*(rate*self.areas)[fid].sum()
+                local_rates=factor*timestep*rate
+                self.local_influx = (local_rates*self.areas)[fid].sum()
+                self.stage_c[:] = self.stage_c[:] + local_rates
+            else:
+                local_rates = factor*timestep*rate
+                self.local_influx=(local_rates*self.areas[fid]).sum()
                 self.stage_c[indices] = self.stage_c[indices] \
-                       + factor*rate*timestep
+                       + local_rates
         else: # Be more careful if rate < 0
             if indices is None:
-                self.local_influx=(num.minimum(factor*timestep*rate, self.stage_c[:]-self.elev_c[:])*self.areas)[fid].sum()
-                self.stage_c[:] = num.maximum(self.stage_c  \
-                       + factor*rate*timestep, self.elev_c )
-            else:
-                self.local_influx=(num.minimum(factor*timestep*rate, self.stage_c[indices]-self.elev_c[indices])*self.areas)[fid].sum()
-                self.stage_c[indices] = num.maximum(self.stage_c[indices] \
-                       + factor*rate*timestep, self.elev_c[indices])
+                #self.local_influx=(num.minimum(factor*timestep*rate, self.stage_c[:]-self.elev_c[:])*self.areas)[fid].sum()
+                #self.stage_c[:] = num.maximum(self.stage_c  \
+                #       + factor*rate*timestep, self.elev_c )
+                local_rates = num.maximum(factor*timestep*rate, self.elev_c[:]-self.stage_c[:])
+                self.local_influx = (local_rates*self.areas)[fid].sum()
+                self.stage_c[:] = self.stage_c + local_rates 
+            else:
+                #self.local_influx=(num.minimum(factor*timestep*rate, self.stage_c[indices]-self.elev_c[indices])*self.areas)[fid].sum()
+                #self.stage_c[indices] = num.maximum(self.stage_c[indices] \
+                #       + factor*rate*timestep, self.elev_c[indices])
+                
+                local_rates = num.maximum(factor*timestep*rate, self.elev_c[indices]-self.stage_c[indices])
+                self.local_influx=(local_rates*self.areas)[fid].sum()
+                self.stage_c[indices] = self.stage_c[indices] + local_rates 
         # Update mass inflows from fractional steps
         self.domain.fractional_step_volume_integral+=self.local_influx

trunk/anuga_core/source/anuga/operators/test_rate_operators.py

@@ -25 +25 @@
 
 
+warnings.simplefilter("ignore")
+
 
 class Test_rate_operators(unittest.TestCase):
@@ -140 +142 @@
 
         stage_ex = [ 0.,  0.,   1.,  0.]
+        step_integral = -6.0
+
+        #print domain.quantities['elevation'].centroid_values
+        #print domain.quantities['stage'].centroid_values
+        #print domain.quantities['xmomentum'].centroid_values
+        #print domain.quantities['ymomentum'].centroid_values
+        #print domain.fractional_step_volume_integral
+        #print factor*domain.timestep*(rate*domain.areas[indices]).sum()
+        
+        #increment = factor*domain.timestep*rate*domain.areas
+        
+        
+
+        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
+        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
+        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
+        assert num.allclose(domain.fractional_step_volume_integral, step_integral)
+
+    def test_rate_operator_negative_rate_full(self):
+        from anuga.config import rho_a, rho_w, eta_w
+        from math import pi, cos, sin
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #             bac,     bce,     ecf,     dbe
+        vertices = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        domain = Domain(points, vertices)
+
+        #Flat surface with 1m of water
+        domain.set_quantity('elevation', 0)
+        domain.set_quantity('stage', 10.0)
+        domain.set_quantity('friction', 0)
+
+        Br = Reflective_boundary(domain)
+        domain.set_boundary({'exterior': Br})
 
 #        print domain.quantities['elevation'].centroid_values
@@ -146 +190 @@
 #        print domain.quantities['ymomentum'].centroid_values
 
-        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
-        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
-        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
-        assert num.allclose(domain.fractional_step_volume_integral, factor*domain.timestep*(rate*domain.areas[indices]).sum())
-
-
+        # Apply operator to these triangles
+        indices = [0,1,3]
+
+
+
+        #Catchment_Rain_Polygon = read_polygon(join('CatchmentBdy.csv'))
+        #rainfall = file_function(join('1y120m.tms'), quantities=['rainfall'])
+        rate = -1.0
+        factor = 10.0
+        default_rate= 0.0
+
+        operator = Rate_operator(domain, rate=rate, factor=factor, \
+                      indices=None, default_rate = default_rate)
+
+
+        # Apply Operator
+        domain.timestep = 2.0
+        operator()
+
+        stage_ex = [ 0.,  0.,   0.,  0.]
+        step_integral = -80.0
+
+        #print domain.quantities['elevation'].centroid_values
+        #print domain.quantities['stage'].centroid_values
+        #print domain.quantities['xmomentum'].centroid_values
+        #print domain.quantities['ymomentum'].centroid_values
+        #print domain.fractional_step_volume_integral
+
+
+        assert num.allclose(domain.quantities['stage'].centroid_values, stage_ex)
+        assert num.allclose(domain.quantities['xmomentum'].centroid_values, 0.0)
+        assert num.allclose(domain.quantities['ymomentum'].centroid_values, 0.0)
+        assert num.allclose(domain.fractional_step_volume_integral, step_integral)
+        
     def test_rate_operator_rate_from_file(self):
         from anuga.config import rho_a, rho_w, eta_w

trunk/anuga_core/validation_tests/case_studies/towradgi/simulation.py

@@ -27 +27 @@
 
 #from project import *
-
 
 

trunk/anuga_core/validation_tests/reports/all_tests_report.tex

@@ -122 +122 @@
 %%======================
 \inputresults{../experimental_data/dam_break_yeh_petroff}
+\inputresults{../experimental_data/okushiri}
+
 \inputresults{../behaviour_only/lid_driven_cavity}
 \inputresults{../behaviour_only/weir_1}
+\inputresults{../behaviour_only/bridge_hecras}
+\inputresults{../behaviour_only/lateral_weir_hecras}
+\inputresults{../behaviour_only/tides_hecras}
 
 \inputresults{../other_references/radial_dam_break_dry}
 \inputresults{../other_references/radial_dam_break_wet}
-
-\inputresults{../case_studies/okushiri}
 
 
@@ -139 +142 @@
 
 
-To setup a new validation test, create a test directory under the
-\textsc{Tests} directory. In that directory there should be the test code, a
+To setup a new validation test, create a test directory under one of 
+the validation directories. In that directory there should be the test code, a
 \TeX{} file \texttt{results.tex} and a python script
 \texttt{produce\_results.py}, which runs the simulation and produces the
-outputs. Copy the format from the other test directory. 
+outputs. Copy the format from one of the other test directories. 
 
 In this \TeX{} file, \texttt{report.tex}, add a line
 \begin{verbatim}
-\inputresults{Tests/Directory/Name}
+\inputresults{../Directory/Name}
 \end{verbatim}
 
@@ -175 +178 @@
 
 \begin{verbatim}
-from anuga.validation_utilities.fabricate import *
-from anuga.validation_utilities import run_validation_script
+import anuga
+from anuga.validation_utilities import produce_report
 
-# Setup the python scripts which produce the output for this
-# validation test
-def build():
-    run_validation_script('run_problem.py')
-    run_validation_script('plot_problem.py')
-    run('python','typeset_report.py')
-    pass
+args = anuga.get_args()
 
-def clean():
-    autoclean()
-
-main()
+produce_report('run_simulation.py', args=args)
 \end{verbatim}
-This script uses \texttt{fabricate} which automatically determines dependences
-and only runs the command if the parameters alg and cfl have changed,
-or input/output or source files have changed. (\texttt{fabricate} is a replacement for
-standard \texttt{make})