Changeset 9659

Timestamp:

Feb 11, 2015, 11:26:16 AM (9 years ago)

Author:

davies

Message:

Experimenting with different smooths for internal boundary functions

Location:

trunk/anuga_core/anuga

Files:

• parallel/parallel_internal_boundary_operator.py (modified) (4 diffs)
• shallow_water/shallow_water_domain.py (modified) (3 diffs)
• structures/internal_boundary_operator.py (modified) (5 diffs)

trunk/anuga_core/anuga/parallel/parallel_internal_boundary_operator.py

@@ -100 +100 @@
         self.smoothing_timescale = 0.
         self.smooth_Q = 0.
+        self.smooth_delta_total_energy = 0.
         # Set them based on a call to the discharge routine with smoothing_timescale=0.
         # [values of self.smooth_* are required in discharge_routine, hence dummy values above]
@@ -106 +107 @@
         # Finally, set the smoothing timescale we actually want
         self.smoothing_timescale = smoothing_timescale
+        self.smooth_delta_total_energy = self.delta_total_energy
 
     def parallel_safe(self):
@@ -162 +164 @@
         # Determine the direction of the flow
         if self.myid == self.master_proc:
-            if self.use_velocity_head:
-                self.delta_total_energy = enq_total_energy0 - enq_total_energy1
-                self.driving_energy = max(enq_total_energy0, enq_total_energy1)
-                # Compute discharge
-                Q = self.internal_boundary_function(enq_total_energy0, enq_total_energy1)
-            else:
-                self.delta_total_energy = enq_stage0 - enq_stage1
-                self.driving_energy = max(enq_stage0, enq_stage1)
-                # Compute discharge
-                Q = self.internal_boundary_function(enq_stage0, enq_stage1)
-
-            # Other variables required by anuga's structure operator are not used
+            # Variables required by anuga's structure operator which are not
+            # used
             barrel_velocity = numpy.nan
             outlet_culvert_depth = numpy.nan
@@ -179 +171 @@
             case = ''
 
-            # Use time-smoothed discharge
+            # 'Timescale' for smoothed discharge and energy
             ts = self.domain.timestep/max(self.domain.timestep, self.smoothing_timescale, 1.0e-30)
+
+            # Energy or stage as head
+            if self.use_velocity_head:
+                E0 = enq_total_energy0
+                E1 = enq_total_energy1
+            else:
+                E0 = enq_stage0
+                E1 = enq_stage1
+
+            self.delta_total_energy = E0 - E1
+            self.driving_energy = max(E0, E1)
+
+            # Compute a 'smoothed' delta_total_energy for improved numerical
+            # stability
+            self.smooth_delta_total_energy = self.smooth_delta_total_energy +\
+                ts*(self.delta_total_energy - self.smooth_delta_total_energy)
+
+            if numpy.sign(self.smooth_delta_total_energy) != numpy.sign(self.delta_total_energy):
+                self.smooth_delta_total_energy = 0.
+
+            # Compute the 'tailwater' energy from the 'headwater' energy
+            # and the smooth_delta_total_energy Note if ts = 1 (no
+            # smoothing), then the raw inlet energies are used
+            if E0 >= E1:
+                inlet0_energy = 1.0*E0
+                inlet1_energy = inlet0_energy - self.smooth_delta_total_energy
+
+            else:
+                inlet1_energy = 1.0*E1
+                inlet0_energy = inlet1_energy + self.smooth_delta_total_energy
+
+            # Compute discharge
+            Q = self.internal_boundary_function(inlet0_energy, inlet1_energy)
             self.smooth_Q = self.smooth_Q + ts*(Q - self.smooth_Q)
-            
+
             if numpy.sign(self.smooth_Q) != numpy.sign(Q):
                 # The flow direction of the 'instantaneous Q' based on the

trunk/anuga_core/anuga/shallow_water/shallow_water_domain.py

@@ -1444 +1444 @@
             send(volume,0)
         
-        barrier()
+        #barrier()
     
         if myid == 0:
@@ -1483 +1483 @@
             send(flux_integral,0)
         
-        barrier()
+        #barrier()
     
         if myid == 0:
@@ -1518 +1518 @@
             send(flux_integral,0)
         
-        barrier()
+        #barrier()
     
         if myid == 0:

trunk/anuga_core/anuga/structures/internal_boundary_operator.py

@@ -100 +100 @@
         self.smoothing_timescale = 0.
         self.smooth_Q = 0.
+        self.smooth_delta_total_energy = 0.
         # Set them based on a call to the discharge routine with smoothing_timescale=0.
         # [values of self.smooth_* are required in discharge_routine, hence dummy values above]
@@ -106 +107 @@
         # Finally, set the smoothing timescale we actually want
         self.smoothing_timescale = smoothing_timescale
-
+        self.smooth_delta_total_energy = self.delta_total_energy
+    
+    ###########################################################################
     def discharge_routine(self):
         """Procedure to determine the inflow and outflow inlets.
@@ -131 +134 @@
             self.inlet0_energy = self.inlets[0].get_enquiry_stage()
             self.inlet1_energy = self.inlets[1].get_enquiry_stage()
-
-        # Compute discharge
-        Q = self.internal_boundary_function(self.inlet0_energy, self.inlet1_energy)
+        
+        # Store these variables for anuga's structure output
+        self.driving_energy = max(self.inlet0_energy, self.inlet1_energy)
+        self.delta_total_energy = self.inlet0_energy - self.inlet1_energy
 
         # Other variables required by anuga's structure operator are not used
@@ -141 +145 @@
         case = ''
 
+        # ts is used for smoothing discharge and delta_total_energy
+        if self.domain.timestep > 0.:
+            ts = self.domain.timestep/max(self.domain.timestep, self.smoothing_timescale, 1.0e-30)
+        else:
+            ts = 1.0
+
+        # Compute a 'smoothed' delta_total_energy
+        self.smooth_delta_total_energy = self.smooth_delta_total_energy +\
+            ts*(self.delta_total_energy - self.smooth_delta_total_energy)
+
+        if numpy.sign(self.smooth_delta_total_energy) != numpy.sign(self.delta_total_energy):
+            self.smooth_delta_total_energy = 0.
+
+        # Compute the 'tailwater' energy from the 'headwater' energy and the smooth_delta_total_energy
+        # Note if ts = 1 (no smoothing), then the raw inlet energies are used
+        if self.inlet0_energy >= self.inlet1_energy:
+            inlet0_energy = 1.0*self.inlet0_energy
+            inlet1_energy = inlet0_energy - self.smooth_delta_total_energy
+
+            # Compute discharge
+            Q = self.internal_boundary_function(inlet0_energy, inlet1_energy)
+        else:
+            inlet1_energy = 1.0*self.inlet1_energy
+            inlet0_energy = inlet1_energy + self.smooth_delta_total_energy
+
+            # Compute discharge
+            Q = self.internal_boundary_function(inlet0_energy, inlet1_energy)
+
         # Use time-smoothed discharge
-        ts = self.domain.timestep/max(self.domain.timestep, self.smoothing_timescale, 1.0e-30)
         self.smooth_Q = self.smooth_Q + ts*(Q - self.smooth_Q)
 
@@ -163 +194 @@
             Q = min( abs(self.smooth_Q), abs(Q) )
 
-        # Store these variables for anuga's structure output
-        self.driving_energy = max(self.inlet0_energy, self.inlet1_energy)
-        self.delta_total_energy = self.inlet0_energy - self.inlet1_energy
-
         return Q, barrel_velocity, outlet_culvert_depth