Changeset 9681 for trunk/anuga_core/anuga/structures

Timestamp:

Feb 23, 2015, 9:21:46 AM (10 years ago)

Author:

davies

Message:

Major changes + bugfix for internal_boundary_operator and related

Location:

trunk/anuga_core/anuga/structures

Files:

• boyd_box_operator_Amended3.py (modified) (1 diff)
• internal_boundary_operator.py (modified) (4 diffs)
• structure_operator.py (modified) (11 diffs)

trunk/anuga_core/anuga/structures/boyd_box_operator_Amended3.py

@@ -32 +32 @@
                  verbose=False):
 
-                     
+        raise Exception('(Feb 2015) This code seems broken -- the structure operator call is incorrect -- retire?')
+
         anuga.Structure_operator.__init__(self,
                                           domain,

trunk/anuga_core/anuga/structures/internal_boundary_operator.py

@@ -3 +3 @@
 import numpy
 from numpy.linalg import solve
+import scipy
+import scipy.optimize as sco
 
 
@@ -75 +77 @@
                                           zero_outflow_momentum=zero_outflow_momentum,
                                           use_old_momentum_method=False,
+                                          always_use_Q_wetdry_adjustment=False,
                                           force_constant_inlet_elevations=force_constant_inlet_elevations,
                                           description=description,
@@ -227 +230 @@
     def discharge_routine_implicit(self):
         """Uses semi-implicit discharge estimation:
-          Discharge = 0.5*(Q(H0, T0) + Q(H0 + delta_H, T0+delta_T))
+          Discharge = (1-theta)*Q(H0, T0) + theta*Q(H0 + delta_H, T0+delta_T))
         where H0 = headwater stage, T0 = tailwater stage, delta_H = change in
         headwater stage over a timestep, delta_T = change in tailwater stage over a
-        timestep, and Q is the discharge function
-
-        We can estimate delta_H, delta_T by solving the following system (based on mass conservation):
-          A0*delta_H = -dt/2*( Q(H0, T0) + Q(H0 + delta_H, T0 + delta_T))
-          A1*delta_T =  dt/2*( Q(H0, T0) + Q(H0 + delta_H, T0 + delta_T))
-        where A0, A1 are the inlet areas, and dt is the timestep
-
-        We linearise the system with the approximation:
-          Q(H0 + delta_H, T0 + delta_T) ~= Q(H0, T0) + delQ/delH * delta_H + delQ/delT*delta_T
-        where delQ/delH and delQ/delT are evaluated with numerical finite differences 
-
+        timestep, and Q is the discharge function, and theta is a constant in
+        [0,1] determining the degree of implicitness (currently hard-coded).
+
+        Note this is effectively assuming:
+        1) Q is a function of stage, not energy (so we can relate mass change directly to delta_H, delta_T). We
+           could generalise it to the energy case ok.
+        2) The stage is computed on the exchange line (or the change in
+            stage at the enquiry point is effectively the same as that on the exchange
+            line)
 
         """
@@ -263 +264 @@
             E0 = self.inlet0_energy
             E1 = self.inlet1_energy
-
-            # Numerical derivatives
-            dE0 = 0.01
-            dE1 = 0.01
-            dQ_dE0 = (self.internal_boundary_function(E0+dE0, E1) - self.internal_boundary_function(E0-dE0, E1))/(2*dE0)
-            dQ_dE1 = (self.internal_boundary_function(E0, E1+dE1) - self.internal_boundary_function(E0, E1-dE1))/(2*dE1)
-
-            # Assemble and solve matrix
             A0 = self.inlets[0].area
             A1 = self.inlets[1].area
-            hdt = 0.5*dt
-
-            M11 = A0 + dQ_dE0*hdt
-            M12 = dQ_dE1*hdt
-            M21 = -dQ_dE0*hdt
-            M22 = A1 - dQ_dE1*hdt
-
-            lhs = numpy.array([ [M11, M12], [M21, M22]])
-            rhs = numpy.array([ -Q0*dt, Q0*dt])
-            # sol contains delta_E0, delta_E1
-            sol = solve(lhs, rhs)
-           
-            #Q = 0.5*(Q0 + ( Q0 + sol[0]*dQ_dE0 + sol[1]*dQ_dE1))
-            Q1 =  self.internal_boundary_function(E0 + sol[0], E1 + sol[1])
-            Q = 0.5*(Q0 + Q1)
+            theta = 1.0
+
+            sol = numpy.array([0., 0.]) # estimate of (delta_H, delta_T)
+            areas = numpy.array([A0, A1])
+
+            # Use scipy root finding
+            def F_to_solve(sol):
+                Q1 =  self.internal_boundary_function(E0 + sol[0], E1 + sol[1])
+                discharge = (1.0-theta)*Q0 + theta*Q1
+                output = sol*areas - discharge*dt*numpy.array([-1., 1.])
+                return(output) 
+
+            final_sol = sco.root(F_to_solve, sol, method='lm').x
+            Q1 =  self.internal_boundary_function(E0 + final_sol[0], E1 + final_sol[1])
+            Q = (1.0-theta)*Q0 + theta*Q1
+
         else:
             Q = Q0

trunk/anuga_core/anuga/structures/structure_operator.py

@@ -39 +39 @@
                  zero_outflow_momentum=True,
                  use_old_momentum_method=True,
+                 always_use_Q_wetdry_adjustment=True,
                  force_constant_inlet_elevations=False,
                  description=None,
@@ -94 +95 @@
             raise Exception(msg)
         self.use_old_momentum_method = use_old_momentum_method
+        self.always_use_Q_wetdry_adjustment = always_use_Q_wetdry_adjustment
 
 
@@ -119 +121 @@
         self.accumulated_flow = 0.0
         self.discharge = 0.0
+        self.discharge_function_value = 0.0
         self.velocity = 0.0
         self.outlet_depth = 0.0
@@ -218 +221 @@
         old_inflow_ymom = self.inflow.get_average_ymom()
 
-        #semi_implicit = True
-        #if semi_implicit:   
-
-        # FIXME: This update replaces Q with Q*new_inflow_depth/old_inflow_depth
-        # This has good wetting and drying properties but means that
-        # discharge != Q.
-        # We should be able to check if old_inflow_depth*old_inflow_area > Q*dt,
-        # and in that case we don't need this implicit trick, and could
-        # have the discharge = Q (whereas in the nearly-dry case we want
-        # the trick).
-
         # Implement the update of flow over a timestep by
         # using a semi-implict update. This ensures that
@@ -237 +229 @@
             dt_Q_on_d = 0.0
 
-        # The depth update is: 
+        always_use_Q_wetdry_adjustment = self.always_use_Q_wetdry_adjustment
+        use_Q_wetdry_adjustment = ((always_use_Q_wetdry_adjustment) |\
+            (old_inflow_depth*self.inflow.get_area() <= Q*timestep))
+        # If we use the 'Q_adjustment', then the discharge is rescaled so that
+        # the depth update is: 
         #    new_inflow_depth*inflow_area = 
         #    old_inflow_depth*inflow_area - 
         #    timestep*Q*(new_inflow_depth/old_inflow_depth)
-        # The last term in () is a wet-dry improvement trick
-        # Note inflow_area is the area of all triangles in the inflow
-        # region -- so this is a volumetric balance equation
+        # The last term in () is a wet-dry improvement trick (which rescales Q)
+        #
+        # Before Feb 2015 this rescaling was always done (even if the flow was
+        # not near wet-dry). Now if use_old_Q_adjustment=False, the rescaling
+        # is only done if required to avoid drying
+        #
         #
         factor = 1.0/(1.0 + dt_Q_on_d/self.inflow.get_area())
-        new_inflow_depth = old_inflow_depth*factor
+
+
+        if use_Q_wetdry_adjustment:
+            # FIXME:
+            new_inflow_depth = old_inflow_depth*factor
+
+            if(old_inflow_depth>0.):
+                timestep_star = timestep*new_inflow_depth/old_inflow_depth
+            else:
+                timestep_star = 0.
+
+        else:
+            new_inflow_depth = old_inflow_depth - timestep*Q/self.inflow.get_area()
+            timestep_star = timestep
 
         if(self.use_old_momentum_method):
@@ -256 +268 @@
         else:
             # For the momentum balance, note that Q also advects the momentum,
-            # The volumetric momentum flux should be Q*momentum, where
+            # The volumetric momentum flux should be Q*momentum/depth, where
             # momentum has an average value of new_inflow_mom (or
-            # old_inflow_mom). For consistency we keep using the
-            # (new_inflow_depth/old_inflow_depth) factor for discharge:
+            # old_inflow_mom).  We use old_inflow_depth for depth
             #
             #     new_inflow_xmom*inflow_area = 
             #     old_inflow_xmom*inflow_area - 
-            #     [timestep*Q*(new_inflow_depth/old_inflow_depth)]*new_inflow_xmom
+            #     [timestep*Q]*new_inflow_xmom/old_inflow_depth
             # and:
             #     new_inflow_ymom*inflow_area = 
             #     old_inflow_ymom*inflow_area - 
-            #     [timestep*Q*(new_inflow_depth/old_inflow_depth)]*new_inflow_ymom
+            #     [timestep*Q]*new_inflow_ymom/old_inflow_depth
             #
-            # The choice of new_inflow_mom in the final term at the end might be
-            # replaced with old_inflow_mom
+            # The units balance: m^2/s*m^2 = m^2/s*m^2 - s*m^3/s*m^2/s *m^(-1)
             #
-            factor2 = 1.0/(1.0 + dt_Q_on_d*new_inflow_depth/self.inflow.get_area())
+            if old_inflow_depth > 0.:
+                if use_Q_wetdry_adjustment:
+                    # Replace dt*Q with dt*Q*new_inflow_depth/old_inflow_depth = dt_Q_on_d*new_inflow_depth
+                    factor2 = 1.0/(1.0 + dt_Q_on_d*new_inflow_depth/(old_inflow_depth*self.inflow.get_area()))
+                else:
+                    factor2 = 1.0/(1.0 + timestep*Q/(old_inflow_depth*self.inflow.get_area()))
+            else:
+                factor2 = 0.
+
             new_inflow_xmom = old_inflow_xmom*factor2
             new_inflow_ymom = old_inflow_ymom*factor2
@@ -289 +307 @@
 
         # set outflow
-        if old_inflow_depth > 0.0 :
-            timestep_star = timestep*new_inflow_depth/old_inflow_depth
-        else:
-            timestep_star = 0.0
-
         outflow_extra_depth = Q*timestep_star/self.outflow.get_area()
         outflow_direction = - self.outflow.outward_culvert_vector
@@ -306 +319 @@
         self.accumulated_flow += gain
         self.discharge  = Q*timestep_star/timestep 
+        self.discharge_function_value = Q
         self.velocity =   barrel_speed
         self.outlet_depth = outlet_depth
@@ -541 +555 @@
         
         message += 'Discharge [m^3/s]: %.2f\n' % self.discharge
+        message += 'Discharge_function_value [m^3/s]: %.2f\n' % self.discharge_function_value
         message += 'Velocity  [m/s]: %.2f\n' % self.velocity
         message += 'Outlet Depth  [m]: %.2f\n' % self.outlet_depth
@@ -562 +577 @@
             self.log_filename = self.domain.get_datadir() + '/' + self.label + '.log'
             log_to_file(self.log_filename, self.statistics(), mode='w')
-            log_to_file(self.log_filename, 'time, discharge, velocity, accumulated_flow, driving_energy, delta_total_energy')
+            log_to_file(self.log_filename, 'time, discharge, discharge_function_value, velocity, accumulated_flow, driving_energy, delta_total_energy')
 
             #log_to_file(self.log_filename, self.culvert_type)
@@ -571 +586 @@
         message  = '%.5f, ' % self.domain.get_time()
         message += '%.5f, ' % self.discharge
+        message += '%.5f, ' % self.discharge_function_value
         message += '%.5f, ' % self.velocity
         message += '%.5f, ' % self.accumulated_flow