Changeset 8254

Timestamp:

Nov 22, 2011, 1:59:47 PM (13 years ago)

Author:

paul

Message:

Cleaned up sqpipe code

Location:

trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe

Files:

• parabolic_canal.py (modified) (1 diff)
• solver.h (deleted)
• sqpipe_comp_flux.c (added)
• sqpipe_default.h (added)
• sqpipe_domain.py (modified) (3 diffs)
• sqpipe_dual_fixed_width_a_d.h (deleted)
• sqpipe_dual_fixed_width_a_d_comp_flux.c (deleted)
• sqpipe_dual_fixed_width_a_d_domain.py (deleted)
• sqpipe_free_surface_fixed_width_a_d.h (deleted)
• sqpipe_free_surface_fixed_width_a_d_comp_flux.c (deleted)
• sqpipe_free_surface_fixed_width_a_d_domain.py (deleted)
• sqpipe_free_surface_fixed_width_w_d.h (deleted)
• sqpipe_free_surface_fixed_width_w_d_comp_flux.c (deleted)
• sqpipe_free_surface_fixed_width_w_d_domain.py (deleted)
• sqpipe_pressurised_fixed_width_a_d.h (deleted)
• sqpipe_pressurised_fixed_width_a_d_comp_flux.c (deleted)
• sqpipe_pressurised_fixed_width_a_d_domain.py (deleted)
• sqpipe_template.h (deleted)
• sqpipe_template_comp_flux.c (deleted)
• sqpipe_template_domain.py (deleted)
• sqpipe_test_domain.py (deleted)
• test_sqpipe_domain.py (modified) (1 diff)
• test_sqpipe_dual_fixed_width_a_d.py (deleted)
• test_sqpipe_free_surface_fixed_width_a_d.py (deleted)
• test_sqpipe_free_surface_fixed_width_w_d.py (deleted)
• test_sqpipe_pressurised_fixed_width_a_d.py (deleted)
• test_sqpipe_template.py (deleted)

trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/parabolic_canal.py

@@ -107 +107 @@
     ztplot, = pylab.plot(x, z+t)
 
-    plot1.set_ylim([-1,40])
+    plot1.set_ylim([-1,50])
     pylab.xlabel('Position')
     pylab.ylabel('Stage')

trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_domain.py

@@ -46 +46 @@
 
 from anuga_1d.base.generic_domain import Generic_domain
+from anuga_1d.sqpipe.sqpipe_forcing_terms import *
+from anuga_1d.sqpipe.sqpipe_boundary_conditions import *
 
 #Shallow water domain
-class Sqpipe_domain(Generic_domain):
-
-    def __init__(self, coordinates, conserved_quantities, forcing_terms = [], boundary = None, state = None, update_state_flag = True, bulk_modulus = 1.0, tagged_elements = None):
-
+class Domain(Generic_domain):
+
+    def __init__(self, coordinates, boundary = None, forcing_terms = [], state = None, update_state_flag = True, bulk_modulus = 1.0, tagged_elements = None):
+        conserved_quantities = ['area', 'discharge']
         evolved_quantities = ['area', 'discharge', 'elevation', 'height', 'velocity','width','top','stage']
         other_quantities = ['friction']
@@ -149 +151 @@
 
         # Import flux method and call it
-        from anuga_1d.sqpipe.sqpipe_dual_fixed_width_a_d_comp_flux import compute_fluxes_ext
+        from anuga_1d.sqpipe.sqpipe_comp_flux import compute_fluxes_ext
         self.flux_timestep = compute_fluxes_ext(timestep,self)        
 
@@ -171 +173 @@
 
     def distribute_to_vertices_and_edges(self):
-        # Should be implemented by the derived class
-        raise Exception("Not implemented")
+        # Shortcuts
+        h0 = self.h0
+        epsilon = self.epsilon
+
+        area      = self.quantities['area']
+        discharge = self.quantities['discharge']
+        bed       = self.quantities['elevation']
+        height    = self.quantities['height']
+        velocity  = self.quantities['velocity']
+        width     = self.quantities['width']
+        top       = self.quantities['top']
+        stage     = self.quantities['stage']
+
+        #Arrays   
+        a_C   = area.centroid_values
+        d_C   = discharge.centroid_values
+        z_C   = bed.centroid_values
+        h_C   = height.centroid_values
+        u_C   = velocity.centroid_values
+        b_C   = width.centroid_values
+        t_C   = top.centroid_values
+        w_C   = stage.centroid_values
+
+        #Calculate height (and fix negatives)better be non-negative!
+        h_C[:] = a_C/(b_C + h0/(b_C + h0)) # Do we need to protect against small b? Make a b0 instead of h0 or use epsilon?
+        w_C[:] = z_C + h_C
+        u_C[:]  = d_C/(h_C + h0/(h_C + h0))
+
+        for name in ['velocity', 'stage' ]:
+            Q = self.quantities[name]
+            if self.order == 1:
+                Q.extrapolate_first_order()
+            elif self.order == 2:
+                Q.extrapolate_second_order()
+            else:
+                raise 'Unknown order'
+
+        # These have been extrapolated
+        w_V = self.quantities['stage'].vertex_values
+        u_V = self.quantities['velocity'].vertex_values
+
+        # These are the given geometry and remain fixed
+        z_V  = bed.vertex_values
+        b_V  = width.vertex_values
+        t_V  = top.vertex_values
+
+        # These need to be updated
+        a_V  = area.vertex_values
+        d_V  = discharge.vertex_values
+        h_V  = height.vertex_values
+
+        h_V[:,:]    = w_V - z_V
+
+        # Fix up small heights
+        h_0 = numpy.where(h_V[:,0] < 0.0, 0.0, h_V[:,0])
+        h_1 = numpy.where(h_V[:,0] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,1])
+
+        h_V[:,0] = h_0
+        h_V[:,1] = h_1
+
+        h_0 = numpy.where(h_V[:,1] < 0.0, h_V[:,1]+h_V[:,0], h_V[:,0])
+        h_1 = numpy.where(h_V[:,1] < 0.0, 0.0, h_V[:,1])
+
+        h_V[:,0] = h_0
+        h_V[:,1] = h_1
+
+        # It may still be possible that h is small 
+        # If we set h to zero, we should also set u to 0
+        h_V[:,:] = numpy.where(h_V < epsilon, 0.0, h_V)
+        u_V[:,:] = numpy.where(h_V < epsilon, 0.0, u_V)
+
+        # Reconstruct conserved quantities and make everything consistent
+        # with new h values
+        w_V[:,:] = z_V + h_V
+        a_V[:,:] = h_V * b_V
+        d_V[:,:] = u_V * a_V
+
+
+        return
         
     def evolve(self, yieldstep = None, finaltime = None, duration = None,

trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/test_sqpipe_domain.py

@@ -1 +1 @@
 import anuga_1d.sqpipe.parabolic_canal
 import time
-import anuga_1d.sqpipe.sqpipe_test_domain as dom
+import anuga_1d.sqpipe.sqpipe_domain as dom